Bug: gamma_inc

[Hans Ekkehard Plesser <hans dot plesser at itf dot nlh dot no>]

Hi!

The incomplete Gamma function implemented in GSL 1.0 fails for certain
arguments, e.g.

	   gsl_sf_gamma_inc_P(a==34, x==32)

The failure occurs in routine gamma_inc_Q_CF in gamma_inc.c:

For any x == a-2, the denominator in rhok (l. 206) vanishes for k==2.
For x==32, 64, this makes rhok=inf, so that tk=inf and the iterations
runs through 5000 times operating on NaNs, before throwing a
GSL_MAXITER error.  For x not a power of two, one is saved by
numerical inaccuracy, but 

 run_gamma_inc(37, 35-20*eps)

yields

k      tk                sum

1   -3.60000e+01      -3.50000e+01
2   -5.06705e+15      -5.06705e+15
3    5.06705e+15       1.00049e+00
4    2.05716e+00       3.05765e+00

where the sum after step k=3 is mostly numerical noise.

It appears that one can find a condition on x and a for each k that
will make the denominator in the expression for rhok dissapear for
each k, although these expression become quite complex for k>2.  

Thus, the Gautschi equivalent method is probably not ideal for the
computation of the given continued fraction.

Best regards,
Hans

ps: I am not subscribed to the mailing list---please reply directly!

---------------------------------------------------------------------
Dr. Hans Ekkehard Plesser             Tel.  :           +47 6494 8832
Physics Section / ITF                 Fax   :           +47 6494 8810
Agricultural University of Norway     e-mail: hans.plesser@itf.nlh.no
N-1432 Ås, Norway		      WWW   :    arken.nlh.no/~itfhep
---------------------------------------------------------------------