cl-gmath

• Results without allocation
• Conventions
• What there is
• 2D transforms
• Doubles
• Batches
• SIMD
• Threads
• Status
• Tests

Graphics math for Common Lisp: vectors, 3x3 and 4x4 matrices, 2D affine transforms and quaternions of single floats, and of doubles, with the transforms, cameras and projections a renderer needs.

(asdf:load-system :cl-gmath)

(defpackage #:example
  (:use #:cl)
  (:local-nicknames (#:gm #:cl-gmath)))
(in-package #:example)

;; A model turned a quarter turn about Y and moved back, seen by a camera.
(defvar *model* (gm:mat4-compose (gm:make-vec3 :z -2)
                                 (gm:quat-from-axis-angle (gm:make-vec3 :y 1) (gm:radians 90))
                                 (gm:make-vec3 :x 1 :y 1 :z 1)))
(defvar *view* (gm:mat4-look-at (gm:make-vec3 :y 1 :z 3)      ; the eye
                                (gm:make-vec3 :z -2)          ; looks at the model
                                (gm:make-vec3 :y 1)))         ; with +Y up
(defvar *projection* (gm:mat4-perspective (gm:radians 60) 16/9 0.1 100))

(defvar *model-view-projection* (gm:mat4* *projection* (gm:mat4* *view* *model*)))

;; Where the model's front, (0, 0, 1), lands: turned to the right of the
;; centre of the screen, level with it, and far away, at a depth near 0.
(gm:mat4-project-point *model-view-projection* (gm:make-vec3 :z 1))
;; => #(0.19107175 1.4094049e-8 0.018630244)

Values are plain specialised arrays: a vec3 is a (simple-array single-float (3)), which can go to the GPU as it is. Functions are named for the type they take: vec3+, mat4*, quat-slerp. Scalar arguments may be any real. Each type has a counterpart of doubles, such as dvec3, with the same functions (see Doubles).

A vector or quaternion is also a typed structure over its array. Its constructor takes the components by keyword, any reals, those not given being 0, and it has an accessor for each component, which setf writes with a float of the array's precision:

(let ((velocity (gm:make-vec3 :x 3 :y -2)))
  (setf (gm:vec3-y velocity) 0.0)           ; landed: a single float, not 0
  (gm:vec3-length velocity))
;; => 3.0

An accessor checks for its own type, so vec3-x refuses a vec4.

Results without allocation

A function that returns a vector, matrix or quaternion takes an optional last argument, result. Without it, you get a new value. With it, the function writes into result, which may be one of its arguments, and allocates nothing. Code that runs every frame passes it:

(let ((world (gm:mat4-identity))
      (position (gm:make-vec3)))
  (gm:mat4-compose (gm:make-vec3 :x 1 :y 2 :z 3) (gm:make-quat) (gm:make-vec3 :x 1 :y 1 :z 1)
                   world)                                     ; into WORLD
  (gm:mat4-transform-point world (gm:make-vec3) position)     ; into POSITION
  (gm:vec3-normalize position position))                      ; in place
;; => #(0.26726124 0.5345225 0.8017837)

The tests check that each such call allocates nothing.

Conventions

• Storage: matrices are stored column by column, as GLSL stores them, and multiply column vectors: (gm:mat4*vec4 m v) is M v, and (gm:mat4* a b) applies B first. (gm:mref m row column) reads an element. gm:make-mat4 and gm:make-mat3 take their elements column by column, as GLSL's constructors do.
• Space: right-handed, +Y up, and a camera looks down −Z, as in Godot, Bevy, glTF and Blender. A positive angle turns counterclockwise when its axis points at you.
• Clip space: Vulkan's, which cl-gfx uses on every backend. X and Y run from −1 to 1 with Y down, so the top of the view is at Y −1, and depth runs from 0 to 1.
• Reversed depth: mat4-perspective, mat4-perspective-infinite and mat4-orthographic put the near plane at depth 1 and the far plane at 0. Floating-point depth then keeps its precision where perspective needs it, far away. Draw with the depth test :greater and clear depth to 0.
• Angles are in radians. radians and degrees convert.

What there is

Values	Functions
Making values	make-vec2 make-vec3 make-vec4 make-quat, by keyword: (make-vec3 :x 1 :z -1); make-mat3 make-mat3x2 make-mat4, by position, column by column
vec2 vec3 vec4	+ - * / (per component), -scale -negate -dot -length -length-squared -distance -distance-squared -normalize -lerp -min -max -copy -set -approx=, vec3-cross; -angle -reflect -project -project-on-plane -clamp-length -move-toward; components vec2-x vec2-y vec3-x to vec4-w
mat3 mat4	* *vec3/*vec4 -identity -transpose -determinant -inverse -copy -approx=; mref; mat3-from-mat4 mat3-normal-matrix
Points	mat4-transform-point mat4-transform-direction mat4-project-point; in 2D, mat3x2-transform-point mat3x2-transform-direction, and mat3-transform-point-2d mat3-transform-direction-2d of homogeneous points
mat3x2, 2D transforms	mat3x2* -identity -compose -translation -rotation -scaling -inverse -determinant -copy -approx=; mat3-from-mat3x2 mat3x2-from-mat3
Single numbers	clamp lerp inverse-lerp remap smoothstep move-toward, of any reals
quat	quat-identity quat-from-axis-angle quat* quat-rotate quat-slerp quat-conjugate quat-inverse quat-normalize quat-length quat-dot quat-copy quat-approx=; components quat-x to quat-w
Transforms	mat4-translation mat4-scaling mat4-rotation mat4-rotation-x/-y/-z mat4-from-quat mat4-compose; in 2D, mat3-translation-2d mat3-rotation-2d mat3-scaling-2d
Cameras	mat4-look-at mat4-perspective mat4-perspective-infinite mat4-orthographic
Doubles	dvec2 dvec3 dvec4 dmat3 dmat3x2 dmat4 dquat, with every function above for them: dvec3+ dmat4* dmat4*dvec4 dmat3x2-compose dquat-slerp and so on
Between precisions	vec3-from-dvec3 dvec3-from-vec3, and the same for vec2 vec4 mat3 mat3x2 mat4 quat

vec3-normalize leaves a zero vector at zero, as Godot does. Operations with no answer signal gmath-error with what to change: inverting a singular matrix, a rotation about a zero axis, or a camera whose up is the direction it looks in. Arguments of the wrong type, such as a vec4 given to vec3+, signal type-error. SBCL usually reports them at compile time.

2D transforms

A mat3x2 is a 2D affine transform: a matrix of 3 columns and 2 rows, as GLSL names it, holding the X axis, the Y axis and the origin. Its six elements, in order, are the a b c d e f of a canvas transform: a point x, y goes to a x + c y + e, b x + d y + f. mat3x2-compose makes a 2D node's transform, scaled, then turned, then moved, and mat3x2* puts a parent's before it:

(let* ((parent (gm:mat3x2-compose (gm:make-vec2 :x 100 :y 50) 0 (gm:make-vec2 :x 2 :y 2)))
       (child (gm:mat3x2-compose (gm:make-vec2 :x 10) (/ pi 2) (gm:make-vec2 :x 1 :y 1)))
       (global (gm:mat3x2* parent child)))
  (values (gm:mat3x2-transform-point global (gm:make-vec2 :x 1))
          (gm:mat3x2-transform-point (gm:mat3x2-inverse global)
                                     (gm:make-vec2 :x 120 :y 52))))
;; => #(120.0 52.0), #(1.0 0.0)

A positive angle turns +X toward +Y: counterclockwise where Y points up, and clockwise on a screen, where it points down. A mat3x2 computes exactly what the mat3 of homogeneous points with the bottom row 0 0 1 does, which mat3-from-mat3x2 makes of it.

Doubles

Every type has a counterpart of doubles, dvec2 to dvec4, dmat3, dmat3x2, dmat4 and dquat, specialised arrays of double floats, and every function one for them, named with each type's d: dvec3+, dmat4*, dmat4*dvec4, dmat4-from-dquat, dquat-slerp, and the batches dmat4*-batch and dmat4-transform-points over a double-array. Both are defined from one source, so they cannot drift apart.

dvec3-from-vec3 and its kin, one for each type, convert exactly, and vec3-from-dvec3 and its kin round each component to the nearest single float. mref takes either precision, and radians and degrees give a double for a double. A dvec3's components, dvec3-x and so on, are doubles, so setf writes them with doubles: (setf (gm:dvec3-x position) 1d0), not 1.0.

Doubles are for positions far from the origin, where single floats are a metre apart 10,000 km out. A renderer keeps positions in doubles and gives the GPU floats relative to the camera, subtracted in doubles first, as Unreal's translated world space does:

(let ((position (gm:make-dvec3 :x 1d7))               ; 10,000 km out
      (camera (gm:make-dvec3 :x (- 1d7 0.25d0))))     ; a quarter of a metre away
  (gm:vec3-from-dvec3 (gm:dvec3- position camera)))   ; what the GPU gets
;; => #(0.25 0.0 0.0)

Done in single floats, the same subtraction gives #(0.0 0.0 0.0).

The functions of doubles allocate nothing, as those of single floats do.

Batches

Batch routines work on many points or matrices at once, stored in a flat float-array: points as x, y and z one after another, matrices as 16 floats each.

• mat4-transform-points transforms the points of one array into another, or into the same one.
• mat4-transform-directions does the same for directions, without the translation.
• mat4*-batch multiplies one matrix by each matrix of an array, such as a parent's transform applied to its children's.

Each takes :start and :end, a range of items, which is exactly what a chunk of cl-threads' parallel-for gives it:

(th:parallel-for pool 0 point-count
                 (lambda (start end)
                   (gm:mat4-transform-points model points world-points
                                             :start start :end end)))

They allocate nothing.

SIMD

Where SBCL has sb-simd, on x86-64 and on ARM64 from SBCL 2.6.7, the 4x4 operations and the batches use SIMD instructions. Each has a kernel for each instruction set that measured faster than the one below it, and runs the best one the processor has. sb-simd's instruction-set-case makes the choice when cl-gmath is loaded and again whenever an image starts, so an image saved on a machine with AVX still runs on one without it.

Every kernel computes exactly what the scalar code computes: the same operations in the same order, sums added in pairs, and no fused multiply-add, which rounds once where a multiplication and an addition round twice. So every processor gets the same bits, and the tests compare each kernel with the scalar code exactly. That is what a game needs to compute the same on every player's machine, as lockstep multiplayer and replays do.

The kernels, and the time each takes on a Ryzen 9 9950X3D with SBCL 2.6.9, per call or per item of a batch:

Operation	Scalar	SSE, SSE2	AVX
mat4*	24 ns	3.7 ns	1.5 ns
mat4*vec4	3.6 ns	1.3 ns
mat4-transpose	2.3 ns	1.3 ns
mat4*-batch, a matrix	6.8 ns	4.6 ns	2.4 ns
mat4-transform-points, a point	1.1 ns	0.93 ns	0.44 ns
mat4-transform-directions, a direction	1.0 ns	0.88 ns	0.41 ns
dmat4*	24 ns	7.3 ns	2.7 ns
dmat4*dvec4	3.6 ns	2.2 ns	1.1 ns
dmat4-transpose	2.3 ns	1.2 ns	0.92 ns
dmat4*-batch, a matrix	6.9 ns		4.7 ns

A call through a function, such as mat4*, adds about 2 ns of checking its arguments and choosing its kernel; a batch adds it once. Batches of points work on four points at a time with SSE and eight with AVX, their x, y and z gathered into registers of their own. The rest is scalar.

Some instruction sets are left out, because they would not help:

• SSE4.1 adds blends, rounding, and inserting and extracting lanes, of which only the extraction could serve, for transforming points one at a time. That took 1.78 ns a point, slower than SSE's four at a time, and than the scalar code.
• AVX2 adds, of what cl-gmath could use, a permute of doubles across a whole register, which made dmat4* 7% faster and dmat4*dvec4 a third slower. Its broadcast from a register is AVX2's too, which is why the AVX kernels spread values with AVX's own permutes.
• AVX-512, which sb-simd has from SBCL 2.6.9, has no permutes there yet, which a 4x4 product needs to spread over 16 lanes, and AVX already holds a whole column of either precision.
• FMA would change the results, and in the AVX kernels made dmat4* and mat4*-batch about 9% faster, dmat4*dvec4 and dmat4*-batch slower, and the rest no faster.

On ARM64, NEON kernels do mat4*, mat4*vec4, mat4-transpose, mat4*-batch, dmat4*, dmat4*dvec4 and dmat4-transpose, the counterparts of the x86-64 kernels that were faster. They are tested under emulation, which computes exactly but cannot be timed (see Tests).

To build without SIMD, push :cl-gmath-scalar onto *features* before cl-gmath is compiled, then rebuild it with (asdf:load-system :cl-gmath :force '("cl-gmath")). Compiled files built the other way are not rebuilt on their own.

Threads

Functions keep no state. A value may be read from several threads at once. A value being written, as a result or with setf, belongs to one thread.

Status

Built in the engine's step 6, the scalar code, and step 12, its first SSE kernels. Doubles, the choice of kernels when an image starts, the AVX, SSE2 and NEON kernels, 2D affine transforms, the functions of single numbers, and the structures' constructors and accessors followed on 2026-10-14.

Tests

stift test

The tests run every kernel the processor has, so a machine with AVX tests the SSE and AVX kernels alike. The NEON kernels need ARM64. Without an ARM64 machine, QEMU's user-mode emulator runs the suite on x86-64 with an ARM64 SBCL of 2.6.7 or later, such as Alpine's, unpacked with Alpine's minimal root file system into a directory, here sysroot, as CI does on every push:

qemu-aarch64 -L sysroot sysroot/usr/bin/sbcl --non-interactive \
  --eval "(load \"$(stift setup)\")" --eval '(asdf:test-system :cl-gmath)'