[Bug-gsl] gsl_ran_gamma wrong for certain parameters

[Tino Kluge]

gsl_ran_gamma doesn't cope with parameters a,b producing a mean of 1
and very small standard deviations. See the programme below. I'm
using fedora core precompiled packages (gsl-1.6-2 and
gsl-devel-1.6-2).

---------------------------------------------------------------------

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main(char** argv, int args) {
 int    i;
 double a,b;
 gsl_rng * rgen = gsl_rng_alloc(gsl_rng_taus);
 
 a=5E9;
 b=1.0/a;
 printf("a=%e, b=%e, mean=%f, std dev=%f\n",a,b,a*b,sqrt(a)*b);
 for(i=1;i<10;i++){
   printf("%f\n",gsl_ran_gamma(rgen,a,b));
 }
 return 1;
}

---------------------------------------------------------------------

It produces seemingly correct samples till about a=4E9.

$ gcc gsl_gamma_test.c -lgsl -lgslcblas -lm -o gsl_gamma_test
$ ./gsl_gamma_test
a=4.000000e+09, b=2.500000e-10, mean=1.000000, std dev=0.000016
1.000020
1.000007
1.000000
1.000017
1.000023
0.999999
0.999995
0.999998
0.999994

---------------------------------------------------------------------

But completely changes behaviour at about a=5E9.

$ gcc gsl_gamma_test.c -lgsl -lgslcblas -lm -o gsl_gamma_test
$ ./gsl_gamma_test
a=5.000000e+09, b=2.000000e-10, mean=1.000000, std dev=0.000014
0.141013
0.141001
0.141012
0.141002
0.141006
0.141009
0.141012
0.141005
0.141002