1 // random number generation (out of line) -*- C++ -*-
3 // Copyright (C) 2009-2018 Free Software Foundation, Inc.
5 // This file is part of the GNU ISO C++ Library. This library is free
6 // software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the
8 // Free Software Foundation; either version 3, or (at your option)
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
16 // 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.
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/>.
25 /** @file bits/random.tcc
26 * This is an internal header file, included by other library headers.
27 * Do not attempt to use it directly. @headername{random}
33 #include <numeric> // std::accumulate and std::partial_sum
35 namespace std _GLIBCXX_VISIBILITY(default)
37 _GLIBCXX_BEGIN_NAMESPACE_VERSION
40 * (Further) implementation-space details.
44 // General case for x = (ax + c) mod m -- use Schrage's algorithm
45 // to avoid integer overflow.
47 // Preconditions: a > 0, m > 0.
49 // Note: only works correctly for __m % __a < __m / __a.
50 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
52 _Mod<_Tp, __m, __a, __c, false, true>::
59 static const _Tp __q = __m / __a;
60 static const _Tp __r = __m % __a;
62 _Tp __t1 = __a * (__x % __q);
63 _Tp __t2 = __r * (__x / __q);
67 __x = __m - __t2 + __t1;
72 const _Tp __d = __m - __x;
81 template<typename _InputIterator, typename _OutputIterator,
84 __normalize(_InputIterator __first, _InputIterator __last,
85 _OutputIterator __result, const _Tp& __factor)
87 for (; __first != __last; ++__first, ++__result)
88 *__result = *__first / __factor;
92 } // namespace __detail
94 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
96 linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
98 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
100 linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
102 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
104 linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
106 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
108 linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
111 * Seeds the LCR with integral value @p __s, adjusted so that the
112 * ring identity is never a member of the convergence set.
114 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
116 linear_congruential_engine<_UIntType, __a, __c, __m>::
117 seed(result_type __s)
119 if ((__detail::__mod<_UIntType, __m>(__c) == 0)
120 && (__detail::__mod<_UIntType, __m>(__s) == 0))
123 _M_x = __detail::__mod<_UIntType, __m>(__s);
127 * Seeds the LCR engine with a value generated by @p __q.
129 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
130 template<typename _Sseq>
131 typename std::enable_if<std::is_class<_Sseq>::value>::type
132 linear_congruential_engine<_UIntType, __a, __c, __m>::
135 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
137 const _UIntType __k = (__k0 + 31) / 32;
138 uint_least32_t __arr[__k + 3];
139 __q.generate(__arr + 0, __arr + __k + 3);
140 _UIntType __factor = 1u;
141 _UIntType __sum = 0u;
142 for (size_t __j = 0; __j < __k; ++__j)
144 __sum += __arr[__j + 3] * __factor;
145 __factor *= __detail::_Shift<_UIntType, 32>::__value;
150 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
151 typename _CharT, typename _Traits>
152 std::basic_ostream<_CharT, _Traits>&
153 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
154 const linear_congruential_engine<_UIntType,
155 __a, __c, __m>& __lcr)
157 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
158 typedef typename __ostream_type::ios_base __ios_base;
160 const typename __ios_base::fmtflags __flags = __os.flags();
161 const _CharT __fill = __os.fill();
162 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
163 __os.fill(__os.widen(' '));
172 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
173 typename _CharT, typename _Traits>
174 std::basic_istream<_CharT, _Traits>&
175 operator>>(std::basic_istream<_CharT, _Traits>& __is,
176 linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
178 typedef std::basic_istream<_CharT, _Traits> __istream_type;
179 typedef typename __istream_type::ios_base __ios_base;
181 const typename __ios_base::fmtflags __flags = __is.flags();
182 __is.flags(__ios_base::dec);
191 template<typename _UIntType,
192 size_t __w, size_t __n, size_t __m, size_t __r,
193 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
194 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
197 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
198 __s, __b, __t, __c, __l, __f>::word_size;
200 template<typename _UIntType,
201 size_t __w, size_t __n, size_t __m, size_t __r,
202 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
203 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
206 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
207 __s, __b, __t, __c, __l, __f>::state_size;
209 template<typename _UIntType,
210 size_t __w, size_t __n, size_t __m, size_t __r,
211 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
212 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
215 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
216 __s, __b, __t, __c, __l, __f>::shift_size;
218 template<typename _UIntType,
219 size_t __w, size_t __n, size_t __m, size_t __r,
220 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
221 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
224 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
225 __s, __b, __t, __c, __l, __f>::mask_bits;
227 template<typename _UIntType,
228 size_t __w, size_t __n, size_t __m, size_t __r,
229 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
230 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
233 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
234 __s, __b, __t, __c, __l, __f>::xor_mask;
236 template<typename _UIntType,
237 size_t __w, size_t __n, size_t __m, size_t __r,
238 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
239 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
242 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
243 __s, __b, __t, __c, __l, __f>::tempering_u;
245 template<typename _UIntType,
246 size_t __w, size_t __n, size_t __m, size_t __r,
247 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
248 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
251 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
252 __s, __b, __t, __c, __l, __f>::tempering_d;
254 template<typename _UIntType,
255 size_t __w, size_t __n, size_t __m, size_t __r,
256 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
257 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
260 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
261 __s, __b, __t, __c, __l, __f>::tempering_s;
263 template<typename _UIntType,
264 size_t __w, size_t __n, size_t __m, size_t __r,
265 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
266 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
269 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
270 __s, __b, __t, __c, __l, __f>::tempering_b;
272 template<typename _UIntType,
273 size_t __w, size_t __n, size_t __m, size_t __r,
274 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
275 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
278 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
279 __s, __b, __t, __c, __l, __f>::tempering_t;
281 template<typename _UIntType,
282 size_t __w, size_t __n, size_t __m, size_t __r,
283 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
284 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
287 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
288 __s, __b, __t, __c, __l, __f>::tempering_c;
290 template<typename _UIntType,
291 size_t __w, size_t __n, size_t __m, size_t __r,
292 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
293 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
296 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
297 __s, __b, __t, __c, __l, __f>::tempering_l;
299 template<typename _UIntType,
300 size_t __w, size_t __n, size_t __m, size_t __r,
301 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
302 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
305 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
306 __s, __b, __t, __c, __l, __f>::
307 initialization_multiplier;
309 template<typename _UIntType,
310 size_t __w, size_t __n, size_t __m, size_t __r,
311 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
312 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
315 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
316 __s, __b, __t, __c, __l, __f>::default_seed;
318 template<typename _UIntType,
319 size_t __w, size_t __n, size_t __m, size_t __r,
320 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
321 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
324 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
325 __s, __b, __t, __c, __l, __f>::
326 seed(result_type __sd)
328 _M_x[0] = __detail::__mod<_UIntType,
329 __detail::_Shift<_UIntType, __w>::__value>(__sd);
331 for (size_t __i = 1; __i < state_size; ++__i)
333 _UIntType __x = _M_x[__i - 1];
334 __x ^= __x >> (__w - 2);
336 __x += __detail::__mod<_UIntType, __n>(__i);
337 _M_x[__i] = __detail::__mod<_UIntType,
338 __detail::_Shift<_UIntType, __w>::__value>(__x);
343 template<typename _UIntType,
344 size_t __w, size_t __n, size_t __m, size_t __r,
345 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
346 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
348 template<typename _Sseq>
349 typename std::enable_if<std::is_class<_Sseq>::value>::type
350 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
351 __s, __b, __t, __c, __l, __f>::
354 const _UIntType __upper_mask = (~_UIntType()) << __r;
355 const size_t __k = (__w + 31) / 32;
356 uint_least32_t __arr[__n * __k];
357 __q.generate(__arr + 0, __arr + __n * __k);
360 for (size_t __i = 0; __i < state_size; ++__i)
362 _UIntType __factor = 1u;
363 _UIntType __sum = 0u;
364 for (size_t __j = 0; __j < __k; ++__j)
366 __sum += __arr[__k * __i + __j] * __factor;
367 __factor *= __detail::_Shift<_UIntType, 32>::__value;
369 _M_x[__i] = __detail::__mod<_UIntType,
370 __detail::_Shift<_UIntType, __w>::__value>(__sum);
376 if ((_M_x[0] & __upper_mask) != 0u)
379 else if (_M_x[__i] != 0u)
384 _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
388 template<typename _UIntType, size_t __w,
389 size_t __n, size_t __m, size_t __r,
390 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
391 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
394 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
395 __s, __b, __t, __c, __l, __f>::
398 const _UIntType __upper_mask = (~_UIntType()) << __r;
399 const _UIntType __lower_mask = ~__upper_mask;
401 for (size_t __k = 0; __k < (__n - __m); ++__k)
403 _UIntType __y = ((_M_x[__k] & __upper_mask)
404 | (_M_x[__k + 1] & __lower_mask));
405 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
406 ^ ((__y & 0x01) ? __a : 0));
409 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
411 _UIntType __y = ((_M_x[__k] & __upper_mask)
412 | (_M_x[__k + 1] & __lower_mask));
413 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
414 ^ ((__y & 0x01) ? __a : 0));
417 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
418 | (_M_x[0] & __lower_mask));
419 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
420 ^ ((__y & 0x01) ? __a : 0));
424 template<typename _UIntType, size_t __w,
425 size_t __n, size_t __m, size_t __r,
426 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
427 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
430 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
431 __s, __b, __t, __c, __l, __f>::
432 discard(unsigned long long __z)
434 while (__z > state_size - _M_p)
436 __z -= state_size - _M_p;
442 template<typename _UIntType, size_t __w,
443 size_t __n, size_t __m, size_t __r,
444 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
445 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
448 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
449 __s, __b, __t, __c, __l, __f>::result_type
450 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
451 __s, __b, __t, __c, __l, __f>::
454 // Reload the vector - cost is O(n) amortized over n calls.
455 if (_M_p >= state_size)
458 // Calculate o(x(i)).
459 result_type __z = _M_x[_M_p++];
460 __z ^= (__z >> __u) & __d;
461 __z ^= (__z << __s) & __b;
462 __z ^= (__z << __t) & __c;
468 template<typename _UIntType, size_t __w,
469 size_t __n, size_t __m, size_t __r,
470 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
471 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
472 _UIntType __f, typename _CharT, typename _Traits>
473 std::basic_ostream<_CharT, _Traits>&
474 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
475 const mersenne_twister_engine<_UIntType, __w, __n, __m,
476 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
478 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
479 typedef typename __ostream_type::ios_base __ios_base;
481 const typename __ios_base::fmtflags __flags = __os.flags();
482 const _CharT __fill = __os.fill();
483 const _CharT __space = __os.widen(' ');
484 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
487 for (size_t __i = 0; __i < __n; ++__i)
488 __os << __x._M_x[__i] << __space;
496 template<typename _UIntType, size_t __w,
497 size_t __n, size_t __m, size_t __r,
498 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
499 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
500 _UIntType __f, typename _CharT, typename _Traits>
501 std::basic_istream<_CharT, _Traits>&
502 operator>>(std::basic_istream<_CharT, _Traits>& __is,
503 mersenne_twister_engine<_UIntType, __w, __n, __m,
504 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
506 typedef std::basic_istream<_CharT, _Traits> __istream_type;
507 typedef typename __istream_type::ios_base __ios_base;
509 const typename __ios_base::fmtflags __flags = __is.flags();
510 __is.flags(__ios_base::dec | __ios_base::skipws);
512 for (size_t __i = 0; __i < __n; ++__i)
513 __is >> __x._M_x[__i];
521 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
523 subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
525 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
527 subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
529 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
531 subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
533 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
535 subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
537 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
539 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
540 seed(result_type __value)
542 std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
543 __lcg(__value == 0u ? default_seed : __value);
545 const size_t __n = (__w + 31) / 32;
547 for (size_t __i = 0; __i < long_lag; ++__i)
549 _UIntType __sum = 0u;
550 _UIntType __factor = 1u;
551 for (size_t __j = 0; __j < __n; ++__j)
553 __sum += __detail::__mod<uint_least32_t,
554 __detail::_Shift<uint_least32_t, 32>::__value>
555 (__lcg()) * __factor;
556 __factor *= __detail::_Shift<_UIntType, 32>::__value;
558 _M_x[__i] = __detail::__mod<_UIntType,
559 __detail::_Shift<_UIntType, __w>::__value>(__sum);
561 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
565 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
566 template<typename _Sseq>
567 typename std::enable_if<std::is_class<_Sseq>::value>::type
568 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
571 const size_t __k = (__w + 31) / 32;
572 uint_least32_t __arr[__r * __k];
573 __q.generate(__arr + 0, __arr + __r * __k);
575 for (size_t __i = 0; __i < long_lag; ++__i)
577 _UIntType __sum = 0u;
578 _UIntType __factor = 1u;
579 for (size_t __j = 0; __j < __k; ++__j)
581 __sum += __arr[__k * __i + __j] * __factor;
582 __factor *= __detail::_Shift<_UIntType, 32>::__value;
584 _M_x[__i] = __detail::__mod<_UIntType,
585 __detail::_Shift<_UIntType, __w>::__value>(__sum);
587 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
591 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
592 typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
594 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
597 // Derive short lag index from current index.
598 long __ps = _M_p - short_lag;
602 // Calculate new x(i) without overflow or division.
603 // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
606 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
608 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
613 __xi = (__detail::_Shift<_UIntType, __w>::__value
614 - _M_x[_M_p] - _M_carry + _M_x[__ps]);
619 // Adjust current index to loop around in ring buffer.
620 if (++_M_p >= long_lag)
626 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
627 typename _CharT, typename _Traits>
628 std::basic_ostream<_CharT, _Traits>&
629 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
630 const subtract_with_carry_engine<_UIntType,
633 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
634 typedef typename __ostream_type::ios_base __ios_base;
636 const typename __ios_base::fmtflags __flags = __os.flags();
637 const _CharT __fill = __os.fill();
638 const _CharT __space = __os.widen(' ');
639 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
642 for (size_t __i = 0; __i < __r; ++__i)
643 __os << __x._M_x[__i] << __space;
644 __os << __x._M_carry << __space << __x._M_p;
651 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
652 typename _CharT, typename _Traits>
653 std::basic_istream<_CharT, _Traits>&
654 operator>>(std::basic_istream<_CharT, _Traits>& __is,
655 subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
657 typedef std::basic_ostream<_CharT, _Traits> __istream_type;
658 typedef typename __istream_type::ios_base __ios_base;
660 const typename __ios_base::fmtflags __flags = __is.flags();
661 __is.flags(__ios_base::dec | __ios_base::skipws);
663 for (size_t __i = 0; __i < __r; ++__i)
664 __is >> __x._M_x[__i];
665 __is >> __x._M_carry;
673 template<typename _RandomNumberEngine, size_t __p, size_t __r>
675 discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
677 template<typename _RandomNumberEngine, size_t __p, size_t __r>
679 discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
681 template<typename _RandomNumberEngine, size_t __p, size_t __r>
682 typename discard_block_engine<_RandomNumberEngine,
683 __p, __r>::result_type
684 discard_block_engine<_RandomNumberEngine, __p, __r>::
687 if (_M_n >= used_block)
689 _M_b.discard(block_size - _M_n);
696 template<typename _RandomNumberEngine, size_t __p, size_t __r,
697 typename _CharT, typename _Traits>
698 std::basic_ostream<_CharT, _Traits>&
699 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
700 const discard_block_engine<_RandomNumberEngine,
703 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
704 typedef typename __ostream_type::ios_base __ios_base;
706 const typename __ios_base::fmtflags __flags = __os.flags();
707 const _CharT __fill = __os.fill();
708 const _CharT __space = __os.widen(' ');
709 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
712 __os << __x.base() << __space << __x._M_n;
719 template<typename _RandomNumberEngine, size_t __p, size_t __r,
720 typename _CharT, typename _Traits>
721 std::basic_istream<_CharT, _Traits>&
722 operator>>(std::basic_istream<_CharT, _Traits>& __is,
723 discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
725 typedef std::basic_istream<_CharT, _Traits> __istream_type;
726 typedef typename __istream_type::ios_base __ios_base;
728 const typename __ios_base::fmtflags __flags = __is.flags();
729 __is.flags(__ios_base::dec | __ios_base::skipws);
731 __is >> __x._M_b >> __x._M_n;
738 template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
739 typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
741 independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
744 typedef typename _RandomNumberEngine::result_type _Eresult_type;
745 const _Eresult_type __r
746 = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
747 ? _M_b.max() - _M_b.min() + 1 : 0);
748 const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
749 const unsigned __m = __r ? std::__lg(__r) : __edig;
751 typedef typename std::common_type<_Eresult_type, result_type>::type
753 const unsigned __cdig = std::numeric_limits<__ctype>::digits;
756 __ctype __s0, __s1, __y0, __y1;
758 for (size_t __i = 0; __i < 2; ++__i)
760 __n = (__w + __m - 1) / __m + __i;
761 __n0 = __n - __w % __n;
762 const unsigned __w0 = __w / __n; // __w0 <= __m
768 __s0 = __ctype(1) << __w0;
776 __y0 = __s0 * (__r / __s0);
778 __y1 = __s1 * (__r / __s1);
780 if (__r - __y0 <= __y0 / __n)
787 result_type __sum = 0;
788 for (size_t __k = 0; __k < __n0; ++__k)
792 __u = _M_b() - _M_b.min();
793 while (__y0 && __u >= __y0);
794 __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
796 for (size_t __k = __n0; __k < __n; ++__k)
800 __u = _M_b() - _M_b.min();
801 while (__y1 && __u >= __y1);
802 __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
808 template<typename _RandomNumberEngine, size_t __k>
810 shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
812 template<typename _RandomNumberEngine, size_t __k>
813 typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
814 shuffle_order_engine<_RandomNumberEngine, __k>::
817 size_t __j = __k * ((_M_y - _M_b.min())
818 / (_M_b.max() - _M_b.min() + 1.0L));
825 template<typename _RandomNumberEngine, size_t __k,
826 typename _CharT, typename _Traits>
827 std::basic_ostream<_CharT, _Traits>&
828 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
829 const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
831 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
832 typedef typename __ostream_type::ios_base __ios_base;
834 const typename __ios_base::fmtflags __flags = __os.flags();
835 const _CharT __fill = __os.fill();
836 const _CharT __space = __os.widen(' ');
837 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
841 for (size_t __i = 0; __i < __k; ++__i)
842 __os << __space << __x._M_v[__i];
843 __os << __space << __x._M_y;
850 template<typename _RandomNumberEngine, size_t __k,
851 typename _CharT, typename _Traits>
852 std::basic_istream<_CharT, _Traits>&
853 operator>>(std::basic_istream<_CharT, _Traits>& __is,
854 shuffle_order_engine<_RandomNumberEngine, __k>& __x)
856 typedef std::basic_istream<_CharT, _Traits> __istream_type;
857 typedef typename __istream_type::ios_base __ios_base;
859 const typename __ios_base::fmtflags __flags = __is.flags();
860 __is.flags(__ios_base::dec | __ios_base::skipws);
863 for (size_t __i = 0; __i < __k; ++__i)
864 __is >> __x._M_v[__i];
872 template<typename _IntType, typename _CharT, typename _Traits>
873 std::basic_ostream<_CharT, _Traits>&
874 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
875 const uniform_int_distribution<_IntType>& __x)
877 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
878 typedef typename __ostream_type::ios_base __ios_base;
880 const typename __ios_base::fmtflags __flags = __os.flags();
881 const _CharT __fill = __os.fill();
882 const _CharT __space = __os.widen(' ');
883 __os.flags(__ios_base::scientific | __ios_base::left);
886 __os << __x.a() << __space << __x.b();
893 template<typename _IntType, typename _CharT, typename _Traits>
894 std::basic_istream<_CharT, _Traits>&
895 operator>>(std::basic_istream<_CharT, _Traits>& __is,
896 uniform_int_distribution<_IntType>& __x)
898 typedef std::basic_istream<_CharT, _Traits> __istream_type;
899 typedef typename __istream_type::ios_base __ios_base;
901 const typename __ios_base::fmtflags __flags = __is.flags();
902 __is.flags(__ios_base::dec | __ios_base::skipws);
906 __x.param(typename uniform_int_distribution<_IntType>::
907 param_type(__a, __b));
914 template<typename _RealType>
915 template<typename _ForwardIterator,
916 typename _UniformRandomNumberGenerator>
918 uniform_real_distribution<_RealType>::
919 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
920 _UniformRandomNumberGenerator& __urng,
921 const param_type& __p)
923 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
924 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
926 auto __range = __p.b() - __p.a();
928 *__f++ = __aurng() * __range + __p.a();
931 template<typename _RealType, typename _CharT, typename _Traits>
932 std::basic_ostream<_CharT, _Traits>&
933 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
934 const uniform_real_distribution<_RealType>& __x)
936 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
937 typedef typename __ostream_type::ios_base __ios_base;
939 const typename __ios_base::fmtflags __flags = __os.flags();
940 const _CharT __fill = __os.fill();
941 const std::streamsize __precision = __os.precision();
942 const _CharT __space = __os.widen(' ');
943 __os.flags(__ios_base::scientific | __ios_base::left);
945 __os.precision(std::numeric_limits<_RealType>::max_digits10);
947 __os << __x.a() << __space << __x.b();
951 __os.precision(__precision);
955 template<typename _RealType, typename _CharT, typename _Traits>
956 std::basic_istream<_CharT, _Traits>&
957 operator>>(std::basic_istream<_CharT, _Traits>& __is,
958 uniform_real_distribution<_RealType>& __x)
960 typedef std::basic_istream<_CharT, _Traits> __istream_type;
961 typedef typename __istream_type::ios_base __ios_base;
963 const typename __ios_base::fmtflags __flags = __is.flags();
964 __is.flags(__ios_base::skipws);
968 __x.param(typename uniform_real_distribution<_RealType>::
969 param_type(__a, __b));
976 template<typename _ForwardIterator,
977 typename _UniformRandomNumberGenerator>
979 std::bernoulli_distribution::
980 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
981 _UniformRandomNumberGenerator& __urng,
982 const param_type& __p)
984 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
985 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
987 auto __limit = __p.p() * (__aurng.max() - __aurng.min());
990 *__f++ = (__aurng() - __aurng.min()) < __limit;
993 template<typename _CharT, typename _Traits>
994 std::basic_ostream<_CharT, _Traits>&
995 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
996 const bernoulli_distribution& __x)
998 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
999 typedef typename __ostream_type::ios_base __ios_base;
1001 const typename __ios_base::fmtflags __flags = __os.flags();
1002 const _CharT __fill = __os.fill();
1003 const std::streamsize __precision = __os.precision();
1004 __os.flags(__ios_base::scientific | __ios_base::left);
1005 __os.fill(__os.widen(' '));
1006 __os.precision(std::numeric_limits<double>::max_digits10);
1010 __os.flags(__flags);
1012 __os.precision(__precision);
1017 template<typename _IntType>
1018 template<typename _UniformRandomNumberGenerator>
1019 typename geometric_distribution<_IntType>::result_type
1020 geometric_distribution<_IntType>::
1021 operator()(_UniformRandomNumberGenerator& __urng,
1022 const param_type& __param)
1024 // About the epsilon thing see this thread:
1025 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1026 const double __naf =
1027 (1 - std::numeric_limits<double>::epsilon()) / 2;
1028 // The largest _RealType convertible to _IntType.
1029 const double __thr =
1030 std::numeric_limits<_IntType>::max() + __naf;
1031 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1036 __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1037 while (__cand >= __thr);
1039 return result_type(__cand + __naf);
1042 template<typename _IntType>
1043 template<typename _ForwardIterator,
1044 typename _UniformRandomNumberGenerator>
1046 geometric_distribution<_IntType>::
1047 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1048 _UniformRandomNumberGenerator& __urng,
1049 const param_type& __param)
1051 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1052 // About the epsilon thing see this thread:
1053 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1054 const double __naf =
1055 (1 - std::numeric_limits<double>::epsilon()) / 2;
1056 // The largest _RealType convertible to _IntType.
1057 const double __thr =
1058 std::numeric_limits<_IntType>::max() + __naf;
1059 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1066 __cand = std::floor(std::log(1.0 - __aurng())
1067 / __param._M_log_1_p);
1068 while (__cand >= __thr);
1070 *__f++ = __cand + __naf;
1074 template<typename _IntType,
1075 typename _CharT, typename _Traits>
1076 std::basic_ostream<_CharT, _Traits>&
1077 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1078 const geometric_distribution<_IntType>& __x)
1080 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1081 typedef typename __ostream_type::ios_base __ios_base;
1083 const typename __ios_base::fmtflags __flags = __os.flags();
1084 const _CharT __fill = __os.fill();
1085 const std::streamsize __precision = __os.precision();
1086 __os.flags(__ios_base::scientific | __ios_base::left);
1087 __os.fill(__os.widen(' '));
1088 __os.precision(std::numeric_limits<double>::max_digits10);
1092 __os.flags(__flags);
1094 __os.precision(__precision);
1098 template<typename _IntType,
1099 typename _CharT, typename _Traits>
1100 std::basic_istream<_CharT, _Traits>&
1101 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1102 geometric_distribution<_IntType>& __x)
1104 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1105 typedef typename __istream_type::ios_base __ios_base;
1107 const typename __ios_base::fmtflags __flags = __is.flags();
1108 __is.flags(__ios_base::skipws);
1112 __x.param(typename geometric_distribution<_IntType>::param_type(__p));
1114 __is.flags(__flags);
1118 // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1119 template<typename _IntType>
1120 template<typename _UniformRandomNumberGenerator>
1121 typename negative_binomial_distribution<_IntType>::result_type
1122 negative_binomial_distribution<_IntType>::
1123 operator()(_UniformRandomNumberGenerator& __urng)
1125 const double __y = _M_gd(__urng);
1127 // XXX Is the constructor too slow?
1128 std::poisson_distribution<result_type> __poisson(__y);
1129 return __poisson(__urng);
1132 template<typename _IntType>
1133 template<typename _UniformRandomNumberGenerator>
1134 typename negative_binomial_distribution<_IntType>::result_type
1135 negative_binomial_distribution<_IntType>::
1136 operator()(_UniformRandomNumberGenerator& __urng,
1137 const param_type& __p)
1139 typedef typename std::gamma_distribution<double>::param_type
1143 _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1145 std::poisson_distribution<result_type> __poisson(__y);
1146 return __poisson(__urng);
1149 template<typename _IntType>
1150 template<typename _ForwardIterator,
1151 typename _UniformRandomNumberGenerator>
1153 negative_binomial_distribution<_IntType>::
1154 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1155 _UniformRandomNumberGenerator& __urng)
1157 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1160 const double __y = _M_gd(__urng);
1162 // XXX Is the constructor too slow?
1163 std::poisson_distribution<result_type> __poisson(__y);
1164 *__f++ = __poisson(__urng);
1168 template<typename _IntType>
1169 template<typename _ForwardIterator,
1170 typename _UniformRandomNumberGenerator>
1172 negative_binomial_distribution<_IntType>::
1173 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1174 _UniformRandomNumberGenerator& __urng,
1175 const param_type& __p)
1177 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1178 typename std::gamma_distribution<result_type>::param_type
1179 __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1183 const double __y = _M_gd(__urng, __p2);
1185 std::poisson_distribution<result_type> __poisson(__y);
1186 *__f++ = __poisson(__urng);
1190 template<typename _IntType, typename _CharT, typename _Traits>
1191 std::basic_ostream<_CharT, _Traits>&
1192 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1193 const negative_binomial_distribution<_IntType>& __x)
1195 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1196 typedef typename __ostream_type::ios_base __ios_base;
1198 const typename __ios_base::fmtflags __flags = __os.flags();
1199 const _CharT __fill = __os.fill();
1200 const std::streamsize __precision = __os.precision();
1201 const _CharT __space = __os.widen(' ');
1202 __os.flags(__ios_base::scientific | __ios_base::left);
1203 __os.fill(__os.widen(' '));
1204 __os.precision(std::numeric_limits<double>::max_digits10);
1206 __os << __x.k() << __space << __x.p()
1207 << __space << __x._M_gd;
1209 __os.flags(__flags);
1211 __os.precision(__precision);
1215 template<typename _IntType, typename _CharT, typename _Traits>
1216 std::basic_istream<_CharT, _Traits>&
1217 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1218 negative_binomial_distribution<_IntType>& __x)
1220 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1221 typedef typename __istream_type::ios_base __ios_base;
1223 const typename __ios_base::fmtflags __flags = __is.flags();
1224 __is.flags(__ios_base::skipws);
1228 __is >> __k >> __p >> __x._M_gd;
1229 __x.param(typename negative_binomial_distribution<_IntType>::
1230 param_type(__k, __p));
1232 __is.flags(__flags);
1237 template<typename _IntType>
1239 poisson_distribution<_IntType>::param_type::
1242 #if _GLIBCXX_USE_C99_MATH_TR1
1245 const double __m = std::floor(_M_mean);
1246 _M_lm_thr = std::log(_M_mean);
1247 _M_lfm = std::lgamma(__m + 1);
1248 _M_sm = std::sqrt(__m);
1250 const double __pi_4 = 0.7853981633974483096156608458198757L;
1251 const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1253 _M_d = std::round(std::max<double>(6.0, std::min(__m, __dx)));
1254 const double __cx = 2 * __m + _M_d;
1255 _M_scx = std::sqrt(__cx / 2);
1258 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1259 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1264 _M_lm_thr = std::exp(-_M_mean);
1268 * A rejection algorithm when mean >= 12 and a simple method based
1269 * upon the multiplication of uniform random variates otherwise.
1270 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1274 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1275 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1277 template<typename _IntType>
1278 template<typename _UniformRandomNumberGenerator>
1279 typename poisson_distribution<_IntType>::result_type
1280 poisson_distribution<_IntType>::
1281 operator()(_UniformRandomNumberGenerator& __urng,
1282 const param_type& __param)
1284 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1286 #if _GLIBCXX_USE_C99_MATH_TR1
1287 if (__param.mean() >= 12)
1291 // See comments above...
1292 const double __naf =
1293 (1 - std::numeric_limits<double>::epsilon()) / 2;
1294 const double __thr =
1295 std::numeric_limits<_IntType>::max() + __naf;
1297 const double __m = std::floor(__param.mean());
1299 const double __spi_2 = 1.2533141373155002512078826424055226L;
1300 const double __c1 = __param._M_sm * __spi_2;
1301 const double __c2 = __param._M_c2b + __c1;
1302 const double __c3 = __c2 + 1;
1303 const double __c4 = __c3 + 1;
1305 const double __178 = 0.0128205128205128205128205128205128L;
1307 const double __e178 = 1.0129030479320018583185514777512983L;
1308 const double __c5 = __c4 + __e178;
1309 const double __c = __param._M_cb + __c5;
1310 const double __2cx = 2 * (2 * __m + __param._M_d);
1312 bool __reject = true;
1315 const double __u = __c * __aurng();
1316 const double __e = -std::log(1.0 - __aurng());
1322 const double __n = _M_nd(__urng);
1323 const double __y = -std::abs(__n) * __param._M_sm - 1;
1324 __x = std::floor(__y);
1325 __w = -__n * __n / 2;
1329 else if (__u <= __c2)
1331 const double __n = _M_nd(__urng);
1332 const double __y = 1 + std::abs(__n) * __param._M_scx;
1333 __x = std::ceil(__y);
1334 __w = __y * (2 - __y) * __param._M_1cx;
1335 if (__x > __param._M_d)
1338 else if (__u <= __c3)
1339 // NB: This case not in the book, nor in the Errata,
1340 // but should be ok...
1342 else if (__u <= __c4)
1344 else if (__u <= __c5)
1347 // Only in the Errata, see libstdc++/83237.
1352 const double __v = -std::log(1.0 - __aurng());
1353 const double __y = __param._M_d
1354 + __v * __2cx / __param._M_d;
1355 __x = std::ceil(__y);
1356 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1359 __reject = (__w - __e - __x * __param._M_lm_thr
1360 > __param._M_lfm - std::lgamma(__x + __m + 1));
1362 __reject |= __x + __m >= __thr;
1366 return result_type(__x + __m + __naf);
1372 double __prod = 1.0;
1376 __prod *= __aurng();
1379 while (__prod > __param._M_lm_thr);
1385 template<typename _IntType>
1386 template<typename _ForwardIterator,
1387 typename _UniformRandomNumberGenerator>
1389 poisson_distribution<_IntType>::
1390 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1391 _UniformRandomNumberGenerator& __urng,
1392 const param_type& __param)
1394 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1395 // We could duplicate everything from operator()...
1397 *__f++ = this->operator()(__urng, __param);
1400 template<typename _IntType,
1401 typename _CharT, typename _Traits>
1402 std::basic_ostream<_CharT, _Traits>&
1403 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1404 const poisson_distribution<_IntType>& __x)
1406 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1407 typedef typename __ostream_type::ios_base __ios_base;
1409 const typename __ios_base::fmtflags __flags = __os.flags();
1410 const _CharT __fill = __os.fill();
1411 const std::streamsize __precision = __os.precision();
1412 const _CharT __space = __os.widen(' ');
1413 __os.flags(__ios_base::scientific | __ios_base::left);
1415 __os.precision(std::numeric_limits<double>::max_digits10);
1417 __os << __x.mean() << __space << __x._M_nd;
1419 __os.flags(__flags);
1421 __os.precision(__precision);
1425 template<typename _IntType,
1426 typename _CharT, typename _Traits>
1427 std::basic_istream<_CharT, _Traits>&
1428 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1429 poisson_distribution<_IntType>& __x)
1431 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1432 typedef typename __istream_type::ios_base __ios_base;
1434 const typename __ios_base::fmtflags __flags = __is.flags();
1435 __is.flags(__ios_base::skipws);
1438 __is >> __mean >> __x._M_nd;
1439 __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1441 __is.flags(__flags);
1446 template<typename _IntType>
1448 binomial_distribution<_IntType>::param_type::
1451 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1455 #if _GLIBCXX_USE_C99_MATH_TR1
1456 if (_M_t * __p12 >= 8)
1459 const double __np = std::floor(_M_t * __p12);
1460 const double __pa = __np / _M_t;
1461 const double __1p = 1 - __pa;
1463 const double __pi_4 = 0.7853981633974483096156608458198757L;
1464 const double __d1x =
1465 std::sqrt(__np * __1p * std::log(32 * __np
1466 / (81 * __pi_4 * __1p)));
1467 _M_d1 = std::round(std::max<double>(1.0, __d1x));
1468 const double __d2x =
1469 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1470 / (__pi_4 * __pa)));
1471 _M_d2 = std::round(std::max<double>(1.0, __d2x));
1474 const double __spi_2 = 1.2533141373155002512078826424055226L;
1475 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1476 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1477 _M_c = 2 * _M_d1 / __np;
1478 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1479 const double __a12 = _M_a1 + _M_s2 * __spi_2;
1480 const double __s1s = _M_s1 * _M_s1;
1481 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1483 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1484 const double __s2s = _M_s2 * _M_s2;
1485 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1486 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1487 _M_lf = (std::lgamma(__np + 1)
1488 + std::lgamma(_M_t - __np + 1));
1489 _M_lp1p = std::log(__pa / __1p);
1491 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1495 _M_q = -std::log(1 - __p12);
1498 template<typename _IntType>
1499 template<typename _UniformRandomNumberGenerator>
1500 typename binomial_distribution<_IntType>::result_type
1501 binomial_distribution<_IntType>::
1502 _M_waiting(_UniformRandomNumberGenerator& __urng,
1503 _IntType __t, double __q)
1507 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1514 const double __e = -std::log(1.0 - __aurng());
1515 __sum += __e / (__t - __x);
1518 while (__sum <= __q);
1524 * A rejection algorithm when t * p >= 8 and a simple waiting time
1525 * method - the second in the referenced book - otherwise.
1526 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1530 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1531 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1533 template<typename _IntType>
1534 template<typename _UniformRandomNumberGenerator>
1535 typename binomial_distribution<_IntType>::result_type
1536 binomial_distribution<_IntType>::
1537 operator()(_UniformRandomNumberGenerator& __urng,
1538 const param_type& __param)
1541 const _IntType __t = __param.t();
1542 const double __p = __param.p();
1543 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1544 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1547 #if _GLIBCXX_USE_C99_MATH_TR1
1548 if (!__param._M_easy)
1552 // See comments above...
1553 const double __naf =
1554 (1 - std::numeric_limits<double>::epsilon()) / 2;
1555 const double __thr =
1556 std::numeric_limits<_IntType>::max() + __naf;
1558 const double __np = std::floor(__t * __p12);
1561 const double __spi_2 = 1.2533141373155002512078826424055226L;
1562 const double __a1 = __param._M_a1;
1563 const double __a12 = __a1 + __param._M_s2 * __spi_2;
1564 const double __a123 = __param._M_a123;
1565 const double __s1s = __param._M_s1 * __param._M_s1;
1566 const double __s2s = __param._M_s2 * __param._M_s2;
1571 const double __u = __param._M_s * __aurng();
1577 const double __n = _M_nd(__urng);
1578 const double __y = __param._M_s1 * std::abs(__n);
1579 __reject = __y >= __param._M_d1;
1582 const double __e = -std::log(1.0 - __aurng());
1583 __x = std::floor(__y);
1584 __v = -__e - __n * __n / 2 + __param._M_c;
1587 else if (__u <= __a12)
1589 const double __n = _M_nd(__urng);
1590 const double __y = __param._M_s2 * std::abs(__n);
1591 __reject = __y >= __param._M_d2;
1594 const double __e = -std::log(1.0 - __aurng());
1595 __x = std::floor(-__y);
1596 __v = -__e - __n * __n / 2;
1599 else if (__u <= __a123)
1601 const double __e1 = -std::log(1.0 - __aurng());
1602 const double __e2 = -std::log(1.0 - __aurng());
1604 const double __y = __param._M_d1
1605 + 2 * __s1s * __e1 / __param._M_d1;
1606 __x = std::floor(__y);
1607 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1608 -__y / (2 * __s1s)));
1613 const double __e1 = -std::log(1.0 - __aurng());
1614 const double __e2 = -std::log(1.0 - __aurng());
1616 const double __y = __param._M_d2
1617 + 2 * __s2s * __e1 / __param._M_d2;
1618 __x = std::floor(-__y);
1619 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1623 __reject = __reject || __x < -__np || __x > __t - __np;
1626 const double __lfx =
1627 std::lgamma(__np + __x + 1)
1628 + std::lgamma(__t - (__np + __x) + 1);
1629 __reject = __v > __param._M_lf - __lfx
1630 + __x * __param._M_lp1p;
1633 __reject |= __x + __np >= __thr;
1637 __x += __np + __naf;
1639 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1641 __ret = _IntType(__x) + __z;
1645 __ret = _M_waiting(__urng, __t, __param._M_q);
1648 __ret = __t - __ret;
1652 template<typename _IntType>
1653 template<typename _ForwardIterator,
1654 typename _UniformRandomNumberGenerator>
1656 binomial_distribution<_IntType>::
1657 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1658 _UniformRandomNumberGenerator& __urng,
1659 const param_type& __param)
1661 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1662 // We could duplicate everything from operator()...
1664 *__f++ = this->operator()(__urng, __param);
1667 template<typename _IntType,
1668 typename _CharT, typename _Traits>
1669 std::basic_ostream<_CharT, _Traits>&
1670 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1671 const binomial_distribution<_IntType>& __x)
1673 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1674 typedef typename __ostream_type::ios_base __ios_base;
1676 const typename __ios_base::fmtflags __flags = __os.flags();
1677 const _CharT __fill = __os.fill();
1678 const std::streamsize __precision = __os.precision();
1679 const _CharT __space = __os.widen(' ');
1680 __os.flags(__ios_base::scientific | __ios_base::left);
1682 __os.precision(std::numeric_limits<double>::max_digits10);
1684 __os << __x.t() << __space << __x.p()
1685 << __space << __x._M_nd;
1687 __os.flags(__flags);
1689 __os.precision(__precision);
1693 template<typename _IntType,
1694 typename _CharT, typename _Traits>
1695 std::basic_istream<_CharT, _Traits>&
1696 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1697 binomial_distribution<_IntType>& __x)
1699 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1700 typedef typename __istream_type::ios_base __ios_base;
1702 const typename __ios_base::fmtflags __flags = __is.flags();
1703 __is.flags(__ios_base::dec | __ios_base::skipws);
1707 __is >> __t >> __p >> __x._M_nd;
1708 __x.param(typename binomial_distribution<_IntType>::
1709 param_type(__t, __p));
1711 __is.flags(__flags);
1716 template<typename _RealType>
1717 template<typename _ForwardIterator,
1718 typename _UniformRandomNumberGenerator>
1720 std::exponential_distribution<_RealType>::
1721 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1722 _UniformRandomNumberGenerator& __urng,
1723 const param_type& __p)
1725 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1726 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1729 *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1732 template<typename _RealType, typename _CharT, typename _Traits>
1733 std::basic_ostream<_CharT, _Traits>&
1734 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1735 const exponential_distribution<_RealType>& __x)
1737 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1738 typedef typename __ostream_type::ios_base __ios_base;
1740 const typename __ios_base::fmtflags __flags = __os.flags();
1741 const _CharT __fill = __os.fill();
1742 const std::streamsize __precision = __os.precision();
1743 __os.flags(__ios_base::scientific | __ios_base::left);
1744 __os.fill(__os.widen(' '));
1745 __os.precision(std::numeric_limits<_RealType>::max_digits10);
1747 __os << __x.lambda();
1749 __os.flags(__flags);
1751 __os.precision(__precision);
1755 template<typename _RealType, typename _CharT, typename _Traits>
1756 std::basic_istream<_CharT, _Traits>&
1757 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1758 exponential_distribution<_RealType>& __x)
1760 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1761 typedef typename __istream_type::ios_base __ios_base;
1763 const typename __ios_base::fmtflags __flags = __is.flags();
1764 __is.flags(__ios_base::dec | __ios_base::skipws);
1768 __x.param(typename exponential_distribution<_RealType>::
1769 param_type(__lambda));
1771 __is.flags(__flags);
1777 * Polar method due to Marsaglia.
1779 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1780 * New York, 1986, Ch. V, Sect. 4.4.
1782 template<typename _RealType>
1783 template<typename _UniformRandomNumberGenerator>
1784 typename normal_distribution<_RealType>::result_type
1785 normal_distribution<_RealType>::
1786 operator()(_UniformRandomNumberGenerator& __urng,
1787 const param_type& __param)
1790 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1793 if (_M_saved_available)
1795 _M_saved_available = false;
1800 result_type __x, __y, __r2;
1803 __x = result_type(2.0) * __aurng() - 1.0;
1804 __y = result_type(2.0) * __aurng() - 1.0;
1805 __r2 = __x * __x + __y * __y;
1807 while (__r2 > 1.0 || __r2 == 0.0);
1809 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1810 _M_saved = __x * __mult;
1811 _M_saved_available = true;
1812 __ret = __y * __mult;
1815 __ret = __ret * __param.stddev() + __param.mean();
1819 template<typename _RealType>
1820 template<typename _ForwardIterator,
1821 typename _UniformRandomNumberGenerator>
1823 normal_distribution<_RealType>::
1824 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1825 _UniformRandomNumberGenerator& __urng,
1826 const param_type& __param)
1828 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1833 if (_M_saved_available)
1835 _M_saved_available = false;
1836 *__f++ = _M_saved * __param.stddev() + __param.mean();
1842 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1845 while (__f + 1 < __t)
1847 result_type __x, __y, __r2;
1850 __x = result_type(2.0) * __aurng() - 1.0;
1851 __y = result_type(2.0) * __aurng() - 1.0;
1852 __r2 = __x * __x + __y * __y;
1854 while (__r2 > 1.0 || __r2 == 0.0);
1856 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1857 *__f++ = __y * __mult * __param.stddev() + __param.mean();
1858 *__f++ = __x * __mult * __param.stddev() + __param.mean();
1863 result_type __x, __y, __r2;
1866 __x = result_type(2.0) * __aurng() - 1.0;
1867 __y = result_type(2.0) * __aurng() - 1.0;
1868 __r2 = __x * __x + __y * __y;
1870 while (__r2 > 1.0 || __r2 == 0.0);
1872 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1873 _M_saved = __x * __mult;
1874 _M_saved_available = true;
1875 *__f = __y * __mult * __param.stddev() + __param.mean();
1879 template<typename _RealType>
1881 operator==(const std::normal_distribution<_RealType>& __d1,
1882 const std::normal_distribution<_RealType>& __d2)
1884 if (__d1._M_param == __d2._M_param
1885 && __d1._M_saved_available == __d2._M_saved_available)
1887 if (__d1._M_saved_available
1888 && __d1._M_saved == __d2._M_saved)
1890 else if(!__d1._M_saved_available)
1899 template<typename _RealType, typename _CharT, typename _Traits>
1900 std::basic_ostream<_CharT, _Traits>&
1901 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1902 const normal_distribution<_RealType>& __x)
1904 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1905 typedef typename __ostream_type::ios_base __ios_base;
1907 const typename __ios_base::fmtflags __flags = __os.flags();
1908 const _CharT __fill = __os.fill();
1909 const std::streamsize __precision = __os.precision();
1910 const _CharT __space = __os.widen(' ');
1911 __os.flags(__ios_base::scientific | __ios_base::left);
1913 __os.precision(std::numeric_limits<_RealType>::max_digits10);
1915 __os << __x.mean() << __space << __x.stddev()
1916 << __space << __x._M_saved_available;
1917 if (__x._M_saved_available)
1918 __os << __space << __x._M_saved;
1920 __os.flags(__flags);
1922 __os.precision(__precision);
1926 template<typename _RealType, typename _CharT, typename _Traits>
1927 std::basic_istream<_CharT, _Traits>&
1928 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1929 normal_distribution<_RealType>& __x)
1931 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1932 typedef typename __istream_type::ios_base __ios_base;
1934 const typename __ios_base::fmtflags __flags = __is.flags();
1935 __is.flags(__ios_base::dec | __ios_base::skipws);
1937 double __mean, __stddev;
1938 __is >> __mean >> __stddev
1939 >> __x._M_saved_available;
1940 if (__x._M_saved_available)
1941 __is >> __x._M_saved;
1942 __x.param(typename normal_distribution<_RealType>::
1943 param_type(__mean, __stddev));
1945 __is.flags(__flags);
1950 template<typename _RealType>
1951 template<typename _ForwardIterator,
1952 typename _UniformRandomNumberGenerator>
1954 lognormal_distribution<_RealType>::
1955 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1956 _UniformRandomNumberGenerator& __urng,
1957 const param_type& __p)
1959 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1961 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
1964 template<typename _RealType, typename _CharT, typename _Traits>
1965 std::basic_ostream<_CharT, _Traits>&
1966 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1967 const lognormal_distribution<_RealType>& __x)
1969 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1970 typedef typename __ostream_type::ios_base __ios_base;
1972 const typename __ios_base::fmtflags __flags = __os.flags();
1973 const _CharT __fill = __os.fill();
1974 const std::streamsize __precision = __os.precision();
1975 const _CharT __space = __os.widen(' ');
1976 __os.flags(__ios_base::scientific | __ios_base::left);
1978 __os.precision(std::numeric_limits<_RealType>::max_digits10);
1980 __os << __x.m() << __space << __x.s()
1981 << __space << __x._M_nd;
1983 __os.flags(__flags);
1985 __os.precision(__precision);
1989 template<typename _RealType, typename _CharT, typename _Traits>
1990 std::basic_istream<_CharT, _Traits>&
1991 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1992 lognormal_distribution<_RealType>& __x)
1994 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1995 typedef typename __istream_type::ios_base __ios_base;
1997 const typename __ios_base::fmtflags __flags = __is.flags();
1998 __is.flags(__ios_base::dec | __ios_base::skipws);
2001 __is >> __m >> __s >> __x._M_nd;
2002 __x.param(typename lognormal_distribution<_RealType>::
2003 param_type(__m, __s));
2005 __is.flags(__flags);
2009 template<typename _RealType>
2010 template<typename _ForwardIterator,
2011 typename _UniformRandomNumberGenerator>
2013 std::chi_squared_distribution<_RealType>::
2014 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2015 _UniformRandomNumberGenerator& __urng)
2017 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2019 *__f++ = 2 * _M_gd(__urng);
2022 template<typename _RealType>
2023 template<typename _ForwardIterator,
2024 typename _UniformRandomNumberGenerator>
2026 std::chi_squared_distribution<_RealType>::
2027 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2028 _UniformRandomNumberGenerator& __urng,
2030 std::gamma_distribution<result_type>::param_type& __p)
2032 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2034 *__f++ = 2 * _M_gd(__urng, __p);
2037 template<typename _RealType, typename _CharT, typename _Traits>
2038 std::basic_ostream<_CharT, _Traits>&
2039 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2040 const chi_squared_distribution<_RealType>& __x)
2042 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2043 typedef typename __ostream_type::ios_base __ios_base;
2045 const typename __ios_base::fmtflags __flags = __os.flags();
2046 const _CharT __fill = __os.fill();
2047 const std::streamsize __precision = __os.precision();
2048 const _CharT __space = __os.widen(' ');
2049 __os.flags(__ios_base::scientific | __ios_base::left);
2051 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2053 __os << __x.n() << __space << __x._M_gd;
2055 __os.flags(__flags);
2057 __os.precision(__precision);
2061 template<typename _RealType, typename _CharT, typename _Traits>
2062 std::basic_istream<_CharT, _Traits>&
2063 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2064 chi_squared_distribution<_RealType>& __x)
2066 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2067 typedef typename __istream_type::ios_base __ios_base;
2069 const typename __ios_base::fmtflags __flags = __is.flags();
2070 __is.flags(__ios_base::dec | __ios_base::skipws);
2073 __is >> __n >> __x._M_gd;
2074 __x.param(typename chi_squared_distribution<_RealType>::
2077 __is.flags(__flags);
2082 template<typename _RealType>
2083 template<typename _UniformRandomNumberGenerator>
2084 typename cauchy_distribution<_RealType>::result_type
2085 cauchy_distribution<_RealType>::
2086 operator()(_UniformRandomNumberGenerator& __urng,
2087 const param_type& __p)
2089 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2096 const _RealType __pi = 3.1415926535897932384626433832795029L;
2097 return __p.a() + __p.b() * std::tan(__pi * __u);
2100 template<typename _RealType>
2101 template<typename _ForwardIterator,
2102 typename _UniformRandomNumberGenerator>
2104 cauchy_distribution<_RealType>::
2105 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2106 _UniformRandomNumberGenerator& __urng,
2107 const param_type& __p)
2109 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2110 const _RealType __pi = 3.1415926535897932384626433832795029L;
2111 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2120 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2124 template<typename _RealType, typename _CharT, typename _Traits>
2125 std::basic_ostream<_CharT, _Traits>&
2126 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2127 const cauchy_distribution<_RealType>& __x)
2129 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2130 typedef typename __ostream_type::ios_base __ios_base;
2132 const typename __ios_base::fmtflags __flags = __os.flags();
2133 const _CharT __fill = __os.fill();
2134 const std::streamsize __precision = __os.precision();
2135 const _CharT __space = __os.widen(' ');
2136 __os.flags(__ios_base::scientific | __ios_base::left);
2138 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2140 __os << __x.a() << __space << __x.b();
2142 __os.flags(__flags);
2144 __os.precision(__precision);
2148 template<typename _RealType, typename _CharT, typename _Traits>
2149 std::basic_istream<_CharT, _Traits>&
2150 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2151 cauchy_distribution<_RealType>& __x)
2153 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2154 typedef typename __istream_type::ios_base __ios_base;
2156 const typename __ios_base::fmtflags __flags = __is.flags();
2157 __is.flags(__ios_base::dec | __ios_base::skipws);
2161 __x.param(typename cauchy_distribution<_RealType>::
2162 param_type(__a, __b));
2164 __is.flags(__flags);
2169 template<typename _RealType>
2170 template<typename _ForwardIterator,
2171 typename _UniformRandomNumberGenerator>
2173 std::fisher_f_distribution<_RealType>::
2174 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2175 _UniformRandomNumberGenerator& __urng)
2177 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2179 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2182 template<typename _RealType>
2183 template<typename _ForwardIterator,
2184 typename _UniformRandomNumberGenerator>
2186 std::fisher_f_distribution<_RealType>::
2187 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2188 _UniformRandomNumberGenerator& __urng,
2189 const param_type& __p)
2191 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2192 typedef typename std::gamma_distribution<result_type>::param_type
2194 param_type __p1(__p.m() / 2);
2195 param_type __p2(__p.n() / 2);
2197 *__f++ = ((_M_gd_x(__urng, __p1) * n())
2198 / (_M_gd_y(__urng, __p2) * m()));
2201 template<typename _RealType, typename _CharT, typename _Traits>
2202 std::basic_ostream<_CharT, _Traits>&
2203 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2204 const fisher_f_distribution<_RealType>& __x)
2206 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2207 typedef typename __ostream_type::ios_base __ios_base;
2209 const typename __ios_base::fmtflags __flags = __os.flags();
2210 const _CharT __fill = __os.fill();
2211 const std::streamsize __precision = __os.precision();
2212 const _CharT __space = __os.widen(' ');
2213 __os.flags(__ios_base::scientific | __ios_base::left);
2215 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2217 __os << __x.m() << __space << __x.n()
2218 << __space << __x._M_gd_x << __space << __x._M_gd_y;
2220 __os.flags(__flags);
2222 __os.precision(__precision);
2226 template<typename _RealType, typename _CharT, typename _Traits>
2227 std::basic_istream<_CharT, _Traits>&
2228 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2229 fisher_f_distribution<_RealType>& __x)
2231 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2232 typedef typename __istream_type::ios_base __ios_base;
2234 const typename __ios_base::fmtflags __flags = __is.flags();
2235 __is.flags(__ios_base::dec | __ios_base::skipws);
2238 __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
2239 __x.param(typename fisher_f_distribution<_RealType>::
2240 param_type(__m, __n));
2242 __is.flags(__flags);
2247 template<typename _RealType>
2248 template<typename _ForwardIterator,
2249 typename _UniformRandomNumberGenerator>
2251 std::student_t_distribution<_RealType>::
2252 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2253 _UniformRandomNumberGenerator& __urng)
2255 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2257 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2260 template<typename _RealType>
2261 template<typename _ForwardIterator,
2262 typename _UniformRandomNumberGenerator>
2264 std::student_t_distribution<_RealType>::
2265 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2266 _UniformRandomNumberGenerator& __urng,
2267 const param_type& __p)
2269 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2270 typename std::gamma_distribution<result_type>::param_type
2271 __p2(__p.n() / 2, 2);
2273 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2276 template<typename _RealType, typename _CharT, typename _Traits>
2277 std::basic_ostream<_CharT, _Traits>&
2278 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2279 const student_t_distribution<_RealType>& __x)
2281 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2282 typedef typename __ostream_type::ios_base __ios_base;
2284 const typename __ios_base::fmtflags __flags = __os.flags();
2285 const _CharT __fill = __os.fill();
2286 const std::streamsize __precision = __os.precision();
2287 const _CharT __space = __os.widen(' ');
2288 __os.flags(__ios_base::scientific | __ios_base::left);
2290 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2292 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2294 __os.flags(__flags);
2296 __os.precision(__precision);
2300 template<typename _RealType, typename _CharT, typename _Traits>
2301 std::basic_istream<_CharT, _Traits>&
2302 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2303 student_t_distribution<_RealType>& __x)
2305 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2306 typedef typename __istream_type::ios_base __ios_base;
2308 const typename __ios_base::fmtflags __flags = __is.flags();
2309 __is.flags(__ios_base::dec | __ios_base::skipws);
2312 __is >> __n >> __x._M_nd >> __x._M_gd;
2313 __x.param(typename student_t_distribution<_RealType>::param_type(__n));
2315 __is.flags(__flags);
2320 template<typename _RealType>
2322 gamma_distribution<_RealType>::param_type::
2325 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2327 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2328 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2332 * Marsaglia, G. and Tsang, W. W.
2333 * "A Simple Method for Generating Gamma Variables"
2334 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2336 template<typename _RealType>
2337 template<typename _UniformRandomNumberGenerator>
2338 typename gamma_distribution<_RealType>::result_type
2339 gamma_distribution<_RealType>::
2340 operator()(_UniformRandomNumberGenerator& __urng,
2341 const param_type& __param)
2343 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2346 result_type __u, __v, __n;
2347 const result_type __a1 = (__param._M_malpha
2348 - _RealType(1.0) / _RealType(3.0));
2354 __n = _M_nd(__urng);
2355 __v = result_type(1.0) + __param._M_a2 * __n;
2359 __v = __v * __v * __v;
2362 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2363 && (std::log(__u) > (0.5 * __n * __n + __a1
2364 * (1.0 - __v + std::log(__v)))));
2366 if (__param.alpha() == __param._M_malpha)
2367 return __a1 * __v * __param.beta();
2374 return (std::pow(__u, result_type(1.0) / __param.alpha())
2375 * __a1 * __v * __param.beta());
2379 template<typename _RealType>
2380 template<typename _ForwardIterator,
2381 typename _UniformRandomNumberGenerator>
2383 gamma_distribution<_RealType>::
2384 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2385 _UniformRandomNumberGenerator& __urng,
2386 const param_type& __param)
2388 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2389 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2392 result_type __u, __v, __n;
2393 const result_type __a1 = (__param._M_malpha
2394 - _RealType(1.0) / _RealType(3.0));
2396 if (__param.alpha() == __param._M_malpha)
2403 __n = _M_nd(__urng);
2404 __v = result_type(1.0) + __param._M_a2 * __n;
2408 __v = __v * __v * __v;
2411 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2412 && (std::log(__u) > (0.5 * __n * __n + __a1
2413 * (1.0 - __v + std::log(__v)))));
2415 *__f++ = __a1 * __v * __param.beta();
2424 __n = _M_nd(__urng);
2425 __v = result_type(1.0) + __param._M_a2 * __n;
2429 __v = __v * __v * __v;
2432 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2433 && (std::log(__u) > (0.5 * __n * __n + __a1
2434 * (1.0 - __v + std::log(__v)))));
2440 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2441 * __a1 * __v * __param.beta());
2445 template<typename _RealType, typename _CharT, typename _Traits>
2446 std::basic_ostream<_CharT, _Traits>&
2447 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2448 const gamma_distribution<_RealType>& __x)
2450 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2451 typedef typename __ostream_type::ios_base __ios_base;
2453 const typename __ios_base::fmtflags __flags = __os.flags();
2454 const _CharT __fill = __os.fill();
2455 const std::streamsize __precision = __os.precision();
2456 const _CharT __space = __os.widen(' ');
2457 __os.flags(__ios_base::scientific | __ios_base::left);
2459 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2461 __os << __x.alpha() << __space << __x.beta()
2462 << __space << __x._M_nd;
2464 __os.flags(__flags);
2466 __os.precision(__precision);
2470 template<typename _RealType, typename _CharT, typename _Traits>
2471 std::basic_istream<_CharT, _Traits>&
2472 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2473 gamma_distribution<_RealType>& __x)
2475 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2476 typedef typename __istream_type::ios_base __ios_base;
2478 const typename __ios_base::fmtflags __flags = __is.flags();
2479 __is.flags(__ios_base::dec | __ios_base::skipws);
2481 _RealType __alpha_val, __beta_val;
2482 __is >> __alpha_val >> __beta_val >> __x._M_nd;
2483 __x.param(typename gamma_distribution<_RealType>::
2484 param_type(__alpha_val, __beta_val));
2486 __is.flags(__flags);
2491 template<typename _RealType>
2492 template<typename _UniformRandomNumberGenerator>
2493 typename weibull_distribution<_RealType>::result_type
2494 weibull_distribution<_RealType>::
2495 operator()(_UniformRandomNumberGenerator& __urng,
2496 const param_type& __p)
2498 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2500 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2501 result_type(1) / __p.a());
2504 template<typename _RealType>
2505 template<typename _ForwardIterator,
2506 typename _UniformRandomNumberGenerator>
2508 weibull_distribution<_RealType>::
2509 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2510 _UniformRandomNumberGenerator& __urng,
2511 const param_type& __p)
2513 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2514 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2516 auto __inv_a = result_type(1) / __p.a();
2519 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2523 template<typename _RealType, typename _CharT, typename _Traits>
2524 std::basic_ostream<_CharT, _Traits>&
2525 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2526 const weibull_distribution<_RealType>& __x)
2528 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2529 typedef typename __ostream_type::ios_base __ios_base;
2531 const typename __ios_base::fmtflags __flags = __os.flags();
2532 const _CharT __fill = __os.fill();
2533 const std::streamsize __precision = __os.precision();
2534 const _CharT __space = __os.widen(' ');
2535 __os.flags(__ios_base::scientific | __ios_base::left);
2537 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2539 __os << __x.a() << __space << __x.b();
2541 __os.flags(__flags);
2543 __os.precision(__precision);
2547 template<typename _RealType, typename _CharT, typename _Traits>
2548 std::basic_istream<_CharT, _Traits>&
2549 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2550 weibull_distribution<_RealType>& __x)
2552 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2553 typedef typename __istream_type::ios_base __ios_base;
2555 const typename __ios_base::fmtflags __flags = __is.flags();
2556 __is.flags(__ios_base::dec | __ios_base::skipws);
2560 __x.param(typename weibull_distribution<_RealType>::
2561 param_type(__a, __b));
2563 __is.flags(__flags);
2568 template<typename _RealType>
2569 template<typename _UniformRandomNumberGenerator>
2570 typename extreme_value_distribution<_RealType>::result_type
2571 extreme_value_distribution<_RealType>::
2572 operator()(_UniformRandomNumberGenerator& __urng,
2573 const param_type& __p)
2575 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2577 return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2581 template<typename _RealType>
2582 template<typename _ForwardIterator,
2583 typename _UniformRandomNumberGenerator>
2585 extreme_value_distribution<_RealType>::
2586 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2587 _UniformRandomNumberGenerator& __urng,
2588 const param_type& __p)
2590 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2591 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2595 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2599 template<typename _RealType, typename _CharT, typename _Traits>
2600 std::basic_ostream<_CharT, _Traits>&
2601 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2602 const extreme_value_distribution<_RealType>& __x)
2604 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2605 typedef typename __ostream_type::ios_base __ios_base;
2607 const typename __ios_base::fmtflags __flags = __os.flags();
2608 const _CharT __fill = __os.fill();
2609 const std::streamsize __precision = __os.precision();
2610 const _CharT __space = __os.widen(' ');
2611 __os.flags(__ios_base::scientific | __ios_base::left);
2613 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2615 __os << __x.a() << __space << __x.b();
2617 __os.flags(__flags);
2619 __os.precision(__precision);
2623 template<typename _RealType, typename _CharT, typename _Traits>
2624 std::basic_istream<_CharT, _Traits>&
2625 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2626 extreme_value_distribution<_RealType>& __x)
2628 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2629 typedef typename __istream_type::ios_base __ios_base;
2631 const typename __ios_base::fmtflags __flags = __is.flags();
2632 __is.flags(__ios_base::dec | __ios_base::skipws);
2636 __x.param(typename extreme_value_distribution<_RealType>::
2637 param_type(__a, __b));
2639 __is.flags(__flags);
2644 template<typename _IntType>
2646 discrete_distribution<_IntType>::param_type::
2649 if (_M_prob.size() < 2)
2655 const double __sum = std::accumulate(_M_prob.begin(),
2656 _M_prob.end(), 0.0);
2657 // Now normalize the probabilites.
2658 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2660 // Accumulate partial sums.
2661 _M_cp.reserve(_M_prob.size());
2662 std::partial_sum(_M_prob.begin(), _M_prob.end(),
2663 std::back_inserter(_M_cp));
2664 // Make sure the last cumulative probability is one.
2665 _M_cp[_M_cp.size() - 1] = 1.0;
2668 template<typename _IntType>
2669 template<typename _Func>
2670 discrete_distribution<_IntType>::param_type::
2671 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2672 : _M_prob(), _M_cp()
2674 const size_t __n = __nw == 0 ? 1 : __nw;
2675 const double __delta = (__xmax - __xmin) / __n;
2677 _M_prob.reserve(__n);
2678 for (size_t __k = 0; __k < __nw; ++__k)
2679 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2684 template<typename _IntType>
2685 template<typename _UniformRandomNumberGenerator>
2686 typename discrete_distribution<_IntType>::result_type
2687 discrete_distribution<_IntType>::
2688 operator()(_UniformRandomNumberGenerator& __urng,
2689 const param_type& __param)
2691 if (__param._M_cp.empty())
2692 return result_type(0);
2694 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2697 const double __p = __aurng();
2698 auto __pos = std::lower_bound(__param._M_cp.begin(),
2699 __param._M_cp.end(), __p);
2701 return __pos - __param._M_cp.begin();
2704 template<typename _IntType>
2705 template<typename _ForwardIterator,
2706 typename _UniformRandomNumberGenerator>
2708 discrete_distribution<_IntType>::
2709 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2710 _UniformRandomNumberGenerator& __urng,
2711 const param_type& __param)
2713 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2715 if (__param._M_cp.empty())
2718 *__f++ = result_type(0);
2722 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2727 const double __p = __aurng();
2728 auto __pos = std::lower_bound(__param._M_cp.begin(),
2729 __param._M_cp.end(), __p);
2731 *__f++ = __pos - __param._M_cp.begin();
2735 template<typename _IntType, typename _CharT, typename _Traits>
2736 std::basic_ostream<_CharT, _Traits>&
2737 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2738 const discrete_distribution<_IntType>& __x)
2740 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2741 typedef typename __ostream_type::ios_base __ios_base;
2743 const typename __ios_base::fmtflags __flags = __os.flags();
2744 const _CharT __fill = __os.fill();
2745 const std::streamsize __precision = __os.precision();
2746 const _CharT __space = __os.widen(' ');
2747 __os.flags(__ios_base::scientific | __ios_base::left);
2749 __os.precision(std::numeric_limits<double>::max_digits10);
2751 std::vector<double> __prob = __x.probabilities();
2752 __os << __prob.size();
2753 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2754 __os << __space << *__dit;
2756 __os.flags(__flags);
2758 __os.precision(__precision);
2762 template<typename _IntType, typename _CharT, typename _Traits>
2763 std::basic_istream<_CharT, _Traits>&
2764 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2765 discrete_distribution<_IntType>& __x)
2767 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2768 typedef typename __istream_type::ios_base __ios_base;
2770 const typename __ios_base::fmtflags __flags = __is.flags();
2771 __is.flags(__ios_base::dec | __ios_base::skipws);
2776 std::vector<double> __prob_vec;
2777 __prob_vec.reserve(__n);
2778 for (; __n != 0; --__n)
2782 __prob_vec.push_back(__prob);
2785 __x.param(typename discrete_distribution<_IntType>::
2786 param_type(__prob_vec.begin(), __prob_vec.end()));
2788 __is.flags(__flags);
2793 template<typename _RealType>
2795 piecewise_constant_distribution<_RealType>::param_type::
2798 if (_M_int.size() < 2
2799 || (_M_int.size() == 2
2800 && _M_int[0] == _RealType(0)
2801 && _M_int[1] == _RealType(1)))
2808 const double __sum = std::accumulate(_M_den.begin(),
2811 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2814 _M_cp.reserve(_M_den.size());
2815 std::partial_sum(_M_den.begin(), _M_den.end(),
2816 std::back_inserter(_M_cp));
2818 // Make sure the last cumulative probability is one.
2819 _M_cp[_M_cp.size() - 1] = 1.0;
2821 for (size_t __k = 0; __k < _M_den.size(); ++__k)
2822 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2825 template<typename _RealType>
2826 template<typename _InputIteratorB, typename _InputIteratorW>
2827 piecewise_constant_distribution<_RealType>::param_type::
2828 param_type(_InputIteratorB __bbegin,
2829 _InputIteratorB __bend,
2830 _InputIteratorW __wbegin)
2831 : _M_int(), _M_den(), _M_cp()
2833 if (__bbegin != __bend)
2837 _M_int.push_back(*__bbegin);
2839 if (__bbegin == __bend)
2842 _M_den.push_back(*__wbegin);
2850 template<typename _RealType>
2851 template<typename _Func>
2852 piecewise_constant_distribution<_RealType>::param_type::
2853 param_type(initializer_list<_RealType> __bl, _Func __fw)
2854 : _M_int(), _M_den(), _M_cp()
2856 _M_int.reserve(__bl.size());
2857 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2858 _M_int.push_back(*__biter);
2860 _M_den.reserve(_M_int.size() - 1);
2861 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2862 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2867 template<typename _RealType>
2868 template<typename _Func>
2869 piecewise_constant_distribution<_RealType>::param_type::
2870 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2871 : _M_int(), _M_den(), _M_cp()
2873 const size_t __n = __nw == 0 ? 1 : __nw;
2874 const _RealType __delta = (__xmax - __xmin) / __n;
2876 _M_int.reserve(__n + 1);
2877 for (size_t __k = 0; __k <= __nw; ++__k)
2878 _M_int.push_back(__xmin + __k * __delta);
2880 _M_den.reserve(__n);
2881 for (size_t __k = 0; __k < __nw; ++__k)
2882 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2887 template<typename _RealType>
2888 template<typename _UniformRandomNumberGenerator>
2889 typename piecewise_constant_distribution<_RealType>::result_type
2890 piecewise_constant_distribution<_RealType>::
2891 operator()(_UniformRandomNumberGenerator& __urng,
2892 const param_type& __param)
2894 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2897 const double __p = __aurng();
2898 if (__param._M_cp.empty())
2901 auto __pos = std::lower_bound(__param._M_cp.begin(),
2902 __param._M_cp.end(), __p);
2903 const size_t __i = __pos - __param._M_cp.begin();
2905 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2907 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2910 template<typename _RealType>
2911 template<typename _ForwardIterator,
2912 typename _UniformRandomNumberGenerator>
2914 piecewise_constant_distribution<_RealType>::
2915 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2916 _UniformRandomNumberGenerator& __urng,
2917 const param_type& __param)
2919 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2920 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2923 if (__param._M_cp.empty())
2932 const double __p = __aurng();
2934 auto __pos = std::lower_bound(__param._M_cp.begin(),
2935 __param._M_cp.end(), __p);
2936 const size_t __i = __pos - __param._M_cp.begin();
2938 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2940 *__f++ = (__param._M_int[__i]
2941 + (__p - __pref) / __param._M_den[__i]);
2945 template<typename _RealType, typename _CharT, typename _Traits>
2946 std::basic_ostream<_CharT, _Traits>&
2947 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2948 const piecewise_constant_distribution<_RealType>& __x)
2950 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2951 typedef typename __ostream_type::ios_base __ios_base;
2953 const typename __ios_base::fmtflags __flags = __os.flags();
2954 const _CharT __fill = __os.fill();
2955 const std::streamsize __precision = __os.precision();
2956 const _CharT __space = __os.widen(' ');
2957 __os.flags(__ios_base::scientific | __ios_base::left);
2959 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2961 std::vector<_RealType> __int = __x.intervals();
2962 __os << __int.size() - 1;
2964 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2965 __os << __space << *__xit;
2967 std::vector<double> __den = __x.densities();
2968 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2969 __os << __space << *__dit;
2971 __os.flags(__flags);
2973 __os.precision(__precision);
2977 template<typename _RealType, typename _CharT, typename _Traits>
2978 std::basic_istream<_CharT, _Traits>&
2979 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2980 piecewise_constant_distribution<_RealType>& __x)
2982 typedef std::basic_istream<_CharT, _Traits> __istream_type;
2983 typedef typename __istream_type::ios_base __ios_base;
2985 const typename __ios_base::fmtflags __flags = __is.flags();
2986 __is.flags(__ios_base::dec | __ios_base::skipws);
2991 std::vector<_RealType> __int_vec;
2992 __int_vec.reserve(__n + 1);
2993 for (size_t __i = 0; __i <= __n; ++__i)
2997 __int_vec.push_back(__int);
3000 std::vector<double> __den_vec;
3001 __den_vec.reserve(__n);
3002 for (size_t __i = 0; __i < __n; ++__i)
3006 __den_vec.push_back(__den);
3009 __x.param(typename piecewise_constant_distribution<_RealType>::
3010 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3012 __is.flags(__flags);
3017 template<typename _RealType>
3019 piecewise_linear_distribution<_RealType>::param_type::
3022 if (_M_int.size() < 2
3023 || (_M_int.size() == 2
3024 && _M_int[0] == _RealType(0)
3025 && _M_int[1] == _RealType(1)
3026 && _M_den[0] == _M_den[1]))
3034 _M_cp.reserve(_M_int.size() - 1);
3035 _M_m.reserve(_M_int.size() - 1);
3036 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3038 const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3039 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3040 _M_cp.push_back(__sum);
3041 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3044 // Now normalize the densities...
3045 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3047 // ... and partial sums...
3048 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3050 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3052 // Make sure the last cumulative probablility is one.
3053 _M_cp[_M_cp.size() - 1] = 1.0;
3056 template<typename _RealType>
3057 template<typename _InputIteratorB, typename _InputIteratorW>
3058 piecewise_linear_distribution<_RealType>::param_type::
3059 param_type(_InputIteratorB __bbegin,
3060 _InputIteratorB __bend,
3061 _InputIteratorW __wbegin)
3062 : _M_int(), _M_den(), _M_cp(), _M_m()
3064 for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
3066 _M_int.push_back(*__bbegin);
3067 _M_den.push_back(*__wbegin);
3073 template<typename _RealType>
3074 template<typename _Func>
3075 piecewise_linear_distribution<_RealType>::param_type::
3076 param_type(initializer_list<_RealType> __bl, _Func __fw)
3077 : _M_int(), _M_den(), _M_cp(), _M_m()
3079 _M_int.reserve(__bl.size());
3080 _M_den.reserve(__bl.size());
3081 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3083 _M_int.push_back(*__biter);
3084 _M_den.push_back(__fw(*__biter));
3090 template<typename _RealType>
3091 template<typename _Func>
3092 piecewise_linear_distribution<_RealType>::param_type::
3093 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3094 : _M_int(), _M_den(), _M_cp(), _M_m()
3096 const size_t __n = __nw == 0 ? 1 : __nw;
3097 const _RealType __delta = (__xmax - __xmin) / __n;
3099 _M_int.reserve(__n + 1);
3100 _M_den.reserve(__n + 1);
3101 for (size_t __k = 0; __k <= __nw; ++__k)
3103 _M_int.push_back(__xmin + __k * __delta);
3104 _M_den.push_back(__fw(_M_int[__k] + __delta));
3110 template<typename _RealType>
3111 template<typename _UniformRandomNumberGenerator>
3112 typename piecewise_linear_distribution<_RealType>::result_type
3113 piecewise_linear_distribution<_RealType>::
3114 operator()(_UniformRandomNumberGenerator& __urng,
3115 const param_type& __param)
3117 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3120 const double __p = __aurng();
3121 if (__param._M_cp.empty())
3124 auto __pos = std::lower_bound(__param._M_cp.begin(),
3125 __param._M_cp.end(), __p);
3126 const size_t __i = __pos - __param._M_cp.begin();
3128 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3130 const double __a = 0.5 * __param._M_m[__i];
3131 const double __b = __param._M_den[__i];
3132 const double __cm = __p - __pref;
3134 _RealType __x = __param._M_int[__i];
3139 const double __d = __b * __b + 4.0 * __a * __cm;
3140 __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3146 template<typename _RealType>
3147 template<typename _ForwardIterator,
3148 typename _UniformRandomNumberGenerator>
3150 piecewise_linear_distribution<_RealType>::
3151 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3152 _UniformRandomNumberGenerator& __urng,
3153 const param_type& __param)
3155 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3156 // We could duplicate everything from operator()...
3158 *__f++ = this->operator()(__urng, __param);
3161 template<typename _RealType, typename _CharT, typename _Traits>
3162 std::basic_ostream<_CharT, _Traits>&
3163 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3164 const piecewise_linear_distribution<_RealType>& __x)
3166 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
3167 typedef typename __ostream_type::ios_base __ios_base;
3169 const typename __ios_base::fmtflags __flags = __os.flags();
3170 const _CharT __fill = __os.fill();
3171 const std::streamsize __precision = __os.precision();
3172 const _CharT __space = __os.widen(' ');
3173 __os.flags(__ios_base::scientific | __ios_base::left);
3175 __os.precision(std::numeric_limits<_RealType>::max_digits10);
3177 std::vector<_RealType> __int = __x.intervals();
3178 __os << __int.size() - 1;
3180 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3181 __os << __space << *__xit;
3183 std::vector<double> __den = __x.densities();
3184 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3185 __os << __space << *__dit;
3187 __os.flags(__flags);
3189 __os.precision(__precision);
3193 template<typename _RealType, typename _CharT, typename _Traits>
3194 std::basic_istream<_CharT, _Traits>&
3195 operator>>(std::basic_istream<_CharT, _Traits>& __is,
3196 piecewise_linear_distribution<_RealType>& __x)
3198 typedef std::basic_istream<_CharT, _Traits> __istream_type;
3199 typedef typename __istream_type::ios_base __ios_base;
3201 const typename __ios_base::fmtflags __flags = __is.flags();
3202 __is.flags(__ios_base::dec | __ios_base::skipws);
3207 std::vector<_RealType> __int_vec;
3208 __int_vec.reserve(__n + 1);
3209 for (size_t __i = 0; __i <= __n; ++__i)
3213 __int_vec.push_back(__int);
3216 std::vector<double> __den_vec;
3217 __den_vec.reserve(__n + 1);
3218 for (size_t __i = 0; __i <= __n; ++__i)
3222 __den_vec.push_back(__den);
3225 __x.param(typename piecewise_linear_distribution<_RealType>::
3226 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3228 __is.flags(__flags);
3233 template<typename _IntType>
3234 seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3236 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3237 _M_v.push_back(__detail::__mod<result_type,
3238 __detail::_Shift<result_type, 32>::__value>(*__iter));
3241 template<typename _InputIterator>
3242 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3244 for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3245 _M_v.push_back(__detail::__mod<result_type,
3246 __detail::_Shift<result_type, 32>::__value>(*__iter));
3249 template<typename _RandomAccessIterator>
3251 seed_seq::generate(_RandomAccessIterator __begin,
3252 _RandomAccessIterator __end)
3254 typedef typename iterator_traits<_RandomAccessIterator>::value_type
3257 if (__begin == __end)
3260 std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3262 const size_t __n = __end - __begin;
3263 const size_t __s = _M_v.size();
3264 const size_t __t = (__n >= 623) ? 11
3269 const size_t __p = (__n - __t) / 2;
3270 const size_t __q = __p + __t;
3271 const size_t __m = std::max(size_t(__s + 1), __n);
3273 for (size_t __k = 0; __k < __m; ++__k)
3275 _Type __arg = (__begin[__k % __n]
3276 ^ __begin[(__k + __p) % __n]
3277 ^ __begin[(__k - 1) % __n]);
3278 _Type __r1 = __arg ^ (__arg >> 27);
3279 __r1 = __detail::__mod<_Type,
3280 __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
3284 else if (__k <= __s)
3285 __r2 += __k % __n + _M_v[__k - 1];
3288 __r2 = __detail::__mod<_Type,
3289 __detail::_Shift<_Type, 32>::__value>(__r2);
3290 __begin[(__k + __p) % __n] += __r1;
3291 __begin[(__k + __q) % __n] += __r2;
3292 __begin[__k % __n] = __r2;
3295 for (size_t __k = __m; __k < __m + __n; ++__k)
3297 _Type __arg = (__begin[__k % __n]
3298 + __begin[(__k + __p) % __n]
3299 + __begin[(__k - 1) % __n]);
3300 _Type __r3 = __arg ^ (__arg >> 27);
3301 __r3 = __detail::__mod<_Type,
3302 __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
3303 _Type __r4 = __r3 - __k % __n;
3304 __r4 = __detail::__mod<_Type,
3305 __detail::_Shift<_Type, 32>::__value>(__r4);
3306 __begin[(__k + __p) % __n] ^= __r3;
3307 __begin[(__k + __q) % __n] ^= __r4;
3308 __begin[__k % __n] = __r4;
3312 template<typename _RealType, size_t __bits,
3313 typename _UniformRandomNumberGenerator>
3315 generate_canonical(_UniformRandomNumberGenerator& __urng)
3317 static_assert(std::is_floating_point<_RealType>::value,
3318 "template argument must be a floating point type");
3321 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3323 const long double __r = static_cast<long double>(__urng.max())
3324 - static_cast<long double>(__urng.min()) + 1.0L;
3325 const size_t __log2r = std::log(__r) / std::log(2.0L);
3326 const size_t __m = std::max<size_t>(1UL,
3327 (__b + __log2r - 1UL) / __log2r);
3329 _RealType __sum = _RealType(0);
3330 _RealType __tmp = _RealType(1);
3331 for (size_t __k = __m; __k != 0; --__k)
3333 __sum += _RealType(__urng() - __urng.min()) * __tmp;
3336 __ret = __sum / __tmp;
3337 if (__builtin_expect(__ret >= _RealType(1), 0))
3339 #if _GLIBCXX_USE_C99_MATH_TR1
3340 __ret = std::nextafter(_RealType(1), _RealType(0));
3342 __ret = _RealType(1)
3343 - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3349 _GLIBCXX_END_NAMESPACE_VERSION