hypgeom.h – support for hypergeometric series¶
This module provides functions for highprecision evaluation of series of the form
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 highprecision 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(k1), \tilde Q(k) = Q(k) A(k1) 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
for given polynomials
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
We can rewrite the hypergeometric term ratio as
where
and where \(\tilde a_i = a_{pi} / a_p\), \(\tilde b_i = b_{qi} / b_q\). Next, we define
Now, if \(k > C\), the magnitude of the numerator of \(F(k)\) is bounded from above by
and if \(k > 2D\), the magnitude of the denominator of \(F(k)\) is bounded from below by
Putting the inequalities together gives the following bound, valid for \(k > K = \max(C, 2D)\):
Let \(r = qp\) and \(\tilde z = z a_p / b_q\). Assuming \(k > \max(C, 2D, {\tilde z}^{1/r})\), we have
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 floatingpoint 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.
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 firstorder 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(1z)  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 Wfunction 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(kD)}{(kC)(k2D)}.\]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.
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}^{n1} 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.