The Ziggurat Method for Generating Gaussian Random Numbers
By Jochen Voss, last updated .
Introduction
The Ziggurat method for generating random numbers is a special case of the rejection method, where the density of the proposal distribution looks like the silhouette of a Ziggurat. The method can be used to produce very fast random number generators for distributions like the normal distribution or the exponential distribution. The main advantage comes from the fact that for a high percentage of the generated numbers no slow floating point operations are necessary.
The technical details of the Ziggurat method are, for example, contained in the article The Ziggurat Method for Generating Random Variables by G. Marsaglia and W. W. Tsang.
On this page I provide an implementation of the Ziggurat method to generate Gaussian random numbers. I hope to see this function to be included in the GSL library at some point of time. Therefore I compare the speed of my implementation with the two functions for Gaussian random number generation, which currently are included in the GSL library.
Implementation
My implementation of the Ziggurat method can be downloaded below. The file gauss.c contains the implementation, the other files are only there for illustration.
There are some differences between my implementation and the one in the article cited above: 1) I use only 128 steps instead of 256 to decrease the amount of required tabulated data. 2) I use an exponential distribution together with rejection for the tails of the base strip to simplify the implementation. Both changes seem to have no significant performance impact.
Speed Measurements
Since the Ziggurat method needs 32-bit random integers, it works best
for random number generators with gsl_rng_min=0 and
gsl_rng_max=232-1. The following table shows the
speed improvement for the different random number generators. For the
default generator mt19937 the speed up factor is 3.9.
| RNG | Ziggurat | Box-Muller | ratio method | improvement |
| borosh13 | 8.89 | 2.31 | 2.39 | 3.7 |
| cmrg | 2.84 | 1.34 | 1.3 | 2.1 |
| coveyou | 8.88 | 2.39 | 2.07 | 3.7 |
| fishman18 | 1.23 | 0.882 | 0.863 | 1.4 |
| fishman20 | 5.88 | 2.21 | 2.13 | 2.7 |
| fishman2x | 4.35 | 1.69 | 1.51 | 2.6 |
| gfsr4 | 9.03 | 2.19 | 2.29 | 4.0 |
| knuthran | 6.47 | 2.21 | 2.27 | 2.8 |
| knuthran2 | 0.662 | 0.527 | 0.499 | 1.3 |
| lecuyer21 | 6.61 | 2.22 | 2.17 | 3.0 |
| minstd | 5.61 | 2.16 | 2.1 | 2.6 |
| mrg | 4.38 | 1.73 | 1.66 | 2.5 |
| mt19937 | 8.28 | 2.01 | 2.14 | 3.9 |
| mt19937_1999 | 8.28 | 2.02 | 2.14 | 3.9 |
| mt19937_1998 | 8.27 | 2.01 | 2.12 | 3.9 |
| r250 | 9.1 | 2.2 | 2.05 | 4.1 |
| ran0 | 5.61 | 2.17 | 2.1 | 2.6 |
| ran1 | 5.94 | 2.1 | 2.04 | 2.8 |
| ran2 | 4.4 | 1.68 | 1.63 | 2.6 |
| ran3 | 6.69 | 2.05 | 2.12 | 3.2 |
| rand | 8.47 | 2.37 | 2.37 | 3.6 |
| rand48 | 8.14 | 0.887 | 0.84 | 9.2 |
| random128-bsd | 7.91 | 2.2 | 2.3 | 3.4 |
| random128-glibc2 | 7.72 | 2.18 | 2.31 | 3.3 |
| random128-libc5 | 7.94 | 2.2 | 2.31 | 3.4 |
| random256-bsd | 8.02 | 2.18 | 2.29 | 3.5 |
| random256-glibc2 | 8 | 2.2 | 2.28 | 3.5 |
| random256-libc5 | 7.72 | 2.18 | 2.29 | 3.4 |
| random32-bsd | 8.14 | 2.24 | 2.34 | 3.5 |
| random32-glibc2 | 7.93 | 2.25 | 2.34 | 3.4 |
| random32-libc5 | 8.19 | 2.24 | 2.34 | 3.5 |
| random64-bsd | 7.71 | 2.2 | 2.32 | 3.3 |
| random64-glibc2 | 7.71 | 2.21 | 2.33 | 3.3 |
| random64-libc5 | 7.68 | 2.2 | 2.32 | 3.3 |
| random8-bsd | 8.45 | 2.37 | 2.36 | 3.6 |
| random8-glibc2 | 8.45 | 2.37 | 2.36 | 3.6 |
| random8-libc5 | 8.47 | 2.37 | 2.37 | 3.6 |
Table 1. This table compares the speed of
Gaussian random number generation for different methods and different
random number generators. The first column gives the GSL name for the
underlying random number generator. Columns 2, 3, and 4 give the speed (in
million generated numbers per second) of the Ziggurat method
gsl_ran_gaussian_ziggurat, the Box-Muller method
gsl_ran_gaussian and the ratio method
gsl_ran_gaussian_ratio_method. All measurements were done on
a Pentium 4 with 2 GHz, the code was compiled with
gcc -O2. The last column gives the ratio of the values for
the Ziggurat method and for the better one of the other two methods.
Statistical Tests
In order to validate my implementation I performed several statistical tests on my random number generator. The program to perform these tests can be downloaded below. The (a little bit arbitrary) list of tests used is:
- the mean of
gsl_ran_gaussian_ziggurat(rng, M_SQRT2)is between -0.003 and +0.003 - the probability that
gsl_ran_gaussian_ziggurat(rng, M_SQRT2)is greater than 1 is 0.2397500611±0.0005 - the probability that
gsl_ran_gaussian_ziggurat(rng, M_SQRT2)is greater than 3 is 0.0169474268±0.0001 - the difference of
gsl_ran_gaussian_ziggurat(rng, sqrt(1.2))andgsl_ran_gaussian_ziggurat(rng, sqrt(2.8))has mean 0±0.005 - the probability that the difference between
gsl_ran_gaussian_ziggurat(rng, sqrt(1.2))andgsl_ran_gaussian_ziggurat(rng, sqrt(2.8))is smaller than -1 is 0.3085375387±0.001
All three methods, gsl_ran_gaussian_ziggurat,
gsl_ran_gaussian and
gsl_ran_gaussian_ratio_method, pass these tests.
Download
- gauss.c (implementation of the Ziggurat method)
- ziggurat.c (program to generate the auxiliary tables, only for illustration)
- timegauss.c (speed measurement, used to generate Table 1 above)
- testgauss.c (statistical tests)
References
- the GNU Scientific Library
- G. Marsaglia and W.W. Tsang:
The Ziggurat Method for Generating Random Variables.
Journal of Statistical Software, vol. 5, no. 8, pp. 1–7,
2000.
link