Класс Matrix во время компиляции

Задуманный как небольшой проект для тестирования различных функций C ++ 20, а также для того, чтобы узнать немного больше о матрицах и их использовании, я решил реализовать относительно простой матричный класс.

После его реализации я понял, что каждая его часть должна быть выполнена в контексте времени компиляции, и повторно реализовал ее (и немного очистил). Я предоставлю здесь весь код, а затем выделю некоторые части, в которых я не уверен, есть ли лучший способ. Рядом с ним есть очень простой векторный класс для использования с оператором умножения матриц.

#include <array>
#include <stdexcept>
#include <functional>

namespace ctm
{
    typedef std::size_t index_t;
    typedef long double real;

    template <index_t N>
    requires (N > 0)
    struct vector
    {
        std::array<real, N> values;

        constexpr vector(const std::array<real, N>& vals)
        {
            for (index_t i = 0; i < N; ++i)
            {
                values[i] = vals[i];
            }
        };

        constexpr real dot(const vector<N>& other) const
        {
            real sum = 0.L;

            for (index_t i = 0; i < N; ++i)
            {
                sum += values[i] * other.values[i];
            }

            return sum;
        };
    };

    template <index_t M, index_t N>
    requires (M > 0 && N > 0)
    struct matrix
    {
        std::array<std::array<real, N>, M> elements = {};

        constexpr matrix() = default;

        constexpr matrix(const std::array<real, M* N>& elems)
        {
            for (index_t i = 0; i < M; ++i)
            {
                for (index_t j = 0; j < N; ++j)
                {
                    elements[i][j] = elems[i * N + j];
                }
            }
        };

        constexpr matrix(const std::function<real(index_t, index_t)>& func)
        {
            for (index_t i = 1; i <= M; ++i)
            {
                for (index_t j = 1; j <= N; ++j)
                {
                    elements[i - 1][j - 1] = func(i, j);
                }
            }
        };

        // accessor functions

        // 0-based indexing
        constexpr std::array<real, N>& operator [](index_t index)
        {
            return elements[index];
        };

        // 0-based indexing
        constexpr const std::array<real, N>& operator [](index_t index) const
        {
            return elements[index];
        };

        // 1-based indexing
        constexpr real& operator ()(index_t index_m, index_t index_n)
        {
            if (index_m == 0 || index_n == 0)
            {
                throw std::out_of_range("compound access is 1-indexed");
            }
            else if (index_m > M || index_n > N)
            {
                throw std::out_of_range("compound access out of range");
            }

            return elements[index_m - 1][index_n - 1];
        };

        // 1-based indexing
        constexpr const real& operator ()(index_t index_m, index_t index_n) const
        {
            if (index_m == 0 || index_n == 0)
            {
                throw std::out_of_range("compound access is 1-indexed");
            }
            else if (index_m > M || index_n > N)
            {
                throw std::out_of_range("compound access out of range");
            }

            return elements[index_m - 1][index_n - 1];
        };

        // 1-based indexing
        constexpr vector<N> row(index_t index) const
        {
            if (index == 0)
            {
                throw std::out_of_range{ "row access is 1-indexed " };
            }
            else if (index > M)
            {
                throw std::out_of_range("row access out of range");
            }

            return vector<N>{ elements[index - 1] };
        };

        // 1-based indexing
        constexpr vector<M> column(index_t index) const
        {
            if (index == 0)
            {
                throw std::out_of_range{ "column access is 1-indexed " };
            }
            else if (index > N)
            {
                throw std::out_of_range("column access out of range");
            }

            std::array<real, M> column = {};

            for (index_t j = 0; j < M; ++j)
            {
                column[j] = elements[j][index - 1];
            }

            return vector<M>{ column };
        };

        constexpr matrix<N, M> transpose() const
        {
            matrix<N, M> result = {};

            for (index_t j = 0; j < N; ++j)
            {
                for (index_t i = 0; i < M; ++i)
                {
                    result.elements[j][i] = elements[i][j];
                }
            }

            return result;
        };

        consteval std::pair<index_t, index_t> size() const
        {
            return { M, N };
        };

        // arithmetic operators

        constexpr matrix& operator +=(const matrix& other)
        {
            for (index_t i = 0; i < M; ++i)
            {
                for (index_t j = 0; j < N; ++j)
                {
                    elements[i][j] += other.elements[i][j];
                }
            }

            return *this;
        };

        constexpr matrix& operator -=(const matrix& other)
        {
            for (index_t i = 0; i < M; ++i)
            {
                for (index_t j = 0; j < N; ++j)
                {
                    elements[i][j] -= other.elements[i][j];
                }
            }

            return *this;
        };

        constexpr matrix& operator *=(real scalar)
        {
            for (index_t i = 0; i < M; ++i)
            {
                for (index_t j = 0; j < N; ++j)
                {
                    elements[i][j] *= scalar;
                }
            }

            return *this;
        };

        constexpr matrix& operator /=(real scalar)
        {
            for (index_t i = 0; i < M; ++i)
            {
                for (index_t j = 0; j < N; ++j)
                {
                    elements[i][j] /= scalar;
                }
            }

            return *this;
        };

        constexpr friend matrix operator +(matrix left, const matrix& right)
        {
            return left += right;
        };

        constexpr friend matrix operator -(matrix left, const matrix& right)
        {
            return left -= right;
        };

        constexpr friend matrix operator *(matrix mat, real scalar)
        {
            return mat *= scalar;
        };

        constexpr friend matrix operator *(real scalar, matrix mat)
        {
            return mat *= scalar;
        };

        constexpr friend matrix operator /(matrix mat, real scalar)
        {
            return mat /= scalar;
        };

        template <index_t P>
        constexpr friend matrix<M, P> operator *(const matrix<M, N>& left, const matrix<N, P>& right)
        {
            matrix<M, P> result;

            for (index_t i = 1; i <= M; ++i)
            {
                for (index_t j = 1; j <= P; ++j)
                {
                    vector<N> left_row = left.row(i);
                    vector<N> right_column = right.column(j);
                    result(i, j) = left_row.dot(right_column);
                }
            }

            return result;
        };

        // comparison operator

        constexpr friend bool operator ==(const matrix&, const matrix&) = default;

        // the following functions are only valid for square matrices (M == N)

        // 1-based indexing
        template <index_t m = M, index_t n = N>
        requires (m == n && m > 1)
        constexpr matrix<m - 1, n - 1> first_minor(index_t i, index_t j) const
        {
            if (i == 0 || j == 0)
            {
                throw std::out_of_range("minor row and column indices are 1-indexed");
            }
            else if (i > M || j > N)
            {
                throw std::out_of_range("minor row and column indices out of range");
            }

            matrix<m - 1, n - 1> minor = {};

            index_t k = 0;
            index_t l = 0;
            for (index_t p = 0; p < m; ++p)
            {
                // skip over the row specified by i
                if (p == i - 1)
                {
                    continue;
                }

                for (index_t q = 0; q < n; ++q)
                {
                    // skip over the column specified by j
                    if (q == j - 1)
                    {
                        continue;
                    }

                    minor.elements[k][l] = elements[p][q];

                    ++l;
                }

                ++k;
                l = 0;
            }

            return minor;
        };

        constexpr real trace() const requires(M == N)
        {
            real sum = 0.L;

            for (index_t i = 0, j = 0; i < M; ++i, ++j)
            {
                sum += elements[i][j];
            }

            return sum;
        };

        constexpr real determinant() const requires(M == N && M == 2)
        {
            return elements[0][0] * elements[1][1] - elements[0][1] * elements[1][0];
        };

        constexpr real determinant() const requires(M == N && M > 2)
        {
            real sum = 0.L;
            real parity = 1.L;

            for (index_t j = 0; j < N; ++j)
            {
                auto submatrix = first_minor(1, j + 1);

                sum += elements[0][j] * parity * submatrix.determinant();
                parity = -parity;
            }

            return sum;
        };

        constexpr bool symmetric() const requires(M == N)
        {
            return *this == transpose();
        };

        constexpr bool skew_symmetric() const requires(M == N)
        {
            return *this == (transpose() *= -1.L);
        };

        // static square matrix creators

        static constexpr matrix<M, N> identity() requires(M == N)
        {
            matrix<M, N> id = {};

            for (index_t i = 0, j = 0; i < M; ++i, ++j)
            {
                id.elements[i][j] = 1.L;
            }

            return id;
        };

        static constexpr matrix<M, N> diagonal(const std::array<real, M>& elems) requires(M == N)
        {
            matrix<M, N> diag = {};

            for (index_t i = 0, j = 0; i < M; ++i, ++j)
            {
                diag.elements[i][j] = elems[i];
            }

            return diag;
        };
    };
};

Одна область, которая вызвала проблемы в обеих моих реализациях, — это конструктор, который принимает объект std :: function. Матрицы могут иметь свои элементы, определяемые результатом функции, которая принимает индексы элемента в качестве входных данных, и я хотел воспроизвести это. Это самый эффективный способ? Мне пришлось сделать это явным, потому что matrix<M, N> может быть преобразован в std::function<real(index_t, index_t)>, поскольку он определяет real operator ()(index_t, index_t) для доступа к 1-индексированной матрице.

Кроме того, есть ли способ потенциально очистить first_minorдекларация? Мне не нравится то, что я должен делать, чтобы эта одна работала, но если я попытаюсь следовать тому же шаблону, что и другие функции, requires выражение, он жалуется.

Общие советы также были бы признательны, спасибо!

Стилистически я бы предпочел увидеть requires-предложения (вроде noexcept-clauses и конечные возвращаемые типы) с отступом. То есть где у вас

template <class T>
requires foo<T>
auto foo()
noexcept(bar)
-> baz

Я бы предпочел увидеть

template<class T>
    requires foo<T>
auto foo()
    noexcept(bar)
    -> baz

просто чтобы лучше выделить важные части.

Ты используешь size_t для вашего типа индекса, что означает, что вам нужно в особом случае 0 во многих местах. При работе с «целыми числами» предпочитайте типы со знаком; Беззнаковые типы лучше зарезервировать для манипуляций с битовой маской. Видеть «The ‘unsigned для диапазона значений «антипаттерн» (13.03.2018). Где у вас есть

typedef std::size_t index_t;
[...]
    constexpr real& operator ()(index_t index_m, index_t index_n)
    {
        if (index_m == 0 || index_n == 0)
        {
            throw std::out_of_range("compound access is 1-indexed");
        }
        else if (index_m > M || index_n > N)
        {
            throw std::out_of_range("compound access out of range");
        }

Я бы предпочел увидеть

using index_t = int;  // actually you don't even need this typedef for anything
[...]
    constexpr real& operator()(int m, int n)
    {
        if (!(1 <= m && m <= M) || !(1 <= n && n <= N)) {
            throw std::out_of_range("compound access out of range");
        }

Если пользователь пишет mymatrix(-1, -1), вы не хотите, чтобы ваш код молча преобразовывал это в mymatrix(ULONG_MAX, ULONG_MAX) — это просто сбивает с толку!

Также обратите внимание на мои стилистические исправления: подтянутые скобы для экономии вертикальной недвижимости; идиоматические сравнения диапазонов (и, кстати, если вы думали, что мы могли бы использовать C ++ 20 std::in_range здесь вы ошибаетесь, он не делает то, что вы думаете). И записывать имена операторов одним словом без пробелов: operator() это имя этой функции, так же как и column это имя функции под ним.

    constexpr friend matrix operator +(matrix left, const matrix& right)
    {
        return left += right;
    };

Завершающая точка с запятой не нужна. Это не В самом деле имеет значение, потому что ваш matrix легко копируется и в любом случае не может извлечь выгоду из семантики перемещения; но обратите внимание, что вы просите здесь

• Сделайте объект с именем left
• Добавлять right в left, давая matrix& (это относится к left)
• Скопируйте матрицу, на которую ссылается это matrix& в слот возврата

В типичном, удобном для перемещения коде (например, если ваш matrix были std::string или же std::vector), вам лучше сказать

• Сделайте объект с именем left
• Добавлять right в left
• Переместить (или переместить-убрать) left сам в слот возврата

который будет написан

    friend constexpr matrix operator+(matrix a, const matrix& b)
    {
        a += b;
        return a;  // implicit move happens here
    }

Обратите внимание также на «порядок прилагательных»: так же, как в английском языке у нас может быть «милый маленький старый французский нож для строгания», в C ++ мы можем иметь

friend static inline constexpr virtual explicit operator T()

(friend предшествует обычному объявлению и, таким образом, «его легче всего удалить»; static и inline спецификаторы хранилища; constexpr еще не влияет на то, как система типов видит сущность; virtual влияет на это; и explicit по сути является ключевым словом, которое гласит: «Осторожно, здесь идет конструктор или оператор преобразования», синтаксически заменяющее тип возвращаемого значения нормальной функции.)

РЕДАКТИРОВАТЬ: это теперь сообщение в блоге. Видеть «static constexpr unsigned long это «прекрасный маленький старый французский ножик для резки» C ++ « (2021-04-03).

constexpr matrix(const std::array<real, M* N>& elems)
                                          ^

Это очень неудачное заблудшее пространство!

constexpr matrix(const std::function<real(index_t, index_t)>& func)

Поскольку вы не можете позвонить std::function во время компиляции нет смысла ставить constexpr на этом конструкторе. Но ты это знал. Вероятно, вместо этого вам следовало написать на C ++ 20 либо

template<std::invocable<int, int> F>
constexpr matrix(const F& func)
{

или же

template<class F>
    requires std::is_invocable_r<real, const F&, int, int>
constexpr matrix(const F& func)
{

Первый красивее. Последний более строгий: он проверяет не только F неприкосновенен, но это const F& вызывается (так как это то, что вы на самом деле будете вызывать в теле функции), а также то, что результат этого вызова может быть преобразован в real. Другой вариант:

template<class F>
constexpr matrix(const F& func)
    requires requires { { func(1,1) } -> std::convertible_to<real>; }
{

Но с точки зрения синтаксиса это выглядит гораздо неприятнее, и в конечном итоге проверяется то же самое, что и is_invocable_r версия.

В любом случае, смысл всего этого в том, чтобы избегать стирания типа. Вы хотите, чтобы этот конструктор принимал любой вызываемый объект? — итак делать он принимает любой вызываемый объект! Избавьтесь от посредников. Это также улучшит генерацию кода, поскольку оптимизатор, вероятно, сможет встроить все это целиком и даже развернуть эти циклы.

Тем не мение. Это, вероятно, столкнет вас с обычной проблемой с жадными шаблонами конструкторов, потому что угадайте, какой тип вы сделали вызываемым с подписью real(int, int)? Верно: matrix сам! Решение здесь не делай этого. Объедините обычный явный шаблон конструктора с обычным .at() функция-член (как и vector и map иметь), пока вы ждете мульти-операнда operator[] приземлиться на C ++ 2b.

template<std::invocable<int, int> F>
constexpr explicit matrix(const F& func) { ... }

constexpr real& at(int i, int j) { ... }
constexpr const real& at(int i, int j) const { ... }

Обратите внимание, что я пометил этот конструктор explicit, чтобы предотвратить неявные преобразования из std::function к ctm::matrix. Вы хотите, чтобы эти преобразования были возможный, но вы никогда не хотите, чтобы они произошли неявно! Как правило, каждый конструктор, который вы пишете, должен быть explicit, за исключением конструкторов копирования и перемещения (потому что вы хочу копии должны происходить неявно), и за исключением очень редких случаев, таких как std::stringконструктор из const char *.

Вы спросили о first_minor:

    template <index_t m = M, index_t n = N>
    requires (m == n && m > 1)
    constexpr matrix<m - 1, n - 1> first_minor(index_t i, index_t j) const

Поначалу я даже не могу понять, что вы здесь пытаетесь сделать. Эти параметры шаблона m и n фактические параметры, которые вы ожидаете от вызывающего абонента? или просто копии M и N что вы используете для метапрограммирования? Если вы добавляете метапрограммирование в свой список параметров шаблона в любом случае, то сначала вы должны попытаться сделать это так, как мы это делали с C ++ 11:

    template<class = std::enable_if_t<(M == N) && (M > 1)>
    constexpr matrix<M-1, N-1> first_minor(int i, int j) const

Однако это все равно не удается, когда либо M или же N равно 1, потому что в этот момент возвращаемый тип matrix<M-1, N-1> становится плохо сформированным (из-за того, что на matrix). Поэтому мы должны делать что-то хитрое, чтобы поддерживать его в правильном виде. Например:

    using SmallerMatrix = matrix<(M > 1) ? M-1 : 1, (N > 1) ? N-1 : 1>;

    constexpr SmallerMatrix first_minor(int i, int j) const
        requires (M == N && M > 1)

Более «умный» способ написания «C ++ 20-ish» был бы

    constexpr auto first_minor(int i, int j) const
        -> matrix<std::clamp(M-1, 1, M), std::clamp(N-1, 1, N)>
        requires (M == N && M > 1)

(Если вы уже прислушались к моему совету использовать int для index_t, то есть. В противном случае компилятор будет жаловаться, что std::clamp вывел конфликтующие типы для своего параметра типа шаблона. Все это как бы держится в одной системе. :))

Наконец, напишите несколько модульных тестов! Вы могли бы улучшить этот вопрос, сославшись на полный компилируемый пример на Godbolt. Так было бы проще поиграть с first_minor метапрограммирование, например.