Compile-time Matrix Class

Question

Intended as a small project to test out various C++20 features, as well as learn a little bit more about matrices and their uses, I decided to implement a relatively simple matrix class.

After implementing it, I figured out that every single part of it should be able to be done in a compile time context, and re-implemented it (and cleaned it up a little). I'll provide the whole code here, and then highlight some parts where I'm unsure of whether there's a better way. There's a very barebones vector class alongside it for use with the matrix multiplication operator.

#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;
        };
    };
};

One area that has caused trouble in both of my implementations is the constructor that takes a std::function object. Matrices can have their elements defined by the result of a function that takes the indices of the element as an input, and I wanted to replicate that. Is this the most effective way? I had to make it explicit because a matrix<M, N> is able to be converted into a std::function<real(index_t, index_t)>, since it defines real operator ()(index_t, index_t) for doing 1-indexed matrix access.

Also, is there any way to potentially clean up the first_minor's declaration? I dislike what I have to do to make that one work, but if I try to follow the same pattern as the other functions with a requires expression, it complains.

General advice would also be appreciated, thanks!

Answer 1

Stylistically, I'd prefer to see requires-clauses (like noexcept-clauses and trailing return types) indented. That is, where you have

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

I'd prefer to see

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

just to highlight the important parts better.

---

You use size_t for your index type, which means you have to special-case 0 in a lot of places. Prefer signed types when dealing with "integers"; unsigned types are better reserved for bitmask manipulations. See "The 'unsigned for value range' antipattern" (2018-03-13). Where you have

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");
        }

I'd rather see

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");
        }

If the user writes mymatrix(-1, -1), you don't want your code silently converting that into mymatrix(ULONG_MAX, ULONG_MAX) — that's just confusing!

Also, notice my stylistic cleanups: cuddled braces to save vertical real estate; idiomatic range-comparisons (and btw if you thought we could maybe use C++20 std::in_range here, you'd be wrong, it doesn't do what you'd think it does). And writing the names of operators as a single word without whitespace: operator() is the name of this function, just as much as column is the name of the function below it.

---

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

The trailing semicolon is unnecessary. It doesn't really matter because your matrix is trivially copyable and can't benefit from move semantics anyway; but notice that what you're asking for here is

• Make an object named left
• Add right into left, yielding a matrix& (that happens to refer to left)
• Copy the matrix referred to by that matrix& into the return slot

In typical move-friendly code (like, if your matrix were a std::string or std::vector), you'd be much better off saying

• Make an object named left
• Add right into left
• Move (or move-elide) left itself into the return slot

which would be spelled

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

Notice also the "adjective order": just as in English we can have a "lovely little old French whittling knife," in C++ we can have a

friend static inline constexpr virtual explicit operator T()

(friend precedes an otherwise ordinary declaration and is thus the "easiest to strip off"; static and inline are storage specifiers; constexpr doesn't yet affect how the type system sees the entity; virtual affects it; and explicit is basically the keyword that says "look out, here comes a constructor or conversion operator," going syntactically in place of a normal function's return type.)

EDIT: This is now a blog post. See "static constexpr unsigned long is C++'s 'lovely little old French whittling knife'" (2021-04-03).

---

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

That's a very unfortunate stray space!

---

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

Since you can't call a std::function at compile time, there's no point to putting constexpr on this constructor. But you knew that. Probably what you should have written instead — in C++20 — is either

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

or

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

The former is prettier. The latter is more rigorous: it checks not only that F is invocable, but that const F& is invocable (since that's what you'll actually be invoking within the body of the function) and also that the result of that invocation is convertible to real. Another option is:

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

But that's much ickier-looking syntactically, and ends up checking basically the same thing as the is_invocable_r version.

Anyway, the point of all of these is to avoid type erasure. You want this constructor to take any callable object?— so make it take any callable object! Cut out the middleman. This will also produce better codegen, since the optimizer can probably inline the whole thing and even unroll those loops.

However. This probably runs you smack into the usual problem with greedy constructor templates, because guess what type you've made invocable with the signature real(int, int)? That's right: matrix itself! The solution here is don't do that. Combine a normal explicit constructor template with a normal .at() member function (just like vector and map have), while you wait for multi-operand operator[] to land in 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 { ... }

Notice that I've marked this constructor explicit, in order to prevent implicit conversions from std::function to ctm::matrix. You want those conversions to be possible, but you never want them to happen implicitly! As a general rule, every constructor you write should be explicit, except for the copy and move constructors (because you want copies to happen implicitly), and except for extremely rare cases such as std::string's constructor from const char *.

---

You asked about 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

At first, I can't even really tell what you're trying to do here. Are those template parameters m and n actual parameters you expect the caller to pass in? or just copies of M and N that you're using for metaprogramming? If you're putting metaprogramming into your template parameter list anyway, then you should first try to do it in the way we've done it since 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

However, this still fails when either M or N is 1, because at that point the return type matrix<M-1, N-1> becomes ill-formed (thanks to the requires-clause on matrix). So we have to do something sneaky to keep it well-formed. For example:

    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)

A more "C++20-ish" clever way to write that would be

    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)

(As long as you have already taken my advice to use int for index_t, that is. Otherwise the compiler will complain that std::clamp has deduced conflicting types for its template type parameter. All this stuff kind of hangs together as one system. :))

---

Finally, write some unit tests! You could have improved this question by linking to a complete compilable example on Godbolt. That would have made it easier to play around with the first_minor metaprogramming, for example.