# pmndrs/math and TypeGPU: out-parameters on the CPU, value types on the GPU
> Satyajit Ghana — Head of Engineering @ Inkers Technology
> canonical: https://ai.thesatyajit.com/articles/pmndrs-math-typegpu
> date: 2026-09-26
> tags: webgpu, typescript, gpu, math, geometry, 3d, performance, benchmarks, open-source, procedural-generation, explainer
On 25 September [pmndrs](https://x.com/pmndrs/status/2103487633502474300) announced
[`math`](https://github.com/pmndrs/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](https://x.com/iwoplaza/status/2103562357301805443),
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:
```ts
// 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.
```ts
// 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:
| | math 0.1.0 | gl-matrix 3.4.4 | wgpu-matrix 3.4.2 | three.js r186 |
|---|---|---|---|---|
| A vector is | plain array of doubles | `Float32Array` (switchable to `Array`) | `Float32Array` (`Float64Array` variants too) | a class with `x`, `y`, `z` fields |
| Add | `add(out, a, b)` | `add(out, a, b)` | `add(a, b, dst?)` | `a.add(b)`, mutates `a` |
| Without an output | not optional | not optional | allocates a new array | clone first, which allocates |
| `mat3` | 9 numbers | 9 numbers | 12, padded like WGSL | 9 numbers |
| Beyond linear algebra | shapes, hulls, noise, springs, IK | nothing | nothing | boxes, frustums, rays, spheres |
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`](https://github.com/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":
| bench | change | bench | change |
|---|---|---|---|
| `perlin2d` | −42% | `simplex4d` | −13% |
| `perlin3d` | −39% | `worley2d` | −41% |
| `simplex2d` | −17% | `worley3d` | −54% |
| `simplex3d` | −20% | `curl2` | −20% |
| `fbm` x5 | −6% | | |
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:
| import | bytes | gzipped |
|---|---|---|
| `vec3.add` via `import { vec3 } from 'math'` | 6,928 | 2,414 |
| `vec3.add` from `core/vec3` directly | 122 | 107 |
| `mat4`, `quat` and `vec3`, three calls, via `math` | 31,346 | 9,303 |
| the same three calls, from the modules directly | 760 | 433 |
| everything `math` exports | 63,578 | 15,921 |
`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:
```ts
// 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:
```ts
// 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:
```ts
// 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.
## 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:
$$
\cos\beta = \frac{d^2 - l_1^2 - l_2^2}{2\,l_1 l_2},
\qquad
\theta_1 = \operatorname{atan2}(t_y, t_x) - \operatorname{atan2}(l_2 \sin\beta,\; l_1 + l_2\cos\beta)
$$
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:
$$
a = \frac{l_1^2 - l_2^2 + d^2}{2d}, \qquad h = \pm\sqrt{l_1^2 - a^2}, \qquad
\mathbf{e} = a\,\hat{\mathbf{u}} + h\,\hat{\mathbf{u}}_\perp
$$
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.
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
$$
\mathbf{n} \propto (-h_x,\; 1,\; -h_z).
$$
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:
```ts
// 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
$$
\frac{h(x+\Delta) - h(x-\Delta)}{2\Delta} = h'(x)\,\frac{\sin k\Delta}{k\Delta}.
$$
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.
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](/articles/voxel-musou#three-hundred-is-a-pool-size) is the same
allocate-once discipline on three.js. [OpenDLSS-NR's WebGPU port](/articles/opendlss-nr#what-same-output-is-measured-to-mean)
checks a JavaScript oracle against its WGSL, the CPU/GPU agreement problem from the other end. And
[the tiny browser models](/articles/tiny-browser-models) are WebGPU compute with a neural network as
the payload.