From 093416ceee27970b5aac29be54ab494d1f159624 Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sat, 27 Feb 2021 13:19:46 +0700
Subject: [PATCH] update singles

---
 singles/vmath.hpp/vmath.hpp | 427 ++++++++++++++++++++++++------------
 untests/CMakeLists.txt      |   2 +-
 2 files changed, 291 insertions(+), 138 deletions(-)

diff --git a/singles/vmath.hpp/vmath.hpp b/singles/vmath.hpp/vmath.hpp
index b85f903..c68a0eb 100644
--- a/singles/vmath.hpp/vmath.hpp
+++ b/singles/vmath.hpp/vmath.hpp
@@ -10,6 +10,7 @@
 #include <cstddef>
 #include <cstdlib>
 
+#include <functional>
 #include <initializer_list>
 #include <iterator>
 #include <limits>
@@ -590,6 +591,9 @@ namespace vmath_hpp::detail
         template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
         constexpr explicit vec_base(const vec_base<U, 4>& other): vec_base(other[0], other[1]) {}
 
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit vec_base(const U* p): vec_base(p[0], p[1]) {}
+
         [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
             switch ( index ) {
             default:
@@ -631,6 +635,9 @@ namespace vmath_hpp::detail
         template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
         constexpr explicit vec_base(const vec_base<U, 4>& other): vec_base(other[0], other[1], other[2]) {}
 
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit vec_base(const U* p): vec_base(p[0], p[1], p[2]) {}
+
         [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
             switch ( index ) {
             default:
@@ -676,6 +683,9 @@ namespace vmath_hpp::detail
         template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
         constexpr vec_base(const vec_base<U, 4>& other): vec_base(other[0], other[1], other[2], other[3]) {}
 
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit vec_base(const U* p): vec_base(p[0], p[1], p[2], p[3]) {}
+
         [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
             switch ( index ) {
             default:
@@ -1276,6 +1286,64 @@ namespace vmath_hpp
         return (xs = (xs ^ ys));
     }
 
+    // operator<<
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator<<(const vec<T, Size>& xs, T y) {
+        return map_join([y](T x){ return x << y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator<<(T x, const vec<T, Size>& ys) {
+        return map_join([x](T y){ return x << y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator<<(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return x << y; }, xs, ys);
+    }
+
+    // operator<<=
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size>& operator<<=(vec<T, Size>& xs, T y) {
+        return (xs = (xs << y));
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size>& operator<<=(vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return (xs = (xs << ys));
+    }
+
+    // operator>>
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator>>(const vec<T, Size>& xs, T y) {
+        return map_join([y](T x){ return x >> y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator>>(T x, const vec<T, Size>& ys) {
+        return map_join([x](T y){ return x >> y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator>>(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return x >> y; }, xs, ys);
+    }
+
+    // operator>>=
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size>& operator>>=(vec<T, Size>& xs, T y) {
+        return (xs = (xs >> y));
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size>& operator>>=(vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return (xs = (xs >> ys));
+    }
+
     // operator&&
 
     template < typename T, std::size_t Size >
@@ -1884,6 +1952,11 @@ namespace vmath_hpp::detail
         constexpr explicit mat_base(const mat_base<U, 4>& other): mat_base(
             row_type{other.rows[0]},
             row_type{other.rows[1]}) {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit mat_base(const U* p): mat_base(
+            row_type{p + 0u * row_type::size},
+            row_type{p + 1u * row_type::size}) {}
     };
 
     template < typename T >
@@ -1952,6 +2025,12 @@ namespace vmath_hpp::detail
             row_type{other.rows[0]},
             row_type{other.rows[1]},
             row_type{other.rows[2]}) {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit mat_base(const U* p): mat_base(
+            row_type{p + 0u * row_type::size},
+            row_type{p + 1u * row_type::size},
+            row_type{p + 2u * row_type::size}) {}
     };
 
     template < typename T >
@@ -2029,6 +2108,13 @@ namespace vmath_hpp::detail
             row_type{other.rows[1], T{0}},
             row_type{other.rows[2], T{0}},
             row_type{T{0}, T{0}, T{0}, T{1}}) {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit mat_base(const U* p): mat_base(
+            row_type{p + 0u * row_type::size},
+            row_type{p + 1u * row_type::size},
+            row_type{p + 2u * row_type::size},
+            row_type{p + 3u * row_type::size}) {}
     };
 }
 
@@ -2566,6 +2652,64 @@ namespace vmath_hpp
         return (xs = (xs ^ ys));
     }
 
+    // operator<<
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator<<(const mat<T, Size>& xs, T y) {
+        return map_join([y](const vec<T, Size>& x){ return x << y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator<<(T x, const mat<T, Size>& ys) {
+        return map_join([x](const vec<T, Size>& y){ return x << y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator<<(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return x << y; }, xs, ys);
+    }
+
+    // operator<<=
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size>& operator<<=(mat<T, Size>& xs, T y) {
+        return (xs = (xs << y));
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size>& operator<<=(mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return (xs = (xs << ys));
+    }
+
+    // operator>>
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator>>(const mat<T, Size>& xs, T y) {
+        return map_join([y](const vec<T, Size>& x){ return x >> y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator>>(T x, const mat<T, Size>& ys) {
+        return map_join([x](const vec<T, Size>& y){ return x >> y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator>>(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return x >> y; }, xs, ys);
+    }
+
+    // operator>>=
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size>& operator>>=(mat<T, Size>& xs, T y) {
+        return (xs = (xs >> y));
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size>& operator>>=(mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return (xs = (xs >> ys));
+    }
+
     // operator&&
 
     template < typename T, std::size_t Size >
@@ -2720,6 +2864,10 @@ namespace vmath_hpp
 
 namespace vmath_hpp
 {
+    //
+    // transpose
+    //
+
     namespace impl
     {
         template < typename T >
@@ -2786,6 +2934,98 @@ namespace vmath_hpp
             m[3][0], m[3][1], m[3][2], m[3][3]);
     }
 
+    //
+    // adjugate
+    //
+
+    namespace impl
+    {
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        mat<T, 2> adjugate_2x2_impl(
+            T a, T b,
+            T c, T d)
+        {
+            return {
+                +d, -b,
+                -c, +a};
+        }
+
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        mat<T, 3> adjugate_3x3_impl(
+            T a, T b, T c,
+            T d, T e, T f,
+            T g, T h, T i)
+        {
+            return {
+                e * i - f * h,
+                c * h - b * i,
+                b * f - c * e,
+                f * g - d * i,
+                a * i - c * g,
+                c * d - a * f,
+                d * h - e * g,
+                b * g - a * h,
+                a * e - b * d};
+        }
+
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        mat<T, 4> adjugate_4x4_impl(
+            T a, T b, T c, T d,
+            T e, T f, T g, T h,
+            T i, T j, T k, T l,
+            T m, T n, T o, T p)
+        {
+            return {
+                f * (k * p - l * o) + g * (l * n - j * p) + h * (j * o - k * n),
+                j * (c * p - d * o) + k * (d * n - b * p) + l * (b * o - c * n),
+                n * (c * h - d * g) + o * (d * f - b * h) + p * (b * g - c * f),
+                b * (h * k - g * l) + c * (f * l - h * j) + d * (g * j - f * k),
+                g * (i * p - l * m) + h * (k * m - i * o) + e * (l * o - k * p),
+                k * (a * p - d * m) + l * (c * m - a * o) + i * (d * o - c * p),
+                o * (a * h - d * e) + p * (c * e - a * g) + m * (d * g - c * h),
+                c * (h * i - e * l) + d * (e * k - g * i) + a * (g * l - h * k),
+                h * (i * n - j * m) + e * (j * p - l * n) + f * (l * m - i * p),
+                l * (a * n - b * m) + i * (b * p - d * n) + j * (d * m - a * p),
+                p * (a * f - b * e) + m * (b * h - d * f) + n * (d * e - a * h),
+                d * (f * i - e * j) + a * (h * j - f * l) + b * (e * l - h * i),
+                e * (k * n - j * o) + f * (i * o - k * m) + g * (j * m - i * n),
+                i * (c * n - b * o) + j * (a * o - c * m) + k * (b * m - a * n),
+                m * (c * f - b * g) + n * (a * g - c * e) + o * (b * e - a * f),
+                a * (f * k - g * j) + b * (g * i - e * k) + c * (e * j - f * i)};
+        }
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 2> adjugate(const mat<T, 2>& m) {
+        return impl::adjugate_2x2_impl(
+            m[0][0], m[0][1],
+            m[1][0], m[1][1]);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> adjugate(const mat<T, 3>& m) {
+        return impl::adjugate_3x3_impl(
+            m[0][0], m[0][1], m[0][2],
+            m[1][0], m[1][1], m[1][2],
+            m[2][0], m[2][1], m[2][2]);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> adjugate(const mat<T, 4>& m) {
+        return impl::adjugate_4x4_impl(
+            m[0][0], m[0][1], m[0][2], m[0][3],
+            m[1][0], m[1][1], m[1][2], m[1][3],
+            m[2][0], m[2][1], m[2][2], m[2][3],
+            m[3][0], m[3][1], m[3][2], m[3][3]);
+    }
+
+    //
+    // determinant
+    //
+
     namespace impl
     {
         template < typename T >
@@ -2794,12 +3034,7 @@ namespace vmath_hpp
             T a, T b,
             T c, T d)
         {
-            /// REFERENCE:
-            /// http://www.euclideanspace.com/maths/algebra/matrix/functions/determinant/twoD/
-
-            return
-                + a * d
-                - b * c;
+            return a * d - b * c;
         }
 
         template < typename T >
@@ -2809,14 +3044,10 @@ namespace vmath_hpp
             T d, T e, T f,
             T g, T h, T i)
         {
-            /// REFERENCE:
-            /// http://www.euclideanspace.com/maths/algebra/matrix/functions/determinant/threeD/
-
             return
-                + a * determinant_2x2_impl(e, f, h, i)
-                - b * determinant_2x2_impl(d, f, g, i)
-                + c * determinant_2x2_impl(d, e, g, h);
-
+                + a * (e * i - f * h)
+                - b * (d * i - f * g)
+                + c * (d * h - e * g);
         }
 
         template < typename T >
@@ -2827,14 +3058,11 @@ namespace vmath_hpp
             T i, T j, T k, T l,
             T m, T n, T o, T p)
         {
-            /// REFERENCE:
-            /// http://www.euclideanspace.com/maths/algebra/matrix/functions/determinant/fourD/
-
             return
-                + a * determinant_3x3_impl(f, g, h, j, k, l, n, o, p)
-                - b * determinant_3x3_impl(e, g, h, i, k, l, m, o, p)
-                + c * determinant_3x3_impl(e, f, h, i, j, l, m, n, p)
-                - d * determinant_3x3_impl(e, f, g, i, j, k, m, n, o);
+                + a * (f * (k * p - l * o) - (j * (g * p - h * o)) + (n * (g * l - h * k)))
+                - b * (e * (k * p - l * o) - (i * (g * p - h * o)) + (m * (g * l - h * k)))
+                + c * (e * (j * p - l * n) - (i * (f * p - h * n)) + (m * (f * l - h * j)))
+                - d * (e * (j * o - k * n) - (i * (f * o - g * n)) + (m * (f * k - g * j)));
         }
     }
 
@@ -2862,118 +3090,13 @@ namespace vmath_hpp
             m[3][0], m[3][1], m[3][2], m[3][3]);
     }
 
-    namespace impl
-    {
-        template < typename T >
-        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
-        mat<T, 2> inverse_2x2_impl(
-            T a, T b,
-            T c, T d)
-        {
-            /// REFERENCE:
-            /// http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/twoD/
+    //
+    // inverse
+    //
 
-            const T inv_det = rcp(determinant_2x2_impl(
-                a, b,
-                c, d));
-
-            const mat inv_m{
-                d, -b,
-                -c, a};
-
-            return inv_m * inv_det;
-        }
-
-        template < typename T >
-        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
-        mat<T, 3> inverse_3x3_impl(
-            T a, T b, T c,
-            T d, T e, T f,
-            T g, T h, T i)
-        {
-            /// REFERENCE:
-            /// http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/threeD/
-
-            const T inv_det = rcp(determinant_3x3_impl(
-                a, b, c,
-                d, e, f,
-                g, h, i));
-
-            const mat inv_m{
-                e * i - f * h,
-                c * h - b * i,
-                b * f - c * e,
-                f * g - d * i,
-                a * i - c * g,
-                c * d - a * f,
-                d * h - e * g,
-                b * g - a * h,
-                a * e - b * d};
-
-            return inv_m * inv_det;
-        }
-
-        template < typename T >
-        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
-        mat<T, 4> inverse_4x4_impl(
-            T a, T b, T c, T d,
-            T e, T f, T g, T h,
-            T i, T j, T k, T l,
-            T m, T n, T o, T p)
-        {
-            /// REFERENCE:
-            /// http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/
-
-            const T inv_det = rcp(determinant_4x4_impl(
-                a, b, c, d,
-                e, f, g, h,
-                i, j, k, l,
-                m, n, o, p));
-
-            const mat inv_m{
-                (f * (k * p - l * o) + g * (l * n - j * p) + h * (j * o - k * n)),
-                (j * (c * p - d * o) + k * (d * n - b * p) + l * (b * o - c * n)),
-                (n * (c * h - d * g) + o * (d * f - b * h) + p * (b * g - c * f)),
-                (b * (h * k - g * l) + c * (f * l - h * j) + d * (g * j - f * k)),
-                (g * (i * p - l * m) + h * (k * m - i * o) + e * (l * o - k * p)),
-                (k * (a * p - d * m) + l * (c * m - a * o) + i * (d * o - c * p)),
-                (o * (a * h - d * e) + p * (c * e - a * g) + m * (d * g - c * h)),
-                (c * (h * i - e * l) + d * (e * k - g * i) + a * (g * l - h * k)),
-                (h * (i * n - j * m) + e * (j * p - l * n) + f * (l * m - i * p)),
-                (l * (a * n - b * m) + i * (b * p - d * n) + j * (d * m - a * p)),
-                (p * (a * f - b * e) + m * (b * h - d * f) + n * (d * e - a * h)),
-                (d * (f * i - e * j) + a * (h * j - f * l) + b * (e * l - h * i)),
-                (e * (k * n - j * o) + f * (i * o - k * m) + g * (j * m - i * n)),
-                (i * (c * n - b * o) + j * (a * o - c * m) + k * (b * m - a * n)),
-                (m * (c * f - b * g) + n * (a * g - c * e) + o * (b * e - a * f)),
-                (a * (f * k - g * j) + b * (g * i - e * k) + c * (e * j - f * i))};
-
-            return inv_m * inv_det;
-        }
-    }
-
-    template < typename T >
-    [[nodiscard]] constexpr mat<T, 2> inverse(const mat<T, 2>& m) {
-        return impl::inverse_2x2_impl(
-            m[0][0], m[0][1],
-            m[1][0], m[1][1]);
-    }
-
-    template < typename T >
-    [[nodiscard]] constexpr mat<T, 3> inverse(const mat<T, 3>& m) {
-        return impl::inverse_3x3_impl(
-            m[0][0], m[0][1], m[0][2],
-            m[1][0], m[1][1], m[1][2],
-            m[2][0], m[2][1], m[2][2]);
-    }
-
-    template < typename T >
-    [[nodiscard]] constexpr mat<T, 4> inverse(const mat<T, 4>& m) {
-        return impl::inverse_4x4_impl(
-            m[0][0], m[0][1], m[0][2], m[0][3],
-            m[1][0], m[1][1], m[1][2], m[1][3],
-            m[2][0], m[2][1], m[2][2], m[2][3],
-            m[3][0], m[3][1], m[3][2], m[3][3]);
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr mat<T, Size> inverse(const mat<T, Size>& m) {
+        return adjugate(m) * rcp(determinant(m));
     }
 }
 
@@ -3002,6 +3125,9 @@ namespace vmath_hpp::detail
         template < typename U, std::enable_if_t<std::is_convertible_v<T, U>, int> = 0 >
         constexpr explicit operator vec<U, 4>() const { return vec<U, 4>(v, s); }
 
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit qua_base(const U* p): qua_base(p[0], p[1], p[2], p[3]) {}
+
         [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
             switch ( index ) {
             default:
@@ -3612,9 +3738,10 @@ namespace vmath_hpp
     }
 
     template < typename T, std::size_t Size >
-    [[nodiscard]] constexpr vec<T, Size> component(vec<T, Size> v, std::size_t index, T component) {
-        v[index] = component;
-        return v;
+    [[nodiscard]] constexpr vec<T, Size> component(const vec<T, Size>& v, std::size_t index, T component) {
+        vec vv = v;
+        vv[index] = component;
+        return vv;
     }
 
     // row
@@ -3625,9 +3752,10 @@ namespace vmath_hpp
     }
 
     template < typename T, std::size_t Size >
-    [[nodiscard]] constexpr mat<T, Size> row(mat<T, Size> m, std::size_t index, const vec<T, Size>& row) {
-        m.rows[index] = row;
-        return m;
+    [[nodiscard]] constexpr mat<T, Size> row(const mat<T, Size>& m, std::size_t index, const vec<T, Size>& row) {
+        mat mm = m;
+        mm[index] = row;
+        return mm;
     }
 
     // column
@@ -3657,6 +3785,31 @@ namespace vmath_hpp
         return impl::column_impl(m, index, column, std::make_index_sequence<Size>{});
     }
 
+    // diagonal
+
+    namespace impl
+    {
+        template < typename T, std::size_t Size, std::size_t... Is >
+        [[nodiscard]] constexpr vec<T, Size> diagonal_impl(const mat<T, Size>& m, std::index_sequence<Is...>) {
+            return { m[Is][Is]... };
+        }
+
+        template < typename T, std::size_t Size, std::size_t... Is >
+        [[nodiscard]] constexpr mat<T, Size> diagonal_impl(const mat<T, Size>& m, const vec<T, Size>& diagonal, std::index_sequence<Is...>) {
+            return { component(m[Is], Is, diagonal[Is])... };
+        }
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> diagonal(const mat<T, Size>& m) {
+        return impl::diagonal_impl(m, std::make_index_sequence<Size>{});
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr mat<T, Size> diagonal(const mat<T, Size>& m, const vec<T, Size>& diagonal) {
+        return impl::diagonal_impl(m, diagonal, std::make_index_sequence<Size>{});
+    }
+
     // real
 
     template < typename T >
diff --git a/untests/CMakeLists.txt b/untests/CMakeLists.txt
index edc5959..0220203 100644
--- a/untests/CMakeLists.txt
+++ b/untests/CMakeLists.txt
@@ -18,7 +18,7 @@ endif()
 
 file(GLOB_RECURSE UNTESTS_SOURCES "*.cpp" "*.hpp")
 add_executable(${PROJECT_NAME} ${UNTESTS_SOURCES})
-target_link_libraries(${PROJECT_NAME} vmath.hpp)
+target_link_libraries(${PROJECT_NAME} vmath.hpp.singles)
 
 target_compile_options(${PROJECT_NAME}
     PRIVATE