libstdc++
bits/random.tcc
1 // random number generation (out of line) -*- C++ -*-
2 
3 // Copyright (C) 2009-2014 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 
35 namespace std _GLIBCXX_VISIBILITY(default)
36 {
37  /*
38  * (Further) implementation-space details.
39  */
40  namespace __detail
41  {
42  _GLIBCXX_BEGIN_NAMESPACE_VERSION
43 
44  // General case for x = (ax + c) mod m -- use Schrage's algorithm
45  // to avoid integer overflow.
46  //
47  // Preconditions: a > 0, m > 0.
48  //
49  // Note: only works correctly for __m % __a < __m / __a.
50  template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
51  _Tp
52  _Mod<_Tp, __m, __a, __c, false, true>::
53  __calc(_Tp __x)
54  {
55  if (__a == 1)
56  __x %= __m;
57  else
58  {
59  static const _Tp __q = __m / __a;
60  static const _Tp __r = __m % __a;
61 
62  _Tp __t1 = __a * (__x % __q);
63  _Tp __t2 = __r * (__x / __q);
64  if (__t1 >= __t2)
65  __x = __t1 - __t2;
66  else
67  __x = __m - __t2 + __t1;
68  }
69 
70  if (__c != 0)
71  {
72  const _Tp __d = __m - __x;
73  if (__d > __c)
74  __x += __c;
75  else
76  __x = __c - __d;
77  }
78  return __x;
79  }
80 
81  template<typename _InputIterator, typename _OutputIterator,
82  typename _Tp>
83  _OutputIterator
84  __normalize(_InputIterator __first, _InputIterator __last,
85  _OutputIterator __result, const _Tp& __factor)
86  {
87  for (; __first != __last; ++__first, ++__result)
88  *__result = *__first / __factor;
89  return __result;
90  }
91 
92  _GLIBCXX_END_NAMESPACE_VERSION
93  } // namespace __detail
94 
95 _GLIBCXX_BEGIN_NAMESPACE_VERSION
96 
97  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
98  constexpr _UIntType
99  linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
100 
101  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
102  constexpr _UIntType
103  linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
104 
105  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
106  constexpr _UIntType
107  linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
108 
109  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
110  constexpr _UIntType
111  linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
112 
113  /**
114  * Seeds the LCR with integral value @p __s, adjusted so that the
115  * ring identity is never a member of the convergence set.
116  */
117  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
118  void
119  linear_congruential_engine<_UIntType, __a, __c, __m>::
120  seed(result_type __s)
121  {
122  if ((__detail::__mod<_UIntType, __m>(__c) == 0)
123  && (__detail::__mod<_UIntType, __m>(__s) == 0))
124  _M_x = 1;
125  else
126  _M_x = __detail::__mod<_UIntType, __m>(__s);
127  }
128 
129  /**
130  * Seeds the LCR engine with a value generated by @p __q.
131  */
132  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
133  template<typename _Sseq>
134  typename std::enable_if<std::is_class<_Sseq>::value>::type
135  linear_congruential_engine<_UIntType, __a, __c, __m>::
136  seed(_Sseq& __q)
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  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
161  typedef typename __ostream_type::ios_base __ios_base;
162 
163  const typename __ios_base::fmtflags __flags = __os.flags();
164  const _CharT __fill = __os.fill();
165  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
166  __os.fill(__os.widen(' '));
167 
168  __os << __lcr._M_x;
169 
170  __os.flags(__flags);
171  __os.fill(__fill);
172  return __os;
173  }
174 
175  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
176  typename _CharT, typename _Traits>
177  std::basic_istream<_CharT, _Traits>&
178  operator>>(std::basic_istream<_CharT, _Traits>& __is,
179  linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
180  {
181  typedef std::basic_istream<_CharT, _Traits> __istream_type;
182  typedef typename __istream_type::ios_base __ios_base;
183 
184  const typename __ios_base::fmtflags __flags = __is.flags();
185  __is.flags(__ios_base::dec);
186 
187  __is >> __lcr._M_x;
188 
189  __is.flags(__flags);
190  return __is;
191  }
192 
193 
194  template<typename _UIntType,
195  size_t __w, size_t __n, size_t __m, size_t __r,
196  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
197  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
198  _UIntType __f>
199  constexpr size_t
200  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
201  __s, __b, __t, __c, __l, __f>::word_size;
202 
203  template<typename _UIntType,
204  size_t __w, size_t __n, size_t __m, size_t __r,
205  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
206  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
207  _UIntType __f>
208  constexpr size_t
209  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
210  __s, __b, __t, __c, __l, __f>::state_size;
211 
212  template<typename _UIntType,
213  size_t __w, size_t __n, size_t __m, size_t __r,
214  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
215  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
216  _UIntType __f>
217  constexpr size_t
218  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
219  __s, __b, __t, __c, __l, __f>::shift_size;
220 
221  template<typename _UIntType,
222  size_t __w, size_t __n, size_t __m, size_t __r,
223  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
224  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
225  _UIntType __f>
226  constexpr size_t
227  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
228  __s, __b, __t, __c, __l, __f>::mask_bits;
229 
230  template<typename _UIntType,
231  size_t __w, size_t __n, size_t __m, size_t __r,
232  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
233  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
234  _UIntType __f>
235  constexpr _UIntType
236  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
237  __s, __b, __t, __c, __l, __f>::xor_mask;
238 
239  template<typename _UIntType,
240  size_t __w, size_t __n, size_t __m, size_t __r,
241  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
242  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
243  _UIntType __f>
244  constexpr size_t
245  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
246  __s, __b, __t, __c, __l, __f>::tempering_u;
247 
248  template<typename _UIntType,
249  size_t __w, size_t __n, size_t __m, size_t __r,
250  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
251  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
252  _UIntType __f>
253  constexpr _UIntType
254  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
255  __s, __b, __t, __c, __l, __f>::tempering_d;
256 
257  template<typename _UIntType,
258  size_t __w, size_t __n, size_t __m, size_t __r,
259  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
260  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
261  _UIntType __f>
262  constexpr size_t
263  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
264  __s, __b, __t, __c, __l, __f>::tempering_s;
265 
266  template<typename _UIntType,
267  size_t __w, size_t __n, size_t __m, size_t __r,
268  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
269  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
270  _UIntType __f>
271  constexpr _UIntType
272  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
273  __s, __b, __t, __c, __l, __f>::tempering_b;
274 
275  template<typename _UIntType,
276  size_t __w, size_t __n, size_t __m, size_t __r,
277  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
278  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
279  _UIntType __f>
280  constexpr size_t
281  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
282  __s, __b, __t, __c, __l, __f>::tempering_t;
283 
284  template<typename _UIntType,
285  size_t __w, size_t __n, size_t __m, size_t __r,
286  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
287  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
288  _UIntType __f>
289  constexpr _UIntType
290  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
291  __s, __b, __t, __c, __l, __f>::tempering_c;
292 
293  template<typename _UIntType,
294  size_t __w, size_t __n, size_t __m, size_t __r,
295  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
296  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
297  _UIntType __f>
298  constexpr size_t
299  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
300  __s, __b, __t, __c, __l, __f>::tempering_l;
301 
302  template<typename _UIntType,
303  size_t __w, size_t __n, size_t __m, size_t __r,
304  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
305  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
306  _UIntType __f>
307  constexpr _UIntType
308  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
309  __s, __b, __t, __c, __l, __f>::
310  initialization_multiplier;
311 
312  template<typename _UIntType,
313  size_t __w, size_t __n, size_t __m, size_t __r,
314  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
315  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
316  _UIntType __f>
317  constexpr _UIntType
318  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
319  __s, __b, __t, __c, __l, __f>::default_seed;
320 
321  template<typename _UIntType,
322  size_t __w, size_t __n, size_t __m, size_t __r,
323  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
324  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
325  _UIntType __f>
326  void
327  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
328  __s, __b, __t, __c, __l, __f>::
329  seed(result_type __sd)
330  {
331  _M_x[0] = __detail::__mod<_UIntType,
332  __detail::_Shift<_UIntType, __w>::__value>(__sd);
333 
334  for (size_t __i = 1; __i < state_size; ++__i)
335  {
336  _UIntType __x = _M_x[__i - 1];
337  __x ^= __x >> (__w - 2);
338  __x *= __f;
339  __x += __detail::__mod<_UIntType, __n>(__i);
340  _M_x[__i] = __detail::__mod<_UIntType,
341  __detail::_Shift<_UIntType, __w>::__value>(__x);
342  }
343  _M_p = state_size;
344  }
345 
346  template<typename _UIntType,
347  size_t __w, size_t __n, size_t __m, size_t __r,
348  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
349  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
350  _UIntType __f>
351  template<typename _Sseq>
352  typename std::enable_if<std::is_class<_Sseq>::value>::type
353  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
354  __s, __b, __t, __c, __l, __f>::
355  seed(_Sseq& __q)
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  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
482  typedef typename __ostream_type::ios_base __ios_base;
483 
484  const typename __ios_base::fmtflags __flags = __os.flags();
485  const _CharT __fill = __os.fill();
486  const _CharT __space = __os.widen(' ');
487  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
488  __os.fill(__space);
489 
490  for (size_t __i = 0; __i < __n; ++__i)
491  __os << __x._M_x[__i] << __space;
492  __os << __x._M_p;
493 
494  __os.flags(__flags);
495  __os.fill(__fill);
496  return __os;
497  }
498 
499  template<typename _UIntType, size_t __w,
500  size_t __n, size_t __m, size_t __r,
501  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
502  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
503  _UIntType __f, typename _CharT, typename _Traits>
504  std::basic_istream<_CharT, _Traits>&
505  operator>>(std::basic_istream<_CharT, _Traits>& __is,
506  mersenne_twister_engine<_UIntType, __w, __n, __m,
507  __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
508  {
509  typedef std::basic_istream<_CharT, _Traits> __istream_type;
510  typedef typename __istream_type::ios_base __ios_base;
511 
512  const typename __ios_base::fmtflags __flags = __is.flags();
513  __is.flags(__ios_base::dec | __ios_base::skipws);
514 
515  for (size_t __i = 0; __i < __n; ++__i)
516  __is >> __x._M_x[__i];
517  __is >> __x._M_p;
518 
519  __is.flags(__flags);
520  return __is;
521  }
522 
523 
524  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
525  constexpr size_t
526  subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
527 
528  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
529  constexpr size_t
530  subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
531 
532  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
533  constexpr size_t
534  subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
535 
536  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
537  constexpr _UIntType
538  subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
539 
540  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
541  void
542  subtract_with_carry_engine<_UIntType, __w, __s, __r>::
543  seed(result_type __value)
544  {
545  std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
546  __lcg(__value == 0u ? default_seed : __value);
547 
548  const size_t __n = (__w + 31) / 32;
549 
550  for (size_t __i = 0; __i < long_lag; ++__i)
551  {
552  _UIntType __sum = 0u;
553  _UIntType __factor = 1u;
554  for (size_t __j = 0; __j < __n; ++__j)
555  {
556  __sum += __detail::__mod<uint_least32_t,
557  __detail::_Shift<uint_least32_t, 32>::__value>
558  (__lcg()) * __factor;
559  __factor *= __detail::_Shift<_UIntType, 32>::__value;
560  }
561  _M_x[__i] = __detail::__mod<_UIntType,
562  __detail::_Shift<_UIntType, __w>::__value>(__sum);
563  }
564  _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
565  _M_p = 0;
566  }
567 
568  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
569  template<typename _Sseq>
570  typename std::enable_if<std::is_class<_Sseq>::value>::type
571  subtract_with_carry_engine<_UIntType, __w, __s, __r>::
572  seed(_Sseq& __q)
573  {
574  const size_t __k = (__w + 31) / 32;
575  uint_least32_t __arr[__r * __k];
576  __q.generate(__arr + 0, __arr + __r * __k);
577 
578  for (size_t __i = 0; __i < long_lag; ++__i)
579  {
580  _UIntType __sum = 0u;
581  _UIntType __factor = 1u;
582  for (size_t __j = 0; __j < __k; ++__j)
583  {
584  __sum += __arr[__k * __i + __j] * __factor;
585  __factor *= __detail::_Shift<_UIntType, 32>::__value;
586  }
587  _M_x[__i] = __detail::__mod<_UIntType,
588  __detail::_Shift<_UIntType, __w>::__value>(__sum);
589  }
590  _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
591  _M_p = 0;
592  }
593 
594  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
595  typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
596  result_type
597  subtract_with_carry_engine<_UIntType, __w, __s, __r>::
598  operator()()
599  {
600  // Derive short lag index from current index.
601  long __ps = _M_p - short_lag;
602  if (__ps < 0)
603  __ps += long_lag;
604 
605  // Calculate new x(i) without overflow or division.
606  // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
607  // cannot overflow.
608  _UIntType __xi;
609  if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
610  {
611  __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
612  _M_carry = 0;
613  }
614  else
615  {
616  __xi = (__detail::_Shift<_UIntType, __w>::__value
617  - _M_x[_M_p] - _M_carry + _M_x[__ps]);
618  _M_carry = 1;
619  }
620  _M_x[_M_p] = __xi;
621 
622  // Adjust current index to loop around in ring buffer.
623  if (++_M_p >= long_lag)
624  _M_p = 0;
625 
626  return __xi;
627  }
628 
629  template<typename _UIntType, size_t __w, size_t __s, size_t __r,
630  typename _CharT, typename _Traits>
631  std::basic_ostream<_CharT, _Traits>&
632  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
633  const subtract_with_carry_engine<_UIntType,
634  __w, __s, __r>& __x)
635  {
636  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
637  typedef typename __ostream_type::ios_base __ios_base;
638 
639  const typename __ios_base::fmtflags __flags = __os.flags();
640  const _CharT __fill = __os.fill();
641  const _CharT __space = __os.widen(' ');
642  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
643  __os.fill(__space);
644 
645  for (size_t __i = 0; __i < __r; ++__i)
646  __os << __x._M_x[__i] << __space;
647  __os << __x._M_carry << __space << __x._M_p;
648 
649  __os.flags(__flags);
650  __os.fill(__fill);
651  return __os;
652  }
653 
654  template<typename _UIntType, size_t __w, size_t __s, size_t __r,
655  typename _CharT, typename _Traits>
656  std::basic_istream<_CharT, _Traits>&
657  operator>>(std::basic_istream<_CharT, _Traits>& __is,
658  subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
659  {
660  typedef std::basic_ostream<_CharT, _Traits> __istream_type;
661  typedef typename __istream_type::ios_base __ios_base;
662 
663  const typename __ios_base::fmtflags __flags = __is.flags();
664  __is.flags(__ios_base::dec | __ios_base::skipws);
665 
666  for (size_t __i = 0; __i < __r; ++__i)
667  __is >> __x._M_x[__i];
668  __is >> __x._M_carry;
669  __is >> __x._M_p;
670 
671  __is.flags(__flags);
672  return __is;
673  }
674 
675 
676  template<typename _RandomNumberEngine, size_t __p, size_t __r>
677  constexpr size_t
678  discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
679 
680  template<typename _RandomNumberEngine, size_t __p, size_t __r>
681  constexpr size_t
682  discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
683 
684  template<typename _RandomNumberEngine, size_t __p, size_t __r>
685  typename discard_block_engine<_RandomNumberEngine,
686  __p, __r>::result_type
687  discard_block_engine<_RandomNumberEngine, __p, __r>::
688  operator()()
689  {
690  if (_M_n >= used_block)
691  {
692  _M_b.discard(block_size - _M_n);
693  _M_n = 0;
694  }
695  ++_M_n;
696  return _M_b();
697  }
698 
699  template<typename _RandomNumberEngine, size_t __p, size_t __r,
700  typename _CharT, typename _Traits>
701  std::basic_ostream<_CharT, _Traits>&
702  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
703  const discard_block_engine<_RandomNumberEngine,
704  __p, __r>& __x)
705  {
706  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
707  typedef typename __ostream_type::ios_base __ios_base;
708 
709  const typename __ios_base::fmtflags __flags = __os.flags();
710  const _CharT __fill = __os.fill();
711  const _CharT __space = __os.widen(' ');
712  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
713  __os.fill(__space);
714 
715  __os << __x.base() << __space << __x._M_n;
716 
717  __os.flags(__flags);
718  __os.fill(__fill);
719  return __os;
720  }
721 
722  template<typename _RandomNumberEngine, size_t __p, size_t __r,
723  typename _CharT, typename _Traits>
724  std::basic_istream<_CharT, _Traits>&
725  operator>>(std::basic_istream<_CharT, _Traits>& __is,
726  discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
727  {
728  typedef std::basic_istream<_CharT, _Traits> __istream_type;
729  typedef typename __istream_type::ios_base __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);<