Algorithms for hypergeometric functions

The algorithms used to compute hypergeometric functions are described in [Joh2016]. Here, we state the most important error bounds.

Convergent series

Let

\[T(k) = \frac{\prod_{i=0}^{p-1} (a_i)_k}{\prod_{i=0}^{q-1} (b_i)_k} z^k.\]

We compute a factor C such that

\[\left|\sum_{k=n}^{\infty} T(k)\right| \le C |T(n)|.\]

We check that \(\operatorname{Re}(b+n) > 0\) for all lower parameters b. If this does not hold, C is set to infinity. Otherwise, we cancel out pairs of parameters \(a\) and \(b\) against each other. We have

\[\left|\frac{a+k}{b+k}\right| = \left|1 + \frac{a-b}{b+k}\right| \le 1 + \frac{|a-b|}{|b+n|}\]

and

\[\left|\frac{1}{b+k}\right| \le \frac{1}{|b+n|}\]

for all \(k \ge n\). This gives us a constant D such that \(T(k+1) \le D T(k)\) for all \(k \ge n\). If \(D \ge 1\), we set C to infinity. Otherwise, we take \(C = \sum_{k=0}^{\infty} D^k = (1-D)^{-1}\).

Convergent series of power series

The same principle is used to get tail bounds for with \(a_i, b_i, z \in \mathbb{C}[[x]]\), or more precisely, bounds for each coefficient in \(\sum_{k=N}^{\infty} T(k) \in \mathbb{C}[[x]] / \langle x^n \rangle\) given \(a_i, b_i, z \in \mathbb{C}[[x]] / \langle x^n \rangle\). First, we fix some notation, assuming that \(A\) and \(B\) are power series:

• \(A_{[k]}\) denotes the coefficient of \(x^k\) in \(A\), and \(A_{[m:n]}\) denotes the power series \(\sum_{k=m}^{n-1} A_{[k]} x^k\).

• \(|A|\) denotes \(\sum_{k=0}^{\infty} |A_{[k]}| x^k\) (this can be viewed as an element of \(\mathbb{R}_{\ge 0}[[x]]\)).

• \(A \le B\) signifies that \(|A|_{[k]} \le |B|_{[k]}\) holds for all \(k\).

• We define \(\mathcal{R}(B) = |B_{[0]}| - |B_{[1:\infty]}|\).

Using the formulas

\[(A B)_{[k]} = \sum_{j=0}^k A_{[j]} B_{[k-j]}, \quad (1 / B)_{[k]} = \frac{1}{B_{[0]}} \sum_{j=1}^k -B_{[j]} (1/B)_{[k-j]},\]

it is easy prove the following bounds for the coefficients of sums, products and quotients of formal power series:

\[|A + B| \le |A| + |B|, \quad |A B| \le |A| |B|, \quad |A / B| \le |A| / \mathcal{R}(B).\]

If \(p \le q\) and \(\operatorname{Re}({b_i}_{[0]}+N) > 0\) for all \(b_i\), then we may take

\[D = |z| \, \prod_{i=1}^p \left(1 + \frac{|a_i-b_i|}{\mathcal{R}(b_i+N)}\right) \prod_{i=p+1}^{q} \frac{1}{\mathcal{R}(b_i + N)}.\]

If \(D_{[0]} < 1\),then \((1 - D)^{-1} |T(n)|\) gives the error bound.

Note when adding and multiplying power series with (complex) interval coefficients, we can use point-valued upper bounds for the absolute values instead of performing interval arithmetic throughout. For \(\mathcal{R}(B)\), we must then pick a lower bound for \(|B_{[0]}|\) and upper bounds for the coefficients of \(|B_{[1:\infty]}|\).

Asymptotic series for the confluent hypergeometric function

Let \(U(a,b,z)\) denote the confluent hypergeometric function of the second kind with the principal branch cut, and let \(U^{*} = z^a U(a,b,z)\). For all \(z \ne 0\) and \(b \notin \mathbb{Z}\) (but valid for all \(b\) as a limit), we have (DLMF 13.2.42)

\[U(a,b,z) = \frac{\Gamma(1-b)}{\Gamma(a-b+1)} M(a,b,z) + \frac{\Gamma(b-1)}{\Gamma(a)} z^{1-b} M(a-b+1,2-b,z).\]

Moreover, for all \(z \ne 0\) we have

\[\frac{{}_1F_1(a,b,z)}{\Gamma(b)} = \frac{(-z)^{-a}}{\Gamma(b-a)} U^{*}(a,b,z) + \frac{z^{a-b} e^z}{\Gamma(a)} U^{*}(b-a,b,-z)\]

which is equivalent to DLMF 13.2.41 (but simpler in form).

We have the asymptotic expansion

\[U^{*}(a,b,z) \sim {}_2F_0(a, a-b+1, -1/z)\]

where \({}_2F_0(a,b,z)\) denotes a formal hypergeometric series, i.e.

\[U^{*}(a,b,z) = \sum_{k=0}^{n-1} \frac{(a)_k (a-b+1)_k}{k! (-z)^k} + \varepsilon_n(z).\]

The error term \(\varepsilon_n(z)\) is bounded according to DLMF 13.7. A case distinction is made depending on whether \(z\) lies in one of three regions which we index by \(R\). Our formula for the error bound increases with the value of \(R\), so we can always choose the larger out of two indices if \(z\) lies in the union of two regions.

Let \(r = |b-2a|\). If \(\operatorname{Re}(z) \ge r\), set \(R = 1\). Otherwise, if \(\operatorname{Im}(z) \ge r\) or \(\operatorname{Re}(z) \ge 0 \land |z| \ge r\), set \(R = 2\). Otherwise, if \(|z| \ge 2r\), set \(R = 3\). Otherwise, the bound is infinite. If the bound is finite, we have

\[|\varepsilon_n(z)| \le 2 \alpha C_n \left|\frac{(a)_n (a-b+1)_n}{n! z^n} \right| \exp(2 \alpha \rho C_1 / |z|)\]

in terms of the following auxiliary quantities

\[\sigma = |(b-2a)/z|\]

\[\begin{split}C_n = \begin{cases} 1 & \text{if } R = 1 \\ \chi(n) & \text{if } R = 2 \\ (\chi(n) + \sigma \nu^2 n) \nu^n & \text{if } R = 3 \end{cases}\end{split}\]

\[\nu = \left(\tfrac{1}{2} + \tfrac{1}{2}\sqrt{1-4\sigma^2}\right)^{-1/2} \le 1 + 2 \sigma^2\]

\[\chi(n) = \sqrt{\pi} \Gamma(\tfrac{1}{2}n+1) / \Gamma(\tfrac{1}{2} n + \tfrac{1}{2})\]

\[\begin{split}\sigma' = \begin{cases} \sigma & \text{if } R \ne 3 \\ \nu \sigma & \text{if } R = 3 \end{cases}\end{split}\]

\[\alpha = (1 - \sigma')^{-1}\]

\[\rho = \tfrac{1}{2} |2a^2-2ab+b| + \sigma' (1+ \tfrac{1}{4} \sigma') (1-\sigma')^{-2}\]

Asymptotic series for Airy functions

Error bounds are based on Olver (DLMF section 9.7). For \(\arg(z) < \pi\) and \(\zeta = (2/3) z^{3/2}\), we have

\[\operatorname{Ai}(z) = \frac{e^{-\zeta}}{2 \sqrt{\pi} z^{1/4}} \left[S_n(\zeta) + R_n(z)\right], \quad \operatorname{Ai}'(z) = -\frac{z^{1/4} e^{-\zeta}}{2 \sqrt{\pi}} \left[(S'_n(\zeta) + R'_n(z)\right]\]

\[S_n(\zeta) = \sum_{k=0}^{n-1} (-1)^k \frac{u(k)}{\zeta^k}, \quad S'_n(\zeta) = \sum_{k=0}^{n-1} (-1)^k \frac{v(k)}{\zeta^k}\]

\[u(k) = \frac{(1/6)_k (5/6)_k}{2^k k!}, \quad v(k) = \frac{6k+1}{1-6k} u(k).\]

Assuming that n is positive, the error terms are bounded by

\[|R_n(z)| \le C |u(n)| |\zeta|^{-n}, \quad |R'_n(z)| \le C |v(n)| |\zeta|^{-n}\]

where

\[\begin{split}C = \begin{cases} 2 \exp(7 / (36 |\zeta|)) & |\arg(z)| \le \pi/3 \\ 2 \chi(n) \exp(7 \pi / (72 |\zeta|)) & \pi/3 \le |\arg(z)| \le 2\pi/3 \\ 4 \chi(n) \exp(7 \pi / (36 |\operatorname{re}(\zeta)|)) |\cos(\arg(\zeta))|^{-n} & 2\pi/3 \le |\arg(z)| < \pi. \end{cases}\end{split}\]

For computing Bi when z is roughly in the positive half-plane, we use the connection formulas

\[ \begin{align}\begin{aligned}\operatorname{Bi}(z) = -i (2 w^{+1} \operatorname{Ai}(z w^{-2}) - \operatorname{Ai}(z))\\\operatorname{Bi}(z) = +i (2 w^{-1} \operatorname{Ai}(z w^{+2}) - \operatorname{Ai}(z))\end{aligned}\end{align} \]

where \(w = \exp(\pi i/3)\). Combining roots of unity gives

\[\operatorname{Bi}(z) = \frac{1}{2 \sqrt{\pi} z^{1/4}} [2X + iY]\]

\[\operatorname{Bi}(z) = \frac{1}{2 \sqrt{\pi} z^{1/4}} [2X - iY]\]

\[X = \exp(+\zeta) [S_n(-\zeta) + R_n(z w^{\mp 2})], \quad Y = \exp(-\zeta) [S_n(\zeta) + R_n(z)]\]

where the upper formula is valid for \(-\pi/3 < \arg(z) < \pi\) and the lower formula is valid for \(-\pi < \arg(z) < \pi/3\). We proceed analogously for the derivative of Bi.

In the negative half-plane, we use the connection formulas

\[\operatorname{Ai}(z) = e^{+\pi i/3} \operatorname{Ai}(z_1) + e^{-\pi i/3} \operatorname{Ai}(z_2)\]

\[\operatorname{Bi}(z) = e^{-\pi i/6} \operatorname{Ai}(z_1) + e^{+\pi i/6} \operatorname{Ai}(z_2)\]

where \(z_1 = -z e^{+\pi i/3}\), \(z_2 = -z e^{-\pi i/3}\). Provided that \(|\arg(-z)| < 2 \pi / 3\), we have \(|\arg(z_1)|, |\arg(z_2)| < \pi\), and thus the asymptotic expansion for Ai can be used. As before, we collect roots of unity to obtain

\[\operatorname{Ai}(z) = A_1 [S_n(i \zeta) + R_n(z_1)] + A_2 [S_n(-i \zeta) + R_n(z_2)]\]

\[\operatorname{Bi}(z) = A_3 [S_n(i \zeta) + R_n(z_1)] + A_4 [S_n(-i \zeta) + R_n(z_2)]\]

where \(\zeta = (2/3) (-z)^{3/2}\) and

\[A_1 = \frac{\exp(-i (\zeta - \pi/4))}{2 \sqrt{\pi} (-z)^{1/4}}, \quad A_2 = \frac{\exp(+i (\zeta - \pi/4))}{2 \sqrt{\pi} (-z)^{1/4}}, \quad A_3 = -i A_1, \quad A_4 = +i A_2.\]

The differentiated formulas are analogous.

Corner case of the Gauss hypergeometric function

In the corner case where \(z\) is near \(\exp(\pm \pi i / 3)\), none of the linear fractional transformations is effective. In this case, we use Taylor series to analytically continue the solution of the hypergeometric differential equation from the origin. The function \(f(z) = {}_2F_1(a,b,c,z_0+z)\) satisfies

\[f''(z) = -\frac{((z_0+z)(a+b+1)-c)}{(z_0+z)(z_0-1+z)} f'(z) - \frac{a b}{(z_0+z)(z_0-1+z)} f(z).\]

Knowing \(f(0), f'(0)\), we can compute the consecutive derivatives recursively, and evaluating the truncated Taylor series allows us to compute \(f(z), f'(z)\) to high accuracy for sufficiently small \(z\). Some experimentation showed that two continuation steps

\[0 \quad \to \quad 0.375 \pm 0.625i \quad \to \quad 0.5 \pm 0.8125i \quad \to \quad z\]

gives good performance. Error bounds for the truncated Taylor series are obtained using the Cauchy-Kovalevskaya majorant method, following the outline in [Hoe2001]. The differential equation is majorized by

\[g''(z) = \frac{N+1}{2} \left( \frac{\nu}{1-\nu z} \right) g'(z) + \frac{(N+1)N}{2} \left( \frac{\nu}{1-\nu z} \right)^2 g(z)\]

provided that \(N\) and \(\nu \ge \max(1/|z_0|, 1/|z_0-1|)\) are chosen sufficiently large. It follows that we can compute explicit numbers \(A, N, \nu\) such that the simple solution \(g(z) = A (1-\nu z)^{-N}\) of the differential equation provides the bound

\[|f_{[k]}| \le g_{[k]} = A {{N+k} \choose k} \nu^k.\]