A simple method for generating gamma variables

We offer a procedure for generating a gamma variate as the cube of a suitably scaled normal variate. It is fast and simple, assuming one has a fast way to generate normal variables. In brief: generate a normal variate x and a uniform variate U until ln(U) < 0.5x(2) + d - dv + dln(v), then return dv. Here, the gamma parameter is <alpha> greater than or equal to 1, and v = (1 + x/ root 9d)(3), with d = alpha - 1/3. The efficiency is high, exceeding 0.951, 0.981, 0.992, 0.996 at alpha = 1, 2, 4, 8. The procedure can be made to run faster by means of a simple squeeze that avoids the two logarithms most of the time: return dv if U < 1 - 0.0331x(4). We give a short C program for any <alpha> greater than or equal to 1, and show how to boost an alpha < 1 into an <alpha> > 1. The gamma procedure is particularly fast for C implementation if the normal variate is generated in-line, via the #define feature. We include such an inline version, leased on our ziggurat method. With it, and an inline uniform generator, gamma variates can be produced in 400MHz CPUs at better than 1.3 million per second, with the parameter alpha changing from call to call.