1/* boost random/normal_distribution.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 * Revision history
14 * 2001-02-18 moved to individual header files
15 */
16
17#ifndef BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
18#define BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
19
20#include <boost/config/no_tr1/cmath.hpp>
21#include <istream>
22#include <iosfwd>
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>
37
38namespace boost {
39namespace random {
40
41namespace detail {
42
43// tables for the ziggurat algorithm
44template<class RealType>
45struct normal_table {
46 static const RealType table_x[129];
47 static const RealType table_y[129];
48};
49
50template<class RealType>
51const 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,
84 0
85};
86
87template<class RealType>
88const 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,
121 1
122};
123
124template<class Engine>
125inline typename boost::make_unsigned<typename Engine::result_type>::type
126generate_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
142template<class RealType, std::size_t w, class Engine>
143std::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
192template<class RealType, std::size_t w, class Engine>
193inline 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
200template<class RealType, std::size_t w, class Engine>
201inline 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}
207
208template<class RealType = double>
209struct unit_normal_distribution
210{
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;
215 for(;;) {
216 std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
217 int i = vals.second;
218 int sign = (i & 1) * 2 - 1;
219 i = i >> 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;
225 }
226 }
227 static RealType f(RealType x) {
228 using std::exp;
229 return exp(-x*x/2);
230 }
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]);
235 for(;;) {
236 RealType x = exponential(eng)/tail_start;
237 RealType y = exponential(eng);
238 if(2*y > x*x) return x + tail_start;
239 }
240 }
241};
242
243}
244
245// deterministic Box-Muller method, uses trigonometric functions
246
247/**
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}}
253 * \f$,
254 * where mean and sigma are the parameters of the distribution.
255 */
256template<class RealType = double>
257class normal_distribution
258{
259public:
260 typedef RealType input_type;
261 typedef RealType result_type;
262
263 class param_type {
264 public:
265 typedef normal_distribution distribution_type;
266
267 /**
268 * Constructs a @c param_type with a given mean and
269 * standard deviation.
270 *
271 * Requires: sigma >= 0
272 */
273 explicit param_type(RealType mean_arg = RealType(0.0),
274 RealType sigma_arg = RealType(1.0))
275 : _mean(mean_arg),
276 _sigma(sigma_arg)
277 {}
278
279 /** Returns the mean of the distribution. */
280 RealType mean() const { return _mean; }
281
282 /** Returns the standand deviation of the distribution. */
283 RealType sigma() const { return _sigma; }
284
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; }
288
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; }
292
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; }
296
297 /** Returns true if the two sets of parameters are the different. */
298 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
299
300 private:
301 RealType _mean;
302 RealType _sigma;
303 };
304
305 /**
306 * Constructs a @c normal_distribution object. @c mean and @c sigma are
307 * the parameters for the distribution.
308 *
309 * Requires: sigma >= 0
310 */
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)
314 {
315 BOOST_ASSERT(_sigma >= RealType(0));
316 }
317
318 /**
319 * Constructs a @c normal_distribution object from its parameters.
320 */
321 explicit normal_distribution(const param_type& parm)
322 : _mean(parm.mean()), _sigma(parm.sigma())
323 {}
324
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; }
329
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(); }
336
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)
341 {
342 _mean = parm.mean();
343 _sigma = parm.sigma();
344 }
345
346 /**
347 * Effects: Subsequent uses of the distribution do not depend
348 * on values produced by any engine prior to invoking reset.
349 */
350 void reset() { }
351
352 /** Returns a normal variate. */
353 template<class Engine>
354 result_type operator()(Engine& eng)
355 {
356 detail::unit_normal_distribution<RealType> impl;
357 return impl(eng) * _sigma + _mean;
358 }
359
360 /** Returns a normal variate with parameters specified by @c param. */
361 template<class URNG>
362 result_type operator()(URNG& urng, const param_type& parm)
363 {
364 return normal_distribution(parm)(urng);
365 }
366
367 /** Writes a @c normal_distribution to a @c std::ostream. */
368 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, normal_distribution, nd)
369 {
370 os << nd._mean << " " << nd._sigma;
371 return os;
372 }
373
374 /** Reads a @c normal_distribution from a @c std::istream. */
375 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, normal_distribution, nd)
376 {
377 is >> std::ws >> nd._mean >> std::ws >> nd._sigma;
378 return is;
379 }
380
381 /**
382 * Returns true if the two instances of @c normal_distribution will
383 * return identical sequences of values given equal generators.
384 */
385 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(normal_distribution, lhs, rhs)
386 {
387 return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma;
388 }
389
390 /**
391 * Returns true if the two instances of @c normal_distribution will
392 * return different sequences of values given equal generators.
393 */
394 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(normal_distribution)
395
396private:
397 RealType _mean, _sigma;
398
399};
400
401} // namespace random
402
403using random::normal_distribution;
404
405} // namespace boost
406
407#endif // BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
408

source code of boost/boost/random/normal_distribution.hpp