Function fpexp

src/big-decimal/utils.ts:276–326

(x: bigint, bits: number)

 * @returns  exp(x/2^bits) * 2^bits as a bigint
 */
export function fpexp(x: bigint, bits: number): bigint {
  // exp(0) = 1
  if (x === 0n) return 1n << BigInt(bits);

  // Elevated internal scale. `targetHalvings` extra halvings shorten the Taylor
  // series to O(√bits) terms; `guard` guard bits absorb the ~2^targetHalvings
  // error amplification of the squaring-back phase (and the final downshift's
  // rounding).
  const targetHalvings = Math.round(EXP_REDUCE_COEF * Math.sqrt(bits));
  const guard = targetHalvings + 8;
  const bits2 = bits + guard;
  const B = BigInt(bits2);
  const scale = 1n << B;

  // Lift the input to the elevated scale.
  const g = BigInt(guard);
  let r = x << g;

  // Argument reduction: halve until |r/scale| < 2^−targetHalvings. The shift
  // count is derived from the ACTUAL magnitude of r, so a small argument (e.g.
  // fpln's Newton seed near ln 1 = 0) is never over-reduced to 0 — the shift
  // always leaves ~bits2 − targetHalvings significant bits. Sign is carried
  // separately (a right-shift of a negative bigint would round toward −∞).
  const neg = r < 0n;
  let ra = neg ? -r : r;
  const kExtra = Math.max(0, bitLength(ra) - bits2 + targetHalvings);
  ra >>= BigInt(kExtra);
  r = neg ? -ra : ra;

  // Taylor series: exp(r/scale) = 1 + r/scale + r²/(2!·scale²) + ...
  // In fixed-point: sum = scale + r + r²/(2·scale) + r³/(6·scale²) + ...
  // Incremental: term_n = term_{n-1} * r / (n * scale)
  //   base-2: ((term * r) >> bits2) / n   — shift + small-divisor division
  let sum = scale; // 1.0
  let term = r; // r/scale in fixed-point
  sum += term;

  for (let n = 2; ; n++) {
    term = ((term * r) >> B) / BigInt(n);
    if (term === 0n) break;
    sum += term;
  }

  // Squaring phase: exp(x/scale) = exp((x/scale)/2^kExtra)^(2^kExtra).
  for (let i = 0; i < kExtra; i++) {
    sum = (sum * sum) >> B;
  }

  // Downshift from the elevated scale back to `bits`, rounding to nearest.
  return (sum + (1n << (g - 1n))) >> g;
}

// AGM-vs-direct-log crossover, in bits. Below LN_AGM_MIN_BITS the direct-log
// kernel (`fplnDirect`: one √-reduced fpexp + a short log1p series) wins; above