Function fplnDirect

src/big-decimal/utils.ts:391–429

* Direct-log fixed-point natural logarithm (machine-seed + log1p correction). * * ln(v) = y₀ + log1p(v·e^{−y₀} − 1), where y₀ = Math.log(v) is a ~48-bit-accurate * float64 seed. Because v·e^{−y₀} = 1 + ε with ε ≈ 2^{−48}, the correction * log1p(ε) = 2·atanh(ε/(2+ε)) needs only ≈ bits/96 series t

(x: bigint, bits: number)

 * @returns  ln(x/2^bits) * 2^bits as a bigint.
 */
function fplnDirect(x: bigint, bits: number): bigint {
  const B = BigInt(bits);
  const scale = 1n << B;

  // Overflow-safe ~48-bit float seed y₀ ≈ ln(x/scale)·scale. Extract the top
  // ~52 bits of x so Number() stays finite even when 2^bits overflows a double
  // (bits ≥ 1024); the log1p correction restores full precision regardless.
  const topBits = 52;
  const sh = bits - topBits;
  let y0Fp: bigint;
  if (sh > 0) {
    const v = Number(x >> BigInt(sh)) / 2 ** topBits;
    y0Fp = BigInt(Math.round(Math.log(v) * 2 ** topBits)) << BigInt(sh);
  } else {
    const v = Number(x) / Number(scale);
    y0Fp = BigInt(Math.round(Math.log(v) * Number(scale)));
  }

  // p = x·e^{−y₀} = (1 + ε)·scale, ε ≈ 2^{−48}. The single fpexp call.
  const p = (x * fpexp(-y0Fp, bits)) >> B;
  const eps = p - scale; // ε·scale (may be negative if y₀ overshot ln v)

  // log1p(ε) = 2·atanh(u), u = ε/(2+ε). The atanh terms all share u's sign, so
  // sum in magnitude: an arithmetic right-shift of a negative bigint floors
  // toward −∞ and sticks at −1 instead of reaching 0, which would never
  // terminate the loop.
  const u = (eps << B) / (p + scale);
  const uAbs = u < 0n ? -u : u;
  const u2 = (uAbs * uAbs) >> B;
  let term = uAbs;
  let sum = uAbs;
  for (let n = 3n; ; n += 2n) {
    term = (term * u2) >> B;
    if (term === 0n) break;
    sum += term / n;
  }

  return y0Fp + (u < 0n ? -2n : 2n) * sum;
}

/**
 * ln(2) · 2^bits, cached. Computed once via the Newton path (NOT the AGM path,