1 // random number generation (out of line) -*- C++ -*-
3 // Copyright (C) 2009-2022 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
39 /// @cond undocumented
40 // (Further) implementation-space details.
43 // General case for x = (ax + c) mod m -- use Schrage's algorithm
44 // to avoid integer overflow.
46 // Preconditions: a > 0, m > 0.
48 // Note: only works correctly for __m % __a < __m / __a.
49 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
51 _Mod<_Tp, __m, __a, __c, false, true>::
58 static const _Tp __q = __m / __a;
59 static const _Tp __r = __m % __a;
61 _Tp __t1 = __a * (__x % __q);
62 _Tp __t2 = __r * (__x / __q);
66 __x = __m - __t2 + __t1;
71 const _Tp __d = __m - __x;
80 template<typename _InputIterator, typename _OutputIterator,
83 __normalize(_InputIterator __first, _InputIterator __last,
84 _OutputIterator __result, const _Tp& __factor)
86 for (; __first != __last; ++__first, ++__result)
87 *__result = *__first / __factor;
91 } // namespace __detail
94 #if ! __cpp_inline_variables
95 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
97 linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
99 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
101 linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
103 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
105 linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
107 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
109 linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
113 * Seeds the LCR with integral value @p __s, adjusted so that the
114 * ring identity is never a member of the convergence set.
116 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
118 linear_congruential_engine<_UIntType, __a, __c, __m>::
119 seed(result_type __s)
121 if ((__detail::__mod<_UIntType, __m>(__c) == 0)
122 && (__detail::__mod<_UIntType, __m>(__s) == 0))
125 _M_x = __detail::__mod<_UIntType, __m>(__s);
129 * Seeds the LCR engine with a value generated by @p __q.
131 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
132 template<typename _Sseq>
134 linear_congruential_engine<_UIntType, __a, __c, __m>::
136 -> _If_seed_seq<_Sseq>
138 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
140 const _UIntType __k = (__k0 + 31) / 32;
141 uint_least32_t __arr[__k + 3];
142 __q.generate(__arr + 0, __arr + __k + 3);
143 _UIntType __factor = 1u;
144 _UIntType __sum = 0u;
145 for (size_t __j = 0; __j < __k; ++__j)
147 __sum += __arr[__j + 3] * __factor;
148 __factor *= __detail::_Shift<_UIntType, 32>::__value;
153 template<typename _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_engine<_UIntType,
158 __a, __c, __m>& __lcr)
160 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
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(' '));
174 template<typename _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_engine<_UIntType, __a, __c, __m>& __lcr)
180 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
182 const typename __ios_base::fmtflags __flags = __is.flags();
183 __is.flags(__ios_base::dec);
191 #if ! __cpp_inline_variables
192 template<typename _UIntType,
193 size_t __w, size_t __n, size_t __m, size_t __r,
194 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
195 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
198 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
199 __s, __b, __t, __c, __l, __f>::word_size;
201 template<typename _UIntType,
202 size_t __w, size_t __n, size_t __m, size_t __r,
203 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
204 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
207 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
208 __s, __b, __t, __c, __l, __f>::state_size;
210 template<typename _UIntType,
211 size_t __w, size_t __n, size_t __m, size_t __r,
212 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
213 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
216 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
217 __s, __b, __t, __c, __l, __f>::shift_size;
219 template<typename _UIntType,
220 size_t __w, size_t __n, size_t __m, size_t __r,
221 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
222 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
225 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
226 __s, __b, __t, __c, __l, __f>::mask_bits;
228 template<typename _UIntType,
229 size_t __w, size_t __n, size_t __m, size_t __r,
230 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
231 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
234 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
235 __s, __b, __t, __c, __l, __f>::xor_mask;
237 template<typename _UIntType,
238 size_t __w, size_t __n, size_t __m, size_t __r,
239 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
240 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
243 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
244 __s, __b, __t, __c, __l, __f>::tempering_u;
246 template<typename _UIntType,
247 size_t __w, size_t __n, size_t __m, size_t __r,
248 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
249 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
252 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
253 __s, __b, __t, __c, __l, __f>::tempering_d;
255 template<typename _UIntType,
256 size_t __w, size_t __n, size_t __m, size_t __r,
257 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
258 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
261 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
262 __s, __b, __t, __c, __l, __f>::tempering_s;
264 template<typename _UIntType,
265 size_t __w, size_t __n, size_t __m, size_t __r,
266 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
267 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
270 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
271 __s, __b, __t, __c, __l, __f>::tempering_b;
273 template<typename _UIntType,
274 size_t __w, size_t __n, size_t __m, size_t __r,
275 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
276 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
279 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
280 __s, __b, __t, __c, __l, __f>::tempering_t;
282 template<typename _UIntType,
283 size_t __w, size_t __n, size_t __m, size_t __r,
284 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
285 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
288 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
289 __s, __b, __t, __c, __l, __f>::tempering_c;
291 template<typename _UIntType,
292 size_t __w, size_t __n, size_t __m, size_t __r,
293 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
294 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
297 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
298 __s, __b, __t, __c, __l, __f>::tempering_l;
300 template<typename _UIntType,
301 size_t __w, size_t __n, size_t __m, size_t __r,
302 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
303 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
306 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
307 __s, __b, __t, __c, __l, __f>::
308 initialization_multiplier;
310 template<typename _UIntType,
311 size_t __w, size_t __n, size_t __m, size_t __r,
312 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
313 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
316 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
317 __s, __b, __t, __c, __l, __f>::default_seed;
320 template<typename _UIntType,
321 size_t __w, size_t __n, size_t __m, size_t __r,
322 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
323 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
326 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
327 __s, __b, __t, __c, __l, __f>::
328 seed(result_type __sd)
330 _M_x[0] = __detail::__mod<_UIntType,
331 __detail::_Shift<_UIntType, __w>::__value>(__sd);
333 for (size_t __i = 1; __i < state_size; ++__i)
335 _UIntType __x = _M_x[__i - 1];
336 __x ^= __x >> (__w - 2);
338 __x += __detail::__mod<_UIntType, __n>(__i);
339 _M_x[__i] = __detail::__mod<_UIntType,
340 __detail::_Shift<_UIntType, __w>::__value>(__x);
345 template<typename _UIntType,
346 size_t __w, size_t __n, size_t __m, size_t __r,
347 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
348 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
350 template<typename _Sseq>
352 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
353 __s, __b, __t, __c, __l, __f>::
355 -> _If_seed_seq<_Sseq>
357 const _UIntType __upper_mask = (~_UIntType()) << __r;
358 const size_t __k = (__w + 31) / 32;
359 uint_least32_t __arr[__n * __k];
360 __q.generate(__arr + 0, __arr + __n * __k);
363 for (size_t __i = 0; __i < state_size; ++__i)
365 _UIntType __factor = 1u;
366 _UIntType __sum = 0u;
367 for (size_t __j = 0; __j < __k; ++__j)
369 __sum += __arr[__k * __i + __j] * __factor;
370 __factor *= __detail::_Shift<_UIntType, 32>::__value;
372 _M_x[__i] = __detail::__mod<_UIntType,
373 __detail::_Shift<_UIntType, __w>::__value>(__sum);
379 if ((_M_x[0] & __upper_mask) != 0u)
382 else if (_M_x[__i] != 0u)
387 _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
391 template<typename _UIntType, size_t __w,
392 size_t __n, size_t __m, size_t __r,
393 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
394 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
397 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
398 __s, __b, __t, __c, __l, __f>::
401 const _UIntType __upper_mask = (~_UIntType()) << __r;
402 const _UIntType __lower_mask = ~__upper_mask;
404 for (size_t __k = 0; __k < (__n - __m); ++__k)
406 _UIntType __y = ((_M_x[__k] & __upper_mask)
407 | (_M_x[__k + 1] & __lower_mask));
408 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
409 ^ ((__y & 0x01) ? __a : 0));
412 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
414 _UIntType __y = ((_M_x[__k] & __upper_mask)
415 | (_M_x[__k + 1] & __lower_mask));
416 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
417 ^ ((__y & 0x01) ? __a : 0));
420 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
421 | (_M_x[0] & __lower_mask));
422 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
423 ^ ((__y & 0x01) ? __a : 0));
427 template<typename _UIntType, size_t __w,
428 size_t __n, size_t __m, size_t __r,
429 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
430 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
433 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
434 __s, __b, __t, __c, __l, __f>::
435 discard(unsigned long long __z)
437 while (__z > state_size - _M_p)
439 __z -= state_size - _M_p;
445 template<typename _UIntType, size_t __w,
446 size_t __n, size_t __m, size_t __r,
447 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
448 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
451 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
452 __s, __b, __t, __c, __l, __f>::result_type
453 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
454 __s, __b, __t, __c, __l, __f>::
457 // Reload the vector - cost is O(n) amortized over n calls.
458 if (_M_p >= state_size)
461 // Calculate o(x(i)).
462 result_type __z = _M_x[_M_p++];
463 __z ^= (__z >> __u) & __d;
464 __z ^= (__z << __s) & __b;
465 __z ^= (__z << __t) & __c;
471 template<typename _UIntType, size_t __w,
472 size_t __n, size_t __m, size_t __r,
473 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
474 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
475 _UIntType __f, typename _CharT, typename _Traits>
476 std::basic_ostream<_CharT, _Traits>&
477 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
478 const mersenne_twister_engine<_UIntType, __w, __n, __m,
479 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
481 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
483 const typename __ios_base::fmtflags __flags = __os.flags();
484 const _CharT __fill = __os.fill();
485 const _CharT __space = __os.widen(' ');
486 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
489 for (size_t __i = 0; __i < __n; ++__i)
490 __os << __x._M_x[__i] << __space;
498 template<typename _UIntType, size_t __w,
499 size_t __n, size_t __m, size_t __r,
500 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
501 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
502 _UIntType __f, typename _CharT, typename _Traits>
503 std::basic_istream<_CharT, _Traits>&
504 operator>>(std::basic_istream<_CharT, _Traits>& __is,
505 mersenne_twister_engine<_UIntType, __w, __n, __m,
506 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
508 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
510 const typename __ios_base::fmtflags __flags = __is.flags();
511 __is.flags(__ios_base::dec | __ios_base::skipws);
513 for (size_t __i = 0; __i < __n; ++__i)
514 __is >> __x._M_x[__i];
521 #if ! __cpp_inline_variables
522 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
524 subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
526 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
528 subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
530 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
532 subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
534 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
536 subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
539 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
541 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
542 seed(result_type __value)
544 std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
545 __lcg(__value == 0u ? default_seed : __value);
547 const size_t __n = (__w + 31) / 32;
549 for (size_t __i = 0; __i < long_lag; ++__i)
551 _UIntType __sum = 0u;
552 _UIntType __factor = 1u;
553 for (size_t __j = 0; __j < __n; ++__j)
555 __sum += __detail::__mod<uint_least32_t,
556 __detail::_Shift<uint_least32_t, 32>::__value>
557 (__lcg()) * __factor;
558 __factor *= __detail::_Shift<_UIntType, 32>::__value;
560 _M_x[__i] = __detail::__mod<_UIntType,
561 __detail::_Shift<_UIntType, __w>::__value>(__sum);
563 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
567 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
568 template<typename _Sseq>
570 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
572 -> _If_seed_seq<_Sseq>
574 const size_t __k = (__w + 31) / 32;
575 uint_least32_t __arr[__r * __k];
576 __q.generate(__arr + 0, __arr + __r * __k);
578 for (size_t __i = 0; __i < long_lag; ++__i)
580 _UIntType __sum = 0u;
581 _UIntType __factor = 1u;
582 for (size_t __j = 0; __j < __k; ++__j)
584 __sum += __arr[__k * __i + __j] * __factor;
585 __factor *= __detail::_Shift<_UIntType, 32>::__value;
587 _M_x[__i] = __detail::__mod<_UIntType,
588 __detail::_Shift<_UIntType, __w>::__value>(__sum);
590 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
594 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
595 typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
597 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
600 // Derive short lag index from current index.
601 long __ps = _M_p - short_lag;
605 // Calculate new x(i) without overflow or division.
606 // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
609 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
611 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
616 __xi = (__detail::_Shift<_UIntType, __w>::__value
617 - _M_x[_M_p] - _M_carry + _M_x[__ps]);
622 // Adjust current index to loop around in ring buffer.
623 if (++_M_p >= long_lag)
629 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
630 typename _CharT, typename _Traits>
631 std::basic_ostream<_CharT, _Traits>&
632 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
633 const subtract_with_carry_engine<_UIntType,
636 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
638 const typename __ios_base::fmtflags __flags = __os.flags();
639 const _CharT __fill = __os.fill();
640 const _CharT __space = __os.widen(' ');
641 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
644 for (size_t __i = 0; __i < __r; ++__i)
645 __os << __x._M_x[__i] << __space;
646 __os << __x._M_carry << __space << __x._M_p;
653 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
654 typename _CharT, typename _Traits>
655 std::basic_istream<_CharT, _Traits>&
656 operator>>(std::basic_istream<_CharT, _Traits>& __is,
657 subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
659 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
661 const typename __ios_base::fmtflags __flags = __is.flags();
662 __is.flags(__ios_base::dec | __ios_base::skipws);
664 for (size_t __i = 0; __i < __r; ++__i)
665 __is >> __x._M_x[__i];
666 __is >> __x._M_carry;
673 #if ! __cpp_inline_variables
674 template<typename _RandomNumberEngine, size_t __p, size_t __r>
676 discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
678 template<typename _RandomNumberEngine, size_t __p, size_t __r>
680 discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
683 template<typename _RandomNumberEngine, size_t __p, size_t __r>
684 typename discard_block_engine<_RandomNumberEngine,
685 __p, __r>::result_type
686 discard_block_engine<_RandomNumberEngine, __p, __r>::
689 if (_M_n >= used_block)
691 _M_b.discard(block_size - _M_n);
698 template<typename _RandomNumberEngine, size_t __p, size_t __r,
699 typename _CharT, typename _Traits>
700 std::basic_ostream<_CharT, _Traits>&
701 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
702 const discard_block_engine<_RandomNumberEngine,
705 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
707 const typename __ios_base::fmtflags __flags = __os.flags();
708 const _CharT __fill = __os.fill();
709 const _CharT __space = __os.widen(' ');
710 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
713 __os << __x.base() << __space << __x._M_n;
720 template<typename _RandomNumberEngine, size_t __p, size_t __r,
721 typename _CharT, typename _Traits>
722 std::basic_istream<_CharT, _Traits>&
723 operator>>(std::basic_istream<_CharT, _Traits>& __is,
724 discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
726 using __ios_base = typename basic_istream<_CharT, _Traits>::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);
807 #if ! __cpp_inline_variables
808 template<typename _RandomNumberEngine, size_t __k>
810 shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
815 // Determine whether an integer is representable as double.
816 template<typename _Tp>
818 __representable_as_double(_Tp __x) noexcept
820 static_assert(numeric_limits<_Tp>::is_integer, "");
821 static_assert(!numeric_limits<_Tp>::is_signed, "");
822 // All integers <= 2^53 are representable.
823 return (__x <= (1ull << __DBL_MANT_DIG__))
824 // Between 2^53 and 2^54 only even numbers are representable.
825 || (!(__x & 1) && __detail::__representable_as_double(__x >> 1));
828 // Determine whether x+1 is representable as double.
829 template<typename _Tp>
831 __p1_representable_as_double(_Tp __x) noexcept
833 static_assert(numeric_limits<_Tp>::is_integer, "");
834 static_assert(!numeric_limits<_Tp>::is_signed, "");
835 return numeric_limits<_Tp>::digits < __DBL_MANT_DIG__
836 || (bool(__x + 1u) // return false if x+1 wraps around to zero
837 && __detail::__representable_as_double(__x + 1u));
841 template<typename _RandomNumberEngine, size_t __k>
842 typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
843 shuffle_order_engine<_RandomNumberEngine, __k>::
846 constexpr result_type __range = max() - min();
848 const result_type __y = _M_y - min();
849 // Avoid using slower long double arithmetic if possible.
850 if _GLIBCXX17_CONSTEXPR (__detail::__p1_representable_as_double(__range))
851 __j *= __y / (__range + 1.0);
853 __j *= __y / (__range + 1.0L);
860 template<typename _RandomNumberEngine, size_t __k,
861 typename _CharT, typename _Traits>
862 std::basic_ostream<_CharT, _Traits>&
863 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
864 const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
866 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
868 const typename __ios_base::fmtflags __flags = __os.flags();
869 const _CharT __fill = __os.fill();
870 const _CharT __space = __os.widen(' ');
871 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
875 for (size_t __i = 0; __i < __k; ++__i)
876 __os << __space << __x._M_v[__i];
877 __os << __space << __x._M_y;
884 template<typename _RandomNumberEngine, size_t __k,
885 typename _CharT, typename _Traits>
886 std::basic_istream<_CharT, _Traits>&
887 operator>>(std::basic_istream<_CharT, _Traits>& __is,
888 shuffle_order_engine<_RandomNumberEngine, __k>& __x)
890 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
892 const typename __ios_base::fmtflags __flags = __is.flags();
893 __is.flags(__ios_base::dec | __ios_base::skipws);
896 for (size_t __i = 0; __i < __k; ++__i)
897 __is >> __x._M_v[__i];
905 template<typename _IntType, typename _CharT, typename _Traits>
906 std::basic_ostream<_CharT, _Traits>&
907 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
908 const uniform_int_distribution<_IntType>& __x)
910 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
912 const typename __ios_base::fmtflags __flags = __os.flags();
913 const _CharT __fill = __os.fill();
914 const _CharT __space = __os.widen(' ');
915 __os.flags(__ios_base::scientific | __ios_base::left);
918 __os << __x.a() << __space << __x.b();
925 template<typename _IntType, typename _CharT, typename _Traits>
926 std::basic_istream<_CharT, _Traits>&
927 operator>>(std::basic_istream<_CharT, _Traits>& __is,
928 uniform_int_distribution<_IntType>& __x)
931 = typename uniform_int_distribution<_IntType>::param_type;
932 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
934 const typename __ios_base::fmtflags __flags = __is.flags();
935 __is.flags(__ios_base::dec | __ios_base::skipws);
938 if (__is >> __a >> __b)
939 __x.param(param_type(__a, __b));
946 template<typename _RealType>
947 template<typename _ForwardIterator,
948 typename _UniformRandomNumberGenerator>
950 uniform_real_distribution<_RealType>::
951 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
952 _UniformRandomNumberGenerator& __urng,
953 const param_type& __p)
955 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
956 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
958 auto __range = __p.b() - __p.a();
960 *__f++ = __aurng() * __range + __p.a();
963 template<typename _RealType, typename _CharT, typename _Traits>
964 std::basic_ostream<_CharT, _Traits>&
965 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
966 const uniform_real_distribution<_RealType>& __x)
968 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
970 const typename __ios_base::fmtflags __flags = __os.flags();
971 const _CharT __fill = __os.fill();
972 const std::streamsize __precision = __os.precision();
973 const _CharT __space = __os.widen(' ');
974 __os.flags(__ios_base::scientific | __ios_base::left);
976 __os.precision(std::numeric_limits<_RealType>::max_digits10);
978 __os << __x.a() << __space << __x.b();
982 __os.precision(__precision);
986 template<typename _RealType, typename _CharT, typename _Traits>
987 std::basic_istream<_CharT, _Traits>&
988 operator>>(std::basic_istream<_CharT, _Traits>& __is,
989 uniform_real_distribution<_RealType>& __x)
992 = typename uniform_real_distribution<_RealType>::param_type;
993 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
995 const typename __ios_base::fmtflags __flags = __is.flags();
996 __is.flags(__ios_base::skipws);
999 if (__is >> __a >> __b)
1000 __x.param(param_type(__a, __b));
1002 __is.flags(__flags);
1007 template<typename _ForwardIterator,
1008 typename _UniformRandomNumberGenerator>
1010 std::bernoulli_distribution::
1011 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1012 _UniformRandomNumberGenerator& __urng,
1013 const param_type& __p)
1015 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1016 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1018 auto __limit = __p.p() * (__aurng.max() - __aurng.min());
1021 *__f++ = (__aurng() - __aurng.min()) < __limit;
1024 template<typename _CharT, typename _Traits>
1025 std::basic_ostream<_CharT, _Traits>&
1026 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1027 const bernoulli_distribution& __x)
1029 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1031 const typename __ios_base::fmtflags __flags = __os.flags();
1032 const _CharT __fill = __os.fill();
1033 const std::streamsize __precision = __os.precision();
1034 __os.flags(__ios_base::scientific | __ios_base::left);
1035 __os.fill(__os.widen(' '));
1036 __os.precision(std::numeric_limits<double>::max_digits10);
1040 __os.flags(__flags);
1042 __os.precision(__precision);
1047 template<typename _IntType>
1048 template<typename _UniformRandomNumberGenerator>
1049 typename geometric_distribution<_IntType>::result_type
1050 geometric_distribution<_IntType>::
1051 operator()(_UniformRandomNumberGenerator& __urng,
1052 const param_type& __param)
1054 // About the epsilon thing see this thread:
1055 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1056 const double __naf =
1057 (1 - std::numeric_limits<double>::epsilon()) / 2;
1058 // The largest _RealType convertible to _IntType.
1059 const double __thr =
1060 std::numeric_limits<_IntType>::max() + __naf;
1061 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1066 __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1067 while (__cand >= __thr);
1069 return result_type(__cand + __naf);
1072 template<typename _IntType>
1073 template<typename _ForwardIterator,
1074 typename _UniformRandomNumberGenerator>
1076 geometric_distribution<_IntType>::
1077 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1078 _UniformRandomNumberGenerator& __urng,
1079 const param_type& __param)
1081 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1082 // About the epsilon thing see this thread:
1083 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1084 const double __naf =
1085 (1 - std::numeric_limits<double>::epsilon()) / 2;
1086 // The largest _RealType convertible to _IntType.
1087 const double __thr =
1088 std::numeric_limits<_IntType>::max() + __naf;
1089 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1096 __cand = std::floor(std::log(1.0 - __aurng())
1097 / __param._M_log_1_p);
1098 while (__cand >= __thr);
1100 *__f++ = __cand + __naf;
1104 template<typename _IntType,
1105 typename _CharT, typename _Traits>
1106 std::basic_ostream<_CharT, _Traits>&
1107 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1108 const geometric_distribution<_IntType>& __x)
1110 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1112 const typename __ios_base::fmtflags __flags = __os.flags();
1113 const _CharT __fill = __os.fill();
1114 const std::streamsize __precision = __os.precision();
1115 __os.flags(__ios_base::scientific | __ios_base::left);
1116 __os.fill(__os.widen(' '));
1117 __os.precision(std::numeric_limits<double>::max_digits10);
1121 __os.flags(__flags);
1123 __os.precision(__precision);
1127 template<typename _IntType,
1128 typename _CharT, typename _Traits>
1129 std::basic_istream<_CharT, _Traits>&
1130 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1131 geometric_distribution<_IntType>& __x)
1133 using param_type = typename geometric_distribution<_IntType>::param_type;
1134 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1136 const typename __ios_base::fmtflags __flags = __is.flags();
1137 __is.flags(__ios_base::skipws);
1141 __x.param(param_type(__p));
1143 __is.flags(__flags);
1147 // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1148 template<typename _IntType>
1149 template<typename _UniformRandomNumberGenerator>
1150 typename negative_binomial_distribution<_IntType>::result_type
1151 negative_binomial_distribution<_IntType>::
1152 operator()(_UniformRandomNumberGenerator& __urng)
1154 const double __y = _M_gd(__urng);
1156 // XXX Is the constructor too slow?
1157 std::poisson_distribution<result_type> __poisson(__y);
1158 return __poisson(__urng);
1161 template<typename _IntType>
1162 template<typename _UniformRandomNumberGenerator>
1163 typename negative_binomial_distribution<_IntType>::result_type
1164 negative_binomial_distribution<_IntType>::
1165 operator()(_UniformRandomNumberGenerator& __urng,
1166 const param_type& __p)
1168 typedef typename std::gamma_distribution<double>::param_type
1172 _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1174 std::poisson_distribution<result_type> __poisson(__y);
1175 return __poisson(__urng);
1178 template<typename _IntType>
1179 template<typename _ForwardIterator,
1180 typename _UniformRandomNumberGenerator>
1182 negative_binomial_distribution<_IntType>::
1183 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1184 _UniformRandomNumberGenerator& __urng)
1186 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1189 const double __y = _M_gd(__urng);
1191 // XXX Is the constructor too slow?
1192 std::poisson_distribution<result_type> __poisson(__y);
1193 *__f++ = __poisson(__urng);
1197 template<typename _IntType>
1198 template<typename _ForwardIterator,
1199 typename _UniformRandomNumberGenerator>
1201 negative_binomial_distribution<_IntType>::
1202 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1203 _UniformRandomNumberGenerator& __urng,
1204 const param_type& __p)
1206 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1207 typename std::gamma_distribution<result_type>::param_type
1208 __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1212 const double __y = _M_gd(__urng, __p2);
1214 std::poisson_distribution<result_type> __poisson(__y);
1215 *__f++ = __poisson(__urng);
1219 template<typename _IntType, typename _CharT, typename _Traits>
1220 std::basic_ostream<_CharT, _Traits>&
1221 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1222 const negative_binomial_distribution<_IntType>& __x)
1224 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1226 const typename __ios_base::fmtflags __flags = __os.flags();
1227 const _CharT __fill = __os.fill();
1228 const std::streamsize __precision = __os.precision();
1229 const _CharT __space = __os.widen(' ');
1230 __os.flags(__ios_base::scientific | __ios_base::left);
1231 __os.fill(__os.widen(' '));
1232 __os.precision(std::numeric_limits<double>::max_digits10);
1234 __os << __x.k() << __space << __x.p()
1235 << __space << __x._M_gd;
1237 __os.flags(__flags);
1239 __os.precision(__precision);
1243 template<typename _IntType, typename _CharT, typename _Traits>
1244 std::basic_istream<_CharT, _Traits>&
1245 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1246 negative_binomial_distribution<_IntType>& __x)
1249 = typename negative_binomial_distribution<_IntType>::param_type;
1250 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1252 const typename __ios_base::fmtflags __flags = __is.flags();
1253 __is.flags(__ios_base::skipws);
1257 if (__is >> __k >> __p >> __x._M_gd)
1258 __x.param(param_type(__k, __p));
1260 __is.flags(__flags);
1265 template<typename _IntType>
1267 poisson_distribution<_IntType>::param_type::
1270 #if _GLIBCXX_USE_C99_MATH_TR1
1273 const double __m = std::floor(_M_mean);
1274 _M_lm_thr = std::log(_M_mean);
1275 _M_lfm = std::lgamma(__m + 1);
1276 _M_sm = std::sqrt(__m);
1278 const double __pi_4 = 0.7853981633974483096156608458198757L;
1279 const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1281 _M_d = std::round(std::max<double>(6.0, std::min(__m, __dx)));
1282 const double __cx = 2 * __m + _M_d;
1283 _M_scx = std::sqrt(__cx / 2);
1286 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1287 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1292 _M_lm_thr = std::exp(-_M_mean);
1296 * A rejection algorithm when mean >= 12 and a simple method based
1297 * upon the multiplication of uniform random variates otherwise.
1298 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1302 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1303 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1305 template<typename _IntType>
1306 template<typename _UniformRandomNumberGenerator>
1307 typename poisson_distribution<_IntType>::result_type
1308 poisson_distribution<_IntType>::
1309 operator()(_UniformRandomNumberGenerator& __urng,
1310 const param_type& __param)
1312 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1314 #if _GLIBCXX_USE_C99_MATH_TR1
1315 if (__param.mean() >= 12)
1319 // See comments above...
1320 const double __naf =
1321 (1 - std::numeric_limits<double>::epsilon()) / 2;
1322 const double __thr =
1323 std::numeric_limits<_IntType>::max() + __naf;
1325 const double __m = std::floor(__param.mean());
1327 const double __spi_2 = 1.2533141373155002512078826424055226L;
1328 const double __c1 = __param._M_sm * __spi_2;
1329 const double __c2 = __param._M_c2b + __c1;
1330 const double __c3 = __c2 + 1;
1331 const double __c4 = __c3 + 1;
1333 const double __178 = 0.0128205128205128205128205128205128L;
1335 const double __e178 = 1.0129030479320018583185514777512983L;
1336 const double __c5 = __c4 + __e178;
1337 const double __c = __param._M_cb + __c5;
1338 const double __2cx = 2 * (2 * __m + __param._M_d);
1340 bool __reject = true;
1343 const double __u = __c * __aurng();
1344 const double __e = -std::log(1.0 - __aurng());
1350 const double __n = _M_nd(__urng);
1351 const double __y = -std::abs(__n) * __param._M_sm - 1;
1352 __x = std::floor(__y);
1353 __w = -__n * __n / 2;
1357 else if (__u <= __c2)
1359 const double __n = _M_nd(__urng);
1360 const double __y = 1 + std::abs(__n) * __param._M_scx;
1361 __x = std::ceil(__y);
1362 __w = __y * (2 - __y) * __param._M_1cx;
1363 if (__x > __param._M_d)
1366 else if (__u <= __c3)
1367 // NB: This case not in the book, nor in the Errata,
1368 // but should be ok...
1370 else if (__u <= __c4)
1372 else if (__u <= __c5)
1375 // Only in the Errata, see libstdc++/83237.
1380 const double __v = -std::log(1.0 - __aurng());
1381 const double __y = __param._M_d
1382 + __v * __2cx / __param._M_d;
1383 __x = std::ceil(__y);
1384 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1387 __reject = (__w - __e - __x * __param._M_lm_thr
1388 > __param._M_lfm - std::lgamma(__x + __m + 1));
1390 __reject |= __x + __m >= __thr;
1394 return result_type(__x + __m + __naf);
1400 double __prod = 1.0;
1404 __prod *= __aurng();
1407 while (__prod > __param._M_lm_thr);
1413 template<typename _IntType>
1414 template<typename _ForwardIterator,
1415 typename _UniformRandomNumberGenerator>
1417 poisson_distribution<_IntType>::
1418 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1419 _UniformRandomNumberGenerator& __urng,
1420 const param_type& __param)
1422 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1423 // We could duplicate everything from operator()...
1425 *__f++ = this->operator()(__urng, __param);
1428 template<typename _IntType,
1429 typename _CharT, typename _Traits>
1430 std::basic_ostream<_CharT, _Traits>&
1431 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1432 const poisson_distribution<_IntType>& __x)
1434 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1436 const typename __ios_base::fmtflags __flags = __os.flags();
1437 const _CharT __fill = __os.fill();
1438 const std::streamsize __precision = __os.precision();
1439 const _CharT __space = __os.widen(' ');
1440 __os.flags(__ios_base::scientific | __ios_base::left);
1442 __os.precision(std::numeric_limits<double>::max_digits10);
1444 __os << __x.mean() << __space << __x._M_nd;
1446 __os.flags(__flags);
1448 __os.precision(__precision);
1452 template<typename _IntType,
1453 typename _CharT, typename _Traits>
1454 std::basic_istream<_CharT, _Traits>&
1455 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1456 poisson_distribution<_IntType>& __x)
1458 using param_type = typename poisson_distribution<_IntType>::param_type;
1459 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1461 const typename __ios_base::fmtflags __flags = __is.flags();
1462 __is.flags(__ios_base::skipws);
1465 if (__is >> __mean >> __x._M_nd)
1466 __x.param(param_type(__mean));
1468 __is.flags(__flags);
1473 template<typename _IntType>
1475 binomial_distribution<_IntType>::param_type::
1478 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1482 #if _GLIBCXX_USE_C99_MATH_TR1
1483 if (_M_t * __p12 >= 8)
1486 const double __np = std::floor(_M_t * __p12);
1487 const double __pa = __np / _M_t;
1488 const double __1p = 1 - __pa;
1490 const double __pi_4 = 0.7853981633974483096156608458198757L;
1491 const double __d1x =
1492 std::sqrt(__np * __1p * std::log(32 * __np
1493 / (81 * __pi_4 * __1p)));
1494 _M_d1 = std::round(std::max<double>(1.0, __d1x));
1495 const double __d2x =
1496 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1497 / (__pi_4 * __pa)));
1498 _M_d2 = std::round(std::max<double>(1.0, __d2x));
1501 const double __spi_2 = 1.2533141373155002512078826424055226L;
1502 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1503 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1504 _M_c = 2 * _M_d1 / __np;
1505 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1506 const double __a12 = _M_a1 + _M_s2 * __spi_2;
1507 const double __s1s = _M_s1 * _M_s1;
1508 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1510 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1511 const double __s2s = _M_s2 * _M_s2;
1512 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1513 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1514 _M_lf = (std::lgamma(__np + 1)
1515 + std::lgamma(_M_t - __np + 1));
1516 _M_lp1p = std::log(__pa / __1p);
1518 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1522 _M_q = -std::log(1 - __p12);
1525 template<typename _IntType>
1526 template<typename _UniformRandomNumberGenerator>
1527 typename binomial_distribution<_IntType>::result_type
1528 binomial_distribution<_IntType>::
1529 _M_waiting(_UniformRandomNumberGenerator& __urng,
1530 _IntType __t, double __q)
1534 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1541 const double __e = -std::log(1.0 - __aurng());
1542 __sum += __e / (__t - __x);
1545 while (__sum <= __q);
1551 * A rejection algorithm when t * p >= 8 and a simple waiting time
1552 * method - the second in the referenced book - otherwise.
1553 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1557 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1558 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1560 template<typename _IntType>
1561 template<typename _UniformRandomNumberGenerator>
1562 typename binomial_distribution<_IntType>::result_type
1563 binomial_distribution<_IntType>::
1564 operator()(_UniformRandomNumberGenerator& __urng,
1565 const param_type& __param)
1568 const _IntType __t = __param.t();
1569 const double __p = __param.p();
1570 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1571 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1574 #if _GLIBCXX_USE_C99_MATH_TR1
1575 if (!__param._M_easy)
1579 // See comments above...
1580 const double __naf =
1581 (1 - std::numeric_limits<double>::epsilon()) / 2;
1582 const double __thr =
1583 std::numeric_limits<_IntType>::max() + __naf;
1585 const double __np = std::floor(__t * __p12);
1588 const double __spi_2 = 1.2533141373155002512078826424055226L;
1589 const double __a1 = __param._M_a1;
1590 const double __a12 = __a1 + __param._M_s2 * __spi_2;
1591 const double __a123 = __param._M_a123;
1592 const double __s1s = __param._M_s1 * __param._M_s1;
1593 const double __s2s = __param._M_s2 * __param._M_s2;
1598 const double __u = __param._M_s * __aurng();
1604 const double __n = _M_nd(__urng);
1605 const double __y = __param._M_s1 * std::abs(__n);
1606 __reject = __y >= __param._M_d1;
1609 const double __e = -std::log(1.0 - __aurng());
1610 __x = std::floor(__y);
1611 __v = -__e - __n * __n / 2 + __param._M_c;
1614 else if (__u <= __a12)
1616 const double __n = _M_nd(__urng);
1617 const double __y = __param._M_s2 * std::abs(__n);
1618 __reject = __y >= __param._M_d2;
1621 const double __e = -std::log(1.0 - __aurng());
1622 __x = std::floor(-__y);
1623 __v = -__e - __n * __n / 2;
1626 else if (__u <= __a123)
1628 const double __e1 = -std::log(1.0 - __aurng());
1629 const double __e2 = -std::log(1.0 - __aurng());
1631 const double __y = __param._M_d1
1632 + 2 * __s1s * __e1 / __param._M_d1;
1633 __x = std::floor(__y);
1634 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1635 -__y / (2 * __s1s)));
1640 const double __e1 = -std::log(1.0 - __aurng());
1641 const double __e2 = -std::log(1.0 - __aurng());
1643 const double __y = __param._M_d2
1644 + 2 * __s2s * __e1 / __param._M_d2;
1645 __x = std::floor(-__y);
1646 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1650 __reject = __reject || __x < -__np || __x > __t - __np;
1653 const double __lfx =
1654 std::lgamma(__np + __x + 1)
1655 + std::lgamma(__t - (__np + __x) + 1);
1656 __reject = __v > __param._M_lf - __lfx
1657 + __x * __param._M_lp1p;
1660 __reject |= __x + __np >= __thr;
1664 __x += __np + __naf;
1666 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1668 __ret = _IntType(__x) + __z;
1672 __ret = _M_waiting(__urng, __t, __param._M_q);
1675 __ret = __t - __ret;
1679 template<typename _IntType>
1680 template<typename _ForwardIterator,
1681 typename _UniformRandomNumberGenerator>
1683 binomial_distribution<_IntType>::
1684 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1685 _UniformRandomNumberGenerator& __urng,
1686 const param_type& __param)
1688 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1689 // We could duplicate everything from operator()...
1691 *__f++ = this->operator()(__urng, __param);
1694 template<typename _IntType,
1695 typename _CharT, typename _Traits>
1696 std::basic_ostream<_CharT, _Traits>&
1697 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1698 const binomial_distribution<_IntType>& __x)
1700 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1702 const typename __ios_base::fmtflags __flags = __os.flags();
1703 const _CharT __fill = __os.fill();
1704 const std::streamsize __precision = __os.precision();
1705 const _CharT __space = __os.widen(' ');
1706 __os.flags(__ios_base::scientific | __ios_base::left);
1708 __os.precision(std::numeric_limits<double>::max_digits10);
1710 __os << __x.t() << __space << __x.p()
1711 << __space << __x._M_nd;
1713 __os.flags(__flags);
1715 __os.precision(__precision);
1719 template<typename _IntType,
1720 typename _CharT, typename _Traits>
1721 std::basic_istream<_CharT, _Traits>&
1722 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1723 binomial_distribution<_IntType>& __x)
1725 using param_type = typename binomial_distribution<_IntType>::param_type;
1726 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1728 const typename __ios_base::fmtflags __flags = __is.flags();
1729 __is.flags(__ios_base::dec | __ios_base::skipws);
1733 if (__is >> __t >> __p >> __x._M_nd)
1734 __x.param(param_type(__t, __p));
1736 __is.flags(__flags);
1741 template<typename _RealType>
1742 template<typename _ForwardIterator,
1743 typename _UniformRandomNumberGenerator>
1745 std::exponential_distribution<_RealType>::
1746 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1747 _UniformRandomNumberGenerator& __urng,
1748 const param_type& __p)
1750 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1751 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1754 *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1757 template<typename _RealType, typename _CharT, typename _Traits>
1758 std::basic_ostream<_CharT, _Traits>&
1759 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1760 const exponential_distribution<_RealType>& __x)
1762 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1764 const typename __ios_base::fmtflags __flags = __os.flags();
1765 const _CharT __fill = __os.fill();
1766 const std::streamsize __precision = __os.precision();
1767 __os.flags(__ios_base::scientific | __ios_base::left);
1768 __os.fill(__os.widen(' '));
1769 __os.precision(std::numeric_limits<_RealType>::max_digits10);
1771 __os << __x.lambda();
1773 __os.flags(__flags);
1775 __os.precision(__precision);
1779 template<typename _RealType, typename _CharT, typename _Traits>
1780 std::basic_istream<_CharT, _Traits>&
1781 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1782 exponential_distribution<_RealType>& __x)
1785 = typename exponential_distribution<_RealType>::param_type;
1786 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1788 const typename __ios_base::fmtflags __flags = __is.flags();
1789 __is.flags(__ios_base::dec | __ios_base::skipws);
1792 if (__is >> __lambda)
1793 __x.param(param_type(__lambda));
1795 __is.flags(__flags);
1801 * Polar method due to Marsaglia.
1803 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1804 * New York, 1986, Ch. V, Sect. 4.4.
1806 template<typename _RealType>
1807 template<typename _UniformRandomNumberGenerator>
1808 typename normal_distribution<_RealType>::result_type
1809 normal_distribution<_RealType>::
1810 operator()(_UniformRandomNumberGenerator& __urng,
1811 const param_type& __param)
1814 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1817 if (_M_saved_available)
1819 _M_saved_available = false;
1824 result_type __x, __y, __r2;
1827 __x = result_type(2.0) * __aurng() - 1.0;
1828 __y = result_type(2.0) * __aurng() - 1.0;
1829 __r2 = __x * __x + __y * __y;
1831 while (__r2 > 1.0 || __r2 == 0.0);
1833 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1834 _M_saved = __x * __mult;
1835 _M_saved_available = true;
1836 __ret = __y * __mult;
1839 __ret = __ret * __param.stddev() + __param.mean();
1843 template<typename _RealType>
1844 template<typename _ForwardIterator,
1845 typename _UniformRandomNumberGenerator>
1847 normal_distribution<_RealType>::
1848 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1849 _UniformRandomNumberGenerator& __urng,
1850 const param_type& __param)
1852 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1857 if (_M_saved_available)
1859 _M_saved_available = false;
1860 *__f++ = _M_saved * __param.stddev() + __param.mean();
1866 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1869 while (__f + 1 < __t)
1871 result_type __x, __y, __r2;
1874 __x = result_type(2.0) * __aurng() - 1.0;
1875 __y = result_type(2.0) * __aurng() - 1.0;
1876 __r2 = __x * __x + __y * __y;
1878 while (__r2 > 1.0 || __r2 == 0.0);
1880 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1881 *__f++ = __y * __mult * __param.stddev() + __param.mean();
1882 *__f++ = __x * __mult * __param.stddev() + __param.mean();
1887 result_type __x, __y, __r2;
1890 __x = result_type(2.0) * __aurng() - 1.0;
1891 __y = result_type(2.0) * __aurng() - 1.0;
1892 __r2 = __x * __x + __y * __y;
1894 while (__r2 > 1.0 || __r2 == 0.0);
1896 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1897 _M_saved = __x * __mult;
1898 _M_saved_available = true;
1899 *__f = __y * __mult * __param.stddev() + __param.mean();
1903 template<typename _RealType>
1905 operator==(const std::normal_distribution<_RealType>& __d1,
1906 const std::normal_distribution<_RealType>& __d2)
1908 if (__d1._M_param == __d2._M_param
1909 && __d1._M_saved_available == __d2._M_saved_available)
1911 if (__d1._M_saved_available
1912 && __d1._M_saved == __d2._M_saved)
1914 else if(!__d1._M_saved_available)
1923 template<typename _RealType, typename _CharT, typename _Traits>
1924 std::basic_ostream<_CharT, _Traits>&
1925 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1926 const normal_distribution<_RealType>& __x)
1928 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1930 const typename __ios_base::fmtflags __flags = __os.flags();
1931 const _CharT __fill = __os.fill();
1932 const std::streamsize __precision = __os.precision();
1933 const _CharT __space = __os.widen(' ');
1934 __os.flags(__ios_base::scientific | __ios_base::left);
1936 __os.precision(std::numeric_limits<_RealType>::max_digits10);
1938 __os << __x.mean() << __space << __x.stddev()
1939 << __space << __x._M_saved_available;
1940 if (__x._M_saved_available)
1941 __os << __space << __x._M_saved;
1943 __os.flags(__flags);
1945 __os.precision(__precision);
1949 template<typename _RealType, typename _CharT, typename _Traits>
1950 std::basic_istream<_CharT, _Traits>&
1951 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1952 normal_distribution<_RealType>& __x)
1954 using param_type = typename normal_distribution<_RealType>::param_type;
1955 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1957 const typename __ios_base::fmtflags __flags = __is.flags();
1958 __is.flags(__ios_base::dec | __ios_base::skipws);
1960 double __mean, __stddev;
1962 if (__is >> __mean >> __stddev >> __saved_avail)
1964 if (!__saved_avail || (__is >> __x._M_saved))
1966 __x._M_saved_available = __saved_avail;
1967 __x.param(param_type(__mean, __stddev));
1971 __is.flags(__flags);
1976 template<typename _RealType>
1977 template<typename _ForwardIterator,
1978 typename _UniformRandomNumberGenerator>
1980 lognormal_distribution<_RealType>::
1981 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1982 _UniformRandomNumberGenerator& __urng,
1983 const param_type& __p)
1985 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1987 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
1990 template<typename _RealType, typename _CharT, typename _Traits>
1991 std::basic_ostream<_CharT, _Traits>&
1992 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1993 const lognormal_distribution<_RealType>& __x)
1995 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1997 const typename __ios_base::fmtflags __flags = __os.flags();
1998 const _CharT __fill = __os.fill();
1999 const std::streamsize __precision = __os.precision();
2000 const _CharT __space = __os.widen(' ');
2001 __os.flags(__ios_base::scientific | __ios_base::left);
2003 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2005 __os << __x.m() << __space << __x.s()
2006 << __space << __x._M_nd;
2008 __os.flags(__flags);
2010 __os.precision(__precision);
2014 template<typename _RealType, typename _CharT, typename _Traits>
2015 std::basic_istream<_CharT, _Traits>&
2016 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2017 lognormal_distribution<_RealType>& __x)
2020 = typename lognormal_distribution<_RealType>::param_type;
2021 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2023 const typename __ios_base::fmtflags __flags = __is.flags();
2024 __is.flags(__ios_base::dec | __ios_base::skipws);
2027 if (__is >> __m >> __s >> __x._M_nd)
2028 __x.param(param_type(__m, __s));
2030 __is.flags(__flags);
2034 template<typename _RealType>
2035 template<typename _ForwardIterator,
2036 typename _UniformRandomNumberGenerator>
2038 std::chi_squared_distribution<_RealType>::
2039 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2040 _UniformRandomNumberGenerator& __urng)
2042 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2044 *__f++ = 2 * _M_gd(__urng);
2047 template<typename _RealType>
2048 template<typename _ForwardIterator,
2049 typename _UniformRandomNumberGenerator>
2051 std::chi_squared_distribution<_RealType>::
2052 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2053 _UniformRandomNumberGenerator& __urng,
2055 std::gamma_distribution<result_type>::param_type& __p)
2057 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2059 *__f++ = 2 * _M_gd(__urng, __p);
2062 template<typename _RealType, typename _CharT, typename _Traits>
2063 std::basic_ostream<_CharT, _Traits>&
2064 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2065 const chi_squared_distribution<_RealType>& __x)
2067 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2069 const typename __ios_base::fmtflags __flags = __os.flags();
2070 const _CharT __fill = __os.fill();
2071 const std::streamsize __precision = __os.precision();
2072 const _CharT __space = __os.widen(' ');
2073 __os.flags(__ios_base::scientific | __ios_base::left);
2075 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2077 __os << __x.n() << __space << __x._M_gd;
2079 __os.flags(__flags);
2081 __os.precision(__precision);
2085 template<typename _RealType, typename _CharT, typename _Traits>
2086 std::basic_istream<_CharT, _Traits>&
2087 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2088 chi_squared_distribution<_RealType>& __x)
2091 = typename chi_squared_distribution<_RealType>::param_type;
2092 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2094 const typename __ios_base::fmtflags __flags = __is.flags();
2095 __is.flags(__ios_base::dec | __ios_base::skipws);
2098 if (__is >> __n >> __x._M_gd)
2099 __x.param(param_type(__n));
2101 __is.flags(__flags);
2106 template<typename _RealType>
2107 template<typename _UniformRandomNumberGenerator>
2108 typename cauchy_distribution<_RealType>::result_type
2109 cauchy_distribution<_RealType>::
2110 operator()(_UniformRandomNumberGenerator& __urng,
2111 const param_type& __p)
2113 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2120 const _RealType __pi = 3.1415926535897932384626433832795029L;
2121 return __p.a() + __p.b() * std::tan(__pi * __u);
2124 template<typename _RealType>
2125 template<typename _ForwardIterator,
2126 typename _UniformRandomNumberGenerator>
2128 cauchy_distribution<_RealType>::
2129 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2130 _UniformRandomNumberGenerator& __urng,
2131 const param_type& __p)
2133 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2134 const _RealType __pi = 3.1415926535897932384626433832795029L;
2135 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2144 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2148 template<typename _RealType, typename _CharT, typename _Traits>
2149 std::basic_ostream<_CharT, _Traits>&
2150 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2151 const cauchy_distribution<_RealType>& __x)
2153 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2155 const typename __ios_base::fmtflags __flags = __os.flags();
2156 const _CharT __fill = __os.fill();
2157 const std::streamsize __precision = __os.precision();
2158 const _CharT __space = __os.widen(' ');
2159 __os.flags(__ios_base::scientific | __ios_base::left);
2161 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2163 __os << __x.a() << __space << __x.b();
2165 __os.flags(__flags);
2167 __os.precision(__precision);
2171 template<typename _RealType, typename _CharT, typename _Traits>
2172 std::basic_istream<_CharT, _Traits>&
2173 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2174 cauchy_distribution<_RealType>& __x)
2176 using param_type = typename cauchy_distribution<_RealType>::param_type;
2177 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2179 const typename __ios_base::fmtflags __flags = __is.flags();
2180 __is.flags(__ios_base::dec | __ios_base::skipws);
2183 if (__is >> __a >> __b)
2184 __x.param(param_type(__a, __b));
2186 __is.flags(__flags);
2191 template<typename _RealType>
2192 template<typename _ForwardIterator,
2193 typename _UniformRandomNumberGenerator>
2195 std::fisher_f_distribution<_RealType>::
2196 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2197 _UniformRandomNumberGenerator& __urng)
2199 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2201 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2204 template<typename _RealType>
2205 template<typename _ForwardIterator,
2206 typename _UniformRandomNumberGenerator>
2208 std::fisher_f_distribution<_RealType>::
2209 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2210 _UniformRandomNumberGenerator& __urng,
2211 const param_type& __p)
2213 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2214 typedef typename std::gamma_distribution<result_type>::param_type
2216 param_type __p1(__p.m() / 2);
2217 param_type __p2(__p.n() / 2);
2219 *__f++ = ((_M_gd_x(__urng, __p1) * n())
2220 / (_M_gd_y(__urng, __p2) * m()));
2223 template<typename _RealType, typename _CharT, typename _Traits>
2224 std::basic_ostream<_CharT, _Traits>&
2225 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2226 const fisher_f_distribution<_RealType>& __x)
2228 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2230 const typename __ios_base::fmtflags __flags = __os.flags();
2231 const _CharT __fill = __os.fill();
2232 const std::streamsize __precision = __os.precision();
2233 const _CharT __space = __os.widen(' ');
2234 __os.flags(__ios_base::scientific | __ios_base::left);
2236 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2238 __os << __x.m() << __space << __x.n()
2239 << __space << __x._M_gd_x << __space << __x._M_gd_y;
2241 __os.flags(__flags);
2243 __os.precision(__precision);
2247 template<typename _RealType, typename _CharT, typename _Traits>
2248 std::basic_istream<_CharT, _Traits>&
2249 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2250 fisher_f_distribution<_RealType>& __x)
2253 = typename fisher_f_distribution<_RealType>::param_type;
2254 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2256 const typename __ios_base::fmtflags __flags = __is.flags();
2257 __is.flags(__ios_base::dec | __ios_base::skipws);
2260 if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y)
2261 __x.param(param_type(__m, __n));
2263 __is.flags(__flags);
2268 template<typename _RealType>
2269 template<typename _ForwardIterator,
2270 typename _UniformRandomNumberGenerator>
2272 std::student_t_distribution<_RealType>::
2273 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2274 _UniformRandomNumberGenerator& __urng)
2276 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2278 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2281 template<typename _RealType>
2282 template<typename _ForwardIterator,
2283 typename _UniformRandomNumberGenerator>
2285 std::student_t_distribution<_RealType>::
2286 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2287 _UniformRandomNumberGenerator& __urng,
2288 const param_type& __p)
2290 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2291 typename std::gamma_distribution<result_type>::param_type
2292 __p2(__p.n() / 2, 2);
2294 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2297 template<typename _RealType, typename _CharT, typename _Traits>
2298 std::basic_ostream<_CharT, _Traits>&
2299 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2300 const student_t_distribution<_RealType>& __x)
2302 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2304 const typename __ios_base::fmtflags __flags = __os.flags();
2305 const _CharT __fill = __os.fill();
2306 const std::streamsize __precision = __os.precision();
2307 const _CharT __space = __os.widen(' ');
2308 __os.flags(__ios_base::scientific | __ios_base::left);
2310 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2312 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2314 __os.flags(__flags);
2316 __os.precision(__precision);
2320 template<typename _RealType, typename _CharT, typename _Traits>
2321 std::basic_istream<_CharT, _Traits>&
2322 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2323 student_t_distribution<_RealType>& __x)
2326 = typename student_t_distribution<_RealType>::param_type;
2327 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2329 const typename __ios_base::fmtflags __flags = __is.flags();
2330 __is.flags(__ios_base::dec | __ios_base::skipws);
2333 if (__is >> __n >> __x._M_nd >> __x._M_gd)
2334 __x.param(param_type(__n));
2336 __is.flags(__flags);
2341 template<typename _RealType>
2343 gamma_distribution<_RealType>::param_type::
2346 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2348 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2349 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2353 * Marsaglia, G. and Tsang, W. W.
2354 * "A Simple Method for Generating Gamma Variables"
2355 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2357 template<typename _RealType>
2358 template<typename _UniformRandomNumberGenerator>
2359 typename gamma_distribution<_RealType>::result_type
2360 gamma_distribution<_RealType>::
2361 operator()(_UniformRandomNumberGenerator& __urng,
2362 const param_type& __param)
2364 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2367 result_type __u, __v, __n;
2368 const result_type __a1 = (__param._M_malpha
2369 - _RealType(1.0) / _RealType(3.0));
2375 __n = _M_nd(__urng);
2376 __v = result_type(1.0) + __param._M_a2 * __n;
2380 __v = __v * __v * __v;
2383 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2384 && (std::log(__u) > (0.5 * __n * __n + __a1
2385 * (1.0 - __v + std::log(__v)))));
2387 if (__param.alpha() == __param._M_malpha)
2388 return __a1 * __v * __param.beta();
2395 return (std::pow(__u, result_type(1.0) / __param.alpha())
2396 * __a1 * __v * __param.beta());
2400 template<typename _RealType>
2401 template<typename _ForwardIterator,
2402 typename _UniformRandomNumberGenerator>
2404 gamma_distribution<_RealType>::
2405 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2406 _UniformRandomNumberGenerator& __urng,
2407 const param_type& __param)
2409 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2410 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2413 result_type __u, __v, __n;
2414 const result_type __a1 = (__param._M_malpha
2415 - _RealType(1.0) / _RealType(3.0));
2417 if (__param.alpha() == __param._M_malpha)
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)))));
2436 *__f++ = __a1 * __v * __param.beta();
2445 __n = _M_nd(__urng);
2446 __v = result_type(1.0) + __param._M_a2 * __n;
2450 __v = __v * __v * __v;
2453 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2454 && (std::log(__u) > (0.5 * __n * __n + __a1
2455 * (1.0 - __v + std::log(__v)))));
2461 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2462 * __a1 * __v * __param.beta());
2466 template<typename _RealType, typename _CharT, typename _Traits>
2467 std::basic_ostream<_CharT, _Traits>&
2468 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2469 const gamma_distribution<_RealType>& __x)
2471 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2473 const typename __ios_base::fmtflags __flags = __os.flags();
2474 const _CharT __fill = __os.fill();
2475 const std::streamsize __precision = __os.precision();
2476 const _CharT __space = __os.widen(' ');
2477 __os.flags(__ios_base::scientific | __ios_base::left);
2479 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2481 __os << __x.alpha() << __space << __x.beta()
2482 << __space << __x._M_nd;
2484 __os.flags(__flags);
2486 __os.precision(__precision);
2490 template<typename _RealType, typename _CharT, typename _Traits>
2491 std::basic_istream<_CharT, _Traits>&
2492 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2493 gamma_distribution<_RealType>& __x)
2495 using param_type = typename gamma_distribution<_RealType>::param_type;
2496 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2498 const typename __ios_base::fmtflags __flags = __is.flags();
2499 __is.flags(__ios_base::dec | __ios_base::skipws);
2501 _RealType __alpha_val, __beta_val;
2502 if (__is >> __alpha_val >> __beta_val >> __x._M_nd)
2503 __x.param(param_type(__alpha_val, __beta_val));
2505 __is.flags(__flags);
2510 template<typename _RealType>
2511 template<typename _UniformRandomNumberGenerator>
2512 typename weibull_distribution<_RealType>::result_type
2513 weibull_distribution<_RealType>::
2514 operator()(_UniformRandomNumberGenerator& __urng,
2515 const param_type& __p)
2517 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2519 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2520 result_type(1) / __p.a());
2523 template<typename _RealType>
2524 template<typename _ForwardIterator,
2525 typename _UniformRandomNumberGenerator>
2527 weibull_distribution<_RealType>::
2528 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2529 _UniformRandomNumberGenerator& __urng,
2530 const param_type& __p)
2532 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2533 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2535 auto __inv_a = result_type(1) / __p.a();
2538 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2542 template<typename _RealType, typename _CharT, typename _Traits>
2543 std::basic_ostream<_CharT, _Traits>&
2544 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2545 const weibull_distribution<_RealType>& __x)
2547 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2549 const typename __ios_base::fmtflags __flags = __os.flags();
2550 const _CharT __fill = __os.fill();
2551 const std::streamsize __precision = __os.precision();
2552 const _CharT __space = __os.widen(' ');
2553 __os.flags(__ios_base::scientific | __ios_base::left);
2555 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2557 __os << __x.a() << __space << __x.b();
2559 __os.flags(__flags);
2561 __os.precision(__precision);
2565 template<typename _RealType, typename _CharT, typename _Traits>
2566 std::basic_istream<_CharT, _Traits>&
2567 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2568 weibull_distribution<_RealType>& __x)
2570 using param_type = typename weibull_distribution<_RealType>::param_type;
2571 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2573 const typename __ios_base::fmtflags __flags = __is.flags();
2574 __is.flags(__ios_base::dec | __ios_base::skipws);
2577 if (__is >> __a >> __b)
2578 __x.param(param_type(__a, __b));
2580 __is.flags(__flags);
2585 template<typename _RealType>
2586 template<typename _UniformRandomNumberGenerator>
2587 typename extreme_value_distribution<_RealType>::result_type
2588 extreme_value_distribution<_RealType>::
2589 operator()(_UniformRandomNumberGenerator& __urng,
2590 const param_type& __p)
2592 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2594 return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2598 template<typename _RealType>
2599 template<typename _ForwardIterator,
2600 typename _UniformRandomNumberGenerator>
2602 extreme_value_distribution<_RealType>::
2603 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2604 _UniformRandomNumberGenerator& __urng,
2605 const param_type& __p)
2607 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2608 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2612 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2616 template<typename _RealType, typename _CharT, typename _Traits>
2617 std::basic_ostream<_CharT, _Traits>&
2618 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2619 const extreme_value_distribution<_RealType>& __x)
2621 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2623 const typename __ios_base::fmtflags __flags = __os.flags();
2624 const _CharT __fill = __os.fill();
2625 const std::streamsize __precision = __os.precision();
2626 const _CharT __space = __os.widen(' ');
2627 __os.flags(__ios_base::scientific | __ios_base::left);
2629 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2631 __os << __x.a() << __space << __x.b();
2633 __os.flags(__flags);
2635 __os.precision(__precision);
2639 template<typename _RealType, typename _CharT, typename _Traits>
2640 std::basic_istream<_CharT, _Traits>&
2641 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2642 extreme_value_distribution<_RealType>& __x)
2645 = typename extreme_value_distribution<_RealType>::param_type;
2646 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2648 const typename __ios_base::fmtflags __flags = __is.flags();
2649 __is.flags(__ios_base::dec | __ios_base::skipws);
2652 if (__is >> __a >> __b)
2653 __x.param(param_type(__a, __b));
2655 __is.flags(__flags);
2660 template<typename _IntType>
2662 discrete_distribution<_IntType>::param_type::
2665 if (_M_prob.size() < 2)
2671 const double __sum = std::accumulate(_M_prob.begin(),
2672 _M_prob.end(), 0.0);
2673 __glibcxx_assert(__sum > 0);
2674 // Now normalize the probabilites.
2675 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2677 // Accumulate partial sums.
2678 _M_cp.reserve(_M_prob.size());
2679 std::partial_sum(_M_prob.begin(), _M_prob.end(),
2680 std::back_inserter(_M_cp));
2681 // Make sure the last cumulative probability is one.
2682 _M_cp[_M_cp.size() - 1] = 1.0;
2685 template<typename _IntType>
2686 template<typename _Func>
2687 discrete_distribution<_IntType>::param_type::
2688 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2689 : _M_prob(), _M_cp()
2691 const size_t __n = __nw == 0 ? 1 : __nw;
2692 const double __delta = (__xmax - __xmin) / __n;
2694 _M_prob.reserve(__n);
2695 for (size_t __k = 0; __k < __nw; ++__k)
2696 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2701 template<typename _IntType>
2702 template<typename _UniformRandomNumberGenerator>
2703 typename discrete_distribution<_IntType>::result_type
2704 discrete_distribution<_IntType>::
2705 operator()(_UniformRandomNumberGenerator& __urng,
2706 const param_type& __param)
2708 if (__param._M_cp.empty())
2709 return result_type(0);
2711 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2714 const double __p = __aurng();
2715 auto __pos = std::lower_bound(__param._M_cp.begin(),
2716 __param._M_cp.end(), __p);
2718 return __pos - __param._M_cp.begin();
2721 template<typename _IntType>
2722 template<typename _ForwardIterator,
2723 typename _UniformRandomNumberGenerator>
2725 discrete_distribution<_IntType>::
2726 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2727 _UniformRandomNumberGenerator& __urng,
2728 const param_type& __param)
2730 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2732 if (__param._M_cp.empty())
2735 *__f++ = result_type(0);
2739 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2744 const double __p = __aurng();
2745 auto __pos = std::lower_bound(__param._M_cp.begin(),
2746 __param._M_cp.end(), __p);
2748 *__f++ = __pos - __param._M_cp.begin();
2752 template<typename _IntType, typename _CharT, typename _Traits>
2753 std::basic_ostream<_CharT, _Traits>&
2754 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2755 const discrete_distribution<_IntType>& __x)
2757 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2759 const typename __ios_base::fmtflags __flags = __os.flags();
2760 const _CharT __fill = __os.fill();
2761 const std::streamsize __precision = __os.precision();
2762 const _CharT __space = __os.widen(' ');
2763 __os.flags(__ios_base::scientific | __ios_base::left);
2765 __os.precision(std::numeric_limits<double>::max_digits10);
2767 std::vector<double> __prob = __x.probabilities();
2768 __os << __prob.size();
2769 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2770 __os << __space << *__dit;
2772 __os.flags(__flags);
2774 __os.precision(__precision);
2780 template<typename _ValT, typename _CharT, typename _Traits>
2781 basic_istream<_CharT, _Traits>&
2782 __extract_params(basic_istream<_CharT, _Traits>& __is,
2783 vector<_ValT>& __vals, size_t __n)
2785 __vals.reserve(__n);
2790 __vals.push_back(__val);
2796 } // namespace __detail
2798 template<typename _IntType, typename _CharT, typename _Traits>
2799 std::basic_istream<_CharT, _Traits>&
2800 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2801 discrete_distribution<_IntType>& __x)
2803 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2805 const typename __ios_base::fmtflags __flags = __is.flags();
2806 __is.flags(__ios_base::dec | __ios_base::skipws);
2811 std::vector<double> __prob_vec;
2812 if (__detail::__extract_params(__is, __prob_vec, __n))
2813 __x.param({__prob_vec.begin(), __prob_vec.end()});
2816 __is.flags(__flags);
2821 template<typename _RealType>
2823 piecewise_constant_distribution<_RealType>::param_type::
2826 if (_M_int.size() < 2
2827 || (_M_int.size() == 2
2828 && _M_int[0] == _RealType(0)
2829 && _M_int[1] == _RealType(1)))
2836 const double __sum = std::accumulate(_M_den.begin(),
2838 __glibcxx_assert(__sum > 0);
2840 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2843 _M_cp.reserve(_M_den.size());
2844 std::partial_sum(_M_den.begin(), _M_den.end(),
2845 std::back_inserter(_M_cp));
2847 // Make sure the last cumulative probability is one.
2848 _M_cp[_M_cp.size() - 1] = 1.0;
2850 for (size_t __k = 0; __k < _M_den.size(); ++__k)
2851 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2854 template<typename _RealType>
2855 template<typename _InputIteratorB, typename _InputIteratorW>
2856 piecewise_constant_distribution<_RealType>::param_type::
2857 param_type(_InputIteratorB __bbegin,
2858 _InputIteratorB __bend,
2859 _InputIteratorW __wbegin)
2860 : _M_int(), _M_den(), _M_cp()
2862 if (__bbegin != __bend)
2866 _M_int.push_back(*__bbegin);
2868 if (__bbegin == __bend)
2871 _M_den.push_back(*__wbegin);
2879 template<typename _RealType>
2880 template<typename _Func>
2881 piecewise_constant_distribution<_RealType>::param_type::
2882 param_type(initializer_list<_RealType> __bl, _Func __fw)
2883 : _M_int(), _M_den(), _M_cp()
2885 _M_int.reserve(__bl.size());
2886 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2887 _M_int.push_back(*__biter);
2889 _M_den.reserve(_M_int.size() - 1);
2890 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2891 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2896 template<typename _RealType>
2897 template<typename _Func>
2898 piecewise_constant_distribution<_RealType>::param_type::
2899 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2900 : _M_int(), _M_den(), _M_cp()
2902 const size_t __n = __nw == 0 ? 1 : __nw;
2903 const _RealType __delta = (__xmax - __xmin) / __n;
2905 _M_int.reserve(__n + 1);
2906 for (size_t __k = 0; __k <= __nw; ++__k)
2907 _M_int.push_back(__xmin + __k * __delta);
2909 _M_den.reserve(__n);
2910 for (size_t __k = 0; __k < __nw; ++__k)
2911 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2916 template<typename _RealType>
2917 template<typename _UniformRandomNumberGenerator>
2918 typename piecewise_constant_distribution<_RealType>::result_type
2919 piecewise_constant_distribution<_RealType>::
2920 operator()(_UniformRandomNumberGenerator& __urng,
2921 const param_type& __param)
2923 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2926 const double __p = __aurng();
2927 if (__param._M_cp.empty())
2930 auto __pos = std::lower_bound(__param._M_cp.begin(),
2931 __param._M_cp.end(), __p);
2932 const size_t __i = __pos - __param._M_cp.begin();
2934 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2936 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2939 template<typename _RealType>
2940 template<typename _ForwardIterator,
2941 typename _UniformRandomNumberGenerator>
2943 piecewise_constant_distribution<_RealType>::
2944 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2945 _UniformRandomNumberGenerator& __urng,
2946 const param_type& __param)
2948 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2949 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2952 if (__param._M_cp.empty())
2961 const double __p = __aurng();
2963 auto __pos = std::lower_bound(__param._M_cp.begin(),
2964 __param._M_cp.end(), __p);
2965 const size_t __i = __pos - __param._M_cp.begin();
2967 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2969 *__f++ = (__param._M_int[__i]
2970 + (__p - __pref) / __param._M_den[__i]);
2974 template<typename _RealType, typename _CharT, typename _Traits>
2975 std::basic_ostream<_CharT, _Traits>&
2976 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2977 const piecewise_constant_distribution<_RealType>& __x)
2979 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2981 const typename __ios_base::fmtflags __flags = __os.flags();
2982 const _CharT __fill = __os.fill();
2983 const std::streamsize __precision = __os.precision();
2984 const _CharT __space = __os.widen(' ');
2985 __os.flags(__ios_base::scientific | __ios_base::left);
2987 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2989 std::vector<_RealType> __int = __x.intervals();
2990 __os << __int.size() - 1;
2992 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2993 __os << __space << *__xit;
2995 std::vector<double> __den = __x.densities();
2996 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2997 __os << __space << *__dit;
2999 __os.flags(__flags);
3001 __os.precision(__precision);
3005 template<typename _RealType, typename _CharT, typename _Traits>
3006 std::basic_istream<_CharT, _Traits>&
3007 operator>>(std::basic_istream<_CharT, _Traits>& __is,
3008 piecewise_constant_distribution<_RealType>& __x)
3010 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3012 const typename __ios_base::fmtflags __flags = __is.flags();
3013 __is.flags(__ios_base::dec | __ios_base::skipws);
3018 std::vector<_RealType> __int_vec;
3019 if (__detail::__extract_params(__is, __int_vec, __n + 1))
3021 std::vector<double> __den_vec;
3022 if (__detail::__extract_params(__is, __den_vec, __n))
3024 __x.param({ __int_vec.begin(), __int_vec.end(),
3025 __den_vec.begin() });
3030 __is.flags(__flags);
3035 template<typename _RealType>
3037 piecewise_linear_distribution<_RealType>::param_type::
3040 if (_M_int.size() < 2
3041 || (_M_int.size() == 2
3042 && _M_int[0] == _RealType(0)
3043 && _M_int[1] == _RealType(1)
3044 && _M_den[0] == _M_den[1]))
3052 _M_cp.reserve(_M_int.size() - 1);
3053 _M_m.reserve(_M_int.size() - 1);
3054 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3056 const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3057 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3058 _M_cp.push_back(__sum);
3059 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3061 __glibcxx_assert(__sum > 0);
3063 // Now normalize the densities...
3064 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3066 // ... and partial sums...
3067 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3069 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3071 // Make sure the last cumulative probablility is one.
3072 _M_cp[_M_cp.size() - 1] = 1.0;
3075 template<typename _RealType>
3076 template<typename _InputIteratorB, typename _InputIteratorW>
3077 piecewise_linear_distribution<_RealType>::param_type::
3078 param_type(_InputIteratorB __bbegin,
3079 _InputIteratorB __bend,
3080 _InputIteratorW __wbegin)
3081 : _M_int(), _M_den(), _M_cp(), _M_m()
3083 for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
3085 _M_int.push_back(*__bbegin);
3086 _M_den.push_back(*__wbegin);
3092 template<typename _RealType>
3093 template<typename _Func>
3094 piecewise_linear_distribution<_RealType>::param_type::
3095 param_type(initializer_list<_RealType> __bl, _Func __fw)
3096 : _M_int(), _M_den(), _M_cp(), _M_m()
3098 _M_int.reserve(__bl.size());
3099 _M_den.reserve(__bl.size());
3100 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3102 _M_int.push_back(*__biter);
3103 _M_den.push_back(__fw(*__biter));
3109 template<typename _RealType>
3110 template<typename _Func>
3111 piecewise_linear_distribution<_RealType>::param_type::
3112 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3113 : _M_int(), _M_den(), _M_cp(), _M_m()
3115 const size_t __n = __nw == 0 ? 1 : __nw;
3116 const _RealType __delta = (__xmax - __xmin) / __n;
3118 _M_int.reserve(__n + 1);
3119 _M_den.reserve(__n + 1);
3120 for (size_t __k = 0; __k <= __nw; ++__k)
3122 _M_int.push_back(__xmin + __k * __delta);
3123 _M_den.push_back(__fw(_M_int[__k] + __delta));
3129 template<typename _RealType>
3130 template<typename _UniformRandomNumberGenerator>
3131 typename piecewise_linear_distribution<_RealType>::result_type
3132 piecewise_linear_distribution<_RealType>::
3133 operator()(_UniformRandomNumberGenerator& __urng,
3134 const param_type& __param)
3136 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3139 const double __p = __aurng();
3140 if (__param._M_cp.empty())
3143 auto __pos = std::lower_bound(__param._M_cp.begin(),
3144 __param._M_cp.end(), __p);
3145 const size_t __i = __pos - __param._M_cp.begin();
3147 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3149 const double __a = 0.5 * __param._M_m[__i];
3150 const double __b = __param._M_den[__i];
3151 const double __cm = __p - __pref;
3153 _RealType __x = __param._M_int[__i];
3158 const double __d = __b * __b + 4.0 * __a * __cm;
3159 __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3165 template<typename _RealType>
3166 template<typename _ForwardIterator,
3167 typename _UniformRandomNumberGenerator>
3169 piecewise_linear_distribution<_RealType>::
3170 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3171 _UniformRandomNumberGenerator& __urng,
3172 const param_type& __param)
3174 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3175 // We could duplicate everything from operator()...
3177 *__f++ = this->operator()(__urng, __param);
3180 template<typename _RealType, typename _CharT, typename _Traits>
3181 std::basic_ostream<_CharT, _Traits>&
3182 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3183 const piecewise_linear_distribution<_RealType>& __x)
3185 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3187 const typename __ios_base::fmtflags __flags = __os.flags();
3188 const _CharT __fill = __os.fill();
3189 const std::streamsize __precision = __os.precision();
3190 const _CharT __space = __os.widen(' ');
3191 __os.flags(__ios_base::scientific | __ios_base::left);
3193 __os.precision(std::numeric_limits<_RealType>::max_digits10);
3195 std::vector<_RealType> __int = __x.intervals();
3196 __os << __int.size() - 1;
3198 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3199 __os << __space << *__xit;
3201 std::vector<double> __den = __x.densities();
3202 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3203 __os << __space << *__dit;
3205 __os.flags(__flags);
3207 __os.precision(__precision);
3211 template<typename _RealType, typename _CharT, typename _Traits>
3212 std::basic_istream<_CharT, _Traits>&
3213 operator>>(std::basic_istream<_CharT, _Traits>& __is,
3214 piecewise_linear_distribution<_RealType>& __x)
3216 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3218 const typename __ios_base::fmtflags __flags = __is.flags();
3219 __is.flags(__ios_base::dec | __ios_base::skipws);
3224 vector<_RealType> __int_vec;
3225 if (__detail::__extract_params(__is, __int_vec, __n + 1))
3227 vector<double> __den_vec;
3228 if (__detail::__extract_params(__is, __den_vec, __n + 1))
3230 __x.param({ __int_vec.begin(), __int_vec.end(),
3231 __den_vec.begin() });
3235 __is.flags(__flags);
3240 template<typename _IntType, typename>
3241 seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3243 _M_v.reserve(__il.size());
3244 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3245 _M_v.push_back(__detail::__mod<result_type,
3246 __detail::_Shift<result_type, 32>::__value>(*__iter));
3249 template<typename _InputIterator>
3250 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3252 if _GLIBCXX17_CONSTEXPR (__is_random_access_iter<_InputIterator>::value)
3253 _M_v.reserve(std::distance(__begin, __end));
3255 for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3256 _M_v.push_back(__detail::__mod<result_type,
3257 __detail::_Shift<result_type, 32>::__value>(*__iter));
3260 template<typename _RandomAccessIterator>
3262 seed_seq::generate(_RandomAccessIterator __begin,
3263 _RandomAccessIterator __end)
3265 typedef typename iterator_traits<_RandomAccessIterator>::value_type
3268 if (__begin == __end)
3271 std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3273 const size_t __n = __end - __begin;
3274 const size_t __s = _M_v.size();
3275 const size_t __t = (__n >= 623) ? 11
3280 const size_t __p = (__n - __t) / 2;
3281 const size_t __q = __p + __t;
3282 const size_t __m = std::max(size_t(__s + 1), __n);
3284 #ifndef __UINT32_TYPE__
3287 _Up(uint_least32_t v) : _M_v(v & 0xffffffffu) { }
3289 operator uint_least32_t() const { return _M_v; }
3291 uint_least32_t _M_v;
3293 using uint32_t = _Up;
3296 // k == 0, every element in [begin,end) equals 0x8b8b8b8bu
3298 uint32_t __r1 = 1371501266u;
3299 uint32_t __r2 = __r1 + __s;
3300 __begin[__p] += __r1;
3301 __begin[__q] = (uint32_t)__begin[__q] + __r2;
3305 for (size_t __k = 1; __k <= __s; ++__k)
3307 const size_t __kn = __k % __n;
3308 const size_t __kpn = (__k + __p) % __n;
3309 const size_t __kqn = (__k + __q) % __n;
3310 uint32_t __arg = (__begin[__kn]
3312 ^ __begin[(__k - 1) % __n]);
3313 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3314 uint32_t __r2 = __r1 + (uint32_t)__kn + _M_v[__k - 1];
3315 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3316 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3317 __begin[__kn] = __r2;
3320 for (size_t __k = __s + 1; __k < __m; ++__k)
3322 const size_t __kn = __k % __n;
3323 const size_t __kpn = (__k + __p) % __n;
3324 const size_t __kqn = (__k + __q) % __n;
3325 uint32_t __arg = (__begin[__kn]
3327 ^ __begin[(__k - 1) % __n]);
3328 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3329 uint32_t __r2 = __r1 + (uint32_t)__kn;
3330 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3331 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3332 __begin[__kn] = __r2;
3335 for (size_t __k = __m; __k < __m + __n; ++__k)
3337 const size_t __kn = __k % __n;
3338 const size_t __kpn = (__k + __p) % __n;
3339 const size_t __kqn = (__k + __q) % __n;
3340 uint32_t __arg = (__begin[__kn]
3342 + __begin[(__k - 1) % __n]);
3343 uint32_t __r3 = 1566083941u * (__arg ^ (__arg >> 27));
3344 uint32_t __r4 = __r3 - __kn;
3345 __begin[__kpn] ^= __r3;
3346 __begin[__kqn] ^= __r4;
3347 __begin[__kn] = __r4;
3351 template<typename _RealType, size_t __bits,
3352 typename _UniformRandomNumberGenerator>
3354 generate_canonical(_UniformRandomNumberGenerator& __urng)
3356 static_assert(std::is_floating_point<_RealType>::value,
3357 "template argument must be a floating point type");
3360 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3362 const long double __r = static_cast<long double>(__urng.max())
3363 - static_cast<long double>(__urng.min()) + 1.0L;
3364 const size_t __log2r = std::log(__r) / std::log(2.0L);
3365 const size_t __m = std::max<size_t>(1UL,
3366 (__b + __log2r - 1UL) / __log2r);
3368 _RealType __sum = _RealType(0);
3369 _RealType __tmp = _RealType(1);
3370 for (size_t __k = __m; __k != 0; --__k)
3372 __sum += _RealType(__urng() - __urng.min()) * __tmp;
3375 __ret = __sum / __tmp;
3376 if (__builtin_expect(__ret >= _RealType(1), 0))
3378 #if _GLIBCXX_USE_C99_MATH_TR1
3379 __ret = std::nextafter(_RealType(1), _RealType(0));
3381 __ret = _RealType(1)
3382 - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3388 _GLIBCXX_END_NAMESPACE_VERSION