nmadd, modMul, modMul53, random

f64math.js

function nmadd(a, b, p) {
  const at = (2**27 + 1) * a;
  const ahi = at - (at - a);
  const alo = a - ahi;
  const bt = (2**27 + 1) * b;
  const bhi = bt - (bt - b);
  const blo = b - bhi;
  return ((p - ahi * bhi) - ahi * blo - alo * bhi) - alo * blo;
}

function modMul(a, b, m, mInv) {
  const p = a * b;
  const r1 = nmadd(a, b, p);
  const q = Math.floor(p * mInv);
  const r2 = nmadd(q, m, p);
  return r2 < (r1 + m) ? (r2 < r1 ? r2 - (r1 - m) : r2 - r1) : r2 - (r1 + m);
}

function modMul53(a, b, m, mInv) {
  const p = a * b;
  const r1 = nmadd(a, b, p);
  const q = Math.floor(p * mInv + 0.5);
  const r2 = nmadd(q, m, p);
  const r = r2 - r1;
  return (r - m < 0.0 - r ? (r + m < 0.0 - r ? r + m : r) : r - m);
}

function random(seed = -0.0 + 1.0) {
  seed = -0.0 + seed;
  return function () {
    // the parameters are from https://www.ams.org/journals/mcom/1999-68-225/S0025-5718-99-00996-5/S0025-5718-99-00996-5.pdf
    seed = seed * (-0.0 + 185852.0);
    seed = seed - Math.floor(seed * (1.0 / 34359738337.0)) * 34359738337.0;
    return seed * (1.0 / 34359738337.0);
  };
}