Problems with gsl_cdf_gamma_Pinv(0.01,a,b)

[Benjamin Redelings <benjamin_redelings@ncsu.edu>]

Hi,

I have been using GSL 1.10 on a 64-bit linux system (Debian/AMD64), and 
I've been trying to figure out how to handle problems like the following:

    // works fine
    gsl_cdf_gamma_Pinv(0.01, 11.888, 1.0);

    // fails to converge
    gsl_cdf_gamma_Pinv(0.01,11.887411491530846,1.0);

Since version 1.10, this failure is a hard failure - e.g. it returns 
GSL_NAN or it terminates. 

My explorations indicate that the range on 'a' in this function is not 
very large, in the grand scheme of things.  It seems like it can be 
expected to work if 'a' is in [0.25, 10] and also for some values outside.

Q1: I'm very much interested in advice on how to handle this kind of 
situation.  Should I expect this to work in GSL, or is there some reason 
why GSL does not want this to work (e.g. speed issues)?  Should I write 
some code myself to revert to a gaussian when a > 10?

Q2: I notice that the R code for this function is more intricate and 
seems to handle values of alpha that are very high (e.g. 500) and very 
low (0.01).  At least, qgamma(0.0001,500,1) and qgamma(0.0001,0.01,1) 
give me answers inside the R program.

The R algorithm seems to use the Chi^2 distribution (which is, I 
suppose, like the Normal distribution, effectively scale-independent).  
Is there any interest in expanding gsl_cdf_gamma_Pinv( ) to handle 
similar ranges to the code in R?

https://svn.r-project.org/R/trunk/src/nmath/qgamma.c

Thanks!

-BenRI