Function digammaCore

src/compute-engine/numerics/special-functions.ts:829–880

(ce: ComputeEngine, z: BigNum)

  const table = computeBernoulliEven(buildTerms, ce._deadline);
  bernoulliCache.set(ce, table);

  // Return the full cached table (length ≥ minTerms), NOT a `minTerms` slice.
  // Every caller already clamps its own consumption — the four asymptotic
  // series (gammaln/digamma/trigamma/polygamma) iterate to
  // `Math.min(maxTerms, table.length)` and break at the tolerance well before
  // that, and `zetaCore` indexes a single fixed element `table[k-1]` — so no
  // caller ever over-sums past its requested `minTerms`. The asymptotic-series
  // divergence trap (extra terms make the sum worse) is therefore handled at
  // the caller's clamp, and the per-call `slice` copy is pure overhead. (The
  // clamp, not the table length, is what enforces the asymptotic cutoff.)
  return table;
}

/**
 * Bignum Digamma function ψ(z) = d/dz ln(Γ(z))
 * Same algorithm as machine `digamma`: reflection for negative z,
 * recurrence to shift z > 7, then asymptotic expansion with Bernoulli numbers.
 */
export function bigDigamma(ce: ComputeEngine, z: BigNum): BigNum {
  if (!z.isFinite()) return BigDecimal.NAN;
  return withGuardDigits(SPECIAL_FN_GUARD, () => digammaCore(ce, z));
}

function digammaCore(ce: ComputeEngine, z: BigNum): BigNum {
  // Reflection formula for negative values: ψ(1-z) = ψ(z) + π·cot(πz)
  if (z.isNegative()) {
    if (z.isInteger()) return BigDecimal.NAN; // poles at non-positive integers
    const pi = BigDecimal.PI;
    const piZ = pi.mul(z);
    const cotPiZ = piZ.cos().div(piZ.sin());
    return digammaCore(ce, BigDecimal.ONE.sub(z)).sub(pi.mul(cotPiZ));
  }

  if (z.isZero()) return BigDecimal.NAN; // pole

  const p = BigDecimal.precision;
  const guard = 10;
  // Shift to w ≈ p (not the old 0.37·p, which left the asymptotic series' floor
  // at the target tolerance, forcing it to run its full ≈π·w terms). The larger
  // shift converges the series in ≈0.4·p terms; see the note in `gammalnCore`.
  const shift = Math.max(7, Math.ceil(p));

  // Recurrence: ψ(z+1) = ψ(z) + 1/z — shift z up until z > shift
  let result = new BigDecimal(0);
  let w = z;
  while (w.lt(shift)) {
    result = result.sub(BigDecimal.ONE.div(w));
    w = w.add(BigDecimal.ONE);
  }

  const maxTerms = Math.max(20, Math.ceil(0.6 * p) + 20);
  const bernoulli = getBernoulliRationals(ce, maxTerms);

  // Asymptotic expansion: ψ(w) ~ ln(w) - 1/(2w) - Σ B_{2k}/(2k·w^{2k}).
  // Round the running power w^{2k} each step: `mul` keeps the full product, so an