path: root/include/math/mat4.h

#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.
///
/// The values are given row by row:
///
/// [ m00 m01 m02 m03 ]
/// [ m10 m11 m12 m13 ]
/// [ m20 m21 m22 m23 ]
/// [ m30 m31 m32 m33 ]
static inline mat4 mat4_make(
    R m00, R m01, R m02, R m03, R m10, R m11, R m12, R m13, R m20, R m21, R m22,
    R m23, R m30, R m31, R m32, R m33) {
  // clang-format off
  // We store the matrix in columns, this is why the following looks flipped.
  mat4 m;
  m.val[0][0] = m00; m.val[0][1] = m10; m.val[0][2] = m20; m.val[0][3] = m30;
  m.val[1][0] = m01; m.val[1][1] = m11; m.val[1][2] = m21; m.val[1][3] = m31;
  m.val[2][0] = m02; m.val[2][1] = m12; m.val[2][2] = m22; m.val[2][3] = m32;
  m.val[3][0] = m03; m.val[3][1] = m13; m.val[3][2] = m23; m.val[3][3] = m33;
  return m;
  // clang-format on
}

/// Construct a matrix from a column-major matrix array.
static inline mat4 mat4_from_array(const R M[16]) {
  // clang-format off
  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;
  // clang-format on
}

/// Construct the identity matrix.
static inline mat4 mat4_id() {
  // clang-format off
  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);
  // clang-format on
}

/// Construct a matrix from 4 column vectors.
static inline mat4 mat4_from_vec4(vec4 v0, vec4 v1, vec4 v2, vec4 v3) {
  // clang-format off
  return mat4_make(
      v0.x, v1.x, v2.x, v3.x,
      v0.y, v1.y, v2.y, v3.y,
      v0.z, v1.z, v2.z, v3.z,
      0,0,0,1);
  // clang-format on
}

/// Construct a transformation matrix from a position and vectors forming a
/// coordinate system.
static inline mat4 mat4_from_vec3(
    vec3 right, vec3 up, vec3 forward, vec3 position) {
  // clang-format off
  return mat4_make(
      right.x, up.x, forward.x, position.x,
      right.y, up.y, forward.y, position.y,
      right.z, up.z, forward.z, position.z,
      0, 0, 0, 1);
  // clant-format on
}

/// 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, vec3_neg(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));
}