math: Make vectors have an array instead of different fields.

25 files changed, 391 insertions, 382 deletions

diff --git a/libdimension-python/dimension.pxd b/libdimension-python/dimension.pxd
index cb152be..4edfbf2 100644
--- a/libdimension-python/dimension.pxd
+++ b/libdimension-python/dimension.pxd
@@ -115,9 +115,7 @@ cdef extern from "../libdimension/dimension.h":
   double dmnsn_degrees(double radians)
 
   ctypedef struct dmnsn_vector:
-    double x
-    double y
-    double z
+    double n[3]
 
   dmnsn_vector dmnsn_new_vector(double x, double y, double z)
 
@@ -126,7 +124,6 @@ cdef extern from "../libdimension/dimension.h":
   dmnsn_vector dmnsn_vector_add(dmnsn_vector lhs, dmnsn_vector rhs)
   dmnsn_vector dmnsn_vector_sub(dmnsn_vector lhs, dmnsn_vector rhs)
   dmnsn_vector dmnsn_vector_mul(double lhs, dmnsn_vector rhs)
-  dmnsn_vector dmnsn_vector_div(dmnsn_vector lhs, double rhs)
 
   dmnsn_vector dmnsn_vector_cross(dmnsn_vector lhs, dmnsn_vector rhs)
   double dmnsn_vector_dot(dmnsn_vector lhs, dmnsn_vector rhs)
diff --git a/libdimension-python/dimension.pyx b/libdimension-python/dimension.pyx
index 34d720e..77c9cd8 100644
--- a/libdimension-python/dimension.pyx
+++ b/libdimension-python/dimension.pyx
@@ -234,15 +234,15 @@ cdef class Vector:
   property x:
     """The x coordinate."""
     def __get__(self):
-      return self._v.x
+      return self._v.n[0]
   property y:
     """The y coordinate."""
     def __get__(self):
-      return self._v.y
+      return self._v.n[1]
   property z:
     """The z coordinate."""
     def __get__(self):
-      return self._v.z
+      return self._v.n[2]
 
   def __pos__(self):
     return self
@@ -261,7 +261,7 @@ cdef class Vector:
     else:
       return _Vector(dmnsn_vector_mul(lhs, (<Vector?>rhs)._v))
   def __truediv__(Vector lhs not None, double rhs):
-    return _Vector(dmnsn_vector_div(lhs._v, rhs))
+    return _Vector(dmnsn_vector_mul(1.0/rhs, lhs._v))
 
   def __richcmp__(lhs, rhs, int op):
     equal = (Vector(lhs) - Vector(rhs)).norm() < dmnsn_epsilon
diff --git a/libdimension/bench/geometry.c b/libdimension/bench/geometry.c
index 59d27e3..ae30d17 100644
--- a/libdimension/bench/geometry.c
+++ b/libdimension/bench/geometry.c
@@ -107,12 +107,6 @@ main(void)
   });
   printf("dmnsn_vector_mul(): %ld\n", sandglass.grains);
 
-  // dmnsn_vector_div()
-  sandglass_bench_fine(&sandglass, {
-    vector = dmnsn_vector_div(vector, 2.0);
-  });
-  printf("dmnsn_vector_div(): %ld\n", sandglass.grains);
-
   // dmnsn_vector_cross()
   sandglass_bench_fine(&sandglass, {
     vector = dmnsn_vector_cross(vector, vector2);
diff --git a/libdimension/bench/prtree.c b/libdimension/bench/prtree.c
index 6b5e7c1..aa29317 100644
--- a/libdimension/bench/prtree.c
+++ b/libdimension/bench/prtree.c
@@ -29,7 +29,7 @@ static bool
 dmnsn_fake_intersection_fn(const dmnsn_object *object, dmnsn_ray ray,
                            dmnsn_intersection *intersection)
 {
-  intersection->t = (object->aabb.min.z - ray.x0.z)/ray.n.z;
+  intersection->t = (object->aabb.min.Z - ray.x0.Z)/ray.n.Z;
   intersection->normal = dmnsn_x;
   return true;
 }
@@ -45,13 +45,10 @@ dmnsn_fake_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
 {
   dmnsn_vector a, b;
 
-  a.x = 2.0*((double)rand())/RAND_MAX - 1.0;
-  a.y = 2.0*((double)rand())/RAND_MAX - 1.0;
-  a.z = 2.0*((double)rand())/RAND_MAX - 1.0;
-
-  b.x = 2.0*((double)rand())/RAND_MAX - 1.0;
-  b.y = 2.0*((double)rand())/RAND_MAX - 1.0;
-  b.z = 2.0*((double)rand())/RAND_MAX - 1.0;
+  for (unsigned int i = 0; i < 3; ++i) {
+    a.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0;
+    b.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0;
+  }
 
   return dmnsn_new_aabb(dmnsn_vector_min(a, b), dmnsn_vector_max(a, b));
 }
diff --git a/libdimension/bvh/bvh.c b/libdimension/bvh/bvh.c
index eab2c28..9f460a2 100644
--- a/libdimension/bvh/bvh.c
+++ b/libdimension/bvh/bvh.c
@@ -186,7 +186,7 @@ dmnsn_optimize_ray(dmnsn_ray ray)
 {
   dmnsn_optimized_ray optray = {
     .x0    = ray.x0,
-    .n_inv = dmnsn_new_vector(1.0/ray.n.x, 1.0/ray.n.y, 1.0/ray.n.z)
+    .n_inv = dmnsn_new_vector(1.0/ray.n.X, 1.0/ray.n.Y, 1.0/ray.n.Z)
   };
   return optray;
 }
@@ -202,20 +202,20 @@ dmnsn_ray_box_intersection(dmnsn_optimized_ray optray, dmnsn_aabb box, double t)
   // or tmax == -inf, while rays inside the box will have tmin and tmax
   // unchanged.
 
-  double tx1 = (box.min.x - optray.x0.x)*optray.n_inv.x;
-  double tx2 = (box.max.x - optray.x0.x)*optray.n_inv.x;
+  double tx1 = (box.min.X - optray.x0.X)*optray.n_inv.X;
+  double tx2 = (box.max.X - optray.x0.X)*optray.n_inv.X;
 
   double tmin = dmnsn_min(tx1, tx2);
   double tmax = dmnsn_max(tx1, tx2);
 
-  double ty1 = (box.min.y - optray.x0.y)*optray.n_inv.y;
-  double ty2 = (box.max.y - optray.x0.y)*optray.n_inv.y;
+  double ty1 = (box.min.Y - optray.x0.Y)*optray.n_inv.Y;
+  double ty2 = (box.max.Y - optray.x0.Y)*optray.n_inv.Y;
 
   tmin = dmnsn_max(tmin, dmnsn_min(ty1, ty2));
   tmax = dmnsn_min(tmax, dmnsn_max(ty1, ty2));
 
-  double tz1 = (box.min.z - optray.x0.z)*optray.n_inv.z;
-  double tz2 = (box.max.z - optray.x0.z)*optray.n_inv.z;
+  double tz1 = (box.min.Z - optray.x0.Z)*optray.n_inv.Z;
+  double tz2 = (box.max.Z - optray.x0.Z)*optray.n_inv.Z;
 
   tmin = dmnsn_max(tmin, dmnsn_min(tz1, tz2));
   tmax = dmnsn_min(tmax, dmnsn_max(tz1, tz2));
diff --git a/libdimension/bvh/prtree.c b/libdimension/bvh/prtree.c
index c8e4e54..84fca58 100644
--- a/libdimension/bvh/prtree.c
+++ b/libdimension/bvh/prtree.c
@@ -74,48 +74,48 @@ enum {
 static int
 dmnsn_xmin_comp(const void *l, const void *r)
 {
-  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.x;
-  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.x;
+  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.X;
+  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.X;
   return (lval > rval) - (lval < rval);
 }
 
 static int
 dmnsn_ymin_comp(const void *l, const void *r)
 {
-  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.y;
-  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.y;
+  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.Y;
+  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.Y;
   return (lval > rval) - (lval < rval);
 }
 
 static int
 dmnsn_zmin_comp(const void *l, const void *r)
 {
-  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.z;
-  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.z;
+  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.Z;
+  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.Z;
   return (lval > rval) - (lval < rval);
 }
 
 static int
 dmnsn_xmax_comp(const void *l, const void *r)
 {
-  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.x;
-  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.x;
+  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.X;
+  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.X;
   return (lval < rval) - (lval > rval);
 }
 
 static int
 dmnsn_ymax_comp(const void *l, const void *r)
 {
-  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.y;
-  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.y;
+  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.Y;
+  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.Y;
   return (lval < rval) - (lval > rval);
 }
 
 static int
 dmnsn_zmax_comp(const void *l, const void *r)
 {
-  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.z;
-  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.z;
+  double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.Z;
+  double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.Z;
   return (lval < rval) - (lval > rval);
 }
 
diff --git a/libdimension/dimension.h b/libdimension/dimension.h
index db67a48..2489669 100644
--- a/libdimension/dimension.h
+++ b/libdimension/dimension.h
@@ -34,9 +34,6 @@
  * all rendering-related tasks for Dimension.
  */
 
-#ifndef DMNSN_H
-#define DMNSN_H
-
 /* Include all modules. */
 #include <dimension/base.h>
 #include <dimension/platform.h>
@@ -47,5 +44,3 @@
 #include <dimension/pattern.h>
 #include <dimension/model.h>
 #include <dimension/render.h>
-
-#endif /* DMNSN_H */
diff --git a/libdimension/dimension/math.h b/libdimension/dimension/math.h
index 603373f..e110ffb 100644
--- a/libdimension/dimension/math.h
+++ b/libdimension/dimension/math.h
@@ -43,3 +43,14 @@ extern "C" {
 #endif
 
 #endif /* DMNSN_MATH_H */
+
+#if defined(DMNSN_SHORT_NAMES) && !defined(DMNSN_SHORT_NAMES_DEFINED)
+  #define DMNSN_SHORT_NAMES_DEFINED
+
+  /** Short name for the \a x component of a \c dmnsn_vector. */
+  #define X n[0]
+  /** Short name for the \a y component of a \c dmnsn_vector. */
+  #define Y n[1]
+  /** Short name for the \a z component of a \c dmnsn_vector. */
+  #define Z n[2]
+#endif
diff --git a/libdimension/dimension/math/aabb.h b/libdimension/dimension/math/aabb.h
index 14cc575..4eb7870 100644
--- a/libdimension/dimension/math/aabb.h
+++ b/libdimension/dimension/math/aabb.h
@@ -57,8 +57,8 @@ DMNSN_INLINE dmnsn_aabb
 dmnsn_zero_aabb(void)
 {
   dmnsn_aabb box = {
-    {  DMNSN_INFINITY,  DMNSN_INFINITY,  DMNSN_INFINITY },
-    { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY }
+    { {  DMNSN_INFINITY,  DMNSN_INFINITY,  DMNSN_INFINITY } },
+    { { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY } }
   };
   return box;
 }
@@ -68,8 +68,8 @@ DMNSN_INLINE dmnsn_aabb
 dmnsn_infinite_aabb(void)
 {
   dmnsn_aabb box = {
-    { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY },
-    {  DMNSN_INFINITY,  DMNSN_INFINITY,  DMNSN_INFINITY }
+    { { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY } },
+    { {  DMNSN_INFINITY,  DMNSN_INFINITY,  DMNSN_INFINITY } }
   };
   return box;
 }
@@ -94,15 +94,26 @@ dmnsn_symmetric_aabb(dmnsn_vector r)
 DMNSN_INLINE bool
 dmnsn_aabb_contains(dmnsn_aabb box, dmnsn_vector p)
 {
-  return (p.x >= box.min.x && p.y >= box.min.y && p.z >= box.min.z)
-      && (p.x <= box.max.x && p.y <= box.max.y && p.z <= box.max.z);
+  for (unsigned int i = 0; i < 3; ++i) {
+    if (p.n[i] < box.min.n[i]) {
+      return false;
+    }
+  }
+
+  for (unsigned int i = 0; i < 3; ++i) {
+    if (p.n[i] > box.max.n[i]) {
+      return false;
+    }
+  }
+
+  return true;
 }
 
 /** Return whether a bounding box is infinite. */
 DMNSN_INLINE bool
 dmnsn_aabb_is_infinite(dmnsn_aabb box)
 {
-  return box.min.x == -DMNSN_INFINITY;
+  return box.min.n[0] == -DMNSN_INFINITY;
 }
 
 /**
diff --git a/libdimension/dimension/math/matrix.h b/libdimension/dimension/math/matrix.h
index 7471bf5..e67121e 100644
--- a/libdimension/dimension/math/matrix.h
+++ b/libdimension/dimension/math/matrix.h
@@ -39,152 +39,218 @@ typedef struct dmnsn_matrix {
   "[%g\t%g\t%g\t%g]\n"                          \
   "[%g\t%g\t%g\t%g]"
 /** The appropriate arguements to printf() a matrix. */
-#define DMNSN_MATRIX_PRINTF(m)                          \
-  (m).n[0][0], (m).n[0][1], (m).n[0][2], (m).n[0][3],   \
-  (m).n[1][0], (m).n[1][1], (m).n[1][2], (m).n[1][3],   \
-  (m).n[2][0], (m).n[2][1], (m).n[2][2], (m).n[2][3],   \
+#define DMNSN_MATRIX_PRINTF(M)                              \
+  (M).n[0][0], (M).n[0][1], (M).n[0][2], (M).n[0][3],   \
+  (M).n[1][0], (M).n[1][1], (M).n[1][2], (M).n[1][3],   \
+  (M).n[2][0], (M).n[2][1], (M).n[2][2], (M).n[2][3],   \
   0.0, 0.0, 0.0, 1.0
 
-/** Construct a new transformation matrix. */
+/** Create a transformation matrix. */
 DMNSN_INLINE dmnsn_matrix
-dmnsn_new_matrix(double a0, double a1, double a2, double a3,
-                 double b0, double b1, double b2, double b3,
-                 double c0, double c1, double c2, double c3)
+dmnsn_new_matrix(double a0, double b0, double c0, double d0,
+                 double a1, double b1, double c1, double d1,
+                 double a2, double b2, double c2, double d2)
 {
-  dmnsn_matrix m = { { { a0, a1, a2, a3 },
-                       { b0, b1, b2, b3 },
-                       { c0, c1, c2, c3 } } };
-  return m;
+  dmnsn_matrix M = {
+    {
+      { a0, b0, c0, d0 },
+      { a1, b1, c1, d1 },
+      { a2, b2, c2, d2 }
+    }
+  };
+  return M;
 }
 
-/** Construct a new transformation matrix from column vectors. */
+/** Create a transformation matrix from column vectors. */
 DMNSN_INLINE dmnsn_matrix
-dmnsn_new_matrix4(dmnsn_vector a, dmnsn_vector b, dmnsn_vector c,
-                  dmnsn_vector d)
+dmnsn_new_matrix4(dmnsn_vector a, dmnsn_vector b, dmnsn_vector c, dmnsn_vector d)
 {
-  dmnsn_matrix m = { { { a.x, b.x, c.x, d.x },
-                       { a.y, b.y, c.y, d.y },
-                       { a.z, b.z, c.z, d.z } } };
-  return m;
+  dmnsn_matrix M;
+
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    M.n[i][0] = a.n[i];
+  }
+  for (i = 0; i < 3; ++i) {
+    M.n[i][1] = b.n[i];
+  }
+  for (i = 0; i < 3; ++i) {
+    M.n[i][2] = c.n[i];
+  }
+  for (i = 0; i < 3; ++i) {
+    M.n[i][3] = d.n[i];
+  }
+
+  return M;
 }
 
 /** Extract column vectors from a matrix. */
 DMNSN_INLINE dmnsn_vector
 dmnsn_matrix_column(dmnsn_matrix M, unsigned int i)
 {
-  return dmnsn_new_vector(M.n[0][i], M.n[1][i], M.n[2][i]);
+  dmnsn_vector v;
+  unsigned int j;
+  for (j = 0; j < 3; ++j) {
+    v.n[j] = M.n[j][i];
+  }
+  return v;
 }
 
 /** Return the identity matrix. */
-dmnsn_matrix dmnsn_identity_matrix(void);
+DMNSN_INLINE dmnsn_matrix
+dmnsn_identity_matrix(void)
+{
+  return dmnsn_new_matrix(
+    1.0, 0.0, 0.0, 0.0,
+    0.0, 1.0, 0.0, 0.0,
+    0.0, 0.0, 1.0, 0.0
+  );
+}
 
 /**
- * A scale transformation.
- * @param[in] s  A vector with components representing the scaling factor in
- *               each axis.
- * @return The transformation matrix.
+ * Return a scale transformation.
+ * @param[in]  s  A vector with components representing the scaling factor in
+ *                each axis.
  */
-dmnsn_matrix dmnsn_scale_matrix(dmnsn_vector s);
+DMNSN_INLINE dmnsn_matrix
+dmnsn_scale_matrix(dmnsn_vector s)
+{
+  return dmnsn_new_matrix(
+    s.n[0], 0.0,    0.0,    0.0,
+    0.0,    s.n[1], 0.0,    0.0,
+    0.0,    0.0,    s.n[2], 0.0
+  );
+}
+
 /**
- * A translation.
- * @param[in] d  The vector to translate by.
- * @return The transformation matrix.
+ * Set \p M to a translation matrix.
+ * @param[in]  d  The vector to translate by.
  */
-dmnsn_matrix dmnsn_translation_matrix(dmnsn_vector d);
+DMNSN_INLINE dmnsn_matrix
+dmnsn_translation_matrix(dmnsn_vector d)
+{
+  return dmnsn_new_matrix(
+    1.0, 0.0, 0.0, d.n[0],
+    0.0, 1.0, 0.0, d.n[1],
+    0.0, 0.0, 1.0, d.n[2]
+  );
+}
+
 /**
- * A left-handed rotation.
- * @param[in] theta  A vector representing an axis and angle.
- *                   @f$ axis = \vec{\theta}/|\vec{\theta}| @f$,
- *                   @f$ angle = |\vec{\theta}| @f$
- * @return The transformation matrix.
+ * Return a rotation matrix.
+ * @param[in]  theta  A vector representing an axis and angle.
+ *                    @f$ axis = \vec{\theta}/|\vec{\theta}| @f$,
+ *                    @f$ angle = |\vec{\theta}| @f$
  */
 dmnsn_matrix dmnsn_rotation_matrix(dmnsn_vector theta);
 
 /**
- * An alignment matrix.
+ * Return an alignment matrix.
  * @param[in] from   The initial vector.
  * @param[in] to     The desired direction.
  * @param[in] axis1  The first axis about which to rotate.
  * @param[in] axis2  The second axis about which to rotate.
- * @return A transformation matrix that will rotate \p from to \p to.
  */
-dmnsn_matrix dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to,
-                                    dmnsn_vector axis1, dmnsn_vector axis2);
+dmnsn_matrix dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis1, dmnsn_vector axis2);
 
 /** Invert a matrix. */
-dmnsn_matrix dmnsn_matrix_inverse(dmnsn_matrix A);
+dmnsn_matrix dmnsn_matrix_inverse(dmnsn_matrix M);
 
 /** Multiply two matricies. */
-dmnsn_matrix dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs);
+DMNSN_INLINE dmnsn_matrix dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs)
+{
+  dmnsn_matrix M;
+
+  unsigned int i, j, k;
+  for (i = 0; i < 3; ++i) {
+    for (j = 0; j < 3; ++j) {
+      M.n[i][j] = 0.0;
+      for (k = 0; k < 3; ++k) {
+        M.n[i][j] += lhs.n[i][k]*rhs.n[k][j];
+      }
+    }
+
+    M.n[i][3] = lhs.n[i][3];
+    for (k = 0; k < 3; ++k) {
+      M.n[i][3] += lhs.n[i][k]*rhs.n[k][3];
+    }
+  }
+
+  return M;
+}
 
 /** Transform a point by a matrix. */
 DMNSN_INLINE dmnsn_vector
-dmnsn_transform_point(dmnsn_matrix T, dmnsn_vector v)
+dmnsn_transform_point(dmnsn_matrix M, dmnsn_vector v)
 {
-  /* 9 multiplications, 9 additions */
   dmnsn_vector r;
-  r.x = T.n[0][0]*v.x + T.n[0][1]*v.y + T.n[0][2]*v.z + T.n[0][3];
-  r.y = T.n[1][0]*v.x + T.n[1][1]*v.y + T.n[1][2]*v.z + T.n[1][3];
-  r.z = T.n[2][0]*v.x + T.n[2][1]*v.y + T.n[2][2]*v.z + T.n[2][3];
+  unsigned int i, j;
+  for (i = 0; i < 3; ++i) {
+    r.n[i] = M.n[i][3];
+    for (j = 0; j < 3; ++j) {
+      r.n[i] += M.n[i][j]*v.n[j];
+    }
+  }
   return r;
 }
 
 /** Transform a direction by a matrix. */
 DMNSN_INLINE dmnsn_vector
-dmnsn_transform_direction(dmnsn_matrix T, dmnsn_vector v)
+dmnsn_transform_direction(dmnsn_matrix M, dmnsn_vector v)
 {
-  /* 9 multiplications, 6 additions */
   dmnsn_vector r;
-  r.x = T.n[0][0]*v.x + T.n[0][1]*v.y + T.n[0][2]*v.z;
-  r.y = T.n[1][0]*v.x + T.n[1][1]*v.y + T.n[1][2]*v.z;
-  r.z = T.n[2][0]*v.x + T.n[2][1]*v.y + T.n[2][2]*v.z;
+  unsigned int i, j;
+  for (i = 0; i < 3; ++i) {
+    r.n[i] = 0.0;
+    for (j = 0; j < 3; ++j) {
+      r.n[i] += M.n[i][j]*v.n[j];
+    }
+  }
   return r;
 }
 
 /**
  * Transform a pseudovector by a matrix.
- * @param[in] Tinv  The inverse of the transformation matrix.
- * @param[in] v     The pseudovector to transform
+ * @param[in] Minv  The inverse of the transformation matrix.
+ * @param[in] v     The pseudovector to transform.
  * @return The transformed pseudovector.
  */
 DMNSN_INLINE dmnsn_vector
-dmnsn_transform_normal(dmnsn_matrix Tinv, dmnsn_vector v)
+dmnsn_transform_normal(dmnsn_matrix Minv, dmnsn_vector v)
 {
-  /* Multiply by the transpose of the inverse
-     (9 multiplications, 6 additions) */
+  /* Multiply by the transpose of the inverse */
   dmnsn_vector r;
-  r.x = Tinv.n[0][0]*v.x + Tinv.n[1][0]*v.y + Tinv.n[2][0]*v.z;
-  r.y = Tinv.n[0][1]*v.x + Tinv.n[1][1]*v.y + Tinv.n[2][1]*v.z;
-  r.z = Tinv.n[0][2]*v.x + Tinv.n[1][2]*v.y + Tinv.n[2][2]*v.z;
+  unsigned int i, j;
+  for (i = 0; i < 3; ++i) {
+    r.n[i] = 0.0;
+    for (j = 0; j < 3; ++j) {
+      r.n[i] += Minv.n[j][i]*v.n[j];
+    }
+  }
   return r;
 }
 
-/**
- * Transform a ray by a matrix.
- * \f$ n' = T(l.\vec{x_0} + l.\vec{n}) - T(l.\vec{x_0}) \f$,
- * \f$ \vec{x_0}' = T(l.\vec{x_0}) \f$
- */
+/** Transform a ray by a matrix. */
 DMNSN_INLINE dmnsn_ray
-dmnsn_transform_ray(dmnsn_matrix T, dmnsn_ray l)
+dmnsn_transform_ray(dmnsn_matrix M, dmnsn_ray l)
 {
-  /* 18 multiplications, 15 additions */
   dmnsn_ray ret;
-  ret.x0 = dmnsn_transform_point(T, l.x0);
-  ret.n  = dmnsn_transform_direction(T, l.n);
+  ret.x0 = dmnsn_transform_point(M, l.x0);
+  ret.n = dmnsn_transform_direction(M, l.n);
   return ret;
 }
 
 /** Transform a bounding box by a matrix. */
-dmnsn_aabb dmnsn_transform_aabb(dmnsn_matrix T, dmnsn_aabb box);
+dmnsn_aabb dmnsn_transform_aabb(dmnsn_matrix M, dmnsn_aabb box);
 
 /** Return whether a matrix contains any NaN components. */
 DMNSN_INLINE bool
-dmnsn_matrix_isnan(dmnsn_matrix m)
+dmnsn_matrix_isnan(dmnsn_matrix M)
 {
-  size_t i, j;
+  unsigned int i, j;
   for (i = 0; i < 3; ++i) {
     for (j = 0; j < 4; ++j) {
-      if (dmnsn_isnan(m.n[i][j])) {
+      if (dmnsn_isnan(M.n[i][j])) {
         return true;
       }
     }
diff --git a/libdimension/dimension/math/vector.h b/libdimension/dimension/math/vector.h
index 8eacee9..90b9a3d 100644
--- a/libdimension/dimension/math/vector.h
+++ b/libdimension/dimension/math/vector.h
@@ -32,26 +32,24 @@
 
 /** A vector in 3 dimensions. */
 typedef struct dmnsn_vector {
-  double x; /**< The x component. */
-  double y; /**< The y component. */
-  double z; /**< The z component. */
+  double n[3]; /**< The components. */
 } dmnsn_vector;
 
 /** A standard format string for vectors. */
 #define DMNSN_VECTOR_FORMAT "<%g, %g, %g>"
 /** The appropriate arguements to printf() a vector. */
-#define DMNSN_VECTOR_PRINTF(v) (v).x, (v).y, (v).z
+#define DMNSN_VECTOR_PRINTF(v) (v).n[0], (v).n[1], (v).n[2]
 
 /* Constants */
 
 /** The zero vector. */
-static const dmnsn_vector dmnsn_zero = { 0.0, 0.0, 0.0 };
+static const dmnsn_vector dmnsn_zero = { { 0.0, 0.0, 0.0 } };
 /** The x vector. */
-static const dmnsn_vector dmnsn_x = { 1.0, 0.0, 0.0 };
+static const dmnsn_vector dmnsn_x = { { 1.0, 0.0, 0.0 } };
 /** The y vector. */
-static const dmnsn_vector dmnsn_y = { 0.0, 1.0, 0.0 };
+static const dmnsn_vector dmnsn_y = { { 0.0, 1.0, 0.0 } };
 /** The z vector. */
-static const dmnsn_vector dmnsn_z = { 0.0, 0.0, 1.0 };
+static const dmnsn_vector dmnsn_z = { { 0.0, 0.0, 1.0 } };
 
 /* Shorthand for vector construction */
 
@@ -59,7 +57,7 @@ static const dmnsn_vector dmnsn_z = { 0.0, 0.0, 1.0 };
 DMNSN_INLINE dmnsn_vector
 dmnsn_new_vector(double x, double y, double z)
 {
-  dmnsn_vector v = { x, y, z };
+  dmnsn_vector v = { { x, y, z } };
   return v;
 }
 
@@ -69,8 +67,11 @@ dmnsn_new_vector(double x, double y, double z)
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_negate(dmnsn_vector rhs)
 {
-  /* 3 negations */
-  dmnsn_vector v = { -rhs.x, -rhs.y, -rhs.z };
+  dmnsn_vector v;
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    v.n[i] = -rhs.n[i];
+  }
   return v;
 }
 
@@ -78,8 +79,11 @@ dmnsn_vector_negate(dmnsn_vector rhs)
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_add(dmnsn_vector lhs, dmnsn_vector rhs)
 {
-  /* 3 additions */
-  dmnsn_vector v = { lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z };
+  dmnsn_vector v;
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    v.n[i] = lhs.n[i] + rhs.n[i];
+  }
   return v;
 }
 
@@ -87,8 +91,11 @@ dmnsn_vector_add(dmnsn_vector lhs, dmnsn_vector rhs)
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_sub(dmnsn_vector lhs, dmnsn_vector rhs)
 {
-  /* 3 additions */
-  dmnsn_vector v = { lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z };
+  dmnsn_vector v;
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    v.n[i] = lhs.n[i] - rhs.n[i];
+  }
   return v;
 }
 
@@ -96,17 +103,11 @@ dmnsn_vector_sub(dmnsn_vector lhs, dmnsn_vector rhs)
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_mul(double lhs, dmnsn_vector rhs)
 {
-  /* 3 multiplications */
-  dmnsn_vector v = { lhs*rhs.x, lhs*rhs.y, lhs*rhs.z };
-  return v;
-}
-
-/** Divide a vector by a scalar. */
-DMNSN_INLINE dmnsn_vector
-dmnsn_vector_div(dmnsn_vector lhs, double rhs)
-{
-  /* 3 divisions */
-  dmnsn_vector v = { lhs.x/rhs, lhs.y/rhs, lhs.z/rhs };
+  dmnsn_vector v;
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    v.n[i] = lhs*rhs.n[i];
+  }
   return v;
 }
 
@@ -114,18 +115,22 @@ dmnsn_vector_div(dmnsn_vector lhs, double rhs)
 DMNSN_INLINE double
 dmnsn_vector_dot(dmnsn_vector lhs, dmnsn_vector rhs)
 {
-  /* 3 multiplications, 2 additions */
-  return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
+  double result = 0.0;
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    result += lhs.n[i]*rhs.n[i];
+  }
+  return result;
 }
 
 /** Return the cross product of two vectors. */
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_cross(dmnsn_vector lhs, dmnsn_vector rhs)
 {
-  /* 6 multiplications, 3 additions */
-  dmnsn_vector v = { lhs.y*rhs.z - lhs.z*rhs.y,
-                     lhs.z*rhs.x - lhs.x*rhs.z,
-                     lhs.x*rhs.y - lhs.y*rhs.x };
+  dmnsn_vector v;
+  v.n[0] = lhs.n[1]*rhs.n[2] - lhs.n[2]*rhs.n[1];
+  v.n[1] = lhs.n[2]*rhs.n[0] - lhs.n[0]*rhs.n[2];
+  v.n[2] = lhs.n[0]*rhs.n[1] - lhs.n[1]*rhs.n[0];
   return v;
 }
 
@@ -133,7 +138,6 @@ dmnsn_vector_cross(dmnsn_vector lhs, dmnsn_vector rhs)
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_proj(dmnsn_vector u, dmnsn_vector d)
 {
-  /* 1 division, 9 multiplications, 4 additions */
   return dmnsn_vector_mul(dmnsn_vector_dot(u, d)/dmnsn_vector_dot(d, d), d);
 }
 
@@ -141,7 +145,6 @@ dmnsn_vector_proj(dmnsn_vector u, dmnsn_vector d)
 DMNSN_INLINE double
 dmnsn_vector_norm(dmnsn_vector n)
 {
-  /* 1 sqrt, 3 multiplications, 2 additions */
   return sqrt(dmnsn_vector_dot(n, n));
 }
 
@@ -149,35 +152,42 @@ dmnsn_vector_norm(dmnsn_vector n)
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_normalized(dmnsn_vector n)
 {
-  /* 1 sqrt, 3 divisions, 3 multiplications, 2 additions */
-  return dmnsn_vector_div(n, dmnsn_vector_norm(n));
+  return dmnsn_vector_mul(1.0/dmnsn_vector_norm(n), n);
 }
 
 /** Return the component-wise minimum of two vectors. */
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_min(dmnsn_vector a, dmnsn_vector b)
 {
-  return dmnsn_new_vector(
-    dmnsn_min(a.x, b.x),
-    dmnsn_min(a.y, b.y),
-    dmnsn_min(a.z, b.z)
-  );
+  dmnsn_vector v;
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    v.n[i] = dmnsn_min(a.n[i], b.n[i]);
+  }
+  return v;
 }
 
 /** Return the component-wise maximum of two vectors. */
 DMNSN_INLINE dmnsn_vector
 dmnsn_vector_max(dmnsn_vector a, dmnsn_vector b)
 {
-  return dmnsn_new_vector(
-    dmnsn_max(a.x, b.x),
-    dmnsn_max(a.y, b.y),
-    dmnsn_max(a.z, b.z)
-  );
+  dmnsn_vector v;
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    v.n[i] = dmnsn_max(a.n[i], b.n[i]);
+  }
+  return v;
 }
 
 /** Return whether a vector contains any NaN components. */
 DMNSN_INLINE bool
 dmnsn_vector_isnan(dmnsn_vector v)
 {
-  return dmnsn_isnan(v.x) || dmnsn_isnan(v.y) || dmnsn_isnan(v.z);
+  unsigned int i;
+  for (i = 0; i < 3; ++i) {
+    if (dmnsn_isnan(v.n[i])) {
+      return true;
+    }
+  }
+  return false;
 }
diff --git a/libdimension/internal.h b/libdimension/internal.h
index 3db2612..4b1910f 100644
--- a/libdimension/internal.h
+++ b/libdimension/internal.h
@@ -27,6 +27,8 @@
 #ifndef DMNSN_INTERNAL_H
 #define DMNSN_INTERNAL_H
 
+#define DMNSN_SHORT_NAMES
+
 #include "dimension/base.h"
 
 #include "internal/compiler.h"
diff --git a/libdimension/math/matrix.c b/libdimension/math/matrix.c
index 25590d8..f0050aa 100644
--- a/libdimension/math/matrix.c
+++ b/libdimension/math/matrix.c
@@ -27,53 +27,22 @@
 #include "dimension/math.h"
 #include <math.h>
 
-// Identity matrix
-dmnsn_matrix
-dmnsn_identity_matrix(void)
-{
-  return dmnsn_new_matrix(1.0, 0.0, 0.0, 0.0,
-                          0.0, 1.0, 0.0, 0.0,
-                          0.0, 0.0, 1.0, 0.0);
-}
-
-// Scaling matrix
-dmnsn_matrix
-dmnsn_scale_matrix(dmnsn_vector s)
-{
-  return dmnsn_new_matrix(s.x, 0.0, 0.0, 0.0,
-                          0.0, s.y, 0.0, 0.0,
-                          0.0, 0.0, s.z, 0.0);
-}
-
-// Translation matrix
-dmnsn_matrix
-dmnsn_translation_matrix(dmnsn_vector d)
-{
-  return dmnsn_new_matrix(1.0, 0.0, 0.0, d.x,
-                          0.0, 1.0, 0.0, d.y,
-                          0.0, 0.0, 1.0, d.z);
-}
-
 // Left-handed rotation matrix; theta/|theta| = axis, |theta| = angle
 dmnsn_matrix
 dmnsn_rotation_matrix(dmnsn_vector theta)
 {
-  // Two trig calls, 25 multiplications, 13 additions
-
   double angle = dmnsn_vector_norm(theta);
   if (fabs(angle) < dmnsn_epsilon) {
     return dmnsn_identity_matrix();
   }
-  dmnsn_vector axis = dmnsn_vector_div(theta, angle);
+  dmnsn_vector axis = dmnsn_vector_mul(1.0/angle, theta);
 
   // Shorthand to make dmnsn_new_matrix() call legible
-
   double s = sin(angle);
   double t = 1.0 - cos(angle);
-
-  double x = axis.x;
-  double y = axis.y;
-  double z = axis.z;
+  double x = axis.X;
+  double y = axis.Y;
+  double z = axis.Z;
 
   return dmnsn_new_matrix(
     1.0 + t*(x*x - 1.0), -z*s + t*x*y,        y*s + t*x*z,         0.0,
@@ -87,7 +56,7 @@ static double
 dmnsn_axis_angle(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis)
 {
   from = dmnsn_vector_sub(from, dmnsn_vector_proj(from, axis));
-  to   = dmnsn_vector_sub(to,   dmnsn_vector_proj(to, axis));
+  to = dmnsn_vector_sub(to, dmnsn_vector_proj(to, axis));
 
   double fromnorm = dmnsn_vector_norm(from);
   double tonorm   = dmnsn_vector_norm(to);
@@ -95,8 +64,8 @@ dmnsn_axis_angle(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis)
     return 0.0;
   }
 
-  from = dmnsn_vector_div(from, fromnorm);
-  to   = dmnsn_vector_div(to,   tonorm);
+  from = dmnsn_vector_mul(1.0/fromnorm, from);
+  to = dmnsn_vector_mul(1.0/tonorm, to);
 
   double angle = acos(dmnsn_vector_dot(from, to));
 
@@ -109,8 +78,7 @@ dmnsn_axis_angle(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis)
 
 // Alignment matrix
 dmnsn_matrix
-dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to,
-                       dmnsn_vector axis1, dmnsn_vector axis2)
+dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis1, dmnsn_vector axis2)
 {
   double theta1 = dmnsn_axis_angle(from, to, axis1);
   dmnsn_matrix align1 = dmnsn_rotation_matrix(dmnsn_vector_mul(theta1, axis1));
@@ -141,13 +109,13 @@ static dmnsn_matrix2 dmnsn_matrix2_sub(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs);
 static dmnsn_matrix2 dmnsn_matrix2_mul(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs);
 
 /// Invert a matrix with the slower cofactor algorithm, if partitioning failed.
-static dmnsn_matrix dmnsn_matrix_inverse_generic(dmnsn_matrix A);
+static dmnsn_matrix dmnsn_matrix_inverse_generic(const dmnsn_matrix M);
 /// Get the [\p row, \p col] cofactor of A.
-static double dmnsn_matrix_cofactor(dmnsn_matrix A, size_t row, size_t col);
+static double dmnsn_matrix_cofactor(dmnsn_matrix M, size_t row, size_t col);
 
 // Invert a matrix, by partitioning
 dmnsn_matrix
-dmnsn_matrix_inverse(dmnsn_matrix A)
+dmnsn_matrix_inverse(dmnsn_matrix M)
 {
   // Use partitioning to invert a matrix:
   //
@@ -162,11 +130,8 @@ dmnsn_matrix_inverse(dmnsn_matrix A)
   //      RR = -SS*R*inv(P), and
   //      SS = inv(S - R*inv(P)*Q).
 
-  // The algorithm uses 2 inversions, 6 multiplications, and 2 subtractions,
-  // giving 52 multiplications, 34 additions, and 8 divisions.
-
   dmnsn_matrix2 P, Q, R, S, Pi, RPi, PiQ, RPiQ, PP, QQ, RR, SS;
-  double Pdet = A.n[0][0]*A.n[1][1] - A.n[0][1]*A.n[1][0];
+  double Pdet = M.n[0][0]*M.n[1][1] - M.n[0][1]*M.n[1][0];
 
   if (dmnsn_unlikely(fabs(Pdet) < dmnsn_epsilon)) {
     // If P is close to singular, try a more generic algorithm; this is very
@@ -175,22 +140,24 @@ dmnsn_matrix_inverse(dmnsn_matrix A)
     //   [ 1 1 1 0 ]
     //   [ 0 1 1 0 ]
     //   [ 0 0 0 1 ]
-    return dmnsn_matrix_inverse_generic(A);
+    return dmnsn_matrix_inverse_generic(M);
   }
 
+  double Pdet_inv = 1.0/Pdet;
+
   // Partition the matrix
-  P = dmnsn_new_matrix2(A.n[0][0], A.n[0][1],
-                        A.n[1][0], A.n[1][1]);
-  Q = dmnsn_new_matrix2(A.n[0][2], A.n[0][3],
-                        A.n[1][2], A.n[1][3]);
-  R = dmnsn_new_matrix2(A.n[2][0], A.n[2][1],
+  P = dmnsn_new_matrix2(M.n[0][0], M.n[0][1],
+                        M.n[1][0], M.n[1][1]);
+  Q = dmnsn_new_matrix2(M.n[0][2], M.n[0][3],
+                        M.n[1][2], M.n[1][3]);
+  R = dmnsn_new_matrix2(M.n[2][0], M.n[2][1],
                         0.0,       0.0);
-  S = dmnsn_new_matrix2(A.n[2][2], A.n[2][3],
+  S = dmnsn_new_matrix2(M.n[2][2], M.n[2][3],
                         0.0,       1.0);
 
   // Do this inversion ourselves, since we already have the determinant
-  Pi = dmnsn_new_matrix2( P.n[1][1]/Pdet, -P.n[0][1]/Pdet,
-                         -P.n[1][0]/Pdet,  P.n[0][0]/Pdet);
+  Pi = dmnsn_new_matrix2( P.n[1][1]*Pdet_inv, -P.n[0][1]*Pdet_inv,
+                         -P.n[1][0]*Pdet_inv,  P.n[0][0]*Pdet_inv);
 
   // Calculate R*inv(P), inv(P)*Q, and R*inv(P)*Q
   RPi  = dmnsn_matrix2_mul(R, Pi);
@@ -204,9 +171,11 @@ dmnsn_matrix_inverse(dmnsn_matrix A)
   PP = dmnsn_matrix2_sub(Pi, dmnsn_matrix2_mul(PiQ, RR));
 
   // Reconstruct the matrix
-  return dmnsn_new_matrix(PP.n[0][0], PP.n[0][1], QQ.n[0][0], QQ.n[0][1],
-                          PP.n[1][0], PP.n[1][1], QQ.n[1][0], QQ.n[1][1],
-                          RR.n[0][0], RR.n[0][1], SS.n[0][0], SS.n[0][1]);
+  return dmnsn_new_matrix(
+    PP.n[0][0], PP.n[0][1], QQ.n[0][0], QQ.n[0][1],
+    PP.n[1][0], PP.n[1][1], QQ.n[1][0], QQ.n[1][1],
+    RR.n[0][0], RR.n[0][1], SS.n[0][0], SS.n[0][1]
+  );
 }
 
 // For nice shorthand
@@ -222,10 +191,9 @@ dmnsn_new_matrix2(double a1, double a2, double b1, double b2)
 static dmnsn_matrix2
 dmnsn_matrix2_inverse(dmnsn_matrix2 A)
 {
-  // 4 divisions, 2 multiplications, 1 addition
-  double det = A.n[0][0]*A.n[1][1] - A.n[0][1]*A.n[1][0];
-  return dmnsn_new_matrix2( A.n[1][1]/det, -A.n[0][1]/det,
-                           -A.n[1][0]/det,  A.n[0][0]/det);
+  double inv_det = 1.0/(A.n[0][0]*A.n[1][1] - A.n[0][1]*A.n[1][0]);
+  return dmnsn_new_matrix2( A.n[1][1]*inv_det, -A.n[0][1]*inv_det,
+                           -A.n[1][0]*inv_det,  A.n[0][0]*inv_det);
 }
 
 // Also basically a shorthand
@@ -240,7 +208,6 @@ dmnsn_matrix2_negate(dmnsn_matrix2 A)
 static dmnsn_matrix2
 dmnsn_matrix2_sub(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs)
 {
-  // 4 additions
   return dmnsn_new_matrix2(
     lhs.n[0][0] - rhs.n[0][0], lhs.n[0][1] - rhs.n[0][1],
     lhs.n[1][0] - rhs.n[1][0], lhs.n[1][1] - rhs.n[1][1]
@@ -251,7 +218,6 @@ dmnsn_matrix2_sub(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs)
 static dmnsn_matrix2
 dmnsn_matrix2_mul(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs)
 {
-  // 8 multiplications, 4 additions
   return dmnsn_new_matrix2(
     lhs.n[0][0]*rhs.n[0][0] + lhs.n[0][1]*rhs.n[1][0],
       lhs.n[0][0]*rhs.n[0][1] + lhs.n[0][1]*rhs.n[1][1],
@@ -262,7 +228,7 @@ dmnsn_matrix2_mul(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs)
 
 // Invert a matrix, if partitioning failed (|P| == 0)
 static dmnsn_matrix
-dmnsn_matrix_inverse_generic(dmnsn_matrix A)
+dmnsn_matrix_inverse_generic(dmnsn_matrix M)
 {
   // For A = [ A'      b ]  A^-1 = [ A'^-1   -(A'^-1)*b ]
   //         [ 0 ... 0 1 ],        [ 0 ... 0      1     ].
@@ -274,20 +240,22 @@ dmnsn_matrix_inverse_generic(dmnsn_matrix A)
   // Perform a Laplace expansion along the first row to give us the adjugate's
   // first column and the determinant
   for (size_t j = 0; j < 3; ++j) {
-    C = dmnsn_matrix_cofactor(A, 0, j);
-    det += A.n[0][j]*C;
+    C = dmnsn_matrix_cofactor(M, 0, j);
+    det += M.n[0][j]*C;
     inv.n[j][0] = C;
   }
 
+  double inv_det = 1.0/det;
+
   // Divide the first column by the determinant
   for (size_t j = 0; j < 3; ++j) {
-    inv.n[j][0] /= det;
+    inv.n[j][0] *= inv_det;
   }
 
   // Find the rest of A'
   for (size_t j = 0; j < 3; ++j) {
     for (size_t i = 1; i < 3; ++i) {
-      inv.n[j][i] = dmnsn_matrix_cofactor(A, i, j)/det;
+      inv.n[j][i] = dmnsn_matrix_cofactor(M, i, j)*inv_det;
     }
     inv.n[j][3] = 0.0;
   }
@@ -295,7 +263,7 @@ dmnsn_matrix_inverse_generic(dmnsn_matrix A)
   // Find the translational component of the inverse
   for (size_t i = 0; i < 3; ++i) {
     for (size_t j = 0; j < 3; ++j) {
-      inv.n[i][3] -= inv.n[i][j]*A.n[j][3];
+      inv.n[i][3] -= inv.n[i][j]*M.n[j][3];
     }
   }
 
@@ -306,15 +274,14 @@ dmnsn_matrix_inverse_generic(dmnsn_matrix A)
 // upper-left 3x3 corner of A by ignoring row `row' and column `col',
 // times (-1)^(row + col)
 static double
-dmnsn_matrix_cofactor(dmnsn_matrix A, size_t row, size_t col)
+dmnsn_matrix_cofactor(dmnsn_matrix M, size_t row, size_t col)
 {
-  // 2 multiplications, 1 addition
   double n[4];
   size_t k = 0;
   for (size_t i = 0; i < 3; ++i) {
     for (size_t j = 0; j < 3; ++j) {
       if (i != row && j != col) {
-        n[k] = A.n[i][j];
+        n[k] = M.n[i][j];
         ++k;
       }
     }
@@ -328,61 +295,30 @@ dmnsn_matrix_cofactor(dmnsn_matrix A, size_t row, size_t col)
   }
 }
 
-// 4x4 matrix multiplication
-dmnsn_matrix
-dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs)
-{
-  // 36 multiplications, 27 additions
-  dmnsn_matrix r;
-
-  r.n[0][0] = lhs.n[0][0]*rhs.n[0][0] + lhs.n[0][1]*rhs.n[1][0] + lhs.n[0][2]*rhs.n[2][0];
-  r.n[0][1] = lhs.n[0][0]*rhs.n[0][1] + lhs.n[0][1]*rhs.n[1][1] + lhs.n[0][2]*rhs.n[2][1];
-  r.n[0][2] = lhs.n[0][0]*rhs.n[0][2] + lhs.n[0][1]*rhs.n[1][2] + lhs.n[0][2]*rhs.n[2][2];
-  r.n[0][3] = lhs.n[0][0]*rhs.n[0][3] + lhs.n[0][1]*rhs.n[1][3] + lhs.n[0][2]*rhs.n[2][3] + lhs.n[0][3];
-
-  r.n[1][0] = lhs.n[1][0]*rhs.n[0][0] + lhs.n[1][1]*rhs.n[1][0] + lhs.n[1][2]*rhs.n[2][0];
-  r.n[1][1] = lhs.n[1][0]*rhs.n[0][1] + lhs.n[1][1]*rhs.n[1][1] + lhs.n[1][2]*rhs.n[2][1];
-  r.n[1][2] = lhs.n[1][0]*rhs.n[0][2] + lhs.n[1][1]*rhs.n[1][2] + lhs.n[1][2]*rhs.n[2][2];
-  r.n[1][3] = lhs.n[1][0]*rhs.n[0][3] + lhs.n[1][1]*rhs.n[1][3] + lhs.n[1][2]*rhs.n[2][3] + lhs.n[1][3];
-
-  r.n[2][0] = lhs.n[2][0]*rhs.n[0][0] + lhs.n[2][1]*rhs.n[1][0] + lhs.n[2][2]*rhs.n[2][0];
-  r.n[2][1] = lhs.n[2][0]*rhs.n[0][1] + lhs.n[2][1]*rhs.n[1][1] + lhs.n[2][2]*rhs.n[2][1];
-  r.n[2][2] = lhs.n[2][0]*rhs.n[0][2] + lhs.n[2][1]*rhs.n[1][2] + lhs.n[2][2]*rhs.n[2][2];
-  r.n[2][3] = lhs.n[2][0]*rhs.n[0][3] + lhs.n[2][1]*rhs.n[1][3] + lhs.n[2][2]*rhs.n[2][3] + lhs.n[2][3];
-
-  return r;
-}
-
 // Give an axis-aligned box that contains the given box transformed by `lhs'
 dmnsn_aabb
-dmnsn_transform_aabb(dmnsn_matrix trans, dmnsn_aabb box)
+dmnsn_transform_aabb(dmnsn_matrix M, dmnsn_aabb box)
 {
   // Infinite/zero bounding box support
-  if (isinf(box.min.x)) {
+  if (isinf(box.min.X)) {
     return box;
   }
 
   // Taking the "absolute value" of the matrix saves some min/max calculations
-  for (int i = 0; i < 3; ++i) {
-    for (int j = 0; j < 3; ++j) {
-      trans.n[i][j] = fabs(trans.n[i][j]);
+  dmnsn_vector Mabs[3];
+  for (unsigned int i = 0; i < 3; ++i) {
+    for (unsigned int j = 0; j < 3; ++j) {
+      Mabs[i].n[j] = fabs(M.n[j][i]);
     }
   }
 
-  dmnsn_vector Mt = dmnsn_matrix_column(trans, 3);
+  dmnsn_vector Mt = dmnsn_matrix_column(M, 3);
   dmnsn_aabb ret = { Mt, Mt };
 
-  dmnsn_vector Mz = dmnsn_matrix_column(trans, 2);
-  ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.z, Mz));
-  ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.z, Mz));
-
-  dmnsn_vector My = dmnsn_matrix_column(trans, 1);
-  ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.y, My));
-  ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.y, My));
-
-  dmnsn_vector Mx = dmnsn_matrix_column(trans, 0);
-  ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.x, Mx));
-  ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.x, Mx));
+  for (unsigned int i = 0; i < 3; ++i) {
+    ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.n[i], Mabs[i]));
+    ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.n[i], Mabs[i]));
+  }
 
   return ret;
 }
diff --git a/libdimension/model/objects/cone.c b/libdimension/model/objects/cone.c
index 26e59ca..60317bb 100644
--- a/libdimension/model/objects/cone.c
+++ b/libdimension/model/objects/cone.c
@@ -43,11 +43,11 @@ dmnsn_cone_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
 
   // Solve (x0 + nx*t)^2 + (z0 + nz*t)^2 == (((r2 - r1)*(y0 + ny*t) + r1 + r2)/2)^2
   double poly[3], x[2];
-  poly[2] = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
-  poly[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
-            - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
-  poly[0] = l.x0.x*l.x0.x + l.x0.z*l.x0.z
-            - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4.0;
+  poly[2] = l.n.X*l.n.X + l.n.Z*l.n.Z - l.n.Y*l.n.Y*(r2 - r1)*(r2 - r1)/4.0;
+  poly[1] = 2.0*(l.n.X*l.x0.X + l.n.Z*l.x0.Z)
+            - l.n.Y*(r2 - r1)*(l.x0.Y*(r2 - r1) + r2 + r1)/2.0;
+  poly[0] = l.x0.X*l.x0.X + l.x0.Z*l.x0.Z
+            - (l.x0.Y*(r2 - r1) + r2 + r1)*(l.x0.Y*(r2 - r1) + r2 + r1)/4.0;
 
   size_t n = dmnsn_polynomial_solve(poly, 2, x);
 
@@ -58,7 +58,7 @@ dmnsn_cone_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
       t = dmnsn_min(t, x[1]);
       p = dmnsn_ray_point(l, t);
 
-      if (p.y <= -1.0 || p.y >= 1.0) {
+      if (p.Y <= -1.0 || p.Y >= 1.0) {
         t = dmnsn_max(x[0], x[1]);
         p = dmnsn_ray_point(l, t);
       }
@@ -66,9 +66,9 @@ dmnsn_cone_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
       p = dmnsn_ray_point(l, t);
     }
 
-    if (t >= 0.0 && p.y >= -1.0 && p.y <= 1.0) {
-      double r = ((r2 - r1)*p.y + r1 + r2)/2.0;
-      dmnsn_vector norm = dmnsn_new_vector(p.x, -r*(r2 - r1)/2.0, p.z);
+    if (t >= 0.0 && p.Y >= -1.0 && p.Y <= 1.0) {
+      double r = ((r2 - r1)*p.Y + r1 + r2)/2.0;
+      dmnsn_vector norm = dmnsn_new_vector(p.X, -r*(r2 - r1)/2.0, p.Z);
       intersection->t      = t;
       intersection->normal = norm;
       return true;
@@ -84,9 +84,9 @@ dmnsn_cone_inside_fn(const dmnsn_object *object, dmnsn_vector point)
 {
   const dmnsn_cone *cone = (const dmnsn_cone *)object;
   double r1 = cone->r1, r2 = cone->r2;
-  double r = (point.y*(r2 - r1) + r1 + r2)/2.0;
-  return point.x*point.x + point.z*point.z < r*r
-         && point.y > -1.0 && point.y < 1.0;
+  double r = (point.Y*(r2 - r1) + r1 + r2)/2.0;
+  return point.X*point.X + point.Z*point.Z < r*r
+         && point.Y > -1.0 && point.Y < 1.0;
 }
 
 /// Cone bounding callback.
@@ -118,12 +118,12 @@ static bool
 dmnsn_cone_cap_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
                                dmnsn_intersection *intersection)
 {
-  if (l.n.y != 0.0) {
+  if (l.n.Y != 0.0) {
     const dmnsn_cone_cap *cap = (const dmnsn_cone_cap *)object;
     double r = cap->r;
-    double t = -l.x0.y/l.n.y;
+    double t = -l.x0.Y/l.n.Y;
     dmnsn_vector p = dmnsn_ray_point(l, t);
-    if (t >= 0.0 && p.x*p.x + p.z*p.z <= r*r) {
+    if (t >= 0.0 && p.X*p.X + p.Z*p.Z <= r*r) {
       intersection->t      = t;
       intersection->normal = dmnsn_new_vector(0.0, -1.0, 0.0);
       return true;
@@ -188,12 +188,8 @@ dmnsn_new_cone(dmnsn_pool *pool, double r1, double r2, bool open)
   // Implement closed cones as a union with the caps
   dmnsn_object *cap1 = dmnsn_new_cone_cap(pool, r1);
   dmnsn_object *cap2 = dmnsn_new_cone_cap(pool, r2);
-  cap1->intrinsic_trans = dmnsn_translation_matrix(
-    dmnsn_new_vector(0.0, -1.0, 0.0)
-  );
-  cap2->intrinsic_trans = dmnsn_translation_matrix(
-    dmnsn_new_vector(0.0, +1.0, 0.0)
-  );
+  cap1->intrinsic_trans = dmnsn_translation_matrix(dmnsn_new_vector(0.0, -1.0, 0.0));
+  cap2->intrinsic_trans = dmnsn_translation_matrix(dmnsn_new_vector(0.0, +1.0, 0.0));
   // Flip the normal around for the top cap
   cap2->intrinsic_trans.n[1][1] = -1.0;
 
diff --git a/libdimension/model/objects/cube.c b/libdimension/model/objects/cube.c
index 7d6fe0f..0434d89 100644
--- a/libdimension/model/objects/cube.c
+++ b/libdimension/model/objects/cube.c
@@ -23,6 +23,7 @@
  * Cubes.
  */
 
+#include "internal.h"
 #include "dimension/model.h"
 #include <math.h>
 
@@ -36,8 +37,8 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray,
   dmnsn_vector nmin, nmax;
   double tmin, tmax;
 
-  double tx1 = (-1.0 - ray.x0.x)/ray.n.x;
-  double tx2 = (+1.0 - ray.x0.x)/ray.n.x;
+  double tx1 = (-1.0 - ray.x0.X)/ray.n.X;
+  double tx2 = (+1.0 - ray.x0.X)/ray.n.X;
 
   if (tx1 < tx2) {
     tmin = tx1;
@@ -54,8 +55,8 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray,
   if (tmin > tmax)
     return false;
 
-  double ty1 = (-1.0 - ray.x0.y)/ray.n.y;
-  double ty2 = (+1.0 - ray.x0.y)/ray.n.y;
+  double ty1 = (-1.0 - ray.x0.Y)/ray.n.Y;
+  double ty2 = (+1.0 - ray.x0.Y)/ray.n.Y;
 
   if (ty1 < ty2) {
     if (ty1 > tmin) {
@@ -80,8 +81,8 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray,
   if (tmin > tmax)
     return false;
 
-  double tz1 = (-1.0 - ray.x0.z)/ray.n.z;
-  double tz2 = (+1.0 - ray.x0.z)/ray.n.z;
+  double tz1 = (-1.0 - ray.x0.Z)/ray.n.Z;
+  double tz2 = (+1.0 - ray.x0.Z)/ray.n.Z;
 
   if (tz1 < tz2) {
     if (tz1 > tmin) {
@@ -124,9 +125,9 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray,
 static bool
 dmnsn_cube_inside_fn(const dmnsn_object *cube, dmnsn_vector point)
 {
-  return point.x > -1.0 && point.x < 1.0
-      && point.y > -1.0 && point.y < 1.0
-      && point.z > -1.0 && point.z < 1.0;
+  return point.X > -1.0 && point.X < 1.0
+      && point.Y > -1.0 && point.Y < 1.0
+      && point.Z > -1.0 && point.Z < 1.0;
 }
 
 /// Boundary callback for a cube.
diff --git a/libdimension/model/objects/sphere.c b/libdimension/model/objects/sphere.c
index e1ca784..2c8a8ef 100644
--- a/libdimension/model/objects/sphere.c
+++ b/libdimension/model/objects/sphere.c
@@ -29,8 +29,7 @@
 
 /// Sphere intersection callback.
 static bool
-dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_ray l,
-                             dmnsn_intersection *intersection)
+dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_ray l, dmnsn_intersection *intersection)
 {
   // Solve (x0 + nx*t)^2 + (y0 + ny*t)^2 + (z0 + nz*t)^2 == 1
   double poly[3], x[2];
@@ -58,12 +57,12 @@ dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_ray l,
 static bool
 dmnsn_sphere_inside_fn(const dmnsn_object *sphere, dmnsn_vector point)
 {
-  return point.x*point.x + point.y*point.y + point.z*point.z < 1.0;
+  return point.X*point.X + point.Y*point.Y + point.Z*point.Z < 1.0;
 }
 
 /// Helper for sphere bounding box calculation.
 static inline double
-dmnsn_implicit_dot(double row[4])
+dmnsn_implicit_dot(const double row[4])
 {
   double ret = 0.0;
   for (int i = 0; i < 3; ++i) {
@@ -84,18 +83,18 @@ dmnsn_sphere_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
 
   double cx = trans.n[0][3];
   double dx = sqrt(dmnsn_implicit_dot(trans.n[0]));
-  box.min.x = cx - dx;
-  box.max.x = cx + dx;
+  box.min.X = cx - dx;
+  box.max.X = cx + dx;
 
   double cy = trans.n[1][3];
   double dy = sqrt(dmnsn_implicit_dot(trans.n[1]));
-  box.min.y = cy - dy;
-  box.max.y = cy + dy;
+  box.min.Y = cy - dy;
+  box.max.Y = cy + dy;
 
   double cz = trans.n[2][3];
   double dz = sqrt(dmnsn_implicit_dot(trans.n[2]));
-  box.min.z = cz - dz;
-  box.max.z = cz + dz;
+  box.min.Z = cz - dz;
+  box.max.Z = cz + dz;
 
   return box;
 }
diff --git a/libdimension/model/objects/torus.c b/libdimension/model/objects/torus.c
index b4baebd..e4894b3 100644
--- a/libdimension/model/objects/torus.c
+++ b/libdimension/model/objects/torus.c
@@ -41,18 +41,18 @@ dmnsn_torus_bound_intersection(const dmnsn_torus *torus, dmnsn_ray l)
   double rmax2 = rmax*rmax, rmin2 = rmin*rmin;
 
   // Try the caps first
-  double tlower = (-r - l.x0.y)/l.n.y;
-  double tupper = (+r - l.x0.y)/l.n.y;
+  double tlower = (-r - l.x0.Y)/l.n.Y;
+  double tupper = (+r - l.x0.Y)/l.n.Y;
   dmnsn_vector lower = dmnsn_ray_point(l, tlower);
   dmnsn_vector upper = dmnsn_ray_point(l, tupper);
-  double ldist2 = lower.x*lower.x + lower.z*lower.z;
-  double udist2 = upper.x*upper.x + upper.z*upper.z;
+  double ldist2 = lower.X*lower.X + lower.Z*lower.Z;
+  double udist2 = upper.X*upper.X + upper.Z*upper.Z;
   if ((ldist2 < rmin2 || ldist2 > rmax2) && (udist2 < rmin2 || udist2 > rmax2)) {
     // No valid intersection with the caps, try the cylinder walls
-    double dist2 = l.x0.x*l.x0.x + l.x0.z*l.x0.z;
+    double dist2 = l.x0.X*l.x0.X + l.x0.Z*l.x0.Z;
     double bigcyl[3], smallcyl[3];
-    bigcyl[2]   = smallcyl[2] = l.n.x*l.n.x + l.n.z*l.n.z;
-    bigcyl[1]   = smallcyl[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z);
+    bigcyl[2]   = smallcyl[2] = l.n.X*l.n.X + l.n.Z*l.n.Z;
+    bigcyl[1]   = smallcyl[1] = 2.0*(l.n.X*l.x0.X + l.n.Z*l.x0.Z);
     bigcyl[0]   = dist2 - rmax2;
     smallcyl[0] = dist2 - rmin2;
 
@@ -63,7 +63,7 @@ dmnsn_torus_bound_intersection(const dmnsn_torus *torus, dmnsn_ray l)
     size_t i;
     for (i = 0; i < n; ++i) {
       dmnsn_vector p = dmnsn_ray_point(l, x[i]);
-      if (p.y >= -r && p.y <= r)
+      if (p.Y >= -r && p.Y <= r)
         break;
     }
 
@@ -90,8 +90,8 @@ dmnsn_torus_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
   }
 
   // This bit of algebra here is correct
-  dmnsn_vector x0mod = dmnsn_new_vector(l.x0.x, -l.x0.y, l.x0.z);
-  dmnsn_vector nmod  = dmnsn_new_vector(l.n.x,  -l.n.y,  l.n.z);
+  dmnsn_vector x0mod = dmnsn_new_vector(l.x0.X, -l.x0.Y, l.x0.Z);
+  dmnsn_vector nmod  = dmnsn_new_vector(l.n.X,  -l.n.Y,  l.n.Z);
   double nn      = dmnsn_vector_dot(l.n, l.n);
   double nx0     = dmnsn_vector_dot(l.n, l.x0);
   double x0x0    = dmnsn_vector_dot(l.x0, l.x0);
@@ -123,7 +123,7 @@ dmnsn_torus_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
   dmnsn_vector p = dmnsn_ray_point(l, t);
   dmnsn_vector center = dmnsn_vector_mul(
     R,
-    dmnsn_vector_normalized(dmnsn_new_vector(p.x, 0.0, p.z))
+    dmnsn_vector_normalized(dmnsn_new_vector(p.X, 0.0, p.Z))
   );
   dmnsn_vector normal = dmnsn_vector_sub(p, center);
 
@@ -137,8 +137,8 @@ static bool
 dmnsn_torus_inside_fn(const dmnsn_object *object, dmnsn_vector point)
 {
   const dmnsn_torus *torus = (const dmnsn_torus *)object;
-  double dmajor = torus->major - sqrt(point.x*point.x + point.z*point.z);
-  return dmajor*dmajor + point.y*point.y < torus->minor*torus->minor;
+  double dmajor = torus->major - sqrt(point.X*point.X + point.Z*point.Z);
+  return dmajor*dmajor + point.Y*point.Y < torus->minor*torus->minor;
 }
 
 /// Torus bounding callback.
diff --git a/libdimension/model/objects/triangle.c b/libdimension/model/objects/triangle.c
index 5af3301..87372d7 100644
--- a/libdimension/model/objects/triangle.c
+++ b/libdimension/model/objects/triangle.c
@@ -33,16 +33,15 @@ static inline bool
 dmnsn_ray_triangle_intersection(dmnsn_ray l, double *t, double *u, double *v)
 {
   // See the change of basis in dmnsn_triangle_basis()
-  *t = -l.x0.z/l.n.z;
-  *u = l.x0.x + (*t)*l.n.x;
-  *v = l.x0.y + (*t)*l.n.y;
+  *t = -l.x0.Z/l.n.Z;
+  *u = l.x0.X + (*t)*l.n.X;
+  *v = l.x0.Y + (*t)*l.n.Y;
   return *t >= 0.0 && *u >= 0.0 && *v >= 0.0 && *u + *v <= 1.0;
 }
 
 /// Triangle intersection callback.
 DMNSN_HOT static bool
-dmnsn_triangle_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
-                               dmnsn_intersection *intersection)
+dmnsn_triangle_intersection_fn(const dmnsn_object *object, dmnsn_ray l, dmnsn_intersection *intersection)
 {
   double t, u, v;
   if (dmnsn_ray_triangle_intersection(l, &t, &u, &v)) {
diff --git a/libdimension/model/objects/triangle_fan.c b/libdimension/model/objects/triangle_fan.c
index 93768a9..ddca581 100644
--- a/libdimension/model/objects/triangle_fan.c
+++ b/libdimension/model/objects/triangle_fan.c
@@ -40,9 +40,9 @@ static inline dmnsn_vector
 dmnsn_change_basis(const double coeffs[6], dmnsn_vector v)
 {
   return dmnsn_new_vector(
-    coeffs[0]*v.x + coeffs[1]*v.z + v.y,
-    coeffs[2]*v.x + coeffs[3]*v.z,
-    coeffs[4]*v.x + coeffs[5]*v.z
+    coeffs[0]*v.X + coeffs[1]*v.Z + v.Y,
+    coeffs[2]*v.X + coeffs[3]*v.Z,
+    coeffs[4]*v.X + coeffs[5]*v.Z
   );
 }
 
@@ -51,9 +51,9 @@ static inline dmnsn_vector
 dmnsn_change_normal_basis(const double coeffs[6], dmnsn_vector n)
 {
   return dmnsn_new_vector(
-    coeffs[0]*n.x + coeffs[2]*n.y + coeffs[4]*n.z,
-    n.x,
-    coeffs[1]*n.x + coeffs[3]*n.y + coeffs[5]*n.z
+    coeffs[0]*n.X + coeffs[2]*n.Y + coeffs[4]*n.Z,
+    n.X,
+    coeffs[1]*n.X + coeffs[3]*n.Y + coeffs[5]*n.Z
   );
 }
 
@@ -80,12 +80,11 @@ dmnsn_compress_coeffs(double coeffs[6], dmnsn_matrix incremental)
 static inline dmnsn_matrix
 dmnsn_decompress_coeffs(const double coeffs[6])
 {
-  dmnsn_matrix incremental = dmnsn_new_matrix(
+  return dmnsn_new_matrix(
     coeffs[0], 1.0, coeffs[1], 0.0,
     coeffs[2], 0.0, coeffs[3], 0.0,
     coeffs[4], 0.0, coeffs[5], 0.0
   );
-  return incremental;
 }
 
 /// Make a change-of-basis matrix for a triangle.
@@ -100,9 +99,9 @@ dmnsn_triangle_basis(dmnsn_vector a, dmnsn_vector ab, dmnsn_vector ac)
 static inline bool
 dmnsn_ray_triangle_intersection(dmnsn_ray l, double *t, double *u, double *v)
 {
-  *t = -l.x0.z/l.n.z;
-  *u = l.x0.x + (*t)*l.n.x;
-  *v = l.x0.y + (*t)*l.n.y;
+  *t = -l.x0.Z/l.n.Z;
+  *u = l.x0.X + (*t)*l.n.X;
+  *v = l.x0.Y + (*t)*l.n.Y;
   return *t >= 0.0 && *u >= 0.0 && *v >= 0.0 && *u + *v <= 1.0;
 }
 
@@ -149,7 +148,7 @@ dmnsn_triangle_fan_inside_fn(const dmnsn_object *object, dmnsn_vector point)
   return false;
 }
 
-/// Computes the bounding box for the first triangle
+/// Computes the bounding box for the first triangle.
 static inline dmnsn_aabb
 dmnsn_bound_first_triangle(dmnsn_matrix trans)
 {
@@ -336,9 +335,9 @@ dmnsn_new_smooth_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], dmnsn_v
 
     nc = dmnsn_vector_normalized(dmnsn_transform_normal(Pabc, normals[i + 3]));
     dmnsn_vector nac = dmnsn_vector_sub(nc, na);
-    coeffs[6] = nac.x;
-    coeffs[7] = nac.y;
-    coeffs[8] = nac.z;
+    coeffs[6] = nac.X;
+    coeffs[7] = nac.Y;
+    coeffs[8] = nac.Z;
 
     P = newP;
   }
diff --git a/libdimension/model/pigments/canvas_pigment.c b/libdimension/model/pigments/canvas_pigment.c
index bb83b0a..c44ffa6 100644
--- a/libdimension/model/pigments/canvas_pigment.c
+++ b/libdimension/model/pigments/canvas_pigment.c
@@ -23,6 +23,7 @@
  * Image maps.
  */
 
+#include "internal.h"
 #include "dimension/model.h"
 
 /// Canvas pigment type.
@@ -38,8 +39,8 @@ dmnsn_canvas_pigment_fn(const dmnsn_pigment *pigment, dmnsn_vector v)
   const dmnsn_canvas_pigment *canvas_pigment = (const dmnsn_canvas_pigment *)pigment;
   dmnsn_canvas *canvas = canvas_pigment->canvas;
 
-  size_t x = llround((fmod(v.x, 1.0) + 1.0)*(canvas->width  - 1));
-  size_t y = llround((fmod(v.y, 1.0) + 1.0)*(canvas->height - 1));
+  size_t x = llround((fmod(v.X, 1.0) + 1.0)*(canvas->width  - 1));
+  size_t y = llround((fmod(v.Y, 1.0) + 1.0)*(canvas->height - 1));
   return dmnsn_canvas_get_pixel(canvas, x%canvas->width, y%canvas->height);
 }
 
diff --git a/libdimension/model/texture.c b/libdimension/model/texture.c
index b7eb7ef..dc25ac9 100644
--- a/libdimension/model/texture.c
+++ b/libdimension/model/texture.c
@@ -45,8 +45,7 @@ dmnsn_texture_initialize(dmnsn_texture *texture)
   texture->trans_inv = dmnsn_matrix_inverse(texture->trans);
 
   if (!texture->pigment->initialized) {
-    texture->pigment->trans = dmnsn_matrix_mul(texture->trans,
-                                               texture->pigment->trans);
+    texture->pigment->trans = dmnsn_matrix_mul(texture->trans, texture->pigment->trans);
     dmnsn_pigment_initialize(texture->pigment);
   }
 }
diff --git a/libdimension/pattern/checker.c b/libdimension/pattern/checker.c
index cce9623..9f378ac 100644
--- a/libdimension/pattern/checker.c
+++ b/libdimension/pattern/checker.c
@@ -23,15 +23,16 @@
  * Checker pattern.
  */
 
+#include "internal.h"
 #include "dimension/pattern.h"
 
 /// Checker pattern callback.
 static double
 dmnsn_checker_pattern_fn(const dmnsn_pattern *checker, dmnsn_vector v)
 {
-  double xmod = fmod(v.x, 2.0);
-  double ymod = fmod(v.y, 2.0);
-  double zmod = fmod(v.z, 2.0);
+  double xmod = fmod(v.X, 2.0);
+  double ymod = fmod(v.Y, 2.0);
+  double zmod = fmod(v.Z, 2.0);
 
   if (xmod < -dmnsn_epsilon)
     xmod += 2.0;
diff --git a/libdimension/pattern/leopard.c b/libdimension/pattern/leopard.c
index 1a7bce0..f17bb19 100644
--- a/libdimension/pattern/leopard.c
+++ b/libdimension/pattern/leopard.c
@@ -23,6 +23,7 @@
  * Leopard pattern.
  */
 
+#include "internal.h"
 #include "dimension/pattern.h"
 #include <math.h>
 
@@ -30,7 +31,7 @@
 static double
 dmnsn_leopard_pattern_fn(const dmnsn_pattern *leopard, dmnsn_vector v)
 {
-  double val = (sin(v.x) + sin(v.y) + sin(v.z))/3.0;
+  double val = (sin(v.X) + sin(v.Y) + sin(v.Z))/3.0;
   return val*val;
 }
 
diff --git a/libdimension/render/render.c b/libdimension/render/render.c
index 842b41e..c0878a0 100644
--- a/libdimension/render/render.c
+++ b/libdimension/render/render.c
@@ -189,10 +189,7 @@ dmnsn_rtstate_initialize(dmnsn_rtstate *state,
   state->interior     = intersection->object->interior;
 
   state->r = dmnsn_ray_point(intersection->ray, intersection->t);
-  state->pigment_r = dmnsn_transform_point(
-    intersection->object->pigment_trans,
-    state->r
-  );
+  state->pigment_r = dmnsn_transform_point(intersection->object->pigment_trans, state->r);
   state->viewer = dmnsn_vector_normalized(
     dmnsn_vector_negate(intersection->ray.n)
   );
diff --git a/libdimension/tests/bvh/prtree.c b/libdimension/tests/bvh/prtree.c
index fb0cea6..31f1490 100644
--- a/libdimension/tests/bvh/prtree.c
+++ b/libdimension/tests/bvh/prtree.c
@@ -30,13 +30,13 @@
 #include <stdio.h>
 #include <stdlib.h>
 
-unsigned int calls = 0;
+static unsigned int calls = 0;
 
 static bool
 dmnsn_fake_intersection_fn(const dmnsn_object *object, dmnsn_ray ray,
                            dmnsn_intersection *intersection)
 {
-  intersection->t = (object->aabb.min.z - ray.x0.z)/ray.n.z;
+  intersection->t = (object->aabb.min.Z - ray.x0.Z)/ray.n.Z;
   intersection->normal = dmnsn_x;
   ++calls;
   return true;
@@ -47,13 +47,10 @@ dmnsn_randomize_aabb(dmnsn_object *object)
 {
   dmnsn_vector a, b;
 
-  a.x = 2.0*((double)rand())/RAND_MAX - 1.0;
-  a.y = 2.0*((double)rand())/RAND_MAX - 1.0;
-  a.z = 2.0*((double)rand())/RAND_MAX - 1.0;
-
-  b.x = 2.0*((double)rand())/RAND_MAX - 1.0;
-  b.y = 2.0*((double)rand())/RAND_MAX - 1.0;
-  b.z = 2.0*((double)rand())/RAND_MAX - 1.0;
+  for (unsigned int i = 0; i < 3; ++i) {
+    a.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0;
+    b.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0;
+  }
 
   object->aabb.min = dmnsn_vector_min(a, b);
   object->aabb.max = dmnsn_vector_max(a, b);