libstdc++
bits/random.tcc
Go to the documentation of this file.
1 // random number generation (out of line) -*- C++ -*-
2 
3 // Copyright (C) 2009-2013 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
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>
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>
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();
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>
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>
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>
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(' ');
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>
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();
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  {
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>
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>
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(' ');
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>
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();
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
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>
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(' ');
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>
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();
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
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 
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 
811  template<typename _RandomNumberEngine, size_t __k>
812  constexpr size_t
814 
815  template<typename _RandomNumberEngine, size_t __k>
819  {
820  size_t __j = __k * ((_M_y - _M_b.min())
821  / (_M_b.max() - _M_b.min() + 1.0L));
822  _M_y = _M_v[__j];
823  _M_v[__j] = _M_b();
824 
825  return _M_y;
826  }
827 
828  template<typename _RandomNumberEngine, size_t __k,
829  typename _CharT, typename _Traits>
831  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
833  {
834  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
835  typedef typename __ostream_type::ios_base __ios_base;
836 
837  const typename __ios_base::fmtflags __flags = __os.flags();
838  const _CharT __fill = __os.fill();
839  const _CharT __space = __os.widen(' ');
841  __os.fill(__space);
842 
843  __os << __x.base();
844  for (size_t __i = 0; __i < __k; ++__i)
845  __os << __space << __x._M_v[__i];
846  __os << __space << __x._M_y;
847 
848  __os.flags(__flags);
849  __os.fill(__fill);
850  return __os;
851  }
852 
853  template<typename _RandomNumberEngine, size_t __k,
854  typename _CharT, typename _Traits>
857  shuffle_order_engine<_RandomNumberEngine, __k>& __x)
858  {
859  typedef std::basic_istream<_CharT, _Traits> __istream_type;
860  typedef typename __istream_type::ios_base __ios_base;
861 
862  const typename __ios_base::fmtflags __flags = __is.flags();
864 
865  __is >> __x._M_b;
866  for (size_t __i = 0; __i < __k; ++__i)
867  __is >> __x._M_v[__i];
868  __is >> __x._M_y;
869 
870  __is.flags(__flags);
871  return __is;
872  }
873 
874 
875  template<typename _IntType>
876  template<typename _UniformRandomNumberGenerator>
877  typename uniform_int_distribution<_IntType>::result_type
879  operator()(_UniformRandomNumberGenerator& __urng,
880  const param_type& __param)
881  {
882  typedef typename _UniformRandomNumberGenerator::result_type
883  _Gresult_type;
884  typedef typename std::make_unsigned<result_type>::type __utype;
886  __uctype;
887 
888  const __uctype __urngmin = __urng.min();
889  const __uctype __urngmax = __urng.max();
890  const __uctype __urngrange = __urngmax - __urngmin;
891  const __uctype __urange
892  = __uctype(__param.b()) - __uctype(__param.a());
893 
894  __uctype __ret;
895 
896  if (__urngrange > __urange)
897  {
898  // downscaling
899  const __uctype __uerange = __urange + 1; // __urange can be zero
900  const __uctype __scaling = __urngrange / __uerange;
901  const __uctype __past = __uerange * __scaling;
902  do
903  __ret = __uctype(__urng()) - __urngmin;
904  while (__ret >= __past);
905  __ret /= __scaling;
906  }
907  else if (__urngrange < __urange)
908  {
909  // upscaling
910  /*
911  Note that every value in [0, urange]
912  can be written uniquely as
913 
914  (urngrange + 1) * high + low
915 
916  where
917 
918  high in [0, urange / (urngrange + 1)]
919 
920  and
921 
922  low in [0, urngrange].
923  */
924  __uctype __tmp; // wraparound control
925  do
926  {
927  const __uctype __uerngrange = __urngrange + 1;
928  __tmp = (__uerngrange * operator()
929  (__urng, param_type(0, __urange / __uerngrange)));
930  __ret = __tmp + (__uctype(__urng()) - __urngmin);
931  }
932  while (__ret > __urange || __ret < __tmp);
933  }
934  else
935  __ret = __uctype(__urng()) - __urngmin;
936 
937  return __ret + __param.a();
938  }
939 
940 
941  template<typename _IntType>
942  template<typename _ForwardIterator,
943  typename _UniformRandomNumberGenerator>
944  void
945  uniform_int_distribution<_IntType>::
946  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
947  _UniformRandomNumberGenerator& __urng,
948  const param_type& __param)
949  {
950  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
951  typedef typename _UniformRandomNumberGenerator::result_type
952  _Gresult_type;
953  typedef typename std::make_unsigned<result_type>::type __utype;
954  typedef typename std::common_type<_Gresult_type, __utype>::type
955  __uctype;
956 
957  const __uctype __urngmin = __urng.min();
958  const __uctype __urngmax = __urng.max();
959  const __uctype __urngrange = __urngmax - __urngmin;
960  const __uctype __urange
961  = __uctype(__param.b()) - __uctype(__param.a());
962 
963  __uctype __ret;
964 
965  if (__urngrange > __urange)
966  {
967  if (__detail::_Power_of_2(__urngrange + 1)
968  && __detail::_Power_of_2(__urange + 1))
969  {
970  while (__f != __t)
971  {
972  __ret = __uctype(__urng()) - __urngmin;
973  *__f++ = (__ret & __urange) + __param.a();
974  }
975  }
976  else
977  {
978  // downscaling
979  const __uctype __uerange = __urange + 1; // __urange can be zero
980  const __uctype __scaling = __urngrange / __uerange;
981  const __uctype __past = __uerange * __scaling;
982  while (__f != __t)
983  {
984  do
985  __ret = __uctype(__urng()) - __urngmin;
986  while (__ret >= __past);
987  *__f++ = __ret / __scaling + __param.a();
988  }
989  }
990  }
991  else if (__urngrange < __urange)
992  {
993  // upscaling
994  /*
995  Note that every value in [0, urange]
996  can be written uniquely as
997 
998  (urngrange + 1) * high + low
999 
1000  where
1001 
1002  high in [0, urange / (urngrange + 1)]
1003 
1004  and
1005 
1006  low in [0, urngrange].
1007  */
1008  __uctype __tmp; // wraparound control
1009  while (__f != __t)
1010  {
1011  do
1012  {
1013  const __uctype __uerngrange = __urngrange + 1;
1014  __tmp = (__uerngrange * operator()
1015  (__urng, param_type(0, __urange / __uerngrange)));
1016  __ret = __tmp + (__uctype(__urng()) - __urngmin);
1017  }
1018  while (__ret > __urange || __ret < __tmp);
1019  *__f++ = __ret;
1020  }
1021  }
1022  else
1023  while (__f != __t)
1024  *__f++ = __uctype(__urng()) - __urngmin + __param.a();
1025  }
1026 
1027  template<typename _IntType, typename _CharT, typename _Traits>
1029  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1031  {
1032  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1033  typedef typename __ostream_type::ios_base __ios_base;
1034 
1035  const typename __ios_base::fmtflags __flags = __os.flags();
1036  const _CharT __fill = __os.fill();
1037  const _CharT __space = __os.widen(' ');
1039  __os.fill(__space);
1040 
1041  __os << __x.a() << __space << __x.b();
1042 
1043  __os.flags(__flags);
1044  __os.fill(__fill);
1045  return __os;
1046  }
1047 
1048  template<typename _IntType, typename _CharT, typename _Traits>
1052  {
1053  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1054  typedef typename __istream_type::ios_base __ios_base;
1055 
1056  const typename __ios_base::fmtflags __flags = __is.flags();
1058 
1059  _IntType __a, __b;
1060  __is >> __a >> __b;
1062  param_type(__a, __b));
1063 
1064  __is.flags(__flags);
1065  return __is;
1066  }
1067 
1068 
1069  template<typename _RealType>
1070  template<typename _ForwardIterator,
1071  typename _UniformRandomNumberGenerator>
1072  void
1073  uniform_real_distribution<_RealType>::
1074  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1075  _UniformRandomNumberGenerator& __urng,
1076  const param_type& __p)
1077  {
1078  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1079  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1080  __aurng(__urng);
1081  auto __range = __p.b() - __p.a();
1082  while (__f != __t)
1083  *__f++ = __aurng() * __range + __p.a();
1084  }
1085 
1086  template<typename _RealType, typename _CharT, typename _Traits>
1087  std::basic_ostream<_CharT, _Traits>&
1088  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1089  const uniform_real_distribution<_RealType>& __x)
1090  {
1091  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1092  typedef typename __ostream_type::ios_base __ios_base;
1093 
1094  const typename __ios_base::fmtflags __flags = __os.flags();
1095  const _CharT __fill = __os.fill();
1096  const std::streamsize __precision = __os.precision();
1097  const _CharT __space = __os.widen(' ');
1099  __os.fill(__space);
1101 
1102  __os << __x.a() << __space << __x.b();
1103 
1104  __os.flags(__flags);
1105  __os.fill(__fill);
1106  __os.precision(__precision);
1107  return __os;
1108  }
1109 
1110  template<typename _RealType, typename _CharT, typename _Traits>
1114  {
1115  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1116  typedef typename __istream_type::ios_base __ios_base;
1117 
1118  const typename __ios_base::fmtflags __flags = __is.flags();
1119  __is.flags(__ios_base::skipws);
1120 
1121  _RealType __a, __b;
1122  __is >> __a >> __b;
1124  param_type(__a, __b));
1125 
1126  __is.flags(__flags);
1127  return __is;
1128  }
1129 
1130 
1131  template<typename _ForwardIterator,
1132  typename _UniformRandomNumberGenerator>
1133  void
1134  std::bernoulli_distribution::
1135  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1136  _UniformRandomNumberGenerator& __urng,
1137  const param_type& __p)
1138  {
1139  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1140  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1141  __aurng(__urng);
1142  auto __limit = __p.p() * (__aurng.max() - __aurng.min());
1143 
1144  while (__f != __t)
1145  *__f++ = (__aurng() - __aurng.min()) < __limit;
1146  }
1147 
1148  template<typename _CharT, typename _Traits>
1149  std::basic_ostream<_CharT, _Traits>&
1150  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1151  const bernoulli_distribution& __x)
1152  {
1153  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1154  typedef typename __ostream_type::ios_base __ios_base;
1155 
1156  const typename __ios_base::fmtflags __flags = __os.flags();
1157  const _CharT __fill = __os.fill();
1158  const std::streamsize __precision = __os.precision();
1160  __os.fill(__os.widen(' '));
1162 
1163  __os << __x.p();
1164 
1165  __os.flags(__flags);
1166  __os.fill(__fill);
1167  __os.precision(__precision);
1168  return __os;
1169  }
1170 
1171 
1172  template<typename _IntType>
1173  template<typename _UniformRandomNumberGenerator>
1174  typename geometric_distribution<_IntType>::result_type
1176  operator()(_UniformRandomNumberGenerator& __urng,
1177  const param_type& __param)
1178  {
1179  // About the epsilon thing see this thread:
1180  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1181  const double __naf =
1183  // The largest _RealType convertible to _IntType.
1184  const double __thr =
1186  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1187  __aurng(__urng);
1188 
1189  double __cand;
1190  do
1191  __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1192  while (__cand >= __thr);
1193 
1194  return result_type(__cand + __naf);
1195  }
1196 
1197  template<typename _IntType>
1198  template<typename _ForwardIterator,
1199  typename _UniformRandomNumberGenerator>
1200  void
1201  geometric_distribution<_IntType>::
1202  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1203  _UniformRandomNumberGenerator& __urng,
1204  const param_type& __param)
1205  {
1206  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1207  // About the epsilon thing see this thread:
1208  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1209  const double __naf =
1210  (1 - std::numeric_limits<double>::epsilon()) / 2;
1211  // The largest _RealType convertible to _IntType.
1212  const double __thr =
1213  std::numeric_limits<_IntType>::max() + __naf;
1214  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1215  __aurng(__urng);
1216 
1217  while (__f != __t)
1218  {
1219  double __cand;
1220  do
1221  __cand = std::floor(std::log(1.0 - __aurng())
1222  / __param._M_log_1_p);
1223  while (__cand >= __thr);
1224 
1225  *__f++ = __cand + __naf;
1226  }
1227  }
1228 
1229  template<typename _IntType,
1230  typename _CharT, typename _Traits>
1232  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1234  {
1235  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1236  typedef typename __ostream_type::ios_base __ios_base;
1237 
1238  const typename __ios_base::fmtflags __flags = __os.flags();
1239  const _CharT __fill = __os.fill();
1240  const std::streamsize __precision = __os.precision();
1242  __os.fill(__os.widen(' '));
1244 
1245  __os << __x.p();
1246 
1247  __os.flags(__flags);
1248  __os.fill(__fill);
1249  __os.precision(__precision);
1250  return __os;
1251  }
1252 
1253  template<typename _IntType,
1254  typename _CharT, typename _Traits>
1258  {
1259  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1260  typedef typename __istream_type::ios_base __ios_base;
1261 
1262  const typename __ios_base::fmtflags __flags = __is.flags();
1263  __is.flags(__ios_base::skipws);
1264 
1265  double __p;
1266  __is >> __p;
1268 
1269  __is.flags(__flags);
1270  return __is;
1271  }
1272 
1273  // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1274  template<typename _IntType>
1275  template<typename _UniformRandomNumberGenerator>
1276  typename negative_binomial_distribution<_IntType>::result_type
1278  operator()(_UniformRandomNumberGenerator& __urng)
1279  {
1280  const double __y = _M_gd(__urng);
1281 
1282  // XXX Is the constructor too slow?
1284  return __poisson(__urng);
1285  }
1286 
1287  template<typename _IntType>
1288  template<typename _UniformRandomNumberGenerator>
1289  typename negative_binomial_distribution<_IntType>::result_type
1291  operator()(_UniformRandomNumberGenerator& __urng,
1292  const param_type& __p)
1293  {
1295  param_type;
1296 
1297  const double __y =
1298  _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1299 
1301  return __poisson(__urng);
1302  }
1303 
1304  template<typename _IntType>
1305  template<typename _ForwardIterator,
1306  typename _UniformRandomNumberGenerator>
1307  void
1308  negative_binomial_distribution<_IntType>::
1309  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1310  _UniformRandomNumberGenerator& __urng)
1311  {
1312  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1313  while (__f != __t)
1314  {
1315  const double __y = _M_gd(__urng);
1316 
1317  // XXX Is the constructor too slow?
1319  *__f++ = __poisson(__urng);
1320  }
1321  }
1322 
1323  template<typename _IntType>
1324  template<typename _ForwardIterator,
1325  typename _UniformRandomNumberGenerator>
1326  void
1327  negative_binomial_distribution<_IntType>::
1328  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1329  _UniformRandomNumberGenerator& __urng,
1330  const param_type& __p)
1331  {
1332  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1333  typename std::gamma_distribution<result_type>::param_type
1334  __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1335 
1336  while (__f != __t)
1337  {
1338  const double __y = _M_gd(__urng, __p2);
1339 
1341  *__f++ = __poisson(__urng);
1342  }
1343  }
1344 
1345  template<typename _IntType, typename _CharT, typename _Traits>
1347  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1348  const negative_binomial_distribution<_IntType>& __x)
1349  {
1350  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1351  typedef typename __ostream_type::ios_base __ios_base;
1352 
1353  const typename __ios_base::fmtflags __flags = __os.flags();
1354  const _CharT __fill = __os.fill();
1355  const std::streamsize __precision = __os.precision();
1356  const _CharT __space = __os.widen(' ');
1358  __os.fill(__os.widen(' '));
1360 
1361  __os << __x.k() << __space << __x.p()
1362  << __space << __x._M_gd;
1363 
1364  __os.flags(__flags);
1365  __os.fill(__fill);
1366  __os.precision(__precision);
1367  return __os;
1368  }
1369 
1370  template<typename _IntType, typename _CharT, typename _Traits>
1373  negative_binomial_distribution<_IntType>& __x)
1374  {
1375  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1376  typedef typename __istream_type::ios_base __ios_base;
1377 
1378  const typename __ios_base::fmtflags __flags = __is.flags();
1379  __is.flags(__ios_base::skipws);
1380 
1381  _IntType __k;
1382  double __p;
1383  __is >> __k >> __p >> __x._M_gd;
1384  __x.param(typename negative_binomial_distribution<_IntType>::
1385  param_type(__k, __p));
1386 
1387  __is.flags(__flags);
1388  return __is;
1389  }
1390 
1391 
1392  template<typename _IntType>
1393  void
1394  poisson_distribution<_IntType>::param_type::
1395  _M_initialize()
1396  {
1397 #if _GLIBCXX_USE_C99_MATH_TR1
1398  if (_M_mean >= 12)
1399  {
1400  const double __m = std::floor(_M_mean);
1401  _M_lm_thr = std::log(_M_mean);
1402  _M_lfm = std::lgamma(__m + 1);
1403  _M_sm = std::sqrt(__m);
1404 
1405  const double __pi_4 = 0.7853981633974483096156608458198757L;
1406  const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1407  / __pi_4));
1408  _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
1409  const double __cx = 2 * __m + _M_d;
1410  _M_scx = std::sqrt(__cx / 2);
1411  _M_1cx = 1 / __cx;
1412 
1413  _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1414  _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1415  / _M_d;
1416  }
1417  else
1418 #endif
1419  _M_lm_thr = std::exp(-_M_mean);
1420  }
1421 
1422  /**
1423  * A rejection algorithm when mean >= 12 and a simple method based
1424  * upon the multiplication of uniform random variates otherwise.
1425  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1426  * is defined.
1427  *
1428  * Reference:
1429  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1430  * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1431  */
1432  template<typename _IntType>
1433  template<typename _UniformRandomNumberGenerator>
1434  typename poisson_distribution<_IntType>::result_type
1436  operator()(_UniformRandomNumberGenerator& __urng,
1437  const param_type& __param)
1438  {
1439  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1440  __aurng(__urng);
1441 #if _GLIBCXX_USE_C99_MATH_TR1
1442  if (__param.mean() >= 12)
1443  {
1444  double __x;
1445 
1446  // See comments above...
1447  const double __naf =
1449  const double __thr =
1451 
1452  const double __m = std::floor(__param.mean());
1453  // sqrt(pi / 2)
1454  const double __spi_2 = 1.2533141373155002512078826424055226L;
1455  const double __c1 = __param._M_sm * __spi_2;
1456  const double __c2 = __param._M_c2b + __c1;
1457  const double __c3 = __c2 + 1;
1458  const double __c4 = __c3 + 1;
1459  // e^(1 / 78)
1460  const double __e178 = 1.0129030479320018583185514777512983L;
1461  const double __c5 = __c4 + __e178;
1462  const double __c = __param._M_cb + __c5;
1463  const double __2cx = 2 * (2 * __m + __param._M_d);
1464 
1465  bool __reject = true;
1466  do
1467  {
1468  const double __u = __c * __aurng();
1469  const double __e = -std::log(1.0 - __aurng());
1470 
1471  double __w = 0.0;
1472 
1473  if (__u <= __c1)
1474  {
1475  const double __n = _M_nd(__urng);
1476  const double __y = -std::abs(__n) * __param._M_sm - 1;
1477  __x = std::floor(__y);
1478  __w = -__n * __n / 2;
1479  if (__x < -__m)
1480  continue;
1481  }
1482  else if (__u <= __c2)
1483  {
1484  const double __n = _M_nd(__urng);
1485  const double __y = 1 + std::abs(__n) * __param._M_scx;
1486  __x = std::ceil(__y);
1487  __w = __y * (2 - __y) * __param._M_1cx;
1488  if (__x > __param._M_d)
1489  continue;
1490  }
1491  else if (__u <= __c3)
1492  // NB: This case not in the book, nor in the Errata,
1493  // but should be ok...
1494  __x = -1;
1495  else if (__u <= __c4)
1496  __x = 0;
1497  else if (__u <= __c5)
1498  __x = 1;
1499  else
1500  {
1501  const double __v = -std::log(1.0 - __aurng());
1502  const double __y = __param._M_d
1503  + __v * __2cx / __param._M_d;
1504  __x = std::ceil(__y);
1505  __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1506  }
1507 
1508  __reject = (__w - __e - __x * __param._M_lm_thr
1509  > __param._M_lfm - std::lgamma(__x + __m + 1));
1510 
1511  __reject |= __x + __m >= __thr;
1512 
1513  } while (__reject);
1514 
1515  return result_type(__x + __m + __naf);
1516  }
1517  else
1518 #endif
1519  {
1520  _IntType __x = 0;
1521  double __prod = 1.0;
1522 
1523  do
1524  {
1525  __prod *= __aurng();
1526  __x += 1;
1527  }
1528  while (__prod > __param._M_lm_thr);
1529 
1530  return __x - 1;
1531  }
1532  }
1533 
1534  template<typename _IntType>
1535  template<typename _ForwardIterator,
1536  typename _UniformRandomNumberGenerator>
1537  void
1539  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1540  _UniformRandomNumberGenerator& __urng,
1541  const param_type& __param)
1542  {
1543  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1544  // We could duplicate everything from operator()...
1545  while (__f != __t)
1546  *__f++ = this->operator()(__urng, __param);
1547  }
1548 
1549  template<typename _IntType,
1550  typename _CharT, typename _Traits>
1551  std::basic_ostream<_CharT, _Traits>&
1552  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1553  const poisson_distribution<_IntType>& __x)
1554  {
1555  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1556  typedef typename __ostream_type::ios_base __ios_base;
1557 
1558  const typename __ios_base::fmtflags __flags = __os.flags();
1559  const _CharT __fill = __os.fill();
1560  const std::streamsize __precision = __os.precision();
1561  const _CharT __space = __os.widen(' ');
1563  __os.fill(__space);
1565 
1566  __os << __x.mean() << __space << __x._M_nd;
1567 
1568  __os.flags(__flags);
1569  __os.fill(__fill);
1570  __os.precision(__precision);
1571  return __os;
1572  }
1573 
1574  template<typename _IntType,
1575  typename _CharT, typename _Traits>
1578  poisson_distribution<_IntType>& __x)
1579  {
1580  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1581  typedef typename __istream_type::ios_base __ios_base;
1582 
1583  const typename __ios_base::fmtflags __flags = __is.flags();
1584  __is.flags(__ios_base::skipws);
1585 
1586  double __mean;
1587  __is >> __mean >> __x._M_nd;
1588  __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1589 
1590  __is.flags(__flags);
1591  return __is;
1592  }
1593 
1594 
1595  template<typename _IntType>
1596  void
1597  binomial_distribution<_IntType>::param_type::
1598  _M_initialize()
1599  {
1600  const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1601 
1602  _M_easy = true;
1603 
1604 #if _GLIBCXX_USE_C99_MATH_TR1
1605  if (_M_t * __p12 >= 8)
1606  {
1607  _M_easy = false;
1608  const double __np = std::floor(_M_t * __p12);
1609  const double __pa = __np / _M_t;
1610  const double __1p = 1 - __pa;
1611 
1612  const double __pi_4 = 0.7853981633974483096156608458198757L;
1613  const double __d1x =
1614  std::sqrt(__np * __1p * std::log(32 * __np
1615  / (81 * __pi_4 * __1p)));
1616  _M_d1 = std::round(std::max(1.0, __d1x));
1617  const double __d2x =
1618  std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1619  / (__pi_4 * __pa)));
1620  _M_d2 = std::round(std::max(1.0, __d2x));
1621 
1622  // sqrt(pi / 2)
1623  const double __spi_2 = 1.2533141373155002512078826424055226L;
1624  _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1625  _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1626  _M_c = 2 * _M_d1 / __np;
1627  _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1628  const double __a12 = _M_a1 + _M_s2 * __spi_2;
1629  const double __s1s = _M_s1 * _M_s1;
1630  _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1631  * 2 * __s1s / _M_d1
1632  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1633  const double __s2s = _M_s2 * _M_s2;
1634  _M_s = (_M_a123 + 2 * __s2s / _M_d2
1635  * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1636  _M_lf = (std::lgamma(__np + 1)
1637  + std::lgamma(_M_t - __np + 1));
1638  _M_lp1p = std::log(__pa / __1p);
1639 
1640  _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1641  }
1642  else
1643 #endif
1644  _M_q = -std::log(1 - __p12);
1645  }
1646 
1647  template<typename _IntType>
1648  template<typename _UniformRandomNumberGenerator>
1649  typename binomial_distribution<_IntType>::result_type
1650  binomial_distribution<_IntType>::
1651  _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1652  {
1653  _IntType __x = 0;
1654  double __sum = 0.0;
1655  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1656  __aurng(__urng);
1657 
1658  do
1659  {
1660  if (__t == __x)
1661  return __x;
1662  const double __e = -std::log(1.0 - __aurng());
1663  __sum += __e / (__t - __x);
1664  __x += 1;
1665  }
1666  while (__sum <= _M_param._M_q);
1667 
1668  return __x - 1;
1669  }
1670 
1671  /**
1672  * A rejection algorithm when t * p >= 8 and a simple waiting time
1673  * method - the second in the referenced book - otherwise.
1674  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1675  * is defined.
1676  *
1677  * Reference:
1678  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1679  * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1680  */
1681  template<typename _IntType>
1682  template<typename _UniformRandomNumberGenerator>
1683  typename binomial_distribution<_IntType>::result_type
1685  operator()(_UniformRandomNumberGenerator& __urng,
1686  const param_type& __param)
1687  {
1688  result_type __ret;
1689  const _IntType __t = __param.t();
1690  const double __p = __param.p();
1691  const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1692  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1693  __aurng(__urng);
1694 
1695 #if _GLIBCXX_USE_C99_MATH_TR1
1696  if (!__param._M_easy)
1697  {
1698  double __x;
1699 
1700  // See comments above...
1701  const double __naf =
1703  const double __thr =
1705 
1706  const double __np = std::floor(__t * __p12);
1707 
1708  // sqrt(pi / 2)
1709  const double __spi_2 = 1.2533141373155002512078826424055226L;
1710  const double __a1 = __param._M_a1;
1711  const double __a12 = __a1 + __param._M_s2 * __spi_2;
1712  const double __a123 = __param._M_a123;
1713  const double __s1s = __param._M_s1 * __param._M_s1;
1714  const double __s2s = __param._M_s2 * __param._M_s2;
1715 
1716  bool __reject;
1717  do
1718  {
1719  const double __u = __param._M_s * __aurng();
1720 
1721  double __v;
1722 
1723  if (__u <= __a1)
1724  {
1725  const double __n = _M_nd(__urng);
1726  const double __y = __param._M_s1 * std::abs(__n);
1727  __reject = __y >= __param._M_d1;
1728  if (!__reject)
1729  {
1730  const double __e = -std::log(1.0 - __aurng());
1731  __x = std::floor(__y);
1732  __v = -__e - __n * __n / 2 + __param._M_c;
1733  }
1734  }
1735  else if (__u <= __a12)
1736  {
1737  const double __n = _M_nd(__urng);
1738  const double __y = __param._M_s2 * std::abs(__n);
1739  __reject = __y >= __param._M_d2;
1740  if (!__reject)
1741  {
1742  const double __e = -std::log(1.0 - __aurng());
1743  __x = std::floor(-__y);
1744  __v = -__e - __n * __n / 2;
1745  }
1746  }
1747  else if (__u <= __a123)
1748  {
1749  const double __e1 = -std::log(1.0 - __aurng());
1750  const double __e2 = -std::log(1.0 - __aurng());
1751 
1752  const double __y = __param._M_d1
1753  + 2 * __s1s * __e1 / __param._M_d1;
1754  __x = std::floor(__y);
1755  __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1756  -__y / (2 * __s1s)));
1757  __reject = false;
1758  }
1759  else
1760  {
1761  const double __e1 = -std::log(1.0 - __aurng());
1762  const double __e2 = -std::log(1.0 - __aurng());
1763 
1764  const double __y = __param._M_d2
1765  + 2 * __s2s * __e1 / __param._M_d2;
1766  __x = std::floor(-__y);
1767  __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1768  __reject = false;
1769  }
1770 
1771  __reject = __reject || __x < -__np || __x > __t - __np;
1772  if (!__reject)
1773  {
1774  const double __lfx =
1775  std::lgamma(__np + __x + 1)
1776  + std::lgamma(__t - (__np + __x) + 1);
1777  __reject = __v > __param._M_lf - __lfx
1778  + __x * __param._M_lp1p;
1779  }
1780 
1781  __reject |= __x + __np >= __thr;
1782  }
1783  while (__reject);
1784 
1785  __x += __np + __naf;
1786 
1787  const _IntType __z = _M_waiting(__urng, __t - _IntType(__x));
1788  __ret = _IntType(__x) + __z;
1789  }
1790  else
1791 #endif
1792  __ret = _M_waiting(__urng, __t);
1793 
1794  if (__p12 != __p)
1795  __ret = __t - __ret;
1796  return __ret;
1797  }
1798 
1799  template<typename _IntType>
1800  template<typename _ForwardIterator,
1801  typename _UniformRandomNumberGenerator>
1802  void
1804  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1805  _UniformRandomNumberGenerator& __urng,
1806  const param_type& __param)
1807  {
1808  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1809  // We could duplicate everything from operator()...
1810  while (__f != __t)
1811  *__f++ = this->operator()(__urng, __param);
1812  }
1813 
1814  template<typename _IntType,
1815  typename _CharT, typename _Traits>
1816  std::basic_ostream<_CharT, _Traits>&
1817  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1818  const binomial_distribution<_IntType>& __x)
1819  {
1820  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1821  typedef typename __ostream_type::ios_base __ios_base;
1822 
1823  const typename __ios_base::fmtflags __flags = __os.flags();
1824  const _CharT __fill = __os.fill();
1825  const std::streamsize __precision = __os.precision();
1826  const _CharT __space = __os.widen(' ');
1828  __os.fill(__space);
1830 
1831  __os << __x.t() << __space << __x.p()
1832  << __space << __x._M_nd;
1833 
1834  __os.flags(__flags);
1835  __os.fill(__fill);
1836  __os.precision(__precision);
1837  return __os;
1838  }
1839 
1840  template<typename _IntType,
1841  typename _CharT, typename _Traits>
1844  binomial_distribution<_IntType>& __x)
1845  {
1846  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1847  typedef typename __istream_type::ios_base __ios_base;
1848 
1849  const typename __ios_base::fmtflags __flags = __is.flags();
1851 
1852  _IntType __t;
1853  double __p;
1854  __is >> __t >> __p >> __x._M_nd;
1855  __x.param(typename binomial_distribution<_IntType>::
1856  param_type(__t, __p));
1857 
1858  __is.flags(__flags);
1859  return __is;
1860  }
1861 
1862 
1863  template<typename _RealType>
1864  template<typename _ForwardIterator,
1865  typename _UniformRandomNumberGenerator>
1866  void
1868  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1869  _UniformRandomNumberGenerator& __urng,
1870  const param_type& __p)
1871  {
1872  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1873  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1874  __aurng(__urng);
1875  while (__f != __t)
1876  *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1877  }
1878 
1879  template<typename _RealType, typename _CharT, typename _Traits>
1880  std::basic_ostream<_CharT, _Traits>&
1881  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1882  const exponential_distribution<_RealType>& __x)
1883  {
1884  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1885  typedef typename __ostream_type::ios_base __ios_base;
1886 
1887  const typename __ios_base::fmtflags __flags = __os.flags();
1888  const _CharT __fill = __os.fill();
1889  const std::streamsize __precision = __os.precision();
1891  __os.fill(__os.widen(' '));
1893 
1894  __os << __x.lambda();
1895 
1896  __os.flags(__flags);
1897  __os.fill(__fill);
1898  __os.precision(__precision);
1899  return __os;
1900  }
1901 
1902  template<typename _RealType, typename _CharT, typename _Traits>
1906  {
1907  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1908  typedef typename __istream_type::ios_base __ios_base;
1909 
1910  const typename __ios_base::fmtflags __flags = __is.flags();
1912 
1913  _RealType __lambda;
1914  __is >> __lambda;
1916  param_type(__lambda));
1917 
1918  __is.flags(__flags);
1919  return __is;
1920  }
1921 
1922 
1923  /**
1924  * Polar method due to Marsaglia.
1925  *
1926  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1927  * New York, 1986, Ch. V, Sect. 4.4.
1928  */
1929  template<typename _RealType>
1930  template<typename _UniformRandomNumberGenerator>
1931  typename normal_distribution<_RealType>::result_type
1933  operator()(_UniformRandomNumberGenerator& __urng,
1934  const param_type& __param)
1935  {
1936  result_type __ret;
1937  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1938  __aurng(__urng);
1939 
1940  if (_M_saved_available)
1941  {
1942  _M_saved_available = false;
1943  __ret = _M_saved;
1944  }
1945  else
1946  {
1947  result_type __x, __y, __r2;
1948  do
1949  {
1950  __x = result_type(2.0) * __aurng() - 1.0;
1951  __y = result_type(2.0) * __aurng() - 1.0;
1952  __r2 = __x * __x + __y * __y;
1953  }
1954  while (__r2 > 1.0 || __r2 == 0.0);
1955 
1956  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1957  _M_saved = __x * __mult;
1958  _M_saved_available = true;
1959  __ret = __y * __mult;
1960  }
1961 
1962  __ret = __ret * __param.stddev() + __param.mean();
1963  return __ret;
1964  }
1965 
1966  template<typename _RealType>
1967  template<typename _ForwardIterator,
1968  typename _UniformRandomNumberGenerator>
1969  void
1971  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1972  _UniformRandomNumberGenerator& __urng,
1973  const param_type& __param)
1974  {
1975  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1976 
1977  if (__f == __t)
1978  return;
1979 
1980  if (_M_saved_available)
1981  {
1982  _M_saved_available = false;
1983  *__f++ = _M_saved * __param.stddev() + __param.mean();
1984 
1985  if (__f == __t)
1986  return;
1987  }
1988 
1989  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1990  __aurng(__urng);
1991 
1992  while (__f + 1 < __t)
1993  {
1994  result_type __x, __y, __r2;
1995  do
1996  {
1997  __x = result_type(2.0) * __aurng() - 1.0;
1998  __y = result_type(2.0) * __aurng() - 1.0;
1999  __r2 = __x * __x + __y * __y;
2000  }
2001  while (__r2 > 1.0 || __r2 == 0.0);
2002 
2003  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2004  *__f++ = __y * __mult * __param.stddev() + __param.mean();
2005  *__f++ = __x * __mult * __param.stddev() + __param.mean();
2006  }
2007 
2008  if (__f != __t)
2009  {
2010  result_type __x, __y, __r2;
2011  do
2012  {
2013  __x = result_type(2.0) * __aurng() - 1.0;
2014  __y = result_type(2.0) * __aurng() - 1.0;
2015  __r2 = __x * __x + __y * __y;
2016  }
2017  while (__r2 > 1.0 || __r2 == 0.0);
2018 
2019  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2020  _M_saved = __x * __mult;
2021  _M_saved_available = true;
2022  *__f = __y * __mult * __param.stddev() + __param.mean();
2023  }
2024  }
2025 
2026  template<typename _RealType>
2027  bool
2030  {
2031  if (__d1._M_param == __d2._M_param
2032  && __d1._M_saved_available == __d2._M_saved_available)
2033  {
2034  if (__d1._M_saved_available
2035  && __d1._M_saved == __d2._M_saved)
2036  return true;
2037  else if(!__d1._M_saved_available)
2038  return true;
2039  else
2040  return false;
2041  }
2042  else
2043  return false;
2044  }
2045 
2046  template<typename _RealType, typename _CharT, typename _Traits>
2048  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2049  const normal_distribution<_RealType>& __x)
2050  {
2051  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2052  typedef typename __ostream_type::ios_base __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(' ');
2059  __os.fill(__space);
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>
2076  normal_distribution<_RealType>& __x)
2077  {
2078  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2079  typedef typename __istream_type::ios_base __ios_base;
2080 
2081  const typename __ios_base::fmtflags __flags = __is.flags();
2083 
2084  double __mean, __stddev;
2085  __is >> __mean >> __stddev
2086  >> __x._M_saved_available;
2087  if (__x._M_saved_available)
2088  __is >> __x._M_saved;
2089  __x.param(typename normal_distribution<_RealType>::
2090  param_type(__mean, __stddev));
2091 
2092  __is.flags(__flags);
2093  return __is;
2094  }
2095 
2096 
2097  template<typename _RealType>
2098  template<typename _ForwardIterator,
2099  typename _UniformRandomNumberGenerator>
2100  void
2101  lognormal_distribution<_RealType>::
2102  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2103  _UniformRandomNumberGenerator& __urng,
2104  const param_type& __p)
2105  {
2106  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2107  while (__f != __t)
2108  *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
2109  }
2110 
2111  template<typename _RealType, typename _CharT, typename _Traits>
2112  std::basic_ostream<_CharT, _Traits>&
2113  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2114  const lognormal_distribution<_RealType>& __x)
2115  {
2116  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2117  typedef typename __ostream_type::ios_base __ios_base;
2118 
2119  const typename __ios_base::fmtflags __flags = __os.flags();
2120  const _CharT __fill = __os.fill();
2121  const std::streamsize __precision = __os.precision();
2122  const _CharT __space = __os.widen(' ');
2124  __os.fill(__space);
2126 
2127  __os << __x.m() << __space << __x.s()
2128  << __space << __x._M_nd;
2129 
2130  __os.flags(__flags);
2131  __os.fill(__fill);
2132  __os.precision(__precision);
2133  return __os;
2134  }
2135 
2136  template<typename _RealType, typename _CharT, typename _Traits>
2139  lognormal_distribution<_RealType>& __x)
2140  {
2141  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2142  typedef typename __istream_type::ios_base __ios_base;
2143 
2144  const typename __ios_base::fmtflags __flags = __is.flags();
2146 
2147  _RealType __m, __s;
2148  __is >> __m >> __s >> __x._M_nd;
2149  __x.param(typename lognormal_distribution<_RealType>::
2150  param_type(__m, __s));
2151 
2152  __is.flags(__flags);
2153  return __is;
2154  }
2155 
2156  template<typename _RealType>
2157  template<typename _ForwardIterator,
2158  typename _UniformRandomNumberGenerator>
2159  void
2161  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2162  _UniformRandomNumberGenerator& __urng)
2163  {
2164  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2165  while (__f != __t)
2166  *__f++ = 2 * _M_gd(__urng);
2167  }
2168 
2169  template<typename _RealType>
2170  template<typename _ForwardIterator,
2171  typename _UniformRandomNumberGenerator>
2172  void
2173  std::chi_squared_distribution<_RealType>::
2174  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2175  _UniformRandomNumberGenerator& __urng,
2176  const typename
2177  std::gamma_distribution<result_type>::param_type& __p)
2178  {
2179  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2180  while (__f != __t)
2181  *__f++ = 2 * _M_gd(__urng, __p);
2182  }
2183 
2184  template<typename _RealType, typename _CharT, typename _Traits>
2185  std::basic_ostream<_CharT, _Traits>&
2186  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2187  const chi_squared_distribution<_RealType>& __x)
2188  {
2189  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2190  typedef typename __ostream_type::ios_base __ios_base;
2191 
2192  const typename __ios_base::fmtflags __flags = __os.flags();
2193  const _CharT __fill = __os.fill();
2194  const std::streamsize __precision = __os.precision();
2195  const _CharT __space = __os.widen(' ');
2197  __os.fill(__space);
2199 
2200  __os << __x.n() << __space << __x._M_gd;
2201 
2202  __os.flags(__flags);
2203  __os.fill(__fill);
2204  __os.precision(__precision);
2205  return __os;
2206  }
2207 
2208  template<typename _RealType, typename _CharT, typename _Traits>
2211  chi_squared_distribution<_RealType>& __x)
2212  {
2213  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2214  typedef typename __istream_type::ios_base __ios_base;
2215 
2216  const typename __ios_base::fmtflags __flags = __is.flags();
2218 
2219  _RealType __n;
2220  __is >> __n >> __x._M_gd;
2221  __x.param(typename chi_squared_distribution<_RealType>::
2222  param_type(__n));
2223 
2224  __is.flags(__flags);
2225  return __is;
2226  }
2227 
2228 
2229  template<typename _RealType>
2230  template<typename _UniformRandomNumberGenerator>
2231  typename cauchy_distribution<_RealType>::result_type
2233  operator()(_UniformRandomNumberGenerator& __urng,
2234  const param_type& __p)
2235  {
2236  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2237  __aurng(__urng);
2238  _RealType __u;
2239  do
2240  __u = __aurng();
2241  while (__u == 0.5);
2242 
2243  const _RealType __pi = 3.1415926535897932384626433832795029L;
2244  return __p.a() + __p.b() * std::tan(__pi * __u);
2245  }
2246 
2247  template<typename _RealType>
2248  template<typename _ForwardIterator,
2249  typename _UniformRandomNumberGenerator>
2250  void
2251  cauchy_distribution<_RealType>::
2252  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2253  _UniformRandomNumberGenerator& __urng,
2254  const param_type& __p)
2255  {
2256  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2257  const _RealType __pi = 3.1415926535897932384626433832795029L;
2258  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2259  __aurng(__urng);
2260  while (__f != __t)
2261  {
2262  _RealType __u;
2263  do
2264  __u = __aurng();
2265  while (__u == 0.5);
2266 
2267  *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2268  }
2269  }
2270 
2271  template<typename _RealType, typename _CharT, typename _Traits>
2273  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2274  const cauchy_distribution<_RealType>& __x)
2275  {
2276  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2277  typedef typename __ostream_type::ios_base __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(' ');
2284  __os.fill(__space);
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>
2299  {
2300  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2301  typedef typename __istream_type::ios_base __ios_base;
2302 
2303  const typename __ios_base::fmtflags __flags = __is.flags();
2305 
2306  _RealType __a, __b;
2307  __is >> __a >> __b;
2308  __x.param(typename cauchy_distribution<_RealType>::
2309  param_type(__a, __b));
2310 
2311  __is.flags(__flags);
2312  return __is;
2313  }
2314 
2315 
2316  template<typename _RealType>
2317  template<typename _ForwardIterator,
2318  typename _UniformRandomNumberGenerator>
2319  void
2321  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2322  _UniformRandomNumberGenerator& __urng)
2323  {
2324  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2325  while (__f != __t)
2326  *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2327  }
2328 
2329  template<typename _RealType>
2330  template<typename _ForwardIterator,
2331  typename _UniformRandomNumberGenerator>
2332  void
2333  std::fisher_f_distribution<_RealType>::
2334  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2335  _UniformRandomNumberGenerator& __urng,
2336  const param_type& __p)
2337  {
2338  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2339  typedef typename std::gamma_distribution<result_type>::param_type
2340  param_type;
2341  param_type __p1(__p.m() / 2);
2342  param_type __p2(__p.n() / 2);
2343  while (__f != __t)
2344  *__f++ = ((_M_gd_x(__urng, __p1) * n())
2345  / (_M_gd_y(__urng, __p2) * m()));
2346  }
2347 
2348  template<typename _RealType, typename _CharT, typename _Traits>
2349  std::basic_ostream<_CharT, _Traits>&
2350  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2351  const fisher_f_distribution<_RealType>& __x)
2352  {
2353  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2354  typedef typename __ostream_type::ios_base __ios_base;
2355 
2356  const typename __ios_base::fmtflags __flags = __os.flags();
2357  const _CharT __fill = __os.fill();
2358  const std::streamsize __precision = __os.precision();
2359  const _CharT __space = __os.widen(' ');
2361  __os.fill(__space);
2363 
2364  __os << __x.m() << __space << __x.n()
2365  << __space << __x._M_gd_x << __space << __x._M_gd_y;
2366 
2367  __os.flags(__flags);
2368  __os.fill(__fill);
2369  __os.precision(__precision);
2370  return __os;
2371  }
2372 
2373  template<typename _RealType, typename _CharT, typename _Traits>
2376  fisher_f_distribution<_RealType>& __x)
2377  {
2378  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2379  typedef typename __istream_type::ios_base __ios_base;
2380 
2381  const typename __ios_base::fmtflags __flags = __is.flags();
2383 
2384  _RealType __m, __n;
2385  __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
2386  __x.param(typename fisher_f_distribution<_RealType>::
2387  param_type(__m, __n));
2388 
2389  __is.flags(__flags);
2390  return __is;
2391  }
2392 
2393 
2394  template<typename _RealType>
2395  template<typename _ForwardIterator,
2396  typename _UniformRandomNumberGenerator>
2397  void
2399  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2400  _UniformRandomNumberGenerator& __urng)
2401  {
2402  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2403  while (__f != __t)
2404  *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2405  }
2406 
2407  template<typename _RealType>
2408  template<typename _ForwardIterator,
2409  typename _UniformRandomNumberGenerator>
2410  void
2411  std::student_t_distribution<_RealType>::
2412  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2413  _UniformRandomNumberGenerator& __urng,
2414  const param_type& __p)
2415  {
2416  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2417  typename std::gamma_distribution<result_type>::param_type
2418  __p2(__p.n() / 2, 2);
2419  while (__f != __t)
2420  *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2421  }
2422 
2423  template<typename _RealType, typename _CharT, typename _Traits>
2424  std::basic_ostream<_CharT, _Traits>&
2425  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2426  const student_t_distribution<_RealType>& __x)
2427  {
2428  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2429  typedef typename __ostream_type::ios_base __ios_base;
2430 
2431  const typename __ios_base::fmtflags __flags = __os.flags();
2432  const _CharT __fill = __os.fill();
2433  const std::streamsize __precision = __os.precision();
2434  const _CharT __space = __os.widen(' ');
2436  __os.fill(__space);
2438 
2439  __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2440 
2441  __os.flags(__flags);
2442  __os.fill(__fill);
2443  __os.precision(__precision);
2444  return __os;
2445  }
2446 
2447  template<typename _RealType, typename _CharT, typename _Traits>
2450  student_t_distribution<_RealType>& __x)
2451  {
2452  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2453  typedef typename __istream_type::ios_base __ios_base;
2454 
2455  const typename __ios_base::fmtflags __flags = __is.flags();
2457 
2458  _RealType __n;
2459  __is >> __n >> __x._M_nd >> __x._M_gd;
2460  __x.param(typename student_t_distribution<_RealType>::param_type(__n));
2461 
2462  __is.flags(__flags);
2463  return __is;
2464  }
2465 
2466 
2467  template<typename _RealType>
2468  void
2469  gamma_distribution<_RealType>::param_type::
2470  _M_initialize()
2471  {
2472  _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2473 
2474  const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2475  _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2476  }
2477 
2478  /**
2479  * Marsaglia, G. and Tsang, W. W.
2480  * "A Simple Method for Generating Gamma Variables"
2481  * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2482  */
2483  template<typename _RealType>
2484  template<typename _UniformRandomNumberGenerator>
2485  typename gamma_distribution<_RealType>::result_type
2487  operator()(_UniformRandomNumberGenerator& __urng,
2488  const param_type& __param)
2489  {
2490  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2491  __aurng(__urng);
2492 
2493  result_type __u, __v, __n;
2494  const result_type __a1 = (__param._M_malpha
2495  - _RealType(1.0) / _RealType(3.0));
2496 
2497  do
2498  {
2499  do
2500  {
2501  __n = _M_nd(__urng);
2502  __v = result_type(1.0) + __param._M_a2 * __n;
2503  }
2504  while (__v <= 0.0);
2505 
2506  __v = __v * __v * __v;
2507  __u = __aurng();
2508  }
2509  while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2510  && (std::log(__u) > (0.5 * __n * __n + __a1
2511  * (1.0 - __v + std::log(__v)))));
2512 
2513  if (__param.alpha() == __param._M_malpha)
2514  return __a1 * __v * __param.beta();
2515  else
2516  {
2517  do
2518  __u = __aurng();
2519  while (__u == 0.0);
2520 
2521  return (std::pow(__u, result_type(1.0) / __param.alpha())
2522  * __a1 * __v * __param.beta());
2523  }
2524  }
2525 
2526  template<typename _RealType>
2527  template<typename _ForwardIterator,
2528  typename _UniformRandomNumberGenerator>
2529  void
2531  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2532  _UniformRandomNumberGenerator& __urng,
2533  const param_type& __param)
2534  {
2535  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2536  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2537  __aurng(__urng);
2538 
2539  result_type __u, __v, __n;
2540  const result_type __a1 = (__param._M_malpha
2541  - _RealType(1.0) / _RealType(3.0));
2542 
2543  if (__param.alpha() == __param._M_malpha)
2544  while (__f != __t)
2545  {
2546  do
2547  {
2548  do
2549  {
2550  __n = _M_nd(__urng);
2551  __v = result_type(1.0) + __param._M_a2 * __n;
2552  }
2553  while (__v <= 0.0);
2554 
2555  __v = __v * __v * __v;
2556  __u = __aurng();
2557  }
2558  while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2559  && (std::log(__u) > (0.5 * __n * __n + __a1
2560  * (1.0 - __v + std::log(__v)))));
2561 
2562  *__f++ = __a1 * __v * __param.beta();
2563  }
2564  else
2565  while (__f != __t)
2566  {
2567  do
2568  {
2569  do
2570  {
2571  __n = _M_nd(__urng);
2572  __v = result_type(1.0) + __param._M_a2 * __n;
2573  }
2574  while (__v <= 0.0);
2575 
2576  __v = __v * __v * __v;
2577  __u = __aurng();
2578  }
2579  while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2580  && (std::log(__u) > (0.5 * __n * __n + __a1
2581  * (1.0 - __v + std::log(__v)))));
2582 
2583  do
2584  __u = __aurng();
2585  while (__u == 0.0);
2586 
2587  *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2588  * __a1 * __v * __param.beta());
2589  }
2590  }
2591 
2592  template<typename _RealType, typename _CharT, typename _Traits>
2594  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2595  const gamma_distribution<_RealType>& __x)
2596  {
2597  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2598  typedef typename __ostream_type::ios_base __ios_base;
2599 
2600  const typename __ios_base::fmtflags __flags = __os.flags();
2601  const _CharT __fill = __os.fill();
2602  const std::streamsize __precision = __os.precision();
2603  const _CharT __space = __os.widen(' ');
2605  __os.fill(__space);
2607 
2608  __os << __x.alpha() << __space << __x.beta()
2609  << __space << __x._M_nd;
2610 
2611  __os.flags(__flags);
2612  __os.fill(__fill);
2613  __os.precision(__precision);
2614  return __os;
2615  }
2616 
2617  template<typename _RealType, typename _CharT, typename _Traits>
2620  gamma_distribution<_RealType>& __x)
2621  {
2622  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2623  typedef typename __istream_type::ios_base __ios_base;
2624 
2625  const typename __ios_base::fmtflags __flags = __is.flags();
2627 
2628  _RealType __alpha_val, __beta_val;
2629  __is >> __alpha_val >> __beta_val >> __x._M_nd;
2630  __x.param(typename gamma_distribution<_RealType>::
2631  param_type(__alpha_val, __beta_val));
2632 
2633  __is.flags(__flags);
2634  return __is;
2635  }
2636 
2637 
2638  template<typename _RealType>
2639  template<typename _UniformRandomNumberGenerator>
2640  typename weibull_distribution<_RealType>::result_type
2642  operator()(_UniformRandomNumberGenerator& __urng,
2643  const param_type& __p)
2644  {
2645  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2646  __aurng(__urng);
2647  return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2648  result_type(1) / __p.a());
2649  }
2650 
2651  template<typename _RealType>
2652  template<typename _ForwardIterator,
2653  typename _UniformRandomNumberGenerator>
2654  void
2655  weibull_distribution<_RealType>::
2656  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2657  _UniformRandomNumberGenerator& __urng,
2658  const param_type& __p)
2659  {
2660  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2661  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2662  __aurng(__urng);
2663  auto __inv_a = result_type(1) / __p.a();
2664 
2665  while (__f != __t)
2666  *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2667  __inv_a);
2668  }
2669 
2670  template<typename _RealType, typename _CharT, typename _Traits>
2671  std::basic_ostream<_CharT, _Traits>&
2672  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2673  const weibull_distribution<_RealType>& __x)
2674  {
2675  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2676  typedef typename __ostream_type::ios_base __ios_base;
2677 
2678  const typename __ios_base::fmtflags __flags = __os.flags();
2679  const _CharT __fill = __os.fill();
2680  const std::streamsize __precision = __os.precision();
2681  const _CharT __space = __os.widen(' ');
2683  __os.fill(__space);
2685 
2686  __os << __x.a() << __space << __x.b();
2687 
2688  __os.flags(__flags);
2689  __os.fill(__fill);
2690  __os.precision(__precision);
2691  return __os;
2692  }
2693 
2694  template<typename _RealType, typename _CharT, typename _Traits>
2698  {
2699  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2700  typedef typename __istream_type::ios_base __ios_base;
2701 
2702  const typename __ios_base::fmtflags __flags = __is.flags();
2704 
2705  _RealType __a, __b;
2706  __is >> __a >> __b;
2707  __x.param(typename weibull_distribution<_RealType>::
2708  param_type(__a, __b));
2709 
2710  __is.flags(__flags);
2711  return __is;
2712  }
2713 
2714 
2715  template<typename _RealType>
2716  template<typename _UniformRandomNumberGenerator>
2717  typename extreme_value_distribution<_RealType>::result_type
2719  operator()(_UniformRandomNumberGenerator& __urng,
2720  const param_type& __p)
2721  {
2722  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2723  __aurng(__urng);
2724  return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2725  - __aurng()));
2726  }
2727 
2728  template<typename _RealType>
2729  template<typename _ForwardIterator,
2730  typename _UniformRandomNumberGenerator>
2731  void
2732  extreme_value_distribution<_RealType>::
2733  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2734  _UniformRandomNumberGenerator& __urng,
2735  const param_type& __p)
2736  {
2737  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2738  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2739  __aurng(__urng);
2740 
2741  while (__f != __t)
2742  *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2743  - __aurng()));
2744  }
2745 
2746  template<typename _RealType, typename _CharT, typename _Traits>
2747  std::basic_ostream<_CharT, _Traits>&
2748  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2749  const extreme_value_distribution<_RealType>& __x)
2750  {
2751  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2752  typedef typename __ostream_type::ios_base __ios_base;
2753 
2754  const typename __ios_base::fmtflags __flags = __os.flags();
2755  const _CharT __fill = __os.fill();
2756  const std::streamsize __precision = __os.precision();
2757  const _CharT __space = __os.widen(' ');
2759  __os.fill(__space);
2761 
2762  __os << __x.a() << __space << __x.b();
2763 
2764  __os.flags(__flags);
2765  __os.fill(__fill);
2766  __os.precision(__precision);
2767  return __os;
2768  }
2769 
2770  template<typename _RealType, typename _CharT, typename _Traits>
2774  {
2775  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2776  typedef typename __istream_type::ios_base __ios_base;
2777 
2778  const typename __ios_base::fmtflags __flags = __is.flags();
2780 
2781  _RealType __a, __b;
2782  __is >> __a >> __b;
2784  param_type(__a, __b));
2785 
2786  __is.flags(__flags);
2787  return __is;
2788  }
2789 
2790 
2791  template<typename _IntType>
2792  void
2793  discrete_distribution<_IntType>::param_type::
2794  _M_initialize()
2795  {
2796  if (_M_prob.size() < 2)
2797  {
2798  _M_prob.clear();
2799  return;
2800  }
2801 
2802  const double __sum = std::accumulate(_M_prob.begin(),
2803  _M_prob.end(), 0.0);
2804  // Now normalize the probabilites.
2805  __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2806  __sum);
2807  // Accumulate partial sums.
2808  _M_cp.reserve(_M_prob.size());
2809  std::partial_sum(_M_prob.begin(), _M_prob.end(),
2810  std::back_inserter(_M_cp));
2811  // Make sure the last cumulative probability is one.
2812  _M_cp[_M_cp.size() - 1] = 1.0;
2813  }
2814 
2815  template<typename _IntType>
2816  template<typename _Func>
2817  discrete_distribution<_IntType>::param_type::
2818  param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2819  : _M_prob(), _M_cp()
2820  {
2821  const size_t __n = __nw == 0 ? 1 : __nw;
2822  const double __delta = (__xmax - __xmin) / __n;
2823 
2824  _M_prob.reserve(__n);
2825  for (size_t __k = 0; __k < __nw; ++__k)
2826  _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2827 
2828  _M_initialize();
2829  }
2830 
2831  template<typename _IntType>
2832  template<typename _UniformRandomNumberGenerator>
2833  typename discrete_distribution<_IntType>::result_type
2834  discrete_distribution<_IntType>::
2835  operator()(_UniformRandomNumberGenerator& __urng,
2836  const param_type& __param)
2837  {
2838  if (__param._M_cp.empty())
2839  return result_type(0);
2840 
2841  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2842  __aurng(__urng);
2843 
2844  const double __p = __aurng();
2845  auto __pos = std::lower_bound(__param._M_cp.begin(),
2846  __param._M_cp.end(), __p);
2847 
2848  return __pos - __param._M_cp.begin();
2849  }
2850 
2851  template<typename _IntType>
2852  template<typename _ForwardIterator,
2853  typename _UniformRandomNumberGenerator>
2854  void
2855  discrete_distribution<_IntType>::
2856  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2857  _UniformRandomNumberGenerator& __urng,
2858  const param_type& __param)
2859  {
2860  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2861 
2862  if (__param._M_cp.empty())
2863  {
2864  while (__f != __t)
2865  *__f++ = result_type(0);
2866  return;
2867  }
2868 
2869  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2870  __aurng(__urng);
2871 
2872  while (__f != __t)
2873  {
2874  const double __p = __aurng();
2875  auto __pos = std::lower_bound(__param._M_cp.begin(),
2876  __param._M_cp.end(), __p);
2877 
2878  *__f++ = __pos - __param._M_cp.begin();
2879  }
2880  }
2881 
2882  template<typename _IntType, typename _CharT, typename _Traits>
2884  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2885  const discrete_distribution<_IntType>& __x)
2886  {
2887  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2888  typedef typename __ostream_type::ios_base __ios_base;
2889 
2890  const typename __ios_base::fmtflags __flags = __os.flags();
2891  const _CharT __fill = __os.fill();
2892  const std::streamsize __precision = __os.precision();
2893  const _CharT __space = __os.widen(' ');
2895  __os.fill(__space);
2897 
2898  std::vector<double> __prob = __x.probabilities();
2899  __os << __prob.size();
2900  for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2901  __os << __space << *__dit;
2902 
2903  __os.flags(__flags);
2904  __os.fill(__fill);
2905  __os.precision(__precision);
2906  return __os;
2907  }
2908 
2909  template<typename _IntType, typename _CharT, typename _Traits>
2912  discrete_distribution<_IntType>& __x)
2913  {
2914  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2915  typedef typename __istream_type::ios_base __ios_base;
2916 
2917  const typename __ios_base::fmtflags __flags = __is.flags();
2919 
2920  size_t __n;
2921  __is >> __n;
2922 
2923  std::vector<double> __prob_vec;
2924  __prob_vec.reserve(__n);
2925  for (; __n != 0; --__n)
2926  {
2927  double __prob;
2928  __is >> __prob;
2929  __prob_vec.push_back(__prob);
2930  }
2931 
2932  __x.param(typename discrete_distribution<_IntType>::
2933  param_type(__prob_vec.begin(), __prob_vec.end()));
2934 
2935  __is.flags(__flags);
2936  return __is;
2937  }
2938 
2939 
2940  template<typename _RealType>
2941  void
2942  piecewise_constant_distribution<_RealType>::param_type::
2943  _M_initialize()
2944  {
2945  if (_M_int.size() < 2
2946  || (_M_int.size() == 2
2947  && _M_int[0] == _RealType(0)
2948  && _M_int[1] == _RealType(1)))
2949  {
2950  _M_int.clear();
2951  _M_den.clear();
2952  return;
2953  }
2954 
2955  const double __sum = std::accumulate(_M_den.begin(),
2956  _M_den.end(), 0.0);
2957 
2958  __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2959  __sum);
2960 
2961  _M_cp.reserve(_M_den.size());
2962  std::partial_sum(_M_den.begin(), _M_den.end(),
2963  std::back_inserter(_M_cp));
2964 
2965  // Make sure the last cumulative probability is one.
2966  _M_cp[_M_cp.size() - 1] = 1.0;
2967 
2968  for (size_t __k = 0; __k < _M_den.size(); ++__k)
2969  _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2970  }
2971 
2972  template<typename _RealType>
2973  template<typename _InputIteratorB, typename _InputIteratorW>
2974  piecewise_constant_distribution<_RealType>::param_type::
2975  param_type(_InputIteratorB __bbegin,
2976  _InputIteratorB __bend,
2977  _InputIteratorW __wbegin)
2978  : _M_int(), _M_den(), _M_cp()
2979  {
2980  if (__bbegin != __bend)
2981  {
2982  for (;;)
2983  {
2984  _M_int.push_back(*__bbegin);
2985  ++__bbegin;
2986  if (__bbegin == __bend)
2987  break;
2988 
2989  _M_den.push_back(*__wbegin);
2990  ++__wbegin;
2991  }
2992  }
2993 
2994  _M_initialize();
2995  }
2996 
2997  template<typename _RealType>
2998  template<typename _Func>
2999  piecewise_constant_distribution<_RealType>::param_type::
3000  param_type(initializer_list<_RealType> __bl, _Func __fw)
3001  : _M_int(), _M_den(), _M_cp()
3002  {
3003  _M_int.reserve(__bl.size());
3004  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3005  _M_int.push_back(*__biter);
3006 
3007  _M_den.reserve(_M_int.size() - 1);
3008  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3009  _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
3010 
3011  _M_initialize();
3012  }
3013 
3014  template<typename _RealType>
3015  template<typename _Func>
3016  piecewise_constant_distribution<_RealType>::param_type::
3017  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3018  : _M_int(), _M_den(), _M_cp()
3019  {
3020  const size_t __n = __nw == 0 ? 1 : __nw;
3021  const _RealType __delta = (__xmax - __xmin) / __n;
3022 
3023  _M_int.reserve(__n + 1);
3024  for (size_t __k = 0; __k <= __nw; ++__k)
3025  _M_int.push_back(__xmin + __k * __delta);
3026 
3027  _M_den.reserve(__n);
3028  for (size_t __k = 0; __k < __nw; ++__k)
3029  _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
3030 
3031  _M_initialize();
3032  }
3033 
3034  template<typename _RealType>
3035  template<typename _UniformRandomNumberGenerator>
3036  typename piecewise_constant_distribution<_RealType>::result_type
3037  piecewise_constant_distribution<_RealType>::
3038  operator()(_UniformRandomNumberGenerator& __urng,
3039  const param_type& __param)
3040  {
3041  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3042  __aurng(__urng);
3043 
3044  const double __p = __aurng();
3045  if (__param._M_cp.empty())
3046  return __p;
3047 
3048  auto __pos = std::lower_bound(__param._M_cp.begin(),
3049  __param._M_cp.end(), __p);
3050  const size_t __i = __pos - __param._M_cp.begin();
3051 
3052  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3053 
3054  return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
3055  }
3056 
3057  template<typename _RealType>
3058  template<typename _ForwardIterator,
3059  typename _UniformRandomNumberGenerator>
3060  void
3061  piecewise_constant_distribution<_RealType>::
3062  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3063  _UniformRandomNumberGenerator& __urng,
3064  const param_type& __param)
3065  {
3066  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3067  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3068  __aurng(__urng);
3069 
3070  if (__param._M_cp.empty())
3071  {
3072  while (__f != __t)
3073  *__f++ = __aurng();
3074  return;
3075  }
3076 
3077  while (__f != __t)
3078  {
3079  const double __p = __aurng();
3080 
3081  auto __pos = std::lower_bound(__param._M_cp.begin(),
3082  __param._M_cp.end(), __p);
3083  const size_t __i = __pos - __param._M_cp.begin();
3084 
3085  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3086 
3087  *__f++ = (__param._M_int[__i]
3088  + (__p - __pref) / __param._M_den[__i]);
3089  }
3090  }
3091 
3092  template<typename _RealType, typename _CharT, typename _Traits>
3094  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3095  const piecewise_constant_distribution<_RealType>& __x)
3096  {
3097  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
3098  typedef typename __ostream_type::ios_base __ios_base;
3099 
3100  const typename __ios_base::fmtflags __flags = __os.flags();
3101  const _CharT __fill = __os.fill();
3102  const std::streamsize __precision = __os.precision();
3103  const _CharT __space = __os.widen(' ');
3105  __os.fill(__space);
3107 
3108  std::vector<_RealType> __int = __x.intervals();
3109  __os << __int.size() - 1;
3110 
3111  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3112  __os << __space << *__xit;
3113 
3114  std::vector<double> __den = __x.densities();
3115  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3116  __os << __space << *__dit;
3117 
3118  __os.flags(__flags);
3119  __os.fill(__fill);
3120  __os.precision(__precision);
3121  return __os;
3122  }
3123 
3124  template<typename _RealType, typename _CharT, typename _Traits>
3127  piecewise_constant_distribution<_RealType>& __x)
3128  {
3129  typedef std::basic_istream<_CharT, _Traits> __istream_type;
3130  typedef typename __istream_type::ios_base __ios_base;
3131 
3132  const typename __ios_base::fmtflags __flags = __is.flags();
3134 
3135  size_t __n;
3136  __is >> __n;
3137 
3138  std::vector<_RealType> __int_vec;
3139  __int_vec.reserve(__n + 1);
3140  for (size_t __i = 0; __i <= __n; ++__i)
3141  {
3142  _RealType __int;
3143  __is >> __int;
3144  __int_vec.push_back(__int);
3145  }
3146 
3147  std::vector<double> __den_vec;
3148  __den_vec.reserve(__n);
3149  for (size_t __i = 0; __i < __n; ++__i)
3150  {
3151  double __den;
3152  __is >> __den;
3153  __den_vec.push_back(__den);
3154  }
3155 
3156  __x.param(typename piecewise_constant_distribution<_RealType>::
3157  param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3158 
3159  __is.flags(__flags);
3160  return __is;
3161  }
3162 
3163 
3164  template<typename _RealType>
3165  void
3166  piecewise_linear_distribution<_RealType>::param_type::
3167  _M_initialize()
3168  {
3169  if (_M_int.size() < 2
3170  || (_M_int.size() == 2
3171  && _M_int[0] == _RealType(0)
3172  && _M_int[1] == _RealType(1)
3173  && _M_den[0] == _M_den[1]))
3174  {
3175  _M_int.clear();
3176  _M_den.clear();
3177  return;
3178  }
3179 
3180  double __sum = 0.0;
3181  _M_cp.reserve(_M_int.size() - 1);
3182  _M_m.reserve(_M_int.size() - 1);
3183  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3184  {
3185  const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3186  __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3187  _M_cp.push_back(__sum);
3188  _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3189  }
3190 
3191  // Now normalize the densities...
3192  __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3193  __sum);
3194  // ... and partial sums...
3195  __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3196  // ... and slopes.
3197  __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3198 
3199  // Make sure the last cumulative probablility is one.
3200  _M_cp[_M_cp.size() - 1] = 1.0;
3201  }
3202 
3203  template<typename _RealType>
3204  template<typename _InputIteratorB, typename _InputIteratorW>
3205  piecewise_linear_distribution<_RealType>::param_type::
3206  param_type(_InputIteratorB __bbegin,
3207  _InputIteratorB __bend,
3208  _InputIteratorW __wbegin)
3209  : _M_int(), _M_den(), _M_cp(), _M_m()
3210  {
3211  for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
3212  {
3213  _M_int.push_back(*__bbegin);
3214  _M_den.push_back(*__wbegin);
3215  }
3216 
3217  _M_initialize();
3218  }
3219 
3220  template<typename _RealType>
3221  template<typename _Func>
3222  piecewise_linear_distribution<_RealType>::param_type::
3223  param_type(initializer_list<_RealType> __bl, _Func __fw)
3224  : _M_int(), _M_den(), _M_cp(), _M_m()
3225  {
3226  _M_int.reserve(__bl.size());
3227  _M_den.reserve(__bl.size());
3228  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3229  {
3230  _M_int.push_back(*__biter);
3231  _M_den.push_back(__fw(*__biter));
3232  }
3233 
3234  _M_initialize();
3235  }
3236 
3237  template<typename _RealType>
3238  template<typename _Func>
3239  piecewise_linear_distribution<_RealType>::param_type::
3240  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3241  : _M_int(), _M_den(), _M_cp(), _M_m()
3242  {
3243  const size_t __n = __nw == 0 ? 1 : __nw;
3244  const _RealType __delta = (__xmax - __xmin) / __n;
3245 
3246  _M_int.reserve(__n + 1);
3247  _M_den.reserve(__n + 1);
3248  for (size_t __k = 0; __k <= __nw; ++__k)
3249  {
3250  _M_int.push_back(__xmin + __k * __delta);
3251  _M_den.push_back(__fw(_M_int[__k] + __delta));
3252  }
3253 
3254  _M_initialize();
3255  }
3256 
3257  template<typename _RealType>
3258  template<typename _UniformRandomNumberGenerator>
3259  typename piecewise_linear_distribution<_RealType>::result_type
3260  piecewise_linear_distribution<_RealType>::
3261  operator()(_UniformRandomNumberGenerator& __urng,
3262  const param_type& __param)
3263  {
3264  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3265  __aurng(__urng);
3266 
3267  const double __p = __aurng();
3268  if (__param._M_cp.empty())
3269  return __p;
3270 
3271  auto __pos = std::lower_bound(__param._M_cp.begin(),
3272  __param._M_cp.end(), __p);
3273  const size_t __i = __pos - __param._M_cp.begin();
3274 
3275  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3276 
3277  const double __a = 0.5 * __param._M_m[__i];
3278  const double __b = __param._M_den[__i];
3279  const double __cm = __p - __pref;
3280 
3281  _RealType __x = __param._M_int[__i];
3282  if (__a == 0)
3283  __x += __cm / __b;
3284  else
3285  {
3286  const double __d = __b * __b + 4.0 * __a * __cm;
3287  __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3288  }
3289 
3290  return __x;
3291  }
3292 
3293  template<typename _RealType>
3294  template<typename _ForwardIterator,
3295  typename _UniformRandomNumberGenerator>
3296  void
3297  piecewise_linear_distribution<_RealType>::
3298  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3299  _UniformRandomNumberGenerator& __urng,
3300  const param_type& __param)
3301  {
3302  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3303  // We could duplicate everything from operator()...
3304  while (__f != __t)
3305  *__f++ = this->operator()(__urng, __param);
3306  }
3307 
3308  template<typename _RealType, typename _CharT, typename _Traits>
3309  std::basic_ostream<_CharT, _Traits>&
3310  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3311  const piecewise_linear_distribution<_RealType>& __x)
3312  {
3313  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
3314  typedef typename __ostream_type::ios_base __ios_base;
3315 
3316  const typename __ios_base::fmtflags __flags = __os.flags();
3317  const _CharT __fill = __os.fill();
3318  const std::streamsize __precision = __os.precision();
3319  const _CharT __space = __os.widen(' ');
3321  __os.fill(__space);
3323 
3324  std::vector<_RealType> __int = __x.intervals();
3325  __os << __int.size() - 1;
3326 
3327  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3328  __os << __space << *__xit;
3329 
3330  std::vector<double> __den = __x.densities();
3331  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3332  __os << __space << *__dit;
3333 
3334  __os.flags(__flags);
3335  __os.fill(__fill);
3336  __os.precision(__precision);
3337  return __os;
3338  }
3339 
3340  template<typename _RealType, typename _CharT, typename _Traits>
3343  piecewise_linear_distribution<_RealType>& __x)
3344  {
3345  typedef std::basic_istream<_CharT, _Traits> __istream_type;
3346  typedef typename __istream_type::ios_base __ios_base;
3347 
3348  const typename __ios_base::fmtflags __flags = __is.flags();
3350 
3351  size_t __n;
3352  __is >> __n;
3353 
3354  std::vector<_RealType> __int_vec;
3355  __int_vec.reserve(__n + 1);
3356  for (size_t __i = 0; __i <= __n; ++__i)
3357  {
3358  _RealType __int;
3359  __is >> __int;
3360  __int_vec.push_back(__int);
3361  }
3362 
3363  std::vector<double> __den_vec;
3364  __den_vec.reserve(__n + 1);
3365  for (size_t __i = 0; __i <= __n; ++__i)
3366  {
3367  double __den;
3368  __is >> __den;
3369  __den_vec.push_back(__den);
3370  }
3371 
3372  __x.param(typename piecewise_linear_distribution<_RealType>::
3373  param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3374 
3375  __is.flags(__flags);
3376  return __is;
3377  }
3378 
3379 
3380  template<typename _IntType>
3381  seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3382  {
3383  for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3384  _M_v.push_back(__detail::__mod<result_type,
3385  __detail::_Shift<result_type, 32>::__value>(*__iter));
3386  }
3387 
3388  template<typename _InputIterator>
3389  seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3390  {
3391  for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3392  _M_v.push_back(__detail::__mod<result_type,
3393  __detail::_Shift<result_type, 32>::__value>(*__iter));
3394  }
3395 
3396  template<typename _RandomAccessIterator>
3397  void
3398  seed_seq::generate(_RandomAccessIterator __begin,
3399  _RandomAccessIterator __end)
3400  {
3401  typedef typename iterator_traits<_RandomAccessIterator>::value_type
3402  _Type;
3403 
3404  if (__begin == __end)
3405  return;
3406 
3407  std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3408 
3409  const size_t __n = __end - __begin;
3410  const size_t __s = _M_v.size();
3411  const size_t __t = (__n >= 623) ? 11
3412  : (__n >= 68) ? 7
3413  : (__n >= 39) ? 5
3414  : (__n >= 7) ? 3
3415  : (__n - 1) / 2;
3416  const size_t __p = (__n - __t) / 2;
3417  const size_t __q = __p + __t;
3418  const size_t __m = std::max(size_t(__s + 1), __n);
3419 
3420  for (size_t __k = 0; __k < __m; ++__k)
3421  {
3422  _Type __arg = (__begin[__k % __n]
3423  ^ __begin[(__k + __p) % __n]
3424  ^ __begin[(__k - 1) % __n]);
3425  _Type __r1 = __arg ^ (__arg >> 27);
3426  __r1 = __detail::__mod<_Type,
3427  __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
3428  _Type __r2 = __r1;
3429  if (__k == 0)
3430  __r2 += __s;
3431  else if (__k <= __s)
3432  __r2 += __k % __n + _M_v[__k - 1];
3433  else
3434  __r2 += __k % __n;
3435  __r2 = __detail::__mod<_Type,
3436  __detail::_Shift<_Type, 32>::__value>(__r2);
3437  __begin[(__k + __p) % __n] += __r1;
3438  __begin[(__k + __q) % __n] += __r2;
3439  __begin[__k % __n] = __r2;
3440  }
3441 
3442  for (size_t __k = __m; __k < __m + __n; ++__k)
3443  {
3444  _Type __arg = (__begin[__k % __n]
3445  + __begin[(__k + __p) % __n]
3446  + __begin[(__k - 1) % __n]);
3447  _Type __r3 = __arg ^ (__arg >> 27);
3448  __r3 = __detail::__mod<_Type,
3449  __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
3450  _Type __r4 = __r3 - __k % __n;
3451  __r4 = __detail::__mod<_Type,
3452  __detail::_Shift<_Type, 32>::__value>(__r4);
3453  __begin[(__k + __p) % __n] ^= __r3;
3454  __begin[(__k + __q) % __n] ^= __r4;
3455  __begin[__k % __n] = __r4;
3456  }
3457  }
3458 
3459  template<typename _RealType, size_t __bits,
3460  typename _UniformRandomNumberGenerator>
3461  _RealType
3462  generate_canonical(_UniformRandomNumberGenerator& __urng)
3463  {
3464  const size_t __b
3465  = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3466  __bits);
3467  const long double __r = static_cast<long double>(__urng.max())
3468  - static_cast<long double>(__urng.min()) + 1.0L;
3469  const size_t __log2r = std::log(__r) / std::log(2.0L);
3470  size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
3471  _RealType __sum = _RealType(0);
3472  _RealType __tmp = _RealType(1);
3473  for (; __k != 0; --__k)
3474  {
3475  __sum += _RealType(__urng() - __urng.min()) * __tmp;
3476  __tmp *= __r;
3477  }
3478  return __sum / __tmp;
3479  }
3480 
3481 _GLIBCXX_END_NAMESPACE_VERSION
3482 } // namespace
3483 
3484 #endif