13 files changed, 1296 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..b9102d0
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,23 @@
+cmake_minimum_required(VERSION 3.16)
+
+project(math)
+
+# Library.
+
+add_library(math INTERFACE)
+
+target_include_directories(math INTERFACE
+  include)
+
+target_link_libraries(math INTERFACE
+  -lm)
+
+# Test.
+
+add_executable(math_test
+  test/mat4_test.c)
+
+target_link_libraries(math_test
+  math)
+
+target_compile_options(math_test PRIVATE -DUNIT_TEST -DNDEBUG -Wall -Wextra)
diff --git a/include/math/camera.h b/include/math/camera.h
new file mode 100644
index 0000000..1621927
--- /dev/null
+++ b/include/math/camera.h
@@ -0,0 +1,42 @@
+#pragma once
+
+#include "mat4.h"
+#include "spatial3.h"
+
+typedef struct Camera {
+  Spatial3 spatial;
+  mat4 projection;
+} Camera;
+
+/// Create an orthographic camera.
+///
+/// The camera is positioned at the origin with canonical right/up/forward
+/// vectors.
+///
+/// \param left   The coordinate for the left vertical clipping plane.
+/// \param right  The coordinate for the right vertical clipping plane.
+/// \param bottom The coordinate for the bottom horizontal clipping plane.
+/// \param top    The coordinate for the top horizontal clipping plane.
+/// \param near   The distance to the near clipping plane.
+/// \param far    The distance to the far clipping plane.
+static inline Camera camera_orthographic(R left, R right, R bottom, R top,
+                                         R near, R far) {
+  return (Camera){.spatial = spatial3_make(),
+                  .projection =
+                      mat4_ortho(left, right, bottom, top, near, far)};
+}
+
+/// Create a perspective camera.
+///
+/// The camera is positioned at the origin with canonical right/up/forward
+/// vectors.
+///
+/// \param fovy   The vertical field of view angle in degrees.
+/// \param aspect The aspect ratio that determines the field of view in the
+///               x-direction.
+/// \param near   Distance to the near clipping plane.
+/// \param far    Distance to the far clipping plane.
+static inline Camera camera_perspective(R fovy, R aspect, R near, R far) {
+  return (Camera){.spatial = spatial3_make(),
+                  .projection = mat4_perspective(fovy, aspect, near, far)};
+}
diff --git a/include/math/defs.h b/include/math/defs.h
new file mode 100644
index 0000000..acfb037
--- /dev/null
+++ b/include/math/defs.h
@@ -0,0 +1,51 @@
+#pragma once
+
+//
+// Configuration macros:
+//
+// MATH_USE_FLOAT
+//   - Use floats instead of doubles for scalar values.
+
+#include <float.h>
+#include <math.h>
+
+#ifdef MATH_USE_DOUBLE
+typedef double R;
+#define R_MIN DBL_MIN
+#define R_MAX DBL_MAX
+#else // floats
+typedef float R;
+#define R_MIN FLT_MIN
+#define R_MAX FLT_MAX
+#endif
+
+#define PI 3.14159265359
+#define INV_PI 0.31830988618
+
+/// Radians per degree.
+#define TO_RAD (PI / 180.0)
+
+/// Degrees per radian.
+#define TO_DEG (180.0 / PI)
+
+#ifdef MATH_USE_DOUBLE
+static inline double min(R a, R b) { return fmin(a, b); }
+static inline double max(R a, R b) { return fmax(a, b); }
+#else // floats
+static inline float min(R a, R b) { return fminf(a, b); }
+static inline float max(R a, R b) { return fmaxf(a, b); }
+#endif
+
+static inline R rabs(R x) { return x >= 0.0 ? x : -x; }
+static inline R clamp(R x, R low, R high) { return max(low, min(high, x)); }
+static inline R sq(R x) { return x * x; }
+static inline R sign(R x) {
+  if (x < 0) {
+    return -1;
+  } else if (x > 0) {
+    return 1;
+  } else {
+    return 0;
+  }
+}
+static inline R lerp(R a, R b, R t) { return a + (b - a) * t; }
diff --git a/include/math/float.h b/include/math/float.h
new file mode 100644
index 0000000..9d289b9
--- /dev/null
+++ b/include/math/float.h
@@ -0,0 +1,9 @@
+#pragma once
+
+#include <math.h>
+#include <stdbool.h>
+
+/// Compare two floats for equality.
+static inline bool float_eq(float a, float b, float eps) {
+  return fabsf(a - b) <= eps;
+}
diff --git a/include/math/fwd.h b/include/math/fwd.h
new file mode 100644
index 0000000..6c72ded
--- /dev/null
+++ b/include/math/fwd.h
@@ -0,0 +1,9 @@
+/// Forward declarations for all math objects.
+#pragma once
+
+typedef struct Camera Camera;
+typedef struct mat4 mat4;
+typedef struct spatial3 spatial3;
+typedef struct vec2 vec2;
+typedef struct vec3 vec3;
+typedef struct vec4 vec4;
diff --git a/include/math/mat4.h b/include/math/mat4.h
new file mode 100644
index 0000000..e6c707a
--- /dev/null
+++ b/include/math/mat4.h
@@ -0,0 +1,488 @@
+#pragma once
+
+#include "defs.h"
+#include "float.h"
+#include "vec3.h"
+#include "vec4.h"
+
+#include <stdbool.h>
+
+/// A 4x4 column-major matrix.
+typedef struct mat4 {
+  R val[4][4];
+} mat4;
+
+/// Construct a matrix from 16 values.
+static inline mat4 mat4_make(R m00, R m10, R m20, R m30, R m01, R m11, R m21,
+                             R m31, R m02, R m12, R m22, R m32, R m03, R m13,
+                             R m23, R m33) {
+  mat4 m;
+  m.val[0][0] = m00;
+  m.val[0][1] = m01;
+  m.val[0][2] = m02;
+  m.val[0][3] = m03;
+
+  m.val[1][0] = m10;
+  m.val[1][1] = m11;
+  m.val[1][2] = m12;
+  m.val[1][3] = m13;
+
+  m.val[2][0] = m20;
+  m.val[2][1] = m21;
+  m.val[2][2] = m22;
+  m.val[2][3] = m23;
+
+  m.val[3][0] = m30;
+  m.val[3][1] = m31;
+  m.val[3][2] = m32;
+  m.val[3][3] = m33;
+  return m;
+}
+
+/// Construct a matrix from a column-major matrix array.
+static inline mat4 mat4_from_array(const R M[16]) {
+  mat4 m;
+  m.val[0][0] = M[0];
+  m.val[0][1] = M[1];
+  m.val[0][2] = M[2];
+  m.val[0][3] = M[3];
+
+  m.val[1][0] = M[4];
+  m.val[1][1] = M[5];
+  m.val[1][2] = M[6];
+  m.val[1][3] = M[7];
+
+  m.val[2][0] = M[8];
+  m.val[2][1] = M[9];
+  m.val[2][2] = M[10];
+  m.val[2][3] = M[11];
+
+  m.val[3][0] = M[12];
+  m.val[3][1] = M[13];
+  m.val[3][2] = M[14];
+  m.val[3][3] = M[15];
+  return m;
+}
+
+/// Construct the identity matrix.
+static inline mat4 mat4_id() {
+  return mat4_make(1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
+                   0.0, 0.0, 0.0, 1.0);
+}
+
+/// Construct a matrix from 4 column vectors.
+static inline mat4 mat4_from_vec4(vec4 v0, vec4 v1, vec4 v2, vec4 v3) {
+  return mat4_make(v0.x, v0.y, v0.z, v0.w, v1.x, v1.y, v1.z, v1.w, v2.x, v2.y,
+                   v2.z, v2.w, v3.x, v3.y, v3.z, v3.w);
+}
+
+/// Construct a transformation matrix from 4 vectors.
+static inline mat4 mat4_from_vec3(vec3 right, vec3 up, vec3 forward,
+                                  vec3 position) {
+  return mat4_make(right.x, right.y, right.z, 0.0, up.x, up.y, up.z, 0.0,
+                   forward.x, forward.y, forward.z, 0.0, position.x, position.y,
+                   position.z, 1.0);
+}
+
+/// Return the value at the specified position.
+static inline R mat4_at(mat4 m, int row, int col) { return m.val[col][row]; }
+
+/// Return the matrix's first column.
+static inline vec3 mat4_v0(mat4 m) { return *((vec3*)m.val[0]); }
+
+/// Return the matrix's second column.
+static inline vec3 mat4_v1(mat4 m) { return *((vec3*)m.val[1]); }
+
+/// Return the matrix's third column.
+static inline vec3 mat4_v2(mat4 m) { return *((vec3*)m.val[2]); }
+
+/// Return the matrix's fourth column.
+static inline vec3 mat4_v3(mat4 m) { return *((vec3*)m.val[3]); }
+
+/// Set the matrix's first column.
+static inline void mat4_set_v0(mat4* m, vec3 v) { *((vec3*)m->val[0]) = v; }
+
+/// Set the matrix's second column.
+static inline void mat4_set_v1(mat4* m, vec3 v) { *((vec3*)m->val[1]) = v; }
+
+/// Set the matrix's third column.
+static inline void mat4_set_v2(mat4* m, vec3 v) { *((vec3*)m->val[2]) = v; }
+
+/// Set the matrix's fourth column.
+static inline void mat4_set_v3(mat4* m, vec3 v) { *((vec3*)m->val[3]) = v; }
+
+/// Set the 3x3 portion of the first matrix equal to the 3x3 portion of the
+/// second matrix.
+static inline void mat4_set_3x3(mat4* m, mat4 n) {
+  m->val[0][0] = n.val[0][0];
+  m->val[0][1] = n.val[0][1];
+  m->val[0][2] = n.val[0][2];
+
+  m->val[1][0] = n.val[1][0];
+  m->val[1][1] = n.val[1][1];
+  m->val[1][2] = n.val[1][2];
+
+  m->val[2][0] = n.val[2][0];
+  m->val[2][1] = n.val[2][1];
+  m->val[2][2] = n.val[2][2];
+}
+
+/// Multiply two matrices.
+/// A * B = AB.
+static inline mat4 mat4_mul(mat4 A, mat4 B) {
+  R m00 = mat4_at(A, 0, 0) * mat4_at(B, 0, 0) +
+          mat4_at(A, 0, 1) * mat4_at(B, 1, 0) +
+          mat4_at(A, 0, 2) * mat4_at(B, 2, 0) +
+          mat4_at(A, 0, 3) * mat4_at(B, 3, 0);
+  R m01 = mat4_at(A, 0, 0) * mat4_at(B, 0, 1) +
+          mat4_at(A, 0, 1) * mat4_at(B, 1, 1) +
+          mat4_at(A, 0, 2) * mat4_at(B, 2, 1) +
+          mat4_at(A, 0, 3) * mat4_at(B, 3, 1);
+  R m02 = mat4_at(A, 0, 0) * mat4_at(B, 0, 2) +
+          mat4_at(A, 0, 1) * mat4_at(B, 1, 2) +
+          mat4_at(A, 0, 2) * mat4_at(B, 2, 2) +
+          mat4_at(A, 0, 3) * mat4_at(B, 3, 2);
+  R m03 = mat4_at(A, 0, 0) * mat4_at(B, 0, 3) +
+          mat4_at(A, 0, 1) * mat4_at(B, 1, 3) +
+          mat4_at(A, 0, 2) * mat4_at(B, 2, 3) +
+          mat4_at(A, 0, 3) * mat4_at(B, 3, 3);
+
+  R m10 = mat4_at(A, 1, 0) * mat4_at(B, 0, 0) +
+          mat4_at(A, 1, 1) * mat4_at(B, 1, 0) +
+          mat4_at(A, 1, 2) * mat4_at(B, 2, 0) +
+          mat4_at(A, 1, 3) * mat4_at(B, 3, 0);
+  R m11 = mat4_at(A, 1, 0) * mat4_at(B, 0, 1) +
+          mat4_at(A, 1, 1) * mat4_at(B, 1, 1) +
+          mat4_at(A, 1, 2) * mat4_at(B, 2, 1) +
+          mat4_at(A, 1, 3) * mat4_at(B, 3, 1);
+  R m12 = mat4_at(A, 1, 0) * mat4_at(B, 0, 2) +
+          mat4_at(A, 1, 1) * mat4_at(B, 1, 2) +
+          mat4_at(A, 1, 2) * mat4_at(B, 2, 2) +
+          mat4_at(A, 1, 3) * mat4_at(B, 3, 2);
+  R m13 = mat4_at(A, 1, 0) * mat4_at(B, 0, 3) +
+          mat4_at(A, 1, 1) * mat4_at(B, 1, 3) +
+          mat4_at(A, 1, 2) * mat4_at(B, 2, 3) +
+          mat4_at(A, 1, 3) * mat4_at(B, 3, 3);
+
+  R m20 = mat4_at(A, 2, 0) * mat4_at(B, 0, 0) +
+          mat4_at(A, 2, 1) * mat4_at(B, 1, 0) +
+          mat4_at(A, 2, 2) * mat4_at(B, 2, 0) +
+          mat4_at(A, 2, 3) * mat4_at(B, 3, 0);
+  R m21 = mat4_at(A, 2, 0) * mat4_at(B, 0, 1) +
+          mat4_at(A, 2, 1) * mat4_at(B, 1, 1) +
+          mat4_at(A, 2, 2) * mat4_at(B, 2, 1) +
+          mat4_at(A, 2, 3) * mat4_at(B, 3, 1);
+  R m22 = mat4_at(A, 2, 0) * mat4_at(B, 0, 2) +
+          mat4_at(A, 2, 1) * mat4_at(B, 1, 2) +
+          mat4_at(A, 2, 2) * mat4_at(B, 2, 2) +
+          mat4_at(A, 2, 3) * mat4_at(B, 3, 2);
+  R m23 = mat4_at(A, 2, 0) * mat4_at(B, 0, 3) +
+          mat4_at(A, 2, 1) * mat4_at(B, 1, 3) +
+          mat4_at(A, 2, 2) * mat4_at(B, 2, 3) +
+          mat4_at(A, 2, 3) * mat4_at(B, 3, 3);
+
+  R m30 = mat4_at(A, 3, 0) * mat4_at(B, 0, 0) +
+          mat4_at(A, 3, 1) * mat4_at(B, 1, 0) +
+          mat4_at(A, 3, 2) * mat4_at(B, 2, 0) +
+          mat4_at(A, 3, 3) * mat4_at(B, 3, 0);
+  R m31 = mat4_at(A, 3, 0) * mat4_at(B, 0, 1) +
+          mat4_at(A, 3, 1) * mat4_at(B, 1, 1) +
+          mat4_at(A, 3, 2) * mat4_at(B, 2, 1) +
+          mat4_at(A, 3, 3) * mat4_at(B, 3, 1);
+  R m32 = mat4_at(A, 3, 0) * mat4_at(B, 0, 2) +
+          mat4_at(A, 3, 1) * mat4_at(B, 1, 2) +
+          mat4_at(A, 3, 2) * mat4_at(B, 2, 2) +
+          mat4_at(A, 3, 3) * mat4_at(B, 3, 2);
+  R m33 = mat4_at(A, 3, 0) * mat4_at(B, 0, 3) +
+          mat4_at(A, 3, 1) * mat4_at(B, 1, 3) +
+          mat4_at(A, 3, 2) * mat4_at(B, 2, 3) +
+          mat4_at(A, 3, 3) * mat4_at(B, 3, 3);
+
+  return mat4_make(m00, m01, m02, m03, m10, m11, m12, m13, m20, m21, m22, m23,
+                   m30, m31, m32, m33);
+}
+
+/// Return the translation component of the matrix.
+static inline mat4 mat4_translation(mat4 m) {
+  return mat4_make(1.0, 0.0, 0.0, mat4_at(m, 0, 3), 0.0, 1.0, 0.0,
+                   mat4_at(m, 1, 3), 0.0, 0.0, 1.0, mat4_at(m, 2, 3), 0.0, 0.0,
+                   0.0, 1.0);
+}
+
+/// Return the rotation component of the matrix.
+static inline mat4 mat4_rotation(mat4 m) {
+  return mat4_make(mat4_at(m, 0, 0), mat4_at(m, 0, 1), mat4_at(m, 0, 2), 0.0,
+                   mat4_at(m, 1, 0), mat4_at(m, 1, 1), mat4_at(m, 1, 2), 0.0,
+                   mat4_at(m, 2, 0), mat4_at(m, 2, 1), mat4_at(m, 2, 2), 0.0,
+                   0.0, 0.0, 0.0, 1.0);
+}
+
+/// Create an X-axis rotation matrix.
+static inline mat4 mat4_rotx(R angle) {
+  const R s = sin(angle);
+  const R c = cos(angle);
+  return mat4_make(1, 0, 0, 0, 0, c, -s, 0, 0, s, c, 0, 0, 0, 0, 1);
+}
+
+/// Create a Y-axis rotation matrix.
+static inline mat4 mat4_roty(R angle) {
+  const R s = sin(angle);
+  const R c = cos(angle);
+  return mat4_make(c, 0, s, 0, 0, 1, 0, 0, -s, 0, c, 0, 0, 0, 0, 1);
+}
+
+/// Create a Z-axis rotation matrix.
+static inline mat4 mat4_rotz(R angle) {
+  const R s = sin(angle);
+  const R c = cos(angle);
+  return mat4_make(c, -s, 0, 0, s, c, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
+}
+
+/// Create a rotation matrix.
+static inline mat4 mat4_rot(vec3 axis, R angle) {
+  const R s = sin(angle);
+  const R c = cos(angle);
+  const R x = axis.x;
+  const R y = axis.y;
+  const R z = axis.z;
+  const R xy = x * y;
+  const R xz = x * z;
+  const R yz = y * z;
+  const R sx = s * x;
+  const R sy = s * y;
+  const R sz = s * z;
+  const R omc = 1.0 - c;
+  return mat4_make(c + omc * x * x, omc * xy - sz, omc * xz + sy, 0,
+                   omc * xy + sz, c + omc * y * y, omc * yz - sx, 0,
+                   omc * xz - sy, omc * yz + sx, c + omc * z * z, 0, 0, 0, 0,
+                   1);
+}
+
+/// Create a scaling matrix.
+static inline mat4 mat4_scale(vec3 s) {
+  return mat4_make(s.x, 0, 0, 0, 0, s.y, 0, 0, 0, 0, s.z, 0, 0, 0, 0, 1);
+}
+
+/// Create a translation matrix.
+static inline mat4 mat4_translate(vec3 v) {
+  return mat4_make(1, 0, 0, v.x, 0, 1, 0, v.y, 0, 0, 1, v.z, 0, 0, 0, 1);
+}
+
+/// The X-axis reflection matrix.
+static inline mat4 mat4_reflectx() {
+  return mat4_make(-1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
+}
+
+/// The Y-axis reflection matrix.
+static inline mat4 mat4_reflecty() {
+  return mat4_make(1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
+}
+
+/// The Z-axis reflection matrix.
+static inline mat4 mat4_reflectz() {
+  return mat4_make(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1);
+}
+
+/// Create a transformation matrix from the given forward vector.
+static inline mat4 mat4_from_forward(vec3 forward) {
+  const vec3 f = vec3_normalize(forward);
+  const vec3 r = vec3_normalize(vec3_cross(f, up3()));
+  const vec3 u = vec3_normalize(vec3_cross(r, f));
+  return mat4_make(r.x, u.x, -f.x, 0.0, r.y, u.y, -f.y, 0.0, r.z, u.z, -f.z,
+                   0.0, 0.0, 0.0, 0.0, 1.0);
+}
+
+/// Create a transformation matrix.
+static inline mat4 mat4_lookat(vec3 position, vec3 target, vec3 up) {
+  const vec3 fwd = vec3_normalize(vec3_sub(target, position));
+  const vec3 right = vec3_normalize(vec3_cross(fwd, up));
+  up = vec3_normalize(vec3_cross(right, fwd));
+  return mat4_from_vec3(right, up, fwd, position);
+}
+
+/// Create an orthographic projection matrix.
+/// \param left   The coordinate for the left vertical clipping plane.
+/// \param right  The coordinate for the right vertical clipping plane.
+/// \param bottom The coordinate for the bottom horizontal clipping plane.
+/// \param top    The coordinate for the top horizontal clipping plane.
+/// \param near   The distance to the near clipping plane.
+/// \param far    The distance to the far clipping plane.
+static inline mat4 mat4_ortho(R left, R right, R bottom, R top, R near, R far) {
+  const R tx = -(right + left) / (right - left);
+  const R ty = -(top + bottom) / (top - bottom);
+  const R tz = -(far + near) / (far - near);
+  return mat4_make(2 / (right - left), 0, 0, tx, 0, 2 / (top - bottom), 0, ty,
+                   0, 0, -2 / (far - near), tz, 0, 0, 0, 1);
+}
+
+/// Create a perspective projection matrix.
+/// \param fovy   The vertical field of view angle in degrees.
+/// \param aspect The aspect ratio that determines the field of view in the
+///               x-direction.
+/// \param near   Distance to the near clipping plane.
+/// \param far    Distance to the far clipping plane.
+static inline mat4 mat4_perspective(R fovy, R aspect, R near, R far) {
+  R f = tan(fovy / 2.0);
+  assert(f > 0.0);
+  f = 1.0 / f;
+  const R a = near - far;
+  return mat4_make(f / aspect, 0, 0, 0, 0, f, 0, 0, 0, 0, (far + near) / a,
+                   (2 * far * near / a), 0, 0, -1, 0);
+}
+
+/// Create the inverse of a perspective projection matrix.
+/// \param fovy   The vertical field of view angle in degrees.
+/// \param aspect The aspect ratio that determines the field of view in the
+///               x-direction.
+/// \param near   Distance to the near clipping plane.
+/// \param far    Distance to the far clipping plane.
+static inline mat4 mat4_perspective_inverse(R fovy, R aspect, R near, R far) {
+  R f = tan(fovy / 2.0);
+  assert(f > 0.0);
+  f = 1.0 / f;
+  const R a = far * near;
+  const R P32 = 0.5 * (near - far) / a;
+  const R P33 = 0.5 * (far + near) / a;
+  return mat4_make(aspect / f, 0, 0, 0, 0, 1.0f / f, 0, 0, 0, 0, 0, -1, 0, 0,
+                   P32, P33);
+}
+
+/// Return the matrix's determinant.
+static inline R mat4_det(mat4 m) {
+  const R* M = (const R*)(m.val);
+  const R inv0 = M[5] * M[10] * M[15] - M[5] * M[11] * M[14] -
+                 M[9] * M[6] * M[15] + M[9] * M[7] * M[14] +
+                 M[13] * M[6] * M[11] - M[13] * M[7] * M[10];
+  const R inv1 = -M[4] * M[10] * M[15] + M[4] * M[11] * M[14] +
+                 M[8] * M[6] * M[15] - M[8] * M[7] * M[14] -
+                 M[12] * M[6] * M[11] + M[12] * M[7] * M[10];
+  const R inv2 = M[4] * M[9] * M[15] - M[4] * M[11] * M[13] -
+                 M[8] * M[5] * M[15] + M[8] * M[7] * M[13] +
+                 M[12] * M[5] * M[11] - M[12] * M[7] * M[9];
+  const R inv3 = -M[4] * M[9] * M[14] + M[4] * M[10] * M[13] +
+                 M[8] * M[5] * M[14] - M[8] * M[6] * M[13] -
+                 M[12] * M[5] * M[10] + M[12] * M[6] * M[9];
+  return M[0] * inv0 + M[1] * inv1 + M[2] * inv2 + M[3] * inv3;
+}
+
+/// Invert the matrix.
+static inline mat4 mat4_inverse(mat4 m) {
+  const R* M = (const R*)(m.val);
+
+  R inv[16];
+  inv[0] = M[5] * M[10] * M[15] - M[5] * M[11] * M[14] - M[9] * M[6] * M[15] +
+           M[9] * M[7] * M[14] + M[13] * M[6] * M[11] - M[13] * M[7] * M[10];
+  inv[4] = -M[4] * M[10] * M[15] + M[4] * M[11] * M[14] + M[8] * M[6] * M[15] -
+           M[8] * M[7] * M[14] - M[12] * M[6] * M[11] + M[12] * M[7] * M[10];
+  inv[8] = M[4] * M[9] * M[15] - M[4] * M[11] * M[13] - M[8] * M[5] * M[15] +
+           M[8] * M[7] * M[13] + M[12] * M[5] * M[11] - M[12] * M[7] * M[9];
+  inv[12] = -M[4] * M[9] * M[14] + M[4] * M[10] * M[13] + M[8] * M[5] * M[14] -
+            M[8] * M[6] * M[13] - M[12] * M[5] * M[10] + M[12] * M[6] * M[9];
+  inv[1] = -M[1] * M[10] * M[15] + M[1] * M[11] * M[14] + M[9] * M[2] * M[15] -
+           M[9] * M[3] * M[14] - M[13] * M[2] * M[11] + M[13] * M[3] * M[10];
+  inv[5] = M[0] * M[10] * M[15] - M[0] * M[11] * M[14] - M[8] * M[2] * M[15] +
+           M[8] * M[3] * M[14] + M[12] * M[2] * M[11] - M[12] * M[3] * M[10];
+  inv[9] = -M[0] * M[9] * M[15] + M[0] * M[11] * M[13] + M[8] * M[1] * M[15] -
+           M[8] * M[3] * M[13] - M[12] * M[1] * M[11] + M[12] * M[3] * M[9];
+  inv[13] = M[0] * M[9] * M[14] - M[0] * M[10] * M[13] - M[8] * M[1] * M[14] +
+            M[8] * M[2] * M[13] + M[12] * M[1] * M[10] - M[12] * M[2] * M[9];
+  inv[2] = M[1] * M[6] * M[15] - M[1] * M[7] * M[14] - M[5] * M[2] * M[15] +
+           M[5] * M[3] * M[14] + M[13] * M[2] * M[7] - M[13] * M[3] * M[6];
+  inv[6] = -M[0] * M[6] * M[15] + M[0] * M[7] * M[14] + M[4] * M[2] * M[15] -
+           M[4] * M[3] * M[14] - M[12] * M[2] * M[7] + M[12] * M[3] * M[6];
+  inv[10] = M[0] * M[5] * M[15] - M[0] * M[7] * M[13] - M[4] * M[1] * M[15] +
+            M[4] * M[3] * M[13] + M[12] * M[1] * M[7] - M[12] * M[3] * M[5];
+  inv[14] = -M[0] * M[5] * M[14] + M[0] * M[6] * M[13] + M[4] * M[1] * M[14] -
+            M[4] * M[2] * M[13] - M[12] * M[1] * M[6] + M[12] * M[2] * M[5];
+  inv[3] = -M[1] * M[6] * M[11] + M[1] * M[7] * M[10] + M[5] * M[2] * M[11] -
+           M[5] * M[3] * M[10] - M[9] * M[2] * M[7] + M[9] * M[3] * M[6];
+  inv[7] = M[0] * M[6] * M[11] - M[0] * M[7] * M[10] - M[4] * M[2] * M[11] +
+           M[4] * M[3] * M[10] + M[8] * M[2] * M[7] - M[8] * M[3] * M[6];
+  inv[11] = -M[0] * M[5] * M[11] + M[0] * M[7] * M[9] + M[4] * M[1] * M[11] -
+            M[4] * M[3] * M[9] - M[8] * M[1] * M[7] + M[8] * M[3] * M[5];
+  inv[15] = M[0] * M[5] * M[10] - M[0] * M[6] * M[9] - M[4] * M[1] * M[10] +
+            M[4] * M[2] * M[9] + M[8] * M[1] * M[6] - M[8] * M[2] * M[5];
+
+  R det = M[0] * inv[0] + M[1] * inv[4] + M[2] * inv[8] + M[3] * inv[12];
+  assert(det != 0.0);
+  det = 1.0 / det;
+  return mat4_make(inv[0] * det, inv[4] * det, inv[8] * det, inv[12] * det,
+                   inv[1] * det, inv[5] * det, inv[9] * det, inv[13] * det,
+                   inv[2] * det, inv[6] * det, inv[10] * det, inv[14] * det,
+                   inv[3] * det, inv[7] * det, inv[11] * det, inv[15] * det);
+}
+
+/// Invert the transformation matrix.
+/// This is much faster than the more general inverse() function, but assumes
+/// that the matrix is of the form TR, where T is a translation and R a
+/// rotation.
+static inline mat4 mat4_inverse_transform(mat4 m) {
+  const vec3 r = mat4_v0(m);
+  const vec3 u = mat4_v1(m);
+  const vec3 f = mat4_v2(m);
+  const vec3 t = mat4_v3(m);
+  return mat4_make(r.x, r.y, r.z, -vec3_dot(r, t), u.x, u.y, u.z,
+                   -vec3_dot(u, t), f.x, f.y, f.z, -vec3_dot(f, t), 0.0, 0.0,
+                   0.0, 1.0);
+}
+
+/// Transpose the matrix.
+static inline mat4 mat4_transpose(mat4 m) {
+  return mat4_make(
+      mat4_at(m, 0, 0), mat4_at(m, 1, 0), mat4_at(m, 2, 0), mat4_at(m, 3, 0),
+      mat4_at(m, 0, 1), mat4_at(m, 1, 1), mat4_at(m, 2, 1), mat4_at(m, 3, 1),
+      mat4_at(m, 0, 2), mat4_at(m, 1, 2), mat4_at(m, 2, 2), mat4_at(m, 3, 2),
+      mat4_at(m, 0, 3), mat4_at(m, 1, 3), mat4_at(m, 2, 3), mat4_at(m, 3, 3));
+}
+
+/// Transform the vector with the matrix.
+static inline vec3 mat4_mul_vec3(mat4 m, vec3 v, R w) {
+  vec3 u;
+  u.x = mat4_at(m, 0, 0) * v.x + mat4_at(m, 0, 1) * v.y +
+        mat4_at(m, 0, 2) * v.z + mat4_at(m, 0, 3) * w;
+  u.y = mat4_at(m, 1, 0) * v.x + mat4_at(m, 1, 1) * v.y +
+        mat4_at(m, 1, 2) * v.z + mat4_at(m, 1, 3) * w;
+  u.z = mat4_at(m, 2, 0) * v.x + mat4_at(m, 2, 1) * v.y +
+        mat4_at(m, 2, 2) * v.z + mat4_at(m, 2, 3) * w;
+  return u;
+}
+
+/// Return the vector multiplied by the matrix.
+static inline vec4 mat4_mul_vec4(mat4 m, vec4 v) {
+  vec4 u;
+  u.x = mat4_at(m, 0, 0) * v.x + mat4_at(m, 0, 1) * v.y +
+        mat4_at(m, 0, 2) * v.z + mat4_at(m, 0, 3) * v.w;
+  u.y = mat4_at(m, 1, 0) * v.x + mat4_at(m, 1, 1) * v.y +
+        mat4_at(m, 1, 2) * v.z + mat4_at(m, 1, 3) * v.w;
+  u.z = mat4_at(m, 2, 0) * v.x + mat4_at(m, 2, 1) * v.y +
+        mat4_at(m, 2, 2) * v.z + mat4_at(m, 2, 3) * v.w;
+  u.w = mat4_at(m, 3, 0) * v.x + mat4_at(m, 3, 1) * v.y +
+        mat4_at(m, 3, 2) * v.z + mat4_at(m, 3, 3) * v.w;
+  return u;
+}
+
+/// Compare two matrices for equality.
+/// Returns true if the difference between each ij-value across matrices is
+/// within |eps|, false if there is at least one ij-value difference that is
+/// greater than eps.
+static inline bool mat4_eq(mat4 m, mat4 w, float eps) {
+  return (float_eq(mat4_at(m, 0, 0), mat4_at(w, 0, 0), eps) &&
+          float_eq(mat4_at(m, 0, 1), mat4_at(w, 0, 1), eps) &&
+          float_eq(mat4_at(m, 0, 2), mat4_at(w, 0, 2), eps) &&
+          float_eq(mat4_at(m, 0, 3), mat4_at(w, 0, 3), eps) &&
+
+          float_eq(mat4_at(m, 1, 0), mat4_at(w, 1, 0), eps) &&
+          float_eq(mat4_at(m, 1, 1), mat4_at(w, 1, 1), eps) &&
+          float_eq(mat4_at(m, 1, 2), mat4_at(w, 1, 2), eps) &&
+          float_eq(mat4_at(m, 1, 3), mat4_at(w, 1, 3), eps) &&
+
+          float_eq(mat4_at(m, 2, 0), mat4_at(w, 2, 0), eps) &&
+          float_eq(mat4_at(m, 2, 1), mat4_at(w, 2, 1), eps) &&
+          float_eq(mat4_at(m, 2, 2), mat4_at(w, 2, 2), eps) &&
+          float_eq(mat4_at(m, 2, 3), mat4_at(w, 2, 3), eps) &&
+
+          float_eq(mat4_at(m, 3, 0), mat4_at(w, 3, 0), eps) &&
+          float_eq(mat4_at(m, 3, 1), mat4_at(w, 3, 1), eps) &&
+          float_eq(mat4_at(m, 3, 2), mat4_at(w, 3, 2), eps) &&
+          float_eq(mat4_at(m, 3, 3), mat4_at(w, 3, 3), eps));
+}
diff --git a/include/math/spatial3.h b/include/math/spatial3.h
new file mode 100644
index 0000000..f8caf5d
--- /dev/null
+++ b/include/math/spatial3.h
@@ -0,0 +1,194 @@
+#pragma once
+
+#include "mat4.h"
+#include "vec3.h"
+
+/// An object in 3D space.
+typedef struct Spatial3 {
+  vec3 p; // Position.
+  vec3 r; // Right vector.
+  vec3 u; // Up vector.
+  vec3 f; // Forward vector.
+} Spatial3;
+
+/// Construct a spatial with position 0 and canonical direction vectors.
+static inline Spatial3 spatial3_make() {
+  return (Spatial3){.p = zero3(), .r = right3(), .u = up3(), .f = forward3()};
+}
+
+/// Return the spatial's transformation matrix (from local to world
+/// coordinates).
+static inline mat4 spatial3_transform(const Spatial3* spatial) {
+  const vec3 p = spatial->p;
+  const vec3 r = spatial->r;
+  const vec3 u = spatial->u;
+  const vec3 f = spatial->f;
+  return mat4_make(r.x, u.x, -f.x, p.x, r.y, u.y, -f.y, p.y, r.z, u.z, -f.z,
+                   p.z, 0.0, 0.0, 0.0, 1.0f);
+}
+
+/// Return the spatial's inverse transformation matrix (from world to local
+/// coordinates).
+static inline mat4 spatial3_inverse_transform(const Spatial3* spatial) {
+  return mat4_inverse_transform(spatial3_transform(spatial));
+}
+
+/// Move the spatial by the given vector.
+static inline void spatial3_move(Spatial3* spatial, vec3 v) {
+  spatial->p = vec3_add(spatial->p, v);
+}
+
+/// Move the spatial along its forward vector.
+static inline void spatial3_move_forwards(Spatial3* spatial, R speed) {
+  spatial->p = vec3_add(spatial->p, vec3_scale(spatial->f, speed));
+}
+
+/// Move the spatial along its backwards vector.
+static inline void spatial3_move_backwards(Spatial3* spatial, R speed) {
+  spatial->p = vec3_add(spatial->p, vec3_scale(spatial->f, -speed));
+}
+
+/// Move the spatial along its left vector.
+static inline void spatial3_move_left(Spatial3* spatial, R speed) {
+  spatial->p = vec3_add(spatial->p, vec3_scale(spatial->r, -speed));
+}
+
+/// Move the spatial along its right vector.
+static inline void spatial3_move_right(Spatial3* spatial, R speed) {
+  spatial->p = vec3_add(spatial->p, vec3_scale(spatial->r, speed));
+}
+
+/// Move the spatial along the global up vector.
+static inline void spatial3_move_up(Spatial3* spatial, R speed) {
+  spatial->p = vec3_add(spatial->p, vec3_scale(up3(), speed));
+}
+
+/// Move the spatial along the global down vector.
+static inline void spatial3_move_down(Spatial3* spatial, R speed) {
+  spatial->p = vec3_add(spatial->p, vec3_scale(up3(), -speed));
+}
+
+/// Rotate the spatial about the given axis by the given angle.
+static inline void spatial3_rotate(Spatial3* spatial, vec3 axis, R angle) {
+  mat4 transf = spatial3_transform(spatial);
+  const vec3 local_axis =
+      vec3_normalize(mat4_mul_vec3(mat4_inverse_transform(transf), axis, 0.0));
+  transf = mat4_mul(transf, mat4_rot(local_axis, angle));
+  spatial->r = vec3_normalize(mat4_v0(transf));
+  spatial->u = vec3_normalize(mat4_v1(transf));
+  spatial->f = vec3_normalize(vec3_neg(mat4_v2(transf)));
+}
+
+/// Rotate the spatial about its local y axis.
+static inline void spatial3_yaw(Spatial3* spatial, const R angle) {
+  const R sa = sin(angle);
+  const R ca = cos(angle);
+  spatial->f = vec3_normalize(
+      vec3_sub(vec3_scale(spatial->f, ca), vec3_scale(spatial->r, sa)));
+  spatial->r = vec3_normalize(vec3_cross(spatial->f, spatial->u));
+}
+
+/// Rotate the spatial about its local x axis.
+static inline void spatial3_pitch(Spatial3* spatial, const R angle) {
+  const R sa = sin(angle);
+  const R ca = cos(angle);
+  spatial->f = vec3_normalize(
+      vec3_add(vec3_scale(spatial->f, ca), vec3_scale(spatial->u, sa)));
+  spatial->u = vec3_normalize(vec3_cross(spatial->r, spatial->f));
+}
+
+/// Rotate the spatial about its local z axis.
+static inline void spatial3_roll(Spatial3* spatial, const R angle) {
+  const R sa = sin(angle);
+  const R ca = cos(angle);
+  spatial->u = vec3_normalize(
+      vec3_sub(vec3_scale(spatial->u, ca), vec3_scale(spatial->r, sa)));
+  spatial->r = vec3_normalize(vec3_cross(spatial->f, spatial->u));
+}
+
+/// Set the spatial's transformation matrix.
+static inline void spatial3_set_transform(Spatial3* spatial, mat4 transform) {
+  spatial->r = mat4_v0(transform);
+  spatial->u = mat4_v1(transform);
+  spatial->f = mat4_v2(transform);
+  spatial->p = mat4_v3(transform);
+}
+
+static inline void spatial3_set_forward(Spatial3* spatial, vec3 forward) {
+  spatial->f = vec3_normalize(forward);
+  // Use aux vector to define right vector orthogonal to forward.
+  if (vec3_eq(vec3_abs(spatial->f), up3())) {
+    spatial->r = vec3_normalize(vec3_cross(spatial->f, forward3()));
+  } else {
+    spatial->r = vec3_normalize(vec3_cross(spatial->f, up3()));
+  }
+  spatial->u = vec3_normalize(vec3_cross(spatial->r, spatial->f));
+}
+
+/// Make the spatial look at the given target.
+static inline void spatial3_lookat(Spatial3* spatial, vec3 target) {
+  spatial3_set_forward(spatial, vec3_sub(target, spatial->p));
+}
+
+/// Make the spatial look at the given target.
+static inline void spatial3_lookat_spatial(Spatial3* spatial,
+                                           const Spatial3* target) {
+  spatial3_set_forward(spatial, vec3_sub(target->p, spatial->p));
+}
+
+/// Make the spatial orbit around the given target.
+/// \param target  Target position.
+/// \param radius  Radial distance.
+/// \param azimuth Azimuthal (horizontal) angle.
+/// \param zenith  Polar (vertical) angle.
+static inline void spatial3_orbit(Spatial3* spatial, vec3 target, R radius,
+                                  R azimuth, R zenith) {
+  const R sx = sin(azimuth);
+  const R sy = sin(zenith);
+  const R cx = cos(azimuth);
+  const R cy = cos(zenith);
+  spatial->p = (vec3){target.x + radius * cy * sx, target.y + radius * sy,
+                      target.z + radius * cx * cy};
+}
+
+/// Make the spatial orbit around the given target.
+/// \param target  Target spatial.
+/// \param radius  Radial distance.
+/// \param azimuth Azimuthal (horizontal) angle.
+/// \param zenith  Polar (vertical) angle.
+static inline void spatial3_orbit_spatial(Spatial3* spatial,
+                                          const Spatial3* target, R radius,
+                                          R azimuth, R zenith) {
+  spatial3_orbit(spatial, target->p, radius, azimuth, zenith);
+}
+
+// The functions below are provided so that client code can work with a Spatial3
+// with no assumptions on the data type's definition.
+
+/// Return the spatial's position.
+static inline vec3 spatial3_position(const Spatial3* spatial) {
+  return spatial->p;
+}
+
+/// Return the spatial's right vector.
+static inline vec3 spatial3_right(const Spatial3* spatial) {
+  return spatial->r;
+}
+
+/// Return the spatial's up vector.
+static inline vec3 spatial3_up(const Spatial3* spatial) { return spatial->u; }
+
+/// Return the spatial's forward vector.
+static inline vec3 spatial3_forward(const Spatial3* spatial) {
+  return spatial->f;
+}
+
+static inline void spatial3_set_position(Spatial3* spatial, vec3 p) {
+  spatial->p = p;
+}
+
+static inline void spatial3_setx(Spatial3* spatial, R x) { spatial->p.x = x; }
+
+static inline void spatial3_sety(Spatial3* spatial, R y) { spatial->p.y = y; }
+
+static inline void spatial3_setz(Spatial3* spatial, R z) { spatial->p.z = z; }
diff --git a/include/math/vec.h b/include/math/vec.h
new file mode 100644
index 0000000..7b870e1
--- /dev/null
+++ b/include/math/vec.h
@@ -0,0 +1,16 @@
+// Common functions for vectors.
+//
+// This header file contains functions that operate with different kinds of
+// vectors that would result in circular dependencies if defined in any of the
+// vec2/3/4 header files.
+#pragma once
+
+#include "vec2.h"
+#include "vec3.h"
+#include "vec4.h"
+
+/// Construct a 3D vector from a 2D one with z = 0.
+static inline vec3 vec3_from_vec2(vec2 v) { return vec3{v.x, v.y, 0, 0}; }
+
+/// Project a 4D vector onto w=0.
+static inline vec3 vec3_from_vec4(vec4 v) { return vec3{v.x, v.y, v.z, 0}; }
diff --git a/include/math/vec2.h b/include/math/vec2.h
new file mode 100644
index 0000000..0a3f692
--- /dev/null
+++ b/include/math/vec2.h
@@ -0,0 +1,10 @@
+#pragma once
+
+#include "defs.h"
+
+typedef struct vec2 {
+  R x, y;
+} vec2;
+
+/// Construct a vector from 2 coordinates.
+static inline vec2 vec2_make(R x, R y) { return (vec2){x, y}; }
diff --git a/include/math/vec3.h b/include/math/vec3.h
new file mode 100644
index 0000000..3c3b053
--- /dev/null
+++ b/include/math/vec3.h
@@ -0,0 +1,143 @@
+#pragma once
+
+#include "defs.h"
+
+#include <assert.h>
+#include <stdbool.h>
+
+/// A 3D vector.
+typedef struct vec3 {
+  R x, y, z;
+} vec3;
+
+/// Construct a vector from 3 coordinates.
+static inline vec3 vec3_make(R x, R y, R z) { return (vec3){x, y, z}; }
+
+/// Construct a vector from an array.
+static inline vec3 vec3_from_array(const R xyz[3]) {
+  return (vec3){xyz[0], xyz[1], xyz[2]};
+}
+
+/// Construct a vector from a single scalar value.
+/// x = y = z = val.
+static inline vec3 vec3_from_scalar(R val) { return (vec3){val, val, val}; }
+
+/// Normalize the vector.
+static inline vec3 vec3_normalize(vec3 v) {
+  R n = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
+  assert(n > 0);
+  return (vec3){v.x / n, v.y / n, v.z / n};
+}
+
+/// Return the vector's ith coordinate.
+static inline R vec3_ith(vec3 v, int i) {
+  assert(i >= 0 && i < 3);
+  return ((const R*)&v)[i];
+}
+
+/// Negate the given vector.
+static inline vec3 vec3_neg(vec3 v) { return (vec3){-v.x, -v.y, -v.z}; }
+
+/// Add two vectors.
+static inline vec3 vec3_add(vec3 a, vec3 b) {
+  return (vec3){a.x + b.x, a.y + b.y, a.z + b.z};
+}
+
+/// Subtract two vectors.
+static inline vec3 vec3_sub(vec3 a, vec3 b) {
+  return (vec3){a.x - b.x, a.y - b.y, a.z - b.z};
+}
+
+/// Modulate two vectors (component-wise multiplication).
+static inline vec3 vec3_mul(vec3 a, vec3 b) {
+  return (vec3){a.x * b.x, a.y * b.y, a.z * b.z};
+}
+
+/// Divide two vectors component-wise.
+static inline vec3 vec3_div(vec3 a, vec3 b) {
+  return (vec3){a.x / b.x, a.y / b.y, a.z / b.z};
+}
+
+/// Scale a vector by a scalar value.
+static inline vec3 vec3_scale(vec3 v, R s) {
+  return (vec3){v.x * s, v.y * s, v.z * s};
+}
+
+/// Compare two vectors for equality.
+static inline bool vec3_eq(vec3 a, vec3 b) {
+  return a.x == b.x && a.y == b.y && a.z == b.z;
+}
+
+/// Return the absolute value of the vector.
+static inline vec3 vec3_abs(vec3 v) {
+  return (vec3){rabs(v.x), rabs(v.y), rabs(v.z)};
+}
+
+/// Compare two vectors for inequality.
+static inline bool vec3_ne(vec3 a, vec3 b) { return !(vec3_eq(a, b)); }
+
+/// Return the vector's squared magnitude.
+static inline R vec3_norm2(vec3 v) { return v.x * v.x + v.y * v.y + v.z * v.z; }
+
+/// Return the vector's magnitude.
+static inline R vec3_norm(vec3 v) { return sqrt(vec3_norm2(v)); }
+
+/// Return the squared distance between two points.
+static inline R vec3_dist2(vec3 a, vec3 b) {
+  const vec3 v = vec3_sub(b, a);
+  return vec3_norm2(v);
+}
+
+/// Return the distance between two points.
+static inline R vec3_dist(vec3 a, vec3 b) { return sqrt(vec3_dist2(a, b)); }
+
+/// Return the given vector divided by its magnitude.
+static inline vec3 normalize(vec3 v) {
+  const R n = vec3_norm(v);
+  assert(n > 0);
+  return (vec3){v.x / n, v.y / n, v.z / n};
+}
+
+/// Return the dot product of two vectors.
+static inline R vec3_dot(vec3 a, vec3 b) {
+  return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+/// Return the cross product of two vectors.
+static inline vec3 vec3_cross(vec3 a, vec3 b) {
+  return (vec3){a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z,
+                a.x * b.y - a.y * b.x};
+}
+
+/// Reflect the vector about the normal.
+static inline vec3 vec3_reflect(vec3 v, vec3 n) {
+  // r = v - 2 * dot(v, n) * n
+  return vec3_sub(v, vec3_scale(n, 2 * vec3_dot(v, n)));
+}
+
+/// Refract the vector about the normal.
+static inline vec3 refract(vec3 v, vec3 n, R e) {
+  // k = 1 - e^2(1 - dot(n,v) * dot(n,v))
+  const R k = 1.0 - e * e * (1.0 - vec3_dot(n, v) * vec3_dot(n, v));
+  assert(k >= 0);
+  // r = e*v - (e * dot(n,v) + sqrt(k)) * n
+  return vec3_sub(vec3_scale(v, e),
+                  vec3_scale(n, e * vec3_dot(n, v) * sqrt(k)));
+}
+
+/// Elevate the vector to a power.
+static inline vec3 vec3_pow(vec3 v, R p) {
+  return (vec3){pow(v.x, p), pow(v.y, p), pow(v.z, p)};
+}
+
+/// The (1, 0, 0) vector.
+static inline vec3 right3() { return (vec3){1.0, 0.0, 0.0}; }
+
+/// The (0, 1, 0) vector.
+static inline vec3 up3() { return (const vec3){0.0, 1.0, 0.0}; }
+
+/// The (0, 0, -1) vector.
+static inline vec3 forward3() { return (const vec3){0.0, 0.0, -1.0}; }
+
+/// The (0, 0, 0) vector.
+static inline vec3 zero3() { return (const vec3){0.0, 0.0, 0.0}; }
diff --git a/include/math/vec4.h b/include/math/vec4.h
new file mode 100644
index 0000000..4ab843b
--- /dev/null
+++ b/include/math/vec4.h
@@ -0,0 +1,15 @@
+#pragma once
+
+#include "defs.h"
+
+typedef struct vec4 {
+  R x, y, w, z;
+} vec4;
+
+/// Construct a vector from 4 coordinates.
+static inline vec4 vec4_make(R x, R y, R z, R w) { return (vec4){x, y, z, w}; }
+
+/// Construct a vector from an array.
+static inline vec4 vec4_from_array(const R xyzw[4]) {
+  return (vec4){xyzw[0], xyzw[1], xyzw[2], xyzw[3]};
+}
diff --git a/test/mat4_test.c b/test/mat4_test.c
new file mode 100644
index 0000000..d4f61bc
--- /dev/null
+++ b/test/mat4_test.c
@@ -0,0 +1,111 @@
+#include <math/mat4.h>
+
+#include <math/float.h>
+
+#include "test.h"
+
+#include <stdio.h>
+
+static const float eps = 1e-9;
+
+static inline void print_mat4(mat4 m) {
+  printf("\n");
+  for (int i = 0; i < 4; ++i) {
+    for (int j = 0; j < 4; ++j) {
+      printf("%6.3f  ", mat4_at(m, i, j));
+    }
+    printf("\n");
+  }
+}
+
+TEST_CASE(lookat_right) {
+  const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
+                             /*target=*/vec3_make(1, 0, 0),
+                             /*up=*/vec3_make(0, 1, 0));
+
+  const mat4 expected = mat4_from_array((const float*)(float[4][4]){
+      {0, 0, 1, 0},
+      {0, 1, 0, 0},
+      {1, 0, 0, 0},
+      {0, 0, 0, 1},
+  });
+
+  TEST_TRUE(mat4_eq(m, expected, eps));
+}
+
+TEST_CASE(lookat_left) {
+  const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
+                             /*target=*/vec3_make(-1, 0, 0),
+                             /*up=*/vec3_make(0, 1, 0));
+
+  const mat4 expected = mat4_from_array((const float*)(float[4][4]){
+      {0, 0, -1, 0},
+      {0, 1, 0, 0},
+      {-1, 0, 0, 0},
+      {0, 0, 0, 1},
+  });
+
+  TEST_TRUE(mat4_eq(m, expected, eps));
+}
+
+TEST_CASE(lookat_up) {
+  const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
+                             /*target=*/vec3_make(0, 1, 0),
+                             /*up=*/vec3_make(0, 0, 1));
+
+  const mat4 expected = mat4_from_array((const float*)(float[4][4]){
+      {1, 0, 0, 0},
+      {0, 0, 1, 0},
+      {0, 1, 0, 0},
+      {0, 0, 0, 1},
+  });
+
+  TEST_TRUE(mat4_eq(m, expected, eps));
+}
+
+TEST_CASE(lookat_down) {
+  const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
+                             /*target=*/vec3_make(0, -1, 0),
+                             /*up=*/vec3_make(0, 0, -1));
+
+  const mat4 expected = mat4_from_array((const float*)(float[4][4]){
+      {1, 0, 0, 0},
+      {0, 0, -1, 0},
+      {0, -1, 0, 0},
+      {0, 0, 0, 1},
+  });
+
+  TEST_TRUE(mat4_eq(m, expected, eps));
+}
+
+TEST_CASE(lookat_forward) {
+  const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
+                             /*target=*/vec3_make(0, 0, -1),
+                             /*up=*/vec3_make(0, 1, 0));
+
+  const mat4 expected = mat4_from_array((const float*)(float[4][4]){
+      {1, 0, 0, 0},
+      {0, 1, 0, 0},
+      {0, 0, -1, 0},
+      {0, 0, 0, 1},
+  });
+
+  TEST_TRUE(mat4_eq(m, expected, eps));
+}
+
+TEST_CASE(lookat_back) {
+  const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
+                             /*target=*/vec3_make(0, 0, 1),
+                             /*up=*/vec3_make(0, 1, 0));
+
+  const mat4 expected = mat4_from_array((const float*)(float[4][4]){
+      {-1, 0, 0, 0},
+      {0, 1, 0, 0},
+      {0, 0, 1, 0},
+      {0, 0, 0, 1},
+  });
+
+  TEST_TRUE(mat4_eq(m, expected, eps));
+}
+
+int main() { return 0; }
diff --git a/test/test.h b/test/test.h
new file mode 100644
index 0000000..fd8dc22
--- /dev/null
+++ b/test/test.h
@@ -0,0 +1,185 @@
+// SPDX-License-Identifier: MIT
+#pragma once
+
+#ifdef UNIT_TEST
+
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#if defined(__DragonFly__) || defined(__FreeBSD__) || defined(__FreeBSD_kernel__) ||     \
+    defined(__NetBSD__) || defined(__OpenBSD__)
+#define USE_SYSCTL_FOR_ARGS 1
+// clang-format off
+#include <sys/types.h>
+#include <sys/sysctl.h>
+// clang-format on
+#include <unistd.h>        // getpid
+#endif
+
+struct test_file_metadata;
+
+struct test_failure {
+        bool present;
+        const char *message;
+        const char *file;
+        int line;
+};
+
+struct test_case_metadata {
+        void (*fn)(struct test_case_metadata *, struct test_file_metadata *);
+        struct test_failure failure;
+        const char *name;
+        struct test_case_metadata *next;
+};
+
+struct test_file_metadata {
+        bool registered;
+        const char *name;
+        struct test_file_metadata *next;
+        struct test_case_metadata *tests;
+};
+
+struct test_file_metadata __attribute__((weak)) * test_file_head;
+
+#define SET_FAILURE(_message)                                                             \
+        metadata->failure = (struct test_failure) {                                       \
+                .message = _message, .file = __FILE__, .line = __LINE__, .present = true, \
+        }
+
+#define TEST_EQUAL(a, b)                                                                 \
+        do {                                                                             \
+                if ((a) != (b)) {                                                        \
+                        SET_FAILURE(#a " != " #b);                                       \
+                        return;                                                          \
+                }                                                                        \
+        } while (0)
+
+#define TEST_TRUE(a)                                                                     \
+        do {                                                                             \
+                if (!(a)) {                                                              \
+                        SET_FAILURE(#a " is not true");                                  \
+                        return;                                                          \
+                }                                                                        \
+        } while (0)
+
+#define TEST_STREQUAL(a, b)                                                              \
+        do {                                                                             \
+                if (strcmp(a, b) != 0) {                                                 \
+                        SET_FAILURE(#a " != " #b);                                       \
+                        return;                                                          \
+                }                                                                        \
+        } while (0)
+
+#define TEST_CASE(_name)                                                                  \
+        static void __test_h_##_name(struct test_case_metadata *,                         \
+                                     struct test_file_metadata *);                        \
+        static struct test_file_metadata __test_h_file;                                   \
+        static struct test_case_metadata __test_h_meta_##_name = {                        \
+            .name = #_name,                                                               \
+            .fn = __test_h_##_name,                                                       \
+        };                                                                                \
+        static void __attribute__((constructor(101))) __test_h_##_name##_register(void) { \
+                __test_h_meta_##_name.next = __test_h_file.tests;                         \
+                __test_h_file.tests = &__test_h_meta_##_name;                             \
+                if (!__test_h_file.registered) {                                          \
+                        __test_h_file.name = __FILE__;                                    \
+                        __test_h_file.next = test_file_head;                              \
+                        test_file_head = &__test_h_file;                                  \
+                        __test_h_file.registered = true;                                  \
+                }                                                                         \
+        }                                                                                 \
+        static void __test_h_##_name(                                                     \
+            struct test_case_metadata *metadata __attribute__((unused)),                  \
+            struct test_file_metadata *file_metadata __attribute__((unused)))
+
+extern void __attribute__((weak)) (*test_h_unittest_setup)(void);
+/// Run defined tests, return true if all tests succeeds
+/// @param[out] tests_run if not NULL, set to whether tests were run
+static inline void __attribute__((constructor(102))) run_tests(void) {
+        bool should_run = false;
+#ifdef USE_SYSCTL_FOR_ARGS
+        int mib[] = {
+                CTL_KERN,
+#if defined(__NetBSD__) || defined(__OpenBSD__)
+                KERN_PROC_ARGS,
+                getpid(),
+                KERN_PROC_ARGV,
+#else
+                KERN_PROC,
+                KERN_PROC_ARGS,
+                getpid(),
+#endif
+        };
+        char *arg = NULL;
+        size_t arglen;
+        sysctl(mib, sizeof(mib) / sizeof(mib[0]), NULL, &arglen, NULL, 0);
+        arg = malloc(arglen);
+        sysctl(mib, sizeof(mib) / sizeof(mib[0]), arg, &arglen, NULL, 0);
+#else
+        FILE *cmdlinef = fopen("/proc/self/cmdline", "r");
+        char *arg = NULL;
+        int arglen;
+        fscanf(cmdlinef, "%ms%n", &arg, &arglen);
+        fclose(cmdlinef);
+#endif
+        for (char *pos = arg; pos < arg + arglen; pos += strlen(pos) + 1) {
+                if (strcmp(pos, "--unittest") == 0) {
+                        should_run = true;
+                        break;
+                }
+        }
+        free(arg);
+
+        if (!should_run) {
+                return;
+        }
+
+        if (&test_h_unittest_setup) {
+                test_h_unittest_setup();
+        }
+
+        struct test_file_metadata *i = test_file_head;
+        int failed = 0, success = 0;
+        while (i) {
+                fprintf(stderr, "Running tests from %s:\n", i->name);
+                struct test_case_metadata *j = i->tests;
+                while (j) {
+                        fprintf(stderr, "\t%s ... ", j->name);
+                        j->failure.present = false;
+                        j->fn(j, i);
+                        if (j->failure.present) {
+                                fprintf(stderr, "failed (%s at %s:%d)\n", j->failure.message,
+                                        j->failure.file, j->failure.line);
+                                failed++;
+                        } else {
+                                fprintf(stderr, "passed\n");
+                                success++;
+                        }
+                        j = j->next;
+                }
+                fprintf(stderr, "\n");
+                i = i->next;
+        }
+        int total = failed + success;
+        fprintf(stderr, "Test results: passed %d/%d, failed %d/%d\n", success, total,
+                failed, total);
+        exit(failed == 0 ? EXIT_SUCCESS : EXIT_FAILURE);
+}
+
+#else
+
+#include <stdbool.h>
+
+#define TEST_CASE(name) static void __attribute__((unused)) __test_h_##name(void)
+
+#define TEST_EQUAL(a, b)                                                                 \
+        (void)(a);                                                                       \
+        (void)(b)
+#define TEST_TRUE(a) (void)(a)
+#define TEST_STREQUAL(a, b)                                                              \
+        (void)(a);                                                                       \
+        (void)(b)
+
+#endif