more of the new matrix vector implementation

src/core/include/pw/core/math.hpp

@@ -25,6 +25,8 @@
 
 #include <cmath>
 #include <algorithm>
+#include <numbers>
+
 //#include <pw/core/matrix.hpp>
 //#include <pw/core/quaternion.hpp>
 

src/core/include/pw/core/matrix.hpp

@@ -26,9 +26,193 @@
 #include <pw/core/globals.hpp>
 #include <pw/core/math.hpp>
 #include <pw/core/vector.hpp>
+#include <utility>
 
 namespace pw {
+template <typename Scalar, unsigned int Rows, unsigned int Cols>
+struct matrix final {
 
+    static_assert(Rows > 0 && Cols > 0, "undefined matrix type");
+
+    static constexpr unsigned int diag_size{std::min(Rows, Cols)};
+    static constexpr unsigned int rows{Rows};
+    static constexpr unsigned int cols{Cols};
+    static constexpr bool is_square{Rows == Cols};
+
+    using diag_type      = vector<Scalar, diag_size>;
+    using transpose_type = matrix<Scalar, Cols, Rows>;
+    using minor_type     = matrix<Scalar, Cols - 1, Rows - 1>;
+    using column_type    = vector<Scalar, Rows>;
+
+    vector<vector<Scalar, Cols>, Rows> m_{};
+
+    auto&& data(this auto&& self) {
+        return std::forward<decltype(self)>(self).m_[0].data();
+    }
+
+    auto&& operator[](this auto&& self, int r) {
+        return std::forward<decltype(self)>(self).m_[r];
+    }
+
+    constexpr auto diagonal() const noexcept -> diag_type {
+        return [this]<std::size_t... Is>(std::index_sequence<Is...>) {
+            return diag_type{(*this)[Is][Is]...};
+        }(std::make_index_sequence<diag_size>{});
+    }
+
+    constexpr auto trace() const noexcept -> Scalar {
+        static_assert(Rows == Cols,
+                      "trace is only defined for square matrices");
+        return [this, sum = Scalar{0}]<std::size_t... Is>(
+                   std::index_sequence<Is...>) mutable {
+            return ((*this)[Is][Is] + ... + sum);
+        }(std::make_index_sequence<diag_size>{});
+    }
+
+    constexpr auto
+    column(std::unsigned_integral auto c) const noexcept -> column_type {
+        return [this, &c]<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            return column_type{(*this)[Rs][c]...};
+        }(std::make_index_sequence<Cols>{});
+    }
+
+    constexpr auto transposed() const noexcept -> transpose_type {
+        return [this]<std::size_t... Cs>(std::index_sequence<Cs...>) {
+            return transpose_type{(*this).column(Cs)...};
+        }(std::make_index_sequence<Cols>{});
+    }
+
+    static constexpr auto all(const Scalar& value) {
+        return [&value]<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            return matrix{vector<Scalar, Cols>::all((value) + Rs - Rs)...};
+        }(std::make_index_sequence<Rows>{});
+    }
+
+    static constexpr auto identity() noexcept -> matrix {
+        static_assert(Rows == Cols,
+                      "identity only defined for square matrices");
+        return []<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            return matrix{vector<Scalar, Cols>::basis(Rs)...};
+        }(std::make_index_sequence<Rows>{});
+    }
+
+    constexpr auto minor(std::unsigned_integral auto r0,
+                         std::unsigned_integral auto c0) const noexcept {
+        static_assert(Rows > 1 && Cols > 1, "cannot create minor matrix");
+        return [this, &r0, &c0]<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            return matrix<Scalar, Rows - 1, Cols - 1>{
+                (*this)[(Rs < r0) ? Rs : Rs + 1].minor(c0)...};
+        }(std::make_index_sequence<Rows - 1>{});
+    }
+
+    constexpr auto determinant() const noexcept -> Scalar {
+        static_assert(Rows == Cols,
+                      "determinant only defined for square matrices");
+        return [this]<std::size_t... Cs>(std::index_sequence<Cs...>) {
+            return ((((Cs % 2 == 0) ? (*this)[0][Cs] : -(*this)[0][Cs]) *
+                     (*this).minor(0u, Cs).determinant()) +
+                    ...);
+        }(std::make_index_sequence<Cols>{});
+    }
+
+    constexpr auto inverse() const noexcept -> matrix {
+        const auto inv_det{Scalar(1) / this->determinant()};
+        matrix<Scalar, Rows, Cols> inv; // no initialization needed!
+        for (auto c = 0u; c < Cols; c++) {
+            for (auto r = 0u; r < Rows; r++) {
+                const auto min_det = this->minor(r, c).determinant();
+                const auto co_fact{((r + c) % 2 == 1) ? -min_det : min_det};
+                inv[r][c] = inv_det * co_fact;
+            }
+        }
+        return inv;
+    }
+
+    constexpr matrix operator*(const Scalar& v) const noexcept {
+        return [this, &v]<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            return matrix{((*this)[Rs] * v)...};
+        }(std::make_index_sequence<Rows>{});
+    }
+    constexpr matrix operator/(const Scalar& v) const noexcept {
+        return operator*(Scalar{1} / v);
+    }
+
+    constexpr matrix operator+(const Scalar& v) const noexcept {
+        return [this, &v]<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            return matrix{((*this)[Rs] + v)...};
+        }(std::make_index_sequence<Rows>{});
+    }
+    constexpr matrix operator-(const Scalar& v) const noexcept {
+        return operator+(-v);
+    }
+
+    constexpr auto operator*=(const Scalar& v) noexcept {
+        [this, &v]<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            (((*this)[Rs] *= v), ...);
+        }(std::make_index_sequence<Rows>{});
+        return *this;
+    }
+
+    constexpr auto operator/=(const Scalar& v) noexcept {
+        return operator*=(Scalar{1} / v);
+    }
+
+    constexpr auto operator+=(const Scalar& v) noexcept {
+        [this, &v]<std::size_t... Rs>(std::index_sequence<Rs...>) {
+            (((*this)[Rs] += v), ...);
+        }(std::make_index_sequence<Rows>{});
+        return *this;
+    }
+
+    constexpr auto operator-=(const Scalar& v) noexcept {
+        return operator+=(-v);
+    }
+};
+
+template <> constexpr float matrix<float, 1, 1>::determinant() const noexcept {
+    return (*this)[0][0];
+}
+
+template <>
+constexpr double matrix<double, 1, 1>::determinant() const noexcept {
+    return (*this)[0][0];
+}
+
+template <typename Scalar, unsigned int Ra, unsigned int CaRb, unsigned int Cb>
+constexpr auto operator*(const matrix<Scalar, Ra, CaRb>& A,
+                         const matrix<Scalar, CaRb, Cb>& B) noexcept {
+    matrix<Scalar, Ra, Cb> result{};
+    [&]<std::size_t... Ia>(std::index_sequence<Ia...>) {
+        (
+            [&]<std::size_t... Ib>(auto&& Ai, std::index_sequence<Ib...>) {
+            (
+                [&]<std::size_t... iI>(auto&& Ai, auto&& Bi,
+                                       std::index_sequence<iI...>) {
+                ((result[Ai][Bi] += A[Ai][iI] * B[iI][Bi]),
+                 ...); // do the actual multiplication
+            }(Ai, Ib, std::make_index_sequence<CaRb>{}),
+                ...); // forward row of A, colum of B and unpack inner loop
+        }(Ia, std::make_index_sequence<Cb>{}),
+            ...); // forward current row, unpack A
+    }(std::make_index_sequence<Ra>{}); // rows of A as input
+    return result;
+}
+
+template <typename ScalarA, typename ScalarB, unsigned int Ra, unsigned int Ca>
+constexpr auto operator*(const matrix<ScalarA, Ra, Ca>& A,
+                         const vector<ScalarB, Ca>& B) noexcept {
+    // first step - should move to concepts to allow for std::array and
+    // std::span to be allowed
+    vector<std::common_type_t<ScalarA, ScalarB>, Ca> result{};
+    [&]<std::size_t... Ia>(std::index_sequence<Ia...>) {
+        (
+            [&]<std::size_t... Ib>(auto&& Ai, std::index_sequence<Ib...>) {
+            ((result[Ib] += A[Ai][Ib] * B[Ib]), ...);
+        }(Ia, std::make_index_sequence<Ca>{}),
+            ...); // forward current row, unpack A
+    }(std::make_index_sequence<Ra>{}); // rows of A as input
+    return result;
+}
 
 } // namespace pw
 

src/core/include/pw/core/point.hpp

@@ -34,9 +34,6 @@ struct point_ {
 
     T_ x {0}, y {0};
 
-    point_() = default;
-	point_(T_ x_,T_ y_) : x(x_), y(y_) {}
-
     template <typename To_>
     point_<To_> cast() const { return point_<To_>(static_cast<To_>(x),static_cast<To_>(y)); }
 

src/core/include/pw/core/vector.hpp

@@ -23,6 +23,8 @@
 #ifndef PW_CORE_VECTOR_HPP
 #define PW_CORE_VECTOR_HPP
 
+#include <pw/core/globals.hpp>
+
 namespace pw {
 
 template <typename Scalar, unsigned int size> struct vector final {
@@ -112,10 +114,56 @@ template <typename Scalar, unsigned int size> struct vector final {
             return vector{{(*this)[Ss] - v[Ss]...}};
         }(std::make_index_sequence<size>{});
     }
+
+    constexpr vector& operator*=(const Scalar& v) noexcept {
+        [this, &v]<std::size_t... Ss>(std::index_sequence<Ss...>) {
+            (((*this)[Ss] *= v), ...);
+        }(std::make_index_sequence<size>{});
+        return *this;
+    }
+
+    constexpr vector& operator/=(const Scalar& v) noexcept {
+        return operator*=(Scalar{1} / v);
+    }
+
+    constexpr vector& operator+=(const Scalar& v) noexcept {
+        [this, &v]<std::size_t... Ss>(std::index_sequence<Ss...>) {
+            (((*this)[Ss] += v), ...);
+        }(std::make_index_sequence<size>{});
+        return *this;
+    }
+
+    constexpr vector& operator-=(const Scalar& v) noexcept {
+        return operator+=(-v);
+    }
+
+    constexpr auto project() const noexcept -> vector<Scalar, size - 1> {
+        return [this]<std::size_t... Ss>(std::index_sequence<Ss...>) {
+            return vector<Scalar, size - 1>{
+                {(*this)[Ss] / (*this)[size - 1]...}};
+        }(std::make_index_sequence<size - 1>{});
+    }
+
+    constexpr auto
+    unproject(auto&& w) const noexcept -> vector<Scalar, size + 1> {
+        return [this, &w]<std::size_t... Ss>(std::index_sequence<Ss...>) {
+            return vector<Scalar, size + 1>{{(Ss < size) ? (*this)[Ss] : w...}};
+        }(std::make_index_sequence<size + 1>{});
+    }
+
+    template <typename ScalarCast>
+    constexpr auto cast() const noexcept {
+        return [this]<std::size_t... Ss>(std::index_sequence<Ss...>) {
+            return vector{ ScalarCast{(*this)[Ss]}...};
+        }(std::make_index_sequence<size>{});
+    }
 };
 
-template <typename Scalar>
+// deduction guide for vector
-struct vector2 : vector<Scalar, 2> {
+template <class... U, class CT = std::common_type_t<U...>>
+vector(U...) -> vector<CT, sizeof...(U)>;
+
+template <typename Scalar> struct vector2 : vector<Scalar, 2> {
 
     constexpr const Scalar& x() const { return (*this)[0]; }
     constexpr Scalar& x() { return (*this)[0]; }
@@ -159,11 +207,9 @@ template <typename Scalar> struct vector4 : vector<Scalar, 4> {
     constexpr const Scalar& w() const { return (*this)[3]; }
     constexpr Scalar& w() { return (*this)[3]; }
 
-    constexpr auto xyz() const { return base_type::swizzle(0,1,2); }
+    constexpr auto xyz() const { return base_type::swizzle(0, 1, 2); }
 };
 
-
-
 using vector2f = vector2<float>;
 using vector2d = vector2<double>;
 using vector2i = vector2<int>;
@@ -171,7 +217,6 @@ using vector2i      = vector2<int>;
 using vector3f = vector3<float>;
 using vector3d = vector3<double>;
 
-
 } // namespace pw
 
 #endif

src/core/tests/CMakeLists.txt

@@ -7,8 +7,8 @@ macro(make_test arg1)
 endmacro()
 
 
-# make_test(pwcore_test_matrix)
+make_test(pwcore_test_matrix)
-# make_test(pwcore_test_vector)
+make_test(pwcore_test_vector)
 # make_test(pwcore_test_axisangle)
 # make_test(pwcore_test_quaternion)
 # make_test(pwcore_test_transform_tools)

tests/CMakeLists.txt

@@ -10,9 +10,6 @@ set(pixwerx_test_files
     pw_test_core.cpp
   )
 
-
-
-
 add_executable(pixwerx_tests ${pixwerx_test_files})
 
 target_link_libraries(pixwerx_tests PRIVATE

tests/pw_test_core.cpp

@@ -3,42 +3,55 @@
 #include <pw/core/matrix.hpp>
 #include <pw/core/vector.hpp>
 
-
-
 TEST_CASE("core", "[matrix]") {
 
-    auto m44 = pw::matrix<float,2,2>{};
+    auto mat = pw::matrix<float,2,2>{};
 
-    REQUIRE(m44[0][0] == 0);
+    REQUIRE(mat[0][0] == 0);
-    REQUIRE(m44[0][1] == 0);
+    REQUIRE(mat[0][1] == 0);
-    REQUIRE(m44[1][0] == 0);
+    REQUIRE(mat[1][0] == 0);
-    REQUIRE(m44[1][1] == 0);
+    REQUIRE(mat[1][1] == 0);
 
-    // auto m44_1 = pw::matrix<2, 2, float>{ 5, 4, 3, 2 };
+    auto mat_1 = pw::matrix<float,2, 2>{
+        pw::vector{5.f, 4},
+        pw::vector{3.f, 2}
+    };
 
-    // REQUIRE(m44_1(0, 0) == 5);
+    REQUIRE(mat_1[0][0] == 5);
-    // REQUIRE(m44_1(0, 1) == 3);
+    REQUIRE(mat_1[0][1] == 4);
-    // REQUIRE(m44_1(1, 0) == 4);
+    REQUIRE(mat_1[1][0] == 3);
-    // REQUIRE(m44_1(1, 1) == 2);
+    REQUIRE(mat_1[1][1] == 2);
 
-    // auto m44_2 = pw::matrix2x2_<float>::all(42.42f);
+    const auto val {42.42f};
 
-    // REQUIRE(m44_2(0, 0) == 42.42f);
+    auto mat_2 = pw::matrix<float,2,2>::all(val);
-    // REQUIRE(m44_2(0, 1) == 42.42f);
-    // REQUIRE(m44_2(1, 0) == 42.42f);
-    // REQUIRE(m44_2(1, 1) == 42.42f);
 
-    // m44 *= 2;
+    REQUIRE(mat_2[0][0] == val);
-    // REQUIRE(m44(0, 0) == 84.84f);
+    REQUIRE(mat_2[0][1] == val);
-    // REQUIRE(m44(0, 1) == 84.84f);
+    REQUIRE(mat_2[1][0] == val);
-    // REQUIRE(m44(1, 0) == 84.84f);
+    REQUIRE(mat_2[1][1] == val);
-    // REQUIRE(m44(1, 1) == 84.84f);
 
-    // m44.set_identity();
+    mat_2 *= 2;
-    // REQUIRE(m44(0, 0) == 1.0f);
+    REQUIRE(mat_2[0][0] == (val * 2));
-    // REQUIRE(m44(0, 1) == 0.0f);
+    REQUIRE(mat_2[0][1] == (val * 2));
-    // REQUIRE(m44(1, 0) == 0.0f);
+    REQUIRE(mat_2[1][0] == (val * 2));
-    // REQUIRE(m44(1, 1) == 1.0f);
+    REQUIRE(mat_2[1][1] == (val * 2));
+
+
+    mat_2 = pw::matrix<float, 2, 2>::identity();
+
+    REQUIRE(mat_2[0][0] == 1.0f);
+    REQUIRE(mat_2[0][1] == 0.0f);
+    REQUIRE(mat_2[1][0] == 0.0f);
+    REQUIRE(mat_2[1][1] == 1.0f);
+
+
+    auto mat_2_inv = mat_2.inverse();
+
+    REQUIRE(mat_2_inv[0][0] == 1.0f);
+    REQUIRE(mat_2_inv[0][1] == 0.0f);
+    REQUIRE(mat_2_inv[1][0] == 0.0f);
+    REQUIRE(mat_2_inv[1][1] == 1.0f);
 
     // auto m44_2 = pw::matrix2x2::make(0, 1, 2, 3);
 
@@ -60,7 +73,14 @@ TEST_CASE("core", "[matrix]") {
 
 TEST_CASE("core", "[vector.matrix]") {
 
-    // auto vec4_1 = pw::vector4_<float>{};
+    auto vec4_1 = pw::vector4<float>{};
+
+    REQUIRE(vec4_1[0] == 0.0f);
+    REQUIRE(vec4_1[1] == 0.0f);
+    REQUIRE(vec4_1[2] == 0.0f);
+    REQUIRE(vec4_1[3] == 0.0f);
+
+
     // auto mat4_1 = pw::matrix_<4,1,float>{};