Orange/include/Orange/Math/Matrix.h

#pragma once

#include <Orange/Math/Vector.h>
#include <Orange/Math/Angle.h>

namespace orange
{
    template <typename T, size_t Rows, size_t Columns>
    struct Matrix
    {
        using RowVector = Vec<T, Columns>;

        constexpr Matrix(T scale = T{ 1 })
        {
            for (size_t i = 0; i < Rows; i++)
            {
                RowVector vector{};
                vector[i] = scale;
                data[i] = vector;
            }
        }

        constexpr Matrix(const RowVector& scale)
        {
            for (size_t i = 0; i < Rows; i++)
            {
                RowVector vector{};
                vector[i] = scale[i];
                data[i] = vector;
            }
        }

        template <typename... Args>
        constexpr Matrix(const Args&... args)
            : data {{ args... }}
        {
            static_assert(sizeof...(Args) == Rows);
        }

        constexpr Matrix(const T components[Columns][Rows])
        {
            std::copy(&components[0], &components[Columns][Rows], data.begin());
        }

        constexpr Matrix(const Matrix& other) = default;


        constexpr       RowVector& operator[](size_t index)       { return data[index]; }
        constexpr const RowVector& operator[](size_t index) const { return data[index]; }


        constexpr const RowVector* begin() const { return data.begin(); }
        constexpr       RowVector* begin()       { return data.begin(); }
        constexpr const RowVector* end()   const { return data.end();   }
        constexpr       RowVector* end()         { return data.end();   }


        constexpr bool operator==(const Matrix& other) const
        {
            return Equal(begin(), end(), other.begin());
        }

        constexpr bool operator!=(const Matrix& other) const
        {
            return !operator==(other);
        }


        template <typename UnaryOperation>
        constexpr Matrix TransformResult(UnaryOperation op) const
        {
            return TransformResult<Matrix>(begin(), end(), op);
        }
        
        template <typename BinaryOperation>
        constexpr Matrix TransformResult(const RowVector *other, BinaryOperation op) const
        {
            return TransformResult<Matrix>(begin(), end(), other, op);
        }

        template <typename BinaryOperation>
        constexpr Matrix TransformResult(const Matrix& other, BinaryOperation op) const
        {
            return TransformResult(other.begin(), op);
        }

        template <typename UnaryOperation>
        constexpr Matrix& TransformInPlace(UnaryOperation op)
        {
            Transform(begin(), end(), begin(), op);
            return *this;
        }
        
        template <typename BinaryOperation>
        constexpr Matrix& TransformInPlace(const RowVector *other, BinaryOperation op)
        {
            Transform(begin(), end(), other, begin(), op);
            return *this;
        }

        template <typename BinaryOperation>
        constexpr Matrix& TransformInPlace(const Matrix& other, BinaryOperation op)
        {
            return TransformInPlace(other.begin(), op);
        }

        // Simple math operations

        constexpr Matrix operator+(const Matrix& other) const
        {
            return TransformResult(other, Math::Add);
        }

        constexpr Matrix operator-(const Matrix& other) const
        {
            return TransformResult(other, Math::Subtract);
        }

        constexpr Matrix operator*(const T& scalar) const
        {
            return TransformResult([scalar](const RowVector& value) { return value * scalar; });
        }

        constexpr Matrix operator/(const T& scalar) const
        {
            return TransformResult([scalar](const RowVector&  value) { return value / scalar; });
        }

        constexpr Matrix operator%(const T& scalar) const
        {
            return TransformResult([scalar](const RowVector&  value) { return value % scalar; });
        }


        constexpr Matrix& operator+=(const Matrix& other)
        {
            return TransformInPlace(Math::Add);
        }

        constexpr Matrix& operator-=(const Matrix& other)
        {
            return TransformInPlace(Math::Subtract);
        }

        constexpr Matrix& operator*=(const T& scalar)
        {
            return TransformInPlace([scalar](const RowVector&  value) { return value * scalar; });
        }

        constexpr Matrix& operator/=(const T& scalar)
        {
            return TransformInPlace([scalar](const RowVector&  value) { return value / scalar; });
        }

        constexpr Matrix& operator%=(const T& scalar)
        {
            return TransformInPlace([scalar](const RowVector&  value) { return value % scalar; });
        }

        // Real matrix operations
        // TODO: Handle mixing-matching column counts
        constexpr Matrix operator*(const Matrix& other) const
        {
            const Matrix &us = *this;

            Matrix out{0.0f};
            for (size_t i = 0; i < Rows; i++)
            {
                for (size_t j = 0; j < Columns; j++)
                {
                    for (size_t k = 0; k < Columns; k++)
                        out[i][j] += other[i][k] * us[k][j];
                }
            }
            return out;
        }

        constexpr Matrix& operator*=(const Matrix& other) const
        {
            return (*this = *this * other);
        }

        constexpr RowVector operator*(const RowVector& v) const
        {
            auto mul = TransformResult(v.begin(), Math::Multiply);
            return Accumulate(mul.begin(), mul.end(), RowVector{});
        }

        Array<RowVector, Rows> data;
    };

    template <typename T, size_t Rows, size_t Columns>
    constexpr Matrix<T, Rows, Columns> operator*(T scalar, const Matrix<T, Rows, Columns>& matrix)
    {
        using J = Matrix<T, Rows, Columns>::RowVector;
        return matrix.TransformResult([scalar](J value) { return scalar * value; });
    }

    template <typename T, size_t Rows, size_t Columns>
    constexpr Matrix<T, Rows, Columns> operator/(T scalar, const Matrix<T, Rows, Columns>& matrix)
    {
        using J = Matrix<T, Rows, Columns>::RowVector;
        return matrix.TransformResult([scalar](J value) { return scalar / value; });
    }

    template <typename T, size_t Rows, size_t Columns>
    constexpr Matrix<T, Rows, Columns> operator%(T scalar, const Matrix<T, Rows, Columns>& matrix)
    {
        using J = Matrix<T, Rows, Columns>::RowVector;
        return matrix.TransformResult([scalar](J value) { return scalar % value; });
    }


    template <typename T, size_t Rows, size_t Columns>
    constexpr Matrix<T, Rows - 1, Columns - 1> minor(const Matrix<T, Rows, Columns>& a, size_t column, size_t row)
    {
        Matrix<T, Rows - 1, Columns - 1> mtx;
        for (size_t y = 0, my = 0; y < Rows; y++) {
            if (y == row) continue;
            for (size_t x = 0, mx = 0; x < Columns; x++) {
                if (x == column) continue;
                mtx[my][mx] = a[y][x];
                mx++;
            }
            my++;
        }
        return mtx;
    }

    template <typename T, size_t Rows, size_t Columns>
    constexpr T determinant(const Matrix<T, Rows, Columns>& a)
    {
        static_assert(Rows == Columns);

        if constexpr (Rows == 1) {
            return a[0][0];
        } else {
            T result = T{};
            for (size_t x = 0; x < Columns; x++) {
                const T sign = (x % 2) ? T{ -1 } : T{ 1 };
                result += sign * a[0][x] * determinant(minor(a, x, 0));
            }
            return result;
        }
    }

    template <typename T, size_t Rows, size_t Columns>
    constexpr Matrix<T, Columns, Rows> transpose(const Matrix<T, Rows, Columns>& a)
    {
        Matrix<T, Columns, Rows> result;
        for (size_t y = 0; y < Rows; y++) {
            for (size_t x = 0; x < Columns; x++)
                result[x][y] = a[y][x];
        }
        return result;
    }

    template <typename T, size_t Rows, size_t Columns>
    constexpr Matrix<T, Rows, Columns> hadamard(const Matrix<T, Rows, Columns>& x, const Matrix<T, Rows, Columns>& y)
    {
        return x.TransformResult(y, Math::Multiply);
    }

    using mat4 = Matrix<float, 4, 4>;

}