1 // random number generation (out of line) -*- C++ -*-
3 // Copyright (C) 2006 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.
32 _GLIBCXX_BEGIN_NAMESPACE(tr1)
35 * Implementation-space details.
39 // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
42 // Because a and c are compile-time integral constants the compiler kindly
43 // elides any unreachable paths.
45 // Preconditions: a > 0, m > 0.
47 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
57 static const _Tp __q = __m / __a;
58 static const _Tp __r = __m % __a;
60 _Tp __t1 = __a * (__x % __q);
61 _Tp __t2 = __r * (__x / __q);
65 __x = __m - __t2 + __t1;
70 const _Tp __d = __m - __x;
80 // Special case for m == 0 -- use unsigned integer overflow as modulo
82 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
83 struct _Mod<_Tp, __a, __c, __m, true>
87 { return __a * __x + __c; }
90 // Dispatch based on modulus value to prevent divide-by-zero compile-time
91 // errors when m == 0.
92 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
95 { return _Mod<_Tp, __a, __c, __m, __m == 0>::__calc(__x); }
98 template<typename _RealType>
101 static const std::streamsize __value =
102 2 + std::numeric_limits<_RealType>::digits * 3010/10000;
105 template<typename _ValueT>
106 struct _To_Unsigned_Type
107 { typedef _ValueT _Type; };
110 struct _To_Unsigned_Type<short>
111 { typedef unsigned short _Type; };
114 struct _To_Unsigned_Type<int>
115 { typedef unsigned int _Type; };
118 struct _To_Unsigned_Type<long>
119 { typedef unsigned long _Type; };
121 #ifdef _GLIBCXX_USE_LONG_LONG
123 struct _To_Unsigned_Type<long long>
124 { typedef unsigned long long _Type; };
127 } // namespace _Private
131 * Seeds the LCR with integral value @p __x0, adjusted so that the
132 * ring identity is never a member of the convergence set.
134 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
136 linear_congruential<_UIntType, __a, __c, __m>::
137 seed(unsigned long __x0)
139 if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
140 && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
141 _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
143 _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
147 * Seeds the LCR engine with a value generated by @p __g.
149 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
152 linear_congruential<_UIntType, __a, __c, __m>::
153 seed(_Gen& __g, false_type)
155 _UIntType __x0 = __g();
156 if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
157 && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
158 _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
160 _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
164 * Returns a value that is less than or equal to all values potentially
165 * returned by operator(). The return value of this function does not
166 * change during the lifetime of the object..
168 * The minumum depends on the @p __c parameter: if it is zero, the
169 * minimum generated must be > 0, otherwise 0 is allowed.
171 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
172 typename linear_congruential<_UIntType, __a, __c, __m>::result_type
173 linear_congruential<_UIntType, __a, __c, __m>::
175 { return (_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0) ? 1 : 0; }
178 * Gets the maximum possible value of the generated range.
180 * For a linear congruential generator, the maximum is always @p __m - 1.
182 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
183 typename linear_congruential<_UIntType, __a, __c, __m>::result_type
184 linear_congruential<_UIntType, __a, __c, __m>::
186 { return (__m == 0) ? std::numeric_limits<_UIntType>::max() : (__m - 1); }
189 * Gets the next generated value in sequence.
191 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
192 typename linear_congruential<_UIntType, __a, __c, __m>::result_type
193 linear_congruential<_UIntType, __a, __c, __m>::
196 _M_x = _Private::__mod<_UIntType, __a, __c, __m>(_M_x);
200 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
201 typename _CharT, typename _Traits>
202 std::basic_ostream<_CharT, _Traits>&
203 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
204 const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
206 const std::ios_base::fmtflags __flags = __os.flags();
207 const _CharT __fill = __os.fill();
208 __os.flags(std::ios_base::dec | std::ios_base::fixed
209 | std::ios_base::left);
210 __os.fill(__os.widen(' '));
219 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
220 typename _CharT, typename _Traits>
221 std::basic_istream<_CharT, _Traits>&
222 operator>>(std::basic_istream<_CharT, _Traits>& __is,
223 linear_congruential<_UIntType, __a, __c, __m>& __lcr)
225 const std::ios_base::fmtflags __flags = __is.flags();
226 __is.flags(std::ios_base::dec);
235 template<class _UIntType, int __w, int __n, int __m, int __r,
236 _UIntType __a, int __u, int __s,
237 _UIntType __b, int __t, _UIntType __c, int __l>
239 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
240 __b, __t, __c, __l>::
241 seed(unsigned long __value)
243 _M_x[0] = _Private::__mod<_UIntType, 1, 0,
244 _Private::_Shift<_UIntType, __w>::__value>(__value);
246 for (int __i = 1; __i < state_size; ++__i)
248 _UIntType __x = _M_x[__i - 1];
249 __x ^= __x >> (__w - 2);
252 _M_x[__i] = _Private::__mod<_UIntType, 1, 0,
253 _Private::_Shift<_UIntType, __w>::__value>(__x);
258 template<class _UIntType, int __w, int __n, int __m, int __r,
259 _UIntType __a, int __u, int __s,
260 _UIntType __b, int __t, _UIntType __c, int __l>
263 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
264 __b, __t, __c, __l>::
265 seed(_Gen& __gen, false_type)
267 for (int __i = 0; __i < state_size; ++__i)
268 _M_x[__i] = _Private::__mod<_UIntType, 1, 0,
269 _Private::_Shift<_UIntType, __w>::__value>(__gen());
273 template<class _UIntType, int __w, int __n, int __m, int __r,
274 _UIntType __a, int __u, int __s,
275 _UIntType __b, int __t, _UIntType __c, int __l>
277 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
278 __b, __t, __c, __l>::result_type
279 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
280 __b, __t, __c, __l>::
283 // Reload the vector - cost is O(n) amortized over n calls.
284 if (_M_p >= state_size)
286 const _UIntType __upper_mask = (~_UIntType()) << __r;
287 const _UIntType __lower_mask = ~__upper_mask;
289 for (int __k = 0; __k < (__n - __m); ++__k)
291 _UIntType __y = ((_M_x[__k] & __upper_mask)
292 | (_M_x[__k + 1] & __lower_mask));
293 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
294 ^ ((__y & 0x01) ? __a : 0));
297 for (int __k = (__n - __m); __k < (__n - 1); ++__k)
299 _UIntType __y = ((_M_x[__k] & __upper_mask)
300 | (_M_x[__k + 1] & __lower_mask));
301 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
302 ^ ((__y & 0x01) ? __a : 0));
305 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
306 | (_M_x[0] & __lower_mask));
307 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
308 ^ ((__y & 0x01) ? __a : 0));
312 // Calculate o(x(i)).
313 result_type __z = _M_x[_M_p++];
315 __z ^= (__z << __s) & __b;
316 __z ^= (__z << __t) & __c;
322 template<class _UIntType, int __w, int __n, int __m, int __r,
323 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
324 _UIntType __c, int __l,
325 typename _CharT, typename _Traits>
326 std::basic_ostream<_CharT, _Traits>&
327 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
328 const mersenne_twister<_UIntType, __w, __n, __m,
329 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
331 const std::ios_base::fmtflags __flags = __os.flags();
332 const _CharT __fill = __os.fill();
333 const _CharT __space = __os.widen(' ');
334 __os.flags(std::ios_base::dec | std::ios_base::fixed
335 | std::ios_base::left);
338 for (int __i = 0; __i < __n - 1; ++__i)
339 __os << __x._M_x[__i] << __space;
340 __os << __x._M_x[__n - 1];
347 template<class _UIntType, int __w, int __n, int __m, int __r,
348 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
349 _UIntType __c, int __l,
350 typename _CharT, typename _Traits>
351 std::basic_istream<_CharT, _Traits>&
352 operator>>(std::basic_istream<_CharT, _Traits>& __is,
353 mersenne_twister<_UIntType, __w, __n, __m,
354 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
356 const std::ios_base::fmtflags __flags = __is.flags();
357 __is.flags(std::ios_base::dec | std::ios_base::skipws);
359 for (int __i = 0; __i < __n; ++__i)
360 __is >> __x._M_x[__i];
367 template<typename _IntType, _IntType __m, int __s, int __r>
369 subtract_with_carry<_IntType, __m, __s, __r>::
370 seed(unsigned long __value)
375 std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
378 for (int __i = 0; __i < long_lag; ++__i)
379 _M_x[__i] = _Private::__mod<_IntType, 1, 0, modulus>(__lcg());
381 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
385 template<typename _IntType, _IntType __m, int __s, int __r>
388 subtract_with_carry<_IntType, __m, __s, __r>::
389 seed(_Gen& __gen, false_type)
391 const int __n = (std::numeric_limits<_IntType>::digits + 31) / 32;
393 typedef typename _Private::_Select<(sizeof(unsigned) == 4),
394 unsigned, unsigned long>::_Type _UInt32Type;
396 typedef typename _Private::_To_Unsigned_Type<_IntType>::_Type
399 for (int __i = 0; __i < long_lag; ++__i)
402 _UIntType __factor = 1;
403 for (int __j = 0; __j < __n; ++__j)
405 __tmp += (_Private::__mod<_UInt32Type, 1, 0, 0>(__gen())
407 __factor *= _Private::_Shift<_UIntType, 32>::__value;
409 _M_x[__i] = _Private::__mod<_UIntType, 1, 0, modulus>(__tmp);
411 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
415 template<typename _IntType, _IntType __m, int __s, int __r>
416 typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
417 subtract_with_carry<_IntType, __m, __s, __r>::
420 // Derive short lag index from current index.
421 int __ps = _M_p - short_lag;
425 // Calculate new x(i) without overflow or division.
427 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
429 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
434 __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
439 // Adjust current index to loop around in ring buffer.
440 if (_M_p >= long_lag)
446 template<typename _IntType, _IntType __m, int __s, int __r,
447 typename _CharT, typename _Traits>
448 std::basic_ostream<_CharT, _Traits>&
449 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
450 const subtract_with_carry<_IntType, __m, __s, __r>& __x)
452 const std::ios_base::fmtflags __flags = __os.flags();
453 const _CharT __fill = __os.fill();
454 const _CharT __space = __os.widen(' ');
455 __os.flags(std::ios_base::dec | std::ios_base::fixed
456 | std::ios_base::left);
459 for (int __i = 0; __i < __r; ++__i)
460 __os << __x._M_x[__i] << __space;
461 __os << __x._M_carry;
468 template<typename _IntType, _IntType __m, int __s, int __r,
469 typename _CharT, typename _Traits>
470 std::basic_istream<_CharT, _Traits>&
471 operator>>(std::basic_istream<_CharT, _Traits>& __is,
472 subtract_with_carry<_IntType, __m, __s, __r>& __x)
474 const std::ios_base::fmtflags __flags = __is.flags();
475 __is.flags(std::ios_base::dec | std::ios_base::skipws);
477 for (int __i = 0; __i < __r; ++__i)
478 __is >> __x._M_x[__i];
479 __is >> __x._M_carry;
486 template<class _UniformRandomNumberGenerator, int __p, int __r>
487 typename discard_block<_UniformRandomNumberGenerator,
488 __p, __r>::result_type
489 discard_block<_UniformRandomNumberGenerator, __p, __r>::
492 if (_M_n >= used_block)
494 while (_M_n < block_size)
505 template<class _UniformRandomNumberGenerator, int __p, int __r,
506 typename _CharT, typename _Traits>
507 std::basic_ostream<_CharT, _Traits>&
508 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
509 const discard_block<_UniformRandomNumberGenerator,
512 const std::ios_base::fmtflags __flags = __os.flags();
513 const _CharT __fill = __os.fill();
514 const _CharT __space = __os.widen(' ');
515 __os.flags(std::ios_base::dec | std::ios_base::fixed
516 | std::ios_base::left);
519 __os << __x._M_b << __space << __x._M_n;
526 template<class _UniformRandomNumberGenerator, int __p, int __r,
527 typename _CharT, typename _Traits>
528 std::basic_istream<_CharT, _Traits>&
529 operator>>(std::basic_istream<_CharT, _Traits>& __is,
530 discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
532 const std::ios_base::fmtflags __flags = __is.flags();
533 __is.flags(std::ios_base::dec | std::ios_base::skipws);
535 __is >> __x._M_b >> __x._M_n;
542 template<class _UniformRandomNumberGenerator1, int __s1,
543 class _UniformRandomNumberGenerator2, int __s2,
544 typename _CharT, typename _Traits>
545 std::basic_ostream<_CharT, _Traits>&
546 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
547 const xor_combine<_UniformRandomNumberGenerator1, __s1,
548 _UniformRandomNumberGenerator2, __s2>& __x)
550 const std::ios_base::fmtflags __flags = __os.flags();
551 const _CharT __fill = __os.fill();
552 const _CharT __space = __os.widen(' ');
553 __os.flags(std::ios_base::dec | std::ios_base::fixed
554 | std::ios_base::left);
557 __os << __x.base1() << __space << __x.base2();
564 template<class _UniformRandomNumberGenerator1, int __s1,
565 class _UniformRandomNumberGenerator2, int __s2,
566 typename _CharT, typename _Traits>
567 std::basic_istream<_CharT, _Traits>&
568 operator>>(std::basic_istream<_CharT, _Traits>& __is,
569 xor_combine<_UniformRandomNumberGenerator1, __s1,
570 _UniformRandomNumberGenerator2, __s2>& __x)
572 const std::ios_base::fmtflags __flags = __is.flags();
573 __is.flags(std::ios_base::skipws);
575 __is >> __x._M_b1 >> __x._M_b2;
582 template<typename _IntType, typename _CharT, typename _Traits>
583 std::basic_ostream<_CharT, _Traits>&
584 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
585 const uniform_int<_IntType>& __x)
587 const std::ios_base::fmtflags __flags = __os.flags();
588 const _CharT __fill = __os.fill();
589 const _CharT __space = __os.widen(' ');
590 __os.flags(std::ios_base::scientific | std::ios_base::left);
593 __os << __x.min() << __space << __x.max();
600 template<typename _IntType, typename _CharT, typename _Traits>
601 std::basic_istream<_CharT, _Traits>&
602 operator>>(std::basic_istream<_CharT, _Traits>& __is,
603 uniform_int<_IntType>& __x)
605 const std::ios_base::fmtflags __flags = __is.flags();
606 __is.flags(std::ios_base::dec | std::ios_base::skipws);
608 __is >> __x._M_min >> __x._M_max;
615 template<typename _CharT, typename _Traits>
616 std::basic_ostream<_CharT, _Traits>&
617 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
618 const bernoulli_distribution& __x)
620 const std::ios_base::fmtflags __flags = __os.flags();
621 const _CharT __fill = __os.fill();
622 const std::streamsize __precision = __os.precision();
623 __os.flags(std::ios_base::scientific | std::ios_base::left);
624 __os.fill(__os.widen(' '));
625 __os.precision(_Private::_Max_digits10<double>::__value);
631 __os.precision(__precision);
636 template<typename _IntType, typename _RealType,
637 typename _CharT, typename _Traits>
638 std::basic_ostream<_CharT, _Traits>&
639 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
640 const geometric_distribution<_IntType, _RealType>& __x)
642 const std::ios_base::fmtflags __flags = __os.flags();
643 const _CharT __fill = __os.fill();
644 const std::streamsize __precision = __os.precision();
645 __os.flags(std::ios_base::scientific | std::ios_base::left);
646 __os.fill(__os.widen(' '));
647 __os.precision(_Private::_Max_digits10<_RealType>::__value);
653 __os.precision(__precision);
658 template<typename _RealType, typename _CharT, typename _Traits>
659 std::basic_ostream<_CharT, _Traits>&
660 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
661 const uniform_real<_RealType>& __x)
663 const std::ios_base::fmtflags __flags = __os.flags();
664 const _CharT __fill = __os.fill();
665 const std::streamsize __precision = __os.precision();
666 const _CharT __space = __os.widen(' ');
667 __os.flags(std::ios_base::scientific | std::ios_base::left);
669 __os.precision(_Private::_Max_digits10<_RealType>::__value);
671 __os << __x.min() << __space << __x.max();
675 __os.precision(__precision);
679 template<typename _RealType, typename _CharT, typename _Traits>
680 std::basic_istream<_CharT, _Traits>&
681 operator>>(std::basic_istream<_CharT, _Traits>& __is,
682 uniform_real<_RealType>& __x)
684 const std::ios_base::fmtflags __flags = __is.flags();
685 __is.flags(std::ios_base::skipws);
687 __is >> __x._M_min >> __x._M_max;
694 template<typename _RealType, typename _CharT, typename _Traits>
695 std::basic_ostream<_CharT, _Traits>&
696 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
697 const exponential_distribution<_RealType>& __x)
699 const std::ios_base::fmtflags __flags = __os.flags();
700 const _CharT __fill = __os.fill();
701 const std::streamsize __precision = __os.precision();
702 __os.flags(std::ios_base::scientific | std::ios_base::left);
703 __os.fill(__os.widen(' '));
704 __os.precision(_Private::_Max_digits10<_RealType>::__value);
706 __os << __x.lambda();
710 __os.precision(__precision);
716 * Classic Box-Muller method.
719 * Box, G. E. P. and Muller, M. E. "A Note on the Generation of
720 * Random Normal Deviates." Ann. Math. Stat. 29, 610-611, 1958.
722 template<typename _RealType>
723 template<class _UniformRandomNumberGenerator>
724 typename normal_distribution<_RealType>::result_type
725 normal_distribution<_RealType>::
726 operator()(_UniformRandomNumberGenerator& __urng)
730 if (_M_saved_available)
732 _M_saved_available = false;
737 result_type __x, __y, __r2;
740 __x = result_type(2.0) * __urng() - result_type(1.0);
741 __y = result_type(2.0) * __urng() - result_type(1.0);
742 __r2 = __x * __x + __y * __y;
744 while (__r2 > result_type(1.0) || __r2 == result_type(0));
746 const result_type __mult = std::sqrt(-result_type(2.0)
747 * std::log(__r2) / __r2);
748 _M_saved = __x * __mult;
749 _M_saved_available = true;
750 __ret = __y * __mult;
753 return __ret * _M_sigma + _M_mean;
756 template<typename _RealType, typename _CharT, typename _Traits>
757 std::basic_ostream<_CharT, _Traits>&
758 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
759 const normal_distribution<_RealType>& __x)
761 const std::ios_base::fmtflags __flags = __os.flags();
762 const _CharT __fill = __os.fill();
763 const std::streamsize __precision = __os.precision();
764 const _CharT __space = __os.widen(' ');
765 __os.flags(std::ios_base::scientific | std::ios_base::left);
767 __os.precision(_Private::_Max_digits10<_RealType>::__value);
769 __os << __x.mean() << __space
770 << __x.sigma() << __space
771 << __x._M_saved << __space
772 << __x._M_saved_available;
776 __os.precision(__precision);
780 template<typename _RealType, typename _CharT, typename _Traits>
781 std::basic_istream<_CharT, _Traits>&
782 operator>>(std::basic_istream<_CharT, _Traits>& __is,
783 normal_distribution<_RealType>& __x)
785 const std::ios_base::fmtflags __flags = __is.flags();
786 __is.flags(std::ios_base::dec | std::ios_base::skipws);
788 __is >> __x._M_mean >> __x._M_sigma
789 >> __x._M_saved >> __x._M_saved_available;
797 * Cheng's rejection algorithm GB for alpha >= 1 and a modification
798 * of Vaduva's rejection from Weibull algorithm due to Devroye for
802 * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
803 * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
805 * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
806 * and Composition Procedures." Math. Operationsforschung and Statistik,
807 * Series in Statistics, 8, 545-576, 1977.
809 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
810 * New York, 1986, Sect. 3.4.
812 template<typename _RealType>
813 template<class _UniformRandomNumberGenerator>
814 typename gamma_distribution<_RealType>::result_type
815 gamma_distribution<_RealType>::
816 operator()(_UniformRandomNumberGenerator& __urng)
823 const result_type __b = _M_alpha
824 - result_type(1.3862943611198906188344642429163531L);
825 const result_type __c = _M_alpha + std::sqrt(2 * _M_alpha - 1);
828 const result_type __k = 2.5040773967762740733732583523868748L;
830 result_type __z, __r;
833 const result_type __u = __urng();
834 const result_type __v = __urng();
836 const result_type __y = _M_alpha * std::log(__v / (1 - __v));
837 __x = _M_alpha * std::exp(__v);
839 __z = __u * __v * __v;
840 __r = __b + __c * __y - __x;
842 while (__r < result_type(4.5) * __z - __k
843 && __r < std::log(__z));
847 const result_type __c = 1 / _M_alpha;
848 const result_type __d =
849 std::pow(_M_alpha, _M_alpha / (1 - _M_alpha)) * (1 - _M_alpha);
851 result_type __z, __e;
854 __z = -std::log(__urng());
855 __e = -std::log(__urng());
857 __x = std::pow(__z, __c);
859 while (__z + __e > __d + __x);
865 template<typename _RealType, typename _CharT, typename _Traits>
866 std::basic_ostream<_CharT, _Traits>&
867 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
868 const gamma_distribution<_RealType>& __x)
870 const std::ios_base::fmtflags __flags = __os.flags();
871 const _CharT __fill = __os.fill();
872 const std::streamsize __precision = __os.precision();
873 __os.flags(std::ios_base::scientific | std::ios_base::left);
874 __os.fill(__os.widen(' '));
875 __os.precision(_Private::_Max_digits10<_RealType>::__value);
881 __os.precision(__precision);
885 _GLIBCXX_END_NAMESPACE