Ziggurat: Fast Normal-Distributed and Exponential-Distributed Pseudo-Random Numbers

In producing my BMPF implementation, I did some rework of David Bateman's BSD-licensed re-implementation of Marsaglia and Tsang's Ziggurat Method of Gaussian (and exponential) pseudo-random number generation. This eventually grew into some fairly custom stuff. I replaced the various PRNG's lying around with the high-quality high-performance ISAAC PRNG. You can clone the GIT repository at git://svcs.cs.pdx.edu/git/ziggurat.git to get the source code and revision history.

I originally included a Ziggurat generator for the function (1-x)n used in the BMPF resampler. Unfortunately, this generator had a subtle mathematical bug, and has been replaced by a correct but slow implementation using uniform() and pow().

The Ziggurat method is a vast performance improvement over the better-known and simpler Box-Müller method in producing normal variates. The basic cost of producing a variate is essentially two table lookups, a floating-point multiply, a floating-point compare, and some amortized stuff; I can generate more than 15M Gaussian variates per second on my fast Intel Core II Duo box. The tables are small, about 6KB, and thus fit easily into L0 cache or onto a microcontroller—the code could be adapted for fixed-point with some effort.

My "improvements" include pre-compiled tables, a library wrapper to make it easier to link against, a header containing GCC inline functions for all the performance-critical bits, a fairly generic API, robust "make check" test code, random() / srandom() emulation, etc.