1 // random number generation (out of line) -*- C++ -*-
3 // Copyright (C) 2009-2019 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 3, 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 // Under Section 7 of GPL version 3, you are granted additional
17 // permissions described in the GCC Runtime Library Exception, version
18 // 3.1, as published by the Free Software Foundation.
20 // You should have received a copy of the GNU General Public License and
21 // a copy of the GCC Runtime Library Exception along with this program;
22 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23 // <http://www.gnu.org/licenses/>.
26 /** @file tr1/random.tcc
27 * This is an internal header file, included by other library headers.
28 * Do not attempt to use it directly. @headername{tr1/random}
31 #ifndef _GLIBCXX_TR1_RANDOM_TCC
32 #define _GLIBCXX_TR1_RANDOM_TCC 1
34 namespace std _GLIBCXX_VISIBILITY(default)
36 _GLIBCXX_BEGIN_NAMESPACE_VERSION
41 * (Further) implementation-space details.
45 // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
48 // Because a and c are compile-time integral constants the compiler kindly
49 // 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
97 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
99 linear_congruential<_UIntType, __a, __c, __m>::multiplier;
101 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
103 linear_congruential<_UIntType, __a, __c, __m>::increment;
105 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
107 linear_congruential<_UIntType, __a, __c, __m>::modulus;
110 * Seeds the LCR with integral value @p __x0, adjusted so that the
111 * ring identity is never a member of the convergence set.
113 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
115 linear_congruential<_UIntType, __a, __c, __m>::
116 seed(unsigned long __x0)
118 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
119 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
120 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
122 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
126 * Seeds the LCR engine with a value generated by @p __g.
128 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
131 linear_congruential<_UIntType, __a, __c, __m>::
132 seed(_Gen& __g, false_type)
134 _UIntType __x0 = __g();
135 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
136 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
137 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
139 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
143 * Gets the next generated value in sequence.
145 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
146 typename linear_congruential<_UIntType, __a, __c, __m>::result_type
147 linear_congruential<_UIntType, __a, __c, __m>::
150 _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
154 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
155 typename _CharT, typename _Traits>
156 std::basic_ostream<_CharT, _Traits>&
157 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
158 const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
160 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
161 typedef typename __ostream_type::ios_base __ios_base;
163 const typename __ios_base::fmtflags __flags = __os.flags();
164 const _CharT __fill = __os.fill();
165 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
166 __os.fill(__os.widen(' '));
175 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
176 typename _CharT, typename _Traits>
177 std::basic_istream<_CharT, _Traits>&
178 operator>>(std::basic_istream<_CharT, _Traits>& __is,
179 linear_congruential<_UIntType, __a, __c, __m>& __lcr)
181 typedef std::basic_istream<_CharT, _Traits> __istream_type;
182 typedef typename __istream_type::ios_base __ios_base;
184 const typename __ios_base::fmtflags __flags = __is.flags();
185 __is.flags(__ios_base::dec);
194 template<class _UIntType, int __w, int __n, int __m, int __r,
195 _UIntType __a, int __u, int __s,
196 _UIntType __b, int __t, _UIntType __c, int __l>
198 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
199 __b, __t, __c, __l>::word_size;
201 template<class _UIntType, int __w, int __n, int __m, int __r,
202 _UIntType __a, int __u, int __s,
203 _UIntType __b, int __t, _UIntType __c, int __l>
205 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
206 __b, __t, __c, __l>::state_size;
208 template<class _UIntType, int __w, int __n, int __m, int __r,
209 _UIntType __a, int __u, int __s,
210 _UIntType __b, int __t, _UIntType __c, int __l>
212 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
213 __b, __t, __c, __l>::shift_size;
215 template<class _UIntType, int __w, int __n, int __m, int __r,
216 _UIntType __a, int __u, int __s,
217 _UIntType __b, int __t, _UIntType __c, int __l>
219 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
220 __b, __t, __c, __l>::mask_bits;
222 template<class _UIntType, int __w, int __n, int __m, int __r,
223 _UIntType __a, int __u, int __s,
224 _UIntType __b, int __t, _UIntType __c, int __l>
226 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
227 __b, __t, __c, __l>::parameter_a;
229 template<class _UIntType, int __w, int __n, int __m, int __r,
230 _UIntType __a, int __u, int __s,
231 _UIntType __b, int __t, _UIntType __c, int __l>
233 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
234 __b, __t, __c, __l>::output_u;
236 template<class _UIntType, int __w, int __n, int __m, int __r,
237 _UIntType __a, int __u, int __s,
238 _UIntType __b, int __t, _UIntType __c, int __l>
240 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
241 __b, __t, __c, __l>::output_s;
243 template<class _UIntType, int __w, int __n, int __m, int __r,
244 _UIntType __a, int __u, int __s,
245 _UIntType __b, int __t, _UIntType __c, int __l>
247 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
248 __b, __t, __c, __l>::output_b;
250 template<class _UIntType, int __w, int __n, int __m, int __r,
251 _UIntType __a, int __u, int __s,
252 _UIntType __b, int __t, _UIntType __c, int __l>
254 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
255 __b, __t, __c, __l>::output_t;
257 template<class _UIntType, int __w, int __n, int __m, int __r,
258 _UIntType __a, int __u, int __s,
259 _UIntType __b, int __t, _UIntType __c, int __l>
261 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
262 __b, __t, __c, __l>::output_c;
264 template<class _UIntType, int __w, int __n, int __m, int __r,
265 _UIntType __a, int __u, int __s,
266 _UIntType __b, int __t, _UIntType __c, int __l>
268 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
269 __b, __t, __c, __l>::output_l;
271 template<class _UIntType, int __w, int __n, int __m, int __r,
272 _UIntType __a, int __u, int __s,
273 _UIntType __b, int __t, _UIntType __c, int __l>
275 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
276 __b, __t, __c, __l>::
277 seed(unsigned long __value)
279 _M_x[0] = __detail::__mod<_UIntType, 1, 0,
280 __detail::_Shift<_UIntType, __w>::__value>(__value);
282 for (int __i = 1; __i < state_size; ++__i)
284 _UIntType __x = _M_x[__i - 1];
285 __x ^= __x >> (__w - 2);
288 _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
289 __detail::_Shift<_UIntType, __w>::__value>(__x);
294 template<class _UIntType, int __w, int __n, int __m, int __r,
295 _UIntType __a, int __u, int __s,
296 _UIntType __b, int __t, _UIntType __c, int __l>
299 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
300 __b, __t, __c, __l>::
301 seed(_Gen& __gen, false_type)
303 for (int __i = 0; __i < state_size; ++__i)
304 _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
305 __detail::_Shift<_UIntType, __w>::__value>(__gen());
309 template<class _UIntType, int __w, int __n, int __m, int __r,
310 _UIntType __a, int __u, int __s,
311 _UIntType __b, int __t, _UIntType __c, int __l>
313 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
314 __b, __t, __c, __l>::result_type
315 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
316 __b, __t, __c, __l>::
319 // Reload the vector - cost is O(n) amortized over n calls.
320 if (_M_p >= state_size)
322 const _UIntType __upper_mask = (~_UIntType()) << __r;
323 const _UIntType __lower_mask = ~__upper_mask;
325 for (int __k = 0; __k < (__n - __m); ++__k)
327 _UIntType __y = ((_M_x[__k] & __upper_mask)
328 | (_M_x[__k + 1] & __lower_mask));
329 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
330 ^ ((__y & 0x01) ? __a : 0));
333 for (int __k = (__n - __m); __k < (__n - 1); ++__k)
335 _UIntType __y = ((_M_x[__k] & __upper_mask)
336 | (_M_x[__k + 1] & __lower_mask));
337 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
338 ^ ((__y & 0x01) ? __a : 0));
341 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
342 | (_M_x[0] & __lower_mask));
343 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
344 ^ ((__y & 0x01) ? __a : 0));
348 // Calculate o(x(i)).
349 result_type __z = _M_x[_M_p++];
351 __z ^= (__z << __s) & __b;
352 __z ^= (__z << __t) & __c;
358 template<class _UIntType, int __w, int __n, int __m, int __r,
359 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
360 _UIntType __c, int __l,
361 typename _CharT, typename _Traits>
362 std::basic_ostream<_CharT, _Traits>&
363 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
364 const mersenne_twister<_UIntType, __w, __n, __m,
365 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
367 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
368 typedef typename __ostream_type::ios_base __ios_base;
370 const typename __ios_base::fmtflags __flags = __os.flags();
371 const _CharT __fill = __os.fill();
372 const _CharT __space = __os.widen(' ');
373 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
376 for (int __i = 0; __i < __n - 1; ++__i)
377 __os << __x._M_x[__i] << __space;
378 __os << __x._M_x[__n - 1];
385 template<class _UIntType, int __w, int __n, int __m, int __r,
386 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
387 _UIntType __c, int __l,
388 typename _CharT, typename _Traits>
389 std::basic_istream<_CharT, _Traits>&
390 operator>>(std::basic_istream<_CharT, _Traits>& __is,
391 mersenne_twister<_UIntType, __w, __n, __m,
392 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
394 typedef std::basic_istream<_CharT, _Traits> __istream_type;
395 typedef typename __istream_type::ios_base __ios_base;
397 const typename __ios_base::fmtflags __flags = __is.flags();
398 __is.flags(__ios_base::dec | __ios_base::skipws);
400 for (int __i = 0; __i < __n; ++__i)
401 __is >> __x._M_x[__i];
408 template<typename _IntType, _IntType __m, int __s, int __r>
410 subtract_with_carry<_IntType, __m, __s, __r>::modulus;
412 template<typename _IntType, _IntType __m, int __s, int __r>
414 subtract_with_carry<_IntType, __m, __s, __r>::long_lag;
416 template<typename _IntType, _IntType __m, int __s, int __r>
418 subtract_with_carry<_IntType, __m, __s, __r>::short_lag;
420 template<typename _IntType, _IntType __m, int __s, int __r>
422 subtract_with_carry<_IntType, __m, __s, __r>::
423 seed(unsigned long __value)
428 std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
431 for (int __i = 0; __i < long_lag; ++__i)
432 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());
434 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
438 template<typename _IntType, _IntType __m, int __s, int __r>
441 subtract_with_carry<_IntType, __m, __s, __r>::
442 seed(_Gen& __gen, false_type)
444 const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;
446 for (int __i = 0; __i < long_lag; ++__i)
449 _UIntType __factor = 1;
450 for (int __j = 0; __j < __n; ++__j)
452 __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
453 (__gen()) * __factor;
454 __factor *= __detail::_Shift<_UIntType, 32>::__value;
456 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
458 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
462 template<typename _IntType, _IntType __m, int __s, int __r>
463 typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
464 subtract_with_carry<_IntType, __m, __s, __r>::
467 // Derive short lag index from current index.
468 int __ps = _M_p - short_lag;
472 // Calculate new x(i) without overflow or division.
473 // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
476 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
478 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
483 __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
488 // Adjust current index to loop around in ring buffer.
489 if (++_M_p >= long_lag)
495 template<typename _IntType, _IntType __m, int __s, int __r,
496 typename _CharT, typename _Traits>
497 std::basic_ostream<_CharT, _Traits>&
498 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
499 const subtract_with_carry<_IntType, __m, __s, __r>& __x)
501 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
502 typedef typename __ostream_type::ios_base __ios_base;
504 const typename __ios_base::fmtflags __flags = __os.flags();
505 const _CharT __fill = __os.fill();
506 const _CharT __space = __os.widen(' ');
507 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
510 for (int __i = 0; __i < __r; ++__i)
511 __os << __x._M_x[__i] << __space;
512 __os << __x._M_carry;
519 template<typename _IntType, _IntType __m, int __s, int __r,
520 typename _CharT, typename _Traits>
521 std::basic_istream<_CharT, _Traits>&
522 operator>>(std::basic_istream<_CharT, _Traits>& __is,
523 subtract_with_carry<_IntType, __m, __s, __r>& __x)
525 typedef std::basic_ostream<_CharT, _Traits> __istream_type;
526 typedef typename __istream_type::ios_base __ios_base;
528 const typename __ios_base::fmtflags __flags = __is.flags();
529 __is.flags(__ios_base::dec | __ios_base::skipws);
531 for (int __i = 0; __i < __r; ++__i)
532 __is >> __x._M_x[__i];
533 __is >> __x._M_carry;
540 template<typename _RealType, int __w, int __s, int __r>
542 subtract_with_carry_01<_RealType, __w, __s, __r>::word_size;
544 template<typename _RealType, int __w, int __s, int __r>
546 subtract_with_carry_01<_RealType, __w, __s, __r>::long_lag;
548 template<typename _RealType, int __w, int __s, int __r>
550 subtract_with_carry_01<_RealType, __w, __s, __r>::short_lag;
552 template<typename _RealType, int __w, int __s, int __r>
554 subtract_with_carry_01<_RealType, __w, __s, __r>::
555 _M_initialize_npows()
557 for (int __j = 0; __j < __n; ++__j)
558 #if _GLIBCXX_USE_C99_MATH_TR1
559 _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32);
561 _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
565 template<typename _RealType, int __w, int __s, int __r>
567 subtract_with_carry_01<_RealType, __w, __s, __r>::
568 seed(unsigned long __value)
573 // _GLIBCXX_RESOLVE_LIB_DEFECTS
574 // 512. Seeding subtract_with_carry_01 from a single unsigned long.
575 std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
581 template<typename _RealType, int __w, int __s, int __r>
584 subtract_with_carry_01<_RealType, __w, __s, __r>::
585 seed(_Gen& __gen, false_type)
587 for (int __i = 0; __i < long_lag; ++__i)
589 for (int __j = 0; __j < __n - 1; ++__j)
590 _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
591 _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
592 __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
596 for (int __j = 0; __j < __n; ++__j)
597 if (_M_x[long_lag - 1][__j] != 0)
606 template<typename _RealType, int __w, int __s, int __r>
607 typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
608 subtract_with_carry_01<_RealType, __w, __s, __r>::
611 // Derive short lag index from current index.
612 int __ps = _M_p - short_lag;
616 _UInt32Type __new_carry;
617 for (int __j = 0; __j < __n - 1; ++__j)
619 if (_M_x[__ps][__j] > _M_x[_M_p][__j]
620 || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
625 _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
626 _M_carry = __new_carry;
629 if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
630 || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
635 _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
636 __detail::_Shift<_UInt32Type, __w % 32>::__value>
637 (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
638 _M_carry = __new_carry;
640 result_type __ret = 0.0;
641 for (int __j = 0; __j < __n; ++__j)
642 __ret += _M_x[_M_p][__j] * _M_npows[__j];
644 // Adjust current index to loop around in ring buffer.
645 if (++_M_p >= long_lag)
651 template<typename _RealType, int __w, int __s, int __r,
652 typename _CharT, typename _Traits>
653 std::basic_ostream<_CharT, _Traits>&
654 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
655 const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
657 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
658 typedef typename __ostream_type::ios_base __ios_base;
660 const typename __ios_base::fmtflags __flags = __os.flags();
661 const _CharT __fill = __os.fill();
662 const _CharT __space = __os.widen(' ');
663 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
666 for (int __i = 0; __i < __r; ++__i)
667 for (int __j = 0; __j < __x.__n; ++__j)
668 __os << __x._M_x[__i][__j] << __space;
669 __os << __x._M_carry;
676 template<typename _RealType, int __w, int __s, int __r,
677 typename _CharT, typename _Traits>
678 std::basic_istream<_CharT, _Traits>&
679 operator>>(std::basic_istream<_CharT, _Traits>& __is,
680 subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
682 typedef std::basic_istream<_CharT, _Traits> __istream_type;
683 typedef typename __istream_type::ios_base __ios_base;
685 const typename __ios_base::fmtflags __flags = __is.flags();
686 __is.flags(__ios_base::dec | __ios_base::skipws);
688 for (int __i = 0; __i < __r; ++__i)
689 for (int __j = 0; __j < __x.__n; ++__j)
690 __is >> __x._M_x[__i][__j];
691 __is >> __x._M_carry;
697 template<class _UniformRandomNumberGenerator, int __p, int __r>
699 discard_block<_UniformRandomNumberGenerator, __p, __r>::block_size;
701 template<class _UniformRandomNumberGenerator, int __p, int __r>
703 discard_block<_UniformRandomNumberGenerator, __p, __r>::used_block;
705 template<class _UniformRandomNumberGenerator, int __p, int __r>
706 typename discard_block<_UniformRandomNumberGenerator,
707 __p, __r>::result_type
708 discard_block<_UniformRandomNumberGenerator, __p, __r>::
711 if (_M_n >= used_block)
713 while (_M_n < block_size)
724 template<class _UniformRandomNumberGenerator, int __p, int __r,
725 typename _CharT, typename _Traits>
726 std::basic_ostream<_CharT, _Traits>&
727 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
728 const discard_block<_UniformRandomNumberGenerator,
731 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
732 typedef typename __ostream_type::ios_base __ios_base;
734 const typename __ios_base::fmtflags __flags = __os.flags();
735 const _CharT __fill = __os.fill();
736 const _CharT __space = __os.widen(' ');
737 __os.flags(__ios_base::dec | __ios_base::fixed
741 __os << __x._M_b << __space << __x._M_n;
748 template<class _UniformRandomNumberGenerator, int __p, int __r,
749 typename _CharT, typename _Traits>
750 std::basic_istream<_CharT, _Traits>&
751 operator>>(std::basic_istream<_CharT, _Traits>& __is,
752 discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
754 typedef std::basic_istream<_CharT, _Traits> __istream_type;
755 typedef typename __istream_type::ios_base __ios_base;
757 const typename __ios_base::fmtflags __flags = __is.flags();
758 __is.flags(__ios_base::dec | __ios_base::skipws);
760 __is >> __x._M_b >> __x._M_n;
767 template<class _UniformRandomNumberGenerator1, int __s1,
768 class _UniformRandomNumberGenerator2, int __s2>
770 xor_combine<_UniformRandomNumberGenerator1, __s1,
771 _UniformRandomNumberGenerator2, __s2>::shift1;
773 template<class _UniformRandomNumberGenerator1, int __s1,
774 class _UniformRandomNumberGenerator2, int __s2>
776 xor_combine<_UniformRandomNumberGenerator1, __s1,
777 _UniformRandomNumberGenerator2, __s2>::shift2;
779 template<class _UniformRandomNumberGenerator1, int __s1,
780 class _UniformRandomNumberGenerator2, int __s2>
782 xor_combine<_UniformRandomNumberGenerator1, __s1,
783 _UniformRandomNumberGenerator2, __s2>::
786 const int __w = std::numeric_limits<result_type>::digits;
788 const result_type __m1 =
789 std::min(result_type(_M_b1.max() - _M_b1.min()),
790 __detail::_Shift<result_type, __w - __s1>::__value - 1);
792 const result_type __m2 =
793 std::min(result_type(_M_b2.max() - _M_b2.min()),
794 __detail::_Shift<result_type, __w - __s2>::__value - 1);
796 // NB: In TR1 s1 is not required to be >= s2.
798 _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
800 _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
803 template<class _UniformRandomNumberGenerator1, int __s1,
804 class _UniformRandomNumberGenerator2, int __s2>
805 typename xor_combine<_UniformRandomNumberGenerator1, __s1,
806 _UniformRandomNumberGenerator2, __s2>::result_type
807 xor_combine<_UniformRandomNumberGenerator1, __s1,
808 _UniformRandomNumberGenerator2, __s2>::
809 _M_initialize_max_aux(result_type __a, result_type __b, int __d)
811 const result_type __two2d = result_type(1) << __d;
812 const result_type __c = __a * __two2d;
814 if (__a == 0 || __b < __two2d)
817 const result_type __t = std::max(__c, __b);
818 const result_type __u = std::min(__c, __b);
820 result_type __ub = __u;
822 for (__p = 0; __ub != 1; __ub >>= 1)
825 const result_type __two2p = result_type(1) << __p;
826 const result_type __k = __t / __two2p;
829 return (__k + 1) * __two2p - 1;
832 return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
836 return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
841 template<class _UniformRandomNumberGenerator1, int __s1,
842 class _UniformRandomNumberGenerator2, int __s2,
843 typename _CharT, typename _Traits>
844 std::basic_ostream<_CharT, _Traits>&
845 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
846 const xor_combine<_UniformRandomNumberGenerator1, __s1,
847 _UniformRandomNumberGenerator2, __s2>& __x)
849 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
850 typedef typename __ostream_type::ios_base __ios_base;
852 const typename __ios_base::fmtflags __flags = __os.flags();
853 const _CharT __fill = __os.fill();
854 const _CharT __space = __os.widen(' ');
855 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
858 __os << __x.base1() << __space << __x.base2();
865 template<class _UniformRandomNumberGenerator1, int __s1,
866 class _UniformRandomNumberGenerator2, int __s2,
867 typename _CharT, typename _Traits>
868 std::basic_istream<_CharT, _Traits>&
869 operator>>(std::basic_istream<_CharT, _Traits>& __is,
870 xor_combine<_UniformRandomNumberGenerator1, __s1,
871 _UniformRandomNumberGenerator2, __s2>& __x)
873 typedef std::basic_istream<_CharT, _Traits> __istream_type;
874 typedef typename __istream_type::ios_base __ios_base;
876 const typename __ios_base::fmtflags __flags = __is.flags();
877 __is.flags(__ios_base::skipws);
879 __is >> __x._M_b1 >> __x._M_b2;
886 template<typename _IntType>
887 template<typename _UniformRandomNumberGenerator>
888 typename uniform_int<_IntType>::result_type
889 uniform_int<_IntType>::
890 _M_call(_UniformRandomNumberGenerator& __urng,
891 result_type __min, result_type __max, true_type)
893 // XXX Must be fixed to work well for *arbitrary* __urng.max(),
894 // __urng.min(), __max, __min. Currently works fine only in the
895 // most common case __urng.max() - __urng.min() >= __max - __min,
896 // with __urng.max() > __urng.min() >= 0.
897 typedef typename __gnu_cxx::__add_unsigned<typename
898 _UniformRandomNumberGenerator::result_type>::__type __urntype;
899 typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
901 typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
903 __urntype, __utype>::__type __uctype;
907 const __urntype __urnmin = __urng.min();
908 const __urntype __urnmax = __urng.max();
909 const __urntype __urnrange = __urnmax - __urnmin;
910 const __uctype __urange = __max - __min;
911 const __uctype __udenom = (__urnrange <= __urange
912 ? 1 : __urnrange / (__urange + 1));
914 __ret = (__urntype(__urng()) - __urnmin) / __udenom;
915 while (__ret > __max - __min);
917 return __ret + __min;
920 template<typename _IntType, typename _CharT, typename _Traits>
921 std::basic_ostream<_CharT, _Traits>&
922 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
923 const uniform_int<_IntType>& __x)
925 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
926 typedef typename __ostream_type::ios_base __ios_base;
928 const typename __ios_base::fmtflags __flags = __os.flags();
929 const _CharT __fill = __os.fill();
930 const _CharT __space = __os.widen(' ');
931 __os.flags(__ios_base::scientific | __ios_base::left);
934 __os << __x.min() << __space << __x.max();
941 template<typename _IntType, typename _CharT, typename _Traits>
942 std::basic_istream<_CharT, _Traits>&
943 operator>>(std::basic_istream<_CharT, _Traits>& __is,
944 uniform_int<_IntType>& __x)
946 typedef std::basic_istream<_CharT, _Traits> __istream_type;
947 typedef typename __istream_type::ios_base __ios_base;
949 const typename __ios_base::fmtflags __flags = __is.flags();
950 __is.flags(__ios_base::dec | __ios_base::skipws);
952 __is >> __x._M_min >> __x._M_max;
959 template<typename _CharT, typename _Traits>
960 std::basic_ostream<_CharT, _Traits>&
961 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
962 const bernoulli_distribution& __x)
964 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
965 typedef typename __ostream_type::ios_base __ios_base;
967 const typename __ios_base::fmtflags __flags = __os.flags();
968 const _CharT __fill = __os.fill();
969 const std::streamsize __precision = __os.precision();
970 __os.flags(__ios_base::scientific | __ios_base::left);
971 __os.fill(__os.widen(' '));
972 __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);
978 __os.precision(__precision);
983 template<typename _IntType, typename _RealType>
984 template<class _UniformRandomNumberGenerator>
985 typename geometric_distribution<_IntType, _RealType>::result_type
986 geometric_distribution<_IntType, _RealType>::
987 operator()(_UniformRandomNumberGenerator& __urng)
989 // About the epsilon thing see this thread:
990 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
991 const _RealType __naf =
992 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
993 // The largest _RealType convertible to _IntType.
994 const _RealType __thr =
995 std::numeric_limits<_IntType>::max() + __naf;
999 __cand = std::ceil(std::log(__urng()) / _M_log_p);
1000 while (__cand >= __thr);
1002 return result_type(__cand + __naf);
1005 template<typename _IntType, typename _RealType,
1006 typename _CharT, typename _Traits>
1007 std::basic_ostream<_CharT, _Traits>&
1008 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1009 const geometric_distribution<_IntType, _RealType>& __x)
1011 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1012 typedef typename __ostream_type::ios_base __ios_base;
1014 const typename __ios_base::fmtflags __flags = __os.flags();
1015 const _CharT __fill = __os.fill();
1016 const std::streamsize __precision = __os.precision();
1017 __os.flags(__ios_base::scientific | __ios_base::left);
1018 __os.fill(__os.widen(' '));
1019 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1023 __os.flags(__flags);
1025 __os.precision(__precision);
1030 template<typename _IntType, typename _RealType>
1032 poisson_distribution<_IntType, _RealType>::
1035 #if _GLIBCXX_USE_C99_MATH_TR1
1038 const _RealType __m = std::floor(_M_mean);
1039 _M_lm_thr = std::log(_M_mean);
1040 _M_lfm = std::tr1::lgamma(__m + 1);
1041 _M_sm = std::sqrt(__m);
1043 const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
1044 const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
1046 _M_d = std::tr1::round(std::max(_RealType(6),
1047 std::min(__m, __dx)));
1048 const _RealType __cx = 2 * __m + _M_d;
1049 _M_scx = std::sqrt(__cx / 2);
1052 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1053 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
1057 _M_lm_thr = std::exp(-_M_mean);
1061 * A rejection algorithm when mean >= 12 and a simple method based
1062 * upon the multiplication of uniform random variates otherwise.
1063 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1067 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1068 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1070 template<typename _IntType, typename _RealType>
1071 template<class _UniformRandomNumberGenerator>
1072 typename poisson_distribution<_IntType, _RealType>::result_type
1073 poisson_distribution<_IntType, _RealType>::
1074 operator()(_UniformRandomNumberGenerator& __urng)
1076 #if _GLIBCXX_USE_C99_MATH_TR1
1081 // See comments above...
1082 const _RealType __naf =
1083 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
1084 const _RealType __thr =
1085 std::numeric_limits<_IntType>::max() + __naf;
1087 const _RealType __m = std::floor(_M_mean);
1089 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1090 const _RealType __c1 = _M_sm * __spi_2;
1091 const _RealType __c2 = _M_c2b + __c1;
1092 const _RealType __c3 = __c2 + 1;
1093 const _RealType __c4 = __c3 + 1;
1095 const _RealType __e178 = 1.0129030479320018583185514777512983L;
1096 const _RealType __c5 = __c4 + __e178;
1097 const _RealType __c = _M_cb + __c5;
1098 const _RealType __2cx = 2 * (2 * __m + _M_d);
1100 bool __reject = true;
1103 const _RealType __u = __c * __urng();
1104 const _RealType __e = -std::log(__urng());
1106 _RealType __w = 0.0;
1110 const _RealType __n = _M_nd(__urng);
1111 const _RealType __y = -std::abs(__n) * _M_sm - 1;
1112 __x = std::floor(__y);
1113 __w = -__n * __n / 2;
1117 else if (__u <= __c2)
1119 const _RealType __n = _M_nd(__urng);
1120 const _RealType __y = 1 + std::abs(__n) * _M_scx;
1121 __x = std::ceil(__y);
1122 __w = __y * (2 - __y) * _M_1cx;
1126 else if (__u <= __c3)
1127 // NB: This case not in the book, nor in the Errata,
1128 // but should be ok...
1130 else if (__u <= __c4)
1132 else if (__u <= __c5)
1136 const _RealType __v = -std::log(__urng());
1137 const _RealType __y = _M_d + __v * __2cx / _M_d;
1138 __x = std::ceil(__y);
1139 __w = -_M_d * _M_1cx * (1 + __y / 2);
1142 __reject = (__w - __e - __x * _M_lm_thr
1143 > _M_lfm - std::tr1::lgamma(__x + __m + 1));
1145 __reject |= __x + __m >= __thr;
1149 return result_type(__x + __m + __naf);
1155 _RealType __prod = 1.0;
1162 while (__prod > _M_lm_thr);
1168 template<typename _IntType, typename _RealType,
1169 typename _CharT, typename _Traits>
1170 std::basic_ostream<_CharT, _Traits>&
1171 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1172 const poisson_distribution<_IntType, _RealType>& __x)
1174 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1175 typedef typename __ostream_type::ios_base __ios_base;
1177 const typename __ios_base::fmtflags __flags = __os.flags();
1178 const _CharT __fill = __os.fill();
1179 const std::streamsize __precision = __os.precision();
1180 const _CharT __space = __os.widen(' ');
1181 __os.flags(__ios_base::scientific | __ios_base::left);
1183 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1185 __os << __x.mean() << __space << __x._M_nd;
1187 __os.flags(__flags);
1189 __os.precision(__precision);
1193 template<typename _IntType, typename _RealType,
1194 typename _CharT, typename _Traits>
1195 std::basic_istream<_CharT, _Traits>&
1196 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1197 poisson_distribution<_IntType, _RealType>& __x)
1199 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1200 typedef typename __istream_type::ios_base __ios_base;
1202 const typename __ios_base::fmtflags __flags = __is.flags();
1203 __is.flags(__ios_base::skipws);
1205 __is >> __x._M_mean >> __x._M_nd;
1206 __x._M_initialize();
1208 __is.flags(__flags);
1213 template<typename _IntType, typename _RealType>
1215 binomial_distribution<_IntType, _RealType>::
1218 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1222 #if _GLIBCXX_USE_C99_MATH_TR1
1223 if (_M_t * __p12 >= 8)
1226 const _RealType __np = std::floor(_M_t * __p12);
1227 const _RealType __pa = __np / _M_t;
1228 const _RealType __1p = 1 - __pa;
1230 const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
1231 const _RealType __d1x =
1232 std::sqrt(__np * __1p * std::log(32 * __np
1233 / (81 * __pi_4 * __1p)));
1234 _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));
1235 const _RealType __d2x =
1236 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1237 / (__pi_4 * __pa)));
1238 _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));
1241 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1242 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1243 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1244 _M_c = 2 * _M_d1 / __np;
1245 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1246 const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
1247 const _RealType __s1s = _M_s1 * _M_s1;
1248 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1250 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1251 const _RealType __s2s = _M_s2 * _M_s2;
1252 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1253 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1254 _M_lf = (std::tr1::lgamma(__np + 1)
1255 + std::tr1::lgamma(_M_t - __np + 1));
1256 _M_lp1p = std::log(__pa / __1p);
1258 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1262 _M_q = -std::log(1 - __p12);
1265 template<typename _IntType, typename _RealType>
1266 template<class _UniformRandomNumberGenerator>
1267 typename binomial_distribution<_IntType, _RealType>::result_type
1268 binomial_distribution<_IntType, _RealType>::
1269 _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1272 _RealType __sum = 0;
1276 const _RealType __e = -std::log(__urng());
1277 __sum += __e / (__t - __x);
1280 while (__sum <= _M_q);
1286 * A rejection algorithm when t * p >= 8 and a simple waiting time
1287 * method - the second in the referenced book - otherwise.
1288 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1292 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1293 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1295 template<typename _IntType, typename _RealType>
1296 template<class _UniformRandomNumberGenerator>
1297 typename binomial_distribution<_IntType, _RealType>::result_type
1298 binomial_distribution<_IntType, _RealType>::
1299 operator()(_UniformRandomNumberGenerator& __urng)
1302 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1304 #if _GLIBCXX_USE_C99_MATH_TR1
1309 // See comments above...
1310 const _RealType __naf =
1311 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
1312 const _RealType __thr =
1313 std::numeric_limits<_IntType>::max() + __naf;
1315 const _RealType __np = std::floor(_M_t * __p12);
1316 const _RealType __pa = __np / _M_t;
1319 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1320 const _RealType __a1 = _M_a1;
1321 const _RealType __a12 = __a1 + _M_s2 * __spi_2;
1322 const _RealType __a123 = _M_a123;
1323 const _RealType __s1s = _M_s1 * _M_s1;
1324 const _RealType __s2s = _M_s2 * _M_s2;
1329 const _RealType __u = _M_s * __urng();
1335 const _RealType __n = _M_nd(__urng);
1336 const _RealType __y = _M_s1 * std::abs(__n);
1337 __reject = __y >= _M_d1;
1340 const _RealType __e = -std::log(__urng());
1341 __x = std::floor(__y);
1342 __v = -__e - __n * __n / 2 + _M_c;
1345 else if (__u <= __a12)
1347 const _RealType __n = _M_nd(__urng);
1348 const _RealType __y = _M_s2 * std::abs(__n);
1349 __reject = __y >= _M_d2;
1352 const _RealType __e = -std::log(__urng());
1353 __x = std::floor(-__y);
1354 __v = -__e - __n * __n / 2;
1357 else if (__u <= __a123)
1359 const _RealType __e1 = -std::log(__urng());
1360 const _RealType __e2 = -std::log(__urng());
1362 const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
1363 __x = std::floor(__y);
1364 __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
1365 -__y / (2 * __s1s)));
1370 const _RealType __e1 = -std::log(__urng());
1371 const _RealType __e2 = -std::log(__urng());
1373 const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
1374 __x = std::floor(-__y);
1375 __v = -__e2 - _M_d2 * __y / (2 * __s2s);
1379 __reject = __reject || __x < -__np || __x > _M_t - __np;
1382 const _RealType __lfx =
1383 std::tr1::lgamma(__np + __x + 1)
1384 + std::tr1::lgamma(_M_t - (__np + __x) + 1);
1385 __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
1388 __reject |= __x + __np >= __thr;
1392 __x += __np + __naf;
1394 const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x));
1395 __ret = _IntType(__x) + __z;
1399 __ret = _M_waiting(__urng, _M_t);
1402 __ret = _M_t - __ret;
1406 template<typename _IntType, typename _RealType,
1407 typename _CharT, typename _Traits>
1408 std::basic_ostream<_CharT, _Traits>&
1409 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1410 const binomial_distribution<_IntType, _RealType>& __x)
1412 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1413 typedef typename __ostream_type::ios_base __ios_base;
1415 const typename __ios_base::fmtflags __flags = __os.flags();
1416 const _CharT __fill = __os.fill();
1417 const std::streamsize __precision = __os.precision();
1418 const _CharT __space = __os.widen(' ');
1419 __os.flags(__ios_base::scientific | __ios_base::left);
1421 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1423 __os << __x.t() << __space << __x.p()
1424 << __space << __x._M_nd;
1426 __os.flags(__flags);
1428 __os.precision(__precision);
1432 template<typename _IntType, typename _RealType,
1433 typename _CharT, typename _Traits>
1434 std::basic_istream<_CharT, _Traits>&
1435 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1436 binomial_distribution<_IntType, _RealType>& __x)
1438 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1439 typedef typename __istream_type::ios_base __ios_base;
1441 const typename __ios_base::fmtflags __flags = __is.flags();
1442 __is.flags(__ios_base::dec | __ios_base::skipws);
1444 __is >> __x._M_t >> __x._M_p >> __x._M_nd;
1445 __x._M_initialize();
1447 __is.flags(__flags);
1452 template<typename _RealType, typename _CharT, typename _Traits>
1453 std::basic_ostream<_CharT, _Traits>&
1454 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1455 const uniform_real<_RealType>& __x)
1457 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1458 typedef typename __ostream_type::ios_base __ios_base;
1460 const typename __ios_base::fmtflags __flags = __os.flags();
1461 const _CharT __fill = __os.fill();
1462 const std::streamsize __precision = __os.precision();
1463 const _CharT __space = __os.widen(' ');
1464 __os.flags(__ios_base::scientific | __ios_base::left);
1466 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1468 __os << __x.min() << __space << __x.max();
1470 __os.flags(__flags);
1472 __os.precision(__precision);
1476 template<typename _RealType, typename _CharT, typename _Traits>
1477 std::basic_istream<_CharT, _Traits>&
1478 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1479 uniform_real<_RealType>& __x)
1481 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1482 typedef typename __istream_type::ios_base __ios_base;
1484 const typename __ios_base::fmtflags __flags = __is.flags();
1485 __is.flags(__ios_base::skipws);
1487 __is >> __x._M_min >> __x._M_max;
1489 __is.flags(__flags);
1494 template<typename _RealType, typename _CharT, typename _Traits>
1495 std::basic_ostream<_CharT, _Traits>&
1496 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1497 const exponential_distribution<_RealType>& __x)
1499 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1500 typedef typename __ostream_type::ios_base __ios_base;
1502 const typename __ios_base::fmtflags __flags = __os.flags();
1503 const _CharT __fill = __os.fill();
1504 const std::streamsize __precision = __os.precision();
1505 __os.flags(__ios_base::scientific | __ios_base::left);
1506 __os.fill(__os.widen(' '));
1507 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1509 __os << __x.lambda();
1511 __os.flags(__flags);
1513 __os.precision(__precision);
1519 * Polar method due to Marsaglia.
1521 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1522 * New York, 1986, Ch. V, Sect. 4.4.
1524 template<typename _RealType>
1525 template<class _UniformRandomNumberGenerator>
1526 typename normal_distribution<_RealType>::result_type
1527 normal_distribution<_RealType>::
1528 operator()(_UniformRandomNumberGenerator& __urng)
1532 if (_M_saved_available)
1534 _M_saved_available = false;
1539 result_type __x, __y, __r2;
1542 __x = result_type(2.0) * __urng() - 1.0;
1543 __y = result_type(2.0) * __urng() - 1.0;
1544 __r2 = __x * __x + __y * __y;
1546 while (__r2 > 1.0 || __r2 == 0.0);
1548 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1549 _M_saved = __x * __mult;
1550 _M_saved_available = true;
1551 __ret = __y * __mult;
1554 __ret = __ret * _M_sigma + _M_mean;
1558 template<typename _RealType, typename _CharT, typename _Traits>
1559 std::basic_ostream<_CharT, _Traits>&
1560 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1561 const normal_distribution<_RealType>& __x)
1563 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1564 typedef typename __ostream_type::ios_base __ios_base;
1566 const typename __ios_base::fmtflags __flags = __os.flags();
1567 const _CharT __fill = __os.fill();
1568 const std::streamsize __precision = __os.precision();
1569 const _CharT __space = __os.widen(' ');
1570 __os.flags(__ios_base::scientific | __ios_base::left);
1572 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1574 __os << __x._M_saved_available << __space
1575 << __x.mean() << __space
1577 if (__x._M_saved_available)
1578 __os << __space << __x._M_saved;
1580 __os.flags(__flags);
1582 __os.precision(__precision);
1586 template<typename _RealType, typename _CharT, typename _Traits>
1587 std::basic_istream<_CharT, _Traits>&
1588 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1589 normal_distribution<_RealType>& __x)
1591 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1592 typedef typename __istream_type::ios_base __ios_base;
1594 const typename __ios_base::fmtflags __flags = __is.flags();
1595 __is.flags(__ios_base::dec | __ios_base::skipws);
1597 __is >> __x._M_saved_available >> __x._M_mean
1599 if (__x._M_saved_available)
1600 __is >> __x._M_saved;
1602 __is.flags(__flags);
1607 template<typename _RealType>
1609 gamma_distribution<_RealType>::
1613 _M_l_d = std::sqrt(2 * _M_alpha - 1);
1615 _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
1620 * Cheng's rejection algorithm GB for alpha >= 1 and a modification
1621 * of Vaduva's rejection from Weibull algorithm due to Devroye for
1625 * Cheng, R. C. The Generation of Gamma Random Variables with Non-integral
1626 * Shape Parameter. Applied Statistics, 26, 71-75, 1977.
1628 * Vaduva, I. Computer Generation of Gamma Gandom Variables by Rejection
1629 * and Composition Procedures. Math. Operationsforschung and Statistik,
1630 * Series in Statistics, 8, 545-576, 1977.
1632 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1633 * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
1635 template<typename _RealType>
1636 template<class _UniformRandomNumberGenerator>
1637 typename gamma_distribution<_RealType>::result_type
1638 gamma_distribution<_RealType>::
1639 operator()(_UniformRandomNumberGenerator& __urng)
1647 const result_type __b = _M_alpha
1648 - result_type(1.3862943611198906188344642429163531L);
1649 const result_type __c = _M_alpha + _M_l_d;
1650 const result_type __1l = 1 / _M_l_d;
1653 const result_type __k = 2.5040773967762740733732583523868748L;
1657 const result_type __u = __urng();
1658 const result_type __v = __urng();
1660 const result_type __y = __1l * std::log(__v / (1 - __v));
1661 __x = _M_alpha * std::exp(__y);
1663 const result_type __z = __u * __v * __v;
1664 const result_type __r = __b + __c * __y - __x;
1666 __reject = __r < result_type(4.5) * __z - __k;
1668 __reject = __r < std::log(__z);
1674 const result_type __c = 1 / _M_alpha;
1678 const result_type __z = -std::log(__urng());
1679 const result_type __e = -std::log(__urng());
1681 __x = std::pow(__z, __c);
1683 __reject = __z + __e < _M_l_d + __x;
1691 template<typename _RealType, typename _CharT, typename _Traits>
1692 std::basic_ostream<_CharT, _Traits>&
1693 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1694 const gamma_distribution<_RealType>& __x)
1696 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1697 typedef typename __ostream_type::ios_base __ios_base;
1699 const typename __ios_base::fmtflags __flags = __os.flags();
1700 const _CharT __fill = __os.fill();
1701 const std::streamsize __precision = __os.precision();
1702 __os.flags(__ios_base::scientific | __ios_base::left);
1703 __os.fill(__os.widen(' '));
1704 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1706 __os << __x.alpha();
1708 __os.flags(__flags);
1710 __os.precision(__precision);
1715 _GLIBCXX_END_NAMESPACE_VERSION