Constexpr determinant (2 dimensional std::array)

Question

I need to write a constexpr function that computes a determinant at compile time. The most obvious solution is to use Laplace expansion. C++14 is supported.

#include <array>
#include <utility>

constexpr int get_cofactor_coef(int i, int j) {
    return (i + j) % 2 == 0 ? 1 : -1;
}

template <int N>
constexpr int determinant(const std::array<std::array<int, N>, N>& a) {
    int det = 0;

    for (size_t i = 0u; i < N; ++i) {
        det += get_cofactor_coef(i, 1) * a[i][0] * determinant<N-1>(GET_SUBMATRIX_OF_A<N-1, I, J>(a);
    }

    return det;
}

template <>
constexpr int determinant<2>(const std::array<std::array<int, 2>, 2>& a) {
    return a[0][0] * a[1][1] - a[0][1] * a[1][0];
}

template <>
constexpr int determinant<1>(const std::array<std::array<int, 1>, 1>& a) {
    return a[0][0];
}

The problem is that I have absolutely no clue how to write the GET_SUBMATRIX_OF_A.

I know that I need to:

1. Generate a sequence (using std::integer_sequence probably);
2. Exclude from this sequence the i-th row;
3. Copy all but the first (0-th) column;

My constexpr skills are next to non-existent. Head on attempts to pass a to another function result in weird errors like error: '* & a' is not a constant expression.

All help is greatly appreciated!

Answer 1

The problem is that the non-const std::array<T, N>::operator[] (returning T&) is not constexpr until C++17, making it difficult to set the elements of the minor.

However, there is an escape hatch, which is that std::get<I>(std::array&) is constexpr, and it is perfectly legitimate to perform pointer arithmetic on the result, so we can rewrite

a[i]  // constexpr since C++17

as

(&std::get<0>(a))[i]  // constexpr in C++14!!

That is, we use std::get to obtain a constexpr reference to the first member of the array, take a pointer to it, and use the built-in [] operator on the pointer and index.

Then a two-level array member access a[i][j] becomes the horrendously ugly but still constexpr (&std::get<0>((&std::get<0>(a))[i]))[j], meaning we can write get_submatrix_of_a as an ordinary constexpr function:

template<std::size_t N>
constexpr std::array<std::array<int, N - 1>, N - 1>
get_submatrix_of_a(const std::array<std::array<int, N>, N>& a, int i, int j) {
    std::array<std::array<int, N - 1>, N - 1> r{};
    for (int ii = 0; ii != N - 1; ++ii)
        for (int jj = 0; jj != N - 1; ++jj)
            (&std::get<0>(((&std::get<0>(r))[ii])))[jj] = a[ii + (ii >= i ? 1 : 0)][jj + (jj >= j ? 1 : 0)];
    return r;
}

Remember that const std::array<T, N>::operator[] is already constexpr in C++14, so we don't need to rewrite the RHS of the minor construction.