Ziggurat implementation of boost::random::exponential_distribution

[Jason Rhinelander <jason@…>]

The Ziggurat algorithm (implemented as of boost 1.56.0 in boost::random::normal_distribution) is suitable for any distribution with a decreasing density (or symmetric distribution with a decreasing density in one half, as in the case of normal_distribution). In particular the Marsaglia and Tsang (2000) paper that described the algorithm in the first place described both the normal distribution case and the exponential distribution.

The attached patch (which I'll shortly also submit as a git pull request) changes boost::random::exponential_distribution to use the Ziggurat algorithm.

In my testing, drawing double values, this ziggurat implementation improves the performance of the exponential distribution by about 3.9x (debian amd64, Intel Sandy Bridge CPU, g++ 5.3.1, using -O3) to 4.4x (debian amd64, Intel Haswell CPU, g++ 6 snapshot 20160313, using -O3). Using -march=native in both cases increases the relative gains further (to about 4.1x and 4.8x, respectively).

Since several other distributions rely (either directly or indirectly) on exponential_distribution, this should have notable speed improvements for drawing from them, as well.

This change has a couple of consequences, which are essentially the same as the consequences that changing normal_distribution to ziggurat had (and that some mailing list posts commented on):

• the values from an random::exponential_distribution object will change, obviously.
• the values from dependent distributions will change (laplace, gamma, normal, and hyperexponential, directly, plus various distributions that make use of those).
• the RNG state can sometimes advance more than it would have before (in particular, whenever we need to get a tail draw).

Some other comments about the details of this patch:

• I moved (without changing) the detail code for drawing an int/float pair from random/normal_distribution.hpp into a new random/detail/int_float_pair.hpp header, since it's needed by both the existing normal and new exponential ziggurat implementations, and has nothing specific to normal.
• I used a table of size 257 (versus normal's 129) so as to keep the uniform int draw as an 8-bit value (which is what normal uses); since exponential doesn't need a sign bit, that leaves the full 8 bits for selecting the table index. This difference (128 vs 256) also follows Marsaglia and Tsang's original reference implementations of normal and exponential.
• I couldn't find an exact source for the tables in normal_distribution.hpp, so I created a program to generate and output them. It can reproduce either the table for exponential or normal for any given table size; the normal output agrees *exactly* with the existing normal_distribution tables, and the table_x[1] value agrees (to 13 significant digits) with the Marsaglia and Tsang reference value. I'm not sure what to do with this program, though: is it suitable for dropping into random/extra?

[Jason Rhinelander <jason@…>]

patch: ziggurat for exponential_distribution

[Jason Rhinelander <jason@…>]

normal and exponential table generator

[Steven Watanabe]

I've tracked down the program that I used to generate the tables.

[Steven Watanabe]

Source for the tables used by boost::random::normal_distribution