Function zetaCore

src/compute-engine/numerics/special-functions.ts:1256–1345

(ce: ComputeEngine, s: BigNum)

  if (x.isNegative() || x.gt(BigDecimal.ONE)) return BigDecimal.NAN;
  if (x.isZero()) return BigDecimal.ZERO;
  if (x.eq(BigDecimal.ONE)) return BigDecimal.ONE;
  return withGuardDigits(SPECIAL_FN_GUARD, () => {
    const bt = gammalnCore(ce, a.add(b))
      .sub(gammalnCore(ce, a))
      .sub(gammalnCore(ce, b))
      .add(a.mul(x.ln()))
      .add(b.mul(BigDecimal.ONE.sub(x).ln()))
      .exp();
    const boundary = a.add(BigDecimal.ONE).div(a.add(b).add(BigDecimal.TWO));
    if (x.lt(boundary))
      return bt.mul(bigBetaContinuedFraction(ce, a, b, x)).div(a);
    return BigDecimal.ONE.sub(
      bt.mul(bigBetaContinuedFraction(ce, b, a, BigDecimal.ONE.sub(x))).div(b)
    );
  });
}

/**
 * Bignum Riemann zeta function ζ(s).
 *
 * The general case uses the Cohen–Villegas–Zagier acceleration of the
 * alternating Dirichlet eta series (error ~(3+√8)^{−n}, so ~1.3 digits of
 * accuracy per term); the kernel runs with working-precision guard digits so
 * the result is accurate to the full requested precision.
 */
export function bigZeta(ce: ComputeEngine, s: BigNum): BigNum {
  if (!s.isFinite()) return BigDecimal.NAN;
  if (s.eq(1)) return new BigDecimal(Infinity); // pole
  return withGuardDigits(SPECIAL_FN_GUARD, () => zetaCore(ce, s));
}

function zetaCore(ce: ComputeEngine, s: BigNum): BigNum {
  const pi = BigDecimal.PI;

  // Special value: ζ(0) = -1/2
  if (s.isZero()) return BigDecimal.HALF.neg();

  // Special values for positive even integers: ζ(2k) = (-1)^{k+1} B_{2k} (2π)^{2k} / (2(2k)!)
  // Capped at MAX_EXACT_ZETA_EVEN_N: the Bernoulli-number computation is
  // O(k²)-ish in the Bernoulli index k = s/2, so this closed form is only
  // worth it for modest s — beyond the cap the general series below (which
  // is precision-scaled and magnitude-independent) computes the same value.
  if (s.isInteger() && s.isPositive()) {
    const sn = s.toNumber();
    if (sn % 2 === 0 && sn >= 2 && sn <= MAX_EXACT_ZETA_EVEN_N) {
      const k = sn / 2;
      const bernoulli = getBernoulliRationals(ce, k);
      const [bNum, bDen] = bernoulli[k - 1];
      const bernAbs = new BigDecimal(absBigint(bNum).toString()).div(
        new BigDecimal(bDen.toString())
      );
      const twoPi = pi.mul(2);
      let factVal: BigNum = BigDecimal.ONE;
      let steps = 0;
      for (let i = 2; i <= sn; i++) {
        if ((++steps & 0xfff) === 0) checkDeadline(ce._deadlineFrame);
        factVal = factVal.mul(i);
      }