From 05cc7f0bf2f6832f08e963a769e28e90be0936dc Mon Sep 17 00:00:00 2001
From: BlackMATov <blackmat@mail.ru>
Date: Sat, 27 Feb 2021 08:03:31 +0700
Subject: [PATCH] add single header vmath variant

---
 CMakeLists.txt              |    4 +
 singles/vmath.hpp/vmath.hpp | 4403 +++++++++++++++++++++++++++++++++++
 untests/CMakeLists.txt      |    2 +-
 3 files changed, 4408 insertions(+), 1 deletion(-)
 create mode 100644 singles/vmath.hpp/vmath.hpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 542fef8..6d64e76 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -11,6 +11,10 @@ add_library(${PROJECT_NAME} INTERFACE)
 target_compile_features(${PROJECT_NAME} INTERFACE cxx_std_17)
 target_include_directories(${PROJECT_NAME} INTERFACE headers)
 
+add_library(${PROJECT_NAME}.singles INTERFACE)
+target_compile_features(${PROJECT_NAME}.singles INTERFACE cxx_std_17)
+target_include_directories(${PROJECT_NAME}.singles INTERFACE singles)
+
 if(BUILD_AS_STANDALONE)
     option(BUILD_WITH_UNTESTS "Build with unit tests" ON)
     if(BUILD_WITH_UNTESTS)
diff --git a/singles/vmath.hpp/vmath.hpp b/singles/vmath.hpp/vmath.hpp
new file mode 100644
index 0000000..b85f903
--- /dev/null
+++ b/singles/vmath.hpp/vmath.hpp
@@ -0,0 +1,4403 @@
+/*******************************************************************************
+ * This file is part of the "https://github.com/blackmatov/vmath.hpp"
+ * For conditions of distribution and use, see copyright notice in LICENSE.md
+ * Copyright (C) 2020-2021, by Matvey Cherevko (blackmatov@gmail.com)
+ ******************************************************************************/
+
+#pragma once
+
+#include <cmath>
+#include <cstddef>
+#include <cstdlib>
+
+#include <initializer_list>
+#include <iterator>
+#include <limits>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+
+#if defined(_MSC_VER)
+#  define VMATH_HPP_FORCE_INLINE __forceinline
+#elif defined(__clang__) || defined(__GNUC__)
+#  define VMATH_HPP_FORCE_INLINE inline __attribute__((__always_inline__))
+#else
+#  define VMATH_HPP_FORCE_INLINE inline
+#endif
+
+#if !defined(__cpp_exceptions) && !defined(__EXCEPTIONS) && !defined(_CPPUNWIND)
+#  define VMATH_HPP_NO_EXCEPTIONS
+#endif
+
+#ifdef VMATH_HPP_NO_EXCEPTIONS
+#  define VMATH_HPP_THROW(...) std::abort()
+#else
+#  define VMATH_HPP_THROW(...) throw __VA_ARGS__
+#endif
+
+#define VMATH_HPP_THROW_IF(pred, ...)\
+    ( (pred) ? VMATH_HPP_THROW(__VA_ARGS__) : (void)0 )
+
+namespace vmath_hpp
+{
+    struct uninit_t { explicit uninit_t() = default; };
+    inline constexpr uninit_t uninit{};
+
+    struct zero_init_t { explicit zero_init_t() = default; };
+    inline constexpr zero_init_t zero_init{};
+
+    struct unit_init_t { explicit unit_init_t() = default; };
+    inline constexpr unit_init_t unit_init{};
+
+    struct identity_init_t { explicit identity_init_t() = default; };
+    inline constexpr identity_init_t identity_init{};
+}
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    class vec;
+
+    using bvec2 = vec<bool, 2>;
+    using bvec3 = vec<bool, 3>;
+    using bvec4 = vec<bool, 4>;
+
+    using ivec2 = vec<int, 2>;
+    using ivec3 = vec<int, 3>;
+    using ivec4 = vec<int, 4>;
+
+    using uvec2 = vec<unsigned, 2>;
+    using uvec3 = vec<unsigned, 3>;
+    using uvec4 = vec<unsigned, 4>;
+
+    using fvec2 = vec<float, 2>;
+    using fvec3 = vec<float, 3>;
+    using fvec4 = vec<float, 4>;
+
+    using dvec2 = vec<double, 2>;
+    using dvec3 = vec<double, 3>;
+    using dvec4 = vec<double, 4>;
+}
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    class mat;
+
+    using bmat2 = mat<bool, 2>;
+    using bmat3 = mat<bool, 3>;
+    using bmat4 = mat<bool, 4>;
+
+    using imat2 = mat<int, 2>;
+    using imat3 = mat<int, 3>;
+    using imat4 = mat<int, 4>;
+
+    using umat2 = mat<unsigned, 2>;
+    using umat3 = mat<unsigned, 3>;
+    using umat4 = mat<unsigned, 4>;
+
+    using fmat2 = mat<float, 2>;
+    using fmat3 = mat<float, 3>;
+    using fmat4 = mat<float, 4>;
+
+    using dmat2 = mat<double, 2>;
+    using dmat3 = mat<double, 3>;
+    using dmat4 = mat<double, 4>;
+}
+
+namespace vmath_hpp
+{
+    template < typename T >
+    class qua;
+
+    using fqua = qua<float>;
+    using dqua = qua<double>;
+}
+
+//
+// Common Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr abs(T x) noexcept {
+        if constexpr ( std::is_signed_v<T> ) {
+            return x < T{0} ? -x : x;
+        } else {
+            return x;
+        }
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr sqr(T x) noexcept {
+        return 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 rcp(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_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]] 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 lerp(T x, T y, T x_a, T y_a) noexcept {
+        return x * x_a + y * 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) * rcp(edge1 - edge0), T{0}, T{1});
+        return t * t * (T{3} - T{2} * t);
+    }
+}
+
+//
+// Angle and Trigonometric Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr radians(T degrees) noexcept {
+        return degrees * T(0.01745329251994329576923690768489);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr degrees(T radians) noexcept {
+        return radians * T(57.295779513082320876798154814105);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    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);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    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);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    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);
+    }
+
+    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);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    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);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    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);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    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);
+    }
+
+    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 >
+    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);
+    }
+}
+
+//
+// Exponential Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    pow(T x, T y) noexcept {
+        return std::pow(x, y);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    exp(T x) noexcept {
+        return std::exp(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    log(T x) noexcept {
+        return std::log(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    exp2(T x) noexcept {
+        return std::exp2(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    log2(T x) noexcept {
+        return std::log2(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    sqrt(T x) noexcept {
+        return std::sqrt(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    rsqrt(T x) noexcept {
+        return rcp(sqrt(x));
+    }
+}
+
+//
+// Geometric Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, typename U
+             , typename V = decltype(std::declval<T>() * std::declval<U>()) >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<V>, V>
+    constexpr dot(T x, U y) noexcept {
+        return { x * y };
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr length(T x) noexcept {
+        return abs(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr rlength(T x) noexcept {
+        return rcp(abs(x));
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr length2(T x) noexcept {
+        return dot(x, x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr rlength2(T x) noexcept {
+        return rcp(dot(x, x));
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr distance(T x, T y) noexcept {
+        return length(y - x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, T>
+    constexpr distance2(T x, T y) noexcept {
+        return length2(y - x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr normalize(T x) noexcept {
+        return x * rlength(x);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr faceforward(T n, T i, T nref) noexcept {
+        return dot(nref, i) < T{0} ? n : -n;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    constexpr reflect(T i, T n) noexcept {
+        return i - T{2} * dot(n, i) * n;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    refract(T i, T n, T eta) noexcept {
+        const T d = dot(n, i);
+        const T k = T{1} - sqr(eta) * (T{1} - sqr(d));
+        return k < T{0} ? T{0} : (eta * i - (eta * d + sqrt(k)) * n);
+    }
+}
+
+//
+// Relational Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    constexpr any(T x) noexcept {
+        return !!x;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    constexpr all(T x) noexcept {
+        return !!x;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    constexpr approx(T x, T y) noexcept {
+        if constexpr ( std::is_floating_point_v<T> ) {
+            /// REFERENCE:
+            /// http://www.realtimecollisiondetection.net/pubs/Tolerances
+            const T epsilon = std::numeric_limits<T>::epsilon();
+            return abs(x - y) <= epsilon * max(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 approx(T x, T y, T epsilon) noexcept {
+        if constexpr ( std::is_floating_point_v<T> ) {
+            /// REFERENCE:
+            /// http://www.realtimecollisiondetection.net/pubs/Tolerances
+            return abs(x - y) <= epsilon * max(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 less(T x, T y) noexcept {
+        return x < y;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    constexpr less_equal(T x, T y) noexcept {
+        return x <= y;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    constexpr greater(T x, T y) noexcept {
+        return x > y;
+    }
+
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
+    constexpr greater_equal(T x, T y) noexcept {
+        return x >= y;
+    }
+
+    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;
+    }
+
+    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;
+    }
+}
+
+namespace vmath_hpp::detail
+{
+    template < typename T, std::size_t Size >
+    class vec_base;
+
+    template < typename T >
+    class vec_base<T, 2> {
+    public:
+        T x, y;
+    public:
+        constexpr vec_base()
+        : vec_base{zero_init} {}
+
+        constexpr vec_base(uninit_t) {}
+        constexpr vec_base(zero_init_t): vec_base{T{0}} {}
+        constexpr vec_base(unit_init_t): vec_base{T{1}} {}
+
+        constexpr explicit vec_base(T v): x{v}, y{v} {}
+        constexpr vec_base(T x, T y): x{x}, y{y} {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr vec_base(const vec_base<U, 2>& 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 vec_base<U, 3>& 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 vec_base<U, 4>& other): vec_base(other[0], other[1]) {}
+
+        [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
+            switch ( index ) {
+            default:
+            case 0: return x;
+            case 1: return y;
+            }
+        }
+
+        [[nodiscard]] constexpr const T& operator[](std::size_t index) const noexcept {
+            switch ( index ) {
+            default:
+            case 0: return x;
+            case 1: return y;
+            }
+        }
+    };
+
+    template < typename T >
+    class vec_base<T, 3> {
+    public:
+        T x, y, z;
+    public:
+        constexpr vec_base()
+        : vec_base{zero_init} {}
+
+        constexpr vec_base(uninit_t) {}
+        constexpr vec_base(zero_init_t): vec_base{T{0}} {}
+        constexpr vec_base(unit_init_t): vec_base{T{1}} {}
+
+        constexpr explicit vec_base(T v): x{v}, y{v}, z{v} {}
+        constexpr vec_base(T x, T y, T z): x{x}, y{y}, z{z} {}
+
+        constexpr vec_base(const vec_base<T, 2>& xy, T z): vec_base(xy[0], xy[1], z) {}
+        constexpr vec_base(T x, const vec_base<T, 2>& yz): vec_base(x, yz[0], yz[1]) {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr vec_base(const vec_base<U, 3>& 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 vec_base<U, 4>& other): vec_base(other[0], other[1], other[2]) {}
+
+        [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
+            switch ( index ) {
+            default:
+            case 0: return x;
+            case 1: return y;
+            case 2: return z;
+            }
+        }
+
+        [[nodiscard]] constexpr const T& operator[](std::size_t index) const noexcept {
+            switch ( index ) {
+            default:
+            case 0: return x;
+            case 1: return y;
+            case 2: return z;
+            }
+        }
+    };
+
+    template < typename T >
+    class vec_base<T, 4> {
+    public:
+        T x, y, z, w;
+    public:
+        constexpr vec_base()
+        : vec_base{zero_init} {}
+
+        constexpr vec_base(uninit_t) {}
+        constexpr vec_base(zero_init_t) : vec_base{T{0}} {}
+        constexpr vec_base(unit_init_t) : vec_base{T{1}} {}
+
+        constexpr explicit vec_base(T v): x{v}, y{v}, z{v}, w{v} {}
+        constexpr vec_base(T x, T y, T z, T w): x{x}, y{y}, z{z}, w{w} {}
+
+        constexpr vec_base(const vec_base<T, 2>& xy, T z, T w): vec_base(xy[0], xy[1], z, w) {}
+        constexpr vec_base(T x, const vec_base<T, 2>& yz, T w): vec_base(x, yz[0], yz[1], w) {}
+        constexpr vec_base(T x, T y, const vec_base<T, 2>& zw): vec_base(x, y, zw[0], zw[1]) {}
+        constexpr vec_base(const vec_base<T, 2>& xy, const vec_base<T, 2>& zw): vec_base(xy[0], xy[1], zw[0], zw[1]) {}
+
+        constexpr vec_base(const vec_base<T, 3>& xyz, T w): vec_base(xyz[0], xyz[1], xyz[2], w) {}
+        constexpr vec_base(T x, const vec_base<T, 3>& yzw): vec_base(x, yzw[0], yzw[1], yzw[2]) {}
+
+        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]) {}
+
+        [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
+            switch ( index ) {
+            default:
+            case 0: return x;
+            case 1: return y;
+            case 2: return z;
+            case 3: return w;
+            }
+        }
+
+        [[nodiscard]] constexpr const T& operator[](std::size_t index) const noexcept {
+            switch ( index ) {
+            default:
+            case 0: return x;
+            case 1: return y;
+            case 2: return z;
+            case 3: return w;
+            }
+        }
+    };
+}
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    class vec final : public detail::vec_base<T, Size> {
+    public:
+        using self_type = vec;
+        using base_type = detail::vec_base<T, Size>;
+        using component_type = T;
+
+        using pointer = component_type*;
+        using const_pointer = const component_type*;
+
+        using reference = component_type&;
+        using const_reference = const component_type&;
+
+        using iterator = pointer;
+        using const_iterator = const_pointer;
+        using reverse_iterator = std::reverse_iterator<iterator>;
+        using const_reverse_iterator = std::reverse_iterator<const_iterator>;
+
+        static inline constexpr std::size_t size = Size;
+    public:
+        using base_type::vec_base;
+        using base_type::operator[];
+
+        void swap(vec& other) noexcept(std::is_nothrow_swappable_v<T>) {
+            for ( std::size_t i = 0; i < Size; ++i ) {
+                using std::swap;
+                swap((*this)[i], other[i]);
+            }
+        }
+
+        [[nodiscard]] iterator begin() noexcept { return iterator(data()); }
+        [[nodiscard]] const_iterator begin() const noexcept { return const_iterator(data()); }
+        [[nodiscard]] iterator end() noexcept { return iterator(data() + Size); }
+        [[nodiscard]] const_iterator end() const noexcept { return const_iterator(data() + Size); }
+
+        [[nodiscard]] reverse_iterator rbegin() noexcept { return reverse_iterator(end()); }
+        [[nodiscard]] const_reverse_iterator rbegin() const noexcept { return const_reverse_iterator(end()); }
+        [[nodiscard]] reverse_iterator rend() noexcept { return reverse_iterator(begin()); }
+        [[nodiscard]] const_reverse_iterator rend() const noexcept { return const_reverse_iterator(begin()); }
+
+        [[nodiscard]] const_iterator cbegin() const noexcept { return begin(); }
+        [[nodiscard]] const_iterator cend() const noexcept { return end(); }
+        [[nodiscard]] const_reverse_iterator crbegin() const noexcept { return rbegin(); }
+        [[nodiscard]] const_reverse_iterator crend() const noexcept { return rend(); }
+
+        [[nodiscard]] pointer data() noexcept {
+            return &(*this)[0];
+        }
+
+        [[nodiscard]] const_pointer data() const noexcept {
+            return &(*this)[0];
+        }
+
+        [[nodiscard]] constexpr reference at(std::size_t index) {
+            VMATH_HPP_THROW_IF(index >= size, std::out_of_range("vec::at"));
+            return (*this)[index];
+        }
+
+        [[nodiscard]] constexpr const_reference at(std::size_t index) const {
+            VMATH_HPP_THROW_IF(index >= size, std::out_of_range("vec::at"));
+            return (*this)[index];
+        }
+    };
+}
+
+namespace vmath_hpp
+{
+    // vec2
+
+    template < typename T >
+    vec(T, T) -> vec<T, 2>;
+
+    // vec3
+
+    template < typename T >
+    vec(T, T, T) -> vec<T, 3>;
+
+    template < typename T >
+    vec(const vec<T, 2>&, T) -> vec<T, 3>;
+
+    template < typename T >
+    vec(T, const vec<T, 2>&) -> vec<T, 3>;
+
+    // vec4
+
+    template < typename T >
+    vec(T, T, T, T) -> vec<T, 4>;
+
+    template < typename T >
+    vec(const vec<T, 2>&, T, T) -> vec<T, 4>;
+
+    template < typename T >
+    vec(T, const vec<T, 2>&, T) -> vec<T, 4>;
+
+    template < typename T >
+    vec(T, T, const vec<T, 2>&) -> vec<T, 4>;
+
+    template < typename T >
+    vec(const vec<T, 2>&, const vec<T, 2>&) -> vec<T, 4>;
+
+    template < typename T >
+    vec(const vec<T, 3>&, T) -> vec<T, 4>;
+
+    template < typename T >
+    vec(T, const vec<T, 3>&) -> vec<T, 4>;
+
+    // swap
+
+    template < typename T, std::size_t Size >
+    void swap(vec<T, Size>& l, vec<T, Size>& r) noexcept(noexcept(l.swap(r))) {
+        l.swap(r);
+    }
+}
+
+namespace vmath_hpp::detail::impl
+{
+    template < typename A, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        std::index_sequence<Is...>)
+    {
+        return vec{ f(a[Is])... };
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        const vec<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return vec{ f(a[Is], b[Is])... };
+    }
+
+    template < typename A, typename B, typename C, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        const vec<B, Size>& b,
+        const vec<C, Size>& c,
+        std::index_sequence<Is...>)
+    {
+        return vec{ f(a[Is], b[Is], c[Is])... };
+    }
+
+    template < typename A, typename B, typename C, typename D, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        const vec<B, Size>& b,
+        const vec<C, Size>& c,
+        const vec<D, Size>& d,
+        std::index_sequence<Is...>)
+    {
+        return vec{ f(a[Is], b[Is], c[Is], d[Is])... };
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold_join_impl(
+        F&& f,
+        A init,
+        const vec<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return ((init = f(std::move(init), b[Is])), ...);
+    }
+
+    template < typename A, std::size_t Size, typename F, std::size_t I, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        std::index_sequence<I, Is...>)
+    {
+        A init = a[I];
+        return ((init = f(std::move(init), a[Is])), ...);
+    }
+
+    template < typename A, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        std::index_sequence<Is...>)
+    {
+        return (... && f(a[Is]));
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        const vec<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return (... && f(a[Is], b[Is]));
+    }
+
+    template < typename A, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        std::index_sequence<Is...>)
+    {
+        return (... || f(a[Is]));
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        const vec<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return (... || f(a[Is], b[Is]));
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_plus_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        const vec<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return (... + f(a[Is], b[Is]));
+    }
+}
+
+namespace vmath_hpp::detail
+{
+    template < typename A, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join(F&& f, const vec<A, Size>& a) {
+        return impl::map_join_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join(F&& f, const vec<A, Size>& a, const vec<B, Size>& b) {
+        return impl::map_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, typename C, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join(F&& f, const vec<A, Size>& a, const vec<B, Size>& b, const vec<C, Size>& c) {
+        return impl::map_join_impl(std::forward<F>(f), a, b, c, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, typename C, typename D, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join(F&& f, const vec<A, Size>& a, const vec<B, Size>& b, const vec<C, Size>& c, const vec<D, Size>& d) {
+        return impl::map_join_impl(std::forward<F>(f), a, b, c, d, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold_join(F&& f, A init, const vec<B, Size>& b) {
+        return impl::fold_join_impl(std::forward<F>(f), std::move(init), b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_join(F&& f, const vec<A, Size>& a) {
+        return impl::fold1_join_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join(F&& f, const vec<A, Size>& a) {
+        return impl::fold1_and_join_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join(F&& f, const vec<A, Size>& a, const vec<B, Size>& b) {
+        return impl::fold1_and_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join(F&& f, const vec<A, Size>& a) {
+        return impl::fold1_or_join_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join(F&& f, const vec<A, Size>& a, const vec<B, Size>& b) {
+        return impl::fold1_or_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_plus_join(F&& f, const vec<A, Size>& a, const vec<B, Size>& b) {
+        return impl::fold1_plus_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+}
+
+//
+// Operators
+//
+
+namespace vmath_hpp
+{
+    // +operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(const vec<T, Size>& xs) {
+        return map_join([](T x){ return +x; }, xs);
+    }
+
+    // -operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(const vec<T, Size>& xs) {
+        return map_join([](T x){ return -x; }, xs);
+    }
+
+    // ~operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator~(const vec<T, Size>& xs) {
+        return map_join([](T x){ return ~x; }, xs);
+    }
+
+    // !operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator!(const vec<T, Size>& xs) {
+        return map_join([](T x){ return !x; }, xs);
+    }
+
+    // ++operator
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size>& operator++(vec<T, Size>& xs) {
+        return (xs = xs + T{1});
+    }
+
+    // --operator
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size>& operator--(vec<T, Size>& xs) {
+        return (xs = xs - T{1});
+    }
+
+    // operator++
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator++(vec<T, Size>& xs, int) {
+        vec<T, Size> ys = xs;
+        ++xs;
+        return ys;
+    }
+
+    // operator--
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator--(vec<T, Size>& xs, int) {
+        vec<T, Size> ys = xs;
+        --xs;
+        return ys;
+    }
+
+    // operator+
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(const vec<T, Size>& xs, U y) {
+        return map_join([y](T x){ return x + y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(T x, const vec<U, Size>& ys) {
+        return map_join([x](U y){ return x + y; }, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(const vec<T, Size>& xs, const vec<U, Size>& ys) {
+        return map_join([](T x, U 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, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(const vec<T, Size>& xs, U y) {
+        return map_join([y](T x){ return x - y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(T x, const vec<U, Size>& ys) {
+        return map_join([x](U y){ return x - y; }, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(const vec<T, Size>& xs, const vec<U, Size>& ys) {
+        return map_join([](T x, U 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, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator*(const vec<T, Size>& xs, U y) {
+        return map_join([y](T x){ return x * y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator*(T x, const vec<U, Size>& ys) {
+        return map_join([x](U y){ return x * y; }, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator*(const vec<T, Size>& xs, const vec<U, Size>& ys) {
+        return map_join([](T x, U 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, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator/(const vec<T, Size>& xs, U y) {
+        return map_join([y](T x){ return x / y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator/(T x, const vec<U, Size>& ys) {
+        return map_join([x](U y){ return x / y; }, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator/(const vec<T, Size>& xs, const vec<U, Size>& ys) {
+        return map_join([](T x, U 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 >
+    [[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 >
+    [[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 >
+    [[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 >
+    [[nodiscard]] constexpr bool operator==(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return fold1_and_join([](T x, T y){ return x == y; }, xs, ys);
+    }
+
+    // operator!=
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr bool operator!=(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return fold1_or_join([](T x, T y){ return x != y; }, xs, ys);
+    }
+
+    // operator<
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr bool operator<(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        for ( std::size_t i = 0; i < Size; ++i ) {
+            if ( xs[i] < ys[i] ) {
+                return true;
+            }
+            if ( ys[i] < xs[i] ) {
+                return false;
+            }
+        }
+        return false;
+    }
+}
+
+//
+// Common Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> abs(const vec<T, Size>& xs) {
+        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);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> rcp(const vec<T, Size>& xs) {
+        return map_join([](T x) { return rcp(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> floor(const vec<T, Size>& xs) {
+        return map_join([](T x) { return floor(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> trunc(const vec<T, Size>& xs) {
+        return map_join([](T x) { return trunc(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> round(const vec<T, Size>& xs) {
+        return map_join([](T x) { return round(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> ceil(const vec<T, Size>& xs) {
+        return map_join([](T x) { return ceil(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> fract(const vec<T, Size>& xs) {
+        return map_join([](T x) { return fract(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> fmod(const vec<T, Size>& xs, T y) {
+        return map_join([y](T x) { return fmod(x, y); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> fmod(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y) { return fmod(x, y); }, xs, ys);
+    }
+
+    namespace impl
+    {
+        template < typename T, std::size_t Size, std::size_t... Is >
+        constexpr VMATH_HPP_FORCE_INLINE
+        vec<T, Size> modf_impl(const vec<T, Size>& xs, vec<T, Size>* is, std::index_sequence<Is...>) {
+            return { modf(xs[Is], &(*is)[Is])... };
+        }
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> modf(const vec<T, Size>& xs, vec<T, Size>* is) {
+        return impl::modf_impl(xs, is, std::make_index_sequence<Size>{});
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr 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]] constexpr 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);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> min(const vec<T, Size>& xs, T y) {
+        return map_join([y](T x) { return min(x, y); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> min(T x, const vec<T, Size>& ys) {
+        return map_join([x](T y) { return min(x, y); }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> min(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y) { return min(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T max(const vec<T, Size>& xs) {
+        return fold1_join([](T acc, T x){ return max(acc, x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> max(const vec<T, Size>& xs, T y) {
+        return map_join([y](T x) { return max(x, y); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> max(T x, const vec<T, Size>& ys) {
+        return map_join([x](T y) { return max(x, y); }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> max(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y) { return max(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> clamp(const vec<T, Size>& xs, T min_x, T max_x) {
+        return map_join([min_x, max_x](T x) { return clamp(x, min_x, max_x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> clamp(const vec<T, Size>& xs, const vec<T, Size>& min_xs, const vec<T, Size>& max_xs) {
+        return map_join([](T x, T min_x, T max_x) { return clamp(x, min_x, max_x); }, xs, min_xs, max_xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> saturate(const vec<T, Size>& xs) {
+        return map_join([](T x) { return saturate(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> lerp(const vec<T, Size>& xs, const vec<T, Size>& ys, T a) {
+        return map_join([a](T x, T y) { return lerp(x, y, a); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> lerp(const vec<T, Size>& xs, const vec<T, Size>& ys, T x_a, T y_a) {
+        return map_join([x_a, y_a](T x, T y) { return lerp(x, y, x_a, y_a); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> lerp(
+        const vec<T, Size>& xs,
+        const vec<T, Size>& ys,
+        const vec<T, Size>& as)
+    {
+        return map_join([](T x, T y, T a) { return lerp(x, y, a); }, xs, ys, as);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> lerp(
+        const vec<T, Size>& xs,
+        const vec<T, Size>& ys,
+        const vec<T, Size>& xs_a,
+        const vec<T, Size>& ys_a)
+    {
+        return map_join([](T x, T y, T x_a, T y_a) { return lerp(x, y, x_a, y_a); }, xs, ys, xs_a, ys_a);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> step(T edge, const vec<T, Size>& xs) {
+        return map_join([edge](T x) { return step(edge, x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> step(const vec<T, Size>& edges, const vec<T, Size>& xs) {
+        return map_join([](T edge, T x) { return step(edge, x); }, edges, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> smoothstep(T edge0, T edge1, const vec<T, Size>& xs) {
+        return map_join([edge0, edge1](T x) { return smoothstep(edge0, edge1, x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> smoothstep(const vec<T, Size>& edges0, const vec<T, Size>& edges1, const vec<T, Size>& xs) {
+        return map_join([](T edge0, T edge1, T x) { return smoothstep(edge0, edge1, x); }, edges0, edges1, xs);
+    }
+}
+
+//
+// Angle and Trigonometric Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> radians(const vec<T, Size>& degrees) {
+        return map_join([](T x) { return radians(x); }, degrees);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> degrees(const vec<T, Size>& radians) {
+        return map_join([](T x) { return degrees(x); }, radians);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> sin(const vec<T, Size>& xs) {
+        return map_join([](T x) { return sin(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> cos(const vec<T, Size>& xs) {
+        return map_join([](T x) { return cos(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> tan(const vec<T, Size>& xs) {
+        return map_join([](T x) { return tan(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> asin(const vec<T, Size>& xs) {
+        return map_join([](T x) { return asin(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> acos(const vec<T, Size>& xs) {
+        return map_join([](T x) { return acos(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> atan(const vec<T, Size>& xs) {
+        return map_join([](T x) { return atan(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> atan2(const vec<T, Size>& ys, const vec<T, Size>& xs) {
+        return map_join([](T y, T x) { return atan2(y, x); }, ys, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> sinh(const vec<T, Size>& xs) {
+        return map_join([](T x) { return sinh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> cosh(const vec<T, Size>& xs) {
+        return map_join([](T x) { return cosh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> tanh(const vec<T, Size>& xs) {
+        return map_join([](T x) { return tanh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> asinh(const vec<T, Size>& xs) {
+        return map_join([](T x) { return asinh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> acosh(const vec<T, Size>& xs) {
+        return map_join([](T x) { return acosh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> atanh(const vec<T, Size>& xs) {
+        return map_join([](T x) { return atanh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr std::pair<vec<T, Size>, vec<T, Size>> sincos(const vec<T, Size>& xs) {
+        return { sin(xs), cos(xs) };
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr void sincos(const vec<T, Size>& xs, vec<T, Size>* ss, vec<T, Size>* cs) {
+        *ss = sin(xs);
+        *cs = cos(xs);
+    }
+}
+
+//
+// Exponential Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> pow(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y) { return pow(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> exp(const vec<T, Size>& xs) {
+        return map_join([](T x) { return exp(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> log(const vec<T, Size>& xs) {
+        return map_join([](T x) { return log(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> exp2(const vec<T, Size>& xs) {
+        return map_join([](T x) { return exp2(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> log2(const vec<T, Size>& xs) {
+        return map_join([](T x) { return log2(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> sqrt(const vec<T, Size>& xs) {
+        return map_join([](T x) { return sqrt(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> rsqrt(const vec<T, Size>& xs) {
+        return map_join([](T x) { return rsqrt(x); }, xs);
+    }
+}
+
+//
+// Geometric Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, typename U, std::size_t Size
+             , typename V = decltype(std::declval<T>() * std::declval<U>()) >
+    [[nodiscard]] constexpr V dot(const vec<T, Size>& xs, const vec<U, Size>& ys) {
+        return fold1_plus_join([](T x, U y){ return x * y; }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T length(const vec<T, Size>& xs) {
+        return sqrt(dot(xs, xs));
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T rlength(const vec<T, Size>& xs) {
+        return rsqrt(dot(xs, xs));
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T length2(const vec<T, Size>& xs) {
+        return dot(xs, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T rlength2(const vec<T, Size>& xs) {
+        return rcp(dot(xs, xs));
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T distance(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return length(ys - xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T distance2(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return length2(ys - xs);
+    }
+
+    template < typename T, typename U
+             , typename V = decltype(std::declval<T>() * std::declval<U>()) >
+    [[nodiscard]] constexpr V cross(const vec<T, 2>& xs, const vec<U, 2>& ys) {
+        return { xs.x * ys.y - xs.y * ys.x };
+    }
+
+    template < typename T, typename U
+             , typename V = decltype(std::declval<T>() * std::declval<U>()) >
+    [[nodiscard]] constexpr vec<V, 3> cross(const vec<T, 3>& xs, const vec<U, 3>& ys) {
+        return {
+            xs.y * ys.z - xs.z * ys.y,
+            xs.z * ys.x - xs.x * ys.z,
+            xs.x * ys.y - xs.y * ys.x};
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> normalize(const vec<T, Size>& xs) {
+        return xs * rlength(xs);
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> faceforward(const vec<T, Size>& n, const vec<T, Size>& i, const vec<T, Size>& nref) {
+        return dot(nref, i) < T{0} ? n : -n;
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> reflect(const vec<T, Size>& i, const vec<T, Size>& n) {
+        return i - T{2} * dot(n, i) * n;
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> refract(const vec<T, Size>& i, const vec<T, Size>& n, T eta) {
+        const T d = dot(n, i);
+        const T k = T{1} - sqr(eta) * (T{1} - sqr(d));
+        return k < T{0} ? vec<T, Size>{T{0}} : (eta * i - (eta * d + sqrt(k)) * n);
+    }
+}
+
+//
+// Relational Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size
+             , typename U = decltype(any(std::declval<T>())) >
+    [[nodiscard]] constexpr U any(const vec<T, Size>& xs) {
+        return fold1_or_join([](T x){ return any(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(all(std::declval<T>())) >
+    [[nodiscard]] constexpr U all(const vec<T, Size>& xs) {
+        return fold1_and_join([](T x){ return all(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(approx(
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> approx(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return approx(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(approx(
+                 std::declval<T>(),
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> approx(const vec<T, Size>& xs, const vec<T, Size>& ys, T epsilon) {
+        return map_join([epsilon](T x, T y){ return approx(x, y, epsilon); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(less(
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> less(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return less(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(less_equal(
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> less_equal(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return less_equal(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(greater(
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> greater(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return greater(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(greater_equal(
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> greater_equal(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return greater_equal(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(equal_to(
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> equal_to(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return equal_to(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(not_equal_to(
+                 std::declval<T>(),
+                 std::declval<T>())) >
+    [[nodiscard]] constexpr vec<U, Size> not_equal_to(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return map_join([](T x, T y){ return not_equal_to(x, y); }, xs, ys);
+    }
+}
+
+namespace vmath_hpp::detail
+{
+    template < typename T, std::size_t Size >
+    class mat_base;
+
+    template < typename T >
+    class mat_base<T, 2> {
+    public:
+        using row_type = vec<T, 2>;
+        row_type rows[2];
+    public:
+        constexpr mat_base()
+        : mat_base(identity_init) {}
+
+        constexpr mat_base(uninit_t) {}
+        constexpr mat_base(zero_init_t): mat_base{zero_init, zero_init} {}
+        constexpr mat_base(unit_init_t): mat_base{unit_init, unit_init} {}
+        constexpr mat_base(identity_init_t): mat_base{T{1}} {}
+
+        constexpr explicit mat_base(T d)
+        : rows{
+            {d, T{0}},
+            {T{0}, d}} {}
+
+        constexpr explicit mat_base(const row_type& d)
+        : rows{
+            {d[0], T{0}},
+            {T{0}, d[1]}} {}
+
+        constexpr mat_base(
+            T m11, T m12,
+            T m21, T m22)
+        : rows{
+            {m11, m12},
+            {m21, m22}} {}
+
+        constexpr mat_base(
+            const row_type& row0,
+            const row_type& row1)
+        : rows{row0, row1} {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr mat_base(const mat_base<U, 2>& 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 mat_base<U, 3>& 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 mat_base<U, 4>& other): mat_base(
+            row_type{other.rows[0]},
+            row_type{other.rows[1]}) {}
+    };
+
+    template < typename T >
+    class mat_base<T, 3> {
+    public:
+        using row_type = vec<T, 3>;
+        row_type rows[3];
+    public:
+        constexpr mat_base()
+        : mat_base(identity_init) {}
+
+        constexpr mat_base(uninit_t) {}
+        constexpr mat_base(zero_init_t): mat_base{zero_init, zero_init, zero_init} {}
+        constexpr mat_base(unit_init_t): mat_base{unit_init, unit_init, unit_init} {}
+        constexpr mat_base(identity_init_t): mat_base{T{1}} {}
+
+        constexpr explicit mat_base(T d)
+        : rows{
+            {d, T{0}, T{0}},
+            {T{0}, d, T{0}},
+            {T{0}, T{0}, d}} {}
+
+        constexpr explicit mat_base(const row_type& d)
+        : rows{
+            {d[0], T{0}, T{0}},
+            {T{0}, d[1], T{0}},
+            {T{0}, T{0}, d[2]}} {}
+
+        constexpr mat_base(
+            T m11, T m12, T m13,
+            T m21, T m22, T m23,
+            T m31, T m32, T m33)
+        : rows{
+            {m11, m12, m13},
+            {m21, m22, m23},
+            {m31, m32, m33}} {}
+
+        constexpr mat_base(
+            const row_type& row0,
+            const row_type& row1,
+            const row_type& row2)
+        : rows{row0, row1, row2} {}
+
+        constexpr mat_base(
+            const mat_base<T, 2>& m,
+            const vec_base<T, 2>& v)
+        : rows{
+            {m.rows[0], T{0}},
+            {m.rows[1], T{0}},
+            {v, T{1}}} {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr mat_base(const mat_base<U, 3>& other): mat_base(
+            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 mat_base<U, 2>& other): mat_base(
+            row_type{other.rows[0], T{0}},
+            row_type{other.rows[1], T{0}},
+            row_type{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 mat_base<U, 4>& other): mat_base(
+            row_type{other.rows[0]},
+            row_type{other.rows[1]},
+            row_type{other.rows[2]}) {}
+    };
+
+    template < typename T >
+    class mat_base<T, 4> {
+    public:
+        using row_type = vec<T, 4>;
+        row_type rows[4];
+    public:
+        constexpr mat_base()
+        : mat_base(identity_init) {}
+
+        constexpr mat_base(uninit_t) {}
+        constexpr mat_base(zero_init_t): mat_base{zero_init, zero_init, zero_init, zero_init} {}
+        constexpr mat_base(unit_init_t): mat_base{unit_init, unit_init, unit_init, unit_init} {}
+        constexpr mat_base(identity_init_t): mat_base{T{1}} {}
+
+        constexpr explicit mat_base(T d)
+        : rows{
+            {d, T{0}, T{0}, T{0}},
+            {T{0}, d, T{0}, T{0}},
+            {T{0}, T{0}, d, T{0}},
+            {T{0}, T{0}, T{0}, d}} {}
+
+        constexpr explicit mat_base(const row_type& d)
+        : rows{
+            {d[0], T{0}, T{0}, T{0}},
+            {T{0}, d[1], T{0}, T{0}},
+            {T{0}, T{0}, d[2], T{0}},
+            {T{0}, T{0}, T{0}, d[3]}} {}
+
+        constexpr mat_base(
+            T m11, T m12, T m13, T m14,
+            T m21, T m22, T m23, T m24,
+            T m31, T m32, T m33, T m34,
+            T m41, T m42, T m43, T m44)
+        : rows{
+            {m11, m12, m13, m14},
+            {m21, m22, m23, m24},
+            {m31, m32, m33, m34},
+            {m41, m42, m43, m44}} {}
+
+        constexpr mat_base(
+            const row_type& row0,
+            const row_type& row1,
+            const row_type& row2,
+            const row_type& row3)
+        : rows{row0, row1, row2, row3} {}
+
+        constexpr mat_base(
+            const mat_base<T, 3>& m,
+            const vec_base<T, 3>& v)
+        : rows{
+            {m.rows[0], T{0}},
+            {m.rows[1], T{0}},
+            {m.rows[2], T{0}},
+            {v, T{1}}} {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr mat_base(const mat_base<U, 4>& other): mat_base(
+            row_type{other.rows[0]},
+            row_type{other.rows[1]},
+            row_type{other.rows[2]},
+            row_type{other.rows[3]}) {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr explicit mat_base(const mat_base<U, 2>& other): mat_base(
+            row_type{other.rows[0], T{0}, T{0}},
+            row_type{other.rows[1], T{0}, T{0}},
+            row_type{T{0}, T{0}, T{1}, 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 mat_base<U, 3>& other): mat_base(
+            row_type{other.rows[0], T{0}},
+            row_type{other.rows[1], T{0}},
+            row_type{other.rows[2], T{0}},
+            row_type{T{0}, T{0}, T{0}, T{1}}) {}
+    };
+}
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    class mat final : public detail::mat_base<T, Size> {
+    public:
+        using self_type = mat;
+        using base_type = detail::mat_base<T, Size>;
+        using component_type = T;
+
+        using row_type = vec<T, Size>;
+
+        using pointer = row_type*;
+        using const_pointer = const row_type*;
+
+        using reference = row_type&;
+        using const_reference = const row_type&;
+
+        using iterator = pointer;
+        using const_iterator = const_pointer;
+        using reverse_iterator = std::reverse_iterator<iterator>;
+        using const_reverse_iterator = std::reverse_iterator<const_iterator>;
+
+        static inline constexpr std::size_t size = Size;
+    public:
+        using base_type::mat_base;
+        using base_type::rows;
+
+        void swap(mat& other) noexcept(std::is_nothrow_swappable_v<T>) {
+            for ( std::size_t i = 0; i < Size; ++i ) {
+                using std::swap;
+                swap(rows[i], other.rows[i]);
+            }
+        }
+
+        [[nodiscard]] iterator begin() noexcept { return iterator(data()); }
+        [[nodiscard]] const_iterator begin() const noexcept { return const_iterator(data()); }
+        [[nodiscard]] iterator end() noexcept { return iterator(data() + Size); }
+        [[nodiscard]] const_iterator end() const noexcept { return const_iterator(data() + Size); }
+
+        [[nodiscard]] reverse_iterator rbegin() noexcept { return reverse_iterator(end()); }
+        [[nodiscard]] const_reverse_iterator rbegin() const noexcept { return const_reverse_iterator(end()); }
+        [[nodiscard]] reverse_iterator rend() noexcept { return reverse_iterator(begin()); }
+        [[nodiscard]] const_reverse_iterator rend() const noexcept { return const_reverse_iterator(begin()); }
+
+        [[nodiscard]] const_iterator cbegin() const noexcept { return begin(); }
+        [[nodiscard]] const_iterator cend() const noexcept { return end(); }
+        [[nodiscard]] const_reverse_iterator crbegin() const noexcept { return rbegin(); }
+        [[nodiscard]] const_reverse_iterator crend() const noexcept { return rend(); }
+
+        [[nodiscard]] pointer data() noexcept {
+            return &rows[0];
+        }
+
+        [[nodiscard]] const_pointer data() const noexcept {
+            return &rows[0];
+        }
+
+        [[nodiscard]] constexpr reference at(std::size_t index) {
+            VMATH_HPP_THROW_IF(index >= size, std::out_of_range("mat::at"));
+            return rows[index];
+        }
+
+        [[nodiscard]] constexpr const_reference at(std::size_t index) const {
+            VMATH_HPP_THROW_IF(index >= size, std::out_of_range("mat::at"));
+            return rows[index];
+        }
+
+        [[nodiscard]] constexpr reference operator[](std::size_t index) noexcept {
+            return rows[index];
+        }
+
+        [[nodiscard]] constexpr const_reference operator[](std::size_t index) const noexcept {
+            return rows[index];
+        }
+    };
+}
+
+namespace vmath_hpp
+{
+    // mat2
+
+    template < typename T >
+    mat(T, T, T, T) -> mat<T, 2>;
+
+    template < typename T >
+    mat(const vec<T, 2>&, const vec<T, 2>&) -> mat<T, 2>;
+
+    template < typename T >
+    mat(std::initializer_list<T>, std::initializer_list<T>) -> mat<T, 2>;
+
+    // mat3
+
+    template < typename T >
+    mat(T, T, T, T, T, T, T, T, T) -> mat<T, 3>;
+
+    template < typename T >
+    mat(const vec<T, 3>&, const vec<T, 3>&, const vec<T, 3>&) -> mat<T, 3>;
+
+    template < typename T >
+    mat(const mat<T, 2>&, const vec<T, 2>&) -> mat<T, 3>;
+
+    template < typename T >
+    mat(std::initializer_list<T>, std::initializer_list<T>, std::initializer_list<T>) -> mat<T, 3>;
+
+    // mat4
+
+    template < typename T >
+    mat(T, T, T, T, T, T, T, T, T, T, T, T, T, T, T, T) -> mat<T, 4>;
+
+    template < typename T >
+    mat(const vec<T, 4>&, const vec<T, 4>&, const vec<T, 4>&, const vec<T, 4>&) -> mat<T, 4>;
+
+    template < typename T >
+    mat(const mat<T, 3>&, const vec<T, 3>&) -> mat<T, 4>;
+
+    template < typename T >
+    mat(std::initializer_list<T>, std::initializer_list<T>, std::initializer_list<T>, std::initializer_list<T>) -> mat<T, 4>;
+
+    // swap
+
+    template < typename T, std::size_t Size >
+    void swap(mat<T, Size>& l, mat<T, Size>& r) noexcept(noexcept(l.swap(r))) {
+        l.swap(r);
+    }
+}
+
+namespace vmath_hpp::detail::impl
+{
+    template < typename A, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join_impl(
+        F&& f,
+        const mat<A, Size>& a,
+        std::index_sequence<Is...>)
+    {
+        return mat{ f(a[Is])... };
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join_impl(
+        F&& f,
+        const mat<A, Size>& a,
+        const mat<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return mat{ f(a[Is], b[Is])... };
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold_join_impl(
+        F&& f,
+        A init,
+        const mat<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return ((init = f(std::move(init), b[Is])), ...);
+    }
+
+    template < typename A, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join_impl(
+        F&& f,
+        const mat<A, Size>& a,
+        std::index_sequence<Is...>)
+    {
+        return (... && f(a[Is]));
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join_impl(
+        F&& f,
+        const mat<A, Size>& a,
+        const mat<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return (... && f(a[Is], b[Is]));
+    }
+
+    template < typename A, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join_impl(
+        F&& f,
+        const mat<A, Size>& a,
+        std::index_sequence<Is...>)
+    {
+        return (... || f(a[Is]));
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join_impl(
+        F&& f,
+        const mat<A, Size>& a,
+        const mat<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return (... || f(a[Is], b[Is]));
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_plus_join_impl(
+        F&& f,
+        const vec<A, Size>& a,
+        const mat<B, Size>& b,
+        std::index_sequence<Is...>)
+    {
+        return (... + f(a[Is], b[Is]));
+    }
+}
+
+namespace vmath_hpp::detail
+{
+    template < typename A, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join(F&& f, const mat<A, Size>& a) {
+        return impl::map_join_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join(F&& f, const mat<A, Size>& a, const mat<B, Size>& b) {
+        return impl::map_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold_join(F&& f, A init, const mat<B, Size>& b) {
+        return impl::fold_join_impl(std::forward<F>(f), std::move(init), b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join(F&& f, const mat<A, Size>& a) {
+        return impl::fold1_and_join_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_and_join(F&& f, const mat<A, Size>& a, const mat<B, Size>& b) {
+        return impl::fold1_and_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join(F&& f, const mat<A, Size>& a) {
+        return impl::fold1_or_join_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_or_join(F&& f, const mat<A, Size>& a, const mat<B, Size>& b) {
+        return impl::fold1_or_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold1_plus_join(F&& f, const vec<A, Size>& a, const mat<B, Size>& b) {
+        return impl::fold1_plus_join_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+}
+
+//
+// Operators
+//
+
+namespace vmath_hpp
+{
+    // +operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(const mat<T, Size>& xs) {
+        return map_join([](const vec<T, Size>& x){ return +x; }, xs);
+    }
+
+    // -operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(const mat<T, Size>& xs) {
+        return map_join([](const vec<T, Size>& x){ return -x; }, xs);
+    }
+
+    // ~operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator~(const mat<T, Size>& xs) {
+        return map_join([](const vec<T, Size>& x){ return ~x; }, xs);
+    }
+
+    // !operator
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr auto operator!(const mat<T, Size>& xs) {
+        return map_join([](const vec<T, Size>& x){ return !x; }, xs);
+    }
+
+    // ++operator
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size>& operator++(mat<T, Size>& xs) {
+        return (xs = xs + T{1});
+    }
+
+    // --operator
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size>& operator--(mat<T, Size>& xs) {
+        return (xs = xs - T{1});
+    }
+
+    // operator++
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator++(mat<T, Size>& xs, int) {
+        mat<T, Size> ys = xs;
+        ++xs;
+        return ys;
+    }
+
+    // operator--
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator--(mat<T, Size>& xs, int) {
+        mat<T, Size> ys = xs;
+        --xs;
+        return ys;
+    }
+
+    // operator+
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(const mat<T, Size>& xs, U y) {
+        return map_join([y](const vec<T, Size>& x){ return x + y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(T x, const mat<U, Size>& ys) {
+        return map_join([x](const vec<U, Size>& y){ return x + y; }, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator+(const mat<T, Size>& xs, const mat<U, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<U, 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, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(const mat<T, Size>& xs, U y) {
+        return map_join([y](const vec<T, Size>& x){ return x - y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(T x, const mat<U, Size>& ys) {
+        return map_join([x](const vec<U, Size>& y){ return x - y; }, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator-(const mat<T, Size>& xs, const mat<U, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<U, 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, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator*(const mat<T, Size>& xs, U y) {
+        return map_join([y](const vec<T, Size>& x){ return x * y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator*(T x, const mat<U, Size>& ys) {
+        return map_join([x](const vec<U, Size>& y){ return x * y; }, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator*(const vec<T, Size>& xs, const mat<U, Size>& ys) {
+        return fold1_plus_join([](T x, const vec<U, Size>& y){ return x * y; }, xs, ys);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator*(const mat<T, Size>& xs, const mat<U, Size>& ys) {
+        return map_join([&ys](const vec<T, Size>& x){ return x * ys; }, xs);
+    }
+
+    // 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 vec<T, Size>& operator*=(vec<T, Size>& xs, const mat<T, Size>& ys) {
+        return (xs = (xs * ys));
+    }
+
+    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, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator/(const mat<T, Size>& xs, U y) {
+        return map_join([y](const vec<T, Size>& x){ return x / y; }, xs);
+    }
+
+    template < typename T, typename U, std::size_t Size >
+    [[nodiscard]] constexpr auto operator/(T x, const mat<U, Size>& ys) {
+        return map_join([x](const vec<U, Size>& y){ return x / y; }, ys);
+    }
+
+    // operator/=
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size>& operator/=(mat<T, Size>& xs, T y) {
+        return (xs = (xs / y));
+    }
+
+    // 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 >
+    [[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 >
+    [[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 >
+    [[nodiscard]] constexpr bool operator==(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return fold1_and_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return x == y; }, xs, ys);
+    }
+
+    // operator!=
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr bool operator!=(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return fold1_or_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return x != y; }, xs, ys);
+    }
+
+    // operator<
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr bool operator<(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        for ( std::size_t i = 0; i < Size; ++i ) {
+            if ( xs[i] < ys[i] ) {
+                return true;
+            }
+            if ( ys[i] < xs[i] ) {
+                return false;
+            }
+        }
+        return false;
+    }
+}
+
+//
+// Relational Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size
+             , typename U = decltype(any(std::declval<vec<T, Size>>())) >
+    [[nodiscard]] constexpr U any(const mat<T, Size>& xs) {
+        return fold1_or_join([](const vec<T, Size>& x){ return any(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = decltype(all(std::declval<vec<T, Size>>())) >
+    [[nodiscard]] constexpr U all(const mat<T, Size>& xs) {
+        return fold1_and_join([](const vec<T, Size>& x){ return all(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(approx(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> approx(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return approx(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(approx(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>(),
+                 std::declval<T>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> approx(const mat<T, Size>& xs, const mat<T, Size>& ys, T epsilon) {
+        return map_join([epsilon](const vec<T, Size>& x, const vec<T, Size>& y){ return approx(x, y, epsilon); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(less(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> less(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return less(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(less_equal(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> less_equal(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return less_equal(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(greater(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> greater(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return greater(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(greater_equal(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> greater_equal(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return greater_equal(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(equal_to(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> equal_to(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return equal_to(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size
+             , typename U = typename decltype(not_equal_to(
+                 std::declval<vec<T, Size>>(),
+                 std::declval<vec<T, Size>>()))::component_type >
+    [[nodiscard]] constexpr mat<U, Size> not_equal_to(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return map_join([](const vec<T, Size>& x, const vec<T, Size>& y){ return not_equal_to(x, y); }, xs, ys);
+    }
+}
+
+//
+// Matrix Functions
+//
+
+namespace vmath_hpp
+{
+    namespace impl
+    {
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        mat<T, 2> transpose_2x2_impl(
+            T a, T c,
+            T b, T d)
+        {
+            return {
+                a, b,
+                c, d};
+        }
+
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        mat<T, 3> transpose_3x3_impl(
+            T a, T d, T g,
+            T b, T e, T h,
+            T c, T f, T i)
+        {
+            return {
+                a, b, c,
+                d, e, f,
+                g, h, i};
+        }
+
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        mat<T, 4> transpose_4x4_impl(
+            T a, T e, T i, T m,
+            T b, T f, T j, T n,
+            T c, T g, T k, T o,
+            T d, T h, T l, T p)
+        {
+            return {
+                a, b, c, d,
+                e, f, g, h,
+                i, j, k, l,
+                m, n, o, p};
+        }
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 2> transpose(const mat<T, 2>& m) {
+        return impl::transpose_2x2_impl(
+            m[0][0], m[0][1],
+            m[1][0], m[1][1]);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> transpose(const mat<T, 3>& m) {
+        return impl::transpose_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> transpose(const mat<T, 4>& m) {
+        return impl::transpose_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]);
+    }
+
+    namespace impl
+    {
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        T determinant_2x2_impl(
+            T a, T b,
+            T c, T d)
+        {
+            /// REFERENCE:
+            /// http://www.euclideanspace.com/maths/algebra/matrix/functions/determinant/twoD/
+
+            return
+                + a * d
+                - b * c;
+        }
+
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        T determinant_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/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);
+
+        }
+
+        template < typename T >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        T determinant_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/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);
+        }
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr T determinant(const mat<T, 2>& m) {
+        return impl::determinant_2x2_impl(
+            m[0][0], m[0][1],
+            m[1][0], m[1][1]);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr T determinant(const mat<T, 3>& m) {
+        return impl::determinant_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 T determinant(const mat<T, 4>& m) {
+        return impl::determinant_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]);
+    }
+
+    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/
+
+            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]);
+    }
+}
+
+namespace vmath_hpp::detail
+{
+    template < typename T >
+    class qua_base {
+    public:
+        vec<T, 3> v;
+        T s;
+    public:
+        constexpr qua_base()
+        : qua_base(identity_init) {}
+
+        constexpr qua_base(uninit_t) {}
+        constexpr qua_base(zero_init_t) : qua_base{zero_init, T{0}} {}
+        constexpr qua_base(identity_init_t) : qua_base{zero_init, T{1}} {}
+
+        constexpr qua_base(T vx, T vy, T vz, T s): v{vx, vy, vz}, s{s} {}
+        constexpr qua_base(const vec<T, 3>& v, T s): v{v}, s{s} {}
+        constexpr explicit qua_base(const vec<T, 4>& vs): v{vs[0], vs[1], vs[2]}, s{vs[3]} {}
+
+        template < typename U, std::enable_if_t<std::is_convertible_v<U, T>, int> = 0 >
+        constexpr qua_base(const qua_base<U>& other): qua_base(other.v, other.s) {}
+
+        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); }
+
+        [[nodiscard]] constexpr T& operator[](std::size_t index) noexcept {
+            switch ( index ) {
+            default:
+            case 0: return v.x;
+            case 1: return v.y;
+            case 2: return v.z;
+            case 3: return s;
+            }
+        }
+
+        [[nodiscard]] constexpr const T& operator[](std::size_t index) const noexcept {
+            switch ( index ) {
+            default:
+            case 0: return v.x;
+            case 1: return v.y;
+            case 2: return v.z;
+            case 3: return s;
+            }
+        }
+    };
+}
+
+namespace vmath_hpp
+{
+    template < typename T >
+    class qua final : public detail::qua_base<T> {
+    public:
+        using self_type = qua;
+        using base_type = detail::qua_base<T>;
+        using component_type = T;
+
+        using imag_type = vec<T, 3>;
+        using real_type = T;
+
+        using pointer = component_type*;
+        using const_pointer = const component_type*;
+
+        using reference = component_type&;
+        using const_reference = const component_type&;
+
+        using iterator = pointer;
+        using const_iterator = const_pointer;
+        using reverse_iterator = std::reverse_iterator<iterator>;
+        using const_reverse_iterator = std::reverse_iterator<const_iterator>;
+
+        static inline constexpr std::size_t size = 4;
+    public:
+        using base_type::qua_base;
+        using base_type::operator[];
+
+        void swap(qua& other) noexcept(std::is_nothrow_swappable_v<T>) {
+            for ( std::size_t i = 0; i < size; ++i ) {
+                using std::swap;
+                swap((*this)[i], other[i]);
+            }
+        }
+
+        [[nodiscard]] iterator begin() noexcept { return iterator(data()); }
+        [[nodiscard]] const_iterator begin() const noexcept { return const_iterator(data()); }
+        [[nodiscard]] iterator end() noexcept { return iterator(data() + size); }
+        [[nodiscard]] const_iterator end() const noexcept { return const_iterator(data() + size); }
+
+        [[nodiscard]] reverse_iterator rbegin() noexcept { return reverse_iterator(end()); }
+        [[nodiscard]] const_reverse_iterator rbegin() const noexcept { return const_reverse_iterator(end()); }
+        [[nodiscard]] reverse_iterator rend() noexcept { return reverse_iterator(begin()); }
+        [[nodiscard]] const_reverse_iterator rend() const noexcept { return const_reverse_iterator(begin()); }
+
+        [[nodiscard]] const_iterator cbegin() const noexcept { return begin(); }
+        [[nodiscard]] const_iterator cend() const noexcept { return end(); }
+        [[nodiscard]] const_reverse_iterator crbegin() const noexcept { return rbegin(); }
+        [[nodiscard]] const_reverse_iterator crend() const noexcept { return rend(); }
+
+        [[nodiscard]] pointer data() noexcept {
+            return &(*this)[0];
+        }
+
+        [[nodiscard]] const_pointer data() const noexcept {
+            return &(*this)[0];
+        }
+
+        [[nodiscard]] constexpr reference at(std::size_t index) {
+            VMATH_HPP_THROW_IF(index >= size, std::out_of_range("qua::at"));
+            return (*this)[index];
+        }
+
+        [[nodiscard]] constexpr const_reference at(std::size_t index) const {
+            VMATH_HPP_THROW_IF(index >= size, std::out_of_range("qua::at"));
+            return (*this)[index];
+        }
+    };
+}
+
+namespace vmath_hpp
+{
+    // vec
+
+    template < typename T >
+    vec(const qua<T>&) -> vec<T, 4>;
+
+    // qua
+
+    template < typename T >
+    qua(T, T, T, T) -> qua<T>;
+
+    template < typename T >
+    qua(const vec<T, 3>&, T) -> qua<T>;
+
+    template < typename T >
+    qua(const vec<T, 4>&) -> qua<T>;
+
+    template < typename T >
+    qua(std::initializer_list<T>) -> qua<T>;
+
+    // swap
+
+    template < typename T >
+    void swap(qua<T>& l, qua<T>& r) noexcept(noexcept(l.swap(r))) {
+        l.swap(r);
+    }
+}
+
+namespace vmath_hpp::detail
+{
+    template < typename A, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto map_join(F&& f, const qua<A>& a) {
+        return qua(map_join(std::forward<F>(f), vec{a}));
+    }
+
+    template < typename A, typename B, typename F >
+    [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+    auto fold_join(F&& f, A init, const qua<B>& b) {
+        return fold_join(std::forward<F>(f), std::move(init), vec{b});
+    }
+}
+
+//
+// Operators
+//
+
+namespace vmath_hpp
+{
+    // +operator
+
+    template < typename T >
+    [[nodiscard]] constexpr auto operator+(const qua<T>& xs) {
+        return qua(+vec<T, 4>{xs});
+    }
+
+    // -operator
+
+    template < typename T >
+    [[nodiscard]] constexpr auto operator-(const qua<T>& xs) {
+        return qua(-vec<T, 4>{xs});
+    }
+
+    // operator+
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator+(const qua<T>& xs, const qua<U>& ys) {
+        return qua(vec{xs} + vec{ys});
+    }
+
+    // operator+=
+
+    template < typename T >
+    constexpr qua<T>& operator+=(qua<T>& xs, const qua<T>& ys) {
+        return (xs = (xs + ys));
+    }
+
+    // operator-
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator-(const qua<T>& xs, const qua<U>& ys) {
+        return qua(vec{xs} - vec{ys});
+    }
+
+    // operator-=
+
+    template < typename T >
+    constexpr qua<T>& operator-=(qua<T>& xs, const qua<T>& ys) {
+        return (xs = (xs - ys));
+    }
+
+    // operator*
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator*(const qua<T>& xs, U y) {
+        return qua(vec{xs} * y);
+    }
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator*(T x, const qua<U>& ys) {
+        return qua(x * vec{ys});
+    }
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator*(const vec<T, 3>& xs, const qua<U>& ys) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/transforms/
+
+        const vec qv2 = cross(ys.v, xs) * T{2};
+        return xs + qv2 * ys.s + cross(ys.v, qv2);
+    }
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator*(const qua<T>& xs, const qua<U>& ys) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/arithmetic/
+
+        return qua{
+            cross(ys.v, xs.v) + ys.s * xs.v + xs.s * ys.v,
+            ys.s * xs.s - dot(ys.v, xs.v)};
+    }
+
+    // operator*=
+
+    template < typename T >
+    constexpr qua<T>& operator*=(qua<T>& xs, T y) {
+        return (xs = (xs * y));
+    }
+
+    template < typename T >
+    constexpr vec<T, 3>& operator*=(vec<T, 3>& xs, const qua<T>& ys) {
+        return (xs = (xs * ys));
+    }
+
+    template < typename T >
+    constexpr qua<T>& operator*=(qua<T>& xs, const qua<T>& ys) {
+        return (xs = (xs * ys));
+    }
+
+    // operator/
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator/(const qua<T>& xs, U y) {
+        return qua(vec{xs} / y);
+    }
+
+    template < typename T, typename U >
+    [[nodiscard]] constexpr auto operator/(T x, const qua<U>& ys) {
+        return qua(x / vec{ys});
+    }
+
+    // operator/=
+
+    template < typename T >
+    constexpr qua<T>& operator/=(qua<T>& xs, T y) {
+        return (xs = (xs / y));
+    }
+
+    // operator==
+
+    template < typename T >
+    [[nodiscard]] constexpr bool operator==(const qua<T>& xs, const qua<T>& ys) {
+        return vec{xs} == vec{ys};
+    }
+
+    // operator!=
+
+    template < typename T >
+    [[nodiscard]] constexpr bool operator!=(const qua<T>& xs, const qua<T>& ys) {
+        return vec{xs} != vec{ys};
+    }
+
+    // operator<
+
+    template < typename T >
+    [[nodiscard]] constexpr bool operator<(const qua<T>& xs, const qua<T>& ys) {
+        return vec{xs} < vec{ys};
+    }
+}
+
+//
+// Common Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> lerp(const qua<T>& xs, const qua<T>& ys, T a) {
+        return qua(lerp(vec{xs}, vec{ys}, a));
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> lerp(const qua<T>& xs, const qua<T>& ys, T x_a, T y_a) {
+        return qua(lerp(vec{xs}, vec{ys}, x_a, y_a));
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> nlerp(const qua<T>& unit_xs, const qua<T>& unit_ys, T a) {
+        const T xs_scale = T{1} - a;
+        const T ys_scale = a * sign(dot(unit_xs, unit_ys));
+        return normalize(lerp(unit_xs, unit_ys, xs_scale, ys_scale));
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> slerp(const qua<T>& unit_xs, const qua<T>& unit_ys, T a) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/slerp/
+
+        const T raw_cos_theta = dot(unit_xs, unit_ys);
+        const T raw_cos_theta_sign = sign(raw_cos_theta);
+
+        // half degree linear threshold: cos((pi / 180) * 0.25)
+        if ( const T cos_theta = raw_cos_theta * raw_cos_theta_sign; cos_theta < T{0.99999f} ) {
+            const T theta = acos(cos_theta);
+            const T rsin_theta = rsqrt(T{1} - sqr(cos_theta));
+            const T xs_scale = sin((T{1} - a) * theta) * rsin_theta;
+            const T ys_scale = sin(a * theta) * raw_cos_theta_sign * rsin_theta;
+            return lerp(unit_xs, unit_ys, xs_scale, ys_scale);
+        } else {
+            // use linear interpolation for small angles
+            const T xs_scale = T{1} - a;
+            const T ys_scale = a * raw_cos_theta_sign;
+            return normalize(lerp(unit_xs, unit_ys, xs_scale, ys_scale));
+        }
+    }
+}
+
+//
+// Geometric Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, typename U
+             , typename V = decltype(dot(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<U, 4>>())) >
+    [[nodiscard]] constexpr V dot(const qua<T>& xs, const qua<U>& ys) {
+        return dot(vec{xs}, vec{ys});
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr T length(const qua<T>& xs) {
+        return length(vec{xs});
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr T rlength(const qua<T>& xs) {
+        return rlength(vec{xs});
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr T length2(const qua<T>& xs) {
+        return length2(vec{xs});
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr T rlength2(const qua<T>& xs) {
+        return rlength2(vec{xs});
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr T distance(const qua<T>& xs, const qua<T>& ys) {
+        const qua zs = xs * conjugate(ys);
+        return T{2} * atan2(length(zs.v), abs(zs.s));
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> normalize(const qua<T>& xs) {
+        return qua(normalize(vec{xs}));
+    }
+}
+
+//
+// Relational Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T
+             , typename U = decltype(any(std::declval<vec<T, 4>>())) >
+    [[nodiscard]] constexpr U any(const qua<T>& xs) {
+        return any(vec{xs});
+    }
+
+    template < typename T
+             , typename U = decltype(all(std::declval<vec<T, 4>>())) >
+    [[nodiscard]] constexpr U all(const qua<T>& xs) {
+        return all(vec{xs});
+    }
+
+    template < typename T
+             , typename U = typename decltype(approx(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> approx(const qua<T>& xs, const qua<T>& ys) {
+        return approx(vec{xs}, vec{ys});
+    }
+
+    template < typename T
+             , typename U = typename decltype(approx(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>(),
+                 std::declval<T>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> approx(const qua<T>& xs, const qua<T>& ys, T epsilon) {
+        return approx(vec{xs}, vec{ys}, epsilon);
+    }
+
+    template < typename T
+             , typename U = typename decltype(less(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> less(const qua<T>& xs, const qua<T>& ys) {
+        return less(vec{xs}, vec{ys});
+    }
+
+    template < typename T
+             , typename U = typename decltype(less_equal(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> less_equal(const qua<T>& xs, const qua<T>& ys) {
+        return less_equal(vec{xs}, vec{ys});
+    }
+
+    template < typename T
+             , typename U = typename decltype(greater(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> greater(const qua<T>& xs, const qua<T>& ys) {
+        return greater(vec{xs}, vec{ys});
+    }
+
+    template < typename T
+             , typename U = typename decltype(greater_equal(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> greater_equal(const qua<T>& xs, const qua<T>& ys) {
+        return greater_equal(vec{xs}, vec{ys});
+    }
+
+    template < typename T
+             , typename U = typename decltype(equal_to(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> equal_to(const qua<T>& xs, const qua<T>& ys) {
+        return equal_to(vec{xs}, vec{ys});
+    }
+
+    template < typename T
+             , typename U = typename decltype(not_equal_to(
+                 std::declval<vec<T, 4>>(),
+                 std::declval<vec<T, 4>>()))::component_type >
+    [[nodiscard]] constexpr vec<U, 4> not_equal_to(const qua<T>& xs, const qua<T>& ys) {
+        return not_equal_to(vec{xs}, vec{ys});
+    }
+}
+
+//
+// Quaternion Functions
+//
+
+namespace vmath_hpp
+{
+    // conjugate
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> conjugate(const qua<T>& q) {
+        return {-q.v, q.s};
+    }
+
+    // inverse
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> inverse(const qua<T>& q) {
+        return conjugate(q) * rlength2(q);
+    }
+}
+
+//
+// Units
+//
+
+namespace vmath_hpp
+{
+    template < typename T > inline constexpr vec<T, 2> zero2{zero_init};
+    template < typename T > inline constexpr vec<T, 3> zero3{zero_init};
+    template < typename T > inline constexpr vec<T, 4> zero4{zero_init};
+
+    template < typename T > inline constexpr vec<T, 2> unit2{unit_init};
+    template < typename T > inline constexpr vec<T, 3> unit3{unit_init};
+    template < typename T > inline constexpr vec<T, 4> unit4{unit_init};
+
+    template < typename T > inline constexpr vec<T, 2> unit2_x{T{1}, T{0}};
+    template < typename T > inline constexpr vec<T, 2> unit2_y{T{0}, T{1}};
+
+    template < typename T > inline constexpr vec<T, 3> unit3_x{T{1}, T{0}, T{0}};
+    template < typename T > inline constexpr vec<T, 3> unit3_y{T{0}, T{1}, T{0}};
+    template < typename T > inline constexpr vec<T, 3> unit3_z{T{0}, T{0}, T{1}};
+
+    template < typename T > inline constexpr vec<T, 4> unit4_x{T{1}, T{0}, T{0}, T{0}};
+    template < typename T > inline constexpr vec<T, 4> unit4_y{T{0}, T{1}, T{0}, T{0}};
+    template < typename T > inline constexpr vec<T, 4> unit4_z{T{0}, T{0}, T{1}, T{0}};
+    template < typename T > inline constexpr vec<T, 4> unit4_w{T{0}, T{0}, T{0}, T{1}};
+
+    template < typename T > inline constexpr mat<T, 2> mzero2{zero_init};
+    template < typename T > inline constexpr mat<T, 3> mzero3{zero_init};
+    template < typename T > inline constexpr mat<T, 4> mzero4{zero_init};
+
+    template < typename T > inline constexpr mat<T, 2> munit2{unit_init};
+    template < typename T > inline constexpr mat<T, 3> munit3{unit_init};
+    template < typename T > inline constexpr mat<T, 4> munit4{unit_init};
+
+    template < typename T > inline constexpr mat<T, 2> midentity2{identity_init};
+    template < typename T > inline constexpr mat<T, 3> midentity3{identity_init};
+    template < typename T > inline constexpr mat<T, 4> midentity4{identity_init};
+
+    template < typename T > inline constexpr qua<T> qzero{zero_init};
+    template < typename T > inline constexpr qua<T> qidentity{identity_init};
+}
+
+//
+// Hash
+//
+
+namespace vmath_hpp::detail
+{
+    struct hash_combiner {
+        template < typename T >
+        [[nodiscard]] std::size_t operator()(std::size_t seed, const T& x) noexcept {
+            return (seed ^= std::hash<T>{}(x) + 0x9e3779b9 + (seed << 6) + ( seed >> 2));
+        }
+    };
+
+    template < typename T, size_t Size >
+    [[nodiscard]] std::size_t hash(const vec<T, Size>& v) noexcept {
+        return fold_join(hash_combiner{}, std::size_t{}, v);
+    }
+
+    template < typename T, size_t Size >
+    [[nodiscard]] std::size_t hash(const mat<T, Size>& m) noexcept {
+        return fold_join(hash_combiner{}, std::size_t{}, m);
+    }
+
+    template < typename T >
+    [[nodiscard]] std::size_t hash(const qua<T>& q) noexcept {
+        return fold_join(hash_combiner{}, std::size_t{}, q);
+    }
+}
+
+namespace std
+{
+    template < typename T, size_t Size >
+    struct hash<vmath_hpp::vec<T, Size>> {
+        size_t operator()(const vmath_hpp::vec<T, Size>& v) const noexcept {
+            return vmath_hpp::detail::hash(v);
+        }
+    };
+
+    template < typename T, size_t Size >
+    struct hash<vmath_hpp::mat<T, Size>> {
+        size_t operator()(const vmath_hpp::mat<T, Size>& m) const noexcept {
+            return vmath_hpp::detail::hash(m);
+        }
+    };
+
+    template < typename T >
+    struct hash<vmath_hpp::qua<T>> {
+        size_t operator()(const vmath_hpp::qua<T>& q) const noexcept {
+            return vmath_hpp::detail::hash(q);
+        }
+    };
+}
+
+//
+// Cast
+//
+
+namespace vmath_hpp
+{
+    template < typename To, typename From >
+    [[nodiscard]] std::enable_if_t<
+        std::is_arithmetic_v<To> && std::is_arithmetic_v<From>
+    , To>
+    constexpr cast_to(From x) noexcept {
+        return static_cast<To>(x);
+    }
+
+    template < typename To, typename From, std::size_t Size >
+    [[nodiscard]] constexpr vec<To, Size> cast_to(const vec<From, Size>& v) {
+        return map_join([](From x){ return cast_to<To>(x); }, v);
+    }
+
+    template < typename To, typename From, std::size_t Size >
+    [[nodiscard]] constexpr mat<To, Size> cast_to(const mat<From, Size>& m) {
+        return map_join([](const vec<From, Size>& v){ return cast_to<To>(v); }, m);
+    }
+
+    template < typename To, typename From >
+    [[nodiscard]] constexpr qua<To> cast_to(const qua<From>& q) {
+        return map_join([](From x){ return cast_to<To>(x); }, q);
+    }
+}
+
+//
+// Access
+//
+
+namespace vmath_hpp
+{
+    // component
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T component(const vec<T, Size>& v, std::size_t index) {
+        return v[index];
+    }
+
+    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;
+    }
+
+    // row
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> row(const mat<T, Size>& m, std::size_t index) {
+        return m.rows[index];
+    }
+
+    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;
+    }
+
+    // column
+
+    namespace impl
+    {
+        template < typename T, std::size_t Size, std::size_t... Is >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        vec<T, Size> column_impl(const mat<T, Size>& m, std::size_t index, std::index_sequence<Is...>) {
+            return { m[Is][index]... };
+        }
+
+        template < typename T, std::size_t Size, std::size_t... Is >
+        [[nodiscard]] constexpr VMATH_HPP_FORCE_INLINE
+        mat<T, Size> column_impl(const mat<T, Size>& m, std::size_t index, const vec<T, Size>& column, std::index_sequence<Is...>) {
+            return { component(m[Is], index, column[Is])... };
+        }
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> column(const mat<T, Size>& m, std::size_t index) {
+        return impl::column_impl(m, index, std::make_index_sequence<Size>{});
+    }
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr mat<T, Size> column(const mat<T, Size>& m, std::size_t index, const vec<T, Size>& column) {
+        return impl::column_impl(m, index, column, std::make_index_sequence<Size>{});
+    }
+
+    // real
+
+    template < typename T >
+    [[nodiscard]] constexpr T real(const qua<T>& q) {
+        return q.s;
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> real(qua<T> q, T real) {
+        q.s = real;
+        return q;
+    }
+
+    // imag
+
+    template < typename T >
+    [[nodiscard]] constexpr vec<T, 3> imag(const qua<T>& q) {
+        return q.v;
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> imag(qua<T> q, const vec<T, 3>& imag) {
+        q.v = imag;
+        return q;
+    }
+}
+
+//
+// Matrix Transform 3D
+//
+
+namespace vmath_hpp
+{
+    // trs
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> trs(const vec<T, 3>& t, const mat<T, 3>& r) {
+        return {
+            { r[0], T{0} },
+            { r[1], T{0} },
+            { r[2], T{0} },
+            { t,    T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> trs(const vec<T, 3>& t, const mat<T, 3>& r, const vec<T, 3>& s) {
+        return {
+            { r[0] * s[0], T{0} },
+            { r[1] * s[1], T{0} },
+            { r[2] * s[2], T{0} },
+            { t,           T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> trs(const vec<T, 3>& t, const qua<T>& r) {
+        return trs(t, rotate(r));
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> trs(const vec<T, 3>& t, const qua<T>& r, const vec<T, 3>& s) {
+        return trs(t, rotate(r), s);
+    }
+
+    // translate
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> translate(const vec<T, 3>& v) {
+        /// REFERENCE:
+        /// https://en.wikipedia.org/wiki/Translation_(geometry)
+
+        return {
+            { T{1}, T{0}, T{0}, T{0} },
+            { T{0}, T{1}, T{0}, T{0} },
+            { T{0}, T{0}, T{1}, T{0} },
+            { v.x,  v.y,  v.z,  T{1} }};
+    }
+
+    // rotate
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> rotate(const qua<T>& q) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/
+
+        const auto [qv, qs] = normalize(q);
+
+        const T x2 = qv.x * T{2};
+        const T y2 = qv.y * T{2};
+        const T z2 = qv.z * T{2};
+
+        const T sx2 = qs * x2;
+        const T sy2 = qs * y2;
+        const T sz2 = qs * z2;
+
+        const T xx2 = qv.x * x2;
+        const T xy2 = qv.x * y2;
+        const T xz2 = qv.x * z2;
+
+        const T yy2 = qv.y * y2;
+        const T yz2 = qv.y * z2;
+        const T zz2 = qv.z * z2;
+
+        return {
+            { T{1} - (yy2 + zz2), (xy2 + sz2),        (xz2 - sy2)        },
+            { (xy2 - sz2),        T{1} - (xx2 + zz2), (yz2 + sx2)        },
+            { (xz2 + sy2),        (yz2 - sx2),        T{1} - (xx2 + yy2) }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> rotate4(const qua<T>& q) {
+        return mat<T, 4>(rotate(q));
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> rotate(T angle, const vec<T, 3>& axis) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToMatrix/
+
+        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 xxm = xx * ic;
+        const T yym = yy * ic;
+        const T zzm = zz * ic;
+
+        const T xym = x * y * ic;
+        const T xzm = x * z * ic;
+        const T yzm = y * z * ic;
+
+        return {
+            { xxm + c,  xym + zs, xzm - ys },
+            { xym - zs, yym + c,  yzm + xs },
+            { xzm + ys, yzm - xs, zzm + c  }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> rotate4(T angle, const vec<T, 3>& axis) {
+        return mat<T, 4>(rotate(angle, axis));
+    }
+
+    // rotate_x
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> rotate_x(T angle) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/algebra/matrix/orthogonal/rotation/
+
+        const auto [s, c] = sincos(angle);
+
+        return {
+            { T{1}, T{0}, T{0} },
+            { T{0},   c,    s  },
+            { T{0},  -s,    c  }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> rotate4_x(T angle) {
+        return mat<T, 4>(rotate_x(angle));
+    }
+
+    // rotate_y
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> rotate_y(T angle) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/algebra/matrix/orthogonal/rotation/
+
+        const auto [s, c] = sincos(angle);
+
+        return {
+            {   c,  T{0},  -s  },
+            { T{0}, T{1}, T{0} },
+            {   s,  T{0},   c  }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> rotate4_y(T angle) {
+        return mat<T, 4>(rotate_y(angle));
+    }
+
+    // rotate_z
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> rotate_z(T angle) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/algebra/matrix/orthogonal/rotation/
+
+        const auto [s, c] = sincos(angle);
+
+        return {
+            {   c,    s,  T{0} },
+            {  -s,    c,  T{0} },
+            { T{0}, T{0}, T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> rotate4_z(T angle) {
+        return mat<T, 4>(rotate_z(angle));
+    }
+
+    // scale
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> scale(const vec<T, 3>& v) {
+        /// REFERENCE:
+        /// https://en.wikipedia.org/wiki/Scaling_(geometry)
+
+        return {
+            { v.x,  T{0}, T{0} },
+            { T{0}, v.y,  T{0} },
+            { T{0}, T{0}, v.z  }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> scale4(const vec<T, 3>& v) {
+        return mat<T, 4>(scale(v));
+    }
+
+    // look_at
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> look_at_lh(
+        const vec<T, 3>& dir,
+        const vec<T, 3>& up)
+    {
+        /// REFERENCE:
+        /// https://www.euclideanspace.com/maths/algebra/vectors/lookat/
+
+        const vec az = normalize(dir);
+        const vec ax = normalize(cross(up, az));
+        const vec ay = cross(az, ax);
+
+        return {
+            { ax.x, ay.x, az.x },
+            { ax.y, ay.y, az.y },
+            { ax.z, ay.z, az.z }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> look_at_lh(
+        const vec<T, 3>& eye,
+        const vec<T, 3>& at,
+        const vec<T, 3>& up)
+    {
+        /// REFERENCE:
+        /// https://www.euclideanspace.com/maths/algebra/vectors/lookat/
+
+        const vec az = normalize(at - eye);
+        const vec ax = normalize(cross(up, az));
+        const vec 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, T{0} },
+            { ax.y, ay.y, az.y, T{0} },
+            { ax.z, ay.z, az.z, T{0} },
+            { -dx,  -dy,  -dz,  T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> look_at_rh(
+        const vec<T, 3>& dir,
+        const vec<T, 3>& up)
+    {
+        /// REFERENCE:
+        /// https://www.euclideanspace.com/maths/algebra/vectors/lookat/
+
+        const vec az = normalize(-dir);
+        const vec ax = normalize(cross(up, az));
+        const vec ay = cross(az, ax);
+
+        return {
+            { ax.x, ay.x, az.x },
+            { ax.y, ay.y, az.y },
+            { ax.z, ay.z, az.z }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> look_at_rh(
+        const vec<T, 3>& eye,
+        const vec<T, 3>& at,
+        const vec<T, 3>& up)
+    {
+        /// REFERENCE:
+        /// https://www.euclideanspace.com/maths/algebra/vectors/lookat/
+
+        const vec az = normalize(eye - at);
+        const vec ax = normalize(cross(up, az));
+        const vec 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, T{0} },
+            { ax.y, ay.y, az.y, T{0} },
+            { ax.z, ay.z, az.z, T{0} },
+            { -dx,  -dy,  -dz,  T{1} }};
+    }
+}
+
+//
+// Matrix Transform 2D
+//
+
+namespace vmath_hpp
+{
+    // trs
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> trs(const vec<T, 2>& t, const mat<T, 2>& r) {
+        return {
+            { r[0], T{0} },
+            { r[1], T{0} },
+            { t,    T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> trs(const vec<T, 2>& t, const mat<T, 2>& r, const vec<T, 2>& s) {
+        return {
+            { r[0] * s[0], T{0} },
+            { r[1] * s[1], T{0} },
+            { t,           T{1} }};
+    }
+
+    // translate
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> translate(const vec<T, 2>& v) {
+        /// REFERENCE:
+        /// https://en.wikipedia.org/wiki/Translation_(geometry)
+
+        return {
+            { T{1}, T{0}, T{0} },
+            { T{0}, T{1}, T{0} },
+            { v.x,  v.y,  T{1} }};
+    }
+
+    // rotate
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 2> rotate(T angle) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/algebra/matrix/orthogonal/rotation/
+
+        const auto [s, c] = sincos(angle);
+
+        return {
+            {  c, s },
+            { -s, c }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> rotate3(T angle) {
+        return mat<T, 3>(rotate(angle));
+    }
+
+    // scale
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 2> scale(const vec<T, 2>& v) {
+        /// REFERENCE:
+        /// https://en.wikipedia.org/wiki/Scaling_(geometry)
+
+        return {
+            { v.x,  T{0} },
+            { T{0}, v.y  }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> scale3(const vec<T, 2>& v) {
+        return mat<T, 3>(scale(v));
+    }
+
+    // shear
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 2> shear(const vec<T, 2>& v) {
+        /// REFERENCE:
+        /// https://en.wikipedia.org/wiki/Shear_matrix
+
+        return {
+            { T{1}, v.y  },
+            { v.x,  T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 3> shear3(const vec<T, 2>& v) {
+        return mat<T, 3>(shear(v));
+    }
+}
+
+//
+// Matrix Projections
+//
+
+namespace vmath_hpp
+{
+    // orthographic
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> orthographic_lh(T width, T height, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixortholh
+
+        const T rwidth = rcp(width);
+        const T rheight = rcp(height);
+        const T frange = rcp(zfar - znear);
+
+        const T sx = T{2} * rwidth;
+        const T sy = T{2} * rheight;
+        const T sz = frange;
+        const T tz = -frange * znear;
+
+        return {
+            {  sx,   T{0},  T{0},  T{0} },
+            { T{0},   sy,   T{0},  T{0} },
+            { T{0},  T{0},   sz,   T{0} },
+            { T{0},  T{0},   tz,   T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> orthographic_rh(T width, T height, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixorthorh
+
+        const T rwidth = rcp(width);
+        const T rheight = rcp(height);
+        const T frange = rcp(znear - zfar);
+
+        const T sx = T{2} * rwidth;
+        const T sy = T{2} * rheight;
+        const T sz = frange;
+        const T tz = frange * znear;
+
+        return {
+            {  sx,  T{0}, T{0}, T{0} },
+            { T{0},  sy,  T{0}, T{0} },
+            { T{0}, T{0},  sz,  T{0} },
+            { T{0}, T{0},  tz,  T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> orthographic_lh(T left, T right, T bottom, T top, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixorthooffcenterlh
+
+        const T rwidth = rcp(right - left);
+        const T rheight = rcp(top - bottom);
+        const T frange = rcp(zfar - znear);
+
+        return {
+            { T{2} * rwidth,            T{0},                      T{0},            T{0} },
+            { T{0},                     T{2} * rheight,            T{0},            T{0} },
+            { T{0},                     T{0},                      frange,          T{0} },
+            { -(left + right) * rwidth, -(top + bottom) * rheight, -frange * znear, T{1} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> orthographic_rh(T left, T right, T bottom, T top, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixorthooffcenterrh
+
+        const T rwidth = rcp(right - left);
+        const T rheight = rcp(top - bottom);
+        const T frange = rcp(znear - zfar);
+
+        return {
+            { T{2} * rwidth,            T{0},                      T{0},           T{0} },
+            { T{0},                     T{2} * rheight,            T{0},           T{0} },
+            { T{0},                     T{0},                      frange,         T{0} },
+            { -(left + right) * rwidth, -(top + bottom) * rheight, frange * znear, T{1} }};
+    }
+
+    // perspective
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> perspective_lh(T width, T height, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixperspectivelh
+
+        const T sx = T{2} * znear * rcp(width);
+        const T sy = T{2} * znear * rcp(height);
+        const T sz = zfar * rcp(zfar - znear);
+        const T tz = (znear * zfar) * rcp(znear - zfar);
+
+        return {
+            { T{sx}, T{0},  T{0},  T{0} },
+            { T{0},  T{sy}, T{0},  T{0} },
+            { T{0},  T{0},  T{sz}, T{1} },
+            { T{0},  T{0},  T{tz}, T{0} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> perspective_rh(T width, T height, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixperspectiverh
+
+        const T sx = T{2} * znear * rcp(width);
+        const T sy = T{2} * znear * rcp(height);
+        const T sz = zfar * rcp(znear - zfar);
+        const T tz = (znear * zfar) * rcp(znear - zfar);
+
+        return {
+            {  sx,  T{0}, T{0},  T{0} },
+            { T{0},  sy,  T{0},  T{0} },
+            { T{0}, T{0},  sz,  -T{1} },
+            { T{0}, T{0},  tz,   T{0} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> perspective_lh(T left, T right, T bottom, T top, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixperspectiveoffcenterlh
+
+        const T znear2 = T{2} * znear;
+        const T rwidth = rcp(right - left);
+        const T rheight = rcp(top - bottom);
+        const T frange = zfar * rcp(zfar - znear);
+
+        return {
+            { znear2 * rwidth,          T{0},                      T{0},            T{0} },
+            { T{0},                     znear2 * rheight,          T{0},            T{0} },
+            { -(left + right) * rwidth, -(top + bottom) * rheight, frange,          T{1} },
+            { T{0},                     T{0},                      -frange * znear, T{0} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> perspective_rh(T left, T right, T bottom, T top, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixperspectiveoffcenterrh
+
+        const T znear2 = T{2} * znear;
+        const T rwidth = rcp(right - left);
+        const T rheight = rcp(top - bottom);
+        const T frange = zfar * rcp(znear - zfar);
+
+        return {
+            { znear2 * rwidth,         T{0},                     T{0},            T{0} },
+            { T{0},                    znear2 * rheight,         T{0},            T{0} },
+            { (left + right) * rwidth, (top + bottom) * rheight, frange,         -T{1} },
+            { T{0},                    T{0},                     frange * znear,  T{0} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> perspective_fov_lh(T fovy, T aspect, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixperspectivefovlh
+
+        const T sy = rcp(tan(fovy * T{0.5f}));
+        const T sx = sy * rcp(aspect);
+        const T sz = zfar * rcp(zfar - znear);
+        const T tz = (znear * zfar) * rcp(znear - zfar);
+
+        return {
+            {  sx,  T{0}, T{0}, T{0} },
+            { T{0},  sy,  T{0}, T{0} },
+            { T{0}, T{0},  sz,  T{1} },
+            { T{0}, T{0},  tz,  T{0} }};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr mat<T, 4> perspective_fov_rh(T fovy, T aspect, T znear, T zfar) {
+        /// REFERENCE:
+        /// https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixperspectivefovrh
+
+        const T sy = rcp(tan(fovy * T{0.5f}));
+        const T sx = sy * rcp(aspect);
+        const T sz = zfar * rcp(znear - zfar);
+        const T tz = (znear * zfar) * rcp(znear - zfar);
+        return  {
+            {  sx,  T{0}, T{0},  T{0} },
+            { T{0},  sy,  T{0},  T{0} },
+            { T{0}, T{0},  sz,  -T{1} },
+            { T{0}, T{0},  tz,   T{0} }};
+    }
+}
+
+//
+// Vector Transform
+//
+
+namespace vmath_hpp
+{
+    // angle
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr T angle(const vec<T, Size>& x, const vec<T, Size>& y) {
+        const T rs = rsqrt(length2(x) * length2(y));
+        return acos(clamp(dot(x, y) * rs, -T{1}, T{1}));
+    }
+
+    // rotate
+
+    template < typename T >
+    [[nodiscard]] constexpr vec<T, 2> rotate(const vec<T, 2>& v, T angle) {
+        return v * rotate(angle);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr vec<T, 3> rotate_x(const vec<T, 3>& v, T angle) {
+        return v * qrotate(angle, unit3_x<T>);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr vec<T, 3> rotate_y(const vec<T, 3>& v, T angle) {
+        return v * qrotate(angle, unit3_y<T>);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr vec<T, 3> rotate_z(const vec<T, 3>& v, T angle) {
+        return v * qrotate(angle, unit3_z<T>);
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr vec<T, 3> rotate(const vec<T, 3>& v, T angle, const vec<T, 3>& axis) {
+        return v * qrotate(angle, axis);
+    }
+
+    // project
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> project(const vec<T, Size>& v, const vec<T, Size>& normal) {
+        return dot(v, normal) * rlength2(normal) * normal;
+    }
+
+    // perpendicular
+
+    template < typename T, std::size_t Size >
+    [[nodiscard]] constexpr vec<T, Size> perpendicular(const vec<T, Size>& v, const vec<T, Size>& normal) {
+        return v - project(v, normal);
+    }
+}
+
+//
+// Quaternion Transform
+//
+
+namespace vmath_hpp
+{
+    // qrotate
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> qrotate(const mat<T, 3>& m) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
+
+        auto xyzw = T{0.5f} * sqrt(max(T{0}, vec{
+            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]}));
+
+        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 >
+    [[nodiscard]] constexpr qua<T> qrotate(const vec<T, 3>& from, const vec<T, 3>& to) {
+        /// REFERENCE:
+        /// http://lolengine.net/blog/2014/02/24/quaternion-from-two-vectors-final
+
+        const T n = sqrt(length2(from) * length2(to));
+        const T s = dot(from, to) + n;
+
+        if ( s < T{0.000001f} * n ) {
+            return abs(from.z) < abs(from.x)
+                ? normalize(qua{vec{-from.y, from.x, T{0}}, T{0}})
+                : normalize(qua{vec{T{0}, -from.z, from.y}, T{0}});
+        }
+
+        return normalize(qua{cross(from, to), s});
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> qrotate(T angle, const vec<T, 3>& axis) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToQuaternion/
+
+        const auto [s, c] = sincos(angle * T{0.5f});
+        const auto [x, y, z] = normalize(axis);
+
+        return {vec{x, y, z} * s, c};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> qrotate_x(T angle) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToQuaternion/
+
+        const auto [s, c] = sincos(angle * T{0.5f});
+
+        return {vec{s, T{0}, T{0}}, c};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> qrotate_y(T angle) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToQuaternion/
+
+        const auto [s, c] = sincos(angle * T{0.5f});
+
+        return {vec{T{0}, s, T{0}}, c};
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> qrotate_z(T angle) {
+        /// REFERENCE:
+        /// http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToQuaternion/
+
+        const auto [s, c] = sincos(angle * T{0.5f});
+
+        return {vec{T{0}, T{0}, s}, c};
+    }
+
+    // qlook_at
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> qlook_at_lh(const vec<T, 3>& dir, const vec<T, 3>& up) {
+        return qrotate(look_at_lh(dir, up));
+    }
+
+    template < typename T >
+    [[nodiscard]] constexpr qua<T> qlook_at_rh(const vec<T, 3>& dir, const vec<T, 3>& up) {
+        return qrotate(look_at_rh(dir, up));
+    }
+}
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