source: Daodan/MSYS2/mingw32/include/c++/11.2.0/tr1/random.tcc@ 1177

Last change on this file since 1177 was 1166, checked in by rossy, 3 years ago

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