moonmath — the arithmetic of space

Every frame a renderer draws is an argument in geometry, carried out several million times a second. This library is that argument written down once: the vectors a vertex shader already has, the matrix a camera is, the rotations that interpolate better than a matrix does, and the predicates that answer whether two things touch. It is not a renderer. It talks to no GPU and loads no mesh; it produces the sixteen numbers you hand to one.

§ I

The page and the buffer

Fig. 1. The unit cube, seen through look_at and perspective at 30° and drawn with the numbers they returned, in a frame of the same proportions. Perspective shortens the far face; the three edges that meet behind it are dashed, as a draughtsman would leave them.

A matrix is four columns. That is not a storage detail to be optimised away later; it is the shape of the thing, and it is the shape a mat4x4<f32> has in a WGSL buffer. So Mat4 holds c0, c1, c2, c3, and to_array writes them in that order, sixteen floats that a single copy puts where the shader expects them.

$$ M = \begin{bmatrix} c_{0x} & c_{1x} & c_{2x} & c_{3x} \\ c_{0y} & c_{1y} & c_{2y} & c_{3y} \\ c_{0z} & c_{1z} & c_{2z} & c_{3z} \\ c_{0w} & c_{1w} & c_{2w} & c_{3w} \end{bmatrix} $$

Fig. 2. Sixteen floats, in the order they sit in the buffer. The fourth group is the translation, which is why it lands at the end of the array.

A vector multiplies against those columns, not against rows: the result is the columns, weighted by the components of the vector and added.

$$ Mv \;=\; \sum_{k=0}^{3} c_k\, v_k $$

Everything else follows from that one line. Two matrices multiply by applying one to each column of the other, so the product applies the right-hand matrix first — $(AB)v = A(Bv)$ — which is the order a scene graph composes in: a model, then a view, then a projection.

listing 1 · the sixteen floats

let moved = @moonmath.Mat4::from_translation(
  @moonmath.Vec3::new(3.0, 0.0, 0.0),
)
let buffer = moved.to_array()

printed

[1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 3, 0, 0, 1]

The translation sits in the twelfth, thirteenth and fourteenth places because it sits in the fourth column. That is the whole of what “column-major” means, and it is why a matrix can be copied into a uniform buffer without being transposed on the way.

§ II

Where the depth goes

WebGPU clips depth to the range zero to one. A projection matrix is the thing that agrees to it, and it is short enough to derive on one page.

Work in view space: the camera at the origin, looking down $-z$. A point in front of it has $z = -d$ with $d > 0$, and $d$ is the distance the depth will be measured against. The vertical field of view is $\theta$, the aspect ratio $a$, the near plane $n$ and the far plane $f$.

$$ x_c = \frac{x}{a\tan(\theta/2)}, \qquad y_c = \frac{y}{\tan(\theta/2)}, \qquad w_c = -z $$
Scale $x$ and $y$ so that the edges of the frustum land on ±1 after the divide, and keep the distance in $w$.
$$ \frac{x_c}{w_c} = \frac{x}{-z\,a\tan(\theta/2)} $$
At the edge of the frustum $x = -z\,a\tan(\theta/2)$, so the quotient is one. The divide is what makes a perspective perspective; nothing else here is nonlinear.
$$ z_c = Az + B, \qquad \frac{z_c}{-z}\Big|_{z=-n} = 0, \qquad \frac{z_c}{-z}\Big|_{z=-f} = 1 $$
Depth is the one row that gets to be affine in $z$: two unknowns, two conditions, and the planes are pinned.
$$ \begin{aligned} \frac{-An + B}{n} = 0 \;&\Rightarrow\; B = An, \\[2pt] \frac{-Af + An}{f} = 1 \;&\Rightarrow\; A = -\frac{f}{f-n} \end{aligned} $$
Substituting back gives $B = -\dfrac{fn}{f-n}$, both negative once $f > n$.
$$ z_{\text{ndc}} = \frac{fz + fn}{(f-n)\,z} $$
At $z=-n$ the numerator is $-fn + fn = 0$; at $z=-f$ it is $-f^2 + fn = -f(f-n)$, over $(f-n)(-f)$, which is one. The library checks both, and prints $0$ and $1$.
$$ \frac{dz_{\text{ndc}}}{dz} = -\frac{fn}{(f-n)\,z^2} $$
The step grows with the inverse square of the distance, so most of the range is spent within a few metres of the eye. That is not an accident of the formula; it is the formula.

Assembled, the four rows become the matrix that Mat4::perspective returns:

$$ P = \begin{bmatrix} \dfrac{1}{a\tan(\theta/2)} & 0 & 0 & 0 \\[4pt] 0 & \dfrac{1}{\tan(\theta/2)} & 0 & 0 \\[4pt] 0 & 0 & \dfrac{f}{n-f} & \dfrac{fn}{n-f} \\[4pt] 0 & 0 & -1 & 0 \end{bmatrix} $$

Fig. 3. The frustum from the side. The projection sends the plane at $n$ to depth 0 and the plane at $f$ to 1; everything between arrives in between, crowded towards the near end.

printed · the two ends

near depth = 0
far depth = 1

the other convention

OpenGL clips to minus one to one. The same derivation with those two conditions gives $A = \dfrac{f+n}{n-f}$ and $B = \dfrac{2fn}{n-f}$, and the _gl variants are exactly that — perspective_gl, orthographic_gl, frustum_gl. Nothing else changes.

§ III

Twelve orders of turning

Three turns about coordinate axes, taken in some order, describe any rotation. There are twelve such orders, and the library names all twelve because a file written by another tool will use one of them.

Six of the orders turn about three different axes — the Tait–Bryan family, which is what a yaw, a pitch and a roll are. The other six turn about the first axis twice, and are called proper Euler angles, though they are no more proper than the others.

Tbl. 1 · the twelve orders, and where each one loses a degree of freedom
order	turns about	locks when
Xyz	x, then y, then z	y = ±90°
Xzy	x, then z, then y	z = ±90°
Yxz	y, then x, then z	x = ±90°
Yzx	y, then z, then x	z = ±90°
Zxy	z, then x, then y	x = ±90°
Zyx	z, then y, then x	y = ±90°
Xyx	x, then y, then x	y = 0° or 180°
Xzx	x, then z, then x	z = 0° or 180°
Yxy	y, then x, then y	x = 0° or 180°
Yzy	y, then z, then y	z = 0° or 180°
Zxz	z, then x, then z	x = 0° or 180°
Zyz	z, then y, then z	y = 0° or 180°

from_euler reads a Vec3 by axis: the component named x is the turn about $x$, wherever in the order that turn comes. For the six Tait–Bryan orders that is exact and complete. For the other six it cannot be, because two turns about $x$ do not fit in the one component named after it — so those six take a Turns record, whose fields are named by their place in the order rather than by an axis.

listing 2 · two readings of the same turn

// By axis: the third component turns about z, which
// Zyx turns first.
let yaw = @moonmath.Quat::from_euler(
  @moonmath.EulerOrder::Zyx,
  @moonmath.Vec3::new(0.0, 0.0, @moonmath.to_radians(90.0)),
)
println(yaw)
println(yaw.to_euler(@moonmath.EulerOrder::Zyx))

// In order: Xyx turns about x twice, so the turns are
// named by place.
let precessing = @moonmath.Quat::from_euler_turns(
  @moonmath.EulerOrder::Xyx,
  @moonmath.Turns::new(0.3, 1.2, -0.4),
)
println(precessing)

printed

Quat(0, 0, 0.7071068, 0.7071067)
Vec3(0, 0, 1.5707965)
Quat(-0.041249573, 0.53040975, 0.19361468, 0.8243042)

§ IV

Two hands, one template

A GPU buffer holds single precision, and a simulation that runs for hours or spans a solar system does not. Past a few kilometres, a Float can no longer tell one metre from the next: at a billion, the gap between two representable numbers is sixty-four.

listing 3 · the same sum in both hands

let far = @f64.Vec3::new(1.0e9, 0.0, 0.0)
let step = @f64.Vec3::new(1.0e-3, 0.0, 0.0)
println("\{far + step}")

let near = @moonmath.Vec3::new(1.0e9, 0.0, 0.0)
let small = @moonmath.Vec3::new(1.0e-3, 0.0, 0.0)
println("\{near + small}")

printed

Vec3(1000000000.001, 0, 0)
Vec3(1000000000, 0, 0)

The millimetre survives in the first line and is gone in the second. The double-precision package is the same library: the same types, the same names, the same operators, generated from one template into two packages, so a program written against one reads the same against the other with a different import.

listing 4 · importing the second hand

import {
  "moonbitstack/moonmath",
  "moonbitstack/moonmath/f64",
}

let single = @moonmath.Mat4::identity()
let double = @f64.Mat4::identity()

The integer vectors — IVec2, UVec3 and their siblings — exist only in single precision, because they count grid cells and pixels, and a pixel has no fractional part.

§ V

What is in the box

Twenty-seven types, grouped by what they are for. Each links to its entry in the lexicon, where every field, method and value has a sentence.

family	types
Vectors	Vec2, Vec3, Vec4 — one, two, three and four components, with dot, cross, normalize, lerp, reflect, refract and the swizzles
Integer vectors	IVec2, IVec3, IVec4, UVec2, UVec3, UVec4 — signed and unsigned, for grid cells, pixels and indices
Matrices	Mat2, Mat3, Mat4 — products, transposes, determinants, inverses, and the cameras: `look_at`, `perspective`, `orthographic`, `frustum`
Rotation	Quat, EulerOrder, Turns — axis and angle, two vectors, slerp and nlerp, twelve Euler orders, and the conversion to and from both matrix sizes
Placement	Trs, Affine2, Affine3 — a rotation, a scale and a translation kept apart, and the two transforms that leave the last row of a matrix alone
Geometry	Rect, Aabb2, Aabb3, Ray2, Ray3, Plane, Sphere, Triangle, Frustum — shape, and whether two shapes meet
Values	`approx`, `lerp`, `smoothstep`, `remap`, `wrap`, `to_radians` and the constants beside them

turn the page From a first vector to a culled scene

The Arithmeticof Space

The page and the buffer

Where the depth goes

Twelve orders of turning

Two hands, one template

What is in the box

The Arithmetic
of Space