pmndrs/math and TypeGPU: out-parameters on the CPU, value types on the GPU

mdjsonmcpcopy↗ai

2026-09-26 · 21 min · webgpu · typescript · gpu · math · geometry · 3d · performance · benchmarks · open-source · procedural-generation · explainer

On 25 September pmndrs announced math v0.1.0: "Our resident expert @isaac_mason_ hand crafted the playful web's most optimized math engine." Four bullets followed: allocation-free, monomorphic and benchmarked for speed; tiny and tree-shakable; works with WebGL, WebGPU, Wasm or your renderer; data-oriented. About five hours later Iwo Plaza, a TypeGPU developer who wrote its 0.11 and 0.12 release posts, quoted it: "Built IK and terrain with math, rendering with typegpu."

That is two libraries with opposite answers to one question: where does a vector live, and who owns it? I cloned pmndrs/math at the math@0.1.0 tag (c6713e3, 11 September) and software-mansion/TypeGPU at b819e00 (26 September) and read both, plus gl-matrix (6f96d57, 3.4.4), wgpu-matrix (3dba901, 3.4.2) and three.js's src/math (r186) at their current heads for the comparisons. The analysis runs none of their code; my own code does the measuring. One widget near the end runs math itself.

A gl-matrix API on plain arrays

The repository is pmndrs's older maath helpers, rewritten over August and renamed on 23 August (8a24199). The README says where the core comes from: the vector, quaternion and matrix code "started life as a port of mathcat, which started as a TypeScript port of glMatrix". It shows. Every function takes the output first and returns it:

// pmndrs/math, src/core/vec3.ts
export type Vec3 = [x: number, y: number, z: number];
 
export function create(): Vec3 {
    return [0, 0, 0];
}
 
export function add(out: Vec3, a: Vec3, b: Vec3): Vec3 {
    out[0] = a[0] + b[0];
    out[1] = a[1] + b[1];
    out[2] = a[2] + b[2];
    return out;
}

The difference from gl-matrix is the storage. The skill file the package ships for coding agents, skills/math/SKILL.md, puts it plainly: "Every type is a plain fixed-length tuple of numbers — no classes, no wrappers, and not a typed array". A Vec3 is a JavaScript array of three doubles; a Mat4 is sixteen, column-major, with the translation in m[12] to m[14]. The same file argues for it: "Don't assume a typed array is faster. A packed plain array is already unboxed and can still grow." Typed arrays are for interop and footprint, it says, and a Float32Array costs "a narrowing conversion on every write".

The rest of the style follows from the caller owning the data:

• Aliasing is allowed. vec3.normalize(v, v) and vec3.cross(a, a, b) work because no function writes a component of out before it has read every input it still needs. cross reads all six inputs into locals first; mat4.multiply caches all sixteen of a, then one column of b at a time.
• Scratch lives at module scope, named _owner_purpose, allocated once. The SKILL file warns that it is not reentrant, and says to pass workspace in for recursive or worker code.
• State objects keep one shape. fabrik2.createChain2() builds every field, including solveDistance: Number.POSITIVE_INFINITY and an empty bestSolution array, in the same order every time, so the hidden class never changes. That is what "monomorphic" means here: the same shapes arrive at every call site.

The whole library is eight entrypoints: math (vectors, quaternions, Euler angles, matrices, spherical and polar coordinates), math/shapes (boxes, OBBs, planes, frustums, raycasts), math/geometry (quickhull in 2-D and 3-D, polygon decomposition and triangulation), math/time (easings and springs), math/random (ISAAC, mulberry32), math/noise (Perlin, simplex, Worley, fBm, curl), math/color and math/ik (FABRIK in 2-D and 3-D).

Is it allocation-free?

The per-frame core is. I searched src/ for arrays and typed arrays created inside functions. The vector, quaternion and matrix operations write into out and create nothing, with one exception that matters because it ships: measured at the tag, vec3.rotateX, rotateY and rotateZ each allocate two arrays per call.

// pmndrs/math at math@0.1.0, src/core/vec3.ts
export function rotateX(out: Vec3, a: Vec3, b: Vec3, rad: number): Vec3 {
    const p: number[] = [];
    const r: number[] = [];
    //Translate point to the origin
    p[0] = a[0] - b[0];
    // ...

That is gl-matrix's code, which still has it. A contributor rewrote all three with locals in 55e8c25 on 14 September, three days after the release; it is on main and in the canary builds, not in 0.1.0, which is what npm i math installs today.

The allocations elsewhere are deliberate. Quickhull and polygon decomposition return variable-length results, so they allocate them. addBone allocates a bone and grows the chain's best-pose scratch by four numbers "so solve never allocates". One pattern is worth knowing about. The fractal helpers take a callback, and the documented call is fbm((f) => simplex2d.sample(gen, x * f, y * f), 5, 2, 0.5). Reasoned: that arrow captures x and y, so written per vertex it is a new closure per sample, unless V8 inlines fbm and escape-analyses the closure away.

The benchmark files show how seriously the zero is taken. A comment in benches/ik/fabrik2.bench.ts explains why the target advances by an integer step: "a module-level let holding a double is boxed into a fresh HeapNumber on every write, which would show up as ~16 bytes/iter of allocation and hide the solver's own zero".

How it differs from gl-matrix, three.js and wgpu-matrix

Measured, reading each library's source:

gl-matrix and math share a signature, so porting between them is mostly deleting new Float32Array. wgpu-matrix's optional destination is convenient, and a hidden allocation every time it is left off. three.js keeps module scratch internally too (const _vector = /*@__PURE__*/ new Vector3()); math's SKILL file says to marshal its objects through toArray(scratch), since toArray() bare allocates.

What "most optimized" rests on

The benchmarks live in benches/ and run on @pmndrs/labs 0.9. Measured: 24 files and 123 benchmarks. Of those, 97 are micro-benchmarks, one function over 10,000 inputs. Six, in five groups, are composite: a small feature built from many calls, namely frustum culling 4,096 spheres or boxes, funnel string-pulling, a 512-sphere physics step, closest-hit raycasting, and a 4,096-node transform hierarchy. The README describes the method as "fresh-process blocks, Mann-Whitney U on block medians, Hodges-Lehmann effect size". Each sample starts with a forced GC, and benchmarks that return their output have it snapshotted, so a refactor that changes a result fails the comparison.

That is careful regression testing. It is not a ranking, for two reasons:

• Nothing else is measured. No benchmark imports gl-matrix, three.js, wgpu-matrix or anything else. Every comparison is against math's own saved baseline.
• No results are committed. Runs go to .labs/, which .gitignore excludes.

I searched all 333 commit messages in the tag's history for timings. One has them. Reported in 37519dc (3 September), after the noise tables were flattened into typed arrays, against the full-suite baseline over "8 fresh-process blocks, p < .001 on every bench":

Negative is time saved, math against math. The closest thing to a cross-representation measurement is on an unmerged branch, claude/wasm-matrix-mult-simd-5pqf4h (d8b68e9), a spike that is not in 0.1.0. Its README reports a 4,096-node tree of matrix multiplies: plain arrays about 138 µs, a flat Float32Array about 147 µs, scalar Wasm about 113 µs and SIMD Wasm about 34 µs, averaged over seven runs on a 2.05GHz Xeon container. It adds that "Labs reports this machine as unstable", with about ±7.5% resolution. Read at one significant figure, it backs the SKILL file's claim that plain arrays are not slower than Float32Array, and says SIMD is the only big lever.

So "most optimized" is a claim about process, and a well-run one: each function is regression-gated, allocation-checked and output-checked. The repository has no number that places it against another library.

Tree-shaking is per namespace

"Only pay for what you use" is checkable, so I checked it. Measured with esbuild 0.28.2 (minified ESM), bundling math's source; rollup.config.mjs uses preserveModules precisely so the published dist/ keeps this module graph:

math re-exports each module as a namespace (export * as vec3), and esbuild keeps a whole namespace once any member is used. The package's exports map only offers the eight entrypoints, so a consumer cannot import core/vec3 directly. Under esbuild you pay per namespace, not per function. I did not test Rollup or Rolldown, which handle namespace objects differently.

TypeGPU: schemas first, then shaders in TypeScript

TypeGPU solves a different problem: the bytes between JavaScript and WGSL. Everything starts from a schema. d.f32, d.vec3f, d.mat4x4f, d.struct({ … }) and d.arrayOf(T, n) describe WGSL types, compute their size and alignment, and type the buffers made from them:

// TypeGPU docs, apis/buffers.mdx (abridged)
const Particle = d.struct({
  position: d.vec3f,
  velocity: d.vec3f,
  health: d.f32,
});
 
const buffer = root
  .createBuffer(d.arrayOf(Particle, 100), initial)
  .$usage('storage');
 
const value = await buffer.read(); // typed as an array of { position, velocity, health }

The layout rules are WGSL's, applied for you. The docs' own example: a struct of a vec3i and an f32 has size 16 and alignment 16, and each element of d.arrayOf(d.vec3f, N) occupies 16 bytes, "12 bytes of data + 4 bytes of padding".

The second half is shaders. A function whose body starts with 'use gpu' is picked up at build time by unplugin-typegpu, which stores its syntax tree in a compact format called tinyest next to the function. At runtime tgpu.resolve walks that tree and emits WGSL. The function stays ordinary JavaScript too:

// TypeGPU README
const neighborhood = (a: number, r: number) => {
  'use gpu';
  return d.vec2f(a - r, a + r);
};
 
// #1) Can be called in JS
const range = neighborhood(1.1, 0.5);
 
// #2) Used to generate WGSL
const main = () => {
  'use gpu';
  return neighborhood(1.1, 0.5);
};
const wgsl = tgpu.resolve([main]);

Types are inferred at each call site, not read from the TypeScript annotations, and each new combination of argument types gets its own WGSL overload. The docs flag the sharp edge: a literal that passes Number.isInteger is an integer, so const bar = 1.0 generates 1i. A shell, tgpu.fn([d.f32, d.f32], d.vec2f)(…), pins the signature. Operators on vectors need a companion tool, tsover; without it you call std.add and friends.

Can the same function run on both processors? That is the stated goal: "Our goal is for all functions to have matching behavior on the CPU and GPU". The std functions are built with dualImpl, a JavaScript body plus a WGSL generator. A few, such as the barriers, throw when called on the CPU. TypeGPU also wraps every division in an f32 cast so the two sides agree.

The CPU side has value semantics, because WGSL does. Measured in std/numeric.ts: std.normalize returns a new vector, and so does std.cross. That is the opposite of math. Here is the same quaternion rotation in both, math's on the CPU and TypeGPU's mesh-skinning example on the GPU:

// pmndrs/math, vec3.transformQuat (abridged): locals in, writes to out
let uvx = qy * z - qz * y;            // uv = q.xyz × v, and uvy, uvz
let uuvx = qy * uvz - qz * uvy;       // uuv = q.xyz × uv, and uuvy, uuvz
const w2 = qw * 2;
uvx *= w2;                            // and uvy, uvz
uuvx *= 2;                            // and uuvy, uuvz
out[0] = x + uvx + uuvx;              // and out[1], out[2]
 
// TypeGPU docs, examples/simple/mesh-skinning
const rotateByUnitQuat = (value: d.v3f, quaternion: d.v4f): d.v3f => {
  'use gpu';
  const tangent = 2 * std.cross(quaternion.xyz, value);
  return value + quaternion.w * tangent + std.cross(quaternion.xyz, tangent);
};

Both compute $v + 2w(q \times v) + 2\,q \times (q \times v)$. One shape suits a JIT that hates garbage, the other suits a shader compiler that has no heap.

expand

Do they interoperate?

Not by design. Neither repository mentions the other. TypeGPU's docs have a page on working with wgpu-matrix, and measured, 14 of its 96 examples import wgpu-matrix directly. math's own examples render through gpucat, Isaac Mason's renderer, not TypeGPU.

What does connect them, reasoned from each side's source and docs:

1. Buffers take plain arrays. TypeGPU's .write() accepts "Plain JS array [1, 2, 3]" with "No TypeGPU wrapper allocated", and the docs recommend it for per-frame writes. A math Vec3 is exactly that. A math Mat4 is sixteen column-major numbers, which is WGSL's mat4x4f order, and a math Mat3 is the nine packed numbers TypeGPU accepts for mat3x3f.
2. Raw typed arrays must match GPU padding. If you pack with math's vec3.toBuffer into a Float32Array for d.arrayOf(d.vec3f, N), the stride is 4 floats, not 3: write at i * 4, or use TypeGPU's common.writeSoA, which inserts the padding.
3. TypeGPU vectors pass as math vectors, at a cost. VecBase extends Array, and d.v3f is declared as extending a three-number tuple, so vec3.add(out, a, b) type-checks with TypeGPU vectors. But each index write goes through an f32 cast: Math.fround, and a throw on anything non-finite. A call site fed both plain arrays and TypeGPU instances also stops being monomorphic. TypeGPU matrices are declared as number views, not tuples, so math's Mat4 functions reject them without a cast.
4. math cannot run in a shader. The docs say "Only functions marked with 'use gpu' can be called from within a shader", and math's functions carry no directive.

So the natural split is the one Iwo Plaza described: solve on the CPU in math, render on the GPU in TypeGPU, and meet at buffer.write. I could not find his demo's source. His public TypeGPU projects, the wayfare engine and phoure, depend on wgpu-matrix, not math.

What IK needs from a math library

The limb case is two bones: shoulder to elbow, elbow to hand, a target for the hand. It has a closed form. Put the shoulder at the origin, let $d$ be the distance to the target, clamped into the ring $[\lvert l_1 - l_2 \rvert,\, l_1 + l_2]$ the arm can reach. The law of cosines gives the elbow's bend $\beta$, and the shoulder angle follows:

That is the angle form: one acos, two atan2. A renderer needs positions, and those need no trig at all. The elbow sits on both circles, radius $l_1$ about the shoulder and $l_2$ about the target, a distance $a$ along the shoulder-target line and $h$ off it:

where $\hat{\mathbf{u}}$ is the unit vector to the target and the sign of $h$ picks the side the elbow bends to. Square roots and division are correctly rounded under IEEE-754, so this form gives the same bits in every JavaScript engine, and it is plain arithmetic a shader runs as written. You need the angles only to drive a rig through joint rotations, and even then math's quat.rotationTo between bone directions avoids trig except when the two directions are opposite.

math does not ship this closed form. math/ik is FABRIK, after Aristidou and Lasenby, with Caliko's constraint model: wedge limits in 2-D; ball and hinge joints in 3-D; multi-chain structures. Its defaults are 20 iterations at most, stop within 0.01 of the target, and give up when an iteration that fails to improve moves less than 1e-4. FABRIK is general. It handles a spine or a tail, where no closed form exists. For two bones it has a known weak spot, and I measured it.

Measured, with my own two-bone FABRIK using math/ik's defaults and the near-straight start pose math's fabrik2 bench uses, both bones of length 1 (not math's code):

• Cold start, 1,368 targets on a polar grid, 19 radii by 72 angles, all reachable: a median of 4 iterations, and 272 targets still more than 0.01 away when it stopped. 238 of the 288 targets within 0.4 of the shoulder missed. A folded arm is where FABRIK crawls.
• Warm start, the bench's sweeping target for 10,000 frames, each solve starting from the last pose: a median of 1 iteration and a 95th percentile of 8. 228 frames ended outside 0.01, every one with the target within 0.34 of the shoulder.
• Closed form on the same 10,000 targets: worst error 5.0e-16.

expand

The widget runs the closed form and, as a dashed ghost, my FABRIK from a cold start, on an arm with an upper bone of 1 and a lower bone you set. It starts on a folded pose where FABRIK needs 11 passes. The toggle switches the closed form between two implementations that return identical positions. Out-parameters write into arrays the widget allocated once. The allocating version returns a new [x, y] from every operation, nine per solve. At the default 400 limbs a frame and 60 frames a second, that is 216,000 short-lived arrays a second. Reasoned: V8's young generation makes each one cheap, and escape analysis can remove some once the solver is inlined. Neither is something the code controls, and that is the whole case for writing into arrays you own.

What terrain needs from a math library

A heightfield is a grid of heights $h(x, z)$. Rendering it needs a normal per vertex. Differentiate the surface $y = h(x, z)$: the tangents are $(1, h_x, 0)$ and $(0, h_z, 1)$, and their cross product is

On a grid with spacing $\Delta$, central differences give $h_x \approx (h_R - h_L)/2\Delta$. Multiply through by $2\Delta$ and you get what math's simplex-2d-noise-terrain example computes for every vertex of its 96 by 96 grid, every frame:

// pmndrs/math, examples/src/example-simplex-2d-noise-terrain.ts
const nx = hL - hR;
const ny = 2 * SPACING;
const nz = hD - hU;
const len = Math.hypot(nx, ny, nz) || 1;

Central differences are biased, and the bias has a closed form. For one sine component of wavenumber $k$, the difference quotient is the true slope times

Fine detail loses its slope first, so shading flattens exactly the octaves that make terrain look rough. The example's grid covers 6 units with 96 points, a spacing of 0.063. The widget exaggerates the spacing so the effect is visible. Its terrain is two sines standing in for two noise octaves, the second at 2.3 times the frequency and 0.4 of the amplitude, the example's ratios. At a spacing of 0.3 the second octave keeps 0.809 of its slope and the worst normal is 6.7° off.

expand

math's noise functions return a height and nothing else, so its terrain has to difference. TypeGPU's own noise package takes the other route: @typegpu/noise has a 'use gpu' sampleWithGradient that returns Perlin noise and its analytic gradient in one vec3f. That gives the normal exactly, with no neighbour lookups.

That matters when IK meets terrain. A foot's IK target is a point on the ground; if a vertex shader displaces the ground, the CPU needs the same height there or the foot floats. Reasoned: one 'use gpu' height function called from both sides guarantees that, to within the f32-versus-f64 difference. So does uploading math's CPU heights and sampling that buffer instead of recomputing.

The real library, running in this page

The widgets above are my own code. This one is not. It imports math 0.1.0 from npm, pinned, built from the same c6713e3 the article reads. It runs the library in your browser on the two problems above.

The left panel is a six-bone chain built with fabrik2.createChain2, addBone and addConsecutiveBone. The joint-limit box adds setLocalJoint at 0.7 radians either way, about 40°. Every animation frame calls fabrik2.solve(chain, target), starting from the previous frame's pose. The target sits on the ground, and the ground is the right panel's terrain, sampled along its dashed row. FABRIK knows nothing about that ground, so a joint can dip through it; keeping limbs out of the terrain is the caller's job. The right panel is a 128 by 96 heightfield: simplex2d.create(seed), then fbm per point, shaded with the central-difference normals derived above. One sampler closure per rebuild reads its point from two variables, so the documented fbm call's per-sample closure never appears. The target sweeps by itself unless your system asks for reduced motion, and the loop stops while off screen.

The timings say "your browser, this frame" because that is all they are. Browsers coarsen performance.now(), so a single six-bone solve often reads as below timer resolution. The widget also shows the mean over 120 frames. Measured, in headless Chromium on the shared container I write on: the timer clamps to 100 µs; the 120-frame mean came out between 23 and 38 µs over four runs; a rebuild at eight octaves, 99,592 simplex2d.sample calls, took between 7.7 and 11 ms over six.

It never touches the vec3 rotations that allocate in 0.1.0; fabrik2 imports eight vec2 functions, and the published vec2.rotate writes through locals. The server renders an empty, fixed-size shell, and the library runs only after hydration, so none of its numbers reach the HTML.

Measured, the cost, with esbuild 0.28.2 on this widget's engine: the library adds 9,651 bytes minified, about 3.6 kB gzipped, on top of 7,164 bytes of the widget's own engine code. It splits as the tree-shaking section predicts. fbm is a named export and costs 105 bytes. The eight vec2 functions that fabrik2 imports by name cost 591. simplex2d and the permutation tables behind its seed cost 559 and 1,791. fabrik2 itself is reachable only as a namespace from math/ik, so all 6,559 bytes of it come along. fabrik3 and the other noise modules drop out.

What it adds up to

math is a disciplined library: a gl-matrix API on plain arrays of doubles that the caller owns, with a benchmark suite that gates regressions statistically. "Allocation-free" holds for the core, except three rotations in 0.1.0. "Tree-shakable" holds per namespace under esbuild. "Most optimized" is not something the repository measures: every benchmark compares math with math.

TypeGPU gives WGSL types a TypeScript shape and lets one function body serve JavaScript and the GPU, with value semantics on both. The two meet at the buffer, which takes math's plain arrays as they are: out-parameters where a garbage collector lives, value types where it does not.

Related: Voxel Musou's crowd is the same allocate-once discipline on three.js. OpenDLSS-NR's WebGPU port checks a JavaScript oracle against its WGSL, the CPU/GPU agreement problem from the other end. And the tiny browser models are WebGPU compute with a neural network as the payload.