libstdc++
bits/random.tcc
Go to the documentation of this file.
1// random number generation (out of line) -*- C++ -*-
2
3// Copyright (C) 2009-2024 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library. This library is free
6// software; you can redistribute it and/or modify it under the
7// terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 3, or (at your option)
9// any later version.
10
11// This library is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14// GNU General Public License for more details.
15
16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.
19
20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23// <http://www.gnu.org/licenses/>.
24
25/** @file bits/random.tcc
26 * This is an internal header file, included by other library headers.
27 * Do not attempt to use it directly. @headername{random}
28 */
29
30#ifndef _RANDOM_TCC
31#define _RANDOM_TCC 1
32
33#include <numeric> // std::accumulate and std::partial_sum
34
35namespace std _GLIBCXX_VISIBILITY(default)
36{
37_GLIBCXX_BEGIN_NAMESPACE_VERSION
38
39 /// @cond undocumented
40 // (Further) implementation-space details.
41 namespace __detail
42 {
43 // General case for x = (ax + c) mod m -- use Schrage's algorithm
44 // to avoid integer overflow.
45 //
46 // Preconditions: a > 0, m > 0.
47 //
48 // Note: only works correctly for __m % __a < __m / __a.
49 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
50 _Tp
51 _Mod<_Tp, __m, __a, __c, false, true>::
52 __calc(_Tp __x)
53 {
54 if (__a == 1)
55 __x %= __m;
56 else
57 {
58 static const _Tp __q = __m / __a;
59 static const _Tp __r = __m % __a;
60
61 _Tp __t1 = __a * (__x % __q);
62 _Tp __t2 = __r * (__x / __q);
63 if (__t1 >= __t2)
64 __x = __t1 - __t2;
65 else
66 __x = __m - __t2 + __t1;
67 }
68
69 if (__c != 0)
70 {
71 const _Tp __d = __m - __x;
72 if (__d > __c)
73 __x += __c;
74 else
75 __x = __c - __d;
76 }
77 return __x;
78 }
79
80 template<typename _InputIterator, typename _OutputIterator,
81 typename _Tp>
82 _OutputIterator
83 __normalize(_InputIterator __first, _InputIterator __last,
84 _OutputIterator __result, const _Tp& __factor)
85 {
86 for (; __first != __last; ++__first, ++__result)
87 *__result = *__first / __factor;
88 return __result;
89 }
90
91 } // namespace __detail
92 /// @endcond
93
94#if ! __cpp_inline_variables
95 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
96 constexpr _UIntType
97 linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
98
99 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
100 constexpr _UIntType
101 linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
102
103 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
104 constexpr _UIntType
105 linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
106
107 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
108 constexpr _UIntType
109 linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
110#endif
111
112 /**
113 * Seeds the LCR with integral value @p __s, adjusted so that the
114 * ring identity is never a member of the convergence set.
115 */
116 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
117 void
118 linear_congruential_engine<_UIntType, __a, __c, __m>::
119 seed(result_type __s)
120 {
121 if ((__detail::__mod<_UIntType, __m>(__c) == 0)
122 && (__detail::__mod<_UIntType, __m>(__s) == 0))
123 _M_x = 1;
124 else
125 _M_x = __detail::__mod<_UIntType, __m>(__s);
126 }
127
128 /**
129 * Seeds the LCR engine with a value generated by @p __q.
130 */
131 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
132 template<typename _Sseq>
133 auto
134 linear_congruential_engine<_UIntType, __a, __c, __m>::
135 seed(_Sseq& __q)
136 -> _If_seed_seq<_Sseq>
137 {
138 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
139 : std::__lg(__m);
140 const _UIntType __k = (__k0 + 31) / 32;
141 uint_least32_t __arr[__k + 3];
142 __q.generate(__arr + 0, __arr + __k + 3);
143 _UIntType __factor = 1u;
144 _UIntType __sum = 0u;
145 for (size_t __j = 0; __j < __k; ++__j)
146 {
147 __sum += __arr[__j + 3] * __factor;
148 __factor *= __detail::_Shift<_UIntType, 32>::__value;
149 }
150 seed(__sum);
151 }
152
153 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
154 typename _CharT, typename _Traits>
155 std::basic_ostream<_CharT, _Traits>&
156 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
157 const linear_congruential_engine<_UIntType,
158 __a, __c, __m>& __lcr)
159 {
160 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
161
162 const typename __ios_base::fmtflags __flags = __os.flags();
163 const _CharT __fill = __os.fill();
164 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
165 __os.fill(__os.widen(' '));
166
167 __os << __lcr._M_x;
168
169 __os.flags(__flags);
170 __os.fill(__fill);
171 return __os;
172 }
173
174 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
175 typename _CharT, typename _Traits>
176 std::basic_istream<_CharT, _Traits>&
177 operator>>(std::basic_istream<_CharT, _Traits>& __is,
178 linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
179 {
180 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
181
182 const typename __ios_base::fmtflags __flags = __is.flags();
183 __is.flags(__ios_base::dec);
184
185 __is >> __lcr._M_x;
186
187 __is.flags(__flags);
188 return __is;
189 }
190
191#if ! __cpp_inline_variables
192 template<typename _UIntType,
193 size_t __w, size_t __n, size_t __m, size_t __r,
194 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
195 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
196 _UIntType __f>
197 constexpr size_t
198 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
199 __s, __b, __t, __c, __l, __f>::word_size;
200
201 template<typename _UIntType,
202 size_t __w, size_t __n, size_t __m, size_t __r,
203 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
204 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
205 _UIntType __f>
206 constexpr size_t
207 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
208 __s, __b, __t, __c, __l, __f>::state_size;
209
210 template<typename _UIntType,
211 size_t __w, size_t __n, size_t __m, size_t __r,
212 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
213 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
214 _UIntType __f>
215 constexpr size_t
216 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
217 __s, __b, __t, __c, __l, __f>::shift_size;
218
219 template<typename _UIntType,
220 size_t __w, size_t __n, size_t __m, size_t __r,
221 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
222 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
223 _UIntType __f>
224 constexpr size_t
225 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
226 __s, __b, __t, __c, __l, __f>::mask_bits;
227
228 template<typename _UIntType,
229 size_t __w, size_t __n, size_t __m, size_t __r,
230 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
231 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
232 _UIntType __f>
233 constexpr _UIntType
234 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
235 __s, __b, __t, __c, __l, __f>::xor_mask;
236
237 template<typename _UIntType,
238 size_t __w, size_t __n, size_t __m, size_t __r,
239 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
240 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
241 _UIntType __f>
242 constexpr size_t
243 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
244 __s, __b, __t, __c, __l, __f>::tempering_u;
245
246 template<typename _UIntType,
247 size_t __w, size_t __n, size_t __m, size_t __r,
248 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
249 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
250 _UIntType __f>
251 constexpr _UIntType
252 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
253 __s, __b, __t, __c, __l, __f>::tempering_d;
254
255 template<typename _UIntType,
256 size_t __w, size_t __n, size_t __m, size_t __r,
257 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
258 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
259 _UIntType __f>
260 constexpr size_t
261 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
262 __s, __b, __t, __c, __l, __f>::tempering_s;
263
264 template<typename _UIntType,
265 size_t __w, size_t __n, size_t __m, size_t __r,
266 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
267 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
268 _UIntType __f>
269 constexpr _UIntType
270 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
271 __s, __b, __t, __c, __l, __f>::tempering_b;
272
273 template<typename _UIntType,
274 size_t __w, size_t __n, size_t __m, size_t __r,
275 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
276 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
277 _UIntType __f>
278 constexpr size_t
279 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
280 __s, __b, __t, __c, __l, __f>::tempering_t;
281
282 template<typename _UIntType,
283 size_t __w, size_t __n, size_t __m, size_t __r,
284 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
285 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
286 _UIntType __f>
287 constexpr _UIntType
288 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
289 __s, __b, __t, __c, __l, __f>::tempering_c;
290
291 template<typename _UIntType,
292 size_t __w, size_t __n, size_t __m, size_t __r,
293 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
294 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
295 _UIntType __f>
296 constexpr size_t
297 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
298 __s, __b, __t, __c, __l, __f>::tempering_l;
299
300 template<typename _UIntType,
301 size_t __w, size_t __n, size_t __m, size_t __r,
302 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
303 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
304 _UIntType __f>
305 constexpr _UIntType
306 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
307 __s, __b, __t, __c, __l, __f>::
308 initialization_multiplier;
309
310 template<typename _UIntType,
311 size_t __w, size_t __n, size_t __m, size_t __r,
312 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
313 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
314 _UIntType __f>
315 constexpr _UIntType
316 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
317 __s, __b, __t, __c, __l, __f>::default_seed;
318#endif
319
320 template<typename _UIntType,
321 size_t __w, size_t __n, size_t __m, size_t __r,
322 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
323 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
324 _UIntType __f>
325 void
326 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
327 __s, __b, __t, __c, __l, __f>::
328 seed(result_type __sd)
329 {
330 _M_x[0] = __detail::__mod<_UIntType,
331 __detail::_Shift<_UIntType, __w>::__value>(__sd);
332
333 for (size_t __i = 1; __i < state_size; ++__i)
334 {
335 _UIntType __x = _M_x[__i - 1];
336 __x ^= __x >> (__w - 2);
337 __x *= __f;
338 __x += __detail::__mod<_UIntType, __n>(__i);
339 _M_x[__i] = __detail::__mod<_UIntType,
340 __detail::_Shift<_UIntType, __w>::__value>(__x);
341 }
342 _M_p = state_size;
343 }
344
345 template<typename _UIntType,
346 size_t __w, size_t __n, size_t __m, size_t __r,
347 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
348 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
349 _UIntType __f>
350 template<typename _Sseq>
351 auto
352 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
353 __s, __b, __t, __c, __l, __f>::
354 seed(_Sseq& __q)
355 -> _If_seed_seq<_Sseq>
356 {
357 const _UIntType __upper_mask = (~_UIntType()) << __r;
358 const size_t __k = (__w + 31) / 32;
359 uint_least32_t __arr[__n * __k];
360 __q.generate(__arr + 0, __arr + __n * __k);
361
362 bool __zero = true;
363 for (size_t __i = 0; __i < state_size; ++__i)
364 {
365 _UIntType __factor = 1u;
366 _UIntType __sum = 0u;
367 for (size_t __j = 0; __j < __k; ++__j)
368 {
369 __sum += __arr[__k * __i + __j] * __factor;
370 __factor *= __detail::_Shift<_UIntType, 32>::__value;
371 }
372 _M_x[__i] = __detail::__mod<_UIntType,
373 __detail::_Shift<_UIntType, __w>::__value>(__sum);
374
375 if (__zero)
376 {
377 if (__i == 0)
378 {
379 if ((_M_x[0] & __upper_mask) != 0u)
380 __zero = false;
381 }
382 else if (_M_x[__i] != 0u)
383 __zero = false;
384 }
385 }
386 if (__zero)
387 _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
388 _M_p = state_size;
389 }
390
391 template<typename _UIntType, size_t __w,
392 size_t __n, size_t __m, size_t __r,
393 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
394 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
395 _UIntType __f>
396 void
397 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
398 __s, __b, __t, __c, __l, __f>::
399 _M_gen_rand(void)
400 {
401 const _UIntType __upper_mask = (~_UIntType()) << __r;
402 const _UIntType __lower_mask = ~__upper_mask;
403
404 for (size_t __k = 0; __k < (__n - __m); ++__k)
405 {
406 _UIntType __y = ((_M_x[__k] & __upper_mask)
407 | (_M_x[__k + 1] & __lower_mask));
408 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
409 ^ ((__y & 0x01) ? __a : 0));
410 }
411
412 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
413 {
414 _UIntType __y = ((_M_x[__k] & __upper_mask)
415 | (_M_x[__k + 1] & __lower_mask));
416 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
417 ^ ((__y & 0x01) ? __a : 0));
418 }
419
420 _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
421 | (_M_x[0] & __lower_mask));
422 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
423 ^ ((__y & 0x01) ? __a : 0));
424 _M_p = 0;
425 }
426
427 template<typename _UIntType, size_t __w,
428 size_t __n, size_t __m, size_t __r,
429 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
430 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
431 _UIntType __f>
432 void
433 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
434 __s, __b, __t, __c, __l, __f>::
435 discard(unsigned long long __z)
436 {
437 while (__z > state_size - _M_p)
438 {
439 __z -= state_size - _M_p;
440 _M_gen_rand();
441 }
442 _M_p += __z;
443 }
444
445 template<typename _UIntType, size_t __w,
446 size_t __n, size_t __m, size_t __r,
447 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
448 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
449 _UIntType __f>
450 typename
451 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
452 __s, __b, __t, __c, __l, __f>::result_type
453 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
454 __s, __b, __t, __c, __l, __f>::
455 operator()()
456 {
457 // Reload the vector - cost is O(n) amortized over n calls.
458 if (_M_p >= state_size)
459 _M_gen_rand();
460
461 // Calculate o(x(i)).
462 result_type __z = _M_x[_M_p++];
463 __z ^= (__z >> __u) & __d;
464 __z ^= (__z << __s) & __b;
465 __z ^= (__z << __t) & __c;
466 __z ^= (__z >> __l);
467
468 return __z;
469 }
470
471 template<typename _UIntType, size_t __w,
472 size_t __n, size_t __m, size_t __r,
473 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
474 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
475 _UIntType __f, typename _CharT, typename _Traits>
476 std::basic_ostream<_CharT, _Traits>&
477 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
478 const mersenne_twister_engine<_UIntType, __w, __n, __m,
479 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
480 {
481 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
482
483 const typename __ios_base::fmtflags __flags = __os.flags();
484 const _CharT __fill = __os.fill();
485 const _CharT __space = __os.widen(' ');
486 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
487 __os.fill(__space);
488
489 for (size_t __i = 0; __i < __n; ++__i)
490 __os << __x._M_x[__i] << __space;
491 __os << __x._M_p;
492
493 __os.flags(__flags);
494 __os.fill(__fill);
495 return __os;
496 }
497
498 template<typename _UIntType, size_t __w,
499 size_t __n, size_t __m, size_t __r,
500 _UIntType __a, size_t __u, _UIntType __d, size_t __s,
501 _UIntType __b, size_t __t, _UIntType __c, size_t __l,
502 _UIntType __f, typename _CharT, typename _Traits>
503 std::basic_istream<_CharT, _Traits>&
504 operator>>(std::basic_istream<_CharT, _Traits>& __is,
505 mersenne_twister_engine<_UIntType, __w, __n, __m,
506 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
507 {
508 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
509
510 const typename __ios_base::fmtflags __flags = __is.flags();
511 __is.flags(__ios_base::dec | __ios_base::skipws);
512
513 for (size_t __i = 0; __i < __n; ++__i)
514 __is >> __x._M_x[__i];
515 __is >> __x._M_p;
516
517 __is.flags(__flags);
518 return __is;
519 }
520
521#if ! __cpp_inline_variables
522 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
523 constexpr size_t
524 subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
525
526 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
527 constexpr size_t
528 subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
529
530 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
531 constexpr size_t
532 subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
533
534 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
535 constexpr uint_least32_t
536 subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
537#endif
538
539 template<typename _UIntType, size_t __w, size_t __s, size_t __r>
540 void
541 subtract_with_carry_engine<_UIntType, __w, __s, __r>::
542 seed(result_type __value)
543 {
544 // _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 // Avoid using slower long double arithmetic if possible.
853 if _GLIBCXX17_CONSTEXPR (__detail::__p1_representable_as_double(__range))
854 __j *= __y / (__range + 1.0);
855 else
856 __j *= __y / (__range + 1.0L);
857 _M_y = _M_v[__j];
858 _M_v[__j] = _M_b();
859
860 return _M_y;
861 }
862
863 template<typename _RandomNumberEngine, size_t __k,
864 typename _CharT, typename _Traits>
865 std::basic_ostream<_CharT, _Traits>&
866 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
867 const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
868 {
869 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
870
871 const typename __ios_base::fmtflags __flags = __os.flags();
872 const _CharT __fill = __os.fill();
873 const _CharT __space = __os.widen(' ');
874 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
875 __os.fill(__space);
876
877 __os << __x.base();
878 for (size_t __i = 0; __i < __k; ++__i)
879 __os << __space << __x._M_v[__i];
880 __os << __space << __x._M_y;
881
882 __os.flags(__flags);
883 __os.fill(__fill);
884 return __os;
885 }
886
887 template<typename _RandomNumberEngine, size_t __k,
888 typename _CharT, typename _Traits>
889 std::basic_istream<_CharT, _Traits>&
890 operator>>(std::basic_istream<_CharT, _Traits>& __is,
891 shuffle_order_engine<_RandomNumberEngine, __k>& __x)
892 {
893 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
894
895 const typename __ios_base::fmtflags __flags = __is.flags();
896 __is.flags(__ios_base::dec | __ios_base::skipws);
897