How do you evaluate the Riemann zeta function ζ(s) for Re(s) < 0?

For Re(s) ≥ 0, ζ(s) is obtained through η(s) = (1-2^{1-s})ζ(s), with η(s) computed by alternating-series acceleration (the Hasse series / Euler transformation). For Re(s) 1, where the series evaluation is performed. At negative even s, sin(πs/2) = 0 yields the trivial zeros ζ(-2n) = 0 directly.

Which algorithm computes the Hurwitz zeta function ζ(s,a)?

Euler-Maclaurin summation. The first N terms are summed directly, then the integral term (a+N)^{1-s}/(s-1), the midpoint correction (1/2)(a+N)^{-s}, and Bernoulli-number correction terms B_{2j} are added. An argument shift N makes a+N large enough that the correction series converges. The double version fixes N=15 and M=8; the arbitrary-precision Float version adjusts N and the number of correction terms to the precision and truncates at a relative ε.

The Gauss hypergeometric ₂F₁ diverges for |z| ≥ 1; how is that handled?

The power series converges only for |z| < 1. In the region z < -0.5, the Pfaff transformation ₂F₁(a,b;c;z) = (1-z)^{-a} ₂F₁(a,c-b;c;z/(z-1)) maps the argument into the convergent region |z/(z-1)| < 1 before the series is summed. When a or b is a non-positive integer the series becomes a finite polynomial and can be evaluated for all z. Other |z| ≥ 1 values return NaN.

How is the Lambert W function obtained?

W(x)·e^{W(x)} = x is solved by Halley iteration (cubic convergence). The principal branch W₀ (x ≥ -1/e) and the second branch W₋₁ (-1/e ≤ x < 0) are handled separately, iterating from an initial guess chosen by the magnitude of x (a series near -1/e, a Fritsch approximation, or the asymptotic form ln(x)-ln(ln(x))). The complex-argument version uses the same Halley iteration, forming the initial guess in double and refining to arbitrary precision.

Special Functions: Zeta, Hypergeometric & More

Q: How do you evaluate the Riemann zeta function ζ(s) for Re(s) < 0?

For Re(s) ≥ 0, ζ(s) is obtained through η(s) = (1-2^{1-s})ζ(s), with η(s) computed by alternating-series acceleration (the Hasse series / Euler transformation). For Re(s) 1, where the series evaluation is performed. At negative even s, sin(πs/2) = 0 yields the trivial zeros ζ(-2n) = 0 directly.

Q: Which algorithm computes the Hurwitz zeta function ζ(s,a)?

Euler-Maclaurin summation. The first N terms are summed directly, then the integral term (a+N)^{1-s}/(s-1), the midpoint correction (1/2)(a+N)^{-s}, and Bernoulli-number correction terms B_{2j} are added. An argument shift N makes a+N large enough that the correction series converges. The double version fixes N=15 and M=8; the arbitrary-precision Float version adjusts N and the number of correction terms to the precision and truncates at a relative ε.

Q: The Gauss hypergeometric ₂F₁ diverges for |z| ≥ 1; how is that handled?

The power series converges only for |z| < 1. In the region z < -0.5, the Pfaff transformation ₂F₁(a,b;c;z) = (1-z)^{-a} ₂F₁(a,c-b;c;z/(z-1)) maps the argument into the convergent region |z/(z-1)| < 1 before the series is summed. When a or b is a non-positive integer the series becomes a finite polynomial and can be evaluated for all z. Other |z| ≥ 1 values return NaN.

Q: How is the Lambert W function obtained?

W(x)·e^{W(x)} = x is solved by Halley iteration (cubic convergence). The principal branch W₀ (x ≥ -1/e) and the second branch W₋₁ (-1/e ≤ x < 0) are handled separately, iterating from an initial guess chosen by the magnitude of x (a series near -1/e, a Fritsch approximation, or the asymptotic form ln(x)-ln(ln(x))). The complex-argument version uses the same Halley iteration, forming the initial guess in double and refining to arbitrary precision.

Overview

Three tools recur throughout special-function evaluation: series summation, asymptotic expansion, and analytic continuation. sangi's special namespace evaluates each function by choosing among these as appropriate.

At the center stand the zeta family and the hypergeometric family, which together form a framework that expresses many special functions in a unified way. The Lerch transcendent $\Phi(z,s,a) = \sum_{k=0}^{\infty} z^k/(k+a)^s$ contains them all: the Hurwitz zeta $\zeta(s,a)$ at $z=1$, the alternating sum (Dirichlet $\beta$) at $z=-1$, and the polylogarithm $\mathrm{Li}_s(z)$ when $z$ is kept free. On the other side, the generalized hypergeometric $_pF_q$ bundles elementary functions, Bessel functions, and confluent functions together as special cases of one power series.

This page works through the implemented evaluation algorithms in order, taking the zeta and hypergeometric families as the axis and proceeding to the exponential, logarithmic, and trigonometric integrals derived from them, and then on to Lambert W, Owen's T, Coulomb wave functions, the Fermi-Dirac/Debye/transport integrals, and the angular-momentum coupling coefficients (Wigner symbols).

Related API: SpecialFunctions. The supported types are generally float / double / long double (native floating point), with some functions also supporting the arbitrary-precision Float and the complex Complex<double> / Complex<Float>.

The Zeta Family

Functions provided: Riemann $\zeta(s)$, Hurwitz $\zeta(s,a)$, Dirichlet $\eta(s)$, the Clausen functions $\mathrm{Cl}_n(\theta)$, and the Lerch transcendent $\Phi(z,s,a)$.

Riemann ζ via the Dirichlet η

The naive defining series $\sum_{n\ge 1} n^{-s}$ for $\zeta(s)$ diverges for $\mathrm{Re}(s) \le 1$ and converges slowly even for $\mathrm{Re}(s) > 1$. We therefore route through the alternating zeta (Dirichlet $\eta$):

$$\eta(s) = \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{k^s} = \bigl(1 - 2^{1-s}\bigr)\,\zeta(s), \qquad \zeta(s) = \frac{\eta(s)}{1 - 2^{1-s}}.$$

Because $\eta(s)$ is an alternating series, acceleration is effective. For complex native types we use the Hasse series (Euler transformation):

$$\eta(s) = \sum_{n=0}^{N} \frac{1}{2^{n+1}} \sum_{k=0}^{n} (-1)^k \binom{n}{k}\,(k+1)^{-s}.$$

The error shrinks by about one bit per term ($2^{-N}$), so for double-equivalent precision $N=50$ delivers roughly 15 digits. The arbitrary-precision Complex<Float> version takes the number of Euler-transformation terms proportional to the bit count of the working precision and truncates at a relative $\varepsilon$. The real double version delegates to std::riemann_zeta / std::pow while returning special values such as $\eta(1)=\ln 2$ directly to avoid the poles.

Analytic continuation via the functional equation

For $\mathrm{Re}(s) < 0$, the functional equation reflects into the convergent region $\mathrm{Re}(1-s) > 1$:

$$\zeta(s) = 2^{s}\,\pi^{s-1}\,\sin\!\left(\frac{\pi s}{2}\right) \Gamma(1-s)\,\zeta(1-s).$$

The $\zeta(1-s)$ on the right lies in the convergent region and can be evaluated via $\eta$, while $\Gamma(1-s)$ calls the gamma-function implementation. At negative even $s=-2n$, $\sin(\pi s/2)=0$ gives the trivial zeros $\zeta(-2n)=0$ directly.

Hurwitz ζ and Euler-Maclaurin summation

The Hurwitz zeta $\zeta(s,a) = \sum_{k=0}^{\infty}(a+k)^{-s}$ sums the first $N$ terms directly and approximates the remainder by Euler-Maclaurin summation:

$$\zeta(s,a) \approx \sum_{k=0}^{N-1}(a+k)^{-s} + \frac{(a+N)^{1-s}}{s-1} + \frac{(a+N)^{-s}}{2} + \sum_{j=1}^{M} \frac{B_{2j}}{(2j)!}\,(s)_{2j-1}\,(a+N)^{-(s+2j-1)},$$

where $B_{2j}$ are the Bernoulli numbers and $(s)_{2j-1}$ is a rising factorial. The integral term $(a+N)^{1-s}/(s-1)$ and the midpoint correction $(a+N)^{-s}/2$ are followed by the Bernoulli correction series. Making $a+N$ large enough through the argument shift $N$ lets the correction series converge quickly.

The double version fixes the empirically optimal $N=15$ and $M=8$ ($M>8$ degrades accuracy because the Bernoulli numbers diverge). The arbitrary-precision Complex<Float> version takes $N$ so that $a+N \gtrsim 0.4\,\mathrm{wp}$ relative to the working precision, caps the number of correction terms at $M \approx \pi(a+N)$, and truncates early once a term falls below the relative $\varepsilon$. The Bernoulli numbers $B_{2k}$ are shared with the gamma implementation. When $a=1$ it reduces to $\zeta(s,1)=\zeta(s)$, and $s=1$ is a pole, so NaN is returned.

The Clausen functions

The Clausen function $\mathrm{Cl}_2(\theta) = \sum_{k=1}^{\infty}\sin(k\theta)/k^2 = -\int_0^\theta \ln\!\bigl|2\sin(t/2)\bigr|\,dt$ is evaluated after folding the argument into $[0,\pi]$ using periodicity and antisymmetry.

For small $\theta$: the small-angle expansion containing a logarithm, $\mathrm{Cl}_2(\theta) = \theta(1-\ln\theta) + \theta^3/36 - \theta^5/3600 + \cdots$
For general $\theta$: the Fourier series $\sum \sin(k\theta)/k^2$ summed directly with Kahan compensated summation

The generalized version $\mathrm{Cl}_n(\theta)$ sums $\sum\cos(k\theta)/k^n$ when $n$ is even and $\sum\sin(k\theta)/k^n$ when $n$ is odd, down to a relative $\varepsilon$.

The Lerch transcendent and the unified representation

The Lerch transcendent $\Phi(z,s,a) = \sum_{k=0}^{\infty} z^k/(k+a)^s$ unifies the polylogarithm, Hurwitz zeta, Riemann zeta, and Dirichlet $\eta$ in a single framework:

$$\mathrm{Li}_s(z) = z\,\Phi(z,s,1),\quad \zeta(s,a) = \Phi(1,s,a),\quad \zeta(s) = \Phi(1,s,1),\quad \eta(s) = -\Phi(-1,s,1).$$

The implementation branches on the argument values: $a^{-s}$ when $z=0$, $1/(1-z)$ when $s=0$, the Hurwitz $\zeta(s,a)$ when $z=1$ (finite only for $s>1$), the Hasse-accelerated alternating sum when $z=-1$, and the direct Taylor series when $|z|<1$. Analytic continuation for $|z|>1$ ($z \neq \pm 1$) is a future extension and currently returns NaN.

The Hypergeometric Family

Functions provided: Gauss $_2F_1(a,b;c;z)$, the confluent $_1F_1(a;b;z)$ (Kummer M), the confluent limit $_0F_1(;b;z)$, the generalized $_pF_q$, Meijer G, and Mittag-Leffler $E_{\alpha,\beta}(z)$.

The power series and its region of convergence

The generalized hypergeometric function is defined by a power series:

$$_pF_q(a_1,\ldots,a_p;\,b_1,\ldots,b_q;\,z) = \sum_{k=0}^{\infty} \frac{(a_1)_k \cdots (a_p)_k}{(b_1)_k \cdots (b_q)_k}\,\frac{z^k}{k!}.$$

Successive terms are built up by recursively multiplying the ratio $t_{k+1}/t_k = z\,\prod_i(a_i+k)\big/\bigl[(k+1)\prod_j(b_j+k)\bigr]$, so the Pochhammer symbols never need to be computed explicitly. The region of convergence is determined as follows:

$p \le q$: converges for all $z$ (entire function). $_0F_1$ and $_1F_1$ fall here
$p = q+1$: converges for $|z| < 1$, with the case $z=1$ parameter-dependent. $_2F_1$ falls here
$p > q+1$: a formal (asymptotic) series — this implementation does not handle it and returns NaN

In every case, if any $a_i$ is a non-positive integer the sum becomes a finite polynomial, which can be truncated and evaluated for all $z$. If a $b_j$ is a non-positive integer (a gamma pole) and the corresponding $a_i$ does not terminate first, it is a pole and NaN is returned.

Extending the convergence region with transformation formulas

Since $_2F_1$ converges as a series only for $|z| < 1$, the argument is mapped in poorly convergent regions:

$$\text{Pfaff transformation:}\quad {}_2F_1(a,b;c;z) = (1-z)^{-a}\,{}_2F_1\!\left(a,\,c-b;\,c;\,\frac{z}{z-1}\right).$$

When $z < -0.5$, $z/(z-1)$ moves close to the origin, so mapping first and then summing the series converges quickly. For the confluent $_1F_1$, the Kummer transformation $M(a,b,z) = e^{z}\,M(b-a,b,-z)$ is used for $z<0$ to avoid cancellation. $_0F_1$ is an entire function converging for all $z$, with relations such as $_0F_1(;1;-z^2/4)=J_0(z)$ and $_0F_1(;1/2;-z^2/4)=\cos z$.

Meijer G (partial implementation)

The Meijer G function $G^{m,n}_{p,q}$ generalizes the hypergeometric family still further. This implementation is partial, using a residue representation, and is restricted to the setting $n=0$ with $a$ empty, $m \le q$, and each $b_j$ a distinct simple pole. In this setting G is the sum of residues at the simple poles at $b_h$, given by a product of a coefficient (a product of gamma functions), $z^{b_h}$, and $_0F_{q-1}$. Resonance where the differences between $b_j$ are integers (producing logarithmic terms) and the general form with $n>0$ are future extensions and return NaN.

The Mittag-Leffler function

The Mittag-Leffler function is a fractional-order generalization of the exponential and appears in the solutions of fractional differential equations:

$$E_{\alpha,\beta}(z) = \sum_{k=0}^{\infty} \frac{z^k}{\Gamma(\alpha k + \beta)}, \qquad E_{\alpha}(z) = E_{\alpha,1}(z).$$

For small $|z|$ the Taylor series builds each term through the reciprocal gamma $1/\Gamma(\alpha k+\beta)$ (gamma poles are handled as $1/\Gamma=0$). Once $|z|$ exceeds a threshold the Taylor series becomes slow and unstable, so the method switches to the Wiman asymptotic for $z>0$ and to a decaying asymptotic with the exponential term removed for $z<0$:

$$E_{\alpha,\beta}(z) \approx \frac{1}{\alpha}\,z^{(1-\beta)/\alpha}\,e^{z^{1/\alpha}} - \sum_{k=1}^{N} \frac{z^{-k}}{\Gamma(\beta-\alpha k)}\quad(z>0).$$

Because the asymptotic series diverges, it is truncated once the term magnitudes start to increase. Special values that reduce to elementary functions are returned directly, such as $E_{1,1}(z)=e^z$ at $\alpha=1,\beta=1$ and $E_{2,1}(z)=\cosh\sqrt{z}$ at $\alpha=2,\beta=1,z\ge 0$. Only $\alpha \in (0,2]$ is supported; otherwise an exception is thrown.

Exponential, Logarithmic, and Trigonometric Integrals

Functions provided: the exponential integral $\mathrm{Ei}(x)$, the logarithmic integral $\mathrm{li}(x)$, the generalized exponential integral $E_n(x)$, the sine integral $\mathrm{Si}(x)$, the cosine integral $\mathrm{Ci}(x)$, the hyperbolic versions $\mathrm{Shi},\mathrm{Chi}$, and the dilogarithm $\mathrm{Li}_2(x)$.

Choosing between series and continued fractions

The exponential integral uses a convergent series near the origin:

$$\mathrm{Ei}(x) = \gamma + \ln|x| + \sum_{n=1}^{\infty} \frac{x^n}{n\cdot n!},$$

where $\gamma$ is the Euler-Mascheroni constant. The logarithmic integral is implemented as $\mathrm{li}(x) = \mathrm{Ei}(\ln x)$.

The generalized exponential integral $E_n(x) = \int_1^\infty e^{-xt}/t^n\,dt$ switches method by the magnitude of the argument:

$x>1$ or $n>10$: a continued fraction by the modified Lentz method, $E_n(x) = e^{-x}\,\mathrm{CF}$. Fast and stable for large $x$
otherwise: first obtain $E_1(x)$ by series, then push upward with the recurrence $n\,E_{n+1}(x) = e^{-x} - x\,E_n(x)$

Trigonometric integrals

The sine and cosine integrals are evaluated by Taylor series:

$$\mathrm{Si}(x) = \sum_{k=0}^{\infty} \frac{(-1)^k\,x^{2k+1}}{(2k+1)\,(2k+1)!}, \qquad \mathrm{Ci}(x) = \gamma + \ln|x| + \sum_{k=1}^{\infty} \frac{(-1)^k\,x^{2k}}{2k\,(2k)!}.$$

Terms are updated recursively and truncated at a relative $\varepsilon$. $\mathrm{Si}$ is an odd function with $\mathrm{Si}(\infty)=\pi/2$, and $\mathrm{Ci}$ takes real values only for $x>0$. The hyperbolic versions $\mathrm{Shi},\mathrm{Chi}$ are evaluated by the series with the alternating signs removed, switching to evaluation via $\mathrm{Ei}$ (e.g. $\mathrm{Shi}=(\mathrm{Ei}(x)-\mathrm{Ei}(-x))/2$) for $|x|>20$.

The dilogarithm

The dilogarithm $\mathrm{Li}_2(x) = \sum_{n\ge 1} x^n/n^2$ uses the direct Taylor series for $|x|\le 0.5$ and the reflection formula for $0.5 < |x| \le 1$ to ensure convergence:

$$\mathrm{Li}_2(x) = -\mathrm{Li}_2(1-x) + \frac{\pi^2}{6} - \ln x\,\ln(1-x).$$

The complex version additionally uses the inversion formula $\mathrm{Li}_2(z) = -\mathrm{Li}_2(1/z) - \pi^2/6 - (\ln(-z))^2/2$ for $|z|>1$, applying the reflection only after confirming $|1-z| < |z|$ so that it does not recurse infinitely.

Other Special Functions

Lambert W (Halley iteration, branches)

Lambert W is the function satisfying $W(x)\,e^{W(x)} = x$, solved by Halley iteration (cubic convergence):

$$w_{n+1} = w_n - \frac{w_n e^{w_n} - x}{e^{w_n}(w_n+1) - \dfrac{(w_n+2)(w_n e^{w_n} - x)}{2w_n+2}}.$$

The principal branch $W_0$ ($x \ge -1/e$, $W_0 \ge -1$) and the second branch $W_{-1}$ ($-1/e \le x < 0$, $W_{-1} \le -1$) are handled separately. To speed convergence, the initial guess is chosen by the region of the argument:

Near $-1/e$: the series $W \approx -1 + p - p^2/3 + \cdots$ using $p=\sqrt{2(ex+1)}$
Small $x$: the Fritsch approximation
Large $x$: the asymptotic form $W \approx \ln x - \ln(\ln x)$

The complex-argument $W_0,W_{-1}$ are evaluated by the same Halley iteration. The arbitrary-precision version forms the initial guess in double and refines to the target precision with Halley iteration.

Owen's T

Owen's T function $T(h,a) = \dfrac{1}{2\pi}\displaystyle\int_0^a \dfrac{e^{-h^2(1+t^2)/2}}{1+t^2}\,dt$ is used for the quadrant probabilities of the bivariate normal distribution and for the CDFs of the non-central $t$ and skew-normal distributions.

Native floating point: 10-point Gauss-Legendre quadrature computes the inner integral directly
Arbitrary-precision Float: a Cauchy-product series $\sum d_n t^{2n}$ ($d_n = a_n - d_{n-1}$) of the integrand $e^{-h^2 t^2/2}/(1+t^2)$ is expanded

For $a>1$, Owen's transformation $T(h,a) = \tfrac{1}{2}[\Phi(h)+\Phi(ah)] - \Phi(h)\Phi(ah) - T(ah,1/a)$ reduces it to the evaluation of $a \le 1$.

Coulomb wave functions

The regular Coulomb wave function $F_L(\eta,\rho)$ is evaluated through the normalization constant $C_L(\eta)$ and a power series:

$$F_L(\eta,\rho) = C_L(\eta)\,\rho^{L+1}\sum_{k\ge 0} a_k\,\rho^k,\qquad a_{k+1} = \frac{2\eta a_k - a_{k-1}}{(k+1)(k+2L+2)}.$$

The normalization constant $C_L(\eta)$ is built recursively from the Sommerfeld parameter $2\pi\eta/(e^{2\pi\eta}-1)$. The irregular function $G_L(\eta,\rho)$ is constructed from the asymptotic phase (the Coulomb phase shift) and the Wronskian relation $F_L G_L' - F_L' G_L = 1$ for moderate or larger $\rho$, and returns NaN for small $\rho$ where it is singular.

Fermi-Dirac integrals

The complete Fermi-Dirac integral $F_s(x) = \dfrac{1}{\Gamma(s+1)}\displaystyle\int_0^\infty \dfrac{t^s}{e^{t-x}+1}\,dt$ appears in the electron density in semiconductor physics.

$x>20$: the Sommerfeld expansion $F_s(x) \approx \dfrac{x^{s+1}}{(s+1)\Gamma(s+1)}\bigl[1 + \tfrac{s(s+1)\pi^2}{6x^2} + \cdots\bigr]$
otherwise: 16-point Gauss-Laguerre quadrature

Dedicated entry points are provided for $F_{1/2},F_{-1/2},F_{3/2}$.

Debye functions and transport integrals

The Debye function $D_n(x) = \dfrac{n}{x^n}\displaystyle\int_0^x \dfrac{t^n}{e^t-1}\,dt$ appears in the specific heat of solids ($D_3$) and elsewhere, and the transport integral $J_n(x) = \displaystyle\int_0^x \dfrac{t^n e^t}{(e^t-1)^2}\,dt$ appears in transport phenomena. Both compute the inner integral by 15-point Gauss-Legendre quadrature, treating the singularity of the integrand near the origin (its asymptotic to $t^{n-1}$ or $t^{n-2}$) separately to preserve integrability. The synchrotron functions $F(x),G(x)$ are routed through asymptotic and integral approximations of the modified Bessel function $K_\nu$.

Angular-momentum coupling coefficients (Wigner symbols)

The Wigner 3j/6j/9j symbols and the Clebsch-Gordan coefficients appear in the coupling of angular momenta. Half-integer angular momenta are stored with $2j$ as an integer ($\texttt{two\_j} = 2j$) to avoid rounding error.

3j: after checking the triangle conditions and selection rules, the Racah formula is evaluated. Factorials are handled with $\ln$ factorials (lgamma), and the sum is taken via $\exp(\ln\text{prefactor} - \ln\text{term})$ to avoid overflow
Clebsch-Gordan: reduced to a 3j multiplied by a phase factor and $\sqrt{2J+1}$
6j: four triangle conditions and the Racah sum formula
9j: expanded into a sum over an intermediate variable $t$ of a triple product of 6j symbols, $\sum_t (-1)^{2t}(2t+1)\,\{\cdots\}\{\cdots\}\{\cdots\}$

Comparison Table

Function	Main evaluation method	Convergence region / applicability
Riemann $\zeta(s)$	alternating-series acceleration via $\eta$ (Hasse / Euler transformation) + functional equation	all $s \neq 1$ (including $\mathrm{Re}(s)<0$ via functional equation)
Hurwitz $\zeta(s,a)$	Euler-Maclaurin summation (Bernoulli correction)	$s \neq 1$, $a > 0$
Dirichlet $\eta(s)$	Hasse series / $(1-2^{1-s})\zeta(s)$	all $s$
Clausen $\mathrm{Cl}_n(\theta)$	small-angle expansion + Fourier series (Kahan summation)	$\theta$ folded by period
Lerch $\Phi(z,s,a)$	branch on value (Hurwitz / Hasse / Taylor)	$\|z\|\le 1$ ($\|z\|>1$ future)
$_2F_1(a,b;c;z)$	power series + Pfaff transformation	$\|z\|<1$ (all $z$ for polynomials)
$_1F_1(a;b;z)$	power series + Kummer transformation	all $z$ (entire)
$_0F_1(;b;z)$	power series	all $z$ (entire)
general $_pF_q$	recursive power series	$p\le q$ all $z$ / $p=q{+}1$ for $\|z\|<1$
Meijer G (partial)	residue sum over simple poles	$n=0$, $a$ empty, distinct $b_j$ (simple poles)
Mittag-Leffler $E_{\alpha,\beta}$	Taylor + Wiman / decaying asymptotic	$\alpha\in(0,2]$, real $z$
$\mathrm{Ei},\mathrm{li},\mathrm{Si},\mathrm{Ci},\mathrm{Li}_2$	series / reflection formula	per-function domain
$E_n(x)$	continued fraction (Lentz) / series + recurrence	continued fraction for $x>1$, series for small $x$
Lambert W	Halley iteration (cubic)	$W_0: x\ge -1/e$, $W_{-1}: -1/e\le x<0$
Owen's T	10-point Gauss-Legendre / Cauchy-product series	reduced to $a\le 1$ by transformation
Coulomb $F_L,G_L$	power series + asymptotic phase + Wronskian	$F$ for $\rho\ge 0$, $G$ for moderate or larger $\rho$
Fermi-Dirac $F_s$	16-point Gauss-Laguerre / Sommerfeld expansion	asymptotic for $x>20$, quadrature otherwise
Debye $D_n$ / transport $J_n$	15-point Gauss-Legendre quadrature	$x\ge 0$, $n\ge 1$ ($J_n$ requires $n\ge 2$)
Wigner 3j/6j/9j, CG	Racah formula ($\ln$ factorials)	when triangle conditions and selection rules hold

Broadly: the zeta family relies on series acceleration + functional equation, the hypergeometric family on power series + transformation formulas, the integral-defined functions (Owen's T, Fermi-Dirac, Debye, transport) on Gauss quadrature, and the equation-defined function (Lambert W) on Halley iteration. That is the basic strategy.

References

NIST Digital Library of Mathematical Functions (DLMF). https://dlmf.nist.gov/ — §16 (generalized hypergeometric, Meijer G), §25 (zeta, polylogarithm, Lerch).
Abramowitz, M. & Stegun, I. A. (1964). Handbook of Mathematical Functions. National Bureau of Standards.
Borwein, P. (2000). "An efficient algorithm for the Riemann zeta function". Constructive, Experimental, and Nonlinear Analysis, CMS Conf. Proc. 27, 29–34.
Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J. & Knuth, D. E. (1996). "On the Lambert W function". Advances in Computational Mathematics, 5(1), 329–359.
Gorenflo, R., Kilbas, A. A., Mainardi, F. & Rogosin, S. V. (2014). Mittag-Leffler Functions, Related Topics and Applications. Springer.
Patefield, M. & Tandy, D. (2000). "Fast and accurate calculation of Owen's T function". Journal of Statistical Software, 5(5), 1–25.
Edmonds, A. R. (1957). Angular Momentum in Quantum Mechanics. Princeton University Press (Wigner symbols, Racah formula).