diff --git a/bench/octave/README.md b/bench/octave/README.md
new file mode 100644
index 0000000..02dbb53
--- /dev/null
+++ b/bench/octave/README.md
@@ -0,0 +1,14 @@
+# Octave cross-check
+
+Independent numerical verification of every public function in
+`Graphmath.{Vec2,Vec3,Vec4,Mat22,Mat33,Mat44,Quatern}` against GNU Octave.
+
+```
+bench/octave/run.sh            # needs mix + octave-cli on PATH
+```
+
+`gen_cases.exs` calls each function with seeded random inputs (plus documented
+error and edge cases) and writes the arguments and results to an Octave data
+file. `check.m` recomputes each result from the documented conventions and
+prints the maximum relative error per function, flagging anything above 1e-9.
+Cases named `err.*` / `edge.*` are listed in the data file for inspection only.
diff --git a/bench/octave/check.m b/bench/octave/check.m
new file mode 100644
index 0000000..390714c
--- /dev/null
+++ b/bench/octave/check.m
@@ -0,0 +1,204 @@
+1;
+% Octave reference checker for graphmath. Usage:
+%   mix run bench/octave/gen_cases.exs /tmp/cases.m
+%   octave-cli --no-gui bench/octave/check.m /tmp/cases.m
+% Recomputes every case from the documented conventions (row-major tuples,
+% row-vector transforms => expected matrix is R^T, right-handed CCW rotations,
+% OpenGL depth, Hamilton product on {w,x,y,z}) and prints per-function max error.
+% Reference implementations. Matrices arrive as flat row-major tuples.
+function M = mat(v) n = sqrt(numel(v)); M = reshape(v, n, n)'; end
+function v = flat(M) v = reshape(M', 1, []); end
+function u = nrm(v) u = v / norm(v); end
+% quaternion {w,x,y,z} Hamilton product
+function r = qmul(p, q)
+  w = p(1)*q(1) - dot(p(2:4), q(2:4));
+  v = p(1)*q(2:4) + q(1)*p(2:4) + cross(p(2:4), q(2:4));
+  r = [w v];
+end
+% column-vector rotation matrix from unit quaternion
+function R = q2R(q)
+  w=q(1);x=q(2);y=q(3);z=q(4);
+  R = [1-2*(y*y+z*z), 2*(x*y-w*z), 2*(x*z+w*y);
+       2*(x*y+w*z), 1-2*(x*x+z*z), 2*(y*z-w*x);
+       2*(x*z-w*y), 2*(y*z+w*x), 1-2*(x*x+y*y)];
+end
+function R = axang(k, th)  % Rodrigues, column convention
+  K = [0 -k(3) k(2); k(3) 0 -k(1); -k(2) k(1) 0];
+  R = eye(3) + sin(th)*K + (1-cos(th))*K*K;
+end
+function M = orient_ref(p, f, up)
+  f = nrm(f); right = nrm(cross(f, nrm(up))); u = cross(right, f);
+  M = [right 0; u 0; -f 0; p 1];
+end
+function e = ref(c)
+  a = c.a; e = [];
+  nm = c.name; p = strsplit(nm, '.'); mod = p{1}; fn = p{2};
+  switch mod
+  case {'v2','v3','v4'}
+    switch fn
+    case 'add', e = a{1}+a{2};
+    case 'subtract', e = a{1}-a{2};
+    case 'multiply', e = a{1}.*a{2};
+    case 'scale', e = a{1}*a{2};
+    case 'dot', e = dot(a{1},a{2});
+    case 'cross', e = cross(a{1},a{2});
+    case 'length', e = norm(a{1});
+    case 'length_squared', e = dot(a{1},a{1});
+    case 'length_manhattan', e = sum(abs(a{1}));
+    case 'chebyshev_distance', e = max(abs(a{1}-a{2}));
+    case 'minkowski_distance', e = sum(abs(a{1}-a{2}).^a{3})^(1/a{3});
+    case 'p_norm', e = sum(abs(a{1}).^a{2})^(1/a{2});
+    case 'lerp', e = (1-a{3})*a{1} + a{3}*a{2};
+    case 'negate', e = -a{1};
+    case 'normalize', e = nrm(a{1});
+    case 'project', e = dot(a{1},a{2})/dot(a{2},a{2}) * a{2};
+    case 'perp', e = [-a{1}(2) a{1}(1)];
+    case 'perp_prod', e = a{1}(1)*a{2}(2) - a{2}(1)*a{1}(2);
+    case 'rotate'
+      if strcmp(mod,'v2'), th=a{2}; e = ([cos(th) -sin(th); sin(th) cos(th)]*a{1}')';
+      else e = (axang(a{2}, a{3})*a{1}')'; end
+    case 'scalar_triple', e = dot(a{1}, cross(a{2}, a{3}));
+    case 'weighted_sum', e = a{1}*a{2} + a{3}*a{4};
+    case 'near', e = norm(a{1}-a{2}) < a{3};
+    case 'equal3', e = all(abs(a{1}-a{2}) <= a{3});
+    case 'equal', e = all(a{1}==a{2});
+    case {'create1','create2','create3','create4'}, e = [a{:}];
+    case 'from_point3', e = [a{1} 1];
+    case 'from_direction3', e = [a{1} 0];
+    case 'random_circle', e = 1; c.out = norm(c.out);
+    case 'random_sphere', e = 1; c.out = norm(c.out);
+    case {'random_disc','random_ball'}, e = 1; c.out = double(norm(c.out) <= 1);
+    case 'random_box', e = 1; c.out = double(all(c.out >= 0 & c.out <= 1));
+    end
+  case {'m22','m33','m44'}
+    n = mod(2) - '0';
+    switch fn
+    case 'identity', e = flat(eye(n));
+    case 'zero', e = zeros(1,n*n);
+    case 'add', e = a{1}+a{2};
+    case 'subtract', e = a{1}-a{2};
+    case 'scale', e = a{1}*a{2};
+    case 'multiply', e = flat(mat(a{1})*mat(a{2}));
+    case 'multiply_transpose', e = flat(mat(a{1})*mat(a{2})');
+    case 'apply', e = (mat(a{1})*a{2}')';
+    case 'apply_left', e = a{1}*mat(a{2});
+    case 'apply_transpose', e = (mat(a{1})'*a{2}')';
+    case 'apply_left_transpose', e = a{1}*mat(a{2})';
+    case 'transform_vector'
+      if n==2, e = a{2}*mat(a{1}); else v = [a{2} 0]*mat(a{1}); e = v(1:n-1); end
+    case 'transform_point', v = [a{2} 1]*mat(a{1}); e = v(1:n-1);
+    case 'at', M = mat(a{1}); e = M(a{2}+1, a{3}+1);
+    case 'submatrix', M = mat(a{1}); M(a{2}+1,:)=[]; M(:,a{3}+1)=[]; e = flat(M);
+    case 'cofactor', M = mat(a{1}); M(a{2}+1,:)=[]; M(:,a{3}+1)=[]; e = (-1)^(a{2}+a{3})*det(M);
+    case {'row0','row1','row2','row3'}, M = mat(a{1}); e = M(fn(4)-'0'+1, :);
+    case {'column0','column1','column2','column3'}, M = mat(a{1}); e = M(:, fn(7)-'0'+1)';
+    case 'diag', e = diag(mat(a{1}))';
+    case 'determinant', e = det(mat(a{1}));
+    case 'trace', e = trace(mat(a{1}));
+    case 'inverse', e = flat(inv(mat(a{1})));
+    case 'round', e = round(a{1}*10^a{2})/10^a{2};
+    % ---- constructors: row-vector convention => M = R' (R column-convention)
+    case 'make_rotate'  % about +Z CCW
+      th=a{1}; R = [cos(th) -sin(th); sin(th) cos(th)];
+      if n==3, R = blkdiag(R,1); end; e = flat(R');
+    case 'make_rotate_x', R = blkdiag(axang([1 0 0], a{1}), 1); e = flat(R');
+    case 'make_rotate_y', R = blkdiag(axang([0 1 0], a{1}), 1); e = flat(R');
+    case 'make_rotate_z', R = blkdiag(axang([0 0 1], a{1}), 1); e = flat(R');
+    case 'make_translate', M = eye(n); M(n, 1:n-1) = [a{:}]; e = flat(M);
+    case 'make_scale1', e = flat(eye(n)*a{1});
+    case {'make_scale2','make_scale3','make_scale4'}, e = flat(diag([a{:}]));
+    case {'make_reflect','make_reflect_3d','make_reflect_2d'}
+      nv = a{1}; d = numel(nv); H = eye(d) - 2*(nv'*nv)/dot(nv,nv);
+      if numel(a) == 2  % affine: p' = p - 2*((n.p - off)/|n|^2) n
+        M = eye(d+1); M(1:d,1:d) = H'; M(d+1,1:d) = 2*a{2}*nv/dot(nv,nv);
+      else M = H'; end
+      e = flat(M);
+    case {'make_shear_x','make_shear_x_2d','make_shear_x_3d'}
+      S = eye(n); if n==2 || strcmp(fn,'make_shear_x_2d'), S(1,2)=a{1}; else S(1,2)=a{1}; S(1,3)=a{2}; end
+      e = flat(S');
+    case {'make_shear_y','make_shear_y_2d','make_shear_y_3d'}
+      S = eye(n); if n==2 || strcmp(fn,'make_shear_y_2d'), S(2,1)=a{1}; else S(2,1)=a{1}; S(2,3)=a{2}; end
+      e = flat(S');
+    case {'make_shear_z','make_shear_z_3d'}
+      S = eye(n); S(3,1)=a{1}; S(3,2)=a{2}; e = flat(S');
+    case 'orthonormalize'
+      A = mat(a{1}); [Q,R] = qr(A'); Q = Q*diag(sign(diag(R))); e = flat(Q');
+    case 'orient', e = flat(orient_ref(a{1}, a{2}, a{3}));
+    case 'look_at', e = flat(inv(orient_ref(a{1}, a{2}-a{1}, a{3})));
+    case 'make_billboard', e = flat(orient_ref(a{1}, a{2}-a{1}, a{3}));
+    case 'make_billboard_axis'
+      y = nrm(a{3}); d = a{2}-a{1}; d = d - dot(d,y)*y; z = -nrm(d); x = cross(y, z);
+      e = flat([x 0; y 0; z 0; a{1} 1]);
+    case 'perspective'  % gluPerspective, column convention, then transpose
+      fovy=a{1}; asp=a{2}; nr=a{3}; fr=a{4}; f = 1/tan(fovy/2);
+      P = [f/asp 0 0 0; 0 f 0 0; 0 0 (fr+nr)/(nr-fr) 2*fr*nr/(nr-fr); 0 0 -1 0];
+      e = flat(P');
+    case 'ortho'  % glOrtho
+      l=a{1};r=a{2};b=a{3};t=a{4};nr=a{5};fr=a{6};
+      P = [2/(r-l) 0 0 -(r+l)/(r-l); 0 2/(t-b) 0 -(t+b)/(t-b); 0 0 -2/(fr-nr) -(fr+nr)/(fr-nr); 0 0 0 1];
+      e = flat(P');
+    end
+  case 'q'
+    switch fn
+    case 'identity', e = [1 0 0 0];
+    case 'zero', e = [0 0 0 0];
+    case {'create','from_list'}, e = [a{:}];
+    case 'add', e = a{1}+a{2};
+    case 'subtract', e = a{1}-a{2};
+    case 'multiply', e = qmul(a{1}, a{2});
+    case 'scale', e = a{1}*a{2};
+    case 'dot', e = dot(a{1},a{2});
+    case 'norm', e = norm(a{1});
+    case {'normalize','normalize_strict'}, e = nrm(a{1});
+    case 'conjugate', e = a{1}.*[1 -1 -1 -1];
+    case 'inverse', e = a{1}.*[1 -1 -1 -1]/dot(a{1},a{1});
+    case 'equal', e = abs(dot(a{1},a{2})) >= 1;
+    case 'equal_elements', e = all(abs(a{1}-a{2}) <= a{3});
+    case 'from_axis_angle', e = [cos(a{1}/2) sin(a{1}/2)*a{2}];
+    case 'to_rotation_matrix_33', e = flat(q2R(a{1})');
+    case 'to_rotation_matrix_44', e = flat(blkdiag(q2R(a{1})', 1));
+    case 'transform_vector', e = (q2R(a{1})*a{2}')';
+    case 'from_rotation_matrix'
+      % compare as orientation: q2R(out)' must equal input matrix
+      e = a{1}; c.out = flat(q2R(nrm(c.out))');
+    case 'slerp'
+      p = a{1}; q = a{2}; t = a{3}; d = dot(p,q);
+      if d < 0, q = -q; d = -d; end
+      if d > 1 - 1e-9, e = nrm((1-t)*p + t*q);
+      else th = acos(min(d,1)); e = (sin((1-t)*th)*p + sin(t*th)*q)/sin(th); end
+      if dot(e, c.out) < 0, c.out = -c.out; end  % same orientation, either sign
+    case 'integrate'
+      q = a{1}; om = a{2}; dt = a{3}; th = norm(om)*dt;
+      if th == 0, dq = [1 0 0 0]; else dq = [cos(th/2) sin(th/2)*om/norm(om)]; end
+      e = nrm(qmul(dq, q));
+    case 'random', e = 1; c.out = norm(c.out);
+    case {'get_roll_z','get_pitch_x','get_yaw_y'}, e = a{1};
+    case 'euler'
+      % Ogre convention: roll about Z, pitch about X, yaw about Y; check each
+      % extractor against the quaternion's rotation matrix.
+      R = q2R(nrm(a{1}));
+      e = [atan2(R(2,1), R(1,1)), atan2(R(3,2), R(2,2)), atan2(R(1,3), R(3,3))];
+    end
+  end
+  e = double(e); c.out = double(c.out);
+  if isempty(e), error('no reference for %s', nm); end
+  if numel(e) ~= numel(c.out), error('shape mismatch for %s: %d vs %d', nm, numel(e), numel(c.out)); end
+  e = max(abs(e(:) - c.out(:)) ./ max(1, abs(e(:))));  % relative-ish error
+end
+
+function s = ifelse_str(b) if b, s = "  <-- FAIL"; else s = ""; end; end
+run(argv(){1});
+names = cellfun(@(c) c.name, C, 'UniformOutput', false);
+[u, ~, idx] = unique(names);
+worst = zeros(size(u));
+for k = 1:numel(C)
+  if strncmp(C{k}.name,'err.',4) || strncmp(C{k}.name,'edge.',5), continue; end
+  if ~isequal(C{k}.err, 0), fprintf('RAISED %s %s\n', C{k}.name, C{k}.err); continue; end
+  err = ref(C{k});
+  worst(idx(k)) = max(worst(idx(k)), err);
+  if err > 1e-9
+    fprintf('MISMATCH %-28s err=%.3e case=%d\n', C{k}.name, err, k);
+  end
+end
+fprintf('\n%-28s %s\n', 'function', 'max err');
+for i = 1:numel(u), fprintf('%-28s %.2e %s\n', u{i}, worst(i), ifelse_str(worst(i) > 1e-9)); end
diff --git a/bench/octave/gen_cases.exs b/bench/octave/gen_cases.exs
new file mode 100644
index 0000000..c77ef5f
--- /dev/null
+++ b/bench/octave/gen_cases.exs
@@ -0,0 +1,553 @@
+# Usage: mix run bench/octave/gen_cases.exs /tmp/cases.m  (then run check.m on it)
+# Drives every public graphmath function with seeded random inputs and writes
+# an Octave data file: C{k} = struct(name, args{...}, out, err)
+alias Graphmath.{Vec2, Vec3, Vec4, Mat22, Mat33, Mat44, Quatern}
+
+defmodule Gen do
+  def start(path) do
+    :rand.seed(:exsss, {20_260_916, 1, 2})
+    {:ok, io} = File.open(path, [:write])
+    Process.put(:io, io)
+    Process.put(:k, 0)
+    IO.write(io, "C = {};\n")
+  end
+
+  def finish, do: File.close(Process.get(:io))
+
+  def r(), do: :rand.uniform() * 10.0 - 5.0
+  def rp(), do: :rand.uniform() * 4.0 + 0.5
+  def ang(), do: :rand.uniform() * 2 * :math.pi() - :math.pi()
+  def v2(), do: {r(), r()}
+  def v3(), do: {r(), r(), r()}
+  def v4(), do: {r(), r(), r(), r()}
+  def m2(), do: List.to_tuple(for _ <- 1..4, do: r())
+  def m3(), do: List.to_tuple(for _ <- 1..9, do: r())
+  def m4(), do: List.to_tuple(for _ <- 1..16, do: r())
+  def unit3(), do: Vec3.normalize(v3())
+  def uq(), do: Quatern.normalize(v4())
+
+  def fmt(x) when is_float(x), do: Float.to_string(x)
+  def fmt(x) when is_integer(x), do: Integer.to_string(x)
+  def fmt(true), do: "1"
+  def fmt(false), do: "0"
+  def fmt(x) when is_tuple(x), do: "[" <> Enum.map_join(Tuple.to_list(x), " ", &fmt/1) <> "]"
+  def fmt(x) when is_list(x), do: "[" <> Enum.map_join(x, " ", &fmt/1) <> "]"
+
+  def emit(name, args, fun) do
+    k = Process.get(:k) + 1
+    Process.put(:k, k)
+    io = Process.get(:io)
+
+    {out, err} =
+      try do
+        {apply(fun, args), 0}
+      rescue
+        e -> {[], "#{inspect(e.__struct__)}"}
+      end
+
+    a = Enum.map_join(args, ", ", &fmt/1)
+    IO.write(io, "C{#{k}} = struct('name','#{name}','a',{{#{a}}},'out',#{fmt(out)},'err'")
+    if err == 0, do: IO.write(io, ",0);\n"), else: IO.write(io, ",'#{err}');\n")
+  end
+
+  def times(n, f), do: for(_ <- 1..n, do: f.())
+end
+
+defmodule Main do
+  import Gen
+
+  def run() do
+    start(System.argv() |> hd())
+    n = 6
+
+    # ---------------- Vec2
+    times(n, fn -> emit("v2.add", [v2(), v2()], &Vec2.add/2) end)
+    times(n, fn -> emit("v2.subtract", [v2(), v2()], &Vec2.subtract/2) end)
+    times(n, fn -> emit("v2.multiply", [v2(), v2()], &Vec2.multiply/2) end)
+    times(n, fn -> emit("v2.scale", [v2(), r()], &Vec2.scale/2) end)
+    times(n, fn -> emit("v2.dot", [v2(), v2()], &Vec2.dot/2) end)
+    times(n, fn -> emit("v2.length", [v2()], &Vec2.length/1) end)
+    times(n, fn -> emit("v2.length_squared", [v2()], &Vec2.length_squared/1) end)
+    times(n, fn -> emit("v2.length_manhattan", [v2()], &Vec2.length_manhattan/1) end)
+    times(n, fn -> emit("v2.chebyshev_distance", [v2(), v2()], &Vec2.chebyshev_distance/2) end)
+
+    times(n, fn ->
+      emit("v2.minkowski_distance", [v2(), v2(), rp()], &Vec2.minkowski_distance/3)
+    end)
+
+    times(n, fn -> emit("v2.p_norm", [v2(), rp()], &Vec2.p_norm/2) end)
+    times(n, fn -> emit("v2.lerp", [v2(), v2(), :rand.uniform()], &Vec2.lerp/3) end)
+    times(n, fn -> emit("v2.negate", [v2()], &Vec2.negate/1) end)
+    times(n, fn -> emit("v2.normalize", [v2()], &Vec2.normalize/1) end)
+    times(n, fn -> emit("v2.project", [v2(), v2()], &Vec2.project/2) end)
+    times(n, fn -> emit("v2.perp", [v2()], &Vec2.perp/1) end)
+    times(n, fn -> emit("v2.perp_prod", [v2(), v2()], &Vec2.perp_prod/2) end)
+    times(n, fn -> emit("v2.rotate", [v2(), ang()], &Vec2.rotate/2) end)
+    times(n, fn -> emit("v2.weighted_sum", [r(), v2(), r(), v2()], &Vec2.weighted_sum/4) end)
+    times(n, fn -> emit("v2.near", [v2(), v2(), rp()], &Vec2.near/3) end)
+    times(n, fn -> emit("v2.equal3", [v2(), v2(), rp()], &Vec2.equal/3) end)
+
+    times(n, fn ->
+      a = v2()
+      emit("v2.equal", [a, a], &Vec2.equal/2)
+    end)
+
+    times(n, fn -> emit("v2.equal", [v2(), v2()], &Vec2.equal/2) end)
+    times(n, fn -> emit("v2.create2", [r(), r()], &Vec2.create/2) end)
+    times(n, fn -> emit("v2.create1", [[r(), r()]], &Vec2.create/1) end)
+    times(n, fn -> emit("v2.random_circle", [], &Vec2.random_circle/0) end)
+    times(n, fn -> emit("v2.random_disc", [], &Vec2.random_disc/0) end)
+    times(n, fn -> emit("v2.random_box", [], &Vec2.random_box/0) end)
+
+    # ---------------- Vec3
+    times(n, fn -> emit("v3.add", [v3(), v3()], &Vec3.add/2) end)
+    times(n, fn -> emit("v3.subtract", [v3(), v3()], &Vec3.subtract/2) end)
+    times(n, fn -> emit("v3.multiply", [v3(), v3()], &Vec3.multiply/2) end)
+    times(n, fn -> emit("v3.scale", [v3(), r()], &Vec3.scale/2) end)
+    times(n, fn -> emit("v3.dot", [v3(), v3()], &Vec3.dot/2) end)
+    times(n, fn -> emit("v3.cross", [v3(), v3()], &Vec3.cross/2) end)
+    times(n, fn -> emit("v3.length", [v3()], &Vec3.length/1) end)
+    times(n, fn -> emit("v3.length_squared", [v3()], &Vec3.length_squared/1) end)
+    times(n, fn -> emit("v3.length_manhattan", [v3()], &Vec3.length_manhattan/1) end)
+    times(n, fn -> emit("v3.chebyshev_distance", [v3(), v3()], &Vec3.chebyshev_distance/2) end)
+
+    times(n, fn ->
+      emit("v3.minkowski_distance", [v3(), v3(), rp()], &Vec3.minkowski_distance/3)
+    end)
+
+    times(n, fn -> emit("v3.p_norm", [v3(), rp()], &Vec3.p_norm/2) end)
+    times(n, fn -> emit("v3.lerp", [v3(), v3(), :rand.uniform()], &Vec3.lerp/3) end)
+    times(n, fn -> emit("v3.negate", [v3()], &Vec3.negate/1) end)
+    times(n, fn -> emit("v3.normalize", [v3()], &Vec3.normalize/1) end)
+    times(n, fn -> emit("v3.rotate", [v3(), unit3(), ang()], &Vec3.rotate/3) end)
+    times(n, fn -> emit("v3.scalar_triple", [v3(), v3(), v3()], &Vec3.scalar_triple/3) end)
+    times(n, fn -> emit("v3.weighted_sum", [r(), v3(), r(), v3()], &Vec3.weighted_sum/4) end)
+    times(n, fn -> emit("v3.near", [v3(), v3(), rp()], &Vec3.near/3) end)
+    times(n, fn -> emit("v3.equal3", [v3(), v3(), rp()], &Vec3.equal/3) end)
+    times(n, fn -> emit("v3.create3", [r(), r(), r()], &Vec3.create/3) end)
+    times(n, fn -> emit("v3.create1", [[r(), r(), r()]], &Vec3.create/1) end)
+    times(n, fn -> emit("v3.random_sphere", [], &Vec3.random_sphere/0) end)
+    times(n, fn -> emit("v3.random_ball", [], &Vec3.random_ball/0) end)
+    times(n, fn -> emit("v3.random_box", [], &Vec3.random_box/0) end)
+
+    # ---------------- Vec4
+    times(n, fn -> emit("v4.add", [v4(), v4()], &Vec4.add/2) end)
+    times(n, fn -> emit("v4.subtract", [v4(), v4()], &Vec4.subtract/2) end)
+    times(n, fn -> emit("v4.multiply", [v4(), v4()], &Vec4.multiply/2) end)
+    times(n, fn -> emit("v4.scale", [v4(), r()], &Vec4.scale/2) end)
+    times(n, fn -> emit("v4.dot", [v4(), v4()], &Vec4.dot/2) end)
+    times(n, fn -> emit("v4.length", [v4()], &Vec4.length/1) end)
+    times(n, fn -> emit("v4.length_squared", [v4()], &Vec4.length_squared/1) end)
+    times(n, fn -> emit("v4.length_manhattan", [v4()], &Vec4.length_manhattan/1) end)
+    times(n, fn -> emit("v4.chebyshev_distance", [v4(), v4()], &Vec4.chebyshev_distance/2) end)
+
+    times(n, fn ->
+      emit("v4.minkowski_distance", [v4(), v4(), rp()], &Vec4.minkowski_distance/3)
+    end)
+
+    times(n, fn -> emit("v4.p_norm", [v4(), rp()], &Vec4.p_norm/2) end)
+    times(n, fn -> emit("v4.lerp", [v4(), v4(), :rand.uniform()], &Vec4.lerp/3) end)
+    times(n, fn -> emit("v4.negate", [v4()], &Vec4.negate/1) end)
+    times(n, fn -> emit("v4.normalize", [v4()], &Vec4.normalize/1) end)
+    times(n, fn -> emit("v4.project", [v4(), v4()], &Vec4.project/2) end)
+    times(n, fn -> emit("v4.weighted_sum", [r(), v4(), r(), v4()], &Vec4.weighted_sum/4) end)
+    times(n, fn -> emit("v4.near", [v4(), v4(), rp()], &Vec4.near/3) end)
+    times(n, fn -> emit("v4.equal3", [v4(), v4(), rp()], &Vec4.equal/3) end)
+    times(n, fn -> emit("v4.create4", [r(), r(), r(), r()], &Vec4.create/4) end)
+    times(n, fn -> emit("v4.create1", [[r(), r(), r(), r()]], &Vec4.create/1) end)
+    times(n, fn -> emit("v4.from_point3", [v3()], &Vec4.from_point3/1) end)
+    times(n, fn -> emit("v4.from_direction3", [v3()], &Vec4.from_direction3/1) end)
+
+    # ---------------- Mat22
+    emit("m22.identity", [], &Mat22.identity/0)
+    emit("m22.zero", [], &Mat22.zero/0)
+    times(n, fn -> emit("m22.add", [m2(), m2()], &Mat22.add/2) end)
+    times(n, fn -> emit("m22.subtract", [m2(), m2()], &Mat22.subtract/2) end)
+    times(n, fn -> emit("m22.scale", [m2(), r()], &Mat22.scale/2) end)
+    times(n, fn -> emit("m22.multiply", [m2(), m2()], &Mat22.multiply/2) end)
+    times(n, fn -> emit("m22.multiply_transpose", [m2(), m2()], &Mat22.multiply_transpose/2) end)
+    times(n, fn -> emit("m22.apply", [m2(), v2()], &Mat22.apply/2) end)
+    times(n, fn -> emit("m22.apply_left", [v2(), m2()], &Mat22.apply_left/2) end)
+    times(n, fn -> emit("m22.apply_transpose", [m2(), v2()], &Mat22.apply_transpose/2) end)
+
+    times(n, fn ->
+      emit("m22.apply_left_transpose", [v2(), m2()], &Mat22.apply_left_transpose/2)
+    end)
+
+    times(n, fn -> emit("m22.transform_vector", [m2(), v2()], &Mat22.transform_vector/2) end)
+    for i <- 0..1, j <- 0..1, do: emit("m22.at", [m2(), i, j], &Mat22.at/3)
+    times(n, fn -> emit("m22.row0", [m2()], &Mat22.row0/1) end)
+    times(n, fn -> emit("m22.row1", [m2()], &Mat22.row1/1) end)
+    times(n, fn -> emit("m22.column0", [m2()], &Mat22.column0/1) end)
+    times(n, fn -> emit("m22.column1", [m2()], &Mat22.column1/1) end)
+    times(n, fn -> emit("m22.diag", [m2()], &Mat22.diag/1) end)
+    times(n, fn -> emit("m22.determinant", [m2()], &Mat22.determinant/1) end)
+    times(n, fn -> emit("m22.trace", [m2()], &Mat22.trace/1) end)
+    times(n, fn -> emit("m22.inverse", [m2()], &Mat22.inverse/1) end)
+    times(n, fn -> emit("m22.round", [m2(), 2], &Mat22.round/2) end)
+    times(n, fn -> emit("m22.make_rotate", [ang()], &Mat22.make_rotate/1) end)
+    times(n, fn -> emit("m22.make_scale1", [r()], &Mat22.make_scale/1) end)
+    times(n, fn -> emit("m22.make_scale2", [r(), r()], &Mat22.make_scale/2) end)
+    times(n, fn -> emit("m22.make_reflect", [v2()], &Mat22.make_reflect/1) end)
+    times(n, fn -> emit("m22.make_shear_x", [r()], &Mat22.make_shear_x/1) end)
+    times(n, fn -> emit("m22.make_shear_y", [r()], &Mat22.make_shear_y/1) end)
+
+    # ---------------- Mat33
+    emit("m33.identity", [], &Mat33.identity/0)
+    emit("m33.zero", [], &Mat33.zero/0)
+    times(n, fn -> emit("m33.add", [m3(), m3()], &Mat33.add/2) end)
+    times(n, fn -> emit("m33.subtract", [m3(), m3()], &Mat33.subtract/2) end)
+    times(n, fn -> emit("m33.scale", [m3(), r()], &Mat33.scale/2) end)
+    times(n, fn -> emit("m33.multiply", [m3(), m3()], &Mat33.multiply/2) end)
+    times(n, fn -> emit("m33.multiply_transpose", [m3(), m3()], &Mat33.multiply_transpose/2) end)
+    times(n, fn -> emit("m33.apply", [m3(), v3()], &Mat33.apply/2) end)
+    times(n, fn -> emit("m33.apply_left", [v3(), m3()], &Mat33.apply_left/2) end)
+    times(n, fn -> emit("m33.apply_transpose", [m3(), v3()], &Mat33.apply_transpose/2) end)
+
+    times(n, fn ->
+      emit("m33.apply_left_transpose", [v3(), m3()], &Mat33.apply_left_transpose/2)
+    end)
+
+    times(n, fn -> emit("m33.transform_point", [m3(), v2()], &Mat33.transform_point/2) end)
+    times(n, fn -> emit("m33.transform_vector", [m3(), v2()], &Mat33.transform_vector/2) end)
+    for i <- 0..2, j <- 0..2, do: emit("m33.at", [m3(), i, j], &Mat33.at/3)
+    for i <- 0..2, j <- 0..2, do: emit("m33.submatrix", [m3(), i, j], &Mat33.submatrix/3)
+    for i <- 0..2, j <- 0..2, do: emit("m33.cofactor", [m3(), i, j], &Mat33.cofactor/3)
+    times(n, fn -> emit("m33.row0", [m3()], &Mat33.row0/1) end)
+    times(n, fn -> emit("m33.row1", [m3()], &Mat33.row1/1) end)
+    times(n, fn -> emit("m33.row2", [m3()], &Mat33.row2/1) end)
+    times(n, fn -> emit("m33.column0", [m3()], &Mat33.column0/1) end)
+    times(n, fn -> emit("m33.column1", [m3()], &Mat33.column1/1) end)
+    times(n, fn -> emit("m33.column2", [m3()], &Mat33.column2/1) end)
+    times(n, fn -> emit("m33.diag", [m3()], &Mat33.diag/1) end)
+    times(n, fn -> emit("m33.determinant", [m3()], &Mat33.determinant/1) end)
+    times(n, fn -> emit("m33.trace", [m3()], &Mat33.trace/1) end)
+    times(n, fn -> emit("m33.inverse", [m3()], &Mat33.inverse/1) end)
+    times(n, fn -> emit("m33.round", [m3(), 3], &Mat33.round/2) end)
+    times(n, fn -> emit("m33.make_rotate", [ang()], &Mat33.make_rotate/1) end)
+    times(n, fn -> emit("m33.make_translate", [r(), r()], &Mat33.make_translate/2) end)
+    times(n, fn -> emit("m33.make_scale1", [r()], &Mat33.make_scale/1) end)
+    times(n, fn -> emit("m33.make_scale3", [r(), r(), r()], &Mat33.make_scale/3) end)
+    times(n, fn -> emit("m33.make_reflect_2d", [v2(), r()], &Mat33.make_reflect_2d/2) end)
+    times(n, fn -> emit("m33.make_reflect_3d", [v3()], &Mat33.make_reflect_3d/1) end)
+    times(n, fn -> emit("m33.make_shear_x_2d", [r()], &Mat33.make_shear_x_2d/1) end)
+    times(n, fn -> emit("m33.make_shear_y_2d", [r()], &Mat33.make_shear_y_2d/1) end)
+    times(n, fn -> emit("m33.make_shear_x_3d", [r(), r()], &Mat33.make_shear_x_3d/2) end)
+    times(n, fn -> emit("m33.make_shear_y_3d", [r(), r()], &Mat33.make_shear_y_3d/2) end)
+    times(n, fn -> emit("m33.make_shear_z_3d", [r(), r()], &Mat33.make_shear_z_3d/2) end)
+    times(n, fn -> emit("m33.orthonormalize", [m3()], &Mat33.orthonormalize/1) end)
+
+    # ---------------- Mat44
+    emit("m44.identity", [], &Mat44.identity/0)
+    emit("m44.zero", [], &Mat44.zero/0)
+    times(n, fn -> emit("m44.add", [m4(), m4()], &Mat44.add/2) end)
+    times(n, fn -> emit("m44.subtract", [m4(), m4()], &Mat44.subtract/2) end)
+    times(n, fn -> emit("m44.scale", [m4(), r()], &Mat44.scale/2) end)
+    times(n, fn -> emit("m44.multiply", [m4(), m4()], &Mat44.multiply/2) end)
+    times(n, fn -> emit("m44.multiply_transpose", [m4(), m4()], &Mat44.multiply_transpose/2) end)
+    times(n, fn -> emit("m44.apply", [m4(), v4()], &Mat44.apply/2) end)
+    times(n, fn -> emit("m44.apply_left", [v4(), m4()], &Mat44.apply_left/2) end)
+    times(n, fn -> emit("m44.apply_transpose", [m4(), v4()], &Mat44.apply_transpose/2) end)
+
+    times(n, fn ->
+      emit("m44.apply_left_transpose", [v4(), m4()], &Mat44.apply_left_transpose/2)
+    end)
+
+    times(n, fn -> emit("m44.transform_point", [m4(), v3()], &Mat44.transform_point/2) end)
+    times(n, fn -> emit("m44.transform_vector", [m4(), v3()], &Mat44.transform_vector/2) end)
+    for i <- 0..3, j <- 0..3, do: emit("m44.at", [m4(), i, j], &Mat44.at/3)
+    for i <- 0..3, j <- 0..3, do: emit("m44.submatrix", [m4(), i, j], &Mat44.submatrix/3)
+    for i <- 0..3, j <- 0..3, do: emit("m44.cofactor", [m4(), i, j], &Mat44.cofactor/3)
+    times(n, fn -> emit("m44.row0", [m4()], &Mat44.row0/1) end)
+    times(n, fn -> emit("m44.row1", [m4()], &Mat44.row1/1) end)
+    times(n, fn -> emit("m44.row2", [m4()], &Mat44.row2/1) end)
+    times(n, fn -> emit("m44.row3", [m4()], &Mat44.row3/1) end)
+    times(n, fn -> emit("m44.column0", [m4()], &Mat44.column0/1) end)
+    times(n, fn -> emit("m44.column1", [m4()], &Mat44.column1/1) end)
+    times(n, fn -> emit("m44.column2", [m4()], &Mat44.column2/1) end)
+    times(n, fn -> emit("m44.column3", [m4()], &Mat44.column3/1) end)
+    times(n, fn -> emit("m44.diag", [m4()], &Mat44.diag/1) end)
+    times(n, fn -> emit("m44.determinant", [m4()], &Mat44.determinant/1) end)
+    times(n, fn -> emit("m44.trace", [m4()], &Mat44.trace/1) end)
+    times(n, fn -> emit("m44.inverse", [m4()], &Mat44.inverse/1) end)
+    times(n, fn -> emit("m44.round", [m4(), 3], &Mat44.round/2) end)
+    times(n, fn -> emit("m44.make_rotate_x", [ang()], &Mat44.make_rotate_x/1) end)
+    times(n, fn -> emit("m44.make_rotate_y", [ang()], &Mat44.make_rotate_y/1) end)
+    times(n, fn -> emit("m44.make_rotate_z", [ang()], &Mat44.make_rotate_z/1) end)
+    times(n, fn -> emit("m44.make_translate", [r(), r(), r()], &Mat44.make_translate/3) end)
+    times(n, fn -> emit("m44.make_scale1", [r()], &Mat44.make_scale/1) end)
+    times(n, fn -> emit("m44.make_scale4", [r(), r(), r(), r()], &Mat44.make_scale/4) end)
+    times(n, fn -> emit("m44.make_reflect", [v3(), r()], &Mat44.make_reflect/2) end)
+    times(n, fn -> emit("m44.make_shear_x", [r(), r()], &Mat44.make_shear_x/2) end)
+    times(n, fn -> emit("m44.make_shear_y", [r(), r()], &Mat44.make_shear_y/2) end)
+    times(n, fn -> emit("m44.make_shear_z", [r(), r()], &Mat44.make_shear_z/2) end)
+    times(n, fn -> emit("m44.orient", [v3(), v3(), v3()], &Mat44.orient/3) end)
+    times(n, fn -> emit("m44.look_at", [v3(), v3(), v3()], &Mat44.look_at/3) end)
+    times(n, fn -> emit("m44.make_billboard", [v3(), v3(), v3()], &Mat44.make_billboard/3) end)
+
+    times(n, fn ->
+      emit("m44.make_billboard_axis", [v3(), v3(), v3()], &Mat44.make_billboard_axis/3)
+    end)
+
+    times(n, fn ->
+      emit(
+        "m44.perspective",
+        [:rand.uniform() * 2.5 + 0.3, rp(), rp() * 0.2, rp() * 10 + 5],
+        &Mat44.perspective/4
+      )
+    end)
+
+    times(n, fn ->
+      {a, b} = {r(), r()}
+      {c, d} = {r(), r()}
+      {e, f} = {r(), r()}
+      emit("m44.ortho", [a, b, c, d, e, f], &Mat44.ortho/6)
+    end)
+
+    # ---------------- Quatern
+    emit("q.identity", [], &Quatern.identity/0)
+    emit("q.zero", [], &Quatern.zero/0)
+    times(n, fn -> emit("q.create", [r(), r(), r(), r()], &Quatern.create/4) end)
+    times(n, fn -> emit("q.from_list", [[r(), r(), r(), r()]], &Quatern.from_list/1) end)
+    times(n, fn -> emit("q.add", [v4(), v4()], &Quatern.add/2) end)
+    times(n, fn -> emit("q.subtract", [v4(), v4()], &Quatern.subtract/2) end)
+    times(n, fn -> emit("q.multiply", [v4(), v4()], &Quatern.multiply/2) end)
+    times(n, fn -> emit("q.scale", [v4(), r()], &Quatern.scale/2) end)
+    times(n, fn -> emit("q.dot", [v4(), v4()], &Quatern.dot/2) end)
+    times(n, fn -> emit("q.norm", [v4()], &Quatern.norm/1) end)
+    times(n, fn -> emit("q.normalize", [v4()], &Quatern.normalize/1) end)
+    times(n, fn -> emit("q.normalize_strict", [v4()], &Quatern.normalize_strict/1) end)
+    times(n, fn -> emit("q.conjugate", [v4()], &Quatern.conjugate/1) end)
+    times(n, fn -> emit("q.inverse", [v4()], &Quatern.inverse/1) end)
+    times(n, fn -> emit("q.equal", [v4(), v4()], &Quatern.equal/2) end)
+
+    times(n, fn ->
+      a = v4()
+      emit("q.equal", [a, a], &Quatern.equal/2)
+    end)
+
+    times(n, fn -> emit("q.equal_elements", [v4(), v4(), rp()], &Quatern.equal_elements/3) end)
+    times(n, fn -> emit("q.from_axis_angle", [ang(), unit3()], &Quatern.from_axis_angle/2) end)
+    times(n, fn -> emit("q.to_rotation_matrix_33", [uq()], &Quatern.to_rotation_matrix_33/1) end)
+    times(n, fn -> emit("q.to_rotation_matrix_44", [uq()], &Quatern.to_rotation_matrix_44/1) end)
+    times(n, fn -> emit("q.transform_vector", [uq(), v3()], &Quatern.transform_vector/2) end)
+    # from_rotation_matrix: feed a matrix from a known quaternion; also feed negative-trace cases
+    times(n, fn ->
+      q = uq()
+
+      emit(
+        "q.from_rotation_matrix",
+        [Quatern.to_rotation_matrix_33(q)],
+        &Quatern.from_rotation_matrix/1
+      )
+    end)
+
+    for ax <- [{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}],
+        th <- [3.0, -3.0, :math.pi()] do
+      q = Quatern.from_axis_angle(th, ax)
+
+      emit(
+        "q.from_rotation_matrix",
+        [Quatern.to_rotation_matrix_33(q)],
+        &Quatern.from_rotation_matrix/1
+      )
+    end
+
+    times(n, fn -> emit("q.slerp", [uq(), uq(), :rand.uniform()], &Quatern.slerp/3) end)
+
+    times(3, fn ->
+      q = uq()
+      emit("q.slerp", [q, Quatern.scale(q, -1.0), :rand.uniform()], &Quatern.slerp/3)
+    end)
+
+    times(3, fn ->
+      q = uq()
+      emit("q.slerp", [q, q, :rand.uniform()], &Quatern.slerp/3)
+    end)
+
+    times(n, fn ->
+      emit("q.integrate", [uq(), v3(), :rand.uniform() * 0.5], &Quatern.integrate/3)
+    end)
+
+    times(3, fn -> emit("q.integrate", [uq(), v3(), 1.0e-9], &Quatern.integrate/3) end)
+    times(n, fn -> emit("q.random", [], &Quatern.random/0) end)
+    # euler extraction against single-axis rotations
+    times(n, fn ->
+      th = ang()
+
+      emit("q.get_roll_z", [th], fn t ->
+        Quatern.get_roll(Quatern.from_axis_angle(t, {0.0, 0.0, 1.0}))
+      end)
+    end)
+
+    times(n, fn ->
+      th = ang()
+
+      emit("q.get_pitch_x", [th], fn t ->
+        Quatern.get_pitch(Quatern.from_axis_angle(t, {1.0, 0.0, 0.0}))
+      end)
+    end)
+
+    times(n, fn ->
+      th = ang()
+
+      emit("q.get_yaw_y", [th], fn t ->
+        Quatern.get_yaw(Quatern.from_axis_angle(t, {0.0, 1.0, 0.0}))
+      end)
+    end)
+
+    times(n, fn ->
+      q = uq()
+
+      emit("q.euler", [q], fn q ->
+        {Quatern.get_roll(q), Quatern.get_pitch(q), Quatern.get_yaw(q)}
+      end)
+    end)
+
+    # --- documented error cases and edge cases
+    times(2, fn -> emit("q.equal", [uq(), uq()], &Quatern.equal/2) end)
+
+    times(3, fn ->
+      q = uq()
+      emit("q.equal", [q, Quatern.scale(q, -1.0)], &Quatern.equal/2)
+    end)
+
+    times(3, fn ->
+      q = uq()
+      emit("q.equal", [q, q], &Quatern.equal/2)
+    end)
+
+    emit("err.v2.project_zero", [v2(), {0.0, 0.0}], &Vec2.project/2)
+    emit("err.v4.project_zero", [v4(), {0.0, 0.0, 0.0, 0.0}], &Vec4.project/2)
+    emit("err.m22.make_reflect_zero", [{0.0, 0.0}], &Mat22.make_reflect/1)
+    emit("err.m33.make_reflect_2d_zero", [{0.0, 0.0}, 1.0], &Mat33.make_reflect_2d/2)
+    emit("err.m33.make_reflect_3d_zero", [{0.0, 0.0, 0.0}], &Mat33.make_reflect_3d/1)
+    emit("err.m44.make_reflect_zero", [{0.0, 0.0, 0.0}, 1.0], &Mat44.make_reflect/2)
+
+    emit(
+      "err.m33.orthonormalize_zero_row",
+      [{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0}],
+      &Mat33.orthonormalize/1
+    )
+
+    emit(
+      "err.m33.orthonormalize_dependent",
+      [{1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 1.0, 0.0}],
+      &Mat33.orthonormalize/1
+    )
+
+    emit(
+      "err.m33.orthonormalize_coplanar",
+      [{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0}],
+      &Mat33.orthonormalize/1
+    )
+
+    emit("err.m44.orient_zero_fwd", [v3(), {0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}], &Mat44.orient/3)
+    emit("err.m44.orient_parallel_up", [v3(), {0.0, 2.0, 0.0}, {0.0, -1.0, 0.0}], &Mat44.orient/3)
+
+    emit(
+      "err.m44.look_at_coincident",
+      [{1.0, 2.0, 3.0}, {1.0, 2.0, 3.0}, {0.0, 1.0, 0.0}],
+      &Mat44.look_at/3
+    )
+
+    emit(
+      "err.m44.look_at_zero_up",
+      [{1.0, 2.0, 3.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}],
+      &Mat44.look_at/3
+    )
+
+    emit(
+      "err.m44.billboard_coincident",
+      [{1.0, 2.0, 3.0}, {1.0, 2.0, 3.0}, {0.0, 1.0, 0.0}],
+      &Mat44.make_billboard/3
+    )
+
+    emit(
+      "err.m44.billboard_axis_on_axis",
+      [{0.0, 0.0, 0.0}, {0.0, 5.0, 0.0}, {0.0, 1.0, 0.0}],
+      &Mat44.make_billboard_axis/3
+    )
+
+    emit(
+      "err.m44.billboard_axis_zero_axis",
+      [{0.0, 0.0, 0.0}, {0.0, 5.0, 0.0}, {0.0, 0.0, 0.0}],
+      &Mat44.make_billboard_axis/3
+    )
+
+    emit("err.m44.perspective_fov0", [0.0, 1.0, 0.1, 10.0], &Mat44.perspective/4)
+    emit("err.m44.perspective_fovpi", [:math.pi(), 1.0, 0.1, 10.0], &Mat44.perspective/4)
+    emit("err.m44.perspective_aspect0", [1.0, 0.0, 0.1, 10.0], &Mat44.perspective/4)
+    emit("err.m44.perspective_near0", [1.0, 1.0, 0.0, 10.0], &Mat44.perspective/4)
+    emit("err.m44.perspective_near_ge_far", [1.0, 1.0, 10.0, 10.0], &Mat44.perspective/4)
+    emit("err.m44.ortho_x_equal", [1.0, 1.0, 0.0, 1.0, 0.0, 1.0], &Mat44.ortho/6)
+    emit("err.m44.ortho_y_equal", [0.0, 1.0, 1.0, 1.0, 0.0, 1.0], &Mat44.ortho/6)
+    emit("err.m44.ortho_z_equal", [0.0, 1.0, 0.0, 1.0, 2.0, 2.0], &Mat44.ortho/6)
+    emit("err.m22.inverse_singular", [{1.0, 2.0, 2.0, 4.0}], &Mat22.inverse/1)
+
+    emit(
+      "err.m33.inverse_singular",
+      [{1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 1.0}],
+      &Mat33.inverse/1
+    )
+
+    emit("err.m44.inverse_singular", [Mat44.zero()], &Mat44.inverse/1)
+    emit("err.q.normalize_strict_zero", [{0.0, 0.0, 0.0, 0.0}], &Quatern.normalize_strict/1)
+    emit("edge.q.normalize_zero", [{0.0, 0.0, 0.0, 0.0}], &Quatern.normalize/1)
+    emit("edge.q.inverse_zero", [{0.0, 0.0, 0.0, 0.0}], &Quatern.inverse/1)
+    emit("edge.q.integrate_zero_q", [{0.0, 0.0, 0.0, 0.0}, v3(), 0.1], &Quatern.integrate/3)
+    emit("edge.q.integrate_zero_omega", [uq(), {0.0, 0.0, 0.0}, 0.1], &Quatern.integrate/3)
+    emit("edge.q.integrate_zero_dt", [uq(), v3(), 0.0], &Quatern.integrate/3)
+    emit("edge.v3.normalize_zero", [{0.0, 0.0, 0.0}], &Vec3.normalize/1)
+    emit("edge.v2.normalize_zero", [{0.0, 0.0}], &Vec2.normalize/1)
+    emit("edge.v4.normalize_zero", [{0.0, 0.0, 0.0, 0.0}], &Vec4.normalize/1)
+    emit("edge.v2.near_exact", [{0.0, 0.0}, {3.0, 4.0}, 5.0], &Vec2.near/3)
+    emit("edge.v3.create_ints", [[1, 2, 3]], &Vec3.create/1)
+    emit("edge.q.from_list_ints", [[1, 0, 0, 0]], &Quatern.from_list/1)
+    # doc examples
+    emit("edge.m44.orient_doc", [{2.0, 3.0, 4.0}, {0.0, 0.0, -1.0}, {0.0, 1.0, 0.0}], fn p,
+                                                                                         f,
+                                                                                         u ->
+      Mat44.transform_point(Mat44.orient(p, f, u), {1.0, 2.0, -3.0})
+    end)
+
+    emit("edge.m44.look_at_doc", [{0.0, 0.0, 5.0}, {0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}], fn e,
+                                                                                         c,
+                                                                                         u ->
+      Mat44.transform_point(Mat44.look_at(e, c, u), {2.0, 3.0, 0.0})
+    end)
+
+    emit("edge.m44.perspective_doc", [], fn ->
+      {x, y, z, w} =
+        Mat44.apply_left({0.0, 0.0, -1.0, 1.0}, Mat44.perspective(:math.pi() / 2, 1.0, 1.0, 3.0))
+
+      {x / w, y / w, z / w}
+    end)
+
+    emit("edge.m44.perspective_far", [], fn ->
+      {x, y, z, w} =
+        Mat44.apply_left({0.0, 0.0, -3.0, 1.0}, Mat44.perspective(:math.pi() / 2, 1.0, 1.0, 3.0))
+
+      {x / w, y / w, z / w}
+    end)
+
+    emit("edge.m44.ortho_doc", [], fn ->
+      Mat44.transform_point(Mat44.ortho(-2.0, 6.0, -4.0, 4.0, 1.0, 3.0), {2.0, 0.0, -2.0})
+    end)
+
+    emit("edge.m44.rotate_y_doc", [], fn ->
+      Mat44.transform_vector(Mat44.make_rotate_y(:math.pi() / 2), {1.0, 0.0, 0.0})
+    end)
+
+    emit(
+      "edge.m33.orthonormalize_doc",
+      [{2.0, 0.0, 0.0, 1.0, 3.0, 0.0, 4.0, 5.0, 6.0}],
+      &Mat33.orthonormalize/1
+    )
+
+    finish()
+    IO.puts("wrote #{Process.get(:k)} cases")
+  end
+end
+
+Main.run()
diff --git a/bench/octave/run.sh b/bench/octave/run.sh
new file mode 100755
index 0000000..171bfbb
--- /dev/null
+++ b/bench/octave/run.sh
@@ -0,0 +1,7 @@
+#!/usr/bin/env bash
+# Cross-check every public graphmath function against Octave reference math.
+set -euo pipefail
+cd "$(dirname "$0")/../.."
+out="${1:-$(mktemp -t graphmath-cases-XXXX.m)}"
+mix run bench/octave/gen_cases.exs "$out"
+octave-cli --no-gui bench/octave/check.m "$out"