Eight steps, from a vector to a culled scene

The package is one module with two packages inside: the default single precision one, and moonmath/f64. Adding the module is enough for both.

moon add moonbitstack/moonmath

A vector is a value: immutable, small, and returned by everything that takes one. The first useful thing to know is the handedness, because it decides a sign you will meet on the first day.

let forward = @moonmath.Vec3::new(0.0, 0.0, -1.0)  // the way the eye looks
let up = @moonmath.Vec3::unit_y()
let right = forward.cross(up)                       // +x, not -x

println(right)
println(forward.dot(up))

Vec3(1, 0, 0)
0

x runs east, y north and z towards the viewer, so cross(x, y) is z and a camera looking down −z has its right hand along +x. A left-handed system would give the opposite sign here, and a mirrored world at the end of it.

The operators work component by component, between two values of the same type. a * b is not a dot product and not a matrix product: it is the product of the components, which is what a shader means by a * b too. The dot product has a name because it is a different thing.

let a = @moonmath.Vec3::new(1.0, 2.0, 3.0)
let b = @moonmath.Vec3::new(4.0, -5.0, 6.0)

println(a + b)          // component by component
println(a * b)          // the same, multiplied
println(a.dot(b))       // one number
println(a.length())
println(a.normalize())
println(a.lerp(b, 0.25))  // a quarter of the way to b
println(a.min(b))        // component by component, the smaller

Vec3(5, -3, 9)
Vec3(4, -10, 18)
12
3.7416575
Vec3(0.26726124, 0.5345225, 0.8017837)
Vec3(1.75, 0.25, 3.75)
Vec3(1, -5, 3)

The zero vector has no direction, and normalizing it would divide by zero. The library returns the zero vector rather than a NaN, and try_normalize is there for callers who need to know which case they are in.

let nothing = @moonmath.Vec3::zero()

println(nothing.normalize())       // no NaN escapes
println(nothing.try_normalize())   // the caller can tell

Vec3(0, 0, 0)
None

The scalars beside them

Angles are radians everywhere in the library, and the conversion happens at the edge, where the number came from a file or a person. The easing and wrapping helpers are plain functions, not methods, because they belong to no type in particular.

println(@moonmath.smoothstep(0.0, 1.0, 0.25))  // 3t² − 2t³, at t = 0.25
println(@moonmath.wrap(7.5, 5.0))              // into 0 … 5
println(@moonmath.to_degrees(@moonmath.half_pi))

0.15625
2.5
90

A matrix is built from its columns, because that is the order it is stored in. Three columns, three vectors, one matrix.

let m = @moonmath.Mat3::from_cols(
  @moonmath.Vec3::new(1.0, 2.0, 3.0),
  @moonmath.Vec3::new(4.0, 5.0, 6.0),
  @moonmath.Vec3::new(7.0, 8.0, 10.0),
)

println(m.determinant())
println(m * m.inverse().unwrap())
println((m * m.inverse().unwrap()).approx(@moonmath.Mat3::identity()))

-3
Mat3[Vec3(1, 0, 0), Vec3(0, 1, 0), Vec3(0, 0, 1)]
true

inverse returns an Option, because a matrix can be singular and the caller has to decide what that means. When the case cannot happen, unwrap says so at the point where the knowledge is. The / operator aborts instead, since it has nowhere to report anything.

A matrix times a vector is not *, because the two sides would have different types. It is a method, and there is one for each size:

let v3 = @moonmath.Vec3::new(1.0, 0.0, 0.0)
println(m.mul_vec3(v3))

let v4 = @moonmath.Vec4::new(1.0, 0.0, 0.0, 1.0)
let big = @moonmath.Mat4::from_mat3(m)
println(big.mul_vec4(v4))

Vec3(1, 2, 3)
Vec4(1, 2, 3, 1)

A camera is a view matrix and a projection matrix. The first moves the world so that the eye is at the origin looking down −z; the second squeezes what it sees into the unit cube the clipper expects. Multiplying them gives the matrix a vertex shader takes.

let eye = @moonmath.Vec3::new(0.0, 0.0, 6.0)
let view = @moonmath.Mat4::look_at(
  eye, @moonmath.Vec3::zero(), @moonmath.Vec3::unit_y(),
)
let projection = @moonmath.Mat4::perspective(
  @moonmath.to_radians(60.0), 1.5, 0.1, 100.0,
)
let view_projection = projection * view

println(view_projection.transform_point(@moonmath.Vec3::zero()))

Vec3(0, 0, 0.9843177)

The origin of the world is dead centre of the screen, so its x and y are zero, and its depth is 0.98 — nearly at the far plane, because the eye is six metres away and the far plane is a hundred. Point the same matrix at the two planes and the convention shows itself:

// 0.1 in front of the eye, and 100 in front of the eye.
let on_near = @moonmath.Vec3::new(0.0, 0.0, 5.9)
let on_far = @moonmath.Vec3::new(0.0, 0.0, -94.0)
println(view_projection.transform_point(on_near))
println(view_projection.transform_point(on_far))

Vec3(0, 0, 0)
Vec3(0, 0, 1)

Zero and one, exactly, as WebGPU clips. The _gl variants of the same four projections map the planes to minus one and one instead, for a pipeline that wants that.

A model matrix is a scale, a rotation and a translation, and keeping the three apart is what makes a transform editable. Trs holds them; Mat4::from_trs turns them into the matrix the shader wants; to_trs reads them back.

let placement = @moonmath.Trs::new(
  @moonmath.Vec3::new(0.0, 2.0, 0.0),
  @moonmath.Quat::from_axis_angle(
    @moonmath.Vec3::unit_y(),
    @moonmath.half_pi,
  ),
  @moonmath.Vec3::splat(0.5),
)
let model = @moonmath.Mat4::from_trs(placement)

// A local +x, halved, turned a quarter turn about y, then lifted two.
println(model.transform_point(@moonmath.Vec3::new(1.0, 0.0, 0.0)))
println(model.transform_vector(@moonmath.Vec3::new(1.0, 0.0, 0.0)))

Vec3(-5.9604645e-8, 2, -0.49999997)
Vec3(-5.9604645e-8, 0, -0.49999997)

The point is lifted by the translation and the direction is not, which is the whole difference between the fourth component being one and zero. The two tiny values are a quarter turn done in single precision: \(\cos 90°\) is not quite zero, and the matrix says so rather than rounding it away.

let back = model.to_trs().unwrap()

println(back.translation)
println(back.scale)

Vec3(0, 2, 0)
Vec3(0.49999997, 0.5, 0.49999997)

Not every matrix can be taken apart this way. A shear leans the three columns towards each other, and no rotation and scale produces that, so to_trs answers None — the same answer a singular matrix gets.

A quaternion is four numbers that interpolate, and a matrix is what a shader wants. Both directions are here, and a rotation survives the round trip.

let turn = @moonmath.Quat::from_axis_angle(
  @moonmath.Vec3::unit_y(),
  @moonmath.quarter_pi,
)
let half_way = @moonmath.Quat::identity().slerp(turn, 0.5)

println(turn.rotate(@moonmath.Vec3::unit_x()))
println(half_way.rotate(@moonmath.Vec3::unit_x()))
println(turn.to_mat3())
println(@moonmath.Quat::from_mat3(turn.to_mat3()).approx(turn))

Vec3(0.7071067, 0, -0.7071068)
Vec3(0.9238795, 0, -0.3826835)
Mat3[Vec3(0.7071067, 0, -0.7071068), Vec3(0, 1, 0), Vec3(0.7071068, 0, 0.7071067)]
true

A quarter turn about y sends +x towards −z, and half of that turn sends it almost nowhere: slerp walks the arc, so the half-way rotation is a rotation by 22.5°, not a linear blend of the four components.

let yaw = @moonmath.Quat::from_euler(
  @moonmath.EulerOrder::Zyx,
  @moonmath.Vec3::new(0.0, 0.0, @moonmath.to_radians(90.0)),
)
println(yaw.to_euler(@moonmath.EulerOrder::Zyx))

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

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

The first component of the vector is the turn about x, wherever in the order that turn comes. All twelve orders are in the lexicon, with the angle at which each one stops being able to say what it means.

A bone in a skinned mesh is a turn and a move at once, and the two have to be weighted together. DualQuat is that pair: two quaternions, one carrying the rotation and one carrying half the translation turned by it.

let bone = @moonmath.DualQuat::from_rotation_translation(
  @moonmath.Quat::from_axis_angle(@moonmath.Vec3::unit_x(), 0.4),
  @moonmath.Vec3::new(0.0, 0.0, 1.0),
)
let rest = @moonmath.DualQuat::from_rotation_translation(
  @moonmath.Quat::from_axis_angle(@moonmath.Vec3::unit_x(), 1.2),
  @moonmath.Vec3::new(0.0, 2.0, 3.0),
)

let half = bone.slerp(rest, 0.5)
println(half.rotation())
println(half.translation())
println(half.transform_point(@moonmath.Vec3::unit_y()))

let blended = (bone.scale(0.5) + rest.scale(0.5)).normalize()
println(blended.to_mat4().transform_point(@moonmath.Vec3::unit_y()))
println(blended.approx(half))

Quat(0.38941836, 0, 0, 0.92106104)
Vec3(0, 1, 2.0000002)
Vec3(0, 1.6967067, 2.7173564)
Vec3(0, 1.8994168, 2.5146463)
false

The two lines that end in a point are two different curves. A slerp takes the rotation along its arc and the translation along the straight line between the two; a weighted sum is what a skinning loop has to do with the bones it is given, and normalizing it divides both halves by the length of the summed rotation, so the translation comes out longer whenever the two rotations differ. Both are here because an animation clip wants the first and a skin wants the second.

A ray is a start and a direction. Asking where it meets a box gives a distance along it, or None.

let room = @moonmath.Aabb3::new(
  @moonmath.Vec3::new(-4.0, 0.0, -4.0),
  @moonmath.Vec3::new(4.0, 3.0, 4.0),
)
let down = @moonmath.Ray3::new(
  @moonmath.Vec3::new(0.0, 2.0, 0.0),
  @moonmath.Vec3::new(0.0, -1.0, 0.0),
)
println(down.intersect_aabb(room))   // the floor is two below
println(down.at(2.0))

let ball = @moonmath.Sphere::new(@moonmath.Vec3::zero(), 1.0)
let ahead = @moonmath.Ray3::new(
  @moonmath.Vec3::new(0.0, 0.0, 5.0),
  @moonmath.Vec3::new(0.0, 0.0, -1.0),
)
println(ahead.intersect_sphere(ball))

Some(2)
Vec3(0, 0, 0)
Some(4)

The last piece a renderer needs is the cheapest: throw away what the camera cannot see before drawing it. A frustum built from the same view-projection matrix answers for a box in one call per plane.

let frustum = @moonmath.Frustum::from_matrix(view_projection).unwrap()
println(frustum.intersects_aabb(room))       // in front of the camera

let far_room = @moonmath.Aabb3::new(
  @moonmath.Vec3::splat(50.0),
  @moonmath.Vec3::splat(60.0),
)
println(frustum.intersects_aabb(far_room))   // behind it

true
false

Five more shapes sit beside the vectors and the matrices, each answering a question the fundamentals do not. A Half is a sixteen-bit float, the storage a packed vertex uses. A Shear is the off-diagonal linear map a skew is. A Color3 carries the arithmetic of a blend, and rgb_to_hsv turns it into the hue, saturation and value an editor shows. An Interval is a one-dimensional range, and a Line3 is the infinite line a construction axis or a normal is.

println(@moonmath.Half::from_float(1.5))
println(@moonmath.Shear::new(0.3, 0.0, 0.0).to_mat3().get(0, 1))
println(@moonmath.rgb_to_hsv(@moonmath.Color3::new(1.0, 0.0, 0.0)))
println(
  @moonmath.Interval::new(0.0, 1.0).union(@moonmath.Interval::new(2.0, 3.0)),
)
let axis = @moonmath.Line3::new(@moonmath.Vec3::zero(), @moonmath.Vec3::unit_x())
println(axis.distance_to_point(@moonmath.Vec3::new(3.0, 4.0, 0.0)))
println(@moonmath.Vec3::random_on_sphere(0.5, 0.0))

Half(1.5)
0.3
Color3(0, 1, 1)
Interval(0, 3)
4
Vec3(1, 0, 0)

The sampling functions take the two or three uniform values in [0, 1) a caller's random source provides, and return the point with the right distribution. Nothing here owns a generator, so the same call with the same values always returns the same point.

moonbitstack/moonmath/special is a second package for the functions a probability, a physics or a numerical caller reaches for: the gamma and beta families with their incomplete and logarithmic forms, the error function and its inverse, the four Bessel functions, the elliptic integrals in Carlson form, the hypergeometric functions and the Riemann zeta.

println(@special.gamma(0.5))              // the square root of pi
println(@special.bessel_j(0.0, 2.4))      // a struck drum
println(@special.gammainc(2.0, 1.0))      // the incomplete gamma
println(@special.incbeta(2.0, 3.0, 0.5))  // the incomplete beta
println(@special.hyper2f1(1.0, 1.0, 2.0, 0.5))
println(@special.zeta(2.0))               // pi squared over six

1.7724538509055165
0.002507683297243856
0.2642411176571134
0.6875000000000018
1.3862943611198901
1.6449340668482266

Every body is written once. It lives in special/generic, against an Ops record that names the operations and constants a number type has to supply — addition, a square root, a logarithm, pi. The package you call is a face over that body: special hands it the double algebra, bigfloat_ops hands it a real at whatever width you ask for, and neither can drift from the other because there is only one implementation to drift.

let ops = @generic.bigfloat_ops(50)
let half = (ops.from_decimal)("0.5")
let value = @generic.gamma(half, ops)

Which method answers is decided by the argument, not by taste. The Bessel functions are a power series up to x = 1, an asymptotic expansion once x is several times the order squared, and Steed's method in between — a continued fraction for the ratio of neighbouring orders, a second one for the ratio of the two kinds, and the Wronskian to fix the scale. K takes the Thompson and Barnett fraction past x = 2 for the same reason: its series and its expansion both degrade, and no choice of crossover closes the gap between them. The incomplete gamma and beta are a series on one side of the mean and a continued fraction on the other.

The two hypergeometric functions are where double precision stops being enough on its own, and both say how far short they fall. A series multiplies its way up to terms far above the answer and adds them back down, and what survives is the digits the two ends had in common. hyper1f1 reports the largest term it formed; hyper2f1 reports the same, and past a half its transformation multiplies a gamma ratio of 10^60 by series whose own terms reach 10^55, so those scales are folded into the first term of the series they multiply rather than applied to the sum afterwards. When the report comes back short, the double face re-runs the whole thing on moonnum’s BigFloat at a width that holds the digits the cancellation took, and rounds once. A 2F1 whose c - a - b is a whole number is the one case the transformation cannot take: its two powers are linearly dependent there, and the value is read off a logarithm instead, which stays bounded however close the argument gets to one.

The numbers are checked three ways, and only one of them is this package's own arithmetic. mpmath supplies references at double, at fifty and at a hundred digits; Boost.Math's test tables supply thirty-three more, fourteen thousand rows, across the whole range — the three large-parameter tables among them, 2,402 rows with parameters past a million, replayed row for row, and those are what pin the routing down — and the two agree with each other, which is what makes either worth believing.

Import moonbitstack/moonmath/f64 and every type appears again with Double where Float was. The names, the operators and the documentation are the same; only the numbers change.

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)

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

At a billion, single precision counts in steps of sixty-four, so a millimetre is not a small distance — it is no distance at all. A simulation that spans a solar system wants the second package; what it sends to the GPU still wants the first.

turn the page The lexicon, every public name