.. _hypgeom:
**hypgeom.h** -- support for hypergeometric series
===============================================================================
This module provides functions for high-precision evaluation of series
of the form
.. math ::
\sum_{k=0}^{n-1} \frac{A(k)}{B(k)} \prod_{j=1}^k \frac{P(k)}{Q(k)} z^k
where `A, B, P, Q` are polynomials. The present version only supports
`A, B, P, Q \in \mathbb{Z}[k]` (represented using the
FLINT *fmpz_poly_t* type). This module also provides functions
for high-precision evaluation of infinite series (`n \to \infty`),
with automatic, rigorous error bounding.
Note that we can standardize to `A = B = 1` by
setting `\tilde P(k) = P(k) A(k) B(k-1), \tilde Q(k) = Q(k) A(k-1) B(k)`.
However, separating out `A` and `B` is convenient and improves
efficiency during evaluation.
Strategy for error bounding
-------------------------------------------------------------------------------
We wish to evaluate `S(z) = \sum_{k=0}^{\infty} T(k) z^k` where `T(k)`
satisfies `T(0) = 1` and
.. math ::
T(k) = R(k) T(k-1) = \left( \frac{P(k)}{Q(k)} \right) T(k-1)
for given polynomials
.. math ::
P(k) = a_p k^p + a_{p-1} k^{p-1} + \ldots a_0
Q(k) = b_q k^q + b_{q-1} k^{q-1} + \ldots b_0.
For convergence, we require `p < q`, or `p = q` with `|z| |a_p| < |b_q|`.
We also assume that `P(k)` and `Q(k)` have no roots among the positive
integers (if there are positive integer roots, the sum is either finite
or undefined). With these conditions satisfied, our goal is to find a
parameter `n \ge 0` such that
.. math ::
\left\lvert \sum_{k=n}^{\infty} T(k) z^k \right\rvert \le 2^{-d}.
We can rewrite the hypergeometric term ratio as
.. math ::
z R(k) = z \frac{P(k)}{Q(k)} =
z \left( \frac{a_p}{b_q} \right) \frac{1}{k^{q-p}} F(k)
where
.. math ::
F(k) = \frac{
1 + \tilde a_{1} / k + \tilde a_{2} / k^2 + \ldots + \tilde a_q / k^p
}{
1 + \tilde b_{1} / k + \tilde b_{2} / k^2 + \ldots + \tilde b_q / k^q
} = 1 + O(1/k)
and where `\tilde a_i = a_{p-i} / a_p`, `\tilde b_i = b_{q-i} / b_q`.
Next, we define
.. math ::
C = \max_{1 \le i \le p} |\tilde a_i|^{(1/i)},
\quad D = \max_{1 \le i \le q} |\tilde b_i|^{(1/i)}.
Now, if `k > C`, the magnitude of the numerator of `F(k)` is
bounded from above by
.. math ::
1 + \sum_{i=1}^p \left(\frac{C}{k}\right)^i \le 1 + \frac{C}{k-C}
and if `k > 2D`, the magnitude of the denominator of `F(k)` is bounded
from below by
.. math ::
1 - \sum_{i=1}^q \left(\frac{D}{k}\right)^i \ge 1 + \frac{D}{D-k}.
Putting the inequalities together gives the following bound,
valid for `k > K = \max(C, 2D)`:
.. math ::
|F(k)| \le \frac{k (k-D)}{(k-C)(k-2D)} = \left(1 + \frac{C}{k-C} \right)
\left(1 + \frac{D}{k-2D} \right).
Let `r = q-p` and `\tilde z = |z a_p / b_q|`. Assuming
`k > \max(C, 2D, {\tilde z}^{1/r})`, we have
.. math ::
|z R(k)| \le G(k) = \frac{\tilde z F(k)}{k^r}
where `G(k)` is monotonically decreasing. Now we just need to find an
`n` such that `G(n) < 1` and for which `|T(n)| / (1 - G(n)) \le 2^{-d}`.
This can be done by computing a floating-point guess for `n` then
trying successively larger values.
This strategy leaves room for some improvement. For example, if
`\tilde b_1` is positive and large, the bound `B` becomes very pessimistic
(a larger positive `\tilde b_1` causes faster convergence,
not slower convergence).
Types, macros and constants
-------------------------------------------------------------------------------
.. type:: hypgeom_struct
.. type:: hypgeom_t
Stores polynomials *A*, *B*, *P*, *Q* and precomputed bounds,
representing a fixed hypergeometric series.
Memory management
-------------------------------------------------------------------------------
.. function:: void hypgeom_init(hypgeom_t hyp)
.. function:: void hypgeom_clear(hypgeom_t hyp)
Error bounding
-------------------------------------------------------------------------------
.. function:: slong hypgeom_estimate_terms(const mag_t z, int r, slong d)
Computes an approximation of the largest `n` such
that `|z|^n/(n!)^r = 2^{-d}`, giving a first-order estimate of the
number of terms needed to approximate the sum of a hypergeometric
series of weight `r \ge 0` and argument `z` to an absolute
precision of `d \ge 0` bits. If `r = 0`, the direct solution of the
equation is given by `n = (\log(1-z) - d \log 2) / \log z`.
If `r > 0`, using `\log n! \approx n \log n - n` gives an equation
that can be solved in terms of the Lambert *W*-function as
`n = (d \log 2) / (r\,W\!(t))` where
`t = (d \log 2) / (e r z^{1/r})`.
The evaluation is done using double precision arithmetic.
The function aborts if the computed value of `n` is greater
than or equal to LONG_MAX / 2.
.. function:: slong hypgeom_bound(mag_t error, int r, slong C, slong D, slong K, const mag_t TK, const mag_t z, slong prec)
Computes a truncation parameter sufficient to achieve *prec* bits
of absolute accuracy, according to the strategy described above.
The input consists of `r`, `C`, `D`, `K`, precomputed bound for `T(K)`,
and `\tilde z = z (a_p / b_q)`, such that for `k > K`, the hypergeometric
term ratio is bounded by
.. math ::
\frac{\tilde z}{k^r} \frac{k(k-D)}{(k-C)(k-2D)}.
Given this information, we compute a `\varepsilon` and an
integer `n` such that
`\left| \sum_{k=n}^{\infty} T(k) \right| \le \varepsilon \le 2^{-\mathrm{prec}}`.
The output variable *error* is set to the value of `\varepsilon`,
and `n` is returned.
.. function:: void hypgeom_precompute(hypgeom_t hyp)
Precomputes the bounds data `C`, `D`, `K` and an upper bound for `T(K)`.
Summation
-------------------------------------------------------------------------------
.. function:: void arb_hypgeom_sum(arb_t P, arb_t Q, const hypgeom_t hyp, const slong n, slong prec)
Computes `P, Q` such that `P / Q = \sum_{k=0}^{n-1} T(k)` where `T(k)`
is defined by *hyp*,
using binary splitting and a working precision of *prec* bits.
.. function:: void arb_hypgeom_infsum(arb_t P, arb_t Q, hypgeom_t hyp, slong tol, slong prec)
Computes `P, Q` such that `P / Q = \sum_{k=0}^{\infty} T(k)` where `T(k)`
is defined by *hyp*, using binary splitting and
working precision of *prec* bits.
The number of terms is chosen automatically to bound the
truncation error by at most `2^{-\mathrm{tol}}`.
The bound for the truncation error is included in the output
as part of *P*.