# What are the efficient, accurate algorithms for evaluation of hypergeometric functions?

I'm curious to know what good numerical algorithms exist for evaluation of the generalized hypergeometric function (or series), defined as

$${}_pF_q(a_1,\ldots,a_p;b_1,\ldots,b_q;z) = \sum_{k=0}^{\infty} \frac{(a_1)_k\cdots(a_p)_k}{(b_1)_k\cdots(b_q)_k}\frac{z^k}{k!}$$

In general, this series is not necessarily going to converge very fast (or at all), so summing up terms one by one seems less than ideal. Is there some alternate method that works better? To be specific, I'm looking for something that will give 4 or 5 digits of precision with a reasonable number of calculations.

The most common cases that I usually see used are $p=1,q=1$ and $p=2,q=1$, but in the particular project I'm working on, I have a need for $p=1,q=2$. Obviously a general algorithm for any $p$ and $q$ is ideal, but I'll take what I can get.

-
If your case it is not covered in Abramowitz and Stegun's Handbook (people.math.sfu.ca/~cbm/aands/subj.htm), which it is not, you are basically doomed to figure it out on your own, I'm afraid... –  Jaime Dec 10 '12 at 18:52

In a single application, it is rather likely that you only will need a small subset of all the possible extremes of the generalized hypergeometric function. It is a very general function, after all. Having an idea about the range of $z$ and of the parameters $a_i, b_i$ would allow giving more specific advice.

In general, the standard method, assuming $p \le q + 1$, is of course to use the defining power series when $|z|$ is small. If $p < q + 1$, it is best to switch to an asymptotic expansion when $|z|$ is large, either because the Taylor series converges too slowly and/or because it becomes too inaccurate due to catastrophic cancellation. The best cutoff between these algorithms depends on the parameters and the accuracy requirements.

For ${}_1F_2$ the asymptotic series is given by http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F2/06/02/03/ It looks rather horrible, but if your $a_1, b_1, b_2$ are fixed, you can compute numerical values for the coefficients in advance. General formulas are found in the DLMF: http://dlmf.nist.gov/16.11 (Note that some care is required to select the correct branch cuts.)

If there is a range where neither the Taylor series nor the asymptotic series works well enough, "exponentially-improved expansions" might be useful. Another possibility worth mentioning is that you could just plug the hypergeometric differential equation into a general-purpose ODE solver. This should work quite well especially if you only need 4-5 digits. This can be used to do analytic continuation from a small $z$ (where the power series works well) to a larger one, or in reverse from a value obtained via an asymptotic series (you might need to do a bit more work to get all the derivatives needed as initial values).

If you need functions with $p = q + 1$ on the whole complex plane then $1/z$ transformation formulas can be used to map the exterior of the unit disk to the interior. Some convergence acceleration algorithm or other method, like numerical integration of the ODE, must be used close to the unit circle. If $p > q + 1$ the radius is convergence is zero, so if the function you want to evaluate is given by such a divergent series, you might need to apply a Borel transform (numerically or symbolically) to reduce it to a convergent series.

For a complete implementation, there are other issues to consider as well (for example, dealing with parameters which are extremely large or very close to negative integers). For sufficiently bad parameters, it will be very difficult to obtain accurate values with double precision, no matter what you do, so arbitrary-precision arithmetic might be needed.

I should note that I have written a nearly complete numerical implementation of the generalized hypergeometric function for the mpmath library (it is currently missing asymptotic series for functions higher than ${}_2F_3$), which might be useful to study or run tests against (assuming it is not fast enough already for your purposes).

-
Excellent! Unfortunately I can't really get more specific about the parameter values because the function pops up in many places with various values. I will definitely be interested in using and/or looking at your implementation in mpmath at some point. –  David Z Dec 10 '12 at 21:30
Fredrik's answer is correct. I would only point out, that I ended up using a rational approximation (from Mathematica) for special values of the "a" and "b" coefficients, because it is accurate for all real "z" (I split the real axis into intervals and used a different rational approximation on each) and very fast. I used mpmath to check the accuracy of my double precision implementation in Fortran. –  Ondřej Čertík Dec 11 '12 at 20:45
The chapter in A&S about the hypergeometric function is actually about $_2F_1$, not the $_1F_2$ the OP is asking about... –  Jaime Dec 11 '12 at 1:01