Analytic continuation of 3F2, 4F3 and higher functions

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.