[Bug-gsl] gsl_cdf_chisq_Pinv bug/feature

[Bret Beck]

While testing my use of gsl_cdf_chisq_Pinv I ran into the following 
error message.

gsl: gammainv.c:105: ERROR: inverse failed to converge
Default GSL error handler invoked.
Aborted


The following code caused the problem.
-----------------------------------------------
---------------- Begin code --------------
-----------------------------------------------
#include <stdio.h>
#include <stdlib.h>

#include <gsl_rng.h>
#include <gsl_randist.h>
#include <gsl_cdf.h>

int main( int avgc, char **argv ) {

    double nu = 1.5, P, v;

    P = 1.93238e-01;
    v = gsl_cdf_chisq_Pinv(  P, nu );
    printf( "P(%.17e) = %.17e\n", P, v );

    P = 1.932383e-01;
    v = gsl_cdf_chisq_Pinv( P, nu );
    printf( "P(%.17e) = %.17e\n", P, v );

    P = 1.93238145206123590e-01;
    v = gsl_cdf_chisq_Pinv( P, nu);
    printf( "P(%.17e) = %.17e\n", P, v );

    exit( EXIT_SUCCESS );
}
-----------------------------------------------
----------------- End code ---------------
-----------------------------------------------

I looked into this and discovered the issue is in the routine 
cdf/gammainv.c. When nu = 1.93238145206123590e-01 and nu = 0.75, the 
number of allowed iterations at line 78 is too small. That is, I 
change the line from

    if (dP == 0.0 || n++ > 32)

to

    if (dP == 0.0 || n++ > 40)

and it no longer fails. In fact, the problem is most likely in line 
60, which is

      double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;

For nu = 1.5 and around P = 1.93238145206123590e-01, xg is very close 
to -sqrt( a ), and the guess x0 can jump from a = 0.75 to 0.

Thanks,
Bret Beck