diff --git a/headers/vmath.hpp/vmath_ext.hpp b/headers/vmath.hpp/vmath_ext.hpp
index 19fd3cc..ff011fe 100644
--- a/headers/vmath.hpp/vmath_ext.hpp
+++ b/headers/vmath.hpp/vmath_ext.hpp
@@ -866,55 +866,17 @@ namespace vmath_hpp
         /// REFERENCE:
         /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
 
-        const T wv = m[0][0] + m[1][1] + m[2][2];
-        const T xv = m[0][0] - m[1][1] - m[2][2];
-        const T yv = m[1][1] - m[0][0] - m[2][2];
-        const T zv = m[2][2] - m[0][0] - m[1][1];
+        auto xyzw = T(0.5) * sqrt(max(zero4<T>, {
+            T(1) + m[0][0] - m[1][1] - m[2][2],
+            T(1) - m[0][0] + m[1][1] - m[2][2],
+            T(1) - m[0][0] - m[1][1] + m[2][2],
+            T(1) + m[0][0] + m[1][1] + m[2][2]}));
 
-        struct pmv {
-            const T v{};
-            const int i{};
-        };
-
-        const auto [mv,mi] = max<pmv>({
-            {wv, 0},
-            {xv, 1},
-            {yv, 2},
-            {zv, 3},
-        }, [](auto&& l, auto&& r){
-            return l.v < r.v;
-        });
-
-        const T qv = T(0.5) * sqrt(T(1) + mv);
-        const T rqv = T(0.25) * rcp(qv);
-
-        switch ( mi ) {
-        default:
-        case 0:
-            return {
-                (m[1][2] - m[2][1]) * rqv,
-                (m[2][0] - m[0][2]) * rqv,
-                (m[0][1] - m[1][0]) * rqv,
-                qv};
-        case 1:
-            return {
-                qv,
-                (m[1][0] + m[0][1]) * rqv,
-                (m[2][0] + m[0][2]) * rqv,
-                (m[1][2] - m[2][1]) * rqv};
-        case 2:
-            return {
-                (m[1][0] + m[0][1]) * rqv,
-                qv,
-                (m[2][1] + m[1][2]) * rqv,
-                (m[2][0] - m[0][2]) * rqv};
-        case 3:
-            return {
-                (m[2][0] + m[0][2]) * rqv,
-                (m[2][1] + m[1][2]) * rqv,
-                qv,
-                (m[0][1] - m[1][0]) * rqv};
-        }
+        return qua(copysign(xyzw, {
+            m[1][2] - m[2][1],
+            m[2][0] - m[0][2],
+            m[0][1] - m[1][0],
+            T(1)}));
     }
 
     template < typename T >
diff --git a/headers/vmath.hpp/vmath_fun.hpp b/headers/vmath.hpp/vmath_fun.hpp
index 0fae129..56276b6 100644
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -14,18 +14,18 @@
 
 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_unsigned_v<T>, T>
+    constexpr abs(T x) noexcept {
+        return x;
+    }
+
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
     constexpr sqr(T x) noexcept {
@@ -86,66 +86,24 @@ namespace vmath_hpp
         return std::modf(x, y);
     }
 
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    copysign(T x, T s) noexcept {
+        return std::copysign(x, s);
+    }
+
     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]] constexpr T min(std::initializer_list<T> xs) {
-        auto iter = xs.begin();
-        auto smallest = iter++;
-        for ( auto last = xs.end(); iter != last; ++iter ) {
-            if ( *iter < *smallest ) {
-                smallest = iter;
-            }
-        }
-        return *smallest;
-    }
-
-    template < typename T, class Compare >
-    [[nodiscard]] constexpr T min(std::initializer_list<T> xs, Compare comp) {
-        auto iter = xs.begin();
-        auto smallest = iter++;
-        for ( auto last = xs.end(); iter != last; ++iter ) {
-            if ( comp(*iter, *smallest) ) {
-                smallest = iter;
-            }
-        }
-        return *smallest;
-    }
-
     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]] constexpr T max(std::initializer_list<T> xs) {
-        auto iter = xs.begin();
-        auto largest = iter++;
-        for ( auto last = xs.end(); iter != last; ++iter ) {
-            if ( *largest < *iter ) {
-                largest = iter;
-            }
-        }
-        return *largest;
-    }
-
-    template < typename T, class Compare >
-    [[nodiscard]] constexpr T max(std::initializer_list<T> xs, Compare comp) {
-        auto iter = xs.begin();
-        auto largest = iter++;
-        for ( auto last = xs.end(); iter != last; ++iter ) {
-            if ( comp(*largest, *iter) ) {
-                largest = iter;
-            }
-        }
-        return *largest;
-    }
-
     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 {
@@ -240,63 +198,94 @@ namespace vmath_hpp
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    sin(T x) noexcept { return std::sin(x); }
+    sin(T x) noexcept {
+        return std::sin(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    cos(T x) noexcept { return std::cos(x); }
+    cos(T x) noexcept {
+        return std::cos(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    tan(T x) noexcept { return std::tan(x); }
+    tan(T x) noexcept {
+        return std::tan(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    asin(T x) noexcept { return std::asin(x); }
+    asin(T x) noexcept {
+        return std::asin(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    acos(T x) noexcept { return std::acos(x); }
+    acos(T x) noexcept {
+        return std::acos(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    atan(T x) noexcept { return std::atan(x); }
+    atan(T x) noexcept {
+        return std::atan(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    atan2(T y, T x) noexcept { return std::atan2(y, x); }
+    atan2(T y, T x) noexcept {
+        return std::atan2(y, x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    sinh(T x) noexcept { return std::sinh(x); }
+    sinh(T x) noexcept {
+        return std::sinh(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    cosh(T x) noexcept { return std::cosh(x); }
+    cosh(T x) noexcept {
+        return std::cosh(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    tanh(T x) noexcept { return std::tanh(x); }
+    tanh(T x) noexcept {
+        return std::tanh(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    asinh(T x) noexcept { return std::asinh(x); }
+    asinh(T x) noexcept {
+        return std::asinh(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    acosh(T x) noexcept { return std::acosh(x); }
+    acosh(T x) noexcept {
+        return std::acosh(x);
+    }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
-    atanh(T x) noexcept { return std::atanh(x); }
+    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)}; }
+    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); }
+    sincos(T x, T* s, T* c) noexcept {
+        *s = sin(x);
+        *c = cos(x);
+    }
 }
 
 //
@@ -436,7 +425,7 @@ namespace vmath_hpp
             /// REFERENCE:
             /// 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)});
+            return abs(x - y) <= epsilon * max(max(T(1), abs(x)), abs(y));
         } else {
             return x == y;
         }
@@ -448,7 +437,7 @@ namespace vmath_hpp
         if constexpr ( std::is_floating_point_v<T> ) {
             /// REFERENCE:
             /// http://www.realtimecollisiondetection.net/pubs/Tolerances
-            return abs(x - y) <= epsilon * max({T(1), abs(x), abs(y)});
+            return abs(x - y) <= epsilon * max(max(T(1), abs(x)), abs(y));
         } else {
             return abs(x - y) <= epsilon;
         }
diff --git a/headers/vmath.hpp/vmath_vec_fun.hpp b/headers/vmath.hpp/vmath_vec_fun.hpp
index 3ab9274..95251d1 100644
--- a/headers/vmath.hpp/vmath_vec_fun.hpp
+++ b/headers/vmath.hpp/vmath_vec_fun.hpp
@@ -641,6 +641,11 @@ namespace vmath_hpp
         return map_join([](T x) { return abs(x); }, xs);
     }
 
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> sqr(const vec<T, Size>& xs) {
+        return map_join([](T x) { return sqr(x); }, xs);
+    }
+
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr vec<T, Size> sign(const vec<T, Size>& xs) {
         return map_join([](T x) { return sign(x); }, xs);
@@ -700,6 +705,16 @@ namespace vmath_hpp
         return impl::modf_impl(xs, is, std::make_index_sequence<Size>{});
     }
 
+    template < typename T, std::size_t Size >
+    [[nodiscard]] vec<T, Size> copysign(const vec<T, Size>& xs, T s) {
+        return map_join([s](T x) { return copysign(x, s); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] vec<T, Size> copysign(const vec<T, Size>& xs, const vec<T, Size>& ss) {
+        return map_join([](T x, T s) { return copysign(x, s); }, xs, ss);
+    }
+
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr T min(const vec<T, Size>& xs) {
         return fold1_join([](T acc, T x){ return min(acc, x); }, xs);
diff --git a/untests/vmath_fun_tests.cpp b/untests/vmath_fun_tests.cpp
index 6cfd9cc..dfd1616 100644
--- a/untests/vmath_fun_tests.cpp
+++ b/untests/vmath_fun_tests.cpp
@@ -86,31 +86,6 @@ TEST_CASE("vmath/fun") {
         STATIC_REQUIRE(min(0.f, 1.f) == uapprox(0.f));
         STATIC_REQUIRE(max(0.f, 1.f) == uapprox(1.f));
 
-        STATIC_REQUIRE(min({0.f}) == uapprox(0.f));
-        STATIC_REQUIRE(min({0.f, 1.f}) == uapprox(0.f));
-        STATIC_REQUIRE(min({3.f, 2.f, 1.f}) == uapprox(1.f));
-        STATIC_REQUIRE(min({4.f, 3.f, 2.f, 1.f}) == uapprox(1.f));
-
-        STATIC_REQUIRE(min({0.f}, [](auto x, auto y){return x < y;}) == uapprox(0.f));
-        STATIC_REQUIRE(min({0.f, 1.f}, [](auto x, auto y){return x < y;}) == uapprox(0.f));
-        STATIC_REQUIRE(min({3.f, 2.f, 1.f}, [](auto x, auto y){return x < y;}) == uapprox(1.f));
-        STATIC_REQUIRE(min({4.f, 3.f, 2.f, 1.f}, [](auto x, auto y){return x < y;}) == uapprox(1.f));
-
-        STATIC_REQUIRE(min({0.f}, [](auto x, auto y){return x > y;}) == uapprox(0.f));
-        STATIC_REQUIRE(min({0.f, 1.f}, [](auto x, auto y){return x > y;}) == uapprox(1.f));
-        STATIC_REQUIRE(min({3.f, 2.f, 1.f}, [](auto x, auto y){return x > y;}) == uapprox(3.f));
-        STATIC_REQUIRE(min({4.f, 3.f, 2.f, 1.f}, [](auto x, auto y){return x > y;}) == uapprox(4.f));
-
-        STATIC_REQUIRE(max({0.f}, [](auto x, auto y){return x < y;}) == uapprox(0.f));
-        STATIC_REQUIRE(max({0.f, 1.f}, [](auto x, auto y){return x < y;}) == uapprox(1.f));
-        STATIC_REQUIRE(max({3.f, 2.f, 1.f}, [](auto x, auto y){return x < y;}) == uapprox(3.f));
-        STATIC_REQUIRE(max({4.f, 3.f, 2.f, 1.f}, [](auto x, auto y){return x < y;}) == uapprox(4.f));
-
-        STATIC_REQUIRE(max({0.f}, [](auto x, auto y){return x > y;}) == uapprox(0.f));
-        STATIC_REQUIRE(max({0.f, 1.f}, [](auto x, auto y){return x > y;}) == uapprox(0.f));
-        STATIC_REQUIRE(max({3.f, 2.f, 1.f}, [](auto x, auto y){return x > y;}) == uapprox(1.f));
-        STATIC_REQUIRE(max({4.f, 3.f, 2.f, 1.f}, [](auto x, auto y){return x > y;}) == uapprox(1.f));
-
         STATIC_REQUIRE(clamp(1.0f, 2.f, 3.f) == uapprox(2.0f));
         STATIC_REQUIRE(clamp(2.5f, 2.f, 3.f) == uapprox(2.5f));
         STATIC_REQUIRE(clamp(3.5f, 2.f, 3.f) == uapprox(3.0f));
diff --git a/untests/vmath_vec_fun_tests.cpp b/untests/vmath_vec_fun_tests.cpp
index 8f67e1c..7efdfa9 100644
--- a/untests/vmath_vec_fun_tests.cpp
+++ b/untests/vmath_vec_fun_tests.cpp
@@ -173,9 +173,25 @@ TEST_CASE("vmath/vec_fun") {
 
     SUBCASE("Common Functions") {
         STATIC_REQUIRE(abs(float2(1.f, -1.f)) == uapprox2(1.f,1.f));
+        STATIC_REQUIRE(sqr(float2(2.f, -3.f)) == uapprox2(4.f,9.f));
         STATIC_REQUIRE(sign(float3(1.f, -1.f, 0.f)) == uapprox3(1.f,-1.f,0.f));
         STATIC_REQUIRE(rcp(float2(2.f, 4.f)) == uapprox2(0.5f,0.25f));
 
+        REQUIRE(copysign(
+            float4(2.f, -4.f, 2.f, -4.f),
+            5.f)
+        == uapprox4(2.f, 4.f, 2.f, 4.f));
+
+        REQUIRE(copysign(
+            float4(2.f, -4.f, 2.f, -4.f),
+            -5.f)
+        == uapprox4(-2.f, -4.f, -2.f, -4.f));
+
+        REQUIRE(copysign(
+            float4(2.f, -4.f, 2.f, -4.f),
+            float4(10.f, 5.f, -4.f, -0.4))
+        == uapprox4(2.f, 4.f, -2.f, -4.f));
+
         (void)floor(float2(1.f, -1.f));
         (void)trunc(float2(1.f, -1.f));
         (void)round(float2(1.f, -1.f));