This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
gsl_ran_gamma
- From: John Lamb <J dot D dot Lamb at btinternet dot com>
- To: gsl-discuss at sources dot redhat dot com
- Date: Wed, 01 Sep 2004 18:37:02 +0100
- Subject: gsl_ran_gamma
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