Special Functions: Evaluating Gamma & Error Functions

Overview

The gamma function family ($\Gamma$, $\ln\Gamma$, $B$, $\psi$, incomplete gamma and beta) and the error function family ($\operatorname{erf}$, $\operatorname{erfc}$, $\operatorname{erfcx}$, $\operatorname{erfi}$, Fresnel, Dawson) appear everywhere in probability distributions, statistics, physics, and signal processing. None of them can be written in closed form with elementary functions, and the key to their evaluation is switching to the fastest-converging method per argument region.

sangi's special functions implement this "region splitting + method switching" with a consistent design:

• The reflection formula reduces the left half-plane to an evaluation in the right half-plane ($\Gamma$, $\ln\Gamma$, $\psi$).
• Argument shifting lifts a small argument to a large one where the asymptotic expansion is effective ($\psi$, $\psi^{(n)}$, and the arbitrary-precision $\ln\Gamma$).
• Switching between series and continued fraction uses only the side that converges fast (incomplete gamma and beta, $\operatorname{erfcx}$).
• Scaling avoids overflow/underflow ($\ln\Gamma$, $\operatorname{erfcx}$).

For the native floating-point types (float / double / long double), functions present in the standard library are delegated to std::tgamma, std::lgamma, std::beta, std::erf, and std::erfc, while functions absent from the standard (digamma, polygamma, incomplete gamma/beta, $\operatorname{erfcx}$, $\operatorname{erfi}$, $\operatorname{erfInv}$, Fresnel, Dawson) are implemented in this library. The Complex type and the arbitrary-precision Float type have dedicated overloads, switching between the Lanczos approximation (complex, double precision) and the Stirling asymptotic expansion (arbitrary precision).

Related API: SpecialFunctions.

Gamma Function (Lanczos and the Reflection Formula)

The gamma function is the continuous extension of the factorial; it satisfies $\Gamma(n+1) = n!$ and the functional equation $\Gamma(z+1) = z\,\Gamma(z)$:

$$\Gamma(z) = \int_0^\infty t^{z-1} e^{-t}\,dt \qquad (\operatorname{Re}(z) > 0)$$

Lanczos approximation

The complex, double-precision gamma(Complex<R>) uses the Lanczos approximation ($g=7$, 9 coefficients). For $\operatorname{Re}(z) \geq 1/2$ it evaluates the following form:

$$\Gamma(z) = \sqrt{2\pi}\;t^{\,z-\frac12}\,e^{-t}\,A(z),\qquad t = z - 1 + g + \tfrac12$$

Here $A(z)$ is a rational function which, with the constant coefficients $c_0,\dots,c_8$, can be written as

$$A(z) = c_0 + \sum_{k=1}^{8} \frac{c_k}{(z-1)+k}$$

The implementation's coefficients are the standard Lanczos set with $g = 7$, $c_0 = 0.99999999999980993$, $c_1 = 676.5203681218851$, …, giving about 15 digits of accuracy in double precision.

Reflection formula (left half-plane)

For $\operatorname{Re}(z) < 1/2$ Lanczos is not effective, so the reflection formula reduces to the right half-plane:

$$\Gamma(z)\,\Gamma(1-z) = \frac{\pi}{\sin(\pi z)} \quad\Longrightarrow\quad \Gamma(z) = \frac{\pi}{\sin(\pi z)\,\Gamma(1-z)}$$

Since $1-z$ lies in the right half-plane, Lanczos is applied there. The zeros of $\sin(\pi z)$ (i.e., integer $z$) correspond to the poles of $\Gamma$.

Log-gamma $\ln\Gamma$ (avoiding overflow)

$\Gamma(z)$ grows rapidly and exceeds the double-precision range for only moderately large $z$. In computations that take ratios or products of $\Gamma$, such as the beta function $B(a,b)$ or the prefactor of the incomplete gamma, one first evaluates the log-gamma $\ln\Gamma$ and exponentiates at the end to avoid overflow:

$$\ln\Gamma(z) = \left(z-\tfrac12\right)\ln t - t + \tfrac12\ln(2\pi) + \ln A(z)$$

In the implementation, products and quotients are replaced by additions and subtractions, as in $B(a,b) = \exp\!\big(\ln\Gamma(a) + \ln\Gamma(b) - \ln\Gamma(a+b)\big)$. In the left half-plane it uses the logarithmic reflection formula $\ln\Gamma(z) = \ln\pi - \ln\sin(\pi z) - \ln\Gamma(1-z)$.

Stirling asymptotic expansion (arbitrary-precision version)

The arbitrary-precision Float overload lnGamma(Complex<Float>, precision) uses the Stirling asymptotic expansion rather than the fixed-coefficient Lanczos:

$$\ln\Gamma(z) = \left(z-\tfrac12\right)\ln z - z + \tfrac12\ln(2\pi) + \sum_{k=1}^{K} \frac{B_{2k}}{2k\,(2k-1)\,z^{2k-1}}$$

The $B_{2k}$ are the Bernoulli numbers, generated by the Akiyama-Tanigawa method according to the requested precision. The asymptotic expansion converges faster the larger $|z|$ is, so when $|z|$ is small the argument shift $\Gamma(z+1) = z\,\Gamma(z)$ is repeated to make $z$ large before expanding. Because the asymptotic series turns into a divergent one past its best term, it must be truncated once the terms begin to grow (optimal truncation). The shift amount and the number of Bernoulli terms are tuned to the requested number of digits. $\Gamma$ itself is obtained as $\exp(\ln\Gamma)$.

Digamma and Polygamma

The digamma function is the derivative of the log-gamma, and the polygamma functions are its higher-order derivatives:

$$\psi(z) = \frac{d}{dz}\ln\Gamma(z) = \frac{\Gamma'(z)}{\Gamma(z)},\qquad \psi^{(n)}(z) = \frac{d^{\,n+1}}{dz^{\,n+1}}\ln\Gamma(z)$$

$\psi_1 = \psi^{(1)}$ is called the trigamma function. These are not in the standard library, so they are implemented here.

Argument shift + asymptotic expansion

The asymptotic expansion of $\psi$ is more accurate the larger $z$ is:

$$\psi(z) \approx \ln z - \frac{1}{2z} - \sum_{k=1}^{K} \frac{B_{2k}}{2k\,z^{2k}}$$

So the recurrence (the logarithmic derivative of the functional equation)

$$\psi(z+1) = \psi(z) + \frac{1}{z}$$

is used to lift the argument to a size where the asymptotic expansion is effective. The native version lifts to $z \geq 8$ and then evaluates with 7 Bernoulli coefficients. The trigamma is likewise shifted by $\psi_1(z) = \psi_1(z+1) + 1/z^2$ before evaluating

$$\psi_1(z) \approx \frac{1}{z} + \frac{1}{2z^2} + \sum_{k=1}^{K} \frac{B_{2k}}{z^{2k+1}}$$

The general $\psi^{(n)}$ ($n \geq 2$) is shifted to $z \geq 10$ and evaluated with a general asymptotic expansion whose coefficients include the sign $(-1)^{n+1}$ and the factorial factor $\prod_{j=1}^{n-1}(2k+j)$ (truncating when a term falls below the relative precision or when the terms start to grow).

Reflection formula (negative arguments)

The left half-plane of $\psi$ is reduced to the right half-plane by the reflection formula:

$$\psi(1-z) - \psi(z) = \pi\cot(\pi z),\qquad \psi_1(1-z) + \psi_1(z) = \frac{\pi^2}{\sin^2(\pi z)}$$

When $z$ is a non-positive integer it is a pole, so NaN is returned. The arbitrary-precision Float version generates Bernoulli numbers to the requested number of digits and, like $\ln\Gamma$, controls the asymptotic series with optimal truncation.

Incomplete Gamma and Incomplete Beta

The regularized incomplete gamma functions are the gamma integral cut off partway through, divided by $\Gamma(a)$ to normalize it onto $[0,1]$:

$$P(a,x) = \frac{\gamma(a,x)}{\Gamma(a)} = \frac{1}{\Gamma(a)}\int_0^x t^{a-1}e^{-t}\,dt,\qquad Q(a,x) = \frac{\Gamma(a,x)}{\Gamma(a)} = 1 - P(a,x)$$

These connect directly to the cumulative distributions of the chi-squared, gamma, and Poisson distributions.

Switching between series and continued fraction

The lower $P(a,x)$ is suited to a Taylor series, and the upper $Q(a,x)$ to a continued fraction. sangi switches between the two at the boundary $x = a+1$:

• $x < a+1$: evaluate $P$ by the Taylor series $$P(a,x) = \frac{x^a e^{-x}}{\Gamma(a)}\sum_{n=0}^{\infty}\frac{x^n}{a(a+1)\cdots(a+n)}$$ and set $Q = 1 - P$.
• $x \geq a+1$: evaluate $Q$ by Legendre's continued fraction $$Q(a,x) = \frac{x^a e^{-x}}{\Gamma(a)}\cdot \cfrac{1}{x+1-a-\cfrac{1\cdot(1-a)}{x+3-a-\cfrac{2\cdot(2-a)}{x+5-a-\cdots}}}$$ and set $P = 1 - Q$.

The prefactor $x^a e^{-x}/\Gamma(a)$ is computed via logarithms as $\exp(-x + a\ln x - \ln\Gamma(a))$ to avoid overflow. The continued fraction is evaluated with the modified Lentz method: since it breaks down as the denominator approaches 0, it clamps with a small quantity $\texttt{tiny}$ while advancing $C_n = b_n + a_n/C_{n-1}$, $D_n = 1/(b_n + a_n D_{n-1})$, $f_n = f_{n-1}\,C_n D_n$, and stops once $|C_n D_n - 1|$ falls below machine precision. The non-regularized versions $\gamma(a,x) = P(a,x)\,\Gamma(a)$ and $\Gamma(a,x) = Q(a,x)\,\Gamma(a)$ are also provided.

Regularized incomplete beta $I_x(a,b)$

The incomplete beta appears in the cumulative distributions of the beta, Student's $t$, and $F$ distributions:

$$I_x(a,b) = \frac{1}{B(a,b)}\int_0^x t^{a-1}(1-t)^{b-1}\,dt$$

This too is evaluated by a continued fraction. Since the continued fraction converges fast on the small-$x$ side, the symmetry relation $I_x(a,b) = 1 - I_{1-x}(b,a)$ is used: when $x \geq (a+1)/(a+b+2)$, swap $x \to 1-x$, $a \leftrightarrow b$ to bring it onto the fast-converging side. The prefactor is computed via logarithms as $\exp(a\ln x + b\ln(1-x) + \ln\Gamma(a+b) - \ln\Gamma(a) - \ln\Gamma(b))/a$, and the continued fraction body is evaluated with the same modified Lentz method as the incomplete gamma (the standard expansion whose coefficient $d_n$ takes a different form for odd and even terms).

Error Functions (erf / erfc / erfcx / erfi / erfInv)

The error function and its relatives appear in the normal distribution, diffusion, the Voigt function, and more:

$$\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt,\qquad \operatorname{erfc}(x) = 1 - \operatorname{erf}(x)$$

The native $\operatorname{erf}$ and $\operatorname{erfc}$ are delegated to the standard library (std::erf, std::erfc). The functions absent from the standard, $\operatorname{erfcx}$, $\operatorname{erfi}$, and $\operatorname{erfInv}$, together with the complex and arbitrary-precision overloads, are implemented here.

Scaled complementary error function $\operatorname{erfcx}$

For large $x$, $\operatorname{erfc}(x)$ underflows exponentially and $e^{x^2}$ overflows, so their product is evaluated stably as a single function:

$$\operatorname{erfcx}(x) = e^{x^2}\operatorname{erfc}(x) \sim \frac{1}{x\sqrt{\pi}} \sum_{n=0}^{\infty}\frac{(-1)^n (2n-1)!!}{(2x^2)^n} \qquad (x \to +\infty)$$

The native version switches method by region:

• $x > 4$: the asymptotic expansion above (optimal truncation once the terms start to grow).
• $0 < x \leq 4$: since $e^{x^2}$ does not yet overflow, the direct $e^{x^2}\operatorname{erfc}(x)$.
• $x < 0$: the reflection $\operatorname{erfcx}(x) = 2e^{x^2} - \operatorname{erfcx}(-x)$.

The complex version computes directly for $|z| \le$ a threshold, and for $\operatorname{Re}(z) \ge 0$ with large $|z|$ it evaluates the Laplace continued fraction ($a_n = n/2$) with the modified Lentz method and sets $\operatorname{erfcx}(z) = 1/(\sqrt{\pi}\,g)$. For $\operatorname{Re}(z) < 0$ it uses the same reflection formula.

Imaginary error function $\operatorname{erfi}$

The imaginary error function is the odd function defined by

$$\operatorname{erfi}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{t^2}\,dt = -i\,\operatorname{erf}(ix)$$

Since $\operatorname{erf}(ix)$ is purely imaginary, $z = ix$ is passed to the complex $\operatorname{erf}$ and its imaginary part is returned. The complex $\operatorname{erf}$ uses a series (derived from $_1F_1$) that sums only positive terms to avoid cancellation, and converges over the whole plane. Because it diverges on the order of $e^{x^2}$, in double precision it clamps to $\pm\infty$ for $|x| \gtrsim 27$.

Inverse error function $\operatorname{erfInv}$

It finds the $y$ satisfying $\operatorname{erf}(y) = p$ by Newton iteration. Since the derivative is $\dfrac{d}{dy}\operatorname{erf}(y) = \dfrac{2}{\sqrt{\pi}}e^{-y^2}$:

$$y_{n+1} = y_n - \frac{\operatorname{erf}(y_n) - p}{\frac{2}{\sqrt{\pi}}e^{-y_n^2}}$$

The initial guess is split by region:

• $|p| \leq 0.7$: the first-order approximation $y_0 = \frac{\sqrt{\pi}}{2}\,p$.
• $|p| > 0.7$: the tail approximation $y_0 = \operatorname{sign}(p)\sqrt{-\ln\big((1-|p|)(1+|p|)\big)}$.

It converges in a few iterations. It returns $\pm\infty$ at $p = \pm 1$ and NaN for $|p| > 1$.

Fresnel Integrals and the Dawson Function

Fresnel integrals $C(x), S(x)$

The Fresnel integrals appear in diffraction, optics, and the clothoid curve (with the same normalization as the DLMF):

$$C(x) = \int_0^x \cos\!\Big(\frac{\pi t^2}{2}\Big)\,dt,\qquad S(x) = \int_0^x \sin\!\Big(\frac{\pi t^2}{2}\Big)\,dt$$

sangi switches to a Taylor series for $|x| \leq 4$ and to an asymptotic expansion for $|x| > 4$. The series is

$$C(x) = \sum_{k=0}^{\infty}\frac{(-1)^k (\pi/2)^{2k}\,x^{4k+1}}{(2k)!\,(4k+1)},\qquad S(x) = \sum_{k=0}^{\infty}\frac{(-1)^k (\pi/2)^{2k+1}\,x^{4k+3}}{(2k+1)!\,(4k+3)}$$

where $C$ and $S$ are built up together recursively. For $|x| > 4$ the asymptotic form using the auxiliary functions $f(x), g(x)$

$$C(x) = \tfrac12 + f(x)\sin\!\Big(\frac{\pi x^2}{2}\Big) - g(x)\cos\!\Big(\frac{\pi x^2}{2}\Big),\qquad S(x) = \tfrac12 - f(x)\cos\!\Big(\frac{\pi x^2}{2}\Big) - g(x)\sin\!\Big(\frac{\pi x^2}{2}\Big)$$

is evaluated with optimal truncation at the smallest term (the term ratios of $f$ and $g$ are $-(4k+1)(4k+3)/(\pi x^2)^2$ and $-(4k+3)(4k+5)/(\pi x^2)^2$, respectively). Since $C, S$ are odd functions, $x<0$ is handled by sign inversion. As $x \to \pm\infty$, $C, S \to \pm 1/2$.

Dawson function $D(x)$

The Dawson function is the odd function defined by

$$D(x) = e^{-x^2}\int_0^x e^{t^2}\,dt$$

and is related to $\operatorname{erfi}$ by $D(x) = \frac{\sqrt{\pi}}{2}e^{-x^2}\operatorname{erfi}(x)$. Unlike $\operatorname{erfi}$ or $e^{x^2}$, it is bounded (maximum $\approx 0.541$), so it is well suited to tail evaluation. sangi switches by region following Rybicki's algorithm:

• $|x| < 0.2$: the Taylor series $D(x) = \sum_{n=0}^{\infty}\dfrac{(-1)^n 2^n x^{2n+1}}{(2n+1)!!}$.
• $|x| > 5$: the asymptotic expansion $D(x) = \dfrac{1}{2x}\sum_{n=0}^{\infty}\dfrac{(2n-1)!!}{(2x^2)^n}$.
• Intermediate: direct numerical integration of $\displaystyle\int_0^x e^{t^2-x^2}\,dt$ by Simpson's rule.

Negative arguments are handled by the oddness $D(-x) = -D(x)$.

Comparison Table

The shared design principles boil down to four points: (1) reduce the left half-plane to the right with the reflection formula, (2) reduce small arguments to large ones with argument shifting, (3) switch at a boundary between series (small arguments) and continued fractions / asymptotic expansions (large arguments), and (4) avoid overflow via logarithms and scaling. Because asymptotic series are divergent series, they always come with optimal truncation at the best term.

References

• Abramowitz, M. & Stegun, I. A. (eds.) (1972). Handbook of Mathematical Functions. Dover. (§6 gamma function, §7 error function, 6.3.18 digamma asymptotic expansion)
• NIST Digital Library of Mathematical Functions (DLMF). https://dlmf.nist.gov/. (§5 gamma function, §7 error / Fresnel / Dawson, §8 incomplete gamma / beta)
• Lanczos, C. (1964). "A precision approximation of the gamma function". Journal of the SIAM, Series B: Numerical Analysis, 1(1), 86–96.
• Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes. 3rd ed. Cambridge University Press. (series/continued fractions for the incomplete gamma and beta, the modified Lentz method, Rybicki's method for Dawson)
• Cody, W. J. (1969). "Rational Chebyshev approximation for the error function". Mathematics of Computation, 23(107), 631–637.
• Rybicki, G. B. (1989). "Dawson's integral and the sampling theorem". Computers in Physics, 3(2), 85–87.