gsl_ran_gamma

[John Lamb]

I looked at this a while back and noticed gsl still uses Ahrens-Dieter. 
I kew there were better methods and a colleague (possibly on this list) 
pointed out a more recent (exact) method by Marsaglia (as in the normal 
variate generator gsl uses) and Tsang, probably the most efficient known.

I tried replacing gsl_ran_gamma with Marsaglia-Tsang (which is actually 
shorter code) and found the following (b = 1):

a = 1.100000: 1532800 variates/second
a = 2.200000: 1529050 variates/second
a = 4.400000: 765578 variates/second
a = 8.800000: 764233 variates/second
a = 17.600000: 509839 variates/second
a = 35.200000: 508905 variates/second
a = 70.400000: 382204 variates/second
a = 140.800000: 381213 variates/second

c = -935079672.157271.

This compares with current gsl:

a = 1.100000: 794212 variates/second
a = 2.200000: 731158 variates/second
a = 4.400000: 354525 variates/second
a = 8.800000: 306560 variates/second
a = 17.600000: 199043 variates/second
a = 35.200000: 186258 variates/second
a = 70.400000: 138425 variates/second
a = 140.800000: 130590 variates/second

c = -934937048.192435.

I.e. Marsaglia and Tsang is obviously better. There's nothing in either 
code that would make it substantially more or less efficient on 
different architectures (I tested it on a pentium4 1.4 GHz, gcc 
-march=pentium4 -O2 -mfpmath=sse -msse -msse2 -malign-double, the rng 
was mt19937).

The tests don't cover a integer because the same (Erlang) method is used 
for both (though for Marsaglia-Tsang n = 12 would be too big).

Would it be worthwhile rewriting the gsl_ran_gamma code? Gamma is 
probably the most important variate after Gaussian and uniform and the 
basis for Chi-square, Erlang, etc.

The obvious downside of doing this is that gsl would generate a 
different sequence of random variates: hance the two values of c: these 
are signed sums of the variates and so very similar in value, but not 
identical.

-- 
JDL