Interpolating Rotations - Slerp and Nlerp

The three ways to blend between two orientations, and the fact that two of them follow the identical path - so the choice is about speed, not route. How big the difference actually is, in degrees, so “cheap is fine” becomes a rule with a number in it. And where precision goes when you convert between quaternions, matrices and Euler angles, which is not where you might guess.

Say you want a character to turn from one orientation to another over half a second. You have the two quaternions. You need the ones in between.

All three methods below are built on linear interpolation, universally shortened to lerp, so it is worth writing down before anything is built on top of it:

Read it as “start at $a$, then go a fraction $t$ of the way towards $b$”. At $t = 0$ it is exactly $a$, at $t = 1$ it is exactly $b$, and at $t = 0.5$ it is the midpoint. Equivalently it is the weighted average $(1-t)a + tb$, which is the form that generalises.

Two things worth knowing now. It works on anything you can subtract and scale - numbers, positions, colours, and the four components of a quaternion. And nothing stops $t$ leaving $[0, 1]$: at $t = 2$ you get a point twice as far as $b$, which is extrapolation rather than an error, and is occasionally what you want and more often a bug.

lerpQuat in the code panel further down is exactly the formula above, applied four times - once per component. That is the entire function.

The names are acronyms: nlerp is Normalized Linear Interpolation (lerp, then normalize), and slerp is Spherical Linear Interpolation (linear along the surface of the sphere rather than through its interior).

The usual telling is that these take three different routes. That is wrong, and the correction matters because it tells you what you are actually choosing between.

nlerp and slerp trace the same arc. Normalizing a straight line between two points on a sphere lands you on the great circle through them - the chord and the arc pass through the same directions, just at different rates. The build check proves this by showing every sample from both methods lies in the two-dimensional plane spanned by the endpoints, to $10^{-16}$.

So slerp does not fix the path. It fixes the speed.

Below, one 150-degree turn, done twice. The dots mark every tenth of the way through the blend.

/**
 * One arc, two rates: evenly spaced slerp ticks against nlerp ticks bunched at the ends.
 */
import * as THREE from "three";
import {
  applyMat4,
  multiplyMat4,
  point,
  translation4,
  type Mat4,
} from "../matrices.ts";
import {
  angleBetweenQuats,
  fromAxisAngle,
  nlerpQuat,
  quatToMat4,
  rotateVector,
  slerpQuat,
  type Quat,
} from "../quaternions.ts";
import { makeCanvas, addSlider, addReadout, addBoxWire } from "./ui.ts";
import type { MountFn } from "./runner.ts";

const Y = { x: 0, y: 1, z: 0 };
const TOTAL = 150;
const FROM: Quat = fromAxisAngle(Y, 0)!;
const TO: Quat = fromAxisAngle(Y, TOTAL)!;
const NOSE = { x: 0, y: 0, z: -1 };
const REACH = 1.3;
const TICKS = 10;

type Blend = (from: Quat, to: Quat, t: number) => Quat | null;

const mount: MountFn = (el) => {
  const { renderer, width, height, background } = makeCanvas(el, 300);

  const scene = new THREE.Scene();
  scene.background = background;
  scene.add(new THREE.GridHelper(10, 10, 0x30363d, 0x21262d));

  const camera = new THREE.PerspectiveCamera(40, width / height, 0.1, 100);
  camera.position.set(0.2, 5.2, 4.6);
  camera.lookAt(0, 0, 0);

  /** One object, its arc, and a dot at every tenth of the way through the blend. */
  function rig(centre: number, color: number, blend: Blend) {
    const box = addBoxWire(scene, color);
    const arc = new THREE.Line(
      new THREE.BufferGeometry(),
      new THREE.LineBasicMaterial({ color: 0x545d68 }),
    );
    scene.add(arc);

    const dots: THREE.Mesh[] = [];
    for (let i = 0; i <= TICKS; i += 1) {
      const dot = new THREE.Mesh(
        new THREE.SphereGeometry(0.065, 10, 8),
        new THREE.MeshBasicMaterial({ color }),
      );
      scene.add(dot);
      dots.push(dot);
    }

    const tipAt = (t: number) => {
      const f = rotateVector(blend(FROM, TO, t)!, NOSE);
      return new THREE.Vector3(centre + f.x * REACH, f.y * REACH, f.z * REACH);
    };

    // The arc never changes, and neither do the ticks. Draw them once.
    const path: THREE.Vector3[] = [];
    for (let i = 0; i <= 120; i += 1) path.push(tipAt(i / 120));
    arc.geometry.setFromPoints(path);
    dots.forEach((d, i) => d.position.copy(tipAt(i / TICKS)));

    return (t: number) => {
      const q = blend(FROM, TO, t)!;
      const m: Mat4 = multiplyMat4(translation4(centre, 0, 0), quatToMat4(q));
      box((c) => {
        const p = applyMat4(m, point(c[0] * 0.55, c[1] * 0.55, c[2] * 0.55));
        return [p.x, p.y, p.z];
      });
      return angleBetweenQuats(FROM, q);
    };
  }

  const nlerpRig = rig(-1.7, 0xd2a8ff, nlerpQuat);
  const slerpRig = rig(1.7, 0x39d3c3, slerpQuat);

  const show = addReadout(el);
  const t = addSlider(
    el,
    "blend from start to end",
    0,
    1,
    0.23,
    draw,
    "",
    0.01,
  );

  function draw() {
    const turnedN = nlerpRig(t());
    const turnedS = slerpRig(t());
    const behind = turnedS - turnedN;
    show(
      `nlerp ${turnedN.toFixed(1)}\u00B0, slerp ${turnedS.toFixed(1)}\u00B0  \u00B7  ` +
        `nlerp is ${Math.abs(behind).toFixed(1)}\u00B0 ` +
        `${behind >= 0 ? "behind" : "ahead"}  \u00B7  slerp is always t \u00D7 ${TOTAL}\u00B0`,
    );
    renderer.render(scene, camera);
  }

  draw();

  return () => renderer.dispose();
};

export default mount;

Look at the dot spacing. Teal slerp dots are evenly spread along the arc, because equal steps in $t$ mean equal steps in angle. Purple nlerp dots bunch up at both ends and spread out in the middle - it starts slow, overtakes, then slows again.

Drag the slider and the readout gives both. They agree exactly at $t = 0$, $t = 0.5$ and $t = 1$, and disagree everywhere else.

The formula is what you would expect once you know the arc is shared:

Two weights that sum along the arc rather than across the chord.

$\Omega$ is capital omega, and it is nothing more than a name for one angle - the traditional letter for it in this formula, not a new operation.

Which angle? The one between the two quaternions, obtained from their dot product exactly the way Part 1 got the angle between two 3D directions. Same idea, two more components: treat each quaternion as a unit vector in four dimensions, dot them together, take the arccosine.

One catch, and it comes straight from Section 3.2’s half-angle. Because a quaternion stores $\theta/2$, the angle between two quaternions is half the rotation angle you would see on screen:

So the 150-degree turn in the scene above has $\Omega = 75°$. That factor of two is also why angleBetweenQuats multiplies by 2 before reporting anything - it converts from quaternion space back to the angle a person would measure.

Sanity check the formula with it. At $t = 0$ the weights become $\sin\Omega / \sin\Omega = 1$ and $\sin 0 / \sin \Omega = 0$, so you get $q_1$ exactly. At $t = 1$ they swap and you get $q_2$. In between, $\sin\Omega$ in the denominator is what scales the pair so the result lands on the arc rather than short of it - which is precisely the job nlerp hands to a renormalize instead.

The check asserts that slerp’s angle from the start is exactly $t \times \text{total}$, to within $10^{-6}$ degrees, at every step of four different turn sizes.

Here is the useful part, because “always use slerp” is advice that ignores the cost.

nlerp’s error depends entirely on how big the turn is. Measured as the worst gap between where nlerp is and where a constant rate would put it:

Those numbers are asserted by the build check, including that the error grows monotonically with the turn size - which is what makes the rule below trustworthy rather than a guess.

So: use nlerp for small blends, slerp for large ones. And note which case dominates in practice. Animation blending, camera smoothing and per-frame turn-towards-target all interpolate tiny angles many times per second, where nlerp is both cheaper and indistinguishable. Slerp earns its keep on deliberate large moves - a scripted camera swinging around, a cutscene, a long reorientation.

Most engines slerp everywhere anyway, because a couple of trig calls per bone per frame stopped mattering some time ago. The reason to know the difference is that when it does matter - a thousand-bone crowd, a physics substep - you can tell which one you can afford to drop.

The third method is the one people write first: blend the four numbers and use the result.

The problem is that a straight line between two points on a sphere passes inside the sphere, so the result is not a unit quaternion. Over a 150-degree turn its length dips to $0.7934$ at the midpoint.

And an un-normalized quaternion does not just rotate. The sandwich $q v q^{*}$ scales by $|q|^2$, so at that midpoint everything comes out at $0.63$ of its proper size.

/**
 * What skipping the renormalize costs: the object shrinks, because the sandwich scales by |q|.
 */
import * as THREE from "three";
import {
  fromAxisAngle,
  lerpQuat,
  normalizeQuat,
  quatLength,
  rotateVector,
  type Quat,
} from "../quaternions.ts";
import {
  makeCanvas,
  addSlider,
  addCheckbox,
  addReadout,
  addBoxWire,
} from "./ui.ts";
import type { MountFn } from "./runner.ts";

const Y = { x: 0, y: 1, z: 0 };
const FROM: Quat = fromAxisAngle(Y, 0)!;
const TO: Quat = fromAxisAngle(Y, 150)!;

const mount: MountFn = (el) => {
  const { renderer, width, height, background } = makeCanvas(el, 290);

  const scene = new THREE.Scene();
  scene.background = background;
  scene.add(new THREE.GridHelper(8, 8, 0x30363d, 0x21262d));

  const camera = new THREE.PerspectiveCamera(38, width / height, 0.1, 100);
  camera.position.set(2.6, 2.4, 3.6);
  camera.lookAt(0, 0, 0);

  // Dashed grey is the size the object is supposed to be, so the shrink has a reference.
  const reference = addBoxWire(scene, 0x7d8590, { dashed: true });
  const actual = addBoxWire(scene, 0xf0883e);

  const show = addReadout(el);
  const t = addSlider(el, "blend from start to end", 0, 1, 0.5, draw, "", 0.01);
  const fix = addCheckbox(el, "normalize before using it", false, draw);

  function draw() {
    const raw = lerpQuat(FROM, TO, t());
    const q = fix() ? normalizeQuat(raw)! : raw;
    const len = quatLength(q);

    // Rotating through the sandwich, which scales by |q| squared when |q| is not 1.
    actual((c) => {
      const v = rotateVector(q, { x: c[0], y: c[1], z: c[2] });
      return [v.x, v.y, v.z];
    });
    // The same corners under a quaternion that has been normalized, whatever the checkbox says.
    const unit = normalizeQuat(raw)!;
    reference((c) => {
      const v = rotateVector(unit, { x: c[0], y: c[1], z: c[2] });
      return [v.x, v.y, v.z];
    });

    show(
      `|q| = ${len.toFixed(3)}, so the object is scaled by ` +
        `${(len * len).toFixed(3)}` +
        (fix() ? "  \u00B7  fixed" : "  \u00B7  shrinking"),
    );
    renderer.render(scene, camera);
  }

  draw();

  return () => renderer.dispose();
};

export default mount;

The dashed grey box is the size the object should be. Drag the blend and watch the orange one shrink away from it, then tick the checkbox to put the normalize back.

That is the whole difference between lerp and nlerp - one square root. It is also why this is rarely a bug you meet in the wild: the shrinking is so obvious that it gets caught immediately. The subtle bugs in rotation code are the double cover from Section 3.2 and the speed difference above, not this one.

/**
 * Four numbers that hold an orientation, with no poles and nothing to lock.
 *
 * The whole construction hangs off one decision: store **half** the angle. That is what makes
 * multiplication compose rotations, makes the conjugate an inverse, and makes the four numbers
 * cover every orientation exactly twice - the double cover, which is the one place quaternions
 * ask something of you in return.
 */
import { direction, point, type Mat4, type Vec3 } from "./matrices.ts";

/** `x, y, z` are the axis scaled by `sin(angle/2)`. `w` is `cos(angle/2)`. */
export type Quat = { x: number; y: number; z: number; w: number };

/** No rotation: a zero-length axis part and a full-size `w`. */
export const IDENTITY_QUAT: Quat = { x: 0, y: 0, z: 0, w: 1 };

const DEG = Math.PI / 180;

export function quatLength(q: Quat): number {
  return Math.hypot(q.x, q.y, q.z, q.w);
}

/** Scale back to unit length. Returns `null` for the one quaternion that has no direction. */
export function normalizeQuat(q: Quat): Quat | null {
  const len = quatLength(q);
  if (len < 1e-12) return null;
  return { x: q.x / len, y: q.y / len, z: q.z / len, w: q.w / len };
}

export function negateQuat(q: Quat): Quat {
  return { x: -q.x, y: -q.y, z: -q.z, w: -q.w };
}

export function dotQuat(a: Quat, b: Quat): number {
  return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
}

/**
 * Build a rotation from an axis and an angle.
 *
 * Note the halves. Turning by 90 degrees stores `cos(45)` and `sin(45)`, and if you ever see
 * a quaternion whose `w` is `cos` of the angle you meant rather than half of it, this is the
 * line to look at.
 */
export function fromAxisAngle(axis: Vec3, degrees: number): Quat | null {
  const len = Math.hypot(axis.x, axis.y, axis.z);
  if (len < 1e-12) return null;
  const half = degrees * DEG * 0.5;
  const s = Math.sin(half) / len;
  return {
    x: axis.x * s,
    y: axis.y * s,
    z: axis.z * s,
    w: Math.cos(half),
  };
}

/** Read the axis and angle back out. The angle comes back doubled, for the same reason. */
export function toAxisAngle(q: Quat): { axis: Vec3; degrees: number } {
  const n = normalizeQuat(q) ?? IDENTITY_QUAT;
  const w = Math.min(1, Math.max(-1, n.w));
  const sinHalf = Math.sqrt(Math.max(0, 1 - w * w));
  // Below this the rotation is too small to have a meaningful axis, so name one.
  if (sinHalf < 1e-9) return { axis: { x: 0, y: 1, z: 0 }, degrees: 0 };
  return {
    axis: { x: n.x / sinHalf, y: n.y / sinHalf, z: n.z / sinHalf },
    degrees: (2 * Math.acos(w)) / DEG,
  };
}

/**
 * Compose two rotations: apply `second` after `first`.
 *
 * Same right-to-left reading as `multiplyMat4`, and for the same reason - so the two can be
 * swapped for each other without any code changing its meaning. Quaternion multiplication does
 * not commute either, which it had better not, since rotations do not.
 */
export function multiplyQuat(second: Quat, first: Quat): Quat {
  const a = second;
  const b = first;
  return {
    w: a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
    x: a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
    y: a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
    z: a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
  };
}

/**
 * Flip the axis part. For a unit quaternion this **is** the inverse - same axis, opposite
 * angle - which is a division replaced by three sign changes.
 */
export function conjugate(q: Quat): Quat {
  return { x: -q.x, y: -q.y, z: -q.z, w: q.w };
}

/**
 * Rotate a vector, by sandwiching it between the quaternion and its conjugate.
 *
 * Written out longhand because the sandwich is the thing worth seeing. Production code uses an
 * algebraically identical shortcut with fewer multiplies.
 *
 * If `q` is **not** unit length the sandwich also scales the vector, by exactly `|q|` squared.
 * That is the mechanism behind the classic "I forgot to normalize and my model shrank".
 */
export function rotateVector(q: Quat, v: Vec3): Vec3 {
  const asQuat: Quat = { x: v.x, y: v.y, z: v.z, w: 0 };
  const out = multiplyQuat(multiplyQuat(q, asQuat), conjugate(q));
  return { x: out.x, y: out.y, z: out.z };
}

/** The same rotation as a 4x4, so it drops straight into Part 2's machinery. */
export function quatToMat4(q: Quat): Mat4 {
  const n = normalizeQuat(q) ?? IDENTITY_QUAT;
  const { x, y, z, w } = n;
  return {
    i: direction(
      1 - 2 * (y * y + z * z),
      2 * (x * y + z * w),
      2 * (x * z - y * w),
    ),
    j: direction(
      2 * (x * y - z * w),
      1 - 2 * (x * x + z * z),
      2 * (y * z + x * w),
    ),
    k: direction(
      2 * (x * z + y * w),
      2 * (y * z - x * w),
      1 - 2 * (x * x + y * y),
    ),
    t: point(0, 0, 0),
  };
}

/** How far apart two orientations are, in degrees. The absolute value handles double cover. */
export function angleBetweenQuats(a: Quat, b: Quat): number {
  const na = normalizeQuat(a);
  const nb = normalizeQuat(b);
  if (na === null || nb === null) return 0;
  const d = Math.min(1, Math.abs(dotQuat(na, nb)));
  return (2 * Math.acos(d)) / DEG;
}

/**
 * The double cover fix, and the single most load-bearing `if` in rotation code.
 *
 * Every orientation has **two** quaternions: `q` and `-q`. They produce the identical matrix,
 * so nothing you can observe about the object tells them apart. But interpolation is arithmetic
 * on the four numbers, and those differ - so blending towards the wrong one of the pair sends
 * the object the long way round, up to 360 degrees when 0 would have done.
 *
 * A negative dot product is exactly the signal that you have the far one. Flip it.
 */
export function shortWayFrom(from: Quat, to: Quat): Quat {
  return dotQuat(from, to) < 0 ? negateQuat(to) : to;
}

/**
 * Blend two orientations the cheap way: interpolate all four numbers, then renormalize.
 *
 * Deliberately does **not** call `shortWayFrom` itself, so a demo can show what happens
 * without it. Section 3.3 compares this against slerp; the short-way fix is needed either way
 * and is not what distinguishes them.
 */
export function nlerpQuat(from: Quat, to: Quat, t: number): Quat | null {
  return normalizeQuat({
    x: from.x + (to.x - from.x) * t,
    y: from.y + (to.y - from.y) * t,
    z: from.z + (to.z - from.z) * t,
    w: from.w + (to.w - from.w) * t,
  });
}

/**
 * Straight linear interpolation of the four numbers, with **no** renormalizing.
 *
 * Included because it is what people write first, and because the result is not a unit
 * quaternion - so anything using it either renormalizes later or quietly scales the object.
 */
export function lerpQuat(from: Quat, to: Quat, t: number): Quat {
  return {
    x: from.x + (to.x - from.x) * t,
    y: from.y + (to.y - from.y) * t,
    z: from.z + (to.z - from.z) * t,
    w: from.w + (to.w - from.w) * t,
  };
}

/**
 * Spherical interpolation: travel the arc between two orientations at a **constant rate**.
 *
 * `nlerpQuat` already follows the right arc - normalizing a straight line between two points on
 * a sphere lands you on the great circle through them. What it gets wrong is the *speed*, and
 * that is the only thing slerp fixes.
 *
 * Like `nlerpQuat`, this does not apply the double-cover flip. Call `shortWayFrom` first, or
 * you will slerp smoothly along the 320-degree route from Section 3.2.
 */
export function slerpQuat(from: Quat, to: Quat, t: number): Quat | null {
  const a = normalizeQuat(from);
  const b = normalizeQuat(to);
  if (a === null || b === null) return null;

  const cos = Math.min(1, Math.max(-1, dotQuat(a, b)));

  /* When the two are nearly the same orientation, `sin(omega)` heads for zero and the weights
     below turn into 0/0. A straight blend is indistinguishable from the arc at that scale, so
     hand off rather than divide. This guard is not optional. */
  if (Math.abs(cos) > 0.9995) return nlerpQuat(a, b, t);

  const omega = Math.acos(cos);
  const sin = Math.sin(omega);
  const wa = Math.sin((1 - t) * omega) / sin;
  const wb = Math.sin(t * omega) / sin;

  return {
    x: a.x * wa + b.x * wb,
    y: a.y * wa + b.y * wb,
    z: a.z * wa + b.z * wb,
    w: a.w * wa + b.w * wb,
  };
}

The other half of this Section. You now have three representations and six conversions between them, and it is worth knowing which ones lose something.

/** Where precision actually goes when you convert between the three representations. */
import { YAW_PITCH_ROLL, fromEuler, toEulerYXZ, type Euler } from "../euler.ts";
import {
  fromAxisAngle,
  multiplyQuat,
  quatLength,
  quatToMat4,
} from "../quaternions.ts";
import { basisOf } from "../spaces.ts";
import { rowsOf, type Mat4 } from "../matrices.ts";
import type { Demo } from "./runner.ts";

const worstEntry = (a: Mat4, b: Mat4) => {
  const x = rowsOf(a).flat();
  const y = rowsOf(b).flat();
  return Math.max(...x.map((n, i) => Math.abs(n - y[i])));
};

/** Angles out and back, and how far the three numbers moved. */
function angleTrip(pitch: number) {
  const start: Euler = { x: pitch, y: 40, z: 25 };
  const m = fromEuler(start, YAW_PITCH_ROLL);
  const back = toEulerYXZ(m);
  return {
    angles: Math.max(
      Math.abs(back.x - start.x),
      Math.abs(back.y - start.y),
      Math.abs(back.z - start.z),
    ),
    matrix: worstEntry(m, fromEuler(back, YAW_PITCH_ROLL)),
  };
}

const demo: Demo = (log) => {
  const mid = angleTrip(30);
  const near = angleTrip(89.999);
  const pole = angleTrip(90);

  log(
    "angles to matrix and back, pitch 30",
    `${mid.angles.toExponential(1)}\u00B0`,
  );
  log("the same at pitch 89.999", `${near.angles.toExponential(1)}\u00B0`);
  log(
    "the same at pitch 90",
    `${pole.angles.toFixed(1)}\u00B0`,
    "roll folded into yaw",
  );
  log(
    "but the orientation at pitch 90",
    pole.matrix.toExponential(1),
    "intact - only the three numbers were lost",
  );

  // Drift under repeated composition, the other place precision goes.
  const step = fromAxisAngle({ x: 0.2, y: 0.9, z: 0.35 }, 0.7)!;
  let q = fromAxisAngle({ x: 0, y: 1, z: 0 }, 0)!;
  let m = quatToMat4(q);
  const stepM = quatToMat4(step);
  const mul = (A: Mat4, B: Mat4): Mat4 => {
    const ra = rowsOf(A);
    const rb = rowsOf(B);
    const g = (r: number, c: number) =>
      ra[r][0] * rb[0][c] +
      ra[r][1] * rb[1][c] +
      ra[r][2] * rb[2][c] +
      ra[r][3] * rb[3][c];
    return {
      i: { x: g(0, 0), y: g(1, 0), z: g(2, 0), w: g(3, 0) },
      j: { x: g(0, 1), y: g(1, 1), z: g(2, 1), w: g(3, 1) },
      k: { x: g(0, 2), y: g(1, 2), z: g(2, 2), w: g(3, 2) },
      t: { x: g(0, 3), y: g(1, 3), z: g(2, 3), w: g(3, 3) },
    };
  };
  for (let i = 0; i < 10000; i += 1) {
    q = multiplyQuat(step, q);
    m = mul(stepM, m);
  }

  const b = basisOf(m);
  const dot = (p: typeof b.i, r: typeof b.i) =>
    p.x * r.x + p.y * r.y + p.z * r.z;

  log(
    "10000 quaternion products, length off by",
    Math.abs(quatLength(q) - 1).toExponential(1),
    "one number to fix",
  );
  log(
    "10000 matrix products, axes off square by",
    Math.max(
      Math.abs(dot(b.i, b.j)),
      Math.abs(dot(b.j, b.k)),
      Math.abs(dot(b.i, b.k)),
    ).toExponential(1),
    "nine numbers to fix",
  );
};

export default demo;

Read the first four rows together, because they tell one story.

Angles out to a matrix and back are accurate to about $10^{-15}$ degrees at pitch 30 - which is to say exact. At pitch 89.999 that has degraded to $10^{-11}$, and at pitch 90 the angles come back 25 degrees different, because toEulerYXZ hits the degenerate case, sets roll to zero and hands the whole rotation to yaw.

And yet the fourth row shows the matrix rebuilt from those “wrong” angles is identical to the original to $10^{-16}$.

That is the precise statement of what Euler conversion costs: the orientation survives, the three numbers do not. Near the pole there are many triples describing the same orientation, and the conversion has to pick one. If you are storing orientations, that is harmless. If you are storing angles - incrementing them, diffing them, driving an animation curve with them - it is not, because the value you read back is not the value you wrote.

The last two rows cover the other kind of loss. Compose a small rotation ten thousand times and error accumulates in both representations - at a broadly similar rate, around $10^{-13}$.

So the honest comparison is not that one drifts less. It is what the repair costs:

• A quaternion drifts in one quantity, its length. Renormalizing is four divisions.
• A matrix drifts in nine entries, losing both unit axis lengths and their right angles. Repairing it means re-orthonormalizing - Gram-Schmidt or similar - which is more work and more code to get wrong.

That, plus the absence of poles, is the real case for holding orientation as a quaternion and converting to a matrix only when something needs to be transformed.

• Animation blending, which is thousands of small slerps a frame between bone orientations, and the main reason skeletons store quaternions.
• Turn-towards-target, where the target moves every frame so the blend restarts constantly and the angle is always small. nlerp territory.
• Cutscene and scripted cameras, which make deliberate large moves where constant speed is the whole point.
• Network interpolation, smoothing between orientation snapshots that arrive late.
• Anything storing rotation as three angles and incrementing them, which is the precision trap above and the reason engines expose angles as a convenience rather than as storage.