[Boost-bugs] [Boost C++ Libraries] #12099: Ziggurat implementation of boost::random::exponential_distribution

Subject: [Boost-bugs] [Boost C++ Libraries] #12099: Ziggurat implementation of boost::random::exponential_distribution
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2016-03-25 20:19:59

#12099: Ziggurat implementation of boost::random::exponential_distribution
-----------------------------------------+---------------------------
 Reporter:  Jason Rhinelander <jason@…>  |      Owner:  no-maintainer
     Type:  Patches                      |     Status:  new
Milestone:  To Be Determined             |  Component:  random
  Version:  Boost Development Trunk      |   Severity:  Optimization
 Keywords:                               |
-----------------------------------------+---------------------------
 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?

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/12099>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.