Fredrik Johansson

Move along

2012-01-26T13:58:00.003+01:00

I'm moving my blog to http://fredrikj.net/blog. If you're living in the future and reading this, you should go there instead!

Some FLINT 2.2 highlights

2011-06-08T10:58:00.004+02:00

Version 2.2 of FLINT (Fast Library for Number Theory) was released last weekend. Some updated benchmarks are available.

In this blog post, I'm going to talk a bit about features I contributed in this version. With apologies to Sebastian Pancratz who wrote a whole lot of code as well -- in particular, a new module for computing with p-adic numbers, and a module for rational functions! Bill Hart also implemented a faster polynomial GCD, which is a quite important update since GCD is crucial for most polynomial business. Anyhow...

Polynomial interpolation

I've added various functions for polynomial interpolation to the fmpz_poly module. In general, these can be used to speed up various computations involving integer or rational polynomials by mapping a given problem to (Z/nZ)[x], Z or even Z/nZ, taking advantage of fast arithmetic in those rings, and then via interpolation recovering a result in Z[x] or Q[x].

Firstly, there are some new Chinese Remainder Theorem functions for integer polynomials, allowing you to reconstruct an integer polynomial from a bunch of polynomials with coefficients modulo different primes. Straightforward code (the actual work is done by functions in the fmpz module), but useful to have. The CRT functions are used by the new modular GCD code.

There are also functions for evaluating an integer polynomial on a set of points, and forming the interpolating polynomial (generally with rational coefficients) given a set of points and the values at those points.

Finally, user-friendly functions for evaluation and interpolation at a power of two (Kronecker segmentation) have been added. The code for this is actually a very old part of FLINT, and possibly some of the most complicated code in the library (packing bits efficiently is surprisingly hard). The new functions just wrap this functionality, but take care of memory management and various special cases, so you can now just safely do something like:


    fmpz_t z;
    fmpz_init(z);
    long bits = fmpz_poly_max_bits(poly) + 1; /* +1 for signs */
    fmpz_poly_bit_pack(z, poly, bits);
    fmpz_poly_bit_unpack(poly, z, bits);  /* recover poly */
    fmpz_clear(z);

Apart from Kronecker segmentation, these functions are not currently asymptotically fast. Fast multi-modulus CRT for coefficient reconstruction is probably not all that important in most circumstances, because it's more common to use evaluation-interpolation techniques for polynomials of large degree and small coefficients than the other way around. Nonetheless, polynomials with large coefficients do arise as well. For example, the vector Bernoulli number code in FLINT relies on fast CRT, and currently uses custom code for this.

Polynomial interpolation uses Lagrange interpolation with barycentric weights, with a few tricks to avoid fractions. This is all implemented using an O(n^2) algorithm, but the actual time complexity is higher due to the fact that the coefficients when working over integers usually will be large, around n! in magnitude.

Here are some timing examples, evaluating and recovering a length-n polynomial with +/-1 coefficients basically as follows:


    x = _fmpz_vec_init(n);
    y = _fmpz_vec_init(n);
    fmpz_poly_init(P);
    fmpz_poly_init(Q);

    for (i = 0; i < n; i++)
        x[i] = -n/2 + i;

    fmpz_poly_randtest(P, state, n, 1);

    fmpz_poly_evaluate_fmpz_vec(y, P, x, n);
    fmpz_poly_interpolate_fmpz_vec(Q, x, y, n);

The bits column below measures the largest value in y, which grows quite large despite the input polynomial having small coefficients:


    n=8  eval=762 ns  interp=13 us  bits=8  ok=1
    n=16  eval=3662 ns  interp=61 us  bits=42  ok=1
    n=32  eval=29 us  interp=673 us  bits=-113  ok=1
    n=64  eval=136 us  interp=4951 us  bits=-316  ok=1
    n=128  eval=625 us  interp=45 ms  bits=-762  ok=1
    n=256  eval=2500 us  interp=792 ms  bits=-1779  ok=1
    n=512  eval=12 ms  interp=10 s  bits=-4089  ok=1

As you can see, the interpolation speed is not too bad for small n, but eventually grows out of control. How to do better?

Naive Lagrange interpolation is not optimal: it is possible to do n-point evaluation and interpolation in essentially O(n log² n) operations. Such algorithms do not necessarily lead to an improvement over the integers (you still have to deal with coefficient explosion), but they should win over finite fields. So the right solution will perhaps be to add polynomial evaluation/interpolation functions based on modular arithmetic.

Rational numbers and matrices

A new module fmpq is provided for computing with arbitrary-precision rational numbers. For the user, the fmpq_t type essentially behaves identically to the MPIR mpq_t type. However, an fmpq_t only takes up two words of memory when the numerator and denominator are small (less than 2⁶²), whereas an mpq_t always requires six words plus additional heap-allocated space for the actual number data.

The fmpq functions are a bit faster than mpq functions in many cases when the numerator and/or denominator is small. But the main improvement should come for vectors, matrices or polynomials of rational numbers, due to the significantly reduced memory usage and memory management overhead (especially when many entries are zero or integer-valued).

Some higher-level functionality is also provided in the fmpq module, e.g. for rational reconstruction. The functions for computing special rational numbers (like Bernoulli numbers) have also been switched over to the fmpq type. Another supported feature is enumeration of the rationals (using the Calkin-Wilf sequence or by height). Generating the 100 million "first" positive rational numbers takes 9.6 seconds done in order of height, or 2.6 seconds in Calkin-Wilf order.

FLINT actually does not use fmpq's to represent polynomials over Q (fmpq_poly), and probably never will. The fmpq_poly module represents a polynomial over Q as an integer polynomial with a single common denominator, which is usually faster. The reason for adding the fmpq_t type is that it enabled developing the new fmpq_mat module, which implements dense matrices of rational numbers. For matrices, a common-denominator representation would be less convenient and in many cases completely impractical.

The new FLINT fmpq_mat module is very fast, or at least very non-slow. It is easy to find examples where it does a simple computation a thousand times faster than the rational matrices in Sage.

There's not actually much code in the fmpq_mat module itself; it does almost all "level 3" linear algebra (computations requiring matrix multiplication or Gaussian elimination) by clearing denominators and computing over the integers. This approach is in fact stolen shamelessly from Sage, but the functions in Sage are highly unoptimized in many cases. The code in Sage still wins for many sufficiently large problems as it has asymptotically fast algorithms for many things we do not (like computing null spaces). See the benchmarks page for more details.

I should not forget to mention that I've implemented Dixon's p-adic algorithm for solving Ax = b for nonsingular square A. (I wish I had a good link for Dixon's algorithm here, but sadly it doesn't appear to be described conveniently anywhere on the web. The original paper is "Exact solution of linear equations using P-adic expansions", if you have the means to get through the Springer paywall.)

This is now used both for solving over both Z and Q. The solver in FLINT is competitive with Sage (which uses IML+ATLAS) up to systems of dimension somewhere between perhaps 100 and 1000 (depending greatly on the size of the entries in the inputs and in the solution!). There's much to do here -- we should eventually have BLAS support in FLINT, which will speed up core matrix arithmetic, but there's room for a lot of algorithmic tuning as well.

There are some other minor new matrix features as well... they can be found in the changelog.

Polynomial matrices

A new module (fmpz_poly_mat) is provided for dense matrices over Z[x], i.e. matrices whose entries are polynomials with integer coefficients. The available functionality includes matrix multiplication, row reduction, and determinants. Matrix multiplication is particularly fast, as it uses the Kronecker segmentation interpolation/evaluation technique described above. (A similar algorithm is provided for determinants, but it's not really optimal as this point.)

The benchmarks page has detailed some detailed timings, so I won't repeat them here -- but generally speaking, the FLINT implementation is an order of magnitude faster than Sage or Magma for matrices of manageable size.

There's much more to be done for polynomial matrices. Row reduction is implemented quite efficiently, but it's too slow as an algorithm for many tasks such as computing null spaces of very large matrices. A future goal is to implement asymptotically fast algorithms (see the paper on x-adic lifting by Burçin Eröcal and Arne Storjohann for example).

100 mpmath one-liners for pi

2011-03-14T16:37:00.006+01:00

Since it's pi day today, I thought I'd share a list of mpmath one-liners for computing the value of pi to high precision using various representations in terms of special functions, infinite series, integrals, etc. Most of them can already be found as doctest examples in some form in the mpmath documentation.

A few of the formulas explicitly involve pi. Using those to calculate pi is rather circular (!), though a few of them could still be used for computing pi using numerical root-finding. In any case, most of the formulas are circular even when pi doesn't appear explicitly since mpmath is likely using its value internally. In any further case, the majority of the formulas are not efficient for computing pi to very high precision (at least as written). Still, ~50 digits is no problem. Enjoy!


from mpmath import *
mp.dps = 50; mp.pretty = True

+pi
180*degree
4*atan(1)
16*acot(5)-4*acot(239)
48*acot(49)+128*acot(57)-20*acot(239)+48*acot(110443)
chop(2*j*log((1-j)/(1+j)))
chop(-2j*asinh(1j))
chop(ci(-inf)/1j)
gamma(0.5)**2
beta(0.5,0.5)
(2/diff(erf, 0))**2
findroot(sin, 3)
findroot(cos, 1)*2
chop(-2j*lambertw(-pi/2))
besseljzero(0.5,1)
3*sqrt(3)/2/hyp2f1((-1,3),(1,3),1,1)
8/(hyp2f1(0.5,0.5,1,0.5)*gamma(0.75)/gamma(1.25))**2
4*(hyp1f2(1,1.5,1,1) / struvel(-0.5, 2))**2
1/meijerg([[],[]], [[0],[0.5]], 0)**2
(meijerg([[],[2]], [[1,1.5],[]], 1, 0.5) / erfc(1))**2
(1-e) / meijerg([[1],[0.5]], [[1],[0.5,0]], 1)
sqrt(psi(1,0.25)-8*catalan)
elliprc(1,2)*4
elliprg(0,1,1)*4
2*agm(1,0.5)*ellipk(0.75)
(gamma(0.75)*jtheta(3,0,exp(-pi)))**4
cbrt(gamma(0.25)**4*agm(1,sqrt(2))**2/8)
sqrt(6*zeta(2))
sqrt(6*(zeta(2,3)+5./4))
sqrt(zeta(2,(3,4))+8*catalan)
exp(-2*zeta(0,1,1))/2
sqrt(12*altzeta(2))
4*dirichlet(1,[0,1,0,-1])
2*catalan/dirichlet(-1,[0,1,0,-1],1)
exp(-dirichlet(0,[0,1,0,-1],1))*gamma(0.25)**2/(2*sqrt(2))
sqrt(7*zeta(3)/(4*diff(lerchphi, (-1,-2,1), (0,1,0))))
sqrt(-12*polylog(2,-1))
sqrt(6*log(2)**2+12*polylog(2,0.5))
chop(root(-81j*(polylog(3,root(1,3,1))+4*zeta(3)/9)/2,3))
2*clsin(1,1)+1
(3+sqrt(3)*sqrt(1+8*clcos(2,1)))/2
root(2,6)*sqrt(e)/(glaisher**6*barnesg(0.5)**4)
nsum(lambda k: 4*(-1)**(k+1)/(2*k-1), [1,inf])
nsum(lambda k: (3**k-1)/4**k*zeta(k+1), [1,inf])
nsum(lambda k: 8/(2*k-1)**2, [1,inf])**0.5
nsum(lambda k: 2*fac(k)/fac2(2*k+1), [0,inf])
nsum(lambda k: fac(k)**2/fac(2*k+1), [0,inf])*3*sqrt(3)/2
nsum(lambda k: fac(k)**2/(phi**(2*k+1)*fac(2*k+1)), [0,inf])*(5*sqrt(phi+2))/2
nsum(lambda k: (4/(8*k+1)-2/(8*k+4)-1/(8*k+5)-1/(8*k+6))/16**k, [0,inf])
2/nsum(lambda k: (-1)**k*(4*k+1)*(fac2(2*k-1)/fac2(2*k))**3, [0,inf])
nsum(lambda k: 72/(k*expm1(k*pi))-96/(k*expm1(2*pi*k))+24/(k*expm1(4*pi*k)), [1,inf])
1/nsum(lambda k: binomial(2*k,k)**3*(42*k+5)/2**(12*k+4), [0,inf])
4/nsum(lambda k: (-1)**k*(1123+21460*k)*fac2(2*k-1)*fac2(4*k-1)/(882**(2*k+1)*32**k*fac(k)**3), [0,inf])
9801/sqrt(8)/nsum(lambda k: fac(4*k)*(1103+26390*k)/(fac(k)**4*396**(4*k)), [0,inf])
426880*sqrt(10005)/nsum(lambda k: (-1)**k*fac(6*k)*(13591409+545140134*k)/(fac(k)**3*fac(3*k)*(640320**3)**k), [0,inf])
4/nsum(lambda k: (6*k+1)*rf(0.5,k)**3/(4**k*fac(k)**3), [0,inf])
(ln(8)+sqrt(48*nsum(lambda m,n: (-1)**(m+n)/(m**2+n**2), [1,inf],[1,inf]) + 9*log(2)**2))/2
-nsum(lambda x,y: (-1)**(x+y)/(x**2+y**2), [-inf,inf], [-inf,inf], ignore=True)/ln2
2*nsum(lambda k: sin(k)/k, [1,inf])+1
quad(lambda x: 2/(x**2+1), [0,inf])
quad(lambda x: exp(-x**2), [-inf,inf])**2
2*quad(lambda x: sqrt(1-x**2), [-1,1])
chop(quad(lambda z: 1/(2j*z), [1,j,-1,-j,1]))
3*(4*log(2+sqrt(3))-quad(lambda x,y: 1/sqrt(1+x**2+y**2), [-1,1],[-1,1]))/2
sqrt(8*quad(lambda x,y: 1/(1-(x*y)**2), [0,1],[0,1]))
sqrt(6*quad(lambda x,y: 1/(1-x*y), [0,1],[0,1]))
sqrt(6*quad(lambda x: x/expm1(x), [0,inf]))
quad(lambda x: (16*x-16)/(x**4-2*x**3+4*x-4), [0,1])
quad(lambda x: sqrt(x-x**2), [0,0.25])*24+3*sqrt(3)/4
mpf(22)/7 - quad(lambda x: x**4*(1-x)**4/(1+x**2), [0,1])
mpf(355)/113 - quad(lambda x: x**8*(1-x)**8*(25+816*x**2)/(1+x**2), [0,1])/3164
2*quadosc(lambda x: sin(x)/x, [0,inf], omega=1)
40*quadosc(lambda x: sin(x)**6/x**6, [0,inf], omega=1)/11
e*quadosc(lambda x: cos(x)/(1+x**2), [-inf,inf], omega=1)
8*quadosc(lambda x: cos(x**2), [0,inf], zeros=lambda n: sqrt(n))**2
2*quadosc(lambda x: sin(exp(x)), [1,inf], zeros=ln)+2*si(e)
exp(2*quad(loggamma, [0,1]))/2
2*nprod(lambda k: sec(pi/2**k), [2,inf])
s=lambda k: sqrt(0.5+s(k-1)/2) if k else 0; 2/nprod(s, [1,inf])
s=lambda k: sqrt(2+s(k-1)) if k else 0; limit(lambda k: sqrt(2-s(k))*2**(k+1), inf)
2*nprod(lambda k: (2*k)**2/((2*k-1)*(2*k+1)), [1,inf])
2*nprod(lambda k: (4*k**2)/(4*k**2-1), [1, inf])
sqrt(6*ln(nprod(lambda k: exp(1/k**2), [1,inf])))
nprod(lambda k: (k**2-1)/(k**2+1), [2,inf])/csch(pi)
nprod(lambda k: (k**2-1)/(k**2+1), [2,inf])*sinh(pi)
nprod(lambda k: (k**4-1)/(k**4+1), [2, inf])*(cosh(sqrt(2)*pi)-cos(sqrt(2)*pi))/sinh(pi)
sinh(pi)/nprod(lambda k: (1-1/k**4), [2, inf])/4
sinh(pi)/nprod(lambda k: (1+1/k**2), [2, inf])/2
(exp(1+euler/2)/nprod(lambda n: (1+1/n)**n * exp(1/(2*n)-1), [1, inf]))**2/2
3*sqrt(2)*cosh(pi*sqrt(3)/2)**2*csch(pi*sqrt(2))/nprod(lambda k: (1+1/k+1/k**2)**2/(1+2/k+3/k**2), [1, inf])
2/e*nprod(lambda k: (1+2/k)**((-1)**(k+1)*k), [1,inf])
limit(lambda k: 16**k/(k*binomial(2*k,k)**2), inf)
limit(lambda x: 4*x*hyp1f2(0.5,1.5,1.5,-x**2), inf)
1/log(limit(lambda n: nprod(lambda k: pi/(2*atan(k)), [n,2*n]), inf),4)
limit(lambda k: 2**(4*k+1)*fac(k)**4/(2*k+1)/fac(2*k)**2, inf)
limit(lambda k: fac(k) / (sqrt(k)*(k/e)**k), inf)**2/2
limit(lambda k: (-(-1)**k*bernoulli(2*k)*2**(2*k-1)/fac(2*k))**(-1/(2*k)), inf)
limit(lambda k: besseljzero(1,k)/k, inf)
1/limit(lambda x: airyai(x)*2*x**0.25*exp(2*x**1.5/3), inf, exp=True)**2
1/limit(lambda x: airybi(x)*x**0.25*exp(-2*x**1.5/3), inf, exp=True)**2

A FLINT example: Lambert W function power series

2011-03-11T13:28:00.010+01:00

Two days ago, a new version of the Fast Library for Number Theory (FLINT) was released. I contributed a lot of new code to this release, including linear algebra speed improvements and new functionality for fast power series arithmetic and computation of special numbers and polynomials (see the release announcement and some of my benchmarking results).

In this blog post I'll demonstrate how to do power series arithmetic with FLINT, using its fmpq_poly module which implements polynomials over the rational numbers Q. Standard operations (addition, multiplication and division) were available before; the functions I've added include square root, log, exp, sin, tan, atan, etc. (all the usual elementary functions). The same functions are also available for power series over a finite field Z/nZ (with word-size n). Everything is asymptotically fast (the running time is linear in the size of the output, up to logarithmic factors).

Of course, transcendental functions are a bit restricted when considered over Q or Z/nZ, since it's only possible to obtain power series expansions at specific rational points (in most cases just x = 0). So at present, some very interesting numerical applications of fast power series arithmetic are not supported. But some time in the future, we'll probably add support for numerical power series over the reals and complexes as well.

As today's example, let us implement the Lambert W function for the power series ring Q[[x]]. The Lambert W function is defined implicitly by the equation x = W(x) exp(W(z)), which can be solved using Newton iteration with the update step w = w - (w exp(w) - x) / ((w+1) exp(w)).

Power series Newton iteration is just like numerical Newton iteration, except that the convergence behavior is much simpler: starting with a correct first-order expansion, each iteration at least doubles the number of correct coefficients.

A simple recursive implementation with asymptotically optimal performance (up to constant factors) looks as follows:


#include <stdio.h>
#include "flint.h"
#include "fmpq_poly.h"

void lambertw(fmpq_poly_t w, fmpq_poly_t x, long n)
{
    if (n == 1)
    {
        fmpq_poly_zero(w);
    }
    else
    {
        fmpq_poly_t t, u, v;

        lambertw(w, x, (n + 1) / 2);

        fmpq_poly_init(t);
        fmpq_poly_init(u);
        fmpq_poly_init(v);

        fmpq_poly_exp_series(t, w, n);
        fmpq_poly_mullow(u, t, w, n);
        fmpq_poly_sub(v, u, x);
        fmpq_poly_add(t, u, t);
        fmpq_poly_div_series(u, v, t, n);
        fmpq_poly_sub(w, w, u);

        fmpq_poly_clear(t);
        fmpq_poly_clear(u);
        fmpq_poly_clear(v);
    }
}

Beyond the base case W(x) = 0 + O(x), the function just computes w to accuracy ceil(n/2), and then extends it to accuracy n using a single Newton step. As we can see, C code directly using the FLINT library interface gets a bit verbose, but this style has the advantage of giving precise control over temporary memory allocation, polynomial lengths, etc. (it is very similar to the interface of GMP/MPIR).

We add a simple test main routine:


int main()
{
   fmpq_poly_t x;
   fmpq_poly_t w;

   fmpq_poly_init(x);
   fmpq_poly_init(w);

   fmpq_poly_set_coeff_ui(x, 1, 1);
   lambertw(w, x, 10);

   fmpq_poly_print_pretty(w, "x");
   printf("\n");

   fmpq_poly_clear(x);
   fmpq_poly_clear(w);
}

The output of the program is:


531441/4480*x^9 - 16384/315*x^8 + 16807/720*x^7 - 54/5*x^6 + 125/24*x^5 - 8/3*x^4 + 3/2*x^3 - 1*x^2 + 1*x

It is well known that the coefficients in this series are given in closed form by (-k)^k-1 / k!, so we can check that the output is correct.

Computing 1000 terms takes just a few seconds. If this sounds like much, remember that the coefficients grow rapidly: together, the computed numerators and denominators have over 2 million digits!

So far this is perhaps not so interesting, as we could compute the coefficients faster using a direct formula. But the nice thing is that arbitrary compositions are allowed, i.e we can compute W(f(x)) for any given power series f, and this will still be just as fast.

Let's consider a nontrivial example: the infinite "power tower" T(z) = z^{z^{z^{z^{.^.}}}}. A moment's reflection shows that this is an analytic function with a rational power series expansion around z = 1. In fact, we have explicitly T(z) = W(-log(z))/(-log(z)). We can compute this series expansion (in the shifted variable x = z - 1) as follows:


int main()
{
    fmpq_poly_t x;
    fmpq_poly_t w;

    long n = 10;

    fmpq_poly_init(x);
    fmpq_poly_init(w);

    fmpq_poly_set_coeff_ui(x, 0, 1);
    fmpq_poly_set_coeff_ui(x, 1, 1);
    fmpq_poly_log_series(x, x, n + 1);
    fmpq_poly_neg(x, x);
    lambertw(w, x, n + 1);
    fmpq_poly_shift_right(w, w, 1);
    fmpq_poly_shift_right(x, x, 1);
    fmpq_poly_div_series(w, w, x, n);

    fmpq_poly_print_pretty(w, "x");
    printf("\n");

    fmpq_poly_clear(x);
    fmpq_poly_clear(w);
}

The only complication is that fmpq_poly_div_series requires a nonzero leading coefficient in the denominator, so we must shift both series down one power.

The program outputs:


118001/2520*x^9 + 123101/5040*x^8 + 4681/360*x^7 + 283/40*x^6 + 4*x^5 + 7/3*x^4 + 3/2*x^3 + 1*x^2 + 1*x + 1

To make things nicer, we assume that the coefficients have the form a_k / k! (i.e. that T(z) is the exponential generating function for a_k) and change the output code to something like the following:


    long k;
    mpq_t t;
    mpz_t u;
    mpq_init(t);
    mpz_init(u);

    for (k = 0; k < n; k++)
    {
        fmpq_poly_get_coeff_mpq(t, w, k);
        mpz_fac_ui(u, k);
        mpz_mul(mpq_numref(t), mpq_numref(t), u);
        mpq_canonicalize(t);
        gmp_printf("%Qd ", t);
    }

    mpq_clear(t);
    mpz_clear(u);

This indeed gives us an integer sequence:


1 1 2 9 56 480 5094 65534 984808 16992144

Now what is the value of the 1000th coefficient (to be precise, a₁₀₀₀, the initial one being the 0th!) in this sequence? After a simple modification of the program, 2.9 seconds of computation gives:


11608872341636087705816513947297167830568265588875720061704
01832880235304267566817912141661469362953389062004053809005
65797054717998071778437757582562676432270594729770831984037
17901116787718293231769568392734610884078152929278291961741
58010228897635319982035567487202368704727403137478203768363
54056589570878404139562784693762331122998711070595645913436
44753733499423283972136827590268687580725109528808039530647
10910254098110789162443473433433758060122558659258182027755
69656436509351036076228649393400187670469063215003559774586
49501015173633083100668758800804388616363320813332492596835
40185987183963214465225072970422690115905543500507650640978
08856685726892919091844572545581642428942983342505179168857
61923031601434642410137173087273453449219217659949560840949
29145910407919393564145312029717057693032572341514569188719
42207889248610196901459400077483577940763454422516589494589
38697976290832628091067571489853751119661925805775760182956
07151657547554699411688610841404991952520564137242651305186
19966880917401902668151574186675809680229260294868082194497
63338464294487320831362657576767926588975644587806316363928
21662453081804476234328933125206970873131871382852201414093
31942812710129867491990841736391939490562342870154316209797
95555638177793757660689621198594912024704112203014400855204
04879191040818216462884689447945725483793082854991264186114
00713712447555062853630274495412279277142852027491666742488
18689076794537156576609645279481454870296442864829766014978
76385015229773871193575960430599394232421616401025152808967
97542967829629757402705726445239053261557399630212654678115
91948563122399554735529747742515102962530483866618795187470
92568029262248891738821070847168914030430887617489382116571
31479578425767585519331805968937010542495567221591600504522
70151935685333213987251220404383044513120115761331175072544
91881860724844683157343078083901966247367831930705346651165
57731933519958498663270193078704185994119446629783305199163
25824443621182783667024174595493553934149891052564101562124
66082538519787858297949190033471879555319648142879656530503
22140399695072998272983889906823049155302053273484019653833
08158019685729676988160041114485564188896445502120959889736
26684734069125268167350474483728161637188322446040542612820
83620649731423678182582137133666912162187578149277916758677
65932622140692260754343559763758688544180440952477345437585
88260535486569816885029406514351482276962081562798684604230
27051552771077659399889469617306015354335528530235916712574
33756257973655927835185354982512983428012895270181767297060
61394636504681554763302758450669487653360858511886083023090
56603401440047692698200295529572915618836122163118770906896
63441094011689868848158568518095899683719854486361541380832
18026233272569661209672552513531416295218659379214599386577
71439492527626159018195922050167504883881038997644963556212
95634222871269535245013411241216112695705600000000000000000
0000000000000000000000000

In fact, if we look up the first 10 coefficients in the On-Line Encyclopedia of Integer Sequences, we find http://oeis.org/A033917. This OEIS entry lists the representation

a(n) = Sum_{k=0..n} Stirling1(n, k)*(k+1)^(k-1)

Since FLINT supports fast vector computation of Stirling numbers, this formula can be implemented efficiently:


#include "fmpz.h"
#include "fmpz_vec.h"
#include "arith.h"

void coefficient(fmpz_t a, long n)
{
    long k;
    fmpz * s;
    fmpz_t t;

    s = _fmpz_vec_init(n + 1);
    fmpz_stirling1_vec(s, n, n + 1);

    fmpz_init(t);
    fmpz_zero(a);
    for (k = 1; k <= n; k++)
    {
        fmpz_set_ui(t, k + 1);
        fmpz_pow_ui(t, t, k - 1);
        fmpz_addmul(a, s + k, t);
    }

    _fmpz_vec_clear(s, n + 1);
    fmpz_clear(t);
}

int main()
{
    fmpz_t a;
    fmpz_init(a);
    coefficient(a, 1000);
    fmpz_print(a);
    printf("\n");
    fmpz_clear(a);
}

And indeed, the output turns out to be the same!

This program is faster, taking only 0.1 seconds to run. But of course, it only gives us a single coefficient, and would be slower for computing a range of values by making repeated calls.

Similar ideas to those presented here (basically, reducing a problem to fast polynomial multiplication using generating functions, Newton iteration, etc.) are used internally by FLINT for computation of the standard elementary functions themselves as well as various special numbers and polynomials (Bernoulli numbers and polynomials, partitions, Stirling numbers, Bell numbers, etc). The internal code uses a lot of tricks to reduce overhead and handle special cases faster, however. (See the previous blog post Fast combinatorial and number-theoretic functions with FLINT 2, and for more recent information the release announcement and benchmarks page linked at the top of this post.)

In other news, I haven't written a lot of code for mpmath or Sage recently. Of course, my hope is that FLINT (2) will make it into Sage in the not too distant future. The fast polynomial and power series arithmetic support in FLINT will also be very useful for future special functions applications (in mpmath and elsewhere).

Announcing mpmath 0.16

2010-09-24T17:27:00.003+02:00

I'm happy to announce the release of mpmath 0.16, which contains the usual bugfixes as well as a slew of new features!

The main focus has been to improve coverage of special functions. Additions include inhomogeneous Bessel functions, Bessel function zeros, incomplete elliptic integrals, and parabolic cylinder functions. As of 0.16, mpmath implements essentially everything listed in the NIST Digital Library of Mathematical Functions chapters 1-20, as well as 21,24,27 and 33. (For 25 and 26 -- combinatorial and number-theoretic functions, see also my post about FLINT 2.)

Another major change is that mpmath 0.16 running in Sage will be much faster thanks to new extension code (currently awaiting review for inclusion in Sage). I've clocked speedups between 1.3x and 2x for various nontrivial pieces of code (such as the mpmath test suite and the torture test programs).

Thanks to William Stein, my work on mpmath during the summer was funded using resources from NSF grant DMS-0757627. This support is gratefully acknowledged.

Most of the new features are described in previous posts on this blog. For convenience, here is a short summary:

Assorted special functions update

The documentation now includes plots to illustrate several of the special functions.

Airy functions have been rewritten for improved speed and accuracy and to support evaluation of derivatives.

Functions airyaizero(), airybizero() for computation of Airy function zeros have been implemented.

Inhomogeneous Airy (Scorer) functions scorergi() and scorerhi() have been implemented.

Four inhomogeneous Bessel functions have been added (lommels1(), lommels2(), angerj(), webere()).

The Lambert W function has been rewritten to fix various bugs and numerical issues

Incomplete elliptic integrals complete

The Legendre and Carlson incomplete elliptic integrals for real and complex arguments have been implemented (ellipf(), ellipe(), ellippi(), elliprf(), elliprc(), elliprj(), elliprd(), elliprg()).

Sage Days 23, and Bessel function zeros

Functions besseljzero() and besselyzero() have been implemented for computing the m-th zero of J_ν(z), J'_ν(z) Y_ν(z), or Y'_ν(z) for any positive integer index m and real order ν ≥ 0.

Post Sage Days 24 report

The Parabolic cylinder functions pcfd(), pcfu(), pcfv(), pcfw() have been implemented.

Euler-Maclaurin summation of hypergeometric series

Hypergeometric functions _pF_p-1(...; ...; z) now support accurate evaluation close to the singularity at z = 1.

A function sumap() has been added for summation of infinite series using the Abel-Plana formula.

Functions diffs_prod() and diffs_prod() have been added for generating high-order derivatives of products or exponentials of functions with known derivatives.

Again, mpmath in Sage is about to get faster

New Cython extension code has been written for Sage to speed up various operations in mpmath, including elementary functions and hypergeometric series.

There are various other changes as well, such as support for matrix slice indexing (contributed by Ioannis Tziakos -- thanks!). As usual, details are available in the changelog and the Changes page on the Google Code project site.

Again, mpmath in Sage is about to get faster

2010-09-22T15:20:00.006+02:00

My summer project on special functions in mpmath and Sage, generously supported by William Stein with funds from NSF grant DMS-0757627, is nearing completion. I will soon release mpmath-0.16, which contains lots of new special functions and bugfixes. Sage users will also benefit from ~1500 lines of new Cython code (preliminary patch here) that speeds up various basic operations. Executing mpmath.runtests() in Sage on my laptop now takes 10.47 seconds (8.60 from a warm cache), compared to 14.21 (11.84) seconds with the new extensions disabled -- a global speedup of 30%.

For comparison, pure-Python mpmath with gmpy as the backend takes 21.46 (18.72) seconds to execute the unit tests and pure-Python mpmath with the pure-Python backend takes 52.33 (45.92) seconds.

Specifically, the new extension code implements exp for real and complex arguments, cos, sin and ln for real arguments, complex exponentiation in some cases, and summation of hypergeometric series, entirely in Cython.

Timings before (new extensions disabled):


sage: import mpmath
sage: x = mpmath.mpf(0.37)
sage: y = mpmath.mpf(0.49)
sage: %timeit mpmath.exp(x)
625 loops, best of 3: 14.5 µs per loop
sage: %timeit mpmath.ln(x)
625 loops, best of 3: 23.2 µs per loop
sage: %timeit mpmath.cos(x)
625 loops, best of 3: 17.2 µs per loop
sage: %timeit x ^ y
625 loops, best of 3: 39.9 µs per loop
sage: %timeit mpmath.hyp1f1(2r,3r,4r)
625 loops, best of 3: 90.3 µs per loop
sage: %timeit mpmath.hyp1f1(x,y,x)
625 loops, best of 3: 83.6 µs per loop
sage: %timeit mpmath.hyp1f1(x,y,mpmath.mpc(x,y))
625 loops, best of 3: 136 µs per loop

Timings after (new extensions enabled):


sage: import mpmath
sage: x = mpmath.mpf(0.37)
sage: y = mpmath.mpf(0.49)
sage: %timeit mpmath.exp(x)
625 loops, best of 3: 2.72 µs per loop
sage: %timeit mpmath.ln(x)
625 loops, best of 3: 7.25 µs per loop
sage: %timeit mpmath.cos(x)
625 loops, best of 3: 4.13 µs per loop
sage: %timeit x ^ y
625 loops, best of 3: 10.5 µs per loop
sage: %timeit mpmath.hyp1f1(2r,3r,4r)
625 loops, best of 3: 47.1 µs per loop
sage: %timeit mpmath.hyp1f1(x,y,x)
625 loops, best of 3: 59.4 µs per loop
sage: %timeit mpmath.hyp1f1(x,y,mpmath.mpc(x,y))
625 loops, best of 3: 83.1 µs per loop

The new elementary functions use a combination of custom algorithms and straightforward MPFR wrappers. Why not just wrap MPFR for everything? There are two primary reasons:

Firstly, because MPFR numbers have a limited range, custom code still needs to be used in the overflowing cases, and this is almost as much work as an implementation-from-scratch. (There are also some more minor incompatibilities, like lack of round-away-from-zero in MPFR, that result in a lot of extra work.)

Secondly, MPFR is not always fast (or as fast as it could be), so it pays off to write custom code. In fact, some of the ordinary Python implementations of functions in mpmath are faster than their MPFR counterparts in various cases, although that is rather exceptional (atan is an example). But generally, at low-mid precisions, it is possible to be perhaps 2-4x faster than MPFR with carefully optimized C code (see fastfunlib). This is a longer-term goal.

Already now, with the new extension code, the mpmath exponential function becomes faster than the Sage RealNumber version (based on MPFR) at low precision:


sage: %timeit mpmath.exp(x)
625 loops, best of 3: 2.75 µs per loop
sage: w = RealField(53)(x)
sage: %timeit w.exp()
625 loops, best of 3: 5.57 µs per loop

As the timings above indicate, hypergeometric series have gotten up to 2x faster. The speedup of the actual summation is much larger, but much of that gain is lost in various Python overheads (more work can be done on this). There should be a noticeable speedup for some hypergeometric function computations, while others will not benefit as much, for the moment.

Another benchmark is the extratest_zeta.py script in mpmath, which exercises the mpmath implementation of the Riemann-Siegel formula for evaluation of ζ(s) for complex s with large imaginary part. Such computations largely depend on elementary function performance (cos, sin, exp, log).

Here are the new timings for mpmath in Sage:


fredrik@scv:~/sage$ ./sage /home/fredrik/mp/mpmath/tests/extratest_zeta.py
399999999 156762524.675 ok = True (time = 1.144)
241389216 97490234.2277 ok = True (time = 9.271)
526196239 202950727.691 ok = True (time = 1.671)
542964976 209039046.579 ok = True (time = 1.189)
1048449112 388858885.231 ok = True (time = 1.774)
1048449113 388858885.384 ok = True (time = 1.604)
1048449114 388858886.002 ok = True (time = 2.096)
1048449115 388858886.002 ok = True (time = 2.587)
1048449116 388858886.691 ok = True (time = 1.546)

This is mpmath in Sage with the new extension code disabled:


fredrik@scv:~/sage$ ./sage /home/fredrik/mp/mpmath/tests/extratest_zeta.py
399999999 156762524.675 ok = True (time = 2.352)
241389216 97490234.2277 ok = True (time = 14.088)
526196239 202950727.691 ok = True (time = 3.036)
542964976 209039046.579 ok = True (time = 2.104)
1048449112 388858885.231 ok = True (time = 3.707)
1048449113 388858885.384 ok = True (time = 3.283)
1048449114 388858886.002 ok = True (time = 4.444)
1048449115 388858886.002 ok = True (time = 5.592)
1048449116 388858886.691 ok = True (time = 3.101)

This is mpmath in ordinary Python mode, using gmpy:


fredrik@scv:~/sage$ python /home/fredrik/mp/mpmath/tests/extratest_zeta.py
399999999 156762524.675 ok = True (time = 2.741)
241389216 97490234.2277 ok = True (time = 13.842)
526196239 202950727.691 ok = True (time = 3.124)
542964976 209039046.579 ok = True (time = 2.143)
1048449112 388858885.231 ok = True (time = 3.257)
1048449113 388858885.384 ok = True (time = 2.912)
1048449114 388858886.002 ok = True (time = 3.953)
1048449115 388858886.002 ok = True (time = 4.964)
1048449116 388858886.691 ok = True (time = 2.762)

With the new extension code, it appears that zeta computations are up to about twice as fast. This speedup could be made much larger as there still is a significant amount of Python overhead left to remove -- also a project for the future.

Fast combinatorial and number-theoretic functions with FLINT 2

2010-09-05T16:50:00.004+02:00

Time for a development update! Recently, I've done only a limited amount of work on mpmath (I have a some almost-finished Cython code for sage.libs.mpmath and new code for numerical integration in mpmath, both to be committed fairly soon -- within a couple of weeks, hopefully).

The last few weeks, I've mostly been contributing to FLINT 2. For those unfamiliar with it, FLINT is a fast C library for computational number theory developed by Bill Hart and others (the other active developers right now are Sebastian Pancratz and Andy Novocin). In particular, FLINT implements ridiculously fast multiprecision integer vectors and polynomials. It also provides very fast primality testing and factorization for word-size integers (32 or 64 bits), among other things. FLINT 2 is an in-progress rewrite of FLINT 1.x, a current standard component in Sage.

What does this have to do with numerical evaluation of special functions (the usual theme of this blog)? In short, my goal is to add code to FLINT 2 for exact special function computations -- combinatorial and number-theoretic functions, special polynomials and the like. Such functions benefit tremendously from the fast integer and polynomial arithmetic available in FLINT 2.

All my code can be found in my public GitHub repository (the most recent commits as of this writing are in the 'factor' branch).

Functions I've implemented so far include:

Möbius μ and Euler φ (totient) functions for word-size and arbitrary-size integers

Divisor sum function σ_k for arbitrary-size integers

Ramanujan τ function (Δ-function q-expansion)

Harmonic numbers 1 + 1/2 + 1/3 + ... + 1/n

Primorials 2 · 3 · 5 · ... · p_n

Stirling numbers (1st and 2nd kind)

The versions in FLINT 2 of these functions should now be faster than all other implementations I've tried (GAP, Pari, Mathematica, the Sage library) for all ranges of arguments, except for those requiring factorization of large integers.

Some of these functions depend fundamentally on the ability to factorize integers efficiently. So far I've only implemented trial division for large integers in FLINT 2, with some clever code to extract large powers of small factors quickly. Sufficiently small cofactors are handled by calling Bill Hart's single-word factoring routines. The resulting code is very fast for "artificial" numbers like factorials, and will eventually be complemented with prime and perfect power detection code, plus fast implementations of Brent's algorithm and other methods. Later on the quadratic sieve from FLINT 1 will probably be ported to FLINT 2, so that FLINT 2 will be able to factor any reasonable number reasonably quickly.

Below, I've posted some benchmark results. A word of caution: all Mathematica timings were done on a different system, which is faster than my own laptop (typically by 30% or so). So in reality, Mathematica performs slightly worse relatively than indicated below. Everything else is timed on my laptop. I have not included test code for the FLINT2 functions (but it's just straightforward C code -- a function call or two between timeit_start and timeit_stop using FLINT 2's profiler module).

Möbius function (the following is basically a raw exercise of the small-integer factoring code):


Pari:
sage: %time pari('sum(n=1,10^6,moebius(n))');
CPU times: user 1.04 s, sys: 0.00 s, total: 1.04 s
Wall time: 1.04 s

Mathematica:
In[1]:= Timing[Sum[MoebiusMu[n], {n,1,10^6}];]
Out[1]= {0.71, Null}

flint2:
650 ms

Divisor sum:


Sage (uses Cython code):
sage: %time sigma(factorial(1000),1000);
CPU times: user 0.47 s, sys: 0.00 s, total: 0.47 s
Wall time: 0.46 s

Mathematica:
In[1]:= Timing[DivisorSigma[1000,1000!];]
Out[1]= {3.01, Null}

flint2:
350 ms

Ramanujan τ function:


Sage (uses FLINT 1):
sage: %time delta_qexp(100000);
CPU times: user 0.42 s, sys: 0.01 s, total: 0.43 s
Wall time: 0.42 s
sage: %time delta_qexp(1000000);
CPU times: user 6.02 s, sys: 0.37 s, total: 6.39 s
Wall time: 6.40 s

flint2:
100000: 230 ms
1000000: 4500 ms

An isolated value (Mathematica seems to be the only other software that knows how to compute this):


Mathematica:
In[1]:= Timing[RamanujanTau[10000!];]
Out[1]= {8.74, Null}

flint2:
280 ms

Harmonic numbers (again, only Mathematica seems to implement these). See also my old blog post How (not) to compute harmonic numbers. I've included the fastest version from there, harmonic5:


Mathematica:
In[1]:= Timing[HarmonicNumber[100000];]
Out[1]= {0.22, Null}
In[2]:= Timing[HarmonicNumber[1000000];]
Out[2]= {6.25, Null}
In[3]:= Timing[HarmonicNumber[10000000];]
Out[3]= {129.13, Null}

harmonic5: (100000):
100000: 0.471 s
1000000: 8.259 s
10000000: 143.639 s

flint2:
100000: 100 ms
1000000: 2560 ms
10000000: 49400 ms

The FLINT 2 function benefits from an improved algorithm that eliminates terms and reduces the size of the temporary numerators and denominators, as well as low-level optimization (the basecase summation directly uses the MPIR mpn interface).

Isolated Stirling numbers of the first kind:


Mathematica:
In[1]:= Timing[StirlingS1[1000,500];]
Out[1]= {0.24, Null}
In[2]:= Timing[StirlingS1[2000,1000];]
Out[2]= {1.79, Null}
In[3]:= Timing[StirlingS1[3000,1500];]
Out[3]= {5.13, Null}

flint 2:
100,500: 100 ms
2000,1000: 740 ms
3000,1500: 1520 ms

Isolated Stirling numbers of the second kind:


Mathematica:
In[1]:= Timing[StirlingS2[1000,500];]
Out11]= {0.21, Null}
In[2]:= Timing[StirlingS2[2000,1000];]
Out[2]= {1.54, Null}
In[3]:= Timing[StirlingS2[3000,1500];]
Out[3]= {4.55, Null}
In[4]:= Timing[StirlingS2[5000,2500];]
Out[4]= {29.25, Null}

flint2:
1000,500: 2 ms
2000,1000: 17 ms
3000,1500: 50 ms
5000,2500: 240 ms

In addition, fast functions are provided for computing a whole row or matrix of Stirling numbers. For example, computing the triangular matrix of ~1.1 million Stirling numbers of the first kind up to S(1500,1500) takes only 1.3 seconds. In Mathematica (again, on the faster system):


In[1]:= Timing[Table[StirlingS1[n,k], {n,0,1500}, {k,0,n}];]
Out[1]= {2.13, Null}

The benchmarks above mostly demonstrate performance for large inputs. Another nice aspect of the FLINT 2 functions is that there typically is very little overhead for small inputs. The high performance is due to a combination of algorithms, low-level optimization, and (most importantly) the fast underlying arithmetic in FLINT 2. I will perhaps write some more about the algorithms (for e.g. Stirling numbers) in a later post.

Euler-Maclaurin summation of hypergeometric series

2010-07-28T21:23:00.009+02:00

I recently fixed another corner case (it's always the corner cases that are hard!) in the hypergeometric function evaluation code in mpmath.

The difficulty in question concerns the functions _pF_p-1(...; ...; z), of which the Gauss hypergeometric function ₂F₁ is an important special case. These functions can be thought of as generalizations of the geometric series

and the natural logarithm

and share the property that the hypergeometric (Maclaurin) series has radius of convergence 1, with a singularity (either algebraic or logarithmic, as above) at z = 1.

Numerical evaluation is fairly easy for z in most of the complex plane. The series for _pF_p-1(...; ...; z) converges like the geometric series, so direct summation works well for, say, |z| ≤ r 0.9. Likewise, there is an analytic continuation formula which transforms z to 1/z, allowing evaluation for (say) |z| ≥ 1.1.

The remaining shell around the unit circle 0.9 < |z| < 1.1 is the difficult region. In fact, ₂F₁ can still (essentially) be computed using exact transformations here, but the higher cases remain difficult. As mentioned in an earlier post, mpmath handles this by calling nsum() which applies convergence acceleration to the slowly converging series. In particular, nsum implements the (generalized) Shanks transformation which is extremely efficient for this particular type of series – in fact, it works outside the radius of convergence (Shanks transformation effectively computes a Padé approximant), so it can be used anywhere in the difficult shell.

Still, the Shanks transformation is most efficient when arg(z) = ±π and unfortunately degenerates as arg(z) → 0, i.e. close to z = 1. nsum also implements Richardson extrapolation, which sometimes works at z = 1, but only for special parameter combinations in the hypergeometric series (although those special combinations happen to be quite important – they include the case where the hypergeometric series reduces to a polygamma function or polylogarithm).

Thus, we have the following evaluation strategy for _pF_p-1:

Blue: direct summation (with respect to z or 1/z)
Yellow: Shanks transformation
Red: ???

A commenter suggested trying some other convergence acceleration techniques (like the Levin transformation). But it seems that they still will fail in some cases.

Thus I've finally implemented the Euler-Maclaurin summation formula for the remaining case. The E-M formula differs from other convergence acceleration techniques in that it is based on analyticity of the integrand and does not require a particular rate of convergence. It is implemented in mpmath as sumem(), but requires a little work to apply efficiently.

To apply the E-M formula to a hypergeometric series Σ T(k), the term T(k) must first of all obviously be extended to an analytic function of k, which is done the straightforward way by considering the rising factorials as quotients of gamma functions a(a+1)...(a+k−1) = Γ(a+k) / Γ(a).

We are left with the difficulty of having to integrate T(k) and also to compute n-th order derivatives of this function for the tail expansion in the E-M formula (neither of which can be done in closed form). Fortunately, numerical integration with quad() works very well, but the derivatives remain a problem.

By default, sumem() computes the derivatives using finite differences. Although this works, it gets extremely expensive at high precision for hypergeometric series due to the rapidly increasing cost of evaluating the gamma function as the precision increases. But T(k) is a product of functions that individually can be differentiated logarithmically in closed form (in terms of polygamma functions). I therefore implemented the equivalent of symbolic high-order differentiation for products and the exponential function using generators.

diffs_prod generates the high-order derivatives of a product, given generators for products of individual factors:


>>> f = lambda x: exp(x)*cos(x)*sin(x)
>>> u = diffs(f, 1)
>>> v = diffs_prod([diffs(exp,1), diffs(cos,1), diffs(sin,1)])
>>> next(u); next(v)
1.23586333600241
1.23586333600241
>>> next(u); next(v)
0.104658952245596
0.104658952245596
>>> next(u); next(v)
-5.96999877552086
-5.96999877552086
>>> next(u); next(v)
-12.4632923122697
-12.4632923122697

diffs_exp generates the high-order derivatives of exp(f(x)), given a generator for the derivatives of f(x):


>>> def diffs_loggamma(x):
...     yield loggamma(x)
...     i = 0
...     while 1:
...         yield psi(i,x)
...         i += 1
...
>>> u = diffs_exp(diffs_loggamma(3))
>>> v = diffs(gamma, 3)
>>> next(u); next(v)
2.0
2.0
>>> next(u); next(v)
1.84556867019693
1.84556867019693
>>> next(u); next(v)
2.49292999190269
2.49292999190269
>>> next(u); next(v)
3.44996501352367
3.44996501352367

Actually only differentiation of the exponential function is necessary, which I didn't realize at first, but the product differentiation is still a nice feature to have.

Exponential differentiation amounts to constructing a polynomial of n+1 variables (the derivatives of f), which has P(n) (partition function) ≈ exp(√n) terms. Despite the rapid growth, it is indeed faster to use this approach to differentiate gamma function products at high precision, since the precision can be kept constant. In fact, this should even work in double precision (fp arithmetic) although I think it's broken right now due to some minor bug.

As a result of all this, here are two examples that now work (the results have full accuracy):


>>> mp.dps = 25
>>> hyper(['1/3',1,'3/2',2], ['1/5','11/6','41/8'], 1)
2.219433352235586121250027

>>> hyper(['1/3',1,'3/2',2], ['1/5','11/6','5/4'], 1)
+inf
>>> eps1 = extradps(6)(lambda: 1 - mpf('1e-6'))()
>>> hyper(['1/3',1,'3/2',2], ['1/5','11/6','5/4'], eps1)
2923978034.412973409330956

Application: roots of polynomials

A neat application is to compute the roots of quintic polynomials. Such roots can be expressed in closed form using hypergeometric functions, as outlined in the Wikipedia article Bring radical. This is more interesting in theory than practice, since polynomial roots can be calculated very efficiently using iterative methods.

In particular, the roots of x⁵ - x + t are given (in some order) by

where

Just verifying one case:


>>> mp.dps = 25; mp.pretty = True
>>> t = mpf(3)
>>> for r in polyroots([1,0,0,0,-1,t]):
...     print r
... 
-1.341293531690699250693537
(0.9790618691253385511579325 - 0.6252190807348055061205923j)
(0.9790618691253385511579325 + 0.6252190807348055061205923j)
(-0.3084151032799889258111639 + 1.249926913731033699905407j)
(-0.3084151032799889258111639 - 1.249926913731033699905407j)
>>> -t*hyper(['1/5','2/5','3/5','4/5'], ['1/2','3/4','5/4'], 3125*t**4/256)
(-0.9790618691253385511579325 + 0.6252190807348055061205923j)

Now, we can choose t such that the hypergeometric argument becomes unity:


>>> t = root(mpf(256)/3125,4)
>>> for r in polyroots([1,0,0,0,-1,t]):
...     print r
... 
-1.103842268886161371052433
(0.6687403049764220240032331 + 3.167963942396665205386494e-14j)
(0.6687403049764220240032331 - 3.172324181973539074536523e-14j)
(-0.1168191705333413384770166 + 1.034453571405662514459057j)
(-0.1168191705333413384770166 - 1.034453571405662514459057j)

>>> -t*hyper(['1/5','2/5','3/5','4/5'], ['1/2','3/4','5/4'], 1)
-0.6687403049764220240032331

The last line would previously fail.

Interesting to note is that the value at the singular point z = 1 corresponds to a double root. This also makes polyroots return inaccurate values (note the imaginary parts), a known deficiency that should be fixed...

Better methods

In fact, there may be better methods than the Euler-Maclaurin formula. One such approach is to use the Abel-Plana formula, which I've also recently implemented as sumap(). The Abel-Plana formula gives the exact value for an infinite series (subject to some special growth conditions on the analytically extended summand) as a sum of two integrals.

The Abel-Plana formula is particularly useful when the summand decreases like a power of k; for example when the sum is a pure zeta function:


>>> sumap(lambda k: 1/k**2.5, [1,inf])
1.34148725725091717975677
>>> zeta(2.5)
1.34148725725091717975677

>>> sumap(lambda k: 1/(k+1j)**(2.5+2.5j), [1,inf])
(-3.385361068546473342286084 - 0.7432082105196321803869551j)
>>> zeta(2.5+2.5j, 1+1j)
(-3.385361068546473342286084 - 0.7432082105196321803869551j)

It actually works very well for the generalized hypergeometric series at the z = 1 singularity (which, as mentioned earlier, is a generalized polygamma function, polylogarithm, or zeta function of integer argument), and near the singularity as well. It should even be slightly faster than the Euler-Maclaurin formula since no derivatives are required. Yet another possibility is to integrate a Mellin-Barnes type contour for the hypergeometric function directly.

But at present, I don't want to modify the existing hypergeometric code because it works and would only get more complicated. Rather, I want to improve nsum so it can handle all of this more easily without external hacks. The numerical integration code should also be improved first, because there are still certain parameter combinations where the hypergeometric function evaluation fails due to slow convergence in the numerical integration (this is due to an implementation issue and not an inherent limitation of the integration algorithm).

Post Sage Days 24 report

2010-07-27T18:32:00.006+02:00

Now that I've overcome the immediate effects of SDWS (Sage Days Withdrawal Syndrome), I should write a brief report about Sage Days 24 which took place last week at the Research Institute for Symbolic Computation (RISC), located in a small town called Hagenberg just outside Linz in Austria. Many thanks for the organizers for inviting me!

Since I'm starting as a Ph.D. student at RISC in October, it was nice to get an early view of the site and meet some current students and staff. The location is quite beautiful, my only complaint being the extremely high temperature last week (they say the weather is more normal most of the year!).

SD24 was titled "Symbolic Computation in Differential Algebra and Special Functions", which pretty accurately describes the topics covered. Some points of interest:

William Stein talked about his vision for Sage (short and long term), including goals for special functions support.

Frédéric Chyzak gave a talk about the Dynamic Dictionary of Mathematical Functions. I really like the approach of starting from a minimal definition of a special function (a differential equation with initial conditions) and generating series expansions, Chebyshev approximations (etc.) algorithmically.

Nico Temme talked about numerical evaluation of special functions, which was familiar grounds to me, but with much food for thought.

I met Manuel Kauers, my future supervisor, who gave a neat presentation of some recent research.

There were many other interesting talks and discussions as well, and I gave a talk about special function evaluation in mpmath (based on my SD23 talk). And of course, I met various other Sage developers and got some coding done as well.

I had hoped to get generalized hypergeometric functions into Sage during SD24. There is now a patch, but it still needs some more work.

I also implemented parabolic cylinder functions (PCFs) in mpmath, the last remaining chapter of hypergeometric functions listed in the DLMF (mpmath now covers chapters 1-20, 22, 24-25, partially 26-27, and 33). Code commit here.

Technically, parabolic cylinder functions are just scaled (linear combinations of) confluent hypergeometric functions or Whittaker functions, but they are a bit tricky to compute due to cancellation and branch cuts.

An example from Nico Temme's talk was to evaluate D_ν(z) with ν = -300.14 and z = 300.15 (Maple 13 takes 5 seconds to give an answer that has the wrong sign and is off by 692 orders of magnitude). The new pcfd function in mpmath instantaneously (in about a millisecond) gives the correct result:


>>> pcfd(-300.14, 300.15)
6.83322814925713e-10526
>>> mp.dps = 30
>>> pcfd('-300.14', '300.15')
6.83322814923312844480669009487e-10526

All the curve plots from DLMF 12.3 can be reproduced as follows:


f1 = lambda x: pcfu(0.5,x)
f2 = lambda x: pcfu(2,x)
f3 = lambda x: pcfu(3.5,x)
f4 = lambda x: pcfu(5,x)
f5 = lambda x: pcfu(8,x)
plot([f1,f2,f3,f4,f5], [-3,3], [0,3])


f1 = lambda x: pcfv(0.5,x)
f2 = lambda x: pcfv(2,x)
f3 = lambda x: pcfv(3.5,x)
f4 = lambda x: pcfv(5,x)
f5 = lambda x: pcfv(8,x)
plot([f1,f2,f3,f4,f5], [-3,3], [-3,3])


f1 = lambda x: pcfu(-0.5,x)
f2 = lambda x: pcfu(-2,x)
f3 = lambda x: pcfu(-3.5,x)
f4 = lambda x: pcfu(-5,x)
plot([f1,f2,f3,f4], [-8,8], [-6,6])


f1 = lambda x: pcfv(-0.5,x)
f2 = lambda x: pcfv(-2,x)
f3 = lambda x: pcfv(-3.5,x)
f4 = lambda x: pcfv(-5,x)
plot([f1,f2,f3,f4], [-8,8], [-6,6])


f1 = lambda x: pcfu(-8,x)
f2 = lambda x: pcfv(-8,x)*gamma(0.5-(-8))
f3 = lambda x: hypot(f1(x), f2(x))
plot([f1,f2,f3], [-6,6], [-120,120])


f1 = lambda x: diff(pcfu,(-8,x),(0,1))
f2 = lambda x: diff(pcfv,(-8,x),(0,1))*gamma(0.5-(-8))
f3 = lambda x: hypot(f1(x), f2(x))
plot([f1,f2,f3], [-6,6], [-220,220])

A little more information can be found in the mpmath documentation section for PCFs.

Thanks again to the NSF grant enabling this work.

Sage Days 23, and Bessel function zeros

2010-07-10T23:18:00.010+02:00

Yesterday I came back home from Sage Days 23 in Leiden. I don't really have time to write a detailed travel report right now, but overall it was very, very nice. Those interested can check out William Stein's photos (1, 2, 3, 4) of the event (I can be seen in a few of the photos, wearing a blue t-shirt during days 1-3 and a yellow t-shirt during day 4).

I didn't get quite as much coding done as I probably had time for (in part because I had a talk to prepare, and in part because I spent a bunch of time playing with MPIR and FLINT 2 and gaining some insights but otherwise not accomplishing much productive).

The code I did write includes a complete Cython implementation of real exp (giving up to a 7x speedup) and complex sqrt for mpmath, as well as code for fast evaluation of the gamma function at rational arguments. I will submit the Cython code to Sage soon; possibly after rewriting the remaining power code in Cython as well.

Anyway, prompted by a recent question on sage-support (and because I finally got inspiration to do it while at the airport), I've now implemented Bessel function zeros in mpmath.

Bessel function zeros are very important in applied mathematics because they are needed for Fourier–Bessel series, used to solve partial differential equations in cylindrical coordinates. I'd like to write down a concrete computational example, but it'd be slightly involved and I don't really have time right now. Perhaps in a later blog post.

The code permits computing the m-th nonnegative simple zero of J_ν(z), J'_ν(z) Y_ν(z), or Y'_ν(z) for any positive integer index m and real order ν ≥ 0. One can use it for large m, large ν, or both:


>>> besseljzero(1,100)
314.9434728377671624580656
>>> besseljzero(100,1)
108.836165898409774363098
>>> besseljzero(100,100)
459.5295465754674695983379

A word of caution: I've chosen the convention Abramowitz & Stegun of including z = 0 for J'₀, which is incompatible with SciPy and the Mathematica BesselZeros package. The A&S convention makes much more sense to me because it makes the zeros a continuous function of the order. This case (ν=0, derivative=1) should be the only incompatibility.

It's a bit tricky to compute the zeros of Bessel functions, because (as far as I know) no globally valid approximations are known; not any reasonably simple ones anyway.

Asymptotic expansions are known with respect to either the order or the index. In my implementation, I've chosen to simply use the asymptotic expansion with respect to the index n. If the n is too small, the code simply chooses a larger index m for which the asymptotic expansion does work, and then isolates all the zeros between 0 and the m-th by counting sign changes. By caching this result, all the consecutive zeros for the same function order can then be computed quickly as well.

This approach is simple and quite robust (up to any silly implementation errors on my part), although the code could certainly be optimized more for large orders.

The following graphic shows the derivative of Y₅₀(z) with its zeros marked:

A neat property of Bessel functions is that they can be expanded as infinite products over their zeros:

These products converge slowly, but work with the aid of convergence acceleration:


>>> v, z = mpf(3), mpf(0.75)
>>> (z/2)**v/gamma(v+1)*
...     nprod(lambda k: 1-(z/besseljzero(v,k))**2, [1,inf])
...
0.008484383423274108843927552
>>> besselj(v,z)
0.008484383423274108843927552

>>> (z/2)**(v-1)/2/gamma(v)*
...     nprod(lambda k: 1-(z/besseljzero(v,k,1))**2, [1,inf])
...
0.0331364636065541211306414
>>> besselj(v,z,1)
0.0331364636065541211306414

I will do more work on Bessel functions in mpmath and Sage, so watch for future updates.

Symbolic infinite series

2010-07-02T15:37:00.007+02:00

I've written a little module for symbolically evaluating hypergeometric infinite series. Technically, a series is hypergeometric if the ratio between successive terms is a rational function of the index k; more intuitively, most "nice-looking" series (such as power series of most special functions) involving powers, factorials, gamma functions, binomial coefficients (etc.) are hypergeometric.

The code currently uses SymPy, but should be possible to convert to Sage's symbolics (it mostly requires basic pattern matching, and evaluation of standard special functions). Perhaps "rewrite" is a better word than "convert" because the current code is a messy hack. It's buggy, and only handles special cases in places where a general approach could be used (such as for rational functions).

The basic idea behind the code is to perform the evaluation in two steps. Firstly, the given series is converted into canonical form. I've used the standard generalized hypergeometric series _pF_q. This is not always an optimal representation, but usually does the job.

In the second step, parameter transformations and table lookups are used to reduce the _pF_q series into elementary (or simpler) functions when possible. If this step fails, the _pF_q function is returned, giving a closed-form solution that can be manipulated or evaluated numerically.

The second step is the most difficult, because the transformation algorithms that need to be used for hypergeometric functions are very complicated in general (and there are all kinds of special cases). So far I've only implemented some simple transformations. For the final lookup stage, it would be straightforward to also implement expansion into Bessel functions, incomplete gamma functions, etc.

Below is a sample of some working test cases, with numerical evaluation of both the evaluated result and the original series as verification. Some of the sums don't quite simplify fully, i.e. there is an unevaluated _pF_q function left. The formulas were mostly generated with SymPy's latex function, so they are not always as pretty as could be.

>>> hypsum(1 / fac(k), k)
E
>>> _.evalf()
2.71828182845905
>>> mpmath.nsum(lambda k: 1 / mpmath.fac(k), [0,mpmath.inf])
2.71828182845905

>>> hypsum(k**6 / fac(k), k)
203*E
>>> _.evalf()
551.811211177186
>>> mpmath.nsum(lambda k: k**6 / mpmath.fac(k), [0,mpmath.inf])
551.811211177186

>>> hypsum(4 * d * k * z**(k+2) / fac(k), k)
4*d*z**3*exp(z)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.0401257942714919
>>> mpmath.nsum(lambda k: 4 * 0.5 * k * 0.25**(k+2) / mpmath.fac(k), [0,mpmath.inf])
0.0401257942714919

>>> hypsum(1 / rf(2,k), k)
-1 + E
>>> _.evalf()
1.71828182845905
>>> mpmath.nsum(lambda k: 1 / mpmath.rf(2,k), [0,mpmath.inf])
1.71828182845905

>>> simplify(hypsum(k**3 * z**(k) / fac(k), k))
z*exp(z) + z**3*exp(z) + 3*z**2*exp(z)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.581824016936633
>>> mpmath.nsum(lambda k: k**3 * 0.25**(k) / mpmath.fac(k), [0,mpmath.inf])
0.581824016936633

>>> simplify(hypsum(1 / binomial(2*k,k), k))
4/3 + 2*pi*3**(1/2)/27  >>> _.evalf()
1.73639985871872
>>> mpmath.nsum(lambda k: 1 / mpmath.binomial(2*k,k), [0,mpmath.inf])
1.73639985871872

>>> hypsum(1 / binomial(3*k,k), k)
2*pi*3**(1/2)*3F2([1/2, 1, 1], [1/3, 2/3], 4/27)/(3*gamma(1/3)*gamma(2/3))
>>> _.evalf()
1.41432204432182
>>> mpmath.nsum(lambda k: 1 / mpmath.binomial(3*k,k), [0,mpmath.inf])
1.41432204432182

>>> hypsum(k**4 / binomial(2*k,k), k)
32/3 + 238*pi*3**(1/2)/243
>>> _.evalf()
15.9961018356512
>>> mpmath.nsum(lambda k: k**4 / mpmath.binomial(2*k,k), [0,mpmath.inf])
15.9961018356512

>>> hypsum(fac(k+1)*fac(k+2)/fac(k+3)/fac(k+4), k)
3F2([1, 2, 3], [4, 5], 1)/72
>>> _.evalf()
0.0217325998184402
>>> mpmath.nsum(lambda k: mpmath.fac(k+1)*mpmath.fac(k+2)/mpmath.fac(k+3)/mpmath.fac(k+4), [0,mpmath.inf])
0.0217325998184402

>>> hypsum(1/((3+k)**2*(8+6*k+k**2)), k)
3F2([1, 2, 3], [4, 5], 1)/72
>>> _.evalf()
0.0217325998184402
>>> mpmath.nsum(lambda k: 1/((3+k)**2*(8+6*k+k**2)), [0,mpmath.inf])
0.0217325998184402

>>> hypsum(fac(2*k+1)/(fac(2*k+2)*(k+1)), k)
pi**2/12
>>> _.evalf()
0.822467033424113
>>> mpmath.nsum(lambda k: mpmath.fac(2*k+1)/(mpmath.fac(2*k+2)*(k+1)), [0,mpmath.inf])
0.822467033424113

>>> hypsum(fac(2*k+1)/(fac(3*k+2)), k)
2*pi*3**(1/2)*2F2([1, 3/2], [4/3, 5/3], 4/27)/(27*gamma(4/3)*gamma(5/3))
>>> _.evalf()
0.553106730441975
>>> mpmath.nsum(lambda k: mpmath.fac(2*k+1)/(mpmath.fac(3*k+2)), [0,mpmath.inf])
0.553106730441975

>>> hypsum(4**k*fac(k+R32)*fac(k+R52)/(fac(2*k+6)), k)
3*pi/256
>>> _.evalf()
0.0368155389092554
>>> mpmath.nsum(lambda k: 4**k*mpmath.fac(k+1.5)*mpmath.fac(k+2.5)/(mpmath.fac(2*k+6)), [0,mpmath.inf], method='e')
0.0368155389092358

>>> hypsum((-1)**k * z**(2*k+1) / fac(k) / (2*k+1), k)
pi**(1/2)*erf(z)/2
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.244887887180256
>>> mpmath.nsum(lambda k: (-1)**k * 0.25**(2*k+1) / mpmath.fac(k) / (2*k+1), [0,mpmath.inf])
0.244887887180256

>>> hypsum(z**(2*k+1) / (2*k+1), k)
atanh(z)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.255412811882995
>>> mpmath.nsum(lambda k: 0.25**(2*k+1) / (2*k+1), [0,mpmath.inf])
0.255412811882995

>>> hypsum(z**k / (2*k+1), k)
atanh(z**(1/2))/z**(1/2)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
1.09861228866811
>>> mpmath.nsum(lambda k: 0.25**k / (2*k+1), [0,mpmath.inf])
1.09861228866811

>>> hypsum((-z)**k / (k+1), k)
log(1 + z)/z
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.892574205256839
>>> mpmath.nsum(lambda k: (-0.25)**k / (k+1), [0,mpmath.inf])
0.892574205256839

>>> hypsum(fac(k-R12)/((1+2*k)*fac(k))*z**(2*k), k)
pi**(1/2)*asin(z)/z
>>> _.evalf(subs={d:'1/2',z:'1/4'})
1.79145636509746
>>> mpmath.nsum(lambda k: mpmath.fac(k-0.5)/((1+2*k)*mpmath.fac(k))*0.25**(2*k), [0,mpmath.inf])
1.79145636509746

>>> hypsum(z**k / fac(k+2), k)
-(1 + z - exp(z))/z**2
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.544406667003864
>>> mpmath.nsum(lambda k: 0.25**k / mpmath.fac(k+2), [0,mpmath.inf])
0.544406667003864

>>> hypsum((3*z**2-5)**(-4*k+3), k)
-(5 - 3*z**2)**3/(1 - 1/(5 - 3*z**2)**4)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
-111.666432272586
>>> mpmath.nsum(lambda k: (3*0.25**2-5)**(-4*k+3), [0,mpmath.inf])
-111.666432272586

>>> hypsum(z**(2*k+1) / fac(2*k+1), k)
sinh(z)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.252612316808168
>>> mpmath.nsum(lambda k: 0.25**(2*k+1) / mpmath.fac(2*k+1), [0,mpmath.inf])
0.252612316808168

>>> hypsum((-1)**k * z**(2*k+1) / fac(2*k+1), k)
sin(z)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.247403959254523
>>> mpmath.nsum(lambda k: (-1)**k * 0.25**(2*k+1) / mpmath.fac(2*k+1), [0,mpmath.inf])
0.247403959254523

>>> hypsum((-1)**k * d**k * (1-z)**(2*k+1) / fac(2*k), k)
(1 - z)*cos(d**(1/2)*(1 - z))
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.646980115568008
>>> mpmath.nsum(lambda k: (-1)**k * 0.5**k * (1-0.25)**(2*k+1) / mpmath.fac(2*k), [0,mpmath.inf])
0.646980115568008

>>> hypsum(k * z**(2*k) / fac(2*k+1), k)
cosh(z)/2 - sinh(z)/(2*z)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
0.0104819163234500
>>> mpmath.nsum(lambda k: k * 0.25**(2*k) / mpmath.fac(2*k+1), [0,mpmath.inf])
0.01048191632345

>>> hypsum(gamma(k-R12)**2/(fac(k)**2), k)
16
>>> _.evalf()
16.0000000000000
>>> mpmath.nsum(lambda k: mpmath.gamma(k-0.5)**2/(mpmath.fac(k)**2), [0,mpmath.inf])
16.0

>>> hypsum(1/(k+1)**2, k)
pi**2/6
>>> _.evalf()
1.64493406684823
>>> mpmath.nsum(lambda k: 1/(k+1)**2, [0,mpmath.inf])
1.64493406684823

>>> hypsum(k/(k+1)**3, k)
-zeta(3) + pi**2/6
>>> _.evalf()
0.442877163688632
>>> mpmath.nsum(lambda k: k/(k+1)**3, [0,mpmath.inf])
0.442877163688632

>>> simplify(hypsum((3+k)/(3+3*k)*z**k/fac(k), k))
-(2 - 2*exp(z) - z*exp(z))/(3*z)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
1.18540958339656
>>> mpmath.nsum(lambda k: (3+k)/(3+3*k)*0.25**k/mpmath.fac(k), [0,mpmath.inf])
1.18540958339656

>>> hypsum(k**2 * z**(2*k+1) / fac(2*k+9), k)
z**3*(-1F2([2], [6, 13/2], z**2/4) + 2*1F2([3], [6, 13/2], z**2/4))/39916800
>>> _.evalf(subs={d:'1/2',z:'1/4'})
3.92066920165299e-10
>>> mpmath.nsum(lambda k: k**2 * 0.25**(2*k+1) / mpmath.fac(2*k+9), [0,mpmath.inf])
3.92066920165299e-10

>>> simplify(hypsum((k+1)*(k+2)*(k+3)*(1+2*k)*z**k, k))
(6 + 42*z)/(1 - 5*z + 10*z**2 - 10*z**3 + 5*z**4 - z**5)
>>> _.evalf(subs={d:'1/2',z:'1/4'})
69.5308641975309
>>> mpmath.nsum(lambda k: (k+1)*(k+2)*(k+3)*(1+2*k)*0.25**k, [0,mpmath.inf])
69.5308641975309

>>> hypsum(k**3 * (-z)**k / (k+1), k)
-z*(2*(1 - z)/(1 + z)**3 - ((2 + 2*z)/(1 + z) - (2 + 2*z)*log(1 + z)/z)/(z*(1 + z)) - 2/(1 + z)**2)/2
>>> _.evalf(subs={d:'1/2',z:'1/4'})
-0.0285742052568390
>>> mpmath.nsum(lambda k: k**3 * (-0.25)**k / (k+1), [0,mpmath.inf])
-0.028574205256839

>>> hypsum(E**k / fac(k), k)
exp(E)
>>> _.evalf()
15.1542622414793
>>> mpmath.nsum(lambda k: mpmath.e**k / mpmath.fac(k), [0,mpmath.inf])
15.1542622414793

>>> hypsum(cos(k) / fac(k), k)
exp(exp(I))/2 + exp(exp(-I))/2
>>> _.evalf()
1.14383564379164 + .0e-20*I
>>> mpmath.nsum(lambda k: mpmath.cos(k) / mpmath.fac(k), [0,mpmath.inf])
1.14383564379164

>>> hypsum(cos(k) / (2*k+1), k)
atanh(exp(-I/2))*exp(I/2)/2 + atanh(exp(I/2))*exp(-I/2)/2
>>> _.evalf()
0.975556628913311
>>> mpmath.nsum(lambda k: mpmath.cos(k) / (2*k+1), [0,mpmath.inf])
0.975556628913311

>>> simplify(hypsum(cos(k) * k / fac(k), k))
exp(I + exp(I))/2 + exp(-I + exp(-I))/2
>>> _.evalf()
-0.458967373729452 + .0e-20*I
>>> mpmath.nsum(lambda k: mpmath.cos(k) * k / mpmath.fac(k), [0,mpmath.inf])
-0.458967373729452

>>> hypsum((-1)**k / (2*k+1), k)
pi/4
>>> _.evalf()
0.785398163397448
>>> mpmath.nsum(lambda k: (-1)**k / (2*k+1), [0,mpmath.inf])
0.785398163397448

>>> hypsum((-1)**k / (2*k+6), k)
-1/4 + log(2)/2
>>> _.evalf()
0.0965735902799727
>>> mpmath.nsum(lambda k: (-1)**k / (2*k+6), [0,mpmath.inf])
0.0965735902799727

>>> hypsum((-1)**k / (2*k+1)**2, k)
Catalan
>>> _.evalf()
0.915965594177219
>>> mpmath.nsum(lambda k: (-1)**k / (2*k+1)**2, [0,mpmath.inf])
0.915965594177219

>>> hypsum((-1)**k / (2*k+2)**2, k)
pi**2/48
>>> _.evalf()
0.205616758356028
>>> mpmath.nsum(lambda k: (-1)**k / (2*k+2)**2, [0,mpmath.inf])
0.205616758356028

>>> hypsum((-1)**k * (k+1) / (2*k+3)**2, k)
3F2([3/2, 3/2, 2], [5/2, 5/2], -1)/9
>>> _.evalf()
0.0652837153898853
>>> mpmath.nsum(lambda k: (-1)**k * (k+1) / (2*k+3)**2, [0,mpmath.inf])
0.0652837153898854

>>> hypsum(1/(k**2+2*k+4), k)
I*3**(1/2)*(-polygamma(0, 1 + I*3**(1/2)) + polygamma(0, 1 - I*3**(1/2)))/6
>>> _.evalf()
0.740267076581851
>>> mpmath.nsum(lambda k: 1/(k**2+2*k+4), [0,mpmath.inf])
0.740267076581851

This work was possible thanks to the support of NSF grant DMS-0757627, which is gratefully acknowledged.

Incomplete elliptic integrals complete

2010-06-21T17:53:00.009+02:00

I'm very happy to finally have implemented incomplete elliptic integrals in mpmath (they were probably the most-requested missing special functions). Now mpmath can compute the three classical (Legendre) elliptic integrals F(φ, m), E(φ, m), Π(n, φ, m) as well Carlson's symmetric integrals R_F, R_C, R_J, R_D, R_G. Previously only the complete integrals K(m) and E(m) were available.

The functions are called, respectively, ellipf, ellipe, ellippi, elliprf, elliprc, elliprj, elliprd, elliprg, although this could change before the release. For the code, see r1162, r1166 and adjacent commits; documentation is available in the elliptic functions section.

This work was possible thanks to the support of NSF grant DMS-0757627, gratefully acknowledged.

Elliptic integral basics

What are elliptic integrals and why are they important? As the name suggests, they are related to ellipses. Given an ellipse of width 2a and height 2b, i.e. satisfying

or using an angular coordinate θ

the arc length along the ellipse from θ = 0 to θ = φ is given by

where m = 1-(a/b)^2 is the so-called elliptic parameter.

Here are some plots of half-ellipses of various proportions, and the corresponding arc lengths:


a = 1
b1 = 0.1; f1 = lambda x: b1*sqrt(1 - (x/a)**2)
b2 = 0.6; f2 = lambda x: b2*sqrt(1 - (x/a)**2)
b3 = 1.0; f3 = lambda x: b3*sqrt(1 - (x/a)**2)
b4 = 2.0; f4 = lambda x: b4*sqrt(1 - (x/a)**2)
plot([f1,f2,f3,f4], [-1,1])

g1 = lambda phi: b1*ellipe(phi, 1-(a/b1)**2)
g2 = lambda phi: b2*ellipe(phi, 1-(a/b2)**2)
g3 = lambda phi: b3*ellipe(phi, 1-(a/b3)**2)
g4 = lambda phi: b4*ellipe(phi, 1-(a/b4)**2)
plot([g1,g2,g3,g4], [0,pi])

The preceding formulas are sometimes seen with a and b switched, which simply corresponds to transposing x and y as the reference axis (above, the integration is assumed to start on the x axis). One can readily confirm that E(φ, m) satisfies obvious geometric symmetries:


>>> a = 0.5
>>> b = 0.25
>>> m1 = 1-(a/b)**2
>>> m2 = 1-(b/a)**2
>>> b*ellipe(pi/2,m1)           # Arc length of quarter-ellipse
0.6055280137842297624017815
>>> a*ellipe(pi/2,m2)
0.6055280137842297624017815
>>> a*ellipe(pi/4,m2) + b*ellipe(pi/4,m1)
0.6055280137842297624017815
>>> a*ellipe(pi/3,m2) + b*ellipe(pi/6,m1)
0.6055280137842297624017815

Elliptic integrals have countless uses in geometry and physics, and even in pure mathematics. They are related to Jacobi elliptic functions (already available in mpmath), which essentially are inverse functions of elliptic integrals:


>>> m = 0.5
>>> sn = ellipfun('sn')
>>> ellipf(asin(sn('0.65', m)), m)
0.65

More information about elliptic integrals can be found in DLMF Chapter 19.

Implementation

Implementing incomplete elliptic integrals properly is fairly complicated. When people have asked for them in the past, I just suggested using numerical integration. This works quite well most of the time, but it's inefficient (and possibly gives poor accuracy) at high precision, for large arguments, or near singularities. Another approach — computing incomplete integrals from Appell hypergeometric functions (available in mpmath) — basically has the same drawbacks.

The numerically conscious way to compute elliptic integrals is to firstly exploit symmetries to obtain a small standard domain, and then using transformation formulas to reduce the magnitude of the arguments until the first few terms of the hypergeometric series give full accuracy (the arithmetic-geometric mean for the complete integrals, which leads to some of the fastest ways to compute π, is a special case of this process).

The transformations for Legendre's integrals are known as Landen's transformations. A more modern alternative, that I've chosen to use, is to represent Legendre's integrals in terms of Carlson's symmetric integrals which have much simpler structure. Although Carlson's integrals are mainly intended for internal use, I've exposed them as top-level functions since it wasn't much extra work and some users may well find them useful. (Carlson's forms seem to become increasingly standard.)

Carlson's paper is very readable and provides almost complete descriptions of the algorithms, including error bounds, which helped my implementation work tremendously. The error bounds are not rigorous in all cases, but good enough most of the time. Unfortunately, there are many special cases to consider, and Carlson does not address all of them explicitly, so putting it all together nevertheless involved a bit of work (and there are still some minor issues to fix).

Complex branch structure

One major issue is how to define the elliptic integrals for complex (or in some places, negative) parameters. Because the integrands defining Legendre's elliptic integrals involve square roots of periodic functions, they have very complicated branch structure. I have chosen the periodic extension of the branches resulting naturally from use of the Carlson forms on -π/2 < Re(z) < π/2. I believe this is equivalent to choosing the principal-branch square root in the integrand and integrating along a path where the principal branch is continuous.

For example, consider E(φ, 3+4i). The integrand (with the principal square root, as computed by sqrt) and the elliptic integral computed by ellipe are plotted below:

>>> cplot(lambda z: sqrt(1-(3+4j)*sin(z)**2), points=50000)

>>> cplot(lambda z: ellipe(z, 3+4j), points=50000)

It can be seen that the cuts have the same shape in both images (the branch points are necessarily the same). As the following figure shows, if we are to integrate from z = 0 to z = -2+4i (X), a straight path (black) will be bad as it crosses a branch cut (O). However, the modified path (blue) avoiding the branch cut by passing through z = -2 is fine:

We can confirm this with numerical integration:


>>> ellipe(-2+4j, 3+4j)
(36.45544164923053561755261 + 49.28080768743760310823023j)
>>> quad(lambda z: sqrt(1-(3+4j)*sin(z)**2), [0,-2+4j])
(19.65109997738606076596704 + 46.63216059660102110020723j)
>>> quad(lambda z: sqrt(1-(3+4j)*sin(z)**2), [0,-2,-2+4j])
(36.45544164923053561755261 + 49.28080768743760310823023j)

Unfortunately, Mathematica doesn't use quite the same branch cuts for its elliptic integrals. I don't have Mathematica available to create a plot for comparison purposes, but I think it basically uses vertical cuts instead of the curvilinear cuts seen in the images above. There are probably good reasons for doing this, presumably that it simplifies symbolic definite integration, although it seems to complicate the evaluation and formulaic representation of the functions. It would perhaps be good to support both conventions.

As far as the symmetric functions are concerned, they have the nice property that although they involve non-principal square roots (chosen so as to be continuous along the real line), Carlson's algorithm gives the continuous branch automatically. For example, if we consider the first function

with x = i-1, y = i, z = 0, then as the following plot shows, the principal branch of the integrand has a discontinuity at t = 1/2:


>>> x,y,z = j-1,j,0
>>> f = lambda t: 1/sqrt((t+x)*(t+y)*(t+z))
>>> plot(f, [0,3])

To obtain the correct integrand, we switch to the negative square root on the left of the discontinuity (the positive sign should be used on the right because, for continuity with real parameters, we want the principal square root as t → +∞):


>>> g = lambda t: f(t) if t >= 0.5 else -f(t)
>>> plot(g, [0,3])

A quick test shows that f is wrong and g is right:


>>> extradps(25)(quad)(f, [0,inf])/2
(1.110136183412179775397387 - 0.04278194644898641314159222j)
>>> extradps(25)(quad)(g, [0,inf])/2
(0.7961258658423391329305694 - 1.213856669836495986430094j)
>>> elliprf(x,y,z)
(0.7961258658423391329305694 - 1.213856669836495986430094j)

In summary, branch cuts of special functions are a tricky business, and standards don't always exist, so both developers and users need to be careful working with them.

Assorted special functions update

2010-06-15T22:56:00.008+02:00

Over the week-and-a-half since the last blog update, I've gotten a bunch more work done on special functions in mpmath.

Function plots in the documentation

I have started adding graphics of special functions to the documentation. So far, I've done most of the Bessel-type functions and orthogonal polynomials. Many more to come!

Example screenshot (see the Bessel functions page):

New inhomogeneous Bessel functions

There exists a large number of lesser-known special functions which are essentially variations of Bessel functions. These include functions which solve the generalized (inhomogeneous) Bessel differential equation

with some specific right-hand side g(z). New additions to mpmath in this category are the Anger function (angerj()), Weber function (webere()), and Lommel functions (lommels1(), lommels2()). See commits here and here.

More information about Anger-Weber functions and Lommel functions can be found in the DLMF.

In the near future, I will probably further improve the implementations of the main Bessel functions. The Bessel functions are mostly implemented as generic hypergeometric functions, but the standard cases can be tuned a great deal with special-purpose code.

Airy functions and related functions

The Airy functions Ai and Bi have been present for quite some time in mpmath. In a recent commit, I have rewritten them for improved rigor and better performance at high precision. There are also some new features, such as the ability to evaluate derivatives or iterated integrals of arbitrary order.

Derivatives:


>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> airyai(1.5, derivative=5)
0.211387453153454489799743
>>> diff(airyai, 1.5, 5)
0.211387453153454489799743
>>> airyai(1.5, derivative=100)
-6.480220187791312407132043e+49
>>> airybi(0, derivative=1000)
3.754976097101270163249629e+854
>>> airybi(0, derivative=1001)
0.0
>>> airybi(0, derivative=1002)
3.756228172694934424506624e+856

Integrals:


>>> airyai(5, derivative=-1)
0.3332875903059178794866562
>>> quad(airyai, [0,5])
0.3332875903059178794866562
>>> airyai(-100000, derivative=-1)
-0.6665753658794626398413214

Also, functions for computing the zeros of Ai and Bi (and the first derivatives) have been added:


>>> airyaizero(1)
-2.338107410459767038489197
>>> airyaizero(2)
-4.087949444130970616636989
>>> airybizero(1)
-1.17371322270912792491998
>>> airybizero(1, derivative=1)
-2.294439682614123246622459
>>> airybizero(1, derivative=1, complex=True)
(0.2149470745374305676088329 + 1.100600143302797880647194j)
>>> airybizero(10000)
-1304.584974702601410702964
>>> airybizero(10000, complex=True)
(652.3059222438076432024695 + 1129.846189716375208308414j)

I have also implemented two new functions related to Airy functions: the Scorer functions Gi and Hi. These are available as scorergi() and scorerhi() respectively.

Here are two plots of the Gi-function, which can also be seen in the documentation:

Interval gamma functions

The interval arithmetic context now implements gamma, rgamma (reciprocal gamma function), factorial as well as loggamma for real as well as complex arguments (commit). For example:


>>> iv.dps = 10
>>> iv.gamma('50.3')
[1.96282982095908e+63, 1.96282982457481e+63]
>>> iv.gamma(iv.mpc('2.7','5.9'))
([0.00269836072064322, 0.00269836072271801] + 
 [0.0120124287790304, 0.0120124287810768]*j)

As a "practical" example, consider evaluating the Riemann-Siegel theta function which involves computing the difference of two log-gamma functions. For input with a large real part, the imaginary part in the result suffers from massive cancellation and may end up with the wrong sign:


>>> mp.dps = 15
>>> mp.siegeltheta(10**50 + 0.25j)
(5.61456887916465e+51 - 0.143091235731175j)
>>> mp.dps = 10; nprint(mp.siegeltheta(10**50 + 0.25j).imag)
-0.143091
>>> mp.dps = 100; nprint(mp.siegeltheta(10**50 + 0.25j).imag)
14.1614

With interval arithmetic, the sign uncertainty is reflected in the output:


>>> iv.dps = 15
>>> iv.siegeltheta(10**50 + 0.25j)
([5.6145688791646467648e+51, 5.6145688791646474294e+51] +
[-5.0706024009129187319e+30, 5.070602400912917606e+30]*j)
>>> iv.dps = 50
>>> iv.siegeltheta(10**50 + 0.25j)
([5614568879164647368060513633451316140100495086670736.0, 5614568879164647368060
513633451316140100495086670744.0] + 
[14.1613521236438249782320715810808676610440
8814838558203, 14.16147419395632497823207158108086766104408814838560341]*j)

As another example, consider evaluating the gamma function of a huge argument. The digits in the answer may be "wrong" because the input is converted from decimal to binary, and the gamma function is sensitive to the input being perturbed:


>>> mp.dps = 15
>>> mp.gamma('123456789012345.1')
6.11544992055093e+1686076589184486
>>> mp.dps = 30
>>> mp.gamma('123456789012345.1')
7.49032018540342193592769680745e+1686076589184486
>>> mp.dps = 60
>>> mp.gamma('123456789012345.1')
7.49032018540342058679709881225047421518964527875047787194339e+1686076589184486

With interval arithmetic, the uncertainty in the input is propagated correctly:


>>> iv.dps = 15
>>> iv.nprint(iv.gamma('123456789012345.1'), mode='diff')
[6.11545e+1686076589184486, 1.01533e+1686076589184487]
>>> iv.dps = 30
>>> iv.nprint(iv.gamma('123456789012345.1'), 20, mode='diff')
7.4903201854034[185631, 219359]e+1686076589184486
>>> iv.dps = 60
>>> iv.nprint(iv.gamma('123456789012345.1'), 50, mode='diff')
7.49032018540342058679709881225047421518964527[60898, 87505]e+1686076589184486

Rewritten Lambert W function

Lastly, the Lambert W function has received a much-needed rewrite (commit) mainly to improve evaluation very close to the branch cut along the negative axis and particularly near the branch point at -1/e for the k = -1, 0, 1 branches.

With the previous implementation, results were frequently inaccurate or ended up on the wrong branch in this region. Here are some hard cases that now work perfectly:


>>> mp.dps = 1000
>>> x = -1/e + mpf('1e-900')
>>> y = -1/e - mpf('1e-900')
>>> z = -1/e + mpf('1e-900')*1j
>>> w = -1/e - mpf('1e-900')*1j
>>> mp.dps = 25
>>> lambertw(x,0); lambertw(y,0); lambertw(z,0); lambertw(w,0)
-1.0
(-1.0 + 2.331643981597124203363536e-450j)
(-1.0 + 1.648721270700128146848651e-450j)
(-1.0 - 1.648721270700128146848651e-450j)
>>> lambertw(x,1); lambertw(y,1); lambertw(z,1); lambertw(w,1)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-1.0 + 1.648721270700128146848651e-450j)
>>> lambertw(x,-1); lambertw(y,-1); lambertw(z,-1); lambertw(w,-1)
-1.0
(-1.0 - 2.331643981597124203363536e-450j)
(-1.0 - 1.648721270700128146848651e-450j)
(-3.088843015613043855957087 - 7.461489285654254556906117j)

To finish this post, I present the following Mathematica bug:


remote1:frejohl:[~]$ math
Mathematica 7.0 for Linux x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.

In[1]:= Im[LambertW[0,-1/E-10^(-900)]]

Out[1]= 0

Announcing mpmath 0.15

2010-06-06T15:18:00.006+02:00

I'm happy to announce the release of mpmath 0.15!

This should've happened earlier, but I was obstructed by other obligations (in particular, finishing my master's thesis). The good news is that I will be working full-time on special functions in mpmath and Sage again this summer thanks to sponsorship provided by William Stein (with money from an NSF grant). I will also come to Sage Days 23 (July 5-9, Leiden, the Netherlands) and Sage Days 24 (July 17-22, Linz, Austria) which should be as fun as SD15.

What's new in mpmath 0.15? As usual, the details can be found in the CHANGES file and the list of commits. Most major changes were covered in detail in previous blog posts, so I will first of all simply link to those posts along with short summaries:

Speedups of elementary functions - cos, sin, atan, cosh, sinh, exp and all derived functions are now faster. Of particular note, trigonometric functions use an asymptotically faster algorithm (all elementary functions have similar asymptotic performance now).

A new gamma function implementation - the gamma function and log-gamma functions have been rewritten scratch, using faster algorithms and code optimizations. The new versions are uniformly faster than the old ones, and tens or hundreds of times faster in many important situations. They are also more accurate in some special cases (near singularities, for complex arguments with extremely small real part, ...).

Numerical multidimensional infinite series - nsum() can now evaluate series (finite or infinite) in any number of dimensions

Computing large zeta zeros with mpmath - Juan Arias de Reyna has implemented code for computing the nth zero of the Riemann zeta function on the critical line for arbitrarily large n. On a related note, the Riemann zeta function code has been optimized (mostly through elimination of low-level overhead) to speed up such computations.

Hypergeneralization - generalized 2D hypergeometric series, bilateral hypergeometric series and some q-analogs (q-Pochhammer symbol, q-factorial, q-gamma function, q-hypergeometric series) have been implemented.

In addition, there are many changes that I haven't had time to blog about yet. I will now write briefly about some of these:

Elliptic functions

The support for working with elliptic functions has been improved. Instead of just the standard three, all 12 Jacobi elliptic functions are now available via ellipfun, e.g. the cd-function:


>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> ellipfun('cd', 3.5, 0.5)
-0.9891101840595543931308394
>>> cd = ellipfun('cd')
>>> cd(3.5, 0.5)
-0.9891101840595543931308394

Functions qfrom(), qbarfrom(), mfrom(), kfrom(), taufrom() have also been added to convert between the various forms of arguments to elliptic functions (nomes (two definitions), parameters, moduli, half-period ratios). If you're like me and never can remember such formulas, let alone avoid messing up when trying to apply them, this can be quite convenient. Some functions (in particular, ellipfun also support direct choice of convention using keyword arguments):


>>> ellipfun('sn',2,1+1j)    # default is m
(1.333680690847080418060154 - 0.2084318767795699518406447j)
>>> ellipfun('sn',2,q=qfrom(m=1+1j))
(1.333680690847080418060154 - 0.2084318767795699518406447j)
>>> ellipfun('sn',2,k=kfrom(m=1+1j))
(1.333680690847080418060154 - 0.2084318767795699518406447j)
>>> ellipfun('sn',2,tau=taufrom(m=1+1j))
(1.333680690847080418060154 - 0.2084318767795699518406447j)

Another new function is the Klein j-invariant (or rather, the "absolute invariant" or "J-function" which differs by a constant factor, and also is the normalization used by Mathematica). A nice application is to compute the Laurent series expansion in terms of the (number-theoretic) nome, giving a famous integer sequence:


>>> mp.dps = 15
>>> taylor(lambda q: 1728*q*kleinj(qbar=q), 0, 5, singular=True)
[1.0, 744.0, 196884.0, 21493760.0, 864299970.0, 20245856256.0]

The J-function also looks pretty when plotted, here as a function of the number-theoretic nome within the unit circle:


fp.cplot(lambda q: fp.kleinj(qbar=q), [-1,1], [-1,1], points=400000)

As a function of the half-period ratio τ, defined in the upper half-plane:

fp.cplot(lambda t: fp.kleinj(tau=t), [-1,2], [0,1.5], points=100000)

Complex interval arithmetic

The support for interval arithmetic has been extended; there is also a new interface for intervals. It is no longer possible to call mpmath.exp etc. directly with interval arguments; instead all interval arithmetic has been moved to a separate context object called iv (similar to the separate interface for working with Python float/complex in mpmath). For example:


>>> iv.dps = 5; iv.pretty = True
>>> iv.mpf(1)/3
[0.3333330154, 0.3333334923]
>>> iv.sin(1)
[0.8414707184, 0.8414716721]
>>> iv.sin([1,2])
[0.8414707184, 1.0]
>>> iv.sin([1,3])
[0.141119957, 1.0]

The interface change was done to simplify the implementation; separating out interval arithmetic also helps avoid inadvertent mixing of interval and non-interval arithmetic.

There is also some support for complex interval arithmetic, including some elementary functions (exp, log, cos, sin, power).


>>> iv.mpc('1.3','1.4') ** 2
([-0.2700066566, -0.2699956894] + [3.639995575, 3.640010834]*j)
>>> iv.mpc('1.3','1.4') ** iv.mpc('1.2','-0.7')
([3.329280853, 3.329311371] + [1.967472076, 1.967498779]*j)
>>> iv.dps = 25
>>> iv.mpc('1.3','1.4') ** 2
([-0.270000000000000000000000061005, -0.269999999999999999999999912371] +
[3.63999999999999999999999988419, 3.64000000000000000000000009099]*j)
>>> mp.cos(2+3j)
(-4.189625690968807230132555 - 9.109227893755336597979197j)
>>> iv.cos(2+3j)
([-4.18962569096880723013255508975, -4.18962569096880723013255498635] + 
[-9.1092278937553365979791973006, -9.10922789375533659797919709381]*j)

Verifying Euler's identity:


>>> iv.dps = 5
>>> iv.exp(1j * iv.pi)
([-1.0, -0.9999990463] + [-1.279517164e-6, 2.535183739e-6]*j)
>>> iv.dps = 25
>>> iv.exp(1j * iv.pi)
([-1.0, -0.999999999999999999999999987075] +
[-1.8806274411915714063022730148e-26, 3.28925138726485156164403137746e-26]*j)

Rigorous complex linear algebra:


>>> iv.dps = 15
>>> A = iv.matrix([[2,-2,1],[1,0,2],[1,1,2]])
>>> b = iv.matrix([1,2+2j,5])
>>> iv.lu_solve(A,b)
[ ([4.0, 4.0] + [-3.3333333333333333339, -3.333333333333333333]*j)]
[([3.0, 3.0] + [-2.0000000000000000004, -1.9999999999999999998]*j)]
[ ([-1.0, -1.0] + [2.6666666666666666665, 2.666666666666666667]*j)]
>>> A*iv.lu_solve(A,b)
[([1.0, 1.0] + [-8.8817841970012523234e-16, 1.7763568394002504647e-15]*j)]
[         ([2.0, 2.0] + [1.9999999999999995559, 2.0000000000000008882]*j)]
[([5.0, 5.0] + [-8.8817841970012523234e-16, 1.7763568394002504647e-15]*j)]

Partial derivatives

This is rather simple, but very convenient. diff can now compute partial derivatives of arbitrary order for multivariable functions. For example, a first derivative with respect to the second argument:


>>> diff(lambda x,y: 3*x*y + 2*y - x, (0.25, 0.5), (0,1))
2.75

The first derivative with respect to both x and y:


>>> diff(lambda x,y: 3*x*y + 2*y - x, (0.25, 0.5), (1,1))
3.0

Second derivative with respect to a parameter of a hypergeometric function ₂F₁(a,b,c,z); tenth derivative with respect to the argument; first derivative with respect to two parameters and fourth derivative with respect to the argument:


>>> diff(hyp2f1, (1,2,3,0.5), (0,0,2,0))
0.2306514802159766624657329
>>> diff(hyp2f1, (1,2,3,0.5), (0,0,0,10))
672869392.1810426015397572
>>> diff(hyp2f1, (1,2,3,0.5), (1,0,1,4))
-585.1335939248098118753126

That's it for now. Stay tuned for many new features in the near future!

Hypergeneralization

2010-03-17T10:21:00.011+01:00

The generalized hypergeometric function can itself be generalized endlessly. I have recently implemented three such extensions in mpmath: bilateral series, two-dimensional hypergeometric series, and q-analog (or "basic") hypergeometric series.

Bilateral series

The bilateral hypergeometric series is the simplest extension, and consists of taking the usual hypergeometric series and extending the range of summation from [0,∞) to (-∞,∞):

This series only converges when |z| = 1 and A = B, but interpreted as a sum of two ordinary hypergeometric series, it can be assigned a value for arbitrary z through the analytic continuation (or Borel regularization) of the ordinary hypergeometric series, which mpmath implements. Anyway, the convergent case is the interesting one. Here one obtains, for instance, Dougall's identity for 2H2:


>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> a,b,c,d = 0.5, 1.5, 2.25, 3.25
>>> bihyper([a,b],[c,d],1)
-14.49118026212345786148847
>>> gammaprod([c,d,1-a,1-b,c+d-a-b-1],[c-a,d-a,c-b,d-b])
-14.49118026212345786148847

As an example of regularization, the divergent 1H0 series can be expressed as the sum of one 2F0 function and one 1F1 function:


>>> a = mpf(0.25)
>>> z = mpf(0.75)
>>> bihyper([a], [], z)
(0.2454393389657273841385582 + 0.2454393389657273841385582j)
>>> hyper([a,1],[],z) + (hyper([1],[1-a],-1/z)-1)
(0.2454393389657273841385582 + 0.2454393389657273841385582j)
>>> hyper([a,1],[],z) + hyper([1],[2-a],-1/z)/z/(a-1)
(0.2454393389657273841385582 + 0.2454393389657273841385582j)

Two-dimensional series

The most common hypergeometric series of two variables, i.e. a twodimensional series whose summand is a hypergeometric expression with respect to both indices separately, is the Appell F1 function, previously available in mpmath as appellf1:

However, much more general functions are possible. There are three other Appell functions: F2, F3, F4. The Horn functions are 34 distinct functions of order two, containing the Appell functions as special cases. The Kampé de Fériet function provides a generalization of Appell functions to arbitrary orders.

The new hyper2d function in mpmath can evaluate all these named functions, and more general functions still. The trick for speed is to write the series as a nested series, where the inner series is a generalized hypergeometric series that can be evaluated efficiently with hyper, and where the outer series has a rational recurrence formula. This rewriting also permits evaluating the analytic continuation with respect to the inner variable (as implemented by hyper).

The user specifies the format of the series in quasi-symbolic form, and the rewriting to nested form is done automatically by mpmath. For example, the Appell F1 function can be computed as


hyper2d({'m+n':[a], 'm':[b1], 'n':[b2]}, {'m+n':[c]}, x, y)

and indeed, this is essentially what appellf1 now does internally. The Appell F2-F4 functions have also been added explicitly as appellf2, appellf3, appellf4.

Hypergeometric functions of two (or more) variables have numerous applications, such as solving high-order algebraic equations, expressing various derivatives and integrals in closed form, and solving differential equations, but I have not yet found any simple examples that make good demonstrations except for F1. (I have mostly found examples of that take half a page to write down.) Any such examples for the documentation would be a welcome contribution!

Some trivial examples from the documentation are:


>>> x, y = mpf(0.25), mpf(0.5)
>>> hyper2d({'m':1,'n':1}, {}, x,y)
2.666666666666666666666667
>>> 1/(1-x)/(1-y)
2.666666666666666666666667
>>> hyper2d({'m':[1,2],'n':[3,4]}, {'m':[5],'n':[6]}, x,y)
4.164358531238938319669856
>>> hyp2f1(1,2,5,x)*hyp2f1(3,4,6,y)
4.164358531238938319669856
>>> hyper2d({'m':1,'n':1},{'m+n':1},x,y)
2.013417124712514809623881
>>> (exp(x)*x-exp(y)*y)/(x-y)
2.013417124712514809623881

An example of a Horn function, H3:


>>> x, y = 0.0625, 0.125
>>> a,b,c = 0.5,0.75,0.625
>>> hyper2d({'2m+n':a,'n':b},{'m+n':c},x,y)
1.190003093972956004227425
>>> nsum(lambda m,n: rf(a,2*m+n)*rf(b,n)/rf(c,m+n)*\
...     x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1.190003093972956004227425

This also demonstrates the recently added generic support for multidimensional infinite series in mpmath. But of course, nsum is much slower than hyper2d.

Hypergeometric q-series

Before introducing the q-analog of the hypergeometric series, I should introduce the q-Pochhammer symbol,

This itself is a new function in mpmath, implemented as qp(a,q,n) (with two- and one-argument forms qp(a,q) and qp(q) also permitted) and is the basis for more general computation involving q-analogs. The q-factorial and q-gamma function have also been added (as qfac and qgamma), but are not yet documented.

The q-analogs have important applications in number theory. As a very neat example, numerically computing the Taylor series of 1/(q, q)_∞ with mpmath gives


>>> taylor(lambda q: 1/qp(q), 0, 10)
[1.0, 1.0, 2.0, 3.0, 5.0, 7.0, 11.0, 15.0, 22.0, 30.0, 42.0]

These are the values of the partition function P(n) for n = 0,1,2,..., i.e. the number of ways of writing n as a sum of positive integers.

Replacing the rising factorials (Pochhammer symbols) in the generalized hypergeometric series with their q-analogs gives the hypergeometric q-series or basic hypergeometric series

This function is implemented as qhyper. Like hyper, it supports arbitrary combinations of real and complex arguments (assuming |q| < 1). Some examples from the documentation:


>>> qhyper([0.5], [2.25], 0.25, 4)
-0.1975849091263356009534385
>>> qhyper([0.5], [2.25], 0.25-0.25j, 4)
(2.806330244925716649839237 + 3.568997623337943121769938j)
>>> qhyper([1+j], [2,3+0.5j], 0.25, 3+4j)
(9.112885171773400017270226 - 1.272756997166375050700388j)

Like hyper, it automatically ensures accurate evaluation for alternating series:


>>> q = 0.998046875
>>> mp.dps=5; qhyper([2],[0.5], q, -0.5)
6.6738e-69
>>> mp.dps=15; qhyper([2],[0.5], q, -0.5)
6.67376764851253e-69
>>> mp.dps=25; qhyper([2],[0.5], q, -0.5)
6.673767648512527695718826e-69
>>> mp.dps=100; qhyper([2],[0.5], q, -0.5)
6.673767648512527695718826106778799352769798151218768443717704076836963752188876
561888441662933081804e-69

With the q-analog of the generalized hypergeometric function implemented, it becomes possible to compute q-exponentials, q-sines, q-orthogonal polynomials, q-Bessel functions, and pretty much anything else. If there is interest, such function could be added explicitly to mpmath.

For more information and examples of the functions discussed in this post, see the sections on hypergeometric functions and q-functions (a little terse at the moment) in the mpmath documentation.

Computing large zeta zeros with mpmath

2010-03-12T14:04:00.006+01:00

Juan Arias de Reyna, who wrote the code used in mpmath for evaluating the Riemann zeta function near the critical strip, has contributed a new implementation of the zetazero function which computes the nth zero on the critical line.

The old version (written by myself) used a lookup table for initial approximations, so it was limited to computing the first few zeros. Juan's code calculates the position of arbitrary zeros using Gram's law, with a sophisticated algorithm (about 300 lines of code) to find the correct zero when Gram's law (which is actually just a heuristic) fails.

A bunch of zeros, from small to large:


from mpmath import *
from timeit import default_timer as clock
mp.dps = 20
for n in range(14):
    t1 = clock()
    v1 = zetazero(10**n)
    t2 = clock()
    v2 = zetazero(10**n+1)
    t3 = clock()
    print "10<sup>%i</sup>  "%n, v1, "%.2f" % (t2-t1)
    print "10<sup>%i</sup>+1"%n, v2, "%.2f" % (t3-t2)


n      value                          time (s)
10⁰    (0.5 + 14.13472514173469379j)  0.12
10⁰+1  (0.5 + 21.022039638771554993j) 0.05
10¹    (0.5 + 49.773832477672302182j) 0.05
10¹+1  (0.5 + 52.970321477714460644j) 0.04
10²    (0.5 + 236.5242296658162058j)  0.32
10²+1  (0.5 + 237.769820480925204j)   0.30
10³    (0.5 + 1419.4224809459956865j) 0.93
10³+1  (0.5 + 1420.416526323751136j)  0.84
10⁴    (0.5 + 9877.7826540055011428j) 3.75
10⁴+1  (0.5 + 9878.6547723856922882j) 3.71
10⁵    (0.5 + 74920.827498994186794j) 2.53
10⁵+1  (0.5 + 74921.929793958414308j) 2.58
10⁶    (0.5 + 600269.67701244495552j) 2.43
10⁶+1  (0.5 + 600270.30109071169866j) 2.89
10⁷    (0.5 + 4992381.014003178666j)  1.63
10⁷+1  (0.5 + 4992381.2627112065366j) 2.30
10⁸    (0.5 + 42653549.760951553903j) 2.36
10⁸+1  (0.5 + 42653550.046758478876j) 2.93
10⁹    (0.5 + 371870203.83702805273j) 6.98
10⁹+1  (0.5 + 371870204.36631304458j) 6.72
10¹⁰   (0.5 + 3293531632.3971367042j) 11.48
10¹⁰+1 (0.5 + 3293531632.6869557853j) 15.92
10¹¹   (0.5 + 29538618431.613072811j) 53.67
10¹¹+1 (0.5 + 29538618432.07777426j)  43.00
10¹²   (0.5 + 267653395648.62594824j) 162.01
10¹²+1 (0.5 + 267653395648.84752313j) 174.67
10¹³   (0.5 + 2445999556030.2468814j) 1262.55
10¹³+1 (0.5 + 2445999556030.6222451j) 1456.38

The function correctly separates close zeros, and can optionally return information about the separation. For example:


>>> mp.dps = 20; mp.pretty = True
>>> zetazero(542964976,info=True)
((0.5 + 209039046.578535272j), [542964969, 542964978], 6, '(013111110)')

The extra information is explained in the docstring for zetazero:

[This means that the] zero is between Gram points 542964969 and 542964978, it is the 6-th zero between them. Finally (01311110) is the pattern of zeros in this interval. The numbers indicate the number of zeros in each Gram interval, (Rosser Blocks between parenthesis). In this case only one Rosser Block of length nine.

Juan reports having computed the 10¹⁵th zero, which is


208514052006405.46942460229754774510611

and verifying that it satisfies the Riemann hypothesis. And fortunately, the values of large zeros computed with mpmath seem to agree with those reported by Andrew Odlyzko and Xavier Gourdon.

The new zetazero function in mpmath is indeed able to separate the closest known pair of zeros, found in Xavier Gourdon's exhaustive search up to 10¹³:


>>> mp.dps = 25
>>> v1 = zetazero(8637740722917, info=True)
>>> v2 = zetazero(8637740722918, info=True)
>>> v1
((0.5 + 2124447368584.392964661515j), [8637740722909L, 8637740722925L], 7L,
'(1)(1)(1)(1)(1)(1)(02)(1)(02)(1)(1)(20)(1)')
>>> v2
((0.5 + 2124447368584.392981706039j), [8637740722909L, 8637740722925L], 8L,
'(1)(1)(1)(1)(1)(1)(02)(1)(02)(1)(1)(20)(1)')
>>> nprint(abs(v1[0]-v2[0]))
1.70445e-5

This computation takes over 2 hours on my laptop.

Of course, zeros can be computed to any desired precision:


>>> mp.dps = 1000
>>> print zetazero(1)
(0.5 + 14.1347251417346937904572519835624702707842571156992431756855674601499634
29809256764949010393171561012779202971548797436766142691469882254582505363239447
13778041338123720597054962195586586020055556672583601077370020541098266150754278
05174425913062544819786510723049387256297383215774203952157256748093321400349904
68034346267314420920377385487141378317356396995365428113079680531491688529067820
82298049264338666734623320078758761792005604868054356801444424651065597568665903
22868651054485944432062407272703209427452221304874872092412385141835146054279015
24478338354254533440044879368067616973008190007313938549837362150130451672696838
92003917628512321285422052396913342583227533516406016976352756375896953767492033
61272092599917304270756830879511844534891800863008264831251691127106829105237596
17977431815170713545316775495153828937849036474709727019948485532209253574357909
22612524773659551801697523346121397731600535412592674745572587780147260983080897
860071253208750939599796666067537838121489190886j)

The only real limitation of the code, as far as I can tell, is that it's impractical for computing a large number of consecutive zeros. For that one needs to use something like the Odlyzko–Schönhage algorithm for multi-evaluation of the zeta function.

A nice exercise would be to parallelize the zeta code using multiprocessing. The speedup for large zeros should be essentially linear with the number of processors.

Speedups of elementary functions

2010-03-11T23:26:00.005+01:00

I have a fairly large backlog of changes in mpmath to blog about. For today, I'll just briefly mention the most recent one, which is a set of optimizations to elementary functions. The commit is here. I intend to provide much faster Cython versions of the elementary functions some time in the future, but since there was room for optimizing the Python versions, I decided to invest a little time in that.

Most importantly, the asymptotic performance of trigonometric function has been improved greatly, due to a more optimized complexity reduction strategy. The faster strategy was previously used only for exp, but cos/sin are virtually identical in theory so fixing this was mostly a coding problem. In fact, the high-precision code for exp/cosh/sinh/cos/sin is largely shared now, which is nice because only one implementation needs to be optimized.

The following image shows performance for computing cos(3.7). The red graph is the old implementation; blue is new. The top subplot shows evaluations/second at low precision, i.e. higher is better, and the bottom subplot shows seconds/evaluation at higher precision, i.e. lower is better.

Already at a few thousand bits, the new code is more than twice as fast, and then it just gets better and better. There is also a little less overhead at low precision now, so the new code is uniformly faster.

Also exp, cosh and sinh have been optimized slightly. Most importantly, there is a little less overhead at low precision, but asymptotic speed has also improved a bit.

Performance for computing cosh(3.7):

Performance for computing exp(3.7):

I should note that the performance depends on some tuning parameters which naturally are system-specific. The results above are on a 32-bit system with gmpy as the backend. I hope to later implement optimized tuning parameters for other systems as well.

Elementary function performance is of course important globally, but it's especially important for some specific functions. For example, the speed of the Riemann zeta function depends almost proportionally on the total speed of evaluating one exp, one log, one sine, and one cosine, since the main part of the work consists of summing the truncated L-series

With the improved elementary functions, and some overhead removals to the zetasum code in the same commit, the Riemann zeta function is up to 2x faster now, and about 50% faster on the critical line.

Numerical multidimensional infinite series

2010-02-12T17:50:00.006+01:00

Another new feature in mpmath: support for multidimensional series in nsum. Some very simple examples (finite and infinite summations can be combined in any desired order):


>>> nsum(lambda i,j: 2**(-i-j), [0,inf], [0,inf])
4.0
>>> nsum(lambda i,j,k: 2**(-i-j-k), [0,inf], [0,inf], [0,inf])
8.0
>>> nsum(lambda i,j,k,l: i/2**j*(i+l)/2**k, [1,2], [0,inf], [1,inf], [2,3])
50.0
>>> nsum((lambda i,j,k,l: 1/(2**(i**2+j**2+k**2+l**2))), *([[-inf,inf]]*4))
20.5423960756379
>>> jtheta(3,0,0.5)**4
20.5423960756379

One could of course also make nested calls to nsum, but having a direct syntax is much more convenient. Nested evaluation is also usually inefficient for convergence acceleration, so this is not what nsum does internally. Instead, it combines all the infinite summations to a single summation over growing hypercubes. The distinction is very important for conditionally convergent series.

For example, nsum can now directly evaluate the Madelung constant in three dimensions (ignore=True is to ignore the singular term at the origin):


>>> nsum(lambda i,j,k: (-1)**(i+j+k)/(i**2+j**2+k**2)**0.5,
...     [-inf,inf], [-inf,inf], [-inf,inf], ignore=True)
-1.74756459463318

While this evaluation takes several seconds, it is somewhat remarkable that a precise value can be obtained at all. A better way to compute the Madelung constant is using the following rapidly convergent 2D series:


>>> f = lambda i,j: -12*pi*sech(0.5*pi*sqrt((2*i+1)**2+(2*j+1)**2))**2
>>> nsum(f, [0,inf], [0,inf])
-1.74756459463318
>>> mp.dps = 100
>>> nsum(f, [0,inf], [0,inf])
-1.74756459463318219063621203554439740348516143662474175815282535076504062353276
117989075836269460789

Another nice application is to evaluate Eisenstein series directly from the definition:


>>> tau = 1j
>>> q = qfrom(tau=tau)
>>> nsum(lambda m,n: (m+n*tau)**(-4), [-inf,inf], [-inf,inf], ignore=True)
(3.1512120021539 + 0.0j)
>>> 0.5*sum(jtheta(n,0,q)**8 for n in [2,3,4])*(2*zeta(4))
(3.1512120021539 + 0.0j)

A new gamma function implementation

2010-02-11T23:09:00.014+01:00

The gamma function is probably the most important nonelementary special function. The gamma function or ratios of gamma functions appear in everything from the functional equation of the Riemann zeta function to the normalization factors of hypergeometric functions. Even expressions that are usually rational functions (such as binomial coefficients) are often most conveniently expressed (or in software, implemented) in terms of the gamma function since this automatically provides correct asymptotics and extension to noninteger parameters. Needless to say, a well-implemented gamma function is fundamental for a special functions library.

The gamma function implementation used in mpmath until now is one of the oldest parts of the code (I wrote it around three years ago), and although it has done its job, it has some flaws. It mainly uses Spouge's approximation, which is fairly efficient but nonoptimal for asymptotically large arguments. The implementation also has some subtle precision issues.

The new gamma function added in this commit is the result of over a year of intermittent work. I did most of the work during last summer but didn't have time to finish it until recently. The improvements are essentially the following:

Removal of various overheads

Special-casing half-integers

Using Taylor series for small real arguments

Using Stirling's series for complex and/or large arguments

Optimizing for computing loggamma

Handling values near poles and near the real axis extremely carefully

The difference in performance is quite dramatic. The smallest speedup occurs for small complex arguments, but it turns out that even there, the Stirling series is faster than Spouge's approximation (at least with the removal of overheads in the new implementation) despite the need to perform argument transformations. For real arguments, and large complex arguments, the new implementation is consistently faster by a large factor, especially for loggamma.

I'll show some benchmarks. The times are in milliseconds; the first value in each cell is the time for the old implementation and the second is the time for the new one, followed by the speedup. The results are on my laptop, with gmpy (I'll try with Sage later):

Firstly, results for gamma(z):

z	dps=15	dps=50	dps=250	dps=1000
3.0	0.0086 0.0077 (1.13x)	0.0089 0.0076 (1.16x)	0.0087 0.0075 (1.16x)	0.0085 0.0077 (1.12x)
5.27	0.1200 0.0342 (3.51x)	0.2513 0.0545 (4.61x)	2.5964 0.4220 (6.15x)	80.7089 11.3011 (7.14x)
2.5	0.1255 0.0158 (7.94x)	0.2398 0.0138 (17.41x)	2.5519 0.0144 (177.72x)	79.0651 0.0149 (5319.87x)
-3.7	0.2518 0.0457 (5.51x)	0.4006 0.0615 (6.52x)	3.1073 0.4283 (7.26x)	87.4200 11.3981 (7.67x)
(-3.0+4.0j)	0.6205 0.3942 (1.57x)	1.0095 0.6093 (1.66x)	9.0896 5.9099 (1.54x)	253.9219 191.3358 (1.33x)
(3.25-4.125j)	0.3548 0.2864 (1.24x)	0.7456 0.4791 (1.56x)	8.6651 5.6367 (1.54x)	252.7598 189.2187 (1.34x)
(15.0+20.0j)	0.3980 0.1973 (2.02x)	0.7526 0.4496 (1.67x)	8.6583 5.4604 (1.59x)	250.8646 188.5849 (1.33x)
52.1	0.1661 0.0586 (2.83x)	0.3050 0.0894 (3.41x)	2.9862 0.5685 (5.25x)	84.8334 12.4169 (6.83x)
123.25	0.1643 0.0638 (2.57x)	0.3071 0.0963 (3.19x)	2.8908 0.9194 (3.14x)	83.8168 14.6588 (5.72x)
-200.7	0.3011 0.1485 (2.03x)	0.4897 0.1914 (2.56x)	3.5145 1.1359 (3.09x)	89.8617 17.5500 (5.12x)
(100.25+100.0j)	0.4238 0.1857 (2.28x)	0.7641 0.2757 (2.77x)	8.6807 3.1403 (2.76x)	249.6731 171.4369 (1.46x)
12345678.9	0.3032 0.0701 (4.32x)	0.4382 0.0857 (5.11x)	3.4078 0.3257 (10.46x)	88.4082 5.2828 (16.74x)
(0.0+12345678.9j)	0.7814 0.1409 (5.54x)	1.1671 0.1802 (6.48x)	10.3985 0.6169 (16.86x)	265.6192 8.3569 (31.78x)
1.0e+20	0.8833 0.0798 (11.06x)	1.3019 0.0899 (14.48x)	6.4010 0.2939 (21.78x)	107.0924 3.5141 (30.47x)
(0.0+1.0e+20j)	3.5025 0.1515 (23.11x)	3.6163 0.1834 (19.72x)	19.5676 0.5308 (36.86x)	322.9033 6.2148 (51.96x)

Secondly, results for loggamma(z):

z	dps=15	dps=50	dps=250	dps=1000
3.0	0.0700 0.0115 (6.07x)	0.0677 0.0116 (5.84x)	0.0683 0.0119 (5.73x)	0.0683 0.0121 (5.65x)
5.27	0.2207 0.0564 (3.91x)	0.3548 0.0796 (4.45x)	2.8925 0.5221 (5.54x)	88.3556 12.6452 (6.99x)
2.5	0.2210 0.0319 (6.93x)	0.3565 0.0374 (9.53x)	2.8595 0.1011 (28.28x)	81.8333 1.3546 (60.41x)
-3.7	0.4392 0.0782 (5.62x)	0.6418 0.0987 (6.50x)	3.5767 0.5429 (6.59x)	88.8484 12.7363 (6.98x)
(-3.0+4.0j)	1.1235 0.4942 (2.27x)	1.5685 0.7224 (2.17x)	10.0879 6.0459 (1.67x)	267.2893 198.7613 (1.34x)
(3.25-4.125j)	0.8789 0.2704 (3.25x)	1.2326 0.4752 (2.59x)	9.7501 5.4572 (1.79x)	264.7649 190.4007 (1.39x)
(15.0+20.0j)	0.8983 0.1248 (7.20x)	1.2923 0.4221 (3.06x)	9.7917 5.2766 (1.86x)	265.2711 187.9775 (1.41x)
52.1	0.2604 0.0812 (3.21x)	0.4186 0.1140 (3.67x)	3.2493 0.6890 (4.72x)	87.0111 13.9280 (6.25x)
123.25	0.2629 0.0478 (5.50x)	0.4227 0.0692 (6.11x)	3.1110 1.0216 (3.05x)	85.9028 15.9917 (5.37x)
-200.7	0.5439 0.1631 (3.34x)	0.7465 0.2221 (3.36x)	3.9305 1.2537 (3.14x)	92.5380 18.8309 (4.91x)
(100.25+100.0j)	0.9103 0.1179 (7.72x)	1.3332 0.1849 (7.21x)	10.1655 3.1083 (3.27x)	266.3565 170.3099 (1.56x)
12345678.9	0.4026 0.0491 (8.21x)	0.5439 0.0525 (10.36x)	3.7183 0.2051 (18.13x)	89.9402 4.0044 (22.46x)
(0.0+12345678.9j)	1.2628 0.0765 (16.52x)	1.7085 0.0858 (19.91x)	11.5023 0.2749 (41.84x)	275.6174 4.4548 (61.87x)
1.0e+20	0.9820 0.0495 (19.84x)	1.3882 0.0540 (25.72x)	6.6983 0.1481 (45.21x)	110.1757 2.1954 (50.19x)
(0.0+1.0e+20j)	4.0385 0.0788 (51.24x)	4.3758 0.0833 (52.54x)	20.9478 0.1936 (108.20x)	335.9407 2.3911 (140.50x)

All times are from a warm cache. The new implementation has slightly longer precomputation time (for the Taylor series coefficients) at high precision, and the Taylor series therefore isn't used above a certain precision (around 1000 digits) despite being much faster on subsequent calls. I will probably add some user-visible function for controlling this.

Finally, several optimizations are possible still. The most important would be to do the implementation in Cython. (Improving the underlying elementary functions will also speed up the gamma function.) Potential algorithmic improvements include faster handling of near-integers (using Taylor series with a reduced number of terms), avoiding the reflection formula for complex arguments in the left half-plane but not too close to the negative half-axis, and performing the argument transformations with a reduced number of multiplications. But I doubt I will have time (or interest -- the present code was demanding enough to finish) to work of any of these things, at least any time soon.

The speedups naturally also affect performance of many other functions, e.g. hypergeometric functions; I might post some benchmarks of that later.

Mpmath 0.14 released

2010-02-05T18:08:00.013+01:00

I've just released version 0.14 of mpmath! Some highlights in this release have been covered in previous development updates on this blog:

mpmath in Sage to become 3x faster -- a Cython extension module will be added to Sage that greatly speeds up mpmath (the Sage patch is currently being reviewed; if all goes to plan, the patch will be accepted soon and the mpmath version in Sage will be updated to 0.14 at the same time).

Zeta evaluation with the Riemann-Siegel expansion -- Juan Arias de Reyna contributed code for very efficient and robust evaluation of the Riemann zeta function for arguments with large imaginary part.

Various features discussed in this update: 3D surface plotting (wrapping matplotlib, matrix calculus (transcendental functions of a matrix argument), Borel summation of divergent hypergeometric series, and options for whether to use Mathematica's conventions for hypergeometric functions have been added.

Improved accuracy for hypergeometric functions with large parameters, and analytic continuation implemented for 3F2, 4F3, and higher hypergeometric functions.

Support for using Python floats/complexes for faster low-precision math. This is very handy for plotting, multidimensional integration, and other things requiring many function evaluations.

There are many other changes as well. See the CHANGES file and the list of commits for more details. Thanks to Juan Arias de Reyna, Vinzent Steinberg, Jorn Baayen and Chris Smith who contributed patches, and various people who tested and submitted bug reports (thanks also to anyone else I forgot to mention).

I'm sorry that it took several months to finish this release! This was partly due to large internal reorganizations done in order to support floats and the Cython backend in Sage. I've also had to take some time off to focus on my studies and other things.

An example: spherical harmonics

Here I will demonstrate three new features with a single example: 3D plotting, one of the newly added special functions (spherical harmonics), and low-precision evaluation with the fp context. Of course, everything here can be done in arbitrary precision with mp as well; fp is just faster for plotting.

The following code visualizes spherical harmonics using a 3D polar plot. I've chosen to visualize the real part; the code is easily edited to show the absolute value (for example) instead.


from mpmath import fp

def Y(m,n):
    def g(theta,phi):
        R = fp.re(fp.spherharm(m,n,theta,phi))**2
        x = R*fp.cos(phi)*fp.sin(theta)
        y = R*fp.sin(phi)*fp.sin(theta)
        z = R*fp.cos(theta)
        return [x,y,z]
    return g

fp.splot(Y(0,0), [0,fp.pi], [0,2*fp.pi])

The first few spherical harmonics are thus:

Y(0,0)
Y(1,0)	Y(1,1)
Y(2,0)	Y(2,1)	Y(2,2)
Y(3,0)	Y(3,1)	Y(3,2)
		Y(3,3)

Some numerical examples

A selection of evaluations that were not implemented, failed, gave inaccurate results or were extremely slow in previous versions of mpmath:


>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True

# reflection formula for Barnes G-function
>>> barnesg(-100000+10000j)
(-7.332534002219787256775675e+22857579769 +
  1.304469872717249495403882e+22857579770j)

# Riemann zeta function, large height
>>> zeta(0.5+100000000j)
(-3.362839487530727943146807 + 1.407234559646447885979583j)

# Accurate evaluation of large-degree Bernoulli polynomial
>>> bernpoly(1000,100)
4.360903799140670486890619e+1996

# Computation of Euler numbers
>>> eulernum(20)
370371188237525.0
>>> eulernum(50, exact=True)
-6053285248188621896314383785111649088103498225146815121L
>>> eulernum(2000000000000000000)
3.19651108713502662532039e+35341231273461153426

# Accurate evaluation of hypergeometric functions with large parameters
>>> hyp2f1(1000,50,25,-1)
-2.886992344761576738138864e-274
>>> legenp(0.5, 300, 0.25)
-1.717508549497252387888262e+578
>>> besseli(2, 10, derivative=100)
821.2386210064833486609255

# Evaluation of 4F3 and 4F2
>>> hyper([1,2.125,3.25,4.5], [5.25,6.5,7], 1+2j)
(1.006737556189825231039221 + 0.3058770612264986390003873j)
>>> hyper([1,2.125,3.25,4.5], [5.25,6.5], 1+2j)
(0.3220085070641791195103201 + 0.5251231752161637627314328j)

# Matrix functions
>>> A = matrix([[2,3,1+j],[1,0,-1],[2,1,5]])
>>> mnorm(A - logm(expm(A)))
9.540212722670391941827867e-25
>>> mnorm(A - expm(logm(A)))
4.375286550134565812513542e-26
>>> nprint(expm(2*A))
[ (-2855.4 + 10733.0j)   (-1871.95 + 8597.9j)  (-7721.42 + 12906.7j)]
[(-3910.68 - 1749.82j)  (-3152.69 - 1272.11j)  (-4460.18 - 3558.54j)]
[ (17039.0 + 22621.8j)   (14320.3 + 17405.9j)   (13587.7 + 35087.8j)]
>>> nprint(expm(A)**2)
[ (-2855.4 + 10733.0j)   (-1871.95 + 8597.9j)  (-7721.42 + 12906.7j)]
[(-3910.68 - 1749.82j)  (-3152.69 - 1272.11j)  (-4460.18 - 3558.54j)]
[ (17039.0 + 22621.8j)   (14320.3 + 17405.9j)   (13587.7 + 35087.8j)]
>>> nprint(chop(sqrtm(A)*sqrtm(A)))
[2.0  3.0  (1.0 + 1.0j)]
[1.0  0.0          -1.0]
[2.0  1.0           5.0]

# A convenience function for precise evaluation of exp(j*pi*x)
>>> expjpi(36); exp(j*pi*36)
(1.0 + 0.0j)
(1.0 + 4.702307256907087875509661e-25j)

# Some fp functions have specialized implementations
# for improved speed and accuracy
>>> fp.erfc(5); float(mp.erfc(5))
1.5374597944280347e-12
1.5374597944280349e-12
>>> fp.erf(5); float(mp.erf(5))
0.99999999999846256
0.99999999999846256
>>> fp.ei(-12-3j); complex(mp.ei(-12-3j))
(4.6085637890261843e-07-3.141592693919709j)
(4.6085637890261901e-07-3.141592693919709j)
>>> fp.zeta(-0.5); float(fp.zeta(-0.5))  # real-argument fp.zeta
-0.20788622497735468
-0.20788622497735468
>>> fp.ci(40); float(mp.ci(40))
0.019020007896208769
0.019020007896208765
>>> fp.gamma(4+5j); complex(mp.gamma(4+5j))
(0.14965532796078551+0.31460331072967695j)
(0.14965532796078521+0.31460331072967596j)

For more examples and images, see the posts linked at the top of this post!

mpmath in Sage to become 3x faster

2010-02-01T01:00:00.013+01:00

I've blogged previously about the potential speedups attainable by rewriting core parts of mpmath in Cython (here and here), but all I had back then was some standalone classes and functions that couldn't easily be integrated. Well, I now finally have a fully working, fully integrated Cython implementation of the mpf and mpc types (plus related utility functions) -- and it successfully runs the entire mpmath test suite!

On my computer (Ubuntu 64 bit, AMD Athlon 64 3700+), running import mpmath; mpmath.runtests() in Sage takes this much time:

Before [mpmath (svn trunk) + Sage 4.3.1]: 62.75 seconds

After [mpmath (svn trunk with modifications) + Cython mpmath types + Sage 4.3.1 + a Sage performance patch by Robert Bradshaw]: 19.87 seconds

I have essentially only implemented arithmetic operations and conversions, so core transcendental functions (as well as square roots and powers) still fall back to the Python code, and they account for the majority of the running time. According to cProfile, about 20% of the total time is spent in exp() and gamma() alone. Hypergeometric series evaluation is also still done in Python. Replacing these functions with Cython versions should give a significant boost, and will be very easy to do in the future (in terms of code infrastructure). So realistically, the running time for the unit tests can be cut much further, probably in half.

The code isn't public yet. I will soon commit the changes to the mpmath svn repository and create a patch for Sage with the Cython extensions. The code just needs a little more cleanup, and there are no doubt a couple of subtle issues (such as corner-case memory leaks) that need fixing... but essentially, it's working, and the tests pass, so I expect it to be releasable quite soon.

Update: after adding just a little more Cythonized wrapper code, the time is now down to 17.64 seconds. I have still not touched exp, log, cos, sin, gamma, or any other core transcendental functions. Expect further improvements.

Update 2 (2010-02-03): preliminary version of the Cython code here

Using Sage numbers in mpmath

2010-01-27T19:14:00.005+01:00

I've written a very basic mpmath context for computing "natively" with Sage's real and complex numbers. It can be tried out by upgrading mpmath in Sage to the svn trunk, applying this patch this patch to fix the helper code in Sage, and then importing this file.

Some examples:


----------------------------------------------------------------------
| Sage Version 4.3.1, Release Date: 2010-01-20                       |
| Type notebook() for the GUI, and license() for information.        |
----------------------------------------------------------------------
sage: from mpsage import SageContext
sage: from mpmath import mp, timing
sage: mp.prec = 150
sage: sc = SageContext(150)
sage: 
sage: print sc.quad(lambda t: sc.exp(-t**2), [0,sc.inf])
0.88622692545275801364908374167057259139877471
sage: print mp.quad(lambda t: mp.exp(-t**2), [0,mp.inf])
0.88622692545275801364908374167057259139877473
sage: timing(sc.quad, lambda t: sc.exp(-t**2), [0,sc.inf])
0.40903496742248535
sage: timing(mp.quad, lambda t: mp.exp(-t**2), [0,mp.inf])
0.43610286712646484
sage: 
sage: print sc.zeta(4+5*I)
0.95121830700949569861514024576822624312072806 + 0.024932824507357722259855155244740571939365897*I
sage: print mp.zeta(4+5*I, method='euler-maclaurin')
(0.95121830700949569861514024576822624312072806 + 0.024932824507357722259855155244740571939365897j)
sage: timing(sc.zeta, 4+5*I)
0.015111517906188966
sage: timing(mp.zeta, 4+5*I, method='euler-maclaurin')
0.028736209869384764
sage: 
sage: print sc.zeta(1+100000000*I)
1.8349308953278466065175598876062914152382527 + 1.0691847904236128969174476736978194200591565*I
sage: print mp.zeta(1+100000000*I)
(1.8349308953278466065175598876062914152563156 + 1.069184790423612896917447673697819420203942j)
sage: timing(sc.zeta, 1+100000000*I)
1.0952680110931396
sage: timing(mp.zeta, 1+100000000*I)
2.0537140369415283

Unfortunately, there are some issues (besides the fact that some methods are missing, so not everything works yet).

This context doesn't provide variable precision, so the user has to manually allocate sufficient extra precision to compensate for rounding errors.

I first tried to do it, but variable precision is very inconvenient to implement using Sage's way of managing precision. There is no direct way to perform operations with a given target precision (independent of the inputs), and the second best option (which is to perform coercions to the target precision everywhere) is very slow, besides losing accuracy. The only way to implement variable precision in a reasonable way is to bypass Sage's RealNumber and ComplexNumber (at least their public interfaces) and wrap MPFR directly using a precision model similar to what MPFR and mpmath uses, where the precision of the result is always specified and independent of the inputs.

Secondly, there is not that much of a speedup (in the examples above, the speedup is at most about 2x). This is mainly due to the fact that the context uses wrapper classes for RealNumber and ComplexNumber, with all interface and conversion code written in Python. So the overhead is about the same as for the corresponding code in vanilla mpmath (where it accounts for about 50% of the total overhead). The reason RealNumber and ComplexNumber can't be used directly, even in a fixed-precision setting, is that mpmath in many places multiplies numbers by floats (usually exact float literals like 0.5), and Sage always coerces to the number with less precision. This could presumably be fixed by replacing all float and complex constants in mpmath, but I'm not in the mood to do that right now.

There should be a significant performance benefit if direct-from-Cython classes and conversion methods were used. It's also very important to optimize certain core functions; for example, quadrature would benefit greatly from a fast fdot implementation.

All things considered, I'm probably not going to continue much more down this particular path. It's better to write a fully compatible Cython context with classes designed directly for mpmath. Trying to wrap Sage's numbers did however help identify a few problems with the existing interfaces, so I might extend the work on this context a little more just to find more such issues.

Zeta evaluation with the Riemann-Siegel expansion

2010-01-19T14:35:00.007+01:00

I'm very grateful to Juan Arias de Reyna who has contributed a module implementing zeta function evaluation using the Riemann-Siegel formula in mpmath. This finally permits rapid evaluation of the Riemann zeta function for arguments with large imaginary part.

To follow tradition on this blog, pictorial proof shall be given. Here is a plot of a segment of the critical line, ζ(1/2+ti) with t between 10¹⁰+50 and 10¹⁰+55:

A complex plot of showing the critical strip, t ranging between 10⁸+40 and 10⁸+45 (note: the y scale is reversed):

Juan Arias de Reyna, who is a professor of mathematics at the University of Seville, has done a thorough job with this code. He has even proved rigorous error bounds for his algorithm (subject to assumptions that the underlying functions in mpmath being well-implemented). The news is that the bounds -- documented in an as-yet unpublished paper -- are valid off the critical line. The code also computes derivatives (up to 4th derivatives), although not as rigorously but still very robustly.

I integrated the code (and added a few optimizations) during the last few days. The zeta function in mpmath now automatically switches between Borwein's algorithm (close to the real line), Euler-Maclaurin summation, and the Riemann-Siegel expansion.

Some example values and timings on my laptop, computing ζ(1/2+10ⁿi) for n up to 12:


>>> from timeit import default_timer as clock
>>> mp.dps = 25
>>>
>>> for n in range(13):
...     t1 = clock(); y = zeta(0.5 + 10**n*j); t2 = clock()
...     print n, y, round(t2-t1,3)
...
0 (0.1439364270771890603243897 - 0.7220997435316730891261751j) 0.02
1 (1.544895220296752766921496 - 0.1153364652712733754365914j) 0.003
2 (2.692619885681324090476096 - 0.02038602960259816177072685j) 0.042
3 (0.3563343671943960550744025 + 0.9319978312329936651150604j) 0.073
4 (-0.3393738026388344575674711 - 0.03709150597320603147434421j) 0.434
5 (1.073032014857753132114076 + 5.780848544363503984261041j) 0.167
6 (0.07608906973822710000556456 + 2.805102101019298955393837j) 0.146
7 (11.45804061057709254500227 - 8.643437226836021723818215j) 0.181
8 (-3.362839487530727943146807 + 1.407234559646447885979583j) 0.336
9 (-2.761748029838060942376813 - 1.677512240989459839213205j) 0.877
10 (0.3568002308560733825395879 + 0.286505849095836103292093j) 2.695
11 (0.6436639255801185727194357 + 0.1168615914616853418448829j) 8.583
12 (2.877961809278403355251079 - 3.206771071318398923493412j) 26.934

YAMDU (yet another mpmath development update)

2010-01-13T11:32:00.009+01:00

At the moment, I'm able to get quite a bit of work done on mpmath. This includes general code cleanup, gardening the issue tracker, and finishing features that were works-in-progress for a long while.

Today's topics are:

3D plotting
Matrix calculus
Borel summation of divergent hypergeometric series
Optional consistency with Mathematica
Internal reorganization

3D plotting

Jorn Baayen submitted a patch that adds 3D surface plotting (along with some other minor changes to the visualization module). Thanks a lot! I previously blogged about 3D plotting using matplotlib.

You can now easily produce plots of functions of the form z = f(x,y):


splot(lambda x, y: 10*sinc(hypot(x,y)), [-10,10], [-10,10])

You can just as easily plot more general parametric functions x,y,z = f(u,v). The following plots a Möbius strip:


def f(u,v):
    d = (1+0.5*v*cos(0.5*u))
    x = d*cos(u)
    y = d*sin(u)
    z = 0.5*v*sin(0.5*u)
    return x,y,z

splot(f, [0,2*pi], [-1,1])

(I'm not sure why there are black spots in the image. This seems to be a matplotlib bug.)

Matrix calculus

It's now possible to compute exponentials, sines, cosines, square roots, logarithms, and complex powers (A^z) of square matrices:


>>> mp.dps = 25; mp.pretty = True
>>> A = matrix([[2,3,1+j],[1,0,-1],[2,1,5]])
>>> nprint(expm(A))
[ (27.6034 + 41.2728j)   (24.9924 + 31.1189j)     (8.67256 + 68.271j)]
[(-15.0962 + 1.51713j)  (-10.1713 + 1.30035j)  (-28.7003 - 0.401789j)]
[ (125.368 + 37.7417j)   (98.9812 + 26.9693j)    (158.616 + 85.1593j)]
>>> nprint(logm(expm(A)))
[(2.0 - 5.35398e-26j)           (3.0 + 6.9172e-26j)           (1.0 + 1.0j)]
[(1.0 + 4.94422e-26j)  (4.42102e-25 - 1.14411e-25j)  (-1.0 - 1.95948e-27j)]
[(2.0 + 1.58966e-26j)          (1.0 - 8.57028e-27j)   (5.0 + 3.23653e-27j)]
>>> nprint(sqrtm(A)**2)
[(2.0 - 8.25839e-30j)          (3.0 - 8.40399e-30j)           (1.0 + 1.0j)]
[(1.0 + 4.63764e-30j)  (4.84564e-30 - 3.27721e-31j)  (-1.0 + 7.47261e-31j)]
[(2.0 - 4.31871e-30j)          (1.0 + 1.78726e-31j)   (5.0 + 2.54582e-29j)]
>>> nprint(powm(powm(A,1+j),1/(1+j)))
[(2.0 - 1.12995e-26j)          (3.0 - 8.93158e-27j)           (1.0 + 1.0j)]
[(1.0 - 1.54352e-26j)  (9.23906e-27 - 1.67262e-26j)  (-1.0 - 2.62243e-27j)]
[(2.0 + 9.97431e-27j)          (1.0 + 1.56341e-26j)   (5.0 - 1.79194e-26j)]

The code also works with the double precision fp context, which obviously is much faster for large matrices. When computing square roots and logarithms, most of the time is spent on matrix inversion, which can be accelerated by substituting in numpy:


>>> from mpmath import fp
>>> A = fp.randmatrix(20)
>>> B = fp.sqrtm(A)
>>> timing(fp.sqrtm, A)
3.6470590535948006
>>>
>>> import numpy
>>> fp.inverse = lambda A: fp.matrix(numpy.linalg.inv(A.tolist()))
>>>
>>> timing(fp.sqrtm, A)
1.0815329881311735
>>> C = fp.sqrtm(A)
>>>
>>> fp.mnorm(A-B**2)
3.5001514923930131e-14
>>> fp.mnorm(A-C**2)
2.6462503380426079e-14

It could be much faster still by doing everything in numpy. Probably in the future fp should be improved to seamlessly use numpy internally for all available matrix operations. Patches are welcome!

Borel summation of divergent hypergeometric series

The hypergeometric series of degree p > q+1 are divergent for |z| > 0. Nevertheless, they define asymptotic expansions for analytic functions which exist in the sense of Borel summation. Previously, mpmath knew only how to evaluate ₂F₀ (whose Borel sum has a closed form), but now it can evaluate the functions of higher degree as well.

To illustrate, the cosine integral Ci(z) has an asymptotic series representation in terms of ₃F₀:

This representation is efficient for very large values of |z|, where the argument of the hypergeometric series is small. But with z = 0.5, say, the series does not yield a single digit. With a value of z = 10 or so, the series only gives about 3 digits with an optimal truncation:


>>> z = 10.0
>>> for n in range(1,16):
...     print sum(rf(0.5,k)*rf(1,k)*rf(1,k)*(-4/(z**2))**k/fac(k) for k in range(n))
...
1.0
0.98
0.9824
0.98168
0.9820832
0.98172032
0.9821993216
0.981327538688
0.9834198176768
0.977017443971072
1.00134646405284
0.888946391275078
1.50939479300832
-2.52351981825774
27.9653146429137

Of course, there are ways to compute the cosine integral using convergent series. But if we pretend that we only know about the asymptotic series, then the Borel regularization means that we can evaluate the function to high precision even for small z:


>>> mp.dps = 25; mp.pretty = True
>>> z = mpf(0.5)
>>> u = -4/z**2
>>> H1 = hyper([1,1,1.5],[],u)
>>> H2 = hyper([0.5,1,1],[],u)
>>>
>>> print log(z)-log(z**2)/2-cos(z)/z**2*H1+sin(z)/z*H2
-0.1777840788066129013358103
>>> print ci(z)
-0.1777840788066129013358103

The z = 10 series above, to high precision:


>>> hyper([0.5,1,1],[],-4/(10.0**2))
0.9819103501017016869905255
>>> mp.dps = 40
>>> hyper([0.5,1,1],[],-4/(10.0**2))
0.9819103501017016869905255431829554224704

It's instructive to visualize the optimal truncations of an asymptotic series compared to the exact solution. Notice how, for an alternating series like this, the truncations alternate between over- and undershooting:


def optimal_asymp(z):
    u = -4/(z**2)
    term_prev = inf
    s = 0.0
    for k in range(30):
        term = rf(0.5,k)*rf(1,k)*rf(1,k)*u**k/fac(k)
        if abs(term) < abs(term_prev):
            s += term
            term_prev = term
        else:
            break
    return s

def exact(z):
    u = -4/(z**2)
    return hyper([0.5,1,1],[],u)

mp.dps = 10
plot([optimal_asymp, exact], [0,10])

The catch with this feature? Computing the Borel regularization involves evaluating (nested) numerical integrals with a hypergeometric function in the integrand. Except for special parameters (where degree reductions happen internally -- the Ci series above happens to be such a case), it's not very fast.

Optional consistency with Mathematica

I often try to follow Mathematica's conventions regarding limits, special values and branch cuts of special functions. This simplifies testing (I can directly compare values with Mathematica), and usually Mathematica's conventions seem well-reasoned.

There are exceptions, however. One such exception concerns Mathematica's interpretation of ₂F₁(a,b,c,z) for b = c a negative integer. Mathematica interprets this as a polynomial, which doesn't make much sense since the denominator of the zero term is zero. It's not even consistent with how Mathematica evaluates this function for a symbolic variable:


In[3]:= Hypergeometric2F1[3,-2,-2,2]

Out[3]= 31

In[4]:= Hypergeometric2F1[3,b,b,2]

Out[4]= -1

Maybe there is a good reason for doing what Mathematica does, but it's at the very least not documented anywhere. I've now changed mpmath to instead interpret the two parameters as eliminating each other and giving a ₁F₀ function (which is what Mathematica does for a generic value). The Mathematica-compatible value can be recovered by explicitly disabling parameter elimination:


>>> mp.pretty = True
>>>
>>> hyp2f1(3,-2,-2,2)
-1.0
>>> hyp2f1(3,-2,-2,2,eliminate=False)
31.0

On a related note, I've fixed the Meijer G-function to switch between its z and 1/z forms automatically to follow Mathematica's definition of the function. This introduces discontinuities in the function for certain orders. The user can explicitly choose which form to use (so as to obtain a continuous function) with series=1 or series=2.

Internal reorganization

In a long overdue change, I've moved many of the modules in mpmath into subdirectories. The multiprecision library code is now located in mpmath/libmp; special functions are in mpmath/functions, calculus routines are in mpmath/calculus, and routines for matrices and linear algebra are in mpmath/matrices.

This shouldn't affect any documented interfaces, but it will break external code that directly uses mpmath internal functions. The mpmath interfaces in SymPy and Sage will need some editing for the next version update. The breakage should be straightforward to fix (mostly just a matter of changing imports).

Since the next version of mpmath is definitely going to break some code, I might use the opportunity to do some other cosmetic interface changes as well. Ideally, after the next version, the interface will be stable until and beyond mpmath 1.0 (whenever that happens). But the version number is still at 0.x for a reason -- no compatibility guarantees.

Accurate hypergeometric functions for large parameters

2010-01-07T13:43:00.005+01:00

Holiday season means more free time for coding (yay!). In this commit, I've fixed the remaining major accuracy issues with hypergeometric functions in mpmath. The effect, generally speaking, is that mpmath will now deal much more robustly with large parameters of hypergeometric-type functions.

Previously, you would obtain an accurate value with parameters of magnitude ≈ 10⁻¹ - 10¹ (say), and often for large parameters as well, but for some functions the accuracy would quickly deteriorate if parameters were increased (or moved closer to singularities). You could still obtain any accuracy by simply increasing the working precision, but you had to manually figure out the amount of extra precision required; the news is that mpmath now automatically gives a value with full accuracy even for large parameters (or screams out loud if it fails to compute an accurate value).

This doesn't mean that you can't trick mpmath into computing a wrong value by choosing sufficiently evil parameters, but it's much harder now than simply plugging in large values. Abnormally large values will now mainly accomplish abnormal slowness (while mpmath implements asymptotic expansions with respect to the argument of hypergeometric functions, evaluation for asymptotically large parameters is a much harder problem as far as I'm aware -- so mpmath in effect accomplishes it by brute force.)

The most trivial change was to change the series summation to aim for control of the relative error instead of the absolute error. This affects alternating series, where large parameters lead to very small function values. Previously, something like the following would happen:


    >>> mp.dps=5; hyp1f1(-5000, 1500, 100)
    1.6353e-14
    >>> mp.dps=15; hyp1f1(-5000, 1500, 100)
    8.09813050863682e-25
    >>> mp.dps=30; hyp1f1(-5000, 1500, 100)
    -1.38318247777802583583082760724e-39

This isn't disastrously bad since the absolute error is controlled (summing this series naively in floating-point arithmetic might give something like 1e+234 due to catastrophic cancellation). But now, you get the relatively accurate answer right away, which is much nicer:


    >>> mp.dps = 5; hyp1f1(-5000, 1500, 100)
    9.1501e-185
    >>> mp.dps = 15; hyp1f1(-5000, 1500, 100)
    9.15012134245639e-185
    >>> mp.dps = 30; hyp1f1(-5000, 1500, 100)
    9.15012134245639443114220499541e-185

Of course, if the value is too small, the evaluation loop must abort eventually. The default behavior now is to raise an exception if relative accuracy cannot be obtained. The user can force a return value by either permitting a higher internal precision before aborting, or specifying a size 2^−zeroprec below which the value is small enough to be considered zero:


    >>> hyp1f1(-8000, 4000, 1500)
    Traceback (most recent call last):
      ...
    ValueError: hypsum() failed to converge to the requested 53 bits of accuracy
    using a working precision of 3568 bits. Try with a higher maxprec,
    maxterms, or set zeroprec.
    >>>
    >>>
    >>> hyp1f1(-8000, 4000, 1500, maxprec=10000)
    4.36754212717293e-1350
    >>> hyp1f1(-8000, 4000, 1500, accurate_small=False)
    -1.99286380611911e-25
    >>> hyp1f1(-8000, 4000, 1500, zeroprec=1000)
    0.0

Since exceptions are inconvenient, it will be necessary to add some more symbolic tests for exact zeros for certain functions (particularly orthogonal polynomials) where exact zeros arise as important special cases.

Also, cancellation detection in hypercomb is now the default for all functions. This fixes, among other things, a bug in the Meijer G-function reported by a user via email. Before:


    >>> meijerg([[],[]],[[0],[]],0.1)
    0.90483741803596
    >>> meijerg([[],[]],[[0,0],[]],0.1)
    1.47347419562219
    >>> meijerg([[],[]],[[0,0,0],[]],0.1)
    1.61746025085449
    >>> meijerg([[],[]],[[0,0,0,0],[]],0.1)   # suspicious
    13194139533312.0

After:


    >>> meijerg([[],[]],[[0],[]],0.1)
    0.90483741803596
    >>> meijerg([[],[]],[[0,0],[]],0.1)
    1.47347419562219
    >>> meijerg([[],[]],[[0,0,0],[]],0.1)
    1.61745972140487
    >>> meijerg([[],[]],[[0,0,0,0],[]],0.1)
    1.56808223438324

Another important change is more correct handling parameters very close to negative integers, particularly those appearing in denominators of series. Previously, unless the integer n was small enough or the precision was set high enough, this situation would yield a bogus value (that of a prematurely truncated series):


    >>> n = -20
    >>> nh = fadd(n, 1e-100, exact=True)
    >>> hyp0f1(n, 0.25)    # nonsense
    0.987581857511351
    >>> hyp0f1(nh, 0.25)   # same nonsense
    0.987581857511351

Now, correctly:


    >>> n = -20
    >>> nh = fadd(n, 1e-100, exact=True)
    >>> hyp0f1(n, 0.25)
    Traceback (most recent call last):
      ...
    ZeroDivisionError: pole in hypergeometric series
    >>> hyp0f1(nh, 0.25)
    1.85014429040103e+49

Finally, and probably most importantly, the limit evaluation code has been changed to adaptively decrease the size of perturbation until convergence. Under fairly general assumptions, the maximum accurate perturbation at precision p can easily be shown to be 2^−(p+C); the problem is that the parameter-dependent constant C isn't generally known. Previously C was just set to a small value (10), and naturally this would break down for some functions when parameters were increased beyond roughly that magnitude.

I briefly considered the idea of estimating C analytically in terms of the parameters, and while I think this can be done rigorously, it seems difficult -- especially to do it tightly (grossly overestimating C would murder performance). The algorithm implemented now is quasi-rigorous, and although there is some slowdown (sometimes by a fair amount), the improved reliability is definitely worth it. A user can also manually supply a perturbation size of their own taste, thereby overriding adaptivity. If the user supplies an intelligent value, this gives both the speed of the old code and full rigor. Probably some intelligent choices for particular functions could be made automatically by mpmath too, to recover the old speed in common cases.

An example of a situation where this becomes necessary is for Legendre functions with certain large parameter combinations. With the old code:


    >>> mp.dps = 15; legenp(0.5, 300, 0.25)   # bad
    4.19317029788723e+626
    >>> mp.dps = 30; legenp(0.5, 300, 0.25)   # bad
    3.72428336871098039162972441784e+611
    >>> mp.dps = 60; legenp(0.5, 300, 0.25)   # bad
    2.93794154954090326636196697693611381787845107728298382876544e+581
    >>> mp.dps = 100; legenp(0.5, 300, 0.25)   # only accurate to a few digits
    -1.71750854949725238788826203712778687036438365374945625996246145924802366061559
    6579520831362887006032e+578

With the new code:


    >>> mp.dps = 15; legenp(0.5, 300, 0.25)
    -1.71750854949725e+578
    >>> mp.dps = 30; legenp(0.5, 300, 0.25)
    -1.71750854949725238788826203713e+578
    >>> mp.dps = 60; legenp(0.5, 300, 0.25)
    -1.71750854949725238788826203712778687063419097363472262860567e+578
    >>> mp.dps = 15; legenp(0.5, 300, 0.25, hmag=300)   # faster, with manual perturbation
    -1.71750854949725e+578

Another example is for Bessel functions differentiated to high order. With the old code:


    >>> mp.dps = 15; besseli(2, 10, derivative=100)   # bad
    2.63560662943646e+34
    >>> mp.dps = 30; besseli(2, 10, derivative=100)   # bad
    23408889310840499424.9813614712
    >>> mp.dps = 60; besseli(2, 10, derivative=100)   # only accurate to a few digits
    821.238621006501815018537509753672810563116338269219460709828

With the new code:


    >>> mp.dps = 15; besseli(2, 10, derivative=100)
    821.238621006483
    >>> mp.dps = 30; besseli(2, 10, derivative=100)
    821.238621006483348660925541651
    >>> mp.dps = 60; besseli(2, 10, derivative=100)
    821.238621006483348660925541650744338655411830854860813048862

In related news, I've given all hypergeometric-type functions the ability to accept kwargs and pass them on to hypercomb and hypsum. This means that tuning and error control options will be available for a large number of functions (Bessel functions, etc), which could be useful for some users who need to do advanced calculations. I have yet to document these features thoroughly (the interface will also perhaps be tweaked before the next release).

Analytic continuation of 3F2, 4F3 and higher functions

2009-12-22T17:01:00.008+01:00

As of a recent commit, mpmath can evaluate the analytic continuation of the generalized hypergeometric function _p+1F_p for any p. Previously ₂F₁ (the Gaussian hypergeometric function) was supported -- see earlier blog posts -- but not ₃F₂ and higher. This addition means that the generalized hypergeometric function is finally supported essentially everywhere where it is "well-posed" (in the sense that the series has nonzero radius of convergence), so it is a rather significant improvement. Unfortunately, the implementation is still not perfect, but I decided to commit the existing code since it is quite useful already (and long overdue).

As proof of operation, I deliver plots of ₃F₂ and ₄F₃, requiring both |z| < 1 and |z| > 1:


from mpmath import *
f1 = lambda z: hyp3f2(1,2,3,4,5,z)
f2 = lambda z: hyper([1,2,3,4],[5,6,7],z)
plot([f1,f2], [-2,2])

A portrait of ₃F₂ restricted to the unit circle:


plot(lambda x: hyp3f2(1,2,3,4,5,exp(2*pi*j*x)), [-1,1])

A numerical value of ₅F₄ with z on the unit circle, with general complex parameters:


>>> mp.dps = 50
>>> print hyper([1+j,0.5,-2j,1,0.5+0.5j],[0.5j,0.25,-j,-1-j],-1)
(-1.8419705729324526212110109087877199070037836117341 - 
0.57799141426978758652682670969879368618437786161612j)

High precision values of ₃F₂ at z = 1 and z = 1 + ε:


>>> print hyp3f2(1,1,2,3.5,1.5,1)
2.2603241503205748814116375624327119385797950825522
>>> print hyp3f2(1,1,2,3.5,1.5,'1.0001')
(2.2626812356987790952469291649495098300894035980837 -
0.00092505569713084902188662960467216683326695892509251j)

A complex plot of ₃F₂:


>>> mp.dps = 5
>>> cplot(lambda z: hyp3f2(2.5,3,4,1,2.25,z), [-2,2], [-2,2],
...     points=50000)

A bit of theoretical background: the hypergeometric series _pF_q has an infinite radius of convergence when p ≤ q, so it can in principle be evaluated then by adding sufficiently many terms at sufficiently high precision (although in practice asymptotic expansions must be used for large arguments, as mpmath does). In the balanced case when p = q+1, i.e. for ₁F₀(a,z), ₂F₁(a,b,c,z), ₃F₂(...), ₄F₃(...), ..., the series converges only when |z| < 1 (plus in some special instances). This is due to the fact that the hypergeometric function has a singularity (a pole or branch point, depending on parameters) at z = 1, in turn owing to the fact that the hypergeometric differential equation is singular at z = 1. (The reason that the p = q+1 and not p = q case is "balanced" is that there is an extra, implicit factorial in the hypergeometric series.)

The main ingredient of the analytic continuation of _p+1F_p is the inversion formula which replaces z with 1/z and thus handles |z| > 1. This was easy to implement -- the only complication is that integer parameters result in singular gamma factors, but the mechanisms to handle those automatically were already in place. There was no particular reason why I hadn't added that code already.

The tricky part is the unit circle, and the close vicinity thereof, where neither series converges (quickly). I've been looking for good ways handle this, with mixed results.

When I posed the problem on this blog several months back, a reader suggested this paper by Wolfgang Bühring which gives a series for the analytic continuation around z = 1. My finding from trying to implement it is that the rate of convergence of the series generally is poor, and therefore it is not immediately effective for computation. However, convergence acceleration with nsum improves the situation considerably in many cases. Some parameter combinations render the convergence acceleration useless, but even then, it can give a few correct digits, so it is better than nothing (although the implementation should probably warn the user when the result probably is inaccurate). I'm unfortunately not aware of any parameter transformations that would substantially improve convergence. The current implementation uses this method for ₃F₂ when |z-1| is small; it should work for 4F3 and higher too, but the series coefficients are much more complicated (involving multiply nested sums), so that's yet to be done.

For the rest of the unit circle, I've settled for simply using convergence acceleration directly. This essentially just amounts to passing the hypergeometric series to nsum, which applies Richardson extrapolation and iterated Shanks transformations. The Shanks transformation is actually perfect for this -- it's almost tailor made for convergence acceleration of the balanced hypergeometric series -- and is able to sum the series even outside the circle of convergence, using just a few terms. This covers most of the unit circle -- the catch is just that the acceleration asymptotically deteriorates and ultimately becomes useless close to z = 1, so complementary method (such as the Bühring's is still required there).

Mathematica seems to support unit-circle evaluation of hypergeometric functions quite well. Unfortunately, I don't know how it does it. According to Wolfram's Some Notes On Internal Implementation page,

The hypergeometric functions use functional equations, stable recurrence relations, series expansions, asymptotic series and Padé approximants. Methods from NSum and NIntegrate are also sometimes used.

This looks similar to what I'm doing -- high order Padé approximants and Mathematica's NSum should be equivalent to the nsum-based series acceleration in mpmath. But probably Mathematica employs some more tricks as well.

I've also tested direct integration of the hypergeometric differential equation and of Mellin-Barnes integral representations, but these approaches don't seem to work much better (at least not without many improvments), and at best seem to be relatively slow.

Python floats and other unusual things spotted in mpmath

2009-09-10T19:29:00.007+02:00

I've just put up a branch (named "mp4") containing changes to mpmath that I've been working on for several days now. This branch includes some rather large changes, including support for fixed-precision (machine-precision) arithmetic, and a new implementation of the multiprecision type. The code is a bit rough at the moment, and not all changes are final, so don't expect this in a release in the immediate future.

New mpf implementation

The mp4 branch uses a new mpf type that handles both real and complex values. It makes complex arithmetic quite a bit faster (up to twice the speed), although real arithmetic is a little slower. Some of the slowdown is due to wrapping old-format functions instead of providing new ones, so there is a bit of unnecessary overhead. The biggest advantage is not that the implementation becomes faster, but that it becomes simpler. For example, a nice feature is that complex calculations that yield real results don't need to be converted back to a real type manually. Unfortunately, writing code that also supports Python float and complex (see below) somewhat reduce the usefulness of such features.

Imaginary infinities are a little more well-behaved (i.e. j*(j*inf) gives -inf and not a complex NaN), and the representation is flexible enough that it might be possible to support unsigned infinities, infinities with an arbitrary complex sign, and possibly even zero with an arbitrary complex sign -- hey, perhaps infinities and zeros with order and residue information attached to them so you can evaluate gamma(-3)**2 / gamma(-5) / gamma(-2) directly? (Well, I'm not sure if this is actually possible, and if it is, it's probably way overkill :-)

This new implementation will not necessarily make it -- it depends on whether the pros outweigh the cons. I would very much appreciate review by others. Regardless of the outcome, writing it has been very useful for identifying and fixing problems in the interface between "low level" and "high level" code in mpmath.

This should particularly simplify the future support for a C or Cython-based backend. In fact, the current fixed-precision backend (see below) only uses about 200 lines of code -- and although it's not complete yet, it handles a large set of calculations. Adding a fast multiprecision backend based on Sage or GMPY can probably be done in an evening now. I'm not going to do such a thing right now, as I want to fix up the existing code before adding anything new; should someone else be interested though, then by all means go ahead.

I've also inserted a few other optimizations. As noted before on this blog, many of the special functions in mpmath are evaluated as linear combinations of hypergeometric series. All those functions basically fall back to a single routine which sums hypergeometric series. The mp4 branch contains a new implementation of this routine that is up to about twice as fast as before. The speed is achieved by code generation: for every different "type" (degree, parameter types) of hypergeometric series, a specialized version is generated with various bookkeeping done in advance so it doesn't have to be repeated in every iteration of the inner loop.

As an example, the following runs at 1.7x the speed:


>>> 1/timing(mpmath.hyp1f1, 1.5, 2.2, 1.35)
5620.8844813722862
>>> 1/timing(mp4.hyp1f1, 1.5, 2.2, 1.35)
9653.1737629459149

Nearly all tests pass with the new mpf type used instead of the old one -- the only significant missing piece is interval arithmetic.

Fixed-precision arithmetic

Several people have requested support for working with regular Python floats and complexes in mpmath. Often you only need a few digits of accuracy, and mpmath's multiprecision arithmetic is unnecessarily slow. Many (though not all) of the algorithms in mpmath work well in fixed precision; the main problem with supporting this feature has been that of providing an appropriate interface that avoids cluttering the code.

In the mp4 branch, this is solved by adding a fixed-precision context. This context uses the same "high-level" methods for special functions, calculus, linear algebra, etc, as the multiprecision context, while low-level functions are replaced with fixed-precision versions. For example exp is just a wrapper around math.exp and cmath.exp.

While the default multiprecision context instance is called mp, the default fixed-precision context instance is called fp. So:


>>> from mp4 import mp, fp
>>> mp.pretty = True
>>> mp.sqrt(-5); type(_)
2.23606797749979j
<class 'mp4.ctx_mp.mpf'>
>>> fp.sqrt(-5); type(_)
2.2360679774997898j
<type 'complex'>

The fixed-precision context is still missing a lot of low-level functions, so many things don't work yet. Let's try a couple of calculations that do work and see how they compare.

Lambert W function

The Lambert W function was the first nontrivial function I tried to get working, since it's very useful in fixed precision and its calculation only requires simple functions. The lambertw method is the same for mp as for fp; the latter just uses float or complex for the arithmetic.


>>> mp.lambertw(7.5)
1.56623095378239
>>> mp.lambertw(3+4j)
(1.28156180612378 + 0.533095222020971j)
>>> fp.lambertw(7.5)
1.5662309537823875
>>> fp.lambertw(3+4j)
(1.2815618061237759+0.53309522202097104j)

The fixed-precision results are very accurate, which is not surprising since the Lambert W function is implemented using a "self-correcting" convergent iteration. In fact, the multiprecision implementation could be sped up by using the fixed-precision version to generate the initial value. The speed difference is quite striking:


>>> 1/timing(mp.lambertw, 7.5)
1249.5319808144905
>>> 1/timing(mp.lambertw, 3+4j)
603.49697841726618
>>> 1/timing(fp.lambertw, 7.5)
36095.559380378654
>>> 1/timing(fp.lambertw, 3+4j)
20281.934235976787

Both the real and complex versions are about 30x faster in fixed precision.

Hurwitz zeta function

The Hurwitz zeta function is implemented mainly using Euler-Maclaurin summation. The main ingredients are Bernoulli numbers and the powers in the truncated L-series. Bernoulli numbers are cached, and can be computed relatively quickly to begin with, so they're not much to worry about. In fact, I implemented fixed-precision Bernoulli numbers by wrapping the arbitrary-precision routine for them, so they are available to full 53-bit accuracy. As it turns out, the fixed-precision evaluation achieves nearly full accuracy:


>>> mp.hurwitz(3.5, 2.25); fp.hurwitz(3.5, 2.25)
0.0890122424923889
0.089012242492388816
>>> t1=timing(mp.hurwitz,3.5,2.25); t2=timing(fp.hurwitz,3.5,2.25); 1/t1; 1/t2; t1/t2
473.66504799548278
7008.0267335004173
14.7953216374269

There is a nice 15x speedup, and it gets even better if we try complex values. Let's evaluate the the zeta function on the critical line 0.5+ti for increasing values of t:


>>> s=0.5+10j
>>> mp.hurwitz(s); fp.hurwitz(s)
(1.54489522029675 - 0.115336465271273j)
(1.5448952202967554-0.11533646527127067j)
>>> t1=timing(mp.hurwitz,s); t2=timing(fp.hurwitz,s); 1/t1; 1/t2; t1/t2
213.37891598750548
4391.4815202596583
20.580672180923461

>>> s=0.5+1000j
>>> mp.hurwitz(s); fp.hurwitz(s)
(0.356334367194396 + 0.931997831232994j)
(0.35633436719476846+0.93199783123336344j)
>>> t1=timing(mp.hurwitz,s); t2=timing(fp.hurwitz,s); 1/t1; 1/t2; t1/t2
8.1703453931669383
352.54843617352128
43.149759185011469

>>> s = 0.5+100000j
>>> mp.hurwitz(s); fp.hurwitz(s)
(1.07303201485775 + 5.7808485443635j)
(1.0730320148426455+5.7808485443942352j)
>>> t1=timing(mp.hurwitz,s); t2=timing(fp.hurwitz,s); 1/t1; 1/t2; t1/t2
0.17125208455924443
9.2754935956407891
54.162806949260492

It's definitely possible to go up to slightly larger heights still. The Euler-Maclaurin truncation in the Hurwitz zeta implementation is not really tuned, and certainly has not been tuned for fixed precision, so the speed can probably be improved.

Hypergeometric functions

Implementing hypergeometric series in fixed precision was trivial. That the multi-precision implementation is actually quite fast can be seen by comparing direct evaluations:


>>> mp.hyp1f1(2.5, 1.2, -4.5)
-0.0164674858506064
>>> fp.hyp1f1(2.5, 1.2, -4.5)
-0.016467485850591584
>>> 1/timing(mp.hyp1f1, 2.5, 1.2, -4.5)
6707.6667199744124
>>> 1/timing(fp.hyp1f1, 2.5, 1.2, -4.5)
12049.135305946565

The float version is only about twice as fast. Unfortunately, hypergeometric series suffer from catastrophic cancellation in fixed precision, as can be seen by trying a larger argument:


>>> mp.hyp1f1(2.5, 1.2, -30.5)
6.62762709628679e-5
>>> fp.hyp1f1(2.5, 1.2, -30.5)
-0.012819333651375751

Potentially, checks could be added so that the fixed-precision series raises an exception or falls back to arbitrary-precision arithmetic internally when catastrophic cancellation occurs. However, it turns out that this evaluation works for even larger arguments when the numerically stable asymptotic series kicks in:


>>> mp.hyp1f1(2.5, 1.2, -300.5)
1.79669541078302e-7
>>> fp.hyp1f1(2.5, 1.2, -300.5+0j)
1.7966954107830195e-07

The reason I used a complex argument is that the asymptotic series uses complex arithmetic internally, and Python has an annoying habit of raising exceptions when performing complex-valued operations involving floats (a proper workaround will have to be added). In this case the speedup is close to an order of magnitude:


>>> 1/timing(mp.hyp1f1, 2.5, 1.2, -300.5)
532.36666412814452
>>> 1/timing(fp.hyp1f1, 2.5, 1.2, -300.5+0j)
4629.4746136865342

Another example, a Bessel function. This function is calculated using a hypergeometric series and also a couple of additional factors, so the speedup is quite large:


>>> mp.besselk(2.5, 7.5)
0.000367862846522012
>>> fp.besselk(2.5, 7.5)
0.00036786284652201188
>>> 1/timing(mp.besselk, 2.5, 7.5)
3410.5578142787444
>>> 1/timing(fp.besselk, 2.5, 7.5)
21879.520083463747

The fixed-precision context does not yet handle cancellation of singularities in hypergeometric functions. This can be implemented as in the multiprecision case by perturbing values, although accuracy will be quite low (at best 5-6 digits; sometimes an accurate result will be impossible to obtain).

Numerical calculus

Root-finding works, as does numerical integration, with nice speedups:


>>> fp.quad(lambda x: fp.cos(x), [2, 5])
-1.8682217014888205
>>> mp.quad(lambda x: mp.cos(x), [2, 5])
-1.86822170148882
>>> 1/timing(fp.quad, lambda x: fp.cos(x), [2, 5])
1368.3622602114056
>>> 1/timing(mp.quad, lambda x: mp.cos(x), [2, 5])
55.480651962846274

>>> fp.findroot(lambda x: fp.cos(x), 2)
1.5707963267948966
>>> mp.findroot(lambda x: mp.cos(x), 2)
1.5707963267949
>>> 1/timing(fp.findroot, lambda x: fp.cos(x), 2)
19160.822293284604
>>> 1/timing(mp.findroot, lambda x: mp.cos(x), 2)
1039.4805452292442

The tests above use well-behaved object functions; some corner cases are likely fragile at this point. I also know, without having tried, that many other calculus functions utterly don't work in fixed precision (not by algorithm, nor by implementation). Some work will be needed to support them even partially. At minimum, several functions will have to be changed to use an epsilon of 10^-5 or so since full 15-16-digit accuracy requires extra working precision which just isn't available.

Plotting

The plot methods work, and are of course way faster in fixed-precision mode. For your enjoyment, here is the gamma function in high resolution; the following image took only a few seconds to generate:


>>> fp.cplot(fp.gamma, points=400000)

With the default number of points (fp.cplot(fp.gamma)), it's of course instantaneous.

Released: mpmath 0.13

2009-08-13T21:49:00.004+02:00

I released mpmath 0.13 today. See the announcement for details. Soon coming to a Sage or SymPy near you!

I've blogged extensively about the new features that went into this version over the last two months. The short version of the changelog is "lots of special functions". For convenience, here a list of all posts:

Torture testing special functions

3D visualization of complex functions with matplotlib

Coulomb wave functions and orthogonal polynomials

Hurwitz zeta function, Dirichlet L-series, Appell F1

Another Mathematica bug

Improved incomplete gamma and exponential integrals; Clausen functions

Meijer G, more hypergeometric functions, fractional differentiation

Massive hypergeometric update

Hypergeometric 2F1, incomplete beta, exponential integrals

I also suggest a look at the documentation and particularly the special functions section.

Torture testing special functions

2009-08-11T02:27:00.007+02:00

Today I implemented a set of torture test for special functions in mpmath.

The function torturer takes a function, say f = airyai or f = lambda z: hyp1f1(1.5, -2.25, z), and attempts to evaluate it with z = a · 10ⁿ as n ranges between -20 and 20, with a ≈ 1, 1+i, i, -1+i, -1. (The "≈" indicates multiplication by an additional quasi-random factor near unity.)

Thus it tests a wide range of magnitudes (both tiny and huge), with pure positive, negative, imaginary, and complex arguments; i.e. covering most distinct regions of the complex plane. (It doesn't cover the lower half-plane, but that's generally just a waste of time since most function algorithms are agnostic about the sign of the imaginary part -- this should be done for some functions in the future, however.)

Each evaluation is also performed at a range of precisions, between 5 and 150 digits by default (for fast functions the max precision is set much higher). The results at two successive precisions are compared to verify that they agree relatively to the lesser of the precisions (with a single-digit error allowed for roundoff error).

This doesn't guarantee correctness, of course -- I might have implemented the wrong formula for the function -- but it's a strong check that whatever formula has been implemented is being evaluated accurately. It does provide some amount of correctness testing for those functions that use more than one algorithm depending on precision and/or the magnitude of the argument (and this includes the hypergeometric functions). If one algorithm is used at low precision and another at high precision, then an error in either will be revealed when comparing results.

My original intention with was just to check robustness of the asymptotic expansions of hypergeometric functions. Testing with a continuous range of magnitudes ensures that there aren't any regions where convergence gets impossibly slow, a loop fails to terminate, an error estimate is insufficient, etc. (And this needs to be done at several different levels of precision.) I then ended up including many other functions in the list of torture tests as well.

The results? I discovered (and fixed!) some 10-20 bugs in mpmath. Most were of the kind "only 43 accurate digits were returned where 50 were expected" and due to inadequate error estimates and/or insufficient extra precision being allocated (mostly in argument transformations). I actually knew about many of the problems, and had just neglected to fix them because they wouldn't show up for "normal" arguments (and in particular, there weren't any existing tests that failed).

Lesson to be learned: if it isn't tested, it probably doesn't work (fully). Semi-automatic exhaustive or randomized testing (the more "evil" inputs, the better), is a necessary complement to testing hand-picked values.

There is more to be done:

Some functions don't yet pass the torture tests

Rather than just large/small, the relevant stress measure for some functions is closeness to a negative integer, closeness to the unit circle, etc

Functions of more than one argument should be tested with randomly selected values for all parameters

But already now, the torture tests have done a good job. Tellingly, functions I implemented recently are much less buggy than functions I implemented long ago (even just three months ago). Periodically rewriting old code as one learns more is bound to improve quality :-)

Speaking of bugs, Mathematica doesn't pass the mpmath torture test suite (far from it). An example of something that will fail:


In[47]:= N[HypergeometricPFQ[{1},{-1/4,1/3},10^10],15]

General::ovfl: Overflow occurred in computation.

General::ovfl: Overflow occurred in computation.

General::ovfl: Overflow occurred in computation.

General::stop: Further output of General::ovfl
     will be suppressed during this calculation.

(... nothing for 5 minutes)

Interrupt> a

Out[47]= $Aborted

And mpmath (instantaneously):


>>> mp.dps = 15; mp.pretty = True
>>> hyp1f2(1,(-1,4),(1,3),10**10)
-3.53521222026601e+86866

By the way, the current set of torture tests takes about 11 minutes to run on my computer with psyco and gmpy (much longer without either). It's possibly a better benchmark than the standard tests.

Update: running the torture tests with multiprocessing on sage.math takes only 28 seconds (which is the time for the slowest individual test). Go parallel!

3D visualization of complex functions with matplotlib

2009-08-09T19:33:00.011+02:00

Hooray! matplotlib 0.99 is out and it has 3D plotting, finally!

I've shown a lot of color plots of complex functions on this blog to demonstrate complex functions in mpmath. These plots are informative, but sometimes a 3D plot (typically of the function's absolute value) gives a much better view. A big advantage of 3D plots over 2D color plots is that far fewer evaluation points are required for a good high-resolution image, and this helps when visualizing the slower functions in mpmath.

There will probably be a 3D plot function in a future version of mpmath (or two functions; for two-variable real, and complex functions), similar in style to the existing matplotlib wrappers plot and cplot. Until I've figured out the details, I'll share a couple of test plots.

Coulomb wave function of a complex argument, F(6,4,z):

Principal-branch logarithmic gamma function:

Gamma function:

Imaginary part of Lambert W function, 0th branch:

Riemann zeta function in the critical strip:

Those are all wireframe plots. It's also possible to do surface plots. Using the Riemann zeta function again, a surface plot looks as follows:

Surface + wireframe plot:

None of these is perfect. I'd like to be able to do a surface plot with a widely spaced wireframe mesh to pronounce the geometry, and smooth colored surface between the meshes. This doesn't seem possible with matplotlib because the surface plot doesn't do smoothing (even with a stride selected); overlaying a wireframe works decently, but the wireframe doesn't render with occlusion and this looks bad for some functions. Since the pure wireframe plot is much faster, I think I prefer it for now.

For complex functions, it'd also be nice with a color function separate from the geometry function -- then you could plot the phase as color in the same picture.

Color helps for visualizing complicated structure, especially for phase plots. Below, I've plotted the phase (actually the sine of the phase, to make it continuous) of a Jacobi theta function θ₃(1+4j/3,q) (restricted to |q| < 0.8, because it gets extremely messy for larger |q|):

The phase of the Barnes G-function:

I could do many more, but that will probably enough for this blog :-)

To reproduce any of these, use the following script with trivial modifications. It's a straightforward extension of the example scripts in the matplotlib documentation.


from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import pylab
import numpy as np
import mpmath
mpmath.dps = 5

# Use instead of arg for a continuous phase
def arg2(x):
    return mpmath.sin(mpmath.arg(x))

#f = lambda z: abs(mpmath.loggamma(z))
#f = lambda z: arg2(mpmath.exp(z))
#f = lambda z: abs(mpmath.besselj(3,z))
f = lambda z: arg2(mpmath.cos(z))

fig = pylab.figure()
ax = Axes3D(fig)
X = np.arange(-5, 5, 0.125)
Y = np.arange(-5, 5, 0.125)
X, Y = np.meshgrid(X, Y)
xn, yn = X.shape
W = X*0
for xk in range(xn):
    for yk in range(yn):
        try:
            z = complex(X[xk,yk],Y[xk,yk])
            w = float(f(z))
            if w != w:
                raise ValueError
            W[xk,yk] = w
        except (ValueError, TypeError, ZeroDivisionError):
            # can handle special values here
            pass
    print xk, xn

# can comment out one of these
ax.plot_surface(X, Y, W, rstride=1, cstride=1, cmap=cm.jet)
ax.plot_wireframe(X, Y, W, rstride=5, cstride=5)

pylab.show()

Coulomb wave functions and orthogonal polynomials

2009-08-04T22:38:00.008+02:00

Work has been a bit slow over the last two weeks, partially due to the fact that I was away during one of them. Nevertheless, I've been able to add a couple of more functions to mpmath, as described in the following post below. With these functions committed, I'm probably going to stop adding features and just do a new release as soon as possible (some documentation and general cleanup will be required first).

Coulomb wave functions

I have, due to a request, implemented the Coulomb wave functions (commit). The Coulomb wave functions are used in quantum physics; they solve the radial part of the Schrödinger equation for a particle in a 1/r potential.

The versions in mpmath are fully general. They work not only for an arbitrarily large or small radius, but for complex values of all parameters. Some example evaluations:


>>> mp.dps = 25; mp.pretty = True
>>> coulombf(4, 2.5, 1000000)
-0.9713831935409625197404938
>>> coulombf(2+3j, 1+j, 100j)
(2.161557009297068729111076e+37 + 1.217455850269101085622437e+39j)
>>> coulombg(-3.5, 2, '1e-100')
2.746460651456730226435998e+252

The definitions of the Coulomb wave functions for complex parameters are based mostly on this paper by N. Michel, which describes a special-purpose C++ implementation. I'm not aware of any other software that implements Coulomb wave functions, except for GSL which only supports limited domains in fixed precision.

The Coulomb wave functions are numerically difficult to calculate: the canonical representation requires complex arithmetic, and the terms vary by many orders of magnitude even for small parameter values. Arbitrary-precision arithmetic, therefore, is almost necessary even to obtain low-precision values, unless one uses much more specialized evaluation methods. The current mpmath versions are simple and robust, although they can be fairly slow in worst cases (i.e. when they need to use a high internal precision). Timings for average/good cases are in the millisecond range,


>>> timing(coulombf, 3, 2, 10)
0.0012244852348553437
>>> timing(coulombg, 3, 2, 10)
0.0068572681723514609

but they can climb up towards a second or so when large cancellations occur. Someone who needs faster Coulomb wave functions, say for a limited range of parameters, could perhaps find mpmath useful for testing a more specialized implementation against.

The speed is plenty for basic visualization purposes. Below I've recreated graphs 14.1, 14.2 and 14.5 from Abramowitz & Stegun.

mp.dps = 5

# Figure 14.1
F1 = lambda x: coulombf(x,1,10)
F2 = lambda x: coulombg(x,1,10)
plot([F1,F2], [0,15], [-1.2, 1.7])

# Figure 14.3
F1 = lambda x: coulombf(0,0,x)
F2 = lambda x: coulombf(0,1,x)
F3 = lambda x: coulombf(0,5,x)
F4 = lambda x: coulombf(0,10,x)
F5 = lambda x: coulombf(0,x/2,x)
plot([F1,F2,F3,F4,F5], [0,25], [-1.2,1.6])

# Figure 14.5
F1 = lambda x: coulombg(0,0,x)
F2 = lambda x: coulombg(0,1,x)
F3 = lambda x: coulombg(0,5,x)
F4 = lambda x: coulombg(0,10,x)
F5 = lambda x: coulombg(0,x/2,x)
plot([F1,F2,F3,F4,F5], [0,30], [-2,2])

Also, a plot for complex values, cplot(lambda z: coulombg(1+j, -0.5, z)):

Orthogonal polynomials

I've added the (associated) Legendre functions of the first and second kind, Gegenbauer, (associated) Laguerre and Hermite functions (commits 1, 2, 3, 4). This means that mpmath finally supports all the classical orthogonal polynomials. As usual, they work in generalized form, for complex values of all parameters.

A fun thing to do is to verify orthogonality of the orthogonal polynomials using numerical quadrature:


>>> mp.dps = 15
>>> chop(quad(lambda t: exp(-t)*t**2*laguerre(3,2,t)*laguerre(4,2,t), [0,inf]))
0.0
>>> chop(quad(lambda t: exp(-t**2)*hermite(3,t)*hermite(4,t), [-inf,inf]))
0.0

As another basic demonstration, here are some Legendre functions of the second kind, visualized:

F1 = lambda x: legenq(0,0,x)
F2 = lambda x: legenq(1,0,x)
F3 = lambda x: legenq(2,0,x)
F4 = lambda x: legenq(3,0,x)
F5 = lambda x: legenq(4,0,x)
plot([F1,F2,F3,F4,F5], [-1,1], [-1.3,1.3])

Implementing these functions is a lot of work, I've found. The general case is simple (just fall back to generic hypergeometric code), but the special cases (singularities, limits, asymptotically large arguments) are easy to get wrong and there are lots of them to test and document.

My general methodology is to implement the special functions as the Wolfram Functions site defines them (if they are listed there), test my implementation against exact formulas, and then test numerical values against Mathematica. Unfortunately, Wolfram Functions often leaves limits undefined, and is not always consistent with Mathematica. In fact, Mathematica is not even consistent with itself. Consider the following:


In[277]:= LaguerreL[-1, b, z] // FunctionExpand

Out[277]= 0

In[278]:= LaguerreL[-1, -2, z] // FunctionExpand

           z  2
          E  z
Out[278]= -----
            2

It's tempting to just leave such a case undefined in mpmath and move on (and perhaps fix it later if it turns out to be used). Ideally mpmath should handle all singular cases correctly and consistently, and it should state in the documentation precisely what it will calculate for any given input so that the user doesn't have to guess. Testing and documenting special cases is very time-consuming, however, and although I'm making progress, this work is far from complete.

Hurwitz zeta function, Dirichlet L-series, Appell F1

2009-07-15T00:02:00.013+02:00

I've added three more functions to mpmath since the last blog update: the Hurwitz zeta function, Dirichlet L-series and the Appell hypergeometric function.

Hurwitz zeta function

The Hurwitz zeta function is available with the syntax hurwitz(s, a=1, derivative=0). It's a separate function from zeta for various reasons (one being to make the plain zeta function as simple as possible). With the optional third argument, it not only computes values, but arbitrarily high derivatives with respect to s (more on which later).

It's quite fast at high precision; e.g. the following evaluation takes just 0.25 seconds (0.6 seconds from cold cache) on my laptop:


>>> mp.dps = 1000
>>> hurwitz(3, (1,3))        # (1,3) specifies the exact parameter 1/3
27.56106119970080377622787797740750928454209531301488108286446031855695818017193
31971447814492522557242833928839896347960163079006102651117997856881812399406784
05731893349475372409101251211346682545314192272823745384665266049022078198043979
99746337871597509019667252833751410875882452970087673701166818534638337151352084
77353765613742227179887954515352387981214398169985388076469409432688912355506092
00433349370785407936365607196038965500743780377965599471761097326271922616743812
30503180321801486524844073873393780084553765789276561564723000726671702263215399
29351207878018547946361756911396554507918933344824138599697851666075594487957278
34823450478752821675163024463415647205608959739100071026336912826930683261155867
06265136825191025873785920661492516925137079914816949085911457084632783280741306
98897458344746469272828267484989756595259917290800848768074029766417786815732919
07917349747600320626141004448429987450854391175868138793015691343081417209625034
88644344900903488585429308918210445984048

With s = π instead, it takes 3.5 seconds. For comparison, Mathematica 6 took 0.5 and 8.3 seconds respectively (on an unloaded remote system which is faster than my laptop, but I'm not going to guess by how much).

The function is a bit slow at low precision, relatively speaking, but fast enough (in the 1-10 millisecond range) for plotting and such. Of course, it works for complex s, and in particular on the critical line 1/2+it. Below, I've plotted |ζ(1/2+it,1)|², |ζ(1/2+it,1/2)|², |ζ(1/2+it,1/3)|²; the first image shows 0 ≤ t < 30 and the second is with 1000 ≤ t < 1030.

To show some more complex values, here are plots of ζ(s, 1), ζ(s, 1/3) and ζ(s, 24/25) for 100 ≤ Im(s) ≤ 110. I picked the range [-2, 3] for the real part to show that the reflection formula for Re(s) < 0 works:

In fact, the Hurwitz zeta implementation works for complex values of both s and a. Here is a plot of hurwitz(3+4j, a) with respect to a:

The evaluation can be slow and/or inaccurate for nonrational values of a, though (in nonconvergent cases where the functional equation isn't used).

Zeta function performance

Right now the implementations of the Riemann zeta function and the Hurwitz zeta function in mpmath are entirely separate. In fact, they even use entirely different algorithms (hurwitz uses Euler-Maclaurin summation, while zeta uses the Borwein approximation). This is useful for verifying correctness of either function.

Regarding performance, it appears that the Borwein approximation is slightly faster for small |Im(s)| while Euler-Maclaurin is massively faster for large |s|.

The following benchmark might give an idea:


>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> timing(zeta, 0.5+10j); timing(hurwitz, 0.5+10j)
0.0027386162207766627
0.0082901877193889192
>>> timing(zeta, 0.5+100j); timing(hurwitz, 0.5+100j)
0.0075319349246901089
0.041037198861866478
>>> timing(zeta, 0.5+1000j); timing(hurwitz, 0.5+1000j)
0.13468499488063479
0.18062463246027072
>>> timing(zeta, 0.5+10000j); timing(hurwitz, 0.5+10000j)
7.2977593520964241
1.0488757649366072
>>> timing(zeta, 0.5+10000j); timing(hurwitz, 0.5+10000j)
0.84306636737350971
1.0601629536714867

As the last evaluation shows, the Borwein algorithm is still not doing too poorly at 0.5+10000j from a warm cache, but at this point the cache is 50 MB large. At 0.5+100000j my laptop ran out of RAM and Firefox crashed (fortunately Blogger autosaves!) so it clearly doesn't scale. The Euler-Maclaurin algorithm caches a lot of Bernoulli numbers, but the memory consumption of this appears to be negligible. With hurwitz (and the function value, for those interested),


>>> timing(hurwitz, 0.5+100000j)
8.0361576680835149
>>> hurwitz(0.5+100000j)
(1.07303201485775 + 5.7808485443635j)

and the memory usage of the Python interpreter process only climbs to 12 MB.

Just by optimizing the pure Python code, I think the Hurwitz zeta function (at low precision) could be improved by about a factor 2 generally and perhaps by a factor 3-4 on the critical line (where square roots can be used instead of logarithms -- zeta uses this optimization but hurwitz doesn't yet); with Cython, perhaps by as much as a factor 10.

The real way to performance is not to use arbitrary-precision arithmetic, though. Euler-Maclaurin summation for zeta functions is remarkably stable in fixed-precision arithmetic, so there is no problem using doubles for most applications. As I wrote a while back on sage-devel, a preliminary version of my Hurwitz zeta code for Python complex was 5x faster than Sage's CDF zeta (in a single test, mind you). If there is interest, I could add such a version, perhaps writing it in Cython for Sage (for even greater speed).

Derivatives

The Hurwitz zeta function can be differentiated easily like so:


>>> mp.dps = 25
>>> hurwitz(2+3j, 1.25, 3)
(-0.01266157985314340398338543 - 0.06298953579777517606036907j)
>>> diff(lambda s: hurwitz(s, 1.25), 2+3j, 3)
(-0.01266157985314340398338543 - 0.06298953579777517606036907j)

For simple applications one can just as well use the numerical diff function as above with reasonable performance, but this isn't really feasible at high precision and/or high orders (numerically computing an nth derivative at precision p requires n+1 function evaluations, each at precision (n+1)×p).

The specialized code in hurwitz requires no extra precision (although it still needs to perform extra work) and therefore scales up to very high precision and high orders. Here for example is a derivative of order 100 (taking 0.5 seconds) and one of order 200 (taking 4 seconds):


>>> hurwitz(2+3j, 1.25, 100)
(2.604086240825183469410664e+107 - 1.388710675207202633247271e+107j)
>>> hurwitz(2+3j, 1.25, 200)
(2.404124816484789309285653e+274 + 6.633220818921777955999455e+273j)

It should be possible to calculate a sequence of derivatives much more quickly than with separate calls, although this isn't implemented yet. One use for this is to produce high-degree Taylor series. This could potentially be used in a future version of mpmath, for example to provide very fast moderate-precision (up to 50 digits, say) function for multi-evaluation of the Riemann zeta function on the real line or not too far away from the real line. This in turn would speed up other functions, such as the prime zeta function and perhaps polylogarithms in some cases, which are computed using series over many Riemann zeta values.

Dirichlet L-series

With the Hurwitz zeta function in place, it was a piece of cake to also supply an evaluation routine for arbitrary Dirichlet L-series.

The way it works it that you provide the character explicitly as a list (so it also works for other periodic sequences than Dirichlet characters), and it evaluates it as a sum of Hurwitz zeta values.

For example (copypasting from the docstring), the following defines and verifies some values of the Dirichlet beta function, which is defined by a Dirichlet character modulo 4:


>>> B = lambda s, d=0: dirichlet(s, [0, 1, 0, -1], d)
>>> B(0); 1./2
0.5
0.5
>>> B(1); pi/4
0.7853981633974483096156609
0.7853981633974483096156609
>>> B(2); +catalan
0.9159655941772190150546035
0.9159655941772190150546035
>>> B(2,1); diff(B, 2)
0.08158073611659279510291217
0.08158073611659279510291217
>>> B(-1,1); 2*catalan/pi
0.5831218080616375602767689
0.5831218080616375602767689
>>> B(0,1); log(gamma(0.25)**2/(2*pi*sqrt(2)))
0.3915943927068367764719453
0.3915943927068367764719454
>>> B(1,1); 0.25*pi*(euler+2*ln2+3*ln(pi)-4*ln(gamma(0.25)))
0.1929013167969124293631898
0.1929013167969124293631898

Appell hypergeometric function

The function appellf1 computes the Appell F1 function which is a hypergeometric series in two variables. I made the implementation reasonably fast by rewriting it as a single series in hyp2f1 (with the side effect of supporting the analytic continuation with that for 2F1).

A useful feature of the Appell F1 function is that it provides a closed form for some integrals, including those of the form

with arbitrary parameter values, for arbitrary endpoints (and that's a quite general integral indeed). Comparing with numerical quadrature:


>>> def integral(a,b,p,q,r,x1,x2):
...     a,b,p,q,r,x1,x2 = map(mpmathify, [a,b,p,q,r,x1,x2])
...     f = lambda x: x**r * (x+a)**p * (x+b)**q
...     def F(x):
...         v = x**(r+1)/(r+1) * (a+x)**p * (b+x)**q
...         v *= (1+x/a)**(-p)
...         v *= (1+x/b)**(-q)
...         v *= appellf1(r+1,-p,-q,2+r,-x/a,-x/b)
...         return v
...     print "Num. quad:", quad(f, [x1,x2])
...     print "Appell F1:", F(x2)-F(x1)
...
>>> integral('1/5','4/3','-2','3','1/2',0,1)
Num. quad: 9.073335358785776206576981
Appell F1: 9.073335358785776206576981
>>> integral('3/2','4/3','-2','3','1/2',0,1)
Num. quad: 1.092829171999626454344678
Appell F1: 1.092829171999626454344678
>>> integral('3/2','4/3','-2','3','1/2',12,25)
Num. quad: 1106.323225040235116498927
Appell F1: 1106.323225040235116498927

Also incomplete elliptic integrals are covered, so you can for example define one like this:


>>> def E(z, m):
...     if (pi/2).ae(z):
...         return ellipe(m)
...     return 2*round(re(z)/pi)*ellipe(m) + mpf(-1)**round(re(z)/pi)*\
...         sin(z)*appellf1(0.5,0.5,-0.5,1.5,sin(z)**2,m*sin(z)**2)
...
>>> z, m = 1, 0.5
>>> E(z,m); quad(lambda t: sqrt(1-m*sin(t)**2), [0,pi/4,3*pi/4,z])
0.9273298836244400669659042
0.9273298836244400669659042
>>> z, m = 3, 2
>>> E(z,m); quad(lambda t: sqrt(1-m*sin(t)**2), [0,pi/4,3*pi/4,z])
(1.057495752337234229715836 + 1.198140234735592207439922j)
(1.057495752337234229715836 + 1.198140234735592207439922j)

The Appell series isn't really faster than numerical quadrature, but it's possibly more robust (the singular points need to be given to quad as above to obtain any accuracy, for example). Mpmath doesn't yet have incomplete elliptic integrals; the version above could be a start, although I'll probably want to try a more canonical approach.

Another Mathematica bug

2009-07-12T20:32:00.004+02:00

Mathematica is great for cross-checking numerical values, but it's not unusual to run into bugs, so triple checking is a good habit.

Here is mpmath's evaluation for two Struve function values:


>>> struvel(1+j, 100j)
(0.1745249349140313153158106 + 0.08029354364282519308755724j)
>>> struvel(1+j, 700j)
(-0.1721150049480079451246076 + 0.1240770953126831093464055j)

The same values in Mathematica:


In[52]:= N[StruveL[1+I, 100I], 25]
Out[52]= 0.1745249349140313153158107 + 0.0802935436428251930875572 I

In[53]:= N[StruveL[1+I, 700I], 25]
Out[53]= -0.2056171312291138282112197 + 0.0509264284065420772723951 I

I'm almost certain that the second value returned by Mathematica is wrong. The value from mpmath agrees with a high-precision direct summation of the series defining the Struve L function, and even Mathematica gives the expected value if one rewrites the L function in terms of the H function:


In[59]:= n=1+I; z=700I
Out[59]= 700 I

In[60]:= N[-I Exp[-n Pi I/2] StruveH[n, I z], 25]
Out[60]= -0.1721150049480079451246076 + 0.1240770953126831093464055 I

Maple also agrees with mpmath:


> evalf(StruveL(1+I, 700*I), 25);
        -0.1721150049480079451246076 + 0.1240770953126831093464055 I

So unless Mathematica uses some nonstandard definition of Struve functions, unannounced, this very much looks like a bug in their implementation.

Wolfram Alpha reproduces the faulty value, so this still appears to be broken in Mathematica 7.

Improved incomplete gamma and exponential integrals; Clausen functions

2009-07-10T21:18:00.015+02:00

The SVN trunk of mpmath now contains much improved implementations of the incomplete gamma function (gammainc()) as well as the exponential integrals (ei(), e1(), expint()). Although the code is not quite perfect yet, this was a rather tedious undertaking, so I'm probably going to work on something entirely different for a while and give these functions another iteration later.

The incomplete gamma function comes in three flavors: the lower incomplete gamma function, the upper incomplete gamma function, and the generalized (two-endpoints) incomplete gamma function. The generalized incomplete gamma function is defined as

which reduces to the lower function when a = 0, and the upper version when b = +∞. A huge number of integrals occurring in pure and applied mathematics have this form (even Gaussian integrals, with a change of variables) so a solid incomplete gamma is quite important. It's especially important to ensure both speed in correctness in asymptotic cases.

The lower incomplete gamma function is easiest to implement, because it's essentially just a rescaled version of the hypergeometric 1F1 function and the 1F1 implementation already works well. Not much change was made here, so I'm going to write some about the other cases instead.

The exponential integrals are essentially the same as the upper incomplete gamma function and mostly share code. Remarks about the upper gamma function therefore also apply to the exponential integrals.

Upper gamma performance

The upper incomplete gamma function is hard because the two main series representations are an asymptotic series involving 2F0 that doesn't always converge, and an 1F1 series that suffers badly from cancellation for even moderately large arguments. The problem is to decide when to use which series.

The 2F0 series is very fast when it converges, while the 1F1 series is quite slow (due to the need for extra precision) just below the point where 2F0 starts to converge. After some experimentation, I decided to change the implementation of 2F0. Instead of performing a heuristic, conservative test to determine whether the series will converge (sometimes claiming falsely that it won't), it now always goes ahead to sum the series and raises an error only when it actually doesn't converge.

Thus the asymptotic series will always be used when possible, and although this leads to a slight slowdown for smaller arguments, it avoids worst-case slowness. The most important use for the incomplete gamma function, I believe, is in asymptotics, so I think this is a correct priority.

As a result, you can now do this:


>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> gammainc(10, 100)
4.083660630910611272288592e-26
>>> gammainc(10, 10000000000000000)
5.290402449901174752972486e-4342944819032375
>>> gammainc(3+4j, 1000000+1000000j)
(-1.257913707524362408877881e-434284 + 2.556691003883483531962095e-434284j)

The following graph compares old and new performance. The y axis shows the reciprocal time (higher is better) for computing gammainc(3.5, x) as x ranges between 0 and 100, at the standard precision of 15 digits. Red is the old implementation and blue is the new. The code also works with complex numbers, of course; replacing x with j*x gives a virtually identical graph (slightly scaled down due to the general overhead of complex arithmetic).

It's very visible where the asymptotic series kicks in, and the speed from then on is about 2000 evaluations per second which is relatively good. The new implementation is regrettably up to 3x slower than the old one for smaller x, although the slowdown is a bit misleading since the old version was broken and gave inaccurate results. The big dip in the blue graph at x = 10 is due to the automatic cancellation correction which the old code didn't use.

The gap between the asymptotic and non-asymptotic cases could be closed by using specialized series code for the lower incomplete gamma function, or using the Legendre continued fraction for intermediate cases (this comes with some problems however, such as accurately estimating the rate of convergence, and the higher overhead for evaluating a continued fraction than a series). This will certainly be worth doing, but I'm not going to pursue those optimizations right now for reasons already stated.

Some good news is that the graph above shows worst-case behavior, where the generic code is used, due the parameter 3.5. I've also implemented fast special-purpose code for the case when the first parameter is a (reasonably small) integer. This also means that the exponential integrals E1(x), Ei(x) as well as En(x) for integer n can be evaluated efficiently.

Here is a speed comparison between the old and new implementations of the ei(x) function, again at standard precision. There is actually no change in algorithm here: the old implementation used a Taylor series for small arguments and an asymptotic series for large arguments. The difference is due to using only low-level code; this turned out to buy a factor 2 in the Taylor case and more than an order of magnitude (!) in the asymptotic case.

The results are similar for the E1 function and with a complex argument. It is similar (only a bit slower) for gammainc(n,x) and expint(n,x) with a small integer value for n, although so far fast code is only implemented for real x in those cases.

Accurate generalized incomplete gamma function

The generalized incomplete gamma function can be written either as the difference of two upper gammas, or as the difference of two lower gammas. Which representation is better depends on the arguments. In general, one will work while the other will lead to total cancellation. gammainc is now clever enough to switch representations.

This uses a difference of lower gamma functions behind the scenes:


>>> gammainc(10000000, 3) - gammainc(10000000, 2)   # Bad
0.0
>>> gammainc(10000000, 2, 3)   # Good
1.755146243738946045873491e+4771204

This uses a difference of upper gamma functions behind the scenes:


>>> gammainc(2, 0, 100000001) - gammainc(2, 0, 100000000)   # Bad
0.0
>>> gammainc(2, 100000000, 100000001)   # Good
4.078258353474186729184421e-43429441

Some demo plots

Here are two plots of the upper gamma functions and exponential integrals (for various values of the first parameter). A lot of time went into getting the correct branch cuts in the low-level code (and writing tests for them), so please appreciate the view of the imaginary parts.


T1 = lambda x: gammainc(-2,x)
T2 = lambda x: gammainc(-1,x)
T3 = lambda x: gammainc(0,x)
T4 = lambda x: gammainc(1,x)
T5 = lambda x: gammainc(2,x)
plot([T1,T2,T3,T4,T5],[-5,5],[-10,10])


T1 = lambda x: expint(-2,x)
T2 = lambda x: expint(-1,x)
T3 = lambda x: expint(0,x)
T4 = lambda x: expint(1,x)
T5 = lambda x: expint(2,x)
plot([T1,T2,T3,T4,T5],[-5,5],[-10,10])

And a complex plot of gammainc(3+4j, 1/z):

A plot of gammainc(1/z, -1/z, 1/z); a rather nonsensical function (but that is besides the point):

Clausen functions

Unrelated to the gamma functions, I've also implemented Clausen functions:

These functions are just polylogarithms in disguise, but convenient as standalone functions. With them one can evaluate certain divergent Fourier series for example:


>>> clsin(-2, 3)
-0.01781786725924798006896962
>>> nsum(lambda k: k**2 * sin(3*k), [1,inf])
-0.01781786725924798006896962

They also work for complex arguments (and are related to zeta functions):


>>> clsin(2+3j, 1+2j)
(1.352229437254898401329125 + 1.401881614736751520048876j)
>>> clcos(2+3j, pi)
(-1.042010539574581174661637 - 0.2070574989958949174656102j)
>>> altzeta(2+3j)
(1.042010539574581174661637 + 0.2070574989958949174656102j)
>>> chop(clcos(zetazero(2), pi/2))
0.0

Meijer G, more hypergeometric functions, fractional differentiation

2009-06-29T23:04:00.011+02:00

My last update and the post before it detailed substantial improvements to the code for hypergeometric functions in mpmath, specifically the support for asymptotic expansions for 0F1, 1F1, 2F1, 2F0, plus the ability to evaluate hypergeometric-type formulas with singular parameters.

Over the past week-and-a-half I've done more work along the same lines. Importantly, I've implemented asymptotic expansions also for 1F2, 2F2 and 2F3 (commits 1, 2), so all hypergeometric functions of degree up to (2,3) now support fast evaluation for |z| → ∞ (1F0 also works, if anyone wonders -- it just happens to be a trivial case).

The next major remaining case is 3F2. It has a 1/z transformation, but this leaves |z| ≈ 1 which I don't know how to deal with. Does anyone who happens to be reading this know methods for evaluating 3F2 close to the unit circle? Taylor expansion around some point other than 0 works to some extent, but it's slow and in particular asymptotically slow close to z = 1, so not much help.

Bessel functions, etc.

As expected, the hypercomb function leads to very simple implementations of a large class of functions. I've now implemented the Whittaker, Struve, and Kelvin functions (they are called whitm, whitw, struveh, struvel, ber, bei, ker, kei). I've yet to update the orthogonal polynomials, but that shouldn't be much work. With this, I will have covered most of Bessel-type and hypergeometric-type functions listed on the Wolfram Functions site.

Speaking of Bessel functions, I also addressed most of the problems with their implementation in this commit. In particular, they can now be evaluated for huge arguments:


>>> mp.dps = 30
>>> print besselj(1, 10**20)
-7.95068198242545016504555020084e-11
>>> print chop(besselj(1, 10**20 * j))
(0.0 + 5.17370851996688078482437428203e+43429448190325182754j)
>>> print bessely(1,10**20)
-6.69800904070342428527377044712e-12
>>> print besselk(1,10**20)
9.66424757155048856421325779143e-43429448190325182776

This wasn't trivial, mainly because although hypercomb generally works, it sometimes becomes impossibly slow when computing functions straight from the definition. Basically: the function of interest might decrease exponentially, but internally it is computed by adding two nearly identical terms that grow exponentially, so the working precision and computation time increases exponentially. It's therefore still necessary to switch between different representations in different parts of the complex plane, and figuring that out involves some work. Some cases probably remain to be fixed.

As a followup to last week, I'll attach plots of J₀(1/z³), Y₀(1/z³) and K₀(1/z³):

The K plot unfortunately took very long time to finish -- almost an hour for 200,000 complex evaluations (J and Y were both much faster), so I'll probably have to optimize besselk a bit further.

Fractional derivatives

Another useful enhancement to the Bessel functions is that they can now be differentiated and integrated directly:


>>> mp.dps = 30
>>> print besselj(0, 3.5, derivative=2)
0.419378462090785430440501275011
>>> print diff(lambda x: besselj(0,x), 3.5, 2)
0.419378462090785430440501275011
>>> print besselj(0,3.5,derivative=-1) - besselj(0,2.5,derivative=-1)
-0.244675206320579138611991019242
>>> print quad(lambda x: besselj(0,x), [2.5, 3.5])
-0.244675206320579138611991019242

In fact, the representation for the derivatives works not just for integer orders (including negative values -- giving iterated integrals), but also for fractional or complex values. This led me to implement a general function differint for computing fractional order derivatives / integrals of arbitrary functions. (Commit.) It works essentially directly with the definition of the Riemann-Liouville differintegral.

It gives the same result as the special-purpose implementation for the Bessel function:


>>> print besselj(0, 3.5, derivative=0.5)
-0.436609427860836504473775239357
>>> print differint(lambda x: besselj(0,x), 3.5, 0.5)
-0.436609427860836504473775239357

One neat application is iterated integration. The following gives a 5-fold integral of f(x) = exp(π x), along with the symbolic evaluation of the same as a check:


>>> print differint(lambda x: exp(pi*x), 3.5, -5, x0=-inf)
194.790546022218468869246881408
>>> print exp(pi*3.5) / pi**5
194.790546022218468869246881408

Does anyone have any other interesting applications for fractional differentiation? I'd be interested in more examples to test with and possibly add to the documentation.

The Meijer G-function

At last, the Meijer G-function is implemented! (Commit.) Personally, I think this is something of a milestone for the mpmath project.

The Meijer G-function is very important because of its role in symbolic definite integration. Basically, definite integrals of Meijer G-functions (and even products of Meijer G-functions) just yield new Meijer G-functions; Mathematica and Maple therefore do many integrals by rewriting the input in terms of Meijer G-functions, applying Meijer G transformations, and converting the result back to simpler functions if possible.

Having the Meijer G-function in mpmath should be useful for anyone who wishes to implement a more powerful definite integrator in SymPy for example. It could also be useful for obtaining numerical values from integrals done by hand.

Looking around for examples to do stress testing with, I found a web page by Viktor Toth: Maple and Meijer's G-function: a numerical instability and a cure. His problem is to accurately evaluate G(-; 0; -1/2,-1,-3/2; -; x) for large real values of x. With my Meijer G-function implementation, I get:


>>> mp.dps = 15
>>> print meijerg([[],[0]],[[-0.5,-1,-1.5],[]],10)
4.89717497704114e-5
>>> print meijerg([[],[0]],[[-0.5,-1,-1.5],[]],100)
1.09696661341118e-12
>>> print meijerg([[],[0]],[[-0.5,-1,-1.5],[]],1000)
0.0
>>> print meijerg([[],[0]],[[-0.5,-1,-1.5],[]],10000)
1.53249554086589e+54

The third value should probably be small but not quite zero, and the last value is clearly bogus. Without looking at the details, the cause is almost certainly catastrophic cancellation of two huge terms. Fortunately, there is a cure for this:


>>> print meijerg([[],[0]],[[-0.5,-1,-1.5],[]],1000, check_cancellation=True)
3.34093555343418e-33
>>> print meijerg([[],[0]],[[-0.5,-1,-1.5],[]],10000, check_cancellation=True)
2.43925769071996e-94

The cancellation check should probably be enabled by default, either to automatically redo the computation as above or at least to issue a warning. The only catch with this is that it might lead to unnecessary slowdown and/or annoyance for the user in some cases, so I'll have to investigate how common those cases are.

Incidentally, I also tried plugging the above two calculations into Mathematica 6.0. It takes a long time to finish (mpmath gives the result instantaneously) -- and returns values that are wrong!


In[1]:= u1 = MeijerG[{{},{0}},{{-1/2,-1,-3/2},{}},1000]

                           3         1
Out[1]= MeijerG[{{}, {0}}, {{-(-), -1, -(-)}, {}}, 1000]
                           2         2

In[2]:= u2 = MeijerG[{{},{0}},{{-1/2,-1,-3/2},{}},10000]

                           3         1
Out[2]= MeijerG[{{}, {0}}, {{-(-), -1, -(-)}, {}}, 10000]
                           2         2

In[3]:= Timing[N[u1,20]]

Out[3]= {8.90265, 0.0017597930166135139087}

In[4]:= Timing[N[u1,50]]

                                                                    -33
Out[4]= {12.8231, 3.3409355534341801158987353523397047765918571151576 10   }

In[5]:=  Timing[N[u2,20]]

                                      50
Out[5]= {59.017, -2.0782671663885270791 10  }

In[6]:= Timing[N[u2,50]]

                                                                     22
Out[6]= {83.3753, -2.8700325450226332558088281915945389986057044454640 10  }

In[7]:= Timing[N[u2,120]]

Out[7]= {451.365, 2.439257690719956395903324691434088756714300374716395499173\

                                                                    -94
>      70196218529840153673260714339051464703903148052541923961351654 10   }

That's several minutes for something mpmath did in less than a second, and with less intervention. So maybe Maple and Mathematica should copy my Meijer G-function code ;-)

Complex roots

A small and trivial, but quite convenient, new feature: the nthroot function can now compute any of the roots of a given number and not just the principal root. I also added root as an alias since I have lazy fingers. Like so:


>>> for k in range(5):
...     r = root(-12, 5, k)
...     print chop(r**5), r
...
-12.0 (1.32982316461435 + 0.966173083818997j)
-12.0 (-0.507947249855734 + 1.56330088863444j)
-12.0 -1.64375182951723
-12.0 (-0.507947249855734 - 1.56330088863444j)
-12.0 (1.32982316461435 - 0.966173083818997j)

While I was at it, I couldn't resist also implementing a function unitroots for computing all the nth roots of unity, optionally just the primitive roots, as well as a function cyclotomic for evaluating the nth cyclotomic polynomial. As it turns out, cyclotomic polynomials are not entirely trivial to evaluate both efficiently and in a numerically stable way -- see the source code for my solution. Commit here. A side effect is that I also implemented a powm1 function that accurately gives x^y - 1 (possibly a useful complement to expm1 for other uses as well):


>>> print power(2,1e-100)-1
0.0
>>> print powm1(2, 1e-100)
6.93147180559945e-101

That will be all for now.

Massive hypergeometric update

2009-06-19T14:55:00.014+02:00

[Update: as further proof that the asymptotic expansions are working, I've plotted 1/erfi(1/z³), Ai(1/z³) and Bi(1/z³) around 0:

End of update.]

Today I committed a large patch to mpmath that significantly improves the state of hypergeometric functions. It's the result of about a week of work (plus some earlier research).

Perhaps most importantly, I've implemented the asymptotic expansions for 0F1 and 1F1 (the expansions for 2F1 were discussed in the previous post). They should now work for arbitrarily large arguments, as such:


>>> from mpmath import *
>>> mp.dps = 25
>>> print hyp0f1(3,100)
786255.7044208151250793134
>>> print hyp0f1(3,100000000)
4.375446848722142947128962e+8675
>>> print hyp0f1(3,10**50)
4.263410645749930620781402e+8685889638065036553022515
>>> print hyp0f1(3,-10**50)
-2.231890729584050840600415e-63
>>> print hyp0f1(1000,1000+10**8*j)
(-1.101783528465991973738237e+4700 - 1.520418042892352360143472e+4700j)
>>> print hyp1f1(2,3,10**10)
2.15550121570157969883678e+4342944809
>>> print hyp1f1(2,3,-10**10)
2.0e-20
>>> print hyp1f1(2,3,10**10*j)
(-9.750120502003974585202174e-11 - 1.746239245451213207369885e-10j)

I also implemented 2F0 and U (Kummer's second function), mostly as a byproduct of the fact that 2F0 is needed for the asymptotic expansions of both 0F1 and 1F1. 2F0 is an interesting function: it is given by a divergent series (it converges only in special cases where it terminates after finitely many steps):

However, it can be assigned a finite value for all z by expressing it in terms of U. Hence, mpmath can now compute the regularized sum of 2F0(a,b,z) for any arguments, say these ones:


>>> print hyp2f0(5, -1.5, 4)
(0.0000005877300438912428637649737 + 89.51091139854661783977495j)
>>> print hyp2f0(5, -1.5, -4)
102.594435262256516621777

(This ought to be a novel feature; SciPy and GSL implement 2F0, but only special cases thereof, and Mathematica doesn't have a direct way to evaluate 2F0.)

It's the asymptotic case where z → 0, where a truncation of the 2F0 series can be used, that is used for the expansions at infinity of 0F1 and 1F1.

Back to 0F1 and 1F1, the asymptotic expansions of these functions are important because they permit many special functions to be evaluated efficiently for large arguments. So far I've fixed erf, erfc, erfi, airyai and airybi to take advantage of this fact (except for erf and erfc of a real variable, all these functions were previously slow even for a moderately large argument, say |z| > 100).

Examples that now work (well, some of them possibly theoretically worked before too, but probably required hours or ages of universe to finish):


>>> print erf(10000+10000j)
(1.000001659143196966967784 - 0.00003985971242709750831972313j)
>>> print erfi(1000000000)
2.526701758277367028229496e+434294481903251818
>>> print erfc(1000-5j)
(-1.27316023652348267063187e-434287 - 4.156805871732993710222905e-434288j)
>>> print airyai(10**10)
1.162235978298741779953693e-289529654602171
>>> print airybi(10**10)
1.369385787943539818688433e+289529654602165
>>> print airyai(-10**10)
0.0001736206448152818510510181
>>> print airybi(-10**10)
0.001775656141692932747610973
>>> print airyai(10**10 * (1+j))
(5.711508683721355528322567e-186339621747698 + 1.867245506962312577848166e-186339621747697j)
>>> print airybi(10**10 * (1+j))
(-6.559955931096196875845858e+186339621747689 - 6.822462726981357180929024e+186339621747690j)

An essential addition in this patch is the function hypercomb which evaluates a linear combination of hypergeometric series, with gamma function and power weights:

This is an extremely general function. Here is a partial list of functions that can be represented more or less directly by means of it:

Regularized hypergeometric series

The generalized incomplete gamma and beta functions (and their regularizations)

Bessel functions

Airy, Whittaker, Kelvin, Struve functions, etc

Error functions

Exponential, trigonometric and hyperbolic integrals

Legendre, Chebyshev, Jacobi, Laguerre, Gegenbauer polynomials

The Meijer G-function

That's most of Abramowitz & Stegun, and means that the remaining hypergeometric-type functions available in Mathematica or Maxima but absent in mpmath will be easy to implement in the near future. With these additions, mpmath will have the most comprehensive support for numerical hypergeometric-type functions of any open source software, and should be very close to Mathematica.

The most important virtue of hypercomb is not that it allows for more concise implementations of various hypergeometric-type functions, although that is a big advantage too. The main idea is that hypercomb can deal with singular subexpressions, and particularly with gamma function poles that cancel against singularities in the hypergeometric series. These cases are almost more common than the nonsingular cases in practice, and hypercomb saves the trouble of handling them in every separate function.

Thus, in principle, a numerically correct implementation of hypercomb leads to correct implementations of all the functions in the list above. It's not a silver bullet, of course. For example, if a particular but very common case of some common function triggers an expensive limit evaluation in hypercomb, then it's probably better to handle that case with special-purpose code. There are also most likely some bugs left in hypercomb, although by now I have tested a rather large set of examples; large enough to be confident that it works soundly.

To show an example of how it works, an implementation of the Bessel J function might look like this:


def besj(n,z):
 z = mpmathify(z)
 h = lambda n: [([z/2],[n],[],[n+1],[],[n+1],-(z/2)**2)]
 return hypercomb(h, [n])

>>> mp.dps = 30
>>> print besselj(3,10.5)
0.1632801643733625756462387
>>> print besselj(-3,10.5)
-0.1632801643733625756462387
>>> print besj(3,10.5)
0.1632801643733625756462387
>>> print besj(-3,10.5)
-0.1632801643733625756462387

It gives the same value as the current Bessel J implementation in mpmath. Note that it works even when n is a negative integer, whereas naively evaluating the equation defining J(n,z) hits a gamma function pole and a division by zero in the hypergeometric series.

Thanks to the support for asymptotic expansions, this implementation (unlike the current besselj will also work happily with large arguments):


>>> print besj(3,1e9)
0.00000521042254280399021290721662036

I'm soon going to fix all the Bessel functions in mpmath along these lines, as I already did with erf, airyai, etc.

Here is another, more involved example (quoting directly from the docstring for hypercomb). The following evaluates

with a=1, z=3. There is a zero factor, two gamma function poles, and the 1F1 function is singular; all singularities cancel out to give a finite value:


>>> from mpmath import *
>>> mp.dps = 15
>>> print hypercomb(lambda a: [([a-1],[1],[a-3],[a-4],[a],[a-1],3)], [1])
-180.769832308689
>>> print -9*exp(3)
-180.769832308689

With some tweaks, the same code perhaps with a few tweaks could be used to symbolically evaluate hypergeometric-type functions (in the sense of rewriting them as pure hypergeometric series, and then possibly evaluating those symbolically as a second step). Symbolic support for hypergeometric functions is a very interesting (and hard) problem, and extremely important for computer algebra, but unfortunately I don't have time to work on that at the moment (there is more than enough to do just on the numerical side).

Hypergeometric 2F1, incomplete beta, exponential integrals

2009-06-11T18:39:00.008+02:00

One of the classes of functions I'm currently looking to improve in mpmath is the hypergeometric functions; particularly 1F1 (equivalently the incomplete gamma function) and the Gauss hypergeometric function 2F1.

For example, the classical orthogonal polynomials (Legendre, Chebyshev, Jacobi) are instances of 2F1 with certain integer parameters, and 2F1 with noninteger parameters allows for generalization of these functions to noninteger orders. Other functions that can be reduced to 2F1 include elliptic integrals (though mpmath uses AGM for these). With a good implementation of 2F1, these functions can be implemented very straightforwardly without a lot of special-purpose code to handle all their corner cases.

Numerical evaluation of 2F1 is far from straightforward, and the hyp2f1 function in mpmath used to be quite fragile. The hypergeometric series only converges for |z| < 1, and rapidly only for |z| << 1. There is a transformation that replaces z with 1/z, but this leaves arguments close to the unit circle which must be handled using further transformations. As if things weren't complicated enough, the transformations involve gamma function factors that often become singular even when the value of 2F1 is actually finite, and obtaining the correct finite value involves appropriately cancelling the singularities against each other.

After about two days of work, I've patched the 2F1 function in mpmath to the point where it should finally work for all complex values of a, b, c, z (see commits here). I'm not going to bet money that there isn't some problematic case left unhandled, but I've done tests for many of the special cases now.

The following is a very simple example that previously triggered a division by zero but now works:


>>> print hyp2f1(3,-1,-1,0.5)
2.5

The following previously returned something like -inf + nan*j, due to incorrect handling of gamma function poles, but now works:


>>> print hyp2f1(1,1,4,3+4j)
(0.492343840009635 + 0.60513406166124j)
>>> print (717./1250-378j/625)-(6324./15625-4032j/15625)*log(-2-4j)  # Exact
(0.492343840009635 + 0.60513406166124j)

Evaluation close to the unit circle used to be completely broken, but should be fine now. A simple test is to integrate along the unit circle:


>>> mp.dps = 25
>>> a, b, c = 1.5, 2, -4.25
>>> print quad(lambda z: hyp2f1(a,b,c,exp(j*z)), [pi/2, 3*pi/2])
(14.97223917917104676241015 + 1.70735170126956043188265e-24j)

Mathematica gives the same value:


In[17]:= NIntegrate[Hypergeometric2F1[3/2,2,-17/4,Exp[I z]],
    {z, Pi/2, 3Pi/2}, WorkingPrecision->25]
                                                                   -26
Out[17]= 14.97223917917104676241014 - 3.514976640925973851950882 10    I

Finally, evaluation at the singular point z = 1 now works and knows whether the result is finite or infinite:


>>> print hyp2f1(1, 0.5, 3, 1)
1.333333333333333333333333
>>> print hyp2f1(1, 4.5, 3, 1)
+inf

As a consequence of these improvements, several mpmath functions (such as the orthogonal polynomials) should now work for almost all complex parameters as well.

The improvements to 2F1 also pave the way for some new functions. One of the many functions that can be reduced to 2F1 is the generalized incomplete beta function:

An implementation of this function (betainc(a,b,x1,x2)) is now available in mpmath trunk. I wrote the basics of this implementation a while back, but it was nearly useless without the recent upgrades to 2F1. Evaluating the incomplete beta function with various choices of parameters proved useful to identify and fix some corner cases in 2F1.

One important application of the incomplete beta integral is that, when regularized, it is the cumulative distribution function of the beta distribution. As a sanity check, the following code successfully reproduces the plot of several beta CDF:s on the Wikipedia page for the beta distribution (I even got the same colors!):


def B(a,b):
    return lambda t: betainc(a,b,0,t,regularized=True)

plot([B(1,3),B(0.5,0.5),B(5,1),B(2,2),B(2,5)], [0,1])

The betainc function is superior to manual numerical integration because of the numerically hairy singularities that occur at x = 0 and x = 1 for some choices of parameters. Thanks to having a good 2F1 implementation, betainc gives accurate results even in those cases.

The betainc function also provides an appropriate analytic continuation of the beta integral, internally via the analytic continuation of 2F1. Thus the beta integral can be evaluated outside of the standard interval [0,1]; for parameters where the integrand is singular at 0 or 1, this is in the sense of a contour that avoids the singularity.

It is interesting to observe how the integration introduces branch cuts; for example, in the following plot, you can see that 0 is a branch point when the first parameter is fractional and 1 is a branch point when the second parameter is fractional (when both are positive integers, the beta integral is just a polynomial, so it then behaves nicely):


# blue, red, green
plot([B(2.5,2), B(3,1.5), B(3,2)], [-0.5,1.5], [-0.5,1.5])

To check which integration path betainc "uses", we can compare with numerical integration. For example, to integrate from 0 to 1.5, we can choose a contour that passes through +i (in the upper half plane) or -i (in the lower half plane):


>>> mp.dps = 25
>>> print betainc(3, 1.5, 0, 1.5)
(0.152380952380952380952381 + 0.4023774302466306150757186j)
>>> print quad(lambda x: x**2*(1-x)**0.5, [0, j, 1.5])
(0.152380952380952380952381 - 0.4023774302466306150757186j)
>>> print quad(lambda x: x**2*(1-x)**0.5, [0, -j, 1.5])
(0.152380952380952380952381 + 0.4023774302466306150757186j)

The sign of the imaginary part shows that betainc gives the equivalent of a contour through the lower half plane. The convention turns out to agree with that used by Mathematica:


In[10]:= Beta[0, 1.5, 3, 1.5]
Out[10]= 0.152381 + 0.402377 I

I'll round things up by noting that I've also implemented the generalized exponential integral (the En-function) in mpmath as expint(n,z). A sample:


>>> print expint(2, 3.5)
0.005801893920899125522331056
>>> print quad(lambda t: exp(-3.5*t)/t**2, [1,inf])
0.005801893920899125522331056

The En-function is based on the incomplete gamma function, which is based on the hypergeometric series 1F1. These functions are still slow and/or inaccurate for certain arguments (in particular, for large ones), so they will require improvements along the lines of those for 2F1. Stay tuned for progress.

In other news, mpmath 0.12 should be in both SymPy and Sage soon. With this announcement I'm just looking for an excuse to tag this post with both 'sympy' and 'sage' so it will show up on both Planet SymPy and Planet Sage :-) Posts purely about mpmath development should be relevant to both audiences though, I hope.

Mpmath 0.12 released

2009-06-09T21:20:00.004+02:00

I'll quote the mailing list announcement:

Mpmath version 0.12 is now available from the website:
http://code.google.com/p/mpmath/

It can also be downloaded from the Python Package Index:
http://pypi.python.org/pypi/mpmath/0.12

Mpmath is a pure-Python library for arbitrary-precision floating-point arithmetic that implements an extensive set of mathematical functions. It can be used as a standalone library or via SymPy (http://code.google.com/p/sympy/).

Version 0.12 mostly contains bug fixes and speed improvements. New features include various special functions from analytic number theory, Newton's method as an option for root-finding, and more versatile printing of intervals. It is now also possible to create multiple working contexts each with its own precision. Finally, mpmath now recognizes being installed in Sage and will automatically wrap Sage's fast integer arithmetic if available.

For more details, see the changelog:
http://mpmath.googlecode.com/svn/tags/0.12/CHANGES

Bug reports and other comments are welcome at the issue tracker at
http://code.google.com/p/mpmath/issues/list or the mpmath mailing list:
http://groups.google.com/group/mpmath

I've previously blogged about some of the new features:

Computing generalized Bell numbers
Approximate prime counting
Fun with zeta functions

See those posts or the documentation for many examples.

Cython mpmath performance, part 2

2009-06-06T17:22:00.004+02:00

A followup to the previous post. Turns out I had left a temporary integer variable undeclared, which slowed things down unnecessarily. After adding a simple cdef, the code runs a bit faster: at 53 bits, an addition now takes 850 ns, a multiplication 920 ns; a complex addition takes 1.25 µs and a complex multiplication takes 2.15 µs.

I've also implemented division and square root. To continue the benchmarking started in the previous post:


Real division, x/y

       mpmath    cython    sage

53     23.7 µs   1.49 µs   854 ns
100    24.6 µs   1.84 µs   1.17 µs
333    27 µs     2.69 µs   2.1 µs
3333   70.3 µs   44 µs     42.8 µs

Real square root, x.sqrt()

53     21.1 µs   1.75 µs   1.48 µs
100    22.5 µs   2.23 µs   1.68 µs
333    32 µs     3.81 µs   3.11 µs
3333   65.2 µs   35.3 µs   33.1 µs

Division is a bit slow compared Sage, but it might be possible to improve. As with the other arithmetic operations, the speedup vs ordinary mpmath is still more than a factor 10 at low precision.

Further, here is a benchmark of the new gamma function algorithm (mentioned in my talk at SD15), which is an order of magnitude faster than the algorithms currently used by both mpmath and MPFR. I had already written a pure Python implementation, which was quite fast, but it's insanely fast in Cython -- at really low precision it's even faster than the exponential function (although this fact probably indicates that the exponential function could be streamlined further). It's still half as fast as exp at 333-bit precision:


Real gamma function, x.gamma()   (gamma(x) in mpmath)

       mpmath    cython    sage

53     332 µs    8.6 µs    195 µs
100    434 µs    12.9 µs   257 µs
333    1.63 ms   52 µs     975 µs
3333   90.9 ms   9.95 ms   347 ms

Yes, that is a whopping 20,000 hundred-digit gamma functions per second!

Making the Cython code fully usable as an mpmath backend might require another week or two of work. This estimate includes rewriting utility code in Cython, but not porting transcendental functions. (The main reason why this will take longer is that I don't just want to translate the existing code, but redo the implementations from scratch to make everything uniformly faster.)

I'm not working full time on Cython mpmath backend; in parallel, I'm working on the mpmath wrapper code for Sage (some of which is ready and has been posted to the Sage trac) and on high-level implementations of various special functions. The new functions I have written partial implementations for so far include the generalized exponential integral (E_n), incomplete beta function, Hurwitz zeta function with derivatives, Appell hypergeometric functions of two variables, and matrix sqrt/exp/log/cos/sin. If anyone feels that their favorite special function is missing from Sage and/or mpmath, I'm taking requests :-)

I've also done little bit of work fixing hypergeometric functions for the hard cases, but really haven't made much progress on this so far (this is a substantial project which will take some time).

By the way, there should be a new mpmath release any day now. I mostly need to review the documentation and do some minor cleanup.

Cython mpmath performance

2009-05-27T12:56:00.004+02:00

I'm presently working on a Cython-based backend for mpmath, with the immediate goal to make mpmath competitive as a component of Sage. The Cython code currently depends on utility functions in the Sage library but should eventually be possible to compile for standalone use or together with SymPy's future Cython backend (linking directly against GMP/MPIR).

So far I've implemented a real type (mpf replacement) and a complex type (mpc replacement); addition, multiplication, and the real exponential function; just enough to do some basic benchmarking. Below is a comparison between standard Python mpmath (running on top of sage.Integer), the new Cython-based types, and Sage's RealNumber / ComplexNumber types. I used the inputs x = sqrt(3)-1, y = sqrt(5)-1, and for the complex cases z = x+yi, w = y+xi. The precision ranges between 53 bits (IEEE double compatible) and 3333 bits (≈1000 decimal digits).


Addition, x+y

       mpmath    cython    sage

53     8.72 µs   979 ns    707 ns
100    8.81 µs   1.01 µs   733 ns
333    10.1 µs   1.14 µs   737 ns
3333   10.6 µs   1.59 µs   1.15 µs

Multiplication, x*y

53     9.04 µs   968 ns    728 ns
100    9.56 µs   1.14 µs   901 ns
333    11.3 µs   1.59 µs   1.2 µs
3333   35.7 µs   25.4 µs   17.8 µs

Real exponential function, exp(x)

53     56.2 µs   9.41 µs   16 µs
100    52.8 µs   11.6 µs   11.3 µs
333    122 µs    25.4 µs   26.9 µs
3333   1.38 ms   831 µs    1.14 ms

Complex addition, z+w

53    14.7 µs   1.45 µs   978 ns
100   14.5 µs   1.7 µs    1 µs
333   16.3 µs   1.85 µs   992 ns
3333  17.7 µs   3.2 µs    1.73 µs

Complex multiplication, z*w

53    35.8 µs   2.72 µs   1.68 µs
100   34.5 µs   2.76 µs   2.08 µs
333   43.2 µs   4.38 µs   3.61 µs
3333  139 µs    98.6 µs   70 µs

A few observations can be made. First off, basic arithmetic with number instances is an order of magnitude faster when fully Cythonized. This isn't surprising because something as simple as an addition of two mpmath.mpf instances has to go through some 50 lines of Python code to check types, check for special cases, and round the result (creating lots of temporary objects, etc).

Reimplemented in Cython, arithmetic comes well within half the speed of Sage's numbers. I believe MPFR (which Sage uses) does arithmetic faster because it works directly with the limb data and uses some tricks to speed up the rounding, whereas my implementation uses the mpz interface straightforwardly. Some difference also comes from the fact that MPFR uses machine-precision integers for exponents, whereas I'm using mpz_t to allow arbitrary-size exponents. In any case, the results are very good.

At very low precision, memory allocations account for probably half the time (for both my Cython-based types and Sage's RealNumber). Thus there is some room for improvement later on by implementing a freelist.

The best news (for special functions) is that my exponential function is as fast as MPFR's, so the same should be true for other functions when I get there. The new exp (which uses a slightly different algorithm from the one currently in mpmath) is actually slightly faster than MPFR, although it should be said that performance for transcendental functions depends heavily on tuning, the inputs, and possibly other factors (I have no idea why Sage's RealNumber.exp in my benchmark is much slower at 53 bits than at 100 bits, for example), so this is not a conclusive result.

The code itself is currently too messy and ad-hoc to share publicly, sorry.

Report from Sage Days 15

2009-05-22T06:21:00.007+02:00

Sage Days 15 is now over. I'm still in Seattle as I write this; I depart early tomorrow. It's been a great week in almost every way. I didn't take photos, but see William Stein's photos and video recordings.

Edit: David Joyner also has photos:
http://sage.math.washington.edu/home/wdj/sagedays/sd15/

People

My prejudices were confirmed insofar as that the Sage developers are an exceptionally friendly bunch. Everybody was helpful and had a good sense of humor. And of course, these people are all very good at what they work on. It's both humbling and inspiring to be surrounded by people who are far more knowledgeable than yourself about almost everything.

There were a couple of talks daily (except for the last two days), most of which were very good; not to mention all the interesting informal discussions. Personally, I thought it was particularly interesting to get Bill Hart's (MPIR and FLINT guy) perspective of his work; he also provided useful technical advice. William's presentation of the general direction of the Sage project and the Cython talk by Robert Bradshaw and Craig Citro were also quite informative.

I gave a short talk about mpmath on Monday (slides). I'm not sure if the talk itself was good, but several people seemed to be interested in my work. I've received questions and comments about mpmath all week, and the general attitute towards using it in Sage seems positive.

Cython

On the purely technical side, by far the coolest thing about SD15 was the chance to finally learn Cython. I've seen Cython in action before, but never actually used it myself. My impressions are exclusively positive. I have plenty of experience with Python and a decent amount of experience with C, and Cython truly combines the best (alternatively, removes the worst) of both worlds. The biggest surprise (and of course this is Cython's selling point) is how simple the interfacing between high level and low level code becomes, and the fact that it is all very robust.

It's exiciting to see that there are several active projects around that attempt to speed up Python. The nice thing about Cython is that it doesn't give you "half the speed of C" or "maybe nearly the speed of C, 3 years from now" -- it gives the real deal, -O3 C, and it works right now. I hope Cython will be part of the standard Python distribution within a year or two from now, and thereby make it even easier to build Cython extensions, especially on Windows.

On a tangential note, it's good to see that progress is being made on the native Windows port of Sage. (I dual boot between Windows and Ubuntu, and vmware Sage on Windows is plainly inconvenient.) The pragmatic attitude of various Sage developers towards issues such as support for proprietary platforms gives the impression of a project in good hands.

Coding

I offered right before SD15 to look at speeding up the prime counting function in Sage by implementing the asymptotically fast Meissel-Lehmer-Lagarias-Miller-Odlyzko method. Doing this from scratch turned out to be unnecessary when Victor Miller emailed me the C implementation of the algorithm which he wrote two decades ago! The code compiles, and is amazingly fast on my laptop; unfortunately, it has some issues (such as not working for small n). It would be possible to wrap it in Sage as-is, and just use a different algorithm for small n. But what's probably an even better solution is that Victor himself (he showed up at SD) is now working on reimplementing his algorithm in Sage.

Most of my productive coding time (I also spent quite a lot of time just poking around at the Sage codebase -- which I've never really looked at before) went into wrapping mpmath for Sage. The purpose of wrapping mpmath is mainly to improve the support for special functions in Sage, firstly by just making the mpmath function library available and secondly by further improving the algorithms in mpmath. (I am getting funded via William's research grant to work on the latter this summer.)

As a first step, I made mpmath able to use sage.Integer instead of Python long. This was mostly trivial, the only problem being a coercion-related bug in sage.Integer that had to be fixed (William helped me write the patch, my first for Sage!).

When all mpmath tests finally passed, they unfortunately turned out to run 40% times slower than pure Python and 2.6x (!) slower than gmpy mpmath. The slowdown was in fact expected since sage.Integer was known to be substantially slower than gmpy.mpz. I therefore helped Robert Bradshaw identify the performance bottlenecks (they were mostly due either to unnecessary coercions or unnecessary memory reallocation), and he promptly patched them up. A nice side effect is that many parts of Sage should be faster now, since much code depends on the Integer type.

After I also implemented a few exra helper functions in Cython, mpmath on top of Sage now runs just a few percent slower than gmpy mpmath. Much of the remaining difference is probably due to inefficiencies in my own Cython code.

Finally, I also wrote wrapper code in Cython that permits Sage to call mpmath functions with Sage objects as input and getting Sage RealNumber or ComplexNumber back, all with reasonably low conversion overhead. This should now make it trivial to fix the state of numerical special functions in Sage (wrapping an mpmath function is literally a one-liner). I'm just waiting for Sage 4.0 with the new symbolics code to come out before I submit a patch and start adding wrapped functions; I want to ensure that the symbolics and numerics integrate well.

The state of mpmath is that it is still a bit slow compared to some numerical functions in Sage, but it's not at all bad. The long term plan is to speed up mpmath further by providing a pure-Cython backend, which along with certain algorithmic improvements should make it roughly as fast as MPFR and Pari for cases where it's currently slower.

Travel

The flight over here was reasonably convenient considering the circa 20-hour duration with transfers accounted for (the repeated security checks were really the biggest hassle). Hopefully the flight back tomorrow will proceed as smoothly.

Seattle, at least the University District, is a fairly nice city (not least because of the trees everywhere). Except for the trip to Microsoft Research on Tuesday, all activities took place
at the University of Washington campus. The campus alone is huge, so there has been plenty to see. I've lived very conveniently at the College Inn, right next to both the campus and The Ave.

I did make one excursion to downtown Seattle. I arrived early last Friday and with nothing on the schedule that day, I had time to meet up with Huy Pham (fellow Doomer and online acquintance) and we went to see Star Trek in IMAX (at the Pacific Science Center), which was well worth the admission.

As a last remark, I'm very grateful to William and the other Sage Days 15 organizers for their good job and for inviting me (and paying for my trip!). I've had a great time, had the opportunity meet interesting people, learned many things, and now I'm sad it's over.

Computing (generalized) Bell numbers

2009-03-04T19:48:00.006+01:00

Over the last few days, I procrastinated a bit from more important tasks by implementing Bell numbers (and more generally, evaluation of Bell polynomials) in mpmath. The Bell polynomials are defined by B₀(x) = 1 and

In particular, B_n(1) is the nth Bell number, which has combinatorial significance as the number of partitions of a set with n elements.

Instead of using the recurrence formula above, I implemented an approximate formula (discussed below) to provide fast approximate evaluation. So it is for example possible to do:


>>> from mpmath import *
>>> mp.dps = 30
>>> print bell(100000)
1.04339424254293899845402468388e+364471

This algorithms turns out to work well for exact evaluation at integer arguments, by setting the precision large enough, e.g. like this:


>>> from mpmath import *
>>> from mpmath.functions import funcwrapper
>>>
>>> @funcwrapper
... def bellexact(n,x=1):
...     mp.prec = 20
...     size = int(log(bell(n,x),2))+10
...     mp.prec = size
...     return int(bell(n,x)+0.5)
...
>>> bellexact(100)
47585391276764833658790768841387207826363669686825611466616334637559114
622672724044217756306953557882560751L
>>>
>>> bellexact(50)
185724268771078270438257767181908917499221852770L
>>> bellexact(40,40)
5288533501514377257070614176831982123270530650242136500078910824523240L

I did a bit of benchmarking for computing large Bell numbers B_n, and this approach appears to be much faster than what other systems use. Below are my timings for SymPy, GAP, Mathematica, and mpmath. (Note: as I do not have a personal Mathematica license, I ran it remotely on a different system, which based on a comparison of integer multiplication speed is 10-20% faster than my own laptop.)


      n        SymPy      GAP   Mathematica    mpmath
    100      0.007 s      0 s       0.005 s   0.012 s
    1000       8.9 s   0.78 s        0.32 s    0.30 s
    3000           -     67 s          12 s     3.6 s
    10000          -        -         396 s      76 s

From what I've read, Bell numbers have some interesting number-theoretical properties such as Touchard's congruence, so this could potentially be useful to someone. (Interestingly, I wonder if Touchard's congruence could be turned into an even faster modular algorithm for Bell numbers, similar to David Harvey's algorithm for Bernoulli numbers.)

The formula I implemented for the Bell polynomials is the infinite series

which is known by the name Dobinski's formula.

As a side note, I implemented this series using a piece of experimental new summation code designed to ensure full accuracy in spite of cancellation. For example, the following evaluation involves alternating terms of magnitude over 10⁴⁰⁰; the internal precision is automatically set to over 400 digits to provide the requested 15 digits:


>>> mp.dps = 15
>>> print bell(10,-1000)
9.55744142784509e+29

I'll probably be working more this summer on implementing similar methods for several mpmath functions that are currently implemented somewhat less robustly (more on that some other time).

One interesting aspect of Dobinski's formula is that it generalizes to arbitrary complex values for n, not just integers. I always prefer to implement functions for the most general possible arguments (and will sometimes hesitate to implement a function if I can only do it for a special case). Generalization is extremely useful: extending functions defined on the integers to the reals permits the use of differentiation, continuous optimization algorithms, etc. Extending further to the complex numbers allows for contour integration. Curiously, Mathematica does not seem to do this for its BellB function.

Unfortunately, there is a difficulty in implementing Dobinski's formula for general n: the 0th term is singular unless n is a nonnegative real number. Removing the term completely breaks compatibility with the standard definition of the Bell polynomials for n = 0. Keeping it either makes the series undefined, or at best renders the function discontinuous at the single point n = 0 which is certainly not desirable. As a workaround, I changed the formula to

The sinc term is continuous (analytic), and zero when n is an integer except when n = 0, so it solves all problems. The sinc function has no connection to Bell polynomials whatsoever as far as I know, but it is the simplest possible "patch" I can think of.

Voila, we now get a smooth transition between for example B₀(x) = 0 and B₁(x) = x:


>>> f = lambda k: (lambda x: bell(k,x))
>>> plot([f(k) for k in linspace(0,1,5)], [-2,3])

It also becomes possible to e.g. compute a Taylor series of B seen as a function of n, using either a step sum or complex integration for the derivatives:


>>> nprint(chop(taylor(lambda n: bell(n,1), 0, 5, method='step')))
[1.0, 0.222119, -0.504269, 3.32909e-2, 0.307601, 2.09191e-3]
>>> nprint(chop(taylor(lambda n: bell(n,1), 0, 5, method='quad')))
[1.0, 0.222119, -0.504269, 3.32909e-2, 0.307601]

We can for example use the Taylor expansion around n = 0 to evaluate the Bell polynomial with n = 1 (not really practical for anything, but a beautiful example of the power of analysis):


>>> print polyval(taylor(lambda n: bell(n,3), 0, 5)[::-1], 1)
2.9949024096858
>>> print polyval(taylor(lambda n: bell(n,3), 0, 10)[::-1], 1)
2.99998868106305
>>> print polyval(taylor(lambda n: bell(n,3), 0, 15)[::-1], 1)
2.99999998751065
>>> print polyval(taylor(lambda n: bell(n,3), 0, 25)[::-1], 1)
3.0
>>> print bell(1,3)
3.0

Finally, I also implemented the raw series

as a separate function, because it seems like it could be useful in its own right. Unfortunately, this function does not appear to have a standardized notation. For certain special values, it reduces to a Bell polynomial times an exponential function; an incomplete gamma function; or a generalized hypergeometric function. It is also related to the Poisson distribution. Other than that, I have not found much information about it. Does this function have a name? I named it polyexp in mpmath, because it can be seen as an exponential analog of the polylogarithm.

As a function of s, it is an L-series, and shows interesting, highly chaotic behavior in the complex plane:


>>> mp.dps = 5
>>> cplot(lambda s: polyexp(s,1), [-20,20],[0,60], points=50000)
>>> plot(lambda s: polyexp(s*j,1), [0,60])

How (not) to compute harmonic numbers

2009-02-21T12:40:00.007+01:00

The nth harmonic number is the nth partial sum of the divergent harmonic series, H_n = 1 + 1/2 + 1/3 + ... + 1/n. The simplest way to compute this quantity is to add it directly the way it is written: 1, 1+1/2, 1+1/2+1/3, and so on. For n approximately greater than 10 or 100, this is algorithm is not a very good one. Why?

Firstly, in floating-point arithmetic, it is much better to use the asymptotic expansion

which requires fewer terms the larger n is. Combining this with recurrence when n is too small gives a fast numerical algorithm for fixed-precision or arbitrary-precision computation of harmonic numbers (especially large ones). It allows one to determine, for example, that it takes exactly 15092688622113788323693563264538101449859497 terms until the partial sums of the harmonic series exceed 100. (This algorithm is implemented in mpmath, and most other comparable software.)

Secondly, the naive algorithm is especially bad for exact computation of harmonic numbers. This might not be obvious at a glance. We are just adding n numbers: what's to improve? The catch is that rational arithmetic has hidden costs: a single addition requires three multiplications and usually a GCD reduction. In contrast with integer or floating-point addition, rational addition gets very slow when the numerators and denominators get large. Computing harmonic numbers the naive way triggers worst-case behavior of rational arithmetic: the denominators grow like n! (if no GCD is performed, the denominator of H_n literally becomes n!).

It should be well known that computing n! as 1, 1·2, 1·2·3, ... is a poor method when n is large. A much better algorithm, assuming that multiplication is subquadratic, is to recursively split the products in half, i.e. compute (1·2·...·n/2) · ((n/2+1)·...·n) and so on. This algorithm has complexity O(log(n) M(n log(n)) compared to O(n² log(n)) for the naive algorithm, where M(n) is the complexity of multiplying two n-bit integers. This approach is called binary splitting, and has many other applications.

The example of factorials directly suggests an analogous algorithm for computing harmonic numbers with reduced complexity: just recursively split the summation in half.

A few more variations are also possible. Instead of working with rational numbers, which require a GCD reduction after each addition, one can work with the unreduced numerators and denominators as integers, and form a rational number at the end. Another way to obtain a pure-integer algorithm for H_n = p/q is to compute the denominator q = n! and then use the formula

for the numerator.

Below, I have implemented five different algorithms for computing harmonic numbers based on the preceding remarks. The code works with either Python 2.6, using the Fraction type from the standard library, or with gmpy, by adding the following respective header code:


# python
from math import factorial as fac
from fractions import Fraction as mpq
one = 1
mpz = int

# gmpy
from gmpy import mpz, mpq, fac
one = mpz(1)

Algorithm 1: directly add 1, 1/2, 1/3, ..., 1/n:


def harmonic1(n):
    return sum(mpq(1,k) for k in xrange(1,n+1))

Algorithm 2: like algorithm 1, GCD reduction postponed to the end:


def harmonic2(n):
    p, q = mpz(1), mpz(1)
    for k in xrange(2,n+1):
        p, q = p*k+q, k*q
    return mpq(p,q)

Algorithm 3: compute denominator; compute numerator series:


def harmonic3(n):
    q = fac(n)
    p = sum(q//k for k in xrange(1,n+1))
    return mpq(p,q)

Algorithm 4: binary splitting:


def _harmonic4(a, b):
    if b-a == 1:
        return mpq(1,a)
    m = (a+b)//2
    return _harmonic4(a,m) + _harmonic4(m,b)

def harmonic4(n):
    return _harmonic4(1,n+1)

Algorithm 5: binary splitting, GCD reduction postponed to the end:


def _harmonic5(a, b):
    if b-a == 1:
        return one, mpz(a)
    m = (a+b)//2
    p, q = _harmonic5(a,m)
    r, s = _harmonic5(m,b)
    return p*s+q*r, q*s

def harmonic5(n):
    return mpq(*_harmonic5(1,n+1))

I did some benchmarking on my laptop (32-bit, 1.6 GHz Intel Celeron M), with n up to 10⁴ for the Python version and 10⁶ for the gmpy version. Here are the results (times in seconds):


python:
n         harmonic1  harmonic2  harmonic3  harmonic4  harmonic5
10         0.000159   0.000009   0.000010   0.000160   0.000025
100        0.004060   0.000233   0.000234   0.002211   0.000358
1000       0.889937   0.028822   0.029111   0.048106   0.026434
10000    769.109402   3.813702   3.823027   2.710409   3.475280

gmpy:
n         harmonic1  harmonic2  harmonic3  harmonic4  harmonic5
10         0.000028   0.000010   0.000010   0.000033   0.000021
100        0.000284   0.000097   0.000111   0.000365   0.000226
1000       0.003543   0.001870   0.004986   0.004280   0.002651
10000      0.103110   0.142100   0.588723   0.059249   0.052265
100000     7.861546  17.002986  74.712460   1.115333   1.200865
1000000  994.548226      -           -     27.927289  24.425421

Algorithm 1 clearly does not do well asymptotically. For the largest n measured, it is 284 times slower than the fastest algorithm with Python arithmetic and 40 times slower with gmpy arithmetic. Interestingly, algorithms 2 and 3 both substantially improve on algorithm 1 with Python arithmetic at least up to n = 10000, but both are worse than algorithm 1 with gmpy. This is presumably because GMP uses a much better GCD algorithm that reduces the relative cost of rational arithmetic.

Algorithms 4 and 5 are the clear winners for large n, but it is hard to tell which is better. The timings seem to fluctuate a bit when repeated, so small differences in the table above are unreliable. I think algorithm 5 could be optimized a bit if implemented in another language, by accumulating temporary values in machine integers.

It is amusing to note that computing the 10000th harmonic number with the recursive algorithm and an optimized implementation of rational arithmetic is 14,700 times faster than the direct algorithm with a non-optimized implementation of rational arithmetic. Choosing a larger n (say, 10⁵) gives an even more sensational ratio, of course, which I'm not patient enough to try. The moral is: Moore's law is not an excuse for doing it wrong.

Algorithm 1 is not all bad, of course. With memoization (or implemented as a generator), it is good enough for sequential generation of the numbers H₁, H₂, H₃, ..., and this tends to be needed much more commonly than isolated large harmonic numbers. SymPy therefore uses the memoized version of algorithm 1.

The algorithms above can be adapted for the task of computing generalized harmonic numbers, i.e. numbers of the form 1 + 1/2^k + 1/3^k + ... + 1/n^k for some integer k. Also worth noting is that algorithm 5 allows for computation of Stirling numbers of the first kind. Since |S(n+1,2)| = n! · H_n, algorithm 5 can be viewed as an efficient algorithm for simultaneous computation of S(n,2) and n!. Expressing Stirling numbers in terms of generalized harmonic numbers, this extends to an algorithm for S(n,k) for any (large) n and (small) k.

A final remark: binary splitting is not the fastest algorithm for computing factorials. The fastest known method is to decompose n! into a product of prime powers and perform balanced multiplication-exponentiation. Can this idea be applied to harmonic numbers and/or Stirling numbers? It's not obvious to me how it would be done. Are there any other, possibly even faster algorithms for harmonic numbers? I've looked around, but been unable to find anything on the subject.

Python 3.1 twice as fast as 3.0

2009-02-19T08:52:00.004+01:00

For large-integer arithmetic, that is. Mark Dickinson has beeen working on a patch that changes the internal representation of long integers from 15-bit to 30-bit digits, and there is an additional patch that adds a few algorithmic optimizations, written by Mark and Mario Pernici. Likely both patches will pass review and be included in Python 3.1. I asked Mark to run the mpmath unit tests to find out how big a difference the patches make, which he kindly did. Here are the results:


    32-bit build:
    unpatched (15-bit digits): 42.91 seconds
    patched (30-bit digits): 40.57 seconds
    patched+optimized: 30.15 seconds

    64-bit build:
    unpatched: 38.39 seconds
    patched: 22.36 seconds
    patched+optimized: 21.55 seconds
    patched+optimized+bitlength: 20.20 seconds

The last result includes the use of the new int.bit_length() method (which I had a small part in implementing) instead of the pure-Python version in mpmath.

It's not surprising that the patches make high-precision arithmetic much faster, but most of the mpmath tests actually work at low precision and depend mostly on general speed of the Python interpreter. There the speedup is perhaps of order 0-30%. For the tests working at very high precision, the improvement is a factor 4 or 5 with the 64-bit build. Then there are a number of tests in between. With some imagination, all unit tests together provide a plausible estimate of actual performance for a wide range of applications.

An excerpt of before/after results for some particular tests, comparing 64-bit unpatched and 64-bit patched+optimized+bitlength:


    # Only tests 53-bit precision
    double_compatibility      ok        1.3149240 s
    double_compatibility      ok        1.1979949 s

    # Logarithms up to 10,000 digits
    log_hp                    ok        4.0845711 s
    log_hp                    ok        0.7967579 s

    # Large Bernoulli numbers
    bernoulli                 ok        6.8625491 s
    bernoulli                 ok        1.4261260 s

    # Gamma function, high precision
    gamma_huge_2              ok        2.4907949 s
    gamma_huge_2              ok        0.5781031 s

    # Numerical quadrature, 15 to 100 digit precision
    basic_integrals           ok        3.0117619 s
    basic_integrals           ok        1.8687689 s

Mpmath does not yet run on Python 3.x; the benchmarking was made with a version of mpmath hacked slightly for the tests to pass. It would be nice to provide a 3.x compatible version of mpmath soon. The 2to3 tool fortunately handles almost all necessary patching; a handful of fairly simple additional changes need to be made. The most annoying problem is the lack of cmp (and particularly the cmp argument for sorted), which has no trivial workaround, but still should not be too hard to fix. In any case, it seems likely that the 30-bit patch will also be backported to Python 2.7, so most users should be able to take advantage of it.

It would be interesting to compare these benchmarks with the GMPY mode in mpmath. GMPY too provides a ~2x speedup for the unit tests. (This factor would be much larger if tests at even higher precision were included. At a million digits, GMPY is orders of magnitude faster than pure Python.) Originally GMPY only sped up the very high precision tests, and was otherwise slower than pure Python, probably due to Python's internal optimizations for int and long instances. This disadvantage was eliminated by implementing some helper functions for mpmath in C in GMPY. Mario Pernici has recently worked on further optimizations along the same lines, which should substantially improve low-precision performance.