diff --git a/modules/CMakeLists.txt b/modules/CMakeLists.txt
index 541daae..6ab5d51 100644
--- a/modules/CMakeLists.txt
+++ b/modules/CMakeLists.txt
@@ -41,7 +41,11 @@ add_module(
     DEPENDENCIES
     preliminary
     TESTS
+    trig.test.cpp
     vec2.test.cpp
+    vec3.test.cpp
+    vec4.test.cpp
+    mat4.test.cpp
 )
 
 add_module(
diff --git a/modules/math/algebra.cppm b/modules/math/algebra.cppm
index ffe8bd8..67176e0 100644
--- a/modules/math/algebra.cppm
+++ b/modules/math/algebra.cppm
@@ -32,12 +32,12 @@ export namespace lt::math {
  *
  * the 1 at [z][3] is to save the Z axis into the resulting W for perspective division.
  *
- * @ref Thanks to pikuma for explaining the math behind this:
+ * @ref Thanks to @pikuma for explaining the math behind this:
  * https://www.youtube.com/watch?v=EqNcqBdrNyI
  */
 template<typename T>
     requires(std::is_arithmetic_v<T>)
-constexpr auto perspective(T field_of_view, T aspect_ratio, T z_near, T z_far)
+constexpr auto perspective(T field_of_view, T aspect_ratio, T z_near, T z_far) -> mat4_impl<T>
 {
 	const T half_fov_tan = std::tan(field_of_view / static_cast<T>(2));
 
diff --git a/modules/math/mat4.cppm b/modules/math/mat4.cppm
index 7d7ca79..d68e7b3 100644
--- a/modules/math/mat4.cppm
+++ b/modules/math/mat4.cppm
@@ -7,13 +7,20 @@ import math.vec4;
 
 export namespace lt::math {
 
+/** A 4 by 4 matrix, column major order
+ *
+ * @todo(Light): Use std::simd when it's implemented. */
 template<typename T = f32>
     requires(std::is_arithmetic_v<T>)
 struct mat4_impl
 {
 	using Column_T = vec4_impl<T>;
 
-	constexpr explicit mat4_impl(T scalar = 0)
+	using Underlying_T = Column_T::Underlying_T;
+
+	static constexpr auto num_elements = 4u * 4u;
+
+	constexpr explicit mat4_impl(T scalar = T {})
 	    : values(
 	          {
 	              Column_T { scalar },
@@ -25,13 +32,12 @@ struct mat4_impl
 	{
 	}
 
-	// clang-format off
-    constexpr mat4_impl(
-        const T&  x0, const T&  y0, const T&  z0, const T&  w0,
-        const T&  x1, const T&  y1, const T&  z1, const T&  w1,
-        const T&  x2, const T&  y2, const T&  z2, const T&  w2,
-        const T&  x3, const T&  y3, const T&  z3, const T&  w3
-    )
+	constexpr mat4_impl(
+	    // clang-format off
+        const T& x0, const T& x1, const T& x2, const T& x3,
+        const T& y0, const T& y1, const T& y2, const T& y3,
+        const T& z0, const T& z1, const T& z2, const T& z3,
+        const T& w0, const T& w1, const T& w2, const T& w3)
 	    // clang-format on
 	    : values({ { x0, x1, x2, x3 }, { y0, y1, y2, y3 }, { z0, z1, z2, z3 }, { w0, w1, w2, w3 } })
 	{
@@ -57,8 +63,41 @@ struct mat4_impl
 		};
 	}
 
+	[[nodiscard]] constexpr auto operator*(const mat4_impl<T> &other) const -> mat4_impl<T>
+	{
+		const auto &[a_x, a_y, a_z, a_w] = values;
+		const auto &[b_x, b_y, b_z, b_w] = other.values;
+
+		return mat4_impl<T>(
+		    // X column
+		    a_x.x * b_x.x + a_y.x * b_x.y + a_z.x * b_x.z + a_w.x * b_x.w,
+		    a_x.y * b_x.x + a_y.y * b_x.y + a_z.y * b_x.z + a_w.y * b_x.w,
+		    a_x.z * b_x.x + a_y.z * b_x.y + a_z.z * b_x.z + a_w.z * b_x.w,
+		    a_x.w * b_x.x + a_y.w * b_x.y + a_z.w * b_x.z + a_w.w * b_x.w,
+
+		    // Y column
+		    a_x.x * b_y.x + a_y.x * b_y.y + a_z.x * b_y.z + a_w.x * b_y.w,
+		    a_x.y * b_y.x + a_y.y * b_y.y + a_z.y * b_y.z + a_w.y * b_y.w,
+		    a_x.z * b_y.x + a_y.z * b_y.y + a_z.z * b_y.z + a_w.z * b_y.w,
+		    a_x.w * b_y.x + a_y.w * b_y.y + a_z.w * b_y.z + a_w.w * b_y.w,
+
+		    // Z column
+		    a_x.x * b_z.x + a_y.x * b_z.y + a_z.x * b_z.z + a_w.x * b_z.w,
+		    a_x.y * b_z.x + a_y.y * b_z.y + a_z.y * b_z.z + a_w.y * b_z.w,
+		    a_x.z * b_z.x + a_y.z * b_z.y + a_z.z * b_z.z + a_w.z * b_z.w,
+		    a_x.w * b_z.x + a_y.w * b_z.y + a_z.w * b_z.z + a_w.w * b_z.w,
+
+		    // W column
+		    a_x.x * b_w.x + a_y.x * b_w.y + a_z.x * b_w.z + a_w.x * b_w.w,
+		    a_x.y * b_w.x + a_y.y * b_w.y + a_z.y * b_w.z + a_w.y * b_w.w,
+		    a_x.z * b_w.x + a_y.z * b_w.y + a_z.z * b_w.z + a_w.z * b_w.w,
+		    a_x.w * b_w.x + a_y.w * b_w.y + a_z.w * b_w.z + a_w.w * b_w.w
+		);
+	}
+
 	[[nodiscard]] constexpr auto operator[](size_t idx) -> Column_T &
 	{
+		debug_check(idx < num_elements, "mat4 out of bound access: {}", idx);
 		return values[idx];
 	}
 
@@ -67,51 +106,80 @@ struct mat4_impl
 		return values[idx];
 	}
 
-	[[nodiscard]] constexpr auto operator*(const mat4_impl<T> &other) const -> mat4_impl<T>
+	[[nodiscard]] static constexpr auto transpose(const mat4_impl<T> &mat) -> mat4_impl<T>
 	{
-		return mat4_impl<T> {};
+		const auto &[x, y, z, w] = mat.values;
+		return mat4_impl<T> {
+			x.x, y.x, z.x, w.x, x.y, y.y, z.y, w.y, x.z, y.z, z.z, w.z, x.w, y.w, z.w, w.w,
+		};
 	}
 
-	[[nodiscard]] constexpr auto operator*(const vec4_impl<T> &other) const -> vec4_impl<T>
+	[[nodiscard]] static constexpr auto translate(const vec3_impl<T> &vec) -> mat4_impl<T>
 	{
-		return vec4_impl<T> {};
+		return mat4_impl<T>(
+		    T { 1 },
+		    T { 0 },
+		    T { 0 },
+		    T { 0 },
+
+		    T { 0 },
+		    T { 1 },
+		    T { 0 },
+		    T { 0 },
+
+		    T { 0 },
+		    T { 0 },
+		    T { 1 },
+		    T { 0 },
+
+		    vec.x,
+		    vec.y,
+		    vec.z,
+		    T { 1 }
+		);
+	}
+
+	[[nodiscard]] static constexpr auto scale(const vec3_impl<T> &vec) -> mat4_impl<T>
+	{
+		return mat4_impl<T>(
+		    vec.x,
+		    T { 0 },
+		    T { 0 },
+		    T { 0 },
+
+		    T { 0 },
+		    vec.y,
+		    T { 0 },
+		    T { 0 },
+
+		    T { 0 },
+		    T { 0 },
+		    vec.z,
+		    T { 0 },
+
+		    T { 0 },
+		    T { 0 },
+		    T { 0 },
+		    T { 1 }
+		);
 	}
 
 	std::array<Column_T, 4u> values;
 };
 
-/** @todo(Light): Implement */
-template<typename T>
-[[nodiscard]] auto translate(const vec3_impl<T> &value) -> mat4_impl<T>
-{
-	return mat4_impl<T> {};
-}
-
-/** @todo(Light): Implement */
-template<typename T>
-[[nodiscard]] auto rotate(f32 value, const vec3_impl<T> &xyz) -> mat4_impl<T>
-{
-	return mat4_impl<T> {};
-}
-
-/** @todo(Light): Implement */
-template<typename T>
-[[nodiscard]] auto scale(const vec3_impl<T> &value) -> mat4_impl<T>
-{
-	return mat4_impl<T> {};
-}
-
-/** @todo(Light): Implement */
-template<typename T>
-[[nodiscard]] auto inverse(const mat4_impl<T> &value) -> mat4_impl<T>
-{
-	return mat4_impl<T> {};
-}
-
 using mat4 = mat4_impl<f32>;
 
-using imat4 = mat4_impl<i32>;
+using mat4_f32 = mat4;
+using mat4_f64 = mat4_impl<f64>;
 
-using umat4 = mat4_impl<u32>;
+using mat4_i8 = mat4_impl<i8>;
+using mat4_i16 = mat4_impl<i16>;
+using mat4_i32 = mat4_impl<i32>;
+using mat4_i64 = mat4_impl<i64>;
+
+using mat4_u8 = mat4_impl<u8>;
+using mat4_u16 = mat4_impl<u16>;
+using mat4_u32 = mat4_impl<u32>;
+using mat4_u64 = mat4_impl<u64>;
 
 } // namespace lt::math
diff --git a/modules/math/mat4.test.cpp b/modules/math/mat4.test.cpp
new file mode 100644
index 0000000..7809214
--- /dev/null
+++ b/modules/math/mat4.test.cpp
@@ -0,0 +1,413 @@
+import test;
+import math.vec3;
+import math.mat4;
+
+using vec3 = ::lt::math::vec3;
+using mat4 = ::lt::math::mat4;
+
+Suite static_tests = "mat4_static_checks"_suite = [] {
+	constexpr auto num_elements = lt::math::mat4::num_elements;
+
+	static_assert(num_elements == 4u * 4u);
+	static_assert(std::is_same_v<lt::math::mat4, lt::math::mat4_f32>);
+
+	static_assert(sizeof(lt::math::mat4_f32) == sizeof(f32) * num_elements);
+	static_assert(sizeof(lt::math::mat4_f64) == sizeof(f64) * num_elements);
+
+	static_assert(sizeof(lt::math::mat4_i8) == sizeof(i8) * num_elements);
+	static_assert(sizeof(lt::math::mat4_i16) == sizeof(i16) * num_elements);
+	static_assert(sizeof(lt::math::mat4_i32) == sizeof(i32) * num_elements);
+	static_assert(sizeof(lt::math::mat4_i64) == sizeof(i64) * num_elements);
+
+	static_assert(sizeof(lt::math::mat4_u8) == sizeof(u8) * num_elements);
+	static_assert(sizeof(lt::math::mat4_u16) == sizeof(u16) * num_elements);
+	static_assert(sizeof(lt::math::mat4_u32) == sizeof(u32) * num_elements);
+	static_assert(sizeof(lt::math::mat4_u64) == sizeof(u64) * num_elements);
+};
+
+Suite raii = "mat4_raii"_suite = [] {
+	Case { "happy paths" } = [] {
+		ignore = mat4 {};
+		ignore = mat4 { 1.0 };
+		ignore = mat4 {
+			1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
+		};
+		ignore = mat4 {
+			mat4::Column_T { 1.0, 2.0, 3.0, 4.0 },
+			mat4::Column_T { 5.0, 6.0, 7.0, 8.0 },
+			mat4::Column_T { 9.0, 10.0, 11.0, 12.0 },
+			mat4::Column_T { 13.0, 14.0, 15.0, 16.0 },
+		};
+	};
+
+	Case { "unhappy paths" } = [] {
+	};
+
+	Case { "many" } = [] {
+		for (auto idx : std::views::iota(0, 1'000'000))
+		{
+			ignore = idx;
+			ignore = mat4 {};
+			ignore = mat4 { 1.0 };
+			ignore = mat4 {
+				1.0, 2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,
+				9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
+			};
+			ignore = mat4 {
+				mat4::Column_T { 1.0, 2.0, 3.0, 4.0 },
+				mat4::Column_T { 5.0, 6.0, 7.0, 8.0 },
+				mat4::Column_T { 9.0, 10.0, 11.0, 12.0 },
+				mat4::Column_T { 13.0, 14.0, 15.0, 16.0 },
+			};
+		}
+	};
+
+	Case { "post default construct has correct state" } = [] {
+		const auto [x, y, z, w] = mat4 {}.values;
+
+		expect_eq(x[0], mat4::Underlying_T {});
+		expect_eq(x[1], mat4::Underlying_T {});
+		expect_eq(x[2], mat4::Underlying_T {});
+		expect_eq(x[3], mat4::Underlying_T {});
+
+		expect_eq(y[0], mat4::Underlying_T {});
+		expect_eq(y[1], mat4::Underlying_T {});
+		expect_eq(y[2], mat4::Underlying_T {});
+		expect_eq(y[3], mat4::Underlying_T {});
+
+		expect_eq(z[0], mat4::Underlying_T {});
+		expect_eq(z[1], mat4::Underlying_T {});
+		expect_eq(z[2], mat4::Underlying_T {});
+		expect_eq(z[3], mat4::Underlying_T {});
+
+		expect_eq(w[0], mat4::Underlying_T {});
+		expect_eq(w[1], mat4::Underlying_T {});
+		expect_eq(w[2], mat4::Underlying_T {});
+		expect_eq(w[3], mat4::Underlying_T {});
+	};
+
+	Case { "post scalar construct has correct state" } = [] {
+		const auto [x, y, z, w] = mat4 { 69.0 }.values;
+
+		expect_eq(x[0], mat4::Underlying_T { 69.0 });
+		expect_eq(x[1], mat4::Underlying_T { 69.0 });
+		expect_eq(x[2], mat4::Underlying_T { 69.0 });
+		expect_eq(x[3], mat4::Underlying_T { 69.0 });
+
+		expect_eq(y[0], mat4::Underlying_T { 69.0 });
+		expect_eq(y[1], mat4::Underlying_T { 69.0 });
+		expect_eq(y[2], mat4::Underlying_T { 69.0 });
+		expect_eq(y[3], mat4::Underlying_T { 69.0 });
+
+		expect_eq(z[0], mat4::Underlying_T { 69.0 });
+		expect_eq(z[1], mat4::Underlying_T { 69.0 });
+		expect_eq(z[2], mat4::Underlying_T { 69.0 });
+		expect_eq(z[3], mat4::Underlying_T { 69.0 });
+
+		expect_eq(w[0], mat4::Underlying_T { 69.0 });
+		expect_eq(w[1], mat4::Underlying_T { 69.0 });
+		expect_eq(w[2], mat4::Underlying_T { 69.0 });
+		expect_eq(w[3], mat4::Underlying_T { 69.0 });
+	};
+
+	Case { "post construct with all values has correct state" } = [] {
+		const auto [x, y, z, w] = mat4 {
+			1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
+		}.values;
+
+		expect_eq(x[0], mat4::Underlying_T { 1.0 });
+		expect_eq(x[1], mat4::Underlying_T { 2.0 });
+		expect_eq(x[2], mat4::Underlying_T { 3.0 });
+		expect_eq(x[3], mat4::Underlying_T { 4.0 });
+
+		expect_eq(y[0], mat4::Underlying_T { 5.0 });
+		expect_eq(y[1], mat4::Underlying_T { 6.0 });
+		expect_eq(y[2], mat4::Underlying_T { 7.0 });
+		expect_eq(y[3], mat4::Underlying_T { 8.0 });
+
+		expect_eq(z[0], mat4::Underlying_T { 9.0 });
+		expect_eq(z[1], mat4::Underlying_T { 10.0 });
+		expect_eq(z[2], mat4::Underlying_T { 11.0 });
+		expect_eq(z[3], mat4::Underlying_T { 12.0 });
+
+		expect_eq(w[0], mat4::Underlying_T { 13.0 });
+		expect_eq(w[1], mat4::Underlying_T { 14.0 });
+		expect_eq(w[2], mat4::Underlying_T { 15.0 });
+		expect_eq(w[3], mat4::Underlying_T { 16.0 });
+	};
+
+	Case { "post construct with columns has correct state" } = [] {
+		const auto [x, y, z, w] = mat4 {
+			mat4::Column_T { 1.0, 2.0, 3.0, 4.0 },
+			mat4::Column_T { 5.0, 6.0, 7.0, 8.0 },
+			mat4::Column_T { 9.0, 10.0, 11.0, 12.0 },
+			mat4::Column_T { 13.0, 14.0, 15.0, 16.0 },
+		}.values;
+
+		expect_eq(x[0], mat4::Underlying_T { 1.0 });
+		expect_eq(x[1], mat4::Underlying_T { 2.0 });
+		expect_eq(x[2], mat4::Underlying_T { 3.0 });
+		expect_eq(x[3], mat4::Underlying_T { 4.0 });
+
+		expect_eq(y[0], mat4::Underlying_T { 5.0 });
+		expect_eq(y[1], mat4::Underlying_T { 6.0 });
+		expect_eq(y[2], mat4::Underlying_T { 7.0 });
+		expect_eq(y[3], mat4::Underlying_T { 8.0 });
+
+		expect_eq(z[0], mat4::Underlying_T { 9.0 });
+		expect_eq(z[1], mat4::Underlying_T { 10.0 });
+		expect_eq(z[2], mat4::Underlying_T { 11.0 });
+		expect_eq(z[3], mat4::Underlying_T { 12.0 });
+
+		expect_eq(w[0], mat4::Underlying_T { 13.0 });
+		expect_eq(w[1], mat4::Underlying_T { 14.0 });
+		expect_eq(w[2], mat4::Underlying_T { 15.0 });
+		expect_eq(w[3], mat4::Underlying_T { 16.0 });
+	};
+
+	Case { "post construct identity matrix has correct state" } = [] {
+		const auto [x, y, z, w] = mat4::identity().values;
+
+		expect_eq(x[0], mat4::Underlying_T { 1 });
+		expect_eq(x[1], mat4::Underlying_T {});
+		expect_eq(x[2], mat4::Underlying_T {});
+		expect_eq(x[3], mat4::Underlying_T {});
+
+		expect_eq(y[0], mat4::Underlying_T {});
+		expect_eq(y[1], mat4::Underlying_T { 1 });
+		expect_eq(y[2], mat4::Underlying_T {});
+		expect_eq(y[3], mat4::Underlying_T {});
+
+		expect_eq(z[0], mat4::Underlying_T {});
+		expect_eq(z[1], mat4::Underlying_T {});
+		expect_eq(z[2], mat4::Underlying_T { 1 });
+		expect_eq(z[3], mat4::Underlying_T {});
+
+		expect_eq(w[0], mat4::Underlying_T {});
+		expect_eq(w[1], mat4::Underlying_T {});
+		expect_eq(w[2], mat4::Underlying_T {});
+		expect_eq(w[3], mat4::Underlying_T { 1 });
+	};
+};
+
+Suite arithmetic_operators = "mat4_arithmetic_operators"_suite = [] {
+	Case { "operator *" } = [] {
+		const auto lhs = mat4 {
+			mat4::Column_T { 1.0, 2.0, 3.0, 4.0 },
+			mat4::Column_T { 5.0, 6.0, 7.0, 8.0 },
+			mat4::Column_T { 9.0, 10.0, 11.0, 12.0 },
+			mat4::Column_T { 13.0, 14.0, 15.0, 16.0 },
+		};
+
+		const auto rhs = mat4 {
+			mat4::Column_T { 17.0, 18.0, 19.0, 20.0 },
+			mat4::Column_T { 21.0, 22.0, 23.0, 24.0 },
+			mat4::Column_T { 25.0, 26.0, 27.0, 28.0 },
+			mat4::Column_T { 29.0, 30.0, 31.0, 32.0 },
+		};
+
+		const auto [x, y, z, w] = (lhs * rhs).values;
+
+		expect_eq(x[0], 538.0);
+		expect_eq(x[1], 612.0);
+		expect_eq(x[2], 686.0);
+		expect_eq(x[3], 760.0);
+
+		expect_eq(y[0], 650.0);
+		expect_eq(y[1], 740.0);
+		expect_eq(y[2], 830.0);
+		expect_eq(y[3], 920.0);
+
+		expect_eq(z[0], 762.0);
+		expect_eq(z[1], 868.0);
+		expect_eq(z[2], 974.0);
+		expect_eq(z[3], 1080.0);
+
+		expect_eq(w[0], 874.0);
+		expect_eq(w[1], 996.0);
+		expect_eq(w[2], 1118.0);
+		expect_eq(w[3], 1240.0);
+	};
+};
+
+Suite access_operators = "mat4_access_operators"_suite = [] {
+	Case { "operator []" } = [] {
+		auto mat = mat4 {
+			1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
+		};
+
+		expect_eq(mat[0][0], 1.0);
+		expect_eq(mat[0][1], 2.0);
+		expect_eq(mat[0][2], 3.0);
+		expect_eq(mat[0][3], 4.0);
+
+		expect_eq(mat[1][0], 5.0);
+		expect_eq(mat[1][1], 6.0);
+		expect_eq(mat[1][2], 7.0);
+		expect_eq(mat[1][3], 8.0);
+
+		expect_eq(mat[2][0], 9.0);
+		expect_eq(mat[2][1], 10.0);
+		expect_eq(mat[2][2], 11.0);
+		expect_eq(mat[2][3], 12.0);
+
+		expect_eq(mat[3][0], 13.0);
+		expect_eq(mat[3][1], 14.0);
+		expect_eq(mat[3][2], 15.0);
+		expect_eq(mat[3][3], 16.0);
+	};
+
+	Case { "operator [] const" } = [] {
+		const auto mat = mat4 {
+			1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
+		};
+
+		expect_eq(mat[0][0], 1.0);
+		expect_eq(mat[0][1], 2.0);
+		expect_eq(mat[0][2], 3.0);
+		expect_eq(mat[0][3], 4.0);
+
+		expect_eq(mat[1][0], 5.0);
+		expect_eq(mat[1][1], 6.0);
+		expect_eq(mat[1][2], 7.0);
+		expect_eq(mat[1][3], 8.0);
+
+		expect_eq(mat[2][0], 9.0);
+		expect_eq(mat[2][1], 10.0);
+		expect_eq(mat[2][2], 11.0);
+		expect_eq(mat[2][3], 12.0);
+
+		expect_eq(mat[3][0], 13.0);
+		expect_eq(mat[3][1], 14.0);
+		expect_eq(mat[3][2], 15.0);
+		expect_eq(mat[3][3], 16.0);
+	};
+};
+
+Suite transformations = "mat4_transformations"_suite = [] {
+	Case { "translate" } = [] {
+		const auto &[x, y, z, w] = mat4::translate(vec3 { 1, 2, 3 }).values;
+
+		// identity basis
+		expect_eq(x[0], 1);
+		expect_eq(x[1], 0);
+		expect_eq(x[2], 0);
+		expect_eq(x[3], 0);
+		expect_eq(y[0], 0);
+		expect_eq(y[1], 1);
+		expect_eq(y[2], 0);
+		expect_eq(y[3], 0);
+		expect_eq(z[0], 0);
+		expect_eq(z[1], 0);
+		expect_eq(z[2], 1);
+		expect_eq(z[3], 0);
+
+		// translation column
+		expect_eq(w[0], 1);
+		expect_eq(w[1], 2);
+		expect_eq(w[2], 3);
+		expect_eq(w[3], 1);
+	};
+
+	Case { "scale" } = [] {
+		const auto [x, y, z, w] = mat4::scale(vec3 { 2, 3, 4 }).values;
+
+		expect_eq(x[0], 2);
+		expect_eq(x[1], 0);
+		expect_eq(x[2], 0);
+		expect_eq(x[3], 0);
+		expect_eq(y[0], 0);
+		expect_eq(y[1], 3);
+		expect_eq(y[2], 0);
+		expect_eq(y[3], 0);
+		expect_eq(z[0], 0);
+		expect_eq(z[1], 0);
+		expect_eq(z[2], 4);
+		expect_eq(z[3], 0);
+		expect_eq(w[0], 0);
+		expect_eq(w[1], 0);
+		expect_eq(w[2], 0);
+		expect_eq(w[3], 1);
+	};
+
+	Case { "scale -> translate" } = [] {
+		const auto scale = mat4::scale(vec3 { 2, 2, 2 });
+		const auto translate = mat4::translate(vec3 { 1, 2, 3 });
+		const auto [x, y, z, w] = (scale * translate).values;
+
+		// scaled basis
+		expect_eq(x[0], 2);
+		expect_eq(x[1], 0);
+		expect_eq(x[2], 0);
+		expect_eq(x[3], 0);
+		expect_eq(y[0], 0);
+		expect_eq(y[1], 2);
+		expect_eq(y[2], 0);
+		expect_eq(y[3], 0);
+		expect_eq(z[0], 0);
+		expect_eq(z[1], 0);
+		expect_eq(z[2], 2);
+		expect_eq(z[3], 0);
+
+		// translation is scaled (local-space translation)
+		expect_eq(w[0], 2); // 1 * 2
+		expect_eq(w[1], 4); // 2 * 2
+		expect_eq(w[2], 6); // 3 * 2
+		expect_eq(w[3], 1);
+	};
+
+	Case { "transpose" } = [] {
+		const auto mat = mat4 {
+			mat4::Column_T { 1, 2, 3, 4 },
+			mat4::Column_T { 5, 6, 7, 8 },
+			mat4::Column_T { 9, 10, 11, 12 },
+			mat4::Column_T { 13, 14, 15, 16 },
+		};
+		const auto [x, y, z, w] = mat4::transpose(mat).values;
+
+		// rows become columns
+		expect_eq(x[0], 1);
+		expect_eq(x[1], 5);
+		expect_eq(x[2], 9);
+		expect_eq(x[3], 13);
+		expect_eq(y[0], 2);
+		expect_eq(y[1], 6);
+		expect_eq(y[2], 10);
+		expect_eq(y[3], 14);
+		expect_eq(z[0], 3);
+		expect_eq(z[1], 7);
+		expect_eq(z[2], 11);
+		expect_eq(z[3], 15);
+		expect_eq(w[0], 4);
+		expect_eq(w[1], 8);
+		expect_eq(w[2], 12);
+		expect_eq(w[3], 16);
+	};
+
+	Case { "transpose twice" } = [] {
+		const auto mat = mat4 {
+			mat4::Column_T { 1, 2, 3, 4 },
+			mat4::Column_T { 5, 6, 7, 8 },
+			mat4::Column_T { 9, 10, 11, 12 },
+			mat4::Column_T { 13, 14, 15, 16 },
+		};
+		const auto [x, y, z, w] = mat4::transpose(mat4::transpose(mat)).values;
+
+		expect_eq(x[0], 1);
+		expect_eq(x[1], 2);
+		expect_eq(x[2], 3);
+		expect_eq(x[3], 4);
+		expect_eq(y[0], 5);
+		expect_eq(y[1], 6);
+		expect_eq(y[2], 7);
+		expect_eq(y[3], 8);
+		expect_eq(z[0], 9);
+		expect_eq(z[1], 10);
+		expect_eq(z[2], 11);
+		expect_eq(z[3], 12);
+		expect_eq(w[0], 13);
+		expect_eq(w[1], 14);
+		expect_eq(w[2], 15);
+		expect_eq(w[3], 16);
+	};
+};
diff --git a/modules/math/trig.cppm b/modules/math/trig.cppm
index 33571b1..93607d9 100644
--- a/modules/math/trig.cppm
+++ b/modules/math/trig.cppm
@@ -4,22 +4,22 @@ import preliminary;
 
 export namespace lt::math {
 
-[[nodiscard]] constexpr auto radians(f32 degrees) -> f32
+[[nodiscard]] constexpr auto to_radians(f32 degrees) -> f32
 {
 	return degrees * 0.01745329251994329576923690768489f;
 }
 
-[[nodiscard]] constexpr auto radians(f64 degrees) -> f64
+[[nodiscard]] constexpr auto to_radians(f64 degrees) -> f64
 {
 	return degrees * 0.01745329251994329576923690768489;
 }
 
-[[nodiscard]] constexpr auto degrees(f32 radians) -> f32
+[[nodiscard]] constexpr auto to_degrees(f32 radians) -> f32
 {
 	return radians * 57.295779513082320876798154814105f;
 }
 
-[[nodiscard]] constexpr auto degrees(f64 radians) -> f64
+[[nodiscard]] constexpr auto to_degrees(f64 radians) -> f64
 {
 	return radians * 57.295779513082320876798154814105;
 }
diff --git a/modules/math/trig.test.cpp b/modules/math/trig.test.cpp
new file mode 100644
index 0000000..a30498e
--- /dev/null
+++ b/modules/math/trig.test.cpp
@@ -0,0 +1,43 @@
+import test;
+import math.trig;
+
+Suite conversions = "trig_conversions"_suite = [] {
+	using ::lt::math::to_degrees;
+	using ::lt::math::to_radians;
+
+	Case { "to_radians <f32>" } = [] {
+		expect_eq(to_radians(f32 { 0.0f }), f32 { 0.0f });
+		expect_eq(to_radians(f32 { 90.0f }), f32 { 1.5707963267948966f });
+		expect_eq(to_radians(f32 { 180.0f }), f32 { 3.1415926535897932f });
+		expect_eq(to_radians(f32 { 360.0f }), f32 { 6.2831853071795864f });
+	};
+
+	Case { "to_radians <f64>" } = [] {
+		expect_eq(to_radians(f64 { 0.0 }), f64 { 0.0 });
+		expect_eq(to_radians(f64 { 90.0 }), f64 { 1.5707963267948966 });
+		expect_eq(to_radians(f64 { 180.0 }), f64 { 3.1415926535897932 });
+		expect_eq(to_radians(f64 { 360.0 }), f64 { 6.2831853071795864 });
+	};
+
+	Case { "to_degrees <f32>" } = [] {
+		expect_eq(to_degrees(f32 { 0.0f }), f32 { 0.0f });
+		expect_eq(to_degrees(f32 { 1.5707963267948966f }), f32 { 90.0f });
+		expect_eq(to_degrees(f32 { 3.1415926535897932f }), f32 { 180.0f });
+		expect_eq(to_degrees(f32 { 6.2831853071795864f }), f32 { 360.0f });
+	};
+
+	Case { "to_degrees <f64>" } = [] {
+		expect_eq(to_degrees(f64 { 0.0 }), f64 { 0.0 });
+		expect_eq(to_degrees(f64 { 1.5707963267948966 }), f64 { 90.0 });
+		expect_eq(to_degrees(f64 { 3.1415926535897932 }), f64 { 180.0 });
+		expect_eq(to_degrees(f64 { 6.2831853071795864 }), f64 { 360.0 });
+	};
+
+	Case { "to_degrees -> to_radians -> to_degrees <f32>" } = [] {
+		expect_eq(to_degrees(to_radians(f32 { 45.0f })), f32 { 45.0f });
+	};
+
+	Case { "to_degrees -> to_radians -> to_degrees <f64>" } = [] {
+		expect_eq(to_degrees(to_radians(f64 { 45.0 })), f64 { 45.0 });
+	};
+};
diff --git a/modules/math/vec2.cppm b/modules/math/vec2.cppm
index bb59e33..f8e2dd4 100644
--- a/modules/math/vec2.cppm
+++ b/modules/math/vec2.cppm
@@ -8,6 +8,8 @@ template<typename T = f32>
     requires(std::is_arithmetic_v<T>)
 struct vec2_impl
 {
+	using Underlying_T = T;
+
 	static constexpr auto num_elements = 2u;
 
 	constexpr vec2_impl(): x(), y()
@@ -66,33 +68,35 @@ struct vec2_impl
 
 	[[nodiscard]] constexpr auto operator[](u8 idx) -> T &
 	{
-		debug_check(idx <= num_elements, "vec2 out of bound: {}", idx);
+		debug_check(idx < num_elements, "vec2 out of bound access: {}", idx);
 		return ((T *)this)[idx];
 	}
 
 	[[nodiscard]] constexpr auto operator[](u8 idx) const -> const T &
 	{
-		debug_check(idx < num_elements, "vec2 out of bound: {}", idx);
+		debug_check(idx < num_elements, "vec2 out of bound access: {}", idx);
 		return ((T *)this)[idx];
 	}
 
-	friend auto operator<<(std::ostream &stream, vec2_impl<T> value) -> std::ostream &
-	{
-		stream << value.x << ", " << value.y;
-		return stream;
-	}
-
 	T x;
 
 	T y;
 };
 
-
 using vec2 = vec2_impl<f32>;
 
-using ivec2 = vec2_impl<i32>;
+using vec2_f32 = vec2;
+using vec2_f64 = vec2_impl<f64>;
 
-using uvec2 = vec2_impl<u32>;
+using vec2_i8 = vec2_impl<i8>;
+using vec2_i16 = vec2_impl<i16>;
+using vec2_i32 = vec2_impl<i32>;
+using vec2_i64 = vec2_impl<i64>;
+
+using vec2_u8 = vec2_impl<u8>;
+using vec2_u16 = vec2_impl<u16>;
+using vec2_u32 = vec2_impl<u32>;
+using vec2_u64 = vec2_impl<u64>;
 
 } // namespace lt::math
 
diff --git a/modules/math/vec2.test.cpp b/modules/math/vec2.test.cpp
index 08a4de0..9364747 100644
--- a/modules/math/vec2.test.cpp
+++ b/modules/math/vec2.test.cpp
@@ -1,7 +1,130 @@
 import test;
 import math.vec2;
 
-Suite raii = "raii"_suite = [] {
-	Case { "happy path" } = [] {
+using vec2 = ::lt::math::vec2;
+using ivec2 = ::lt::math::vec2_i32;
+
+Suite static_tests = "vec3_static_checks"_suite = [] {
+	constexpr auto num_elements = lt::math::vec2::num_elements;
+
+	static_assert(num_elements == 2u);
+	static_assert(std::is_same_v<lt::math::vec2, lt::math::vec2_f32>);
+
+	static_assert(sizeof(lt::math::vec2_f32) == sizeof(f32) * num_elements);
+	static_assert(sizeof(lt::math::vec2_f64) == sizeof(f64) * num_elements);
+
+	static_assert(sizeof(lt::math::vec2_i8) == sizeof(i8) * num_elements);
+	static_assert(sizeof(lt::math::vec2_i16) == sizeof(i16) * num_elements);
+	static_assert(sizeof(lt::math::vec2_i32) == sizeof(i32) * num_elements);
+	static_assert(sizeof(lt::math::vec2_i64) == sizeof(i64) * num_elements);
+
+	static_assert(sizeof(lt::math::vec2_u8) == sizeof(u8) * num_elements);
+	static_assert(sizeof(lt::math::vec2_u16) == sizeof(u16) * num_elements);
+	static_assert(sizeof(lt::math::vec2_u32) == sizeof(u32) * num_elements);
+	static_assert(sizeof(lt::math::vec2_u64) == sizeof(u64) * num_elements);
+};
+
+Suite raii = "vec2_raii"_suite = [] {
+	Case { "happy paths" } = [] {
+		ignore = vec2 {};
+		ignore = vec2 { 2.0 };
+		ignore = vec2 { 2.0, 4.0 };
+	};
+
+	Case { "unhappy paths" } = [] {
+	};
+
+	Case { "many" } = [] {
+		for (auto idx : std::views::iota(0, 1'000'000))
+		{
+			ignore = idx;
+			ignore = vec2 {};
+			ignore = vec2 { 2.0 };
+			ignore = vec2 { 2.0, 4.0 };
+		}
+	};
+
+	Case { "post default construct has correct state" } = [] {
+		const auto vec = vec2 {};
+		expect_eq(vec.x, 0.0);
+		expect_eq(vec.y, 0.0);
+	};
+
+	Case { "post scalar construct has correct state" } = [] {
+		const auto vec = vec2 { 2.0 };
+		expect_eq(vec.x, 2.0);
+		expect_eq(vec.y, 2.0);
+	};
+
+	Case { "post construct with x,y has correct state" } = [] {
+		const auto vec = vec2 { 2.0, 3.0 };
+		expect_eq(vec.x, 2.0);
+		expect_eq(vec.y, 3.0);
+	};
+};
+
+Suite arithmetic_operators = "vec2_operators"_suite = [] {
+	Case { "operator ==" } = [] {
+		const auto lhs = vec2 { 1.0, 2.0 };
+
+		expect_false(lhs == vec2 { {}, 2.0 });
+		expect_false(lhs == vec2 { 1.0, {} });
+		expect_true(lhs == vec2 { 1.0, 2.0 });
+	};
+
+	Case { "operator !=" } = [] {
+		const auto lhs = vec2 { 1.0, 2.0 };
+
+		expect_true(lhs != vec2 { {}, 2.0 });
+		expect_true(lhs != vec2 { 1.0, {} });
+		expect_false(lhs != vec2 { 1.0, 2.0 });
+	};
+
+	Case { "operator +" } = [] {
+		const auto lhs = vec2 { 2.0, 3.0 };
+		const auto rhs = vec2 { 4.0, 5.0 };
+		expect_eq(lhs + rhs, vec2 { 6.0, 8.0 });
+	};
+
+	Case { "operator -" } = [] {
+		const auto lhs = vec2 { 2.0, 3.0 };
+		const auto rhs = vec2 { 4.0, 6.0 };
+		expect_eq(lhs - rhs, vec2 { -2.0, -3.0 });
+	};
+
+	Case { "operator *" } = [] {
+		const auto lhs = vec2 { 2.0, 3.0 };
+		const auto rhs = vec2 { 10.0, 20.0 };
+		expect_eq(lhs * rhs, vec2 { 20.0, 60.0 });
+	};
+
+	Case { "operator /" } = [] {
+		const auto lhs = vec2 { 10.0, 20.0 };
+		const auto rhs = vec2 { 2.0, 20.0 };
+		expect_eq(lhs / rhs, vec2 { 5.0, 1.0 });
+	};
+
+	Case { "operator []" } = [] {
+		auto vec = vec2 { 0.0, 1.0 };
+		expect_eq(vec[0], 0.0);
+		expect_eq(vec[1], 1.0);
+	};
+
+	Case { "operator [] const" } = [] {
+		const auto vec = vec2 { 0.0, 1.0 };
+		expect_eq(vec[0], 0.0);
+		expect_eq(vec[1], 1.0);
+	};
+};
+
+Suite utilities = "vec2_utilities"_suite = [] {
+	Case { "std::format float" } = [] {
+		auto str = std::format("{}", vec2 { 10.0000f, 30.0005f });
+		expect_eq(str, "10, 30.0005");
+	};
+
+	Case { "std::format int" } = [] {
+		auto str = std::format("{}", ivec2 { 10, 30 });
+		expect_eq(str, "10, 30");
 	};
 };
diff --git a/modules/math/vec3.cppm b/modules/math/vec3.cppm
index db13702..2038214 100644
--- a/modules/math/vec3.cppm
+++ b/modules/math/vec3.cppm
@@ -9,6 +9,8 @@ template<typename T = f32>
     requires(std::is_arithmetic_v<T>)
 struct vec3_impl
 {
+	using Underlying_T = T;
+
 	static constexpr auto num_elements = 3u;
 
 	constexpr vec3_impl(): x(), y(), z()
@@ -27,7 +29,7 @@ struct vec3_impl
 	{
 	}
 
-	constexpr vec3_impl(T x, vec2_impl<T> yz): x(x), y(yz.y), z(yz.z)
+	constexpr vec3_impl(T x, vec2_impl<T> yz): x(x), y(yz.x), z(yz.y)
 	{
 	}
 
@@ -79,13 +81,13 @@ struct vec3_impl
 
 	[[nodiscard]] constexpr auto operator[](u8 idx) -> T &
 	{
-		debug_check(idx <= num_elements, "vec3 out of bound: {}", idx);
+		debug_check(idx < num_elements, "vec3 out of bound access: {}", idx);
 		return ((T *)this)[idx];
 	}
 
 	[[nodiscard]] constexpr auto operator[](u8 idx) const -> const T &
 	{
-		debug_check(idx < num_elements, "vec3 out of bound: {}", idx);
+		debug_check(idx < num_elements, "vec3 out of bound access: {}", idx);
 		return ((T *)this)[idx];
 	}
 
@@ -104,9 +106,18 @@ struct vec3_impl
 
 using vec3 = vec3_impl<f32>;
 
-using ivec3 = vec3_impl<i32>;
+using vec3_f32 = vec3;
+using vec3_f64 = vec3_impl<f64>;
 
-using uvec3 = vec3_impl<u32>;
+using vec3_i8 = vec3_impl<i8>;
+using vec3_i16 = vec3_impl<i16>;
+using vec3_i32 = vec3_impl<i32>;
+using vec3_i64 = vec3_impl<i64>;
+
+using vec3_u8 = vec3_impl<u8>;
+using vec3_u16 = vec3_impl<u16>;
+using vec3_u32 = vec3_impl<u32>;
+using vec3_u64 = vec3_impl<u64>;
 
 } // namespace lt::math
 
diff --git a/modules/math/vec3.test.cpp b/modules/math/vec3.test.cpp
new file mode 100644
index 0000000..ca7782a
--- /dev/null
+++ b/modules/math/vec3.test.cpp
@@ -0,0 +1,157 @@
+import test;
+import math.vec2;
+import math.vec3;
+
+using vec2 = ::lt::math::vec2;
+using vec3 = ::lt::math::vec3;
+using ivec3 = ::lt::math::vec3_i32;
+
+Suite static_tests = "vec3_static_checks"_suite = [] {
+	constexpr auto num_elements = lt::math::vec3::num_elements;
+
+	static_assert(num_elements == 3u);
+	static_assert(std::is_same_v<lt::math::vec3, lt::math::vec3_f32>);
+
+	static_assert(sizeof(lt::math::vec3_f32) == sizeof(f32) * num_elements);
+	static_assert(sizeof(lt::math::vec3_f64) == sizeof(f64) * num_elements);
+
+	static_assert(sizeof(lt::math::vec3_i8) == sizeof(i8) * num_elements);
+	static_assert(sizeof(lt::math::vec3_i16) == sizeof(i16) * num_elements);
+	static_assert(sizeof(lt::math::vec3_i32) == sizeof(i32) * num_elements);
+	static_assert(sizeof(lt::math::vec3_i64) == sizeof(i64) * num_elements);
+
+	static_assert(sizeof(lt::math::vec3_u8) == sizeof(u8) * num_elements);
+	static_assert(sizeof(lt::math::vec3_u16) == sizeof(u16) * num_elements);
+	static_assert(sizeof(lt::math::vec3_u32) == sizeof(u32) * num_elements);
+	static_assert(sizeof(lt::math::vec3_u64) == sizeof(u64) * num_elements);
+};
+
+Suite raii = "vec3_raii"_suite = [] {
+	Case { "happy paths" } = [] {
+		ignore = vec3 {};
+		ignore = vec3 { 2.0 };
+		ignore = vec3 { 2.0, 4.0, 6.0 };
+		ignore = vec3 { vec2 { 2.0, 4.0 }, 6.0 };
+		ignore = vec3 { 2.0, vec2 { 4.0, 6.0 } };
+	};
+
+	Case { "unhappy paths" } = [] {
+	};
+
+	Case { "many" } = [] {
+		for (auto idx : std::views::iota(0, 1'000'000))
+		{
+			ignore = idx;
+			ignore = vec3 {};
+			ignore = vec3 { 2.0 };
+			ignore = vec3 { 2.0, 4.0, 6.0 };
+			ignore = vec3 { vec2 { 2.0, 4.0 }, 6.0 };
+			ignore = vec3 { 2.0, vec2 { 4.0, 6.0 } };
+		}
+	};
+
+	Case { "post default construct has correct state" } = [] {
+		const auto vec = vec3 {};
+		expect_eq(vec.x, 0.0);
+		expect_eq(vec.y, 0.0);
+		expect_eq(vec.z, 0.0);
+	};
+
+	Case { "post scalar construct has correct state" } = [] {
+		const auto vec = vec3 { 2.0 };
+		expect_eq(vec.x, 2.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 2.0);
+	};
+
+	Case { "post construct with x,y,z has correct state" } = [] {
+		const auto vec = vec3 { 1.0, 2.0, 3.0 };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.y, 3.0);
+	};
+
+	Case { "post construct with xy,z has correct state" } = [] {
+		const auto vec = vec3 { vec2 { 1.0, 2.0 }, 3.0 };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+	};
+
+	Case { "post construct with x,yz has correct state" } = [] {
+		const auto vec = vec3 { 1.0, vec2 { 2.0, 3.0 } };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+	};
+};
+
+Suite arithmetic_operators = "vec3_operators"_suite = [] {
+	Case { "operator ==" } = [] {
+		const auto lhs = vec3 { 1.0, 2.0, 3.0 };
+
+		expect_false(lhs == vec3 { {}, 2.0, 3.0 });
+		expect_false(lhs == vec3 { 1.0, {}, 3.0 });
+		expect_false(lhs == vec3 { 1.0, 2.0, {} });
+		expect_true(lhs == vec3 { 1.0, 2.0, 3.0 });
+	};
+
+	Case { "operator !=" } = [] {
+		const auto lhs = vec3 { 1.0, 2.0, 3.0 };
+
+		expect_true(lhs != vec3 { {}, 2.0, 3.0 });
+		expect_true(lhs != vec3 { 1.0, {}, 3.0 });
+		expect_true(lhs != vec3 { 1.0, 2.0, {} });
+		expect_false(lhs != vec3 { 1.0, 2.0, 3.0 });
+	};
+
+	Case { "operator +" } = [] {
+		const auto lhs = vec3 { 1.0, 2.0, 3.0 };
+		const auto rhs = vec3 { 4.0, 5.0, 6.0 };
+		expect_eq(lhs + rhs, vec3 { 5.0, 7.0, 9.0 });
+	};
+
+	Case { "operator -" } = [] {
+		const auto lhs = vec3 { 1.0, 2.0, 3.0 };
+		const auto rhs = vec3 { 4.0, 5.0, 7.0 };
+		expect_eq(lhs - rhs, vec3 { -2.0, -3.0, -4.0 });
+	};
+
+	Case { "operator *" } = [] {
+		const auto lhs = vec3 { 1.0, 2.0, 3.0 };
+		const auto rhs = vec3 { 4.0, 5.0, 6.0 };
+		expect_eq(lhs * rhs, vec3 { 4.0, 10.0, 18.0 });
+	};
+
+	Case { "operator /" } = [] {
+		const auto lhs = vec3 { 4.0, 10.0, 30.0 };
+		const auto rhs = vec3 { 1.0, 2.0, 5.0 };
+		expect_eq(lhs / rhs, vec3 { 4.0, 5.0, 6.0 });
+	};
+
+	Case { "operator []" } = [] {
+		auto vec = vec3 { 0.0, 1.0, 2.0 };
+		expect_eq(vec[0], 0.0);
+		expect_eq(vec[1], 1.0);
+		expect_eq(vec[2], 2.0);
+	};
+
+	Case { "operator [] const" } = [] {
+		const auto vec = vec3 { 0.0, 1.0, 2.0 };
+		expect_eq(vec[0], 0.0);
+		expect_eq(vec[1], 1.0);
+		expect_eq(vec[2], 2.0);
+	};
+};
+
+Suite utilities = "vec3_utilities"_suite = [] {
+	Case { "std::format float" } = [] {
+		auto str = std::format("{}", vec3 { 10.0000f, 30.0005f, 40.00005f });
+		expect_eq(str, "10, 30.0005, 40.00005");
+	};
+
+	Case { "std::format int" } = [] {
+		auto str = std::format("{}", ivec3 { 10, 30, 3'000'000 });
+		expect_eq(str, "10, 30, 3000000");
+	};
+};
diff --git a/modules/math/vec4.cppm b/modules/math/vec4.cppm
index db5fdb6..e8a98ec 100644
--- a/modules/math/vec4.cppm
+++ b/modules/math/vec4.cppm
@@ -10,6 +10,8 @@ template<typename T = f32>
     requires(std::is_arithmetic_v<T>)
 struct vec4_impl
 {
+	using Underlying_T = T;
+
 	static constexpr auto num_elements = 4u;
 
 	constexpr vec4_impl(): x(), y(), z(), w()
@@ -28,11 +30,15 @@ struct vec4_impl
 	{
 	}
 
-	constexpr vec4_impl(T x, T y, vec2_impl<T> zw): x(x), y(y), z(zw.z), w(zw.w)
+	constexpr vec4_impl(T x, vec2_impl<T> yz, T w): x(x), y(yz.x), z(yz.y), w(w)
 	{
 	}
 
-	constexpr vec4_impl(vec2_impl<T> xy, vec2_impl<T> zw): x(xy.x), y(xy.y), z(zw.z), w(zw.w)
+	constexpr vec4_impl(T x, T y, vec2_impl<T> zw): x(x), y(y), z(zw.x), w(zw.y)
+	{
+	}
+
+	constexpr vec4_impl(vec2_impl<T> xy, vec2_impl<T> zw): x(xy.x), y(xy.y), z(zw.x), w(zw.y)
 	{
 	}
 
@@ -40,7 +46,7 @@ struct vec4_impl
 	{
 	}
 
-	constexpr vec4_impl(T x, vec3_impl<T> yzw): x(x), y(yzw.y), z(yzw.z), w(yzw.w)
+	constexpr vec4_impl(T x, vec3_impl<T> yzw): x(x), y(yzw.x), z(yzw.y), w(yzw.z)
 	{
 	}
 
@@ -96,13 +102,13 @@ struct vec4_impl
 
 	[[nodiscard]] constexpr auto operator[](u8 idx) -> T &
 	{
-		debug_check(idx <= num_elements, "vec4 out of bound: {}", idx);
+		debug_check(idx < num_elements, "vec4 out of bound access: {}", idx);
 		return ((T *)this)[idx];
 	}
 
 	[[nodiscard]] constexpr auto operator[](u8 idx) const -> const T &
 	{
-		debug_check(idx < num_elements, "vec4 out of bound: {}", idx);
+		debug_check(idx < num_elements, "vec4 out of bound access: {}", idx);
 		return ((T *)this)[idx];
 	}
 
@@ -123,9 +129,18 @@ struct vec4_impl
 
 using vec4 = vec4_impl<f32>;
 
-using ivec4 = vec4_impl<i32>;
+using vec4_f32 = vec4;
+using vec4_f64 = vec4_impl<f64>;
 
-using uvec4 = vec4_impl<u32>;
+using vec4_i8 = vec4_impl<i8>;
+using vec4_i16 = vec4_impl<i16>;
+using vec4_i32 = vec4_impl<i32>;
+using vec4_i64 = vec4_impl<i64>;
+
+using vec4_u8 = vec4_impl<u8>;
+using vec4_u16 = vec4_impl<u16>;
+using vec4_u32 = vec4_impl<u32>;
+using vec4_u64 = vec4_impl<u64>;
 
 } // namespace lt::math
 
diff --git a/modules/math/vec4.test.cpp b/modules/math/vec4.test.cpp
new file mode 100644
index 0000000..041c9d2
--- /dev/null
+++ b/modules/math/vec4.test.cpp
@@ -0,0 +1,215 @@
+import test;
+import math.vec2;
+import math.vec3;
+import math.vec4;
+import logger;
+
+using vec2 = ::lt::math::vec2;
+using vec3 = ::lt::math::vec3;
+using vec4 = ::lt::math::vec4;
+using ivec4 = ::lt::math::vec4_i32;
+
+Suite static_tests = "vec4_static_checks"_suite = [] {
+	constexpr auto num_elements = lt::math::vec4::num_elements;
+
+	static_assert(num_elements == 4u);
+	static_assert(std::is_same_v<lt::math::vec4, lt::math::vec4_f32>);
+
+	static_assert(sizeof(lt::math::vec4_f32) == sizeof(f32) * num_elements);
+	static_assert(sizeof(lt::math::vec4_f64) == sizeof(f64) * num_elements);
+
+	static_assert(sizeof(lt::math::vec4_i8) == sizeof(i8) * num_elements);
+	static_assert(sizeof(lt::math::vec4_i16) == sizeof(i16) * num_elements);
+	static_assert(sizeof(lt::math::vec4_i32) == sizeof(i32) * num_elements);
+	static_assert(sizeof(lt::math::vec4_i64) == sizeof(i64) * num_elements);
+
+	static_assert(sizeof(lt::math::vec4_u8) == sizeof(u8) * num_elements);
+	static_assert(sizeof(lt::math::vec4_u16) == sizeof(u16) * num_elements);
+	static_assert(sizeof(lt::math::vec4_u32) == sizeof(u32) * num_elements);
+	static_assert(sizeof(lt::math::vec4_u64) == sizeof(u64) * num_elements);
+};
+
+Suite raii = "vec4_raii"_suite = [] {
+	Case { "happy paths" } = [] {
+		ignore = vec4 {};
+		ignore = vec4 { 2.0 };
+		ignore = vec4 { 2.0, 4.0, 6.0, 8.0 };
+		ignore = vec4 { vec2 { 2.0, 4.0 }, 6.0, 8.0 };
+		ignore = vec4 { 2.0, 4.0, vec2 { 6.0, 8.0 } };
+		ignore = vec4 { vec2 { 2.0, 4.0 }, vec2 { 6.0, 8.0 } };
+		ignore = vec4 { vec3 { 2.0, 4.0, 6.0 }, 8.0 };
+		ignore = vec4 { 2.0, vec3 { 4.0, 6.0, 8.0 } };
+	};
+
+	Case { "unhappy paths" } = [] {
+	};
+
+	Case { "many" } = [] {
+		for (auto idx : std::views::iota(0, 1'000'000))
+		{
+			ignore = idx;
+			ignore = vec4 {};
+			ignore = vec4 { 2.0 };
+			ignore = vec4 { 2.0, 4.0, 6.0, 8.0 };
+			ignore = vec4 { vec2 { 2.0, 4.0 }, 6.0, 8.0 };
+			ignore = vec4 { 2.0, 4.0, vec2 { 6.0, 8.0 } };
+			ignore = vec4 { vec2 { 2.0, 4.0 }, vec2 { 6.0, 8.0 } };
+			ignore = vec4 { vec3 { 2.0, 4.0, 6.0 }, 8.0 };
+			ignore = vec4 { 2.0, vec3 { 4.0, 6.0, 8.0 } };
+		}
+	};
+
+	Case { "post default construct has correct state" } = [] {
+		const auto vec = vec4 {};
+		expect_eq(vec.x, 0.0);
+		expect_eq(vec.y, 0.0);
+		expect_eq(vec.z, 0.0);
+		expect_eq(vec.w, 0.0);
+	};
+
+	Case { "post scalar construct has correct state" } = [] {
+		const auto vec = vec4 { 2.0 };
+		expect_eq(vec.x, 2.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 2.0);
+		expect_eq(vec.w, 2.0);
+	};
+
+	Case { "post construct with x,y,z,w has correct state" } = [] {
+		const auto vec = vec4 { 1.0, 2.0, 3.0, 4.0 };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.y, 3.0);
+		expect_eq(vec.z, 4.0);
+	};
+
+	Case { "post construct with xy,z,w has correct state" } = [] {
+		const auto vec = vec4 { vec2 { 1.0, 2.0 }, 3.0, 4.0 };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+		expect_eq(vec.w, 4.0);
+	};
+
+	Case { "post construct with x,y,zw has correct state" } = [] {
+		const auto vec = vec4 { 1.0, 2.0, vec2 { 3.0, 4.0 } };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+		expect_eq(vec.w, 4.0);
+	};
+
+	Case { "post construct with x,yz,w has correct state" } = [] {
+		const auto vec = vec4 { 1.0, vec2 { 2.0, 3.0 }, 4.0 };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+		expect_eq(vec.w, 4.0);
+	};
+
+	Case { "post construct with x,y,zw has correct state" } = [] {
+		const auto vec = vec4 { 1.0, 2.0, vec2 { 3.0, 4.0 } };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+		expect_eq(vec.w, 4.0);
+	};
+
+	Case { "post construct with xy,zw has correct state" } = [] {
+		const auto vec = vec4 { vec2 { 1.0, 2.0 }, vec2 { 3.0, 4.0 } };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+		expect_eq(vec.w, 4.0);
+	};
+
+	Case { "post construct with xyz,w has correct state" } = [] {
+		const auto vec = vec4 { vec3 { 1.0, 2.0, 3.0 }, 4.0 };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+		expect_eq(vec.w, 4.0);
+	};
+
+	Case { "post construct with x,yzw has correct state" } = [] {
+		const auto vec = vec4 { 1.0, vec3 { 2.0, 3.0, 4.0 } };
+		expect_eq(vec.x, 1.0);
+		expect_eq(vec.y, 2.0);
+		expect_eq(vec.z, 3.0);
+		expect_eq(vec.w, 4.0);
+	};
+};
+
+Suite arithmetic_operators = "vec4_operators"_suite = [] {
+	Case { "operator ==" } = [] {
+		const auto lhs = vec4 { 1.0, 2.0, 3.0, 4.0 };
+
+		expect_false(lhs == vec4 { {}, 2.0, 3.0, 4.0 });
+		expect_false(lhs == vec4 { 1.0, {}, 3.0, 4.0 });
+		expect_false(lhs == vec4 { 1.0, 2.0, {}, 4.0 });
+		expect_false(lhs == vec4 { 1.0, 2.0, 3.0, {} });
+		expect_true(lhs == vec4 { 1.0, 2.0, 3.0, 4.0 });
+	};
+
+	Case { "operator !=" } = [] {
+		const auto lhs = vec4 { 1.0, 2.0, 3.0, 4.0 };
+
+		expect_true(lhs != vec4 { {}, 2.0, 3.0, 4.0 });
+		expect_true(lhs != vec4 { 1.0, {}, 3.0, 4.0 });
+		expect_true(lhs != vec4 { 1.0, 2.0, {}, 4.0 });
+		expect_true(lhs != vec4 { 1.0, 2.0, 3.0, {} });
+		expect_false(lhs != vec4 { 1.0, 2.0, 3.0, 4.0 });
+	};
+
+	Case { "operator +" } = [] {
+		const auto lhs = vec4 { 1.0, 2.0, 3.0, 4.0 };
+		const auto rhs = vec4 { 5.0, 6.0, 7.0, 8.0 };
+		expect_eq(lhs + rhs, vec4 { 6.0, 8.0, 10.0, 12.0 });
+	};
+
+	Case { "operator -" } = [] {
+		const auto lhs = vec4 { 1.0, 2.0, 3.0, 4.0 };
+		const auto rhs = vec4 { 5.0, 10.0, 15.0, 20.0 };
+		expect_eq(lhs - rhs, vec4 { -4.0, -8.0, -12.0, -16.0 });
+	};
+
+	Case { "operator *" } = [] {
+		const auto lhs = vec4 { 1.0, 2.0, 3.0, 4.0 };
+		const auto rhs = vec4 { 5.0, 6.0, 7.0, 8.0 };
+		expect_eq(lhs * rhs, vec4 { 5.0, 12.0, 21.0, 32.0 });
+	};
+
+	Case { "operator /" } = [] {
+		const auto lhs = vec4 { 5.0, 6.0, 30.0, 8.0 };
+		const auto rhs = vec4 { 1.0, 2.0, 3.0, 4.0 };
+		expect_eq(lhs / rhs, vec4 { 5.0, 3.0, 10.0, 2.0 });
+	};
+
+	Case { "operator []" } = [] {
+		auto vec = vec4 { 0.0, 1.0, 2.0, 3.0 };
+		expect_eq(vec[0], 0.0);
+		expect_eq(vec[1], 1.0);
+		expect_eq(vec[2], 2.0);
+		expect_eq(vec[3], 3.0);
+	};
+
+	Case { "operator [] const" } = [] {
+		const auto vec = vec4 { 0.0, 1.0, 2.0, 3.0 };
+		expect_eq(vec[0], 0.0);
+		expect_eq(vec[1], 1.0);
+		expect_eq(vec[2], 2.0);
+		expect_eq(vec[3], 3.0);
+	};
+};
+
+Suite utilities = "vec4_utilities"_suite = [] {
+	Case { "std::format float" } = [] {
+		auto str = std::format("{}", vec4 { 10.0000f, 30.0005f, 40.00005f, 0.0 });
+		expect_eq(str, "10, 30.0005, 40.00005, 0");
+	};
+
+	Case { "std::format int" } = [] {
+		auto str = std::format("{}", ivec4 { 10, 30, 3'000'000, 13 });
+		expect_eq(str, "10, 30, 3000000, 13");
+	};
+};