Merge pull request #5 from BlackMATov/dev

Dev
2025-12-14 20:31:25 +07:00 · 2020-12-06 23:32:14 +07:00
parent 2917a19cbc 5fe1aab5b6
commit 942e074c63
9 changed files with 393 additions and 216 deletions
--- a/README.md
+++ b/README.md
@@ -617,6 +617,12 @@ 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;
+
 // Vector

 template < typename T, size_t Size >
@@ -663,6 +669,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
@@ -753,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;

@@ -1147,6 +1162,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);

@@ -1268,6 +1301,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)
--- 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 >
@@ -240,6 +241,51 @@ namespace vmath_hpp
        return m * rotate(angle, axis);
    }

+    template < typename T >
+    [[nodiscard]] mat<T, 4> rotate_x(T angle) {
+        const auto [s, c] = sincos(angle);
+        return {
+            1,  0, 0, 0,
+            0,  c, s, 0,
+            0, -s, c, 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 auto [s, c] = sincos(angle);
+        return {
+            c, 0, -s, 0,
+            0, 1,  0, 0,
+            s, 0,  c, 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 auto [s, c] = sincos(angle);
+        return {
+             c, s, 0, 0,
+            -s, c, 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 >
@@ -273,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,
@@ -288,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,
@@ -334,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 >
@@ -434,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);
@@ -451,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);
@@ -468,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);
@@ -485,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);
@@ -504,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);
@@ -517,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);
@@ -530,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};
    }
 }

@@ -572,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 >
@@ -588,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;
+    }
 }
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -8,6 +8,174 @@

 #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, 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 {
+        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
 //
@@ -103,6 +271,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);
+    }
 }

 //
@@ -150,159 +331,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));
    }
 }

@@ -402,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 >
--- 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>
--- 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<
--- 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<
@@ -400,6 +398,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);
+    }
 }

 //
--- 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));
@@ -212,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));
+    }
 }
--- 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") {
@@ -78,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));
@@ -153,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));
--- 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") {