.. _algorithms_hypergeometric: Algorithms for hypergeometric functions =============================================================================== The algorithms used to compute hypergeometric functions are described in [Joh2016]_. Here, we state the most important error bounds. .. _algorithms_hypergeometric_convergent: Convergent series ------------------------------------------------------------------------------- Let .. math :: 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 .. math :: \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 .. math :: \left|\frac{a+k}{b+k}\right| = \left|1 + \frac{a-b}{b+k}\right| \le 1 + \frac{|a-b|}{|b+n|} and .. math :: \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 .. math :: (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: .. math :: |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 .. math :: 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]}|`. .. _algorithms_hypergeometric_asymptotic_confluent: 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) .. math :: 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 .. math :: \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 .. math :: 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. .. math :: 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 .. math :: |\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 .. math :: \sigma = |(b-2a)/z| .. math :: 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} .. math :: \nu = \left(\tfrac{1}{2} + \tfrac{1}{2}\sqrt{1-4\sigma^2}\right)^{-1/2} \le 1 + 2 \sigma^2 .. math :: \chi(n) = \sqrt{\pi} \Gamma(\tfrac{1}{2}n+1) / \Gamma(\tfrac{1}{2} n + \tfrac{1}{2}) .. math :: \sigma' = \begin{cases} \sigma & \text{if } R \ne 3 \\ \nu \sigma & \text{if } R = 3 \end{cases} .. math :: \alpha = (1 - \sigma')^{-1} .. math :: \rho = \tfrac{1}{2} |2a^2-2ab+b| + \sigma' (1+ \tfrac{1}{4} \sigma') (1-\sigma')^{-2} .. _algorithms_hypergeometric_asymptotic_airy: 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 .. math :: \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] .. math :: 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} .. math :: 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 .. math :: |R_n(z)| \le C |u(n)| |\zeta|^{-n}, \quad |R'_n(z)| \le C |v(n)| |\zeta|^{-n} where .. math :: 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} For computing Bi when *z* is roughly in the positive half-plane, we use the connection formulas .. math :: \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)) where `w = \exp(\pi i/3)`. Combining roots of unity gives .. math :: \operatorname{Bi}(z) = \frac{1}{2 \sqrt{\pi} z^{1/4}} [2X + iY] .. math :: \operatorname{Bi}(z) = \frac{1}{2 \sqrt{\pi} z^{1/4}} [2X - iY] .. math :: 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 .. math :: \operatorname{Ai}(z) = e^{+\pi i/3} \operatorname{Ai}(z_1) + e^{-\pi i/3} \operatorname{Ai}(z_2) .. math :: \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 .. math :: \operatorname{Ai}(z) = A_1 [S_n(i \zeta) + R_n(z_1)] + A_2 [S_n(-i \zeta) + R_n(z_2)] .. math :: \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 .. math :: 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 .. math :: 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 .. math :: 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 .. math :: 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 .. math :: |f_{[k]}| \le g_{[k]} = A {{N+k} \choose k} \nu^k.