libstdc++
bits/random.tcc
Go to the documentation of this file.
1// random number generation (out of line) -*- C++ -*-
2
3// Copyright (C) 2009-2020 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library. This library is free
6// software; you can redistribute it and/or modify it under the
7// terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 3, or (at your option)
9// any later version.
10
11// This library is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14// GNU General Public License for more details.
15
16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.
19
20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23// <http://www.gnu.org/licenses/>.
24
25/** @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}
28 */
29
30#ifndef _RANDOM_TCC
31#define _RANDOM_TCC 1
32
33#include <numeric> // std::accumulate and std::partial_sum
34
35namespace std _GLIBCXX_VISIBILITY(default)
36{
37_GLIBCXX_BEGIN_NAMESPACE_VERSION
38
39 /*
40 * (Further) implementation-space details.
41 */
42 namespace __detail
43 {
44 // General case for x = (ax + c) mod m -- use Schrage's algorithm
45 // to avoid integer overflow.
46 //
47 // Preconditions: a > 0, m > 0.
48 //
49 // Note: only works correctly for __m % __a < __m / __a.
50 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
51 _Tp
52 _Mod<_Tp, __m, __a, __c, false, true>::
53 __calc(_Tp __x)
54 {
55 if (__a == 1)
56 __x %= __m;
57 else
58 {
59 static const _Tp __q = __m / __a;
60 static const _Tp __r = __m % __a;
61
62 _Tp __t1 = __a * (__x % __q);
63 _Tp __t2 = __r * (__x / __q);
64 if (__t1 >= __t2)
65 __x = __t1 - __t2;
66 else
67 __x = __m - __t2 + __t1;
68 }
69
70 if (__c != 0)
71 {
72 const _Tp __d = __m - __x;
73 if (__d > __c)
74 __x += __c;
75 else
76 __x = __c - __d;
77 }
78 return __x;
79 }
80
81 template<typename _InputIterator, typename _OutputIterator,
82 typename _Tp>
83 _OutputIterator
84 __normalize(_InputIterator __first, _InputIterator __last,
85 _OutputIterator __result, const _Tp& __factor)
86 {
87 for (; __first != __last; ++__first, ++__result)
88 *__result = *__first / __factor;
89 return __result;
90 }
91
92 } // namespace __detail
93
94 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
95 constexpr _UIntType
97
98 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
99 constexpr _UIntType
101
102 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
103 constexpr _UIntType
105
106 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
107 constexpr _UIntType
108 linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
109
110 /**
111 * Seeds the LCR with integral value @p __s, adjusted so that the
112 * ring identity is never a member of the convergence set.
113 */
114 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
115 void
118 {
119 if ((__detail::__mod<_UIntType, __m>(__c) == 0)
120 && (__detail::__mod<_UIntType, __m>(__s) == 0))
121 _M_x = 1;
122 else
123 _M_x = __detail::__mod<_UIntType, __m>(__s);
124 }
125
126 /**
127 * Seeds the LCR engine with a value generated by @p __q.
128 */
129 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
130 template<typename _Sseq>
131 auto
133 seed(_Sseq& __q)
134 -> _If_seed_seq<_Sseq>
135 {
136 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
137 : std::__lg(__m);
138 const _UIntType __k = (__k0 + 31) / 32;
139 uint_least32_t __arr[__k + 3];
140 __q.generate(__arr + 0, __arr + __k + 3);
141 _UIntType __factor = 1u;
142 _UIntType __sum = 0u;
143 for (size_t __j = 0; __j < __k; ++__j)
144 {
145 __sum += __arr[__j + 3] * __factor;
146 __factor *= __detail::_Shift<_UIntType, 32>::__value;
147 }
148 seed(__sum);
149 }
150
151 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
152 typename _CharT, typename _Traits>
155 const linear_congruential_engine<_UIntType,
156 __a, __c, __m>& __lcr)
157 {
158 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
159
160 const typename __ios_base::fmtflags __flags = __os.flags();
161 const _CharT __fill = __os.fill();
163 __os.fill(__os.widen(' '));
164
165 __os << __lcr._M_x;
166
167 __os.flags(__flags);
168 __os.fill(__fill);
169 return __os;
170 }
171
172 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
173 typename _CharT, typename _Traits>
176 linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
177 {
178 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
179
180 const typename __ios_base::fmtflags __flags = __is.flags();
182
183 __is >> __lcr._M_x;
184
185 __is.flags(__flags);
186 return __is;
187 }
188
189
190 template<typename _UIntType,
191 size_t __w, size_t __n, size_t __m, size_t __r,
192 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
193 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
194 _UIntType __f>
195 constexpr size_t
196 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
197 __s, __b, __t, __c, __l, __f>::word_size;
198
199 template<typename _UIntType,
200 size_t __w, size_t __n, size_t __m, size_t __r,
201 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
202 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
203 _UIntType __f>
204 constexpr size_t
205 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
206 __s, __b, __t, __c, __l, __f>::state_size;
207
208 template<typename _UIntType,
209 size_t __w, size_t __n, size_t __m, size_t __r,
210 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
211 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
212 _UIntType __f>
213 constexpr size_t
214 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
215 __s, __b, __t, __c, __l, __f>::shift_size;
216
217 template<typename _UIntType,
218 size_t __w, size_t __n, size_t __m, size_t __r,
219 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
220 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
221 _UIntType __f>
222 constexpr size_t
223 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
224 __s, __b, __t, __c, __l, __f>::mask_bits;
225
226 template<typename _UIntType,
227 size_t __w, size_t __n, size_t __m, size_t __r,
228 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
229 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
230 _UIntType __f>
231 constexpr _UIntType
232 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
233 __s, __b, __t, __c, __l, __f>::xor_mask;
234
235 template<typename _UIntType,
236 size_t __w, size_t __n, size_t __m, size_t __r,
237 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
238 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
239 _UIntType __f>
240 constexpr size_t
241 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
242 __s, __b, __t, __c, __l, __f>::tempering_u;
243
244 template<typename _UIntType,
245 size_t __w, size_t __n, size_t __m, size_t __r,
246 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
247 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
248 _UIntType __f>
249 constexpr _UIntType
250 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
251 __s, __b, __t, __c, __l, __f>::tempering_d;
252
253 template<typename _UIntType,
254 size_t __w, size_t __n, size_t __m, size_t __r,
255 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
256 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
257 _UIntType __f>
258 constexpr size_t
259 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
260 __s, __b, __t, __c, __l, __f>::tempering_s;
261
262 template<typename _UIntType,
263 size_t __w, size_t __n, size_t __m, size_t __r,
264 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
265 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
266 _UIntType __f>
267 constexpr _UIntType
268 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
269 __s, __b, __t, __c, __l, __f>::tempering_b;
270
271 template<typename _UIntType,
272 size_t __w, size_t __n, size_t __m, size_t __r,
273 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
274 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
275 _UIntType __f>
276 constexpr size_t
277 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
278 __s, __b, __t, __c, __l, __f>::tempering_t;
279
280 template<typename _UIntType,
281 size_t __w, size_t __n, size_t __m, size_t __r,
282 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
283 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
284 _UIntType __f>
285 constexpr _UIntType
286 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
287 __s, __b, __t, __c, __l, __f>::tempering_c;
288
289 template<typename _UIntType,
290 size_t __w, size_t __n, size_t __m, size_t __r,
291 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
292 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
293 _UIntType __f>
294 constexpr size_t
295 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
296 __s, __b, __t, __c, __l, __f>::tempering_l;
297
298 template<typename _UIntType,
299 size_t __w, size_t __n, size_t __m, size_t __r,
300 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
301 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
302 _UIntType __f>
303 constexpr _UIntType
304 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
305 __s, __b, __t, __c, __l, __f>::
306 initialization_multiplier;
307
308 template<typename _UIntType,
309 size_t __w, size_t __n, size_t __m, size_t __r,
310 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
311 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
312 _UIntType __f>
313 constexpr _UIntType
314 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
315 __s, __b, __t, __c, __l, __f>::default_seed;
316
317 template<typename _UIntType,
318 size_t __w, size_t __n, size_t __m, size_t __r,
319 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
320 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
321 _UIntType __f>
322 void
323 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
324 __s, __b, __t, __c, __l, __f>::
325 seed(result_type __sd)
326 {
327 _M_x[0] = __detail::__mod<_UIntType,
328 __detail::_Shift<_UIntType, __w>::__value>(__sd);
329
330 for (size_t __i = 1; __i < state_size; ++__i)
331 {
332 _UIntType __x = _M_x[__i - 1];
333 __x ^= __x >> (__w - 2);
334 __x *= __f;
335 __x += __detail::__mod<_UIntType, __n>(__i);
336 _M_x[__i] = __detail::__mod<_UIntType,
337 __detail::_Shift<_UIntType, __w>::__value>(__x);
338 }
339 _M_p = state_size;
340 }
341
342 template<typename _UIntType,
343 size_t __w, size_t __n, size_t __m, size_t __r,
344 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
345 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
346 _UIntType __f>
347 template<typename _Sseq>
348 auto
349 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
350 __s, __b, __t, __c, __l, __f>::
351 seed(_Sseq& __q)
352 -> _If_seed_seq<_Sseq>
353 {
354 const _UIntType __upper_mask = (~_UIntType()) << __r;
355 const size_t __k = (__w + 31) / 32;
356 uint_least32_t __arr[__n * __k];
357 __q.generate(__arr + 0, __arr + __n * __k);
358
359 bool __zero = true;
360 for (size_t __i = 0; __i < state_size; ++__i)
361 {
362 _UIntType __factor = 1u;
363 _UIntType __sum = 0u;
364 for (size_t __j = 0; __j < __k; ++__j)
365 {
366 __sum += __arr[__k * __i + __j] * __factor;
367 __factor *= __detail::_Shift<_UIntType, 32>::__value;
368 }
369 _M_x[__i] = __detail::__mod<_UIntType,
370 __detail::_Shift<_UIntType, __w>::__value>(__sum);
371
372 if (__zero)
373 {
374 if (__i == 0)
375 {
376 if ((_M_x[0] & __upper_mask) != 0u)
377 __zero = false;
378 }
379 else if (_M_x[__i] != 0u)
380 __zero = false;
381 }
382 }
383 if (__zero)
384 _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
385 _M_p = state_size;
386 }
387
388 template<typename _UIntType, size_t __w,
389 size_t __n, size_t __m, size_t __r,
390 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
391 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
392 _UIntType __f>
393 void
394 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
395 __s, __b, __t, __c, __l, __f>::
396 _M_gen_rand(void)
397 {
398 const _UIntType __upper_mask = (~_UIntType()) << __r;
399 const _UIntType __lower_mask = ~__upper_mask;
400
401 for (size_t __k = 0; __k < (__n - __m); ++__k)
402 {
403 _UIntType __y = ((_M_x[__k] & __upper_mask)
404 | (_M_x[__k + 1] & __lower_mask));
405 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
406 ^ ((__y & 0x01) ? __a : 0));
407 }
408
409 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
410 {
411 _UIntType __y = ((_M_x[__k] & __upper_mask)
412 | (_M_x[__k + 1] & __lower_mask));
413 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
414 ^ ((__y & 0x01) ? __a : 0));
415 }
416
417 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
418 | (_M_x[0] & __lower_mask));
419 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
420 ^ ((__y & 0x01) ? __a : 0));
421 _M_p = 0;
422 }
423
424 template<typename _UIntType, size_t __w,
425 size_t __n, size_t __m, size_t __r,
426 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
427 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
428 _UIntType __f>
429 void
430 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
431 __s, __b, __t, __c, __l, __f>::
432 discard(unsigned long long __z)
433 {
434 while (__z > state_size - _M_p)
435 {
436 __z -= state_size - _M_p;
437 _M_gen_rand();
438 }
439 _M_p += __z;
440 }
441
442 template<typename _UIntType, size_t __w,
443 size_t __n, size_t __m, size_t __r,
444 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
445 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
446 _UIntType __f>
447 typename
448 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
449 __s, __b, __t, __c, __l, __f>::result_type
450 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
451 __s, __b, __t, __c, __l, __f>::
452 operator()()
453 {
454 // Reload the vector - cost is O(n) amortized over n calls.
455 if (_M_p >= state_size)
456 _M_gen_rand();
457
458 // Calculate o(x(i)).
459 result_type __z = _M_x[_M_p++];
460 __z ^= (__z >> __u) & __d;
461 __z ^= (__z << __s) & __b;
462 __z ^= (__z << __t) & __c;
463 __z ^= (__z >> __l);
464
465 return __z;
466 }
467
468 template<typename _UIntType, size_t __w,
469 size_t __n, size_t __m, size_t __r,
470 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
471 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
472 _UIntType __f, typename _CharT, typename _Traits>
475 const mersenne_twister_engine<_UIntType, __w, __n, __m,
476 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
477 {
478 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
479
480 const typename __ios_base::fmtflags __flags = __os.flags();
481 const _CharT __fill = __os.fill();
482 const _CharT __space = __os.widen(' ');
484 __os.fill(__space);
485
486 for (size_t __i = 0; __i < __n; ++__i)
487 __os << __x._M_x[__i] << __space;
488 __os << __x._M_p;
489
490 __os.flags(__flags);
491 __os.fill(__fill);
492 return __os;
493 }
494
495 template<typename _UIntType, size_t __w,
496 size_t __n, size_t __m, size_t __r,
497 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
498 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
499 _UIntType __f, typename _CharT, typename _Traits>
502 mersenne_twister_engine<_UIntType, __w, __n, __m,
503 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
504 {
505 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
506
507 const typename __ios_base::fmtflags __flags = __is.flags();
509
510 for (size_t __i = 0; __i < __n; ++__i)
511 __is >> __x._M_x[__i];
512 __is >> __x._M_p;
513
514 __is.flags(__flags);
515 return __is;
516 }
517
518
519 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
520 constexpr size_t
521 subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
522
523 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
524 constexpr size_t
525 subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
526
527 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
528 constexpr size_t
529 subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
530
531 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
532 constexpr _UIntType
533 subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
534
535 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
536 void
538 seed(result_type __value)
539 {
541 __lcg(__value == 0u ? default_seed : __value);
542
543 const size_t __n = (__w + 31) / 32;
544
545 for (size_t __i = 0; __i < long_lag; ++__i)
546 {
547 _UIntType __sum = 0u;
548 _UIntType __factor = 1u;
549 for (size_t __j = 0; __j < __n; ++__j)
550 {
551 __sum += __detail::__mod<uint_least32_t,
552 __detail::_Shift<uint_least32_t, 32>::__value>
553 (__lcg()) * __factor;
554 __factor *= __detail::_Shift<_UIntType, 32>::__value;
555 }
556 _M_x[__i] = __detail::__mod<_UIntType,
557 __detail::_Shift<_UIntType, __w>::__value>(__sum);
558 }
559 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
560 _M_p = 0;
561 }
562
563 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
564 template<typename _Sseq>
565 auto
567 seed(_Sseq& __q)
568 -> _If_seed_seq<_Sseq>
569 {
570 const size_t __k = (__w + 31) / 32;
571 uint_least32_t __arr[__r * __k];
572 __q.generate(__arr + 0, __arr + __r * __k);
573
574 for (size_t __i = 0; __i < long_lag; ++__i)
575 {
576 _UIntType __sum = 0u;
577 _UIntType __factor = 1u;
578 for (size_t __j = 0; __j < __k; ++__j)
579 {
580 __sum += __arr[__k * __i + __j] * __factor;
581 __factor *= __detail::_Shift<_UIntType, 32>::__value;
582 }
583 _M_x[__i] = __detail::__mod<_UIntType,
584 __detail::_Shift<_UIntType, __w>::__value>(__sum);
585 }
586 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
587 _M_p = 0;
588 }
589
590 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
591 typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
592 result_type
595 {
596 // Derive short lag index from current index.
597 long __ps = _M_p - short_lag;
598 if (__ps < 0)
599 __ps += long_lag;
600
601 // Calculate new x(i) without overflow or division.
602 // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
603 // cannot overflow.
604 _UIntType __xi;
605 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
606 {
607 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
608 _M_carry = 0;
609 }
610 else
611 {
612 __xi = (__detail::_Shift<_UIntType, __w>::__value
613 - _M_x[_M_p] - _M_carry + _M_x[__ps]);
614 _M_carry = 1;
615 }
616 _M_x[_M_p] = __xi;
617
618 // Adjust current index to loop around in ring buffer.
619 if (++_M_p >= long_lag)
620 _M_p = 0;
621
622 return __xi;
623 }
624
625 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
626 typename _CharT, typename _Traits>
629 const subtract_with_carry_engine<_UIntType,
630 __w, __s, __r>& __x)
631 {
632 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
633
634 const typename __ios_base::fmtflags __flags = __os.flags();
635 const _CharT __fill = __os.fill();
636 const _CharT __space = __os.widen(' ');
638 __os.fill(__space);
639
640 for (size_t __i = 0; __i < __r; ++__i)
641 __os << __x._M_x[__i] << __space;
642 __os << __x._M_carry << __space << __x._M_p;
643
644 __os.flags(__flags);
645 __os.fill(__fill);
646 return __os;
647 }
648
649 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
650 typename _CharT, typename _Traits>
653 subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
654 {
655 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
656
657 const typename __ios_base::fmtflags __flags = __is.flags();
659
660 for (size_t __i = 0; __i < __r; ++__i)
661 __is >> __x._M_x[__i];
662 __is >> __x._M_carry;
663 __is >> __x._M_p;
664
665 __is.flags(__flags);
666 return __is;
667 }
668
669
670 template<typename _RandomNumberEngine, size_t __p, size_t __r>
671 constexpr size_t
672 discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
673
674 template<typename _RandomNumberEngine, size_t __p, size_t __r>
675 constexpr size_t
676 discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
677
678 template<typename _RandomNumberEngine, size_t __p, size_t __r>
679 typename discard_block_engine<_RandomNumberEngine,
680 __p, __r>::result_type
683 {
684 if (_M_n >= used_block)
685 {
686 _M_b.discard(block_size - _M_n);
687 _M_n = 0;
688 }
689 ++_M_n;
690 return _M_b();
691 }
692
693 template<typename _RandomNumberEngine, size_t __p, size_t __r,
694 typename _CharT, typename _Traits>
697 const discard_block_engine<_RandomNumberEngine,
698 __p, __r>& __x)
699 {
700 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
701
702 const typename __ios_base::fmtflags __flags = __os.flags();
703 const _CharT __fill = __os.fill();
704 const _CharT __space = __os.widen(' ');
706 __os.fill(__space);
707
708 __os << __x.base() << __space << __x._M_n;
709
710 __os.flags(__flags);
711 __os.fill(__fill);
712 return __os;
713 }
714
715 template<typename _RandomNumberEngine, size_t __p, size_t __r,
716 typename _CharT, typename _Traits>
719 discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
720 {
721 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
722
723 const typename __ios_base::fmtflags __flags = __is.flags();
725
726 __is >> __x._M_b >> __x._M_n;
727
728 __is.flags(__flags);
729 return __is;
730 }
731
732
733 template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
734 typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
735 result_type
738 {
739 typedef typename _RandomNumberEngine::result_type _Eresult_type;
740 const _Eresult_type __r
741 = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
742 ? _M_b.max() - _M_b.min() + 1 : 0);
743 const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
744 const unsigned __m = __r ? std::__lg(__r) : __edig;
745
747 __ctype;
748 const unsigned __cdig = std::numeric_limits<__ctype>::digits;
749
750 unsigned __n, __n0;
751 __ctype __s0, __s1, __y0, __y1;
752
753 for (size_t __i = 0; __i < 2; ++__i)
754 {
755 __n = (__w + __m - 1) / __m + __i;
756 __n0 = __n - __w % __n;
757 const unsigned __w0 = __w / __n; // __w0 <= __m
758
759 __s0 = 0;
760 __s1 = 0;
761 if (__w0 < __cdig)
762 {
763 __s0 = __ctype(1) << __w0;
764 __s1 = __s0 << 1;
765 }
766
767 __y0 = 0;
768 __y1 = 0;
769 if (__r)
770 {
771 __y0 = __s0 * (__r / __s0);
772 if (__s1)
773 __y1 = __s1 * (__r / __s1);
774
775 if (__r - __y0 <= __y0 / __n)
776 break;
777 }
778 else
779 break;
780 }
781
782 result_type __sum = 0;
783 for (size_t __k = 0; __k < __n0; ++__k)
784 {
785 __ctype __u;
786 do
787 __u = _M_b() - _M_b.min();
788 while (__y0 && __u >= __y0);
789 __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
790 }
791 for (size_t __k = __n0; __k < __n; ++__k)
792 {
793 __ctype __u;
794 do
795 __u = _M_b() - _M_b.min();
796 while (__y1 && __u >= __y1);
797 __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
798 }
799 return __sum;
800 }
801
802
803 template<typename _RandomNumberEngine, size_t __k>
804 constexpr size_t
806
807 template<typename _RandomNumberEngine, size_t __k>
811 {
812 size_t __j = __k * ((_M_y - _M_b.min())
813 / (_M_b.max() - _M_b.min() + 1.0L));
814 _M_y = _M_v[__j];
815 _M_v[__j] = _M_b();
816
817 return _M_y;
818 }
819
820 template<typename _RandomNumberEngine, size_t __k,
821 typename _CharT, typename _Traits>
825 {
826 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
827
828 const typename __ios_base::fmtflags __flags = __os.flags();
829 const _CharT __fill = __os.fill();
830 const _CharT __space = __os.widen(' ');
832 __os.fill(__space);
833
834 __os << __x.base();
835 for (size_t __i = 0; __i < __k; ++__i)
836 __os << __space << __x._M_v[__i];
837 __os << __space << __x._M_y;
838
839 __os.flags(__flags);
840 __os.fill(__fill);
841 return __os;
842 }
843
844 template<typename _RandomNumberEngine, size_t __k,
845 typename _CharT, typename _Traits>
849 {
850 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
851
852 const typename __ios_base::fmtflags __flags = __is.flags();
854
855 __is >> __x._M_b;
856 for (size_t __i = 0; __i < __k; ++__i)
857 __is >> __x._M_v[__i];
858 __is >> __x._M_y;
859
860 __is.flags(__flags);
861 return __is;
862 }
863
864
865 template<typename _IntType, typename _CharT, typename _Traits>
868 const uniform_int_distribution<_IntType>& __x)
869 {
870 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
871
872 const typename __ios_base::fmtflags __flags = __os.flags();
873 const _CharT __fill = __os.fill();
874 const _CharT __space = __os.widen(' ');
876 __os.fill(__space);
877
878 __os << __x.a() << __space << __x.b();
879
880 __os.flags(__flags);
881 __os.fill(__fill);
882 return __os;
883 }
884
885 template<typename _IntType, typename _CharT, typename _Traits>
889 {
890 using param_type
892 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
893
894 const typename __ios_base::fmtflags __flags = __is.flags();
896
897 _IntType __a, __b;
898 if (__is >> __a >> __b)
899 __x.param(param_type(__a, __b));
900
901 __is.flags(__flags);
902 return __is;
903 }
904
905
906 template<typename _RealType>
907 template<typename _ForwardIterator,
908 typename _UniformRandomNumberGenerator>
909 void
911 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
912 _UniformRandomNumberGenerator& __urng,
913 const param_type& __p)
914 {
915 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
916 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
917 __aurng(__urng);
918 auto __range = __p.b() - __p.a();
919 while (__f != __t)
920 *__f++ = __aurng() * __range + __p.a();
921 }
922
923 template<typename _RealType, typename _CharT, typename _Traits>
926 const uniform_real_distribution<_RealType>& __x)
927 {
928 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
929
930 const typename __ios_base::fmtflags __flags = __os.flags();
931 const _CharT __fill = __os.fill();
932 const std::streamsize __precision = __os.precision();
933 const _CharT __space = __os.widen(' ');
935 __os.fill(__space);
937
938 __os << __x.a() << __space << __x.b();
939
940 __os.flags(__flags);
941 __os.fill(__fill);
942 __os.precision(__precision);
943 return __os;
944 }
945
946 template<typename _RealType, typename _CharT, typename _Traits>
950 {
951 using param_type
953 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
954
955 const typename __ios_base::fmtflags __flags = __is.flags();
957
958 _RealType __a, __b;
959 if (__is >> __a >> __b)
960 __x.param(param_type(__a, __b));
961
962 __is.flags(__flags);
963 return __is;
964 }
965
966
967 template<typename _ForwardIterator,
968 typename _UniformRandomNumberGenerator>
969 void
970 std::bernoulli_distribution::
971 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
972 _UniformRandomNumberGenerator& __urng,
973 const param_type& __p)
974 {
975 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
976 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
977 __aurng(__urng);
978 auto __limit = __p.p() * (__aurng.max() - __aurng.min());
979
980 while (__f != __t)
981 *__f++ = (__aurng() - __aurng.min()) < __limit;
982 }
983
984 template<typename _CharT, typename _Traits>
987 const bernoulli_distribution& __x)
988 {
989 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
990
991 const typename __ios_base::fmtflags __flags = __os.flags();
992 const _CharT __fill = __os.fill();
993 const std::streamsize __precision = __os.precision();
995 __os.fill(__os.widen(' '));
997
998 __os << __x.p();
999
1000 __os.flags(__flags);
1001 __os.fill(__fill);
1002 __os.precision(__precision);
1003 return __os;
1004 }
1005
1006
1007 template<typename _IntType>
1008 template<typename _UniformRandomNumberGenerator>
1011 operator()(_UniformRandomNumberGenerator& __urng,
1012 const param_type& __param)
1013 {
1014 // About the epsilon thing see this thread:
1015 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1016 const double __naf =
1018 // The largest _RealType convertible to _IntType.
1019 const double __thr =
1021 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1022 __aurng(__urng);
1023
1024 double __cand;
1025 do
1026 __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1027 while (__cand >= __thr);
1028
1029 return result_type(__cand + __naf);
1030 }
1031
1032 template<typename _IntType>
1033 template<typename _ForwardIterator,
1034 typename _UniformRandomNumberGenerator>
1035 void
1037 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1038 _UniformRandomNumberGenerator& __urng,
1039 const param_type& __param)
1040 {
1041 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1042 // About the epsilon thing see this thread:
1043 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1044 const double __naf =
1046 // The largest _RealType convertible to _IntType.
1047 const double __thr =
1049 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1050 __aurng(__urng);
1051
1052 while (__f != __t)
1053 {
1054 double __cand;
1055 do
1056 __cand = std::floor(std::log(1.0 - __aurng())
1057 / __param._M_log_1_p);
1058 while (__cand >= __thr);
1059
1060 *__f++ = __cand + __naf;
1061 }
1062 }
1063
1064 template<typename _IntType,
1065 typename _CharT, typename _Traits>
1067 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1068 const geometric_distribution<_IntType>& __x)
1069 {
1070 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1071
1072 const typename __ios_base::fmtflags __flags = __os.flags();
1073 const _CharT __fill = __os.fill();
1074 const std::streamsize __precision = __os.precision();
1076 __os.fill(__os.widen(' '));
1078
1079 __os << __x.p();
1080
1081 __os.flags(__flags);
1082 __os.fill(__fill);
1083 __os.precision(__precision);
1084 return __os;
1085 }
1086
1087 template<typename _IntType,
1088 typename _CharT, typename _Traits>
1092 {
1093 using param_type = typename geometric_distribution<_IntType>::param_type;
1094 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1095
1096 const typename __ios_base::fmtflags __flags = __is.flags();
1098
1099 double __p;
1100 if (__is >> __p)
1101 __x.param(param_type(__p));
1102
1103 __is.flags(__flags);
1104 return __is;
1105 }
1106
1107 // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1108 template<typename _IntType>
1109 template<typename _UniformRandomNumberGenerator>
1112 operator()(_UniformRandomNumberGenerator& __urng)
1113 {
1114 const double __y = _M_gd(__urng);
1115
1116 // XXX Is the constructor too slow?
1118 return __poisson(__urng);
1119 }
1120
1121 template<typename _IntType>
1122 template<typename _UniformRandomNumberGenerator>
1125 operator()(_UniformRandomNumberGenerator& __urng,
1126 const param_type& __p)
1127 {
1129 param_type;
1130
1131 const double __y =
1132 _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1133
1135 return __poisson(__urng);
1136 }
1137
1138 template<typename _IntType>
1139 template<typename _ForwardIterator,
1140 typename _UniformRandomNumberGenerator>
1141 void
1142 negative_binomial_distribution<_IntType>::
1143 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1144 _UniformRandomNumberGenerator& __urng)
1145 {
1146 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1147 while (__f != __t)
1148 {
1149 const double __y = _M_gd(__urng);
1150
1151 // XXX Is the constructor too slow?
1153 *__f++ = __poisson(__urng);
1154 }
1155 }
1156
1157 template<typename _IntType>
1158 template<typename _ForwardIterator,
1159 typename _UniformRandomNumberGenerator>
1160 void
1161 negative_binomial_distribution<_IntType>::
1162 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1163 _UniformRandomNumberGenerator& __urng,
1164 const param_type& __p)
1165 {
1166 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1168 __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1169
1170 while (__f != __t)
1171 {
1172 const double __y = _M_gd(__urng, __p2);
1173
1175 *__f++ = __poisson(__urng);
1176 }
1177 }
1178
1179 template<typename _IntType, typename _CharT, typename _Traits>
1181 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1182 const negative_binomial_distribution<_IntType>& __x)
1183 {
1184 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1185
1186 const typename __ios_base::fmtflags __flags = __os.flags();
1187 const _CharT __fill = __os.fill();
1188 const std::streamsize __precision = __os.precision();
1189 const _CharT __space = __os.widen(' ');
1191 __os.fill(__os.widen(' '));
1193
1194 __os << __x.k() << __space << __x.p()
1195 << __space << __x._M_gd;
1196
1197 __os.flags(__flags);
1198 __os.fill(__fill);
1199 __os.precision(__precision);
1200 return __os;
1201 }
1202
1203 template<typename _IntType, typename _CharT, typename _Traits>
1206 negative_binomial_distribution<_IntType>& __x)
1207 {
1208 using param_type
1209 = typename negative_binomial_distribution<_IntType>::param_type;
1210 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1211
1212 const typename __ios_base::fmtflags __flags = __is.flags();
1214
1215 _IntType __k;
1216 double __p;
1217 if (__is >> __k >> __p >> __x._M_gd)
1218 __x.param(param_type(__k, __p));
1219
1220 __is.flags(__flags);
1221 return __is;
1222 }
1223
1224
1225 template<typename _IntType>
1226 void
1227 poisson_distribution<_IntType>::param_type::
1228 _M_initialize()
1229 {
1230#if _GLIBCXX_USE_C99_MATH_TR1
1231 if (_M_mean >= 12)
1232 {
1233 const double __m = std::floor(_M_mean);
1234 _M_lm_thr = std::log(_M_mean);
1235 _M_lfm = std::lgamma(__m + 1);
1236 _M_sm = std::sqrt(__m);
1237
1238 const double __pi_4 = 0.7853981633974483096156608458198757L;
1239 const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1240 / __pi_4));
1241 _M_d = std::round(std::max<double>(6.0, std::min(__m, __dx)));
1242 const double __cx = 2 * __m + _M_d;
1243 _M_scx = std::sqrt(__cx / 2);
1244 _M_1cx = 1 / __cx;
1245
1246 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1247 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1248 / _M_d;
1249 }
1250 else
1251#endif
1252 _M_lm_thr = std::exp(-_M_mean);
1253 }
1254
1255 /**
1256 * A rejection algorithm when mean >= 12 and a simple method based
1257 * upon the multiplication of uniform random variates otherwise.
1258 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1259 * is defined.
1260 *
1261 * Reference:
1262 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1263 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1264 */
1265 template<typename _IntType>
1266 template<typename _UniformRandomNumberGenerator>
1269 operator()(_UniformRandomNumberGenerator& __urng,
1270 const param_type& __param)
1271 {
1272 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1273 __aurng(__urng);
1274#if _GLIBCXX_USE_C99_MATH_TR1
1275 if (__param.mean() >= 12)
1276 {
1277 double __x;
1278
1279 // See comments above...
1280 const double __naf =
1282 const double __thr =
1284
1285 const double __m = std::floor(__param.mean());
1286 // sqrt(pi / 2)
1287 const double __spi_2 = 1.2533141373155002512078826424055226L;
1288 const double __c1 = __param._M_sm * __spi_2;
1289 const double __c2 = __param._M_c2b + __c1;
1290 const double __c3 = __c2 + 1;
1291 const double __c4 = __c3 + 1;
1292 // 1 / 78
1293 const double __178 = 0.0128205128205128205128205128205128L;
1294 // e^(1 / 78)
1295 const double __e178 = 1.0129030479320018583185514777512983L;
1296 const double __c5 = __c4 + __e178;
1297 const double __c = __param._M_cb + __c5;
1298 const double __2cx = 2 * (2 * __m + __param._M_d);
1299
1300 bool __reject = true;
1301 do
1302 {
1303 const double __u = __c * __aurng();
1304 const double __e = -std::log(1.0 - __aurng());
1305
1306 double __w = 0.0;
1307
1308 if (__u <= __c1)
1309 {
1310 const double __n = _M_nd(__urng);
1311 const double __y = -std::abs(__n) * __param._M_sm - 1;
1312 __x = std::floor(__y);
1313 __w = -__n * __n / 2;
1314 if (__x < -__m)
1315 continue;
1316 }
1317 else if (__u <= __c2)
1318 {
1319 const double __n = _M_nd(__urng);
1320 const double __y = 1 + std::abs(__n) * __param._M_scx;
1321 __x = std::ceil(__y);
1322 __w = __y * (2 - __y) * __param._M_1cx;
1323 if (__x > __param._M_d)
1324 continue;
1325 }
1326 else if (__u <= __c3)
1327 // NB: This case not in the book, nor in the Errata,
1328 // but should be ok...
1329 __x = -1;
1330 else if (__u <= __c4)
1331 __x = 0;
1332 else if (__u <= __c5)
1333 {
1334 __x = 1;
1335 // Only in the Errata, see libstdc++/83237.
1336 __w = __178;
1337 }
1338 else
1339 {
1340 const double __v = -std::log(1.0 - __aurng());
1341 const double __y = __param._M_d
1342 + __v * __2cx / __param._M_d;
1343 __x = std::ceil(__y);
1344 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1345 }
1346
1347 __reject = (__w - __e - __x * __param._M_lm_thr
1348 > __param._M_lfm - std::lgamma(__x + __m + 1));
1349
1350 __reject |= __x + __m >= __thr;
1351
1352 } while (__reject);
1353
1354 return result_type(__x + __m + __naf);
1355 }
1356 else
1357#endif
1358 {
1359 _IntType __x = 0;
1360 double __prod = 1.0;
1361
1362 do
1363 {
1364 __prod *= __aurng();
1365 __x += 1;
1366 }
1367 while (__prod > __param._M_lm_thr);
1368
1369 return __x - 1;
1370 }
1371 }
1372
1373 template<typename _IntType>
1374 template<typename _ForwardIterator,
1375 typename _UniformRandomNumberGenerator>
1376 void
1378 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1379 _UniformRandomNumberGenerator& __urng,
1380 const param_type& __param)
1381 {
1382 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1383 // We could duplicate everything from operator()...
1384 while (__f != __t)
1385 *__f++ = this->operator()(__urng, __param);
1386 }
1387
1388 template<typename _IntType,
1389 typename _CharT, typename _Traits>
1392 const poisson_distribution<_IntType>& __x)
1393 {
1394 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1395
1396 const typename __ios_base::fmtflags __flags = __os.flags();
1397 const _CharT __fill = __os.fill();
1398 const std::streamsize __precision = __os.precision();
1399 const _CharT __space = __os.widen(' ');
1401 __os.fill(__space);
1403
1404 __os << __x.mean() << __space << __x._M_nd;
1405
1406 __os.flags(__flags);
1407 __os.fill(__fill);
1408 __os.precision(__precision);
1409 return __os;
1410 }
1411
1412 template<typename _IntType,
1413 typename _CharT, typename _Traits>
1416 poisson_distribution<_IntType>& __x)
1417 {
1418 using param_type = typename poisson_distribution<_IntType>::param_type;
1419 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1420
1421 const typename __ios_base::fmtflags __flags = __is.flags();
1423
1424 double __mean;
1425 if (__is >> __mean >> __x._M_nd)
1426 __x.param(param_type(__mean));
1427
1428 __is.flags(__flags);
1429 return __is;
1430 }
1431
1432
1433 template<typename _IntType>
1434 void
1435 binomial_distribution<_IntType>::param_type::
1436 _M_initialize()
1437 {
1438 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1439
1440 _M_easy = true;
1441
1442#if _GLIBCXX_USE_C99_MATH_TR1
1443 if (_M_t * __p12 >= 8)
1444 {
1445 _M_easy = false;
1446 const double __np = std::floor(_M_t * __p12);
1447 const double __pa = __np / _M_t;
1448 const double __1p = 1 - __pa;
1449
1450 const double __pi_4 = 0.7853981633974483096156608458198757L;
1451 const double __d1x =
1452 std::sqrt(__np * __1p * std::log(32 * __np
1453 / (81 * __pi_4 * __1p)));
1454 _M_d1 = std::round(std::max<double>(1.0, __d1x));
1455 const double __d2x =
1456 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1457 / (__pi_4 * __pa)));
1458 _M_d2 = std::round(std::max<double>(1.0, __d2x));
1459
1460 // sqrt(pi / 2)
1461 const double __spi_2 = 1.2533141373155002512078826424055226L;
1462 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1463 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1464 _M_c = 2 * _M_d1 / __np;
1465 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1466 const double __a12 = _M_a1 + _M_s2 * __spi_2;
1467 const double __s1s = _M_s1 * _M_s1;
1468 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1469 * 2 * __s1s / _M_d1
1470 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1471 const double __s2s = _M_s2 * _M_s2;
1472 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1473 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1474 _M_lf = (std::lgamma(__np + 1)
1475 + std::lgamma(_M_t - __np + 1));
1476 _M_lp1p = std::log(__pa / __1p);
1477
1478 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1479 }
1480 else
1481#endif
1482 _M_q = -std::log(1 - __p12);
1483 }
1484
1485 template<typename _IntType>
1486 template<typename _UniformRandomNumberGenerator>
1488 binomial_distribution<_IntType>::
1489 _M_waiting(_UniformRandomNumberGenerator& __urng,
1490 _IntType __t, double __q)
1491 {
1492 _IntType __x = 0;
1493 double __sum = 0.0;
1494 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1495 __aurng(__urng);
1496
1497 do
1498 {
1499 if (__t == __x)
1500 return __x;
1501 const double __e = -std::log(1.0 - __aurng());
1502 __sum += __e / (__t - __x);
1503 __x += 1;
1504 }
1505 while (__sum <= __q);
1506
1507 return __x - 1;
1508 }
1509
1510 /**
1511 * A rejection algorithm when t * p >= 8 and a simple waiting time
1512 * method - the second in the referenced book - otherwise.
1513 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1514 * is defined.
1515 *
1516 * Reference:
1517 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1518 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1519 */
1520 template<typename _IntType>
1521 template<typename _UniformRandomNumberGenerator>
1524 operator()(_UniformRandomNumberGenerator& __urng,
1525 const param_type& __param)
1526 {
1527 result_type __ret;
1528 const _IntType __t = __param.t();
1529 const double __p = __param.p();
1530 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1531 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1532 __aurng(__urng);
1533
1534#if _GLIBCXX_USE_C99_MATH_TR1
1535 if (!__param._M_easy)
1536 {
1537 double __x;
1538
1539 // See comments above...
1540 const double __naf =
1542 const double __thr =
1544
1545 const double __np = std::floor(__t * __p12);
1546
1547 // sqrt(pi / 2)
1548 const double __spi_2 = 1.2533141373155002512078826424055226L;
1549 const double __a1 = __param._M_a1;
1550 const double __a12 = __a1 + __param._M_s2 * __spi_2;
1551 const double __a123 = __param._M_a123;
1552 const double __s1s = __param._M_s1 * __param._M_s1;
1553 const double __s2s = __param._M_s2 * __param._M_s2;
1554
1555 bool __reject;
1556 do
1557 {
1558 const double __u = __param._M_s * __aurng();
1559
1560 double __v;
1561
1562 if (__u <= __a1)
1563 {
1564 const double __n = _M_nd(__urng);
1565 const double __y = __param._M_s1 * std::abs(__n);
1566 __reject = __y >= __param._M_d1;
1567 if (!__reject)
1568 {
1569 const double __e = -std::log(1.0 - __aurng());
1570 __x = std::floor(__y);
1571 __v = -__e - __n * __n / 2 + __param._M_c;
1572 }
1573 }
1574 else if (__u <= __a12)
1575 {
1576 const double __n = _M_nd(__urng);
1577 const double __y = __param._M_s2 * std::abs(__n);
1578 __reject = __y >= __param._M_d2;
1579 if (!__reject)
1580 {
1581 const double __e = -std::log(1.0 - __aurng());
1582 __x = std::floor(-__y);
1583 __v = -__e - __n * __n / 2;
1584 }
1585 }
1586 else if (__u <= __a123)
1587 {
1588 const double __e1 = -std::log(1.0 - __aurng());
1589 const double __e2 = -std::log(1.0 - __aurng());
1590
1591 const double __y = __param._M_d1
1592 + 2 * __s1s * __e1 / __param._M_d1;
1593 __x = std::floor(__y);
1594 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1595 -__y / (2 * __s1s)));
1596 __reject = false;
1597 }
1598 else
1599 {
1600 const double __e1 = -std::log(1.0 - __aurng());
1601 const double __e2 = -std::log(1.0 - __aurng());
1602
1603 const double __y = __param._M_d2
1604 + 2 * __s2s * __e1 / __param._M_d2;
1605 __x = std::floor(-__y);
1606 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1607 __reject = false;
1608 }
1609
1610 __reject = __reject || __x < -__np || __x > __t - __np;
1611 if (!__reject)
1612 {
1613 const double __lfx =
1614 std::lgamma(__np + __x + 1)
1615 + std::lgamma(__t - (__np + __x) + 1);
1616 __reject = __v > __param._M_lf - __lfx
1617 + __x * __param._M_lp1p;
1618 }
1619
1620 __reject |= __x + __np >= __thr;
1621 }
1622 while (__reject);
1623
1624 __x += __np + __naf;
1625
1626 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1627 __param._M_q);
1628 __ret = _IntType(__x) + __z;
1629 }
1630 else
1631#endif
1632 __ret = _M_waiting(__urng, __t, __param._M_q);
1633
1634 if (__p12 != __p)
1635 __ret = __t - __ret;
1636 return __ret;
1637 }
1638
1639 template<typename _IntType>
1640 template<typename _ForwardIterator,
1641 typename _UniformRandomNumberGenerator>
1642 void
1644 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1645 _UniformRandomNumberGenerator& __urng,
1646 const param_type& __param)
1647 {
1648 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1649 // We could duplicate everything from operator()...
1650 while (__f != __t)
1651 *__f++ = this->operator()(__urng, __param);
1652 }
1653
1654 template<typename _IntType,
1655 typename _CharT, typename _Traits>
1658 const binomial_distribution<_IntType>& __x)
1659 {
1660 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1661
1662 const typename __ios_base::fmtflags __flags = __os.flags();
1663 const _CharT __fill = __os.fill();
1664 const std::streamsize __precision = __os.precision();
1665 const _CharT __space = __os.widen(' ');
1667 __os.fill(__space);
1669
1670 __os << __x.t() << __space << __x.p()
1671 << __space << __x._M_nd;
1672
1673 __os.flags(__flags);
1674 __os.fill(__fill);
1675 __os.precision(__precision);
1676 return __os;
1677 }
1678
1679 template<typename _IntType,
1680 typename _CharT, typename _Traits>
1683 binomial_distribution<_IntType>& __x)
1684 {
1685 using param_type = typename binomial_distribution<_IntType>::param_type;
1686 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1687
1688 const typename __ios_base::fmtflags __flags = __is.flags();
1690
1691 _IntType __t;
1692 double __p;
1693 if (__is >> __t >> __p >> __x._M_nd)
1694 __x.param(param_type(__t, __p));
1695
1696 __is.flags(__flags);
1697 return __is;
1698 }
1699
1700
1701 template<typename _RealType>
1702 template<typename _ForwardIterator,
1703 typename _UniformRandomNumberGenerator>
1704 void
1706 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1707 _UniformRandomNumberGenerator& __urng,
1708 const param_type& __p)
1709 {
1710 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1711 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1712 __aurng(__urng);
1713 while (__f != __t)
1714 *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1715 }
1716
1717 template<typename _RealType, typename _CharT, typename _Traits>
1720 const exponential_distribution<_RealType>& __x)
1721 {
1722 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1723
1724 const typename __ios_base::fmtflags __flags = __os.flags();
1725 const _CharT __fill = __os.fill();
1726 const std::streamsize __precision = __os.precision();
1728 __os.fill(__os.widen(' '));
1730
1731 __os << __x.lambda();
1732
1733 __os.flags(__flags);
1734 __os.fill(__fill);
1735 __os.precision(__precision);
1736 return __os;
1737 }
1738
1739 template<typename _RealType, typename _CharT, typename _Traits>
1743 {
1744 using param_type
1746 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1747
1748 const typename __ios_base::fmtflags __flags = __is.flags();
1750
1751 _RealType __lambda;
1752 if (__is >> __lambda)
1753 __x.param(param_type(__lambda));
1754
1755 __is.flags(__flags);
1756 return __is;
1757 }
1758
1759
1760 /**
1761 * Polar method due to Marsaglia.
1762 *
1763 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1764 * New York, 1986, Ch. V, Sect. 4.4.
1765 */
1766 template<typename _RealType>
1767 template<typename _UniformRandomNumberGenerator>
1770 operator()(_UniformRandomNumberGenerator& __urng,
1771 const param_type& __param)
1772 {
1773 result_type __ret;
1774 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1775 __aurng(__urng);
1776
1777 if (_M_saved_available)
1778 {
1779 _M_saved_available = false;
1780 __ret = _M_saved;
1781 }
1782 else
1783 {
1784 result_type __x, __y, __r2;
1785 do
1786 {
1787 __x = result_type(2.0) * __aurng() - 1.0;
1788 __y = result_type(2.0) * __aurng() - 1.0;
1789 __r2 = __x * __x + __y * __y;
1790 }
1791 while (__r2 > 1.0 || __r2 == 0.0);
1792
1793 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1794 _M_saved = __x * __mult;
1795 _M_saved_available = true;
1796 __ret = __y * __mult;
1797 }
1798
1799 __ret = __ret * __param.stddev() + __param.mean();
1800 return __ret;
1801 }
1802
1803 template<typename _RealType>
1804 template<typename _ForwardIterator,
1805 typename _UniformRandomNumberGenerator>
1806 void
1808 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1809 _UniformRandomNumberGenerator& __urng,
1810 const param_type& __param)
1811 {
1812 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1813
1814 if (__f == __t)
1815 return;
1816
1817 if (_M_saved_available)
1818 {
1819 _M_saved_available = false;
1820 *__f++ = _M_saved * __param.stddev() + __param.mean();
1821
1822 if (__f == __t)
1823 return;
1824 }
1825
1826 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1827 __aurng(__urng);
1828
1829 while (__f + 1 < __t)
1830 {
1831 result_type __x, __y, __r2;
1832 do
1833 {
1834 __x = result_type(2.0) * __aurng() - 1.0;
1835 __y = result_type(2.0) * __aurng() - 1.0;
1836 __r2 = __x * __x + __y * __y;
1837 }
1838 while (__r2 > 1.0 || __r2 == 0.0);
1839
1840 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1841 *__f++ = __y * __mult * __param.stddev() + __param.mean();
1842 *__f++ = __x * __mult * __param.stddev() + __param.mean();
1843 }
1844
1845 if (__f != __t)
1846 {
1847 result_type __x, __y, __r2;
1848 do
1849 {
1850 __x = result_type(2.0) * __aurng() - 1.0;
1851 __y = result_type(2.0) * __aurng() - 1.0;
1852 __r2 = __x * __x + __y * __y;
1853 }
1854 while (__r2 > 1.0 || __r2 == 0.0);
1855
1856 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1857 _M_saved = __x * __mult;
1858 _M_saved_available = true;
1859 *__f = __y * __mult * __param.stddev() + __param.mean();
1860 }
1861 }
1862
1863 template<typename _RealType>
1864 bool
1867 {
1868 if (__d1._M_param == __d2._M_param
1869 && __d1._M_saved_available == __d2._M_saved_available)
1870 {
1871 if (__d1._M_saved_available
1872 && __d1._M_saved == __d2._M_saved)
1873 return true;
1874 else if(!__d1._M_saved_available)
1875 return true;
1876 else
1877 return false;
1878 }
1879 else
1880 return false;
1881 }
1882
1883 template<typename _RealType, typename _CharT, typename _Traits>
1886 const normal_distribution<_RealType>& __x)
1887 {
1888 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1889
1890 const typename __ios_base::fmtflags __flags = __os.flags();
1891 const _CharT __fill = __os.fill();
1892 const std::streamsize __precision = __os.precision();
1893 const _CharT __space = __os.widen(' ');
1895 __os.fill(__space);
1897
1898 __os << __x.mean() << __space << __x.stddev()
1899 << __space << __x._M_saved_available;
1900 if (__x._M_saved_available)
1901 __os << __space << __x._M_saved;
1902
1903 __os.flags(__flags);
1904 __os.fill(__fill);
1905 __os.precision(__precision);
1906 return __os;
1907 }
1908
1909 template<typename _RealType, typename _CharT, typename _Traits>
1912 normal_distribution<_RealType>& __x)
1913 {
1914 using param_type = typename normal_distribution<_RealType>::param_type;
1915 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1916
1917 const typename __ios_base::fmtflags __flags = __is.flags();
1919
1920 double __mean, __stddev;
1921 bool __saved_avail;
1922 if (__is >> __mean >> __stddev >> __saved_avail)
1923 {
1924 if (__saved_avail && (__is >> __x._M_saved))
1925 {
1926 __x._M_saved_available = __saved_avail;
1927 __x.param(param_type(__mean, __stddev));
1928 }
1929 }
1930
1931 __is.flags(__flags);
1932 return __is;
1933 }
1934
1935
1936 template<typename _RealType>
1937 template<typename _ForwardIterator,
1938 typename _UniformRandomNumberGenerator>
1939 void
1940 lognormal_distribution<_RealType>::
1941 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1942 _UniformRandomNumberGenerator& __urng,
1943 const param_type& __p)
1944 {
1945 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1946 while (__f != __t)
1947 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
1948 }
1949
1950 template<typename _RealType, typename _CharT, typename _Traits>
1953 const lognormal_distribution<_RealType>& __x)
1954 {
1955 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1956
1957 const typename __ios_base::fmtflags __flags = __os.flags();
1958 const _CharT __fill = __os.fill();
1959 const std::streamsize __precision = __os.precision();
1960 const _CharT __space = __os.widen(' ');
1962 __os.fill(__space);
1964
1965 __os << __x.m() << __space << __x.s()
1966 << __space << __x._M_nd;
1967
1968 __os.flags(__flags);
1969 __os.fill(__fill);
1970 __os.precision(__precision);
1971 return __os;
1972 }
1973
1974 template<typename _RealType, typename _CharT, typename _Traits>
1977 lognormal_distribution<_RealType>& __x)
1978 {
1979 using param_type
1980 = typename lognormal_distribution<_RealType>::param_type;
1981 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1982
1983 const typename __ios_base::fmtflags __flags = __is.flags();
1985
1986 _RealType __m, __s;
1987 if (__is >> __m >> __s >> __x._M_nd)
1988 __x.param(param_type(__m, __s));
1989
1990 __is.flags(__flags);
1991 return __is;
1992 }
1993
1994 template<typename _RealType>
1995 template<typename _ForwardIterator,
1996 typename _UniformRandomNumberGenerator>
1997 void
1999 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2000 _UniformRandomNumberGenerator& __urng)
2001 {
2002 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2003 while (__f != __t)
2004 *__f++ = 2 * _M_gd(__urng);
2005 }
2006
2007 template<typename _RealType>
2008 template<typename _ForwardIterator,
2009 typename _UniformRandomNumberGenerator>
2010 void
2012 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2013 _UniformRandomNumberGenerator& __urng,
2014 const typename
2016 {
2017 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2018 while (__f != __t)
2019 *__f++ = 2 * _M_gd(__urng, __p);
2020 }
2021
2022 template<typename _RealType, typename _CharT, typename _Traits>
2025 const chi_squared_distribution<_RealType>& __x)
2026 {
2027 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2028
2029 const typename __ios_base::fmtflags __flags = __os.flags();
2030 const _CharT __fill = __os.fill();
2031 const std::streamsize __precision = __os.precision();
2032 const _CharT __space = __os.widen(' ');
2034 __os.fill(__space);
2036
2037 __os << __x.n() << __space << __x._M_gd;
2038
2039 __os.flags(__flags);
2040 __os.fill(__fill);
2041 __os.precision(__precision);
2042 return __os;
2043 }
2044
2045 template<typename _RealType, typename _CharT, typename _Traits>
2048 chi_squared_distribution<_RealType>& __x)
2049 {
2050 using param_type
2051 = typename chi_squared_distribution<_RealType>::param_type;
2052 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2053
2054 const typename __ios_base::fmtflags __flags = __is.flags();
2056
2057 _RealType __n;
2058 if (__is >> __n >> __x._M_gd)
2059 __x.param(param_type(__n));
2060
2061 __is.flags(__flags);
2062 return __is;
2063 }
2064
2065
2066 template<typename _RealType>
2067 template<typename _UniformRandomNumberGenerator>
2070 operator()(_UniformRandomNumberGenerator& __urng,
2071 const param_type& __p)
2072 {
2073 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2074 __aurng(__urng);
2075 _RealType __u;
2076 do
2077 __u = __aurng();
2078 while (__u == 0.5);
2079
2080 const _RealType __pi = 3.1415926535897932384626433832795029L;
2081 return __p.a() + __p.b() * std::tan(__pi * __u);
2082 }
2083
2084 template<typename _RealType>
2085 template<typename _ForwardIterator,
2086 typename _UniformRandomNumberGenerator>
2087 void
2089 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2090 _UniformRandomNumberGenerator& __urng,
2091 const param_type& __p)
2092 {
2093 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2094 const _RealType __pi = 3.1415926535897932384626433832795029L;
2095 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2096 __aurng(__urng);
2097 while (__f != __t)
2098 {
2099 _RealType __u;
2100 do
2101 __u = __aurng();
2102 while (__u == 0.5);
2103
2104 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2105 }
2106 }
2107
2108 template<typename _RealType, typename _CharT, typename _Traits>
2111 const cauchy_distribution<_RealType>& __x)
2112 {
2113 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2114
2115 const typename __ios_base::fmtflags __flags = __os.flags();
2116 const _CharT __fill = __os.fill();
2117 const std::streamsize __precision = __os.precision();
2118 const _CharT __space = __os.widen(' ');
2120 __os.fill(__space);
2122
2123 __os << __x.a() << __space << __x.b();
2124
2125 __os.flags(__flags);
2126 __os.fill(__fill);
2127 __os.precision(__precision);
2128 return __os;
2129 }
2130
2131 template<typename _RealType, typename _CharT, typename _Traits>
2135 {
2136 using param_type = typename cauchy_distribution<_RealType>::param_type;
2137 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2138
2139 const typename __ios_base::fmtflags __flags = __is.flags();
2141
2142 _RealType __a, __b;
2143 if (__is >> __a >> __b)
2144 __x.param(param_type(__a, __b));
2145
2146 __is.flags(__flags);
2147 return __is;
2148 }
2149
2150
2151 template<typename _RealType>
2152 template<typename _ForwardIterator,
2153 typename _UniformRandomNumberGenerator>
2154 void
2156 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2157 _UniformRandomNumberGenerator& __urng)
2158 {
2159 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2160 while (__f != __t)
2161 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2162 }
2163
2164 template<typename _RealType>
2165 template<typename _ForwardIterator,
2166 typename _UniformRandomNumberGenerator>
2167 void
2169 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2170 _UniformRandomNumberGenerator& __urng,
2171 const param_type& __p)
2172 {
2173 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2175 param_type;
2176 param_type __p1(__p.m() / 2);
2177 param_type __p2(__p.n() / 2);
2178 while (__f != __t)
2179 *__f++ = ((_M_gd_x(__urng, __p1) * n())
2180 / (_M_gd_y(__urng, __p2) * m()));
2181 }
2182
2183 template<typename _RealType, typename _CharT, typename _Traits>
2186 const fisher_f_distribution<_RealType>& __x)
2187 {
2188 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2189
2190 const typename __ios_base::fmtflags __flags = __os.flags();
2191 const _CharT __fill = __os.fill();
2192 const std::streamsize __precision = __os.precision();
2193 const _CharT __space = __os.widen(' ');
2195 __os.fill(__space);
2197
2198 __os << __x.m() << __space << __x.n()
2199 << __space << __x._M_gd_x << __space << __x._M_gd_y;
2200
2201 __os.flags(__flags);
2202 __os.fill(__fill);
2203 __os.precision(__precision);
2204 return __os;
2205 }
2206
2207 template<typename _RealType, typename _CharT, typename _Traits>
2210 fisher_f_distribution<_RealType>& __x)
2211 {
2212 using param_type
2213 = typename fisher_f_distribution<_RealType>::param_type;
2214 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2215
2216 const typename __ios_base::fmtflags __flags = __is.flags();
2218
2219 _RealType __m, __n;
2220 if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y)
2221 __x.param(param_type(__m, __n));
2222
2223 __is.flags(__flags);
2224 return __is;
2225 }
2226
2227
2228 template<typename _RealType>
2229 template<typename _ForwardIterator,
2230 typename _UniformRandomNumberGenerator>
2231 void
2233 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2234 _UniformRandomNumberGenerator& __urng)
2235 {
2236 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2237 while (__f != __t)
2238 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2239 }
2240
2241 template<typename _RealType>
2242 template<typename _ForwardIterator,
2243 typename _UniformRandomNumberGenerator>
2244 void
2246 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2247 _UniformRandomNumberGenerator& __urng,
2248 const param_type& __p)
2249 {
2250 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2252 __p2(__p.n() / 2, 2);
2253 while (__f != __t)
2254 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2255 }
2256
2257 template<typename _RealType, typename _CharT, typename _Traits>
2260 const student_t_distribution<_RealType>& __x)
2261 {
2262 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2263
2264 const typename __ios_base::fmtflags __flags = __os.flags();
2265 const _CharT __fill = __os.fill();
2266 const std::streamsize __precision = __os.precision();
2267 const _CharT __space = __os.widen(' ');
2269 __os.fill(__space);
2271
2272 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2273
2274 __os.flags(__flags);
2275 __os.fill(__fill);
2276 __os.precision(__precision);
2277 return __os;
2278 }
2279
2280 template<typename _RealType, typename _CharT, typename _Traits>
2283 student_t_distribution<_RealType>& __x)
2284 {
2285 using param_type
2286 = typename student_t_distribution<_RealType>::param_type;
2287 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2288
2289 const typename __ios_base::fmtflags __flags = __is.flags();
2291
2292 _RealType __n;
2293 if (__is >> __n >> __x._M_nd >> __x._M_gd)
2294 __x.param(param_type(__n));
2295
2296 __is.flags(__flags);
2297 return __is;
2298 }
2299
2300
2301 template<typename _RealType>
2302 void
2303 gamma_distribution<_RealType>::param_type::
2304 _M_initialize()
2305 {
2306 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2307
2308 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2309 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2310 }
2311
2312 /**
2313 * Marsaglia, G. and Tsang, W. W.
2314 * "A Simple Method for Generating Gamma Variables"
2315 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2316 */
2317 template<typename _RealType>
2318 template<typename _UniformRandomNumberGenerator>
2321 operator()(_UniformRandomNumberGenerator& __urng,
2322 const param_type& __param)
2323 {
2324 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2325 __aurng(__urng);
2326
2327 result_type __u, __v, __n;
2328 const result_type __a1 = (__param._M_malpha
2329 - _RealType(1.0) / _RealType(3.0));
2330
2331 do
2332 {
2333 do
2334 {
2335 __n = _M_nd(__urng);
2336 __v = result_type(1.0) + __param._M_a2 * __n;
2337 }
2338 while (__v <= 0.0);
2339
2340 __v = __v * __v * __v;
2341 __u = __aurng();
2342 }
2343 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2344 && (std::log(__u) > (0.5 * __n * __n + __a1
2345 * (1.0 - __v + std::log(__v)))));
2346
2347 if (__param.alpha() == __param._M_malpha)
2348 return __a1 * __v * __param.beta();
2349 else
2350 {
2351 do
2352 __u = __aurng();
2353 while (__u == 0.0);
2354
2355 return (std::pow(__u, result_type(1.0) / __param.alpha())
2356 * __a1 * __v * __param.beta());
2357 }
2358 }
2359
2360 template<typename _RealType>
2361 template<typename _ForwardIterator,
2362 typename _UniformRandomNumberGenerator>
2363 void
2365 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2366 _UniformRandomNumberGenerator& __urng,
2367 const param_type& __param)
2368 {
2369 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2370 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2371 __aurng(__urng);
2372
2373 result_type __u, __v, __n;
2374 const result_type __a1 = (__param._M_malpha
2375 - _RealType(1.0) / _RealType(3.0));
2376
2377 if (__param.alpha() == __param._M_malpha)
2378 while (__f != __t)
2379 {
2380 do
2381 {
2382 do
2383 {
2384 __n = _M_nd(__urng);
2385 __v = result_type(1.0) + __param._M_a2 * __n;
2386 }
2387 while (__v <= 0.0);
2388
2389 __v = __v * __v * __v;
2390 __u = __aurng();
2391 }
2392 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2393 && (std::log(__u) > (0.5 * __n * __n + __a1
2394 * (1.0 - __v + std::log(__v)))));
2395
2396 *__f++ = __a1 * __v * __param.beta();
2397 }
2398 else
2399 while (__f != __t)
2400 {
2401 do
2402 {
2403 do
2404 {
2405 __n = _M_nd(__urng);
2406 __v = result_type(1.0) + __param._M_a2 * __n;
2407 }
2408 while (__v <= 0.0);
2409
2410 __v = __v * __v * __v;
2411 __u = __aurng();
2412 }
2413 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2414 && (std::log(__u) > (0.5 * __n * __n + __a1
2415 * (1.0 - __v + std::log(__v)))));
2416
2417 do
2418 __u = __aurng();
2419 while (__u == 0.0);
2420
2421 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2422 * __a1 * __v * __param.beta());
2423 }
2424 }
2425
2426 template<typename _RealType, typename _CharT, typename _Traits>
2429 const gamma_distribution<_RealType>& __x)
2430 {
2431 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2432
2433 const typename __ios_base::fmtflags __flags = __os.flags();
2434 const _CharT __fill = __os.fill();
2435 const std::streamsize __precision = __os.precision();
2436 const _CharT __space = __os.widen(' ');
2438 __os.fill(__space);
2440
2441 __os << __x.alpha() << __space << __x.beta()
2442 << __space << __x._M_nd;
2443
2444 __os.flags(__flags);
2445 __os.fill(__fill);
2446 __os.precision(__precision);
2447 return __os;
2448 }
2449
2450 template<typename _RealType, typename _CharT, typename _Traits>
2453 gamma_distribution<_RealType>& __x)
2454 {
2455 using param_type = typename gamma_distribution<_RealType>::param_type;
2456 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2457
2458 const typename __ios_base::fmtflags __flags = __is.flags();
2460
2461 _RealType __alpha_val, __beta_val;
2462 if (__is >> __alpha_val >> __beta_val >> __x._M_nd)
2463 __x.param(param_type(__alpha_val, __beta_val));
2464
2465 __is.flags(__flags);
2466 return __is;
2467 }
2468
2469
2470 template<typename _RealType>
2471 template<typename _UniformRandomNumberGenerator>
2474 operator()(_UniformRandomNumberGenerator& __urng,
2475 const param_type& __p)
2476 {
2477 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2478 __aurng(__urng);
2479 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2480 result_type(1) / __p.a());
2481 }
2482
2483 template<typename _RealType>
2484 template<typename _ForwardIterator,
2485 typename _UniformRandomNumberGenerator>
2486 void
2488 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2489 _UniformRandomNumberGenerator& __urng,
2490 const param_type& __p)
2491 {
2492 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2493 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2494 __aurng(__urng);
2495 auto __inv_a = result_type(1) / __p.a();
2496
2497 while (__f != __t)
2498 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2499 __inv_a);
2500 }
2501
2502 template<typename _RealType, typename _CharT, typename _Traits>
2505 const weibull_distribution<_RealType>& __x)
2506 {
2507 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2508
2509 const typename __ios_base::fmtflags __flags = __os.flags();
2510 const _CharT __fill = __os.fill();
2511 const std::streamsize __precision = __os.precision();
2512 const _CharT __space = __os.widen(' ');
2514 __os.fill(__space);
2516
2517 __os << __x.a() << __space << __x.b();
2518
2519 __os.flags(__flags);
2520 __os.fill(__fill);
2521 __os.precision(__precision);
2522 return __os;
2523 }
2524
2525 template<typename _RealType, typename _CharT, typename _Traits>
2529 {
2530 using param_type = typename weibull_distribution<_RealType>::param_type;
2531 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2532
2533 const typename __ios_base::fmtflags __flags = __is.flags();
2535
2536 _RealType __a, __b;
2537 if (__is >> __a >> __b)
2538 __x.param(param_type(__a, __b));
2539
2540 __is.flags(__flags);
2541 return __is;
2542 }
2543
2544
2545 template<typename _RealType>
2546 template<typename _UniformRandomNumberGenerator>
2549 operator()(_UniformRandomNumberGenerator& __urng,
2550 const param_type& __p)
2551 {
2552 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2553 __aurng(__urng);
2554 return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2555 - __aurng()));
2556 }
2557
2558 template<typename _RealType>
2559 template<typename _ForwardIterator,
2560 typename _UniformRandomNumberGenerator>
2561 void
2563 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2564 _UniformRandomNumberGenerator& __urng,
2565 const param_type& __p)
2566 {
2567 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2568 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2569 __aurng(__urng);
2570
2571 while (__f != __t)
2572 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2573 - __aurng()));
2574 }
2575
2576 template<typename _RealType, typename _CharT, typename _Traits>
2579 const extreme_value_distribution<_RealType>& __x)
2580 {
2581 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2582
2583 const typename __ios_base::fmtflags __flags = __os.flags();
2584 const _CharT __fill = __os.fill();
2585 const std::streamsize __precision = __os.precision();
2586 const _CharT __space = __os.widen(' ');
2588 __os.fill(__space);
2590
2591 __os << __x.a() << __space << __x.b();
2592
2593 __os.flags(__flags);
2594 __os.fill(__fill);
2595 __os.precision(__precision);
2596 return __os;
2597 }
2598
2599 template<typename _RealType, typename _CharT, typename _Traits>
2603 {
2604 using param_type
2606 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2607
2608 const typename __ios_base::fmtflags __flags = __is.flags();
2610
2611 _RealType __a, __b;
2612 if (__is >> __a >> __b)
2613 __x.param(param_type(__a, __b));
2614
2615 __is.flags(__flags);
2616 return __is;
2617 }
2618
2619
2620 template<typename _IntType>
2621 void
2622 discrete_distribution<_IntType>::param_type::
2623 _M_initialize()
2624 {
2625 if (_M_prob.size() < 2)
2626 {
2627 _M_prob.clear();
2628 return;
2629 }
2630
2631 const double __sum = std::accumulate(_M_prob.begin(),
2632 _M_prob.end(), 0.0);
2633 // Now normalize the probabilites.
2634 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2635 __sum);
2636 // Accumulate partial sums.
2637 _M_cp.reserve(_M_prob.size());
2638 std::partial_sum(_M_prob.begin(), _M_prob.end(),
2639 std::back_inserter(_M_cp));
2640 // Make sure the last cumulative probability is one.
2641 _M_cp[_M_cp.size() - 1] = 1.0;
2642 }
2643
2644 template<typename _IntType>
2645 template<typename _Func>
2646 discrete_distribution<_IntType>::param_type::
2647 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2648 : _M_prob(), _M_cp()
2649 {
2650 const size_t __n = __nw == 0 ? 1 : __nw;
2651 const double __delta = (__xmax - __xmin) / __n;
2652
2653 _M_prob.reserve(__n);
2654 for (size_t __k = 0; __k < __nw; ++__k)
2655 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2656
2657 _M_initialize();
2658 }
2659
2660 template<typename _IntType>
2661 template<typename _UniformRandomNumberGenerator>
2662 typename discrete_distribution<_IntType>::result_type
2663 discrete_distribution<_IntType>::
2664 operator()(_UniformRandomNumberGenerator& __urng,
2665 const param_type& __param)
2666 {
2667 if (__param._M_cp.empty())
2668 return result_type(0);
2669
2670 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2671 __aurng(__urng);
2672
2673 const double __p = __aurng();
2674 auto __pos = std::lower_bound(__param._M_cp.begin(),
2675 __param._M_cp.end(), __p);
2676
2677 return __pos - __param._M_cp.begin();
2678 }
2679
2680 template<typename _IntType>
2681 template<typename _ForwardIterator,
2682 typename _UniformRandomNumberGenerator>
2683 void
2684 discrete_distribution<_IntType>::
2685 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2686 _UniformRandomNumberGenerator& __urng,
2687 const param_type& __param)
2688 {
2689 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2690
2691 if (__param._M_cp.empty())
2692 {
2693 while (__f != __t)
2694 *__f++ = result_type(0);
2695 return;
2696 }
2697
2698 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2699 __aurng(__urng);
2700
2701 while (__f != __t)
2702 {
2703 const double __p = __aurng();
2704 auto __pos = std::lower_bound(__param._M_cp.begin(),
2705 __param._M_cp.end(), __p);
2706
2707 *__f++ = __pos - __param._M_cp.begin();
2708 }
2709 }
2710
2711 template<typename _IntType, typename _CharT, typename _Traits>
2714 const discrete_distribution<_IntType>& __x)
2715 {
2716 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2717
2718 const typename __ios_base::fmtflags __flags = __os.flags();
2719 const _CharT __fill = __os.fill();
2720 const std::streamsize __precision = __os.precision();
2721 const _CharT __space = __os.widen(' ');
2723 __os.fill(__space);
2725
2726 std::vector<double> __prob = __x.probabilities();
2727 __os << __prob.size();
2728 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2729 __os << __space << *__dit;
2730
2731 __os.flags(__flags);
2732 __os.fill(__fill);
2733 __os.precision(__precision);
2734 return __os;
2735 }
2736
2737namespace __detail
2738{
2739 template<typename _ValT, typename _CharT, typename _Traits>
2740 basic_istream<_CharT, _Traits>&
2741 __extract_params(basic_istream<_CharT, _Traits>& __is,
2742 vector<_ValT>& __vals, size_t __n)
2743 {
2744 __vals.reserve(__n);
2745 while (__n--)
2746 {
2747 _ValT __val;
2748 if (__is >> __val)
2749 __vals.push_back(__val);
2750 else
2751 break;
2752 }
2753 return __is;
2754 }
2755} // namespace __detail
2756
2757 template<typename _IntType, typename _CharT, typename _Traits>
2760 discrete_distribution<_IntType>& __x)
2761 {
2762 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2763
2764 const typename __ios_base::fmtflags __flags = __is.flags();
2766
2767 size_t __n;
2768 if (__is >> __n)
2769 {
2770 std::vector<double> __prob_vec;
2771 if (__detail::__extract_params(__is, __prob_vec, __n))
2772 __x.param({__prob_vec.begin(), __prob_vec.end()});
2773 }
2774
2775 __is.flags(__flags);
2776 return __is;
2777 }
2778
2779
2780 template<typename _RealType>
2781 void
2782 piecewise_constant_distribution<_RealType>::param_type::
2783 _M_initialize()
2784 {
2785 if (_M_int.size() < 2
2786 || (_M_int.size() == 2
2787 && _M_int[0] == _RealType(0)
2788 && _M_int[1] == _RealType(1)))
2789 {
2790 _M_int.clear();
2791 _M_den.clear();
2792 return;
2793 }
2794
2795 const double __sum = std::accumulate(_M_den.begin(),
2796 _M_den.end(), 0.0);
2797
2798 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2799 __sum);
2800
2801 _M_cp.reserve(_M_den.size());
2802 std::partial_sum(_M_den.begin(), _M_den.end(),
2803 std::back_inserter(_M_cp));
2804
2805 // Make sure the last cumulative probability is one.
2806 _M_cp[_M_cp.size() - 1] = 1.0;
2807
2808 for (size_t __k = 0; __k < _M_den.size(); ++__k)
2809 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2810 }
2811
2812 template<typename _RealType>
2813 template<typename _InputIteratorB, typename _InputIteratorW>
2814 piecewise_constant_distribution<_RealType>::param_type::
2815 param_type(_InputIteratorB __bbegin,
2816 _InputIteratorB __bend,
2817 _InputIteratorW __wbegin)
2818 : _M_int(), _M_den(), _M_cp()
2819 {
2820 if (__bbegin != __bend)
2821 {
2822 for (;;)
2823 {
2824 _M_int.push_back(*__bbegin);
2825 ++__bbegin;
2826 if (__bbegin == __bend)
2827 break;
2828
2829 _M_den.push_back(*__wbegin);
2830 ++__wbegin;
2831 }
2832 }
2833
2834 _M_initialize();
2835 }
2836
2837 template<typename _RealType>
2838 template<typename _Func>
2839 piecewise_constant_distribution<_RealType>::param_type::
2840 param_type(initializer_list<_RealType> __bl, _Func __fw)
2841 : _M_int(), _M_den(), _M_cp()
2842 {
2843 _M_int.reserve(__bl.size());
2844 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2845 _M_int.push_back(*__biter);
2846
2847 _M_den.reserve(_M_int.size() - 1);
2848 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2849 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2850
2851 _M_initialize();
2852 }
2853
2854 template<typename _RealType>
2855 template<typename _Func>
2856 piecewise_constant_distribution<_RealType>::param_type::
2857 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2858 : _M_int(), _M_den(), _M_cp()
2859 {
2860 const size_t __n = __nw == 0 ? 1 : __nw;
2861 const _RealType __delta = (__xmax - __xmin) / __n;
2862
2863 _M_int.reserve(__n + 1);
2864 for (size_t __k = 0; __k <= __nw; ++__k)
2865 _M_int.push_back(__xmin + __k * __delta);
2866
2867 _M_den.reserve(__n);
2868 for (size_t __k = 0; __k < __nw; ++__k)
2869 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2870
2871 _M_initialize();
2872 }
2873
2874 template<typename _RealType>
2875 template<typename _UniformRandomNumberGenerator>
2876 typename piecewise_constant_distribution<_RealType>::result_type
2877 piecewise_constant_distribution<_RealType>::
2878 operator()(_UniformRandomNumberGenerator& __urng,
2879 const param_type& __param)
2880 {
2881 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2882 __aurng(__urng);
2883
2884 const double __p = __aurng();
2885 if (__param._M_cp.empty())
2886 return __p;
2887
2888 auto __pos = std::lower_bound(__param._M_cp.begin(),
2889 __param._M_cp.end(), __p);
2890 const size_t __i = __pos - __param._M_cp.begin();
2891
2892 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2893
2894 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2895 }
2896
2897 template<typename _RealType>
2898 template<typename _ForwardIterator,
2899 typename _UniformRandomNumberGenerator>
2900 void
2901 piecewise_constant_distribution<_RealType>::
2902 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2903 _UniformRandomNumberGenerator& __urng,
2904 const param_type& __param)
2905 {
2906 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2907 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2908 __aurng(__urng);
2909
2910 if (__param._M_cp.empty())
2911 {
2912 while (__f != __t)
2913 *__f++ = __aurng();
2914 return;
2915 }
2916
2917 while (__f != __t)
2918 {
2919 const double __p = __aurng();
2920
2921 auto __pos = std::lower_bound(__param._M_cp.begin(),
2922 __param._M_cp.end(), __p);
2923 const size_t __i = __pos - __param._M_cp.begin();
2924
2925 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2926
2927 *__f++ = (__param._M_int[__i]
2928 + (__p - __pref) / __param._M_den[__i]);
2929 }
2930 }
2931
2932 template<typename _RealType, typename _CharT, typename _Traits>
2935 const piecewise_constant_distribution<_RealType>& __x)
2936 {
2937 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2938
2939 const typename __ios_base::fmtflags __flags = __os.flags();
2940 const _CharT __fill = __os.fill();
2941 const std::streamsize __precision = __os.precision();
2942 const _CharT __space = __os.widen(' ');
2944 __os.fill(__space);
2946
2947 std::vector<_RealType> __int = __x.intervals();
2948 __os << __int.size() - 1;
2949
2950 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2951 __os << __space << *__xit;
2952
2953 std::vector<double> __den = __x.densities();
2954 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2955 __os << __space << *__dit;
2956
2957 __os.flags(__flags);
2958 __os.fill(__fill);
2959 __os.precision(__precision);
2960 return __os;
2961 }
2962
2963 template<typename _RealType, typename _CharT, typename _Traits>
2966 piecewise_constant_distribution<_RealType>& __x)
2967 {
2968 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2969
2970 const typename __ios_base::fmtflags __flags = __is.flags();
2972
2973 size_t __n;
2974 if (__is >> __n)
2975 {
2976 std::vector<_RealType> __int_vec;
2977 if (__detail::__extract_params(__is, __int_vec, __n + 1))
2978 {
2979 std::vector<double> __den_vec;
2980 if (__detail::__extract_params(__is, __den_vec, __n))
2981 {
2982 __x.param({ __int_vec.begin(), __int_vec.end(),
2983 __den_vec.begin() });
2984 }
2985 }
2986 }
2987
2988 __is.flags(__flags);
2989 return __is;
2990 }
2991
2992
2993 template<typename _RealType>
2994 void
2995 piecewise_linear_distribution<_RealType>::param_type::
2996 _M_initialize()
2997 {
2998 if (_M_int.size() < 2
2999 || (_M_int.size() == 2
3000 && _M_int[0] == _RealType(0)
3001 && _M_int[1] == _RealType(1)
3002 && _M_den[0] == _M_den[1]))
3003 {
3004 _M_int.clear();
3005 _M_den.clear();
3006 return;
3007 }
3008
3009 double __sum = 0.0;
3010 _M_cp.reserve(_M_int.size() - 1);
3011 _M_m.reserve(_M_int.size() - 1);
3012 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3013 {
3014 const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3015 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3016 _M_cp.push_back(__sum);
3017 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3018 }
3019
3020 // Now normalize the densities...
3021 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3022 __sum);
3023 // ... and partial sums...
3024 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3025 // ... and slopes.
3026 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3027
3028 // Make sure the last cumulative probablility is one.
3029 _M_cp[_M_cp.size() - 1] = 1.0;
3030 }
3031
3032 template<typename _RealType>
3033 template<typename _InputIteratorB, typename _InputIteratorW>
3034 piecewise_linear_distribution<_RealType>::param_type::
3035 param_type(_InputIteratorB __bbegin,
3036 _InputIteratorB __bend,
3037 _InputIteratorW __wbegin)
3038 : _M_int(), _M_den(), _M_cp(), _M_m()
3039 {
3040 for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
3041 {
3042 _M_int.push_back(*__bbegin);
3043 _M_den.push_back(*__wbegin);
3044 }
3045
3046 _M_initialize();
3047 }
3048
3049 template<typename _RealType>
3050 template<typename _Func>
3051 piecewise_linear_distribution<_RealType>::param_type::
3052 param_type(initializer_list<_RealType> __bl, _Func __fw)
3053 : _M_int(), _M_den(), _M_cp(), _M_m()
3054 {
3055 _M_int.reserve(__bl.size());
3056 _M_den.reserve(__bl.size());
3057 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3058 {
3059 _M_int.push_back(*__biter);
3060 _M_den.push_back(__fw(*__biter));
3061 }
3062
3063 _M_initialize();
3064 }
3065
3066 template<typename _RealType>
3067 template<typename _Func>
3068 piecewise_linear_distribution<_RealType>::param_type::
3069 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3070 : _M_int(), _M_den(), _M_cp(), _M_m()
3071 {
3072 const size_t __n = __nw == 0 ? 1 : __nw;
3073 const _RealType __delta = (__xmax - __xmin) / __n;
3074
3075 _M_int.reserve(__n + 1);
3076 _M_den.reserve(__n + 1);
3077 for (size_t __k = 0; __k <= __nw; ++__k)
3078 {
3079 _M_int.push_back(__xmin + __k * __delta);
3080 _M_den.push_back(__fw(_M_int[__k] + __delta));
3081 }
3082
3083 _M_initialize();
3084 }
3085
3086 template<typename _RealType>
3087 template<typename _UniformRandomNumberGenerator>
3088 typename piecewise_linear_distribution<_RealType>::result_type
3089 piecewise_linear_distribution<_RealType>::
3090 operator()(_UniformRandomNumberGenerator& __urng,
3091 const param_type& __param)
3092 {
3093 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3094 __aurng(__urng);
3095
3096 const double __p = __aurng();
3097 if (__param._M_cp.empty())
3098 return __p;
3099
3100 auto __pos = std::lower_bound(__param._M_cp.begin(),
3101 __param._M_cp.end(), __p);
3102 const size_t __i = __pos - __param._M_cp.begin();
3103
3104 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3105
3106 const double __a = 0.5 * __param._M_m[__i];
3107 const double __b = __param._M_den[__i];
3108 const double __cm = __p - __pref;
3109
3110 _RealType __x = __param._M_int[__i];
3111 if (__a == 0)
3112 __x += __cm / __b;
3113 else
3114 {
3115 const double __d = __b * __b + 4.0 * __a * __cm;
3116 __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3117 }
3118
3119 return __x;
3120 }
3121
3122 template<typename _RealType>
3123 template<typename _ForwardIterator,
3124 typename _UniformRandomNumberGenerator>
3125 void
3126 piecewise_linear_distribution<_RealType>::
3127 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3128 _UniformRandomNumberGenerator& __urng,
3129 const param_type& __param)
3130 {
3131 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3132 // We could duplicate everything from operator()...
3133 while (__f != __t)
3134 *__f++ = this->operator()(__urng, __param);
3135 }
3136
3137 template<typename _RealType, typename _CharT, typename _Traits>
3140 const piecewise_linear_distribution<_RealType>& __x)
3141 {
3142 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3143
3144 const typename __ios_base::fmtflags __flags = __os.flags();
3145 const _CharT __fill = __os.fill();
3146 const std::streamsize __precision = __os.precision();
3147 const _CharT __space = __os.widen(' ');
3149 __os.fill(__space);
3151
3152 std::vector<_RealType> __int = __x.intervals();
3153 __os << __int.size() - 1;
3154
3155 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3156 __os << __space << *__xit;
3157
3158 std::vector<double> __den = __x.densities();
3159 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3160 __os << __space << *__dit;
3161
3162 __os.flags(__flags);
3163 __os.fill(__fill);
3164 __os.precision(__precision);
3165 return __os;
3166 }
3167
3168 template<typename _RealType, typename _CharT, typename _Traits>
3171 piecewise_linear_distribution<_RealType>& __x)
3172 {
3173 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3174
3175 const typename __ios_base::fmtflags __flags = __is.flags();
3177
3178 size_t __n;
3179 if (__is >> __n)
3180 {
3181 vector<_RealType> __int_vec;
3182 if (__detail::__extract_params(__is, __int_vec, __n + 1))
3183 {
3184 vector<double> __den_vec;
3185 if (__detail::__extract_params(__is, __den_vec, __n + 1))
3186 {
3187 __x.param({ __int_vec.begin(), __int_vec.end(),
3188 __den_vec.begin() });
3189 }
3190 }
3191 }
3192 __is.flags(__flags);
3193 return __is;
3194 }
3195
3196
3197 template<typename _IntType>
3198 seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3199 {
3200 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3201 _M_v.push_back(__detail::__mod<result_type,
3202 __detail::_Shift<result_type, 32>::__value>(*__iter));
3203 }
3204
3205 template<typename _InputIterator>
3206 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3207 {
3208 for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3209 _M_v.push_back(__detail::__mod<result_type,
3210 __detail::_Shift<result_type, 32>::__value>(*__iter));
3211 }
3212
3213 template<typename _RandomAccessIterator>
3214 void
3215 seed_seq::generate(_RandomAccessIterator __begin,
3216 _RandomAccessIterator __end)
3217 {
3218 typedef typename iterator_traits<_RandomAccessIterator>::value_type
3219 _Type;
3220
3221 if (__begin == __end)
3222 return;
3223
3224 std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3225
3226 const size_t __n = __end - __begin;
3227 const size_t __s = _M_v.size();
3228 const size_t __t = (__n >= 623) ? 11
3229 : (__n >= 68) ? 7
3230 : (__n >= 39) ? 5
3231 : (__n >= 7) ? 3
3232 : (__n - 1) / 2;
3233 const size_t __p = (__n - __t) / 2;
3234 const size_t __q = __p + __t;
3235 const size_t __m = std::max(size_t(__s + 1), __n);
3236
3237 for (size_t __k = 0; __k < __m; ++__k)
3238 {
3239 _Type __arg = (__begin[__k % __n]
3240 ^ __begin[(__k + __p) % __n]
3241 ^ __begin[(__k - 1) % __n]);
3242 _Type __r1 = __arg ^ (__arg >> 27);
3243 __r1 = __detail::__mod<_Type,
3244 __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
3245 _Type __r2 = __r1;
3246 if (__k == 0)
3247 __r2 += __s;
3248 else if (__k <= __s)
3249 __r2 += __k % __n + _M_v[__k - 1];
3250 else
3251 __r2 += __k % __n;
3252 __r2 = __detail::__mod<_Type,
3253 __detail::_Shift<_Type, 32>::__value>(__r2);
3254 __begin[(__k + __p) % __n] += __r1;
3255 __begin[(__k + __q) % __n] += __r2;
3256 __begin[__k % __n] = __r2;
3257 }
3258
3259 for (size_t __k = __m; __k < __m + __n; ++__k)
3260 {
3261 _Type __arg = (__begin[__k % __n]
3262 + __begin[(__k + __p) % __n]
3263 + __begin[(__k - 1) % __n]);
3264 _Type __r3 = __arg ^ (__arg >> 27);
3265 __r3 = __detail::__mod<_Type,
3266 __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
3267 _Type __r4 = __r3 - __k % __n;
3268 __r4 = __detail::__mod<_Type,
3269 __detail::_Shift<_Type, 32>::__value>(__r4);
3270 __begin[(__k + __p) % __n] ^= __r3;
3271 __begin[(__k + __q) % __n] ^= __r4;
3272 __begin[__k % __n] = __r4;
3273 }
3274 }
3275
3276 template<typename _RealType, size_t __bits,
3277 typename _UniformRandomNumberGenerator>
3278 _RealType
3279 generate_canonical(_UniformRandomNumberGenerator& __urng)
3280 {
3282 "template argument must be a floating point type");
3283
3284 const size_t __b
3285 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3286 __bits);
3287 const long double __r = static_cast<long double>(__urng.max())
3288 - static_cast<long double>(__urng.min()) + 1.0L;
3289 const size_t __log2r = std::log(__r) / std::log(2.0L);
3290 const size_t __m = std::max<size_t>(1UL,
3291 (__b + __log2r - 1UL) / __log2r);
3292 _RealType __ret;
3293 _RealType __sum = _RealType(0);
3294 _RealType __tmp = _RealType(1);
3295 for (size_t __k = __m; __k != 0; --__k)
3296 {
3297 __sum += _RealType(__urng() - __urng.min()) * __tmp;
3298 __tmp *= __r;
3299 }
3300 __ret = __sum / __tmp;
3301 if (__builtin_expect(__ret >= _RealType(1), 0))
3302 {
3303#if _GLIBCXX_USE_C99_MATH_TR1
3304 __ret = std::nextafter(_RealType(1), _RealType(0));
3305#else
3306 __ret = _RealType(1)
3307 - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3308#endif
3309 }
3310 return __ret;
3311 }
3312
3313_GLIBCXX_END_NAMESPACE_VERSION
3314} // namespace
3315
3316#endif
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
Definition: complex:823
complex< _Tp > tan(const complex< _Tp > &)
Return complex tangent of z.
Definition: complex:959
_Tp abs(const complex< _Tp > &)
Return magnitude of z.
Definition: complex:629
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.
Definition: complex:796
complex< _Tp > pow(const complex< _Tp > &, int)
Return x to the y'th power.
Definition: complex:1018
complex< _Tp > sqrt(const complex< _Tp > &)
Return complex square root of z.
Definition: complex:932
constexpr const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:254
constexpr const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:230
_RealType generate_canonical(_UniformRandomNumberGenerator &__g)
A function template for converting the output of a (integral) uniform random number generator to a fl...
basic_ostream< _Ch_type, _Ch_traits > & operator<<(basic_ostream< _Ch_type, _Ch_traits > &__os, const sub_match< _Bi_iter > &__m)
Inserts a matched string into an output stream.
Definition: regex.h:1563
constexpr back_insert_iterator< _Container > back_inserter(_Container &__x)
constexpr _Tp accumulate(_InputIterator __first, _InputIterator __last, _Tp __init)
Accumulate values in a range.
Definition: stl_numeric.h:134
constexpr _OutputIterator partial_sum(_InputIterator __first, _InputIterator __last, _OutputIterator __result)
Return list of partial sums.
Definition: stl_numeric.h:256
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4930
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
const _RandomNumberEngine & base() const noexcept
Definition: random.h:1422
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:1819
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4079
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4522
void seed(result_type __sd=default_seed)
Seeds the initial state of the random number generator.
static constexpr result_type multiplier
Definition: random.h:262
result_type operator()()
Gets the next value in the generated random number sequence.
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4049
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2521
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4715
result_type operator()()
Gets the next random number in the sequence.
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4960
_RandomNumberEngine::result_type result_type
Definition: random.h:1322
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:3856
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:5170
friend bool operator==(const poisson_distribution &__d1, const poisson_distribution &__d2)
Return true if two Poisson distributions have the same parameters and the sequences that would be gen...
Definition: random.h:4558
static constexpr result_type modulus
Definition: random.h:266
void seed(result_type __s=default_seed)
Reseeds the linear_congruential_engine random number generator engine sequence to the seed __s.
_RealType result_type
Definition: random.h:2400
result_type operator()()
Gets the next value in the generated random number sequence.
friend std::basic_ostream< _CharT, _Traits > & operator<<(std::basic_ostream< _CharT, _Traits > &__os, const std::poisson_distribution< _IntType1 > &__x)
Inserts a poisson_distribution random number distribution __x into the output stream __os.
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:2921
static constexpr result_type increment
Definition: random.h:264
friend std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, std::poisson_distribution< _IntType1 > &__x)
Extracts a poisson_distribution random number distribution __x from the input stream __is.
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2951
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2079
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:5140
ISO C++ entities toplevel namespace is std.
ios_base & scientific(ios_base &__base)
Calls base.setf(ios_base::scientific, ios_base::floatfield).
Definition: ios_base.h:1056
ptrdiff_t streamsize
Integral type for I/O operation counts and buffer sizes.
Definition: postypes.h:98
ios_base & left(ios_base &__base)
Calls base.setf(ios_base::left, ios_base::adjustfield).
Definition: ios_base.h:1006
ios_base & dec(ios_base &__base)
Calls base.setf(ios_base::dec, ios_base::basefield).
Definition: ios_base.h:1023
ios_base & skipws(ios_base &__base)
Calls base.setf(ios_base::skipws).
Definition: ios_base.h:949
ios_base & fixed(ios_base &__base)
Calls base.setf(ios_base::fixed, ios_base::floatfield).
Definition: ios_base.h:1048
constexpr int __lg(int __n)
This is a helper function for the sort routines and for random.tcc.
std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, bitset< _Nb > &__x)
Global I/O operators for bitsets.
Definition: bitset:1470
initializer_list
void clear(iostate __state=goodbit)
[Re]sets the error state.
Definition: basic_ios.tcc:41
char_type widen(char __c) const
Widens characters.
Definition: basic_ios.h:449
char_type fill() const
Retrieves the empty character.
Definition: basic_ios.h:370
Template class basic_istream.
Definition: istream:59
Template class basic_ostream.
Definition: ostream:59
Properties of fundamental types.
Definition: limits:313
static constexpr _Tp max() noexcept
Definition: limits:321
static constexpr _Tp epsilon() noexcept
Definition: limits:333
static constexpr _Tp min() noexcept
Definition: limits:317
is_floating_point
Definition: type_traits:395
common_type
Definition: type_traits:2202
streamsize precision() const
Flags access.
Definition: ios_base.h:696
fmtflags flags() const
Access to format flags.
Definition: ios_base.h:626
A model of a linear congruential random number generator.
Definition: random.h:247
The Marsaglia-Zaman generator.
Definition: random.h:684
Produces random numbers by combining random numbers from some base engine to produce random numbers w...
Definition: random.h:1316
Uniform continuous distribution for random numbers.
Definition: random.h:1732
A normal continuous distribution for random numbers.
Definition: random.h:1962
A gamma continuous distribution for random numbers.
Definition: random.h:2394
A chi_squared_distribution random number distribution.
Definition: random.h:2622
A cauchy_distribution random number distribution.
Definition: random.h:2846
A fisher_f_distribution random number distribution.
Definition: random.h:3054
A student_t_distribution random number distribution.
Definition: random.h:3286
A discrete binomial random number distribution.
Definition: random.h:3730
A discrete geometric random number distribution.
Definition: random.h:3970
A discrete Poisson random number distribution.
Definition: random.h:4411
An exponential continuous distribution for random numbers.
Definition: random.h:4637
A weibull_distribution random number distribution.
Definition: random.h:4852
A extreme_value_distribution random number distribution.
Definition: random.h:5062
iterator begin() noexcept
Definition: stl_vector.h:808
iterator end() noexcept
Definition: stl_vector.h:826
size_type size() const noexcept
Definition: stl_vector.h:915
Uniform discrete distribution for random numbers. A discrete random distribution on the range with e...
param_type param() const
Returns the parameter set of the distribution.
Parallel STL function calls corresponding to stl_numeric.h. The functions defined here mainly do case...