multiply bigint by float in JS (spoiler alert: BigInt has no useful functions)

multiplyFloat.js

/**
 * For a finite normal 64-bit float `f`, extracts integers `sgn`,
 * `exponent`, and `mantissa` such that:
 *
 *   - `sgn` is -1 or +1
 *   - `exponent` is between -1023 and 1024, inclusive
 *   - `mantissa` is between 0 and 2^51 - 1, inclusive
 *   - the number given by `f` equals `sgn * 2^exponent * (1 + mantissa / 2^52)`
 *
 * The results are all bigints within the range of safe integers for
 * 64-bit floats (i.e., converting them to `Number` is lossless).
 *
 * Throws an error if `f` is subnormal (biased exponent is 0).
 */
function decomposeFloat(f: number) {
  if (!isFinite(f)) {
    throw new Error("Input must be finite: " + f);
  }
  const union = new DataView(new ArrayBuffer(8));
  const littleEndian = true; // arbitrary, but faster when matches native arch
  union.setFloat64(0, f, littleEndian);
  const bytes = union.getBigUint64(0, littleEndian);
  const sgn = (-1n) ** (bytes >> 63n);
  const biasedExponent = (bytes & ~(1n << 63n)) >> 52n;
  if (biasedExponent === 0n) {
    throw new Error("Subnormal floats not supported: " + f);
  }
  const exponent = biasedExponent - 1023n;
  const mantissa = bytes & ((1n << 52n) - 1n);
  return {sgn, exponent, mantissa};
}

/**
 * Multiply an exact bigint by a floating point number.
 *
 * This function is exact in the first argument and subject to the usual
 * floating point considerations in the second. Thus, it is the case
 * that
 *
 *      multiplyFloat(g, 1) === g
 *      multiplyFloat(g, x) + multiplyFloat(h, f) === multiplyFloat(g + h, x)
 *      multiplyFloat(k * g, x) === k * multiplyFloat(g, x)
 *
 * for all `BigInt`s `k`, `g`, and `h` and all floats `x`. But it is not
 * necessarily the case that
 *
 *      multiplyFloat(g, x) + multiplyFloat(g, y) === multiplyFloat(g, x + y)
 *
 * for all `BigInt`s `g` and floats `x` and `y`: e.g., when `x === 1`
 * and `y === 1e-16`, we have `x + y === x` even though `y !== 0`.
 */
export function multiplyFloat(g: BigInt, fac: number): BigInt {
  if (fac === 0) {
    // Special case, as 0 is subnormal.
    return 0n;
  }
  const {sgn, exponent, mantissa} = decomposeFloat(fac);
  // from `decomposeFloat` contract, `fac = numerator / denominator`
  // exactly (in arbitrary-precision arithmetic)
  const numerator = sgn * 2n ** (exponent + 1023n) * (2n ** 52n + mantissa);
  const denominator = 2n ** (1023n + 52n);
  // round to nearest, biasing toward zero on exact tie
  return (2n * numerator * g + sgn * denominator) / (2n * denominator);
}

// Tests and discussion at:
// <https://github.com/sourcecred/sourcecred/pull/1715#issuecomment-603354146>

Would you be willing to give this file an open-source license?

Certainly. Please feel free to use it under either MIT or Apache2, at your option. Glad you've found it useful.

Thank you! Here's my version, which primarily fixes the Typescript types by specifying bigint instead of BigInt: https://gist.github.com/NoahAndrews/eb883c5a8a6d5805e76227b02cc84a92