255 lines
12 KiB
JavaScript
255 lines
12 KiB
JavaScript
// server/sim/detmath.js — engine-independent math for the sim. Every function here gives the same bits on every
|
||
// conforming JavaScript engine.
|
||
//
|
||
// ECMA-262 lets Math.hypot, sin, cos, atan2, pow (and **) return "implementation-approximated" values, and the engines
|
||
// differ in the last bits: V8's hypot is not the one of SpiderMonkey (fdlibm) or JavaScriptCore (sqrt(x*x + y*y));
|
||
// Chrome 154's sin / cos / atan2 are correctly rounded while Node 24–26's are fdlibm's; JavaScriptCore computes an
|
||
// integer power by repeated squaring. + − × ÷ and Math.sqrt are correctly rounded IEEE 754 operations, and the language
|
||
// fuses nothing (no FMA, no extended precision), so code built only from them is reproducible bit for bit. That is
|
||
// what these functions are. The sim must not call the approximated Math functions (eslint.config.js, and
|
||
// test/sim/detmath.test.js scans server/sim/).
|
||
//
|
||
// hypot(x, y) √(x² + y²) as Math.sqrt(x*x + y*y), rescaled by exact powers of two outside 2^±500
|
||
// powi(x, n) xⁿ for an integer n: square-and-multiply in double-double, one final rounding
|
||
// sin(x), cos(x) fdlibm 5.3 (__kernel_sin / __kernel_cos / __ieee754_rem_pio2)
|
||
// atan2(y, x) fdlibm (FreeBSD msun e_atan2.c / s_atan.c)
|
||
//
|
||
// sin / cos / atan2 are the algorithms V8 runs in Node, so they return what the server always computed (its golden
|
||
// digests stay); hypot and powi are at least as accurate as V8's. ±∞ / NaN arguments take the exactly specified
|
||
// results of the Math functions.
|
||
|
||
const F64 = new Float64Array(1);
|
||
const I32 = new Int32Array(F64.buffer);
|
||
const HW = new Uint8Array(new Uint16Array([1]).buffer)[0] === 1 ? 1 : 0; // index of the high word
|
||
const LW = 1 - HW;
|
||
/** @param {number} x */
|
||
const hi = (x) => { F64[0] = x; return I32[HW]; };
|
||
/** @param {number} x */
|
||
const lo = (x) => { F64[0] = x; return I32[LW] >>> 0; };
|
||
/** @param {number} h @param {number} l */
|
||
const fromWords = (h, l) => { I32[HW] = h; I32[LW] = l; return F64[0]; };
|
||
/** exact 2^e for a normal e @param {number} e */
|
||
const pow2 = (e) => fromWords((e + 1023) << 20, 0);
|
||
const P500 = pow2(500), PM500 = pow2(-500), P600 = pow2(600), PM600 = pow2(-600);
|
||
|
||
/**
|
||
* √(x² + y²).
|
||
* @param {number} x
|
||
* @param {number} y
|
||
*/
|
||
export function hypot(x, y) {
|
||
// eslint-disable-next-line no-restricted-properties -- ±∞ / NaN only: those results are exact by the spec
|
||
if (!Number.isFinite(x) || !Number.isFinite(y)) return Math.hypot(x, y);
|
||
const m = Math.max(Math.abs(x), Math.abs(y));
|
||
if (m > P500) { const a = x * PM600, b = y * PM600; return Math.sqrt(a * a + b * b) * P600; }
|
||
if (m < PM500 && m !== 0) { const a = x * P600, b = y * P600; return Math.sqrt(a * a + b * b) * PM600; }
|
||
return Math.sqrt(x * x + y * y);
|
||
}
|
||
|
||
// ---- integer powers ----------------------------------------------------------------------------------------------
|
||
|
||
const SPLIT = 134217729; // 2^27 + 1 (Dekker)
|
||
const ACC = [0, 0];
|
||
const BASE = [0, 0];
|
||
/** (ah + al)·(bh + bl) in double-double, into out (TwoProduct by Dekker's split: no FMA) */
|
||
function mul2(ah, al, bh, bl, out) {
|
||
const p = ah * bh;
|
||
let t = SPLIT * ah;
|
||
const ahh = t - (t - ah), ahl = ah - ahh;
|
||
t = SPLIT * bh;
|
||
const bhh = t - (t - bh), bhl = bh - bhh;
|
||
const e = ((ahh * bhh - p) + ahh * bhl + ahl * bhh) + ahl * bhl + (ah * bl + al * bh);
|
||
const s = p + e;
|
||
out[0] = s;
|
||
out[1] = e - (s - p);
|
||
}
|
||
|
||
/**
|
||
* xⁿ for an integer n (the sim's powers are stack counts, chain bounces …). The square-and-multiply runs in
|
||
* double-double and rounds once, so the result is correctly rounded in practice; a non-integer exponent or a huge
|
||
* base falls back to Math.pow. A zero base takes Math.pow's exact results: an odd power keeps the sign of −0, which
|
||
* the double-double sum would lose ((−0)¹ = −0, (−0)⁻¹ = −∞).
|
||
* @param {number} x
|
||
* @param {number} n
|
||
*/
|
||
export function powi(x, n) {
|
||
// eslint-disable-next-line no-restricted-properties -- outside the integer domain only (never in the sim)
|
||
if (!Number.isInteger(n) || !Number.isFinite(x) || Math.abs(x) > P500) return Math.pow(x, n);
|
||
if (x === 0) { const odd = n % 2 !== 0; return n === 0 ? 1 : n > 0 ? (odd ? x : 0) : (odd ? 1 / x : Infinity); }
|
||
if (n < 0) return 1 / powi(x, -n);
|
||
ACC[0] = 1; ACC[1] = 0; BASE[0] = x; BASE[1] = 0;
|
||
let k = n;
|
||
while (k > 0) {
|
||
if (k % 2 === 1) mul2(ACC[0], ACC[1], BASE[0], BASE[1], ACC);
|
||
k = Math.floor(k / 2);
|
||
if (k > 0) mul2(BASE[0], BASE[1], BASE[0], BASE[1], BASE);
|
||
}
|
||
return ACC[0] + ACC[1];
|
||
}
|
||
|
||
// ---- sin / cos (fdlibm 5.3) --------------------------------------------------------------------------------------
|
||
|
||
const S1 = -1.66666666666666324348e-01, S2 = 8.33333333332248946124e-03, S3 = -1.98412698298579493134e-04;
|
||
const S4 = 2.75573137070700676789e-06, S5 = -2.50507602534068634195e-08, S6 = 1.58969099521155010221e-10;
|
||
/** sin on [-π/4, π/4]; y is the tail of x, iy = 0 when y is 0 */
|
||
function kSin(x, y, iy) {
|
||
if ((hi(x) & 0x7fffffff) < 0x3e400000 && Math.trunc(x) === 0) return x;
|
||
const z = x * x, v = z * x;
|
||
const r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
|
||
if (iy === 0) return x + v * (S1 + z * r);
|
||
return x - ((z * (0.5 * y - v * r) - y) - v * S1);
|
||
}
|
||
|
||
const C1 = 4.16666666666666019037e-02, C2 = -1.38888888888741095749e-03, C3 = 2.48015872894767294178e-05;
|
||
const C4 = -2.75573143513906633035e-07, C5 = 2.08757232129817482790e-09, C6 = -1.13596475577881948265e-11;
|
||
/** cos on [-π/4, π/4]; y is the tail of x */
|
||
function kCos(x, y) {
|
||
const ix = hi(x) & 0x7fffffff;
|
||
if (ix < 0x3e400000 && Math.trunc(x) === 0) return 1;
|
||
const z = x * x;
|
||
const r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
|
||
if (ix < 0x3fd33333) return 1 - (0.5 * z - (z * r - x * y));
|
||
const qx = ix > 0x3fe90000 ? 0.28125 : fromWords(ix - 0x00200000, 0);
|
||
return (1 - qx) - ((0.5 * z - qx) - (z * r - x * y));
|
||
}
|
||
|
||
const NPIO2_HW = [
|
||
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A,
|
||
0x4031475C, 0x4032D97C, 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 0x403DD85A, 0x403F6A7A,
|
||
0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
|
||
0x404858EB, 0x404921FB,
|
||
];
|
||
const INV_PIO2 = 6.36619772367581382433e-01;
|
||
const PIO2_1 = 1.57079632673412561417e+00, PIO2_1T = 6.07710050650619224932e-11;
|
||
const PIO2_2 = 6.07710050630396597660e-11, PIO2_2T = 2.02226624879595063154e-21;
|
||
const PIO2_3 = 2.02226624871116645580e-21, PIO2_3T = 8.47842766036889956997e-32;
|
||
const TWO_PI = 6.283185307179586;
|
||
const Y = [0, 0];
|
||
|
||
/** x − n·π/2 into Y (head, tail); returns n. |x| ≤ 2^19·π/2: fdlibm's medium path; beyond that (never in the sim)
|
||
* x is first taken modulo the double nearest 2π (exact %), which stays deterministic though not accurate. */
|
||
function remPio2(x0) {
|
||
let x = x0;
|
||
if ((hi(x) & 0x7fffffff) > 0x413921fb) x = x % TWO_PI;
|
||
const hx = hi(x), ix = hx & 0x7fffffff;
|
||
if (ix <= 0x3fe921fb) { Y[0] = x; Y[1] = 0; return 0; }
|
||
if (ix < 0x4002d97c) { // |x| < 3π/4: n = ±1
|
||
if (hx > 0) {
|
||
let z = x - PIO2_1;
|
||
if (ix !== 0x3ff921fb) { Y[0] = z - PIO2_1T; Y[1] = (z - Y[0]) - PIO2_1T; } else { z -= PIO2_2; Y[0] = z - PIO2_2T; Y[1] = (z - Y[0]) - PIO2_2T; }
|
||
return 1;
|
||
}
|
||
let z = x + PIO2_1;
|
||
if (ix !== 0x3ff921fb) { Y[0] = z + PIO2_1T; Y[1] = (z - Y[0]) + PIO2_1T; } else { z += PIO2_2; Y[0] = z + PIO2_2T; Y[1] = (z - Y[0]) + PIO2_2T; }
|
||
return -1;
|
||
}
|
||
let t = Math.abs(x);
|
||
const n = Math.trunc(t * INV_PIO2 + 0.5);
|
||
let r = t - n * PIO2_1, w = n * PIO2_1T;
|
||
Y[0] = r - w;
|
||
if (!(n < 32 && ix !== NPIO2_HW[n - 1])) {
|
||
const j = ix >> 20;
|
||
if (j - ((hi(Y[0]) >> 20) & 0x7ff) > 16) { // a 2nd step, good to 118 bits
|
||
t = r; w = n * PIO2_2; r = t - w; w = n * PIO2_2T - ((t - r) - w); Y[0] = r - w;
|
||
if (j - ((hi(Y[0]) >> 20) & 0x7ff) > 49) { // a 3rd, 151 bits
|
||
t = r; w = n * PIO2_3; r = t - w; w = n * PIO2_3T - ((t - r) - w); Y[0] = r - w;
|
||
}
|
||
}
|
||
}
|
||
Y[1] = (r - Y[0]) - w;
|
||
if (hx < 0) { Y[0] = -Y[0]; Y[1] = -Y[1]; return -n; }
|
||
return n;
|
||
}
|
||
|
||
/** @param {number} x radians */
|
||
export function sin(x) {
|
||
if (!Number.isFinite(x)) return NaN;
|
||
if ((hi(x) & 0x7fffffff) <= 0x3fe921fb) return kSin(x, 0, 0);
|
||
switch (remPio2(x) & 3) {
|
||
case 0: return kSin(Y[0], Y[1], 1);
|
||
case 1: return kCos(Y[0], Y[1]);
|
||
case 2: return -kSin(Y[0], Y[1], 1);
|
||
default: return -kCos(Y[0], Y[1]);
|
||
}
|
||
}
|
||
|
||
/** @param {number} x radians */
|
||
export function cos(x) {
|
||
if (!Number.isFinite(x)) return NaN;
|
||
if ((hi(x) & 0x7fffffff) <= 0x3fe921fb) return kCos(x, 0);
|
||
switch (remPio2(x) & 3) {
|
||
case 0: return kCos(Y[0], Y[1]);
|
||
case 1: return -kSin(Y[0], Y[1], 1);
|
||
case 2: return -kCos(Y[0], Y[1]);
|
||
default: return kSin(Y[0], Y[1], 1);
|
||
}
|
||
}
|
||
|
||
// ---- atan2 (fdlibm, FreeBSD msun) --------------------------------------------------------------------------------
|
||
|
||
const ATAN_HI = [4.63647609000806093515e-01, 7.85398163397448278999e-01, 9.82793723247329054082e-01, 1.57079632679489655800e+00];
|
||
const ATAN_LO = [2.26987774529616870924e-17, 3.06161699786838301793e-17, 1.39033110312309984516e-17, 6.12323399573676603587e-17];
|
||
const AT = [
|
||
3.33333333333329318027e-01, -1.99999999998764832476e-01, 1.42857142725034663711e-01, -1.11111104054623557880e-01,
|
||
9.09088713343650656196e-02, -7.69187620504482999495e-02, 6.66107313738753120669e-02, -5.83357013379057348645e-02,
|
||
4.97687799461593236017e-02, -3.65315727442169155270e-02, 1.62858201153657823623e-02,
|
||
];
|
||
/** @param {number} x0 */
|
||
function atan(x0) {
|
||
let x = x0;
|
||
const hx = hi(x), ix = hx & 0x7fffffff;
|
||
let id;
|
||
if (ix >= 0x44100000) { // |x| ≥ 2^66
|
||
if (ix > 0x7ff00000 || (ix === 0x7ff00000 && lo(x) !== 0)) return x + x;
|
||
return hx > 0 ? ATAN_HI[3] + ATAN_LO[3] : -ATAN_HI[3] - ATAN_LO[3];
|
||
}
|
||
if (ix < 0x3fdc0000) { // |x| < 0.4375
|
||
if (ix < 0x3e200000) return x;
|
||
id = -1;
|
||
} else {
|
||
x = Math.abs(x);
|
||
if (ix < 0x3ff30000) {
|
||
if (ix < 0x3fe60000) { id = 0; x = (2 * x - 1) / (2 + x); } else { id = 1; x = (x - 1) / (x + 1); }
|
||
} else if (ix < 0x40038000) { id = 2; x = (x - 1.5) / (1 + 1.5 * x); } else { id = 3; x = -1 / x; }
|
||
}
|
||
const z = x * x, w = z * z;
|
||
const s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10])))));
|
||
const s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9]))));
|
||
if (id < 0) return x - x * (s1 + s2);
|
||
const r = ATAN_HI[id] - ((x * (s1 + s2) - ATAN_LO[id]) - x);
|
||
return hx < 0 ? -r : r;
|
||
}
|
||
|
||
const PI_O_4 = 7.8539816339744827900E-01, PI_O_2 = 1.5707963267948965580E+00, PI = 3.1415926535897931160E+00;
|
||
const PI_LO = 1.2246467991473532e-16;
|
||
/**
|
||
* The angle of (x, y), in (−π, π].
|
||
* @param {number} y
|
||
* @param {number} x
|
||
*/
|
||
export function atan2(y, x) {
|
||
if (Number.isNaN(x) || Number.isNaN(y)) return NaN;
|
||
const hx = hi(x), ix = hx & 0x7fffffff, lx = lo(x);
|
||
const hy = hi(y), iy = hy & 0x7fffffff, ly = lo(y);
|
||
if (((hx - 0x3ff00000) | lx) === 0) return atan(y); // x = 1
|
||
let m = ((hy >>> 31) & 1) | ((hx >>> 30) & 2); // 2·sign(x) + sign(y)
|
||
if ((iy | ly) === 0) return m === 2 ? PI : m === 3 ? -PI : y; // y = ±0
|
||
if ((ix | lx) === 0) return hy < 0 ? -PI_O_2 : PI_O_2; // x = ±0
|
||
if (ix === 0x7ff00000) { // x = ±∞
|
||
if (iy === 0x7ff00000) return [PI_O_4, -PI_O_4, 3 * PI_O_4, -3 * PI_O_4][m];
|
||
return [0, -0, PI, -PI][m];
|
||
}
|
||
if (iy === 0x7ff00000) return hy < 0 ? -PI_O_2 : PI_O_2; // y = ±∞
|
||
const k = (iy - ix) >> 20;
|
||
let z;
|
||
if (k > 60) { z = PI_O_2 + 0.5 * PI_LO; m &= 1; } // |y/x| > 2^60
|
||
else if (hx < 0 && k < -60) z = 0; // |y|/x < −2^60
|
||
else z = atan(Math.abs(y / x));
|
||
switch (m) {
|
||
case 0: return z;
|
||
case 1: return -z;
|
||
case 2: return PI - (z - PI_LO);
|
||
default: return (z - PI_LO) - PI;
|
||
}
|
||
}
|