Special Functions

Overview

The sangi special-functions module provides roughly 30 families — over 80 special functions that recur throughout mathematical physics, statistics, and numerical analysis — under a single, unified API. All are free functions in namespace sangi::special.

Three families of scalar type are supported.

• Native floating point — float / double / long double. The template<IsNativeFloat T> form is the base implementation for most functions
• Complex — Complex<double> (and, for some, Complex<Float>). Gamma, error function, Bessel, zeta, elliptic, hypergeometric, and more accept complex arguments
• Arbitrary precision — Float. The multi-precision overloads take a precision argument int precision and are exposed either directly under the parent namespace sangi:: as in sangi::gamma(x, precision), or in-header as Complex<Float> overloads

Every function is marked [[nodiscard]] and returns NaN or ±∞ outside its domain, at poles, or on divergence (see each table for the per-function convention).

Build

// Aggregate header: includes the major families together
#include <math/special/special.hpp>

Header-only. No library linkage required. All functions belong to namespace sangi::special.

special.hpp pulls in the following families: gamma (gamma.hpp), factorial utilities (factorial_utils.hpp), error function (error_function.hpp), Dawson (dawson.hpp), Bessel (bessel.hpp), Airy (airy.hpp), Anger-Weber (anger_weber.hpp), Legendre (legendre.hpp), elliptic integrals (elliptic.hpp), theta (theta.hpp), Weierstrass (weierstrass.hpp), zeta (zeta.hpp), Clausen (clausen.hpp), exponential integral (exponential_integral.hpp), hypergeometric (hypergeometric.hpp), Lambert W (lambert_w.hpp), Owen's T (owens_t.hpp), Coulomb (coulomb.hpp), Fermi-Dirac (fermi_dirac.hpp), Debye (debye.hpp), transport (transport.hpp), and angular-momentum coupling (coupling.hpp).

The following families are not part of special.hpp and must be included individually.

#include <math/special/orthogonal_classical.hpp>  // Hermite, Laguerre, Chebyshev
#include <math/special/fresnel.hpp>               // Fresnel integrals C(x), S(x)
#include <math/special/kelvin.hpp>                // Kelvin functions ber, bei, ker, kei
#include <math/special/struve.hpp>                // Struve functions H_v, L_v
#include <math/special/hypergeometric_pq.hpp>     // generalized hypergeometric pFq, Meijer G
#include <math/special/mittag_leffler.hpp>        // Mittag-Leffler function
#include <math/special/lerch.hpp>                 // Lerch transcendent Phi
#include <math/special/spheroidal.hpp>            // spheroidal wave functions

Gamma Function Family

Provides the gamma function $\Gamma(z)$ together with its logarithm and derivatives (digamma $\psi$, polygamma $\psi^{(n)}$), the beta function $B(a,b)$, and the regularized incomplete gamma and beta functions. gamma.hpp handles these analytic functions, while factorial_utils.hpp covers the combinatorial factorial-type utilities. In addition to native floating-point versions, many functions have Complex<double> and Complex<Float> (precision-argument) overloads.

Function reference

gamma

// native floating point (delegates to std::tgamma)
template<IsNativeFloat T>
[[nodiscard]] T gamma(T x);

// complex argument (Lanczos approximation g=7, n=9)
template<IsNativeFloat R>
[[nodiscard]] Complex<R> gamma(Complex<R> z);

// arbitrary-precision complex (any precision, Stirling asymptotic expansion)
[[nodiscard]] Complex<Float> gamma(Complex<Float> z, int precision);
[[nodiscard]] Complex<Float> gamma(Complex<Float> z);  // default precision

Definition: $\Gamma(z) = \int_0^\infty t^{z-1} e^{-t}\,dt$ ($\operatorname{Re} z > 0$), and the analytic continuation in general. For $\operatorname{Re}(z) < 0.5$ the reflection formula $\Gamma(z) = \pi / (\sin(\pi z)\,\Gamma(1-z))$ is used. Poles at the non-positive integers $0, -1, -2, \dots$. The arbitrary-precision Complex<Float> version takes a precision argument precision (number of decimal digits).

lnGamma

template<IsNativeFloat T>
[[nodiscard]] T lnGamma(T x);                              // delegates to std::lgamma

template<IsNativeFloat R>
[[nodiscard]] Complex<R> lnGamma(Complex<R> z);            // Lanczos + log

[[nodiscard]] Complex<Float> lnGamma(Complex<Float> z, int precision);
[[nodiscard]] Complex<Float> lnGamma(Complex<Float> z);

Description: $\ln\Gamma(z)$. The native version returns the real part only (std::lgamma, corresponding to the principal value $\ln|\Gamma|$); the complex version returns the principal branch of the logarithm. Useful for large arguments where $\Gamma$ would overflow, such as normalizing constants of statistical distributions.

beta

template<IsNativeFloat T>
[[nodiscard]] T beta(T a, T b);

template<IsNativeFloat R>
[[nodiscard]] Complex<R> beta(Complex<R> a, Complex<R> b);

[[nodiscard]] Complex<Float> beta(Complex<Float> a, Complex<Float> b, int precision);
[[nodiscard]] Complex<Float> beta(Complex<Float> a, Complex<Float> b);

Definition: $B(a,b) = \dfrac{\Gamma(a)\,\Gamma(b)}{\Gamma(a+b)} = \displaystyle\int_0^1 t^{a-1}(1-t)^{b-1}\,dt$. The complex version computes $\exp(\ln\Gamma(a) + \ln\Gamma(b) - \ln\Gamma(a+b))$ to suppress overflow.

betaRegularized

template<IsNativeFloat T>
[[nodiscard]] T betaRegularized(T x, T a, T b);            // x in [0,1], a,b > 0

template<IsNativeFloat R>
[[nodiscard]] Complex<R> betaRegularized(Complex<R> z, Complex<R> a, Complex<R> b);

[[nodiscard]] Complex<Float> betaRegularized(Complex<Float> z, Complex<Float> a,
                                             Complex<Float> b, int precision);

Definition: the regularized incomplete beta $I_x(a,b) = \dfrac{1}{B(a,b)}\displaystyle\int_0^x t^{a-1}(1-t)^{b-1}\,dt$. The first argument is the upper limit $x$, followed by the shape parameters $a, b$. $I_0=0,\ I_1=1$. Stabilized with a continued fraction (modified Lentz) and the symmetry transform $I_x(a,b) = 1 - I_{1-x}(b,a)$. Appears in the cumulative distribution functions of the binomial, $F$, and $t$ distributions.

Error Function, Fresnel, Dawson

Provides the error function $\operatorname{erf}$ and its relatives (complementary $\operatorname{erfc}$, scaled complementary $\operatorname{erfcx}$, imaginary $\operatorname{erfi}$, inverse $\operatorname{erf}^{-1}$), the Fresnel integrals $C(x), S(x)$, and the Dawson function $F(x)$. error_function.hpp has complex overloads, whereas fresnel.hpp / dawson.hpp are native floating point only.

Function reference

erf

template<IsNativeFloat T>
[[nodiscard]] T erf(T x);                                  // delegates to std::erf

template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> erf(Complex<R> z);                // series (1F1 form)

Description: the error function that appears in the standard normal distribution. The complex version uses the everywhere-convergent series $\operatorname{erf}(z) = \frac{2z}{\sqrt\pi}e^{-z^2}\sum_{n\ge0}\frac{(2z^2)^n}{1\cdot3\cdots(2n+1)}$. SpecialFunctionScalar R admits double or Float (for Complex<Float> the convergence test uses the working precision of the argument).

erfcx

template<IsNativeFloat T>
[[nodiscard]] T erfcx(T x);

template<IsNativeFloat R>
[[nodiscard]] Complex<R> erfcx(Complex<R> z);

[[nodiscard]] Complex<Float> erfcx(Complex<Float> z, int precision);
[[nodiscard]] Complex<Float> erfcx(Complex<Float> z);

Description: $\operatorname{erfcx}(x) = e^{x^2}\operatorname{erfc}(x)$. Because it preserves significant digits even for large $x$ where $\operatorname{erfc}$ becomes exponentially small, it is used for evaluating Voigt profiles and diffusion equations. For large $x$ it switches to an asymptotic expansion, for small $|z|$ to direct evaluation, and for $\operatorname{Re}(z)\ge0$ to a Laplace continued fraction.

erfInv

template<IsNativeFloat T>
[[nodiscard]] T erfInv(T p);                               // p in (-1, 1)

Description: the inverse error function. Solves $\operatorname{erf}(y) = p$ for the real $y$ by Newton's method (machine precision in 3–4 iterations). Returns $\pm\infty$ at $p = \pm1$ and $\mathrm{NaN}$ for $|p| > 1$. Used to compute quantiles of the normal distribution (probability to value).

fresnelC / fresnelS

template<IsNativeFloat T>
[[nodiscard]] T fresnelC(T x);

template<IsNativeFloat T>
[[nodiscard]] T fresnelS(T x);

Definition: $C(x) = \displaystyle\int_0^x \cos\!\left(\tfrac{\pi}{2}t^2\right)dt$, $S(x) = \displaystyle\int_0^x \sin\!\left(\tfrac{\pi}{2}t^2\right)dt$ (normalized form). Both are odd and tend to $C, S \to 1/2$ as $x\to\infty$. They appear in optical diffraction and the Cornu spiral. A Taylor series is used for small $x$ and an asymptotic expansion for large $x$.

fresnelCS

template<IsNativeFloat T>
[[nodiscard]] std::pair<T, T> fresnelCS(T x);              // { C(x), S(x) }

Description: computes $C(x)$ and $S(x)$ together and returns them as a std::pair<T,T> ($\{C, S\}$, i.e. .first = C(x), .second = S(x)). Avoids redundant computation when both are needed.

dawson

template<IsNativeFloat T>
[[nodiscard]] T dawson(T x);

Definition: the Dawson function $F(x) = e^{-x^2}\displaystyle\int_0^x e^{t^2}\,dt$. Odd, with $F(x)\approx x$ near $x=0$ and $F(x)\approx 1/(2x)$ as $x\to\infty$. Related to $\operatorname{erfi}$ by $F(x) = \frac{\sqrt\pi}{2}e^{-x^2}\operatorname{erfi}(x)$. Implemented with the Rybicki (1989) algorithm.

Exponential & Logarithmic Integrals

exponential_integral.hpp provides the exponential integral $\operatorname{Ei}$ and its generalization $E_n$, the logarithmic integral $\operatorname{li}$, the trigonometric and hyperbolic integrals ($\operatorname{Si}, \operatorname{Ci}, \operatorname{Shi}, \operatorname{Chi}$), and the dilogarithm $\operatorname{Li}_2$. expint (Ei) and dilog ($\operatorname{Li}_2$) have complex overloads.

Function reference

expint / expintN

// exponential integral Ei(x)
template<IsNativeFloat T>
[[nodiscard]] T expint(T x);

template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> expint(Complex<R> z);

// generalized exponential integral E_n(x)
template<IsNativeFloat T>
[[nodiscard]] T expintN(int n, T x);                       // n >= 0

Description: expint is the exponential integral $\operatorname{Ei}(x)$ (logarithmically divergent at $x=0$ on the real axis, principal value). expintN is the generalized exponential integral $E_n(x)$ of order $n$, related by $E_1(x) = -\operatorname{Ei}(-x)$. Here $n$ is a non-negative integer order.

dilog

template<IsNativeFloat T>
[[nodiscard]] T dilog(T x);

template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> dilog(Complex<R> z);

Definition: the dilogarithm (Spence's function) $\operatorname{Li}_2(z) = \displaystyle\sum_{k=1}^\infty \frac{z^k}{k^2} = -\!\!\int_0^z \frac{\ln(1-t)}{t}\,dt$. $\operatorname{Li}_2(1) = \pi^2/6$. The complex version is the principal branch with a branch cut along $[1,\infty)$. Recurs in loop integrals of statistical mechanics and quantum field theory.

Bessel Function Family

Collects the Bessel functions $J_\nu, Y_\nu$ (first and second kind), the modified Bessel functions $I_\nu, K_\nu$, spherical Bessel functions, Hankel $H^{(1,2)}$, Airy $\operatorname{Ai}, \operatorname{Bi}$, Kelvin ($\operatorname{ber}, \operatorname{bei}, \operatorname{ker}, \operatorname{kei}$), Struve $\mathbf{H}_\nu, \mathbf{L}_\nu$, and Anger-Weber. The Bessel core and Airy accept complex arguments (and arbitrary-precision Complex<Float>). Kelvin, Struve, and Anger-Weber are native floating point only.

Argument conventions: a real order is nu (type T), an integer order is n (type int; unsigned int for the spherical Bessel functions), the real argument is x, and the complex argument is z.

Bessel functions (bessel.hpp)

besselJ

// real order nu (real)
template<IsNativeFloat T>
[[nodiscard]] T besselJ(T nu, T x);

// convenience overload for integer order n
template<IsNativeFloat T>
[[nodiscard]] T besselJ(int n, T x);

// complex argument (integer order)
[[nodiscard]] Complex<double> besselJ(int n, const Complex<double>& z);
[[nodiscard]] Complex<Float>  besselJ(int n, const Complex<Float>& z, int precision);
[[nodiscard]] Complex<Float>  besselJ(int n, const Complex<Float>& z);

Definition: $J_\nu(x)$ is the bounded solution of the Bessel differential equation $x^2 y'' + x y' + (x^2 - \nu^2)y = 0$. The real order nu may be any real number; the integer-order version takes int n. The complex overloads take an integer order n and a complex argument z. besselY / besselI / besselK share the same overload structure (real order, integer order, complex).

Airy functions (airy.hpp)

airyAi

template<IsNativeFloat T>
[[nodiscard]] T airyAi(T x);

[[nodiscard]] Complex<double> airyAi(const Complex<double>& z);
[[nodiscard]] Complex<Float>  airyAi(const Complex<Float>& z, int precision);

Definition: solutions of the Airy equation $y'' = x y$. $\operatorname{Ai}(x)$ decays exponentially as $x\to+\infty$, while $\operatorname{Bi}(x)$ diverges. They appear near turning points in quantum mechanics (WKB connection) and in focal diffraction in optics. The arbitrary-precision complex version requires a precision argument precision (there is no default-precision overload).

Kelvin, Struve, Anger-Weber

These are all native floating point only and have no complex or arbitrary-precision overloads. Of kelvin.hpp / struve.hpp / anger_weber.hpp, the headers kelvin.hpp and struve.hpp are not part of special.hpp and must be included individually (anger_weber.hpp is part of the aggregate).

Orthogonal Polynomials

Classical orthogonal polynomials: legendre.hpp provides Legendre $P_n$ and its relatives (associated Legendre $P_n^m$, the radial part of spherical harmonics $Y_n^m$, Jacobi $P_n^{(\alpha,\beta)}$, Gegenbauer $C_n^\lambda$), and orthogonal_classical.hpp provides Hermite, Laguerre, and Chebyshev. Part of the Legendre family accepts complex arguments (and Complex<Float>). orthogonal_classical.hpp is native floating point only and, being outside special.hpp, must be included individually.

Legendre family (legendre.hpp)

legendreP

template<IsNativeFloat T>
[[nodiscard]] T legendreP(int n, T x);

[[nodiscard]] Complex<double> legendreP(int n, const Complex<double>& z);
[[nodiscard]] Complex<Float>  legendreP(int n, const Complex<Float>& z, int precision);

Definition: $P_n(x)$ is the degree-$n$ polynomial orthogonal on $[-1,1]$ with weight $1$ (the polynomial solution of Legendre's equation). The argument n is a non-negative integer degree. assocLegendreP adds an order m to give $P_n^m$, and jacobiP is the generalization with parameters alpha, beta; both have complex overloads.

Hermite, Laguerre, Chebyshev (orthogonal_classical.hpp)

All have the form template<IsNativeFloat T> [[nodiscard]] T f(int n, ...), with non-negative integer degree n. Used for Gauss quadrature nodes in numerical integration, for spectroscopy, and as the basis of spectral methods.

Elliptic Integrals & Functions

elliptic.hpp provides the Legendre-form elliptic integrals ($K, E, F, \Pi$), the Carlson symmetric forms ($R_F, R_D, R_J, R_C$), and the Jacobi elliptic functions ($\operatorname{sn}, \operatorname{cn}, \operatorname{dn}$); theta.hpp provides the Jacobi theta functions $\vartheta_i$; and weierstrass.hpp provides the Weierstrass $\wp, \zeta_W, \sigma$. The argument names used are k for the modulus, phi for the amplitude, and n for the characteristic.

Elliptic integrals & Jacobi elliptic functions (elliptic.hpp)

ellipticK

template<IsNativeFloat T>
[[nodiscard]] T ellipticK(T k);                            // modulus k, |k| < 1

[[nodiscard]] Complex<double> ellipticK(const Complex<double>& k);
[[nodiscard]] Complex<Float>  ellipticK(const Complex<Float>& k, int precision);

Definition: the complete elliptic integral of the first kind $K(k) = \displaystyle\int_0^{\pi/2}\frac{d\theta}{\sqrt{1 - k^2\sin^2\theta}}$. The argument is the modulus $k$ (not the parameter $m = k^2$). Logarithmically divergent as $k\to1$, with $K(0) = \pi/2$. Closely related to the pendulum period and the AGM (arithmetic-geometric mean).

carlsonRF

template<IsNativeFloat T>
[[nodiscard]] T carlsonRF(T x, T y, T z);

[[nodiscard]] Complex<double> carlsonRF(const Complex<double>& x,
                                        const Complex<double>& y,
                                        const Complex<double>& z);

Definition: $R_F(x,y,z) = \tfrac{1}{2}\displaystyle\int_0^\infty \frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}$. The symmetric forms are numerically more stable than the Legendre standard forms, and $K, E, F$ can all be written in Carlson form (e.g. $K(k) = R_F(0, 1-k^2, 1)$). carlsonRJ takes the characteristic p as a fourth argument. The complex overloads of the Carlson family are Complex<double> only.

Jacobi theta functions (theta.hpp)

The Jacobi theta functions $\vartheta_1, \vartheta_2, \vartheta_3, \vartheta_4$. The arguments are the variable z and the nome q (or the half-period-ratio $\tau$ version). Each has native, Float real, Complex<double>, and Complex<Float> overloads.

jacobiTheta1 (representative)

// native floating point (variable z, nome q)
template<IsNativeFloat T>
[[nodiscard]] T jacobiTheta1(T z, T q);

// arbitrary-precision Float real
[[nodiscard]] Float jacobiTheta1(const Float& z, const Float& q, int precision);
[[nodiscard]] Float jacobiTheta1(const Float& z, const Float& q);

// complex (real nome / complex nome)
[[nodiscard]] Complex<double> jacobiTheta1(const Complex<double>& z, double q);
[[nodiscard]] Complex<double> jacobiTheta1(const Complex<double>& z, const Complex<double>& q);

// arbitrary-precision complex (real nome / complex nome, with precision and default precision)
[[nodiscard]] Complex<Float> jacobiTheta1(const Complex<Float>& z, const Float& q, int precision);
[[nodiscard]] Complex<Float> jacobiTheta1(const Complex<Float>& z, const Complex<Float>& q, int precision);

Description: theta series parametrized by the nome $q = e^{i\pi\tau}$ ($|q| < 1$). $\vartheta_2, \vartheta_3, \vartheta_4$ have the same overload structure. Used for constructing elliptic functions, the Jacobi triple product, and heat-kernel representations.

Weierstrass elliptic functions (weierstrass.hpp)

The lattice is specified by the half-periods $\omega_1, \omega_2$. The double overload, for rectangular lattices, supplies the second period as a pure imaginary part omega2_imag. There are no Float / Complex<Float> versions.

Zeta Function Family

zeta.hpp provides the Riemann zeta $\zeta(s)$, the Hurwitz zeta $\zeta(s, a)$, and the Dirichlet eta $\eta(s)$; clausen.hpp provides the Clausen function $\operatorname{Cl}$; and lerch.hpp provides the Lerch transcendent $\Phi$. The three functions in zeta.hpp have native, complex, and Complex<Float> overloads. lerch.hpp is not part of special.hpp and must be included individually.

Function reference

riemannZeta

template<IsNativeFloat T>
[[nodiscard]] T riemannZeta(T s);

template<IsNativeFloat R>
[[nodiscard]] Complex<R> riemannZeta(Complex<R> s);

[[nodiscard]] Complex<Float> riemannZeta(Complex<Float> s, int precision);
[[nodiscard]] Complex<Float> riemannZeta(Complex<Float> s);

Definition: $\zeta(s) = \displaystyle\sum_{n=1}^\infty n^{-s}$ ($\operatorname{Re}s > 1$), and the analytic continuation in general. Simple pole at $s = 1$. $\zeta(2) = \pi^2/6$, and the trivial zeros are the negative even integers. The complex version can evaluate within the critical strip $0 < \operatorname{Re}s < 1$ as well.

hurwitzZeta

template<IsNativeFloat T>
[[nodiscard]] T hurwitzZeta(T s, T a);

template<IsNativeFloat R>
[[nodiscard]] Complex<R> hurwitzZeta(Complex<R> s, Complex<R> a);

[[nodiscard]] Complex<Float> hurwitzZeta(Complex<Float> s, Complex<Float> a, int precision);
[[nodiscard]] Complex<Float> hurwitzZeta(Complex<Float> s, Complex<Float> a);

Definition: $\zeta(s, a) = \displaystyle\sum_{n=0}^\infty (n+a)^{-s}$. It contains the Riemann zeta as $\zeta(s, 1) = \zeta(s)$. Related to the polygamma $\psi^{(n)}$ and the Lerch transcendent, and used in constructing Dirichlet $L$-functions.

lerchPhi

template<IsNativeFloat T>
[[nodiscard]] T lerchPhi(T z, T s, T a);

[[nodiscard]] Float lerchPhi(const Float& z, const Float& s, const Float& a, int precision);
[[nodiscard]] Float lerchPhi(const Float& z, const Float& s, const Float& a);

Definition: the Lerch transcendent $\Phi(z, s, a) = \displaystyle\sum_{n=0}^\infty \frac{z^n}{(n+a)^s}$. A unifying function that contains the polylogarithm $\operatorname{Li}_s(z) = z\,\Phi(z, s, 1)$ and the Hurwitz zeta $\zeta(s, a) = \Phi(1, s, a)$ as special cases. The arbitrary-precision Float version takes a precision argument.

Hypergeometric Functions

hypergeometric.hpp provides the Gauss hypergeometric ${}_2F_1$, the confluent (Kummer) ${}_1F_1$, and ${}_0F_1$; hypergeometric_pq.hpp provides the generalized hypergeometric ${}_pF_q$ and the Meijer G function; and mittag_leffler.hpp provides the Mittag-Leffler function $E_{\alpha,\beta}$. hypergeometric_pq.hpp and mittag_leffler.hpp are not part of special.hpp and must be included individually.

Function reference

hyperg / confHyperg / hyperg0F1

// 2F1(a, b; c; z)
template<IsNativeFloat T>
[[nodiscard]] T hyperg(T a, T b, T c, T z);
template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> hyperg(Complex<R> a, Complex<R> b, Complex<R> c, Complex<R> z);

// 1F1(a; b; z) — Kummer confluent
template<IsNativeFloat T>
[[nodiscard]] T confHyperg(T a, T b, T z);
template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> confHyperg(Complex<R> a, Complex<R> b, Complex<R> z);

// 0F1(; b; z)
template<IsNativeFloat T>
[[nodiscard]] T hyperg0F1(T b, T z);
template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> hyperg0F1(Complex<R> b, Complex<R> z);

Definition: ${}_2F_1(a,b;c;z) = \displaystyle\sum_{k=0}^\infty \frac{(a)_k (b)_k}{(c)_k}\frac{z^k}{k!}$ ($(x)_k$ is the Pochhammer symbol). Unifies many special functions (Legendre, Chebyshev, elliptic integrals, and more). confHyperg ${}_1F_1$ is the parent of the Bessel, Laguerre, and error functions, and hyperg0F1 is related to the Bessel functions $J_\nu$. The complex versions' SpecialFunctionScalar R admits double / Float.

pFq

// native (numerator parameters a, denominator parameters b as span or initializer_list)
template<IsNativeFloat T>
[[nodiscard]] T pFq(std::span<const T> a, std::span<const T> b, T z);
template<IsNativeFloat T>
[[nodiscard]] T pFq(std::initializer_list<T> a, std::initializer_list<T> b, T z);

// arbitrary-precision Float
[[nodiscard]] Float pFq(std::span<const Float> a, std::span<const Float> b,
                        const Float& z, int precision);
[[nodiscard]] Float pFq(std::initializer_list<Float> a, std::initializer_list<Float> b,
                        const Float& z, int precision);

Definition: the generalized hypergeometric ${}_pF_q(a_1,\dots,a_p; b_1,\dots,b_q; z) = \displaystyle\sum_{k=0}^\infty \frac{(a_1)_k\cdots(a_p)_k}{(b_1)_k\cdots(b_q)_k}\frac{z^k}{k!}$. Pass the array of numerator parameters a ($p$ of them) and the array of denominator parameters b ($q$ of them) as a std::span or std::initializer_list. ${}_2F_1$, ${}_1F_1$, and so on are special cases.

// Example: 2F1(1, 2; 3; 0.5) via pFq
double v = sangi::special::pFq<double>({1.0, 2.0}, {3.0}, 0.5);

mittagLeffler

// two-parameter E_{alpha,beta}(z)
template<IsNativeFloat T>
[[nodiscard]] T mittagLeffler(T alpha, T beta, T z);

// single-parameter E_alpha(z) = E_{alpha,1}(z)
template<IsNativeFloat T>
[[nodiscard]] T mittagLeffler(T alpha, T z);

Definition: $E_{\alpha,\beta}(z) = \displaystyle\sum_{k=0}^\infty \frac{z^k}{\Gamma(\alpha k + \beta)}$. Generalizes $E_{1,1}(z) = e^z$ and appears in the solutions of fractional differential equations, anomalous diffusion, and viscoelastic models. With two arguments it is the single-parameter version $E_\alpha(z) = E_{\alpha,1}(z)$; with three, the two-parameter version.

Other Special Functions

Houses Lambert W, Owen's T, the Coulomb wave functions, Fermi-Dirac integrals, Debye functions, transport / synchrotron integrals, angular-momentum coupling coefficients (Wigner), and spheroidal wave functions. spheroidal.hpp is not part of special.hpp and must be included individually.

Lambert W & Owen's T (lambert_w.hpp, owens_t.hpp)

lambertW0 / lambertWm1

template<IsNativeFloat T>
[[nodiscard]] T lambertW0(T x);
[[nodiscard]] Complex<double> lambertW0(const Complex<double>& z);
[[nodiscard]] Complex<Float>  lambertW0(const Complex<Float>& z, int precision);

template<IsNativeFloat T>
[[nodiscard]] T lambertWm1(T x);
[[nodiscard]] Complex<double> lambertWm1(const Complex<double>& z);
[[nodiscard]] Complex<Float>  lambertWm1(const Complex<Float>& z, int precision);

Definition: $W(x)$ is the solution of $W e^{W} = x$. Being multivalued, over the reals it splits into two branches ($W_0$ principal, $W_{-1}$ lower), which merge at $x = -1/e$. owensT(h, a) is $T(h,a) = \frac{1}{2\pi}\int_0^a \frac{e^{-h^2(1+t^2)/2}}{1+t^2}dt$, used in bivariate-normal probability calculations.

Coulomb, Fermi-Dirac, Debye, transport

Angular-momentum coupling (coupling.hpp)

wigner3j / clebschGordan

inline double wigner3j(int two_j1, int two_j2, int two_j3,
                       int two_m1, int two_m2, int two_m3);

inline double clebschGordan(int two_j1, int two_m1,
                            int two_j2, int two_m2,
                            int two_J, int two_M);

inline double wigner6j(int two_j1, int two_j2, int two_j3,
                       int two_j4, int two_j5, int two_j6);

inline double wigner9j(int two_j1, int two_j2, int two_j3,
                       int two_j4, int two_j5, int two_j6,
                       int two_j7, int two_j8, int two_j9);

Argument conventions: angular momenta and magnetic quantum numbers are all passed as doubled integers (two_j = $2j$, two_m = $2m$). This lets half-integer spins be handled exactly using integers only. Used in angular-momentum coupling in quantum mechanics and in selection-rule calculations in spectroscopy. The return type is double (not a template).

Spheroidal wave functions (spheroidal.hpp)

Takes the orders m, n ($n \ge m \ge 0$), the spheroidal parameter c, and the angular variable eta. These are the separated solutions of the wave equation in spheroidal coordinates, arising in antenna and acoustic-scattering analysis.

Example

#include <math/special/special.hpp>
#include <iostream>
using namespace sangi::special;

int main() {
    // Gamma function: Gamma(1/2) = sqrt(pi)
    std::cout << "gamma(0.5)  = " << gamma(0.5) << '\n';
    // Result: gamma(0.5)  = 1.77245385090552   (sqrt(pi))

    // Bessel function of the first kind: J_0(0) = 1
    std::cout << "besselJ(0,0)= " << besselJ(0, 0.0) << '\n';
    // Result: besselJ(0,0)= 1

    // Error function: erf(1)
    std::cout << "erf(1.0)    = " << erf(1.0) << '\n';
    // Result: erf(1.0)    = 0.842700792949715

    // Complete elliptic integral of the first kind: K(0) = pi/2
    std::cout << "ellipticK(0)= " << ellipticK(0.0) << '\n';
    // Result: ellipticK(0)= 1.5707963267949   (pi/2)

    return 0;
}

To call without the using directive, write the fully qualified name, e.g. sangi::special::gamma(0.5). For complex arguments or arbitrary precision (Complex<Float> plus a precision argument), pick a family that has the matching overload (gamma, error function, Bessel, Airy, elliptic, zeta, hypergeometric, and so on).

// Complex argument example: Riemann zeta on the critical line, zeta(1/2 + 14.13i)
#include <math/special/special.hpp>
using namespace sangi;
using namespace sangi::special;

Complex<double> s(0.5, 14.134725);
Complex<double> z = riemannZeta(s);   // near the first nontrivial zero -> |z| ~ 0

// A family that needs an individual include: generalized hypergeometric 2F1(1,2;3;0.5)
#include <math/special/hypergeometric_pq.hpp>
double f = pFq<double>({1.0, 2.0}, {3.0}, 0.5);

Related Mathematical Background

The following articles explain the mathematical concepts underlying the Special Functions module.

• A Tour of Special Functions — A map of the special functions and how they relate
• Gamma Function — Analytic continuation of the factorial and the reflection formula
• Bessel Functions — From the wave equation in cylindrical coordinates
• Error Function — The integral of the normal distribution
• Elliptic Integrals — The pendulum period and Carlson symmetric forms
• Orthogonal Polynomials — The Legendre / Hermite / Laguerre / Chebyshev system
• Hypergeometric Functions — The ${}_pF_q$ that unifies the special functions
• Riemann Zeta Function — Analytic continuation and the functional equation