Sums of hypergeometric series with rational parameters

Many special functions and constants reduce to partial sums

with small rational parameters $a_i = p_i/q_i$, $b_j = u_j/v_j$ and a real argument $w$ of full precision: the error function ($r = 0$, $b = (3/2)$), the sine and cosine integrals (${}_1F_2$ and ${}_2F_3$ series and asymptotic series $\sum_n n!\,(i/z)^n$), Bessel functions of rational order, and the incomplete gamma function. The functions arb_hypgeom_sum_fmpq_arb and arb_hypgeom_sum_fmpq_imag_arb of arb_hypgeom compute such sums exactly as written (no truncation error), with $w = z$ or $w = 1/z$. Their usage is described in Hypergeometric sums with rational parameters. In each term ratio

the numerator and denominator are products of small integers: the cost of the sum is dominated by the operations that involve the full-precision $w$.

Forward recurrence#

arb_hypgeom_sum_fmpq_arb_forward keeps the partial sum as a fraction $R/D$ with an integer-valued ball $D$: with $c_n$ and $d_n$ the integer numerator and denominator of the term ratio,

and a single division $R/D$ at the end. Each step costs one full multiplication by $w$ and multiplications by word-size integers when the products of the parameters fit in a word (otherwise by fmpz). This is the fastest method for a handful of terms.

Rectangular splitting#

The sum is evaluated backwards as a polynomial in $w^m$ with $m \approx \sqrt N$, after computing $w, w^2, \ldots, w^m$:

where each block $B_k$ is a polynomial of degree less than $m$ in $w$ with rational coefficients. Within a block, up to 8 consecutive term ratios are combined over a common denominator, so that the block is computed as one dot product of the powers of $w$ with exact integer coefficients (arb_dot_fmpz, with single-word products when the parameters are small) followed by one division by an integer; the multiplication by $w^m$ between blocks is the only full multiplication per block.

The cost is about $N/m + m$ full multiplications and $O(N)$ multiplications and divisions of balls by small integers, which are linear in $p$. When the terms decrease (for example $r < s$ and $|w|$ moderate), the blocks with small terms are computed at a precision reduced by the logarithm of their magnitude relative to the largest term, estimated in double precision from $\log_2|w|$ and the growth of $(a)_n/(b)_n$. This method is due to Smith [Smi1989] and is analysed with these refinements in [Joh2014b; Joh2014c].

Binary splitting#

When $w = P/Q$ itself has few bits, the matrix product form of Summation algorithms is computed as a product tree of exact (integer-valued) balls; arb_hypgeom_sum_fmpq_arb_bs puts $z$ into the numerators, or into the denominators if reciprocal is set. The cost is quasi-linear in the size of the result [Bre1976; CC1988]. The default choice uses binary splitting from 8192 bits when $z$ has at most $p/1000$ bits.

Real and imaginary parts#

arb_hypgeom_sum_fmpq_imag_arb computes $S_N(iw)$ with real arithmetic: the even and odd terms form the real and imaginary parts, with signs $i^{2k} = (-1)^k$ and $i^{2k+1} = (-1)^k i$. The forward and rectangular splitting versions track the two parts separately (in rectangular splitting, the dot products take every second power of $w$ and coefficient). This saves the complex multiplications that acb arithmetic would perform. Binary splitting, where the leaves are cheap, uses complex balls with the factor $i$ in the leaf matrices.

Use in FLINT#

• arb_hypgeom_sum_fmpq_arb and arb_hypgeom_sum_fmpq_imag_arb with their _forward, _rs and _bs variants.
• arb_hypgeom_si and arb_hypgeom_ci sum their ${}_1F_2$ and ${}_2F_3$ series with arb_hypgeom_sum_fmpq_arb and their asymptotic series with arb_hypgeom_sum_fmpq_imag_arb.
• The special cases $\sum z^k/(a)_k$ and $\sum (a)_k z^{-k}$ with $a = p/q$ are implemented separately by _arb_hypgeom_gamma_lower_sum_rs_1 and _arb_hypgeom_gamma_upper_sum_rs_1 for the error function.
• Multiple-precision floating-point versions of hypergeometric binary splitting for constants are in mp_real (see Floating-point arithmetic in FLINT).