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_min`

=0 and
`gsl_rng_max`

=2^{32}-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.

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))`

and`gsl_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))`

and`gsl_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.

- 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)

- 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

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.