| 1577 | const factN = factorial(n); |
| 1578 | const sign = n % 2 === 0 ? -1 : 1; // (−1)^{n+1} |
| 1579 | const zMin = n + 10; |
| 1580 | while (z < zMin) { |
| 1581 | result += (sign * factN) / Math.pow(z, n + 1); |
| 1582 | z += 1; |
| 1583 | } |
| 1584 | |
| 1585 | // Asymptotic expansion (DLMF 5.15.9, extended to general n — obtained by |
| 1586 | // differentiating DLMF 5.11.2 n times): |
| 1587 | // ψ⁽ⁿ⁾(z) ~ (−1)^{n−1} [ (n−1)!/zⁿ + n!/(2z^{n+1}) |
| 1588 | // + Σ_{k≥1} B₂ₖ·(2k+n−1)!/((2k)!·z^{2k+n}) ] |
| 1589 | const signA = n % 2 === 0 ? -1 : 1; // (−1)^{n−1} |
| 1590 | result += (signA * factorial(n - 1)) / Math.pow(z, n); |
| 1591 | result += (signA * factN) / (2 * Math.pow(z, n + 1)); |
| 1592 | |
| 1593 | // Bernoulli tail. The k-th factor fₖ = (2k+n−1)!/((2k)!·z^{2k+n}) is built |
| 1594 | // incrementally: f₁ = (n+1)!/(2!·z^{n+2}), |
| 1595 | // f_{k+1} = fₖ·(2k+n)(2k+n+1)/((2k+1)(2k+2)·z²). |
| 1596 | // (A previous version used n(n+1)⋯(n+2k−1)/(2k)! — a factor (n−1)! too |
| 1597 | // small; only n ≤ 2 was unaffected since 1! = 0! = 1.) |
| 1598 | let f = (factN * (n + 1)) / (2 * Math.pow(z, n + 2)); |
| 1599 | let prevAbs = Infinity; |
| 1600 | for (let k = 1; k <= BERNOULLI_2K.length; k++) { |
| 1601 | const term = BERNOULLI_2K[k - 1] * f; |
| 1602 | // The expansion is divergent: stop at the smallest term |
| 1603 | if (Math.abs(term) >= prevAbs) break; |
| 1604 | result += signA * term; |
| 1605 | prevAbs = Math.abs(term); |
| 1606 | const m = 2 * k; |
| 1607 | f *= ((m + n) * (m + n + 1)) / ((m + 1) * (m + 2) * z * z); |
| 1608 | } |
| 1609 | |
| 1610 | return result; |
| 1611 | } |
| 1612 | |
| 1613 | function factorial(n: number): number { |
| 1614 | if (n <= 1) return 1; |
| 1615 | let r = 1; |
| 1616 | for (let i = 2; i <= n; i++) r *= i; |
| 1617 | return r; |
| 1618 | } |
| 1619 | |
| 1620 | /** |
| 1621 | * Beta function B(a, b) = Γ(a)Γ(b)/Γ(a+b) |
| 1622 | * Uses gamma directly for small args (more accurate) and gammaln for large. |
| 1623 | */ |
| 1624 | export function beta(a: number, b: number): number { |
| 1625 | // For large arguments, use gammaln to avoid overflow |
| 1626 | if (a > 100 || b > 100 || a + b > 100) { |
| 1627 | return Math.exp(gammaln(a) + gammaln(b) - gammaln(a + b)); |
| 1628 | } |
| 1629 | // Direct computation: more accurate for small arguments |
| 1630 | return (gamma(a) * gamma(b)) / gamma(a + b); |
| 1631 | } |
| 1632 | |
| 1633 | /** |
| 1634 | * Riemann zeta function ζ(s) = Σ_{n=1}^∞ 1/n^s |
| 1635 | * |
| 1636 | * Uses the Cohen–Rodríguez Villegas–Zagier acceleration of the alternating |