Floating point, precision and robust geometry

// The gap between a float32 value and the next one up, read from the bits with a DataView, at the sizes a 3D world
// uses. Math.fround rounds a JavaScript number (a float64) to the nearest float32.
const view = new DataView(new ArrayBuffer(4));
function nextFloat32(x) {
  view.setFloat32(0, x);
  view.setUint32(0, view.getUint32(0) + 1); // for a positive float32, the next bit pattern is the next value
  return view.getFloat32(0);
}
const f32Gap = (x) => nextFloat32(Math.fround(x)) - Math.fround(x);
const f64Gap = (x) => x * Number.EPSILON; // for x a power of two; close to it otherwise
const fmt = (n) => (n >= 0.01 ? String(n) : n.toExponential(3));

console.log('value         float32 gap      float64 gap   in metres, the float32 gap is about');
for (const x of [1, 100, 10000, 1000000, 16777216]) {
  const gap = f32Gap(x);
  const words = gap < 1e-4 ? `${(gap * 1e6).toFixed(2)} micrometres` : gap < 1 ? `${(gap * 1000).toFixed(2)} mm` : `${gap} m`;
  console.log(`${String(x).padEnd(12)}  ${fmt(gap).padEnd(15)}  ${fmt(f64Gap(2 ** Math.floor(Math.log2(x)))).padEnd(12)}  ${words}`);
}
console.log(`\nfloat32 of 16777217 is ${Math.fround(16777217)}: above 2^24 not every integer exists`);
console.log(`float32 of 1000000.03 is ${Math.fround(1000000.03)}: a 3 cm move at a million metres is lost`);
const world = 1000000; // two vertices 3 cm apart, 1,000 km from the world's origin
const apart = Math.fround(world + 0.03) - Math.fround(world);
console.log(`\ntwo vertices at 1000000 and 1000000.03 m, stored as float32: ${apart} m apart`);
const camera = world - 1; // render relative to a camera 1 m away: subtract in float64 first, then store as float32
const rel = Math.fround(world + 0.03 - camera) - Math.fround(world - camera);
console.log(`the same vertices stored relative to a camera at 999999 m: ${rel.toFixed(7)} m apart`);

value         float32 gap      float64 gap   in metres, the float32 gap is about
1             1.192e-7         2.220e-16     0.12 micrometres
100           7.629e-6         1.421e-14     7.63 micrometres
10000         9.766e-4         1.819e-12     0.98 mm
1000000       0.0625           1.164e-10     62.50 mm
16777216      2                3.725e-9      2 m

float32 of 16777217 is 16777216: above 2^24 not every integer exists
float32 of 1000000.03 is 1000000: a 3 cm move at a million metres is lost

two vertices at 1000000 and 1000000.03 m, stored as float32: 0 m apart
the same vertices stored relative to a camera at 999999 m: 0.0300000 m apart

// The orientation test (b - a) x (c - a) on 20,000 nearly collinear triples: c is put on the line through a and b and
// then rounded to float32, as a vertex shader's output would be. The sign is computed three ways and compared.
const f = Math.fround;
// 1. In float32, rounding after every operation, as an f32 shader computes it.
const orient32 = (a, b, c) => f(f(f(b[0] - a[0]) * f(c[1] - a[1])) - f(f(b[1] - a[1]) * f(c[0] - a[0])));
// 2. In float64, JavaScript's own numbers.
const orient64 = (a, b, c) => (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
// 3. Exactly: every float32 from 1024 to 4096 is a whole multiple of 2^-13, so scale to integers and use BigInt.
const exactSign = (a, b, c) => {
  const I = (v) => BigInt(v * 8192);
  const r = (I(b[0]) - I(a[0])) * (I(c[1]) - I(a[1])) - (I(b[1]) - I(a[1])) * (I(c[0]) - I(a[0]));
  return r > 0n ? 1 : r < 0n ? -1 : 0;
};
const sign = (x) => (x > 0 ? 1 : x < 0 ? -1 : 0);

let seed = 2026;
const rand = () => {
  seed ^= seed << 13;
  seed ^= seed >>> 17;
  seed ^= seed << 5;
  return (seed >>> 0) / 2 ** 32;
};
const coord = () => f(1024 + rand() * 3072);
let tried = 0;
let flipped32 = 0;
let zero32 = 0;
let wrong64 = 0;
let example = null;
while (tried < 20000) {
  const a = [coord(), coord()];
  const b = [coord(), coord()];
  const t = rand();
  const c = [f(a[0] + t * (b[0] - a[0])), f(a[1] + t * (b[1] - a[1]))];
  if (c[0] < 1024 || c[1] < 1024) continue;
  tried++;
  const exact = exactSign(a, b, c);
  const s32 = sign(orient32(a, b, c));
  if (s32 !== exact && s32 === 0) zero32++;
  if (s32 !== exact && s32 !== 0) {
    flipped32++;
    if (!example) example = { a, b, c, v32: orient32(a, b, c), exact };
  }
  if (sign(orient64(a, b, c)) !== exact) wrong64++;
}
console.log(`${tried} triples: float32 got the sign backwards ${flipped32} times and returned 0 for ${zero32} triples that are not collinear`);
console.log(`float64 got ${wrong64} wrong: these float32 inputs need at most about 51 bits, and a double holds 53`);
const p = ([x, y]) => `(${x}, ${y})`;
console.log(`\nfor example a = ${p(example.a)}, b = ${p(example.b)}, c = ${p(example.c)}:`);
console.log(`  float32 says ${example.v32} (left), the exact sign is ${example.exact} (right)`);

// The rasteriser's way: snap every coordinate to 1/256 of a pixel and compare whole numbers, exactly.
const snap = (v) => Math.round(v * 256);
const orientSnapped = (a, b, c) => orient64(a.map(snap), b.map(snap), c.map(snap));
const s = orientSnapped(example.a, example.b, example.c);
console.log(`  snapped to 1/256 of a pixel: ${s}, in units of 1/65536 of a square pixel, ${(s / 65536).toFixed(3)} (left)`);
console.log(`  that is exact for the snapped triangle, the one a rasteriser draws: no product exceeds 2^53`);

20000 triples: float32 got the sign backwards 89 times and returned 0 for 2419 triples that are not collinear
float64 got 0 wrong: these float32 inputs need at most about 51 bits, and a double holds 53

for example a = (3054.60986328125, 3474.336669921875), b = (1311.977294921875, 1123.3292236328125), c = (1364.0750732421875, 1193.614990234375):
  float32 says 0.25 (left), the exact sign is -1 (right)
  snapped to 1/256 of a pixel: 50944, in units of 1/65536 of a square pixel, 0.777 (left)
  that is exact for the snapped triangle, the one a rasteriser draws: no product exceeds 2^53

// Shadow rays in float32: 10,000 camera rays hit a sphere of radius 1000 at the origin, and from each hit point that
// faces the light a shadow ray tests whether the sphere itself blocks the light. Every operation rounds to float32
// (Math.fround), as an f32 shader would. A point facing the light can never be shadowed by its own convex sphere.
const f = Math.fround;
const dot = (a, b) => f(f(f(a[0] * b[0]) + f(a[1] * b[1])) + f(a[2] * b[2]));
const sub = (a, b) => [f(a[0] - b[0]), f(a[1] - b[1]), f(a[2] - b[2])];
const add = (a, b) => [f(a[0] + b[0]), f(a[1] + b[1]), f(a[2] + b[2])];
const scale = (a, s) => [f(a[0] * s), f(a[1] * s), f(a[2] * s)];
const normalize = (a) => scale(a, f(1 / f(Math.sqrt(dot(a, a)))));
const RADIUS = f(1000);

// The nearest positive t at which the ray o + t d (d of length 1) meets the sphere, or null.
function hitSphere(o, d) {
  const b = dot(o, d);
  const c = f(dot(o, o) - f(RADIUS * RADIUS));
  const disc = f(f(b * b) - c);
  if (disc < 0) return null;
  const s = f(Math.sqrt(disc));
  if (f(-b - s) > 0) return f(-b - s);
  if (f(-b + s) > 0) return f(-b + s);
  return null;
}

let seed = 1;
const rand = () => {
  seed ^= seed << 13;
  seed ^= seed >>> 17;
  seed ^= seed << 5;
  return (seed >>> 0) / 2 ** 32;
};
const eye = [0, 0, 5000];
const toLight = normalize([1, 1, 1]);
let facing = 0;
let selfShadowed = 0;
let selfShadowedWithOffset = 0;
for (let i = 0; i < 10000; i++) {
  const aim = [f((rand() - 0.5) * 1400), f((rand() - 0.5) * 1400), 0];
  const d = normalize(sub(aim, eye));
  const t = hitSphere(eye, d);
  if (t === null) continue;
  const p = add(eye, scale(d, t)); // the hit point, rounded: a little inside or outside the true surface
  const n = normalize(p);
  if (dot(n, toLight) <= 0) continue; // turned away from the light: dark anyway
  facing++;
  if (hitSphere(p, toLight) !== null) selfShadowed++;
  if (hitSphere(add(p, scale(n, 0.01)), toLight) !== null) selfShadowedWithOffset++; // start the ray 0.01 off the surface
}
console.log(`${facing} hit points face the light, so none of them should be in shadow`);
console.log(`shadow ray from the hit point itself: ${selfShadowed} report the sphere in the way (${((100 * selfShadowed) / facing).toFixed(1)} %)`);
console.log(`shadow ray started 0.01 along the normal: ${selfShadowedWithOffset} do`);

9519 hit points face the light, so none of them should be in shadow
shadow ray from the hit point itself: 4674 report the sphere in the way (49.1 %)
shadow ray started 0.01 along the normal: 0 do

// Two triangles share the edge from a to b, one on each side. A rasteriser tests a pixel centre p against that edge
// with the edge function E(a, b, p) for one triangle and E(b, a, p) for the other; exactly, E(b, a, p) = -E(a, b, p),
// so every pixel belongs to exactly one side. 40 random edges on a 4096 x 2160 target, the pixel centres along each.
const f = Math.fround;
const edge32 = (a, b, p) => f(f(f(b[0] - a[0]) * f(p[1] - a[1])) - f(f(b[1] - a[1]) * f(p[0] - a[0])));
const SUB = 256; // 8 bits of sub-pixel precision: positions snapped to 1/256 of a pixel
const snap = (v) => Math.round(v * SUB);
const edgeFixed = (a, b, p) => {
  const [A, B, P] = [a.map(snap), b.map(snap), p.map(snap)];
  return (B[0] - A[0]) * (P[1] - A[1]) - (B[1] - A[1]) * (P[0] - A[0]); // whole numbers below 2^53: exact
};

let seed = 7;
const rand = () => {
  seed ^= seed << 13;
  seed ^= seed >>> 17;
  seed ^= seed << 5;
  return (seed >>> 0) / 2 ** 32;
};
let samples = 0;
let disagree32 = 0;
let disagreeFixed = 0;
for (let k = 0; k < 40; k++) {
  const a = [f(rand() * 4096), f(rand() * 2160)];
  const b = [f(rand() * 4096), f(rand() * 2160)];
  const steps = Math.ceil(Math.hypot(b[0] - a[0], b[1] - a[1]));
  for (let i = 0; i < steps; i++) {
    const t = i / steps;
    const p = [Math.floor(a[0] + t * (b[0] - a[0])) + 0.5, Math.floor(a[1] + t * (b[1] - a[1])) + 0.5];
    samples++;
    if (edge32(a, b, p) + edge32(b, a, p) !== 0) disagree32++;
    if (edgeFixed(a, b, p) + edgeFixed(b, a, p) !== 0) disagreeFixed++;
  }
}
console.log(`${samples} pixel centres next to 40 shared edges`);
console.log(`float32: E(a, b, p) and E(b, a, p) are not exact opposites at ${disagree32} of them (${((100 * disagree32) / samples).toFixed(1)} %)`);
console.log(`fixed point, 1/${SUB} pixel: they are exact opposites at every one (${disagreeFixed} differ)`);
const group = (n) => String(n).replace(/\B(?=(\d{3})+(?!\d))/g, ',');
console.log(`largest product on this target: (4096 x ${SUB}) x (2160 x ${SUB}) = ${group(4096 * SUB * 2160 * SUB)}, below 2^53 = ${group(2 ** 53)}`);

62590 pixel centres next to 40 shared edges
float32: E(a, b, p) and E(b, a, p) are not exact opposites at 47539 of them (76.0 %)
fixed point, 1/256 pixel: they are exact opposites at every one (0 differ)
largest product on this target: (4096 x 256) x (2160 x 256) = 579,820,584,960, below 2^53 = 9,007,199,254,740,992

// The same test at a larger scale (120 shared edges, every pixel centre in a 2 x 2 block at each step along the edge):
// how often float32 lets both triangles claim a pixel, or neither, and how often 1/256-pixel fixed point does.
// About 740,000 tests: run with Node.js, the browser runner would take too long.
const f = Math.fround;
const edge32 = (a, b, p) => f(f(f(b[0] - a[0]) * f(p[1] - a[1])) - f(f(b[1] - a[1]) * f(p[0] - a[0])));
const snap = (v) => Math.round(v * 256);
const edgeFixed = (a, b, p) => {
  const [A, B, P] = [a.map(snap), b.map(snap), p.map(snap)];
  return (B[0] - A[0]) * (P[1] - A[1]) - (B[1] - A[1]) * (P[0] - A[0]);
};
let seed = 7;
const rand = () => {
  seed ^= seed << 13;
  seed ^= seed >>> 17;
  seed ^= seed << 5;
  return (seed >>> 0) / 2 ** 32;
};
const count = { samples: 0, both32: 0, neither32: 0, bothFixed: 0, neitherFixed: 0 };
for (let k = 0; k < 120; k++) {
  const a = [f(rand() * 4096), f(rand() * 2160)];
  const b = [f(rand() * 4096), f(rand() * 2160)];
  const steps = Math.ceil(Math.hypot(b[0] - a[0], b[1] - a[1]));
  for (let i = 0; i < steps; i++) {
    const t = i / steps;
    const x = a[0] + t * (b[0] - a[0]);
    const y = a[1] + t * (b[1] - a[1]);
    for (const [dx, dy] of [[0, 0], [1, 0], [0, 1], [1, 1]]) {
      const p = [Math.floor(x) + dx + 0.5, Math.floor(y) + dy + 0.5];
      count.samples++;
      const [e1, e2] = [edge32(a, b, p), edge32(b, a, p)];
      if (e1 > 0 && e2 > 0) count.both32++;
      if (e1 < 0 && e2 < 0) count.neither32++;
      const [g1, g2] = [edgeFixed(a, b, p), edgeFixed(b, a, p)];
      if (g1 > 0 && g2 > 0) count.bothFixed++;
      if (g1 < 0 && g2 < 0) count.neitherFixed++;
    }
  }
}
console.log(`${count.samples} pixel centres tested against both sides of a shared edge`);
console.log(`float32:     claimed by both triangles ${count.both32}, by neither ${count.neither32}`);
console.log(`fixed point: claimed by both triangles ${count.bothFixed}, by neither ${count.neitherFixed}`);

739864 pixel centres tested against both sides of a shared edge
float32:     claimed by both triangles 2, by neither 0
fixed point: claimed by both triangles 0, by neither 0

// Six comparisons that an exact === gets wrong, tried with three kinds of tolerance: absolute, relative and a count
// of ULPs (how many representable doubles lie between the two values).
const view = new DataView(new ArrayBuffer(8));
// A double's bits as a signed integer that grows with the value, so subtracting two of them counts the doubles between.
function ordered(x) {
  view.setFloat64(0, x);
  const bits = view.getBigInt64(0);
  return bits < 0n ? -(bits & 0x7fffffffffffffffn) : bits;
}
const ulpDistance = (a, b) => {
  const d = ordered(a) - ordered(b);
  return d < 0n ? -d : d;
};
const absolutely = (a, b, eps = 1e-9) => Math.abs(a - b) <= eps;
const relatively = (a, b, rel = 1e-12) => Math.abs(a - b) <= rel * Math.max(Math.abs(a), Math.abs(b));
const withinUlps = (a, b, n = 4n) => ulpDistance(a, b) <= n;

const f = Math.fround;
const cases = [
  ['0.1 + 0.2 against 0.3', 0.1 + 0.2, 0.3],
  ['sqrt(2) squared against 2', Math.sqrt(2) ** 2, 2],
  ['1e9 + 1e-6 against 1e9', 1e9 + 1e-6, 1e9],
  ['float32 of 0.1 against 0.1', f(0.1), 0.1],
  ['1 - (1 - 1e-12) against 1e-12', 1 - (1 - 1e-12), 1e-12],
  ['1e-20 against 0', 1e-20, 0],
];
console.log('case                                     ===    abs 1e-9  rel 1e-12  4 ULPs   ULPs apart');
for (const [name, a, b] of cases) {
  console.log(`${name.padEnd(40)} ${String(a === b).padEnd(6)} ${String(absolutely(a, b)).padEnd(9)} ${String(relatively(a, b)).padEnd(10)} ${String(withinUlps(a, b)).padEnd(8)} ${ulpDistance(a, b)}`);
}
console.log(`\n0.1 + 0.2 = ${0.1 + 0.2}; sqrt(2) squared = ${Math.sqrt(2) ** 2}`);
console.log(`float32 of 0.1 = ${f(0.1)}; 1 - (1 - 1e-12) = ${1 - (1 - 1e-12)}`);
console.log(`float32 against float64 needs a float32-sized tolerance: relative error ${(Math.abs(f(0.1) - 0.1) / 0.1).toExponential(2)}, 2^-24 = ${(2 ** -24).toExponential(2)}`);

case                                     ===    abs 1e-9  rel 1e-12  4 ULPs   ULPs apart
0.1 + 0.2 against 0.3                    false  true      true       true     1
sqrt(2) squared against 2                false  true      true       true     1
1e9 + 1e-6 against 1e9                   false  false     true       false    8
float32 of 0.1 against 0.1               false  false     false      false    107374182
1 - (1 - 1e-12) against 1e-12            false  true      false      false    109541452305
1e-20 against 0                          false  true      false      false    4307583784117748259

0.1 + 0.2 = 0.30000000000000004; sqrt(2) squared = 2.0000000000000004
float32 of 0.1 = 0.10000000149011612; 1 - (1 - 1e-12) = 9.999778782798785e-13
float32 against float64 needs a float32-sized tolerance: relative error 1.49e-8, 2^-24 = 5.96e-8

// Half floats (IEEE binary16) through ECMAScript's Float16Array and Math.f16round: what colours and positions lose.
const h = new Float16Array(1);
const viaF16 = (x) => ((h[0] = x), h[0]);

// 1. Colours: the 256 levels of an 8-bit channel, stored as k / 255 in a half float, then read back as codes.
let worstColour = 0;
let codesKept = 0;
for (let k = 0; k < 256; k++) {
  const back = viaF16(k / 255);
  worstColour = Math.max(worstColour, Math.abs(back - k / 255));
  if (Math.round(back * 255) === k) codesKept++;
}
console.log(`colours: ${codesKept} of 256 codes survive; largest error ${worstColour.toExponential(2)}, against a step of 1/255 = ${(1 / 255).toExponential(2)}`);

// 2. Positions: the gap to the next half float grows with the value, and the range ends at 65504.
const gap = (x) => {
  const v = viaF16(x);
  h[0] = v;
  const bits = new Uint16Array(h.buffer);
  bits[0] += 1; // the next half float up
  return h[0] - v;
};
for (const x of [1, 10, 100, 1000, 2048, 30000]) console.log(`at ${String(x).padStart(5)}: gap to the next half float ${gap(x)}`);
console.log(`2049 is stored as ${viaF16(2049)}; 65504 as ${viaF16(65504)}; 65520 as ${viaF16(65520)}`);
console.log(`the smallest positive half float, a subnormal: 2^-24 = ${viaF16(2 ** -24).toExponential(4)}; half of it rounds to ${viaF16(2 ** -25)}`);

// 3. Memory, and the double-rounding trap the ECMAScript specification warns about.
console.log(`\n1,000,000 positions of x, y, z: float32 ${(12e6 / 2 ** 20).toFixed(2)} MiB, half floats ${(6e6 / 2 ** 20).toFixed(2)} MiB`);
const k = 1.00048828125000022204;
console.log(`Math.f16round(k) = ${Math.f16round(k)}, but Math.f16round(Math.fround(k)) = ${Math.f16round(Math.fround(k))}`);

colours: 256 of 256 codes survive; largest error 2.43e-4, against a step of 1/255 = 3.92e-3
at     1: gap to the next half float 0.0009765625
at    10: gap to the next half float 0.0078125
at   100: gap to the next half float 0.0625
at  1000: gap to the next half float 0.5
at  2048: gap to the next half float 2
at 30000: gap to the next half float 16
2049 is stored as 2048; 65504 as 65504; 65520 as Infinity
the smallest positive half float, a subnormal: 2^-24 = 5.9605e-8; half of it rounds to 0

1,000,000 positions of x, y, z: float32 11.44 MiB, half floats 5.72 MiB
Math.f16round(k) = 1.0009765625, but Math.f16round(Math.fround(k)) = 1