# hypgeom.h – support for hypergeometric series¶

This module provides functions for high-precision evaluation of series of the form

$\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

$T(k) = R(k) T(k-1) = \left( \frac{P(k)}{Q(k)} \right) T(k-1)$

for given polynomials

\begin{align}\begin{aligned}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.\end{aligned}\end{align}

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

$\left\lvert \sum_{k=n}^{\infty} T(k) z^k \right\rvert \le 2^{-d}.$

We can rewrite the hypergeometric term ratio as

$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

$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

$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

$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

$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)$$:

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

$|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¶

void hypgeom_init(hypgeom_t hyp)
void hypgeom_clear(hypgeom_t hyp)

## Error bounding¶

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.

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

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

void hypgeom_precompute(hypgeom_t hyp)

Precomputes the bounds data $$C$$, $$D$$, $$K$$ and an upper bound for $$T(K)$$.

## Summation¶

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.

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.