matrix.h

Go to the documentation of this file.

// =====================================================================================================================
//  fennec, a free and open source game engine
//  Copyright © 2025  Medusa Slockbower
//
//  This program is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program.  If not, see <https://www.gnu.org/licenses/>.
// =====================================================================================================================
 
 
#ifndef FENNEC_MATH_MATRIX_H
#define FENNEC_MATH_MATRIX_H
 
 
#include <fennec/math/detail/_fwd.h>
#include <fennec/math/detail/_matrix.h>
 
#include <fennec/containers/array.h>
 
#include <fennec/math/geometric.h>
#include <fennec/math/vector_traits.h>
 
namespace fennec
{
 
 
template<typename scalar, size_t rows, size_t...cols>
constexpr vec<scalar, rows> column(const matrix<scalar, rows, cols...>& m, size_t i) noexcept {
    return m[i];
}
 
 
template<typename scalar, size_t rows, size_t...cols>
constexpr vec<scalar, sizeof...(cols)> row(const matrix<scalar, rows, cols...>& m, size_t i) noexcept {
    return vec<scalar, sizeof...(cols)>(m[cols][i]...);
}
 
 
template<typename ScalarT>
using tmat2x2 = mat<ScalarT, 2, 2>; 
template<typename ScalarT> using tmat2x3 = mat<ScalarT, 2, 3>; 
template<typename ScalarT> using tmat2x4 = mat<ScalarT, 2, 4>; 
template<typename ScalarT> using tmat3x2 = mat<ScalarT, 3, 2>; 
template<typename ScalarT> using tmat3x3 = mat<ScalarT, 3, 3>; 
template<typename ScalarT> using tmat3x4 = mat<ScalarT, 3, 4>; 
template<typename ScalarT> using tmat4x2 = mat<ScalarT, 4, 2>; 
template<typename ScalarT> using tmat4x3 = mat<ScalarT, 4, 3>; 
template<typename ScalarT> using tmat4x4 = mat<ScalarT, 4, 4>; 
 
using mat2 = tmat2x2<float_t>; 
using mat3 = tmat3x3<float_t>; 
using mat4 = tmat4x4<float_t>; 
 
using mat2x2 = tmat2x2<float_t>; 
using mat2x3 = tmat2x3<float_t>; 
using mat2x4 = tmat2x4<float_t>; 
using mat3x2 = tmat3x2<float_t>; 
using mat3x3 = tmat3x3<float_t>; 
using mat3x4 = tmat3x4<float_t>; 
using mat4x2 = tmat4x2<float_t>; 
using mat4x3 = tmat4x3<float_t>; 
using mat4x4 = tmat4x4<float_t>; 
 
using dmat2 = tmat2x2<double_t>; 
using dmat3 = tmat3x3<double_t>; 
using dmat4 = tmat3x3<double_t>; 
 
using dmat2x2 = tmat2x2<double_t>; 
using dmat2x3 = tmat2x3<double_t>; 
using dmat2x4 = tmat2x4<double_t>; 
using dmat3x2 = tmat3x2<double_t>; 
using dmat3x3 = tmat3x3<double_t>; 
using dmat3x4 = tmat3x4<double_t>; 
using dmat4x2 = tmat4x2<double_t>; 
using dmat4x3 = tmat4x3<double_t>; 
using dmat4x4 = tmat4x4<double_t>; 
 
 
template<typename scalar, size_t rows, size_t...cols>
constexpr matrix<scalar, rows, cols...> matrixCompMult(const matrix<scalar, rows, cols...>& x, const matrix<scalar, rows, cols...>& y) noexcept {
    return matrix<scalar, rows, cols...>(x[cols] * y[cols] ...);
}
 
template<typename scalar, size_t...s0, size_t...s1>
constexpr matrix<scalar, sizeof...(s0), s1...> outerProduct(const vector<scalar, s0...>& c, const vector<scalar, s1...>& r) noexcept {
    return matrix<scalar, sizeof...(s0), s1...>(
        c * r[s1]...
    );
}
 
template<typename scalar, size_t rows, size_t...cols>
constexpr mat<scalar, rows, sizeof...(cols)> transpose(const matrix<scalar, rows, cols...>& m) noexcept {
    return mat<scalar, rows, sizeof...(cols)>::transpose(m);
}
 
template<typename scalar, size_t rows, size_t...cols>
constexpr scalar determinant(const matrix<scalar, rows, cols...>&) noexcept {
    // ReSharper disable once CppStaticAssertFailure
    static_assert(false, "implementation undefined");
    return 0;
}
 
template<typename scalar, size_t rows, size_t...cols>
constexpr matrix<scalar, rows, cols...> inverse(const matrix<scalar, rows, cols...>&) noexcept {
    // ReSharper disable once CppStaticAssertFailure
    static_assert(false, "implementation undefined");
    return 1;
}
 
template<typename ScalarT, size_t RowsV, size_t...ColIndicesV> requires(is_scalar_v<ScalarT>)
struct matrix
{
// Assertions ==========================================================================================================
 
    static_assert(is_arithmetic_v<ScalarT> or is_bool_v<ScalarT>);
 
 
// Typedefs & Constants ================================================================================================
 
    static constexpr size_t rows = RowsV;
 
    static constexpr size_t columns = sizeof...(ColIndicesV);
 
    static constexpr size_t num_components = rows * columns;
 
    using row_t = vec<ScalarT, columns>;
 
    using column_t = vec<ScalarT, rows>;
 
    using scalar_t = ScalarT;
 
    using matrix_t = matrix;
 
    using transpose_t = mat<ScalarT, rows, columns>;
 
 
    array<column_t, columns> data;
 
 
 
// Constructors ========================================================================================================
 
 
 
    constexpr matrix() = default;
 
 
    constexpr matrix(const matrix_t& mat)
        : data{ mat.data } {
    }
 
    constexpr matrix(matrix_t&& mat) noexcept
        : data{ mat.data } {
    }
 
    template<typename OScalarT>
    constexpr matrix(const matrix<OScalarT, RowsV, ColIndicesV...>& mat)
        : matrix() {
        for (size_t i = 0; i < columns; ++i) {
            for (size_t j = 0; j < rows; ++j) {
                data[i][j] = scalar_t(mat[i][j]);
            }
        }
    }
 
    template<typename OScalarT>
    constexpr matrix(matrix<OScalarT, RowsV, ColIndicesV...>&& mat) noexcept
        : matrix() {
        for (size_t i = 0; i < columns; ++i) {
            for (size_t j = 0; j < rows; ++j) {
                data[i][j] = scalar_t(mat[i][j]);
            }
        }
    }
 
 
    constexpr matrix(scalar_t s) {
        (((ColIndicesV < rows) ? data[ColIndicesV][ColIndicesV] = s : 0), ...);
    }
 
 
    template<typename...ArgsT> requires(total_component_count_v<ArgsT...> == num_components)
    constexpr matrix(ArgsT&&...args) {
        matrix::_construct(fennec::forward<ArgsT>(args)...);
    }
 
 
// Assignment Operators ================================================================================================
 
 
    constexpr matrix_t& operator=(const matrix_t& rhs) {
        data = rhs.data; return *this;
    }
 
    constexpr matrix_t& operator=(matrix_t&& rhs) noexcept {
        data = rhs.data; return *this;
    }
 
 
 
// Access Operators ====================================================================================================
 
 
    constexpr column_t& operator[](size_t i) {
        return data[i];
    }
 
    constexpr const column_t& operator[](size_t i) const {
        return data[i];
    }
 
    constexpr scalar_t& operator[](size_t i, size_t j) {
        return data[i][j];
    }
 
    constexpr scalar_t operator[](size_t i, size_t j) const {
        return data[i][j];
    }
 
 
// Scalar Operators ====================================================================================================
 
 
    constexpr friend matrix_t operator+(const matrix_t& lhs, scalar_t rhs) {
        return matrix_t(lhs[ColIndicesV] + rhs ...);
    }
 
    constexpr friend matrix_t operator+(scalar_t lhs, const matrix_t& rhs) {
        return matrix_t(lhs + rhs[ColIndicesV] ...);
    }
 
    constexpr friend matrix_t operator-(const matrix_t& lhs, scalar_t rhs) {
        return matrix_t(lhs[ColIndicesV] - rhs ...);
    }
 
    constexpr friend matrix_t operator-(scalar_t lhs, const matrix_t& rhs) {
        return matrix_t(lhs - rhs[ColIndicesV] ...);
    }
 
    constexpr friend matrix_t operator*(const matrix_t& lhs, scalar_t rhs) {
        return matrix_t(lhs[ColIndicesV] * rhs ...);
    }
 
    constexpr friend matrix_t operator*(scalar_t lhs, const matrix_t& rhs) {
        return matrix_t(lhs * rhs[ColIndicesV] ...);
    }
 
    constexpr friend matrix_t operator/(const matrix_t& lhs, scalar_t rhs) {
        return matrix_t(lhs[ColIndicesV] / rhs ...);
    }
 
    constexpr friend matrix_t operator/(scalar_t lhs, const matrix_t& rhs) {
        return matrix_t(lhs / rhs[ColIndicesV] ...);
    }
 
 
 
    constexpr friend matrix_t& operator+=(matrix_t& lhs, scalar_t rhs) {
        ((lhs[ColIndicesV] += rhs), ...);
        return lhs;
    }
 
    constexpr friend matrix_t& operator+=(scalar_t lhs, matrix_t& rhs) {
        ((rhs[ColIndicesV] += lhs), ...);
        return rhs;
    }
 
    constexpr friend matrix_t& operator-=(matrix_t& lhs, scalar_t rhs) {
        ((lhs[ColIndicesV] -= rhs), ...);
        return lhs;
    }
 
    constexpr friend matrix_t& operator-=(scalar_t lhs, matrix_t& rhs) {
        ((rhs[ColIndicesV] -= lhs), ...);
        return rhs;
    }
 
    constexpr friend matrix_t& operator*=(matrix_t& lhs, scalar_t rhs) {
        ((lhs[ColIndicesV] *= rhs), ...);
        return lhs;
    }
 
    constexpr friend matrix_t& operator*=(scalar_t lhs, matrix_t& rhs) {
        ((rhs[ColIndicesV] *= lhs), ...);
        return rhs;
    }
 
    constexpr friend matrix_t& operator/=(matrix_t& lhs, scalar_t rhs) {
        ((lhs[ColIndicesV] /= rhs), ...);
        return lhs;
    }
 
    constexpr friend matrix_t& operator/=(scalar_t lhs, matrix_t& rhs) {
        ((rhs[ColIndicesV] /= lhs), ...);
        return rhs;
    }
 
 
 
    constexpr friend matrix_t operator%(const matrix_t& lhs, scalar_t rhs) requires is_integral_v<scalar_t> {
        matrix_t retval = lhs;
        ((retval[ColIndicesV] %= rhs), ...);
        return retval;
    }
 
    constexpr friend matrix_t operator%(scalar_t lhs, const matrix_t& rhs) requires is_integral_v<scalar_t> {
        matrix_t retval = rhs;
        ((retval[ColIndicesV] %= lhs), ...);
        return retval;
    }
 
    constexpr friend matrix_t& operator%=(matrix_t& lhs, scalar_t rhs) requires is_integral_v<scalar_t> {
        ((lhs[ColIndicesV] %= rhs), ...);
        return lhs;
    }
 
    constexpr friend matrix_t& operator%=(scalar_t lhs, matrix_t& rhs) requires is_integral_v<scalar_t> {
        ((rhs[ColIndicesV] %= lhs), ...);
        return rhs;
    }
 
 
// Vector Operators ====================================================================================================
 
 
    constexpr friend column_t operator*(const matrix_t& lhs, const row_t& rhs) {
        return _mul(lhs, rhs);
    }
 
    constexpr friend row_t operator*(const column_t& lhs, const matrix_t& rhs) {
        return row_t(fennec::dot(fennec::column(rhs, ColIndicesV), lhs) ...);
    }
 
 
 
 
// Matrix Operators ====================================================================================================
 
 
    constexpr friend bool operator==(const matrix_t& lhs, const matrix_t& rhs) {
        return lhs.data == rhs.data;
    }
 
    constexpr friend bool operator!=(const matrix_t& lhs, const matrix_t& rhs) {
        return lhs.data != rhs.data;
    }
 
    template<size_t ORowsV, size_t...OColIndicesV> requires(columns == ORowsV)
    constexpr friend matrix<scalar_t, RowsV, OColIndicesV...> operator*(const matrix_t& lhs, const matrix<scalar_t, ORowsV, OColIndicesV...>& rhs) {
        return matrix<scalar_t, RowsV, OColIndicesV...>(
            matrix::_mul(lhs, rhs[OColIndicesV])...
        );
    }
 
    template<size_t ORowsV, size_t...OColIndicesV> requires(columns == ORowsV)
    constexpr matrix<scalar_t, RowsV, OColIndicesV...>& operator*=(const matrix<scalar_t, ORowsV, OColIndicesV...>& rhs) {
        return *this = *this * rhs;
    }
 
 
 
 
// Helpers =============================================================================================================
 
public:
 
    static constexpr matrix_t transpose(const transpose_t& mat) {
        return matrix_t(fennec::row(mat, ColIndicesV)...);
    }
 
private:
 
    // ReSharper disable once CppMemberFunctionMayBeStatic
    template<size_t i0 = 0>
    constexpr void _construct() {
        // base case, does nothing, this will get optimized away
    }
 
    // helper for parsing parameter packs
    template<size_t i0 = 0, typename HeadT, typename...RestT>
    constexpr void _construct(HeadT&& head, RestT&&...rest) {
        matrix::_insert<i0>(head); // insert the head element
        matrix::_construct<i0 + component_count_v<HeadT>>(std::forward<RestT>(rest)...); // propagate to the rest of the arguments
    }
 
    // helper for inserting a scalar value
    template<size_t i0 = 0>
    constexpr void _insert(scalar_t s) {
        data[i0 / rows][i0 % rows] = s;
    }
 
    // helper for inserting a scalar value of differing type
    template<size_t i0 = 0, typename OScalarT> requires(is_arithmetic_v<OScalarT>)
    constexpr void _insert(OScalarT s) {
        data[i0 / rows][i0 % rows] = scalar_t(s);
    }
 
    // helper for inserting a vector
    template<size_t i0 = 0, size_t...i>
    constexpr void _insert(const vector<scalar_t, i...>& v) {
        (matrix::_insert<i0 + i>(v[i]), ...);
    }
 
    // helper for inserting a vector of differing type
    template<size_t i0 = 0, typename OScalarT, size_t...i>
    constexpr void _insert(const vector<OScalarT, i...>& v) {
        (matrix::_insert<i0 + i>(v[i]), ...);
    }
 
    // helper for a linear algebraic multiply
    static constexpr column_t _mul(const matrix_t& lhs, const row_t& rhs) {
        // the compiler will optimize this better than writing out a specific definition
        // when compared to glm or CxxSwizzle, this is faster by a significant margin
        // all implementations provide 7 decimal places of precision
        return ((lhs[ColIndicesV] * rhs[ColIndicesV]) + ...);
    }
};
 
 
// Internal ============================================================================================================
 
}
 
 
#endif // FENNEC_MATH_MATRIX_H