]> git.ipfire.org Git - thirdparty/gcc.git/blame - libstdc++-v3/include/tr1/random.tcc
projects / thirdparty / gcc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update copyright date
[thirdparty/gcc.git] / libstdc++-v3 / include / tr1 / random.tcc
CommitLineData
8e79468d
BK		 1// random number generation (out of line) -*- C++ -*-
2
3// Copyright (C) 2007 Free Software Foundation, Inc.
4//
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)
9// any later version.
10
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.
15
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,
19// USA.
20
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.
29
30/** @file tr1/random.tcc
31 * This is an internal header file, included by other library headers.
32 * You should not attempt to use it directly.
33 */
34
35namespace std
36{
37namespace tr1
38{
39
40 /*
41 * (Further) implementation-space details.
42 */
43 namespace __detail
44 {
45 // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
46 // integer overflow.
47 //
48 // Because a and c are compile-time integral constants the compiler kindly
49 // elides any unreachable paths.
50 //
51 // Preconditions: a > 0, m > 0.
52 //
53 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
54 struct _Mod
55 {
56 static _Tp
57 __calc(_Tp __x)
58 {
59 if (__a == 1)
60 __x %= __m;
61 else
62 {
63 static const _Tp __q = __m / __a;
64 static const _Tp __r = __m % __a;
65
66 _Tp __t1 = __a * (__x % __q);
67 _Tp __t2 = __r * (__x / __q);
68 if (__t1 >= __t2)
69 __x = __t1 - __t2;
70 else
71 __x = __m - __t2 + __t1;
72 }
73
74 if (__c != 0)
75 {
76 const _Tp __d = __m - __x;
77 if (__d > __c)
78 __x += __c;
79 else
80 __x = __c - __d;
81 }
82 return __x;
83 }
84 };
85
86 // Special case for m == 0 -- use unsigned integer overflow as modulo
87 // operator.
88 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
89 struct _Mod<_Tp, __a, __c, __m, true>
90 {
91 static _Tp
92 __calc(_Tp __x)
93 { return __a * __x + __c; }
94 };
95 } // namespace __detail
96
97 /**
98 * Seeds the LCR with integral value @p __x0, adjusted so that the
99 * ring identity is never a member of the convergence set.
100 */
101 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
102 void
103 linear_congruential<_UIntType, __a, __c, __m>::
104 seed(unsigned long __x0)
105 {
106 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
107 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
108 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
109 else
110 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
111 }
112
113 /**
114 * Seeds the LCR engine with a value generated by @p __g.
115 */
116 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
117 template<class _Gen>
118 void
119 linear_congruential<_UIntType, __a, __c, __m>::
120 seed(_Gen& __g, false_type)
121 {
122 _UIntType __x0 = __g();
123 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
124 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
125 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
126 else
127 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
128 }
129
130 /**
131 * Gets the next generated value in sequence.
132 */
133 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
134 typename linear_congruential<_UIntType, __a, __c, __m>::result_type
135 linear_congruential<_UIntType, __a, __c, __m>::
136 operator()()
137 {
138 _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
139 return _M_x;
140 }
141
142 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
143 typename _CharT, typename _Traits>
144 std::basic_ostream<_CharT, _Traits>&
145 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
146 const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
147 {
148 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
149 typedef typename __ostream_type::ios_base __ios_base;
150
151 const typename __ios_base::fmtflags __flags = __os.flags();
152 const _CharT __fill = __os.fill();
153 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
154 __os.fill(__os.widen(' '));
155
156 __os << __lcr._M_x;
157
158 __os.flags(__flags);
159 __os.fill(__fill);
160 return __os;
161 }
162
163 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
164 typename _CharT, typename _Traits>
165 std::basic_istream<_CharT, _Traits>&
166 operator>>(std::basic_istream<_CharT, _Traits>& __is,
167 linear_congruential<_UIntType, __a, __c, __m>& __lcr)
168 {
169 typedef std::basic_istream<_CharT, _Traits> __istream_type;
170 typedef typename __istream_type::ios_base __ios_base;
171
172 const typename __ios_base::fmtflags __flags = __is.flags();
173 __is.flags(__ios_base::dec);
174
175 __is >> __lcr._M_x;
176
177 __is.flags(__flags);
178 return __is;
179 }
180
181
182 template<class _UIntType, int __w, int __n, int __m, int __r,
183 _UIntType __a, int __u, int __s,
184 _UIntType __b, int __t, _UIntType __c, int __l>
185 void
186 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
187 __b, __t, __c, __l>::
188 seed(unsigned long __value)
189 {
190 _M_x[0] = __detail::__mod<_UIntType, 1, 0,
191 __detail::_Shift<_UIntType, __w>::__value>(__value);
192
193 for (int __i = 1; __i < state_size; ++__i)
194 {
195 _UIntType __x = _M_x[__i - 1];
196 __x ^= __x >> (__w - 2);
197 __x *= 1812433253ul;
198 __x += __i;
199 _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
200 __detail::_Shift<_UIntType, __w>::__value>(__x);
201 }
202 _M_p = state_size;
203 }
204
205 template<class _UIntType, int __w, int __n, int __m, int __r,
206 _UIntType __a, int __u, int __s,
207 _UIntType __b, int __t, _UIntType __c, int __l>
208 template<class _Gen>
209 void
210 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
211 __b, __t, __c, __l>::
212 seed(_Gen& __gen, false_type)
213 {
214 for (int __i = 0; __i < state_size; ++__i)
215 _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
216 __detail::_Shift<_UIntType, __w>::__value>(__gen());
217 _M_p = state_size;
218 }
219
220 template<class _UIntType, int __w, int __n, int __m, int __r,
221 _UIntType __a, int __u, int __s,
222 _UIntType __b, int __t, _UIntType __c, int __l>
223 typename
224 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
225 __b, __t, __c, __l>::result_type
226 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
227 __b, __t, __c, __l>::
228 operator()()
229 {
230 // Reload the vector - cost is O(n) amortized over n calls.
231 if (_M_p >= state_size)
232 {
233 const _UIntType __upper_mask = (~_UIntType()) << __r;
234 const _UIntType __lower_mask = ~__upper_mask;
235
236 for (int __k = 0; __k < (__n - __m); ++__k)
237 {
238 _UIntType __y = ((_M_x[__k] & __upper_mask)
239 | (_M_x[__k + 1] & __lower_mask));
240 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
241 ^ ((__y & 0x01) ? __a : 0));
242 }
243
244 for (int __k = (__n - __m); __k < (__n - 1); ++__k)
245 {
246 _UIntType __y = ((_M_x[__k] & __upper_mask)
247 | (_M_x[__k + 1] & __lower_mask));
248 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
249 ^ ((__y & 0x01) ? __a : 0));
250 }
251
252 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
253 | (_M_x[0] & __lower_mask));
254 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
255 ^ ((__y & 0x01) ? __a : 0));
256 _M_p = 0;
257 }
258
259 // Calculate o(x(i)).
260 result_type __z = _M_x[_M_p++];
261 __z ^= (__z >> __u);
262 __z ^= (__z << __s) & __b;
263 __z ^= (__z << __t) & __c;
264 __z ^= (__z >> __l);
265
266 return __z;
267 }
268
269 template<class _UIntType, int __w, int __n, int __m, int __r,
270 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
271 _UIntType __c, int __l,
272 typename _CharT, typename _Traits>
273 std::basic_ostream<_CharT, _Traits>&
274 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
275 const mersenne_twister<_UIntType, __w, __n, __m,
276 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
277 {
278 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
279 typedef typename __ostream_type::ios_base __ios_base;
280
281 const typename __ios_base::fmtflags __flags = __os.flags();
282 const _CharT __fill = __os.fill();
283 const _CharT __space = __os.widen(' ');
284 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
285 __os.fill(__space);
286
287 for (int __i = 0; __i < __n - 1; ++__i)
288 __os << __x._M_x[__i] << __space;
289 __os << __x._M_x[__n - 1];
290
291 __os.flags(__flags);
292 __os.fill(__fill);
293 return __os;
294 }
295
296 template<class _UIntType, int __w, int __n, int __m, int __r,
297 _UIntType __a, int __u, int __s, _UIntType __b, int __t,
298 _UIntType __c, int __l,
299 typename _CharT, typename _Traits>
300 std::basic_istream<_CharT, _Traits>&
301 operator>>(std::basic_istream<_CharT, _Traits>& __is,
302 mersenne_twister<_UIntType, __w, __n, __m,
303 __r, __a, __u, __s, __b, __t, __c, __l>& __x)
304 {
305 typedef std::basic_istream<_CharT, _Traits> __istream_type;
306 typedef typename __istream_type::ios_base __ios_base;
307
308 const typename __ios_base::fmtflags __flags = __is.flags();
309 __is.flags(__ios_base::dec | __ios_base::skipws);
310
311 for (int __i = 0; __i < __n; ++__i)
312 __is >> __x._M_x[__i];
313
314 __is.flags(__flags);
315 return __is;
316 }
317
318
319 template<typename _IntType, _IntType __m, int __s, int __r>
320 void
321 subtract_with_carry<_IntType, __m, __s, __r>::
322 seed(unsigned long __value)
323 {
324 if (__value == 0)
325 __value = 19780503;
326
327 std::_GLIBCXX_TR1 linear_congruential<unsigned long, 40014, 0, 2147483563>
328 __lcg(__value);
329
330 for (int __i = 0; __i < long_lag; ++__i)
331 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());
332
333 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
334 _M_p = 0;
335 }
336
337 template<typename _IntType, _IntType __m, int __s, int __r>
338 template<class _Gen>
339 void
340 subtract_with_carry<_IntType, __m, __s, __r>::
341 seed(_Gen& __gen, false_type)
342 {
343 const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;
344
345 for (int __i = 0; __i < long_lag; ++__i)
346 {
347 _UIntType __tmp = 0;
348 _UIntType __factor = 1;
349 for (int __j = 0; __j < __n; ++__j)
350 {
351 __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
352 (__gen()) * __factor;
353 __factor *= __detail::_Shift<_UIntType, 32>::__value;
354 }
355 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
356 }
357 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
358 _M_p = 0;
359 }
360
361 template<typename _IntType, _IntType __m, int __s, int __r>
362 typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
363 subtract_with_carry<_IntType, __m, __s, __r>::
364 operator()()
365 {
366 // Derive short lag index from current index.
367 int __ps = _M_p - short_lag;
368 if (__ps < 0)
369 __ps += long_lag;
370
371 // Calculate new x(i) without overflow or division.
372 // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
373 // cannot overflow.
374 _UIntType __xi;
375 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
376 {
377 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
378 _M_carry = 0;
379 }
380 else
381 {
382 __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
383 _M_carry = 1;
384 }
385 _M_x[_M_p] = __xi;
386
387 // Adjust current index to loop around in ring buffer.
388 if (++_M_p >= long_lag)
389 _M_p = 0;
390
391 return __xi;
392 }
393
394 template<typename _IntType, _IntType __m, int __s, int __r,
395 typename _CharT, typename _Traits>
396 std::basic_ostream<_CharT, _Traits>&
397 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
398 const subtract_with_carry<_IntType, __m, __s, __r>& __x)
399 {
400 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
401 typedef typename __ostream_type::ios_base __ios_base;
402
403 const typename __ios_base::fmtflags __flags = __os.flags();
404 const _CharT __fill = __os.fill();
405 const _CharT __space = __os.widen(' ');
406 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
407 __os.fill(__space);
408
409 for (int __i = 0; __i < __r; ++__i)
410 __os << __x._M_x[__i] << __space;
411 __os << __x._M_carry;
412
413 __os.flags(__flags);
414 __os.fill(__fill);
415 return __os;
416 }
417
418 template<typename _IntType, _IntType __m, int __s, int __r,
419 typename _CharT, typename _Traits>
420 std::basic_istream<_CharT, _Traits>&
421 operator>>(std::basic_istream<_CharT, _Traits>& __is,
422 subtract_with_carry<_IntType, __m, __s, __r>& __x)
423 {
424 typedef std::basic_ostream<_CharT, _Traits> __istream_type;
425 typedef typename __istream_type::ios_base __ios_base;
426
427 const typename __ios_base::fmtflags __flags = __is.flags();
428 __is.flags(__ios_base::dec | __ios_base::skipws);
429
430 for (int __i = 0; __i < __r; ++__i)
431 __is >> __x._M_x[__i];
432 __is >> __x._M_carry;
433
434 __is.flags(__flags);
435 return __is;
436 }
437
438
439 template<typename _RealType, int __w, int __s, int __r>
440 void
441 subtract_with_carry_01<_RealType, __w, __s, __r>::
442 _M_initialize_npows()
443 {
444 for (int __j = 0; __j < __n; ++__j)
445#if _GLIBCXX_USE_C99_MATH_TR1
446 _M_npows[__j] = std::_GLIBCXX_TR1 ldexp(_RealType(1), -__w + __j * 32);
447#else
448 _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
449#endif
450 }
451
452 template<typename _RealType, int __w, int __s, int __r>
453 void
454 subtract_with_carry_01<_RealType, __w, __s, __r>::
455 seed(unsigned long __value)
456 {
457 if (__value == 0)
458 __value = 19780503;
459
460 // _GLIBCXX_RESOLVE_LIB_DEFECTS
461 // 512. Seeding subtract_with_carry_01 from a single unsigned long.
462 std::_GLIBCXX_TR1 linear_congruential<unsigned long, 40014, 0, 2147483563>
463 __lcg(__value);
464
465 this->seed(__lcg);
466 }
467
468 template<typename _RealType, int __w, int __s, int __r>
469 template<class _Gen>
470 void
471 subtract_with_carry_01<_RealType, __w, __s, __r>::
472 seed(_Gen& __gen, false_type)
473 {
474 for (int __i = 0; __i < long_lag; ++__i)
475 {
476 for (int __j = 0; __j < __n - 1; ++__j)
477 _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
478 _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
479 __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
480 }
481
482 _M_carry = 1;
483 for (int __j = 0; __j < __n; ++__j)
484 if (_M_x[long_lag - 1][__j] != 0)
485 {
486 _M_carry = 0;
487 break;
488 }
489
490 _M_p = 0;
491 }
492
493 template<typename _RealType, int __w, int __s, int __r>
494 typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
495 subtract_with_carry_01<_RealType, __w, __s, __r>::
496 operator()()
497 {
498 // Derive short lag index from current index.
499 int __ps = _M_p - short_lag;
500 if (__ps < 0)
501 __ps += long_lag;
502
503 _UInt32Type __new_carry;
504 for (int __j = 0; __j < __n - 1; ++__j)
505 {
506 if (_M_x[__ps][__j] > _M_x[_M_p][__j]
507 || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
508 __new_carry = 0;
509 else
510 __new_carry = 1;
511
512 _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
513 _M_carry = __new_carry;
514 }
515
516 if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
517 || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
518 __new_carry = 0;
519 else
520 __new_carry = 1;
521
522 _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
523 __detail::_Shift<_UInt32Type, __w % 32>::__value>
524 (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
525 _M_carry = __new_carry;
526
527 result_type __ret = 0.0;
528 for (int __j = 0; __j < __n; ++__j)
529 __ret += _M_x[_M_p][__j] * _M_npows[__j];
530
531 // Adjust current index to loop around in ring buffer.
532 if (++_M_p >= long_lag)
533 _M_p = 0;
534
535 return __ret;
536 }
537
538 template<typename _RealType, int __w, int __s, int __r,
539 typename _CharT, typename _Traits>
540 std::basic_ostream<_CharT, _Traits>&
541 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
542 const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
543 {
544 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
545 typedef typename __ostream_type::ios_base __ios_base;
546
547 const typename __ios_base::fmtflags __flags = __os.flags();
548 const _CharT __fill = __os.fill();
549 const _CharT __space = __os.widen(' ');
550 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
551 __os.fill(__space);
552
553 for (int __i = 0; __i < __r; ++__i)
554 for (int __j = 0; __j < __x.__n; ++__j)
555 __os << __x._M_x[__i][__j] << __space;
556 __os << __x._M_carry;
557
558 __os.flags(__flags);
559 __os.fill(__fill);
560 return __os;
561 }
562
563 template<typename _RealType, int __w, int __s, int __r,
564 typename _CharT, typename _Traits>
565 std::basic_istream<_CharT, _Traits>&
566 operator>>(std::basic_istream<_CharT, _Traits>& __is,
567 subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
568 {
569 typedef std::basic_istream<_CharT, _Traits> __istream_type;
570 typedef typename __istream_type::ios_base __ios_base;
571
572 const typename __ios_base::fmtflags __flags = __is.flags();
573 __is.flags(__ios_base::dec | __ios_base::skipws);
574
575 for (int __i = 0; __i < __r; ++__i)
576 for (int __j = 0; __j < __x.__n; ++__j)
577 __is >> __x._M_x[__i][__j];
578 __is >> __x._M_carry;
579
580 __is.flags(__flags);
581 return __is;
582 }
583
584
585 template<class _UniformRandomNumberGenerator, int __p, int __r>
586 typename discard_block<_UniformRandomNumberGenerator,
587 __p, __r>::result_type
588 discard_block<_UniformRandomNumberGenerator, __p, __r>::
589 operator()()
590 {
591 if (_M_n >= used_block)
592 {
593 while (_M_n < block_size)
594 {
595 _M_b();
596 ++_M_n;
597 }
598 _M_n = 0;
599 }
600 ++_M_n;
601 return _M_b();
602 }
603
604 template<class _UniformRandomNumberGenerator, int __p, int __r,
605 typename _CharT, typename _Traits>
606 std::basic_ostream<_CharT, _Traits>&
607 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
608 const discard_block<_UniformRandomNumberGenerator,
609 __p, __r>& __x)
610 {
611 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
612 typedef typename __ostream_type::ios_base __ios_base;
613
614 const typename __ios_base::fmtflags __flags = __os.flags();
615 const _CharT __fill = __os.fill();
616 const _CharT __space = __os.widen(' ');
617 __os.flags(__ios_base::dec | __ios_base::fixed
618 | __ios_base::left);
619 __os.fill(__space);
620
621 __os << __x._M_b << __space << __x._M_n;
622
623 __os.flags(__flags);
624 __os.fill(__fill);
625 return __os;
626 }
627
628 template<class _UniformRandomNumberGenerator, int __p, int __r,
629 typename _CharT, typename _Traits>
630 std::basic_istream<_CharT, _Traits>&
631 operator>>(std::basic_istream<_CharT, _Traits>& __is,
632 discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
633 {
634 typedef std::basic_istream<_CharT, _Traits> __istream_type;
635 typedef typename __istream_type::ios_base __ios_base;
636
637 const typename __ios_base::fmtflags __flags = __is.flags();
638 __is.flags(__ios_base::dec | __ios_base::skipws);
639
640 __is >> __x._M_b >> __x._M_n;
641
642 __is.flags(__flags);
643 return __is;
644 }
645
646
647 template<class _UniformRandomNumberGenerator1, int __s1,
648 class _UniformRandomNumberGenerator2, int __s2>
649 void
650 xor_combine<_UniformRandomNumberGenerator1, __s1,
651 _UniformRandomNumberGenerator2, __s2>::
652 _M_initialize_max()
653 {
654 const int __w = std::numeric_limits<result_type>::digits;
655
656 const result_type __m1 =
657 std::min(result_type(_M_b1.max() - _M_b1.min()),
658 __detail::_Shift<result_type, __w - __s1>::__value - 1);
659
660 const result_type __m2 =
661 std::min(result_type(_M_b2.max() - _M_b2.min()),
662 __detail::_Shift<result_type, __w - __s2>::__value - 1);
663
664 // NB: In TR1 s1 is not required to be >= s2.
665 if (__s1 < __s2)
666 _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
667 else
668 _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
669 }
670
671 template<class _UniformRandomNumberGenerator1, int __s1,
672 class _UniformRandomNumberGenerator2, int __s2>
673 typename xor_combine<_UniformRandomNumberGenerator1, __s1,
674 _UniformRandomNumberGenerator2, __s2>::result_type
675 xor_combine<_UniformRandomNumberGenerator1, __s1,
676 _UniformRandomNumberGenerator2, __s2>::
677 _M_initialize_max_aux(result_type __a, result_type __b, int __d)
678 {
679 const result_type __two2d = result_type(1) << __d;
680 const result_type __c = __a * __two2d;
681
682 if (__a == 0 || __b < __two2d)
683 return __c + __b;
684
685 const result_type __t = std::max(__c, __b);
686 const result_type __u = std::min(__c, __b);
687
688 result_type __ub = __u;
689 result_type __p;
690 for (__p = 0; __ub != 1; __ub >>= 1)
691 ++__p;
692
693 const result_type __two2p = result_type(1) << __p;
694 const result_type __k = __t / __two2p;
695
696 if (__k & 1)
697 return (__k + 1) * __two2p - 1;
698
699 if (__c >= __b)
700 return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
701 / __two2d,
702 __u % __two2p, __d);
703 else
704 return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
705 / __two2d,
706 __t % __two2p, __d);
707 }
708
709 template<class _UniformRandomNumberGenerator1, int __s1,
710 class _UniformRandomNumberGenerator2, int __s2,
711 typename _CharT, typename _Traits>
712 std::basic_ostream<_CharT, _Traits>&
713 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
714 const xor_combine<_UniformRandomNumberGenerator1, __s1,
715 _UniformRandomNumberGenerator2, __s2>& __x)
716 {
717 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
718 typedef typename __ostream_type::ios_base __ios_base;
719
720 const typename __ios_base::fmtflags __flags = __os.flags();
721 const _CharT __fill = __os.fill();
722 const _CharT __space = __os.widen(' ');
723 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
724 __os.fill(__space);
725
726 __os << __x.base1() << __space << __x.base2();
727
728 __os.flags(__flags);
729 __os.fill(__fill);
730 return __os;
731 }
732
733 template<class _UniformRandomNumberGenerator1, int __s1,
734 class _UniformRandomNumberGenerator2, int __s2,
735 typename _CharT, typename _Traits>
736 std::basic_istream<_CharT, _Traits>&
737 operator>>(std::basic_istream<_CharT, _Traits>& __is,
738 xor_combine<_UniformRandomNumberGenerator1, __s1,
739 _UniformRandomNumberGenerator2, __s2>& __x)
740 {
741 typedef std::basic_istream<_CharT, _Traits> __istream_type;
742 typedef typename __istream_type::ios_base __ios_base;
743
744 const typename __ios_base::fmtflags __flags = __is.flags();
745 __is.flags(__ios_base::skipws);
746
747 __is >> __x._M_b1 >> __x._M_b2;
748
749 __is.flags(__flags);
750 return __is;
751 }
752
753
754 template<typename _IntType>
755 template<typename _UniformRandomNumberGenerator>
756 typename uniform_int<_IntType>::result_type
757 uniform_int<_IntType>::
758 _M_call(_UniformRandomNumberGenerator& __urng,
759 result_type __min, result_type __max, true_type)
760 {
761 // XXX Must be fixed to work well for *arbitrary* __urng.max(),
762 // __urng.min(), __max, __min. Currently works fine only in the
763 // most common case __urng.max() - __urng.min() >= __max - __min,
764 // with __urng.max() > __urng.min() >= 0.
765 typedef typename __gnu_cxx::__add_unsigned<typename
766 _UniformRandomNumberGenerator::result_type>::__type __urntype;
767 typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
768 __utype;
769 typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
770 > sizeof(__utype)),
771 __urntype, __utype>::__type __uctype;
772
773 result_type __ret;
774
775 const __urntype __urnmin = __urng.min();
776 const __urntype __urnmax = __urng.max();
777 const __urntype __urnrange = __urnmax - __urnmin;
778 const __uctype __urange = __max - __min;
779 const __uctype __udenom = (__urnrange <= __urange
780 ? 1 : __urnrange / (__urange + 1));
781 do
782 __ret = (__urntype(__urng()) - __urnmin) / __udenom;
783 while (__ret > __max - __min);
784
785 return __ret + __min;
786 }
787
788 template<typename _IntType, typename _CharT, typename _Traits>
789 std::basic_ostream<_CharT, _Traits>&
790 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
791 const uniform_int<_IntType>& __x)
792 {
793 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
794 typedef typename __ostream_type::ios_base __ios_base;
795
796 const typename __ios_base::fmtflags __flags = __os.flags();
797 const _CharT __fill = __os.fill();
798 const _CharT __space = __os.widen(' ');
799 __os.flags(__ios_base::scientific | __ios_base::left);
800 __os.fill(__space);
801
802 __os << __x.min() << __space << __x.max();
803
804 __os.flags(__flags);
805 __os.fill(__fill);
806 return __os;
807 }
808
809 template<typename _IntType, typename _CharT, typename _Traits>
810 std::basic_istream<_CharT, _Traits>&
811 operator>>(std::basic_istream<_CharT, _Traits>& __is,
812 uniform_int<_IntType>& __x)
813 {
814 typedef std::basic_istream<_CharT, _Traits> __istream_type;
815 typedef typename __istream_type::ios_base __ios_base;
816
817 const typename __ios_base::fmtflags __flags = __is.flags();
818 __is.flags(__ios_base::dec | __ios_base::skipws);
819
820 __is >> __x._M_min >> __x._M_max;
821
822 __is.flags(__flags);
823 return __is;
824 }
825
826
827 template<typename _CharT, typename _Traits>
828 std::basic_ostream<_CharT, _Traits>&
829 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
830 const bernoulli_distribution& __x)
831 {
832 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
833 typedef typename __ostream_type::ios_base __ios_base;
834
835 const typename __ios_base::fmtflags __flags = __os.flags();
836 const _CharT __fill = __os.fill();
837 const std::streamsize __precision = __os.precision();
838 __os.flags(__ios_base::scientific | __ios_base::left);
839 __os.fill(__os.widen(' '));
840 __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);
841
842 __os << __x.p();
843
844 __os.flags(__flags);
845 __os.fill(__fill);
846 __os.precision(__precision);
847 return __os;
848 }
849
850
851 template<typename _IntType, typename _RealType>
852 template<class _UniformRandomNumberGenerator>
853 typename geometric_distribution<_IntType, _RealType>::result_type
854 geometric_distribution<_IntType, _RealType>::
855 operator()(_UniformRandomNumberGenerator& __urng)
856 {
857 // About the epsilon thing see this thread:
858 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
859 const _RealType __naf =
860 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
861 // The largest _RealType convertible to _IntType.
862 const _RealType __thr =
863 std::numeric_limits<_IntType>::max() + __naf;
864
865 _RealType __cand;
866 do
867 __cand = std::ceil(std::log(__urng()) / _M_log_p);
868 while (__cand >= __thr);
869
870 return result_type(__cand + __naf);
871 }
872
873 template<typename _IntType, typename _RealType,
874 typename _CharT, typename _Traits>
875 std::basic_ostream<_CharT, _Traits>&
876 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
877 const geometric_distribution<_IntType, _RealType>& __x)
878 {
879 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
880 typedef typename __ostream_type::ios_base __ios_base;
881
882 const typename __ios_base::fmtflags __flags = __os.flags();
883 const _CharT __fill = __os.fill();
884 const std::streamsize __precision = __os.precision();
885 __os.flags(__ios_base::scientific | __ios_base::left);
886 __os.fill(__os.widen(' '));
887 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
888
889 __os << __x.p();
890
891 __os.flags(__flags);
892 __os.fill(__fill);
893 __os.precision(__precision);
894 return __os;
895 }
896
897
898 template<typename _IntType, typename _RealType>
899 void
900 poisson_distribution<_IntType, _RealType>::
901 _M_initialize()
902 {
903#if _GLIBCXX_USE_C99_MATH_TR1
904 if (_M_mean >= 12)
905 {
906 const _RealType __m = std::floor(_M_mean);
907 _M_lm_thr = std::log(_M_mean);
908 _M_lfm = std::_GLIBCXX_TR1 lgamma(__m + 1);
909 _M_sm = std::sqrt(__m);
910
911 const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
912 const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
913 / __pi_4));
914 _M_d = std::_GLIBCXX_TR1 round(std::max(_RealType(6),
915 std::min(__m, __dx)));
916 const _RealType __cx = 2 * __m + _M_d;
917 _M_scx = std::sqrt(__cx / 2);
918 _M_1cx = 1 / __cx;
919
920 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
921 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
922 }
923 else
924#endif
925 _M_lm_thr = std::exp(-_M_mean);
926 }
927
928 /**
929 * A rejection algorithm when mean >= 12 and a simple method based
930 * upon the multiplication of uniform random variates otherwise.
931 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
932 * is defined.
933 *
934 * Reference:
935 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
936 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
937 */
938 template<typename _IntType, typename _RealType>
939 template<class _UniformRandomNumberGenerator>
940 typename poisson_distribution<_IntType, _RealType>::result_type
941 poisson_distribution<_IntType, _RealType>::
942 operator()(_UniformRandomNumberGenerator& __urng)
943 {
944#if _GLIBCXX_USE_C99_MATH_TR1
945 if (_M_mean >= 12)
946 {
947 _RealType __x;
948
949 // See comments above...
950 const _RealType __naf =
951 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
952 const _RealType __thr =
953 std::numeric_limits<_IntType>::max() + __naf;
954
955 const _RealType __m = std::floor(_M_mean);
956 // sqrt(pi / 2)
957 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
958 const _RealType __c1 = _M_sm * __spi_2;
959 const _RealType __c2 = _M_c2b + __c1;
960 const _RealType __c3 = __c2 + 1;
961 const _RealType __c4 = __c3 + 1;
962 // e^(1 / 78)
963 const _RealType __e178 = 1.0129030479320018583185514777512983L;
964 const _RealType __c5 = __c4 + __e178;
965 const _RealType __c = _M_cb + __c5;
966 const _RealType __2cx = 2 * (2 * __m + _M_d);
967
968 bool __reject = true;
969 do
970 {
971 const _RealType __u = __c * __urng();
972 const _RealType __e = -std::log(__urng());
973
974 _RealType __w = 0.0;
975
976 if (__u <= __c1)
977 {
978 const _RealType __n = _M_nd(__urng);
979 const _RealType __y = -std::abs(__n) * _M_sm - 1;
980 __x = std::floor(__y);
981 __w = -__n * __n / 2;
982 if (__x < -__m)
983 continue;
984 }
985 else if (__u <= __c2)
986 {
987 const _RealType __n = _M_nd(__urng);
988 const _RealType __y = 1 + std::abs(__n) * _M_scx;
989 __x = std::ceil(__y);
990 __w = __y * (2 - __y) * _M_1cx;
991 if (__x > _M_d)
992 continue;
993 }
994 else if (__u <= __c3)
995 // NB: This case not in the book, nor in the Errata,
996 // but should be ok...
997 __x = -1;
998 else if (__u <= __c4)
999 __x = 0;
1000 else if (__u <= __c5)
1001 __x = 1;
1002 else
1003 {
1004 const _RealType __v = -std::log(__urng());
1005 const _RealType __y = _M_d + __v * __2cx / _M_d;
1006 __x = std::ceil(__y);
1007 __w = -_M_d * _M_1cx * (1 + __y / 2);
1008 }
1009
1010 __reject = (__w - __e - __x * _M_lm_thr
1011 > _M_lfm - std::_GLIBCXX_TR1 lgamma(__x + __m + 1));
1012
1013 __reject |= __x + __m >= __thr;
1014
1015 } while (__reject);
1016
1017 return result_type(__x + __m + __naf);
1018 }
1019 else
1020#endif
1021 {
1022 _IntType __x = 0;
1023 _RealType __prod = 1.0;
1024
1025 do
1026 {
1027 __prod *= __urng();
1028 __x += 1;
1029 }
1030 while (__prod > _M_lm_thr);
1031
1032 return __x - 1;
1033 }
1034 }
1035
1036 template<typename _IntType, typename _RealType,
1037 typename _CharT, typename _Traits>
1038 std::basic_ostream<_CharT, _Traits>&
1039 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1040 const poisson_distribution<_IntType, _RealType>& __x)
1041 {
1042 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1043 typedef typename __ostream_type::ios_base __ios_base;
1044
1045 const typename __ios_base::fmtflags __flags = __os.flags();
1046 const _CharT __fill = __os.fill();
1047 const std::streamsize __precision = __os.precision();
1048 const _CharT __space = __os.widen(' ');
1049 __os.flags(__ios_base::scientific | __ios_base::left);
1050 __os.fill(__space);
1051 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1052
1053 __os << __x.mean() << __space << __x._M_nd;
1054
1055 __os.flags(__flags);
1056 __os.fill(__fill);
1057 __os.precision(__precision);
1058 return __os;
1059 }
1060
1061 template<typename _IntType, typename _RealType,
1062 typename _CharT, typename _Traits>
1063 std::basic_istream<_CharT, _Traits>&
1064 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1065 poisson_distribution<_IntType, _RealType>& __x)
1066 {
1067 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1068 typedef typename __istream_type::ios_base __ios_base;
1069
1070 const typename __ios_base::fmtflags __flags = __is.flags();
1071 __is.flags(__ios_base::skipws);
1072
1073 __is >> __x._M_mean >> __x._M_nd;
1074 __x._M_initialize();
1075
1076 __is.flags(__flags);
1077 return __is;
1078 }
1079
1080
1081 template<typename _IntType, typename _RealType>
1082 void
1083 binomial_distribution<_IntType, _RealType>::
1084 _M_initialize()
1085 {
1086 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1087
1088 _M_easy = true;
1089
1090#if _GLIBCXX_USE_C99_MATH_TR1
1091 if (_M_t * __p12 >= 8)
1092 {
1093 _M_easy = false;
1094 const _RealType __np = std::floor(_M_t * __p12);
1095 const _RealType __pa = __np / _M_t;
1096 const _RealType __1p = 1 - __pa;
1097
1098 const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
1099 const _RealType __d1x =
1100 std::sqrt(__np * __1p * std::log(32 * __np
1101 / (81 * __pi_4 * __1p)));
1102 _M_d1 = std::_GLIBCXX_TR1 round(std::max(_RealType(1), __d1x));
1103 const _RealType __d2x =
1104 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1105 / (__pi_4 * __pa)));
1106 _M_d2 = std::_GLIBCXX_TR1 round(std::max(_RealType(1), __d2x));
1107
1108 // sqrt(pi / 2)
1109 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1110 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1111 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1112 _M_c = 2 * _M_d1 / __np;
1113 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1114 const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
1115 const _RealType __s1s = _M_s1 * _M_s1;
1116 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1117 * 2 * __s1s / _M_d1
1118 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1119 const _RealType __s2s = _M_s2 * _M_s2;
1120 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1121 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1122 _M_lf = (std::_GLIBCXX_TR1 lgamma(__np + 1)
1123 + std::_GLIBCXX_TR1 lgamma(_M_t - __np + 1));
1124 _M_lp1p = std::log(__pa / __1p);
1125
1126 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1127 }
1128 else
1129#endif
1130 _M_q = -std::log(1 - __p12);
1131 }
1132
1133 template<typename _IntType, typename _RealType>
1134 template<class _UniformRandomNumberGenerator>
1135 typename binomial_distribution<_IntType, _RealType>::result_type
1136 binomial_distribution<_IntType, _RealType>::
1137 _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1138 {
1139 _IntType __x = 0;
1140 _RealType __sum = 0;
1141
1142 do
1143 {
1144 const _RealType __e = -std::log(__urng());
1145 __sum += __e / (__t - __x);
1146 __x += 1;
1147 }
1148 while (__sum <= _M_q);
1149
1150 return __x - 1;
1151 }
1152
1153 /**
1154 * A rejection algorithm when t * p >= 8 and a simple waiting time
1155 * method - the second in the referenced book - otherwise.
1156 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1157 * is defined.
1158 *
1159 * Reference:
1160 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1161 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1162 */
1163 template<typename _IntType, typename _RealType>
1164 template<class _UniformRandomNumberGenerator>
1165 typename binomial_distribution<_IntType, _RealType>::result_type
1166 binomial_distribution<_IntType, _RealType>::
1167 operator()(_UniformRandomNumberGenerator& __urng)
1168 {
1169 result_type __ret;
1170 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1171
1172#if _GLIBCXX_USE_C99_MATH_TR1
1173 if (!_M_easy)
1174 {
1175 _RealType __x;
1176
1177 // See comments above...
1178 const _RealType __naf =
1179 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
1180 const _RealType __thr =
1181 std::numeric_limits<_IntType>::max() + __naf;
1182
1183 const _RealType __np = std::floor(_M_t * __p12);
1184 const _RealType __pa = __np / _M_t;
1185
1186 // sqrt(pi / 2)
1187 const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1188 const _RealType __a1 = _M_a1;
1189 const _RealType __a12 = __a1 + _M_s2 * __spi_2;
1190 const _RealType __a123 = _M_a123;
1191 const _RealType __s1s = _M_s1 * _M_s1;
1192 const _RealType __s2s = _M_s2 * _M_s2;
1193
1194 bool __reject;
1195 do
1196 {
1197 const _RealType __u = _M_s * __urng();
1198
1199 _RealType __v;
1200
1201 if (__u <= __a1)
1202 {
1203 const _RealType __n = _M_nd(__urng);
1204 const _RealType __y = _M_s1 * std::abs(__n);
1205 __reject = __y >= _M_d1;
1206 if (!__reject)
1207 {
1208 const _RealType __e = -std::log(__urng());
1209 __x = std::floor(__y);
1210 __v = -__e - __n * __n / 2 + _M_c;
1211 }
1212 }
1213 else if (__u <= __a12)
1214 {
1215 const _RealType __n = _M_nd(__urng);
1216 const _RealType __y = _M_s2 * std::abs(__n);
1217 __reject = __y >= _M_d2;
1218 if (!__reject)
1219 {
1220 const _RealType __e = -std::log(__urng());
1221 __x = std::floor(-__y);
1222 __v = -__e - __n * __n / 2;
1223 }
1224 }
1225 else if (__u <= __a123)
1226 {
1227 const _RealType __e1 = -std::log(__urng());
1228 const _RealType __e2 = -std::log(__urng());
1229
1230 const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
1231 __x = std::floor(__y);
1232 __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
1233 -__y / (2 * __s1s)));
1234 __reject = false;
1235 }
1236 else
1237 {
1238 const _RealType __e1 = -std::log(__urng());
1239 const _RealType __e2 = -std::log(__urng());
1240
1241 const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
1242 __x = std::floor(-__y);
1243 __v = -__e2 - _M_d2 * __y / (2 * __s2s);
1244 __reject = false;
1245 }
1246
1247 __reject = __reject || __x < -__np || __x > _M_t - __np;
1248 if (!__reject)
1249 {
1250 const _RealType __lfx =
1251 std::_GLIBCXX_TR1 lgamma(__np + __x + 1)
1252 + std::_GLIBCXX_TR1 lgamma(_M_t - (__np + __x) + 1);
1253 __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
1254 }
1255
1256 __reject |= __x + __np >= __thr;
1257 }
1258 while (__reject);
1259
1260 __x += __np + __naf;
1261
1262 const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x));
1263 __ret = _IntType(__x) + __z;
1264 }
1265 else
1266#endif
1267 __ret = _M_waiting(__urng, _M_t);
1268
1269 if (__p12 != _M_p)
1270 __ret = _M_t - __ret;
1271 return __ret;
1272 }
1273
1274 template<typename _IntType, typename _RealType,
1275 typename _CharT, typename _Traits>
1276 std::basic_ostream<_CharT, _Traits>&
1277 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1278 const binomial_distribution<_IntType, _RealType>& __x)
1279 {
1280 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1281 typedef typename __ostream_type::ios_base __ios_base;
1282
1283 const typename __ios_base::fmtflags __flags = __os.flags();
1284 const _CharT __fill = __os.fill();
1285 const std::streamsize __precision = __os.precision();
1286 const _CharT __space = __os.widen(' ');
1287 __os.flags(__ios_base::scientific | __ios_base::left);
1288 __os.fill(__space);
1289 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1290
1291 __os << __x.t() << __space << __x.p()
1292 << __space << __x._M_nd;
1293
1294 __os.flags(__flags);
1295 __os.fill(__fill);
1296 __os.precision(__precision);
1297 return __os;
1298 }
1299
1300 template<typename _IntType, typename _RealType,
1301 typename _CharT, typename _Traits>
1302 std::basic_istream<_CharT, _Traits>&
1303 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1304 binomial_distribution<_IntType, _RealType>& __x)
1305 {
1306 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1307 typedef typename __istream_type::ios_base __ios_base;
1308
1309 const typename __ios_base::fmtflags __flags = __is.flags();
1310 __is.flags(__ios_base::dec | __ios_base::skipws);
1311
1312 __is >> __x._M_t >> __x._M_p >> __x._M_nd;
1313 __x._M_initialize();
1314
1315 __is.flags(__flags);
1316 return __is;
1317 }
1318
1319
1320 template<typename _RealType, typename _CharT, typename _Traits>
1321 std::basic_ostream<_CharT, _Traits>&
1322 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1323 const uniform_real<_RealType>& __x)
1324 {
1325 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1326 typedef typename __ostream_type::ios_base __ios_base;
1327
1328 const typename __ios_base::fmtflags __flags = __os.flags();
1329 const _CharT __fill = __os.fill();
1330 const std::streamsize __precision = __os.precision();
1331 const _CharT __space = __os.widen(' ');
1332 __os.flags(__ios_base::scientific | __ios_base::left);
1333 __os.fill(__space);
1334 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1335
1336 __os << __x.min() << __space << __x.max();
1337
1338 __os.flags(__flags);
1339 __os.fill(__fill);
1340 __os.precision(__precision);
1341 return __os;
1342 }
1343
1344 template<typename _RealType, typename _CharT, typename _Traits>
1345 std::basic_istream<_CharT, _Traits>&
1346 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1347 uniform_real<_RealType>& __x)
1348 {
1349 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1350 typedef typename __istream_type::ios_base __ios_base;
1351
1352 const typename __ios_base::fmtflags __flags = __is.flags();
1353 __is.flags(__ios_base::skipws);
1354
1355 __is >> __x._M_min >> __x._M_max;
1356
1357 __is.flags(__flags);
1358 return __is;
1359 }
1360
1361
1362 template<typename _RealType, typename _CharT, typename _Traits>
1363 std::basic_ostream<_CharT, _Traits>&
1364 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1365 const exponential_distribution<_RealType>& __x)
1366 {
1367 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1368 typedef typename __ostream_type::ios_base __ios_base;
1369
1370 const typename __ios_base::fmtflags __flags = __os.flags();
1371 const _CharT __fill = __os.fill();
1372 const std::streamsize __precision = __os.precision();
1373 __os.flags(__ios_base::scientific | __ios_base::left);
1374 __os.fill(__os.widen(' '));
1375 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1376
1377 __os << __x.lambda();
1378
1379 __os.flags(__flags);
1380 __os.fill(__fill);
1381 __os.precision(__precision);
1382 return __os;
1383 }
1384
1385
1386 /**
1387 * Polar method due to Marsaglia.
1388 *
1389 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1390 * New York, 1986, Ch. V, Sect. 4.4.
1391 */
1392 template<typename _RealType>
1393 template<class _UniformRandomNumberGenerator>
1394 typename normal_distribution<_RealType>::result_type
1395 normal_distribution<_RealType>::
1396 operator()(_UniformRandomNumberGenerator& __urng)
1397 {
1398 result_type __ret;
1399
1400 if (_M_saved_available)
1401 {
1402 _M_saved_available = false;
1403 __ret = _M_saved;
1404 }
1405 else
1406 {
1407 result_type __x, __y, __r2;
1408 do
1409 {
1410 __x = result_type(2.0) * __urng() - 1.0;
1411 __y = result_type(2.0) * __urng() - 1.0;
1412 __r2 = __x * __x + __y * __y;
1413 }
1414 while (__r2 > 1.0 || __r2 == 0.0);
1415
1416 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1417 _M_saved = __x * __mult;
1418 _M_saved_available = true;
1419 __ret = __y * __mult;
1420 }
1421
1422 __ret = __ret * _M_sigma + _M_mean;
1423 return __ret;
1424 }
1425
1426 template<typename _RealType, typename _CharT, typename _Traits>
1427 std::basic_ostream<_CharT, _Traits>&
1428 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1429 const normal_distribution<_RealType>& __x)
1430 {
1431 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1432 typedef typename __ostream_type::ios_base __ios_base;
1433
1434 const typename __ios_base::fmtflags __flags = __os.flags();
1435 const _CharT __fill = __os.fill();
1436 const std::streamsize __precision = __os.precision();
1437 const _CharT __space = __os.widen(' ');
1438 __os.flags(__ios_base::scientific | __ios_base::left);
1439 __os.fill(__space);
1440 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1441
1442 __os << __x._M_saved_available << __space
1443 << __x.mean() << __space
1444 << __x.sigma();
1445 if (__x._M_saved_available)
1446 __os << __space << __x._M_saved;
1447
1448 __os.flags(__flags);
1449 __os.fill(__fill);
1450 __os.precision(__precision);
1451 return __os;
1452 }
1453
1454 template<typename _RealType, typename _CharT, typename _Traits>
1455 std::basic_istream<_CharT, _Traits>&
1456 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1457 normal_distribution<_RealType>& __x)
1458 {
1459 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1460 typedef typename __istream_type::ios_base __ios_base;
1461
1462 const typename __ios_base::fmtflags __flags = __is.flags();
1463 __is.flags(__ios_base::dec | __ios_base::skipws);
1464
1465 __is >> __x._M_saved_available >> __x._M_mean
1466 >> __x._M_sigma;
1467 if (__x._M_saved_available)
1468 __is >> __x._M_saved;
1469
1470 __is.flags(__flags);
1471 return __is;
1472 }
1473
1474
1475 template<typename _RealType>
1476 void
1477 gamma_distribution<_RealType>::
1478 _M_initialize()
1479 {
1480 if (_M_alpha >= 1)
1481 _M_l_d = std::sqrt(2 * _M_alpha - 1);
1482 else
1483 _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
1484 * (1 - _M_alpha));
1485 }
1486
1487 /**
1488 * Cheng's rejection algorithm GB for alpha >= 1 and a modification
1489 * of Vaduva's rejection from Weibull algorithm due to Devroye for
1490 * alpha < 1.
1491 *
1492 * References:
1493 * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
1494 * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
1495 *
1496 * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
1497 * and Composition Procedures." Math. Operationsforschung and Statistik,
1498 * Series in Statistics, 8, 545-576, 1977.
1499 *
1500 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1501 * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
1502 */
1503 template<typename _RealType>
1504 template<class _UniformRandomNumberGenerator>
1505 typename gamma_distribution<_RealType>::result_type
1506 gamma_distribution<_RealType>::
1507 operator()(_UniformRandomNumberGenerator& __urng)
1508 {
1509 result_type __x;
1510
1511 bool __reject;
1512 if (_M_alpha >= 1)
1513 {
1514 // alpha - log(4)
1515 const result_type __b = _M_alpha
1516 - result_type(1.3862943611198906188344642429163531L);
1517 const result_type __c = _M_alpha + _M_l_d;
1518 const result_type __1l = 1 / _M_l_d;
1519
1520 // 1 + log(9 / 2)
1521 const result_type __k = 2.5040773967762740733732583523868748L;
1522
1523 do
1524 {
1525 const result_type __u = __urng();
1526 const result_type __v = __urng();
1527
1528 const result_type __y = __1l * std::log(__v / (1 - __v));
1529 __x = _M_alpha * std::exp(__y);
1530
1531 const result_type __z = __u * __v * __v;
1532 const result_type __r = __b + __c * __y - __x;
1533
1534 __reject = __r < result_type(4.5) * __z - __k;
1535 if (__reject)
1536 __reject = __r < std::log(__z);
1537 }
1538 while (__reject);
1539 }
1540 else
1541 {
1542 const result_type __c = 1 / _M_alpha;
1543
1544 do
1545 {
1546 const result_type __z = -std::log(__urng());
1547 const result_type __e = -std::log(__urng());
1548
1549 __x = std::pow(__z, __c);
1550
1551 __reject = __z + __e < _M_l_d + __x;
1552 }
1553 while (__reject);
1554 }
1555
1556 return __x;
1557 }
1558
1559 template<typename _RealType, typename _CharT, typename _Traits>
1560 std::basic_ostream<_CharT, _Traits>&
1561 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1562 const gamma_distribution<_RealType>& __x)
1563 {
1564 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1565 typedef typename __ostream_type::ios_base __ios_base;
1566
1567 const typename __ios_base::fmtflags __flags = __os.flags();
1568 const _CharT __fill = __os.fill();
1569 const std::streamsize __precision = __os.precision();
1570 __os.flags(__ios_base::scientific | __ios_base::left);
1571 __os.fill(__os.widen(' '));
1572 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1573
1574 __os << __x.alpha();
1575
1576 __os.flags(__flags);
1577 __os.fill(__fill);
1578 __os.precision(__precision);
1579 return __os;
1580 }
1581
1582}
1583}
Mirror of https://gcc.gnu.org/git/gcc.git