1 // random number generation (out of line) -*- C++ -*-
3 // Copyright (C) 2009 Free Software Foundation, Inc.
5 // This file is part of the GNU ISO C++ Library. This library is free
6 // software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the
8 // Free Software Foundation; either version 2, or (at your option)
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
16 // You should have received a copy of the GNU General Public License along
17 // with this library; see the file COPYING. If not, write to the Free
18 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
21 // As a special exception, you may use this file as part of a free software
22 // library without restriction. Specifically, if other files instantiate
23 // templates or use macros or inline functions from this file, or you compile
24 // this file and link it with other files to produce an executable, this
25 // file does not by itself cause the resulting executable to be covered by
26 // the GNU General Public License. This exception does not however
27 // invalidate any other reasons why the executable file might be covered by
28 // the GNU General Public License.
30 /** @file bits/random.tcc
31 * This is an internal header file, included by other library headers.
32 * You should not attempt to use it directly.
41 * (Further) implementation-space details.
45 // General case for x = (ax + c) mod m -- use Schrage's algorithm to
46 // avoid integer overflow.
48 // Because a and c are compile-time integral constants the compiler
49 // kindly elides any unreachable paths.
51 // Preconditions: a > 0, m > 0.
53 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
63 static const _Tp __q = __m / __a;
64 static const _Tp __r = __m % __a;
66 _Tp __t1 = __a * (__x % __q);
67 _Tp __t2 = __r * (__x / __q);
71 __x = __m - __t2 + __t1;
76 const _Tp __d = __m - __x;
86 // Special case for m == 0 -- use unsigned integer overflow as modulo
88 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
89 struct _Mod<_Tp, __a, __c, __m, true>
93 { return __a * __x + __c; }
95 } // namespace __detail
98 * Seeds the LCR with integral value @p __x0, adjusted so that the
99 * ring identity is never a member of the convergence set.
101 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
103 linear_congruential_engine<_UIntType, __a, __c, __m>::
106 if ((__detail::__mod<_UIntType, 1U, 0U, __m>(__c) == 0U)
107 && (__detail::__mod<_UIntType, 1U, 0U, __m>(__x0) == 0U))
108 _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(1U);
110 _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(__x0);
114 * Seeds the LCR engine with a value generated by @p __g.
116 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
118 linear_congruential_engine<_UIntType, __a, __c, __m>::
121 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
122 : (std::__lg(__m) + 31);
123 const _UIntType __k = __k0 / 32;
124 _UIntType __arr[__k + 3];
125 __q.generate(__arr + 0, __arr + 3);
126 _UIntType __factor = 1U;
127 _UIntType __sum = 0U;
128 for (size_t __i = 0; __i < __k; ++__i)
130 __sum += __arr[__i + 3] * __factor;
131 __factor *= __detail::_Shift<_UIntType, 32>::__value;
134 if ((__detail::__mod<_UIntType, 1U, 0U, __m>(__c) == 0U)
135 && (__detail::__mod<_UIntType, 1U, 0U, __m>(__sum) == 0U))
136 _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(1U);
138 _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(__sum);
142 * Seeds the LCR engine with a value generated by @p __g.
144 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
145 template<typename _Gen>
147 linear_congruential_engine<_UIntType, __a, __c, __m>::
148 seed(_Gen& __g, false_type)
150 _UIntType __x0 = __g();
151 if ((__detail::__mod<_UIntType, 1U, 0U, __m>(__c) == 0U)
152 && (__detail::__mod<_UIntType, 1U, 0U, __m>(__x0) == 0U))
153 _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(1U);
155 _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(__x0);
159 * Gets the next generated value in sequence.
161 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
162 typename linear_congruential_engine<_UIntType, __a, __c, __m>::
164 linear_congruential_engine<_UIntType, __a, __c, __m>::
167 _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
171 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
172 typename _CharT, typename _Traits>
173 std::basic_ostream<_CharT, _Traits>&
174 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
175 const linear_congruential_engine<_UIntType,
176 __a, __c, __m>& __lcr)
178 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
179 typedef typename __ostream_type::ios_base __ios_base;
181 const typename __ios_base::fmtflags __flags = __os.flags();
182 const _CharT __fill = __os.fill();
183 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
184 __os.fill(__os.widen(' '));
193 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
194 typename _CharT, typename _Traits>
195 std::basic_istream<_CharT, _Traits>&
196 operator>>(std::basic_istream<_CharT, _Traits>& __is,
197 linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
199 typedef std::basic_istream<_CharT, _Traits> __istream_type;
200 typedef typename __istream_type::ios_base __ios_base;
202 const typename __ios_base::fmtflags __flags = __is.flags();
203 __is.flags(__ios_base::dec);
212 template<typename _UIntType,
213 size_t __w, size_t __n, size_t __m, size_t __r,
214 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
215 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
218 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
219 __s, __b, __t, __c, __l, __f>::
220 seed(result_type __sd)
222 _M_x[0] = __detail::__mod<_UIntType, 1, 0,
223 __detail::_Shift<_UIntType, __w>::__value>(__sd);
225 for (size_t __i = 1; __i < state_size; ++__i)
227 _UIntType __x = _M_x[__i - 1];
228 __x ^= __x >> (__w - 2);
231 _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
232 __detail::_Shift<_UIntType, __w>::__value>(__x);
237 template<typename _UIntType,
238 size_t __w, size_t __n, size_t __m, size_t __r,
239 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
240 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
243 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
244 __s, __b, __t, __c, __l, __f>::
247 const _UIntType __upper_mask = (~_UIntType()) << __r;
248 const size_t __k = (__w + 31) / 32;
249 _UIntType __arr[__k * __n];
250 __q.generate(__arr + 0, __arr + __k * __n);
253 for (size_t __i = 0; __i < state_size; ++__i)
255 _UIntType __factor = 1U;
256 _UIntType __sum = 0U;
257 for (size_t __j = 0; __j < __k; ++__j)
259 __sum += __arr[__i * __k + __j] * __factor;
260 __factor *= __detail::_Shift<_UIntType, 32>::__value;
262 _M_x[__i] = __detail::__mod<_UIntType, 1U, 0U,
263 __detail::_Shift<_UIntType, __w>::__value>(__sum);
269 if ((_M_x[0] & __upper_mask) != 0U)
272 else if (_M_x[__i] != 0U)
277 _M_x[0] = __detail::_Shift<_UIntType, __w - 1U>::__value;
280 template<typename _UIntType, size_t __w,
281 size_t __n, size_t __m, size_t __r,
282 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
283 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
286 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
287 __s, __b, __t, __c, __l, __f>::result_type
288 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
289 __s, __b, __t, __c, __l, __f>::
292 // Reload the vector - cost is O(n) amortized over n calls.
293 if (_M_p >= state_size)
295 const _UIntType __upper_mask = (~_UIntType()) << __r;
296 const _UIntType __lower_mask = ~__upper_mask;
298 for (size_t __k = 0; __k < (__n - __m); ++__k)
300 _UIntType __y = ((_M_x[__k] & __upper_mask)
301 | (_M_x[__k + 1] & __lower_mask));
302 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
303 ^ ((__y & 0x01) ? __a : 0));
306 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
308 _UIntType __y = ((_M_x[__k] & __upper_mask)
309 | (_M_x[__k + 1] & __lower_mask));
310 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
311 ^ ((__y & 0x01) ? __a : 0));
314 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
315 | (_M_x[0] & __lower_mask));
316 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
317 ^ ((__y & 0x01) ? __a : 0));
321 // Calculate o(x(i)).
322 result_type __z = _M_x[_M_p++];
323 __z ^= (__z >> __u) & __d;
324 __z ^= (__z << __s) & __b;
325 __z ^= (__z << __t) & __c;
331 template<typename _UIntType, size_t __w,
332 size_t __n, size_t __m, size_t __r,
333 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
334 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
335 _UIntType __f, typename _CharT, typename _Traits>
336 std::basic_ostream<_CharT, _Traits>&
337 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
338 const mersenne_twister_engine<_UIntType, __w, __n, __m,
339 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
341 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
342 typedef typename __ostream_type::ios_base __ios_base;
344 const typename __ios_base::fmtflags __flags = __os.flags();
345 const _CharT __fill = __os.fill();
346 const _CharT __space = __os.widen(' ');
347 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
350 for (size_t __i = 0; __i < __n - 1; ++__i)
351 __os << __x._M_x[__i] << __space;
352 __os << __x._M_x[__n - 1];
359 template<typename _UIntType, size_t __w,
360 size_t __n, size_t __m, size_t __r,
361 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
362 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
363 _UIntType __f, typename _CharT, typename _Traits>
364 std::basic_istream<_CharT, _Traits>&
365 operator>>(std::basic_istream<_CharT, _Traits>& __is,
366 mersenne_twister_engine<_UIntType, __w, __n, __m,
367 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
369 typedef std::basic_istream<_CharT, _Traits> __istream_type;
370 typedef typename __istream_type::ios_base __ios_base;
372 const typename __ios_base::fmtflags __flags = __is.flags();
373 __is.flags(__ios_base::dec | __ios_base::skipws);
375 for (size_t __i = 0; __i < __n; ++__i)
376 __is >> __x._M_x[__i];
383 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
385 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
386 seed(result_type __value)
389 __value = default_seed;
391 std::linear_congruential_engine<result_type, 40014U, 0U, 2147483563U>
394 // I hope this is right. The "10000" tests work for the ranluxen.
395 const size_t __n = (word_size + 31) / 32;
397 for (size_t __i = 0; __i < long_lag; ++__i)
399 _UIntType __sum = 0U;
400 _UIntType __factor = 1U;
401 for (size_t __j = 0; __j < __n; ++__j)
403 __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
404 (__lcg()) * __factor;
405 __factor *= __detail::_Shift<_UIntType, 32>::__value;
407 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, _S_modulus>(__sum);
409 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
413 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
415 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
418 const size_t __n = (word_size + 31) / 32;
419 unsigned int __arr[long_lag + __n];
420 __q.generate(__arr + 0, __arr + long_lag + __n);
422 for (size_t __i = 0; __i < long_lag; ++__i)
424 _UIntType __sum = 0U;
425 _UIntType __factor = 1U;
426 for (size_t __j = 0; __j < __n; ++__j)
428 __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
429 (__arr[__i * __n + __j]) * __factor;
430 __factor *= __detail::_Shift<_UIntType, 32>::__value;
432 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, _S_modulus>(__sum);
434 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
438 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
439 typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
441 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
444 // Derive short lag index from current index.
445 long __ps = _M_p - short_lag;
449 // Calculate new x(i) without overflow or division.
450 // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
453 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
455 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
460 __xi = _S_modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
465 // Adjust current index to loop around in ring buffer.
466 if (++_M_p >= long_lag)
472 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
473 typename _CharT, typename _Traits>
474 std::basic_ostream<_CharT, _Traits>&
475 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
476 const subtract_with_carry_engine<_UIntType,
479 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
480 typedef typename __ostream_type::ios_base __ios_base;
482 const typename __ios_base::fmtflags __flags = __os.flags();
483 const _CharT __fill = __os.fill();
484 const _CharT __space = __os.widen(' ');
485 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
488 for (size_t __i = 0; __i < __r; ++__i)
489 __os << __x._M_x[__i] << __space;
490 __os << __x._M_carry;
497 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
498 typename _CharT, typename _Traits>
499 std::basic_istream<_CharT, _Traits>&
500 operator>>(std::basic_istream<_CharT, _Traits>& __is,
501 subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
503 typedef std::basic_ostream<_CharT, _Traits> __istream_type;
504 typedef typename __istream_type::ios_base __ios_base;
506 const typename __ios_base::fmtflags __flags = __is.flags();
507 __is.flags(__ios_base::dec | __ios_base::skipws);
509 for (size_t __i = 0; __i < __r; ++__i)
510 __is >> __x._M_x[__i];
511 __is >> __x._M_carry;
518 template<typename _RandomNumberEngine, size_t __p, size_t __r>
519 typename discard_block_engine<_RandomNumberEngine,
520 __p, __r>::result_type
521 discard_block_engine<_RandomNumberEngine, __p, __r>::
524 if (_M_n >= used_block)
526 _M_b.discard(block_size - _M_n);
533 template<typename _RandomNumberEngine, size_t __p, size_t __r,
534 typename _CharT, typename _Traits>
535 std::basic_ostream<_CharT, _Traits>&
536 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
537 const discard_block_engine<_RandomNumberEngine,
540 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
541 typedef typename __ostream_type::ios_base __ios_base;
543 const typename __ios_base::fmtflags __flags = __os.flags();
544 const _CharT __fill = __os.fill();
545 const _CharT __space = __os.widen(' ');
546 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
549 __os << __x.base() << __space << __x._M_n;
556 template<typename _RandomNumberEngine, size_t __p, size_t __r,
557 typename _CharT, typename _Traits>
558 std::basic_istream<_CharT, _Traits>&
559 operator>>(std::basic_istream<_CharT, _Traits>& __is,
560 discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
562 typedef std::basic_istream<_CharT, _Traits> __istream_type;
563 typedef typename __istream_type::ios_base __ios_base;
565 const typename __ios_base::fmtflags __flags = __is.flags();
566 __is.flags(__ios_base::dec | __ios_base::skipws);
568 __is >> __x._M_b >> __x._M_n;
575 template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
576 typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
578 independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
581 const long double __r = static_cast<long double>(this->max())
582 - static_cast<long double>(this->min()) + 1.0L;
583 const result_type __m = std::log10(__r) / std::log10(2.0L);
584 result_type __n, __n0, __y0, __y1, __s0, __s1;
585 for (size_t __i = 0; __i < 2; ++__i)
587 __n = (__w + __m - 1) / __m + __i;
588 __n0 = __n - __w % __n;
589 const result_type __w0 = __w / __n;
590 const result_type __w1 = __w0 + 1;
593 __y0 = __s0 * (__r / __s0);
594 __y1 = __s1 * (__r / __s1);
595 if (__r - __y0 <= __y0 / __n)
599 result_type __sum = 0;
600 for (size_t __k = 0; __k < __n0; ++__k)
604 __u = _M_b() - this->min();
609 for (size_t __k = __n0; __k < __n; ++__k)
613 __u = _M_b() - this->min();
622 template<typename _RandomNumberEngine, size_t __k>
623 typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
624 shuffle_order_engine<_RandomNumberEngine, __k>::
627 size_t __j = (__k * (_M_y - _M_b.min()))
628 / (_M_b.max() - _M_b.min() + 1);
635 template<typename _RandomNumberEngine, size_t __k,
636 typename _CharT, typename _Traits>
637 std::basic_ostream<_CharT, _Traits>&
638 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
639 const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
641 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
642 typedef typename __ostream_type::ios_base __ios_base;
644 const typename __ios_base::fmtflags __flags = __os.flags();
645 const _CharT __fill = __os.fill();
646 const _CharT __space = __os.widen(' ');
647 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
651 for (size_t __i = 0; __i < __k; ++__i)
652 __os << __space << __x._M_v[__i];
653 __os << __space << __x._M_y;
660 template<typename _RandomNumberEngine, size_t __k,
661 typename _CharT, typename _Traits>
662 std::basic_istream<_CharT, _Traits>&
663 operator>>(std::basic_istream<_CharT, _Traits>& __is,
664 shuffle_order_engine<_RandomNumberEngine, __k>& __x)
666 typedef std::basic_istream<_CharT, _Traits> __istream_type;
667 typedef typename __istream_type::ios_base __ios_base;
669 const typename __ios_base::fmtflags __flags = __is.flags();
670 __is.flags(__ios_base::dec | __ios_base::skipws);
673 for (size_t __i = 0; __i < __k; ++__i)
674 __is >> __x._M_v[__i];
682 template<typename _IntType>
683 template<typename _UniformRandomNumberGenerator>
684 typename uniform_int_distribution<_IntType>::result_type
685 uniform_int_distribution<_IntType>::
686 _M_call(_UniformRandomNumberGenerator& __urng,
687 result_type __min, result_type __max, true_type)
689 // XXX Must be fixed to work well for *arbitrary* __urng.max(),
690 // __urng.min(), __max, __min. Currently works fine only in the
691 // most common case __urng.max() - __urng.min() >= __max - __min,
692 // with __urng.max() > __urng.min() >= 0.
693 typedef typename __gnu_cxx::__add_unsigned<typename
694 _UniformRandomNumberGenerator::result_type>::__type __urntype;
695 typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
697 typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
699 __urntype, __utype>::__type __uctype;
703 const __urntype __urnmin = __urng.min();
704 const __urntype __urnmax = __urng.max();
705 const __urntype __urnrange = __urnmax - __urnmin;
706 const __uctype __urange = __max - __min;
707 const __uctype __udenom = (__urnrange <= __urange
708 ? 1 : __urnrange / (__urange + 1));
710 __ret = (__urntype(__urng()) - __urnmin) / __udenom;
711 while (__ret > __max - __min);
713 return __ret + __min;
716 template<typename _IntType, typename _CharT, typename _Traits>
717 std::basic_ostream<_CharT, _Traits>&
718 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
719 const uniform_int_distribution<_IntType>& __x)
721 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
722 typedef typename __ostream_type::ios_base __ios_base;
724 const typename __ios_base::fmtflags __flags = __os.flags();
725 const _CharT __fill = __os.fill();
726 const _CharT __space = __os.widen(' ');
727 __os.flags(__ios_base::scientific | __ios_base::left);
730 __os << __x.a() << __space << __x.b();
737 template<typename _IntType, typename _CharT, typename _Traits>
738 std::basic_istream<_CharT, _Traits>&
739 operator>>(std::basic_istream<_CharT, _Traits>& __is,
740 uniform_int_distribution<_IntType>& __x)
742 typedef std::basic_istream<_CharT, _Traits> __istream_type;
743 typedef typename __istream_type::ios_base __ios_base;
745 const typename __ios_base::fmtflags __flags = __is.flags();
746 __is.flags(__ios_base::dec | __ios_base::skipws);
750 __x.param(typename uniform_int_distribution<_IntType>::
751 param_type(__a, __b));
758 template<typename _RealType, typename _CharT, typename _Traits>
759 std::basic_ostream<_CharT, _Traits>&
760 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
761 const uniform_real_distribution<_RealType>& __x)
763 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
764 typedef typename __ostream_type::ios_base __ios_base;
766 const typename __ios_base::fmtflags __flags = __os.flags();
767 const _CharT __fill = __os.fill();
768 const std::streamsize __precision = __os.precision();
769 const _CharT __space = __os.widen(' ');
770 __os.flags(__ios_base::scientific | __ios_base::left);
772 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
774 __os << __x.a() << __space << __x.b();
778 __os.precision(__precision);
782 template<typename _RealType, typename _CharT, typename _Traits>
783 std::basic_istream<_CharT, _Traits>&
784 operator>>(std::basic_istream<_CharT, _Traits>& __is,
785 uniform_real_distribution<_RealType>& __x)
787 typedef std::basic_istream<_CharT, _Traits> __istream_type;
788 typedef typename __istream_type::ios_base __ios_base;
790 const typename __ios_base::fmtflags __flags = __is.flags();
791 __is.flags(__ios_base::skipws);
795 __x.param(typename uniform_real_distribution<_RealType>::
796 param_type(__a, __b));
803 template<typename _CharT, typename _Traits>
804 std::basic_ostream<_CharT, _Traits>&
805 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
806 const bernoulli_distribution& __x)
808 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
809 typedef typename __ostream_type::ios_base __ios_base;
811 const typename __ios_base::fmtflags __flags = __os.flags();
812 const _CharT __fill = __os.fill();
813 const std::streamsize __precision = __os.precision();
814 __os.flags(__ios_base::scientific | __ios_base::left);
815 __os.fill(__os.widen(' '));
816 __os.precision(std::numeric_limits<double>::digits10 + 1);
822 __os.precision(__precision);
827 template<typename _IntType>
828 template<typename _UniformRandomNumberGenerator>
829 typename geometric_distribution<_IntType>::result_type
830 geometric_distribution<_IntType>::
831 operator()(_UniformRandomNumberGenerator& __urng,
832 const param_type& __param)
834 // About the epsilon thing see this thread:
835 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
837 (1 - std::numeric_limits<double>::epsilon()) / 2;
838 // The largest _RealType convertible to _IntType.
840 std::numeric_limits<_IntType>::max() + __naf;
841 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
846 __cand = std::ceil(std::log(__aurng()) / __param._M_log_p);
847 while (__cand >= __thr);
849 return result_type(__cand + __naf);
852 template<typename _IntType,
853 typename _CharT, typename _Traits>
854 std::basic_ostream<_CharT, _Traits>&
855 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
856 const geometric_distribution<_IntType>& __x)
858 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
859 typedef typename __ostream_type::ios_base __ios_base;
861 const typename __ios_base::fmtflags __flags = __os.flags();
862 const _CharT __fill = __os.fill();
863 const std::streamsize __precision = __os.precision();
864 __os.flags(__ios_base::scientific | __ios_base::left);
865 __os.fill(__os.widen(' '));
866 __os.precision(std::numeric_limits<double>::digits10 + 1);
872 __os.precision(__precision);
876 template<typename _IntType,
877 typename _CharT, typename _Traits>
878 std::basic_istream<_CharT, _Traits>&
879 operator>>(std::basic_istream<_CharT, _Traits>& __is,
880 geometric_distribution<_IntType>& __x)
882 typedef std::basic_istream<_CharT, _Traits> __istream_type;
883 typedef typename __istream_type::ios_base __ios_base;
885 const typename __ios_base::fmtflags __flags = __is.flags();
886 __is.flags(__ios_base::skipws);
890 __x.param(typename geometric_distribution<_IntType>::param_type(__p));
896 template<typename _IntType>
897 template<typename _UniformRandomNumberGenerator>
898 typename negative_binomial_distribution<_IntType>::result_type
899 negative_binomial_distribution<_IntType>::
900 operator()(_UniformRandomNumberGenerator& __urng,
901 const param_type& __p)
903 typename gamma_distribution<>::param_type
904 __gamma_param(__p.k(), 1.0);
905 gamma_distribution<> __gamma(__gamma_param);
906 double __x = __gamma(__urng);
908 typename poisson_distribution<result_type>::param_type
909 __poisson_param(__x * __p.p() / (1.0 - __p.p()));
910 poisson_distribution<result_type> __poisson(__poisson_param);
911 result_type __m = __poisson(__urng);
916 template<typename _IntType, typename _CharT, typename _Traits>
917 std::basic_ostream<_CharT, _Traits>&
918 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
919 const negative_binomial_distribution<_IntType>& __x)
921 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
922 typedef typename __ostream_type::ios_base __ios_base;
924 const typename __ios_base::fmtflags __flags = __os.flags();
925 const _CharT __fill = __os.fill();
926 const std::streamsize __precision = __os.precision();
927 const _CharT __space = __os.widen(' ');
928 __os.flags(__ios_base::scientific | __ios_base::left);
929 __os.fill(__os.widen(' '));
930 __os.precision(std::numeric_limits<double>::digits10 + 1);
932 __os << __x.k() << __space << __x.p();
936 __os.precision(__precision);
940 template<typename _IntType, typename _CharT, typename _Traits>
941 std::basic_istream<_CharT, _Traits>&
942 operator>>(std::basic_istream<_CharT, _Traits>& __is,
943 negative_binomial_distribution<_IntType>& __x)
945 typedef std::basic_istream<_CharT, _Traits> __istream_type;
946 typedef typename __istream_type::ios_base __ios_base;
948 const typename __ios_base::fmtflags __flags = __is.flags();
949 __is.flags(__ios_base::skipws);
954 __x.param(typename negative_binomial_distribution<_IntType>::
955 param_type(__k, __p));
962 template<typename _IntType>
964 poisson_distribution<_IntType>::param_type::
967 #if _GLIBCXX_USE_C99_MATH_TR1
970 const double __m = std::floor(_M_mean);
971 _M_lm_thr = std::log(_M_mean);
972 _M_lfm = std::lgamma(__m + 1);
973 _M_sm = std::sqrt(__m);
975 const double __pi_4 = 0.7853981633974483096156608458198757L;
976 const double __dx = std::sqrt(2 * __m * std::log(32 * __m
978 _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
979 const double __cx = 2 * __m + _M_d;
980 _M_scx = std::sqrt(__cx / 2);
983 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
984 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
989 _M_lm_thr = std::exp(-_M_mean);
993 * A rejection algorithm when mean >= 12 and a simple method based
994 * upon the multiplication of uniform random variates otherwise.
995 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
999 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1000 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1002 template<typename _IntType>
1003 template<typename _UniformRandomNumberGenerator>
1004 typename poisson_distribution<_IntType>::result_type
1005 poisson_distribution<_IntType>::
1006 operator()(_UniformRandomNumberGenerator& __urng,
1007 const param_type& __param)
1009 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1011 #if _GLIBCXX_USE_C99_MATH_TR1
1012 if (__param.mean() >= 12)
1016 // See comments above...
1017 const double __naf =
1018 (1 - std::numeric_limits<double>::epsilon()) / 2;
1019 const double __thr =
1020 std::numeric_limits<_IntType>::max() + __naf;
1022 const double __m = std::floor(__param.mean());
1024 const double __spi_2 = 1.2533141373155002512078826424055226L;
1025 const double __c1 = __param._M_sm * __spi_2;
1026 const double __c2 = __param._M_c2b + __c1;
1027 const double __c3 = __c2 + 1;
1028 const double __c4 = __c3 + 1;
1030 const double __e178 = 1.0129030479320018583185514777512983L;
1031 const double __c5 = __c4 + __e178;
1032 const double __c = __param._M_cb + __c5;
1033 const double __2cx = 2 * (2 * __m + __param._M_d);
1035 bool __reject = true;
1038 const double __u = __c * __aurng();
1039 const double __e = -std::log(__aurng());
1045 const double __n = _M_nd(__urng);
1046 const double __y = -std::abs(__n) * __param._M_sm - 1;
1047 __x = std::floor(__y);
1048 __w = -__n * __n / 2;
1052 else if (__u <= __c2)
1054 const double __n = _M_nd(__urng);
1055 const double __y = 1 + std::abs(__n) * __param._M_scx;
1056 __x = std::ceil(__y);
1057 __w = __y * (2 - __y) * __param._M_1cx;
1058 if (__x > __param._M_d)
1061 else if (__u <= __c3)
1062 // NB: This case not in the book, nor in the Errata,
1063 // but should be ok...
1065 else if (__u <= __c4)
1067 else if (__u <= __c5)
1071 const double __v = -std::log(__aurng());
1072 const double __y = __param._M_d
1073 + __v * __2cx / __param._M_d;
1074 __x = std::ceil(__y);
1075 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1078 __reject = (__w - __e - __x * __param._M_lm_thr
1079 > __param._M_lfm - std::lgamma(__x + __m + 1));
1081 __reject |= __x + __m >= __thr;
1085 return result_type(__x + __m + __naf);
1091 double __prod = 1.0;
1095 __prod *= __aurng();
1098 while (__prod > __param._M_lm_thr);
1104 template<typename _IntType,
1105 typename _CharT, typename _Traits>
1106 std::basic_ostream<_CharT, _Traits>&
1107 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1108 const poisson_distribution<_IntType>& __x)
1110 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1111 typedef typename __ostream_type::ios_base __ios_base;
1113 const typename __ios_base::fmtflags __flags = __os.flags();
1114 const _CharT __fill = __os.fill();
1115 const std::streamsize __precision = __os.precision();
1116 const _CharT __space = __os.widen(' ');
1117 __os.flags(__ios_base::scientific | __ios_base::left);
1119 __os.precision(std::numeric_limits<double>::digits10 + 1);
1121 __os << __x.mean() << __space << __x._M_nd;
1123 __os.flags(__flags);
1125 __os.precision(__precision);
1129 template<typename _IntType,
1130 typename _CharT, typename _Traits>
1131 std::basic_istream<_CharT, _Traits>&
1132 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1133 poisson_distribution<_IntType>& __x)
1135 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1136 typedef typename __istream_type::ios_base __ios_base;
1138 const typename __ios_base::fmtflags __flags = __is.flags();
1139 __is.flags(__ios_base::skipws);
1142 __is >> __mean >> __x._M_nd;
1143 __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1145 __is.flags(__flags);
1150 template<typename _IntType>
1152 binomial_distribution<_IntType>::param_type::
1155 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1159 #if _GLIBCXX_USE_C99_MATH_TR1
1160 if (_M_t * __p12 >= 8)
1163 const double __np = std::floor(_M_t * __p12);
1164 const double __pa = __np / _M_t;
1165 const double __1p = 1 - __pa;
1167 const double __pi_4 = 0.7853981633974483096156608458198757L;
1168 const double __d1x =
1169 std::sqrt(__np * __1p * std::log(32 * __np
1170 / (81 * __pi_4 * __1p)));
1171 _M_d1 = std::round(std::max(1.0, __d1x));
1172 const double __d2x =
1173 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1174 / (__pi_4 * __pa)));
1175 _M_d2 = std::round(std::max(1.0, __d2x));
1178 const double __spi_2 = 1.2533141373155002512078826424055226L;
1179 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1180 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1181 _M_c = 2 * _M_d1 / __np;
1182 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1183 const double __a12 = _M_a1 + _M_s2 * __spi_2;
1184 const double __s1s = _M_s1 * _M_s1;
1185 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1187 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1188 const double __s2s = _M_s2 * _M_s2;
1189 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1190 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1191 _M_lf = (std::lgamma(__np + 1)
1192 + std::lgamma(_M_t - __np + 1));
1193 _M_lp1p = std::log(__pa / __1p);
1195 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1199 _M_q = -std::log(1 - __p12);
1202 template<typename _IntType>
1203 template<typename _UniformRandomNumberGenerator>
1204 typename binomial_distribution<_IntType>::result_type
1205 binomial_distribution<_IntType>::
1206 _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1210 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1215 const double __e = -std::log(__aurng());
1216 __sum += __e / (__t - __x);
1219 while (__sum <= _M_param._M_q);
1225 * A rejection algorithm when t * p >= 8 and a simple waiting time
1226 * method - the second in the referenced book - otherwise.
1227 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1231 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1232 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1234 template<typename _IntType>
1235 template<typename _UniformRandomNumberGenerator>
1236 typename binomial_distribution<_IntType>::result_type
1237 binomial_distribution<_IntType>::
1238 operator()(_UniformRandomNumberGenerator& __urng,
1239 const param_type& __param)
1242 const _IntType __t = __param.t();
1243 const _IntType __p = __param.p();
1244 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1245 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1248 #if _GLIBCXX_USE_C99_MATH_TR1
1249 if (!__param._M_easy)
1253 // See comments above...
1254 const double __naf =
1255 (1 - std::numeric_limits<double>::epsilon()) / 2;
1256 const double __thr =
1257 std::numeric_limits<_IntType>::max() + __naf;
1259 const double __np = std::floor(__t * __p12);
1262 const double __spi_2 = 1.2533141373155002512078826424055226L;
1263 const double __a1 = __param._M_a1;
1264 const double __a12 = __a1 + __param._M_s2 * __spi_2;
1265 const double __a123 = __param._M_a123;
1266 const double __s1s = __param._M_s1 * __param._M_s1;
1267 const double __s2s = __param._M_s2 * __param._M_s2;
1272 const double __u = __param._M_s * __aurng();
1278 const double __n = _M_nd(__urng);
1279 const double __y = __param._M_s1 * std::abs(__n);
1280 __reject = __y >= __param._M_d1;
1283 const double __e = -std::log(__aurng());
1284 __x = std::floor(__y);
1285 __v = -__e - __n * __n / 2 + __param._M_c;
1288 else if (__u <= __a12)
1290 const double __n = _M_nd(__urng);
1291 const double __y = __param._M_s2 * std::abs(__n);
1292 __reject = __y >= __param._M_d2;
1295 const double __e = -std::log(__aurng());
1296 __x = std::floor(-__y);
1297 __v = -__e - __n * __n / 2;
1300 else if (__u <= __a123)
1302 const double __e1 = -std::log(__aurng());
1303 const double __e2 = -std::log(__aurng());
1305 const double __y = __param._M_d1
1306 + 2 * __s1s * __e1 / __param._M_d1;
1307 __x = std::floor(__y);
1308 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1309 -__y / (2 * __s1s)));
1314 const double __e1 = -std::log(__aurng());
1315 const double __e2 = -std::log(__aurng());
1317 const double __y = __param._M_d2
1318 + 2 * __s2s * __e1 / __param._M_d2;
1319 __x = std::floor(-__y);
1320 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1324 __reject = __reject || __x < -__np || __x > __t - __np;
1327 const double __lfx =
1328 std::lgamma(__np + __x + 1)
1329 + std::lgamma(__t - (__np + __x) + 1);
1330 __reject = __v > __param._M_lf - __lfx
1331 + __x * __param._M_lp1p;
1334 __reject |= __x + __np >= __thr;
1338 __x += __np + __naf;
1340 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x));
1341 __ret = _IntType(__x) + __z;
1345 __ret = _M_waiting(__urng, __t);
1348 __ret = __t - __ret;
1352 template<typename _IntType,
1353 typename _CharT, typename _Traits>
1354 std::basic_ostream<_CharT, _Traits>&
1355 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1356 const binomial_distribution<_IntType>& __x)
1358 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1359 typedef typename __ostream_type::ios_base __ios_base;
1361 const typename __ios_base::fmtflags __flags = __os.flags();
1362 const _CharT __fill = __os.fill();
1363 const std::streamsize __precision = __os.precision();
1364 const _CharT __space = __os.widen(' ');
1365 __os.flags(__ios_base::scientific | __ios_base::left);
1367 __os.precision(std::numeric_limits<double>::digits10 + 1);
1369 __os << __x.t() << __space << __x.p()
1370 << __space << __x._M_nd;
1372 __os.flags(__flags);
1374 __os.precision(__precision);
1378 template<typename _IntType,
1379 typename _CharT, typename _Traits>
1380 std::basic_istream<_CharT, _Traits>&
1381 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1382 binomial_distribution<_IntType>& __x)
1384 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1385 typedef typename __istream_type::ios_base __ios_base;
1387 const typename __ios_base::fmtflags __flags = __is.flags();
1388 __is.flags(__ios_base::dec | __ios_base::skipws);
1392 __is >> __t >> __p >> __x._M_nd;
1393 __x.param(typename binomial_distribution<_IntType>::
1394 param_type(__t, __p));
1396 __is.flags(__flags);
1401 template<typename _RealType, typename _CharT, typename _Traits>
1402 std::basic_ostream<_CharT, _Traits>&
1403 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1404 const exponential_distribution<_RealType>& __x)
1406 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1407 typedef typename __ostream_type::ios_base __ios_base;
1409 const typename __ios_base::fmtflags __flags = __os.flags();
1410 const _CharT __fill = __os.fill();
1411 const std::streamsize __precision = __os.precision();
1412 __os.flags(__ios_base::scientific | __ios_base::left);
1413 __os.fill(__os.widen(' '));
1414 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1416 __os << __x.lambda();
1418 __os.flags(__flags);
1420 __os.precision(__precision);
1424 template<typename _RealType, typename _CharT, typename _Traits>
1425 std::basic_istream<_CharT, _Traits>&
1426 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1427 exponential_distribution<_RealType>& __x)
1429 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1430 typedef typename __istream_type::ios_base __ios_base;
1432 const typename __ios_base::fmtflags __flags = __is.flags();
1433 __is.flags(__ios_base::dec | __ios_base::skipws);
1437 __x.param(typename exponential_distribution<_RealType>::
1438 param_type(__lambda));
1440 __is.flags(__flags);
1445 template<typename _RealType>
1447 operator==(const normal_distribution<_RealType>& __d1,
1448 const normal_distribution<_RealType>& __d2)
1450 if (__d1._M_param == __d2._M_param)
1452 if (__d1._M_saved_available == __d2._M_saved_available)
1454 if (__d1._M_saved_available
1455 && __d1._M_saved == __d2._M_saved)
1457 else if(!__d1._M_saved_available)
1470 * Polar method due to Marsaglia.
1472 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1473 * New York, 1986, Ch. V, Sect. 4.4.
1475 template<typename _RealType>
1476 template<typename _UniformRandomNumberGenerator>
1477 typename normal_distribution<_RealType>::result_type
1478 normal_distribution<_RealType>::
1479 operator()(_UniformRandomNumberGenerator& __urng,
1480 const param_type& __param)
1483 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1486 if (_M_saved_available)
1488 _M_saved_available = false;
1493 result_type __x, __y, __r2;
1496 __x = result_type(2.0) * __aurng() - 1.0;
1497 __y = result_type(2.0) * __aurng() - 1.0;
1498 __r2 = __x * __x + __y * __y;
1500 while (__r2 > 1.0 || __r2 == 0.0);
1502 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1503 _M_saved = __x * __mult;
1504 _M_saved_available = true;
1505 __ret = __y * __mult;
1508 __ret = __ret * __param.stddev() + __param.mean();
1512 template<typename _RealType, typename _CharT, typename _Traits>
1513 std::basic_ostream<_CharT, _Traits>&
1514 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1515 const normal_distribution<_RealType>& __x)
1517 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1518 typedef typename __ostream_type::ios_base __ios_base;
1520 const typename __ios_base::fmtflags __flags = __os.flags();
1521 const _CharT __fill = __os.fill();
1522 const std::streamsize __precision = __os.precision();
1523 const _CharT __space = __os.widen(' ');
1524 __os.flags(__ios_base::scientific | __ios_base::left);
1526 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1528 __os << __x.mean() << __space << __x.stddev()
1529 << __space << __x._M_saved_available;
1530 if (__x._M_saved_available)
1531 __os << __space << __x._M_saved;
1533 __os.flags(__flags);
1535 __os.precision(__precision);
1539 template<typename _RealType, typename _CharT, typename _Traits>
1540 std::basic_istream<_CharT, _Traits>&
1541 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1542 normal_distribution<_RealType>& __x)
1544 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1545 typedef typename __istream_type::ios_base __ios_base;
1547 const typename __ios_base::fmtflags __flags = __is.flags();
1548 __is.flags(__ios_base::dec | __ios_base::skipws);
1550 double __mean, __stddev;
1551 __is >> __mean >> __stddev
1552 >> __x._M_saved_available;
1553 if (__x._M_saved_available)
1554 __is >> __x._M_saved;
1555 __x.param(typename normal_distribution<_RealType>::
1556 param_type(__mean, __stddev));
1558 __is.flags(__flags);
1563 template<typename _RealType>
1564 template<typename _UniformRandomNumberGenerator>
1565 typename lognormal_distribution<_RealType>::result_type
1566 lognormal_distribution<_RealType>::
1567 operator()(_UniformRandomNumberGenerator& __urng,
1568 const param_type& __p)
1570 _RealType __u, __v, __r2, __normal;
1571 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1576 // Choose x,y in uniform square (-1,-1) to (+1,+1).
1577 __u = 2 * __aurng() - 1;
1578 __v = 2 * __aurng() - 1;
1580 // See if it is in the unit circle.
1581 __r2 = __u * __u + __v * __v;
1583 while (__r2 > 1 || __r2 == 0);
1585 __normal = __u * std::sqrt(-2 * std::log(__r2) / __r2);
1587 return std::exp(__p.s() * __normal + __p.m());
1590 template<typename _RealType, typename _CharT, typename _Traits>
1591 std::basic_ostream<_CharT, _Traits>&
1592 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1593 const lognormal_distribution<_RealType>& __x)
1595 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1596 typedef typename __ostream_type::ios_base __ios_base;
1598 const typename __ios_base::fmtflags __flags = __os.flags();
1599 const _CharT __fill = __os.fill();
1600 const std::streamsize __precision = __os.precision();
1601 const _CharT __space = __os.widen(' ');
1602 __os.flags(__ios_base::scientific | __ios_base::left);
1604 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1606 __os << __x.m() << __space << __x.s();
1608 __os.flags(__flags);
1610 __os.precision(__precision);
1614 template<typename _RealType, typename _CharT, typename _Traits>
1615 std::basic_istream<_CharT, _Traits>&
1616 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1617 lognormal_distribution<_RealType>& __x)
1619 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1620 typedef typename __istream_type::ios_base __ios_base;
1622 const typename __ios_base::fmtflags __flags = __is.flags();
1623 __is.flags(__ios_base::dec | __ios_base::skipws);
1627 __x.param(typename lognormal_distribution<_RealType>::
1628 param_type(__m, __s));
1630 __is.flags(__flags);
1635 template<typename _RealType>
1636 template<typename _UniformRandomNumberGenerator>
1637 typename chi_squared_distribution<_RealType>::result_type
1638 chi_squared_distribution<_RealType>::
1639 operator()(_UniformRandomNumberGenerator& __urng,
1640 const param_type& __p)
1642 typename gamma_distribution<_RealType>::param_type
1643 __gamma_param(__p.n() / 2, 1.0);
1644 gamma_distribution<_RealType> __gamma(__gamma_param);
1645 return 2 * __gamma(__urng);
1648 template<typename _RealType, typename _CharT, typename _Traits>
1649 std::basic_ostream<_CharT, _Traits>&
1650 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1651 const chi_squared_distribution<_RealType>& __x)
1653 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1654 typedef typename __ostream_type::ios_base __ios_base;
1656 const typename __ios_base::fmtflags __flags = __os.flags();
1657 const _CharT __fill = __os.fill();
1658 const std::streamsize __precision = __os.precision();
1659 const _CharT __space = __os.widen(' ');
1660 __os.flags(__ios_base::scientific | __ios_base::left);
1662 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1666 __os.flags(__flags);
1668 __os.precision(__precision);
1672 template<typename _RealType, typename _CharT, typename _Traits>
1673 std::basic_istream<_CharT, _Traits>&
1674 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1675 chi_squared_distribution<_RealType>& __x)
1677 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1678 typedef typename __istream_type::ios_base __ios_base;
1680 const typename __ios_base::fmtflags __flags = __is.flags();
1681 __is.flags(__ios_base::dec | __ios_base::skipws);
1685 __x.param(typename chi_squared_distribution<_RealType>::
1688 __is.flags(__flags);
1693 template<typename _RealType>
1694 template<typename _UniformRandomNumberGenerator>
1695 typename cauchy_distribution<_RealType>::result_type
1696 cauchy_distribution<_RealType>::
1697 operator()(_UniformRandomNumberGenerator& __urng,
1698 const param_type& __p)
1700 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1709 return __p.a() + __p.b() * std::tan(M_PI * __u);
1712 template<typename _RealType, typename _CharT, typename _Traits>
1713 std::basic_ostream<_CharT, _Traits>&
1714 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1715 const cauchy_distribution<_RealType>& __x)
1717 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1718 typedef typename __ostream_type::ios_base __ios_base;
1720 const typename __ios_base::fmtflags __flags = __os.flags();
1721 const _CharT __fill = __os.fill();
1722 const std::streamsize __precision = __os.precision();
1723 const _CharT __space = __os.widen(' ');
1724 __os.flags(__ios_base::scientific | __ios_base::left);
1726 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1728 __os << __x.a() << __space << __x.b();
1730 __os.flags(__flags);
1732 __os.precision(__precision);
1736 template<typename _RealType, typename _CharT, typename _Traits>
1737 std::basic_istream<_CharT, _Traits>&
1738 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1739 cauchy_distribution<_RealType>& __x)
1741 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1742 typedef typename __istream_type::ios_base __ios_base;
1744 const typename __ios_base::fmtflags __flags = __is.flags();
1745 __is.flags(__ios_base::dec | __ios_base::skipws);
1749 __x.param(typename cauchy_distribution<_RealType>::
1750 param_type(__a, __b));
1752 __is.flags(__flags);
1757 template<typename _RealType>
1758 template<typename _UniformRandomNumberGenerator>
1759 typename fisher_f_distribution<_RealType>::result_type
1760 fisher_f_distribution<_RealType>::
1761 operator()(_UniformRandomNumberGenerator& __urng,
1762 const param_type& __p)
1764 gamma_distribution<_RealType> __gamma;
1765 _RealType __ym = __gamma(__urng,
1766 typename gamma_distribution<_RealType>::param_type(__p.m() / 2, 2));
1768 _RealType __yn = __gamma(__urng,
1769 typename gamma_distribution<_RealType>::param_type(__p.n() / 2, 2));
1771 return (__ym * __p.n()) / (__yn * __p.m());
1774 template<typename _RealType, typename _CharT, typename _Traits>
1775 std::basic_ostream<_CharT, _Traits>&
1776 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1777 const fisher_f_distribution<_RealType>& __x)
1779 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1780 typedef typename __ostream_type::ios_base __ios_base;
1782 const typename __ios_base::fmtflags __flags = __os.flags();
1783 const _CharT __fill = __os.fill();
1784 const std::streamsize __precision = __os.precision();
1785 const _CharT __space = __os.widen(' ');
1786 __os.flags(__ios_base::scientific | __ios_base::left);
1788 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1790 __os << __x.m() << __space << __x.n();
1792 __os.flags(__flags);
1794 __os.precision(__precision);
1798 template<typename _RealType, typename _CharT, typename _Traits>
1799 std::basic_istream<_CharT, _Traits>&
1800 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1801 fisher_f_distribution<_RealType>& __x)
1803 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1804 typedef typename __istream_type::ios_base __ios_base;
1806 const typename __ios_base::fmtflags __flags = __is.flags();
1807 __is.flags(__ios_base::dec | __ios_base::skipws);
1811 __x.param(typename fisher_f_distribution<_RealType>::
1812 param_type(__m, __n));
1814 __is.flags(__flags);
1820 // This could be operator() for a Gaussian distribution.
1822 template<typename _RealType>
1823 template<typename _UniformRandomNumberGenerator>
1824 typename student_t_distribution<_RealType>::result_type
1825 student_t_distribution<_RealType>::
1826 _M_gaussian(_UniformRandomNumberGenerator& __urng,
1827 const result_type __sigma)
1829 _RealType __x, __y, __r2;
1830 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1835 // Choose x,y in uniform square (-1,-1) to (+1,+1).
1836 __x = 2 * __aurng() - 1;
1837 __y = 2 * __aurng() - 1;
1839 // See if it is in the unit circle.
1840 __r2 = __x * __x + __y * __y;
1842 while (__r2 > 1 || __r2 == 0);
1844 // Box-Muller transform.
1845 return __sigma * __y * std::sqrt(-2 * std::log(__r2) / __r2);
1848 template<typename _RealType>
1849 template<typename _UniformRandomNumberGenerator>
1850 typename student_t_distribution<_RealType>::result_type
1851 student_t_distribution<_RealType>::
1852 operator()(_UniformRandomNumberGenerator& __urng,
1853 const param_type& __param)
1855 if (__param.n() <= 2.0)
1857 _RealType __y1 = _M_gaussian(__urng, 1.0);
1858 typename chi_squared_distribution<_RealType>::param_type
1859 __chisq_param(__param.n());
1860 chi_squared_distribution<_RealType> __chisq(__chisq_param);
1861 _RealType __y2 = __chisq(__urng);
1863 return __y1 / std::sqrt(__y2 / __param.n());
1867 _RealType __y1, __y2, __z;
1870 __y1 = _M_gaussian(__urng, 1.0);
1871 typename exponential_distribution<_RealType>::param_type
1872 __exp_param(1.0 / (__param.n() / 2.0 - 1.0));
1873 exponential_distribution<_RealType>
1874 __exponential(__exp_param);
1875 __y2 = __exponential(__urng);
1877 __z = __y1 * __y1 / (__param.n() - 2.0);
1879 while (1.0 - __z < 0.0 || std::exp(-__y2 - __z) > (1.0 - __z));
1881 // Note that there is a typo in Knuth's formula, the line below
1882 // is taken from the original paper of Marsaglia, Mathematics of
1883 // Computation, 34 (1980), p 234-256
1884 return __y1 / std::sqrt((1.0 - 2.0 / __param.n()) * (1.0 - __z));
1888 template<typename _RealType, typename _CharT, typename _Traits>
1889 std::basic_ostream<_CharT, _Traits>&
1890 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1891 const student_t_distribution<_RealType>& __x)
1893 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1894 typedef typename __ostream_type::ios_base __ios_base;
1896 const typename __ios_base::fmtflags __flags = __os.flags();
1897 const _CharT __fill = __os.fill();
1898 const std::streamsize __precision = __os.precision();
1899 const _CharT __space = __os.widen(' ');
1900 __os.flags(__ios_base::scientific | __ios_base::left);
1902 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1906 __os.flags(__flags);
1908 __os.precision(__precision);
1912 template<typename _RealType, typename _CharT, typename _Traits>
1913 std::basic_istream<_CharT, _Traits>&
1914 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1915 student_t_distribution<_RealType>& __x)
1917 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1918 typedef typename __istream_type::ios_base __ios_base;
1920 const typename __ios_base::fmtflags __flags = __is.flags();
1921 __is.flags(__ios_base::dec | __ios_base::skipws);
1925 __x.param(typename student_t_distribution<_RealType>::param_type(__n));
1927 __is.flags(__flags);
1932 template<typename _RealType>
1934 gamma_distribution<_RealType>::param_type::
1938 _M_l_d = std::sqrt(2 * _M_alpha - 1);
1940 _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
1945 * Cheng's rejection algorithm GB for alpha >= 1 and a modification
1946 * of Vaduva's rejection from Weibull algorithm due to Devroye for
1950 * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
1951 * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
1953 * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
1954 * and Composition Procedures." Math. Operationsforschung and Statistik,
1955 * Series in Statistics, 8, 545-576, 1977.
1957 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1958 * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
1960 template<typename _RealType>
1961 template<typename _UniformRandomNumberGenerator>
1962 typename gamma_distribution<_RealType>::result_type
1963 gamma_distribution<_RealType>::
1964 operator()(_UniformRandomNumberGenerator& __urng,
1965 const param_type& __param)
1968 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1972 const _RealType __alpha = __param.alpha();
1973 const _RealType __beta = __param.beta();
1977 const result_type __b = __alpha
1978 - result_type(1.3862943611198906188344642429163531L);
1979 const result_type __c = __alpha + __param._M_l_d;
1980 const result_type __1l = 1 / __param._M_l_d;
1983 const result_type __k = 2.5040773967762740733732583523868748L;
1987 const result_type __u = __aurng() / __beta;
1988 const result_type __v = __aurng() / __beta;
1990 const result_type __y = __1l * std::log(__v / (1 - __v));
1991 __x = __alpha * std::exp(__y);
1993 const result_type __z = __u * __v * __v;
1994 const result_type __r = __b + __c * __y - __x;
1996 __reject = __r < result_type(4.5) * __z - __k;
1998 __reject = __r < std::log(__z);
2004 const result_type __c = 1 / __alpha;
2008 const result_type __z = -std::log(__aurng() / __beta);
2009 const result_type __e = -std::log(__aurng() / __beta);
2011 __x = std::pow(__z, __c);
2013 __reject = __z + __e < __param._M_l_d + __x;
2018 return __beta * __x;
2021 template<typename _RealType, typename _CharT, typename _Traits>
2022 std::basic_ostream<_CharT, _Traits>&
2023 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2024 const gamma_distribution<_RealType>& __x)
2026 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2027 typedef typename __ostream_type::ios_base __ios_base;
2029 const typename __ios_base::fmtflags __flags = __os.flags();
2030 const _CharT __fill = __os.fill();
2031 const std::streamsize __precision = __os.precision();
2032 const _CharT __space = __os.widen(' ');
2033 __os.flags(__ios_base::scientific | __ios_base::left);
2035 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
2037 __os << __x.alpha() << __space << __x.beta();
2039 __os.flags(__flags);
2041 __os.precision(__precision);
2045 template<typename _RealType, typename _CharT, typename _Traits>
2046 std::basic_istream<_CharT, _Traits>&
2047 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2048 gamma_distribution<_RealType>& __x)
2050 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2051 typedef typename __istream_type::ios_base __ios_base;
2053 const typename __ios_base::fmtflags __flags = __is.flags();
2054 __is.flags(__ios_base::dec | __ios_base::skipws);
2056 _RealType __alpha, __beta;
2057 __is >> __alpha >> __beta;
2058 __x.param(typename gamma_distribution<_RealType>::
2059 param_type(__alpha, __beta));
2061 __is.flags(__flags);
2066 template<typename _RealType, typename _CharT, typename _Traits>
2067 std::basic_ostream<_CharT, _Traits>&
2068 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2069 const weibull_distribution<_RealType>& __x)
2071 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2072 typedef typename __ostream_type::ios_base __ios_base;
2074 const typename __ios_base::fmtflags __flags = __os.flags();
2075 const _CharT __fill = __os.fill();
2076 const std::streamsize __precision = __os.precision();
2077 const _CharT __space = __os.widen(' ');
2078 __os.flags(__ios_base::scientific | __ios_base::left);
2080 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
2082 __os << __x.a() << __space << __x.b();
2084 __os.flags(__flags);
2086 __os.precision(__precision);
2090 template<typename _RealType, typename _CharT, typename _Traits>
2091 std::basic_istream<_CharT, _Traits>&
2092 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2093 weibull_distribution<_RealType>& __x)
2095 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2096 typedef typename __istream_type::ios_base __ios_base;
2098 const typename __ios_base::fmtflags __flags = __is.flags();
2099 __is.flags(__ios_base::dec | __ios_base::skipws);
2103 __x.param(typename weibull_distribution<_RealType>::
2104 param_type(__a, __b));
2106 __is.flags(__flags);
2111 template<typename _RealType>
2112 template<typename _UniformRandomNumberGenerator>
2113 typename extreme_value_distribution<_RealType>::result_type
2114 extreme_value_distribution<_RealType>::
2115 operator()(_UniformRandomNumberGenerator& __urng,
2116 const param_type& __p)
2118 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2120 return __p.a() - __p.b() * std::log(-std::log(__aurng()));
2123 template<typename _RealType, typename _CharT, typename _Traits>
2124 std::basic_ostream<_CharT, _Traits>&
2125 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2126 const extreme_value_distribution<_RealType>& __x)
2128 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2129 typedef typename __ostream_type::ios_base __ios_base;
2131 const typename __ios_base::fmtflags __flags = __os.flags();
2132 const _CharT __fill = __os.fill();
2133 const std::streamsize __precision = __os.precision();
2134 const _CharT __space = __os.widen(' ');
2135 __os.flags(__ios_base::scientific | __ios_base::left);
2137 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
2139 __os << __x.a() << __space << __x.b();
2141 __os.flags(__flags);
2143 __os.precision(__precision);
2147 template<typename _RealType, typename _CharT, typename _Traits>
2148 std::basic_istream<_CharT, _Traits>&
2149 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2150 extreme_value_distribution<_RealType>& __x)
2152 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2153 typedef typename __istream_type::ios_base __ios_base;
2155 const typename __ios_base::fmtflags __flags = __is.flags();
2156 __is.flags(__ios_base::dec | __ios_base::skipws);
2160 __x.param(typename extreme_value_distribution<_RealType>::
2161 param_type(__a, __b));
2163 __is.flags(__flags);
2168 template<typename _IntType>
2170 discrete_distribution<_IntType>::param_type::
2173 if (_M_prob.size() < 2)
2176 _M_prob.push_back(1.0);
2180 double __sum = std::accumulate(_M_prob.begin(), _M_prob.end(), 0.0);
2181 // Now normalize the densities.
2182 std::transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2183 std::bind2nd(std::divides<double>(), __sum));
2184 // Accumulate partial sums.
2185 std::partial_sum(_M_prob.begin(), _M_prob.end(),
2186 std::back_inserter(_M_cp));
2187 // Make sure the last cumulative probablility is one.
2188 _M_cp[_M_cp.size() - 1] = 1.0;
2191 template<typename _IntType>
2192 template<typename _Func>
2193 discrete_distribution<_IntType>::param_type::
2194 param_type(size_t __nw, double __xmin, double __xmax,
2196 : _M_prob(), _M_cp()
2198 for (size_t __i = 0; __i < __nw; ++__i)
2200 const double __x = ((__nw - __i - 0.5) * __xmin
2201 + (__i + 0.5) * __xmax) / __nw;
2202 _M_prob.push_back(__fw(__x));
2208 template<typename _IntType>
2209 template<typename _UniformRandomNumberGenerator>
2210 typename discrete_distribution<_IntType>::result_type
2211 discrete_distribution<_IntType>::
2212 operator()(_UniformRandomNumberGenerator& __urng,
2213 const param_type& __param)
2215 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2218 const double __p = __aurng();
2219 auto __pos = std::lower_bound(__param._M_cp.begin(),
2220 __param._M_cp.end(), __p);
2221 if (__pos == __param._M_cp.end())
2223 const size_t __i = __pos - __param._M_cp.begin();
2228 template<typename _IntType, typename _CharT, typename _Traits>
2229 std::basic_ostream<_CharT, _Traits>&
2230 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2231 const discrete_distribution<_IntType>& __x)
2233 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2234 typedef typename __ostream_type::ios_base __ios_base;
2236 const typename __ios_base::fmtflags __flags = __os.flags();
2237 const _CharT __fill = __os.fill();
2238 const std::streamsize __precision = __os.precision();
2239 const _CharT __space = __os.widen(' ');
2240 __os.flags(__ios_base::scientific | __ios_base::left);
2242 __os.precision(std::numeric_limits<double>::digits10 + 1);
2244 std::vector<double> __prob = __x.probabilities();
2245 __os << __prob.size();
2246 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2247 __os << __space << *__dit;
2249 __os.flags(__flags);
2251 __os.precision(__precision);
2255 template<typename _IntType, typename _CharT, typename _Traits>
2256 std::basic_istream<_CharT, _Traits>&
2257 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2258 discrete_distribution<_IntType>& __x)
2260 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2261 typedef typename __istream_type::ios_base __ios_base;
2263 const typename __ios_base::fmtflags __flags = __is.flags();
2264 __is.flags(__ios_base::dec | __ios_base::skipws);
2269 std::vector<double> __prob_vec;
2270 for (; __n != 0; --__n)
2274 __prob_vec.push_back(__prob);
2277 __x.param(typename discrete_distribution<_IntType>::
2278 param_type(__prob_vec.begin(), __prob_vec.end()));
2280 __is.flags(__flags);
2285 template<typename _RealType>
2287 piecewise_constant_distribution<_RealType>::param_type::
2290 if (_M_int.size() < 2)
2293 _M_int.push_back(_RealType(0));
2294 _M_int.push_back(_RealType(1));
2297 _M_den.push_back(1.0);
2303 for (size_t __i = 0; __i < _M_den.size(); ++__i)
2305 __sum += _M_den[__i] * (_M_int[__i + 1] - _M_int[__i]);
2306 _M_cp.push_back(__sum);
2309 // Now normalize the densities...
2310 std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2311 std::bind2nd(std::divides<double>(), __sum));
2312 // ... and partial sums.
2313 std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
2314 std::bind2nd(std::divides<double>(), __sum));
2315 // Make sure the last cumulative probablility is one.
2316 _M_cp[_M_cp.size() - 1] = 1.0;
2319 template<typename _RealType>
2320 piecewise_constant_distribution<_RealType>::param_type::
2322 : _M_int(), _M_den(), _M_cp()
2323 { _M_initialize(); }
2325 template<typename _RealType>
2326 template<typename _InputIteratorB, typename _InputIteratorW>
2327 piecewise_constant_distribution<_RealType>::param_type::
2328 param_type(_InputIteratorB __bbegin,
2329 _InputIteratorB __bend,
2330 _InputIteratorW __wbegin)
2331 : _M_int(), _M_den(), _M_cp()
2335 _M_int.push_back(*__bbegin);
2337 if (__bbegin != __bend)
2339 _M_den.push_back(*__wbegin);
2343 while (__bbegin != __bend);
2348 template<typename _RealType>
2349 template<typename _Func>
2350 piecewise_constant_distribution<_RealType>::param_type::
2351 param_type(initializer_list<_RealType> __bil, _Func __fw)
2352 : _M_int(), _M_den(), _M_cp()
2354 for (auto __biter = __bil.begin(); __biter != __bil.end(); ++__biter)
2355 _M_int.push_back(*__biter);
2357 for (size_t __i = 0; __i < _M_int.size() - 1; ++__i)
2359 _RealType __x = 0.5 * (_M_int[__i] + _M_int[__i + 1]);
2360 _M_den.push_back(__fw(__x));
2366 template<typename _RealType>
2367 template<typename _Func>
2368 piecewise_constant_distribution<_RealType>::param_type::
2369 param_type(size_t __nw, _RealType __xmin, _RealType __xmax,
2371 : _M_int(), _M_den(), _M_cp()
2373 for (size_t __i = 0; __i <= __nw; ++__i)
2375 const _RealType __x = ((__nw - __i) * __xmin
2376 + __i * __xmax) / __nw;
2377 _M_int.push_back(__x);
2379 for (size_t __i = 0; __i < __nw; ++__i)
2381 const _RealType __x = ((__nw - __i - 0.5) * __xmin
2382 + (__i + 0.5) * __xmax) / __nw;
2383 _M_den.push_back(__fw(__x));
2389 template<typename _RealType>
2390 template<typename _UniformRandomNumberGenerator>
2391 typename piecewise_constant_distribution<_RealType>::result_type
2392 piecewise_constant_distribution<_RealType>::
2393 operator()(_UniformRandomNumberGenerator& __urng,
2394 const param_type& __param)
2396 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2399 const double __p = __aurng();
2400 auto __pos = std::lower_bound(__param._M_cp.begin(),
2401 __param._M_cp.end(), __p);
2402 const size_t __i = __pos - __param._M_cp.begin();
2404 return __param._M_int[__i]
2405 + (__p - __param._M_cp[__i]) / __param._M_den[__i];
2408 template<typename _RealType, typename _CharT, typename _Traits>
2409 std::basic_ostream<_CharT, _Traits>&
2410 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2411 const piecewise_constant_distribution<_RealType>& __x)
2413 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2414 typedef typename __ostream_type::ios_base __ios_base;
2416 const typename __ios_base::fmtflags __flags = __os.flags();
2417 const _CharT __fill = __os.fill();
2418 const std::streamsize __precision = __os.precision();
2419 const _CharT __space = __os.widen(' ');
2420 __os.flags(__ios_base::scientific | __ios_base::left);
2422 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
2424 std::vector<_RealType> __int = __x.intervals();
2425 __os << __int.size() - 1;
2427 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2428 __os << __space << *__xit;
2430 std::vector<double> __den = __x.densities();
2431 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2432 __os << __space << *__dit;
2434 __os.flags(__flags);
2436 __os.precision(__precision);
2440 template<typename _RealType, typename _CharT, typename _Traits>
2441 std::basic_istream<_CharT, _Traits>&
2442 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2443 piecewise_constant_distribution<_RealType>& __x)
2445 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2446 typedef typename __istream_type::ios_base __ios_base;
2448 const typename __ios_base::fmtflags __flags = __is.flags();
2449 __is.flags(__ios_base::dec | __ios_base::skipws);
2454 std::vector<_RealType> __int_vec;
2455 for (size_t __i = 0; __i <= __n; ++__i)
2459 __int_vec.push_back(__int);
2462 std::vector<double> __den_vec;
2463 for (size_t __i = 0; __i < __n; ++__i)
2467 __den_vec.push_back(__den);
2470 __x.param(typename piecewise_constant_distribution<_RealType>::
2471 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2473 __is.flags(__flags);
2478 template<typename _RealType>
2480 piecewise_linear_distribution<_RealType>::param_type::
2483 if (_M_int.size() < 2)
2486 _M_int.push_back(_RealType(0));
2487 _M_int.push_back(_RealType(1));
2490 _M_den.push_back(1.0);
2491 _M_den.push_back(1.0);
2497 for (size_t __i = 0; __i < _M_int.size() - 1; ++__i)
2499 const _RealType __delta = _M_int[__i + 1] - _M_int[__i];
2500 __sum += 0.5 * (_M_den[__i + 1] + _M_den[__i]) * __delta;
2501 _M_cp.push_back(__sum);
2502 _M_m.push_back((_M_den[__i + 1] - _M_den[__i]) / __delta);
2505 // Now normalize the densities...
2506 std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2507 std::bind2nd(std::divides<double>(),__sum));
2508 // ... and partial sums...
2509 std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
2510 std::bind2nd(std::divides<double>(), __sum));
2512 std::transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
2513 std::bind2nd(std::divides<double>(), __sum));
2514 // Make sure the last cumulative probablility is one.
2515 _M_cp[_M_cp.size() - 1] = 1.0;
2518 template<typename _RealType>
2519 piecewise_linear_distribution<_RealType>::param_type::
2521 : _M_int(), _M_den(), _M_cp(), _M_m()
2522 { _M_initialize(); }
2524 template<typename _RealType>
2525 template<typename _InputIteratorB, typename _InputIteratorW>
2526 piecewise_linear_distribution<_RealType>::param_type::
2527 param_type(_InputIteratorB __bbegin,
2528 _InputIteratorB __bend,
2529 _InputIteratorW __wbegin)
2530 : _M_int(), _M_den(), _M_cp(), _M_m()
2532 for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
2534 _M_int.push_back(*__bbegin);
2535 _M_den.push_back(*__wbegin);
2541 template<typename _RealType>
2542 template<typename _Func>
2543 piecewise_linear_distribution<_RealType>::param_type::
2544 param_type(initializer_list<_RealType> __bil, _Func __fw)
2545 : _M_int(), _M_den(), _M_cp(), _M_m()
2547 for (auto __biter = __bil.begin(); __biter != __bil.end(); ++__biter)
2549 _M_int.push_back(*__biter);
2550 _M_den.push_back(__fw(*__biter));
2556 template<typename _RealType>
2557 template<typename _Func>
2558 piecewise_linear_distribution<_RealType>::param_type::
2559 param_type(size_t __nw, _RealType __xmin, _RealType __xmax,
2561 : _M_int(), _M_den(), _M_cp(), _M_m()
2563 for (size_t __i = 0; __i <= __nw; ++__i)
2565 const _RealType __x = ((__nw - __i) * __xmin
2566 + __i * __xmax) / __nw;
2567 _M_int.push_back(__x);
2568 _M_den.push_back(__fw(__x));
2574 template<typename _RealType>
2575 template<typename _UniformRandomNumberGenerator>
2576 typename piecewise_linear_distribution<_RealType>::result_type
2577 piecewise_linear_distribution<_RealType>::
2578 operator()(_UniformRandomNumberGenerator& __urng,
2579 const param_type& __param)
2582 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2585 const double __p = __aurng();
2586 auto __pos = std::lower_bound(__param._M_cp.begin(),
2587 __param._M_cp.end(), __p);
2588 const size_t __i = __pos - __param._M_cp.begin();
2589 const double __a = 0.5 * __param._M_m[__i];
2590 const double __b = __param._M_den[__i];
2591 const double __c = __param._M_cp[__i];
2592 const double __q = -0.5 * (__b
2593 #if _GLIBCXX_USE_C99_MATH_TR1
2594 + std::copysign(std::sqrt(__b * __b
2595 - 4.0 * __a * __c), __b));
2597 + (__b < 0.0 ? -1.0 : 1.0)
2598 * std::sqrt(__b * __b - 4.0 * __a * __c)));
2600 const double __x0 = __param._M_int[__i];
2601 const double __x1 = __q / __a;
2602 const double __x2 = __c / __q;
2603 __x = std::max(__x0 + __x1, __x0 + __x2);
2608 template<typename _RealType, typename _CharT, typename _Traits>
2609 std::basic_ostream<_CharT, _Traits>&
2610 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2611 const piecewise_linear_distribution<_RealType>& __x)
2613 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2614 typedef typename __ostream_type::ios_base __ios_base;
2616 const typename __ios_base::fmtflags __flags = __os.flags();
2617 const _CharT __fill = __os.fill();
2618 const std::streamsize __precision = __os.precision();
2619 const _CharT __space = __os.widen(' ');
2620 __os.flags(__ios_base::scientific | __ios_base::left);
2622 __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
2624 std::vector<_RealType> __int = __x.intervals();
2625 __os << __int.size() - 1;
2627 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2628 __os << __space << *__xit;
2630 std::vector<double> __den = __x.densities();
2631 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2632 __os << __space << *__dit;
2634 __os.flags(__flags);
2636 __os.precision(__precision);
2640 template<typename _RealType, typename _CharT, typename _Traits>
2641 std::basic_istream<_CharT, _Traits>&
2642 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2643 piecewise_linear_distribution<_RealType>& __x)
2645 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2646 typedef typename __istream_type::ios_base __ios_base;
2648 const typename __ios_base::fmtflags __flags = __is.flags();
2649 __is.flags(__ios_base::dec | __ios_base::skipws);
2654 std::vector<_RealType> __int_vec;
2655 for (size_t __i = 0; __i <= __n; ++__i)
2659 __int_vec.push_back(__int);
2662 std::vector<double> __den_vec;
2663 for (size_t __i = 0; __i <= __n; ++__i)
2667 __den_vec.push_back(__den);
2670 __x.param(typename piecewise_linear_distribution<_RealType>::
2671 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2673 __is.flags(__flags);
2678 template<typename _IntType>
2679 seed_seq::seed_seq(std::initializer_list<_IntType> __il)
2681 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
2682 _M_v.push_back(__detail::__mod<result_type, 1, 0,
2683 __detail::_Shift<result_type, 32>::__value>(*__iter));
2686 template<typename _InputIterator>
2687 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
2689 for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
2690 _M_v.push_back(__detail::__mod<result_type, 1, 0,
2691 __detail::_Shift<result_type, 32>::__value>(*__iter));
2694 template<typename _RandomAccessIterator>
2696 seed_seq::generate(_RandomAccessIterator __begin,
2697 _RandomAccessIterator __end)
2699 typedef typename iterator_traits<_RandomAccessIterator>::value_type
2702 if (__begin == __end)
2705 std::fill(__begin, __end, __Type(0x8b8b8b8bU));
2707 const size_t __n = __end - __begin;
2708 const size_t __s = _M_v.size();
2709 const size_t __t = (__n >= 623) ? 11
2714 const size_t __p = (__n - __t) / 2;
2715 const size_t __q = __p + __t;
2716 const size_t __m = std::max(__s + 1, __n);
2718 for (size_t __k = 0; __k < __m; ++__k)
2720 __Type __arg = __begin[__k % __n]
2721 ^ __begin[(__k + __p) % __n]
2722 ^ __begin[(__k - 1) % __n];
2723 __Type __r1 = __arg ^ (__arg << 27);
2724 __r1 = __detail::__mod<__Type, 1664525U, 0U,
2725 __detail::_Shift<__Type, 32>::__value>(__r1);
2729 else if (__k <= __s)
2730 __r2 += __k % __n + _M_v[__k - 1];
2733 __r2 = __detail::__mod<__Type, 1U, 0U,
2734 __detail::_Shift<__Type, 32>::__value>(__r2);
2735 __begin[(__k + __p) % __n] += __r1;
2736 __begin[(__k + __q) % __n] += __r2;
2737 __begin[__k % __n] = __r2;
2740 for (size_t __k = __m; __k < __m + __n; ++__k)
2742 __Type __arg = __begin[__k % __n]
2743 + __begin[(__k + __p) % __n]
2744 + __begin[(__k - 1) % __n];
2745 __Type __r3 = __arg ^ (__arg << 27);
2746 __r3 = __detail::__mod<__Type, 1566083941U, 0U,
2747 __detail::_Shift<__Type, 32>::__value>(__r3);
2748 __Type __r4 = __r3 - __k % __n;
2749 __r4 = __detail::__mod<__Type, 1U, 0U,
2750 __detail::_Shift<__Type, 32>::__value>(__r4);
2751 __begin[(__k + __p) % __n] ^= __r4;
2752 __begin[(__k + __q) % __n] ^= __r3;
2753 __begin[__k % __n] = __r4;
2757 template<typename _RealType, size_t __bits,
2758 typename _UniformRandomNumberGenerator>
2760 generate_canonical(_UniformRandomNumberGenerator& __urng)
2763 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
2765 const long double __r = static_cast<long double>(__urng.max())
2766 - static_cast<long double>(__urng.min()) + 1.0L;
2767 const size_t __log2r = std::log10(__r) / std::log10(2.0L);
2768 size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
2769 _RealType __sum = _RealType(0);
2770 _RealType __tmp = _RealType(1);
2771 for (; __k != 0; --__k)
2773 __sum += _RealType(__urng() - __urng.min()) * __tmp;
2776 return __sum / __tmp;