1 /* boost random/normal_distribution.hpp header file
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)
9 * See http://www.boost.org for most recent version including documentation.
14 * 2001-02-18 moved to individual header files
17 #ifndef BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
18 #define BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
20 #include <boost/config/no_tr1/cmath.hpp>
23 #include <boost/assert.hpp>
24 #include <boost/limits.hpp>
25 #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>
30 #include <boost/random/detail/config.hpp>
31 #include <boost/random/detail/operators.hpp>
32 #include <boost/random/detail/integer_log2.hpp>
33 #include <boost/random/uniform_01.hpp>
34 #include <boost/random/uniform_int_distribution.hpp>
35 #include <boost/random/exponential_distribution.hpp>
36 #include <boost/mpl/bool.hpp>
43 // tables for the ziggurat algorithm
44 template<class RealType>
46 static const RealType table_x[129];
47 static const RealType table_y[129];
50 template<class RealType>
51 const RealType normal_table<RealType>::table_x[129] = {
52 3.7130862467403632609, 3.4426198558966521214, 3.2230849845786185446, 3.0832288582142137009,
53 2.9786962526450169606, 2.8943440070186706210, 2.8231253505459664379, 2.7611693723841538514,
54 2.7061135731187223371, 2.6564064112581924999, 2.6109722484286132035, 2.5690336259216391328,
55 2.5300096723854666170, 2.4934545220919507609, 2.4590181774083500943, 2.4264206455302115930,
56 2.3954342780074673425, 2.3658713701139875435, 2.3375752413355307354, 2.3104136836950021558,
57 2.2842740596736568056, 2.2590595738653295251, 2.2346863955870569803, 2.2110814088747278106,
58 2.1881804320720206093, 2.1659267937448407377, 2.1442701823562613518, 2.1231657086697899595,
59 2.1025731351849988838, 2.0824562379877246441, 2.0627822745039633575, 2.0435215366506694976,
60 2.0246469733729338782, 2.0061338699589668403, 1.9879595741230607243, 1.9701032608497132242,
61 1.9525457295488889058, 1.9352692282919002011, 1.9182573008597320303, 1.9014946531003176140,
62 1.8849670357028692380, 1.8686611409895420085, 1.8525645117230870617, 1.8366654602533840447,
63 1.8209529965910050740, 1.8054167642140487420, 1.7900469825946189862, 1.7748343955807692457,
64 1.7597702248942318749, 1.7448461281083765085, 1.7300541605582435350, 1.7153867407081165482,
65 1.7008366185643009437, 1.6863968467734863258, 1.6720607540918522072, 1.6578219209482075462,
66 1.6436741568569826489, 1.6296114794646783962, 1.6156280950371329644, 1.6017183802152770587,
67 1.5878768648844007019, 1.5740982160167497219, 1.5603772223598406870, 1.5467087798535034608,
68 1.5330878776675560787, 1.5195095847593707806, 1.5059690368565502602, 1.4924614237746154081,
69 1.4789819769830978546, 1.4655259573357946276, 1.4520886428822164926, 1.4386653166774613138,
70 1.4252512545068615734, 1.4118417124397602509, 1.3984319141236063517, 1.3850170377251486449,
71 1.3715922024197322698, 1.3581524543224228739, 1.3446927517457130432, 1.3312079496576765017,
72 1.3176927832013429910, 1.3041418501204215390, 1.2905495919178731508, 1.2769102735516997175,
73 1.2632179614460282310, 1.2494664995643337480, 1.2356494832544811749, 1.2217602305309625678,
74 1.2077917504067576028, 1.1937367078237721994, 1.1795873846544607035, 1.1653356361550469083,
75 1.1509728421389760651, 1.1364898520030755352, 1.1218769225722540661, 1.1071236475235353980,
76 1.0922188768965537614, 1.0771506248819376573, 1.0619059636836193998, 1.0464709007525802629,
77 1.0308302360564555907, 1.0149673952392994716, 0.99886423348064351303, 0.98250080350276038481,
78 0.96585507938813059489, 0.94890262549791195381, 0.93161619660135381056, 0.91396525100880177644,
79 0.89591535256623852894, 0.87742742909771569142, 0.85845684317805086354, 0.83895221428120745572,
80 0.81885390668331772331, 0.79809206062627480454, 0.77658398787614838598, 0.75423066443451007146,
81 0.73091191062188128150, 0.70647961131360803456, 0.68074791864590421664, 0.65347863871504238702,
82 0.62435859730908822111, 0.59296294244197797913, 0.55869217837551797140, 0.52065603872514491759,
83 0.47743783725378787681, 0.42654798630330512490, 0.36287143102841830424, 0.27232086470466385065,
87 template<class RealType>
88 const RealType normal_table<RealType>::table_y[129] = {
89 0, 0.0026696290839025035092, 0.0055489952208164705392, 0.0086244844129304709682,
90 0.011839478657982313715, 0.015167298010672042468, 0.018592102737165812650, 0.022103304616111592615,
91 0.025693291936149616572, 0.029356317440253829618, 0.033087886146505155566, 0.036884388786968774128,
92 0.040742868074790604632, 0.044660862200872429800, 0.048636295860284051878, 0.052667401903503169793,
93 0.056752663481538584188, 0.060890770348566375972, 0.065080585213631873753, 0.069321117394180252601,
94 0.073611501884754893389, 0.077950982514654714188, 0.082338898242957408243, 0.086774671895542968998,
95 0.091257800827634710201, 0.09578784912257815216, 0.10036444102954554013, 0.10498725541035453978,
96 0.10965602101581776100, 0.11437051244988827452, 0.11913054670871858767, 0.12393598020398174246,
97 0.12878670619710396109, 0.13368265258464764118, 0.13862377998585103702, 0.14361008009193299469,
98 0.14864157424369696566, 0.15371831220958657066, 0.15884037114093507813, 0.16400785468492774791,
99 0.16922089223892475176, 0.17447963833240232295, 0.17978427212496211424, 0.18513499701071343216,
100 0.19053204032091372112, 0.19597565311811041399, 0.20146611007620324118, 0.20700370944187380064,
101 0.21258877307373610060, 0.21822164655637059599, 0.22390269938713388747, 0.22963232523430270355,
102 0.23541094226572765600, 0.24123899354775131610, 0.24711694751469673582, 0.25304529850976585934,
103 0.25902456739871074263, 0.26505530225816194029, 0.27113807914102527343, 0.27727350292189771153,
104 0.28346220822601251779, 0.28970486044581049771, 0.29600215684985583659, 0.30235482778947976274,
105 0.30876363800925192282, 0.31522938806815752222, 0.32175291587920862031, 0.32833509837615239609,
106 0.33497685331697116147, 0.34167914123501368412, 0.34844296754987246935, 0.35526938485154714435,
107 0.36215949537303321162, 0.36911445366827513952, 0.37613546951445442947, 0.38322381105988364587,
108 0.39038080824138948916, 0.39760785649804255208, 0.40490642081148835099, 0.41227804010702462062,
109 0.41972433205403823467, 0.42724699830956239880, 0.43484783025466189638, 0.44252871528024661483,
110 0.45029164368692696086, 0.45813871627287196483, 0.46607215269457097924, 0.47409430069824960453,
111 0.48220764633483869062, 0.49041482528932163741, 0.49871863547658432422, 0.50712205108130458951,
112 0.51562823824987205196, 0.52424057267899279809, 0.53296265938998758838, 0.54179835503172412311,
113 0.55075179312105527738, 0.55982741271069481791, 0.56902999107472161225, 0.57836468112670231279,
114 0.58783705444182052571, 0.59745315095181228217, 0.60721953663260488551, 0.61714337082656248870,
115 0.62723248525781456578, 0.63749547734314487428, 0.64794182111855080873, 0.65858200005865368016,
116 0.66942766735770616891, 0.68049184100641433355, 0.69178914344603585279, 0.70333609902581741633,
117 0.71515150742047704368, 0.72725691835450587793, 0.73967724368333814856, 0.75244155918570380145,
118 0.76558417390923599480, 0.77914608594170316563, 0.79317701178385921053, 0.80773829469612111340,
119 0.82290721139526200050, 0.83878360531064722379, 0.85550060788506428418, 0.87324304892685358879,
120 0.89228165080230272301, 0.91304364799203805999, 0.93628268170837107547, 0.96359969315576759960,
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)
128 typedef typename Engine::result_type base_result;
129 typedef typename boost::make_unsigned<base_result>::type base_unsigned;
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;
137 u = detail::subtract<base_result>()(eng(), (eng.min)());
138 } while(y0 != 0 && u > base_unsigned(y0 - 1));
142 template<class RealType, std::size_t w, class Engine>
143 std::pair<RealType, int> generate_int_float_pair(Engine& eng, boost::mpl::true_)
145 typedef typename Engine::result_type base_result;
146 typedef typename boost::make_unsigned<base_result>::type base_unsigned;
148 base_unsigned range =
149 detail::subtract<base_result>()((eng.max)(), (eng.min)());
152 (range == (std::numeric_limits<base_unsigned>::max)()) ?
153 std::numeric_limits<base_unsigned>::digits :
154 detail::integer_log2(range + 1);
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;
164 const std::size_t digits = std::numeric_limits<RealType>::digits;
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));
174 r = RealType(u >> (w%m)) * mult;
176 for(std::size_t i = m - w%m; i + m < digits; ++i) {
177 base_unsigned u = generate_one_digit(eng, m);
179 r *= RealType(0.5)/RealType(base_unsigned(1) << (m - 1));
181 if (m - w%m < digits)
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));
189 return std::make_pair(r, bucket);
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_)
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);
200 template<class RealType, std::size_t w, class Engine>
201 inline std::pair<RealType, int> generate_int_float_pair(Engine& eng)
203 typedef typename Engine::result_type base_result;
204 return generate_int_float_pair<RealType, w>(eng,
205 boost::is_integral<base_result>());
208 template<class RealType = double>
209 struct unit_normal_distribution
211 template<class Engine>
212 RealType operator()(Engine& eng) {
213 const double * const table_x = normal_table<double>::table_x;
214 const double * const table_y = normal_table<double>::table_y;
216 std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
218 int sign = (i & 1) * 2 - 1;
220 RealType x = vals.first * RealType(table_x[i]);
221 if(x < table_x[i + 1]) return x * sign;
222 if(i == 0) return generate_tail(eng) * sign;
223 RealType y = RealType(table_y[i]) + uniform_01<RealType>()(eng) * RealType(table_y[i + 1] - table_y[i]);
224 if (y < f(x)) return x * sign;
227 static RealType f(RealType x) {
231 template<class Engine>
232 RealType generate_tail(Engine& eng) {
233 boost::random::exponential_distribution<RealType> exponential;
234 const RealType tail_start = RealType(normal_table<double>::table_x[1]);
236 RealType x = exponential(eng)/tail_start;
237 RealType y = exponential(eng);
238 if(2*y > x*x) return x + tail_start;
245 // deterministic Box-Muller method, uses trigonometric functions
248 * Instantiations of class template normal_distribution model a
249 * \random_distribution. Such a distribution produces random numbers
250 * @c x distributed with probability density function
251 * \f$\displaystyle p(x) =
252 * \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
254 * where mean and sigma are the parameters of the distribution.
256 template<class RealType = double>
257 class normal_distribution
260 typedef RealType input_type;
261 typedef RealType result_type;
265 typedef normal_distribution distribution_type;
268 * Constructs a @c param_type with a given mean and
269 * standard deviation.
271 * Requires: sigma >= 0
273 explicit param_type(RealType mean_arg = RealType(0.0),
274 RealType sigma_arg = RealType(1.0))
279 /** Returns the mean of the distribution. */
280 RealType mean() const { return _mean; }
282 /** Returns the standand deviation of the distribution. */
283 RealType sigma() const { return _sigma; }
285 /** Writes a @c param_type to a @c std::ostream. */
286 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
287 { os << parm._mean << " " << parm._sigma ; return os; }
289 /** Reads a @c param_type from a @c std::istream. */
290 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
291 { is >> parm._mean >> std::ws >> parm._sigma; return is; }
293 /** Returns true if the two sets of parameters are the same. */
294 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
295 { return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma; }
297 /** Returns true if the two sets of parameters are the different. */
298 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
306 * Constructs a @c normal_distribution object. @c mean and @c sigma are
307 * the parameters for the distribution.
309 * Requires: sigma >= 0
311 explicit normal_distribution(const RealType& mean_arg = RealType(0.0),
312 const RealType& sigma_arg = RealType(1.0))
313 : _mean(mean_arg), _sigma(sigma_arg)
315 BOOST_ASSERT(_sigma >= RealType(0));
319 * Constructs a @c normal_distribution object from its parameters.
321 explicit normal_distribution(const param_type& parm)
322 : _mean(parm.mean()), _sigma(parm.sigma())
325 /** Returns the mean of the distribution. */
326 RealType mean() const { return _mean; }
327 /** Returns the standard deviation of the distribution. */
328 RealType sigma() const { return _sigma; }
330 /** Returns the smallest value that the distribution can produce. */
331 RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
332 { return -std::numeric_limits<RealType>::infinity(); }
333 /** Returns the largest value that the distribution can produce. */
334 RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
335 { return std::numeric_limits<RealType>::infinity(); }
337 /** Returns the parameters of the distribution. */
338 param_type param() const { return param_type(_mean, _sigma); }
339 /** Sets the parameters of the distribution. */
340 void param(const param_type& parm)
343 _sigma = parm.sigma();
347 * Effects: Subsequent uses of the distribution do not depend
348 * on values produced by any engine prior to invoking reset.
352 /** Returns a normal variate. */
353 template<class Engine>
354 result_type operator()(Engine& eng)
356 detail::unit_normal_distribution<RealType> impl;
357 return impl(eng) * _sigma + _mean;
360 /** Returns a normal variate with parameters specified by @c param. */
362 result_type operator()(URNG& urng, const param_type& parm)
364 return normal_distribution(parm)(urng);
367 /** Writes a @c normal_distribution to a @c std::ostream. */
368 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, normal_distribution, nd)
370 os << nd._mean << " " << nd._sigma;
374 /** Reads a @c normal_distribution from a @c std::istream. */
375 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, normal_distribution, nd)
377 is >> std::ws >> nd._mean >> std::ws >> nd._sigma;
382 * Returns true if the two instances of @c normal_distribution will
383 * return identical sequences of values given equal generators.
385 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(normal_distribution, lhs, rhs)
387 return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma;
391 * Returns true if the two instances of @c normal_distribution will
392 * return different sequences of values given equal generators.
394 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(normal_distribution)
397 RealType _mean, _sigma;
401 } // namespace random
403 using random::normal_distribution;
407 #endif // BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP