bug in specfunc/gamma_inc.c

bug in specfunc/gamma_inc.c

[Jason Hooper Stover]

I think there is a bug in gamma_inc.c:gamma_inc_Q_CF().

gsl_sf_gamma_inc_Q (185.0, 200.0) gives the value 0.13594954199834325,
whereas a gp script gives 0.13594954213327904. The disparity
between gsl and gp grows as the second argument increases. (I guess
gp could be wrong, but since it has almost-arbitrary precision, I figured
it's correct.)

The function gamma_inc_Q_CF(a,x) was the subject of a previous bug report.
This function replaced an older function that was said to be unstable
for some regions in the plane.

The problem arises within the region 0.3*a<x<2.3*a, a>10.
This is the region for which N.M. Temme proposed another algorithm in 

"On the computation of the incomplete gamma functions for large values
of the parameters", Algorithms for Approximation, J.C. Mason and M.G.
Cox, eds. The Institute for Mathematics and Its Applications Conference
Series no. 10. 1987.

This algorithm might give a better approximation in this region than the
function gamma_inc_Q_CF().

-Jason

Re: bug in specfunc/gamma_inc.c

[Brian Gough]

Jason Hooper Stover writes:
 > gsl_sf_gamma_inc_Q (185.0, 200.0) gives the value 0.13594954199834325,
 > whereas a gp script gives 0.13594954213327904. The disparity
 > between gsl and gp grows as the second argument increases. (I guess
 > gp could be wrong, 

Hmmm... GSL-1.4 produces the correct answer to full double-precision
accuracy for gsl_sf_gamma_inc_Q(185,200) (for integer "a" this can be
verified analytically).

The gp front-end prints too many digits by default:

$ gp
? incgam(185,200)/gamma(185)
%1 = 0.1359495421332790409768941685
? default(realprecision,300)
? incgam(185,200)/gamma(185)
%2 = 0.13594954199834326027261113221577.....

 > but since it has almost-arbitrary precision, I figured
 > it's correct.)

I can't believe you said that ;-)

-- 
Brian Gough