From f9dd3b610510a602d89253506e096102efd69183 Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sun, 6 Dec 2020 02:59:17 +0700
Subject: [PATCH 1/6] add rotate_x, rotate_y, rotate_z functions

---
 README.md                       | 18 +++++++++++++
 headers/vmath.hpp/vmath_ext.hpp | 48 +++++++++++++++++++++++++++++++++
 untests/vmath_ext_tests.cpp     |  9 +++++++
 3 files changed, 75 insertions(+)

diff --git a/README.md b/README.md
index a345c56..e88e20f 100644
--- a/README.md
+++ b/README.md
@@ -1147,6 +1147,24 @@ mat<T, 4> rotate(T angle, const vec<T, 3>& axis);
 template < typename T >
 mat<T, 4> rotate(const mat<T, 4>& m, T angle, const vec<T, 3>& axis);
 
+template < typename T >
+mat<T, 4> rotate_x(T angle);
+
+template < typename T >
+mat<T, 4> rotate_x(const mat<T, 4>& m, T angle);
+
+template < typename T >
+mat<T, 4> rotate_y(T angle);
+
+template < typename T >
+mat<T, 4> rotate_y(const mat<T, 4>& m, T angle);
+
+template < typename T >
+mat<T, 4> rotate_z(T angle);
+
+template < typename T >
+mat<T, 4> rotate_z(const mat<T, 4>& m, T angle);
+
 template < typename T >
 constexpr mat<T, 4> scale(T x, T y, T z);
 
diff --git a/headers/vmath.hpp/vmath_ext.hpp b/headers/vmath.hpp/vmath_ext.hpp
index d150895..f79bf57 100644
--- a/headers/vmath.hpp/vmath_ext.hpp
+++ b/headers/vmath.hpp/vmath_ext.hpp
@@ -240,6 +240,54 @@ namespace vmath_hpp
         return m * rotate(angle, axis);
     }
 
+    template < typename T >
+    [[nodiscard]] mat<T, 4> rotate_x(T angle) {
+        const T cs = cos(angle);
+        const T sn = sin(angle);
+        return {
+            1, 0,   0,  0,
+            0, cs,  sn, 0,
+            0, -sn, cs, 0,
+            0, 0,   0,  1};
+    }
+
+    template < typename T >
+    [[nodiscard]] mat<T, 4> rotate_x(const mat<T, 4>& m, T angle) {
+        return m * rotate_x(angle);
+    }
+
+    template < typename T >
+    [[nodiscard]] mat<T, 4> rotate_y(T angle) {
+        const T cs = cos(angle);
+        const T sn = sin(angle);
+        return {
+            cs, 0, -sn, 0,
+            0,  1, 0,   0,
+            sn, 0, cs,  0,
+            0,  0, 0,   1};
+    }
+
+    template < typename T >
+    [[nodiscard]] mat<T, 4> rotate_y(const mat<T, 4>& m, T angle) {
+        return m * rotate_y(angle);
+    }
+
+    template < typename T >
+    [[nodiscard]] mat<T, 4> rotate_z(T angle) {
+        const T cs = cos(angle);
+        const T sn = sin(angle);
+        return {
+            cs,  sn, 0, 0,
+            -sn, cs, 0, 0,
+            0,   0,  1, 0,
+            0,   0,  0, 1};
+    }
+
+    template < typename T >
+    [[nodiscard]] mat<T, 4> rotate_z(const mat<T, 4>& m, T angle) {
+        return m * rotate_z(angle);
+    }
+
     // scale
 
     template < typename T >
diff --git a/untests/vmath_ext_tests.cpp b/untests/vmath_ext_tests.cpp
index abee61b..4a8516f 100644
--- a/untests/vmath_ext_tests.cpp
+++ b/untests/vmath_ext_tests.cpp
@@ -145,6 +145,15 @@ TEST_CASE("vmath/ext") {
     SUBCASE("matrix rotate") {
         constexpr float pi = radians(180.f);
         constexpr float pi_2 = radians(90.f);
+        constexpr float pi_4 = radians(45.f);
+
+        REQUIRE(float4(0.f,1.f,0.f,1.f) * rotate_x(pi_2) == approx4(0.f,0.f,1.f,1.f));
+        REQUIRE(float4(0.f,0.f,1.f,1.f) * rotate_y(pi_2) == approx4(1.f,0.f,0.f,1.f));
+        REQUIRE(float4(1.f,0.f,0.f,1.f) * rotate_z(pi_2) == approx4(0.f,1.f,0.f,1.f));
+
+        REQUIRE(float4(0.f,1.f,0.f,1.f) * rotate_x(rotate_x(pi_4),pi_4) == approx4(0.f,0.f,1.f,1.f));
+        REQUIRE(float4(0.f,0.f,1.f,1.f) * rotate_y(rotate_y(pi_4),pi_4) == approx4(1.f,0.f,0.f,1.f));
+        REQUIRE(float4(1.f,0.f,0.f,1.f) * rotate_z(rotate_z(pi_4),pi_4) == approx4(0.f,1.f,0.f,1.f));
 
         REQUIRE(float3(2.f,3.f,1.f) * rotate(pi) == approx3(-2.f,-3.f,1.f));
         REQUIRE(float4(2.f,3.f,4.f,1.f) * rotate(pi,{0.f,0.f,1.f}) == approx4(-2.f,-3.f,4.f,1.f));

From 80287ccaf8ffc325f0f68b253b5f0f57dd3b752e Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sun, 6 Dec 2020 03:07:08 +0700
Subject: [PATCH 2/6] vector project function

---
 README.md                       | 3 +++
 headers/vmath.hpp/vmath_ext.hpp | 7 +++++++
 untests/vmath_ext_tests.cpp     | 5 +++++
 3 files changed, 15 insertions(+)

diff --git a/README.md b/README.md
index e88e20f..e795a55 100644
--- a/README.md
+++ b/README.md
@@ -1286,6 +1286,9 @@ vec<T, 3> rotate(const vec<T, 3>& v, T angle, const vec<T, 3>& normal);
 
 template < typename T >
 vec<T, 4> rotate(const vec<T, 4>& v, T angle, const vec<T, 3>& normal);
+
+template < typename T, std::size_t Size >
+vec<T, Size> project(const vec<T, Size>& v, const vec<T, Size>& normal);
 ```
 
 ## [License (MIT)](./LICENSE.md)
diff --git a/headers/vmath.hpp/vmath_ext.hpp b/headers/vmath.hpp/vmath_ext.hpp
index f79bf57..435d2e5 100644
--- a/headers/vmath.hpp/vmath_ext.hpp
+++ b/headers/vmath.hpp/vmath_ext.hpp
@@ -636,4 +636,11 @@ namespace vmath_hpp
     [[nodiscard]] vec<T, 4> rotate(const vec<T, 4>& v, T angle, const vec<T, 3>& normal) {
         return v * rotate(angle, normal);
     }
+
+    // project
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] vec<T, Size> project(const vec<T, Size>& v, const vec<T, Size>& normal) {
+        return dot(v, normal) / length2(normal) * normal;
+    }
 }
diff --git a/untests/vmath_ext_tests.cpp b/untests/vmath_ext_tests.cpp
index 4a8516f..51446e9 100644
--- a/untests/vmath_ext_tests.cpp
+++ b/untests/vmath_ext_tests.cpp
@@ -221,4 +221,9 @@ TEST_CASE("vmath/ext") {
         REQUIRE(rotate(float3(1.5f,0.f,0.f), radians(90.f), float3(0,0,1)) == approx3(0.f,1.5f,0.f));
         REQUIRE(rotate(float4(1.5f,0.f,0.f,1.f), radians(90.f), float3(0,0,1)) == approx4(0.f,1.5f,0.f,1.f));
     }
+
+    SUBCASE("vector project") {
+        REQUIRE(project(float2(2.f, 2.f), float2(0.f, 1.f)) == approx2(0.f, 2.f));
+        REQUIRE(project(float3(2.f, 2.f, 2.f), float3(0.f, 0.f, 1.f)) == approx3(0.f, 0.f, 2.f));
+    }
 }

From b246cad6afef427402841190ce8823b5dbff2780 Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sun, 6 Dec 2020 04:04:19 +0700
Subject: [PATCH 3/6] add sincos

---
 README.md                           |  6 ++++++
 headers/vmath.hpp/vmath_fun.hpp     | 13 +++++++++++++
 headers/vmath.hpp/vmath_vec_fun.hpp |  6 ++++++
 untests/vmath_fun_tests.cpp         |  7 +++++++
 untests/vmath_vec_fun_tests.cpp     |  7 +++++++
 5 files changed, 39 insertions(+)

diff --git a/README.md b/README.md
index e795a55..a87b95c 100644
--- a/README.md
+++ b/README.md
@@ -617,6 +617,9 @@ T acosh(T x) noexcept;
 template < floating_point T >
 T atanh(T x) noexcept;
 
+template < floating_point T >
+void sincos(T x, T* s, T* c) noexcept;
+
 // Vector
 
 template < typename T, size_t Size >
@@ -663,6 +666,9 @@ vec<T, Size> acosh(const vec<T, Size>& xs);
 
 template < typename T, size_t Size >
 vec<T, Size> atanh(const vec<T, Size>& xs);
+
+template < typename T, size_t Size >
+void sincos(const vec<T, Size>& xs, vec<T, Size>* ss, vec<T, Size>* cs);
 ```
 
 ### Exponential Functions
diff --git a/headers/vmath.hpp/vmath_fun.hpp b/headers/vmath.hpp/vmath_fun.hpp
index c753df6..ab0caea 100644
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -103,6 +103,19 @@ namespace vmath_hpp
     atanh(T x) noexcept {
         return std::atanh(x);
     }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, std::pair<T, T>>
+    sincos(T x) noexcept {
+        return {sin(x), cos(x)};
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, void>
+    sincos(T x, T* s, T* c) noexcept {
+        *s = sin(x);
+        *c = cos(x);
+    }
 }
 
 //
diff --git a/headers/vmath.hpp/vmath_vec_fun.hpp b/headers/vmath.hpp/vmath_vec_fun.hpp
index ffbff21..ce780dc 100644
--- a/headers/vmath.hpp/vmath_vec_fun.hpp
+++ b/headers/vmath.hpp/vmath_vec_fun.hpp
@@ -400,6 +400,12 @@ namespace vmath_hpp
     [[nodiscard]] vec<T, Size> atanh(const vec<T, Size>& xs) {
         return map_join([](T x) { return atanh(x); }, xs);
     }
+
+    template < typename T, size_t Size >
+    void sincos(const vec<T, Size>& xs, vec<T, Size>* ss, vec<T, Size>* cs) {
+        *ss = map_join([](T x){ return sin(x); }, xs);
+        *cs = map_join([](T x){ return cos(x); }, xs);
+    }
 }
 
 //
diff --git a/untests/vmath_fun_tests.cpp b/untests/vmath_fun_tests.cpp
index 11eb4af..be74efd 100644
--- a/untests/vmath_fun_tests.cpp
+++ b/untests/vmath_fun_tests.cpp
@@ -34,6 +34,13 @@ TEST_CASE("vmath/fun") {
         (void)asinh(0.f);
         (void)acosh(0.f);
         (void)atanh(0.f);
+
+        {
+            float out_s{}, out_c{};
+            sincos(15.f, &out_s, &out_c);
+            REQUIRE(out_s == approx(sin(15.f)));
+            REQUIRE(out_c == approx(cos(15.f)));
+        }
     }
 
     SUBCASE("Exponential Functions") {
diff --git a/untests/vmath_vec_fun_tests.cpp b/untests/vmath_vec_fun_tests.cpp
index 67faff9..f5f256a 100644
--- a/untests/vmath_vec_fun_tests.cpp
+++ b/untests/vmath_vec_fun_tests.cpp
@@ -108,6 +108,13 @@ TEST_CASE("vmath/vec_fun") {
         (void)asinh(float2(1.f));
         (void)acosh(float2(1.f));
         (void)atanh(float2(1.f));
+
+        {
+            float2 out_ss{}, out_cs{};
+            sincos(float2(10.f,15.f), &out_ss, &out_cs);
+            REQUIRE(out_ss == approx2(sin(10.f), sin(15.f)));
+            REQUIRE(out_cs == approx2(cos(10.f), cos(15.f)));
+        }
     }
 
     SUBCASE("Exponential Functions") {

From 422d251df1d2e009ab00dcaeae3c3774a2f9f975 Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sun, 6 Dec 2020 04:49:18 +0700
Subject: [PATCH 4/6] normalize axis in rotate function, little style fixes

---
 headers/vmath.hpp/vmath_ext.hpp | 136 +++++++-------
 headers/vmath.hpp/vmath_fun.hpp | 306 ++++++++++++++++----------------
 2 files changed, 221 insertions(+), 221 deletions(-)

diff --git a/headers/vmath.hpp/vmath_ext.hpp b/headers/vmath.hpp/vmath_ext.hpp
index 435d2e5..67c82b2 100644
--- a/headers/vmath.hpp/vmath_ext.hpp
+++ b/headers/vmath.hpp/vmath_ext.hpp
@@ -213,26 +213,27 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> rotate(T angle, const vec<T, 3>& axis) {
-        const T x = axis.x;
-        const T y = axis.y;
-        const T z = axis.z;
-        const T px = x * x;
-        const T py = y * y;
-        const T pz = z * z;
-        const T cs = cos(angle);
-        const T sn = sin(angle);
-        const T ics = T(1) - cs;
-        const T xym = x * y * ics;
-        const T xzm = x * z * ics;
-        const T yzm = y * z * ics;
-        const T xsn = x * sn;
-        const T ysn = y * sn;
-        const T zsn = z * sn;
+        const auto [s, c] = sincos(angle);
+        const auto [x, y, z] = normalize(axis);
+
+        const T xx = x * x;
+        const T yy = y * y;
+        const T zz = z * z;
+
+        const T xs = x * s;
+        const T ys = y * s;
+        const T zs = z * s;
+
+        const T ic = T(1) - c;
+        const T xym = x * y * ic;
+        const T xzm = x * z * ic;
+        const T yzm = y * z * ic;
+
         return {
-            px * ics + cs, xym + zsn,     xzm - ysn,     0,
-            xym - zsn,     py * ics + cs, yzm + xsn,     0,
-            xzm + ysn,     yzm - xsn,     pz * ics + cs, 0,
-            0,             0,             0,             1};
+            xx * ic + c, xym + zs,    xzm - ys,    0,
+            xym - zs,    yy * ic + c, yzm + xs,    0,
+            xzm + ys,    yzm - xs,    zz * ic + c, 0,
+            0,           0,           0,           1};
     }
 
     template < typename T >
@@ -242,13 +243,12 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> rotate_x(T angle) {
-        const T cs = cos(angle);
-        const T sn = sin(angle);
+        const auto [s, c] = sincos(angle);
         return {
-            1, 0,   0,  0,
-            0, cs,  sn, 0,
-            0, -sn, cs, 0,
-            0, 0,   0,  1};
+            1,  0, 0, 0,
+            0,  c, s, 0,
+            0, -s, c, 0,
+            0,  0, 0, 1};
     }
 
     template < typename T >
@@ -258,13 +258,12 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> rotate_y(T angle) {
-        const T cs = cos(angle);
-        const T sn = sin(angle);
+        const auto [s, c] = sincos(angle);
         return {
-            cs, 0, -sn, 0,
-            0,  1, 0,   0,
-            sn, 0, cs,  0,
-            0,  0, 0,   1};
+            c, 0, -s, 0,
+            0, 1,  0, 0,
+            s, 0,  c, 0,
+            0, 0,  0, 1};
     }
 
     template < typename T >
@@ -274,13 +273,12 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> rotate_z(T angle) {
-        const T cs = cos(angle);
-        const T sn = sin(angle);
+        const auto [s, c] = sincos(angle);
         return {
-            cs,  sn, 0, 0,
-            -sn, cs, 0, 0,
-            0,   0,  1, 0,
-            0,   0,  0, 1};
+             c, s, 0, 0,
+            -s, c, 0, 0,
+             0, 0, 1, 0,
+             0, 0, 0, 1};
     }
 
     template < typename T >
@@ -321,9 +319,11 @@ namespace vmath_hpp
         const vec<T, 3> az = normalize(at - eye);
         const vec<T, 3> ax = normalize(cross(up, az));
         const vec<T, 3> ay = cross(az, ax);
+
         const T dx = dot(ax, eye);
         const T dy = dot(ay, eye);
         const T dz = dot(az, eye);
+
         return {
             ax.x, ay.x, az.x, 0,
             ax.y, ay.y, az.y, 0,
@@ -336,9 +336,11 @@ namespace vmath_hpp
         const vec<T, 3> az = normalize(eye - at);
         const vec<T, 3> ax = normalize(cross(up, az));
         const vec<T, 3> ay = cross(az, ax);
+
         const T dx = dot(ax, eye);
         const T dy = dot(ay, eye);
         const T dz = dot(az, eye);
+
         return {
             ax.x, ay.x, az.x, 0,
             ax.y, ay.y, az.y, 0,
@@ -382,12 +384,11 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 3> rotate(T angle) {
-        const T cs = cos(angle);
-        const T sn = sin(angle);
+        const auto [s, c] = sincos(angle);
         return {
-            cs,  sn, 0,
-            -sn, cs, 0,
-            0,   0,  1};
+             c, s, 0,
+            -s, c, 0,
+             0, 0, 1};
     }
 
     template < typename T >
@@ -482,9 +483,9 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> orthographic_lh_zo(T left, T right, T bottom, T top, T znear, T zfar) {
-        const T sx = T(2) / (right - left);
-        const T sy = T(2) / (top - bottom);
-        const T sz = T(1) / (zfar - znear);
+        const T sx = T(2) * reciprocal(right - left);
+        const T sy = T(2) * reciprocal(top - bottom);
+        const T sz = T(1) * reciprocal(zfar - znear);
 
         const T tx = - (right + left) / (right - left);
         const T ty = - (top + bottom) / (top - bottom);
@@ -499,9 +500,9 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> orthographic_lh_no(T left, T right, T bottom, T top, T znear, T zfar) {
-        const T sx = T(2) / (right - left);
-        const T sy = T(2) / (top - bottom);
-        const T sz = T(2) / (zfar - znear);
+        const T sx = T(2) * reciprocal(right - left);
+        const T sy = T(2) * reciprocal(top - bottom);
+        const T sz = T(2) * reciprocal(zfar - znear);
 
         const T tx = - (right + left) / (right - left);
         const T ty = - (top + bottom) / (top - bottom);
@@ -516,9 +517,9 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> orthographic_rh_zo(T left, T right, T bottom, T top, T znear, T zfar) {
-        const T sx = T(2) / (right - left);
-        const T sy = T(2) / (top - bottom);
-        const T sz = -T(1) / (zfar - znear);
+        const T sx = T(2) * reciprocal(right - left);
+        const T sy = T(2) * reciprocal(top - bottom);
+        const T sz = -T(1) * reciprocal(zfar - znear);
 
         const T tx = - (right + left) / (right - left);
         const T ty = - (top + bottom) / (top - bottom);
@@ -533,9 +534,9 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> orthographic_rh_no(T left, T right, T bottom, T top, T znear, T zfar) {
-        const T sx = T(2) / (right - left);
-        const T sy = T(2) / (top - bottom);
-        const T sz = -T(2) / (zfar - znear);
+        const T sx = T(2) * reciprocal(right - left);
+        const T sy = T(2) * reciprocal(top - bottom);
+        const T sz = -T(2) * reciprocal(zfar - znear);
 
         const T tx = - (right + left) / (right - left);
         const T ty = - (top + bottom) / (top - bottom);
@@ -552,7 +553,7 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> perspective_lh_zo(T fov, T aspect, T znear, T zfar) {
-        const T sy = T(1) / tan(fov * T(0.5));
+        const T sy = reciprocal(tan(fov * T(0.5)));
         const T sx = sy / aspect;
         const T sz = zfar / (zfar - znear);
         const T tz = (znear * zfar) / (znear - zfar);
@@ -565,7 +566,7 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> perspective_lh_no(T fov, T aspect, T znear, T zfar) {
-        const T sy = T(1) / tan(fov * T(0.5));
+        const T sy = reciprocal(tan(fov * T(0.5)));
         const T sx = sy / aspect;
         const T sz = (zfar + znear) / (zfar - znear);
         const T tz = (T(2) * znear * zfar) / (znear - zfar);
@@ -578,28 +579,28 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] mat<T, 4> perspective_rh_zo(T fov, T aspect, T znear, T zfar) {
-        const T sy = T(1) / tan(fov * T(0.5));
+        const T sy = reciprocal(tan(fov * T(0.5)));
         const T sx = sy / aspect;
         const T sz = zfar / (znear - zfar);
         const T tz = (znear * zfar) / (znear - zfar);
         return  {
-            sx, 0,  0,  0,
-            0,  sy, 0,  0,
+            sx, 0,  0,   0,
+            0,  sy, 0,   0,
             0,  0,  sz, -1,
-            0,  0,  tz, 0};
+            0,  0,  tz,  0};
     }
 
     template < typename T >
     [[nodiscard]] mat<T, 4> perspective_rh_no(T fov, T aspect, T znear, T zfar) {
-        const T sy = T(1) / tan(fov * T(0.5));
+        const T sy = reciprocal(tan(fov * T(0.5)));
         const T sx = sy / aspect;
         const T sz = (zfar + znear) / (znear - zfar);
         const T tz = (T(2) * znear * zfar) / (znear - zfar);
         return  {
-            sx, 0,  0,  0,
-            0,  sy, 0,  0,
+            sx, 0,  0,   0,
+            0,  sy, 0,   0,
             0,  0,  sz, -1,
-            0,  0,  tz, 0};
+            0,  0,  tz,  0};
     }
 }
 
@@ -620,11 +621,10 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] vec<T, 2> rotate(const vec<T, 2>& v, T angle) {
-        const T cs = cos(angle);
-        const T sn = sin(angle);
+        const auto [s, c] = sincos(angle);
         return {
-            v.x * cs - v.y * sn,
-            v.x * sn + v.y * cs};
+            v.x * c - v.y * s,
+            v.x * s + v.y * c};
     }
 
     template < typename T >
diff --git a/headers/vmath.hpp/vmath_fun.hpp b/headers/vmath.hpp/vmath_fun.hpp
index ab0caea..7a201a1 100644
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -8,6 +8,158 @@
 
 #include "vmath_fwd.hpp"
 
+//
+// Common Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_unsigned_v<T>, T>
+    constexpr abs(T x) noexcept {
+        return x;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_signed_v<T>, T>
+    constexpr abs(T x) noexcept {
+        return x >= T(0) ? x : -x;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr sign(T x) noexcept {
+        return static_cast<T>((T(0) < x) - (x < T(0)));
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr reciprocal(T x) noexcept {
+        return T(1) / x;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    floor(T x) noexcept {
+        return std::floor(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    trunc(T x) noexcept {
+        return std::trunc(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    round(T x) noexcept {
+        return std::round(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    ceil(T x) noexcept {
+        return std::ceil(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    fract(T x) noexcept {
+        return x - floor(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    fmod(T x, T y) noexcept {
+        return std::fmod(x, y);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    modf(T x, T* y) noexcept {
+        return std::modf(x, y);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr min(T x, T y) noexcept {
+        return x < y ? x : y;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr max(T x, T y) noexcept {
+        return x < y ? y : x;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr clamp(T x, T min_x, T max_x) noexcept {
+        return min(max(x, min_x), max_x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr saturate(T x) noexcept {
+        return clamp(x, T(0), T(1));
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr lerp(T x, T y, T a) noexcept {
+        return x * (T(1) - a) + y * a;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr step(T edge, T x) noexcept {
+        return x < edge ? T(0) : T(1);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr smoothstep(T edge0, T edge1, T x) noexcept {
+        const T t = clamp((x - edge0) / (edge1 - edge0), T(0), T(1));
+        return t * t * (T(3) - T(2) * t);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    isnan(T x) noexcept {
+        return std::isnan(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    isinf(T x) noexcept {
+        return std::isinf(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    isfinite(T x) noexcept {
+        return std::isfinite(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    fma(T x, T y, T z) noexcept {
+        return std::fma(x, y, z);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    frexp(T x, int* exp) noexcept {
+        return std::frexp(x, exp);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    ldexp(T x, int exp) noexcept {
+        return std::ldexp(x, exp);
+    }
+}
+
 //
 // Angle and Trigonometry Functions
 //
@@ -163,159 +315,7 @@ namespace vmath_hpp
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
     rsqrt(T x) noexcept {
-        return T(1) / sqrt(x);
-    }
-}
-
-//
-// Common Functions
-//
-
-namespace vmath_hpp
-{
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_unsigned_v<T>, T>
-    constexpr abs(T x) noexcept {
-        return x;
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_signed_v<T>, T>
-    constexpr abs(T x) noexcept {
-        return x >= T(0) ? x : -x;
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
-    constexpr sign(T x) noexcept {
-        return static_cast<T>((T(0) < x) - (x < T(0)));
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    constexpr reciprocal(T x) noexcept {
-        return T(1) / x;
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    floor(T x) noexcept {
-        return std::floor(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    trunc(T x) noexcept {
-        return std::trunc(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    round(T x) noexcept {
-        return std::round(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    ceil(T x) noexcept {
-        return std::ceil(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    fract(T x) noexcept {
-        return x - floor(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    fmod(T x, T y) noexcept {
-        return std::fmod(x, y);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    modf(T x, T* y) noexcept {
-        return std::modf(x, y);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
-    constexpr min(T x, T y) noexcept {
-        return x < y ? x : y;
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
-    constexpr max(T x, T y) noexcept {
-        return x < y ? y : x;
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
-    constexpr clamp(T x, T min_x, T max_x) noexcept {
-        return min(max(x, min_x), max_x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
-    constexpr saturate(T x) noexcept {
-        return clamp(x, T(0), T(1));
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    constexpr lerp(T x, T y, T a) noexcept {
-        return x * (T(1) - a) + y * a;
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    constexpr step(T edge, T x) noexcept {
-        return x < edge ? T(0) : T(1);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    constexpr smoothstep(T edge0, T edge1, T x) noexcept {
-        const T t = clamp((x - edge0) / (edge1 - edge0), T(0), T(1));
-        return t * t * (T(3) - T(2) * t);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
-    isnan(T x) noexcept {
-        return std::isnan(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
-    isinf(T x) noexcept {
-        return std::isinf(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
-    isfinite(T x) noexcept {
-        return std::isfinite(x);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    fma(T x, T y, T z) noexcept {
-        return std::fma(x, y, z);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    frexp(T x, int* exp) noexcept {
-        return std::frexp(x, exp);
-    }
-
-    template < typename T >
-    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    ldexp(T x, int exp) noexcept {
-        return std::ldexp(x, exp);
+        return reciprocal(sqrt(x));
     }
 }
 

From b0659e92393130cb1add250f405a2fac6b585510 Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sun, 6 Dec 2020 05:38:00 +0700
Subject: [PATCH 5/6] variadic min/max

---
 README.md                       |  9 +++++++++
 headers/vmath.hpp/vmath_fun.hpp | 16 ++++++++++++++++
 untests/vmath_fun_tests.cpp     |  9 +++++++--
 3 files changed, 32 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index a87b95c..cbb55a8 100644
--- a/README.md
+++ b/README.md
@@ -617,6 +617,9 @@ T acosh(T x) noexcept;
 template < floating_point T >
 T atanh(T x) noexcept;
 
+template < floating_point T >
+std::pair<T, T> sincos(T x) noexcept;
+
 template < floating_point T >
 void sincos(T x, T* s, T* c) noexcept;
 
@@ -759,9 +762,15 @@ T modf(T x, T* y) noexcept;
 template < arithmetic T >
 constexpr T min(T x, T y) noexcept;
 
+template < arithmetic T, arithmetic... Ts >
+constexpr std::common_type_t<T, Ts...> min(T x, T y, Ts... ts) noexcept;
+
 template < arithmetic T >
 constexpr T max(T x, T y) noexcept;
 
+template < arithmetic T, arithmetic... Ts >
+constexpr std::common_type_t<T, Ts...> max(T x, T y, Ts... ts) noexcept;
+
 template < arithmetic T >
 constexpr T clamp(T x, T min_x, T max_x) noexcept;
 
diff --git a/headers/vmath.hpp/vmath_fun.hpp b/headers/vmath.hpp/vmath_fun.hpp
index 7a201a1..e2d0bc1 100644
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -86,12 +86,28 @@ namespace vmath_hpp
         return x < y ? x : y;
     }
 
+    template < typename T, typename... Ts >
+    [[nodiscard]] std::enable_if_t<
+        std::is_arithmetic_v<T>,
+        std::common_type_t<T, Ts...>>
+    constexpr min(T x, T y, Ts... ts) noexcept {
+        return min(min(x, y), ts...);
+    }
+
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
     constexpr max(T x, T y) noexcept {
         return x < y ? y : x;
     }
 
+    template < typename T, typename... Ts >
+    [[nodiscard]] std::enable_if_t<
+        std::is_arithmetic_v<T>,
+        std::common_type_t<T, Ts...>>
+    constexpr max(T x, T y, Ts... ts) noexcept {
+        return max(max(x, y), ts...);
+    }
+
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
     constexpr clamp(T x, T min_x, T max_x) noexcept {
diff --git a/untests/vmath_fun_tests.cpp b/untests/vmath_fun_tests.cpp
index be74efd..2cfcb89 100644
--- a/untests/vmath_fun_tests.cpp
+++ b/untests/vmath_fun_tests.cpp
@@ -85,8 +85,13 @@ TEST_CASE("vmath/fun") {
             REQUIRE(out_i == approx(1.f));
         }
 
-        STATIC_REQUIRE(min(1.f, 2.f) == approx(1.f));
-        STATIC_REQUIRE(max(1.f, 2.f) == approx(2.f));
+        STATIC_REQUIRE(min(0.f, 1.f) == approx(0.f));
+        STATIC_REQUIRE(min(3.f, 2.f, 1.f) == approx(1.f));
+        STATIC_REQUIRE(min(4.f, 3.f, 2.f, 1.f) == approx(1.f));
+
+        STATIC_REQUIRE(max(0.f, 1.f) == approx(1.f));
+        STATIC_REQUIRE(max(3.f, 2.f, 1.f) == approx(3.f));
+        STATIC_REQUIRE(max(4.f, 3.f, 2.f, 1.f) == approx(4.f));
 
         STATIC_REQUIRE(clamp(1.0f, 2.f, 3.f) == approx(2.0f));
         STATIC_REQUIRE(clamp(2.5f, 2.f, 3.f) == approx(2.5f));

From 5fe1aab5b6182150b28b01137ae0503d6496d7d8 Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sun, 6 Dec 2020 23:15:08 +0700
Subject: [PATCH 6/6] equal floats with tolerance

---
 headers/vmath.hpp/vmath_fun.hpp     | 19 +++++++++++++++----
 headers/vmath.hpp/vmath_fwd.hpp     |  1 +
 headers/vmath.hpp/vmath_mat_fun.hpp |  4 +---
 headers/vmath.hpp/vmath_vec_fun.hpp |  6 ++----
 untests/vmath_fun_tests.cpp         | 10 ++++++++++
 5 files changed, 29 insertions(+), 11 deletions(-)

diff --git a/headers/vmath.hpp/vmath_fun.hpp b/headers/vmath.hpp/vmath_fun.hpp
index e2d0bc1..a876bd8 100644
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -431,25 +431,36 @@ namespace vmath_hpp
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr equal_to(T x, T y) noexcept {
-        return x == y;
+        if constexpr ( std::is_floating_point_v<T> ) {
+            // http://www.realtimecollisiondetection.net/pubs/Tolerances
+            const T epsilon = std::numeric_limits<T>::epsilon();
+            return abs(x - y) <= epsilon * max(T(1), abs(x), abs(y));
+        } else {
+            return x == y;
+        }
     }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr equal_to(T x, T y, T epsilon) noexcept {
-        return abs(x - y) <= epsilon;
+        if constexpr ( std::is_floating_point_v<T> ) {
+            // http://www.realtimecollisiondetection.net/pubs/Tolerances
+            return abs(x - y) <= epsilon * max(T(1), abs(x), abs(y));
+        } else {
+            return abs(x - y) <= epsilon;
+        }
     }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr not_equal_to(T x, T y) noexcept {
-        return x != y;
+        return !equal_to(x, y);
     }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr not_equal_to(T x, T y, T epsilon) noexcept {
-        return abs(x - y) > epsilon;
+        return !equal_to(x, y, epsilon);
     }
 
     template < typename T >
diff --git a/headers/vmath.hpp/vmath_fwd.hpp b/headers/vmath.hpp/vmath_fwd.hpp
index 2b330ec..b7a247a 100644
--- a/headers/vmath.hpp/vmath_fwd.hpp
+++ b/headers/vmath.hpp/vmath_fwd.hpp
@@ -11,6 +11,7 @@
 
 #include <initializer_list>
 #include <iterator>
+#include <limits>
 #include <stdexcept>
 #include <type_traits>
 #include <utility>
diff --git a/headers/vmath.hpp/vmath_mat_fun.hpp b/headers/vmath.hpp/vmath_mat_fun.hpp
index 5960e2e..8ffc74e 100644
--- a/headers/vmath.hpp/vmath_mat_fun.hpp
+++ b/headers/vmath.hpp/vmath_mat_fun.hpp
@@ -340,9 +340,7 @@ namespace vmath_hpp
 
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr bool operator!=(const mat<T, Size>& xs, const mat<T, Size>& ys) {
-        return fold_join([](bool acc, const vec<T, Size>& x, const vec<T, Size>& y){
-            return acc || (x != y);
-        }, false, xs, ys);
+        return !(xs == ys);
     }
 
     // operator<
diff --git a/headers/vmath.hpp/vmath_vec_fun.hpp b/headers/vmath.hpp/vmath_vec_fun.hpp
index ce780dc..d9a0cf0 100644
--- a/headers/vmath.hpp/vmath_vec_fun.hpp
+++ b/headers/vmath.hpp/vmath_vec_fun.hpp
@@ -293,15 +293,13 @@ namespace vmath_hpp
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr bool operator==(const vec<T, Size>& xs, const vec<T, Size>& ys) {
         return fold_join([](bool acc, T x, T y){
-            return acc && (x == y);
+            return acc && equal_to(x, y);
         }, true, xs, ys);
     }
 
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr bool operator!=(const vec<T, Size>& xs, const vec<T, Size>& ys) {
-        return fold_join([](bool acc, T x, T y){
-            return acc || (x != y);
-        }, false, xs, ys);
+        return !(xs == ys);
     }
 
     // operator<
diff --git a/untests/vmath_fun_tests.cpp b/untests/vmath_fun_tests.cpp
index 2cfcb89..33787f5 100644
--- a/untests/vmath_fun_tests.cpp
+++ b/untests/vmath_fun_tests.cpp
@@ -165,6 +165,16 @@ TEST_CASE("vmath/fun") {
         STATIC_REQUIRE_FALSE(not_equal_to(1, 1, 0));
         STATIC_REQUIRE_FALSE(not_equal_to(1, 1, 1));
 
+        STATIC_REQUIRE(equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 0.5f));
+        STATIC_REQUIRE_FALSE(equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 1.5f));
+        STATIC_REQUIRE(equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 90.f));
+        STATIC_REQUIRE_FALSE(equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 110.f));
+
+        STATIC_REQUIRE_FALSE(not_equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 0.5f));
+        STATIC_REQUIRE(not_equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 1.5f));
+        STATIC_REQUIRE_FALSE(not_equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 90.f));
+        STATIC_REQUIRE(not_equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 110.f));
+
         STATIC_REQUIRE_FALSE(any(false));
         STATIC_REQUIRE_FALSE(any(0));
         STATIC_REQUIRE(any(true));