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.