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