1 // random number generation (out of line) -*- C++ -*-
3 // Copyright (C) 2007 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 tr1_impl/random.tcc
31 * This is an internal header file, included by other library headers.
32 * You should not attempt to use it directly.
37 _GLIBCXX_BEGIN_NAMESPACE_TR1
40 * (Further) implementation-space details.
44 // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
47 // Because a and c are compile-time integral constants the compiler kindly
48 // elides any unreachable paths.
50 // Preconditions: a > 0, m > 0.
52 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
62 static const _Tp __q = __m / __a;
63 static const _Tp __r = __m % __a;
65 _Tp __t1 = __a * (__x % __q);
66 _Tp __t2 = __r * (__x / __q);
70 __x = __m - __t2 + __t1;
75 const _Tp __d = __m - __x;
85 // Special case for m == 0 -- use unsigned integer overflow as modulo
87 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
88 struct _Mod<_Tp, __a, __c, __m, true>
92 { return __a * __x + __c; }
94 } // namespace __detail
97 * Seeds the LCR with integral value @p __x0, adjusted so that the
98 * ring identity is never a member of the convergence set.
100 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
102 linear_congruential<_UIntType, __a, __c, __m>::
103 seed(unsigned long __x0)
105 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
106 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
107 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
109 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
113 * Seeds the LCR engine with a value generated by @p __g.
115 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
118 linear_congruential<_UIntType, __a, __c, __m>::
119 seed(_Gen& __g, false_type)
121 _UIntType __x0 = __g();
122 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
123 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
124 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
126 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
130 * Gets the next generated value in sequence.
132 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
133 typename linear_congruential<_UIntType, __a, __c, __m>::result_type
134 linear_congruential<_UIntType, __a, __c, __m>::
137 _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
141 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
142 typename _CharT, typename _Traits>
143 std::basic_ostream<_CharT, _Traits>&
144 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
145 const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
147 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
148 typedef typename __ostream_type::ios_base __ios_base;
150 const typename __ios_base::fmtflags __flags = __os.flags();
151 const _CharT __fill = __os.fill();
152 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
153 __os.fill(__os.widen(' '));
162 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
163 typename _CharT, typename _Traits>
164 std::basic_istream<_CharT, _Traits>&
165 operator>>(std::basic_istream<_CharT, _Traits>& __is,
166 linear_congruential<_UIntType, __a, __c, __m>& __lcr)
168 typedef std::basic_istream<_CharT, _Traits> __istream_type;
169 typedef typename __istream_type::ios_base __ios_base;
171 const typename __ios_base::fmtflags __flags = __is.flags();
172 __is.flags(__ios_base::dec);
181 template<class _UIntType, int __w, int __n, int __m, int __r,
182 _UIntType __a, int __u, int __s,
183 _UIntType __b, int __t, _UIntType __c, int __l>
185 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
186 __b, __t, __c, __l>::
187 seed(unsigned long __value)
189 _M_x[0] = __detail::__mod<_UIntType, 1, 0,
190 __detail::_Shift<_UIntType, __w>::__value>(__value);
192 for (int __i = 1; __i < state_size; ++__i)
194 _UIntType __x = _M_x[__i - 1];
195 __x ^= __x >> (__w - 2);
198 _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
199 __detail::_Shift<_UIntType, __w>::__value>(__x);
204 template<class _UIntType, int __w, int __n, int __m, int __r,
205 _UIntType __a, int __u, int __s,
206 _UIntType __b, int __t, _UIntType __c, int __l>
209 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
210 __b, __t, __c, __l>::
211 seed(_Gen& __gen, false_type)
213 for (int __i = 0; __i < state_size; ++__i)
214 _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
215 __detail::_Shift<_UIntType, __w>::__value>(__gen());
219 template<class _UIntType, int __w, int __n, int __m, int __r,
220 _UIntType __a, int __u, int __s,
221 _UIntType __b, int __t, _UIntType __c, int __l>
223 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
224 __b, __t, __c, __l>::result_type
225 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
226 __b, __t, __c, __l>::
229 // Reload the vector - cost is O(n) amortized over n calls.
230 if (_M_p >= state_size)
232 const _UIntType __upper_mask = (~_UIntType()) << __r;
233 const _UIntType __lower_mask = ~__upper_mask;
235 for (int __k = 0; __k < (__n - __m); ++__k)
237 _UIntType __y = ((_M_x[__k] & __upper_mask)
238 | (_M_x[__k + 1] & __lower_mask));
239 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
240 ^ ((__y & 0x01) ? __a : 0));
243 for (int __k = (__n - __m); __k < (__n - 1); ++__k)
245 _UIntType __y = ((_M_x[__k] & __upper_mask)
246 | (_M_x[__k + 1] & __lower_mask));
247 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
248 ^ ((__y & 0x01) ? __a : 0));
251 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
252 | (_M_x[0] & __lower_mask));
253 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
254 ^ ((__y & 0x01) ? __a : 0));
258 // Calculate o(x(i)).
259 result_type __z = _M_x[_M_p++];
261 __z ^= (__z << __s) & __b;
262 __z ^= (__z << __t) & __c;
268 template<class _UIntType, int __w, int __n, int __m, int __r,
269 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
270 _UIntType __c, int __l,
271 typename _CharT, typename _Traits>
272 std::basic_ostream<_CharT, _Traits>&
273 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
274 const mersenne_twister<_UIntType, __w, __n, __m,
275 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
277 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
278 typedef typename __ostream_type::ios_base __ios_base;
280 const typename __ios_base::fmtflags __flags = __os.flags();
281 const _CharT __fill = __os.fill();
282 const _CharT __space = __os.widen(' ');
283 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
286 for (int __i = 0; __i < __n - 1; ++__i)
287 __os << __x._M_x[__i] << __space;
288 __os << __x._M_x[__n - 1];
295 template<class _UIntType, int __w, int __n, int __m, int __r,
296 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
297 _UIntType __c, int __l,
298 typename _CharT, typename _Traits>
299 std::basic_istream<_CharT, _Traits>&
300 operator>>(std::basic_istream<_CharT, _Traits>& __is,
301 mersenne_twister<_UIntType, __w, __n, __m,
302 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
304 typedef std::basic_istream<_CharT, _Traits> __istream_type;
305 typedef typename __istream_type::ios_base __ios_base;
307 const typename __ios_base::fmtflags __flags = __is.flags();
308 __is.flags(__ios_base::dec | __ios_base::skipws);
310 for (int __i = 0; __i < __n; ++__i)
311 __is >> __x._M_x[__i];
318 template<typename _IntType, _IntType __m, int __s, int __r>
320 subtract_with_carry<_IntType, __m, __s, __r>::
321 seed(unsigned long __value)
326 std::_GLIBCXX_TR1 linear_congruential<unsigned long, 40014, 0, 2147483563>
329 for (int __i = 0; __i < long_lag; ++__i)
330 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());
332 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
336 template<typename _IntType, _IntType __m, int __s, int __r>
339 subtract_with_carry<_IntType, __m, __s, __r>::
340 seed(_Gen& __gen, false_type)
342 const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;
344 for (int __i = 0; __i < long_lag; ++__i)
347 _UIntType __factor = 1;
348 for (int __j = 0; __j < __n; ++__j)
350 __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
351 (__gen()) * __factor;
352 __factor *= __detail::_Shift<_UIntType, 32>::__value;
354 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
356 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
360 template<typename _IntType, _IntType __m, int __s, int __r>
361 typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
362 subtract_with_carry<_IntType, __m, __s, __r>::
365 // Derive short lag index from current index.
366 int __ps = _M_p - short_lag;
370 // Calculate new x(i) without overflow or division.
371 // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
374 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
376 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
381 __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
386 // Adjust current index to loop around in ring buffer.
387 if (++_M_p >= long_lag)
393 template<typename _IntType, _IntType __m, int __s, int __r,
394 typename _CharT, typename _Traits>
395 std::basic_ostream<_CharT, _Traits>&
396 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
397 const subtract_with_carry<_IntType, __m, __s, __r>& __x)
399 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
400 typedef typename __ostream_type::ios_base __ios_base;
402 const typename __ios_base::fmtflags __flags = __os.flags();
403 const _CharT __fill = __os.fill();
404 const _CharT __space = __os.widen(' ');
405 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
408 for (int __i = 0; __i < __r; ++__i)
409 __os << __x._M_x[__i] << __space;
410 __os << __x._M_carry;
417 template<typename _IntType, _IntType __m, int __s, int __r,
418 typename _CharT, typename _Traits>
419 std::basic_istream<_CharT, _Traits>&
420 operator>>(std::basic_istream<_CharT, _Traits>& __is,
421 subtract_with_carry<_IntType, __m, __s, __r>& __x)
423 typedef std::basic_ostream<_CharT, _Traits> __istream_type;
424 typedef typename __istream_type::ios_base __ios_base;
426 const typename __ios_base::fmtflags __flags = __is.flags();
427 __is.flags(__ios_base::dec | __ios_base::skipws);
429 for (int __i = 0; __i < __r; ++__i)
430 __is >> __x._M_x[__i];
431 __is >> __x._M_carry;
438 template<typename _RealType, int __w, int __s, int __r>
440 subtract_with_carry_01<_RealType, __w, __s, __r>::
441 _M_initialize_npows()
443 for (int __j = 0; __j < __n; ++__j)
444 #if _GLIBCXX_USE_C99_MATH_TR1
445 _M_npows[__j] = std::_GLIBCXX_TR1 ldexp(_RealType(1), -__w + __j * 32);
447 _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
451 template<typename _RealType, int __w, int __s, int __r>
453 subtract_with_carry_01<_RealType, __w, __s, __r>::
454 seed(unsigned long __value)
459 // _GLIBCXX_RESOLVE_LIB_DEFECTS
460 // 512. Seeding subtract_with_carry_01 from a single unsigned long.
461 std::_GLIBCXX_TR1 linear_congruential<unsigned long, 40014, 0, 2147483563>
467 template<typename _RealType, int __w, int __s, int __r>
470 subtract_with_carry_01<_RealType, __w, __s, __r>::
471 seed(_Gen& __gen, false_type)
473 for (int __i = 0; __i < long_lag; ++__i)
475 for (int __j = 0; __j < __n - 1; ++__j)
476 _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
477 _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
478 __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
482 for (int __j = 0; __j < __n; ++__j)
483 if (_M_x[long_lag - 1][__j] != 0)
492 template<typename _RealType, int __w, int __s, int __r>
493 typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
494 subtract_with_carry_01<_RealType, __w, __s, __r>::
497 // Derive short lag index from current index.
498 int __ps = _M_p - short_lag;
502 _UInt32Type __new_carry;
503 for (int __j = 0; __j < __n - 1; ++__j)
505 if (_M_x[__ps][__j] > _M_x[_M_p][__j]
506 || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
511 _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
512 _M_carry = __new_carry;
515 if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
516 || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
521 _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
522 __detail::_Shift<_UInt32Type, __w % 32>::__value>
523 (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
524 _M_carry = __new_carry;
526 result_type __ret = 0.0;
527 for (int __j = 0; __j < __n; ++__j)
528 __ret += _M_x[_M_p][__j] * _M_npows[__j];
530 // Adjust current index to loop around in ring buffer.
531 if (++_M_p >= long_lag)
537 template<typename _RealType, int __w, int __s, int __r,
538 typename _CharT, typename _Traits>
539 std::basic_ostream<_CharT, _Traits>&
540 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
541 const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
543 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
544 typedef typename __ostream_type::ios_base __ios_base;
546 const typename __ios_base::fmtflags __flags = __os.flags();
547 const _CharT __fill = __os.fill();
548 const _CharT __space = __os.widen(' ');
549 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
552 for (int __i = 0; __i < __r; ++__i)
553 for (int __j = 0; __j < __x.__n; ++__j)
554 __os << __x._M_x[__i][__j] << __space;
555 __os << __x._M_carry;
562 template<typename _RealType, int __w, int __s, int __r,
563 typename _CharT, typename _Traits>
564 std::basic_istream<_CharT, _Traits>&
565 operator>>(std::basic_istream<_CharT, _Traits>& __is,
566 subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
568 typedef std::basic_istream<_CharT, _Traits> __istream_type;
569 typedef typename __istream_type::ios_base __ios_base;
571 const typename __ios_base::fmtflags __flags = __is.flags();
572 __is.flags(__ios_base::dec | __ios_base::skipws);
574 for (int __i = 0; __i < __r; ++__i)
575 for (int __j = 0; __j < __x.__n; ++__j)
576 __is >> __x._M_x[__i][__j];
577 __is >> __x._M_carry;
584 template<class _UniformRandomNumberGenerator, int __p, int __r>
585 typename discard_block<_UniformRandomNumberGenerator,
586 __p, __r>::result_type
587 discard_block<_UniformRandomNumberGenerator, __p, __r>::
590 if (_M_n >= used_block)
592 while (_M_n < block_size)
603 template<class _UniformRandomNumberGenerator, int __p, int __r,
604 typename _CharT, typename _Traits>
605 std::basic_ostream<_CharT, _Traits>&
606 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
607 const discard_block<_UniformRandomNumberGenerator,
610 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
611 typedef typename __ostream_type::ios_base __ios_base;
613 const typename __ios_base::fmtflags __flags = __os.flags();
614 const _CharT __fill = __os.fill();
615 const _CharT __space = __os.widen(' ');
616 __os.flags(__ios_base::dec | __ios_base::fixed
620 __os << __x._M_b << __space << __x._M_n;
627 template<class _UniformRandomNumberGenerator, int __p, int __r,
628 typename _CharT, typename _Traits>
629 std::basic_istream<_CharT, _Traits>&
630 operator>>(std::basic_istream<_CharT, _Traits>& __is,
631 discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
633 typedef std::basic_istream<_CharT, _Traits> __istream_type;
634 typedef typename __istream_type::ios_base __ios_base;
636 const typename __ios_base::fmtflags __flags = __is.flags();
637 __is.flags(__ios_base::dec | __ios_base::skipws);
639 __is >> __x._M_b >> __x._M_n;
646 template<class _UniformRandomNumberGenerator1, int __s1,
647 class _UniformRandomNumberGenerator2, int __s2>
649 xor_combine<_UniformRandomNumberGenerator1, __s1,
650 _UniformRandomNumberGenerator2, __s2>::
653 const int __w = std::numeric_limits<result_type>::digits;
655 const result_type __m1 =
656 std::min(result_type(_M_b1.max() - _M_b1.min()),
657 __detail::_Shift<result_type, __w - __s1>::__value - 1);
659 const result_type __m2 =
660 std::min(result_type(_M_b2.max() - _M_b2.min()),
661 __detail::_Shift<result_type, __w - __s2>::__value - 1);
663 // NB: In TR1 s1 is not required to be >= s2.
665 _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
667 _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
670 template<class _UniformRandomNumberGenerator1, int __s1,
671 class _UniformRandomNumberGenerator2, int __s2>
672 typename xor_combine<_UniformRandomNumberGenerator1, __s1,
673 _UniformRandomNumberGenerator2, __s2>::result_type
674 xor_combine<_UniformRandomNumberGenerator1, __s1,
675 _UniformRandomNumberGenerator2, __s2>::
676 _M_initialize_max_aux(result_type __a, result_type __b, int __d)
678 const result_type __two2d = result_type(1) << __d;
679 const result_type __c = __a * __two2d;
681 if (__a == 0 || __b < __two2d)
684 const result_type __t = std::max(__c, __b);
685 const result_type __u = std::min(__c, __b);
687 result_type __ub = __u;
689 for (__p = 0; __ub != 1; __ub >>= 1)
692 const result_type __two2p = result_type(1) << __p;
693 const result_type __k = __t / __two2p;
696 return (__k + 1) * __two2p - 1;
699 return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
703 return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
708 template<class _UniformRandomNumberGenerator1, int __s1,
709 class _UniformRandomNumberGenerator2, int __s2,
710 typename _CharT, typename _Traits>
711 std::basic_ostream<_CharT, _Traits>&
712 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
713 const xor_combine<_UniformRandomNumberGenerator1, __s1,
714 _UniformRandomNumberGenerator2, __s2>& __x)
716 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
717 typedef typename __ostream_type::ios_base __ios_base;
719 const typename __ios_base::fmtflags __flags = __os.flags();
720 const _CharT __fill = __os.fill();
721 const _CharT __space = __os.widen(' ');
722 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
725 __os << __x.base1() << __space << __x.base2();
732 template<class _UniformRandomNumberGenerator1, int __s1,
733 class _UniformRandomNumberGenerator2, int __s2,
734 typename _CharT, typename _Traits>
735 std::basic_istream<_CharT, _Traits>&
736 operator>>(std::basic_istream<_CharT, _Traits>& __is,
737 xor_combine<_UniformRandomNumberGenerator1, __s1,
738 _UniformRandomNumberGenerator2, __s2>& __x)
740 typedef std::basic_istream<_CharT, _Traits> __istream_type;
741 typedef typename __istream_type::ios_base __ios_base;
743 const typename __ios_base::fmtflags __flags = __is.flags();
744 __is.flags(__ios_base::skipws);
746 __is >> __x._M_b1 >> __x._M_b2;
753 template<typename _IntType>
754 template<typename _UniformRandomNumberGenerator>
755 typename uniform_int<_IntType>::result_type
756 uniform_int<_IntType>::
757 _M_call(_UniformRandomNumberGenerator& __urng,
758 result_type __min, result_type __max, true_type)
760 // XXX Must be fixed to work well for *arbitrary* __urng.max(),
761 // __urng.min(), __max, __min. Currently works fine only in the
762 // most common case __urng.max() - __urng.min() >= __max - __min,
763 // with __urng.max() > __urng.min() >= 0.
764 typedef typename __gnu_cxx::__add_unsigned<typename
765 _UniformRandomNumberGenerator::result_type>::__type __urntype;
766 typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
768 typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
770 __urntype, __utype>::__type __uctype;
774 const __urntype __urnmin = __urng.min();
775 const __urntype __urnmax = __urng.max();
776 const __urntype __urnrange = __urnmax - __urnmin;
777 const __uctype __urange = __max - __min;
778 const __uctype __udenom = (__urnrange <= __urange
779 ? 1 : __urnrange / (__urange + 1));
781 __ret = (__urntype(__urng()) - __urnmin) / __udenom;
782 while (__ret > __max - __min);
784 return __ret + __min;
787 template<typename _IntType, typename _CharT, typename _Traits>
788 std::basic_ostream<_CharT, _Traits>&
789 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
790 const uniform_int<_IntType>& __x)
792 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
793 typedef typename __ostream_type::ios_base __ios_base;
795 const typename __ios_base::fmtflags __flags = __os.flags();
796 const _CharT __fill = __os.fill();
797 const _CharT __space = __os.widen(' ');
798 __os.flags(__ios_base::scientific | __ios_base::left);
801 __os << __x.min() << __space << __x.max();
808 template<typename _IntType, typename _CharT, typename _Traits>
809 std::basic_istream<_CharT, _Traits>&
810 operator>>(std::basic_istream<_CharT, _Traits>& __is,
811 uniform_int<_IntType>& __x)
813 typedef std::basic_istream<_CharT, _Traits> __istream_type;
814 typedef typename __istream_type::ios_base __ios_base;
816 const typename __ios_base::fmtflags __flags = __is.flags();
817 __is.flags(__ios_base::dec | __ios_base::skipws);
819 __is >> __x._M_min >> __x._M_max;
826 template<typename _CharT, typename _Traits>
827 std::basic_ostream<_CharT, _Traits>&
828 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
829 const bernoulli_distribution& __x)
831 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
832 typedef typename __ostream_type::ios_base __ios_base;
834 const typename __ios_base::fmtflags __flags = __os.flags();
835 const _CharT __fill = __os.fill();
836 const std::streamsize __precision = __os.precision();
837 __os.flags(__ios_base::scientific | __ios_base::left);
838 __os.fill(__os.widen(' '));
839 __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);
845 __os.precision(__precision);
850 template<typename _IntType, typename _RealType>
851 template<class _UniformRandomNumberGenerator>
852 typename geometric_distribution<_IntType, _RealType>::result_type
853 geometric_distribution<_IntType, _RealType>::
854 operator()(_UniformRandomNumberGenerator& __urng)
856 // About the epsilon thing see this thread:
857 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
858 const _RealType __naf =
859 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
860 // The largest _RealType convertible to _IntType.
861 const _RealType __thr =
862 std::numeric_limits<_IntType>::max() + __naf;
866 __cand = std::ceil(std::log(__urng()) / _M_log_p);
867 while (__cand >= __thr);
869 return result_type(__cand + __naf);
872 template<typename _IntType, typename _RealType,
873 typename _CharT, typename _Traits>
874 std::basic_ostream<_CharT, _Traits>&
875 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
876 const geometric_distribution<_IntType, _RealType>& __x)
878 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
879 typedef typename __ostream_type::ios_base __ios_base;
881 const typename __ios_base::fmtflags __flags = __os.flags();
882 const _CharT __fill = __os.fill();
883 const std::streamsize __precision = __os.precision();
884 __os.flags(__ios_base::scientific | __ios_base::left);
885 __os.fill(__os.widen(' '));
886 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
892 __os.precision(__precision);
897 template<typename _IntType, typename _RealType>
899 poisson_distribution<_IntType, _RealType>::
902 #if _GLIBCXX_USE_C99_MATH_TR1
905 const _RealType __m = std::floor(_M_mean);
906 _M_lm_thr = std::log(_M_mean);
907 _M_lfm = std::_GLIBCXX_TR1 lgamma(__m + 1);
908 _M_sm = std::sqrt(__m);
910 const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
911 const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
913 _M_d = std::_GLIBCXX_TR1 round(std::max(_RealType(6),
914 std::min(__m, __dx)));
915 const _RealType __cx = 2 * __m + _M_d;
916 _M_scx = std::sqrt(__cx / 2);
919 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
920 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
924 _M_lm_thr = std::exp(-_M_mean);
928 * A rejection algorithm when mean >= 12 and a simple method based
929 * upon the multiplication of uniform random variates otherwise.
930 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
934 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
935 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
937 template<typename _IntType, typename _RealType>
938 template<class _UniformRandomNumberGenerator>
939 typename poisson_distribution<_IntType, _RealType>::result_type
940 poisson_distribution<_IntType, _RealType>::
941 operator()(_UniformRandomNumberGenerator& __urng)
943 #if _GLIBCXX_USE_C99_MATH_TR1
948 // See comments above...
949 const _RealType __naf =
950 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
951 const _RealType __thr =
952 std::numeric_limits<_IntType>::max() + __naf;
954 const _RealType __m = std::floor(_M_mean);
956 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
957 const _RealType __c1 = _M_sm * __spi_2;
958 const _RealType __c2 = _M_c2b + __c1;
959 const _RealType __c3 = __c2 + 1;
960 const _RealType __c4 = __c3 + 1;
962 const _RealType __e178 = 1.0129030479320018583185514777512983L;
963 const _RealType __c5 = __c4 + __e178;
964 const _RealType __c = _M_cb + __c5;
965 const _RealType __2cx = 2 * (2 * __m + _M_d);
967 bool __reject = true;
970 const _RealType __u = __c * __urng();
971 const _RealType __e = -std::log(__urng());
977 const _RealType __n = _M_nd(__urng);
978 const _RealType __y = -std::abs(__n) * _M_sm - 1;
979 __x = std::floor(__y);
980 __w = -__n * __n / 2;
984 else if (__u <= __c2)
986 const _RealType __n = _M_nd(__urng);
987 const _RealType __y = 1 + std::abs(__n) * _M_scx;
988 __x = std::ceil(__y);
989 __w = __y * (2 - __y) * _M_1cx;
993 else if (__u <= __c3)
994 // NB: This case not in the book, nor in the Errata,
995 // but should be ok...
997 else if (__u <= __c4)
999 else if (__u <= __c5)
1003 const _RealType __v = -std::log(__urng());
1004 const _RealType __y = _M_d + __v * __2cx / _M_d;
1005 __x = std::ceil(__y);
1006 __w = -_M_d * _M_1cx * (1 + __y / 2);
1009 __reject = (__w - __e - __x * _M_lm_thr
1010 > _M_lfm - std::_GLIBCXX_TR1 lgamma(__x + __m + 1));
1012 __reject |= __x + __m >= __thr;
1016 return result_type(__x + __m + __naf);
1022 _RealType __prod = 1.0;
1029 while (__prod > _M_lm_thr);
1035 template<typename _IntType, typename _RealType,
1036 typename _CharT, typename _Traits>
1037 std::basic_ostream<_CharT, _Traits>&
1038 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1039 const poisson_distribution<_IntType, _RealType>& __x)
1041 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1042 typedef typename __ostream_type::ios_base __ios_base;
1044 const typename __ios_base::fmtflags __flags = __os.flags();
1045 const _CharT __fill = __os.fill();
1046 const std::streamsize __precision = __os.precision();
1047 const _CharT __space = __os.widen(' ');
1048 __os.flags(__ios_base::scientific | __ios_base::left);
1050 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1052 __os << __x.mean() << __space << __x._M_nd;
1054 __os.flags(__flags);
1056 __os.precision(__precision);
1060 template<typename _IntType, typename _RealType,
1061 typename _CharT, typename _Traits>
1062 std::basic_istream<_CharT, _Traits>&
1063 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1064 poisson_distribution<_IntType, _RealType>& __x)
1066 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1067 typedef typename __istream_type::ios_base __ios_base;
1069 const typename __ios_base::fmtflags __flags = __is.flags();
1070 __is.flags(__ios_base::skipws);
1072 __is >> __x._M_mean >> __x._M_nd;
1073 __x._M_initialize();
1075 __is.flags(__flags);
1080 template<typename _IntType, typename _RealType>
1082 binomial_distribution<_IntType, _RealType>::
1085 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1089 #if _GLIBCXX_USE_C99_MATH_TR1
1090 if (_M_t * __p12 >= 8)
1093 const _RealType __np = std::floor(_M_t * __p12);
1094 const _RealType __pa = __np / _M_t;
1095 const _RealType __1p = 1 - __pa;
1097 const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
1098 const _RealType __d1x =
1099 std::sqrt(__np * __1p * std::log(32 * __np
1100 / (81 * __pi_4 * __1p)));
1101 _M_d1 = std::_GLIBCXX_TR1 round(std::max(_RealType(1), __d1x));
1102 const _RealType __d2x =
1103 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1104 / (__pi_4 * __pa)));
1105 _M_d2 = std::_GLIBCXX_TR1 round(std::max(_RealType(1), __d2x));
1108 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1109 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1110 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1111 _M_c = 2 * _M_d1 / __np;
1112 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1113 const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
1114 const _RealType __s1s = _M_s1 * _M_s1;
1115 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1117 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1118 const _RealType __s2s = _M_s2 * _M_s2;
1119 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1120 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1121 _M_lf = (std::_GLIBCXX_TR1 lgamma(__np + 1)
1122 + std::_GLIBCXX_TR1 lgamma(_M_t - __np + 1));
1123 _M_lp1p = std::log(__pa / __1p);
1125 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1129 _M_q = -std::log(1 - __p12);
1132 template<typename _IntType, typename _RealType>
1133 template<class _UniformRandomNumberGenerator>
1134 typename binomial_distribution<_IntType, _RealType>::result_type
1135 binomial_distribution<_IntType, _RealType>::
1136 _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1139 _RealType __sum = 0;
1143 const _RealType __e = -std::log(__urng());
1144 __sum += __e / (__t - __x);
1147 while (__sum <= _M_q);
1153 * A rejection algorithm when t * p >= 8 and a simple waiting time
1154 * method - the second in the referenced book - otherwise.
1155 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1159 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1160 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1162 template<typename _IntType, typename _RealType>
1163 template<class _UniformRandomNumberGenerator>
1164 typename binomial_distribution<_IntType, _RealType>::result_type
1165 binomial_distribution<_IntType, _RealType>::
1166 operator()(_UniformRandomNumberGenerator& __urng)
1169 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1171 #if _GLIBCXX_USE_C99_MATH_TR1
1176 // See comments above...
1177 const _RealType __naf =
1178 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
1179 const _RealType __thr =
1180 std::numeric_limits<_IntType>::max() + __naf;
1182 const _RealType __np = std::floor(_M_t * __p12);
1183 const _RealType __pa = __np / _M_t;
1186 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1187 const _RealType __a1 = _M_a1;
1188 const _RealType __a12 = __a1 + _M_s2 * __spi_2;
1189 const _RealType __a123 = _M_a123;
1190 const _RealType __s1s = _M_s1 * _M_s1;
1191 const _RealType __s2s = _M_s2 * _M_s2;
1196 const _RealType __u = _M_s * __urng();
1202 const _RealType __n = _M_nd(__urng);
1203 const _RealType __y = _M_s1 * std::abs(__n);
1204 __reject = __y >= _M_d1;
1207 const _RealType __e = -std::log(__urng());
1208 __x = std::floor(__y);
1209 __v = -__e - __n * __n / 2 + _M_c;
1212 else if (__u <= __a12)
1214 const _RealType __n = _M_nd(__urng);
1215 const _RealType __y = _M_s2 * std::abs(__n);
1216 __reject = __y >= _M_d2;
1219 const _RealType __e = -std::log(__urng());
1220 __x = std::floor(-__y);
1221 __v = -__e - __n * __n / 2;
1224 else if (__u <= __a123)
1226 const _RealType __e1 = -std::log(__urng());
1227 const _RealType __e2 = -std::log(__urng());
1229 const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
1230 __x = std::floor(__y);
1231 __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
1232 -__y / (2 * __s1s)));
1237 const _RealType __e1 = -std::log(__urng());
1238 const _RealType __e2 = -std::log(__urng());
1240 const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
1241 __x = std::floor(-__y);
1242 __v = -__e2 - _M_d2 * __y / (2 * __s2s);
1246 __reject = __reject || __x < -__np || __x > _M_t - __np;
1249 const _RealType __lfx =
1250 std::_GLIBCXX_TR1 lgamma(__np + __x + 1)
1251 + std::_GLIBCXX_TR1 lgamma(_M_t - (__np + __x) + 1);
1252 __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
1255 __reject |= __x + __np >= __thr;
1259 __x += __np + __naf;
1261 const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x));
1262 __ret = _IntType(__x) + __z;
1266 __ret = _M_waiting(__urng, _M_t);
1269 __ret = _M_t - __ret;
1273 template<typename _IntType, typename _RealType,
1274 typename _CharT, typename _Traits>
1275 std::basic_ostream<_CharT, _Traits>&
1276 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1277 const binomial_distribution<_IntType, _RealType>& __x)
1279 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1280 typedef typename __ostream_type::ios_base __ios_base;
1282 const typename __ios_base::fmtflags __flags = __os.flags();
1283 const _CharT __fill = __os.fill();
1284 const std::streamsize __precision = __os.precision();
1285 const _CharT __space = __os.widen(' ');
1286 __os.flags(__ios_base::scientific | __ios_base::left);
1288 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1290 __os << __x.t() << __space << __x.p()
1291 << __space << __x._M_nd;
1293 __os.flags(__flags);
1295 __os.precision(__precision);
1299 template<typename _IntType, typename _RealType,
1300 typename _CharT, typename _Traits>
1301 std::basic_istream<_CharT, _Traits>&
1302 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1303 binomial_distribution<_IntType, _RealType>& __x)
1305 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1306 typedef typename __istream_type::ios_base __ios_base;
1308 const typename __ios_base::fmtflags __flags = __is.flags();
1309 __is.flags(__ios_base::dec | __ios_base::skipws);
1311 __is >> __x._M_t >> __x._M_p >> __x._M_nd;
1312 __x._M_initialize();
1314 __is.flags(__flags);
1319 template<typename _RealType, typename _CharT, typename _Traits>
1320 std::basic_ostream<_CharT, _Traits>&
1321 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1322 const uniform_real<_RealType>& __x)
1324 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1325 typedef typename __ostream_type::ios_base __ios_base;
1327 const typename __ios_base::fmtflags __flags = __os.flags();
1328 const _CharT __fill = __os.fill();
1329 const std::streamsize __precision = __os.precision();
1330 const _CharT __space = __os.widen(' ');
1331 __os.flags(__ios_base::scientific | __ios_base::left);
1333 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1335 __os << __x.min() << __space << __x.max();
1337 __os.flags(__flags);
1339 __os.precision(__precision);
1343 template<typename _RealType, typename _CharT, typename _Traits>
1344 std::basic_istream<_CharT, _Traits>&
1345 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1346 uniform_real<_RealType>& __x)
1348 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1349 typedef typename __istream_type::ios_base __ios_base;
1351 const typename __ios_base::fmtflags __flags = __is.flags();
1352 __is.flags(__ios_base::skipws);
1354 __is >> __x._M_min >> __x._M_max;
1356 __is.flags(__flags);
1361 template<typename _RealType, typename _CharT, typename _Traits>
1362 std::basic_ostream<_CharT, _Traits>&
1363 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1364 const exponential_distribution<_RealType>& __x)
1366 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1367 typedef typename __ostream_type::ios_base __ios_base;
1369 const typename __ios_base::fmtflags __flags = __os.flags();
1370 const _CharT __fill = __os.fill();
1371 const std::streamsize __precision = __os.precision();
1372 __os.flags(__ios_base::scientific | __ios_base::left);
1373 __os.fill(__os.widen(' '));
1374 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1376 __os << __x.lambda();
1378 __os.flags(__flags);
1380 __os.precision(__precision);
1386 * Polar method due to Marsaglia.
1388 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1389 * New York, 1986, Ch. V, Sect. 4.4.
1391 template<typename _RealType>
1392 template<class _UniformRandomNumberGenerator>
1393 typename normal_distribution<_RealType>::result_type
1394 normal_distribution<_RealType>::
1395 operator()(_UniformRandomNumberGenerator& __urng)
1399 if (_M_saved_available)
1401 _M_saved_available = false;
1406 result_type __x, __y, __r2;
1409 __x = result_type(2.0) * __urng() - 1.0;
1410 __y = result_type(2.0) * __urng() - 1.0;
1411 __r2 = __x * __x + __y * __y;
1413 while (__r2 > 1.0 || __r2 == 0.0);
1415 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1416 _M_saved = __x * __mult;
1417 _M_saved_available = true;
1418 __ret = __y * __mult;
1421 __ret = __ret * _M_sigma + _M_mean;
1425 template<typename _RealType, typename _CharT, typename _Traits>
1426 std::basic_ostream<_CharT, _Traits>&
1427 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1428 const normal_distribution<_RealType>& __x)
1430 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1431 typedef typename __ostream_type::ios_base __ios_base;
1433 const typename __ios_base::fmtflags __flags = __os.flags();
1434 const _CharT __fill = __os.fill();
1435 const std::streamsize __precision = __os.precision();
1436 const _CharT __space = __os.widen(' ');
1437 __os.flags(__ios_base::scientific | __ios_base::left);
1439 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1441 __os << __x._M_saved_available << __space
1442 << __x.mean() << __space
1444 if (__x._M_saved_available)
1445 __os << __space << __x._M_saved;
1447 __os.flags(__flags);
1449 __os.precision(__precision);
1453 template<typename _RealType, typename _CharT, typename _Traits>
1454 std::basic_istream<_CharT, _Traits>&
1455 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1456 normal_distribution<_RealType>& __x)
1458 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1459 typedef typename __istream_type::ios_base __ios_base;
1461 const typename __ios_base::fmtflags __flags = __is.flags();
1462 __is.flags(__ios_base::dec | __ios_base::skipws);
1464 __is >> __x._M_saved_available >> __x._M_mean
1466 if (__x._M_saved_available)
1467 __is >> __x._M_saved;
1469 __is.flags(__flags);
1474 template<typename _RealType>
1476 gamma_distribution<_RealType>::
1480 _M_l_d = std::sqrt(2 * _M_alpha - 1);
1482 _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
1487 * Cheng's rejection algorithm GB for alpha >= 1 and a modification
1488 * of Vaduva's rejection from Weibull algorithm due to Devroye for
1492 * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
1493 * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
1495 * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
1496 * and Composition Procedures." Math. Operationsforschung and Statistik,
1497 * Series in Statistics, 8, 545-576, 1977.
1499 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1500 * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
1502 template<typename _RealType>
1503 template<class _UniformRandomNumberGenerator>
1504 typename gamma_distribution<_RealType>::result_type
1505 gamma_distribution<_RealType>::
1506 operator()(_UniformRandomNumberGenerator& __urng)
1514 const result_type __b = _M_alpha
1515 - result_type(1.3862943611198906188344642429163531L);
1516 const result_type __c = _M_alpha + _M_l_d;
1517 const result_type __1l = 1 / _M_l_d;
1520 const result_type __k = 2.5040773967762740733732583523868748L;
1524 const result_type __u = __urng();
1525 const result_type __v = __urng();
1527 const result_type __y = __1l * std::log(__v / (1 - __v));
1528 __x = _M_alpha * std::exp(__y);
1530 const result_type __z = __u * __v * __v;
1531 const result_type __r = __b + __c * __y - __x;
1533 __reject = __r < result_type(4.5) * __z - __k;
1535 __reject = __r < std::log(__z);
1541 const result_type __c = 1 / _M_alpha;
1545 const result_type __z = -std::log(__urng());
1546 const result_type __e = -std::log(__urng());
1548 __x = std::pow(__z, __c);
1550 __reject = __z + __e < _M_l_d + __x;
1558 template<typename _RealType, typename _CharT, typename _Traits>
1559 std::basic_ostream<_CharT, _Traits>&
1560 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1561 const gamma_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 __os.flags(__ios_base::scientific | __ios_base::left);
1570 __os.fill(__os.widen(' '));
1571 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1573 __os << __x.alpha();
1575 __os.flags(__flags);
1577 __os.precision(__precision);
1581 _GLIBCXX_END_NAMESPACE_TR1