// g++ -std=c++11 -Wall -Wextra -Wpedantic -o matrix_ops matrix_ops.cpp && ./matrix_ops


/// Minimum type
template<typename T, int N, int M>
struct Mat
{
    T values[N][M];
    T       * operator[](int i)       { return values[i]; }
    T const * operator[](int i) const { return values[i]; }
};


/// Functions
// https://en.wikipedia.org/wiki/Minor_(linear_algebra)
// https://en.wikipedia.org/wiki/Laplace_expansion
// https://en.wikipedia.org/wiki/Adjugate_matrix
template<typename T, int N>
T cofactor(Mat<T, N, N> const & mat, int i, int j)
{
    auto sub = Mat<T, N-1, N-1>{};
    for (int n = 0; n < N-1; ++n)
    for (int m = 0; m < N-1; ++m)
        sub[n][m] = mat[n+(n>=i)][m+(m>=j)];
    return determinant(sub) * ((i+j)%2 ? -1 : 1);
}
template<typename T>
T cofactor(Mat<T, 1, 1> const &, int, int)
{
    return T(1);
}
template<typename T, int N>
T determinant(Mat<T, N, N> const & mat)
{
    auto det = T(0);
    for (int i = 0; i < N; ++i)
        det += mat[i][0] * cofactor(mat, i, 0);
    return det;
}
template<typename T, int N>
Mat<T, N, N> inverse(Mat<T, N, N> const & mat)
{
    auto adj = Mat<T, N, N>{};
    auto det = T(0);
    for (int i = 0; i < N; ++i)
        det += mat[i][0] * (adj[0][i] = cofactor(mat, i, 0));
    for (int j = 0; j < N; ++j)
    for (int i = 0; i < N; ++i)
        adj[j][i] = T(1) / det * (j ? cofactor(mat, i, j) : adj[0][i]);
    return adj;
}


/// Optional specializations
// Helps also for bigger `N` since those cases recursively reduce to these.
template<typename T>
T determinant(Mat<T, 1, 1> const & mat)
{
    return
        + mat[0][0];
}
template<typename T>
T determinant(Mat<T, 2, 2> const & mat)
{
    return
        + mat[0][0] * mat[1][1]
        - mat[0][1] * mat[1][0];
}
template<typename T>
T determinant(Mat<T, 3, 3> const & mat)
{
    return
        + mat[0][0] * (
            + mat[1][1] * mat[2][2]
            - mat[1][2] * mat[2][1]
        )
        - mat[0][1] * (
            + mat[1][2] * mat[2][0]
            - mat[1][0] * mat[2][2]
        )
        + mat[0][2] * (
            + mat[1][0] * mat[2][1]
            - mat[1][1] * mat[2][0]
        );
}


/// Tests
#include <iomanip>
#include <iostream>
#include <string>
int main()
{
    constexpr auto N = 4;
    constexpr auto P = 3;
    std::cout << std::showpos << std::fixed << std::setprecision(P);
    auto m1 = Mat<float, N, N>{{
        {2, 0, 0, 1},
        {0, 3, 0, 1},
        {0, 0, 4, 1},
        {0, 0, 0, 1},
    }};
    std::cout << determinant(m1) << "\n\n";
    // +24.000
    auto m2    = inverse(m1);
    auto none = std::string(3+P, ' ');
    for (int i = 0; i < N; ++i, std::cout << "\n")
    for (int j = 0; j < N; ++j, std::cout << (j == N ? "" : " "))
        m2[i][j]
            ? std::cout << m2[i][j]
            : std::cout << none;
    // +0.500               -0.500
    //        +0.333        -0.333
    //               +0.250 -0.250
    //                      +1.000
}