mirror of
https://github.com/vale981/arb
synced 2025-03-06 01:41:39 -05:00
328 lines
10 KiB
ReStructuredText
328 lines
10 KiB
ReStructuredText
.. _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.
|
|
|