https://svn.boost.org/trac10/ticket/12099 <p> 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. </p> <p> The attached patch (which I'll shortly also submit as a git pull request) changes boost::random::exponential_distribution to use the Ziggurat algorithm. </p> <p> 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). </p> <p> Since several other distributions rely (either directly or indirectly) on exponential_distribution, this should have notable speed improvements for drawing from them, as well. </p> <p> 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): </p> <ul><li>the values from an random::exponential_distribution object will change, obviously. </li><li>the values from dependent distributions will change (laplace, gamma, normal, and hyperexponential, directly, plus various distributions that make use of those). </li><li>the RNG state can sometimes advance more than it would have before (in particular, whenever we need to get a tail draw). </li></ul><p> Some other comments about the details of this patch: </p> <ul><li>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. </li><li>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. </li><li>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 <code>table_x[1]</code> 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? </li></ul> en-us /htdocs/site/boost.png https://svn.boost.org/trac10/ticket/12099 Trac 1.4.3 Jason Rhinelander <jason@…> Fri, 25 Mar 2016 20:22:43 GMT https://svn.boost.org/trac10/ticket/12099 https://svn.boost.org/trac10/ticket/12099 <ul> <li><strong>attachment</strong> → <span class="trac-field-new">0001-Implement-ziggurat-for-exponential-random-dist.patch</span> </li> </ul> <p> patch: ziggurat for exponential_distribution </p> Ticket Jason Rhinelander <jason@…> Fri, 25 Mar 2016 20:24:26 GMT https://svn.boost.org/trac10/ticket/12099 https://svn.boost.org/trac10/ticket/12099 <ul> <li><strong>attachment</strong> → <span class="trac-field-new">rexp_tables.cpp</span> </li> </ul> <p> normal and exponential table generator </p> Ticket Steven Watanabe Fri, 25 Mar 2016 20:47:53 GMT https://svn.boost.org/trac10/ticket/12099 https://svn.boost.org/trac10/ticket/12099 <ul> <li><strong>attachment</strong> → <span class="trac-field-new">scratch.cpp</span> </li> </ul> <p> Source for the tables used by boost::random::normal_distribution </p> Ticket