The Ziggurat Method for Generating Gaussian Random Numbers

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.

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=2³²-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:

All three methods, gsl_ran_gaussian_ziggurat, gsl_ran_gaussian and gsl_ran_gaussian_ratio_method, pass these tests.

Download

References