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 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.002507683297243826
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, moonnum hands it a fixed-point real at whatever width you ask for, and neither can drift from the other because there is only one implementation to drift.

let bits = @moonnum.Fixed::bits_for_digits(50)
let ops = @generic.fixed_ops(bits)
let half = @moonnum.Fixed::of_decimal("0.5", bits)
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 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 twenty-six more across the whole range; 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