Function fpsqrt

src/big-decimal/utils.ts:87–119

(a: bigint, bits: number)

 * Converge until |x_{n+1} - x_n| <= 1 (one ULP in the fixed-point grid).
 */
export function fpsqrt(a: bigint, bits: number): bigint {
  if (a === 0n) return 0n;
  if (a < 0n) throw new RangeError('fpsqrt: negative input');

  // as = a * scale = a << bits; the result is isqrt(as).
  const as = a << BigInt(bits);

  // Seed the final refinement. At low/medium precision the float-seeded Heron
  // converges in a few full-width divisions and wins outright; at high
  // precision a recursive giant-steps isqrt (root the top half, then refine)
  // does ~3× fewer full-width divisions. Dispatch on `bits` (a plain number
  // compare, no bitLength) so the hot low-precision path is unchanged — fpsqrt
  // callers pass a ≈ O(1)·2^bits, so as ≈ 2^(2·bits) and bits ≥ FP_SQRT_GIANT_BITS
  // ⇒ as ≳ SQRT_GIANT_MIN_BITS (and isqrtGiant falls back to flat Heron if not).
  let x: bigint;
  if (bits < FP_SQRT_GIANT_BITS) {
    x = bigSqrtSeed(as);
    let prev: bigint;
    do {
      prev = x;
      x = (x + as / x) / 2n;
    } while (bigintAbs(x - prev) > 1n);
  } else {
    x = isqrtGiant(as, bitLength(as));
  }

  // One more iteration, then pick whichever of {x, next} has x²
  // closest to `as` (the true floor-root or one above it).
  const next = (x + as / x) / 2n;
  const diffX = bigintAbs(x * x - as);
  const diffNext = bigintAbs(next * next - as);
  return diffNext < diffX ? next : x;
}

/**
 * Seed for an integer square root of `n`: a value within a few bits of