Ticket #12099: 0001-Implement-ziggurat-for-exponential-random-dist.patch

File 0001-Implement-ziggurat-for-exponential-random-dist.patch, 24.9 KB (added by Jason Rhinelander <jason@…>, 7 years ago)

patch: ziggurat for exponential_distribution

  • new file include/boost/random/detail/int_float_pair.hpp 

    From 0b8f06ecfc6e5a1e488f121468e37c240471ead1 Mon Sep 17 00:00:00 2001
From: Jason Rhinelander <jason@imaginary.ca>
Date: Fri, 25 Mar 2016 14:48:56 -0400
Subject: [PATCH] Implement ziggurat for exponential random dist

This implementation closely follows the existing ziggurat implementation
in boost::random::normal_distribution; it uses a larger table (257
instead of 129) because, unlike normal, we don't lose a bit to
determining the sign.
---
 include/boost/random/detail/int_float_pair.hpp    | 120 ++++++++++++++
 include/boost/random/exponential_distribution.hpp | 186 +++++++++++++++++++++-
 include/boost/random/normal_distribution.hpp      |  93 +----------
 3 files changed, 305 insertions(+), 94 deletions(-)
 create mode 100644 include/boost/random/detail/int_float_pair.hpp

diff --git a/include/boost/random/detail/int_float_pair.hpp b/include/boost/random/detail/int_float_pair.hpp
new file mode 100644
index 0000000..160f82c
    - +  
     1/* boost random/detail/int_float_pair.hpp header file
     2 *
     3 * Copyright Jens Maurer 2000-2001
     4 * Copyright Steven Watanabe 2010-2011
     5 * Distributed under the Boost Software License, Version 1.0. (See
     6 * accompanying file LICENSE_1_0.txt or copy at
     7 * http://www.boost.org/LICENSE_1_0.txt)
     8 *
     9 * See http://www.boost.org for most recent version including documentation.
     10 *
     11 * $Id$
     12 *
     13 */
     14
     15#ifndef BOOST_RANDOM_DETAIL_INT_FLOAT_PAIR_HPP
     16#define BOOST_RANDOM_DETAIL_INT_FLOAT_PAIR_HPP
     17
     18#include <utility>
     19#include <boost/integer.hpp>
     20#include <boost/integer/integer_mask.hpp>
     21#include <boost/type_traits/make_unsigned.hpp>
     22#include <boost/type_traits/is_integral.hpp>
     23#include <boost/random/uniform_01.hpp>
     24#include <boost/random/uniform_int_distribution.hpp>
     25#include <boost/random/detail/signed_unsigned_tools.hpp>
     26#include <boost/random/detail/integer_log2.hpp>
     27
     28namespace boost {
     29namespace random {
     30namespace detail {
     31
     32template<class Engine>
     33inline typename boost::make_unsigned<typename Engine::result_type>::type
     34generate_one_digit(Engine& eng, std::size_t bits)
     35{
     36    typedef typename Engine::result_type base_result;
     37    typedef typename boost::make_unsigned<base_result>::type base_unsigned;
     38   
     39    base_unsigned range =
     40        detail::subtract<base_result>()((eng.max)(), (eng.min)());
     41    base_unsigned y0_mask = (base_unsigned(2) << (bits - 1)) - 1;
     42    base_unsigned y0 = (range + 1) & ~y0_mask;
     43    base_unsigned u;
     44    do {
     45        u = detail::subtract<base_result>()(eng(), (eng.min)());
     46    } while(y0 != 0 && u > base_unsigned(y0 - 1));
     47    return u & y0_mask;
     48}
     49
     50template<class RealType, std::size_t w, class Engine>
     51std::pair<RealType, int> generate_int_float_pair(Engine& eng, boost::mpl::true_)
     52{
     53    typedef typename Engine::result_type base_result;
     54    typedef typename boost::make_unsigned<base_result>::type base_unsigned;
     55   
     56    base_unsigned range =
     57        detail::subtract<base_result>()((eng.max)(), (eng.min)());
     58   
     59    std::size_t m =
     60        (range == (std::numeric_limits<base_unsigned>::max)()) ?
     61            std::numeric_limits<base_unsigned>::digits :
     62            detail::integer_log2(range + 1);
     63           
     64    int bucket = 0;
     65    // process as many full digits as possible into the int part
     66    for(std::size_t i = 0; i < w/m; ++i) {
     67        base_unsigned u = generate_one_digit(eng, m);
     68        bucket = (bucket << m) | u;
     69    }
     70    RealType r;
     71
     72    const std::size_t digits = std::numeric_limits<RealType>::digits;
     73    {
     74        base_unsigned u = generate_one_digit(eng, m);
     75        base_unsigned mask = (base_unsigned(1) << (w%m)) - 1;
     76        bucket = (bucket << (w%m)) | (mask & u);
     77        const RealType mult = RealType(1)/RealType(base_unsigned(1) << (m - w%m));
     78        // zero out unused bits
     79        if (m - w%m > digits) {
     80            u &= ~(base_unsigned(1) << (m - digits));
     81        }
     82        r = RealType(u >> (w%m)) * mult;
     83    }
     84    for(std::size_t i = m - w%m; i + m < digits; ++i) {
     85        base_unsigned u = generate_one_digit(eng, m);
     86        r += u;
     87        r *= RealType(0.5)/RealType(base_unsigned(1) << (m - 1));
     88    }
     89    if (m - w%m < digits)
     90    {
     91        const std::size_t remaining = (digits - m + w%m) % m;
     92        base_unsigned u = generate_one_digit(eng, m);
     93        r += u & ((base_unsigned(2) << (remaining - 1)) - 1);
     94        const RealType mult = RealType(0.5)/RealType(base_unsigned(1) << (remaining - 1));
     95        r *= mult;
     96    }
     97    return std::make_pair(r, bucket);
     98}
     99
     100template<class RealType, std::size_t w, class Engine>
     101inline std::pair<RealType, int> generate_int_float_pair(Engine& eng, boost::mpl::false_)
     102{
     103    int bucket = uniform_int_distribution<>(0, (1 << w) - 1)(eng);
     104    RealType r = uniform_01<RealType>()(eng);
     105    return std::make_pair(r, bucket);
     106}
     107
     108template<class RealType, std::size_t w, class Engine>
     109inline std::pair<RealType, int> generate_int_float_pair(Engine& eng)
     110{
     111    typedef typename Engine::result_type base_result;
     112    return generate_int_float_pair<RealType, w>(eng,
     113        boost::is_integral<base_result>());
     114}
     115
     116} // namespace detail
     117} // namespace random
     118} // namespace boost
     119
     120#endif // BOOST_RANDOM_DETAIL_INT_FLOAT_PAIR_HPP

  • include/boost/random/exponential_distribution.hpp 

    diff --git a/include/boost/random/exponential_distribution.hpp b/include/boost/random/exponential_distribution.hpp
index 53eaa93..73d96eb 100644
    a b  
    22 *
    33 * Copyright Jens Maurer 2000-2001
    44 * Copyright Steven Watanabe 2011
     5 * Copyright Jason Rhinelander 2016
    56 * Distributed under the Boost Software License, Version 1.0. (See
    67 * accompanying file LICENSE_1_0.txt or copy at
    78 * http://www.boost.org/LICENSE_1_0.txt)
    … …  
    2324#include <boost/limits.hpp>
    2425#include <boost/random/detail/config.hpp>
    2526#include <boost/random/detail/operators.hpp>
     27#include <boost/random/detail/int_float_pair.hpp>
    2628#include <boost/random/uniform_01.hpp>
    2729
    2830namespace boost {
    2931namespace random {
    3032
     33namespace detail {
     34
     35// tables for the ziggurat algorithm
     36template<class RealType>
     37struct exponential_table {
     38    static const RealType table_x[257];
     39    static const RealType table_y[257];
     40};
     41
     42template<class RealType>
     43const RealType exponential_table<RealType>::table_x[257] = {
     44    8.6971174701310497140, 7.6971174701310497140, 6.9410336293772123602, 6.4783784938325698538,
     45    6.1441646657724730491, 5.8821443157953997963, 5.6664101674540337371, 5.4828906275260628694,
     46    5.3230905057543986131, 5.1814872813015010392, 5.0542884899813047117, 4.9387770859012514838,
     47    4.8329397410251125881, 4.7352429966017412526, 4.6444918854200854873, 4.5597370617073515513,
     48    4.4802117465284221949, 4.4052876934735729805, 4.3344436803172730116, 4.2672424802773661873,
     49    4.2033137137351843802, 4.1423408656640511251, 4.0840513104082974638, 4.0282085446479365106,
     50    3.9746060666737884793, 3.9230625001354895926, 3.8734176703995089983, 3.8255294185223367372,
     51    3.7792709924116678992, 3.7345288940397975350, 3.6912010902374189454, 3.6491955157608538478,
     52    3.6084288131289096339, 3.5688252656483374051, 3.5303158891293438633, 3.4928376547740601814,
     53    3.4563328211327607625, 3.4207483572511205323, 3.3860354424603017887, 3.3521490309001100106,
     54    3.3190474709707487166, 3.2866921715990692095, 3.2550473085704501813, 3.2240795652862645207,
     55    3.1937579032122407483, 3.1640533580259734580, 3.1349388580844407393, 3.1063890623398246660,
     56    3.0783802152540905188, 3.0508900166154554479, 3.0238975044556767713, 2.9973829495161306949,
     57    2.9713277599210896472, 2.9457143948950456386, 2.9205262865127406647, 2.8957477686001416838,
     58    2.8713640120155362592, 2.8473609656351888266, 2.8237253024500354905, 2.8004443702507381944,
     59    2.7775061464397572041, 2.7548991965623453650, 2.7326126361947007411, 2.7106360958679293686,
     60    2.6889596887418041593, 2.6675739807732670816, 2.6464699631518093905, 2.6256390267977886123,
     61    2.6050729387408355373, 2.5847638202141406911, 2.5647041263169053687, 2.5448866271118700928,
     62    2.5253043900378279427, 2.5059507635285939648, 2.4868193617402096807, 2.4679040502973649846,
     63    2.4491989329782498908, 2.4306983392644199088, 2.4123968126888708336, 2.3942890999214583288,
     64    2.3763701405361408194, 2.3586350574093374601, 2.3410791477030346875, 2.3236978743901964559,
     65    2.3064868582835798692, 2.2894418705322694265, 2.2725588255531546952, 2.2558337743672190441,
     66    2.2392628983129087111, 2.2228425031110364013, 2.2065690132576635755, 2.1904389667232199235,
     67    2.1744490099377744673, 2.1585958930438856781, 2.1428764653998416425, 2.1272876713173679737,
     68    2.1118265460190418108, 2.0964902118017147637, 2.0812758743932248696, 2.0661808194905755036,
     69    2.0512024094685848641, 2.0363380802487695916, 2.0215853383189260770, 2.0069417578945183144,
     70    1.9924049782135764992, 1.9779727009573602295, 1.9636426877895480401, 1.9494127580071845659,
     71    1.9352807862970511135, 1.9212447005915276767, 1.9073024800183871196, 1.8934521529393077332,
     72    1.8796917950722108462, 1.8660195276928275962, 1.8524335159111751661, 1.8389319670188793980,
     73    1.8255131289035192212, 1.8121752885263901413, 1.7989167704602903934, 1.7857359354841254047,
     74    1.7726311792313049959, 1.7596009308890742369, 1.7466436519460739352, 1.7337578349855711926,
     75    1.7209420025219350428, 1.7081947058780575683, 1.6955145241015377061, 1.6829000629175537544,
     76    1.6703499537164519163, 1.6578628525741725325, 1.6454374393037234057, 1.6330724165359912048,
     77    1.6207665088282577216, 1.6085184617988580769, 1.5963270412864831349, 1.5841910325326886695,
     78    1.5721092393862294810, 1.5600804835278879161, 1.5481036037145133070, 1.5361774550410318943,
     79    1.5243009082192260050, 1.5124728488721167573, 1.5006921768428164936, 1.4889578055167456003,
     80    1.4772686611561334579, 1.4656236822457450411, 1.4540218188487932264, 1.4424620319720121876,
     81    1.4309432929388794104, 1.4194645827699828254, 1.4080248915695353509, 1.3966232179170417110,
     82    1.3852585682631217189, 1.3739299563284902176, 1.3626364025050864742, 1.3513769332583349176,
     83    1.3401505805295045843, 1.3289563811371163220, 1.3177933761763245480, 1.3066606104151739482,
     84    1.2955571316866007210, 1.2844819902750125450, 1.2734342382962410994, 1.2624129290696153434,
     85    1.2514171164808525098, 1.2404458543344064544, 1.2294981956938491599, 1.2185731922087903071,
     86    1.2076698934267612830, 1.1967873460884031665, 1.1859245934042023557, 1.1750806743109117687,
     87    1.1642546227056790397, 1.1534454666557748056, 1.1426522275816728928, 1.1318739194110786733,
     88    1.1211095477013306083, 1.1103581087274114281, 1.0996185885325976575, 1.0888899619385472598,
     89    1.0781711915113727024, 1.0674612264799681530, 1.0567590016025518414, 1.0460634359770445503,
     90    1.0353734317905289496, 1.0246878730026178052, 1.0140056239570971074, 1.0033255279156973717,
     91    0.99264640550727647009, 0.98196705308506317914, 0.97128624098390397896, 0.96060271166866709917,
     92    0.94991517776407659940, 0.93922231995526297952, 0.92852278474721113999, 0.91781518207004493915,
     93    0.90709808271569100600, 0.89637001558989069006, 0.88562946476175228052, 0.87487486629102585352,
     94    0.86410460481100519511, 0.85331700984237406386, 0.84251035181036928333, 0.83168283773427388393,
     95    0.82083260655441252290, 0.80995772405741906620, 0.79905617735548788109, 0.78812586886949324977,
     96    0.77716460975913043936, 0.76617011273543541328, 0.75513998418198289808, 0.74407171550050873971,
     97    0.73296267358436604916, 0.72181009030875689912, 0.71061105090965570413, 0.69936248110323266174,
     98    0.68806113277374858613, 0.67670356802952337911, 0.66528614139267855405, 0.65380497984766565353,
     99    0.64225596042453703448, 0.63063468493349100113, 0.61893645139487678178, 0.60715622162030085137,
     100    0.59528858429150359384, 0.58332771274877027785, 0.57126731653258903915, 0.55910058551154127652,
     101    0.54682012516331112550, 0.53441788123716615385, 0.52188505159213564105, 0.50921198244365495319,
     102    0.49638804551867159754, 0.48340149165346224782, 0.47023927508216945338, 0.45688684093142071279,
     103    0.44332786607355296305, 0.42954394022541129589, 0.41551416960035700100, 0.40121467889627836229,
     104    0.38661797794112021568, 0.37169214532991786118, 0.35639976025839443721, 0.34069648106484979674,
     105    0.32452911701691008547, 0.30783295467493287307, 0.29052795549123115167, 0.27251318547846547924,
     106    0.25365836338591284433, 0.23379048305967553619, 0.21267151063096745264, 0.18995868962243277774,
     107    0.16512762256418831796, 0.13730498094001380420, 0.10483850756582017915, 0.063852163815003480173,
     108    0
     109};
     110
     111template<class RealType>
     112const RealType exponential_table<RealType>::table_y[257] = {
     113    0, 0.00045413435384149675545, 0.00096726928232717452884, 0.0015362997803015723824,
     114    0.0021459677437189061793, 0.0027887987935740759640, 0.0034602647778369039855, 0.0041572951208337952532,
     115    0.0048776559835423925804, 0.0056196422072054831710, 0.0063819059373191794422, 0.0071633531836349841425,
     116    0.0079630774380170392396, 0.0087803149858089752347, 0.0096144136425022094101, 0.010464810181029979488,
     117    0.011331013597834597488, 0.012212592426255380661, 0.013109164931254991070, 0.014020391403181937334,
     118    0.014945968011691148079, 0.015885621839973162490, 0.016839106826039946359, 0.017806200410911360563,
     119    0.018786700744696029497, 0.019780424338009741737, 0.020787204072578117603, 0.021806887504283582125,
     120    0.022839335406385238829, 0.023884420511558170348, 0.024942026419731782971, 0.026012046645134218076,
     121    0.027094383780955798424, 0.028188948763978634421, 0.029295660224637394015, 0.030414443910466605492,
     122    0.031545232172893605499, 0.032687963508959533317, 0.033842582150874329031, 0.035009037697397411067,
     123    0.036187284781931419754, 0.037377282772959360128, 0.038578995503074859626, 0.039792391023374122670,
     124    0.041017441380414820816, 0.042254122413316231413, 0.043502413568888183301, 0.044762297732943280694,
     125    0.046033761076175166762, 0.047316792913181548703, 0.048611385573379494401, 0.049917534282706374944,
     126    0.051235237055126279830, 0.052564494593071689595, 0.053905310196046085104, 0.055257689676697038322,
     127    0.056621641283742874438, 0.057997175631200659098, 0.059384305633420264487, 0.060783046445479636051,
     128    0.062193415408540996150, 0.063615431999807331076, 0.065049117786753755036, 0.066494496385339779043,
     129    0.067951593421936607770, 0.069420436498728751675, 0.070901055162371828426, 0.072393480875708743023,
     130    0.073897746992364746308, 0.075413888734058408453, 0.076941943170480510100, 0.078481949201606426042,
     131    0.080033947542319910023, 0.081597980709237420930, 0.083174093009632380354, 0.084762330532368125386,
     132    0.086362741140756912277, 0.087975374467270219300, 0.089600281910032864534, 0.091237516631040162057,
     133    0.092887133556043546523, 0.094549189376055853718, 0.096223742550432800103, 0.097910853311492199618,
     134    0.099610583670637128826, 0.10132299742595363588, 0.10304816017125771553, 0.10478613930657016928,
     135    0.10653700405000166218, 0.10830082545103379867, 0.11007767640518539026, 0.11186763167005629731,
     136    0.11367076788274431301, 0.11548716357863353664, 0.11731689921155557057, 0.11916005717532768467,
     137    0.12101672182667483729, 0.12288697950954513498, 0.12477091858083096578, 0.12666862943751066518,
     138    0.12858020454522817870, 0.13050573846833078225, 0.13244532790138752023, 0.13439907170221363078,
     139    0.13636707092642885841, 0.13834942886358021406, 0.14034625107486244210, 0.14235764543247220043,
     140    0.14438372216063476473, 0.14642459387834493787, 0.14848037564386679222, 0.15055118500103990354,
     141    0.15263714202744286154, 0.15473836938446807312, 0.15685499236936522013, 0.15898713896931420572,
     142    0.16113493991759203183, 0.16329852875190180795, 0.16547804187493600915, 0.16767361861725019322,
     143    0.16988540130252766513, 0.17211353531532005700, 0.17435816917135348788, 0.17661945459049489581,
     144    0.17889754657247831241, 0.18119260347549629488, 0.18350478709776746150, 0.18583426276219711495,
     145    0.18818119940425430485, 0.19054576966319540013, 0.19292814997677133873, 0.19532852067956322315,
     146    0.19774706610509886464, 0.20018397469191127727, 0.20263943909370901930, 0.20511365629383770880,
     147    0.20760682772422204205, 0.21011915938898825914, 0.21265086199297827522, 0.21520215107537867786,
     148    0.21777324714870053264, 0.22036437584335949720, 0.22297576805812018050, 0.22560766011668406495,
     149    0.22826029393071670664, 0.23093391716962742173, 0.23362878343743333945, 0.23634515245705964715,
     150    0.23908329026244917002, 0.24184346939887722761, 0.24462596913189210901, 0.24743107566532763894,
     151    0.25025908236886230967, 0.25311029001562948171, 0.25598500703041538015, 0.25888354974901621678,
     152    0.26180624268936295243, 0.26475341883506220209, 0.26772541993204481808, 0.27072259679906003167,
     153    0.27374530965280298302, 0.27679392844851734458, 0.27986883323697289920, 0.28297041453878076010,
     154    0.28609907373707684673, 0.28925522348967773308, 0.29243928816189258772, 0.29565170428126120948,
     155    0.29889292101558177099, 0.30216340067569352897, 0.30546361924459023541, 0.30879406693456016794,
     156    0.31215524877417956945, 0.31554768522712893632, 0.31897191284495723773, 0.32242848495608914289,
     157    0.32591797239355619822, 0.32944096426413633091, 0.33299806876180896713, 0.33658991402867758144,
     158    0.34021714906678004560, 0.34388044470450243010, 0.34758049462163698567, 0.35131801643748334681,
     159    0.35509375286678745925, 0.35890847294874976196, 0.36276297335481777335, 0.36665807978151414890,
     160    0.37059464843514599421, 0.37457356761590215193, 0.37859575940958081092, 0.38266218149600982112,
     161    0.38677382908413768115, 0.39093173698479710717, 0.39513698183329015336, 0.39939068447523107877,
     162    0.40369401253053026739, 0.40804818315203238238, 0.41245446599716116772, 0.41691418643300289465,
     163    0.42142872899761659635, 0.42599954114303435739, 0.43062813728845883923, 0.43531610321563659758,
     164    0.44006510084235387501, 0.44487687341454851593, 0.44975325116275498919, 0.45469615747461548049,
     165    0.45970761564213768669, 0.46478975625042618067, 0.46994482528395999841, 0.47517519303737738299,
     166    0.48048336393045423016, 0.48587198734188493564, 0.49134386959403255500, 0.49690198724154955294,
     167    0.50254950184134769289, 0.50828977641064283495, 0.51412639381474855788, 0.52006317736823356823,
     168    0.52610421398361972602, 0.53225388026304326945, 0.53851687200286186590, 0.54489823767243963663,
     169    0.55140341654064131685, 0.55803828226258748140, 0.56480919291240022434, 0.57172304866482579008,
     170    0.57878735860284503057, 0.58601031847726802755, 0.59340090169173341521, 0.60096896636523224742,
     171    0.60872538207962206507, 0.61668218091520762326, 0.62485273870366592605, 0.63325199421436607968,
     172    0.64189671642726607018, 0.65080583341457104881, 0.66000084107899974178, 0.66950631673192477684,
     173    0.67935057226476538741, 0.68956649611707798890, 0.70019265508278816709, 0.71127476080507597882,
     174    0.72286765959357200702, 0.73503809243142351530, 0.74786862198519510742, 0.76146338884989624862,
     175    0.77595685204011559675, 0.79152763697249565519, 0.80842165152300838005, 0.82699329664305033399,
     176    0.84778550062398962096, 0.87170433238120363669, 0.90046992992574643800, 0.93814368086217467916,
     177    1
     178};
     179
     180template<class RealType = double>
     181struct unit_exponential_distribution
     182{
     183    template<class Engine>
     184    RealType operator()(Engine& eng) {
     185        const double * const table_x = exponential_table<double>::table_x;
     186        const double * const table_y = exponential_table<double>::table_y;
     187        RealType shift(0);
     188        for(;;) {
     189            std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
     190            int i = vals.second;
     191            RealType x = vals.first * RealType(table_x[i]);
     192            if(x < table_x[i + 1]) return shift + x;
     193            // For i=0 we need to generate from the tail, but because this is an exponential
     194            // distribution, the tail looks exactly like the body, so we can simply repeat with a
     195            // shift:
     196            if (i == 0) shift += RealType(table_x[1]);
     197            else {
     198                RealType y = RealType(table_y[i]) + uniform_01<RealType>()(eng) * RealType(table_y[i + 1] - table_y[i]);
     199                if (y < f(x)) return shift + x;
     200            }
     201        }
     202    }
     203    static RealType f(RealType x) {
     204        using std::exp;
     205        return exp(-x);
     206    }
     207};
     208
     209} // namespace detail
     210
     211
    31212/**
    32213 * The exponential distribution is a model of \random_distribution with
    33214 * a single parameter lambda.
    … … public:  
    127308    template<class Engine>
    128309    result_type operator()(Engine& eng) const
    129310    {
    130         using std::log;
    131         return -result_type(1) /
    132             _lambda * log(result_type(1)-uniform_01<RealType>()(eng));
     311        detail::unit_exponential_distribution<RealType> impl;
     312        return impl(eng) / _lambda;
    133313    }
    134314
    135315    /**

  • include/boost/random/normal_distribution.hpp 

    diff --git a/include/boost/random/normal_distribution.hpp b/include/boost/random/normal_distribution.hpp
index 6fffc4a..56fe10a 100644
    a b  
    2323#include <boost/assert.hpp>
    2424#include <boost/limits.hpp>
    2525#include <boost/static_assert.hpp>
    26 #include <boost/integer.hpp>
    27 #include <boost/integer/integer_mask.hpp>
    28 #include <boost/type_traits/is_integral.hpp>
    29 #include <boost/type_traits/make_unsigned.hpp>
    3026#include <boost/random/detail/config.hpp>
    3127#include <boost/random/detail/operators.hpp>
    32 #include <boost/random/detail/integer_log2.hpp>
     28#include <boost/random/detail/int_float_pair.hpp>
    3329#include <boost/random/uniform_01.hpp>
    34 #include <boost/random/uniform_int_distribution.hpp>
    3530#include <boost/random/exponential_distribution.hpp>
    3631#include <boost/mpl/bool.hpp>
    3732
    … … const RealType normal_table<RealType>::table_y[129] = {  
    121116    1
    122117};
    123118
    124 template<class Engine>
    125 inline typename boost::make_unsigned<typename Engine::result_type>::type
    126 generate_one_digit(Engine& eng, std::size_t bits)
    127 {
    128     typedef typename Engine::result_type base_result;
    129     typedef typename boost::make_unsigned<base_result>::type base_unsigned;
    130    
    131     base_unsigned range =
    132         detail::subtract<base_result>()((eng.max)(), (eng.min)());
    133     base_unsigned y0_mask = (base_unsigned(2) << (bits - 1)) - 1;
    134     base_unsigned y0 = (range + 1) & ~y0_mask;
    135     base_unsigned u;
    136     do {
    137         u = detail::subtract<base_result>()(eng(), (eng.min)());
    138     } while(y0 != 0 && u > base_unsigned(y0 - 1));
    139     return u & y0_mask;
    140 }
    141 
    142 template<class RealType, std::size_t w, class Engine>
    143 std::pair<RealType, int> generate_int_float_pair(Engine& eng, boost::mpl::true_)
    144 {
    145     typedef typename Engine::result_type base_result;
    146     typedef typename boost::make_unsigned<base_result>::type base_unsigned;
    147    
    148     base_unsigned range =
    149         detail::subtract<base_result>()((eng.max)(), (eng.min)());
    150    
    151     std::size_t m =
    152         (range == (std::numeric_limits<base_unsigned>::max)()) ?
    153             std::numeric_limits<base_unsigned>::digits :
    154             detail::integer_log2(range + 1);
    155            
    156     int bucket = 0;
    157     // process as many full digits as possible into the int part
    158     for(std::size_t i = 0; i < w/m; ++i) {
    159         base_unsigned u = generate_one_digit(eng, m);
    160         bucket = (bucket << m) | u;
    161     }
    162     RealType r;
    163 
    164     const std::size_t digits = std::numeric_limits<RealType>::digits;
    165     {
    166         base_unsigned u = generate_one_digit(eng, m);
    167         base_unsigned mask = (base_unsigned(1) << (w%m)) - 1;
    168         bucket = (bucket << (w%m)) | (mask & u);
    169         const RealType mult = RealType(1)/RealType(base_unsigned(1) << (m - w%m));
    170         // zero out unused bits
    171         if (m - w%m > digits) {
    172             u &= ~(base_unsigned(1) << (m - digits));
    173         }
    174         r = RealType(u >> (w%m)) * mult;
    175     }
    176     for(std::size_t i = m - w%m; i + m < digits; ++i) {
    177         base_unsigned u = generate_one_digit(eng, m);
    178         r += u;
    179         r *= RealType(0.5)/RealType(base_unsigned(1) << (m - 1));
    180     }
    181     if (m - w%m < digits)
    182     {
    183         const std::size_t remaining = (digits - m + w%m) % m;
    184         base_unsigned u = generate_one_digit(eng, m);
    185         r += u & ((base_unsigned(2) << (remaining - 1)) - 1);
    186         const RealType mult = RealType(0.5)/RealType(base_unsigned(1) << (remaining - 1));
    187         r *= mult;
    188     }
    189     return std::make_pair(r, bucket);
    190 }
    191 
    192 template<class RealType, std::size_t w, class Engine>
    193 inline std::pair<RealType, int> generate_int_float_pair(Engine& eng, boost::mpl::false_)
    194 {
    195     int bucket = uniform_int_distribution<>(0, (1 << w) - 1)(eng);
    196     RealType r = uniform_01<RealType>()(eng);
    197     return std::make_pair(r, bucket);
    198 }
    199 
    200 template<class RealType, std::size_t w, class Engine>
    201 inline std::pair<RealType, int> generate_int_float_pair(Engine& eng)
    202 {
    203     typedef typename Engine::result_type base_result;
    204     return generate_int_float_pair<RealType, w>(eng,
    205         boost::is_integral<base_result>());
    206 }
    207119
    208120template<class RealType = double>
    209121struct unit_normal_distribution
    … … struct unit_normal_distribution  
    240152    }
    241153};
    242154
    243 }
     155} // namespace detail
    244156
    245 // deterministic Box-Muller method, uses trigonometric functions
    246157
    247158/**
    248159 * Instantiations of class template normal_distribution model a