By Jochen Voss, last updated 2014-06-11
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.
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.
Since the Ziggurat method needs 32-bit random integers, it works best
for random number generators with
gsl_rng_max=232-1. The following table shows the
speed improvement for the different random number generators. For the
mt19937 the speed up factor is 3.9.
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.
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:
gsl_ran_gaussian_ziggurat(rng, M_SQRT2)is between -0.003 and +0.003
gsl_ran_gaussian_ziggurat(rng, M_SQRT2)is greater than 1 is 0.2397500611±0.0005
gsl_ran_gaussian_ziggurat(rng, M_SQRT2)is greater than 3 is 0.0169474268±0.0001
gsl_ran_gaussian_ziggurat(rng, sqrt(2.8))has mean 0±0.005
gsl_ran_gaussian_ziggurat(rng, sqrt(2.8))is smaller than -1 is 0.3085375387±0.001
All three methods,
gsl_ran_gaussian_ratio_method, pass these tests.
Copyright © 2014, Jochen Voss. All content on this website (including text, pictures, and any other original works), unless otherwise noted, is licensed under a Creative Commons Attribution-Share Alike 3.0 License.