Test exact dtoa output and compact powers

Author: JairusSW

tests/dtoa/format.ts

@@ -3,6 +3,15 @@ import {
 } from "~lib/util/dtoa";
 import { dtoa as numberDtoa, dtoa_buffered as numberBuffered } from "util/number";
 
+// Called by the JS differential harness with exact input bit patterns.
+export function format64Bits(bits: u64, buffer: usize): u32 {
+  return numberBuffered<f64>(buffer, reinterpret<f64>(bits));
+}
+
+export function format32Bits(bits: u32, buffer: usize): u32 {
+  return numberBuffered<f32>(buffer, reinterpret<f32>(bits));
+}
+
 toDigits32(1);
 assert(gDigits == 1);
 assert(gDigHi == 0x3130303030303030); // "00000001"

tests/dtoa/index.js

@@ -1,4 +1,102 @@
 import asc from "../../dist/asc.js";
+import { shortestFloat32, withAssemblyScriptDotZero } from "./oracle.js";
+import { verifyPowers } from "./powers.js";
+
+verifyPowers();
+
+const f64View = new DataView(new ArrayBuffer(8));
+const f32View = new DataView(new ArrayBuffer(4));
+const sign64 = 1n << 63n;
+const maxPositive64 = sign64 - 1n;
+let randomState = 0x12345678;
+
+function random32() {
+  randomState ^= randomState << 13;
+  randomState ^= randomState >>> 17;
+  randomState ^= randomState << 5;
+  return randomState >>> 0;
+}
+
+function bits64(value) {
+  f64View.setFloat64(0, value);
+  return f64View.getBigUint64(0);
+}
+
+function bits32(value) {
+  f32View.setFloat32(0, value);
+  return f32View.getUint32(0);
+}
+
+function makeCases64() {
+  const patterns = new Set();
+  const addSigned = positive => {
+    if (positive < 0n || positive > maxPositive64) return;
+    patterns.add(positive);
+    patterns.add(positive | sign64);
+  };
+  for (let exponent = 0; exponent <= 2047; ++exponent) {
+    const base = BigInt(exponent) << 52n;
+    for (const significand of [0n, 1n, 15n, (1n << 52n) - 1n]) {
+      addSigned(base | significand);
+    }
+  }
+  for (let exponent = -323; exponent <= 308; ++exponent) {
+    const base = bits64(Number(`1e${exponent}`));
+    for (let delta = -2; delta <= 2; ++delta) addSigned(base + BigInt(delta));
+  }
+  for (let i = 0; i < 2048; ++i) {
+    patterns.add((BigInt(random32()) << 32n) | BigInt(random32()));
+  }
+  return [...patterns].map(bits => {
+    f64View.setBigUint64(0, bits);
+    return [bits, withAssemblyScriptDotZero(f64View.getFloat64(0).toString())];
+  });
+}
+
+function makeCases32() {
+  const patterns = new Set();
+  const addSigned = positive => {
+    if (positive < 0 || positive > 0x7fffffff) return;
+    patterns.add(positive >>> 0);
+    patterns.add((positive | 0x80000000) >>> 0);
+  };
+  for (let exponent = 0; exponent <= 255; ++exponent) {
+    const base = exponent << 23;
+    for (const significand of [0, 1, 15, 0x7fffff]) {
+      addSigned(base | significand);
+    }
+  }
+  for (let exponent = -45; exponent <= 38; ++exponent) {
+    const base = bits32(10 ** exponent);
+    for (let delta = -2; delta <= 2; ++delta) addSigned(base + delta);
+  }
+  for (let i = 0; i < 2048; ++i) patterns.add(random32());
+  return [...patterns].map(bits => [bits, withAssemblyScriptDotZero(shortestFloat32(bits))]);
+}
+
+const cases64 = makeCases64();
+const cases32 = makeCases32();
+
+function compareExactOutput(exports, mode) {
+  const { memory, format64Bits, format32Bits } = exports;
+  const buffer = memory.grow(1) * 65536;
+  for (const [bits, expected] of cases64) {
+    const length = format64Bits(bits, buffer);
+    if (length > 32) throw new Error(`${mode} f64 length ${length} exceeds buffer capacity`);
+    const actual = String.fromCharCode(...new Uint16Array(memory.buffer, buffer, length));
+    if (actual !== expected) {
+      throw new Error(`${mode} f64 0x${bits.toString(16)}: ${JSON.stringify(actual)} != ${JSON.stringify(expected)}`);
+    }
+  }
+  for (const [bits, expected] of cases32) {
+    const length = format32Bits(bits, buffer);
+    if (length > 32) throw new Error(`${mode} f32 length ${length} exceeds buffer capacity`);
+    const actual = String.fromCharCode(...new Uint16Array(memory.buffer, buffer, length));
+    if (actual !== expected) {
+      throw new Error(`${mode} f32 0x${bits.toString(16)}: ${JSON.stringify(actual)} != ${JSON.stringify(expected)}`);
+    }
+  }
+}
 
 for (const optimized of [false, true]) {
   for (const simd of [false, true]) {
@@ -31,6 +129,7 @@ for (const optimized of [false, true]) {
       },
     });
     instance.exports.start();
-    console.log(`dtoa ${mode}: passed`);
+    compareExactOutput(instance.exports, mode);
+    console.log(`dtoa ${mode}: ${cases64.length} f64 and ${cases32.length} f32 exact outputs passed`);
   }
 }

tests/dtoa/oracle.js

@@ -0,0 +1,74 @@
+// Independent f32 shortest-decimal oracle, adapted from xjb-as's BigInt
+// verifier. It works from the exact binary32 value, not a widened f64 string.
+const bits32 = new DataView(new ArrayBuffer(4));
+
+export function float32FromBits(bits) {
+  bits32.setUint32(0, bits >>> 0);
+  return bits32.getFloat32(0);
+}
+
+function formatDecimal(negative, digits, point) {
+  const count = digits.length;
+  let result;
+  if (count <= point && point <= 21) {
+    result = digits + "0".repeat(point - count);
+  } else if (0 < point && point <= 21) {
+    result = digits.slice(0, point) + "." + digits.slice(point);
+  } else if (-6 < point && point <= 0) {
+    result = "0." + "0".repeat(-point) + digits;
+  } else {
+    const exponent = point - 1;
+    const suffix = "e" + (exponent >= 0 ? "+" : "-") + Math.abs(exponent);
+    result = count === 1 ? digits + suffix : digits[0] + "." + digits.slice(1) + suffix;
+  }
+  return (negative ? "-" : "") + result;
+}
+
+export function shortestFloat32(bits) {
+  const value = float32FromBits(bits);
+  if (Number.isNaN(value)) return "NaN";
+  if (value === Infinity) return "Infinity";
+  if (value === -Infinity) return "-Infinity";
+  if (value === 0) return "0";
+
+  const negative = (bits >>> 31) !== 0;
+  const exponentBits = (bits >>> 23) & 0xff;
+  const mantissa = bits & 0x7fffff;
+  const significand = BigInt(exponentBits ? mantissa | 0x800000 : mantissa);
+  const binaryExponent = exponentBits ? exponentBits - 150 : -149;
+  const decimalInteger = binaryExponent >= 0
+    ? significand << BigInt(binaryExponent)
+    : significand * 5n ** BigInt(-binaryExponent);
+  const decimalExponent = Math.min(binaryExponent, 0);
+  const exactDigits = decimalInteger.toString();
+  const magnitude = Math.abs(value);
+
+  // At each length, test the two decimals bracketing the exact value. Choose
+  // the closest one that parses back to the binary32 value; ties go to even.
+  for (let count = 1; count <= exactDigits.length; ++count) {
+    const shift = exactDigits.length - count + decimalExponent;
+    const floor = exactDigits.slice(0, count);
+    const ceil = (BigInt(floor) + 1n).toString();
+    const roundTrips = digits => Math.fround(Number(digits + "e" + shift)) === magnitude;
+    const floorOk = roundTrips(floor);
+    const ceilOk = roundTrips(ceil);
+    if (!floorOk && !ceilOk) continue;
+
+    let chosen;
+    if (floorOk && ceilOk) {
+      const unit = 10n ** BigInt(exactDigits.length - count);
+      const distance = 2n * (decimalInteger % unit);
+      chosen = distance < unit ? floor : distance > unit ? ceil
+        : (Number(floor[floor.length - 1]) & 1) === 0 ? floor : ceil;
+    } else {
+      chosen = floorOk ? floor : ceil;
+    }
+    const point = chosen.length + shift;
+    return formatDecimal(negative, chosen.replace(/0+$/, ""), point);
+  }
+  throw new Error(`no shortest f32 decimal for 0x${bits.toString(16)}`);
+}
+
+export function withAssemblyScriptDotZero(text) {
+  return /^-?\d+$/.test(text) ? text + ".0" : text;
+}

tests/dtoa/powers.js

@@ -0,0 +1,74 @@
+import { readFileSync } from "node:fs";
+
+const source = readFileSync(new URL("../../std/assembly/util/dtoa.ts", import.meta.url), "utf8");
+const mask64 = (1n << 64n) - 1n;
+
+function table(name) {
+  const match = new RegExp(`\\b${name}\\s*=\\s*memory\\.data<[^>]+>\\(\\[([\\s\\S]*?)\\]\\)`).exec(source);
+  if (!match) throw new Error(`missing ${name}`);
+  return [...match[1].matchAll(/0x[0-9a-f]+/gi)].map(([hex]) => BigInt(hex));
+}
+
+function exactPower(index) {
+  const exponent = index - 293;
+  const power = 10n ** BigInt(Math.abs(exponent));
+  const bitCount = power.toString(2).length;
+  return exponent >= 0
+    ? bitCount <= 128 ? power << BigInt(128 - bitCount) : power >> BigInt(bitCount - 128)
+    : (1n << BigInt(127 + bitCount)) / power;
+}
+
+export function verifyPowers() {
+  const majorTable = table("POW10_MAJOR");
+  const fixups = table("POW10_FIXUPS");
+  const shared = table("POW10_SHARED");
+  if (majorTable.length !== 46 || fixups.length !== 20 || shared.length !== 53) {
+    throw new Error("power-table lengths changed; update the verifier with the layout");
+  }
+
+  // Shared highs cover powers 10^27..10^-17. Their rounding direction is
+  // important to the f32 path, even where the corresponding low is not stored.
+  for (let t = 17; t <= 61; ++t) {
+    const exact = exactPower(337 - t);
+    const high = exact >> 64n;
+    const low = exact & mask64;
+    if (shared[t - 17] !== high + BigInt(low !== 0n)) {
+      throw new Error(`shared high mismatch at t=${t}`);
+    }
+    if (t >= 45 && t <= 52 && shared[t] !== low) {
+      throw new Error(`shared low mismatch at t=${t}`);
+    }
+  }
+
+  // Reproduce the compact reconstruction using BigInt multiplication, then
+  // compare all 618 entries against exact powers of ten. The reference above
+  // does not use any cached factors or correction bits.
+  for (let index = 0; index < 618; ++index) {
+    const j = index + 10;
+    const major = (j * 293) >>> 13;
+    if (major !== Math.floor(j / 28)) throw new Error(`quotient mismatch at ${j}`);
+    const minor = j - major * 28;
+    const factor = shared[27 - minor];
+    const anchor = (majorTable[major * 2] << 64n) | majorTable[major * 2 + 1];
+    const product = anchor * factor;
+    const c0 = product & mask64;
+    const c1 = (product >> 64n) & mask64;
+    const c2 = product >> 128n;
+    let high;
+    let low;
+    if (c2 >> 63n) {
+      high = c2;
+      low = c1;
+    } else {
+      high = (c2 << 1n) | (c1 >> 63n);
+      low = ((c1 << 1n) | (c0 >> 63n)) & mask64;
+    }
+    const fixup = (fixups[index >> 5] >> BigInt(index & 31)) & 1n;
+    low = (low - fixup) & mask64;
+    const actual = (high << 64n) | low;
+    const expected = exactPower(index);
+    if (actual !== expected) {
+      throw new Error(`power mismatch at ${index}: ${actual.toString(16)} != ${expected.toString(16)}`);
+    }
+  }
+}