Safe Haskell | Safe-Inferred |
---|---|

Language | Haskell2010 |

## 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} 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{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).

## Synopsis

- data Hypgeom = Hypgeom !(ForeignPtr CHypgeom)
- type CHypgeom = CFlint Hypgeom
- newHypgeom :: IO Hypgeom
- withHypgeom :: Hypgeom -> (Ptr CHypgeom -> IO a) -> IO (Hypgeom, a)
- withNewHypgeom :: (Ptr CHypgeom -> IO a) -> IO (Hypgeom, a)
- hypgeom_init :: Ptr CHypgeom -> IO ()
- hypgeom_clear :: Ptr CHypgeom -> IO ()
- hypgeom_estimate_terms :: Ptr CMag -> CInt -> CLong -> IO CLong
- hypgeom_bound :: Ptr CMag -> CInt -> CLong -> CLong -> CLong -> Ptr CMag -> Ptr CMag -> CLong -> IO CLong
- hypgeom_precompute :: Ptr CHypgeom -> IO ()
- arb_hypgeom_sum :: Ptr CArb -> Ptr CArb -> Ptr CHypgeom -> CLong -> CLong -> IO ()
- arb_hypgeom_infsum :: Ptr CArb -> Ptr CArb -> Ptr CHypgeom -> CLong -> CLong -> IO ()

# Support for hypergeometric series

# Types

#### Instances

Storable CHypgeom Source # | |

Defined in Data.Number.Flint.Hypgeom.FFI |

newHypgeom :: IO Hypgeom Source #

# Memory management

# Error bounding

hypgeom_estimate_terms :: Ptr CMag -> CInt -> CLong -> IO CLong Source #

*hypgeom_estimate_terms* *z* *r* *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.

hypgeom_bound :: Ptr CMag -> CInt -> CLong -> CLong -> CLong -> Ptr CMag -> Ptr CMag -> CLong -> IO CLong Source #

*hypgeom_bound* *error* *r* *C* *D* *K* *TK* *z* *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.

hypgeom_precompute :: Ptr CHypgeom -> IO () Source #

*hypgeom_precompute* *hyp*

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

# Summation

arb_hypgeom_sum :: Ptr CArb -> Ptr CArb -> Ptr CHypgeom -> CLong -> CLong -> IO () Source #

*arb_hypgeom_sum* *P* *Q* *hyp* *n* *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.

arb_hypgeom_infsum :: Ptr CArb -> Ptr CArb -> Ptr CHypgeom -> CLong -> CLong -> IO () Source #

*arb_hypgeom_infsum* *P* *Q* *hyp* *tol* *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*.