diff --git a/engine/native/core/math/constants.cppm b/engine/native/core/math/constants.cppm
index aba5794..241eff6 100644
--- a/engine/native/core/math/constants.cppm
+++ b/engine/native/core/math/constants.cppm
@@ -1,9 +1,7 @@
 module;
 
-#include <concepts>
-#include <cstdint>
-#include <limits>
 #include <numbers>
+#include <limits>
 
 export module core.math.constants;
 import core.defs;
@@ -12,25 +10,27 @@ export namespace draco::math {
     // Limit the depth of recursive algorithms
     constexpr int MAX_RECURSIONS = 100;
 
-    constexpr double SQRT2 = std::numbers::sqrt2_v<double>;
-    constexpr double SQRT3 = std::numbers::sqrt3_v<double>;
-    constexpr double SQRT12 = 1. / SQRT2;
-    constexpr double SQRT13 = std::numbers::inv_sqrt3_v<double>;
-    constexpr double LN2 = std::numbers::ln2_v<double>;
-    constexpr double LN10 = std::numbers::ln10_v<double>;
-    constexpr double PI = std::numbers::pi_v<double>;
-    constexpr double TAU = 2. * PI;
-    constexpr double E = std::numbers::e_v<double>;
-    constexpr double INF = std::numeric_limits<double>::infinity();
-    constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
-    constexpr double DB_CONVERSION_GAIN = 8.6858896380650365530225783783321;
-    constexpr double GAIN_CONVERSION_DB = 0.11512925464970228420089957273422;
-    constexpr double UINT32_MAX_D = 1. / std::numeric_limits<uint32_t>::max();
-    constexpr float  UINT32_MAX_F = 1.f / std::numeric_limits<uint32_t>::max();
+    constexpr float SQRT2 = std::numbers::sqrt2_v<float>;
+    constexpr float SQRT3 = std::numbers::sqrt3_v<float>;
+    constexpr float SQRT12 = 1. / SQRT2;
+    constexpr float SQRT13 = std::numbers::inv_sqrt3_v<float>;
+    constexpr float LN2 = std::numbers::ln2_v<float>;
+    constexpr float LN10 = std::numbers::ln10_v<float>;
+    constexpr float PI = std::numbers::pi_v<float>;
+    constexpr float PI2 = PI * .5;
+    constexpr float TAU = 2. * PI;
+    constexpr float E = std::numbers::e_v<float>;
+    constexpr float INF = std::numeric_limits<float>::infinity();
+    constexpr float NaN = std::numeric_limits<float>::quiet_NaN();
+    constexpr float DB_CONVERSION_GAIN = 8.6858896380650365530225783783321;
+    constexpr float GAIN_CONVERSION_DB = 0.11512925464970228420089957273422;
+    constexpr float UINT32_MAX_F = 1.f / std::numeric_limits<uint32_t>::max();
+    constexpr float DECIMAL_LIMIT_F = 8388608.0f;
 
-    template<std::floating_point T> constexpr T CMP_EPSILON = T{0.000001};
-    template<std::floating_point T> constexpr T CMP_EPSILON2 = CMP_EPSILON<T> * CMP_EPSILON<T>;
+    constexpr float CMP_EPSILON = 0.000001f;
+    constexpr float CMP_EPSILON2 = CMP_EPSILON * CMP_EPSILON;
 
-    template<std::floating_point T> constexpr T CMP_NORMALIZE_TOLERANCE = T{0.000001};
-    template<std::floating_point T> constexpr T CMP_POINT_IN_PLANE_EPSILON = T{0.00001};
+    constexpr float CMP_NORMALIZE_TOLERANCE = 0.000001f;
+    constexpr float CMP_NORMALIZE_TOLERANCE2 = CMP_NORMALIZE_TOLERANCE * CMP_NORMALIZE_TOLERANCE;
+    constexpr float CMP_POINT_IN_PLANE_EPSILON = 0.00001f;
 }
diff --git a/engine/native/core/math/functions.cppm b/engine/native/core/math/functions.cppm
new file mode 100644
index 0000000..f0a69fc
--- /dev/null
+++ b/engine/native/core/math/functions.cppm
@@ -0,0 +1,195 @@
+module;
+
+#include <numbers>
+#include <cmath>
+#include <concepts>
+#include <limits>
+
+export module core.math.functions;
+import core.math.constants;
+import core.defs;
+
+export namespace draco::math {
+    template <arithmetic T>
+    constexpr T sqr(T x) noexcept { return x*x; }
+
+    template <std::floating_point T>
+    [[nodiscard]] constexpr bool is_nan(T val) noexcept {
+        // Only NaN does not equal itself.
+        return val != val;
+    }
+
+    template <std::floating_point T>
+    [[nodiscard]] constexpr bool is_inf(T val) noexcept {
+        return std::isinf(val);
+    }
+
+    template <std::floating_point T>
+    [[nodiscard]] constexpr bool is_finite(T val) noexcept {
+        return std::isfinite(val);
+    }
+
+    template <arithmetic T>
+    constexpr T abs(T value) noexcept {
+        // Manually compute abs for signed types.
+        // Also avoids potential int8_t -> int issues.
+        if constexpr (std::floating_point<T>) {
+            return value < T{0} ? -value : value;
+        } else if constexpr (std::signed_integral<T>) {
+            if (value == std::numeric_limits<T>::min()) {
+                return std::numeric_limits<T>::max(); // define saturating behavior explicitly
+            }
+            return value < T{0} ? -value : value;
+        } else {
+            // unsigned is always positive! :^)
+            return value;
+        }
+    }
+
+	template <arithmetic T>
+	constexpr T sign(T value) noexcept {
+        if constexpr (std::floating_point<T>) {
+			if (value != value) {
+				return value;
+			} else if (value) {
+            	return value < T{0} ? T{-1} : T{1};
+			}
+			return T{0};
+        } else if constexpr (std::signed_integral<T>) {
+            if (value) {
+            	return value < T{0} ? T{-1} : T{1};
+            }
+			return T{0};
+        } else {
+            return T{value != T{0}};
+        }
+	}
+
+	constexpr float floor(float value) noexcept {
+		if (value != value || abs(value) >= DECIMAL_LIMIT_F) {
+			return value;
+		}
+		const float truncated = static_cast<int>(value);
+		return truncated - (value < truncated);
+	}
+
+	constexpr float ceil(float value) noexcept {
+		return -floor(-value);
+	}
+
+	constexpr float trunc(float value) noexcept {
+		if (value != value || abs(value) >= DECIMAL_LIMIT_F) {
+			return value;
+		}
+		return static_cast<int>(value);
+	}
+
+	constexpr float round(float value) noexcept {
+		const float s = sign(value);
+		return s * floor(s * value + 0.5f);
+	}
+
+    template <std::floating_point T>
+    constexpr T deg_to_rad(T y) noexcept {
+        return y * (T{PI} / T{180.});
+    }
+
+    template <std::floating_point T>
+    constexpr T rad_to_deg(T y) noexcept {
+        return y * (T{180.} / T{PI});
+    }
+
+    template <std::floating_point T>
+    T pow(T x, T y) {
+        return static_cast<T>(std::pow(x, y));
+    }
+
+    template <std::floating_point T>
+    constexpr T lerp(T from, T to, T weight) noexcept {
+        return std::lerp(from, to, weight);
+    }
+
+    template <std::floating_point T>
+    constexpr T cubic_interpolate(T from, T to, T before, T after, T weight) noexcept {
+        // weight squared.
+        T w2 = weight * weight;
+        // weight cubed.
+        T w3 = weight * w2;
+
+        // calculate coefficients.
+        T a = -before + to;
+        T b = T{2} * before - T{5} * from + T{4} * to - after;
+        T c = -before + T{3} * from - T{3} * to + after;
+
+        // Catmull-Rom Interpolation:
+        // 0.5 * ((2 * p_from) + (a * w) + (b * w^2) + (c * w^3))
+
+        if consteval {
+            // compile time
+            return T{0.5} * (T{2.}*from + a*weight + b*w2 + c*w3);
+        } else {
+            // runtime
+            return T{0.5} * std::fma(c, w3, std::fma(b, w2, std::fma(a, weight, T{2} * from)));
+        }
+    }
+
+    template <std::floating_point T>
+    constexpr T cubic_interpolate_in_time(
+        T from, T to,
+        T before, T after, T weight,
+        T to_t, T before_t, T after_t) noexcept {
+        /* Barry-Goldman method */
+        T t = lerp(T{0.}, to_t, weight);
+
+        // At least try to make this easier to parse for others.
+        T pre_scale = before_t == T{0.} ? T{0.} : (t - before_t) / -before_t;
+        T to_scale = (to_t == T{0.}) ? T{.5} : t / to_t;
+        T post_range = after_t - to_t;
+        T post_scale = (post_range == T{0.}) ? T{1.} : (t - to_t) / post_range;
+
+        // First layer.
+        T a1 = lerp(before, from, pre_scale);
+        T a2 = lerp(from, to, to_scale);
+        T a3 = lerp(to, after, post_scale);
+
+        // More parsing.
+        T mid_range = to_t - before_t;
+        T from_to_scale = (mid_range == T{0.}) ? T{0.} : (t - before_t) / mid_range;
+        T to_post_scale = (after_t == T{0.}) ? T{1.} : t / after_t;
+
+        // Second layer.
+        T b1 = lerp(a1, a2, from_to_scale);
+        T b2 = lerp(a2, a3, to_post_scale);
+
+        // One more for the road.
+        T final_scale = (to_t == T{0.}) ? T{.5} : t / to_t;
+
+        return lerp(b1, b2, final_scale);
+    }
+
+    template <std::floating_point T>
+    constexpr T bezier_interpolate(T start, T control_1, T control_2, T end, T t) noexcept {
+        /* Formula from Wikipedia article on Bezier curves. */
+        // one minus t.
+        T omt = T{1.} - t;
+        T omt2 = omt * omt;
+        T omt3 = omt2 * omt;
+        T t2 = t * t;
+        T t3 = t2 * t;
+
+        // B(t) = (1-t)^3 * P_0 + 3(1 - t)^2 * t * P_1 + 3(1 - t) * t^2 * P_2 + t^3 * P_3
+        T d = start * omt3 + control_1 * omt2 * t * T{3.} + control_2 * omt * t2 * T{3.} + end * t3;
+        return d;
+    }
+
+    template <std::floating_point T>
+    constexpr T bezier_derivative(T start, T control_1, T control_2, T end, T t) noexcept {
+        /* Formula from Wikipedia article on Bezier curves. */
+        T omt = T{1.} - t;
+        T omt2 = omt * omt;
+        T t2 = t * t;
+
+        T d = (control_1 - start) * T{3.} * omt2 + (control_2 - control_1) * T{6.} * omt * t + (end - control_2) * T{3.} * t2;
+        return d;
+    }
+}
\ No newline at end of file
diff --git a/engine/native/core/math/math.cppm b/engine/native/core/math/math.cppm
index 8134cf8..fe3de8f 100644
--- a/engine/native/core/math/math.cppm
+++ b/engine/native/core/math/math.cppm
@@ -1,153 +1,6 @@
-module;
-
-#include <type_traits>
-#include <concepts>
-#include <limits>
-#include <cmath>
-
 export module core.math;
+
 export import core.math.constants;
-export import core.math.vector4;
+export import core.math.functions;
+export import core.math.types;
 export import core.defs;
-
-export namespace draco::math {
-    template <arithmetic T>
-    [[nodiscard]] constexpr T sqr(T x) noexcept { return x*x; }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr bool is_nan(T val) noexcept {
-        // Only NaN does not equal itself.
-        return val != val;
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr bool is_inf(T val) noexcept {
-        return std::isinf(val);
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr bool is_finite(T val) noexcept {
-        return std::isfinite(val);
-    }
-
-    template <arithmetic T>
-    [[nodiscard]] constexpr T abs(T value) noexcept {
-        // Manually compute abs for signed types.
-        // Also avoids potential int8_t -> int issues.
-        if constexpr (std::floating_point<T>) {
-            return value < T{0} ? -value : value;
-        } else if constexpr (std::signed_integral<T>) {
-            if (value == std::numeric_limits<T>::min()) {
-                return std::numeric_limits<T>::max(); // define saturating behavior explicitly
-            }
-            return value < T{0} ? -value : value;
-        } else {
-            // unsigned is always positive! :^)
-            return value;
-        }
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr T deg_to_rad(T y) noexcept {
-        return y * (T{PI} / T{180.});
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr T rad_to_deg(T y) noexcept {
-        return y * (T{180.} / T{PI});
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] T pow(T x, T y) noexcept {
-        return static_cast<T>(std::pow(x, y));
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr T lerp(T from, T to, T weight) noexcept {
-        return std::lerp(from, to, weight);
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr T cubic_interpolate(T from, T to, T before, T after, T weight) noexcept {
-        // weight squared.
-        T w2 = weight * weight;
-        // weight cubed.
-        T w3 = weight * w2;
-
-        // calculate coefficients.
-        T a = -before + to;
-        T b = T{2} * before - T{5} * from + T{4} * to - after;
-        T c = -before + T{3} * from - T{3} * to + after;
-
-        // Catmull-Rom Interpolation:
-        // 0.5 * ((2 * p_from) + (a * w) + (b * w^2) + (c * w^3))
-
-        if consteval {
-            // compile time
-            return T{0.5} * (T{2.}*from + a*weight + b*w2 + c*w3);
-        } else {
-            // runtime
-            return T{0.5} * std::fma(c, w3, std::fma(b, w2, std::fma(a, weight, T{2} * from)));
-        }
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr T cubic_interpolate_in_time(
-        T from, T to,
-        T before, T after, T weight,
-        T to_t, T before_t, T after_t) noexcept {
-        /* Barry-Goldman method */
-        T t = lerp(T{0.}, to_t, weight);
-
-        // At least try to make this easier to parse for others.
-        T pre_scale = before_t == T{0.} ? T{0.} : (t - before_t) / -before_t;
-        T to_scale = (to_t == T{0.}) ? T{.5} : t / to_t;
-        T post_range = after_t - to_t;
-        T post_scale = (post_range == T{0.}) ? T{1.} : (t - to_t) / post_range;
-
-        // First layer.
-        T a1 = lerp(before, from, pre_scale);
-        T a2 = lerp(from, to, to_scale);
-        T a3 = lerp(to, after, post_scale);
-
-        // More parsing.
-        T mid_range = to_t - before_t;
-        T from_to_scale = (mid_range == T{0.}) ? T{0.} : (t - before_t) / mid_range;
-        T to_post_scale = (after_t == T{0.}) ? T{1.} : t / after_t;
-
-        // Second layer.
-        T b1 = lerp(a1, a2, from_to_scale);
-        T b2 = lerp(a2, a3, to_post_scale);
-
-        // One more for the road.
-        T final_scale = (to_t == T{0.}) ? T{.5} : t / to_t;
-
-        return lerp(b1, b2, final_scale);
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr T bezier_interpolate(T start, T control_1, T control_2, T end, T t) noexcept {
-        /* Formula from Wikipedia article on Bezier curves. */
-        // one minus t.
-        T omt = T{1.} - t;
-        T omt2 = omt * omt;
-        T omt3 = omt2 * omt;
-        T t2 = t * t;
-        T t3 = t2 * t;
-
-        // B(t) = (1-t)^3 * P_0 + 3(1 - t)^2 * t * P_1 + 3(1 - t) * t^2 * P_2 + t^3 * P_3
-        T d = start * omt3 + control_1 * omt2 * t * T{3.} + control_2 * omt * t2 * T{3.} + end * t3;
-        return d;
-    }
-
-    template <std::floating_point T>
-    [[nodiscard]] constexpr T bezier_derivative(T start, T control_1, T control_2, T end, T t) noexcept {
-        /* Formula from Wikipedia article on Bezier curves. */
-        T omt = T{1.} - t;
-        T omt2 = omt * omt;
-        T t2 = t * t;
-
-        T d = (control_1 - start) * T{3.} * omt2 + (control_2 - control_1) * T{6.} * omt * t + (end - control_2) * T{3.} * t2;
-        return d;
-    }
-}
\ No newline at end of file
diff --git a/engine/native/core/math/math.test.cpp b/engine/native/core/math/math.test.cpp
index deca5ed..d538b82 100644
--- a/engine/native/core/math/math.test.cpp
+++ b/engine/native/core/math/math.test.cpp
@@ -1,86 +1,1448 @@
 #define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
 #include <doctest.h>
+#include <type_traits>
 
-import core.math;
-
-TEST_CASE("pow") {
-    double result = draco::math::pow(2., .5);
-    constexpr double expected = draco::math::SQRT2;
-    REQUIRE(result == expected);
+#define R_CHECK_EQ(L_expr, R_expr) {                                 \
+    const auto L_run = (L_expr);                                     \
+    const auto R_run = (R_expr);                                     \
+    static_assert(std::is_same_v<decltype(L_run), decltype(R_run)>); \
+    CHECK_EQ(L_run, R_run);                                          \
 }
 
-TEST_CASE("abs") {
-    using draco::math::abs;
-
-    REQUIRE(abs(-1.f) == 1.f);
-    REQUIRE(abs(4.56f) == 4.56f);
-    REQUIRE(abs(-1.) == 1.);
-    REQUIRE(abs(4.56) == 4.56);
-    REQUIRE(abs(-5) == 5);
-    REQUIRE(abs(3L) == 3L);
-    REQUIRE(abs(-32L) == 32L);
-    REQUIRE(abs(5000ULL) == 5000ULL);
+#define RAC_CHECK_EQ(L_expr, R_expr) {                                 \
+    static constexpr auto L_comp = (L_expr);                           \
+    static constexpr auto R_comp = (R_expr);                           \
+    static_assert(std::is_same_v<decltype(L_comp), decltype(R_comp)>); \
+    static_assert(L_comp == R_comp);                                   \
+    R_CHECK_EQ(L_expr, R_expr);                                        \
 }
 
-TEST_CASE("Vector4 ctor and dtor") {
-    using draco::math::Vector4;
+#define BASIC_R_SUBCASE(name, L_expr, R_expr) \
+    SUBCASE(name) { R_CHECK_EQ(L_expr, R_expr); }
 
-    static constexpr Vector4 v{1.0f, 2.0f, 3.0f, 4.0f};
-    static_assert(v[0] == 1.0f);
-    static_assert(v[1] == 2.0f);
-    static_assert(v[2] == 3.0f);
-    static_assert(v[3] == 4.0f);
+#define BASIC_R_SUBCASE_2(name, L_expr1, R_expr1, L_expr2, R_expr2) \
+    SUBCASE(name) { R_CHECK_EQ(L_expr1, R_expr1); R_CHECK_EQ(L_expr2, R_expr2); }
 
-    REQUIRE(v[0] == 1.0f);
-    REQUIRE(v[1] == 2.0f);
-    REQUIRE(v[2] == 3.0f);
-    REQUIRE(v[3] == 4.0f);
-}
+#define BASIC_RAC_SUBCASE(name, L_expr, R_expr) \
+    SUBCASE(name) { RAC_CHECK_EQ(L_expr, R_expr); }
 
-TEST_CASE("Vector4 swap") {
-    using draco::math::Vector4;
+#define BASIC_RAC_SUBCASE_2(name, L_expr1, R_expr1, L_expr2, R_expr2) \
+    SUBCASE(name) { RAC_CHECK_EQ(L_expr1, R_expr1); RAC_CHECK_EQ(L_expr2, R_expr2); }
+
+import core.math;
 
-    Vector4 a{1.f, 2.f, 3.f, 4.f};
-    Vector4 b{4.f, 3.f, 2.f, 1.f};
+TEST_SUITE("math") {
+    TEST_CASE("pow") {
+        float result = draco::math::pow(2.0f, 0.5f);
+        constexpr float expected = draco::math::SQRT2;
+        CHECK_EQ(result, expected);
+    }
 
-    std::swap(a, b);
+    TEST_CASE("abs") {
+        using draco::math::abs;
 
-    REQUIRE(a == Vector4{4.f, 3.f, 2.f, 1.f});
-    REQUIRE(b == Vector4{1.f, 2.f, 3.f, 4.f});
+        RAC_CHECK_EQ(abs(-1.f), 1.f);
+        RAC_CHECK_EQ(abs(4.56f), 4.56f);
+        RAC_CHECK_EQ(abs(-1.), 1.);
+        RAC_CHECK_EQ(abs(4.56), 4.56);
+        RAC_CHECK_EQ(abs(-5), 5);
+        RAC_CHECK_EQ(abs(3L), 3L);
+        RAC_CHECK_EQ(abs(-32L), 32L);
+        RAC_CHECK_EQ(abs(5000ULL), 5000ULL);
+    }
 }
 
-TEST_CASE("Vector4 dot") {
-    using draco::math::Vector4;
-    using draco::math::dot;
+TEST_SUITE("vector2") {
+    TEST_CASE("constructors") {
+        using draco::math::Vector2;
+        using draco::math::Vector3;
+        using draco::math::Vector4;
 
-    static constexpr Vector4 a{1.0f, 2.0f, 3.0f, 4.0f};
-    static constexpr Vector4 b{5.0f, 6.0f, 7.0f, 8.0f};
+        static constexpr Vector3 a{1.0f, 2.0f, 3.0f};
+        static constexpr Vector4 b{4.0f, 5.0f, 6.0f, 7.0f};
+        
+        BASIC_RAC_SUBCASE("float",
+            ( Vector2(1.0f) ),
+            ( Vector2{1.0f, 1.0f} )
+        );
+        
+        BASIC_RAC_SUBCASE("vec3",
+            ( Vector2(a) ),
+            ( Vector2{1.0f, 2.0f} )
+        );
+        
+        BASIC_RAC_SUBCASE("vec4",
+            ( Vector2(b) ),
+            ( Vector2{4.0f, 5.0f} )
+        );
+    }
 
-    const float result = dot(a, b);
-    // 1 * 5 + 2 * 6 + 3 * 7 + 4 * 8
-    static constexpr float expected = 70.0f;
+    TEST_CASE("access") {
+        using draco::math::Vector2;
 
-    REQUIRE(result == expected);
-}
+        static constexpr Vector2 v(1.0f, 2.0f);
+
+        RAC_CHECK_EQ(v[0], 1.0f);
+        RAC_CHECK_EQ(v[1], 2.0f);
+    }
+
+    TEST_CASE("swizzle") {
+        using draco::math::Vector2;
+        using draco::math::Vector3;
+        using draco::math::Vector4;
+        
+        static constexpr Vector2 v{1.0f, 2.0f};
+
+        BASIC_RAC_SUBCASE("vec2",
+            ( v[1, 0] ),
+            ( Vector2{2.0f, 1.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec3",
+            ( v[1, 1, 0] ),
+            ( Vector3{2.0f, 2.0f, 1.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec4",
+            ( v[0, 1, 1, 0] ),
+            ( Vector4{1.0f, 2.0f, 2.0f, 1.0f} )
+        );
+    }
+
+    TEST_CASE("swap") {
+        using draco::math::Vector2;
+
+        Vector2 a{1.f, 2.f};
+        Vector2 b{2.f, 1.f};
+
+        std::swap(a, b);
+
+        CHECK_EQ(a, Vector2{2.f, 1.f});
+        CHECK_EQ(b, Vector2{1.f, 2.f});
+    }
+
+    TEST_CASE("dot") {
+        using draco::math::Vector2;
+        using draco::math::dot;
+
+        static constexpr Vector2 a{1.0f, 2.0f};
+        static constexpr Vector2 b{3.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE("basic",
+            ( dot(a, b) ),
+            ( 11.0f )
+        );
+
+        BASIC_RAC_SUBCASE("self",
+            ( dot(a, a) ),
+            ( 5.0f )
+        );
+
+        BASIC_RAC_SUBCASE("zero",
+            ( dot(a, Vector2()) ),
+            ( 0.0f )
+        );
+    }
+
+    TEST_CASE("length") {
+        using draco::math::Vector2;
+        using draco::math::length;
+        using draco::math::length_sq;
+
+        static constexpr Vector2 v{3.0f, 4.0f};
+
+        BASIC_R_SUBCASE("normal",
+            ( length(v) ),
+            ( 5.0f )
+        );
+
+        BASIC_RAC_SUBCASE("squared",
+            ( length_sq(v) ),
+            ( 25.0f )
+        );
+    }
+
+    TEST_CASE("distance") {
+        using draco::math::Vector2;
+        using draco::math::distance;
+        using draco::math::distance_sq;
+
+        static constexpr Vector2 a{3.0f, 4.0f};
+        static constexpr Vector2 b{-3.0f, 12.0f};
+
+        BASIC_R_SUBCASE("normal",
+            ( distance(a, b) ),
+            ( 10.0f )
+        );
+
+        BASIC_RAC_SUBCASE("squared",
+            ( distance_sq(a, b) ),
+            ( 100.0f )
+        );
+    }
+
+    TEST_CASE("normalize") {
+        using draco::math::Vector2;
+        using draco::math::length;
+        using draco::math::normalize;
+        using draco::math::normalize_fast;
+
+        static constexpr Vector2 a{3.0f, 4.0f};
+        static constexpr Vector2 b(1e-99);
+
+        const Vector2 result = normalize(a);
+        const Vector2 result_fast = normalize_fast(a);
+        const Vector2 result_zero = normalize(b);
+
+        CHECK_EQ(length(result), 1.0f);
+        CHECK_EQ(result, result_fast);
+        CHECK_EQ(result_zero, Vector2());
+    }
+
+    TEST_CASE("project") {
+        using draco::math::Vector2;
+        using draco::math::project;
+
+        static constexpr Vector2 a{4.0f, 6.0f};
+        static constexpr Vector2 b{2.0f, 2.0f};
+
+        RAC_CHECK_EQ(
+            ( project(a, b) ),
+            ( Vector2{5.0f, 5.0f} )
+        );
+    }
+
+    TEST_CASE("reflect") {
+        using draco::math::Vector2;
+        using draco::math::reflect;
+
+        static constexpr Vector2 a{1.0f, 2.0f};
+        static constexpr Vector2 b{3.0f, 4.0f};
+
+        RAC_CHECK_EQ(
+            ( reflect(a, b) ),
+            ( Vector2{-65.0f, -86.0f} )
+        );
+    }
+
+    TEST_CASE("angle") {
+        using draco::math::Vector2;
+        using draco::math::angle;
+        using draco::math::PI2;
+
+        static constexpr Vector2 a{2.0f, 1.0f};
+        static constexpr Vector2 b{-2.0f, 4.0f};
+
+        R_CHECK_EQ(
+            ( angle(a, b) ),
+            ( PI2 )
+        );
+    }
+
+    TEST_CASE("lerp") {
+        using draco::math::Vector2;
+        using draco::math::lerp;
+
+        static constexpr Vector2 a{1.0f, 2.0f};
+        static constexpr Vector2 b{3.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE("weight = -1",
+            ( lerp(a, b, -1.0f) ),
+            ( Vector2{-1.0f, 0.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 0",
+            ( lerp(a, b, 0.0f) ),
+            ( a )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 0.5",
+            ( lerp(a, b, 0.5f) ),
+            ( Vector2{2.0f, 3.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 1",
+            ( lerp(a, b, 1.0f) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 2",
+            ( lerp(a, b, 2.0f) ),
+            ( Vector2{5.0f, 6.0f} )
+        );
+    }
+
+    TEST_CASE("min") {
+        using draco::math::Vector2;
+        using draco::math::min;
+
+        static constexpr Vector2 a{5.0f, 3.0f};
+        static constexpr Vector2 b{1.0f, 7.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( min(a, b) ),
+            ( Vector2{1.0f, 3.0f} )
+        );
+
+        static constexpr Vector2 expected{4.0f, 3.0f};
+
+        BASIC_RAC_SUBCASE_2("float",
+            ( min(a, 4.0f) ),
+            ( expected ),
+            ( min(4.0f, a) ),
+            ( expected )
+        );
+    }
+
+    TEST_CASE("min_length") {
+        using draco::math::Vector2;
+        using draco::math::length;
+        using draco::math::min_length;
+
+        static constexpr Vector2 a{3.0f, 4.0f};  // len: 5
+        static constexpr Vector2 b{5.0f, 12.0f}; // len: 13
+
+        BASIC_RAC_SUBCASE("vector",
+            ( min_length(a, b) ),
+            ( a )
+        );
+
+        SUBCASE("float") {
+            static constexpr float smaller_length = 1.0f;
+            static constexpr float larger_length = 10.0f;
+
+            const Vector2 result_smaller = min_length(a, smaller_length);
+            const Vector2 result_swapped = min_length(smaller_length, a);
+            const Vector2 result_larger = min_length(larger_length, a);
+
+            CHECK_EQ(length(result_smaller), smaller_length);
+            CHECK_EQ(result_smaller, result_swapped);
+            CHECK_EQ(result_larger, a);
+        }
+    }
+
+    TEST_CASE("max") {
+        using draco::math::Vector2;
+        using draco::math::max;
+
+        static constexpr Vector2 a{5.0f, 3.0f};
+        static constexpr Vector2 b{1.0f, 7.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( max(a, b) ),
+            ( Vector2{5.0f, 7.0f} )
+        );
+
+        static constexpr Vector2 expected{5.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE_2("float",
+            ( max(a, 4.0f) ),
+            ( expected ),
+            ( max(4.0f, a) ),
+            ( expected )
+        );
+    }
+
+    TEST_CASE("max_length") {
+        using draco::math::Vector2;
+        using draco::math::length;
+        using draco::math::max_length;
+
+        static constexpr Vector2 a{3.0f, 4.0f};  // len: 5
+        static constexpr Vector2 b{5.0f, 12.0f}; // len: 13
+
+        BASIC_RAC_SUBCASE("vector",
+            ( max_length(a, b) ),
+            ( b )
+        );
+
+        SUBCASE("float") {
+            static constexpr float smaller_length = 1.0f;
+            static constexpr float larger_length = 10.0f;
+
+            const Vector2 result_smaller = max_length(a, smaller_length);
+            const Vector2 result_swapped = max_length(larger_length, a);
+            const Vector2 result_larger = max_length(larger_length, a);
+
+            CHECK_EQ(length(result_larger), larger_length);
+            CHECK_EQ(result_larger, result_swapped);
+            CHECK_EQ(result_smaller, a);
+        }
+    }
 
-TEST_CASE("Vector4 dot zero") {
-    using draco::math::Vector4;
-    using draco::math::dot;
+    TEST_CASE("clamp") {
+        using draco::math::Vector2;
+        using draco::math::clamp;
 
-    Vector4 a{0.0f, 0.0f, 0.0f, 0.0f};
-    Vector4 b{1.0f, 2.0f, 3.0f, 4.0f};
+        static constexpr Vector2 a{5.0f, 3.0f};
+        static constexpr Vector2 b{1.0f, 7.0f};
+        static constexpr Vector2 c{3.0f, 8.0f};
 
-    REQUIRE(dot(a, b) == 0.0f);
+        BASIC_RAC_SUBCASE("vector",
+            ( clamp(a, b, c) ),
+            ( Vector2{3.0f, 7.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("float",
+            ( clamp(a, 4.0f, 5.0f) ),
+            ( Vector2{5.0f, 4.0f} )
+        );
+    }
+
+    TEST_CASE("clamp_length") {
+        using draco::math::Vector2;
+        using draco::math::length;
+        using draco::math::clamp_length;
+
+        static constexpr Vector2 v{3.0f, 4.0f}; // len: 5
+
+        BASIC_R_SUBCASE("length < min",
+            ( clamp_length(v, 10.0f, 15.0f) ),
+            ( Vector2{6.0f, 8.0f} )
+        );
+
+        BASIC_R_SUBCASE("length == min",
+            ( clamp_length(v, 5.0f, 10.0f) ),
+            ( v )
+        );
+
+        BASIC_R_SUBCASE("length == max",
+            ( clamp_length(v, 3.0f, 5.0f) ),
+            ( v )
+        );
+
+        BASIC_R_SUBCASE("length > max",
+            ( clamp_length(v, 1.0f, 2.5f) ),
+            ( Vector2{1.5f, 2.0f} )
+        );
+    }
+
+    TEST_CASE("abs") {
+        using draco::math::Vector2;
+        using draco::math::abs;
+
+        RAC_CHECK_EQ(
+            ( abs(Vector2{1.0f, -2.0f}) ),
+            ( Vector2{1.0f, 2.0f} )
+        );
+    }
+
+    TEST_CASE("rounding") {
+        using draco::math::Vector2;
+        using draco::math::floor;
+        using draco::math::ceil;
+        using draco::math::trunc;
+        using draco::math::round;
+        
+        static constexpr Vector2 a{0.5f, 1.4f};
+        static constexpr Vector2 b{-1.0f, 1.0f};
+
+        BASIC_RAC_SUBCASE_2("floor",
+            ( floor(a) ),
+            ( Vector2{0.0f, 1.0f} ),
+            ( floor(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("ceil",
+            ( ceil(a) ),
+            ( Vector2{1.0f, 2.0f} ),
+            ( ceil(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("trunc",
+            ( trunc(a) ),
+            ( Vector2{0.0f, 1.0f} ),
+            ( trunc(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("round",
+            ( round(a) ),
+            ( Vector2{1.0f, 1.0f} ),
+            ( round(b) ),
+            ( b )
+        );
+    }
+
+    TEST_CASE("sign") {
+        using draco::math::Vector2;
+        using draco::math::sign;
+
+        RAC_CHECK_EQ(
+            ( sign(Vector2{1.0f, -1.0f}) ),
+            ( Vector2{1.0f, -1.0f} )
+        );
+    }
+
+    TEST_CASE("approx_eq") {
+        using draco::math::Vector2;
+        using draco::math::approx_eq;
+        using draco::math::CMP_EPSILON;
+
+        static constexpr Vector2 v{1.0f, 2.0f};
+        static constexpr Vector2 offset = Vector2::x_axis(CMP_EPSILON);
+
+        BASIC_R_SUBCASE("distance < epsilon",
+            ( approx_eq(v, v + offset * 0.5f) ),
+            ( true )
+        );
+
+        BASIC_R_SUBCASE("distance == epsilon",
+            ( approx_eq(v, v + offset) ),
+            ( true )
+        );
+
+        BASIC_R_SUBCASE("distance > epsilon",
+            ( approx_eq(v, v + offset * 2.0f) ),
+            ( false )
+        );
+    }
 }
 
-TEST_CASE("Vector4 dot self") {
-    using draco::math::Vector4;
-    using draco::math::dot;
+TEST_SUITE("vector3") {
+    TEST_CASE("constructors") {
+        using draco::math::Vector2;
+        using draco::math::Vector3;
+        using draco::math::Vector4;
+
+        static constexpr Vector2 a{1.0f, 2.0f};
+        static constexpr Vector4 b{3.0f, 4.0f, 5.0f, 6.0f};
+        
+        BASIC_RAC_SUBCASE("float",
+            ( Vector3(1.0f) ),
+            ( Vector3{1.0f, 1.0f, 1.0f} )
+        );
+        
+        BASIC_RAC_SUBCASE("vec2, float",
+            ( Vector3(a, 7.0f) ),
+            ( Vector3{1.0f, 2.0f, 7.0f} )
+        );
+        
+        BASIC_RAC_SUBCASE("float, vec2",
+            ( Vector3(7.0f, a) ),
+            ( Vector3{7.0f, 1.0f, 2.0f} )
+        );
+        
+        BASIC_RAC_SUBCASE("vec4",
+            ( Vector3(b) ),
+            ( Vector3{3.0f, 4.0f, 5.0f} )
+        );
+    }
+
+    TEST_CASE("access") {
+        using draco::math::Vector3;
+
+        static constexpr Vector3 v(1.0f, 2.0f, 3.0f);
+
+        RAC_CHECK_EQ(v[0], 1.0f);
+        RAC_CHECK_EQ(v[1], 2.0f);
+        RAC_CHECK_EQ(v[2], 3.0f);
+    }
+
+    TEST_CASE("swizzle") {
+        using draco::math::Vector2;
+        using draco::math::Vector3;
+        using draco::math::Vector4;
+        
+        static constexpr Vector3 v{1.0f, 2.0f, 3.0f};
+
+        BASIC_RAC_SUBCASE("vec2",
+            ( v[1, 0] ),
+            ( Vector2{2.0f, 1.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec3",
+            ( v[1, 2, 0] ),
+            ( Vector3{2.0f, 3.0f, 1.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec4",
+            ( v[0, 2, 1, 0] ),
+            ( Vector4{1.0f, 3.0f, 2.0f, 1.0f} )
+        );
+    }
+
+    TEST_CASE("swap") {
+        using draco::math::Vector3;
+
+        Vector3 a{1.f, 2.f, 3.f};
+        Vector3 b{3.f, 2.f, 1.f};
+
+        std::swap(a, b);
+
+        CHECK_EQ(a, Vector3{3.f, 2.f, 1.f});
+        CHECK_EQ(b, Vector3{1.f, 2.f, 3.f});
+    }
+
+    TEST_CASE("dot") {
+        using draco::math::Vector3;
+        using draco::math::dot;
+
+        static constexpr Vector3 a{1.0f, 2.0f, 3.0f};
+        static constexpr Vector3 b{4.0f, 5.0f, 6.0f};
+
+        BASIC_RAC_SUBCASE("basic",
+            ( dot(a, b) ),
+            ( 32.0f )
+        );
+
+        BASIC_RAC_SUBCASE("self",
+            ( dot(a, a) ),
+            ( 14.0f )
+        );
+
+        BASIC_RAC_SUBCASE("zero",
+            ( dot(a, Vector3()) ),
+            ( 0.0f )
+        );
+    }
+
+    TEST_CASE("length") {
+        using draco::math::Vector3;
+        using draco::math::length;
+        using draco::math::length_sq;
+
+        static constexpr Vector3 v{2.0f, 4.0f, 4.0f};
+
+        BASIC_R_SUBCASE("normal",
+            ( length(v) ),
+            ( 6.0f )
+        );
+
+        BASIC_RAC_SUBCASE("squared",
+            ( length_sq(v) ),
+            ( 36.0f )
+        );
+    }
+
+    TEST_CASE("distance") {
+        using draco::math::Vector3;
+        using draco::math::distance;
+        using draco::math::distance_sq;
+
+        static constexpr Vector3 a{2.0f, 4.0f, 4.0f};
+        static constexpr Vector3 b{-1.0f, -2.0f, -2.0f};
+
+        BASIC_R_SUBCASE("normal",
+            ( distance(a, b) ),
+            ( 9.0f )
+        );
+
+        BASIC_RAC_SUBCASE("squared",
+            ( distance_sq(a, b) ),
+            ( 81.0f )
+        );
+    }
+
+    TEST_CASE("normalize") {
+        using draco::math::Vector3;
+        using draco::math::length;
+        using draco::math::normalize;
+        using draco::math::normalize_fast;
+
+        static constexpr Vector3 a{0.0f, 6.4f, 4.8f};
+        static constexpr Vector3 b(1e-99);
+
+        const Vector3 result = normalize(a);
+        const Vector3 result_fast = normalize_fast(a);
+        const Vector3 result_zero = normalize(b);
+
+        CHECK_EQ(length(result), 1.0f);
+        CHECK_EQ(result, result_fast);
+        CHECK_EQ(result_zero, Vector3());
+    }
+
+    TEST_CASE("project") {
+        using draco::math::Vector3;
+        using draco::math::project;
+
+        static constexpr Vector3 a{2.0f, 8.0f, 4.0f};
+        static constexpr Vector3 b{1.0f, 1.0f, 2.0f};
+
+        RAC_CHECK_EQ(
+            ( project(a, b) ),
+            ( Vector3{3.0f, 3.0f, 6.0f} )
+        );
+    }
+
+    TEST_CASE("reflect") {
+        using draco::math::Vector3;
+        using draco::math::reflect;
+
+        static constexpr Vector3 a{1.0f, 2.0f, 3.0f};
+        static constexpr Vector3 b{4.0f, 5.0f, 6.0f};
+
+        RAC_CHECK_EQ(
+            ( reflect(a, b) ),
+            ( Vector3{-255.0f, -318.0f, -381.0f} )
+        );
+    }
+
+    TEST_CASE("angle") {
+        using draco::math::Vector3;
+        using draco::math::angle;
+        using draco::math::PI;
+
+        static constexpr Vector3 a{2.0f, 4.0f, 4.0f};
+        static constexpr Vector3 b{-4.0f, -8.0f, -8.0f};
+
+        R_CHECK_EQ(
+            ( angle(a, b) ),
+            ( PI )
+        );
+    }
+
+    TEST_CASE("lerp") {
+        using draco::math::Vector3;
+        using draco::math::lerp;
+
+        static constexpr Vector3 a{1.0f, 2.0f, 3.0f};
+        static constexpr Vector3 b{4.0f, 5.0f, 6.0f};
+
+        BASIC_RAC_SUBCASE("weight = -1",
+            ( lerp(a, b, -1.0f) ),
+            ( Vector3{-2.0f, -1.0f, -0.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 0",
+            ( lerp(a, b, 0.0f) ),
+            ( a )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 0.5",
+            ( lerp(a, b, 0.5f) ),
+            ( Vector3{2.5f, 3.5f, 4.5f} )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 1",
+            ( lerp(a, b, 1.0f) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 2",
+            ( lerp(a, b, 2.0f) ),
+            ( Vector3{7.0f, 8.0f, 9.0f} )
+        );
+    }
+
+    TEST_CASE("min") {
+        using draco::math::Vector3;
+        using draco::math::min;
+
+        static constexpr Vector3 a{5.0f, 8.0f, 3.0f};
+        static constexpr Vector3 b{1.0f, 6.0f, 7.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( min(a, b) ),
+            ( Vector3{1.0f, 6.0f, 3.0f} )
+        );
+
+        static constexpr Vector3 expected{4.0f, 4.0f, 3.0f};
 
-    Vector4 v{1.0f, 2.0f, 3.0f, 4.0f};
+        BASIC_RAC_SUBCASE_2("float",
+            ( min(a, 4.0f) ),
+            ( expected ),
+            ( min(4.0f, a) ),
+            ( expected )
+        );
+    }
 
-    const float result = dot(v, v);
-    constexpr float expected = 30.0f;
+    TEST_CASE("min_length") {
+        using draco::math::Vector3;
+        using draco::math::length;
+        using draco::math::min_length;
 
-    REQUIRE(result == expected);
+        static constexpr Vector3 a{2.0f, 4.0f, 4.0f};   // len: 6
+        static constexpr Vector3 b{5.0f, 10.0f, 10.0f}; // len: 15
+
+        BASIC_RAC_SUBCASE("vector",
+            ( min_length(a, b) ),
+            ( a )
+        );
+
+        SUBCASE("float") {
+            static constexpr float smaller_length = 1.0f;
+            static constexpr float larger_length = 10.0f;
+
+            const Vector3 result_smaller = min_length(a, smaller_length);
+            const Vector3 result_swapped = min_length(smaller_length, a);
+            const Vector3 result_larger = min_length(larger_length, a);
+
+            CHECK_EQ(length(result_smaller), smaller_length);
+            CHECK_EQ(result_smaller, result_swapped);
+            CHECK_EQ(result_larger, a);
+        }
+    }
+
+    TEST_CASE("max") {
+        using draco::math::Vector3;
+        using draco::math::max;
+
+        static constexpr Vector3 a{5.0f, 8.0f, 3.0f};
+        static constexpr Vector3 b{1.0f, 6.0f, 7.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( max(a, b) ),
+            ( Vector3{5.0f, 8.0f, 7.0f} )
+        );
+
+        static constexpr Vector3 expected{5.0f, 8.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE_2("float",
+            ( max(a, 4.0f) ),
+            ( expected ),
+            ( max(4.0f, a) ),
+            ( expected )
+        );
+    }
+
+    TEST_CASE("max_length") {
+        using draco::math::Vector3;
+        using draco::math::length;
+        using draco::math::max_length;
+
+        static constexpr Vector3 a{2.0f, 4.0f, 4.0f};   // len: 6
+        static constexpr Vector3 b{5.0f, 10.0f, 10.0f}; // len: 15
+
+        BASIC_RAC_SUBCASE("vector",
+            ( max_length(a, b) ),
+            ( b )
+        );
+
+        SUBCASE("float") {
+            static constexpr float smaller_length = 1.0f;
+            static constexpr float larger_length = 10.0f;
+
+            const Vector3 result_smaller = max_length(a, smaller_length);
+            const Vector3 result_swapped = max_length(larger_length, a);
+            const Vector3 result_larger = max_length(larger_length, a);
+
+            CHECK_EQ(length(result_larger), larger_length);
+            CHECK_EQ(result_larger, result_swapped);
+            CHECK_EQ(result_smaller, a);
+        }
+    }
+
+    TEST_CASE("clamp") {
+        using draco::math::Vector3;
+        using draco::math::clamp;
+
+        static constexpr Vector3 a{5.0f, 8.0f, 3.0f};
+        static constexpr Vector3 b{1.0f, 6.0f, 7.0f};
+        static constexpr Vector3 c{3.0f, 9.0f, 8.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( clamp(a, b, c) ),
+            ( Vector3{3.0f, 8.0f, 7.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("float",
+            ( clamp(a, 4.0f, 5.0f) ),
+            ( Vector3{5.0f, 5.0f, 4.0f} )
+        );
+    }
+
+    TEST_CASE("clamp_length") {
+        using draco::math::Vector3;
+        using draco::math::length;
+        using draco::math::clamp_length;
+
+        static constexpr Vector3 v{2.0f, 4.0f, 4.0f}; // len: 6
+
+        BASIC_R_SUBCASE("length < min",
+            ( clamp_length(v, 12.0f, 14.0f) ),
+            ( Vector3{4.0f, 8.0f, 8.0f} )
+        );
+
+        BASIC_R_SUBCASE("length == min",
+            ( clamp_length(v, 6.0f, 9.0f) ),
+            ( v )
+        );
+
+        BASIC_R_SUBCASE("length == max",
+            ( clamp_length(v, 3.0f, 6.0f) ),
+            ( v )
+        );
+
+        BASIC_R_SUBCASE("length > max",
+            ( clamp_length(v, 1.0f, 3.0f) ),
+            ( Vector3{1.0f, 2.0f, 2.0f} )
+        );
+    }
+
+    TEST_CASE("abs") {
+        using draco::math::Vector3;
+        using draco::math::abs;
+
+        RAC_CHECK_EQ(
+            ( abs(Vector3{1.0f, -2.0f, 0.0f}) ),
+            ( Vector3{1.0f, 2.0f, 0.0f} )
+        );
+    }
+
+    TEST_CASE("rounding") {
+        using draco::math::Vector3;
+        using draco::math::floor;
+        using draco::math::ceil;
+        using draco::math::trunc;
+        using draco::math::round;
+        
+        static constexpr Vector3 a{0.5f, -0.5f, 1.4f};
+        static constexpr Vector3 b{-1.0f, 0.0f, 1.0f};
+
+        BASIC_RAC_SUBCASE_2("floor",
+            ( floor(a) ),
+            ( Vector3{0.0f, -1.0f, 1.0f} ),
+            ( floor(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("ceil",
+            ( ceil(a) ),
+            ( Vector3{1.0f, 0.0f, 2.0f} ),
+            ( ceil(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("trunc",
+            ( trunc(a) ),
+            ( Vector3{0.0f, 0.0f, 1.0f} ),
+            ( trunc(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("round",
+            ( round(a) ),
+            ( Vector3{1.0f, -1.0f, 1.0f} ),
+            ( round(b) ),
+            ( b )
+        );
+    }
+
+    TEST_CASE("sign") {
+        using draco::math::Vector3;
+        using draco::math::sign;
+
+        RAC_CHECK_EQ(
+            ( sign(Vector3{1.0f, -1.0f, 0.0f}) ),
+            ( Vector3{1.0f, -1.0f, 0.0f} )
+        );
+    }
+
+    TEST_CASE("approx_eq") {
+        using draco::math::Vector3;
+        using draco::math::approx_eq;
+        using draco::math::CMP_EPSILON;
+
+        static constexpr Vector3 v{1.0f, 2.0f, 3.0f};
+        static constexpr Vector3 offset = Vector3::x_axis(CMP_EPSILON);
+
+        BASIC_R_SUBCASE("distance < epsilon",
+            ( approx_eq(v, v + offset * 0.5f) ),
+            ( true )
+        );
+
+        BASIC_R_SUBCASE("distance == epsilon",
+            ( approx_eq(v, v + offset) ),
+            ( true )
+        );
+
+        BASIC_R_SUBCASE("distance > epsilon",
+            ( approx_eq(v, v + offset * 2.0f) ),
+            ( false )
+        );
+    }
+
+    TEST_CASE("cross") {
+        using draco::math::Vector3;
+        using draco::math::cross;
+
+        RAC_CHECK_EQ(
+            ( cross(Vector3::x_axis(), Vector3::y_axis()) ),
+            ( Vector3::z_axis() )
+        )
+    }
 }
+
+TEST_SUITE("vector4") {
+    TEST_CASE("constructors") {
+        using draco::math::Vector2;
+        using draco::math::Vector3;
+        using draco::math::Vector4;
+
+        static constexpr Vector2 a{1.0f, 2.0f};
+        static constexpr Vector3 b{3.0f, 4.0f, 5.0f};
+        
+        BASIC_RAC_SUBCASE("float",
+            ( Vector4(1.0f) ),
+            ( Vector4{1.0f, 1.0f, 1.0f, 1.0f} )
+        );
+        
+        BASIC_RAC_SUBCASE("vec2, float, float",
+            ( Vector4(a, 6.0f, 7.0f) ),
+            ( Vector4{1.0f, 2.0f, 6.0f, 7.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("float, vec2, float",
+            ( Vector4(6.0f, a, 7.0f) ),
+            ( Vector4{6.0f, 1.0f, 2.0f, 7.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("float, float, vec2",
+            ( Vector4(6.0f, 7.0f, a) ),
+            ( Vector4{6.0f, 7.0f, 1.0f, 2.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec2, vec2",
+            ( Vector4(a, a) ),
+            ( Vector4{1.0f, 2.0f, 1.0f, 2.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec3, float",
+            ( Vector4(b, 6.0f) ),
+            ( Vector4{3.0f, 4.0f, 5.0f, 6.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("float, vec3",
+            ( Vector4(6.0f, b) ),
+            ( Vector4{6.0f, 3.0f, 4.0f, 5.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec2",
+            ( Vector4(a) ),
+            ( Vector4{1.0f, 2.0f, 0.0f, 0.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec3",
+            ( Vector4(b) ),
+            ( Vector4{3.0f, 4.0f, 5.0f, 0.0f} )
+        );
+    }
+
+    TEST_CASE("access") {
+        using draco::math::Vector4;
+
+        static constexpr Vector4 v(1.0f, 2.0f, 3.0f, 4.0f);
+
+        RAC_CHECK_EQ(v[0], 1.0f);
+        RAC_CHECK_EQ(v[1], 2.0f);
+        RAC_CHECK_EQ(v[2], 3.0f);
+        RAC_CHECK_EQ(v[3], 4.0f);
+    }
+
+    TEST_CASE("swizzle") {
+        using draco::math::Vector2;
+        using draco::math::Vector3;
+        using draco::math::Vector4;
+        
+        static constexpr Vector4 v{1.0f, 2.0f, 3.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE("vec2",
+            ( v[1, 0] ),
+            ( Vector2{2.0f, 1.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec3",
+            ( v[1, 2, 0] ),
+            ( Vector3{2.0f, 3.0f, 1.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("vec4",
+            ( v[0, 2, 1, 3] ),
+            ( Vector4{1.0f, 3.0f, 2.0f, 4.0f} )
+        );
+    }
+
+    TEST_CASE("swap") {
+        using draco::math::Vector4;
+
+        Vector4 a{1.f, 2.f, 3.f, 4.f};
+        Vector4 b{4.f, 3.f, 2.f, 1.f};
+
+        std::swap(a, b);
+
+        CHECK_EQ(a, Vector4{4.f, 3.f, 2.f, 1.f});
+        CHECK_EQ(b, Vector4{1.f, 2.f, 3.f, 4.f});
+    }
+
+    TEST_CASE("dot") {
+        using draco::math::Vector4;
+        using draco::math::dot;
+
+        static constexpr Vector4 a{1.0f, 2.0f, 3.0f, 4.0f};
+        static constexpr Vector4 b{5.0f, 6.0f, 7.0f, 8.0f};
+
+        BASIC_RAC_SUBCASE("basic",
+            ( dot(a, b) ),
+            ( 70.0f )
+        );
+
+        BASIC_RAC_SUBCASE("self",
+            ( dot(a, a) ),
+            ( 30.0f )
+        );
+
+        BASIC_RAC_SUBCASE("zero",
+            ( dot(a, Vector4()) ),
+            ( 0.0f )
+        );
+    }
+
+    TEST_CASE("length") {
+        using draco::math::Vector4;
+        using draco::math::length;
+        using draco::math::length_sq;
+
+        static constexpr Vector4 v{1.0f, 2.0f, 2.0f, 4.0f};
+
+        BASIC_R_SUBCASE("normal",
+            ( length(v) ),
+            ( 5.0f )
+        );
+
+        BASIC_RAC_SUBCASE("squared",
+            ( length_sq(v) ),
+            ( 25.0f )
+        );
+    }
+
+    TEST_CASE("distance") {
+        using draco::math::Vector4;
+        using draco::math::distance;
+        using draco::math::distance_sq;
+
+        static constexpr Vector4 a{1.0f, 2.0f, 2.0f, 4.0f};
+        static constexpr Vector4 b{3.0f, 6.0f, 7.0f, 10.0f};
+
+        BASIC_R_SUBCASE("normal",
+            ( distance(a, b) ),
+            ( 9.0f )
+        );
+
+        BASIC_RAC_SUBCASE("squared",
+            ( distance_sq(a, b) ),
+            ( 81.0f )
+        );
+    }
+
+    TEST_CASE("normalize") {
+        using draco::math::Vector4;
+        using draco::math::length;
+        using draco::math::normalize;
+        using draco::math::normalize_fast;
+
+        static constexpr Vector4 a{2.0f, 4.0f, 5.0f, 6.0f};
+        static constexpr Vector4 b(1e-99);
+
+        const Vector4 result = normalize(a);
+        const Vector4 result_fast = normalize_fast(a);
+        const Vector4 result_zero = normalize(b);
+
+        CHECK_EQ(length(result), 1.0f);
+        CHECK_EQ(result, result_fast);
+        CHECK_EQ(result_zero, Vector4());
+    }
+
+    TEST_CASE("project") {
+        using draco::math::Vector4;
+        using draco::math::project;
+
+        static constexpr Vector4 a{8.0f, 2.0f, 6.0f, 8.0f};
+        static constexpr Vector4 b{12.0f, 14.0f, 8.0f, 6.0f};
+
+        RAC_CHECK_EQ(
+            ( project(a, b) ),
+            ( Vector4{6.0f, 7.0f, 4.0f, 3.0f} )
+        );
+    }
+
+    TEST_CASE("reflect") {
+        using draco::math::Vector4;
+        using draco::math::reflect;
+
+        static constexpr Vector4 a{1.0f, 2.0f, 3.0f, 4.0f};
+        static constexpr Vector4 b{5.0f, 6.0f, 7.0f, 8.0f};
+
+        RAC_CHECK_EQ(
+            ( reflect(a, b) ),
+            ( Vector4{-699.0f, -838.0f, -977.0f, -1116.0f} )
+        );
+    }
+
+    TEST_CASE("angle") {
+        using draco::math::Vector4;
+        using draco::math::angle;
+        using draco::math::PI2;
+
+        static constexpr Vector4 a{1.0f, 5.0f, 1.0f, 3.0f};
+        static constexpr Vector4 b{2.0f, -6.0f, -2.0f, 10.0f};
+
+        R_CHECK_EQ(
+            ( angle(a, b) ),
+            ( PI2 )
+        );
+    }
+
+    TEST_CASE("lerp") {
+        using draco::math::Vector4;
+        using draco::math::lerp;
+
+        static constexpr Vector4 a{1.0f, 2.0f, 3.0f, 4.0f};
+        static constexpr Vector4 b{5.0f, 6.0f, 7.0f, 8.0f};
+
+        BASIC_RAC_SUBCASE("weight = -1",
+            ( lerp(a, b, -1.0f) ),
+            ( Vector4{-3.0f, -2.0f, -1.0f, 0.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 0",
+            ( lerp(a, b, 0.0f) ),
+            ( a )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 0.5",
+            ( lerp(a, b, 0.5f) ),
+            ( Vector4{3.0f, 4.0f, 5.0f, 6.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 1",
+            ( lerp(a, b, 1.0f) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE("weight = 2",
+            ( lerp(a, b, 2.0f) ),
+            ( Vector4{9.0f, 10.0f, 11.0f, 12.0f} )
+        );
+    }
+
+    TEST_CASE("min") {
+        using draco::math::Vector4;
+        using draco::math::min;
+
+        static constexpr Vector4 a{5.0f, 8.0f, 3.0f, 4.0f};
+        static constexpr Vector4 b{1.0f, 6.0f, 7.0f, 2.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( min(a, b) ),
+            ( Vector4{1.0f, 6.0f, 3.0f, 2.0f} )
+        );
+
+        static constexpr Vector4 expected{4.0f, 4.0f, 3.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE_2("float",
+            ( min(a, 4.0f) ),
+            ( expected ),
+            ( min(4.0f, a) ),
+            ( expected )
+        );
+    }
+
+    TEST_CASE("min_length") {
+        using draco::math::Vector4;
+        using draco::math::length;
+        using draco::math::min_length;
+
+        static constexpr Vector4 a{1.0f, 2.0f, 2.0f, 4.0f};   // len: 5
+        static constexpr Vector4 b{1.0f, -3.0f, -1.0f, 5.0f}; // len: 6
+
+        BASIC_RAC_SUBCASE("vector",
+            ( min_length(a, b) ),
+            ( a )
+        );
+
+        SUBCASE("float") {
+            static constexpr float smaller_length = 1.0f;
+            static constexpr float larger_length = 10.0f;
+
+            const Vector4 result_smaller = min_length(a, smaller_length);
+            const Vector4 result_swapped = min_length(smaller_length, a);
+            const Vector4 result_larger = min_length(larger_length, a);
+
+            CHECK_EQ(length(result_smaller), smaller_length);
+            CHECK_EQ(result_smaller, result_swapped);
+            CHECK_EQ(result_larger, a);
+        }
+    }
+
+    TEST_CASE("max") {
+        using draco::math::Vector4;
+        using draco::math::max;
+
+        static constexpr Vector4 a{5.0f, 8.0f, 3.0f, 4.0f};
+        static constexpr Vector4 b{1.0f, 6.0f, 7.0f, 2.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( max(a, b) ),
+            ( Vector4{5.0f, 8.0f, 7.0f, 4.0f} )
+        );
+
+        static constexpr Vector4 expected{5.0f, 8.0f, 4.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE_2("float",
+            ( max(a, 4.0f) ),
+            ( expected ),
+            ( max(4.0f, a) ),
+            ( expected )
+        );
+    }
+
+    TEST_CASE("max_length") {
+        using draco::math::Vector4;
+        using draco::math::length;
+        using draco::math::max_length;
+
+        static constexpr Vector4 a{1.0f, 2.0f, 2.0f, 4.0f};   // len: 5
+        static constexpr Vector4 b{1.0f, -3.0f, -1.0f, 5.0f}; // len: 6
+
+        BASIC_RAC_SUBCASE("vector",
+            ( max_length(a, b) ),
+            ( b )
+        );
+
+        SUBCASE("float") {
+            static constexpr float smaller_length = 1.0f;
+            static constexpr float larger_length = 10.0f;
+
+            const Vector4 result_smaller = max_length(a, smaller_length);
+            const Vector4 result_swapped = max_length(larger_length, a);
+            const Vector4 result_larger = max_length(larger_length, a);
+
+            CHECK_EQ(length(result_larger), larger_length);
+            CHECK_EQ(result_larger, result_swapped);
+            CHECK_EQ(result_smaller, a);
+        }
+    }
+
+    TEST_CASE("clamp") {
+        using draco::math::Vector4;
+        using draco::math::clamp;
+
+        static constexpr Vector4 a{5.0f, 8.0f, 3.0f, 4.0f};
+        static constexpr Vector4 b{1.0f, 6.0f, 7.0f, 2.0f};
+        static constexpr Vector4 c{3.0f, 9.0f, 8.0f, 4.0f};
+
+        BASIC_RAC_SUBCASE("vector",
+            ( clamp(a, b, c) ),
+            ( Vector4{3.0f, 8.0f, 7.0f, 4.0f} )
+        );
+
+        BASIC_RAC_SUBCASE("float",
+            ( clamp(a, 4.0f, 5.0f) ),
+            ( Vector4{5.0f, 5.0f, 4.0f, 4.0f} )
+        );
+    }
+
+    TEST_CASE("clamp_length") {
+        using draco::math::Vector4;
+        using draco::math::length;
+        using draco::math::clamp_length;
+
+        static constexpr Vector4 v{1.0f, -3.0f, -1.0f, 5.0f}; // len: 6
+
+        BASIC_R_SUBCASE("length < min",
+            ( clamp_length(v, 12.0f, 14.0f) ),
+            ( Vector4{2.0f, -6.0f, -2.0f, 10.0f} )
+        );
+
+        BASIC_R_SUBCASE("length == min",
+            ( clamp_length(v, 6.0f, 9.0f) ),
+            ( v )
+        );
+
+        BASIC_R_SUBCASE("length == max",
+            ( clamp_length(v, 3.0f, 6.0f) ),
+            ( v )
+        );
+
+        BASIC_R_SUBCASE("length > max",
+            ( clamp_length(v, 1.0f, 3.0f) ),
+            ( Vector4{0.5f, -1.5f, -0.5f, 2.5f} )
+        );
+    }
+
+    TEST_CASE("abs") {
+        using draco::math::Vector4;
+        using draco::math::abs;
+
+        RAC_CHECK_EQ(
+            ( abs(Vector4{1.0f, -2.0f, -3.0f, 0.0f}) ),
+            ( Vector4{1.0f, 2.0f, 3.0f, 0.0f} )
+        );
+    }
+
+    TEST_CASE("rounding") {
+        using draco::math::Vector4;
+        using draco::math::floor;
+        using draco::math::ceil;
+        using draco::math::trunc;
+        using draco::math::round;
+        
+        static constexpr Vector4 a{0.5f, -0.5f, 1.4f, 1.6f};
+        static constexpr Vector4 b{-1.0f, 0.0f, 1.0f, 2.0f};
+
+        BASIC_RAC_SUBCASE_2("floor",
+            ( floor(a) ),
+            ( Vector4{0.0f, -1.0f, 1.0f, 1.0f} ),
+            ( floor(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("ceil",
+            ( ceil(a) ),
+            ( Vector4{1.0f, 0.0f, 2.0f, 2.0f} ),
+            ( ceil(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("trunc",
+            ( trunc(a) ),
+            ( Vector4{0.0f, 0.0f, 1.0f, 1.0f} ),
+            ( trunc(b) ),
+            ( b )
+        );
+
+        BASIC_RAC_SUBCASE_2("round",
+            ( round(a) ),
+            ( Vector4{1.0f, -1.0f, 1.0f, 2.0f} ),
+            ( round(b) ),
+            ( b )
+        );
+    }
+
+    TEST_CASE("sign") {
+        using draco::math::Vector4;
+        using draco::math::sign;
+
+        RAC_CHECK_EQ(
+            ( sign(Vector4{1.0f, -1.0f, 0.0f, -0.0f}) ),
+            ( Vector4{1.0f, -1.0f, 0.0f, 0.0f} )
+        );
+    }
+
+    TEST_CASE("approx_eq") {
+        using draco::math::Vector4;
+        using draco::math::approx_eq;
+        using draco::math::CMP_EPSILON;
+
+        static constexpr Vector4 v{1.0f, 2.0f, 3.0f, 4.0f};
+        static constexpr Vector4 offset = Vector4::x_axis(CMP_EPSILON);
+
+        BASIC_R_SUBCASE("distance < epsilon",
+            ( approx_eq(v, v + offset * 0.5f) ),
+            ( true )
+        );
+
+        BASIC_R_SUBCASE("distance == epsilon",
+            ( approx_eq(v, v + offset) ),
+            ( true )
+        );
+
+        BASIC_R_SUBCASE("distance > epsilon",
+            ( approx_eq(v, v + offset * 2.0f) ),
+            ( false )
+        );
+    }
+}
\ No newline at end of file
diff --git a/engine/native/core/math/types.cppm b/engine/native/core/math/types.cppm
new file mode 100644
index 0000000..0dbef11
--- /dev/null
+++ b/engine/native/core/math/types.cppm
@@ -0,0 +1,6 @@
+export module core.math.types;
+
+export import :common;
+export import :vector2;
+export import :vector3;
+export import :vector4;
\ No newline at end of file
diff --git a/engine/native/core/math/types_common.cppm b/engine/native/core/math/types_common.cppm
new file mode 100644
index 0000000..fe3d77a
--- /dev/null
+++ b/engine/native/core/math/types_common.cppm
@@ -0,0 +1,170 @@
+export module core.math.types:common;
+import core.defs;
+
+export namespace draco::math {
+	struct Vector2;
+	struct Vector3;
+	struct Vector4;
+
+    struct alignas(8) Vector2 {
+        float x, y;
+
+        // constructors
+        [[nodiscard]] constexpr Vector2() noexcept = default;
+        [[nodiscard]] constexpr explicit Vector2(float n) noexcept;
+        [[nodiscard]] constexpr Vector2(float x, float y) noexcept;
+        [[nodiscard]] constexpr explicit Vector2(const Vector3& xy) noexcept;
+        [[nodiscard]] constexpr explicit Vector2(const Vector4& xy) noexcept;
+        
+        // static
+        [[nodiscard]] static constexpr Vector2 x_axis(float x = 1.0f) noexcept;
+        [[nodiscard]] static constexpr Vector2 y_axis(float y = 1.0f) noexcept;
+        [[nodiscard]] static Vector2 polar(float angle, float radius = 1.0f) noexcept;
+
+        // element access
+        [[nodiscard]] constexpr float& operator[](int i) noexcept;
+        [[nodiscard]] constexpr const float& operator[](int i) const noexcept;
+
+        // swizzle
+        [[nodiscard]] constexpr Vector2 operator[](int i0, int i1) noexcept;
+        [[nodiscard]] constexpr Vector2 operator[](int i0, int i1) const noexcept;
+        [[nodiscard]] constexpr Vector3 operator[](int i0, int i1, int i2) noexcept;
+        [[nodiscard]] constexpr Vector3 operator[](int i0, int i1, int i2) const noexcept;
+        [[nodiscard]] constexpr Vector4 operator[](int i0, int i1, int i2, int i3) noexcept;
+        [[nodiscard]] constexpr Vector4 operator[](int i0, int i1, int i2, int i3) const noexcept;
+
+        // operators
+        [[nodiscard]] constexpr Vector2 operator+() const noexcept;
+        [[nodiscard]] constexpr Vector2 operator-() const noexcept;
+        [[nodiscard]] constexpr bool operator==(const Vector2& other) const noexcept = default;
+        constexpr Vector2& operator+=(const Vector2& other) noexcept;
+        constexpr Vector2& operator+=(float other) noexcept;
+        constexpr Vector2& operator-=(const Vector2& other) noexcept;
+        constexpr Vector2& operator-=(float other) noexcept;
+        constexpr Vector2& operator*=(const Vector2& other) noexcept;
+        constexpr Vector2& operator*=(float other) noexcept;
+        constexpr Vector2& operator/=(const Vector2& other) noexcept;
+        constexpr Vector2& operator/=(float other) noexcept;
+        constexpr Vector2& operator=(float other) noexcept;
+    };
+
+    struct alignas(16) Vector3 {
+        float x, y, z;
+
+        // constructors
+        [[nodiscard]] constexpr Vector3() noexcept = default;
+        [[nodiscard]] constexpr explicit Vector3(float n) noexcept;
+        [[nodiscard]] constexpr Vector3(float x, float y, float z) noexcept;
+        [[nodiscard]] constexpr explicit Vector3(const Vector2& xy, float z = 0.0f) noexcept;
+        [[nodiscard]] constexpr Vector3(float x, const Vector2& yz) noexcept;
+        [[nodiscard]] constexpr explicit Vector3(const Vector4& xyz) noexcept;
+        
+        // static
+        [[nodiscard]] static constexpr Vector3 x_axis(float x = 1.0f) noexcept;
+        [[nodiscard]] static constexpr Vector3 y_axis(float y = 1.0f) noexcept;
+        [[nodiscard]] static constexpr Vector3 z_axis(float z = 1.0f) noexcept;
+        [[nodiscard]] static Vector3 spherical(float azimuth, float inclination, float radius = 1.0f) noexcept;
+        [[nodiscard]] static Vector3 cylindrical(float angle, float radius = 1.0f, float height = 0.0f) noexcept;
+
+        // element access
+        [[nodiscard]] constexpr float& operator[](int i) noexcept;
+        [[nodiscard]] constexpr const float& operator[](int i) const noexcept;
+
+        // swizzle
+        [[nodiscard]] constexpr Vector2 operator[](int i0, int i1) noexcept;
+        [[nodiscard]] constexpr Vector2 operator[](int i0, int i1) const noexcept;
+        [[nodiscard]] constexpr Vector3 operator[](int i0, int i1, int i2) noexcept;
+        [[nodiscard]] constexpr Vector3 operator[](int i0, int i1, int i2) const noexcept;
+        [[nodiscard]] constexpr Vector4 operator[](int i0, int i1, int i2, int i3) noexcept;
+        [[nodiscard]] constexpr Vector4 operator[](int i0, int i1, int i2, int i3) const noexcept;
+
+        // operators
+        [[nodiscard]] constexpr Vector3 operator+() const noexcept;
+        [[nodiscard]] constexpr Vector3 operator-() const noexcept;
+        [[nodiscard]] constexpr bool operator==(const Vector3& other) const noexcept = default;
+        constexpr Vector3& operator+=(const Vector3& other) noexcept;
+        constexpr Vector3& operator+=(float other) noexcept;
+        constexpr Vector3& operator-=(const Vector3& other) noexcept;
+        constexpr Vector3& operator-=(float other) noexcept;
+        constexpr Vector3& operator*=(const Vector3& other) noexcept;
+        constexpr Vector3& operator*=(float other) noexcept;
+        constexpr Vector3& operator/=(const Vector3& other) noexcept;
+        constexpr Vector3& operator/=(float other) noexcept;
+        constexpr Vector3& operator=(float other) noexcept;
+    };
+
+	struct alignas(16) Vector4 {
+        float x, y, z, w;
+
+        // constructors
+        [[nodiscard]] constexpr Vector4() noexcept = default;
+        [[nodiscard]] constexpr explicit Vector4(float n) noexcept;
+        [[nodiscard]] constexpr Vector4(float x, float y, float z, float w) noexcept;
+        [[nodiscard]] constexpr explicit Vector4(const Vector2& xy) noexcept;
+        [[nodiscard]] constexpr Vector4(const Vector2& xy, float z, float w) noexcept;
+        [[nodiscard]] constexpr Vector4(float x, const Vector2& yz, float w) noexcept;
+        [[nodiscard]] constexpr Vector4(float x, float y, const Vector2& zw) noexcept;
+        [[nodiscard]] constexpr Vector4(const Vector2& xy, const Vector2& zw) noexcept;
+        [[nodiscard]] constexpr explicit Vector4(const Vector3& xyz, float w = 0.0f) noexcept;
+        [[nodiscard]] constexpr Vector4(float x, const Vector3& yzw) noexcept;
+        
+        // static
+        [[nodiscard]] static constexpr Vector4 x_axis(float x = 1.0f) noexcept;
+        [[nodiscard]] static constexpr Vector4 y_axis(float y = 1.0f) noexcept;
+        [[nodiscard]] static constexpr Vector4 z_axis(float z = 1.0f) noexcept;
+        [[nodiscard]] static constexpr Vector4 w_axis(float w = 1.0f) noexcept;
+
+        // element access
+        [[nodiscard]] constexpr float& operator[](int i) noexcept;
+        [[nodiscard]] constexpr const float& operator[](int i) const noexcept;
+
+        // swizzle
+        [[nodiscard]] constexpr Vector2 operator[](int i0, int i1) noexcept;
+        [[nodiscard]] constexpr Vector2 operator[](int i0, int i1) const noexcept;
+        [[nodiscard]] constexpr Vector3 operator[](int i0, int i1, int i2) noexcept;
+        [[nodiscard]] constexpr Vector3 operator[](int i0, int i1, int i2) const noexcept;
+        [[nodiscard]] constexpr Vector4 operator[](int i0, int i1, int i2, int i3) noexcept;
+        [[nodiscard]] constexpr Vector4 operator[](int i0, int i1, int i2, int i3) const noexcept;
+
+        // member operators
+        [[nodiscard]] constexpr Vector4 operator+() const noexcept;
+        [[nodiscard]] constexpr Vector4 operator-() const noexcept;
+        [[nodiscard]] constexpr bool operator==(const Vector4& other) const noexcept = default;
+        constexpr Vector4& operator+=(const Vector4& other) noexcept;
+        constexpr Vector4& operator+=(float other) noexcept;
+        constexpr Vector4& operator-=(const Vector4& other) noexcept;
+        constexpr Vector4& operator-=(float other) noexcept;
+        constexpr Vector4& operator*=(const Vector4& other) noexcept;
+        constexpr Vector4& operator*=(float other) noexcept;
+        constexpr Vector4& operator/=(const Vector4& other) noexcept;
+        constexpr Vector4& operator/=(float other) noexcept;
+        constexpr Vector4& operator=(float other) noexcept;
+    };
+}
+
+template<typename T> consteval T select(const int i, const T v1, const T v2) {
+    switch (i) {
+        case 0:  return v1;
+        case 1:  return v2;
+        default: throw "Index out of range";
+    }
+}
+
+template<typename T> consteval T select(const int i, const T v1, const T v2, const T v3) {
+    switch (i) {
+        case 0:  return v1;
+        case 1:  return v2;
+        case 2:  return v3;
+        default: throw "Index out of range";
+    }
+}
+
+template<typename T> consteval T select(const int i, const T v1, const T v2, const T v3, const T v4) {
+    switch (i) {
+        case 0:  return v1;
+        case 1:  return v2;
+        case 2:  return v3;
+        case 3:  return v4;
+        default: throw "Index out of range";
+    }
+}
\ No newline at end of file
diff --git a/engine/native/core/math/vector2.cppm b/engine/native/core/math/vector2.cppm
new file mode 100644
index 0000000..5446043
--- /dev/null
+++ b/engine/native/core/math/vector2.cppm
@@ -0,0 +1,463 @@
+module;
+
+#include <cmath>
+#include <algorithm>
+#include <format>
+#include "platform/simd.h"
+
+#if ARCH_X64
+    #include <immintrin.h>
+#elif ARCH_ARM64
+    #include <arm_neon.h>
+#endif
+
+export module core.math.types:vector2;
+export import :common;
+
+import core.math.constants;
+import core.math.functions;
+import core.defs;
+
+export namespace draco::math {
+    // assertions
+    static_assert(sizeof(Vector2) == 8, "Vector2 must be 8 bytes");
+    static_assert(alignof(Vector2) == 8, "Vector2 must be 8-byte aligned");
+    static_assert(trivial<Vector2>, "Vector2 must be trivial");
+    static_assert(std::is_standard_layout_v<Vector2>, "Vector2 must be standard layout");
+
+    // constructors
+    [[nodiscard]] constexpr Vector2::Vector2(const float n) noexcept
+        : x{n}, y{n} { }
+    
+    [[nodiscard]] constexpr Vector2::Vector2(const float x, const float y) noexcept
+        : x{x}, y{y} { }
+        
+    [[nodiscard]] constexpr Vector2::Vector2(const Vector3& xy) noexcept
+        : x{xy.x}, y{xy.y} { }
+
+    [[nodiscard]] constexpr Vector2::Vector2(const Vector4& xy) noexcept
+        : x{xy.x}, y{xy.y} { }
+    
+    // static
+    [[nodiscard]] constexpr Vector2 Vector2::x_axis(const float x) noexcept {
+        return { x, 0.0f };
+    }
+
+    [[nodiscard]] constexpr Vector2 Vector2::y_axis(const float y) noexcept {
+        return { 0.0f, y };
+    }
+
+    [[nodiscard]] Vector2 Vector2::polar(const float angle, const float radius) noexcept {
+        return { radius * std::cos(angle), radius * std::sin(angle) };
+    }
+    
+    // element access
+    [[nodiscard]] constexpr float& Vector2::operator[](const int i) noexcept {
+        if consteval {
+            return i ? y : x;
+        } else {
+            return (&x)[i];
+        }
+    }
+
+    [[nodiscard]] constexpr const float& Vector2::operator[](const int i) const noexcept {
+        if consteval {
+            return i ? y : x;
+        } else {
+            return (&x)[i];
+        }
+    }
+
+    // swizzle
+    [[nodiscard]] constexpr Vector2 Vector2::operator[](const int i0, const int i1) noexcept {
+        if consteval {
+            return { select(i0, x, y), select(i1, x, y) };
+        } else {
+            return { (&x)[i0], (&x)[i1] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector2 Vector2::operator[](const int i0, const int i1) const noexcept {
+        if consteval {
+            return { select(i0, x, y), select(i1, x, y) };
+        } else {
+            return { (&x)[i0], (&x)[i1] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector2::operator[](const int i0, const int i1, const int i2) noexcept {
+        if consteval {
+            return { select(i0, x, y), select(i1, x, y), select(i2, x, y) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector2::operator[](const int i0, const int i1, const int i2) const noexcept {
+        if consteval {
+            return { select(i0, x, y), select(i1, x, y), select(i2, x, y) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector2::operator[](const int i0, const int i1, const int i2, const int i3) noexcept {
+        if consteval {
+            return { select(i0, x, y), select(i1, x, y), select(i2, x, y), select(i3, x, y) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2], (&x)[i3] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector2::operator[](const int i0, const int i1, const int i2, const int i3) const noexcept {
+        if consteval {
+            return { select(i0, x, y), select(i1, x, y), select(i2, x, y), select(i3, x, y) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2], (&x)[i3] };
+        }
+    }
+
+    // operators
+    constexpr Vector2& Vector2::operator+=(const Vector2& other) noexcept {
+        x += other.x;
+        y += other.y;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator+=(const float other) noexcept {
+        x += other;
+        y += other;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator-=(const Vector2& other) noexcept {
+        x -= other.x;
+        y -= other.y;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator-=(const float other) noexcept {
+        x -= other;
+        y -= other;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator*=(const Vector2& other) noexcept {
+        x *= other.x;
+        y *= other.y;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator*=(const float other) noexcept {
+        x *= other;
+        y *= other;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator/=(const Vector2& other) noexcept {
+        x /= other.x;
+        y /= other.y;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator/=(const float other) noexcept {
+        const float inv = 1.0f / other;
+        x *= inv;
+        y *= inv;
+        return *this;
+    }
+
+    constexpr Vector2& Vector2::operator=(const float other) noexcept {
+        x = other;
+        y = other;
+        return *this;
+    }
+
+    [[nodiscard]] constexpr Vector2 Vector2::operator+() const noexcept {
+        return { x, y };
+    }
+
+    [[nodiscard]] constexpr Vector2 Vector2::operator-() const noexcept {
+        return { -x, -y };
+    }
+    
+    [[nodiscard]] constexpr Vector2 operator+(const Vector2& a, const Vector2& b) noexcept {
+        return { a.x+b.x, a.y+b.y };
+    }
+
+    [[nodiscard]] constexpr Vector2 operator+(const Vector2& a, const float b) noexcept {
+        return { a.x+b, a.y+b };
+    }
+
+    [[nodiscard]] constexpr Vector2 operator+(const float a, const Vector2& b) noexcept {
+        return b+a;
+    }
+
+    [[nodiscard]] constexpr Vector2 operator-(const Vector2& a, const Vector2& b) noexcept {
+        return { a.x-b.x, a.y-b.y };
+    }
+
+    [[nodiscard]] constexpr Vector2 operator-(const Vector2& a, const float b) noexcept {
+        return { a.x-b, a.y-b };
+    }
+
+    [[nodiscard]] constexpr Vector2 operator-(const float a, const Vector2& b) noexcept {
+        return { a-b.x, a-b.y };
+    }
+
+    [[nodiscard]] constexpr Vector2 operator*(const Vector2& a, const Vector2& b) noexcept {
+        return { a.x*b.x, a.y*b.y };
+    }
+
+    [[nodiscard]] constexpr Vector2 operator*(const Vector2& a, const float b) noexcept {
+        return { a.x*b, a.y*b };
+    }
+
+    [[nodiscard]] constexpr Vector2 operator*(const float a, const Vector2& b) noexcept {
+        return b*a;
+    }
+
+    [[nodiscard]] constexpr Vector2 operator/(const Vector2& a, const Vector2& b) noexcept {
+        return { a.x/b.x, a.y/b.y };
+    }
+    
+    [[nodiscard]] constexpr Vector2 operator/(const Vector2& a, const float b) noexcept {
+        return a * (1.0f / b);
+    }
+    
+    [[nodiscard]] constexpr Vector2 operator/(const float a, const Vector2& b) noexcept {
+        return { a/b.x, a/b.y };
+    }
+
+    // functions
+
+    // Returns dot product
+    [[nodiscard]] constexpr float dot(const Vector2& a, const Vector2& b) noexcept {
+        return a.x*b.x + a.y*b.y;
+    }
+
+    // Returns squared magnitude
+    [[nodiscard]] constexpr float length_sq(const Vector2& v) noexcept {
+        return dot(v, v);
+    }
+
+    // Returns magnitude
+    [[nodiscard]] float length(const Vector2& v) noexcept {
+        return std::sqrt(length_sq(v));
+    }
+
+    // Return squared distance between two vectors
+    [[nodiscard]] constexpr float distance_sq(const Vector2& a, const Vector2& b) noexcept {
+        return length_sq(a - b);
+    }
+
+    // Returns distance between two vectors
+    [[nodiscard]] float distance(const Vector2& a, const Vector2& b) noexcept {
+        return length(a - b);
+    }
+
+    // Safe normalize, checks length
+    [[nodiscard]] Vector2 normalize(const Vector2& v) noexcept {
+        const float len = length(v);
+
+        return (len > CMP_NORMALIZE_TOLERANCE) ? v / len : Vector2();
+    }
+    
+    // Faster normalize, it presupposes vector has non-zero length
+    // TODO: add check that v is non-zero on debug builds
+    [[nodiscard]] Vector2 normalize_fast(const Vector2& v) noexcept {
+        return v / length(v);
+    }
+
+    // Returns vector projected onto normal
+    [[nodiscard]] constexpr Vector2 project(const Vector2& vector, const Vector2& normal) noexcept {
+        return normal * (dot(vector, normal) / length_sq(normal));
+    }
+
+    // Returns a vector reflected off a plane defined by its normal
+    [[nodiscard]] constexpr Vector2 reflect(const Vector2& incoming, const Vector2& normal) noexcept {
+        return incoming - 2.0f * dot(incoming, normal) * normal;
+    }
+
+    // Returns the angle between two vectors
+    [[nodiscard]] float angle(const Vector2& a, const Vector2& b) noexcept {
+        return std::acos(dot(a, b) / (length(a) * length(b)));
+    }
+
+    // Returns linear interpolation between two vectors
+    [[nodiscard]] constexpr Vector2 lerp(const Vector2& from, const Vector2& to, const float weight) noexcept {
+        return {
+            lerp(from.x, to.x, weight),
+            lerp(from.y, to.y, weight)
+        };
+    }
+
+    // Returns component-wise minimum
+    [[nodiscard]] constexpr Vector2 min(const Vector2& a, const Vector2& b) noexcept {
+        return {
+            std::min(a.x, b.x),
+            std::min(a.y, b.y)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector2 min(const Vector2& a, const float b) noexcept {
+        return {
+            std::min(a.x, b),
+            std::min(a.y, b)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector2 min(const float a, const Vector2& b) noexcept {
+        return min(b, a);
+    }
+
+    // Returns the vector with the smaller length
+    [[nodiscard]] constexpr Vector2 min_length(const Vector2& a, const Vector2& b) noexcept {
+        return length_sq(a) < length_sq(b) ? a : b;
+    }
+
+    // Returns a vector in the same direction whose length is bounded above by the given value.
+    [[nodiscard]] Vector2 min_length(const Vector2& a, const float b) noexcept {
+        const float len_sq = length_sq(a);
+        
+        if (len_sq > b * b) {
+            return a * (b / std::sqrt(len_sq));
+        } else {
+            return a;
+        }
+    }
+
+    [[nodiscard]] Vector2 min_length(const float a, const Vector2& b) noexcept {
+        return min_length(b, a);
+    }
+
+    // Returns component-wise maximum
+    [[nodiscard]] constexpr Vector2 max(const Vector2& a, const Vector2& b) noexcept {
+        return {
+            std::max(a.x, b.x),
+            std::max(a.y, b.y)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector2 max(const Vector2& a, const float b) noexcept {
+        return {
+            std::max(a.x, b),
+            std::max(a.y, b)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector2 max(const float a, const Vector2& b) noexcept {
+        return max(b, a);
+    }
+
+    // Returns the vector with the larger length
+    [[nodiscard]] constexpr Vector2 max_length(const Vector2& a, const Vector2& b) noexcept {
+        return length_sq(a) > length_sq(b) ? a : b;
+    }
+
+    // Returns a vector in the same direction whose length is bounded below by the given value. Returns the 0 vector if the vector is too small to be normalized.
+    [[nodiscard]] Vector2 max_length(const Vector2& a, const float b) noexcept {
+        const float len_sq = length_sq(a);
+        
+        if (len_sq <= CMP_NORMALIZE_TOLERANCE2) {
+            return Vector2();
+        } else if (len_sq < b * b) {
+            return a * (b / std::sqrt(len_sq));
+        } else {
+            return a;
+        }
+    }
+
+    [[nodiscard]] Vector2 max_length(const float a, const Vector2& b) noexcept {
+        return max_length(b, a);
+    }
+
+    // Clamps each component of x to the range [x_min, x_max]. Presupposes x_min <= x_max.
+    [[nodiscard]] constexpr Vector2 clamp(const Vector2& x, const Vector2& x_min, const Vector2& x_max) noexcept {
+        return max(x_min, min(x, x_max));
+    }
+
+    [[nodiscard]] constexpr Vector2 clamp(const Vector2& x, const float x_min, const float x_max) noexcept {
+        return max(x_min, min(x, x_max));
+    }
+
+    // Clamps the length of the vector to the range [x_min, x_max]. Presupposes x_min <= x_max. Returns the 0 vector if the vector is too small to be normalized.
+    [[nodiscard]] Vector2 clamp_length(const Vector2& v, const float x_min, const float x_max) noexcept {
+        const float len_sq = length_sq(v);
+        
+        if (len_sq <= CMP_NORMALIZE_TOLERANCE2) {
+            return Vector2();
+        } else if (len_sq < x_min * x_min) {
+            return v * (x_min / std::sqrt(len_sq));
+        } else if (len_sq > x_max * x_max) {
+            return v * (x_max / std::sqrt(len_sq));
+        } else {
+            return v;
+        }
+    }
+
+    // Returns component-wise absolute value
+    [[nodiscard]] constexpr Vector2 abs(const Vector2& v) noexcept {
+        return {
+            abs(v.x),
+            abs(v.y)
+        };
+    }
+
+    // Returns component-wise floor
+    [[nodiscard]] constexpr Vector2 floor(const Vector2& v) noexcept {
+        return {
+            floor(v.x),
+            floor(v.y)
+        };
+    }
+
+    // Returns component-wise ceiling
+    [[nodiscard]] constexpr Vector2 ceil(const Vector2& v) noexcept {
+        return {
+            ceil(v.x),
+            ceil(v.y)
+        };
+    }
+
+    // Returns component-wise truncation
+    [[nodiscard]] constexpr Vector2 trunc(const Vector2& v) noexcept {
+        return {
+            trunc(v.x),
+            trunc(v.y)
+        };
+    }
+
+    // Returns component-wise round
+    [[nodiscard]] constexpr Vector2 round(const Vector2& v) noexcept {
+        return {
+            round(v.x),
+            round(v.y)
+        };
+    }
+
+    // Returns component-wise sign. Note that -0 still returns 0
+    [[nodiscard]] constexpr Vector2 sign(const Vector2& v) noexcept {
+        return {
+            sign(v.x),
+            sign(v.y)
+        };
+    }
+
+    // Returns true if the vectors are approximately equal
+    [[nodiscard]] constexpr bool approx_eq(const Vector2& a, const Vector2& b) noexcept {
+        return distance_sq(a, b) < CMP_EPSILON2;
+    }
+}
+
+export namespace std {
+    template<> struct formatter<draco::math::Vector2> : formatter<float> {
+        auto format(const draco::math::Vector2& v, format_context& ctx) const {
+            ctx.advance_to(format_to(ctx.out(), "{{"));
+            ctx.advance_to(formatter<float>::format(v.x, ctx));
+            ctx.advance_to(format_to(ctx.out(), ", "));
+            ctx.advance_to(formatter<float>::format(v.y, ctx));
+            return format_to(ctx.out(), "}}");
+        }
+    };
+}
\ No newline at end of file
diff --git a/engine/native/core/math/vector3.cppm b/engine/native/core/math/vector3.cppm
new file mode 100644
index 0000000..5254e90
--- /dev/null
+++ b/engine/native/core/math/vector3.cppm
@@ -0,0 +1,509 @@
+module;
+
+#include <cmath>
+#include <algorithm>
+#include <format>
+#include "platform/simd.h"
+
+#if ARCH_X64
+    #include <immintrin.h>
+#elif ARCH_ARM64
+    #include <arm_neon.h>
+#endif
+
+export module core.math.types:vector3;
+export import :common;
+
+import core.math.constants;
+import core.math.functions;
+import core.defs;
+
+export namespace draco::math {
+    // assertions
+    static_assert(sizeof(Vector3) == 16, "Vector3 must be 16 bytes");
+    static_assert(alignof(Vector3) == 16, "Vector3 must be 16-byte aligned");
+    static_assert(trivial<Vector3>, "Vector3 must be trivial");
+    static_assert(std::is_standard_layout_v<Vector3>, "Vector3 must be standard layout");
+
+    // constructors
+    [[nodiscard]] constexpr Vector3::Vector3(const float n) noexcept
+        : x{n}, y{n}, z{n} { }
+
+    [[nodiscard]] constexpr Vector3::Vector3(const float x, const float y, const float z) noexcept
+        : x{x}, y{y}, z{z} { }
+
+    [[nodiscard]] constexpr Vector3::Vector3(const Vector2& xy, const float z) noexcept
+        : x{xy.x}, y{xy.y}, z{z} { }
+
+    [[nodiscard]] constexpr Vector3::Vector3(const float x, const Vector2& yz) noexcept
+        : x{x}, y{yz.x}, z{yz.y} { }
+
+    [[nodiscard]] constexpr Vector3::Vector3(const Vector4& xyz) noexcept
+        : x{xyz.x}, y{xyz.y}, z{xyz.z} { }
+    
+    // static
+    [[nodiscard]] constexpr Vector3 Vector3::x_axis(const float x) noexcept {
+        return { x, 0.0f, 0.0f };
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector3::y_axis(const float y) noexcept {
+        return { 0.0f, y, 0.0f };
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector3::z_axis(const float z) noexcept {
+        return { 0.0f, 0.0f, z };
+    }
+
+    [[nodiscard]] Vector3 Vector3::spherical(const float azimuth, const float inclination, const float radius) noexcept {
+        const float sin_incl = radius * std::sin(inclination);
+        return { sin_incl * std::cos(azimuth), radius * std::cos(inclination), sin_incl * std::sin(azimuth) };
+    }
+
+    [[nodiscard]] Vector3 Vector3::cylindrical(const float angle, const float radius, const float height) noexcept {
+        return { radius * std::cos(angle), height, radius * std::sin(angle) };
+    }
+
+    // element access
+    [[nodiscard]] constexpr float& Vector3::operator[](const int i) noexcept {
+        if consteval {
+            switch (i) {
+                case 0: return x;
+                case 1: return y;
+                default:
+                case 2: return z;
+            }
+        } else { return (&x)[i]; }
+    }
+
+    [[nodiscard]] constexpr const float& Vector3::operator[](const int i) const noexcept {
+        if consteval {
+            switch (i) {
+                case 0: return x;
+                case 1: return y;
+                default:
+                case 2: return z;
+            }
+        } else { return (&x)[i]; }
+    }
+
+    // swizzle
+    [[nodiscard]] constexpr Vector2 Vector3::operator[](const int i0, const int i1) noexcept {
+        if consteval {
+            return { select(i0, x, y, z), select(i1, x, y, z) };
+        } else {
+            return { (&x)[i0], (&x)[i1] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector2 Vector3::operator[](const int i0, const int i1) const noexcept {
+        if consteval {
+            return { select(i0, x, y, z), select(i1, x, y, z) };
+        } else {
+            return { (&x)[i0], (&x)[i1] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector3::operator[](const int i0, const int i1, const int i2) noexcept {
+        if consteval {
+            return { select(i0, x, y, z), select(i1, x, y, z), select(i2, x, y, z) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector3::operator[](const int i0, const int i1, const int i2) const noexcept {
+        if consteval {
+            return { select(i0, x, y, z), select(i1, x, y, z), select(i2, x, y, z) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector3::operator[](const int i0, const int i1, const int i2, const int i3) noexcept {
+        if consteval {
+            return { select(i0, x, y, z), select(i1, x, y, z), select(i2, x, y, z), select(i3, x, y, z)  };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2], (&x)[i3] };
+        }
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector3::operator[](const int i0, const int i1, const int i2, const int i3) const noexcept {
+        if consteval {
+            return { select(i0, x, y, z), select(i1, x, y, z), select(i2, x, y, z), select(i3, x, y, z)  };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2], (&x)[i3] };
+        }
+    }
+
+    // operators
+    constexpr Vector3& Vector3::operator+=(const Vector3& other) noexcept {
+        x += other.x;
+        y += other.y;
+        z += other.z;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator+=(const float other) noexcept {
+        x += other;
+        y += other;
+        z += other;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator-=(const Vector3& other) noexcept {
+        x -= other.x;
+        y -= other.y;
+        z -= other.z;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator-=(const float other) noexcept {
+        x -= other;
+        y -= other;
+        z -= other;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator*=(const Vector3& other) noexcept {
+        x *= other.x;
+        y *= other.y;
+        z *= other.z;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator*=(const float other) noexcept {
+        x *= other;
+        y *= other;
+        z *= other;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator/=(const Vector3& other) noexcept {
+        x /= other.x;
+        y /= other.y;
+        z /= other.z;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator/=(const float other) noexcept {
+        const float inv = 1.0f / other;
+        x *= inv;
+        y *= inv;
+        z *= inv;
+        return *this;
+    }
+
+    constexpr Vector3& Vector3::operator=(const float other) noexcept {
+        x = other;
+        y = other;
+        z = other;
+        return *this;
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector3::operator+() const noexcept {
+        return { x, y, z };
+    }
+
+    [[nodiscard]] constexpr Vector3 Vector3::operator-() const noexcept {
+        return { -x, -y, -z };
+    }
+    
+    [[nodiscard]] constexpr Vector3 operator+(const Vector3& a, const Vector3& b) noexcept {
+        return { a.x+b.x, a.y+b.y, a.z+b.z };
+    }
+
+    [[nodiscard]] constexpr Vector3 operator+(const Vector3& a, const float b) noexcept {
+        return { a.x+b, a.y+b, a.z+b };
+    }
+
+    [[nodiscard]] constexpr Vector3 operator+(const float a, const Vector3& b) noexcept {
+        return b+a;
+    }
+
+    [[nodiscard]] constexpr Vector3 operator-(const Vector3& a, const Vector3& b) noexcept {
+        return { a.x-b.x, a.y-b.y, a.z-b.z };
+    }
+
+    [[nodiscard]] constexpr Vector3 operator-(const Vector3& a, const float b) noexcept {
+        return { a.x-b, a.y-b, a.z-b };
+    }
+
+    [[nodiscard]] constexpr Vector3 operator-(const float a, const Vector3& b) noexcept {
+        return { a-b.x, a-b.y, a-b.z };
+    }
+
+    [[nodiscard]] constexpr Vector3 operator*(const Vector3& a, const Vector3& b) noexcept {
+        return { a.x*b.x, a.y*b.y, a.z*b.z };
+    }
+
+    [[nodiscard]] constexpr Vector3 operator*(const Vector3& a, const float b) noexcept {
+        return { a.x*b, a.y*b, a.z*b };
+    }
+
+    [[nodiscard]] constexpr Vector3 operator*(const float a, const Vector3& b) noexcept {
+        return b*a;
+    }
+
+    [[nodiscard]] constexpr Vector3 operator/(const Vector3& a, const Vector3& b) noexcept {
+        return { a.x/b.x, a.y/b.y, a.z/b.z };
+    }
+    
+    [[nodiscard]] constexpr Vector3 operator/(const Vector3& a, const float b) noexcept {
+        return a * (1.0f / b);
+    }
+    
+    [[nodiscard]] constexpr Vector3 operator/(const float a, const Vector3& b) noexcept {
+        return { a/b.x, a/b.y, a/b.z };
+    }
+
+    // functions
+
+    // Returns dot product
+    [[nodiscard]] constexpr float dot(const Vector3& a, const Vector3& b) noexcept {
+        return a.x*b.x + a.y*b.y + a.z*b.z;
+    }
+
+    // Returns squared magnitude
+    [[nodiscard]] constexpr float length_sq(const Vector3& v) noexcept {
+        return dot(v, v);
+    }
+
+    // Returns magnitude
+    [[nodiscard]] float length(const Vector3& v) noexcept {
+        return std::sqrt(length_sq(v));
+    }
+
+    // Return squared distance between two vectors
+    [[nodiscard]] constexpr float distance_sq(const Vector3& a, const Vector3& b) noexcept {
+        return length_sq(a - b);
+    }
+
+    // Returns distance between two vectors
+    [[nodiscard]] float distance(const Vector3& a, const Vector3& b) noexcept {
+        return length(a - b);
+    }
+
+    // Safe normalize, checks length
+    [[nodiscard]] Vector3 normalize(const Vector3& v) noexcept {
+        const float len = length(v);
+
+        return (len > CMP_NORMALIZE_TOLERANCE) ? v / len : Vector3();
+    }
+    
+    // Faster normalize, it presupposes vector has non-zero length
+    // TODO: add check that v is non-zero on debug builds
+    [[nodiscard]] Vector3 normalize_fast(const Vector3& v) noexcept {
+        return v / length(v);
+    }
+
+    // Returns vector projected onto normal
+    [[nodiscard]] constexpr Vector3 project(const Vector3& vector, const Vector3& normal) noexcept {
+        return normal * (dot(vector, normal) / length_sq(normal));
+    }
+
+    // Returns a vector reflected off a plane defined by its normal
+    [[nodiscard]] constexpr Vector3 reflect(const Vector3& incoming, const Vector3& normal) noexcept {
+        return incoming - 2.0f * dot(incoming, normal) * normal;
+    }
+
+    // Returns the angle between two vectors
+    [[nodiscard]] float angle(const Vector3& a, const Vector3& b) noexcept {
+        return std::acos(dot(a, b) / (length(a) * length(b)));
+    }
+
+    // Returns linear interpolation between two vectors
+    [[nodiscard]] constexpr Vector3 lerp(const Vector3& from, const Vector3& to, const float weight) noexcept {
+        return {
+            lerp(from.x, to.x, weight),
+            lerp(from.y, to.y, weight),
+            lerp(from.z, to.z, weight)
+        };
+    }
+
+    // Returns component-wise minimum
+    [[nodiscard]] constexpr Vector3 min(const Vector3& a, const Vector3& b) noexcept {
+        return {
+            std::min(a.x, b.x),
+            std::min(a.y, b.y),
+            std::min(a.z, b.z)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector3 min(const Vector3& a, const float b) noexcept {
+        return {
+            std::min(a.x, b),
+            std::min(a.y, b),
+            std::min(a.z, b)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector3 min(const float a, const Vector3& b) noexcept {
+        return min(b, a);
+    }
+
+    // Returns the vector with the smaller length
+    [[nodiscard]] constexpr Vector3 min_length(const Vector3& a, const Vector3& b) noexcept {
+        return length_sq(a) < length_sq(b) ? a : b;
+    }
+
+    // Returns a vector in the same direction whose length is bounded above by the given value.
+    [[nodiscard]] Vector3 min_length(const Vector3& a, const float b) noexcept {
+        const float len_sq = length_sq(a);
+        
+        if (len_sq > b * b) {
+            return a * (b / std::sqrt(len_sq));
+        } else {
+            return a;
+        }
+    }
+
+    [[nodiscard]] Vector3 min_length(const float a, const Vector3& b) noexcept {
+        return min_length(b, a);
+    }
+
+    // Returns component-wise maximum
+    [[nodiscard]] constexpr Vector3 max(const Vector3& a, const Vector3& b) noexcept {
+        return {
+            std::max(a.x, b.x),
+            std::max(a.y, b.y),
+            std::max(a.z, b.z)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector3 max(const Vector3& a, const float b) noexcept {
+        return {
+            std::max(a.x, b),
+            std::max(a.y, b),
+            std::max(a.z, b)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector3 max(const float a, const Vector3& b) noexcept {
+        return max(b, a);
+    }
+
+    // Returns the vector with the larger length
+    [[nodiscard]] constexpr Vector3 max_length(const Vector3& a, const Vector3& b) noexcept {
+        return length_sq(a) > length_sq(b) ? a : b;
+    }
+
+    // Returns a vector in the same direction whose length is bounded below by the given value. Returns the 0 vector if the vector is too small to be normalized.
+    [[nodiscard]] Vector3 max_length(const Vector3& a, const float b) noexcept {
+        const float len_sq = length_sq(a);
+        
+        if (len_sq <= CMP_NORMALIZE_TOLERANCE2) {
+            return Vector3();
+        } else if (len_sq < b * b) {
+            return a * (b / std::sqrt(len_sq));
+        } else {
+            return a;
+        }
+    }
+
+    [[nodiscard]] Vector3 max_length(const float a, const Vector3& b) noexcept {
+        return max_length(b, a);
+    }
+
+    // Clamps each component of x to the range [x_min, x_max]. Presupposes x_min <= x_max.
+    [[nodiscard]] constexpr Vector3 clamp(const Vector3& x, const Vector3& x_min, const Vector3& x_max) noexcept {
+        return max(x_min, min(x, x_max));
+    }
+
+    [[nodiscard]] constexpr Vector3 clamp(const Vector3& x, const float x_min, const float x_max) noexcept {
+        return max(x_min, min(x, x_max));
+    }
+
+    // Clamps the length of the vector to the range [x_min, x_max]. Presupposes x_min <= x_max. Returns the 0 vector if the vector is too small to be normalized.
+    [[nodiscard]] Vector3 clamp_length(const Vector3& v, const float x_min, const float x_max) noexcept {
+        const float len_sq = length_sq(v);
+        
+        if (len_sq <= CMP_NORMALIZE_TOLERANCE2) {
+            return Vector3();
+        } else if (len_sq < x_min * x_min) {
+            return v * (x_min / std::sqrt(len_sq));
+        } else if (len_sq > x_max * x_max) {
+            return v * (x_max / std::sqrt(len_sq));
+        } else {
+            return v;
+        }
+    }
+
+    // Returns component-wise absolute value
+    [[nodiscard]] constexpr Vector3 abs(const Vector3& v) noexcept {
+        return {
+            abs(v.x),
+            abs(v.y),
+            abs(v.z)
+        };
+    }
+
+    // Returns component-wise floor
+    [[nodiscard]] constexpr Vector3 floor(const Vector3& v) noexcept {
+        return {
+            floor(v.x),
+            floor(v.y),
+            floor(v.z)
+        };
+    }
+
+    // Returns component-wise ceiling
+    [[nodiscard]] constexpr Vector3 ceil(const Vector3& v) noexcept {
+        return {
+            ceil(v.x),
+            ceil(v.y),
+            ceil(v.z)
+        };
+    }
+
+    // Returns component-wise truncation
+    [[nodiscard]] constexpr Vector3 trunc(const Vector3& v) noexcept {
+        return {
+            trunc(v.x),
+            trunc(v.y),
+            trunc(v.z)
+        };
+    }
+
+    // Returns component-wise round
+    [[nodiscard]] constexpr Vector3 round(const Vector3& v) noexcept {
+        return {
+            round(v.x),
+            round(v.y),
+            round(v.z)
+        };
+    }
+
+    // Returns component-wise sign. Note that -0 still returns 0
+    [[nodiscard]] constexpr Vector3 sign(const Vector3& v) noexcept {
+        return {
+            sign(v.x),
+            sign(v.y),
+            sign(v.z)
+        };
+    }
+
+    // Returns true if the vectors are approximately equal
+    [[nodiscard]] constexpr bool approx_eq(const Vector3& a, const Vector3& b) noexcept {
+        return distance_sq(a, b) < CMP_EPSILON2;
+    }
+
+    // Returns cross product
+    [[nodiscard]] constexpr Vector3 cross(const Vector3& a, const Vector3& b) noexcept {
+        return { a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x };
+    }
+}
+
+export namespace std {
+    template<> struct formatter<draco::math::Vector3> : formatter<float> {
+        auto format(const draco::math::Vector3& v, format_context& ctx) const {
+            ctx.advance_to(format_to(ctx.out(), "{{"));
+
+            for (int i = 0; i < 3; ++i) {
+                if (i) ctx.advance_to(format_to(ctx.out(), ", "));
+                ctx.advance_to(formatter<float>::format(v[i], ctx));
+            }
+
+            return format_to(ctx.out(), "}}");
+        }
+    };
+}
\ No newline at end of file
diff --git a/engine/native/core/math/vector4.cppm b/engine/native/core/math/vector4.cppm
index df2d8ec..fce1cbc 100644
--- a/engine/native/core/math/vector4.cppm
+++ b/engine/native/core/math/vector4.cppm
@@ -1,5 +1,8 @@
 module;
 
+#include <cmath>
+#include <algorithm>
+#include <format>
 #include "platform/simd.h"
 
 #if ARCH_X64
@@ -8,105 +11,245 @@ module;
     #include <arm_neon.h>
 #endif
 
-#include <cmath>
-#include <type_traits>
+export module core.math.types:vector4;
+export import :common;
 
-export module core.math.vector4;
+import core.math.constants;
+import core.math.functions;
 import core.defs;
 
 export namespace draco::math {
+    // assertions
+    static_assert(sizeof(Vector4) == 16, "Vector4 must be 16 bytes");
+    static_assert(alignof(Vector4) == 16, "Vector4 must be 16-byte aligned");
+    static_assert(trivial<Vector4>, "Vector4 must be trivial");
+    static_assert(std::is_standard_layout_v<Vector4>, "Vector4 must be standard layout");
 
-    struct alignas(16) Vector4 {
-        float x, y, z, w;
+    // constructors
+    [[nodiscard]] constexpr Vector4::Vector4(const float n) noexcept
+        : x{n}, y{n}, z{n}, w{n} { }
 
-        // constructors.
-        constexpr Vector4() noexcept = default;
-        constexpr Vector4(const float x, const float y, const float z, const float w) noexcept
+    [[nodiscard]] constexpr Vector4::Vector4(const float x, const float y, const float z, const float w) noexcept
         : x{x}, y{y}, z{z}, w{w} { }
 
-        // element access.
-        constexpr float& operator[](const int i) noexcept {
-            if consteval {
-                switch (i) {
-                    case 0: return x;
-                    case 1: return y;
-                    case 2: return z;
-                    default:
-                    case 3: return w;
-                }
-            } else { return (&x)[i]; }
-        }
+    [[nodiscard]] constexpr Vector4::Vector4(const Vector2& xy) noexcept
+        : x{xy.x}, y{xy.y}, z{0.0f}, w{0.0f} { }
 
-        constexpr const float& operator[](const int i) const noexcept {
-            if consteval {
-                switch (i) {
-                    case 0: return x;
-                    case 1: return y;
-                    case 2: return z;
-                    default:
-                    case 3: return w;
-                }
-            } else { return (&x)[i]; }
-        }
+    [[nodiscard]] constexpr Vector4::Vector4(const Vector2& xy, const float z, const float w) noexcept
+        : x{xy.x}, y{xy.y}, z{z}, w{w} { }
+
+    [[nodiscard]] constexpr Vector4::Vector4(const float x, const Vector2& yz, const float w) noexcept
+        : x{x}, y{yz.x}, z{yz.y}, w{w} { }
+
+    [[nodiscard]] constexpr Vector4::Vector4(const float x, const float y, const Vector2& zw) noexcept
+        : x{x}, y{y}, z{zw.x}, w{zw.y} { }
+
+    [[nodiscard]] constexpr Vector4::Vector4(const Vector2& xy, const Vector2& zw) noexcept
+        : x{xy.x}, y{xy.y}, z{zw.x}, w{zw.y} { }
+
+    [[nodiscard]] constexpr Vector4::Vector4(const Vector3& xyz, const float w) noexcept
+        : x{xyz.x}, y{xyz.y}, z{xyz.z}, w{w} { }
+
+    [[nodiscard]] constexpr Vector4::Vector4(const float x, const Vector3& yzw) noexcept
+        : x{x}, y{yzw.x}, z{yzw.y}, w{yzw.z} { }
+    
+    // static
+    [[nodiscard]] constexpr Vector4 Vector4::x_axis(const float x) noexcept {
+        return { x, 0.0f, 0.0f, 0.0f };
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector4::y_axis(const float y) noexcept {
+        return { 0.0f, y, 0.0f, 0.0f };
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector4::z_axis(const float z) noexcept {
+        return { 0.0f, 0.0f, z, 0.0f };
+    }
 
-        [[nodiscard]] constexpr bool operator==(const Vector4& other) const noexcept = default;
+    [[nodiscard]] constexpr Vector4 Vector4::w_axis(const float w) noexcept {
+        return { 0.0f, 0.0f, 0.0f, w };
+    }
 
-        constexpr Vector4& operator+=(const Vector4& other) noexcept {
-            x += other.x;
-            y += other.y;
-            z += other.z;
-            w += other.w;
-            return *this;
+    // element access
+    [[nodiscard]] constexpr float& Vector4::operator[](const int i) noexcept {
+        if consteval {
+            switch (i) {
+                case 0: return x;
+                case 1: return y;
+                case 2: return z;
+                default:
+                case 3: return w;
+            }
+        } else { return (&x)[i]; }
+    }
+
+    [[nodiscard]] constexpr const float& Vector4::operator[](const int i) const noexcept {
+        if consteval {
+            switch (i) {
+                case 0: return x;
+                case 1: return y;
+                case 2: return z;
+                default:
+                case 3: return w;
+            }
+        } else { return (&x)[i]; }
+    }
+
+    // swizzle
+    [[nodiscard]] constexpr Vector2 Vector4::operator[](const int i0, const int i1) noexcept {
+        if consteval {
+            return { select(i0, x, y, z, w), select(i1, x, y, z, w) };
+        } else {
+            return { (&x)[i0], (&x)[i1] };
         }
+    }
 
-        constexpr Vector4& operator-=(const Vector4& other) noexcept {
-            x -= other.x;
-            y -= other.y;
-            z -= other.z;
-            w -= other.w;
-            return *this;
+    [[nodiscard]] constexpr Vector2 Vector4::operator[](const int i0, const int i1) const noexcept {
+        if consteval {
+            return { select(i0, x, y, z, w), select(i1, x, y, z, w) };
+        } else {
+            return { (&x)[i0], (&x)[i1] };
         }
+    }
 
-        constexpr Vector4& operator*=(const Vector4& other) noexcept {
-            x *= other.x;
-            y *= other.y;
-            z *= other.z;
-            w *= other.w;
-            return *this;
+    [[nodiscard]] constexpr Vector3 Vector4::operator[](const int i0, const int i1, const int i2) noexcept {
+        if consteval {
+            return { select(i0, x, y, z, w), select(i1, x, y, z, w), select(i2, x, y, z, w) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2] };
         }
+    }
 
-        constexpr Vector4& operator*=(const float s) noexcept {
-            x *= s;
-            y *= s;
-            z *= s;
-            w *= s;
-            return *this;
+    [[nodiscard]] constexpr Vector3 Vector4::operator[](const int i0, const int i1, const int i2) const noexcept {
+        if consteval {
+            return { select(i0, x, y, z, w), select(i1, x, y, z, w), select(i2, x, y, z, w) };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2] };
         }
+    }
 
-        constexpr Vector4& operator/=(const float s) noexcept {
-            const float inv = 1.0f / s;
-            x *= inv;
-            y *= inv;
-            z *= inv;
-            w *= inv;
-            return *this;
+    [[nodiscard]] constexpr Vector4 Vector4::operator[](const int i0, const int i1, const int i2, const int i3) noexcept {
+        if consteval {
+            return { select(i0, x, y, z, w), select(i1, x, y, z, w), select(i2, x, y, z, w), select(i3, x, y, z, w)  };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2], (&x)[i3] };
         }
-    };
+    }
 
-    // safety features go brr
-    static_assert(sizeof(Vector4) == 16, "Vector4 must be 16 bytes");
-    static_assert(alignof(Vector4) == 16, "Vector4 must be 16-byte aligned");
-    static_assert(trivial<Vector4>, "Vector4 must be trivial");
-    static_assert(std::is_standard_layout_v<Vector4>, "Vector4 must be standard layout");
+    [[nodiscard]] constexpr Vector4 Vector4::operator[](const int i0, const int i1, const int i2, const int i3) const noexcept {
+        if consteval {
+            return { select(i0, x, y, z, w), select(i1, x, y, z, w), select(i2, x, y, z, w), select(i3, x, y, z, w)  };
+        } else {
+            return { (&x)[i0], (&x)[i1], (&x)[i2], (&x)[i3] };
+        }
+    }
+
+    // operators
+    constexpr Vector4& Vector4::operator+=(const Vector4& other) noexcept {
+        x += other.x;
+        y += other.y;
+        z += other.z;
+        w += other.w;
+        return *this;
+    }
+
+    constexpr Vector4& Vector4::operator+=(const float other) noexcept {
+        x += other;
+        y += other;
+        z += other;
+        w += other;
+        return *this;
+    }
+
+    constexpr Vector4& Vector4::operator-=(const Vector4& other) noexcept {
+        x -= other.x;
+        y -= other.y;
+        z -= other.z;
+        w -= other.w;
+        return *this;
+    }
+
+    constexpr Vector4& Vector4::operator-=(const float other) noexcept {
+        x -= other;
+        y -= other;
+        z -= other;
+        w -= other;
+        return *this;
+    }
+
+    constexpr Vector4& Vector4::operator*=(const Vector4& other) noexcept {
+        x *= other.x;
+        y *= other.y;
+        z *= other.z;
+        w *= other.w;
+        return *this;
+    }
+
+    constexpr Vector4& Vector4::operator*=(const float other) noexcept {
+        x *= other;
+        y *= other;
+        z *= other;
+        w *= other;
+        return *this;
+    }
+
+    constexpr Vector4& Vector4::operator/=(const Vector4& other) noexcept {
+        x /= other.x;
+        y /= other.y;
+        z /= other.z;
+        w /= other.w;
+        return *this;
+    }
+
+    constexpr Vector4& Vector4::operator/=(const float other) noexcept {
+        const float inv = 1.0f / other;
+        x *= inv;
+        y *= inv;
+        z *= inv;
+        w *= inv;
+        return *this;
+    }
 
+    constexpr Vector4& Vector4::operator=(const float other) noexcept {
+        x = other;
+        y = other;
+        z = other;
+        w = other;
+        return *this;
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector4::operator+() const noexcept {
+        return { x, y, z, w };
+    }
+
+    [[nodiscard]] constexpr Vector4 Vector4::operator-() const noexcept {
+        return { -x, -y, -z, -w };
+    }
+    
     [[nodiscard]] constexpr Vector4 operator+(const Vector4& a, const Vector4& b) noexcept {
         return { a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w };
     }
 
+    [[nodiscard]] constexpr Vector4 operator+(const Vector4& a, const float b) noexcept {
+        return { a.x+b, a.y+b, a.z+b, a.w+b };
+    }
+
+    [[nodiscard]] constexpr Vector4 operator+(const float a, const Vector4& b) noexcept {
+        return b+a;
+    }
+
     [[nodiscard]] constexpr Vector4 operator-(const Vector4& a, const Vector4& b) noexcept {
         return { a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w };
     }
 
+    [[nodiscard]] constexpr Vector4 operator-(const Vector4& a, const float b) noexcept {
+        return { a.x-b, a.y-b, a.z-b, a.w-b };
+    }
+
+    [[nodiscard]] constexpr Vector4 operator-(const float a, const Vector4& b) noexcept {
+        return { a-b.x, a-b.y, a-b.z, a-b.w };
+    }
+
     [[nodiscard]] constexpr Vector4 operator*(const Vector4& a, const Vector4& b) noexcept {
         return { a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w };
     }
@@ -115,62 +258,300 @@ export namespace draco::math {
         return { a.x*b, a.y*b, a.z*b, a.w*b };
     }
 
-    [[nodiscard]] constexpr Vector4 operator*(const float s, const Vector4& v) noexcept {
-        return v*s;
+    [[nodiscard]] constexpr Vector4 operator*(const float a, const Vector4& b) noexcept {
+        return b*a;
     }
 
-    [[nodiscard]] constexpr Vector4 operator/(const Vector4& v, const float s) noexcept {
-        return v * (1.f/s);
+    [[nodiscard]] constexpr Vector4 operator/(const Vector4& a, const Vector4& b) noexcept {
+        return { a.x/b.x, a.y/b.y, a.z/b.z, a.w/b.w };
+    }
+    
+    [[nodiscard]] constexpr Vector4 operator/(const Vector4& a, const float b) noexcept {
+        return a * (1.0f / b);
+    }
+    
+    [[nodiscard]] constexpr Vector4 operator/(const float a, const Vector4& b) noexcept {
+        return { a/b.x, a/b.y, a/b.z, a/b.w };
     }
 
-    [[nodiscard]] FORCEINLINE float dot(const Vector4 &a, const Vector4 &b) noexcept {
-        #if ARCH_X64
-            // There's only 4 floats, so SSE is what we will use.
-            // If there's a situation with multiple dot calls, we can setup a
-            // way to call 8 / 16 / 32 floats, but over-head could upset gains.
-            // Be sure it occurs commonly enough to matter here.
-            // Shuffle-first reduction worked best here.
-            __m128 va = _mm_load_ps(&a.x);
-            __m128 vb = _mm_load_ps(&b.x);
-
-            __m128 m = _mm_mul_ps(va, vb);
-
-            __m128 shuf = _mm_movehdup_ps(m);
-            __m128 sum  = _mm_add_ps(m, shuf);
-
-            shuf = _mm_movehl_ps(shuf, sum);
-            sum = _mm_add_ss(sum, shuf);
+    // functions
+    
+    // Returns dot product
+    [[nodiscard]] constexpr float dot(const Vector4& a, const Vector4& b) noexcept {
+        if !consteval {
+            #if ARCH_X64
+                // There's only 4 floats, so SSE is what we will use.
+                // If there's a situation with multiple dot calls, we can setup a
+                // way to call 8 / 16 / 32 floats, but over-head could upset gains.
+                // Be sure it occurs commonly enough to matter here.
+                // Shuffle-first reduction worked best here.
+                __m128 va = _mm_load_ps(&a.x);
+                __m128 vb = _mm_load_ps(&b.x);
+
+                __m128 m = _mm_mul_ps(va, vb);
+
+                __m128 shuf = _mm_movehdup_ps(m);
+                __m128 sum  = _mm_add_ps(m, shuf);
+
+                shuf = _mm_movehl_ps(shuf, sum);
+                sum = _mm_add_ss(sum, shuf);
+
+                return _mm_cvtss_f32(sum);
+            #elif ARCH_ARM64
+                #error "ARM64 NEON support not yet implemented."
+            #endif
+        }
 
-            return _mm_cvtss_f32(sum);
-        #elif ARCH_ARM64
-            #error "ARM64 NEON support not yet implemented."
-        #else
-            // scalar.
-            return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
-        #endif
+        return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
     }
 
-    // Returns squared magnitude.
-    [[nodiscard]] FORCEINLINE float length_sq(const Vector4& v) noexcept {
+    // Returns squared magnitude
+    [[nodiscard]] constexpr float length_sq(const Vector4& v) noexcept {
         return dot(v, v);
     }
 
     // Returns magnitude
-    [[nodiscard]] FORCEINLINE float length(const Vector4& v) noexcept {
+    [[nodiscard]] float length(const Vector4& v) noexcept {
         return std::sqrt(length_sq(v));
     }
 
-    // Safe normalize, checks length.
-    [[nodiscard]] FORCEINLINE Vector4 normalize(const Vector4& v) noexcept {
-        const float len = length(v);
+    // Return squared distance between two vectors
+    [[nodiscard]] constexpr float distance_sq(const Vector4& a, const Vector4& b) noexcept {
+        return length_sq(a - b);
+    }
 
-        return (len > 0.0f) ? v / len : Vector4{0,0,0,0};
+    // Returns distance between two vectors
+    [[nodiscard]] float distance(const Vector4& a, const Vector4& b) noexcept {
+        return length(a - b);
     }
 
-    // Faster normalize, it presupposes vector has non-zero length.
-    [[nodiscard]] FORCEINLINE Vector4 normalize_fast(const Vector4& v) noexcept {
+    // Safe normalize, checks length
+    [[nodiscard]] Vector4 normalize(const Vector4& v) noexcept {
+        const float len = length(v);
+
+        return (len > CMP_NORMALIZE_TOLERANCE) ? v / len : Vector4();
+    }
+    
+    // Faster normalize, it presupposes vector has non-zero length
+    // TODO: add check that v is non-zero on debug builds
+    [[nodiscard]] Vector4 normalize_fast(const Vector4& v) noexcept {
         return v / length(v);
     }
 
+    // Returns vector projected onto normal
+    [[nodiscard]] constexpr Vector4 project(const Vector4& vector, const Vector4& normal) noexcept {
+        return normal * (dot(vector, normal) / length_sq(normal));
+    }
+
+    // Returns a vector reflected off a plane defined by its normal
+    [[nodiscard]] constexpr Vector4 reflect(const Vector4& incoming, const Vector4& normal) noexcept {
+        return incoming - 2.0f * dot(incoming, normal) * normal;
+    }
+
+    // Returns the angle between two vectors
+    [[nodiscard]] float angle(const Vector4& a, const Vector4& b) noexcept {
+        return std::acos(dot(a, b) / (length(a) * length(b)));
+    }
+
+    // Returns linear interpolation between two vectors
+    [[nodiscard]] constexpr Vector4 lerp(const Vector4& from, const Vector4& to, const float weight) noexcept {
+        return {
+            lerp(from.x, to.x, weight),
+            lerp(from.y, to.y, weight),
+            lerp(from.z, to.z, weight),
+            lerp(from.w, to.w, weight)
+        };
+    }
+
+    // Returns component-wise minimum
+    [[nodiscard]] constexpr Vector4 min(const Vector4& a, const Vector4& b) noexcept {
+        return {
+            std::min(a.x, b.x),
+            std::min(a.y, b.y),
+            std::min(a.z, b.z),
+            std::min(a.w, b.w)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector4 min(const Vector4& a, const float b) noexcept {
+        return {
+            std::min(a.x, b),
+            std::min(a.y, b),
+            std::min(a.z, b),
+            std::min(a.w, b)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector4 min(const float a, const Vector4& b) noexcept {
+        return min(b, a);
+    }
+
+    // Returns the vector with the smaller length
+    [[nodiscard]] constexpr Vector4 min_length(const Vector4& a, const Vector4& b) noexcept {
+        return length_sq(a) < length_sq(b) ? a : b;
+    }
+
+    // Returns a vector in the same direction whose length is bounded above by the given value.
+    [[nodiscard]] Vector4 min_length(const Vector4& a, const float b) noexcept {
+        const float len_sq = length_sq(a);
+        
+        if (len_sq > b * b) {
+            return a * (b / std::sqrt(len_sq));
+        } else {
+            return a;
+        }
+    }
+
+    [[nodiscard]] Vector4 min_length(const float a, const Vector4& b) noexcept {
+        return min_length(b, a);
+    }
+
+    // Returns component-wise maximum
+    [[nodiscard]] constexpr Vector4 max(const Vector4& a, const Vector4& b) noexcept {
+        return {
+            std::max(a.x, b.x),
+            std::max(a.y, b.y),
+            std::max(a.z, b.z),
+            std::max(a.w, b.w)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector4 max(const Vector4& a, const float b) noexcept {
+        return {
+            std::max(a.x, b),
+            std::max(a.y, b),
+            std::max(a.z, b),
+            std::max(a.w, b)
+        };
+    }
+
+    [[nodiscard]] constexpr Vector4 max(const float a, const Vector4& b) noexcept {
+        return max(b, a);
+    }
+
+    // Returns the vector with the larger length
+    [[nodiscard]] constexpr Vector4 max_length(const Vector4& a, const Vector4& b) noexcept {
+        return length_sq(a) > length_sq(b) ? a : b;
+    }
+
+    // Returns a vector in the same direction whose length is bounded below by the given value. Returns the 0 vector if the vector is too small to be normalized.
+    [[nodiscard]] Vector4 max_length(const Vector4& a, const float b) noexcept {
+        const float len_sq = length_sq(a);
+        
+        if (len_sq <= CMP_NORMALIZE_TOLERANCE2) {
+            return Vector4();
+        } else if (len_sq < b * b) {
+            return a * (b / std::sqrt(len_sq));
+        } else {
+            return a;
+        }
+    }
+
+    [[nodiscard]] Vector4 max_length(const float a, const Vector4& b) noexcept {
+        return max_length(b, a);
+    }
+
+    // Clamps each component of x to the range [x_min, x_max]. Presupposes x_min <= x_max.
+    [[nodiscard]] constexpr Vector4 clamp(const Vector4& x, const Vector4& x_min, const Vector4& x_max) noexcept {
+        return max(x_min, min(x, x_max));
+    }
+
+    [[nodiscard]] constexpr Vector4 clamp(const Vector4& x, const float x_min, const float x_max) noexcept {
+        return max(x_min, min(x, x_max));
+    }
+
+    // Clamps the length of the vector to the range [x_min, x_max]. Presupposes x_min <= x_max. Returns the 0 vector if the vector is too small to be normalized.
+    [[nodiscard]] Vector4 clamp_length(const Vector4& v, const float x_min, const float x_max) noexcept {
+        const float len_sq = length_sq(v);
+        
+        if (len_sq <= CMP_NORMALIZE_TOLERANCE2) {
+            return Vector4();
+        } else if (len_sq < x_min * x_min) {
+            return v * (x_min / std::sqrt(len_sq));
+        } else if (len_sq > x_max * x_max) {
+            return v * (x_max / std::sqrt(len_sq));
+        } else {
+            return v;
+        }
+    }
+
+    // Returns component-wise absolute value
+    [[nodiscard]] constexpr Vector4 abs(const Vector4& v) noexcept {
+        return {
+            abs(v.x),
+            abs(v.y),
+            abs(v.z),
+            abs(v.w)
+        };
+    }
 
+    // Returns component-wise floor
+    [[nodiscard]] constexpr Vector4 floor(const Vector4& v) noexcept {
+        return {
+            floor(v.x),
+            floor(v.y),
+            floor(v.z),
+            floor(v.w)
+        };
+    }
+
+    // Returns component-wise ceiling
+    [[nodiscard]] constexpr Vector4 ceil(const Vector4& v) noexcept {
+        return {
+            ceil(v.x),
+            ceil(v.y),
+            ceil(v.z),
+            ceil(v.w)
+        };
+    }
+
+    // Returns component-wise truncation
+    [[nodiscard]] constexpr Vector4 trunc(const Vector4& v) noexcept {
+        return {
+            trunc(v.x),
+            trunc(v.y),
+            trunc(v.z),
+            trunc(v.w)
+        };
+    }
+
+    // Returns component-wise round
+    [[nodiscard]] constexpr Vector4 round(const Vector4& v) noexcept {
+        return {
+            round(v.x),
+            round(v.y),
+            round(v.z),
+            round(v.w)
+        };
+    }
+
+    // Returns component-wise sign. Note that -0 still returns 0
+    [[nodiscard]] constexpr Vector4 sign(const Vector4& v) noexcept {
+        return {
+            sign(v.x),
+            sign(v.y),
+            sign(v.z),
+            sign(v.w)
+        };
+    }
+
+    // Returns true if the vectors are approximately equal
+    [[nodiscard]] constexpr bool approx_eq(const Vector4& a, const Vector4& b) noexcept {
+        return distance_sq(a, b) < CMP_EPSILON2;
+    }
 } // namespace draco::math
+
+export namespace std {
+    template<> struct formatter<draco::math::Vector4> : formatter<float> {
+        auto format(const draco::math::Vector4& v, format_context& ctx) const {
+            ctx.advance_to(format_to(ctx.out(), "{{"));
+
+            for (int i = 0; i < 4; ++i) {
+                if (i) ctx.advance_to(format_to(ctx.out(), ", "));
+                ctx.advance_to(formatter<float>::format(v[i], ctx));
+            }
+
+            return format_to(ctx.out(), "}}");
+        }
+    };
+}