# 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.$