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