Fast and accurate double precision implementation of incomplete gamma function

Question

What is the state of the art way of implementing double precision special functions? I need the following integral: $$ F_m(t) = \int_0^1 u^{2m} e^{-tu^2} d u = {\gamma(m+{1\over 2}, t)\over 2 t^{m+{1\over 2}}} $$ for $m=0, 1, 2, ...$ and $t>0$, which can be written in terms of the lower incomplete gamma function. Here is my Fortran and C implementation:

https://gist.github.com/3764427

which uses series expansion, sums up the terms until the given accuracy, and then uses recursion relations to efficiently obtain values for lower $m$. I tested it well and I get 1e-15 accuracy for all values of parameters that I need, see the Fortran version's comments for details.

Is there a better way to implement it? Here is a gamma function implementation in gfortran:

https://github.com/mirrors/gcc/blob/master/libgfortran/intrinsics/c99_functions.c#L1781

it is using rational function approximation instead of summing up some infinite series that I am doing. I think that's a better approach, because one should obtain uniform accuracy. Is there some canonical way to approach these things, or does one have to figure out a special algorithm for each special function?

Update 1:

Based on the comments, here is the implementation using SLATEC:

https://gist.github.com/3767621

it reproduces values from my own function, roughly on the level of 1e-15 accuracy. However, I noticed a problem that for t=1e-6 and m=50, the $t^{m+{1\over2}}$ term gets equal to 1e-303 and for higher "m" it simply starts giving incorrect answers. My function doesn't have this problem, because I use a series expansion/recurrence relations directly for $F_m$. Here is an example of a correct value:

$F_{100}$(1e-6)=4.97511945200351715E-003,

but I can't get this using SLATEC because the denominator blows up. As you can see, the actual value of $F_m$ is nice and small.

Update 2:

To avoid the above problem, one can use the function dgamit (Tricomi's incomplete Gamma function), then F(m, t) = dgamit(m+0.5_dp, t) * gamma(m+0.5_dp) / 2, so there is no problem with $t$ anymore, but unfortunately the gamma(m+0.5_dp) blows up for $m\approx 172$. This however might be high enough $m$ for my purposes.

Answer 1

The integral in question is also known as the Boys function, after the British chemist Samuel Francis Boys who introduced its use in the early 1950s. A few years ago, I needed to compute this function in double precision, as fast as possible but accurately. I managed to achieve a relative error on the order of $10^{-15}$ across the entire input domain.

It is generally advantageous to use different approximations for small and large arguments, where the optimal switch-over between "large" and "small" is best determined experimentally, and is in general a function of $m$. For my code, I defined "small" arguments as those satisfying condition $a \le m + 1{1\over 2}$.

For large arguments, I compute

$$\mathrm{F}_m(a) = {1\over 2}\gamma\left(m + {1\over 2}, a\right) \times p \times p, \space \space p = a^{-{1\over 2}\left(m+ {1\over 2}\right)}$$

This order of operations avoids premature underflow. As we need only the lower incomplete gamma function of half-integer orders here rather than a fully general lower incomplete gamma function, it is advantageous from a performance perspective to compute

$$\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right) - \Gamma\left(m + {1\over 2}, a\right)$$

using tabulated values of $\Gamma \left(m + {1\over 2}\right)$ and computing $\Gamma\left(m + {1\over 2}, a\right)$ according to this answer, while carefully avoiding the issue of subtractive cancellation through use of a fused multiply-add operation. A potential further optimization is to observe that for sufficiently large $a$, $\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right)$ to within a given floating-point precision.

For small arguments, I started with a series expansion for the lower incomplete gamma function from

A. Erdelyi, W. Magnus, F. Oberhettinger, and F.G. Tricomi, "Higher Transcendental Functions, Vol. 2". New York, NY: McGraw-Hill 1953

and modified it to compute the Boys function $\mathrm{F}_{m}(a)$ as follows (truncating the series when the term is sufficiently small for a given precision):

$$\mathrm{F}_{m}(a) = {1\over 2}\frac{1}{m + {1\over 2}}\exp(-a)\left(1+\sum_{n=1}^{\infty}\frac{a^{n}}{(1 + m + {1\over 2}) \times\space ...\space \times (n + m + {1\over 2})}\right)$$

There are also interesting and potentially important special cases for low orders of the Boys function, particularly $m = 0, 1, 2, 3$. First, we have $\mathrm{F}_{0}(a) = \sqrt{\frac{\pi}{4a}}\mathrm{erf}\left(\sqrt{a}\right)$, where $\mathrm{erf}$ is the error function provided in Fortran 2008 as the elemental function ERF and in C/C++ as the standard library functions erf and erff.

For fast computation when $m = 1, 2, 3$, I use custom minimax polynomial approximations for small arguments, say $a \lt {2{1\over 2}}$, and forward recursion $\mathrm{F}_{m}(a) = \frac{1}{2a}\left(\left(2m-1\right)\mathrm{F}_{m-1}(a) - \exp(-a)\right)$, for large ones, where issues with subtractive cancellation in the latter are mitigated by the use of fused multiply-add operations.

Where function values have to be computed for a given $a$ across multiple orders $m$, one would want to compute the function value for the highest value of $m$ directly, i.e. as discussed above, then use the numerically stable backward recursion $\mathrm{F}_{m-1} = \frac{1}{2m-1} \left(2a \space \mathrm{F}_{m}(a) + \exp\left(-a\right)\right)$ to compute all the other function values.

Answer 2

You might take a look at Numerical Methods for Special Functions by Amparo Gil, Javier Segura, and Nico M. Temme.

Answer 3

I'd take a look at Abramowicz & Stegun's book, or the newer revision that NIST has published a couple of years ago and that's available online I believe. They also discuss ways to implement things in a stable way.

Answer 4

It doesn't seem to be state-of-the-art but SLATEC at Netlib offers "1400 general purpose mathematical and statistical routines." The incomplete Gamma is available under the special functions here.

Implementing such functions is time consuming and error prone so I wouldn't do it myself unless absolutely necessary. SLATEC has been around for quite a while now and is widely used, at least based on the download counts, so I would expect the implementation to be mature.