Automatic hypergeometric series

October 8, 2012

Arb now contains code for evaluating infinite (rapidly converging) hypergeometric series, i.e. sums of the form $\sum_{k=0}^{\infty} T(k)$ where $T(k) = R(k) T(k-1)$ where $R(k)$ is a rational function (currently limited to having rational coefficients).

A series is defined by providing input polynomials $A, B, P, Q$ such that $T(k) = (A(k) / B(k)) \tilde T(k)$ where $\tilde T(k) = (P(k) / Q(k)) \tilde T(k-1)$. This is often less elegant than $\,_pF_q$ notation, but more versatile since it’s possible to stick with rational numbers even when the polynomial roots (essentially the parameters of the $\,_pF_q$) become (complex) algebraic numbers. $A$ and $B$ are technically redundant, but keeping them separate improves efficiency.

No other information is needed — the algorithm automatically bounds the magnitude of the terms, chooses a nearly optimal truncation point to match the target precision, evaluates the sum efficiently using binary splitting, and adds a provably correct (modulo bugs) error bound to the output.

The code can be used for rigorous numerical evaluation of hypergeometric functions (erf, Bessel functions, etc.). It actually probably runs slightly slower than mpmath at low precision right now, but should be faster when I’ve optimized it. At the moment, the selling point is that it’s very efficient for high-precision evaluation of hypergeometric functions at fixed rational points, enabling easy computation of various mathematical constants.

For example, the Chudnovsky formula

$$\frac{1}{\pi} = 12 \sum^\infty_{k=0} \frac{(-1)^k (6k)! (13591409 + 545140134k)}{(3k)!(k!)^3 640320^{3k + 3/2}}$$

now has the following elegant implementation:

void
fmprb_const_pi_chudnovsky(fmprb_t s, long prec)
{
    hypgeom_t series;
    fmprb_t t, u;

    fmprb_init(t);
    fmprb_init(u);
    hypgeom_init(series);

    fmpz_poly_set_str(series->A, "2  13591409 545140134");
    fmpz_poly_set_str(series->B, "1  1");
    fmpz_poly_set_str(series->P, "4  5 -46 108 -72");
    fmpz_poly_set_str(series->Q, "4  0 0 0 10939058860032000");

    prec += FLINT_CLOG2(prec);
    fmprb_hypgeom_infsum(s, t, series, prec, prec);
    fmprb_sqrt_ui(u, 640320, prec);
    fmprb_mul_ui(u, u, 640320 / 12, prec);
    fmprb_mul(u, u, t, prec);
    fmprb_div(s, u, s, prec);

    hypgeom_clear(series);
    fmprb_clear(t);
    fmprb_clear(u);
}

Adding a few guard bits to the precision is not required for correctness, but gives nicer output (potentially avoiding unnecessary re-computations). As usual, the prec parameter is just intended to give approximate control of the error, and code at a higher level should verify that the radius of the output is small enough for the intended application.

In this case, with prec set to $\lceil 10^6 \log_2(10) \rceil = 3321929$, the algorithm decides to truncate the series before term $k = 70515$, and the final output is $3.14159… \pm 4.7005 \times 10^{-1000006}$, where most of the error comes from the final arithmetic operations.

A few more examples (omitting the boilerplate):

$$e = \sum_{k=0}^{\infty} \frac{1}{k!}$$

    fmpz_poly_set_str(series->A, "1  1");
    fmpz_poly_set_str(series->B, "1  1");
    fmpz_poly_set_str(series->P, "1  1");
    fmpz_poly_set_str(series->Q, "2  0 1");
    fmprb_hypgeom_infsum(s, t, series, prec, prec);
    fmprb_div(s, s, t, prec);

$$\log 2 = \frac{3}{4} \sum_{k=0}^{\infty} \frac{(-1)^k (k!)^2}{2^k (2k+1)!}$$

    fmpz_poly_set_str(series->A, "1  1");
    fmpz_poly_set_str(series->B, "1  1");
    fmpz_poly_set_str(series->P, "2  0 -1");
    fmpz_poly_set_str(series->Q, "2  4 8");
    fmprb_hypgeom_infsum(s, t, series, prec, prec);
    fmprb_mul_ui(s, s, 3, prec);
    fmprb_mul_ui(t, t, 4, prec);
    fmprb_div(s, s, t, prec);

$$\zeta(3) = \frac{1}{64} \sum_{k=0}^\infty \frac{(-1)^k (205k^2 + 250k + 77) (k!)^{10}}{((2k+1)!)^5}$$

    fmpz_poly_set_str(series->A, "3  77 250 205");
    fmpz_poly_set_str(series->B, "1  1");
    fmpz_poly_set_str(series->P, "6  0 0 0 0 0 -1");
    fmpz_poly_set_str(series->Q, "6  32 320 1280 2560 2560 1024");
    fmprb_hypgeom_infsum(s, t, series, prec, prec);
    fmprb_mul_ui(t ,t, 64, prec);
    fmprb_div(s, s, t, prec);

Catalan’s constant:

$$C = \sum_{k=0}^{\infty} \frac{(-1)^k 4^{4 k+1} \left(40 k^2+56 k+19\right) ((k+1)!)^2 ((2 (k+1))!)^3}{(k+1)^3 (2 k+1) ((4
(k+1))!)^2}$$

    fmpz_poly_set_str(series->A, "3  19 56 40");
    fmpz_poly_set_str(series->B, "1  1");
    fmpz_poly_set_str(series->P, "5  0 0 0 32 -64");
    fmpz_poly_set_str(series->Q, "5  9 96 352 512 256");
    fmprb_hypgeom_infsum(s, t, series, prec, prec);
    fmprb_mul_ui(t, t, 18, prec);
    fmprb_div(s, s, t, prec);

And finally, a random series:

$$
r = \sum_{k=0}^{\infty} = \frac{\left(-k^2+7 k+3\right) (2 k)! (k+4)!}{\left(9 k^3+8 k^2+5 k+2\right) k! (k+1)! (k+2)! (k+3)!
\binom{2 k}{k}^2} = 3.079715\ldots
$$

    fmpz_poly_set_str(series->A, "3  3 7 -1");
    fmpz_poly_set_str(series->B, "4  2 5 8 9");
    fmpz_poly_set_str(series->P, "4  0 0 4 1");
    fmpz_poly_set_str(series->Q, "5  -12 2 32 22 4");
    fmprb_hypgeom_infsum(s, t, series, prec, prec);
    fmprb_div(s, s, t, prec);
    fmprb_mul_ui(s, s, 2, prec);

The following table shows the results of some benchmarking (done on a 32-bit system). I’ve compared the code for $\pi$ and $\zeta(3)$ with the old code, which did the error bounding and sum evaluation manually for each constant. We find that there is no slowdown (in fact $\pi$ has gotten very slightly faster). I’ve also included timings for the relevant functions in MPFR (only $\log 2$ is directly comparable, as MPFR uses a different algorithm for all the others).

Constant	Digits	Old	New	MPFR
$e$	10⁶		0.51 s	0.84 s
	10⁷		6.9 s	11 s
$\pi$	10⁶	1.4 s	1.2 s	3.8 s
	10⁷	22 s	22 s	68 s
$\log 2$	10⁶		4.9 s	5.1 s
	10⁷		84 s	91 s
$\zeta(3)$	10⁶	6.0	6.0 s
	10⁷	113 s	113 s
$C$	10⁶		28 s	37 s
	10⁷		524 s	648 s
$r$	10⁶		5.0 s
	10⁷		81 s

The part that isn’t yet automatic is converting an explicit product of factorials, binomial coefficients, rising factorials, powers, and polynomials to a recurrence equation in the right input format. This is no fun to do by hand except in simple cases, and I plan to add code later on to support such conversions in an algorithmic way. (It’s easy to do with most general-purpose computer algebra systems.) I also plan to add more optimizations (e.g. algorithmically choosing the most efficient representation of the recurrence, instead of relying on the user to supply it that way), and supporting holonomic sums in a similar fashion. I will also soon add support for evaluating $\sum_{k=0}^{\infty} T(k) z^k$ where $z$ is a real or complex floating-point ball (perhaps even a matrix?).

fredrikj.net | Blog index | RSS feed | Follow me on Mastodon | Become a sponsor