Function polygammaCore

src/compute-engine/numerics/special-functions.ts:957–1033

(ce: ComputeEngine, nNum: number, z: BigNum)

  result = result.add(BigDecimal.ONE.div(w));
  result = result.add(BigDecimal.ONE.div(w.mul(w).mul(2)));
  let w2kp1 = w.mul(w)._mulToPrecision(w, p); // w^3
  const w2 = w._mulToPrecision(w, p);
  const tol = new BigDecimal(10).pow(-(p + guard));
  const nTerms = Math.min(maxTerms, bernoulli.length);
  for (let k = 0; k < nTerms; k++) {
    const [bNum, bDen] = bernoulli[k];
    const term = new BigDecimal(bNum.toString()).div(
      new BigDecimal(bDen.toString()).mul(w2kp1)
    );
    if (k > 0 && term.abs().lt(tol)) break;
    result = result.add(term);
    w2kp1 = w2kp1._mulToPrecision(w2, p);
  }

  return result;
}

/**
 * Bignum Polygamma function ψₙ(z) = dⁿ/dzⁿ ψ(z)
 * Delegates to bigDigamma/bigTrigamma for n=0,1.
 * For n ≥ 2, uses recurrence + asymptotic expansion.
 */
export function bigPolygamma(ce: ComputeEngine, n: BigNum, z: BigNum): BigNum {
  const nNum = n.toNumber();
  if (!Number.isInteger(nNum) || nNum < 0) return BigDecimal.NAN;
  if (nNum === 0) return bigDigamma(ce, z); // already guarded
  if (nNum === 1) return bigTrigamma(ce, z); // already guarded
  if (!z.isFinite() || z.isZero()) return BigDecimal.NAN;
  return withGuardDigits(SPECIAL_FN_GUARD, () => polygammaCore(ce, nNum, z));
}

function polygammaCore(ce: ComputeEngine, nNum: number, z: BigNum): BigNum {
  // Bignum factorial helper (small n, simple loop)
  const bigFactorial = (m: number): BigNum => {
    let r: BigNum = BigDecimal.ONE;
    for (let i = 2; i <= m; i++) r = r.mul(i);
    return r;
  };

  const p = BigDecimal.precision;
  const guard = 10;
  // Shift to w ≈ p so the asymptotic series converges in ≈0.4·p terms rather
  // than running its full ≈π·w (see `gammalnCore`).
  const shift = Math.max(7, Math.ceil(p));

  // Poles at z = 0, −1, −2, …
  let w = z;
  let result = new BigDecimal(0);
  if (w.isNegative() && w.isInteger()) return BigDecimal.NAN;

  // Recurrence (from ψ(z+1) = ψ(z) + 1/z differentiated n times, DLMF 5.15.5):
  //   ψ⁽ⁿ⁾(z) = ψ⁽ⁿ⁾(z+1) + (−1)^{n+1} n!/z^{n+1}
  // The single shift loop also lifts negative non-integer z past the poles.
  const sign = nNum % 2 === 0 ? -1 : 1; // (−1)^{n+1}
  const bigFactN = bigFactorial(nNum);
  while (w.lt(shift)) {
    result = result.add(
      new BigDecimal(sign).mul(bigFactN).div(w.pow(nNum + 1))