libstdc++
bits/random.tcc
Go to the documentation of this file.
1// random number generation (out of line) -*- C++ -*-
2
3// Copyright (C) 2009-2026 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, (void) ++__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 // _GLIBCXX_RESOLVE_LIB_DEFECTS
545 // 3809. Is std::subtract_with_carry_engine<uint16_t> supposed to work?
546 // 4014. LWG 3809 changes behavior of some existing code
547 std::linear_congruential_engine<uint_least32_t, 40014u, 0u, 2147483563u>
548 __lcg(__value == 0u ? default_seed : __value % 2147483563u);
549
550 const size_t __n = (__w + 31) / 32;
551
552 for (size_t __i = 0; __i < long_lag; ++__i)
553 {
554 _UIntType __sum = 0u;
555 _UIntType __factor = 1u;
556 for (size_t __j = 0; __j < __n; ++__j)
557 {
558 __sum += __detail::__mod<uint_least32_t,
559 __detail::_Shift<uint_least32_t, 32>::__value>
560 (__lcg()) * __factor;
561 __factor *= __detail::_Shift<_UIntType, 32>::__value;
562 }
563 _M_x[__i] = __detail::__mod<_UIntType,
564 __detail::_Shift<_UIntType, __w>::__value>(__sum);
565 }
566 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
567 _M_p = 0;
568 }
569
570 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
571 template<typename _Sseq>
572 auto
573 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
574 seed(_Sseq& __q)
575 -> _If_seed_seq<_Sseq>
576 {
577 const size_t __k = (__w + 31) / 32;
578 uint_least32_t __arr[__r * __k];
579 __q.generate(__arr + 0, __arr + __r * __k);
580
581 for (size_t __i = 0; __i < long_lag; ++__i)
582 {
583 _UIntType __sum = 0u;
584 _UIntType __factor = 1u;
585 for (size_t __j = 0; __j < __k; ++__j)
586 {
587 __sum += __arr[__k * __i + __j] * __factor;
588 __factor *= __detail::_Shift<_UIntType, 32>::__value;
589 }
590 _M_x[__i] = __detail::__mod<_UIntType,
591 __detail::_Shift<_UIntType, __w>::__value>(__sum);
592 }
593 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
594 _M_p = 0;
595 }
596
597 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
598 typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
599 result_type
600 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
601 operator()()
602 {
603 // Derive short lag index from current index.
604 long __ps = _M_p - short_lag;
605 if (__ps < 0)
606 __ps += long_lag;
607
608 // Calculate new x(i) without overflow or division.
609 // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
610 // cannot overflow.
611 _UIntType __xi;
612 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
613 {
614 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
615 _M_carry = 0;
616 }
617 else
618 {
619 __xi = (__detail::_Shift<_UIntType, __w>::__value
620 - _M_x[_M_p] - _M_carry + _M_x[__ps]);
621 _M_carry = 1;
622 }
623 _M_x[_M_p] = __xi;
624
625 // Adjust current index to loop around in ring buffer.
626 if (++_M_p >= long_lag)
627 _M_p = 0;
628
629 return __xi;
630 }
631
632 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
633 typename _CharT, typename _Traits>
634 std::basic_ostream<_CharT, _Traits>&
635 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
636 const subtract_with_carry_engine<_UIntType,
637 __w, __s, __r>& __x)
638 {
639 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
640
641 const typename __ios_base::fmtflags __flags = __os.flags();
642 const _CharT __fill = __os.fill();
643 const _CharT __space = __os.widen(' ');
644 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
645 __os.fill(__space);
646
647 for (size_t __i = 0; __i < __r; ++__i)
648 __os << __x._M_x[__i] << __space;
649 __os << __x._M_carry << __space << __x._M_p;
650
651 __os.flags(__flags);
652 __os.fill(__fill);
653 return __os;
654 }
655
656 template<typename _UIntType, size_t __w, size_t __s, size_t __r,
657 typename _CharT, typename _Traits>
658 std::basic_istream<_CharT, _Traits>&
659 operator>>(std::basic_istream<_CharT, _Traits>& __is,
660 subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
661 {
662 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
663
664 const typename __ios_base::fmtflags __flags = __is.flags();
665 __is.flags(__ios_base::dec | __ios_base::skipws);
666
667 for (size_t __i = 0; __i < __r; ++__i)
668 __is >> __x._M_x[__i];
669 __is >> __x._M_carry;
670 __is >> __x._M_p;
671
672 __is.flags(__flags);
673 return __is;
674 }
675
676#if ! __cpp_inline_variables
677 template<typename _RandomNumberEngine, size_t __p, size_t __r>
678 constexpr size_t
679 discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
680
681 template<typename _RandomNumberEngine, size_t __p, size_t __r>
682 constexpr size_t
683 discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
684#endif
685
686 template<typename _RandomNumberEngine, size_t __p, size_t __r>
687 typename discard_block_engine<_RandomNumberEngine,
688 __p, __r>::result_type
689 discard_block_engine<_RandomNumberEngine, __p, __r>::
690 operator()()
691 {
692 if (_M_n >= used_block)
693 {
694 _M_b.discard(block_size - _M_n);
695 _M_n = 0;
696 }
697 ++_M_n;
698 return _M_b();
699 }
700
701 template<typename _RandomNumberEngine, size_t __p, size_t __r,
702 typename _CharT, typename _Traits>
703 std::basic_ostream<_CharT, _Traits>&
704 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
705 const discard_block_engine<_RandomNumberEngine,
706 __p, __r>& __x)
707 {
708 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
709
710 const typename __ios_base::fmtflags __flags = __os.flags();
711 const _CharT __fill = __os.fill();
712 const _CharT __space = __os.widen(' ');
713 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
714 __os.fill(__space);
715
716 __os << __x.base() << __space << __x._M_n;
717
718 __os.flags(__flags);
719 __os.fill(__fill);
720 return __os;
721 }
722
723 template<typename _RandomNumberEngine, size_t __p, size_t __r,
724 typename _CharT, typename _Traits>
725 std::basic_istream<_CharT, _Traits>&
726 operator>>(std::basic_istream<_CharT, _Traits>& __is,
727 discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
728 {
729 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
730
731 const typename __ios_base::fmtflags __flags = __is.flags();
732 __is.flags(__ios_base::dec | __ios_base::skipws);
733
734 __is >> __x._M_b >> __x._M_n;
735
736 __is.flags(__flags);
737 return __is;
738 }
739
740
741 template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
742 typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
743 result_type
744 independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
745 operator()()
746 {
747 typedef typename _RandomNumberEngine::result_type _Eresult_type;
748 const _Eresult_type __r
749 = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
750 ? _M_b.max() - _M_b.min() + 1 : 0);
751 const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
752 const unsigned __m = __r ? std::__lg(__r) : __edig;
753
754 typedef typename std::common_type<_Eresult_type, result_type>::type
755 __ctype;
756 const unsigned __cdig = std::numeric_limits<__ctype>::digits;
757
758 unsigned __n, __n0;
759 __ctype __s0, __s1, __y0, __y1;
760
761 for (size_t __i = 0; __i < 2; ++__i)
762 {
763 __n = (__w + __m - 1) / __m + __i;
764 __n0 = __n - __w % __n;
765 const unsigned __w0 = __w / __n; // __w0 <= __m
766
767 __s0 = 0;
768 __s1 = 0;
769 if (__w0 < __cdig)
770 {
771 __s0 = __ctype(1) << __w0;
772 __s1 = __s0 << 1;
773 }
774
775 __y0 = 0;
776 __y1 = 0;
777 if (__r)
778 {
779 __y0 = __s0 * (__r / __s0);
780 if (__s1)
781 __y1 = __s1 * (__r / __s1);
782
783 if (__r - __y0 <= __y0 / __n)
784 break;
785 }
786 else
787 break;
788 }
789
790 result_type __sum = 0;
791 for (size_t __k = 0; __k < __n0; ++__k)
792 {
793 __ctype __u;
794 do
795 __u = _M_b() - _M_b.min();
796 while (__y0 && __u >= __y0);
797 __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
798 }
799 for (size_t __k = __n0; __k < __n; ++__k)
800 {
801 __ctype __u;
802 do
803 __u = _M_b() - _M_b.min();
804 while (__y1 && __u >= __y1);
805 __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
806 }
807 return __sum;
808 }
809
810#if ! __cpp_inline_variables
811 template<typename _RandomNumberEngine, size_t __k>
812 constexpr size_t
813 shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
814#endif
815
816 namespace __detail
817 {
818 // Determine whether an integer is representable as double.
819 template<typename _Tp>
820 constexpr bool
821 __representable_as_double(_Tp __x) noexcept
822 {
823 static_assert(numeric_limits<_Tp>::is_integer, "");
824 static_assert(!numeric_limits<_Tp>::is_signed, "");
825 // All integers <= 2^53 are representable.
826 return (__x <= (1ull << __DBL_MANT_DIG__))
827 // Between 2^53 and 2^54 only even numbers are representable.
828 || (!(__x & 1) && __detail::__representable_as_double(__x >> 1));
829 }
830
831 // Determine whether x+1 is representable as double.
832 template<typename _Tp>
833 constexpr bool
834 __p1_representable_as_double(_Tp __x) noexcept
835 {
836 static_assert(numeric_limits<_Tp>::is_integer, "");
837 static_assert(!numeric_limits<_Tp>::is_signed, "");
838 return numeric_limits<_Tp>::digits < __DBL_MANT_DIG__
839 || (bool(__x + 1u) // return false if x+1 wraps around to zero
840 && __detail::__representable_as_double(__x + 1u));
841 }
842 }
843
844 template<typename _RandomNumberEngine, size_t __k>
845 typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
846 shuffle_order_engine<_RandomNumberEngine, __k>::
847 operator()()
848 {
849 constexpr result_type __range = max() - min();
850 size_t __j = __k;
851 const result_type __y = _M_y - min();
852#pragma GCC diagnostic push
853#pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
854 // Avoid using slower long double arithmetic if possible.
855 if constexpr (__detail::__p1_representable_as_double(__range))
856 __j *= __y / (__range + 1.0);
857 else
858 __j *= __y / (__range + 1.0L);
859#pragma GCC diagnostic pop
860 _M_y = _M_v[__j];
861 _M_v[__j] = _M_b();
862
863 return _M_y;
864 }
865
866 template<typename _RandomNumberEngine, size_t __k,
867 typename _CharT, typename _Traits>
868 std::basic_ostream<_CharT, _Traits>&
869 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
870 const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
871 {
872 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
873
874 const typename __ios_base::fmtflags __flags = __os.flags();
875 const _CharT __fill = __os.fill();
876 const _CharT __space = __os.widen(' ');
877 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
878 __os.fill(__space);
879
880 __os << __x.base();
881 for (size_t __i = 0; __i < __k; ++__i)
882 __os << __space << __x._M_v[__i];
883 __os << __space << __x._M_y;
884
885 __os.flags(__flags);
886 __os.fill(__fill);
887 return __os;
888 }
889
890 template<typename _RandomNumberEngine, size_t __k,
891 typename _CharT, typename _Traits>
892 std::basic_istream<_CharT, _Traits>&
893 operator>>(std::basic_istream<_CharT, _Traits>& __is,
894 shuffle_order_engine<_RandomNumberEngine, __k>& __x)
895 {
896 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
897
898 const typename __ios_base::fmtflags __flags = __is.flags();
899 __is.flags(__ios_base::dec | __ios_base::skipws);
900
901 __is >> __x._M_b;
902 for (size_t __i = 0; __i < __k; ++__i)
903 __is >> __x._M_v[__i];
904 __is >> __x._M_y;
905
906 __is.flags(__flags);
907 return __is;
908 }
909
910#if __glibcxx_philox_engine // >= C++26
911
912 template<typename _UIntType, size_t __w, size_t __n, size_t __r,
913 _UIntType... __consts>
914 _UIntType
915 philox_engine<_UIntType, __w, __n, __r, __consts...>::
916 _S_mulhi(_UIntType __a, _UIntType __b)
917 {
918 using __type = typename __detail::_Select_uint_least_t<__w * 2>::type;
919 const __type __num = static_cast<__type>(__a) * __b;
920 return static_cast<_UIntType>(__num >> __w) & max();
921 }
922
923 template<typename _UIntType, size_t __w, size_t __n, size_t __r,
924 _UIntType... __consts>
925 _UIntType
926 philox_engine<_UIntType, __w, __n, __r, __consts...>::
927 _S_mullo(_UIntType __a, _UIntType __b)
928 {
929 return static_cast<_UIntType>((__a * __b) & max());
930 }
931
932 template<typename _UIntType, size_t __w, size_t __n, size_t __r,
933 _UIntType... __consts>
934 void
935 philox_engine<_UIntType, __w, __n, __r, __consts...>::_M_transition()
936 {
937 ++_M_i;
938 if (_M_i != __n)
939 return;
940
941 using __type = typename __detail::_Select_uint_least_t<__w * 2>::type;
942
943 _M_philox();
944 if constexpr (__n == 4)
945 {
946 __type __uh
947 = (static_cast<__type>(_M_x[1]) << __w)
948 | static_cast<__type>(_M_x[0]);
949 ++__uh;
950 __type __lh
951 = (static_cast<__type>(_M_x[3]) << __w)
952 | static_cast<__type>(_M_x[2]);
953 __type __bigMask
954 = ~__type(0) >> ((sizeof(__type) * __CHAR_BIT__) - (__w * 2));
955 if ((__uh & __bigMask) == 0)
956 {
957 ++__lh;
958 __uh = 0;
959 }
960 _M_x[0] = static_cast<_UIntType>(__uh & max());
961 _M_x[1] = static_cast<_UIntType>((__uh >> (__w)) & max());
962 _M_x[2] = static_cast<_UIntType>(__lh & max());
963 _M_x[3] = static_cast<_UIntType>((__lh >> (__w)) & max());
964 }
965 else
966 {
967 __type __num
968 = (static_cast<__type>(_M_x[1]) << __w)
969 | static_cast<__type>(_M_x[0]);
970 ++__num;
971 _M_x[0] = static_cast<_UIntType>(__num & max());
972 _M_x[1] = static_cast<_UIntType>((__num >> __w) & max());
973 }
974 _M_i = 0;
975 }
976
977 template<typename _UIntType, size_t __w, size_t __n, size_t __r,
978 _UIntType... __consts>
979 void
980 philox_engine<_UIntType, __w, __n, __r, __consts...>::_M_philox()
981 {
982 array<_UIntType, __n> __outputSeq = _M_x;
983 for (size_t __j = 0; __j < __r; ++__j)
984 {
985 array<_UIntType, __n> __intermedSeq{};
986 if constexpr (__n == 4)
987 {
988 __intermedSeq[0] = __outputSeq[2];
989 __intermedSeq[1] = __outputSeq[1];
990 __intermedSeq[2] = __outputSeq[0];
991 __intermedSeq[3] = __outputSeq[3];
992 }
993 else
994 {
995 __intermedSeq[0] = __outputSeq[0];
996 __intermedSeq[1] = __outputSeq[1];
997 }
998 for (unsigned long __k = 0; __k < (__n/2); ++__k)
999 {
1000 __outputSeq[2*__k]
1001 = _S_mulhi(__intermedSeq[2*__k], multipliers[__k])
1002 ^ (((_M_k[__k] + (__j * round_consts[__k])) & max()))
1003 ^ __intermedSeq[2*__k+1];
1004
1005 __outputSeq[(2*__k)+1]
1006 = _S_mullo(__intermedSeq[2*__k], multipliers[__k]);
1007 }
1008 }
1009 _M_y = __outputSeq;
1010 }
1011
1012 template<typename _UIntType, size_t __w, size_t __n, size_t __r,
1013 _UIntType... __consts>
1014 template<typename _Sseq>
1015 void
1016 philox_engine<_UIntType, __w, __n, __r, __consts...>::seed(_Sseq& __q)
1017 requires __is_seed_seq<_Sseq>
1018 {
1019 seed(0);
1020
1021 const unsigned __p = 1 + ((__w - 1) / 32);
1022 uint_least32_t __tmpArr[(__n/2) * __p];
1023 __q.generate(__tmpArr + 0, __tmpArr + ((__n/2) * __p));
1024 for (unsigned __k = 0; __k < (__n/2); ++__k)
1025 {
1026 unsigned long long __precalc = 0;
1027 for (unsigned __j = 0; __j < __p; ++__j)
1028 {
1029 unsigned long long __multiplicand = (1ull << (32 * __j));
1030 __precalc += (__tmpArr[__k * __p + __j] * __multiplicand) & max();
1031 }
1032 _M_k[__k] = __precalc;
1033 }
1034 }
1035#endif // philox_engine
1036
1037 template<typename _IntType, typename _CharT, typename _Traits>
1038 std::basic_ostream<_CharT, _Traits>&
1039 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1040 const uniform_int_distribution<_IntType>& __x)
1041 {
1042 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1043
1044 const typename __ios_base::fmtflags __flags = __os.flags();
1045 const _CharT __fill = __os.fill();
1046 const _CharT __space = __os.widen(' ');
1047 __os.flags(__ios_base::scientific | __ios_base::left);
1048 __os.fill(__space);
1049
1050 __os << __x.a() << __space << __x.b();
1051
1052 __os.flags(__flags);
1053 __os.fill(__fill);
1054 return __os;
1055 }
1056
1057 template<typename _IntType, typename _CharT, typename _Traits>
1058 std::basic_istream<_CharT, _Traits>&
1059 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1060 uniform_int_distribution<_IntType>& __x)
1061 {
1062 using param_type
1063 = typename uniform_int_distribution<_IntType>::param_type;
1064 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1065
1066 const typename __ios_base::fmtflags __flags = __is.flags();
1067 __is.flags(__ios_base::dec | __ios_base::skipws);
1068
1069 _IntType __a, __b;
1070 if (__is >> __a >> __b)
1071 __x.param(param_type(__a, __b));
1072
1073 __is.flags(__flags);
1074 return __is;
1075 }
1076
1077
1078 template<typename _RealType>
1079 template<typename _ForwardIterator,
1080 typename _UniformRandomNumberGenerator>
1081 void
1082 uniform_real_distribution<_RealType>::
1083 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1084 _UniformRandomNumberGenerator& __urng,
1085 const param_type& __p)
1086 {
1087 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1088 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1089 __aurng(__urng);
1090 auto __range = __p.b() - __p.a();
1091 while (__f != __t)
1092 *__f++ = __aurng() * __range + __p.a();
1093 }
1094
1095 template<typename _RealType, typename _CharT, typename _Traits>
1096 std::basic_ostream<_CharT, _Traits>&
1097 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1098 const uniform_real_distribution<_RealType>& __x)
1099 {
1100 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1101
1102 const typename __ios_base::fmtflags __flags = __os.flags();
1103 const _CharT __fill = __os.fill();
1104 const std::streamsize __precision = __os.precision();
1105 const _CharT __space = __os.widen(' ');
1106 __os.flags(__ios_base::scientific | __ios_base::left);
1107 __os.fill(__space);
1108 __os.precision(std::numeric_limits<_RealType>::max_digits10);
1109
1110 __os << __x.a() << __space << __x.b();
1111
1112 __os.flags(__flags);
1113 __os.fill(__fill);
1114 __os.precision(__precision);
1115 return __os;
1116 }
1117
1118 template<typename _RealType, typename _CharT, typename _Traits>
1119 std::basic_istream<_CharT, _Traits>&
1120 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1121 uniform_real_distribution<_RealType>& __x)
1122 {
1123 using param_type
1124 = typename uniform_real_distribution<_RealType>::param_type;
1125 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1126
1127 const typename __ios_base::fmtflags __flags = __is.flags();
1128 __is.flags(__ios_base::skipws);
1129
1130 _RealType __a, __b;
1131 if (__is >> __a >> __b)
1132 __x.param(param_type(__a, __b));
1133
1134 __is.flags(__flags);
1135 return __is;
1136 }
1137
1138
1139 template<typename _ForwardIterator,
1140 typename _UniformRandomNumberGenerator>
1141 void
1142 std::bernoulli_distribution::
1143 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1144 _UniformRandomNumberGenerator& __urng,
1145 const param_type& __p)
1146 {
1147 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1148 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1149 __aurng(__urng);
1150 auto __limit = __p.p() * (__aurng.max() - __aurng.min());
1151
1152 while (__f != __t)
1153 *__f++ = (__aurng() - __aurng.min()) < __limit;
1154 }
1155
1156 template<typename _CharT, typename _Traits>
1157 std::basic_ostream<_CharT, _Traits>&
1158 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1159 const bernoulli_distribution& __x)
1160 {
1161 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1162
1163 const typename __ios_base::fmtflags __flags = __os.flags();
1164 const _CharT __fill = __os.fill();
1165 const std::streamsize __precision = __os.precision();
1166 __os.flags(__ios_base::scientific | __ios_base::left);
1167 __os.fill(__os.widen(' '));
1168 __os.precision(std::numeric_limits<double>::max_digits10);
1169
1170 __os << __x.p();
1171
1172 __os.flags(__flags);
1173 __os.fill(__fill);
1174 __os.precision(__precision);
1175 return __os;
1176 }
1177
1178
1179 template<typename _IntType>
1180 template<typename _UniformRandomNumberGenerator>
1181 typename geometric_distribution<_IntType>::result_type
1182 geometric_distribution<_IntType>::
1183 operator()(_UniformRandomNumberGenerator& __urng,
1184 const param_type& __param)
1185 {
1186 // About the epsilon thing see this thread:
1187 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1188 const double __naf =
1189 (1 - std::numeric_limits<double>::epsilon()) / 2;
1190 // The largest _RealType convertible to _IntType.
1191 const double __thr =
1192 std::numeric_limits<_IntType>::max() + __naf;
1193 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1194 __aurng(__urng);
1195
1196 double __cand;
1197 do
1198 __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1199 while (__cand >= __thr);
1200
1201 return result_type(__cand + __naf);
1202 }
1203
1204 template<typename _IntType>
1205 template<typename _ForwardIterator,
1206 typename _UniformRandomNumberGenerator>
1207 void
1208 geometric_distribution<_IntType>::
1209 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1210 _UniformRandomNumberGenerator& __urng,
1211 const param_type& __param)
1212 {
1213 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1214 // About the epsilon thing see this thread:
1215 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1216 const double __naf =
1217 (1 - std::numeric_limits<double>::epsilon()) / 2;
1218 // The largest _RealType convertible to _IntType.
1219 const double __thr =
1220 std::numeric_limits<_IntType>::max() + __naf;
1221 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1222 __aurng(__urng);
1223
1224 while (__f != __t)
1225 {
1226 double __cand;
1227 do
1228 __cand = std::floor(std::log(1.0 - __aurng())
1229 / __param._M_log_1_p);
1230 while (__cand >= __thr);
1231
1232 *__f++ = __cand + __naf;
1233 }
1234 }
1235
1236 template<typename _IntType,
1237 typename _CharT, typename _Traits>
1238 std::basic_ostream<_CharT, _Traits>&
1239 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1240 const geometric_distribution<_IntType>& __x)
1241 {
1242 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1243
1244 const typename __ios_base::fmtflags __flags = __os.flags();
1245 const _CharT __fill = __os.fill();
1246 const std::streamsize __precision = __os.precision();
1247 __os.flags(__ios_base::scientific | __ios_base::left);
1248 __os.fill(__os.widen(' '));
1249 __os.precision(std::numeric_limits<double>::max_digits10);
1250
1251 __os << __x.p();
1252
1253 __os.flags(__flags);
1254 __os.fill(__fill);
1255 __os.precision(__precision);
1256 return __os;
1257 }
1258
1259 template<typename _IntType,
1260 typename _CharT, typename _Traits>
1261 std::basic_istream<_CharT, _Traits>&
1262 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1263 geometric_distribution<_IntType>& __x)
1264 {
1265 using param_type = typename geometric_distribution<_IntType>::param_type;
1266 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1267
1268 const typename __ios_base::fmtflags __flags = __is.flags();
1269 __is.flags(__ios_base::skipws);
1270
1271 double __p;
1272 if (__is >> __p)
1273 __x.param(param_type(__p));
1274
1275 __is.flags(__flags);
1276 return __is;
1277 }
1278
1279 // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1280 template<typename _IntType>
1281 template<typename _UniformRandomNumberGenerator>
1282 typename negative_binomial_distribution<_IntType>::result_type
1283 negative_binomial_distribution<_IntType>::
1284 operator()(_UniformRandomNumberGenerator& __urng)
1285 {
1286 const double __y = _M_gd(__urng);
1287
1288 // XXX Is the constructor too slow?
1289 std::poisson_distribution<result_type> __poisson(__y);
1290 return __poisson(__urng);
1291 }
1292
1293 template<typename _IntType>
1294 template<typename _UniformRandomNumberGenerator>
1295 typename negative_binomial_distribution<_IntType>::result_type
1296 negative_binomial_distribution<_IntType>::
1297 operator()(_UniformRandomNumberGenerator& __urng,
1298 const param_type& __p)
1299 {
1300 typedef typename std::gamma_distribution<double>::param_type
1301 param_type;
1302
1303 const double __y =
1304 _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1305
1306 std::poisson_distribution<result_type> __poisson(__y);
1307 return __poisson(__urng);
1308 }
1309
1310 template<typename _IntType>
1311 template<typename _ForwardIterator,
1312 typename _UniformRandomNumberGenerator>
1313 void
1314 negative_binomial_distribution<_IntType>::
1315 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1316 _UniformRandomNumberGenerator& __urng)
1317 {
1318 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1319 while (__f != __t)
1320 {
1321 const double __y = _M_gd(__urng);
1322
1323 // XXX Is the constructor too slow?
1324 std::poisson_distribution<result_type> __poisson(__y);
1325 *__f++ = __poisson(__urng);
1326 }
1327 }
1328
1329 template<typename _IntType>
1330 template<typename _ForwardIterator,
1331 typename _UniformRandomNumberGenerator>
1332 void
1333 negative_binomial_distribution<_IntType>::
1334 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1335 _UniformRandomNumberGenerator& __urng,
1336 const param_type& __p)
1337 {
1338 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1339 typename std::gamma_distribution<result_type>::param_type
1340 __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1341
1342 while (__f != __t)
1343 {
1344 const double __y = _M_gd(__urng, __p2);
1345
1346 std::poisson_distribution<result_type> __poisson(__y);
1347 *__f++ = __poisson(__urng);
1348 }
1349 }
1350
1351 template<typename _IntType, typename _CharT, typename _Traits>
1352 std::basic_ostream<_CharT, _Traits>&
1353 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1354 const negative_binomial_distribution<_IntType>& __x)
1355 {
1356 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1357
1358 const typename __ios_base::fmtflags __flags = __os.flags();
1359 const _CharT __fill = __os.fill();
1360 const std::streamsize __precision = __os.precision();
1361 const _CharT __space = __os.widen(' ');
1362 __os.flags(__ios_base::scientific | __ios_base::left);
1363 __os.fill(__os.widen(' '));
1364 __os.precision(std::numeric_limits<double>::max_digits10);
1365
1366 __os << __x.k() << __space << __x.p()
1367 << __space << __x._M_gd;
1368
1369 __os.flags(__flags);
1370 __os.fill(__fill);
1371 __os.precision(__precision);
1372 return __os;
1373 }
1374
1375 template<typename _IntType, typename _CharT, typename _Traits>
1376 std::basic_istream<_CharT, _Traits>&
1377 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1378 negative_binomial_distribution<_IntType>& __x)
1379 {
1380 using param_type
1381 = typename negative_binomial_distribution<_IntType>::param_type;
1382 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1383
1384 const typename __ios_base::fmtflags __flags = __is.flags();
1385 __is.flags(__ios_base::skipws);
1386
1387 _IntType __k;
1388 double __p;
1389 if (__is >> __k >> __p >> __x._M_gd)
1390 __x.param(param_type(__k, __p));
1391
1392 __is.flags(__flags);
1393 return __is;
1394 }
1395
1396
1397 template<typename _IntType>
1398 void
1399 poisson_distribution<_IntType>::param_type::
1400 _M_initialize()
1401 {
1402#if _GLIBCXX_USE_C99_MATH_FUNCS
1403 if (_M_mean >= 12)
1404 {
1405 const double __m = std::floor(_M_mean);
1406 _M_lm_thr = std::log(_M_mean);
1407 _M_lfm = std::lgamma(__m + 1);
1408 _M_sm = std::sqrt(__m);
1409
1410 const double __pi_4 = 0.7853981633974483096156608458198757L;
1411 const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1412 / __pi_4));
1413 _M_d = std::round(std::max<double>(6.0, std::min(__m, __dx)));
1414 const double __cx = 2 * __m + _M_d;
1415 _M_scx = std::sqrt(__cx / 2);
1416 _M_1cx = 1 / __cx;
1417
1418 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1419 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1420 / _M_d;
1421 }
1422 else
1423#endif
1424 _M_lm_thr = std::exp(-_M_mean);
1425 }
1426
1427 /**
1428 * A rejection algorithm when mean >= 12 and a simple method based
1429 * upon the multiplication of uniform random variates otherwise.
1430 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_FUNCS
1431 * is defined.
1432 *
1433 * Reference:
1434 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1435 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1436 */
1437 template<typename _IntType>
1438 template<typename _UniformRandomNumberGenerator>
1439 typename poisson_distribution<_IntType>::result_type
1440 poisson_distribution<_IntType>::
1441 operator()(_UniformRandomNumberGenerator& __urng,
1442 const param_type& __param)
1443 {
1444 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1445 __aurng(__urng);
1446#if _GLIBCXX_USE_C99_MATH_FUNCS
1447 if (__param.mean() >= 12)
1448 {
1449 double __x;
1450
1451 // See comments above...
1452 const double __naf =
1453 (1 - std::numeric_limits<double>::epsilon()) / 2;
1454 const double __thr =
1455 std::numeric_limits<_IntType>::max() + __naf;
1456
1457 const double __m = std::floor(__param.mean());
1458 // sqrt(pi / 2)
1459 const double __spi_2 = 1.2533141373155002512078826424055226L;
1460 const double __c1 = __param._M_sm * __spi_2;
1461 const double __c2 = __param._M_c2b + __c1;
1462 const double __c3 = __c2 + 1;
1463 const double __c4 = __c3 + 1;
1464 // 1 / 78
1465 const double __178 = 0.0128205128205128205128205128205128L;
1466 // e^(1 / 78)
1467 const double __e178 = 1.0129030479320018583185514777512983L;
1468 const double __c5 = __c4 + __e178;
1469 const double __c = __param._M_cb + __c5;
1470 const double __2cx = 2 * (2 * __m + __param._M_d);
1471
1472 bool __reject = true;
1473 do
1474 {
1475 const double __u = __c * __aurng();
1476 const double __e = -std::log(1.0 - __aurng());
1477
1478 double __w = 0.0;
1479
1480 if (__u <= __c1)
1481 {
1482 const double __n = _M_nd(__urng);
1483 const double __y = -std::abs(__n) * __param._M_sm - 1;
1484 __x = std::floor(__y);
1485 __w = -__n * __n / 2;
1486 if (__x < -__m)
1487 continue;
1488 }
1489 else if (__u <= __c2)
1490 {
1491 const double __n = _M_nd(__urng);
1492 const double __y = 1 + std::abs(__n) * __param._M_scx;
1493 __x = std::ceil(__y);
1494 __w = __y * (2 - __y) * __param._M_1cx;
1495 if (__x > __param._M_d)
1496 continue;
1497 }
1498 else if (__u <= __c3)
1499 // NB: This case not in the book, nor in the Errata,
1500 // but should be ok...
1501 __x = -1;
1502 else if (__u <= __c4)
1503 __x = 0;
1504 else if (__u <= __c5)
1505 {
1506 __x = 1;
1507 // Only in the Errata, see libstdc++/83237.
1508 __w = __178;
1509 }
1510 else
1511 {
1512 const double __v = -std::log(1.0 - __aurng());
1513 const double __y = __param._M_d
1514 + __v * __2cx / __param._M_d;
1515 __x = std::ceil(__y);
1516 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1517 }
1518
1519 __reject = (__w - __e - __x * __param._M_lm_thr
1520 > __param._M_lfm - std::lgamma(__x + __m + 1));
1521
1522 __reject |= __x + __m >= __thr;
1523
1524 } while (__reject);
1525
1526 return result_type(__x + __m + __naf);
1527 }
1528 else
1529#endif
1530 {
1531 _IntType __x = 0;
1532 double __prod = 1.0;
1533
1534 do
1535 {
1536 __prod *= __aurng();
1537 __x += 1;
1538 }
1539 while (__prod > __param._M_lm_thr);
1540
1541 return __x - 1;
1542 }
1543 }
1544
1545 template<typename _IntType>
1546 template<typename _ForwardIterator,
1547 typename _UniformRandomNumberGenerator>
1548 void
1549 poisson_distribution<_IntType>::
1550 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1551 _UniformRandomNumberGenerator& __urng,
1552 const param_type& __param)
1553 {
1554 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1555 // We could duplicate everything from operator()...
1556 while (__f != __t)
1557 *__f++ = this->operator()(__urng, __param);
1558 }
1559
1560 template<typename _IntType,
1561 typename _CharT, typename _Traits>
1562 std::basic_ostream<_CharT, _Traits>&
1563 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1564 const poisson_distribution<_IntType>& __x)
1565 {
1566 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1567
1568 const typename __ios_base::fmtflags __flags = __os.flags();
1569 const _CharT __fill = __os.fill();
1570 const std::streamsize __precision = __os.precision();
1571 const _CharT __space = __os.widen(' ');
1572 __os.flags(__ios_base::scientific | __ios_base::left);
1573 __os.fill(__space);
1574 __os.precision(std::numeric_limits<double>::max_digits10);
1575
1576 __os << __x.mean() << __space << __x._M_nd;
1577
1578 __os.flags(__flags);
1579 __os.fill(__fill);
1580 __os.precision(__precision);
1581 return __os;
1582 }
1583
1584 template<typename _IntType,
1585 typename _CharT, typename _Traits>
1586 std::basic_istream<_CharT, _Traits>&
1587 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1588 poisson_distribution<_IntType>& __x)
1589 {
1590 using param_type = typename poisson_distribution<_IntType>::param_type;
1591 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1592
1593 const typename __ios_base::fmtflags __flags = __is.flags();
1594 __is.flags(__ios_base::skipws);
1595
1596 double __mean;
1597 if (__is >> __mean >> __x._M_nd)
1598 __x.param(param_type(__mean));
1599
1600 __is.flags(__flags);
1601 return __is;
1602 }
1603
1604
1605 template<typename _IntType>
1606 void
1607 binomial_distribution<_IntType>::param_type::
1608 _M_initialize()
1609 {
1610 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1611
1612 _M_easy = true;
1613
1614#if _GLIBCXX_USE_C99_MATH_FUNCS
1615 if (_M_t * __p12 >= 8)
1616 {
1617 _M_easy = false;
1618 const double __np = std::floor(_M_t * __p12);
1619 const double __pa = __np / _M_t;
1620 const double __1p = 1 - __pa;
1621
1622 const double __pi_4 = 0.7853981633974483096156608458198757L;
1623 const double __d1x =
1624 std::sqrt(__np * __1p * std::log(32 * __np
1625 / (81 * __pi_4 * __1p)));
1626 _M_d1 = std::round(std::max<double>(1.0, __d1x));
1627 const double __d2x =
1628 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1629 / (__pi_4 * __pa)));
1630 _M_d2 = std::round(std::max<double>(1.0, __d2x));
1631
1632 // sqrt(pi / 2)
1633 const double __spi_2 = 1.2533141373155002512078826424055226L;
1634 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1635 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * (_M_t * __1p)));
1636 _M_c = 2 * _M_d1 / __np;
1637 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1638 const double __a12 = _M_a1 + _M_s2 * __spi_2;
1639 const double __s1s = _M_s1 * _M_s1;
1640 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1641 * 2 * __s1s / _M_d1
1642 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1643 const double __s2s = _M_s2 * _M_s2;
1644 _M_s = (_M_a123 + 2 * __s2s / _M_d2
1645 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1646 _M_lf = (std::lgamma(__np + 1)
1647 + std::lgamma(_M_t - __np + 1));
1648 _M_lp1p = std::log(__pa / __1p);
1649
1650 _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1651 }
1652 else
1653#endif
1654 _M_q = -std::log(1 - __p12);
1655 }
1656
1657 template<typename _IntType>
1658 template<typename _UniformRandomNumberGenerator>
1659 typename binomial_distribution<_IntType>::result_type
1660 binomial_distribution<_IntType>::
1661 _M_waiting(_UniformRandomNumberGenerator& __urng,
1662 _IntType __t, double __q)
1663 {
1664 _IntType __x = 0;
1665 double __sum = 0.0;
1666 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1667 __aurng(__urng);
1668
1669 do
1670 {
1671 if (__t == __x)
1672 return __x;
1673 const double __e = -std::log(1.0 - __aurng());
1674 __sum += __e / (__t - __x);
1675 __x += 1;
1676 }
1677 while (__sum <= __q);
1678
1679 return __x - 1;
1680 }
1681
1682 /**
1683 * A rejection algorithm when t * p >= 8 and a simple waiting time
1684 * method - the second in the referenced book - otherwise.
1685 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_FUNCS
1686 * is defined.
1687 *
1688 * Reference:
1689 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1690 * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1691 */
1692 template<typename _IntType>
1693 template<typename _UniformRandomNumberGenerator>
1694 typename binomial_distribution<_IntType>::result_type
1695 binomial_distribution<_IntType>::
1696 operator()(_UniformRandomNumberGenerator& __urng,
1697 const param_type& __param)
1698 {
1699 result_type __ret;
1700 const _IntType __t = __param.t();
1701 const double __p = __param.p();
1702 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1703 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1704 __aurng(__urng);
1705
1706#if _GLIBCXX_USE_C99_MATH_FUNCS
1707 if (!__param._M_easy)
1708 {
1709 double __x;
1710
1711 // See comments above...
1712 const double __naf =
1713 (1 - std::numeric_limits<double>::epsilon()) / 2;
1714 const double __thr =
1715 std::numeric_limits<_IntType>::max() + __naf;
1716
1717 const double __np = std::floor(__t * __p12);
1718
1719 // sqrt(pi / 2)
1720 const double __spi_2 = 1.2533141373155002512078826424055226L;
1721 const double __a1 = __param._M_a1;
1722 const double __a12 = __a1 + __param._M_s2 * __spi_2;
1723 const double __a123 = __param._M_a123;
1724 const double __s1s = __param._M_s1 * __param._M_s1;
1725 const double __s2s = __param._M_s2 * __param._M_s2;
1726
1727 bool __reject;
1728 do
1729 {
1730 const double __u = __param._M_s * __aurng();
1731
1732 double __v;
1733
1734 if (__u <= __a1)
1735 {
1736 const double __n = _M_nd(__urng);
1737 const double __y = __param._M_s1 * std::abs(__n);
1738 __reject = __y >= __param._M_d1;
1739 if (!__reject)
1740 {
1741 const double __e = -std::log(1.0 - __aurng());
1742 __x = std::floor(__y);
1743 __v = -__e - __n * __n / 2 + __param._M_c;
1744 }
1745 }
1746 else if (__u <= __a12)
1747 {
1748 const double __n = _M_nd(__urng);
1749 const double __y = __param._M_s2 * std::abs(__n);
1750 __reject = __y >= __param._M_d2;
1751 if (!__reject)
1752 {
1753 const double __e = -std::log(1.0 - __aurng());
1754 __x = std::floor(-__y);
1755 __v = -__e - __n * __n / 2;
1756 }
1757 }
1758 else if (__u <= __a123)
1759 {
1760 const double __e1 = -std::log(1.0 - __aurng());
1761 const double __e2 = -std::log(1.0 - __aurng());
1762
1763 const double __y = __param._M_d1
1764 + 2 * __s1s * __e1 / __param._M_d1;
1765 __x = std::floor(__y);
1766 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1767 -__y / (2 * __s1s)));
1768 __reject = false;
1769 }
1770 else
1771 {
1772 const double __e1 = -std::log(1.0 - __aurng());
1773 const double __e2 = -std::log(1.0 - __aurng());
1774
1775 const double __y = __param._M_d2
1776 + 2 * __s2s * __e1 / __param._M_d2;
1777 __x = std::floor(-__y);
1778 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1779 __reject = false;
1780 }
1781
1782 __reject = __reject || __x < -__np || __x > __t - __np;
1783 if (!__reject)
1784 {
1785 const double __lfx =
1786 std::lgamma(__np + __x + 1)
1787 + std::lgamma(__t - (__np + __x) + 1);
1788 __reject = __v > __param._M_lf - __lfx
1789 + __x * __param._M_lp1p;
1790 }
1791
1792 __reject |= __x + __np >= __thr;
1793 }
1794 while (__reject);
1795
1796 __x += __np + __naf;
1797
1798 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1799 __param._M_q);
1800 __ret = _IntType(__x) + __z;
1801 }
1802 else
1803#endif
1804 __ret = _M_waiting(__urng, __t, __param._M_q);
1805
1806 if (__p12 != __p)
1807 __ret = __t - __ret;
1808 return __ret;
1809 }
1810
1811 template<typename _IntType>
1812 template<typename _ForwardIterator,
1813 typename _UniformRandomNumberGenerator>
1814 void
1815 binomial_distribution<_IntType>::
1816 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1817 _UniformRandomNumberGenerator& __urng,
1818 const param_type& __param)
1819 {
1820 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1821 // We could duplicate everything from operator()...
1822 while (__f != __t)
1823 *__f++ = this->operator()(__urng, __param);
1824 }
1825
1826 template<typename _IntType,
1827 typename _CharT, typename _Traits>
1828 std::basic_ostream<_CharT, _Traits>&
1829 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1830 const binomial_distribution<_IntType>& __x)
1831 {
1832 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1833
1834 const typename __ios_base::fmtflags __flags = __os.flags();
1835 const _CharT __fill = __os.fill();
1836 const std::streamsize __precision = __os.precision();
1837 const _CharT __space = __os.widen(' ');
1838 __os.flags(__ios_base::scientific | __ios_base::left);
1839 __os.fill(__space);
1840 __os.precision(std::numeric_limits<double>::max_digits10);
1841
1842 __os << __x.t() << __space << __x.p()
1843 << __space << __x._M_nd;
1844
1845 __os.flags(__flags);
1846 __os.fill(__fill);
1847 __os.precision(__precision);
1848 return __os;
1849 }
1850
1851 template<typename _IntType,
1852 typename _CharT, typename _Traits>
1853 std::basic_istream<_CharT, _Traits>&
1854 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1855 binomial_distribution<_IntType>& __x)
1856 {
1857 using param_type = typename binomial_distribution<_IntType>::param_type;
1858 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1859
1860 const typename __ios_base::fmtflags __flags = __is.flags();
1861 __is.flags(__ios_base::dec | __ios_base::skipws);
1862
1863 _IntType __t;
1864 double __p;
1865 if (__is >> __t >> __p >> __x._M_nd)
1866 __x.param(param_type(__t, __p));
1867
1868 __is.flags(__flags);
1869 return __is;
1870 }
1871
1872
1873 template<typename _RealType>
1874 template<typename _ForwardIterator,
1875 typename _UniformRandomNumberGenerator>
1876 void
1877 std::exponential_distribution<_RealType>::
1878 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1879 _UniformRandomNumberGenerator& __urng,
1880 const param_type& __p)
1881 {
1882 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1883 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1884 __aurng(__urng);
1885 while (__f != __t)
1886 *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1887 }
1888
1889 template<typename _RealType, typename _CharT, typename _Traits>
1890 std::basic_ostream<_CharT, _Traits>&
1891 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1892 const exponential_distribution<_RealType>& __x)
1893 {
1894 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1895
1896 const typename __ios_base::fmtflags __flags = __os.flags();
1897 const _CharT __fill = __os.fill();
1898 const std::streamsize __precision = __os.precision();
1899 __os.flags(__ios_base::scientific | __ios_base::left);
1900 __os.fill(__os.widen(' '));
1901 __os.precision(std::numeric_limits<_RealType>::max_digits10);
1902
1903 __os << __x.lambda();
1904
1905 __os.flags(__flags);
1906 __os.fill(__fill);
1907 __os.precision(__precision);
1908 return __os;
1909 }
1910
1911 template<typename _RealType, typename _CharT, typename _Traits>
1912 std::basic_istream<_CharT, _Traits>&
1913 operator>>(std::basic_istream<_CharT, _Traits>& __is,
1914 exponential_distribution<_RealType>& __x)
1915 {
1916 using param_type
1917 = typename exponential_distribution<_RealType>::param_type;
1918 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1919
1920 const typename __ios_base::fmtflags __flags = __is.flags();
1921 __is.flags(__ios_base::dec | __ios_base::skipws);
1922
1923 _RealType __lambda;
1924 if (__is >> __lambda)
1925 __x.param(param_type(__lambda));
1926
1927 __is.flags(__flags);
1928 return __is;
1929 }
1930
1931
1932 /**
1933 * Polar method due to Marsaglia.
1934 *
1935 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1936 * New York, 1986, Ch. V, Sect. 4.4.
1937 */
1938 template<typename _RealType>
1939 template<typename _UniformRandomNumberGenerator>
1940 typename normal_distribution<_RealType>::result_type
1941 normal_distribution<_RealType>::
1942 operator()(_UniformRandomNumberGenerator& __urng,
1943 const param_type& __param)
1944 {
1945 result_type __ret;
1946 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1947 __aurng(__urng);
1948
1949 if (_M_saved_available)
1950 {
1951 _M_saved_available = false;
1952 __ret = _M_saved;
1953 }
1954 else
1955 {
1956 result_type __x, __y, __r2;
1957 do
1958 {
1959 __x = result_type(2.0) * __aurng() - 1.0;
1960 __y = result_type(2.0) * __aurng() - 1.0;
1961 __r2 = __x * __x + __y * __y;
1962 }
1963 while (__r2 > 1.0 || __r2 == 0.0);
1964
1965 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1966 _M_saved = __x * __mult;
1967 _M_saved_available = true;
1968 __ret = __y * __mult;
1969 }
1970
1971 __ret = __ret * __param.stddev() + __param.mean();
1972 return __ret;
1973 }
1974
1975 template<typename _RealType>
1976 template<typename _ForwardIterator,
1977 typename _UniformRandomNumberGenerator>
1978 void
1979 normal_distribution<_RealType>::
1980 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1981 _UniformRandomNumberGenerator& __urng,
1982 const param_type& __param)
1983 {
1984 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1985
1986 if (__f == __t)
1987 return;
1988
1989 if (_M_saved_available)
1990 {
1991 _M_saved_available = false;
1992 *__f++ = _M_saved * __param.stddev() + __param.mean();
1993
1994 if (__f == __t)
1995 return;
1996 }
1997
1998 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1999 __aurng(__urng);
2000
2001 while (__f + 1 < __t)
2002 {
2003 result_type __x, __y, __r2;
2004 do
2005 {
2006 __x = result_type(2.0) * __aurng() - 1.0;
2007 __y = result_type(2.0) * __aurng() - 1.0;
2008 __r2 = __x * __x + __y * __y;
2009 }
2010 while (__r2 > 1.0 || __r2 == 0.0);
2011
2012 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2013 *__f++ = __y * __mult * __param.stddev() + __param.mean();
2014 *__f++ = __x * __mult * __param.stddev() + __param.mean();
2015 }
2016
2017 if (__f != __t)
2018 {
2019 result_type __x, __y, __r2;
2020 do
2021 {
2022 __x = result_type(2.0) * __aurng() - 1.0;
2023 __y = result_type(2.0) * __aurng() - 1.0;
2024 __r2 = __x * __x + __y * __y;
2025 }
2026 while (__r2 > 1.0 || __r2 == 0.0);
2027
2028 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2029 _M_saved = __x * __mult;
2030 _M_saved_available = true;
2031 *__f = __y * __mult * __param.stddev() + __param.mean();
2032 }
2033 }
2034
2035 template<typename _RealType>
2036 bool
2037 operator==(const std::normal_distribution<_RealType>& __d1,
2038 const std::normal_distribution<_RealType>& __d2)
2039 {
2040 if (__d1._M_param == __d2._M_param
2041 && __d1._M_saved_available == __d2._M_saved_available)
2042 return __d1._M_saved_available ? __d1._M_saved == __d2._M_saved : true;
2043 else
2044 return false;
2045 }
2046
2047 template<typename _RealType, typename _CharT, typename _Traits>
2048 std::basic_ostream<_CharT, _Traits>&
2049 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2050 const normal_distribution<_RealType>& __x)
2051 {
2052 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2053
2054 const typename __ios_base::fmtflags __flags = __os.flags();
2055 const _CharT __fill = __os.fill();
2056 const std::streamsize __precision = __os.precision();
2057 const _CharT __space = __os.widen(' ');
2058 __os.flags(__ios_base::scientific | __ios_base::left);
2059 __os.fill(__space);
2060 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2061
2062 __os << __x.mean() << __space << __x.stddev()
2063 << __space << __x._M_saved_available;
2064 if (__x._M_saved_available)
2065 __os << __space << __x._M_saved;
2066
2067 __os.flags(__flags);
2068 __os.fill(__fill);
2069 __os.precision(__precision);
2070 return __os;
2071 }
2072
2073 template<typename _RealType, typename _CharT, typename _Traits>
2074 std::basic_istream<_CharT, _Traits>&
2075 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2076 normal_distribution<_RealType>& __x)
2077 {
2078 using param_type = typename normal_distribution<_RealType>::param_type;
2079 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2080
2081 const typename __ios_base::fmtflags __flags = __is.flags();
2082 __is.flags(__ios_base::dec | __ios_base::skipws);
2083
2084 double __mean, __stddev;
2085 bool __saved_avail;
2086 if (__is >> __mean >> __stddev >> __saved_avail)
2087 {
2088 if (!__saved_avail || (__is >> __x._M_saved))
2089 {
2090 __x._M_saved_available = __saved_avail;
2091 __x.param(param_type(__mean, __stddev));
2092 }
2093 }
2094
2095 __is.flags(__flags);
2096 return __is;
2097 }
2098
2099
2100 template<typename _RealType>
2101 template<typename _ForwardIterator,
2102 typename _UniformRandomNumberGenerator>
2103 void
2104 lognormal_distribution<_RealType>::
2105 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2106 _UniformRandomNumberGenerator& __urng,
2107 const param_type& __p)
2108 {
2109 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2110 while (__f != __t)
2111 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
2112 }
2113
2114 template<typename _RealType, typename _CharT, typename _Traits>
2115 std::basic_ostream<_CharT, _Traits>&
2116 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2117 const lognormal_distribution<_RealType>& __x)
2118 {
2119 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2120
2121 const typename __ios_base::fmtflags __flags = __os.flags();
2122 const _CharT __fill = __os.fill();
2123 const std::streamsize __precision = __os.precision();
2124 const _CharT __space = __os.widen(' ');
2125 __os.flags(__ios_base::scientific | __ios_base::left);
2126 __os.fill(__space);
2127 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2128
2129 __os << __x.m() << __space << __x.s()
2130 << __space << __x._M_nd;
2131
2132 __os.flags(__flags);
2133 __os.fill(__fill);
2134 __os.precision(__precision);
2135 return __os;
2136 }
2137
2138 template<typename _RealType, typename _CharT, typename _Traits>
2139 std::basic_istream<_CharT, _Traits>&
2140 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2141 lognormal_distribution<_RealType>& __x)
2142 {
2143 using param_type
2144 = typename lognormal_distribution<_RealType>::param_type;
2145 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2146
2147 const typename __ios_base::fmtflags __flags = __is.flags();
2148 __is.flags(__ios_base::dec | __ios_base::skipws);
2149
2150 _RealType __m, __s;
2151 if (__is >> __m >> __s >> __x._M_nd)
2152 __x.param(param_type(__m, __s));
2153
2154 __is.flags(__flags);
2155 return __is;
2156 }
2157
2158 template<typename _RealType>
2159 template<typename _ForwardIterator,
2160 typename _UniformRandomNumberGenerator>
2161 void
2162 std::chi_squared_distribution<_RealType>::
2163 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2164 _UniformRandomNumberGenerator& __urng)
2165 {
2166 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2167 while (__f != __t)
2168 *__f++ = 2 * _M_gd(__urng);
2169 }
2170
2171 template<typename _RealType>
2172 template<typename _ForwardIterator,
2173 typename _UniformRandomNumberGenerator>
2174 void
2175 std::chi_squared_distribution<_RealType>::
2176 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2177 _UniformRandomNumberGenerator& __urng,
2178 const typename
2179 std::gamma_distribution<result_type>::param_type& __p)
2180 {
2181 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2182 while (__f != __t)
2183 *__f++ = 2 * _M_gd(__urng, __p);
2184 }
2185
2186 template<typename _RealType, typename _CharT, typename _Traits>
2187 std::basic_ostream<_CharT, _Traits>&
2188 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2189 const chi_squared_distribution<_RealType>& __x)
2190 {
2191 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2192
2193 const typename __ios_base::fmtflags __flags = __os.flags();
2194 const _CharT __fill = __os.fill();
2195 const std::streamsize __precision = __os.precision();
2196 const _CharT __space = __os.widen(' ');
2197 __os.flags(__ios_base::scientific | __ios_base::left);
2198 __os.fill(__space);
2199 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2200
2201 __os << __x.n() << __space << __x._M_gd;
2202
2203 __os.flags(__flags);
2204 __os.fill(__fill);
2205 __os.precision(__precision);
2206 return __os;
2207 }
2208
2209 template<typename _RealType, typename _CharT, typename _Traits>
2210 std::basic_istream<_CharT, _Traits>&
2211 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2212 chi_squared_distribution<_RealType>& __x)
2213 {
2214 using param_type
2215 = typename chi_squared_distribution<_RealType>::param_type;
2216 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2217
2218 const typename __ios_base::fmtflags __flags = __is.flags();
2219 __is.flags(__ios_base::dec | __ios_base::skipws);
2220
2221 _RealType __n;
2222 if (__is >> __n >> __x._M_gd)
2223 __x.param(param_type(__n));
2224
2225 __is.flags(__flags);
2226 return __is;
2227 }
2228
2229
2230 template<typename _RealType>
2231 template<typename _UniformRandomNumberGenerator>
2232 typename cauchy_distribution<_RealType>::result_type
2233 cauchy_distribution<_RealType>::
2234 operator()(_UniformRandomNumberGenerator& __urng,
2235 const param_type& __p)
2236 {
2237 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2238 __aurng(__urng);
2239 _RealType __u;
2240 do
2241 __u = __aurng();
2242 while (__u == 0.5);
2243
2244 const _RealType __pi = 3.1415926535897932384626433832795029L;
2245 return __p.a() + __p.b() * std::tan(__pi * __u);
2246 }
2247
2248 template<typename _RealType>
2249 template<typename _ForwardIterator,
2250 typename _UniformRandomNumberGenerator>
2251 void
2252 cauchy_distribution<_RealType>::
2253 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2254 _UniformRandomNumberGenerator& __urng,
2255 const param_type& __p)
2256 {
2257 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2258 const _RealType __pi = 3.1415926535897932384626433832795029L;
2259 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2260 __aurng(__urng);
2261 while (__f != __t)
2262 {
2263 _RealType __u;
2264 do
2265 __u = __aurng();
2266 while (__u == 0.5);
2267
2268 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2269 }
2270 }
2271
2272 template<typename _RealType, typename _CharT, typename _Traits>
2273 std::basic_ostream<_CharT, _Traits>&
2274 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2275 const cauchy_distribution<_RealType>& __x)
2276 {
2277 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2278
2279 const typename __ios_base::fmtflags __flags = __os.flags();
2280 const _CharT __fill = __os.fill();
2281 const std::streamsize __precision = __os.precision();
2282 const _CharT __space = __os.widen(' ');
2283 __os.flags(__ios_base::scientific | __ios_base::left);
2284 __os.fill(__space);
2285 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2286
2287 __os << __x.a() << __space << __x.b();
2288
2289 __os.flags(__flags);
2290 __os.fill(__fill);
2291 __os.precision(__precision);
2292 return __os;
2293 }
2294
2295 template<typename _RealType, typename _CharT, typename _Traits>
2296 std::basic_istream<_CharT, _Traits>&
2297 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2298 cauchy_distribution<_RealType>& __x)
2299 {
2300 using param_type = typename cauchy_distribution<_RealType>::param_type;
2301 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2302
2303 const typename __ios_base::fmtflags __flags = __is.flags();
2304 __is.flags(__ios_base::dec | __ios_base::skipws);
2305
2306 _RealType __a, __b;
2307 if (__is >> __a >> __b)
2308 __x.param(param_type(__a, __b));
2309
2310 __is.flags(__flags);
2311 return __is;
2312 }
2313
2314
2315 template<typename _RealType>
2316 template<typename _ForwardIterator,
2317 typename _UniformRandomNumberGenerator>
2318 void
2319 std::fisher_f_distribution<_RealType>::
2320 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2321 _UniformRandomNumberGenerator& __urng)
2322 {
2323 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2324 while (__f != __t)
2325 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2326 }
2327
2328 template<typename _RealType>
2329 template<typename _ForwardIterator,
2330 typename _UniformRandomNumberGenerator>
2331 void
2332 std::fisher_f_distribution<_RealType>::
2333 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2334 _UniformRandomNumberGenerator& __urng,
2335 const param_type& __p)
2336 {
2337 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2338 typedef typename std::gamma_distribution<result_type>::param_type
2339 param_type;
2340 param_type __p1(__p.m() / 2);
2341 param_type __p2(__p.n() / 2);
2342 while (__f != __t)
2343 *__f++ = ((_M_gd_x(__urng, __p1) * n())
2344 / (_M_gd_y(__urng, __p2) * m()));
2345 }
2346
2347 template<typename _RealType, typename _CharT, typename _Traits>
2348 std::basic_ostream<_CharT, _Traits>&
2349 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2350 const fisher_f_distribution<_RealType>& __x)
2351 {
2352 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2353
2354 const typename __ios_base::fmtflags __flags = __os.flags();
2355 const _CharT __fill = __os.fill();
2356 const std::streamsize __precision = __os.precision();
2357 const _CharT __space = __os.widen(' ');
2358 __os.flags(__ios_base::scientific | __ios_base::left);
2359 __os.fill(__space);
2360 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2361
2362 __os << __x.m() << __space << __x.n()
2363 << __space << __x._M_gd_x << __space << __x._M_gd_y;
2364
2365 __os.flags(__flags);
2366 __os.fill(__fill);
2367 __os.precision(__precision);
2368 return __os;
2369 }
2370
2371 template<typename _RealType, typename _CharT, typename _Traits>
2372 std::basic_istream<_CharT, _Traits>&
2373 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2374 fisher_f_distribution<_RealType>& __x)
2375 {
2376 using param_type
2377 = typename fisher_f_distribution<_RealType>::param_type;
2378 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2379
2380 const typename __ios_base::fmtflags __flags = __is.flags();
2381 __is.flags(__ios_base::dec | __ios_base::skipws);
2382
2383 _RealType __m, __n;
2384 if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y)
2385 __x.param(param_type(__m, __n));
2386
2387 __is.flags(__flags);
2388 return __is;
2389 }
2390
2391
2392 template<typename _RealType>
2393 template<typename _ForwardIterator,
2394 typename _UniformRandomNumberGenerator>
2395 void
2396 std::student_t_distribution<_RealType>::
2397 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2398 _UniformRandomNumberGenerator& __urng)
2399 {
2400 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2401 while (__f != __t)
2402 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2403 }
2404
2405 template<typename _RealType>
2406 template<typename _ForwardIterator,
2407 typename _UniformRandomNumberGenerator>
2408 void
2409 std::student_t_distribution<_RealType>::
2410 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2411 _UniformRandomNumberGenerator& __urng,
2412 const param_type& __p)
2413 {
2414 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2415 typename std::gamma_distribution<result_type>::param_type
2416 __p2(__p.n() / 2, 2);
2417 while (__f != __t)
2418 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2419 }
2420
2421 template<typename _RealType, typename _CharT, typename _Traits>
2422 std::basic_ostream<_CharT, _Traits>&
2423 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2424 const student_t_distribution<_RealType>& __x)
2425 {
2426 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2427
2428 const typename __ios_base::fmtflags __flags = __os.flags();
2429 const _CharT __fill = __os.fill();
2430 const std::streamsize __precision = __os.precision();
2431 const _CharT __space = __os.widen(' ');
2432 __os.flags(__ios_base::scientific | __ios_base::left);
2433 __os.fill(__space);
2434 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2435
2436 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2437
2438 __os.flags(__flags);
2439 __os.fill(__fill);
2440 __os.precision(__precision);
2441 return __os;
2442 }
2443
2444 template<typename _RealType, typename _CharT, typename _Traits>
2445 std::basic_istream<_CharT, _Traits>&
2446 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2447 student_t_distribution<_RealType>& __x)
2448 {
2449 using param_type
2450 = typename student_t_distribution<_RealType>::param_type;
2451 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2452
2453 const typename __ios_base::fmtflags __flags = __is.flags();
2454 __is.flags(__ios_base::dec | __ios_base::skipws);
2455
2456 _RealType __n;
2457 if (__is >> __n >> __x._M_nd >> __x._M_gd)
2458 __x.param(param_type(__n));
2459
2460 __is.flags(__flags);
2461 return __is;
2462 }
2463
2464
2465 template<typename _RealType>
2466 void
2467 gamma_distribution<_RealType>::param_type::
2468 _M_initialize()
2469 {
2470 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2471
2472 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2473 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2474 }
2475
2476 /**
2477 * Marsaglia, G. and Tsang, W. W.
2478 * "A Simple Method for Generating Gamma Variables"
2479 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2480 */
2481 template<typename _RealType>
2482 template<typename _UniformRandomNumberGenerator>
2483 typename gamma_distribution<_RealType>::result_type
2484 gamma_distribution<_RealType>::
2485 operator()(_UniformRandomNumberGenerator& __urng,
2486 const param_type& __param)
2487 {
2488 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2489 __aurng(__urng);
2490
2491 result_type __u, __v, __n;
2492 const result_type __a1 = (__param._M_malpha
2493 - _RealType(1.0) / _RealType(3.0));
2494
2495 do
2496 {
2497 do
2498 {
2499 __n = _M_nd(__urng);
2500 __v = result_type(1.0) + __param._M_a2 * __n;
2501 }
2502 while (__v <= 0.0);
2503
2504 __v = __v * __v * __v;
2505 __u = __aurng();
2506 }
2507 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2508 && (std::log(__u) > (0.5 * __n * __n + __a1
2509 * (1.0 - __v + std::log(__v)))));
2510
2511 if (__param.alpha() == __param._M_malpha)
2512 return __a1 * __v * __param.beta();
2513 else
2514 {
2515 do
2516 __u = __aurng();
2517 while (__u == 0.0);
2518
2519 return (std::pow(__u, result_type(1.0) / __param.alpha())
2520 * __a1 * __v * __param.beta());
2521 }
2522 }
2523
2524 template<typename _RealType>
2525 template<typename _ForwardIterator,
2526 typename _UniformRandomNumberGenerator>
2527 void
2528 gamma_distribution<_RealType>::
2529 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2530 _UniformRandomNumberGenerator& __urng,
2531 const param_type& __param)
2532 {
2533 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2534 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2535 __aurng(__urng);
2536
2537 result_type __u, __v, __n;
2538 const result_type __a1 = (__param._M_malpha
2539 - _RealType(1.0) / _RealType(3.0));
2540
2541 if (__param.alpha() == __param._M_malpha)
2542 while (__f != __t)
2543 {
2544 do
2545 {
2546 do
2547 {
2548 __n = _M_nd(__urng);
2549 __v = result_type(1.0) + __param._M_a2 * __n;
2550 }
2551 while (__v <= 0.0);
2552
2553 __v = __v * __v * __v;
2554 __u = __aurng();
2555 }
2556 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2557 && (std::log(__u) > (0.5 * __n * __n + __a1
2558 * (1.0 - __v + std::log(__v)))));
2559
2560 *__f++ = __a1 * __v * __param.beta();
2561 }
2562 else
2563 while (__f != __t)
2564 {
2565 do
2566 {
2567 do
2568 {
2569 __n = _M_nd(__urng);
2570 __v = result_type(1.0) + __param._M_a2 * __n;
2571 }
2572 while (__v <= 0.0);
2573
2574 __v = __v * __v * __v;
2575 __u = __aurng();
2576 }
2577 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2578 && (std::log(__u) > (0.5 * __n * __n + __a1
2579 * (1.0 - __v + std::log(__v)))));
2580
2581 do
2582 __u = __aurng();
2583 while (__u == 0.0);
2584
2585 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2586 * __a1 * __v * __param.beta());
2587 }
2588 }
2589
2590 template<typename _RealType, typename _CharT, typename _Traits>
2591 std::basic_ostream<_CharT, _Traits>&
2592 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2593 const gamma_distribution<_RealType>& __x)
2594 {
2595 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2596
2597 const typename __ios_base::fmtflags __flags = __os.flags();
2598 const _CharT __fill = __os.fill();
2599 const std::streamsize __precision = __os.precision();
2600 const _CharT __space = __os.widen(' ');
2601 __os.flags(__ios_base::scientific | __ios_base::left);
2602 __os.fill(__space);
2603 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2604
2605 __os << __x.alpha() << __space << __x.beta()
2606 << __space << __x._M_nd;
2607
2608 __os.flags(__flags);
2609 __os.fill(__fill);
2610 __os.precision(__precision);
2611 return __os;
2612 }
2613
2614 template<typename _RealType, typename _CharT, typename _Traits>
2615 std::basic_istream<_CharT, _Traits>&
2616 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2617 gamma_distribution<_RealType>& __x)
2618 {
2619 using param_type = typename gamma_distribution<_RealType>::param_type;
2620 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2621
2622 const typename __ios_base::fmtflags __flags = __is.flags();
2623 __is.flags(__ios_base::dec | __ios_base::skipws);
2624
2625 _RealType __alpha_val, __beta_val;
2626 if (__is >> __alpha_val >> __beta_val >> __x._M_nd)
2627 __x.param(param_type(__alpha_val, __beta_val));
2628
2629 __is.flags(__flags);
2630 return __is;
2631 }
2632
2633
2634 template<typename _RealType>
2635 template<typename _UniformRandomNumberGenerator>
2636 typename weibull_distribution<_RealType>::result_type
2637 weibull_distribution<_RealType>::
2638 operator()(_UniformRandomNumberGenerator& __urng,
2639 const param_type& __p)
2640 {
2641 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2642 __aurng(__urng);
2643 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2644 result_type(1) / __p.a());
2645 }
2646
2647 template<typename _RealType>
2648 template<typename _ForwardIterator,
2649 typename _UniformRandomNumberGenerator>
2650 void
2651 weibull_distribution<_RealType>::
2652 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2653 _UniformRandomNumberGenerator& __urng,
2654 const param_type& __p)
2655 {
2656 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2657 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2658 __aurng(__urng);
2659 auto __inv_a = result_type(1) / __p.a();
2660
2661 while (__f != __t)
2662 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2663 __inv_a);
2664 }
2665
2666 template<typename _RealType, typename _CharT, typename _Traits>
2667 std::basic_ostream<_CharT, _Traits>&
2668 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2669 const weibull_distribution<_RealType>& __x)
2670 {
2671 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2672
2673 const typename __ios_base::fmtflags __flags = __os.flags();
2674 const _CharT __fill = __os.fill();
2675 const std::streamsize __precision = __os.precision();
2676 const _CharT __space = __os.widen(' ');
2677 __os.flags(__ios_base::scientific | __ios_base::left);
2678 __os.fill(__space);
2679 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2680
2681 __os << __x.a() << __space << __x.b();
2682
2683 __os.flags(__flags);
2684 __os.fill(__fill);
2685 __os.precision(__precision);
2686 return __os;
2687 }
2688
2689 template<typename _RealType, typename _CharT, typename _Traits>
2690 std::basic_istream<_CharT, _Traits>&
2691 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2692 weibull_distribution<_RealType>& __x)
2693 {
2694 using param_type = typename weibull_distribution<_RealType>::param_type;
2695 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2696
2697 const typename __ios_base::fmtflags __flags = __is.flags();
2698 __is.flags(__ios_base::dec | __ios_base::skipws);
2699
2700 _RealType __a, __b;
2701 if (__is >> __a >> __b)
2702 __x.param(param_type(__a, __b));
2703
2704 __is.flags(__flags);
2705 return __is;
2706 }
2707
2708
2709 template<typename _RealType>
2710 template<typename _UniformRandomNumberGenerator>
2711 typename extreme_value_distribution<_RealType>::result_type
2712 extreme_value_distribution<_RealType>::
2713 operator()(_UniformRandomNumberGenerator& __urng,
2714 const param_type& __p)
2715 {
2716 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2717 __aurng(__urng);
2718 return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2719 - __aurng()));
2720 }
2721
2722 template<typename _RealType>
2723 template<typename _ForwardIterator,
2724 typename _UniformRandomNumberGenerator>
2725 void
2726 extreme_value_distribution<_RealType>::
2727 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2728 _UniformRandomNumberGenerator& __urng,
2729 const param_type& __p)
2730 {
2731 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2732 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2733 __aurng(__urng);
2734
2735 while (__f != __t)
2736 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2737 - __aurng()));
2738 }
2739
2740 template<typename _RealType, typename _CharT, typename _Traits>
2741 std::basic_ostream<_CharT, _Traits>&
2742 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2743 const extreme_value_distribution<_RealType>& __x)
2744 {
2745 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2746
2747 const typename __ios_base::fmtflags __flags = __os.flags();
2748 const _CharT __fill = __os.fill();
2749 const std::streamsize __precision = __os.precision();
2750 const _CharT __space = __os.widen(' ');
2751 __os.flags(__ios_base::scientific | __ios_base::left);
2752 __os.fill(__space);
2753 __os.precision(std::numeric_limits<_RealType>::max_digits10);
2754
2755 __os << __x.a() << __space << __x.b();
2756
2757 __os.flags(__flags);
2758 __os.fill(__fill);
2759 __os.precision(__precision);
2760 return __os;
2761 }
2762
2763 template<typename _RealType, typename _CharT, typename _Traits>
2764 std::basic_istream<_CharT, _Traits>&
2765 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2766 extreme_value_distribution<_RealType>& __x)
2767 {
2768 using param_type
2769 = typename extreme_value_distribution<_RealType>::param_type;
2770 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2771
2772 const typename __ios_base::fmtflags __flags = __is.flags();
2773 __is.flags(__ios_base::dec | __ios_base::skipws);
2774
2775 _RealType __a, __b;
2776 if (__is >> __a >> __b)
2777 __x.param(param_type(__a, __b));
2778
2779 __is.flags(__flags);
2780 return __is;
2781 }
2782
2783
2784 template<typename _IntType>
2785 void
2786 discrete_distribution<_IntType>::param_type::
2787 _M_initialize()
2788 {
2789 if (_M_prob.size() < 2)
2790 {
2791 _M_prob.clear();
2792 return;
2793 }
2794
2795 const double __sum = std::accumulate(_M_prob.begin(),
2796 _M_prob.end(), 0.0);
2797 __glibcxx_assert(__sum > 0);
2798 // Now normalize the probabilites.
2799 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2800 __sum);
2801 // Accumulate partial sums.
2802 _M_cp.reserve(_M_prob.size());
2803 std::partial_sum(_M_prob.begin(), _M_prob.end(),
2804 std::back_inserter(_M_cp));
2805 // Make sure the last cumulative probability is one.
2806 _M_cp[_M_cp.size() - 1] = 1.0;
2807 }
2808
2809 template<typename _IntType>
2810 template<typename _Func>
2811 discrete_distribution<_IntType>::param_type::
2812 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2813 : _M_prob(), _M_cp()
2814 {
2815 const size_t __n = __nw == 0 ? 1 : __nw;
2816 const double __delta = (__xmax - __xmin) / __n;
2817
2818 _M_prob.reserve(__n);
2819 for (size_t __k = 0; __k < __nw; ++__k)
2820 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2821
2822 _M_initialize();
2823 }
2824
2825 template<typename _IntType>
2826 template<typename _UniformRandomNumberGenerator>
2827 typename discrete_distribution<_IntType>::result_type
2828 discrete_distribution<_IntType>::
2829 operator()(_UniformRandomNumberGenerator& __urng,
2830 const param_type& __param)
2831 {
2832 if (__param._M_cp.empty())
2833 return result_type(0);
2834
2835 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2836 __aurng(__urng);
2837
2838 const double __p = __aurng();
2839 auto __pos = std::lower_bound(__param._M_cp.begin(),
2840 __param._M_cp.end(), __p);
2841
2842 return __pos - __param._M_cp.begin();
2843 }
2844
2845 template<typename _IntType>
2846 template<typename _ForwardIterator,
2847 typename _UniformRandomNumberGenerator>
2848 void
2849 discrete_distribution<_IntType>::
2850 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2851 _UniformRandomNumberGenerator& __urng,
2852 const param_type& __param)
2853 {
2854 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2855
2856 if (__param._M_cp.empty())
2857 {
2858 while (__f != __t)
2859 *__f++ = result_type(0);
2860 return;
2861 }
2862
2863 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2864 __aurng(__urng);
2865
2866 while (__f != __t)
2867 {
2868 const double __p = __aurng();
2869 auto __pos = std::lower_bound(__param._M_cp.begin(),
2870 __param._M_cp.end(), __p);
2871
2872 *__f++ = __pos - __param._M_cp.begin();
2873 }
2874 }
2875
2876 template<typename _IntType, typename _CharT, typename _Traits>
2877 std::basic_ostream<_CharT, _Traits>&
2878 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2879 const discrete_distribution<_IntType>& __x)
2880 {
2881 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2882
2883 const typename __ios_base::fmtflags __flags = __os.flags();
2884 const _CharT __fill = __os.fill();
2885 const std::streamsize __precision = __os.precision();
2886 const _CharT __space = __os.widen(' ');
2887 __os.flags(__ios_base::scientific | __ios_base::left);
2888 __os.fill(__space);
2889 __os.precision(std::numeric_limits<double>::max_digits10);
2890
2891 std::vector<double> __prob = __x.probabilities();
2892 __os << __prob.size();
2893 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2894 __os << __space << *__dit;
2895
2896 __os.flags(__flags);
2897 __os.fill(__fill);
2898 __os.precision(__precision);
2899 return __os;
2900 }
2901
2902namespace __detail
2903{
2904 template<typename _ValT, typename _CharT, typename _Traits>
2905 basic_istream<_CharT, _Traits>&
2906 __extract_params(basic_istream<_CharT, _Traits>& __is,
2907 vector<_ValT>& __vals, size_t __n)
2908 {
2909 __vals.reserve(__n);
2910 while (__n--)
2911 {
2912 _ValT __val;
2913 if (__is >> __val)
2914 __vals.push_back(__val);
2915 else
2916 break;
2917 }
2918 return __is;
2919 }
2920} // namespace __detail
2921
2922 template<typename _IntType, typename _CharT, typename _Traits>
2923 std::basic_istream<_CharT, _Traits>&
2924 operator>>(std::basic_istream<_CharT, _Traits>& __is,
2925 discrete_distribution<_IntType>& __x)
2926 {
2927 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2928
2929 const typename __ios_base::fmtflags __flags = __is.flags();
2930 __is.flags(__ios_base::dec | __ios_base::skipws);
2931
2932 size_t __n;
2933 if (__is >> __n)
2934 {
2935 std::vector<double> __prob_vec;
2936 if (__detail::__extract_params(__is, __prob_vec, __n))
2937 __x.param({__prob_vec.begin(), __prob_vec.end()});
2938 }
2939
2940 __is.flags(__flags);
2941 return __is;
2942 }
2943
2944
2945 template<typename _RealType>
2946 void
2947 piecewise_constant_distribution<_RealType>::param_type::
2948 _M_initialize()
2949 {
2950 if (_M_int.size() < 2
2951 || (_M_int.size() == 2
2952 && _M_int[0] == _RealType(0)
2953 && _M_int[1] == _RealType(1)))
2954 {
2955 _M_int.clear();
2956 _M_den.clear();
2957 return;
2958 }
2959
2960 const double __sum = std::accumulate(_M_den.begin(),
2961 _M_den.end(), 0.0);
2962 __glibcxx_assert(__sum > 0);
2963
2964 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2965 __sum);
2966
2967 _M_cp.reserve(_M_den.size());
2968 std::partial_sum(_M_den.begin(), _M_den.end(),
2969 std::back_inserter(_M_cp));
2970
2971 // Make sure the last cumulative probability is one.
2972 _M_cp[_M_cp.size() - 1] = 1.0;
2973
2974 for (size_t __k = 0; __k < _M_den.size(); ++__k)
2975 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2976 }
2977
2978 template<typename _RealType>
2979 template<typename _InputIteratorB, typename _InputIteratorW>
2980 piecewise_constant_distribution<_RealType>::param_type::
2981 param_type(_InputIteratorB __bbegin,
2982 _InputIteratorB __bend,
2983 _InputIteratorW __wbegin)
2984 : _M_int(), _M_den(), _M_cp()
2985 {
2986 if (__bbegin != __bend)
2987 {
2988 for (;;)
2989 {
2990 _M_int.push_back(*__bbegin);
2991 ++__bbegin;
2992 if (__bbegin == __bend)
2993 break;
2994
2995 _M_den.push_back(*__wbegin);
2996 ++__wbegin;
2997 }
2998 }
2999
3000 _M_initialize();
3001 }
3002
3003 template<typename _RealType>
3004 template<typename _Func>
3005 piecewise_constant_distribution<_RealType>::param_type::
3006 param_type(initializer_list<_RealType> __bl, _Func __fw)
3007 : _M_int(), _M_den(), _M_cp()
3008 {
3009 _M_int.reserve(__bl.size());
3010 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3011 _M_int.push_back(*__biter);
3012
3013 _M_den.reserve(_M_int.size() - 1);
3014 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3015 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
3016
3017 _M_initialize();
3018 }
3019
3020 template<typename _RealType>
3021 template<typename _Func>
3022 piecewise_constant_distribution<_RealType>::param_type::
3023 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3024 : _M_int(), _M_den(), _M_cp()
3025 {
3026 const size_t __n = __nw == 0 ? 1 : __nw;
3027 const _RealType __delta = (__xmax - __xmin) / __n;
3028
3029 _M_int.reserve(__n + 1);
3030 for (size_t __k = 0; __k <= __nw; ++__k)
3031 _M_int.push_back(__xmin + __k * __delta);
3032
3033 _M_den.reserve(__n);
3034 for (size_t __k = 0; __k < __nw; ++__k)
3035 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
3036
3037 _M_initialize();
3038 }
3039
3040 template<typename _RealType>
3041 template<typename _UniformRandomNumberGenerator>
3042 typename piecewise_constant_distribution<_RealType>::result_type
3043 piecewise_constant_distribution<_RealType>::
3044 operator()(_UniformRandomNumberGenerator& __urng,
3045 const param_type& __param)
3046 {
3047 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3048 __aurng(__urng);
3049
3050 const double __p = __aurng();
3051 if (__param._M_cp.empty())
3052 return __p;
3053
3054 auto __pos = std::lower_bound(__param._M_cp.begin(),
3055 __param._M_cp.end(), __p);
3056 const size_t __i = __pos - __param._M_cp.begin();
3057
3058 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3059
3060 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
3061 }
3062
3063 template<typename _RealType>
3064 template<typename _ForwardIterator,
3065 typename _UniformRandomNumberGenerator>
3066 void
3067 piecewise_constant_distribution<_RealType>::
3068 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3069 _UniformRandomNumberGenerator& __urng,
3070 const param_type& __param)
3071 {
3072 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3073 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3074 __aurng(__urng);
3075
3076 if (__param._M_cp.empty())
3077 {
3078 while (__f != __t)
3079 *__f++ = __aurng();
3080 return;
3081 }
3082
3083 while (__f != __t)
3084 {
3085 const double __p = __aurng();
3086
3087 auto __pos = std::lower_bound(__param._M_cp.begin(),
3088 __param._M_cp.end(), __p);
3089 const size_t __i = __pos - __param._M_cp.begin();
3090
3091 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3092
3093 *__f++ = (__param._M_int[__i]
3094 + (__p - __pref) / __param._M_den[__i]);
3095 }
3096 }
3097
3098 template<typename _RealType, typename _CharT, typename _Traits>
3099 std::basic_ostream<_CharT, _Traits>&
3100 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3101 const piecewise_constant_distribution<_RealType>& __x)
3102 {
3103 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3104
3105 const typename __ios_base::fmtflags __flags = __os.flags();
3106 const _CharT __fill = __os.fill();
3107 const std::streamsize __precision = __os.precision();
3108 const _CharT __space = __os.widen(' ');
3109 __os.flags(__ios_base::scientific | __ios_base::left);
3110 __os.fill(__space);
3111 __os.precision(std::numeric_limits<_RealType>::max_digits10);
3112
3113 std::vector<_RealType> __int = __x.intervals();
3114 __os << __int.size() - 1;
3115
3116 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3117 __os << __space << *__xit;
3118
3119 std::vector<double> __den = __x.densities();
3120 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3121 __os << __space << *__dit;
3122
3123 __os.flags(__flags);
3124 __os.fill(__fill);
3125 __os.precision(__precision);
3126 return __os;
3127 }
3128
3129 template<typename _RealType, typename _CharT, typename _Traits>
3130 std::basic_istream<_CharT, _Traits>&
3131 operator>>(std::basic_istream<_CharT, _Traits>& __is,
3132 piecewise_constant_distribution<_RealType>& __x)
3133 {
3134 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3135
3136 const typename __ios_base::fmtflags __flags = __is.flags();
3137 __is.flags(__ios_base::dec | __ios_base::skipws);
3138
3139 size_t __n;
3140 if (__is >> __n)
3141 {
3142 std::vector<_RealType> __int_vec;
3143 if (__detail::__extract_params(__is, __int_vec, __n + 1))
3144 {
3145 std::vector<double> __den_vec;
3146 if (__detail::__extract_params(__is, __den_vec, __n))
3147 {
3148 __x.param({ __int_vec.begin(), __int_vec.end(),
3149 __den_vec.begin() });
3150 }
3151 }
3152 }
3153
3154 __is.flags(__flags);
3155 return __is;
3156 }
3157
3158
3159 template<typename _RealType>
3160 void
3161 piecewise_linear_distribution<_RealType>::param_type::
3162 _M_initialize()
3163 {
3164 if (_M_int.size() < 2
3165 || (_M_int.size() == 2
3166 && _M_int[0] == _RealType(0)
3167 && _M_int[1] == _RealType(1)
3168 && _M_den[0] == _M_den[1]))
3169 {
3170 _M_int.clear();
3171 _M_den.clear();
3172 return;
3173 }
3174
3175 double __sum = 0.0;
3176 _M_cp.reserve(_M_int.size() - 1);
3177 _M_m.reserve(_M_int.size() - 1);
3178 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3179 {
3180 const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3181 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3182 _M_cp.push_back(__sum);
3183 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3184 }
3185 __glibcxx_assert(__sum > 0);
3186
3187 // Now normalize the densities...
3188 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3189 __sum);
3190 // ... and partial sums...
3191 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3192 // ... and slopes.
3193 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3194
3195 // Make sure the last cumulative probablility is one.
3196 _M_cp[_M_cp.size() - 1] = 1.0;
3197 }
3198
3199 template<typename _RealType>
3200 template<typename _InputIteratorB, typename _InputIteratorW>
3201 piecewise_linear_distribution<_RealType>::param_type::
3202 param_type(_InputIteratorB __bbegin,
3203 _InputIteratorB __bend,
3204 _InputIteratorW __wbegin)
3205 : _M_int(), _M_den(), _M_cp(), _M_m()
3206 {
3207 for (; __bbegin != __bend; ++__bbegin, (void) ++__wbegin)
3208 {
3209 _M_int.push_back(*__bbegin);
3210 _M_den.push_back(*__wbegin);
3211 }
3212
3213 _M_initialize();
3214 }
3215
3216 template<typename _RealType>
3217 template<typename _Func>
3218 piecewise_linear_distribution<_RealType>::param_type::
3219 param_type(initializer_list<_RealType> __bl, _Func __fw)
3220 : _M_int(), _M_den(), _M_cp(), _M_m()
3221 {
3222 _M_int.reserve(__bl.size());
3223 _M_den.reserve(__bl.size());
3224 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3225 {
3226 _M_int.push_back(*__biter);
3227 _M_den.push_back(__fw(*__biter));
3228 }
3229
3230 _M_initialize();
3231 }
3232
3233 template<typename _RealType>
3234 template<typename _Func>
3235 piecewise_linear_distribution<_RealType>::param_type::
3236 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3237 : _M_int(), _M_den(), _M_cp(), _M_m()
3238 {
3239 const size_t __n = __nw == 0 ? 1 : __nw;
3240 const _RealType __delta = (__xmax - __xmin) / __n;
3241
3242 _M_int.reserve(__n + 1);
3243 _M_den.reserve(__n + 1);
3244 for (size_t __k = 0; __k <= __nw; ++__k)
3245 {
3246 _M_int.push_back(__xmin + __k * __delta);
3247 _M_den.push_back(__fw(_M_int[__k] + __delta));
3248 }
3249
3250 _M_initialize();
3251 }
3252
3253 template<typename _RealType>
3254 template<typename _UniformRandomNumberGenerator>
3255 typename piecewise_linear_distribution<_RealType>::result_type
3256 piecewise_linear_distribution<_RealType>::
3257 operator()(_UniformRandomNumberGenerator& __urng,
3258 const param_type& __param)
3259 {
3260 __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3261 __aurng(__urng);
3262
3263 const double __p = __aurng();
3264 if (__param._M_cp.empty())
3265 return __p;
3266
3267 auto __pos = std::lower_bound(__param._M_cp.begin(),
3268 __param._M_cp.end(), __p);
3269 const size_t __i = __pos - __param._M_cp.begin();
3270
3271 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3272
3273 const double __a = 0.5 * __param._M_m[__i];
3274 const double __b = __param._M_den[__i];
3275 const double __cm = __p - __pref;
3276
3277 _RealType __x = __param._M_int[__i];
3278 if (__a == 0)
3279 __x += __cm / __b;
3280 else
3281 {
3282 const double __d = __b * __b + 4.0 * __a * __cm;
3283 __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3284 }
3285
3286 return __x;
3287 }
3288
3289 template<typename _RealType>
3290 template<typename _ForwardIterator,
3291 typename _UniformRandomNumberGenerator>
3292 void
3293 piecewise_linear_distribution<_RealType>::
3294 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3295 _UniformRandomNumberGenerator& __urng,
3296 const param_type& __param)
3297 {
3298 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3299 // We could duplicate everything from operator()...
3300 while (__f != __t)
3301 *__f++ = this->operator()(__urng, __param);
3302 }
3303
3304 template<typename _RealType, typename _CharT, typename _Traits>
3305 std::basic_ostream<_CharT, _Traits>&
3306 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3307 const piecewise_linear_distribution<_RealType>& __x)
3308 {
3309 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3310
3311 const typename __ios_base::fmtflags __flags = __os.flags();
3312 const _CharT __fill = __os.fill();
3313 const std::streamsize __precision = __os.precision();
3314 const _CharT __space = __os.widen(' ');
3315 __os.flags(__ios_base::scientific | __ios_base::left);
3316 __os.fill(__space);
3317 __os.precision(std::numeric_limits<_RealType>::max_digits10);
3318
3319 std::vector<_RealType> __int = __x.intervals();
3320 __os << __int.size() - 1;
3321
3322 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3323 __os << __space << *__xit;
3324
3325 std::vector<double> __den = __x.densities();
3326 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3327 __os << __space << *__dit;
3328
3329 __os.flags(__flags);
3330 __os.fill(__fill);
3331 __os.precision(__precision);
3332 return __os;
3333 }
3334
3335 template<typename _RealType, typename _CharT, typename _Traits>
3336 std::basic_istream<_CharT, _Traits>&
3337 operator>>(std::basic_istream<_CharT, _Traits>& __is,
3338 piecewise_linear_distribution<_RealType>& __x)
3339 {
3340 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3341
3342 const typename __ios_base::fmtflags __flags = __is.flags();
3343 __is.flags(__ios_base::dec | __ios_base::skipws);
3344
3345 size_t __n;
3346 if (__is >> __n)
3347 {
3348 vector<_RealType> __int_vec;
3349 if (__detail::__extract_params(__is, __int_vec, __n + 1))
3350 {
3351 vector<double> __den_vec;
3352 if (__detail::__extract_params(__is, __den_vec, __n + 1))
3353 {
3354 __x.param({ __int_vec.begin(), __int_vec.end(),
3355 __den_vec.begin() });
3356 }
3357 }
3358 }
3359 __is.flags(__flags);
3360 return __is;
3361 }
3362
3363
3364 template<typename _IntType, typename>
3365 seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3366 {
3367 _M_v.reserve(__il.size());
3368 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3369 _M_v.push_back(__detail::__mod<result_type,
3370 __detail::_Shift<result_type, 32>::__value>(*__iter));
3371 }
3372
3373 template<typename _InputIterator>
3374 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3375 {
3376#pragma GCC diagnostic push
3377#pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3378 if constexpr (__is_random_access_iter<_InputIterator>::value)
3379 _M_v.reserve(std::distance(__begin, __end));
3380#pragma GCC diagnostic pop
3381
3382 for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3383 _M_v.push_back(__detail::__mod<result_type,
3384 __detail::_Shift<result_type, 32>::__value>(*__iter));
3385 }
3386
3387 template<typename _RandomAccessIterator>
3388 void
3389 seed_seq::generate(_RandomAccessIterator __begin,
3390 _RandomAccessIterator __end)
3391 {
3392 typedef typename iterator_traits<_RandomAccessIterator>::value_type
3393 _Type;
3394
3395 if (__begin == __end)
3396 return;
3397
3398 std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3399
3400 const size_t __n = __end - __begin;
3401 const size_t __s = _M_v.size();
3402 const size_t __t = (__n >= 623) ? 11
3403 : (__n >= 68) ? 7
3404 : (__n >= 39) ? 5
3405 : (__n >= 7) ? 3
3406 : (__n - 1) / 2;
3407 const size_t __p = (__n - __t) / 2;
3408 const size_t __q = __p + __t;
3409 const size_t __m = std::max(size_t(__s + 1), __n);
3410
3411#ifndef __UINT32_TYPE__
3412 struct _Up
3413 {
3414 _Up(uint_least32_t v) : _M_v(v & 0xffffffffu) { }
3415
3416 operator uint_least32_t() const { return _M_v; }
3417
3418 uint_least32_t _M_v;
3419 };
3420 using uint32_t = _Up;
3421#endif
3422
3423 // k == 0, every element in [begin,end) equals 0x8b8b8b8bu
3424 {
3425 uint32_t __r1 = 1371501266u;
3426 uint32_t __r2 = __r1 + __s;
3427 __begin[__p] += __r1;
3428 __begin[__q] = (uint32_t)__begin[__q] + __r2;
3429 __begin[0] = __r2;
3430 }
3431
3432 for (size_t __k = 1; __k <= __s; ++__k)
3433 {
3434 const size_t __kn = __k % __n;
3435 const size_t __kpn = (__k + __p) % __n;
3436 const size_t __kqn = (__k + __q) % __n;
3437 uint32_t __arg = (__begin[__kn]
3438 ^ __begin[__kpn]
3439 ^ __begin[(__k - 1) % __n]);
3440 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3441 uint32_t __r2 = __r1 + (uint32_t)__kn + _M_v[__k - 1];
3442 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3443 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3444 __begin[__kn] = __r2;
3445 }
3446
3447 for (size_t __k = __s + 1; __k < __m; ++__k)
3448 {
3449 const size_t __kn = __k % __n;
3450 const size_t __kpn = (__k + __p) % __n;
3451 const size_t __kqn = (__k + __q) % __n;
3452 uint32_t __arg = (__begin[__kn]
3453 ^ __begin[__kpn]
3454 ^ __begin[(__k - 1) % __n]);
3455 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3456 uint32_t __r2 = __r1 + (uint32_t)__kn;
3457 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3458 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3459 __begin[__kn] = __r2;
3460 }
3461
3462 for (size_t __k = __m; __k < __m + __n; ++__k)
3463 {
3464 const size_t __kn = __k % __n;
3465 const size_t __kpn = (__k + __p) % __n;
3466 const size_t __kqn = (__k + __q) % __n;
3467 uint32_t __arg = (__begin[__kn]
3468 + __begin[__kpn]
3469 + __begin[(__k - 1) % __n]);
3470 uint32_t __r3 = 1566083941u * (__arg ^ (__arg >> 27));
3471 uint32_t __r4 = __r3 - __kn;
3472 __begin[__kpn] ^= __r3;
3473 __begin[__kqn] ^= __r4;
3474 __begin[__kn] = __r4;
3475 }
3476 }
3477
3478// [rand.util.canonical]
3479// generate_canonical(RNG&)
3480
3481#ifndef _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3482
3483#pragma GCC diagnostic push
3484#pragma GCC diagnostic ignored "-Wc++14-extensions" // for variable templates
3485#pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3486
3487 // __generate_canonical_pow2 is used when Urbg::max()-Urbg::min() is
3488 // a power of two less 1. It works by calling urng() as many times as
3489 // needed to fill the target mantissa, accumulating entropy into an
3490 // integer value, converting that to the float type, and then dividing
3491 // by the range of the integer value (a constexpr power of 2,
3492 // so only adjusts the exponent) to produce a result in [0..1].
3493 // In case of an exact 1.0 result, we re-try.
3494 //
3495 // It needs to work even when the integer type used is only as big
3496 // as the float mantissa, such as uint64_t for long double. So,
3497 // commented-out assignments represent computations the Standard
3498 // prescribes but cannot be performed, or are not used. Names are
3499 // chosen to match the description in the Standard.
3500 //
3501 // When the result is close to zero, the strict Standard-prescribed
3502 // calculation may leave more low-order zeros in the mantissa than
3503 // is usually necessary. When spare entropy has been extracted, as
3504 // is usual for float and double, some or all of the spare entropy
3505 // can commonly be pulled into the result for better randomness.
3506 // Defining _GLIBCXX_GENERATE_CANONICAL_STRICT discards it instead.
3507 //
3508 // When k calls to urng() yield more bits of entropy, log2_Rk_max,
3509 // than fit into UInt, we discard some of it by overflowing, which
3510 // is OK. On converting the integer representation of the sample
3511 // to the float value, we must divide out the (possibly-truncated)
3512 // size log2_Rk.
3513 //
3514 // This implementation works with std::bfloat16, which can exactly
3515 // represent 2^32, but not with std::float16_t, limited to 2^15.
3516
3517 template<typename _RealT, size_t __d, typename _Urbg>
3518 _RealT
3519 __generate_canonical_pow2(_Urbg& __urng)
3520 {
3521 using _UInt = typename __detail::_Select_uint_least_t<__d>::type;
3522
3523 // Parameter __d is the actual target number of bits.
3524 // Commented-out assignments below are of values specified in
3525 // the Standard, but not used here for reasons noted.
3526 // r = 2; // Redundant, we only support radix 2.
3527 using _Rng = decltype(_Urbg::max());
3528 const _Rng __rng_range_less_1 = _Urbg::max() - _Urbg::min();
3529 // R = _UInt(__rng_range_less_1) + 1; // May wrap to 0.
3530 const auto __log2_R = __builtin_popcountg(__rng_range_less_1);
3531 const auto __log2_uint_max = sizeof(_UInt) * __CHAR_BIT__;
3532 // rd = _UInt(1) << __d; // Could overflow, UB.
3533 const unsigned __k = (__d + __log2_R - 1) / __log2_R;
3534 const unsigned __log2_Rk_max = __k * __log2_R;
3535 const unsigned __log2_Rk = // Bits of entropy actually obtained:
3536 __log2_uint_max < __log2_Rk_max ? __log2_uint_max : __log2_Rk_max;
3537 // Rk = _UInt(1) << __log2_Rk; // Likely overflows, UB.
3538 _GLIBCXX14_CONSTEXPR const _RealT __Rk
3539 = _RealT(_UInt(1) << (__log2_Rk - 1)) * _RealT(2.0);
3540#if defined(_GLIBCXX_GENERATE_CANONICAL_STRICT)
3541 const unsigned __log2_x = __log2_Rk - __d; // # of spare entropy bits.
3542#else
3543 const unsigned __log2_x = 0;
3544#endif
3545 _GLIBCXX14_CONSTEXPR const _UInt __x = _UInt(1) << __log2_x;
3546 _GLIBCXX14_CONSTEXPR const _RealT __rd = __Rk / _RealT(__x);
3547 // xrd = __x << __d; // Could overflow.
3548
3549 while (true)
3550 {
3551 _UInt __sum = _UInt(__urng() - _Urbg::min());
3552 for (unsigned __i = __k - 1, __shift = 0; __i > 0; --__i)
3553 {
3554 __shift += __log2_R;
3555 __sum |= _UInt(__urng() - _Urbg::min()) << __shift;
3556 }
3557 const _RealT __ret = _RealT(__sum >> __log2_x) / _RealT(__rd);
3558 if (__ret < _RealT(1.0))
3559 return __ret;
3560 }
3561 }
3562
3563
3564 template<typename _UInt>
3565 struct __gen_canon_log_res
3566 {
3567 unsigned __floor_log;
3568 _UInt __floor_pow;
3569
3570 constexpr __gen_canon_log_res
3571 update(_UInt __base) const
3572 { return {__floor_log + 1, __floor_pow * __base}; }
3573 };
3574
3575
3576 template <typename _UInt1, typename _UInt2,
3577 typename _UComm = __conditional_t<(sizeof(_UInt2) > sizeof(_UInt1)),
3578 _UInt2, _UInt1>>
3579 constexpr __gen_canon_log_res<_UInt1>
3580 __gen_canon_log(_UInt1 __val, _UInt2 __base)
3581 {
3582#if __cplusplus >= 201402L
3583 __gen_canon_log_res<_UInt1> __res{0, _UInt1(1)};
3584 if (_UComm(__base) > _UComm(__val))
3585 return __res;
3586
3587 const _UInt1 __base1(__base);
3588 do
3589 {
3590 __val /= __base1;
3591 __res = __res.update(__base1);
3592 }
3593 while (__val >= __base1);
3594 return __res;
3595#else
3596 return (_UComm(__val) >= _UComm(__base))
3597 ? __gen_canon_log(__val / _UInt1(__base), _UInt1(__base))
3598 .update(_UInt1(__base))
3599 : __gen_canon_log_res<_UInt1>{0, _UInt1(1)};
3600#endif
3601 }
3602
3603 // This version must be used when the range of possible RNG results,
3604 // Urbg::max()-Urbg::min(), is not a power of two less one. The UInt
3605 // type passed must be big enough to represent Rk, R^k, a power of R
3606 // (the range of values produced by the rng) up to twice the length
3607 // of the mantissa.
3608
3609 template<typename _RealT, size_t __d, typename _Urbg>
3610 _RealT
3611 __generate_canonical_any(_Urbg& __urng)
3612 {
3613 // Names below are chosen to match the description in the Standard.
3614 // Parameter d is the actual target number of bits.
3615#if (__cplusplus >= 201402L) || defined(__SIZEOF_INT128__)
3616# define _GLIBCXX_GEN_CANON_CONST constexpr
3617#else
3618# define _GLIBCXX_GEN_CANON_CONST const
3619#endif
3620
3621 using _UIntR = typename make_unsigned<decltype(_Urbg::max())>::type;
3622 // Cannot overflow, as _Urbg::max() - _Urbg::min() is not power of
3623 // two minus one
3624 constexpr _UIntR __R = _UIntR(_Urbg::max() - _Urbg::min()) + 1;
3625 constexpr unsigned __log2R
3626 = sizeof(_UIntR) * __CHAR_BIT__ - __builtin_clzg(__R) - 1;
3627 // We overstimate number of required bits, by computing
3628 // r such that l * log2(R) >= d, so:
3629 // R^l >= (2 ^ log2(R)) ^ l == 2 ^ (log2(r) * l) >= 2^d
3630 // And then requiring l * bit_width(R) bits.
3631 constexpr unsigned __l = (__d + __log2R - 1) / __log2R;
3632 constexpr unsigned __bits = (__log2R + 1) * __l;
3633 using _UInt = typename __detail::_Select_uint_least_t<__bits>::type;
3634
3635 _GLIBCXX_GEN_CANON_CONST _UInt __rd = _UInt(1) << __d;
3636 _GLIBCXX_GEN_CANON_CONST auto __logRrd = __gen_canon_log(__rd, __R);
3637 _GLIBCXX_GEN_CANON_CONST unsigned __k
3638 = __logRrd.__floor_log + (__rd > __logRrd.__floor_pow);
3639
3640 _GLIBCXX_GEN_CANON_CONST _UInt __Rk
3641 = (__k > __logRrd.__floor_log)
3642 ? _UInt(__logRrd.__floor_pow) * _UInt(__R)
3643 : _UInt(__logRrd.__floor_pow);
3644 _GLIBCXX_GEN_CANON_CONST _UInt __x = __Rk / __rd;
3645
3646 while (true)
3647 {
3648 _UInt __Ri{1};
3649 _UInt __sum(__urng() - _Urbg::min());
3650 for (int __i = __k - 1; __i > 0; --__i)
3651 {
3652 __Ri *= _UInt(__R);
3653 __sum += _UInt(__urng() - _Urbg::min()) * __Ri;
3654 }
3655 const _RealT __ret = _RealT(__sum / __x) / _RealT(__rd);
3656 if (__ret < _RealT(1.0))
3657 return __ret;
3658 }
3659#undef _GLIBCXX_GEN_CANON_CONST
3660 }
3661
3662#if !defined(_GLIBCXX_GENERATE_CANONICAL_STRICT)
3663 template <typename _Tp>
3664 const bool __is_rand_dist_float_v = is_floating_point<_Tp>::value;
3665#else
3666 template <typename _Tp> const bool __is_rand_dist_float_v = false;
3667 template <> const bool __is_rand_dist_float_v<float> = true;
3668 template <> const bool __is_rand_dist_float_v<double> = true;
3669 template <> const bool __is_rand_dist_float_v<long double> = true;
3670#endif
3671
3672 // Note, this works even when (__range + 1) overflows:
3673 template <typename _Rng>
3674 constexpr bool __is_power_of_2_less_1(_Rng __range)
3675 { return ((__range + 1) & __range) == 0; };
3676
3677_GLIBCXX_BEGIN_INLINE_ABI_NAMESPACE(_V2)
3678 /** Produce a random floating-point value in the range [0..1)
3679 *
3680 * The result of `std::generate_canonical<RealT,digits>(urng)` is a
3681 * random floating-point value of type `RealT` in the range [0..1),
3682 * using entropy provided by the uniform random bit generator `urng`.
3683 * A value for `digits` may be passed to limit the precision of the
3684 * result to so many bits, but normally `-1u` is passed to get the
3685 * native precision of `RealT`. As many `urng()` calls are made as
3686 * needed to obtain the required entropy. On rare occasions, more
3687 * `urng()` calls are used. It is fastest when the value of
3688 * `Urbg::max()` is a power of two less one, such as from a
3689 * `std::philox4x32` (for `float`) or `philox4x64` (for `double`).
3690 *
3691 * @since C++11
3692 */
3693 template<typename _RealT, size_t __digits,
3694 typename _Urbg>
3695 _RealT
3696 generate_canonical(_Urbg& __urng)
3697 {
3698#ifdef __glibcxx_concepts
3699 static_assert(uniform_random_bit_generator<_Urbg>);
3700#endif
3701 static_assert(__is_rand_dist_float_v<_RealT>,
3702 "template argument must be a floating point type");
3703 static_assert(__digits != 0 && _Urbg::max() > _Urbg::min(),
3704 "random samples with 0 bits are not meaningful");
3705 static_assert(std::numeric_limits<_RealT>::radix == 2,
3706 "only base-2 float types are supported");
3707#if defined(__STDCPP_FLOAT16_T__)
3708 static_assert(! is_same_v<_RealT, _Float16>,
3709 "float16_t type is not supported, consider using bfloat16_t");
3710#endif
3711
3712 const unsigned __d_max = std::numeric_limits<_RealT>::digits;
3713 const unsigned __d = __digits > __d_max ? __d_max : __digits;
3714
3715 // If the RNG range is a power of 2 less 1, the float type mantissa
3716 // is enough bits. If not, we need more.
3717 if constexpr (__is_power_of_2_less_1(_Urbg::max() - _Urbg::min()))
3718 return __generate_canonical_pow2<_RealT, __d>(__urng);
3719 else // Need up to 2x bits.
3720 return __generate_canonical_any<_RealT, __d>(__urng);
3721 }
3722_GLIBCXX_END_INLINE_ABI_NAMESPACE(_V2)
3723
3724#pragma GCC diagnostic pop
3725
3726#else // _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3727
3728 // This is the pre-P0952 definition, to reproduce old results.
3729
3730 template<typename _RealType, size_t __bits,
3731 typename _UniformRandomNumberGenerator>
3732 _RealType
3733 generate_canonical(_UniformRandomNumberGenerator& __urng)
3734 {
3735 static_assert(std::is_floating_point<_RealType>::value,
3736 "template argument must be a floating point type");
3737
3738 const size_t __b
3739 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3740 __bits);
3741 const long double __r = static_cast<long double>(__urng.max())
3742 - static_cast<long double>(__urng.min()) + 1.0L;
3743 const size_t __log2r = std::log(__r) / std::log(2.0L);
3744 const size_t __m = std::max<size_t>(1UL,
3745 (__b + __log2r - 1UL) / __log2r);
3746 _RealType __ret;
3747 _RealType __sum = _RealType(0);
3748 _RealType __tmp = _RealType(1);
3749 for (size_t __k = __m; __k != 0; --__k)
3750 {
3751 __sum += _RealType(__urng() - __urng.min()) * __tmp;
3752 __tmp *= __r;
3753 }
3754 __ret = __sum / __tmp;
3755 if (__builtin_expect(__ret >= _RealType(1), 0))
3756 {
3757# if _GLIBCXX_USE_C99_MATH_FUNCS
3758 __ret = std::nextafter(_RealType(1), _RealType(0));
3759# else
3760 __ret = _RealType(1)
3761 - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3762# endif
3763 }
3764 return __ret;
3765 }
3766
3767#endif // _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3768
3769_GLIBCXX_END_NAMESPACE_VERSION
3770} // namespace
3771
3772#endif
std::log
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
Definition complex:1162
std::tan
complex< _Tp > tan(const complex< _Tp > &)
Return complex tangent of z.
Definition complex:1298
std::abs
_Tp abs(const complex< _Tp > &)
Return magnitude of z.
Definition complex:968
std::exp
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.
Definition complex:1135
std::pow
complex< _Tp > pow(const complex< _Tp > &, int)
Return x to the y'th power.
Definition complex:1357
std::sqrt
complex< _Tp > sqrt(const complex< _Tp > &)
Return complex square root of z.
Definition complex:1271
std::max
constexpr const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
Definition stl_algobase.h:256
std::min
constexpr const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition stl_algobase.h:232
std::generate_canonical
_RealType generate_canonical(_UniformRandomNumberGenerator &__g)
A function template for converting the output of a (integral) uniform random number generator to a fl...
std::back_inserter
constexpr back_insert_iterator< _Container > back_inserter(_Container &__x)
Definition bits/stl_iterator.h:764
std::accumulate
constexpr _Tp accumulate(_InputIterator __first, _InputIterator __last, _Tp __init)
Accumulate values in a range.
Definition stl_numeric.h:133
std::partial_sum
constexpr _OutputIterator partial_sum(_InputIterator __first, _InputIterator __last, _OutputIterator __result)
Return list of partial sums.
Definition stl_numeric.h:255
std
ISO C++ entities toplevel namespace is std.
std::streamsize
ptrdiff_t streamsize
Integral type for I/O operation counts and buffer sizes.
Definition postypes.h:73
std::distance
constexpr iterator_traits< _InputIterator >::difference_type distance(_InputIterator __first, _InputIterator __last)
A generalization of pointer arithmetic.
Definition stl_iterator_base_funcs.h:172
std::__lg
constexpr _Tp __lg(_Tp __n)
This is a helper function for the sort routines and for random.tcc.
Definition stl_algobase.h:1558
std::operator>>
std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, bitset< _Nb > &__x)
Global I/O operators for bitsets.
Definition bitset:1658
std::operator<<
std::basic_ostream< _CharT, _Traits > & operator<<(std::basic_ostream< _CharT, _Traits > &__os, const bitset< _Nb > &__x)
Global I/O operators for bitsets.
Definition bitset:1754
std::__detail
Implementation details not part of the namespace std interface.
__gnu_debug::__base
constexpr _Iterator __base(_Iterator __it)
Definition helper_functions.h:339
std::initializer_list
initializer_list
Definition initializer_list:48
std::array
A standard container for storing a fixed size sequence of elements.
Definition array:103
std::basic_ios::widen
char_type widen(char __c) const
Widens characters.
Definition basic_ios.h:465
std::basic_ios::fill
char_type fill() const
Retrieves the empty character.
Definition basic_ios.h:388
std::basic_istream
Template class basic_istream.
Definition istream:67
std::basic_ostream
Template class basic_ostream.
Definition ostream.h:67
std::__numeric_limits_base::is_integer
static constexpr bool is_integer
Definition limits:233
std::__numeric_limits_base::max_digits10
static constexpr int max_digits10
Definition limits:226
std::__numeric_limits_base::digits
static constexpr int digits
Definition limits:218
std::__numeric_limits_base::is_signed
static constexpr bool is_signed
Definition limits:230
std::__numeric_limits_base::radix
static constexpr int radix
Definition limits:242
std::numeric_limits::max
static constexpr _Tp max() noexcept
Definition limits:328
std::numeric_limits::epsilon
static constexpr _Tp epsilon() noexcept
Definition limits:340
std::is_floating_point
is_floating_point
Definition type_traits:601
std::common_type
common_type
Definition type_traits:2573
std::ios_base::precision
streamsize precision() const
Flags access.
Definition ios_base.h:765
std::ios_base::flags
fmtflags flags() const
Access to format flags.
Definition ios_base.h:694
std::linear_congruential_engine
A model of a linear congruential random number generator.
Definition random.h:704
std::linear_congruential_engine::multiplier
static constexpr result_type multiplier
Definition random.h:719
std::linear_congruential_engine::modulus
static constexpr result_type modulus
Definition random.h:723
std::linear_congruential_engine::seed
void seed(result_type __s=default_seed)
Reseeds the linear_congruential_engine random number generator engine sequence to the seed __s.
Definition bits/random.tcc:119
std::linear_congruential_engine::increment
static constexpr result_type increment
Definition random.h:721
std::subtract_with_carry_engine
The Marsaglia-Zaman generator.
Definition random.h:1151
std::subtract_with_carry_engine::seed
void seed(result_type __sd=0u)
Seeds the initial state of the random number generator.
Definition bits/random.tcc:542
std::subtract_with_carry_engine::operator()
result_type operator()()
Gets the next random number in the sequence.
Definition bits/random.tcc:601
std::discard_block_engine::operator()
result_type operator()()
Gets the next value in the generated random number sequence.
Definition bits/random.tcc:690
std::independent_bits_engine::operator()
result_type operator()()
Gets the next value in the generated random number sequence.
Definition bits/random.tcc:745
std::shuffle_order_engine
Produces random numbers by reordering random numbers from some base engine.
Definition random.h:1801
std::shuffle_order_engine::base
const _RandomNumberEngine & base() const noexcept
Definition random.h:1907
std::shuffle_order_engine::result_type
_RandomNumberEngine::result_type result_type
Definition random.h:1807
std::philox_engine
A discrete pseudorandom number generator with weak cryptographic properties.
Definition random.h:2060
std::uniform_real_distribution
Uniform continuous distribution for random numbers.
Definition random.h:2493
std::uniform_real_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:2582
std::normal_distribution
A normal continuous distribution for random numbers.
Definition random.h:2730
std::normal_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:2849
std::gamma_distribution
A gamma continuous distribution for random numbers.
Definition random.h:3182
std::gamma_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:3311
std::cauchy_distribution
A cauchy_distribution random number distribution.
Definition random.h:3654
std::cauchy_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:3731
std::cauchy_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:3761
std::fisher_f_distribution
A fisher_f_distribution random number distribution.
Definition random.h:3869
std::binomial_distribution
A discrete binomial random number distribution.
Definition random.h:4565
std::binomial_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:4693
std::geometric_distribution
A discrete geometric random number distribution.
Definition random.h:4811
std::geometric_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:4922
std::geometric_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:4892
std::negative_binomial_distribution
A negative_binomial_distribution random number distribution.
Definition random.h:5028
std::negative_binomial_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition bits/random.tcc:1284
std::poisson_distribution
A discrete Poisson random number distribution.
Definition random.h:5265
std::poisson_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:5378
std::poisson_distribution::operator==
friend bool operator==(const poisson_distribution &__d1, const poisson_distribution &__d2)
Return true if two Poisson distributions have the same parameters and the sequences that would be gen...
Definition random.h:5414
std::poisson_distribution::operator<<
friend std::basic_ostream< _CharT, _Traits > & operator<<(std::basic_ostream< _CharT, _Traits > &__os, const std::poisson_distribution< _IntType1 > &__x)
Inserts a poisson_distribution random number distribution __x into the output stream __os.
std::poisson_distribution::operator>>
friend std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, std::poisson_distribution< _IntType1 > &__x)
Extracts a poisson_distribution random number distribution __x from the input stream __is.
std::exponential_distribution
An exponential continuous distribution for random numbers.
Definition random.h:5497
std::exponential_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:5577
std::weibull_distribution
A weibull_distribution random number distribution.
Definition random.h:5719
std::weibull_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:5799
std::weibull_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:5829
std::extreme_value_distribution
A extreme_value_distribution random number distribution.
Definition random.h:5936
std::extreme_value_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:6046
std::extreme_value_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:6016
std::vector
A standard container which offers fixed time access to individual elements in any order.
Definition stl_vector.h:461
std::vector::end
constexpr iterator end() noexcept
Definition stl_vector.h:1008
std::vector::begin
constexpr iterator begin() noexcept
Definition stl_vector.h:988
std::vector::size
constexpr size_type size() const noexcept
Definition stl_vector.h:1107
std::uniform_int_distribution
Uniform discrete distribution for random numbers. A discrete random distribution on the range with e...
Definition uniform_int_dist.h:89
std::uniform_int_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition uniform_int_dist.h:173
std::uniform_random_bit_generator
Requirements for a uniform random bit generator.
Definition uniform_int_dist.h:54