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 _RealType>
82 _OutputIterator
83 __normalize(_InputIterator __first, _InputIterator __last,
84 _OutputIterator __result, const _RealType& __factor)
85 {
86 for (; __first != __last; ++__first, (void) ++__result)
87 *__result = _RealType(_RealType(*__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_configure()
2949 {
2950 const _CalcType __sum = std::accumulate(_M_den.begin(), _M_den.end(),
2951 _CalcType(0));
2952 __glibcxx_assert(__sum > 0);
2953
2954 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2955 __sum);
2956
2957 _CalcType __psum(0);
2958 _M_cp.reserve(_M_den.size());
2959 for (_CalcType __den : _M_den)
2960 _M_cp.push_back(__psum += __den);
2961
2962 // Make sure the last cumulative probability is one.
2963 _M_cp[_M_cp.size() - 1] = _CalcType(1);
2964
2965 for (size_t __k = 0; __k < _M_den.size(); ++__k)
2966 _M_den[__k] = _CalcType(_CalcType(_M_den[__k])
2967 / (_M_int[__k + 1] - _M_int[__k]));
2968 }
2969
2970 template<typename _RealType>
2971 void
2972 piecewise_constant_distribution<_RealType>::param_type::
2973 _M_initialize2(const _RealType* __ints, _CalcType __den)
2974 {
2975 if (__ints[0] == _RealType(0) && __ints[1] == _RealType(1))
2976 return;
2977
2978 _M_int.reserve(2);
2979 _M_int.push_back(__ints[0]);
2980 _M_int.push_back(__ints[1]);
2981
2982 _M_den.reserve(1);
2983 _M_den.push_back(__den);
2984 _M_configure();
2985 }
2986
2987namespace __detail
2988{
2989 template<typename _InputIterator, typename _RealType>
2990 bool
2991 __load_first2(_InputIterator& __first, _InputIterator __last,
2992 _RealType* __out)
2993 {
2994 if (__first == __last)
2995 return false;
2996
2997 *__out = *__first;
2998 ++__first;
2999 if (__first == __last)
3000 return false;
3001
3002 ++__out;
3003 *__out = *__first;
3004 ++__first;
3005 return true;
3006 }
3007} // namespace __detail
3008
3009
3010 template<typename _RealType>
3011 template<typename _InputIteratorB, typename _InputIteratorW>
3012 piecewise_constant_distribution<_RealType>::param_type::
3013 param_type(_InputIteratorB __bbegin,
3014 _InputIteratorB __bend,
3015 _InputIteratorW __wbegin)
3016 : _M_int(), _M_den(), _M_cp()
3017 {
3018 _RealType __ints[2];
3019 if (!__detail::__load_first2(__bbegin, __bend, __ints))
3020 return;
3021
3022 if (__bbegin == __bend)
3023 {
3024 _M_initialize2(__ints, *__wbegin);
3025 return;
3026 }
3027
3028#if __glibcxx_concepts // C++ >= C++20
3029 if constexpr (sized_sentinel_for<_InputIteratorB, _InputIteratorB>
3030 || forward_iterator<_InputIteratorB>)
3031 _M_int.reserve(2 + size_t(ranges::distance(__bbegin, __bend)));
3032#else
3033#pragma GCC diagnostic push
3034#pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3035 if constexpr (is_convertible<__iter_category_t<_InputIteratorB>,
3036 forward_iterator_tag>::value)
3037 _M_int.reserve(2 + size_t(std::distance(__bbegin, __bend)));
3038#pragma GCC diagnostic pop
3039#endif
3040
3041 _M_int.push_back(__ints[0]);
3042 _M_int.push_back(__ints[1]);
3043 for (; __bbegin != __bend; ++__bbegin)
3044 _M_int.push_back(*__bbegin);
3045
3046 // _GLIBCXX_RESOLVE_LIB_DEFECTS
3047 // 4052. Bogus requirements for piecewise_linear_distribution
3048 _M_den.reserve(_M_int.size() - 1);
3049 for (size_t __k = 0; __k < _M_int.size() - 1; (void)++__k, ++__wbegin)
3050 _M_den.push_back(_CalcType(*__wbegin));
3051
3052 _M_configure();
3053 }
3054
3055 template<typename _RealType>
3056 template<typename _Func>
3057 piecewise_constant_distribution<_RealType>::param_type::
3058 param_type(initializer_list<_RealType> __bl, _Func __fw)
3059 : _M_int(), _M_den(), _M_cp()
3060 {
3061 if (__bl.size() < 2)
3062 return;
3063
3064 // _GLIBCXX_RESOLVE_LIB_DEFECTS
3065 // 4052. Bogus requirements for piecewise_linear_distribution
3066 auto __cfw = [&__fw](_RealType __n, _RealType __p)
3067 {
3068#ifdef _GLIBCXX_USE_OLD_PIECEWISE_DISTRIBUTIONS
3069 return __fw(0.5 * (__n + __p));
3070#else
3071 return _RealType(__fw(_RealType(0.5) * (__n + __p)));
3072#endif
3073 };
3074
3075 if (__bl.size() == 2)
3076 {
3077 const _RealType *__ints = __bl.begin();
3078 _CalcType __den = __cfw(__ints[1], __ints[0]);
3079 _M_initialize2(__ints, __den);
3080 return;
3081 }
3082
3083 _M_int = __bl;
3084 _M_den.reserve(_M_int.size() - 1);
3085 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3086 _M_den.push_back(__cfw(_M_int[__k + 1], _M_int[__k]));
3087
3088 _M_configure();
3089 }
3090
3091 template<typename _RealType>
3092 template<typename _Func>
3093 piecewise_constant_distribution<_RealType>::param_type::
3094 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3095 : _M_int(), _M_den(), _M_cp()
3096 {
3097 const size_t __n = __nw == 0 ? 1 : __nw;
3098 const _RealType __delta = (__xmax - __xmin) / __n;
3099 auto __cfw = [&__fw, __delta](_RealType __v)
3100 {
3101#ifdef _GLIBCXX_USE_OLD_PIECEWISE_DISTRIBUTIONS
3102 return __fw(__v + 0.5 * __delta);
3103#else
3104 return _RealType(__fw(__v + _RealType(0.5) * __delta));
3105#endif
3106 };
3107
3108 if (__n == 1)
3109 {
3110 _RealType __ints[2] = { __xmin, __xmin + __delta };
3111 _CalcType __den = __cfw(__xmin);
3112 _M_initialize2(__ints, __den);
3113 return;
3114 }
3115
3116 _M_int.reserve(__n + 1);
3117 for (size_t __k = 0; __k <= __nw; ++__k)
3118 _M_int.push_back(__xmin + __k * __delta);
3119
3120 _M_den.reserve(__n);
3121 for (size_t __k = 0; __k < __nw; ++__k)
3122 _M_den.push_back(__cfw(_M_int[__k]));
3123
3124 _M_configure();
3125 }
3126
3127 template<typename _RealType>
3128 template<typename _AdaptedUniformRandomNumberGenerator>
3129 typename piecewise_constant_distribution<_RealType>::result_type
3130 piecewise_constant_distribution<_RealType>::
3131 __generate_one(_AdaptedUniformRandomNumberGenerator& __aurng,
3132 const param_type& __param)
3133 {
3134 const _CalcType __p = __aurng();
3135 auto __pos = std::lower_bound(__param._M_cp.begin(),
3136 __param._M_cp.end(), __p);
3137 const size_t __i = __pos - __param._M_cp.begin();
3138 const _CalcType __pref = __i > 0 ? __param._M_cp[__i - 1] : _CalcType(0);
3139 return __param._M_int[__i] + (__p - __pref) / _CalcType(__param._M_den[__i]);
3140 }
3141
3142 template<typename _RealType>
3143 template<typename _UniformRandomNumberGenerator>
3144 typename piecewise_constant_distribution<_RealType>::result_type
3145 piecewise_constant_distribution<_RealType>::
3146 operator()(_UniformRandomNumberGenerator& __urng,
3147 const param_type& __param)
3148 {
3149 __detail::_Adaptor<_UniformRandomNumberGenerator, _CalcType>
3150 __aurng(__urng);
3151
3152 if (__param._M_cp.empty())
3153 return __aurng();
3154
3155 return __generate_one(__aurng, __param);
3156 }
3157
3158 template<typename _RealType>
3159 template<typename _ForwardIterator,
3160 typename _UniformRandomNumberGenerator>
3161 void
3162 piecewise_constant_distribution<_RealType>::
3163 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3164 _UniformRandomNumberGenerator& __urng,
3165 const param_type& __param)
3166 {
3167 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3168 __detail::_Adaptor<_UniformRandomNumberGenerator, _CalcType>
3169 __aurng(__urng);
3170
3171 if (__param._M_cp.empty())
3172 while (__f != __t)
3173 *__f++ = __aurng();
3174 else
3175 while (__f != __t)
3176 *__f++ = __generate_one(__urng, __param);
3177 }
3178
3179 template<typename _RealType, typename _CharT, typename _Traits>
3180 std::basic_ostream<_CharT, _Traits>&
3181 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3182 const piecewise_constant_distribution<_RealType>& __x)
3183 {
3184 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3185
3186 const typename __ios_base::fmtflags __flags = __os.flags();
3187 const _CharT __fill = __os.fill();
3188 const std::streamsize __precision = __os.precision();
3189 const _CharT __space = __os.widen(' ');
3190 __os.flags(__ios_base::scientific | __ios_base::left);
3191 __os.fill(__space);
3192 __os.precision(std::numeric_limits<_RealType>::max_digits10);
3193
3194 const auto& __int = __x._M_param._M_int;
3195 if (__int.empty())
3196 __os << size_t(1)
3197 << __space << _RealType(0) << __space << _RealType(1)
3198 << __space << _RealType(1);
3199 else
3200 {
3201 __os << __int.size() - 1;
3202 for (auto __xv : __int)
3203 __os << __space << __xv;
3204 for (auto __dv : __x._M_param._M_den)
3205 __os << __space << __dv;
3206 }
3207
3208 __os.flags(__flags);
3209 __os.fill(__fill);
3210 __os.precision(__precision);
3211 return __os;
3212 }
3213
3214 template<typename _RealType, typename _CharT, typename _Traits>
3215 std::basic_istream<_CharT, _Traits>&
3216 operator>>(std::basic_istream<_CharT, _Traits>& __is,
3217 piecewise_constant_distribution<_RealType>& __x)
3218 {
3219 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3220 using _StorageType
3221 = typename piecewise_constant_distribution<_RealType>::_StorageType;
3222
3223 const typename __ios_base::fmtflags __flags = __is.flags();
3224 __is.flags(__ios_base::dec | __ios_base::skipws);
3225
3226 size_t __n;
3227 if (__is >> __n)
3228 {
3229 std::vector<_RealType> __int_vec;
3230 if (__detail::__extract_params(__is, __int_vec, __n + 1))
3231 {
3232 std::vector<_StorageType> __den_vec;
3233 if (__detail::__extract_params(__is, __den_vec, __n))
3234 {
3235 __x.param({ __int_vec.begin(), __int_vec.end(),
3236 __den_vec.begin() });
3237 }
3238 }
3239 }
3240
3241 __is.flags(__flags);
3242 return __is;
3243 }
3244
3245
3246 template<typename _RealType>
3247 void
3248 piecewise_linear_distribution<_RealType>::param_type::
3249 _M_configure()
3250 {
3251 _CalcType __sum = 0.0;
3252 _M_cp.reserve(_M_int.size() - 1);
3253 _M_m.reserve(_M_int.size() - 1);
3254 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3255 {
3256 const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3257 __sum += _CalcType(0.5)
3258 * (_CalcType(_M_den[__k + 1]) + _CalcType(_M_den[__k]))
3259 * __delta;
3260 _M_cp.push_back(__sum);
3261 _M_m.push_back(
3262 (_CalcType(_M_den[__k + 1]) - _CalcType(_M_den[__k])) / __delta);
3263 }
3264 __glibcxx_assert(__sum > 0);
3265
3266 // Now normalize the densities...
3267 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3268 __sum);
3269 // ... and partial sums...
3270 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3271 // ... and slopes.
3272 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3273
3274 // Make sure the last cumulative probablility is one.
3275 _M_cp[_M_cp.size() - 1] = _CalcType(1);
3276 }
3277
3278 template<typename _RealType>
3279 void
3280 piecewise_linear_distribution<_RealType>::param_type::
3281 _M_initialize2(const _RealType* __ints, const _CalcType* __dens)
3282 {
3283 if (__ints[0] == _RealType(0)
3284 && __ints[1] == _RealType(1)
3285 && __dens[0] == __dens[1])
3286 return;
3287
3288 _M_int.reserve(2);
3289 _M_int.push_back(__ints[0]);
3290 _M_int.push_back(__ints[1]);
3291
3292 _M_den.reserve(2);
3293 _M_den.push_back(__dens[0]);
3294 _M_den.push_back(__dens[1]);
3295 _M_configure();
3296 }
3297
3298 template<typename _RealType>
3299 template<typename _InputIteratorB, typename _InputIteratorW>
3300 piecewise_linear_distribution<_RealType>::param_type::
3301 param_type(_InputIteratorB __bbegin,
3302 _InputIteratorB __bend,
3303 _InputIteratorW __wbegin)
3304 : _M_int(), _M_den(), _M_cp(), _M_m()
3305 {
3306 _RealType __ints[2];
3307 if (!__detail::__load_first2(__bbegin, __bend, __ints))
3308 return;
3309
3310 if (__bbegin == __bend)
3311 {
3312 _CalcType __dens[2];
3313 __dens[0] = *__wbegin;
3314 ++__wbegin;
3315 __dens[1] = *__wbegin;
3316 _M_initialize2(__ints, __dens);
3317 return;
3318 }
3319
3320#if __glibcxx_concepts // C++ >= C++20
3321 if constexpr (sized_sentinel_for<_InputIteratorB, _InputIteratorB>
3322 || forward_iterator<_InputIteratorB>)
3323 _M_int.reserve(2 + size_t(ranges::distance(__bbegin, __bend)));
3324#else
3325#pragma GCC diagnostic push
3326#pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3327 if constexpr (is_convertible<__iter_category_t<_InputIteratorB>,
3328 forward_iterator_tag>::value)
3329 _M_int.reserve(2 + size_t(std::distance(__bbegin, __bend)));
3330#pragma GCC diagnostic pop
3331#endif
3332
3333 _M_int.push_back(__ints[0]);
3334 _M_int.push_back(__ints[1]);
3335 for (; __bbegin != __bend; ++__bbegin)
3336 _M_int.push_back(*__bbegin);
3337
3338 // _GLIBCXX_RESOLVE_LIB_DEFECTS
3339 // 4052. Bogus requirements for piecewise_linear_distribution
3340 _M_den.reserve(_M_int.size());
3341 for (size_t __i = 0; __i < _M_int.size(); (void)++__i, ++__wbegin)
3342 _M_den.push_back(_CalcType(*__wbegin));
3343
3344 _M_configure();
3345 }
3346
3347 template<typename _RealType>
3348 template<typename _Func>
3349 piecewise_linear_distribution<_RealType>::param_type::
3350 param_type(initializer_list<_RealType> __bl, _Func __fw)
3351 : _M_int(), _M_den(), _M_cp(), _M_m()
3352 {
3353 if (__bl.size() < 2)
3354 return;
3355
3356 if (__bl.size() == 2)
3357 {
3358 const _RealType *__ints = __bl.begin();
3359 _CalcType __den[2];
3360 __den[0] = __fw(__ints[0]);
3361 __den[1] = __fw(__ints[1]);
3362 _M_initialize2(__ints, __den);
3363 return;
3364 }
3365
3366 _M_int = __bl;
3367 _M_den.reserve(__bl.size());
3368 for (_RealType __b : __bl)
3369 _M_den.push_back(_CalcType(__fw(__b)));
3370 _M_configure();
3371 }
3372
3373 template<typename _RealType>
3374 template<typename _Func>
3375 piecewise_linear_distribution<_RealType>::param_type::
3376 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3377 : _M_int(), _M_den(), _M_cp(), _M_m()
3378 {
3379 const size_t __n = __nw == 0 ? 1 : __nw;
3380 const _RealType __delta = (__xmax - __xmin) / __n;
3381 const auto __cfw = [&] (_RealType __v)
3382 {
3383#ifdef _GLIBCXX_USE_OLD_PIECEWISE_DISTRIBUTIONS
3384 return __fw(__v + __delta);
3385#else
3386 // _GLIBCXX_RESOLVE_LIB_DEFECTS
3387 // 4052. Bogus requirements for piecewise_linear_distribution
3388 return _RealType(__fw(__v));
3389#endif
3390 };
3391
3392 if (__n == 1)
3393 {
3394 _RealType __ints[2] = { __xmin, __xmin + __delta };
3395 _CalcType __dens[2];
3396 __dens[0] = __cfw(__ints[0]);
3397 __dens[1] = __cfw(__ints[1]);
3398 _M_initialize2(__ints, __dens);
3399 return;
3400 }
3401
3402 _M_int.reserve(__n + 1);
3403 _M_den.reserve(__n + 1);
3404 for (size_t __k = 0; __k <= __nw; ++__k)
3405 {
3406 _M_int.push_back(__xmin + __k * __delta);
3407 _M_den.push_back(__cfw(_M_int[__k]));
3408 }
3409
3410 _M_configure();
3411 }
3412
3413 template<typename _RealType>
3414 template<typename _AdaptedUniformRandomNumberGenerator>
3415 typename piecewise_linear_distribution<_RealType>::result_type
3416 piecewise_linear_distribution<_RealType>::
3417 __generate_one(_AdaptedUniformRandomNumberGenerator& __aurng,
3418 const param_type& __param)
3419 {
3420 const _CalcType __p = __aurng();
3421 auto __pos = std::lower_bound(__param._M_cp.begin(),
3422 __param._M_cp.end(), __p);
3423 const size_t __i = __pos - __param._M_cp.begin();
3424
3425 const _CalcType __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3426
3427 const _CalcType __a = _CalcType(0.5) * __param._M_m[__i];
3428 const _CalcType __b = __param._M_den[__i];
3429 const _CalcType __cm = __p - __pref;
3430
3431 _RealType __x = __param._M_int[__i];
3432 if (__a == 0)
3433 __x += __cm / __b;
3434 else
3435 {
3436 const _CalcType __d = __b * __b + 4.0 * __a * __cm;
3437 __x += _CalcType(0.5) * (std::sqrt(__d) - __b) / __a;
3438 }
3439 return __x;
3440 }
3441
3442
3443 template<typename _RealType>
3444 template<typename _UniformRandomNumberGenerator>
3445 typename piecewise_linear_distribution<_RealType>::result_type
3446 piecewise_linear_distribution<_RealType>::
3447 operator()(_UniformRandomNumberGenerator& __urng,
3448 const param_type& __param)
3449 {
3450 __detail::_Adaptor<_UniformRandomNumberGenerator, _CalcType>
3451 __aurng(__urng);
3452
3453 if (__param._M_cp.empty())
3454 return __aurng();
3455
3456 return __generate_one(__aurng, __param);
3457 }
3458
3459 template<typename _RealType>
3460 template<typename _ForwardIterator,
3461 typename _UniformRandomNumberGenerator>
3462 void
3463 piecewise_linear_distribution<_RealType>::
3464 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3465 _UniformRandomNumberGenerator& __urng,
3466 const param_type& __param)
3467 {
3468 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3469 __detail::_Adaptor<_UniformRandomNumberGenerator, _CalcType>
3470 __aurng(__urng);
3471
3472 if (__param._M_cp.empty())
3473 while (__f != __t)
3474 *__f++ = __aurng();
3475 else
3476 while (__f != __t)
3477 *__f++ = __generate_one(__urng, __param);
3478 }
3479
3480 template<typename _RealType, typename _CharT, typename _Traits>
3481 std::basic_ostream<_CharT, _Traits>&
3482 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3483 const piecewise_linear_distribution<_RealType>& __x)
3484 {
3485 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3486
3487 const typename __ios_base::fmtflags __flags = __os.flags();
3488 const _CharT __fill = __os.fill();
3489 const std::streamsize __precision = __os.precision();
3490 const _CharT __space = __os.widen(' ');
3491 __os.flags(__ios_base::scientific | __ios_base::left);
3492 __os.fill(__space);
3493 __os.precision(std::numeric_limits<_RealType>::max_digits10);
3494
3495 auto const& __int = __x._M_param._M_int;
3496 if (__int.empty())
3497 __os << size_t(1)
3498 << __space << _RealType(0) << __space << _RealType(1)
3499 << __space << _RealType(1) << __space << _RealType(1);
3500 else
3501 {
3502 __os << __int.size() - 1;
3503 for (auto __xv : __int)
3504 __os << __space << __xv;
3505 for (auto __dv : __x._M_param._M_den)
3506 __os << __space << __dv;
3507 }
3508
3509 __os.flags(__flags);
3510 __os.fill(__fill);
3511 __os.precision(__precision);
3512 return __os;
3513 }
3514
3515 template<typename _RealType, typename _CharT, typename _Traits>
3516 std::basic_istream<_CharT, _Traits>&
3517 operator>>(std::basic_istream<_CharT, _Traits>& __is,
3518 piecewise_linear_distribution<_RealType>& __x)
3519 {
3520 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3521 using _StorageType
3522 = typename piecewise_linear_distribution<_RealType>::_StorageType;
3523
3524 const typename __ios_base::fmtflags __flags = __is.flags();
3525 __is.flags(__ios_base::dec | __ios_base::skipws);
3526
3527 size_t __n;
3528 if (__is >> __n)
3529 {
3530 vector<_RealType> __int_vec;
3531 if (__detail::__extract_params(__is, __int_vec, __n + 1))
3532 {
3533 vector<_StorageType> __den_vec;
3534 if (__detail::__extract_params(__is, __den_vec, __n + 1))
3535 {
3536 __x.param({ __int_vec.begin(), __int_vec.end(),
3537 __den_vec.begin() });
3538 }
3539 }
3540 }
3541 __is.flags(__flags);
3542 return __is;
3543 }
3544
3545
3546 template<typename _IntType, typename>
3547 seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3548 {
3549 _M_v.reserve(__il.size());
3550 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3551 _M_v.push_back(__detail::__mod<result_type,
3552 __detail::_Shift<result_type, 32>::__value>(*__iter));
3553 }
3554
3555 template<typename _InputIterator>
3556 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3557 {
3558#pragma GCC diagnostic push
3559#pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3560 if constexpr (__is_random_access_iter<_InputIterator>::value)
3561 _M_v.reserve(std::distance(__begin, __end));
3562#pragma GCC diagnostic pop
3563
3564 for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3565 _M_v.push_back(__detail::__mod<result_type,
3566 __detail::_Shift<result_type, 32>::__value>(*__iter));
3567 }
3568
3569 template<typename _RandomAccessIterator>
3570 void
3571 seed_seq::generate(_RandomAccessIterator __begin,
3572 _RandomAccessIterator __end)
3573 {
3574 typedef typename iterator_traits<_RandomAccessIterator>::value_type
3575 _Type;
3576
3577 if (__begin == __end)
3578 return;
3579
3580 std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3581
3582 const size_t __n = __end - __begin;
3583 const size_t __s = _M_v.size();
3584 const size_t __t = (__n >= 623) ? 11
3585 : (__n >= 68) ? 7
3586 : (__n >= 39) ? 5
3587 : (__n >= 7) ? 3
3588 : (__n - 1) / 2;
3589 const size_t __p = (__n - __t) / 2;
3590 const size_t __q = __p + __t;
3591 const size_t __m = std::max(size_t(__s + 1), __n);
3592
3593#ifndef __UINT32_TYPE__
3594 struct _Up
3595 {
3596 _Up(uint_least32_t v) : _M_v(v & 0xffffffffu) { }
3597
3598 operator uint_least32_t() const { return _M_v; }
3599
3600 uint_least32_t _M_v;
3601 };
3602 using uint32_t = _Up;
3603#endif
3604
3605 // k == 0, every element in [begin,end) equals 0x8b8b8b8bu
3606 {
3607 uint32_t __r1 = 1371501266u;
3608 uint32_t __r2 = __r1 + __s;
3609 __begin[__p] += __r1;
3610 __begin[__q] = (uint32_t)__begin[__q] + __r2;
3611 __begin[0] = __r2;
3612 }
3613
3614 for (size_t __k = 1; __k <= __s; ++__k)
3615 {
3616 const size_t __kn = __k % __n;
3617 const size_t __kpn = (__k + __p) % __n;
3618 const size_t __kqn = (__k + __q) % __n;
3619 uint32_t __arg = (__begin[__kn]
3620 ^ __begin[__kpn]
3621 ^ __begin[(__k - 1) % __n]);
3622 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3623 uint32_t __r2 = __r1 + (uint32_t)__kn + _M_v[__k - 1];
3624 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3625 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3626 __begin[__kn] = __r2;
3627 }
3628
3629 for (size_t __k = __s + 1; __k < __m; ++__k)
3630 {
3631 const size_t __kn = __k % __n;
3632 const size_t __kpn = (__k + __p) % __n;
3633 const size_t __kqn = (__k + __q) % __n;
3634 uint32_t __arg = (__begin[__kn]
3635 ^ __begin[__kpn]
3636 ^ __begin[(__k - 1) % __n]);
3637 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3638 uint32_t __r2 = __r1 + (uint32_t)__kn;
3639 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3640 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3641 __begin[__kn] = __r2;
3642 }
3643
3644 for (size_t __k = __m; __k < __m + __n; ++__k)
3645 {
3646 const size_t __kn = __k % __n;
3647 const size_t __kpn = (__k + __p) % __n;
3648 const size_t __kqn = (__k + __q) % __n;
3649 uint32_t __arg = (__begin[__kn]
3650 + __begin[__kpn]
3651 + __begin[(__k - 1) % __n]);
3652 uint32_t __r3 = 1566083941u * (__arg ^ (__arg >> 27));
3653 uint32_t __r4 = __r3 - __kn;
3654 __begin[__kpn] ^= __r3;
3655 __begin[__kqn] ^= __r4;
3656 __begin[__kn] = __r4;
3657 }
3658 }
3659
3660// [rand.util.canonical]
3661// generate_canonical(RNG&)
3662
3663#ifndef _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3664
3665#pragma GCC diagnostic push
3666#pragma GCC diagnostic ignored "-Wc++14-extensions" // for variable templates
3667#pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3668
3669 // __generate_canonical_pow2 is used when Urbg::max()-Urbg::min() is
3670 // a power of two less 1. It works by calling urng() as many times as
3671 // needed to fill the target mantissa, accumulating entropy into an
3672 // integer value, converting that to the float type, and then dividing
3673 // by the range of the integer value (a constexpr power of 2,
3674 // so only adjusts the exponent) to produce a result in [0..1].
3675 // In case of an exact 1.0 result, we re-try.
3676 //
3677 // It needs to work even when the integer type used is only as big
3678 // as the float mantissa, such as uint64_t for long double. So,
3679 // commented-out assignments represent computations the Standard
3680 // prescribes but cannot be performed, or are not used. Names are
3681 // chosen to match the description in the Standard.
3682 //
3683 // When the result is close to zero, the strict Standard-prescribed
3684 // calculation may leave more low-order zeros in the mantissa than
3685 // is usually necessary. When spare entropy has been extracted, as
3686 // is usual for float and double, some or all of the spare entropy
3687 // can commonly be pulled into the result for better randomness.
3688 // Defining _GLIBCXX_GENERATE_CANONICAL_STRICT discards it instead.
3689 //
3690 // When k calls to urng() yield more bits of entropy, log2_Rk_max,
3691 // than fit into UInt, we discard some of it by overflowing, which
3692 // is OK. On converting the integer representation of the sample
3693 // to the float value, we must divide out the (possibly-truncated)
3694 // size log2_Rk.
3695 //
3696 // This implementation works with std::bfloat16, which can exactly
3697 // represent 2^32, but not with std::float16_t, limited to 2^15.
3698
3699 template<typename _RealT, size_t __d, typename _Urbg>
3700 _RealT
3701 __generate_canonical_pow2(_Urbg& __urng)
3702 {
3703 using _UInt = typename __detail::_Select_uint_least_t<__d>::type;
3704
3705 // Parameter __d is the actual target number of bits.
3706 // Commented-out assignments below are of values specified in
3707 // the Standard, but not used here for reasons noted.
3708 // r = 2; // Redundant, we only support radix 2.
3709 using _Rng = decltype(_Urbg::max());
3710 const _Rng __rng_range_less_1 = _Urbg::max() - _Urbg::min();
3711 // R = _UInt(__rng_range_less_1) + 1; // May wrap to 0.
3712 const auto __log2_R = __builtin_popcountg(__rng_range_less_1);
3713 const auto __log2_uint_max = sizeof(_UInt) * __CHAR_BIT__;
3714 // rd = _UInt(1) << __d; // Could overflow, UB.
3715 const unsigned __k = (__d + __log2_R - 1) / __log2_R;
3716 const unsigned __log2_Rk_max = __k * __log2_R;
3717 const unsigned __log2_Rk = // Bits of entropy actually obtained:
3718 __log2_uint_max < __log2_Rk_max ? __log2_uint_max : __log2_Rk_max;
3719 // Rk = _UInt(1) << __log2_Rk; // Likely overflows, UB.
3720 _GLIBCXX14_CONSTEXPR const _RealT __Rk
3721 = _RealT(_UInt(1) << (__log2_Rk - 1)) * _RealT(2.0);
3722#if defined(_GLIBCXX_GENERATE_CANONICAL_STRICT)
3723 const unsigned __log2_x = __log2_Rk - __d; // # of spare entropy bits.
3724#else
3725 const unsigned __log2_x = 0;
3726#endif
3727 _GLIBCXX14_CONSTEXPR const _UInt __x = _UInt(1) << __log2_x;
3728 _GLIBCXX14_CONSTEXPR const _RealT __rd = __Rk / _RealT(__x);
3729 // xrd = __x << __d; // Could overflow.
3730
3731 while (true)
3732 {
3733 _UInt __sum = _UInt(__urng() - _Urbg::min());
3734 for (unsigned __i = __k - 1, __shift = 0; __i > 0; --__i)
3735 {
3736 __shift += __log2_R;
3737 __sum |= _UInt(__urng() - _Urbg::min()) << __shift;
3738 }
3739 const _RealT __ret = _RealT(__sum >> __log2_x) / _RealT(__rd);
3740 if (__ret < _RealT(1.0))
3741 return __ret;
3742 }
3743 }
3744
3745
3746 template<typename _UInt>
3747 struct __gen_canon_log_res
3748 {
3749 unsigned __floor_log;
3750 _UInt __floor_pow;
3751
3752 constexpr __gen_canon_log_res
3753 update(_UInt __base) const
3754 { return {__floor_log + 1, __floor_pow * __base}; }
3755 };
3756
3757
3758 template <typename _UInt1, typename _UInt2,
3759 typename _UComm = __conditional_t<(sizeof(_UInt2) > sizeof(_UInt1)),
3760 _UInt2, _UInt1>>
3761 constexpr __gen_canon_log_res<_UInt1>
3762 __gen_canon_log(_UInt1 __val, _UInt2 __base)
3763 {
3764#if __cplusplus >= 201402L
3765 __gen_canon_log_res<_UInt1> __res{0, _UInt1(1)};
3766 if (_UComm(__base) > _UComm(__val))
3767 return __res;
3768
3769 const _UInt1 __base1(__base);
3770 do
3771 {
3772 __val /= __base1;
3773 __res = __res.update(__base1);
3774 }
3775 while (__val >= __base1);
3776 return __res;
3777#else
3778 return (_UComm(__val) >= _UComm(__base))
3779 ? __gen_canon_log(__val / _UInt1(__base), _UInt1(__base))
3780 .update(_UInt1(__base))
3781 : __gen_canon_log_res<_UInt1>{0, _UInt1(1)};
3782#endif
3783 }
3784
3785 // This version must be used when the range of possible RNG results,
3786 // Urbg::max()-Urbg::min(), is not a power of two less one. The UInt
3787 // type passed must be big enough to represent Rk, R^k, a power of R
3788 // (the range of values produced by the rng) up to twice the length
3789 // of the mantissa.
3790
3791 template<typename _RealT, size_t __d, typename _Urbg>
3792 _RealT
3793 __generate_canonical_any(_Urbg& __urng)
3794 {
3795 // Names below are chosen to match the description in the Standard.
3796 // Parameter d is the actual target number of bits.
3797#if (__cplusplus >= 201402L) || defined(__SIZEOF_INT128__)
3798# define _GLIBCXX_GEN_CANON_CONST constexpr
3799#else
3800# define _GLIBCXX_GEN_CANON_CONST const
3801#endif
3802
3803 using _UIntR = typename make_unsigned<decltype(_Urbg::max())>::type;
3804 // Cannot overflow, as _Urbg::max() - _Urbg::min() is not power of
3805 // two minus one
3806 constexpr _UIntR __R = _UIntR(_Urbg::max() - _Urbg::min()) + 1;
3807 constexpr unsigned __log2R
3808 = sizeof(_UIntR) * __CHAR_BIT__ - __builtin_clzg(__R) - 1;
3809 // We overstimate number of required bits, by computing
3810 // r such that l * log2(R) >= d, so:
3811 // R^l >= (2 ^ log2(R)) ^ l == 2 ^ (log2(r) * l) >= 2^d
3812 // And then requiring l * bit_width(R) bits.
3813 constexpr unsigned __l = (__d + __log2R - 1) / __log2R;
3814 constexpr unsigned __bits = (__log2R + 1) * __l;
3815 using _UInt = typename __detail::_Select_uint_least_t<__bits>::type;
3816
3817 _GLIBCXX_GEN_CANON_CONST _UInt __rd = _UInt(1) << __d;
3818 _GLIBCXX_GEN_CANON_CONST auto __logRrd = __gen_canon_log(__rd, __R);
3819 _GLIBCXX_GEN_CANON_CONST unsigned __k
3820 = __logRrd.__floor_log + (__rd > __logRrd.__floor_pow);
3821
3822 _GLIBCXX_GEN_CANON_CONST _UInt __Rk
3823 = (__k > __logRrd.__floor_log)
3824 ? _UInt(__logRrd.__floor_pow) * _UInt(__R)
3825 : _UInt(__logRrd.__floor_pow);
3826 _GLIBCXX_GEN_CANON_CONST _UInt __x = __Rk / __rd;
3827
3828 while (true)
3829 {
3830 _UInt __Ri{1};
3831 _UInt __sum(__urng() - _Urbg::min());
3832 for (int __i = __k - 1; __i > 0; --__i)
3833 {
3834 __Ri *= _UInt(__R);
3835 __sum += _UInt(__urng() - _Urbg::min()) * __Ri;
3836 }
3837 const _RealT __ret = _RealT(__sum / __x) / _RealT(__rd);
3838 if (__ret < _RealT(1.0))
3839 return __ret;
3840 }
3841#undef _GLIBCXX_GEN_CANON_CONST
3842 }
3843
3844#if !defined(_GLIBCXX_GENERATE_CANONICAL_STRICT)
3845 template <typename _Tp>
3846 const bool __is_rand_dist_float_v = is_floating_point<_Tp>::value;
3847#else
3848 template <typename _Tp> const bool __is_rand_dist_float_v = false;
3849 template <> const bool __is_rand_dist_float_v<float> = true;
3850 template <> const bool __is_rand_dist_float_v<double> = true;
3851 template <> const bool __is_rand_dist_float_v<long double> = true;
3852#endif
3853
3854 // Note, this works even when (__range + 1) overflows:
3855 template <typename _Rng>
3856 constexpr bool __is_power_of_2_less_1(_Rng __range)
3857 { return ((__range + 1) & __range) == 0; };
3858
3859_GLIBCXX_BEGIN_INLINE_ABI_NAMESPACE(_V2)
3860 /** Produce a random floating-point value in the range [0..1)
3861 *
3862 * The result of `std::generate_canonical<RealT,digits>(urng)` is a
3863 * random floating-point value of type `RealT` in the range [0..1),
3864 * using entropy provided by the uniform random bit generator `urng`.
3865 * A value for `digits` may be passed to limit the precision of the
3866 * result to so many bits, but normally `-1u` is passed to get the
3867 * native precision of `RealT`. As many `urng()` calls are made as
3868 * needed to obtain the required entropy. On rare occasions, more
3869 * `urng()` calls are used. It is fastest when the value of
3870 * `Urbg::max()` is a power of two less one, such as from a
3871 * `std::philox4x32` (for `float`) or `philox4x64` (for `double`).
3872 *
3873 * @since C++11
3874 */
3875 template<typename _RealT, size_t __digits,
3876 typename _Urbg>
3877 _RealT
3878 generate_canonical(_Urbg& __urng)
3879 {
3880#ifdef __glibcxx_concepts
3881 static_assert(uniform_random_bit_generator<_Urbg>);
3882#endif
3883 static_assert(__is_rand_dist_float_v<_RealT>,
3884 "template argument must be a floating point type");
3885 static_assert(__digits != 0 && _Urbg::max() > _Urbg::min(),
3886 "random samples with 0 bits are not meaningful");
3887 static_assert(std::numeric_limits<_RealT>::radix == 2,
3888 "only base-2 float types are supported");
3889#if defined(__STDCPP_FLOAT16_T__)
3890 static_assert(! is_same_v<_RealT, _Float16>,
3891 "float16_t type is not supported, consider using bfloat16_t");
3892#endif
3893
3894 const unsigned __d_max = std::numeric_limits<_RealT>::digits;
3895 const unsigned __d = __digits > __d_max ? __d_max : __digits;
3896
3897 // If the RNG range is a power of 2 less 1, the float type mantissa
3898 // is enough bits. If not, we need more.
3899 if constexpr (__is_power_of_2_less_1(_Urbg::max() - _Urbg::min()))
3900 return __generate_canonical_pow2<_RealT, __d>(__urng);
3901 else // Need up to 2x bits.
3902 return __generate_canonical_any<_RealT, __d>(__urng);
3903 }
3904_GLIBCXX_END_INLINE_ABI_NAMESPACE(_V2)
3905
3906#pragma GCC diagnostic pop
3907
3908#else // _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3909
3910 // This is the pre-P0952 definition, to reproduce old results.
3911
3912 template<typename _RealType, size_t __bits,
3913 typename _UniformRandomNumberGenerator>
3914 _RealType
3915 generate_canonical(_UniformRandomNumberGenerator& __urng)
3916 {
3917 static_assert(std::is_floating_point<_RealType>::value,
3918 "template argument must be a floating point type");
3919
3920 const size_t __b
3921 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3922 __bits);
3923 const long double __r = static_cast<long double>(__urng.max())
3924 - static_cast<long double>(__urng.min()) + 1.0L;
3925 const size_t __log2r = std::log(__r) / std::log(2.0L);
3926 const size_t __m = std::max<size_t>(1UL,
3927 (__b + __log2r - 1UL) / __log2r);
3928 _RealType __ret;
3929 _RealType __sum = _RealType(0);
3930 _RealType __tmp = _RealType(1);
3931 for (size_t __k = __m; __k != 0; --__k)
3932 {
3933 __sum += _RealType(__urng() - __urng.min()) * __tmp;
3934 __tmp *= __r;
3935 }
3936 __ret = __sum / __tmp;
3937 if (__builtin_expect(__ret >= _RealType(1), 0))
3938 {
3939# if _GLIBCXX_USE_C99_MATH_FUNCS
3940 __ret = std::nextafter(_RealType(1), _RealType(0));
3941# else
3942 __ret = _RealType(1)
3943 - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3944# endif
3945 }
3946 return __ret;
3947 }
3948
3949#endif // _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3950
3951_GLIBCXX_END_NAMESPACE_VERSION
3952} // namespace
3953
3954#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:177
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:1702
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:1798
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:53
std::array
A standard container for storing a fixed size sequence of elements.
Definition array:104
std::basic_istream
Template class basic_istream.
Definition istream:73
std::__numeric_limits_base::is_integer
static constexpr bool is_integer
Definition limits:239
std::__numeric_limits_base::max_digits10
static constexpr int max_digits10
Definition limits:232
std::__numeric_limits_base::digits
static constexpr int digits
Definition limits:224
std::__numeric_limits_base::is_signed
static constexpr bool is_signed
Definition limits:236
std::__numeric_limits_base::radix
static constexpr int radix
Definition limits:248
std::numeric_limits::max
static constexpr _Tp max() noexcept
Definition limits:336
std::numeric_limits::epsilon
static constexpr _Tp epsilon() noexcept
Definition limits:348
std::is_floating_point
is_floating_point
Definition type_traits:624
std::common_type
common_type
Definition type_traits:2600
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::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::basic_ostream
Template class basic_ostream.
Definition ostream.h:72
std::linear_congruential_engine
A model of a linear congruential random number generator.
Definition random.h:728
std::linear_congruential_engine::multiplier
static constexpr result_type multiplier
Definition random.h:743
std::linear_congruential_engine::modulus
static constexpr result_type modulus
Definition random.h:747
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:745
std::subtract_with_carry_engine
The Marsaglia-Zaman generator.
Definition random.h:1175
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:1825
std::shuffle_order_engine::base
const _RandomNumberEngine & base() const noexcept
Definition random.h:1931
std::shuffle_order_engine::result_type
_RandomNumberEngine::result_type result_type
Definition random.h:1831
std::philox_engine
A discrete pseudorandom number generator with weak cryptographic properties.
Definition random.h:2084
std::uniform_real_distribution
Uniform continuous distribution for random numbers.
Definition random.h:2517
std::uniform_real_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:2606
std::normal_distribution
A normal continuous distribution for random numbers.
Definition random.h:2754
std::normal_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:2873
std::gamma_distribution
A gamma continuous distribution for random numbers.
Definition random.h:3206
std::gamma_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:3335
std::cauchy_distribution
A cauchy_distribution random number distribution.
Definition random.h:3678
std::cauchy_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:3755
std::cauchy_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:3785
std::fisher_f_distribution
A fisher_f_distribution random number distribution.
Definition random.h:3893
std::binomial_distribution
A discrete binomial random number distribution.
Definition random.h:4589
std::binomial_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:4717
std::geometric_distribution
A discrete geometric random number distribution.
Definition random.h:4835
std::geometric_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:4946
std::geometric_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:4916
std::negative_binomial_distribution
A negative_binomial_distribution random number distribution.
Definition random.h:5052
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:5289
std::poisson_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:5402
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:5438
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:5521
std::exponential_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:5601
std::weibull_distribution
A weibull_distribution random number distribution.
Definition random.h:5743
std::weibull_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:5823
std::weibull_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:5853
std::extreme_value_distribution
A extreme_value_distribution random number distribution.
Definition random.h:5960
std::extreme_value_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition random.h:6070
std::extreme_value_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition random.h:6040
std::forward_iterator_tag
Forward iterators support a superset of input iterator operations.
Definition stl_iterator_base_types.h:101
std::vector
A standard container which offers fixed time access to individual elements in any order.
Definition stl_vector.h:511
std::vector::end
constexpr iterator end() noexcept
Definition stl_vector.h:1060
std::vector::begin
constexpr iterator begin() noexcept
Definition stl_vector.h:1040
std::vector::size
constexpr size_type size() const noexcept
Definition stl_vector.h:1159
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