Module special 

Special functions: the error function and its inverse, the Gaussian tail, the Fresnel integrals, Bessel functions, the gamma function and the regularized incomplete gamma function, the binomial law, the regularized incomplete beta function and its inverse, the exponential integral, the generalized hypergeometric function ${1}F{2}$ and the Hermite polynomials (all re-exported at the root).

The power series share one private summation (the terms until the last falls below $10^{-17}$ of the sum), the continued fractions one modified Lentz iteration (the partial denominators floored at $10^{-300}$, until a factor is within $10^{-16}$ of 1), and the periodic integrands one midpoint rule; the order of the operations of each function is the one it had before they were shared (the same results, bit for bit).

References

1. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes, 3rd ed., Cambridge University Press, 2007, §5.2 (the modified Lentz method), §6.2, §6.3 and §6.4.
2. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, 1965.

Functions

bessel_i0e

The scaled modified Bessel function of the first kind of order zero, $e^{-|x|} I_{0}(x)$, with $I_{0}(x) = \frac{1}{\pi} \int_{0}^{\pi} e^{x\cos\theta} , d\theta$: the midpoint rule on the periodic integrand $e^{|x|(\cos\theta - 1)}$ for $|x| \le 60$ (exponentially convergent), and the asymptotic series $\frac{1}{\sqrt{2\pi x}} \sum_{k} \frac{((2k-1)!!)^{2}}{k! , (8x)^{k}}$ beyond; to about $10^{-14}$ relative. The density of the Rice distribution, $y e^{-(y - v)^{2}/2} , e^{-yv} I_{0}(yv)$, needs the scaled form.

bessel_j0

The Bessel function of the first kind of order zero, $J_{0}(x) = \frac{1}{\pi} \int_{0}^{\pi} \cos(x \sin\theta), d\theta$, by the midpoint rule on that periodic integrand (exponentially convergent: $64 + 2|x|$ points give the machine precision), as for the correlation $J_{0}(2\pi f_{D} \tau)$ of Jakes’ fading.

bessel_k

The modified Bessel function of the second kind $K_{\nu}(x)$ of real order: exp(ln_bessel_k).

beta_inc

The regularized incomplete beta function $I_{x}(a, b) = \frac{1}{B(a, b)} \int_{0}^{x} u^{a-1} (1 - u)^{b-1}, du$, $a, b > 0$, $x \in [0, 1]$: the distribution function of a $\mathrm{Beta}(a, b)$ variable, and the binomial tail $P(X \ge k) = I_{p}(k, n - k + 1)$ of $X \sim \mathrm{Bin}(n, p)$. The continued fraction of Press et al. (Numerical Recipes, §6.4) by the modified Lentz method, on $I_{1-x}(b, a) = 1 - I_{x}(a, b)$ when $x > (a + 1)/(a + b + 2)$, where it converges faster; relative error about $10^{-14}$.

beta_inv

The inverse of beta_inc in $x$: the $q$ quantile of $\mathrm{Beta}(a, b)$, by bisection on $[0, 1]$ (0 and 1 at the ends of $q$).

erf

The error function $\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^{2}}, dt$: its Maclaurin series $\frac{2}{\sqrt{\pi}} \sum_{k} \frac{(-1)^{k} x^{2k+1}}{k!,(2k+1)}$ for $|x| < 2$ (relative error about $10^{-15}$, small arguments included), $1 - \operatorname{erfc}(x)$ beyond (absolute error about $10^{-16}$). Odd.

erfc

Complementary error function $\operatorname{erfc}(x) = \frac{2}{\sqrt{\pi}} \int_{x}^{\infty} e^{-t^{2}}, dt$.

erfcinv

Inverse of erfc on (0, 2), by bisection (±∞ at the ends).

expint_e1

The exponential integral $E_{1}(x) = \int_{x}^{\infty} \frac{e^{-t}}{t}, dt$, $x > 0$: the series $-\gamma - \ln x - \sum_{k \ge 1} \frac{(-x)^{k}}{k, k!}$ for $x \le 1$, the continued fraction $\cfrac{e^{-x}}{x + 1 - \cfrac{1}{x + 3 - \cfrac{4}{x + 5 - \dotsb}}}$ beyond (Press et al., Numerical Recipes, §6.3); relative error about $10^{-15}$. $e^{1/\rho} E_{1}(1/\rho)$ is the ergodic capacity (nats) of the Rayleigh channel at the mean SNR $\rho$.

fresnel

The Fresnel integrals $C(x) = \int_{0}^{x} \cos(\pi t^{2}/2), dt$ and $S(x) = \int_{0}^{x} \sin(\pi t^{2}/2), dt$, as (C, S), to about $10^{-15}$.

gamma_p

The regularized lower incomplete gamma function $P(a, x) = \frac{1}{\Gamma(a)} \int_{0}^{x} u^{a-1} e^{-u}, du$, $a > 0$, $x \ge 0$: the distribution function of a gamma variable of shape $a$ and unit scale (a $\chi^{2}$ variable with $2a$ degrees of freedom, halved). The series $e^{-x} x^{a} \sum_{k} x^{k} / \Gamma(a + k + 1)$ for $x < a + 1$, and $1 - Q(a, x)$ with the continued fraction of $Q$ (modified Lentz) beyond (Press et al., Numerical Recipes, §6.2); relative error about $10^{-14}$ on the side computed directly.

hermite_h

The physicists’ Hermite polynomial $H_{n}(x)$, by the recurrence $H_{n+1}(x) = 2x H_{n}(x) - 2n H_{n-1}(x)$ from $H_{0} = 1$, $H_{1} = 2x$ (Abramowitz & Stegun, eq. 22.7.13); the Hermite prototype of FBMC is a sum of $H_{4i}$.

hyp1f2

The generalized hypergeometric function ${1}F{2}(a; b_{1}, b_{2}; z) = \sum_{k \ge 0} \frac{(a){k}}{(b{1}){k} (b{2}){k}} \frac{z^{k}}{k!}$ by its power series, which converges for every $z$ ($(x){k}$ the rising factorial; $b_{1}, b_{2}$ not zero or negative integers). The ICI of OFDM under Jakes’ fading is $1 - {}{1}F{2}(\tfrac{1}{2}; \tfrac{3}{2}, 2; -(\pi f_{D} T_{u})^{2})$ (Robertson & Kaiser 1999; Nissel, Schwarz & Rupp 2017, eq. (42)).

ln_bessel_k

The natural logarithm of the modified Bessel function of the second kind $K_{\nu}(x)$ of real order $\nu$ and $x > 0$, from $K_{\nu}(x) = \int_{0}^{\infty} e^{-x \cosh t} \cosh(\nu t), dt$ by the trapezoidal rule (exponentially convergent on this even, entire integrand), summed in the log domain so that neither the large values at small $x$ and large $\nu$ nor the small ones at large $x$ overflow; to about $10^{-13}$ relative. The limits: $-\infty$ for $x = +\infty$ ($K_{\nu}(x) \to 0$), $+\infty$ for an infinite order, and the leading term $\ln\left(\tfrac{1}{2} \Gamma(\nu)\right) + \nu \ln(2/x)$ when $\nu / x$ overflows (a subnormal $x$; for $\nu \ge 0.05$, where the next term is below $e^{-70}$ relative); NaN for $x \le 0$, a NaN argument, or an overflowing $\nu / x$ with a smaller order.

ln_binomial_pmf

The logarithms of the binomial probabilities $\ln \binom{n}{k} + k \ln p + (n - k) \ln (1 - p)$, $k = 0, \dots, n$: the law of the number of successes among $n$ trials of probability $p \in [0, 1]$ (the coefficients by the running product $\binom{n}{k} = \binom{n}{k-1} (n + 1 - k)/k$; at $p = 0$ or 1 the impossible counts are $-\infty$).

ln_gamma

The logarithm of the gamma function, $\ln \Gamma(x)$ for $x > 0$, by the Lanczos approximation ($g = 7$, nine coefficients; relative error about $10^{-15}$) and the reflection formula $\Gamma(x)\Gamma(1-x) = \pi / \sin(\pi x)$ below $1/2$.

normal_pdf

The standard normal density $\varphi(x) = e^{-x^{2}/2}/\sqrt{2\pi}$ (its distribution function is $1 - Q(x)$, qfunc).

qfunc

Gaussian tail probability Q(x) = P(N(0,1) > x) = ½ erfc(x/√2).

qfuncinv

Inverse of qfunc on $(0, 1)$: $Q^{-1}(p) = \sqrt{2}, \operatorname{erfc}^{-1}(2p)$ (±∞ at the ends).