libstdc++
bits/random.tcc
Go to the documentation of this file.
1 // random number generation (out of line) -*- C++ -*-
2 
3 // Copyright (C) 2009-2017 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
100 
101  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
102  constexpr _UIntType
104 
105  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
106  constexpr _UIntType
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>
134  typename std::enable_if<std::is_class<_Sseq>::value>::type
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>
180  {
181  typedef std::basic_istream<_CharT, _Traits> __istream_type;
182  typedef typename __istream_type::ios_base __ios_base;
183 
184  const typename __ios_base::fmtflags __flags = __is.flags();
185  __is.flags(__ios_base::dec);
186 
187  __is >> __lcr._M_x;
188 
189  __is.flags(__flags);
190  return __is;
191  }
192 
193 
194  template<typename _UIntType,
195  size_t __w, size_t __n, size_t __m, size_t __r,
196  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
197  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
198  _UIntType __f>
199  constexpr size_t
200  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
201  __s, __b, __t, __c, __l, __f>::word_size;
202 
203  template<typename _UIntType,
204  size_t __w, size_t __n, size_t __m, size_t __r,
205  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
206  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
207  _UIntType __f>
208  constexpr size_t
209  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
210  __s, __b, __t, __c, __l, __f>::state_size;
211 
212  template<typename _UIntType,
213  size_t __w, size_t __n, size_t __m, size_t __r,
214  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
215  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
216  _UIntType __f>
217  constexpr size_t
218  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
219  __s, __b, __t, __c, __l, __f>::shift_size;
220 
221  template<typename _UIntType,
222  size_t __w, size_t __n, size_t __m, size_t __r,
223  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
224  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
225  _UIntType __f>
226  constexpr size_t
227  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
228  __s, __b, __t, __c, __l, __f>::mask_bits;
229 
230  template<typename _UIntType,
231  size_t __w, size_t __n, size_t __m, size_t __r,
232  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
233  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
234  _UIntType __f>
235  constexpr _UIntType
236  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
237  __s, __b, __t, __c, __l, __f>::xor_mask;
238 
239  template<typename _UIntType,
240  size_t __w, size_t __n, size_t __m, size_t __r,
241  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
242  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
243  _UIntType __f>
244  constexpr size_t
245  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
246  __s, __b, __t, __c, __l, __f>::tempering_u;
247 
248  template<typename _UIntType,
249  size_t __w, size_t __n, size_t __m, size_t __r,
250  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
251  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
252  _UIntType __f>
253  constexpr _UIntType
254  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
255  __s, __b, __t, __c, __l, __f>::tempering_d;
256 
257  template<typename _UIntType,
258  size_t __w, size_t __n, size_t __m, size_t __r,
259  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
260  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
261  _UIntType __f>
262  constexpr size_t
263  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
264  __s, __b, __t, __c, __l, __f>::tempering_s;
265 
266  template<typename _UIntType,
267  size_t __w, size_t __n, size_t __m, size_t __r,
268  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
269  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
270  _UIntType __f>
271  constexpr _UIntType
272  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
273  __s, __b, __t, __c, __l, __f>::tempering_b;
274 
275  template<typename _UIntType,
276  size_t __w, size_t __n, size_t __m, size_t __r,
277  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
278  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
279  _UIntType __f>
280  constexpr size_t
281  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
282  __s, __b, __t, __c, __l, __f>::tempering_t;
283 
284  template<typename _UIntType,
285  size_t __w, size_t __n, size_t __m, size_t __r,
286  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
287  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
288  _UIntType __f>
289  constexpr _UIntType
290  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
291  __s, __b, __t, __c, __l, __f>::tempering_c;
292 
293  template<typename _UIntType,
294  size_t __w, size_t __n, size_t __m, size_t __r,
295  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
296  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
297  _UIntType __f>
298  constexpr size_t
299  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
300  __s, __b, __t, __c, __l, __f>::tempering_l;
301 
302  template<typename _UIntType,
303  size_t __w, size_t __n, size_t __m, size_t __r,
304  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
305  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
306  _UIntType __f>
307  constexpr _UIntType
308  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
309  __s, __b, __t, __c, __l, __f>::
310  initialization_multiplier;
311 
312  template<typename _UIntType,
313  size_t __w, size_t __n, size_t __m, size_t __r,
314  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
315  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
316  _UIntType __f>
317  constexpr _UIntType
318  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
319  __s, __b, __t, __c, __l, __f>::default_seed;
320 
321  template<typename _UIntType,
322  size_t __w, size_t __n, size_t __m, size_t __r,
323  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
324  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
325  _UIntType __f>
326  void
327  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
328  __s, __b, __t, __c, __l, __f>::
329  seed(result_type __sd)
330  {
331  _M_x[0] = __detail::__mod<_UIntType,
332  __detail::_Shift<_UIntType, __w>::__value>(__sd);
333 
334  for (size_t __i = 1; __i < state_size; ++__i)
335  {
336  _UIntType __x = _M_x[__i - 1];
337  __x ^= __x >> (__w - 2);
338  __x *= __f;
339  __x += __detail::__mod<_UIntType, __n>(__i);
340  _M_x[__i] = __detail::__mod<_UIntType,
341  __detail::_Shift<_UIntType, __w>::__value>(__x);
342  }
343  _M_p = state_size;
344  }
345 
346  template<typename _UIntType,
347  size_t __w, size_t __n, size_t __m, size_t __r,
348  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
349  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
350  _UIntType __f>
351  template<typename _Sseq>
352  typename std::enable_if<std::is_class<_Sseq>::value>::type
353  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
354  __s, __b, __t, __c, __l, __f>::
355  seed(_Sseq& __q)
356  {
357  const _UIntType __upper_mask = (~_UIntType()) << __r;
358  const size_t __k = (__w + 31) / 32;
359  uint_least32_t __arr[__n * __k];
360  __q.generate(__arr + 0, __arr + __n * __k);
361 
362  bool __zero = true;
363  for (size_t __i = 0; __i < state_size; ++__i)
364  {
365  _UIntType __factor = 1u;
366  _UIntType __sum = 0u;
367  for (size_t __j = 0; __j < __k; ++__j)
368  {
369  __sum += __arr[__k * __i + __j] * __factor;
370  __factor *= __detail::_Shift<_UIntType, 32>::__value;
371  }
372  _M_x[__i] = __detail::__mod<_UIntType,
373  __detail::_Shift<_UIntType, __w>::__value>(__sum);
374 
375  if (__zero)
376  {
377  if (__i == 0)
378  {
379  if ((_M_x[0] & __upper_mask) != 0u)
380  __zero = false;
381  }
382  else if (_M_x[__i] != 0u)
383  __zero = false;
384  }
385  }
386  if (__zero)
387  _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
388  _M_p = state_size;
389  }
390 
391  template<typename _UIntType, size_t __w,
392  size_t __n, size_t __m, size_t __r,
393  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
394  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
395  _UIntType __f>
396  void
397  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
398  __s, __b, __t, __c, __l, __f>::
399  _M_gen_rand(void)
400  {
401  const _UIntType __upper_mask = (~_UIntType()) << __r;
402  const _UIntType __lower_mask = ~__upper_mask;
403 
404  for (size_t __k = 0; __k < (__n - __m); ++__k)
405  {
406  _UIntType __y = ((_M_x[__k] & __upper_mask)
407  | (_M_x[__k + 1] & __lower_mask));
408  _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
409  ^ ((__y & 0x01) ? __a : 0));
410  }
411 
412  for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
413  {
414  _UIntType __y = ((_M_x[__k] & __upper_mask)
415  | (_M_x[__k + 1] & __lower_mask));
416  _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
417  ^ ((__y & 0x01) ? __a : 0));
418  }
419 
420  _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
421  | (_M_x[0] & __lower_mask));
422  _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
423  ^ ((__y & 0x01) ? __a : 0));
424  _M_p = 0;
425  }
426 
427  template<typename _UIntType, size_t __w,
428  size_t __n, size_t __m, size_t __r,
429  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
430  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
431  _UIntType __f>
432  void
433  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
434  __s, __b, __t, __c, __l, __f>::
435  discard(unsigned long long __z)
436  {
437  while (__z > state_size - _M_p)
438  {
439  __z -= state_size - _M_p;
440  _M_gen_rand();
441  }
442  _M_p += __z;
443  }
444 
445  template<typename _UIntType, size_t __w,
446  size_t __n, size_t __m, size_t __r,
447  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
448  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
449  _UIntType __f>
450  typename
451  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
452  __s, __b, __t, __c, __l, __f>::result_type
453  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
454  __s, __b, __t, __c, __l, __f>::
455  operator()()
456  {
457  // Reload the vector - cost is O(n) amortized over n calls.
458  if (_M_p >= state_size)
459  _M_gen_rand();
460 
461  // Calculate o(x(i)).
462  result_type __z = _M_x[_M_p++];
463  __z ^= (__z >> __u) & __d;
464  __z ^= (__z << __s) & __b;
465  __z ^= (__z << __t) & __c;
466  __z ^= (__z >> __l);
467 
468  return __z;
469  }
470 
471  template<typename _UIntType, size_t __w,
472  size_t __n, size_t __m, size_t __r,
473  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
474  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
475  _UIntType __f, typename _CharT, typename _Traits>
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
527 
528  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
529  constexpr size_t
531 
532  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
533  constexpr size_t
535 
536  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
537  constexpr _UIntType
539 
540  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
541  void
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>
570  typename std::enable_if<std::is_class<_Sseq>::value>::type
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>
596  result_type
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>
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
679 
680  template<typename _RandomNumberEngine, size_t __p, size_t __r>
681  constexpr size_t
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>
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>
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 
754  typedef typename std::common_type<_Eresult_type, result_type>::type
755  __ctype;
756  const unsigned __cdig = std::numeric_limits<__ctype>::digits;
757 
758  unsigned __n, __n0;
759  __ctype __s0, __s1, __y0, __y1;
760 
761  for (size_t __i = 0; __i < 2; ++__i)
762  {
763  __n = (__w + __m - 1) / __m + __i;
764  __n0 = __n - __w % __n;
765  const unsigned __w0 = __w / __n; // __w0 <= __m
766 
767  __s0 = 0;
768  __s1 = 0;
769  if (__w0 < __cdig)
770  {
771  __s0 = __ctype(1) << __w0;
772  __s1 = __s0 << 1;
773  }
774 
775  __y0 = 0;
776  __y1 = 0;
777  if (__r)
778  {
779  __y0 = __s0 * (__r / __s0);
780  if (__s1)
781  __y1 = __s1 * (__r / __s1);
782 
783  if (__r - __y0 <= __y0 / __n)
784  break;
785  }
786  else
787  break;
788  }
789 
790  result_type __sum = 0;
791  for (size_t __k = 0; __k < __n0; ++__k)
792  {
793  __ctype __u;
794  do
795  __u = _M_b() - _M_b.min();
796  while (__y0 && __u >= __y0);
797  __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
798  }
799  for (size_t __k = __n0; __k < __n; ++__k)
800  {
801  __ctype __u;
802  do
803  __u = _M_b() - _M_b.min();
804  while (__y1 && __u >= __y1);
805  __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
806  }
807  return __sum;
808  }
809 
810 
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>
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, typename _CharT, typename _Traits>
877  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
879  {
880  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
881  typedef typename __ostream_type::ios_base __ios_base;
882 
883  const typename __ios_base::fmtflags __flags = __os.flags();
884  const _CharT __fill = __os.fill();
885  const _CharT __space = __os.widen(' ');
887  __os.fill(__space);
888 
889  __os << __x.a() << __space << __x.b();
890 
891  __os.flags(__flags);
892  __os.fill(__fill);
893  return __os;
894  }
895 
896  template<typename _IntType, typename _CharT, typename _Traits>
900  {
901  typedef std::basic_istream<_CharT, _Traits> __istream_type;
902  typedef typename __istream_type::ios_base __ios_base;
903 
904  const typename __ios_base::fmtflags __flags = __is.flags();
906 
907  _IntType __a, __b;
908  __is >> __a >> __b;
910  param_type(__a, __b));
911 
912  __is.flags(__flags);
913  return __is;
914  }
915 
916 
917  template<typename _RealType>
918  template<typename _ForwardIterator,
919  typename _UniformRandomNumberGenerator>
920  void
922  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
923  _UniformRandomNumberGenerator& __urng,
924  const param_type& __p)
925  {
926  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
927  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
928  __aurng(__urng);
929  auto __range = __p.b() - __p.a();
930  while (__f != __t)
931  *__f++ = __aurng() * __range + __p.a();
932  }
933 
934  template<typename _RealType, typename _CharT, typename _Traits>
936  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
938  {
939  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
940  typedef typename __ostream_type::ios_base __ios_base;
941 
942  const typename __ios_base::fmtflags __flags = __os.flags();
943  const _CharT __fill = __os.fill();
944  const std::streamsize __precision = __os.precision();
945  const _CharT __space = __os.widen(' ');
947  __os.fill(__space);
949 
950  __os << __x.a() << __space << __x.b();
951 
952  __os.flags(__flags);
953  __os.fill(__fill);
954  __os.precision(__precision);
955  return __os;
956  }
957 
958  template<typename _RealType, typename _CharT, typename _Traits>
962  {
963  typedef std::basic_istream<_CharT, _Traits> __istream_type;
964  typedef typename __istream_type::ios_base __ios_base;
965 
966  const typename __ios_base::fmtflags __flags = __is.flags();
968 
969  _RealType __a, __b;
970  __is >> __a >> __b;
972  param_type(__a, __b));
973 
974  __is.flags(__flags);
975  return __is;
976  }
977 
978 
979  template<typename _ForwardIterator,
980  typename _UniformRandomNumberGenerator>
981  void
982  std::bernoulli_distribution::
983  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
984  _UniformRandomNumberGenerator& __urng,
985  const param_type& __p)
986  {
987  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
988  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
989  __aurng(__urng);
990  auto __limit = __p.p() * (__aurng.max() - __aurng.min());
991 
992  while (__f != __t)
993  *__f++ = (__aurng() - __aurng.min()) < __limit;
994  }
995 
996  template<typename _CharT, typename _Traits>
998  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
999  const bernoulli_distribution& __x)
1000  {
1001  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1002  typedef typename __ostream_type::ios_base __ios_base;
1003 
1004  const typename __ios_base::fmtflags __flags = __os.flags();
1005  const _CharT __fill = __os.fill();
1006  const std::streamsize __precision = __os.precision();
1008  __os.fill(__os.widen(' '));
1010 
1011  __os << __x.p();
1012 
1013  __os.flags(__flags);
1014  __os.fill(__fill);
1015  __os.precision(__precision);
1016  return __os;
1017  }
1018 
1019 
1020  template<typename _IntType>
1021  template<typename _UniformRandomNumberGenerator>
1024  operator()(_UniformRandomNumberGenerator& __urng,
1025  const param_type& __param)
1026  {
1027  // About the epsilon thing see this thread:
1028  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1029  const double __naf =
1031  // The largest _RealType convertible to _IntType.
1032  const double __thr =
1034  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1035  __aurng(__urng);
1036 
1037  double __cand;
1038  do
1039  __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1040  while (__cand >= __thr);
1041 
1042  return result_type(__cand + __naf);
1043  }
1044 
1045  template<typename _IntType>
1046  template<typename _ForwardIterator,
1047  typename _UniformRandomNumberGenerator>
1048  void
1050  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1051  _UniformRandomNumberGenerator& __urng,
1052  const param_type& __param)
1053  {
1054  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1055  // About the epsilon thing see this thread:
1056  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1057  const double __naf =
1059  // The largest _RealType convertible to _IntType.
1060  const double __thr =
1062  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1063  __aurng(__urng);
1064 
1065  while (__f != __t)
1066  {
1067  double __cand;
1068  do
1069  __cand = std::floor(std::log(1.0 - __aurng())
1070  / __param._M_log_1_p);
1071  while (__cand >= __thr);
1072 
1073  *__f++ = __cand + __naf;
1074  }
1075  }
1076 
1077  template<typename _IntType,
1078  typename _CharT, typename _Traits>
1080  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1082  {
1083  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1084  typedef typename __ostream_type::ios_base __ios_base;
1085 
1086  const typename __ios_base::fmtflags __flags = __os.flags();
1087  const _CharT __fill = __os.fill();
1088  const std::streamsize __precision = __os.precision();
1090  __os.fill(__os.widen(' '));
1092 
1093  __os << __x.p();
1094 
1095  __os.flags(__flags);
1096  __os.fill(__fill);
1097  __os.precision(__precision);
1098  return __os;
1099  }
1100 
1101  template<typename _IntType,
1102  typename _CharT, typename _Traits>
1106  {
1107  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1108  typedef typename __istream_type::ios_base __ios_base;
1109 
1110  const typename __ios_base::fmtflags __flags = __is.flags();
1111  __is.flags(__ios_base::skipws);
1112 
1113  double __p;
1114  __is >> __p;
1116 
1117  __is.flags(__flags);
1118  return __is;
1119  }
1120 
1121  // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1122  template<typename _IntType>
1123  template<typename _UniformRandomNumberGenerator>
1126  operator()(_UniformRandomNumberGenerator& __urng)
1127  {
1128  const double __y = _M_gd(__urng);
1129 
1130  // XXX Is the constructor too slow?
1132  return __poisson(__urng);
1133  }
1134 
1135  template<typename _IntType>
1136  template<typename _UniformRandomNumberGenerator>
1139  operator()(_UniformRandomNumberGenerator& __urng,
1140  const param_type& __p)
1141  {
1143  param_type;
1144 
1145  const double __y =
1146  _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1147 
1149  return __poisson(__urng);
1150  }
1151 
1152  template<typename _IntType>
1153  template<typename _ForwardIterator,
1154  typename _UniformRandomNumberGenerator>
1155  void
1157  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1158  _UniformRandomNumberGenerator& __urng)
1159  {
1160  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1161  while (__f != __t)
1162  {
1163  const double __y = _M_gd(__urng);
1164 
1165  // XXX Is the constructor too slow?
1167  *__f++ = __poisson(__urng);
1168  }
1169  }
1170 
1171  template<typename _IntType>
1172  template<typename _ForwardIterator,
1173  typename _UniformRandomNumberGenerator>
1174  void
1176  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1177  _UniformRandomNumberGenerator& __urng,
1178  const param_type& __p)
1179  {
1180  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1182  __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1183 
1184  while (__f != __t)
1185  {
1186  const double __y = _M_gd(__urng, __p2);
1187 
1189  *__f++ = __poisson(__urng);
1190  }
1191  }
1192 
1193  template<typename _IntType, typename _CharT, typename _Traits>
1195  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1197  {
1198  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1199  typedef typename __ostream_type::ios_base __ios_base;
1200 
1201  const typename __ios_base::fmtflags __flags = __os.flags();
1202  const _CharT __fill = __os.fill();
1203  const std::streamsize __precision = __os.precision();
1204  const _CharT __space = __os.widen(' ');
1206  __os.fill(__os.widen(' '));
1208 
1209  __os << __x.k() << __space << __x.p()
1210  << __space << __x._M_gd;
1211 
1212  __os.flags(__flags);
1213  __os.fill(__fill);
1214  __os.precision(__precision);
1215  return __os;
1216  }
1217 
1218  template<typename _IntType, typename _CharT, typename _Traits>
1222  {
1223  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1224  typedef typename __istream_type::ios_base __ios_base;
1225 
1226  const typename __ios_base::fmtflags __flags = __is.flags();
1227  __is.flags(__ios_base::skipws);
1228 
1229  _IntType __k;
1230  double __p;
1231  __is >> __k >> __p >> __x._M_gd;
1233  param_type(__k, __p));
1234 
1235  __is.flags(__flags);
1236  return __is;
1237  }
1238 
1239 
1240  template<typename _IntType>
1241  void
1244  {
1245 #if _GLIBCXX_USE_C99_MATH_TR1
1246  if (_M_mean >= 12)
1247  {
1248  const double __m = std::floor(_M_mean);
1249  _M_lm_thr = std::log(_M_mean);
1250  _M_lfm = std::lgamma(__m + 1);
1251  _M_sm = std::sqrt(__m);
1252 
1253  const double __pi_4 = 0.7853981633974483096156608458198757L;
1254  const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1255  / __pi_4));
1256  _M_d = std::round(std::max<double>(6.0, std::min(__m, __dx)));
1257  const double __cx = 2 * __m + _M_d;
1258  _M_scx = std::sqrt(__cx / 2);
1259  _M_1cx = 1 / __cx;
1260 
1261  _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1262  _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1263  / _M_d;
1264  }
1265  else
1266 #endif
1267  _M_lm_thr = std::exp(-_M_mean);
1268  }
1269 
1270  /**
1271  * A rejection algorithm when mean >= 12 and a simple method based
1272  * upon the multiplication of uniform random variates otherwise.
1273  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1274  * is defined.
1275  *
1276  * Reference:
1277  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1278  * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1279  */
1280  template<typename _IntType>
1281  template<typename _UniformRandomNumberGenerator>
1284  operator()(_UniformRandomNumberGenerator& __urng,
1285  const param_type& __param)
1286  {
1287  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1288  __aurng(__urng);
1289 #if _GLIBCXX_USE_C99_MATH_TR1
1290  if (__param.mean() >= 12)
1291  {
1292  double __x;
1293 
1294  // See comments above...
1295  const double __naf =
1297  const double __thr =
1299 
1300  const double __m = std::floor(__param.mean());
1301  // sqrt(pi / 2)
1302  const double __spi_2 = 1.2533141373155002512078826424055226L;
1303  const double __c1 = __param._M_sm * __spi_2;
1304  const double __c2 = __param._M_c2b + __c1;
1305  const double __c3 = __c2 + 1;
1306  const double __c4 = __c3 + 1;
1307  // e^(1 / 78)
1308  const double __e178 = 1.0129030479320018583185514777512983L;
1309  const double __c5 = __c4 + __e178;
1310  const double __c = __param._M_cb + __c5;
1311  const double __2cx = 2 * (2 * __m + __param._M_d);
1312 
1313  bool __reject = true;
1314  do
1315  {
1316  const double __u = __c * __aurng();
1317  const double __e = -std::log(1.0 - __aurng());
1318 
1319  double __w = 0.0;
1320 
1321  if (__u <= __c1)
1322  {
1323  const double __n = _M_nd(__urng);
1324  const double __y = -std::abs(__n) * __param._M_sm - 1;
1325  __x = std::floor(__y);
1326  __w = -__n * __n / 2;
1327  if (__x < -__m)
1328  continue;
1329  }
1330  else if (__u <= __c2)
1331  {
1332  const double __n = _M_nd(__urng);
1333  const double __y = 1 + std::abs(__n) * __param._M_scx;
1334  __x = std::ceil(__y);
1335  __w = __y * (2 - __y) * __param._M_1cx;
1336  if (__x > __param._M_d)
1337  continue;
1338  }
1339  else if (__u <= __c3)
1340  // NB: This case not in the book, nor in the Errata,
1341  // but should be ok...
1342  __x = -1;
1343  else if (__u <= __c4)
1344  __x = 0;
1345  else if (__u <= __c5)
1346  __x = 1;
1347  else
1348  {
1349  const double __v = -std::log(1.0 - __aurng());
1350  const double __y = __param._M_d
1351  + __v * __2cx / __param._M_d;
1352  __x = std::ceil(__y);
1353  __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1354  }
1355 
1356  __reject = (__w - __e - __x * __param._M_lm_thr
1357  > __param._M_lfm - std::lgamma(__x + __m + 1));
1358 
1359  __reject |= __x + __m >= __thr;
1360 
1361  } while (__reject);
1362 
1363  return result_type(__x + __m + __naf);
1364  }
1365  else
1366 #endif
1367  {
1368  _IntType __x = 0;
1369  double __prod = 1.0;
1370 
1371  do
1372  {
1373  __prod *= __aurng();
1374  __x += 1;
1375  }
1376  while (__prod > __param._M_lm_thr);
1377 
1378  return __x - 1;
1379  }
1380  }
1381 
1382  template<typename _IntType>
1383  template<typename _ForwardIterator,
1384  typename _UniformRandomNumberGenerator>
1385  void
1387  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1388  _UniformRandomNumberGenerator& __urng,
1389  const param_type& __param)
1390  {
1391  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1392  // We could duplicate everything from operator()...
1393  while (__f != __t)
1394  *__f++ = this->operator()(__urng, __param);
1395  }
1396 
1397  template<typename _IntType,
1398  typename _CharT, typename _Traits>
1400  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1401  const poisson_distribution<_IntType>& __x)
1402  {
1403  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1404  typedef typename __ostream_type::ios_base __ios_base;
1405 
1406  const typename __ios_base::fmtflags __flags = __os.flags();
1407  const _CharT __fill = __os.fill();
1408  const std::streamsize __precision = __os.precision();
1409  const _CharT __space = __os.widen(' ');
1411  __os.fill(__space);
1413 
1414  __os << __x.mean() << __space << __x._M_nd;
1415 
1416  __os.flags(__flags);
1417  __os.fill(__fill);
1418  __os.precision(__precision);
1419  return __os;
1420  }
1421 
1422  template<typename _IntType,
1423  typename _CharT, typename _Traits>
1427  {
1428  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1429  typedef typename __istream_type::ios_base __ios_base;
1430 
1431  const typename __ios_base::fmtflags __flags = __is.flags();
1432  __is.flags(__ios_base::skipws);
1433 
1434  double __mean;
1435  __is >> __mean >> __x._M_nd;
1436  __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1437 
1438  __is.flags(__flags);
1439  return __is;
1440  }
1441 
1442 
1443  template<typename _IntType>
1444  void
1447  {
1448  const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1449 
1450  _M_easy = true;
1451 
1452 #if _GLIBCXX_USE_C99_MATH_TR1
1453  if (_M_t * __p12 >= 8)
1454  {
1455  _M_easy = false;
1456  const double __np = std::floor(_M_t * __p12);
1457  const double __pa = __np / _M_t;
1458  const double __1p = 1 - __pa;
1459 
1460  const double __pi_4 = 0.7853981633974483096156608458198757L;
1461  const double __d1x =
1462  std::sqrt(__np * __1p * std::log(32 * __np
1463  / (81 * __pi_4 * __1p)));
1464  _M_d1 = std::round(std::max<double>(1.0, __d1x));
1465  const double __d2x =
1466  std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1467  / (__pi_4 * __pa)));
1468  _M_d2 = std::round(std::max<double>(1.0, __d2x));
1469 
1470  // sqrt(pi / 2)
1471  const double __spi_2 = 1.2533141373155002512078826424055226L;
1472  _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1473  _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1474  _M_c = 2 * _M_d1 / __np;
1475  _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1476  const double __a12 = _M_a1 + _M_s2 * __spi_2;
1477  const double __s1s = _M_s1 * _M_s1;
1478  _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1479  * 2 * __s1s / _M_d1
1480  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1481  const double __s2s = _M_s2 * _M_s2;
1482  _M_s = (_M_a123 + 2 * __s2s / _M_d2
1483  * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1484  _M_lf = (std::lgamma(__np + 1)
1485  + std::lgamma(_M_t - __np + 1));
1486  _M_lp1p = std::log(__pa / __1p);
1487 
1488  _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1489  }
1490  else
1491 #endif
1492  _M_q = -std::log(1 - __p12);
1493  }
1494 
1495  template<typename _IntType>
1496  template<typename _UniformRandomNumberGenerator>
1499  _M_waiting(_UniformRandomNumberGenerator& __urng,
1500  _IntType __t, double __q)
1501  {
1502  _IntType __x = 0;
1503  double __sum = 0.0;
1504  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1505  __aurng(__urng);
1506 
1507  do
1508  {
1509  if (__t == __x)
1510  return __x;
1511  const double __e = -std::log(1.0 - __aurng());
1512  __sum += __e / (__t - __x);
1513  __x += 1;
1514  }
1515  while (__sum <= __q);
1516 
1517  return __x - 1;
1518  }
1519 
1520  /**
1521  * A rejection algorithm when t * p >= 8 and a simple waiting time
1522  * method - the second in the referenced book - otherwise.
1523  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1524  * is defined.
1525  *
1526  * Reference:
1527  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1528  * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1529  */
1530  template<typename _IntType>
1531  template<typename _UniformRandomNumberGenerator>
1534  operator()(_UniformRandomNumberGenerator& __urng,
1535  const param_type& __param)
1536  {
1537  result_type __ret;
1538  const _IntType __t = __param.t();
1539  const double __p = __param.p();
1540  const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1541  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1542  __aurng(__urng);
1543 
1544 #if _GLIBCXX_USE_C99_MATH_TR1
1545  if (!__param._M_easy)
1546  {
1547  double __x;
1548 
1549  // See comments above...
1550  const double __naf =
1552  const double __thr =
1554 
1555  const double __np = std::floor(__t * __p12);
1556 
1557  // sqrt(pi / 2)
1558  const double __spi_2 = 1.2533141373155002512078826424055226L;
1559  const double __a1 = __param._M_a1;
1560  const double __a12 = __a1 + __param._M_s2 * __spi_2;
1561  const double __a123 = __param._M_a123;
1562  const double __s1s = __param._M_s1 * __param._M_s1;
1563  const double __s2s = __param._M_s2 * __param._M_s2;
1564 
1565  bool __reject;
1566  do
1567  {
1568  const double __u = __param._M_s * __aurng();
1569 
1570  double __v;
1571 
1572  if (__u <= __a1)
1573  {
1574  const double __n = _M_nd(__urng);
1575  const double __y = __param._M_s1 * std::abs(__n);
1576  __reject = __y >= __param._M_d1;
1577  if (!__reject)
1578  {
1579  const double __e = -std::log(1.0 - __aurng());
1580  __x = std::floor(__y);
1581  __v = -__e - __n * __n / 2 + __param._M_c;
1582  }
1583  }
1584  else if (__u <= __a12)
1585  {
1586  const double __n = _M_nd(__urng);
1587  const double __y = __param._M_s2 * std::abs(__n);
1588  __reject = __y >= __param._M_d2;
1589  if (!__reject)
1590  {
1591  const double __e = -std::log(1.0 - __aurng());
1592  __x = std::floor(-__y);
1593  __v = -__e - __n * __n / 2;
1594  }
1595  }
1596  else if (__u <= __a123)
1597  {
1598  const double __e1 = -std::log(1.0 - __aurng());
1599  const double __e2 = -std::log(1.0 - __aurng());
1600 
1601  const double __y = __param._M_d1
1602  + 2 * __s1s * __e1 / __param._M_d1;
1603  __x = std::floor(__y);
1604  __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1605  -__y / (2 * __s1s)));
1606  __reject = false;
1607  }
1608  else
1609  {
1610  const double __e1 = -std::log(1.0 - __aurng());
1611  const double __e2 = -std::log(1.0 - __aurng());
1612 
1613  const double __y = __param._M_d2
1614  + 2 * __s2s * __e1 / __param._M_d2;
1615  __x = std::floor(-__y);
1616  __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1617  __reject = false;
1618  }
1619 
1620  __reject = __reject || __x < -__np || __x > __t - __np;
1621  if (!__reject)
1622  {
1623  const double __lfx =
1624  std::lgamma(__np + __x + 1)
1625  + std::lgamma(__t - (__np + __x) + 1);
1626  __reject = __v > __param._M_lf - __lfx
1627  + __x * __param._M_lp1p;
1628  }
1629 
1630  __reject |= __x + __np >= __thr;
1631  }
1632  while (__reject);
1633 
1634  __x += __np + __naf;
1635 
1636  const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1637  __param._M_q);
1638  __ret = _IntType(__x) + __z;
1639  }
1640  else
1641 #endif
1642  __ret = _M_waiting(__urng, __t, __param._M_q);
1643 
1644  if (__p12 != __p)
1645  __ret = __t - __ret;
1646  return __ret;
1647  }
1648 
1649  template<typename _IntType>
1650  template<typename _ForwardIterator,
1651  typename _UniformRandomNumberGenerator>
1652  void
1654  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1655  _UniformRandomNumberGenerator& __urng,
1656  const param_type& __param)
1657  {
1658  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1659  // We could duplicate everything from operator()...
1660  while (__f != __t)
1661  *__f++ = this->operator()(__urng, __param);
1662  }
1663 
1664  template<typename _IntType,
1665  typename _CharT, typename _Traits>
1667  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1669  {
1670  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1671  typedef typename __ostream_type::ios_base __ios_base;
1672 
1673  const typename __ios_base::fmtflags __flags = __os.flags();
1674  const _CharT __fill = __os.fill();
1675  const std::streamsize __precision = __os.precision();
1676  const _CharT __space = __os.widen(' ');
1678  __os.fill(__space);
1680 
1681  __os << __x.t() << __space << __x.p()
1682  << __space << __x._M_nd;
1683 
1684  __os.flags(__flags);
1685  __os.fill(__fill);
1686  __os.precision(__precision);
1687  return __os;
1688  }
1689 
1690  template<typename _IntType,
1691  typename _CharT, typename _Traits>
1695  {
1696  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1697  typedef typename __istream_type::ios_base __ios_base;
1698 
1699  const typename __ios_base::fmtflags __flags = __is.flags();
1701 
1702  _IntType __t;
1703  double __p;
1704  __is >> __t >> __p >> __x._M_nd;
1705  __x.param(typename binomial_distribution<_IntType>::
1706  param_type(__t, __p));
1707 
1708  __is.flags(__flags);
1709  return __is;
1710  }
1711 
1712 
1713  template<typename _RealType>
1714  template<typename _ForwardIterator,
1715  typename _UniformRandomNumberGenerator>
1716  void
1718  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1719  _UniformRandomNumberGenerator& __urng,
1720  const param_type& __p)
1721  {
1722  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1723  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1724  __aurng(__urng);
1725  while (__f != __t)
1726  *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1727  }
1728 
1729  template<typename _RealType, typename _CharT, typename _Traits>
1731  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1733  {
1734  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1735  typedef typename __ostream_type::ios_base __ios_base;
1736 
1737  const typename __ios_base::fmtflags __flags = __os.flags();
1738  const _CharT __fill = __os.fill();
1739  const std::streamsize __precision = __os.precision();
1741  __os.fill(__os.widen(' '));
1743 
1744  __os << __x.lambda();
1745 
1746  __os.flags(__flags);
1747  __os.fill(__fill);
1748  __os.precision(__precision);
1749  return __os;
1750  }
1751 
1752  template<typename _RealType, typename _CharT, typename _Traits>
1756  {
1757  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1758  typedef typename __istream_type::ios_base __ios_base;
1759 
1760  const typename __ios_base::fmtflags __flags = __is.flags();
1762 
1763  _RealType __lambda;
1764  __is >> __lambda;
1766  param_type(__lambda));
1767 
1768  __is.flags(__flags);
1769  return __is;
1770  }
1771 
1772 
1773  /**
1774  * Polar method due to Marsaglia.
1775  *
1776  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1777  * New York, 1986, Ch. V, Sect. 4.4.
1778  */
1779  template<typename _RealType>
1780  template<typename _UniformRandomNumberGenerator>
1783  operator()(_UniformRandomNumberGenerator& __urng,
1784  const param_type& __param)
1785  {
1786  result_type __ret;
1787  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1788  __aurng(__urng);
1789 
1790  if (_M_saved_available)
1791  {
1792  _M_saved_available = false;
1793  __ret = _M_saved;
1794  }
1795  else
1796  {
1797  result_type __x, __y, __r2;
1798  do
1799  {
1800  __x = result_type(2.0) * __aurng() - 1.0;
1801  __y = result_type(2.0) * __aurng() - 1.0;
1802  __r2 = __x * __x + __y * __y;
1803  }
1804  while (__r2 > 1.0 || __r2 == 0.0);
1805 
1806  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1807  _M_saved = __x * __mult;
1808  _M_saved_available = true;
1809  __ret = __y * __mult;
1810  }
1811 
1812  __ret = __ret * __param.stddev() + __param.mean();
1813  return __ret;
1814  }
1815 
1816  template<typename _RealType>
1817  template<typename _ForwardIterator,
1818  typename _UniformRandomNumberGenerator>
1819  void
1821  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1822  _UniformRandomNumberGenerator& __urng,
1823  const param_type& __param)
1824  {
1825  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1826 
1827  if (__f == __t)
1828  return;
1829 
1830  if (_M_saved_available)
1831  {
1832  _M_saved_available = false;
1833  *__f++ = _M_saved * __param.stddev() + __param.mean();
1834 
1835  if (__f == __t)
1836  return;
1837  }
1838 
1839  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1840  __aurng(__urng);
1841 
1842  while (__f + 1 < __t)
1843  {
1844  result_type __x, __y, __r2;
1845  do
1846  {
1847  __x = result_type(2.0) * __aurng() - 1.0;
1848  __y = result_type(2.0) * __aurng() - 1.0;
1849  __r2 = __x * __x + __y * __y;
1850  }
1851  while (__r2 > 1.0 || __r2 == 0.0);
1852 
1853  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1854  *__f++ = __y * __mult * __param.stddev() + __param.mean();
1855  *__f++ = __x * __mult * __param.stddev() + __param.mean();
1856  }
1857 
1858  if (__f != __t)
1859  {
1860  result_type __x, __y, __r2;
1861  do
1862  {
1863  __x = result_type(2.0) * __aurng() - 1.0;
1864  __y = result_type(2.0) * __aurng() - 1.0;
1865  __r2 = __x * __x + __y * __y;
1866  }
1867  while (__r2 > 1.0 || __r2 == 0.0);
1868 
1869  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1870  _M_saved = __x * __mult;
1871  _M_saved_available = true;
1872  *__f = __y * __mult * __param.stddev() + __param.mean();
1873  }
1874  }
1875 
1876  template<typename _RealType>
1877  bool
1878  operator==(const std::normal_distribution<_RealType>& __d1,
1880  {
1881  if (__d1._M_param == __d2._M_param
1882  && __d1._M_saved_available == __d2._M_saved_available)
1883  {
1884  if (__d1._M_saved_available
1885  && __d1._M_saved == __d2._M_saved)
1886  return true;
1887  else if(!__d1._M_saved_available)
1888  return true;
1889  else
1890  return false;
1891  }
1892  else
1893  return false;
1894  }
1895 
1896  template<typename _RealType, typename _CharT, typename _Traits>
1898  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1899  const normal_distribution<_RealType>& __x)
1900  {
1901  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1902  typedef typename __ostream_type::ios_base __ios_base;
1903 
1904  const typename __ios_base::fmtflags __flags = __os.flags();
1905  const _CharT __fill = __os.fill();
1906  const std::streamsize __precision = __os.precision();
1907  const _CharT __space = __os.widen(' ');
1909  __os.fill(__space);
1911 
1912  __os << __x.mean() << __space << __x.stddev()
1913  << __space << __x._M_saved_available;
1914  if (__x._M_saved_available)
1915  __os << __space << __x._M_saved;
1916 
1917  __os.flags(__flags);
1918  __os.fill(__fill);
1919  __os.precision(__precision);
1920  return __os;
1921  }
1922 
1923  template<typename _RealType, typename _CharT, typename _Traits>
1927  {
1928  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1929  typedef typename __istream_type::ios_base __ios_base;
1930 
1931  const typename __ios_base::fmtflags __flags = __is.flags();
1933 
1934  double __mean, __stddev;
1935  __is >> __mean >> __stddev
1936  >> __x._M_saved_available;
1937  if (__x._M_saved_available)
1938  __is >> __x._M_saved;
1939  __x.param(typename normal_distribution<_RealType>::
1940  param_type(__mean, __stddev));
1941 
1942  __is.flags(__flags);
1943  return __is;
1944  }
1945 
1946 
1947  template<typename _RealType>
1948  template<typename _ForwardIterator,
1949  typename _UniformRandomNumberGenerator>
1950  void
1952  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1953  _UniformRandomNumberGenerator& __urng,
1954  const param_type& __p)
1955  {
1956  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1957  while (__f != __t)
1958  *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
1959  }
1960 
1961  template<typename _RealType, typename _CharT, typename _Traits>
1963  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1965  {
1966  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1967  typedef typename __ostream_type::ios_base __ios_base;
1968 
1969  const typename __ios_base::fmtflags __flags = __os.flags();
1970  const _CharT __fill = __os.fill();
1971  const std::streamsize __precision = __os.precision();
1972  const _CharT __space = __os.widen(' ');
1974  __os.fill(__space);
1976 
1977  __os << __x.m() << __space << __x.s()
1978  << __space << __x._M_nd;
1979 
1980  __os.flags(__flags);
1981  __os.fill(__fill);
1982  __os.precision(__precision);
1983  return __os;
1984  }
1985 
1986  template<typename _RealType, typename _CharT, typename _Traits>
1990  {
1991  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1992  typedef typename __istream_type::ios_base __ios_base;
1993 
1994  const typename __ios_base::fmtflags __flags = __is.flags();
1996 
1997  _RealType __m, __s;
1998  __is >> __m >> __s >> __x._M_nd;
2000  param_type(__m, __s));
2001 
2002  __is.flags(__flags);
2003  return __is;
2004  }
2005 
2006  template<typename _RealType>
2007  template<typename _ForwardIterator,
2008  typename _UniformRandomNumberGenerator>
2009  void
2011  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2012  _UniformRandomNumberGenerator& __urng)
2013  {
2014  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2015  while (__f != __t)
2016  *__f++ = 2 * _M_gd(__urng);
2017  }
2018 
2019  template<typename _RealType>
2020  template<typename _ForwardIterator,
2021  typename _UniformRandomNumberGenerator>
2022  void
2024  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2025  _UniformRandomNumberGenerator& __urng,
2026  const typename
2028  {
2029  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2030  while (__f != __t)
2031  *__f++ = 2 * _M_gd(__urng, __p);
2032  }
2033 
2034  template<typename _RealType, typename _CharT, typename _Traits>
2036  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2038  {
2039  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2040  typedef typename __ostream_type::ios_base __ios_base;
2041 
2042  const typename __ios_base::fmtflags __flags = __os.flags();
2043  const _CharT __fill = __os.fill();
2044  const std::streamsize __precision = __os.precision();
2045  const _CharT __space = __os.widen(' ');
2047  __os.fill(__space);
2049 
2050  __os << __x.n() << __space << __x._M_gd;
2051 
2052  __os.flags(__flags);
2053  __os.fill(__fill);
2054  __os.precision(__precision);
2055  return __os;
2056  }
2057 
2058  template<typename _RealType, typename _CharT, typename _Traits>
2062  {
2063  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2064  typedef typename __istream_type::ios_base __ios_base;
2065 
2066  const typename __ios_base::fmtflags __flags = __is.flags();
2068 
2069  _RealType __n;
2070  __is >> __n >> __x._M_gd;
2072  param_type(__n));
2073 
2074  __is.flags(__flags);
2075  return __is;
2076  }
2077 
2078 
2079  template<typename _RealType>
2080  template<typename _UniformRandomNumberGenerator>
2083  operator()(_UniformRandomNumberGenerator& __urng,
2084  const param_type& __p)
2085  {
2086  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2087  __aurng(__urng);
2088  _RealType __u;
2089  do
2090  __u = __aurng();
2091  while (__u == 0.5);
2092 
2093  const _RealType __pi = 3.1415926535897932384626433832795029L;
2094  return __p.a() + __p.b() * std::tan(__pi * __u);
2095  }
2096 
2097  template<typename _RealType>
2098  template<typename _ForwardIterator,
2099  typename _UniformRandomNumberGenerator>
2100  void
2102  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2103  _UniformRandomNumberGenerator& __urng,
2104  const param_type& __p)
2105  {
2106  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2107  const _RealType __pi = 3.1415926535897932384626433832795029L;
2108  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2109  __aurng(__urng);
2110  while (__f != __t)
2111  {
2112  _RealType __u;
2113  do
2114  __u = __aurng();
2115  while (__u == 0.5);
2116 
2117  *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2118  }
2119  }
2120 
2121  template<typename _RealType, typename _CharT, typename _Traits>
2123  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2124  const cauchy_distribution<_RealType>& __x)
2125  {
2126  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2127  typedef typename __ostream_type::ios_base __ios_base;
2128 
2129  const typename __ios_base::fmtflags __flags = __os.flags();
2130  const _CharT __fill = __os.fill();
2131  const std::streamsize __precision = __os.precision();
2132  const _CharT __space = __os.widen(' ');
2134  __os.fill(__space);
2136 
2137  __os << __x.a() << __space << __x.b();
2138 
2139  __os.flags(__flags);
2140  __os.fill(__fill);
2141  __os.precision(__precision);
2142  return __os;
2143  }
2144 
2145  template<typename _RealType, typename _CharT, typename _Traits>
2149  {
2150  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2151  typedef typename __istream_type::ios_base __ios_base;
2152 
2153  const typename __ios_base::fmtflags __flags = __is.flags();
2155 
2156  _RealType __a, __b;
2157  __is >> __a >> __b;
2158  __x.param(typename cauchy_distribution<_RealType>::
2159  param_type(__a, __b));
2160 
2161  __is.flags(__flags);
2162  return __is;
2163  }
2164 
2165 
2166  template<typename _RealType>
2167  template<typename _ForwardIterator,
2168  typename _UniformRandomNumberGenerator>
2169  void
2171  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2172  _UniformRandomNumberGenerator& __urng)
2173  {
2174  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2175  while (__f != __t)
2176  *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2177  }
2178 
2179  template<typename _RealType>
2180  template<typename _ForwardIterator,
2181  typename _UniformRandomNumberGenerator>
2182  void
2184  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2185  _UniformRandomNumberGenerator& __urng,
2186  const param_type& __p)
2187  {
2188  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2190  param_type;
2191  param_type __p1(__p.m() / 2);
2192  param_type __p2(__p.n() / 2);
2193  while (__f != __t)
2194  *__f++ = ((_M_gd_x(__urng, __p1) * n())
2195  / (_M_gd_y(__urng, __p2) * m()));
2196  }
2197 
2198  template<typename _RealType, typename _CharT, typename _Traits>
2200  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2202  {
2203  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2204  typedef typename __ostream_type::ios_base __ios_base;
2205 
2206  const typename __ios_base::fmtflags __flags = __os.flags();
2207  const _CharT __fill = __os.fill();
2208  const std::streamsize __precision = __os.precision();
2209  const _CharT __space = __os.widen(' ');
2211  __os.fill(__space);
2213 
2214  __os << __x.m() << __space << __x.n()
2215  << __space << __x._M_gd_x << __space << __x._M_gd_y;
2216 
2217  __os.flags(__flags);
2218  __os.fill(__fill);
2219  __os.precision(__precision);
2220  return __os;
2221  }
2222 
2223  template<typename _RealType, typename _CharT, typename _Traits>
2227  {
2228  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2229  typedef typename __istream_type::ios_base __ios_base;
2230 
2231  const typename __ios_base::fmtflags __flags = __is.flags();
2233 
2234  _RealType __m, __n;
2235  __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
2237  param_type(__m, __n));
2238 
2239  __is.flags(__flags);
2240  return __is;
2241  }
2242 
2243 
2244  template<typename _RealType>
2245  template<typename _ForwardIterator,
2246  typename _UniformRandomNumberGenerator>
2247  void
2249  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2250  _UniformRandomNumberGenerator& __urng)
2251  {
2252  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2253  while (__f != __t)
2254  *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2255  }
2256 
2257  template<typename _RealType>
2258  template<typename _ForwardIterator,
2259  typename _UniformRandomNumberGenerator>
2260  void
2262  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2263  _UniformRandomNumberGenerator& __urng,
2264  const param_type& __p)
2265  {
2266  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2268  __p2(__p.n() / 2, 2);
2269  while (__f != __t)
2270  *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2271  }
2272 
2273  template<typename _RealType, typename _CharT, typename _Traits>
2275  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2277  {
2278  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2279  typedef typename __ostream_type::ios_base __ios_base;
2280 
2281  const typename __ios_base::fmtflags __flags = __os.flags();
2282  const _CharT __fill = __os.fill();
2283  const std::streamsize __precision = __os.precision();
2284  const _CharT __space = __os.widen(' ');
2286  __os.fill(__space);
2288 
2289  __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2290 
2291  __os.flags(__flags);
2292  __os.fill(__fill);
2293  __os.precision(__precision);
2294  return __os;
2295  }
2296 
2297  template<typename _RealType, typename _CharT, typename _Traits>
2301  {
2302  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2303  typedef typename __istream_type::ios_base __ios_base;
2304 
2305  const typename __ios_base::fmtflags __flags = __is.flags();
2307 
2308  _RealType __n;
2309  __is >> __n >> __x._M_nd >> __x._M_gd;
2311 
2312  __is.flags(__flags);
2313  return __is;
2314  }
2315 
2316 
2317  template<typename _RealType>
2318  void
2321  {
2322  _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2323 
2324  const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2325  _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2326  }
2327 
2328  /**
2329  * Marsaglia, G. and Tsang, W. W.
2330  * "A Simple Method for Generating Gamma Variables"
2331  * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2332  */
2333  template<typename _RealType>
2334  template<typename _UniformRandomNumberGenerator>
2337  operator()(_UniformRandomNumberGenerator& __urng,
2338  const param_type& __param)
2339  {
2340  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2341  __aurng(__urng);
2342 
2343  result_type __u, __v, __n;
2344  const result_type __a1 = (__param._M_malpha
2345  - _RealType(1.0) / _RealType(3.0));
2346 
2347  do
2348  {
2349  do
2350  {
2351  __n = _M_nd(__urng);
2352  __v = result_type(1.0) + __param._M_a2 * __n;
2353  }
2354  while (__v <= 0.0);
2355 
2356  __v = __v * __v * __v;
2357  __u = __aurng();
2358  }
2359  while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2360  && (std::log(__u) > (0.5 * __n * __n + __a1
2361  * (1.0 - __v + std::log(__v)))));
2362 
2363  if (__param.alpha() == __param._M_malpha)
2364  return __a1 * __v * __param.beta();
2365  else
2366  {
2367  do
2368  __u = __aurng();
2369  while (__u == 0.0);
2370 
2371  return (std::pow(__u, result_type(1.0) / __param.alpha())
2372  * __a1 * __v * __param.beta());
2373  }
2374  }
2375 
2376  template<typename _RealType>
2377  template<typename _ForwardIterator,
2378  typename _UniformRandomNumberGenerator>
2379  void
2381  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2382  _UniformRandomNumberGenerator& __urng,
2383  const param_type& __param)
2384  {
2385  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2386  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2387  __aurng(__urng);
2388 
2389  result_type __u, __v, __n;
2390  const result_type __a1 = (__param._M_malpha
2391  - _RealType(1.0) / _RealType(3.0));
2392 
2393  if (__param.alpha() == __param._M_malpha)
2394  while (__f != __t)
2395  {
2396  do
2397  {
2398  do
2399  {
2400  __n = _M_nd(__urng);
2401  __v = result_type(1.0) + __param._M_a2 * __n;
2402  }
2403  while (__v <= 0.0);
2404 
2405  __v = __v * __v * __v;
2406  __u = __aurng();
2407  }
2408  while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2409  && (std::log(__u) > (0.5 * __n * __n + __a1
2410  * (1.0 - __v + std::log(__v)))));
2411 
2412  *__f++ = __a1 * __v * __param.beta();
2413  }
2414  else
2415  while (__f != __t)
2416  {
2417  do
2418  {
2419  do
2420  {
2421  __n = _M_nd(__urng);
2422  __v = result_type(1.0) + __param._M_a2 * __n;
2423  }
2424  while (__v <= 0.0);
2425 
2426  __v = __v * __v * __v;
2427  __u = __aurng();
2428  }
2429  while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2430  && (std::log(__u) > (0.5 * __n * __n + __a1
2431  * (1.0 - __v + std::log(__v)))));
2432 
2433  do
2434  __u = __aurng();
2435  while (__u == 0.0);
2436 
2437  *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2438  * __a1 * __v * __param.beta());
2439  }
2440  }
2441 
2442  template<typename _RealType, typename _CharT, typename _Traits>
2444  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2445  const gamma_distribution<_RealType>& __x)
2446  {
2447  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2448  typedef typename __ostream_type::ios_base __ios_base;
2449 
2450  const typename __ios_base::fmtflags __flags = __os.flags();
2451  const _CharT __fill = __os.fill();
2452  const std::streamsize __precision = __os.precision();
2453  const _CharT __space = __os.widen(' ');
2455  __os.fill(__space);
2457 
2458  __os << __x.alpha() << __space << __x.beta()
2459  << __space << __x._M_nd;
2460 
2461  __os.flags(__flags);
2462  __os.fill(__fill);
2463  __os.precision(__precision);
2464  return __os;
2465  }
2466 
2467  template<typename _RealType, typename _CharT, typename _Traits>
2471  {
2472  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2473  typedef typename __istream_type::ios_base __ios_base;
2474 
2475  const typename __ios_base::fmtflags __flags = __is.flags();
2477 
2478  _RealType __alpha_val, __beta_val;
2479  __is >> __alpha_val >> __beta_val >> __x._M_nd;
2480  __x.param(typename gamma_distribution<_RealType>::
2481  param_type(__alpha_val, __beta_val));
2482 
2483  __is.flags(__flags);
2484  return __is;
2485  }
2486 
2487 
2488  template<typename _RealType>
2489  template<typename _UniformRandomNumberGenerator>
2492  operator()(_UniformRandomNumberGenerator& __urng,
2493  const param_type& __p)
2494  {
2495  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2496  __aurng(__urng);
2497  return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2498  result_type(1) / __p.a());
2499  }
2500 
2501  template<typename _RealType>
2502  template<typename _ForwardIterator,
2503  typename _UniformRandomNumberGenerator>
2504  void
2506  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2507  _UniformRandomNumberGenerator& __urng,
2508  const param_type& __p)
2509  {
2510  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2511  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2512  __aurng(__urng);
2513  auto __inv_a = result_type(1) / __p.a();
2514 
2515  while (__f != __t)
2516  *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2517  __inv_a);
2518  }
2519 
2520  template<typename _RealType, typename _CharT, typename _Traits>
2522  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2524  {
2525  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2526  typedef typename __ostream_type::ios_base __ios_base;
2527 
2528  const typename __ios_base::fmtflags __flags = __os.flags();
2529  const _CharT __fill = __os.fill();
2530  const std::streamsize __precision = __os.precision();
2531  const _CharT __space = __os.widen(' ');
2533  __os.fill(__space);
2535 
2536  __os << __x.a() << __space << __x.b();
2537 
2538  __os.flags(__flags);
2539  __os.fill(__fill);
2540  __os.precision(__precision);
2541  return __os;
2542  }
2543 
2544  template<typename _RealType, typename _CharT, typename _Traits>
2548  {
2549  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2550  typedef typename __istream_type::ios_base __ios_base;
2551 
2552  const typename __ios_base::fmtflags __flags = __is.flags();
2554 
2555  _RealType __a, __b;
2556  __is >> __a >> __b;
2557  __x.param(typename weibull_distribution<_RealType>::
2558  param_type(__a, __b));
2559 
2560  __is.flags(__flags);
2561  return __is;
2562  }
2563 
2564 
2565  template<typename _RealType>
2566  template<typename _UniformRandomNumberGenerator>
2569  operator()(_UniformRandomNumberGenerator& __urng,
2570  const param_type& __p)
2571  {
2572  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2573  __aurng(__urng);
2574  return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2575  - __aurng()));
2576  }
2577 
2578  template<typename _RealType>
2579  template<typename _ForwardIterator,
2580  typename _UniformRandomNumberGenerator>
2581  void
2583  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2584  _UniformRandomNumberGenerator& __urng,
2585  const param_type& __p)
2586  {
2587  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2588  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2589  __aurng(__urng);
2590 
2591  while (__f != __t)
2592  *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2593  - __aurng()));
2594  }
2595 
2596  template<typename _RealType, typename _CharT, typename _Traits>
2598  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2600  {
2601  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2602  typedef typename __ostream_type::ios_base __ios_base;
2603 
2604  const typename __ios_base::fmtflags __flags = __os.flags();
2605  const _CharT __fill = __os.fill();
2606  const std::streamsize __precision = __os.precision();
2607  const _CharT __space = __os.widen(' ');
2609  __os.fill(__space);
2611 
2612  __os << __x.a() << __space << __x.b();
2613 
2614  __os.flags(__flags);
2615  __os.fill(__fill);
2616  __os.precision(__precision);
2617  return __os;
2618  }
2619 
2620  template<typename _RealType, typename _CharT, typename _Traits>
2624  {
2625  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2626  typedef typename __istream_type::ios_base __ios_base;
2627 
2628  const typename __ios_base::fmtflags __flags = __is.flags();
2630 
2631  _RealType __a, __b;
2632  __is >> __a >> __b;
2634  param_type(__a, __b));
2635 
2636  __is.flags(__flags);
2637  return __is;
2638  }
2639 
2640 
2641  template<typename _IntType>
2642  void
2645  {
2646  if (_M_prob.size() < 2)
2647  {
2648  _M_prob.clear();
2649  return;
2650  }
2651 
2652  const double __sum = std::accumulate(_M_prob.begin(),
2653  _M_prob.end(), 0.0);
2654  // Now normalize the probabilites.
2655  __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2656  __sum);
2657  // Accumulate partial sums.
2658  _M_cp.reserve(_M_prob.size());
2659  std::partial_sum(_M_prob.begin(), _M_prob.end(),
2660  std::back_inserter(_M_cp));
2661  // Make sure the last cumulative probability is one.
2662  _M_cp[_M_cp.size() - 1] = 1.0;
2663  }
2664 
2665  template<typename _IntType>
2666  template<typename _Func>
2668  param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2669  : _M_prob(), _M_cp()
2670  {
2671  const size_t __n = __nw == 0 ? 1 : __nw;
2672  const double __delta = (__xmax - __xmin) / __n;
2673 
2674  _M_prob.reserve(__n);
2675  for (size_t __k = 0; __k < __nw; ++__k)
2676  _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2677 
2678  _M_initialize();
2679  }
2680 
2681  template<typename _IntType>
2682  template<typename _UniformRandomNumberGenerator>
2685  operator()(_UniformRandomNumberGenerator& __urng,
2686  const param_type& __param)
2687  {
2688  if (__param._M_cp.empty())
2689  return result_type(0);
2690 
2691  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2692  __aurng(__urng);
2693 
2694  const double __p = __aurng();
2695  auto __pos = std::lower_bound(__param._M_cp.begin(),
2696  __param._M_cp.end(), __p);
2697 
2698  return __pos - __param._M_cp.begin();
2699  }
2700 
2701  template<typename _IntType>
2702  template<typename _ForwardIterator,
2703  typename _UniformRandomNumberGenerator>
2704  void
2706  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2707  _UniformRandomNumberGenerator& __urng,
2708  const param_type& __param)
2709  {
2710  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2711 
2712  if (__param._M_cp.empty())
2713  {
2714  while (__f != __t)
2715  *__f++ = result_type(0);
2716  return;
2717  }
2718 
2719  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2720  __aurng(__urng);
2721 
2722  while (__f != __t)
2723  {
2724  const double __p = __aurng();
2725  auto __pos = std::lower_bound(__param._M_cp.begin(),
2726  __param._M_cp.end(), __p);
2727 
2728  *__f++ = __pos - __param._M_cp.begin();
2729  }
2730  }
2731 
2732  template<typename _IntType, typename _CharT, typename _Traits>
2734  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2736  {
2737  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2738  typedef typename __ostream_type::ios_base __ios_base;
2739 
2740  const typename __ios_base::fmtflags __flags = __os.flags();
2741  const _CharT __fill = __os.fill();
2742  const std::streamsize __precision = __os.precision();
2743  const _CharT __space = __os.widen(' ');
2745  __os.fill(__space);
2747 
2748  std::vector<double> __prob = __x.probabilities();
2749  __os << __prob.size();
2750  for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2751  __os << __space << *__dit;
2752 
2753  __os.flags(__flags);
2754  __os.fill(__fill);
2755  __os.precision(__precision);
2756  return __os;
2757  }
2758 
2759  template<typename _IntType, typename _CharT, typename _Traits>
2763  {
2764  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2765  typedef typename __istream_type::ios_base __ios_base;
2766 
2767  const typename __ios_base::fmtflags __flags = __is.flags();
2769 
2770  size_t __n;
2771  __is >> __n;
2772 
2773  std::vector<double> __prob_vec;
2774  __prob_vec.reserve(__n);
2775  for (; __n != 0; --__n)
2776  {
2777  double __prob;
2778  __is >> __prob;
2779  __prob_vec.push_back(__prob);
2780  }
2781 
2782  __x.param(typename discrete_distribution<_IntType>::
2783  param_type(__prob_vec.begin(), __prob_vec.end()));
2784 
2785  __is.flags(__flags);
2786  return __is;
2787  }
2788 
2789 
2790  template<typename _RealType>
2791  void
2794  {
2795  if (_M_int.size() < 2
2796  || (_M_int.size() == 2
2797  && _M_int[0] == _RealType(0)
2798  && _M_int[1] == _RealType(1)))
2799  {
2800  _M_int.clear();
2801  _M_den.clear();
2802  return;
2803  }
2804 
2805  const double __sum = std::accumulate(_M_den.begin(),
2806  _M_den.end(), 0.0);
2807 
2808  __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2809  __sum);
2810 
2811  _M_cp.reserve(_M_den.size());
2812  std::partial_sum(_M_den.begin(), _M_den.end(),
2813  std::back_inserter(_M_cp));
2814 
2815  // Make sure the last cumulative probability is one.
2816  _M_cp[_M_cp.size() - 1] = 1.0;
2817 
2818  for (size_t __k = 0; __k < _M_den.size(); ++__k)
2819  _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2820  }
2821 
2822  template<typename _RealType>
2823  template<typename _InputIteratorB, typename _InputIteratorW>
2825  param_type(_InputIteratorB __bbegin,
2826  _InputIteratorB __bend,
2827  _InputIteratorW __wbegin)
2828  : _M_int(), _M_den(), _M_cp()
2829  {
2830  if (__bbegin != __bend)
2831  {
2832  for (;;)
2833  {
2834  _M_int.push_back(*__bbegin);
2835  ++__bbegin;
2836  if (__bbegin == __bend)
2837  break;
2838 
2839  _M_den.push_back(*__wbegin);
2840  ++__wbegin;
2841  }
2842  }
2843 
2844  _M_initialize();
2845  }
2846 
2847  template<typename _RealType>
2848  template<typename _Func>
2850  param_type(initializer_list<_RealType> __bl, _Func __fw)
2851  : _M_int(), _M_den(), _M_cp()
2852  {
2853  _M_int.reserve(__bl.size());
2854  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2855  _M_int.push_back(*__biter);
2856 
2857  _M_den.reserve(_M_int.size() - 1);
2858  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2859  _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2860 
2861  _M_initialize();
2862  }
2863 
2864  template<typename _RealType>
2865  template<typename _Func>
2867  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2868  : _M_int(), _M_den(), _M_cp()
2869  {
2870  const size_t __n = __nw == 0 ? 1 : __nw;
2871  const _RealType __delta = (__xmax - __xmin) / __n;
2872 
2873  _M_int.reserve(__n + 1);
2874  for (size_t __k = 0; __k <= __nw; ++__k)
2875  _M_int.push_back(__xmin + __k * __delta);
2876 
2877  _M_den.reserve(__n);
2878  for (size_t __k = 0; __k < __nw; ++__k)
2879  _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2880 
2881  _M_initialize();
2882  }
2883 
2884  template<typename _RealType>
2885  template<typename _UniformRandomNumberGenerator>
2888  operator()(_UniformRandomNumberGenerator& __urng,
2889  const param_type& __param)
2890  {
2891  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2892  __aurng(__urng);
2893 
2894  const double __p = __aurng();
2895  if (__param._M_cp.empty())
2896  return __p;
2897 
2898  auto __pos = std::lower_bound(__param._M_cp.begin(),
2899  __param._M_cp.end(), __p);
2900  const size_t __i = __pos - __param._M_cp.begin();
2901 
2902  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2903 
2904  return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2905  }
2906 
2907  template<typename _RealType>
2908  template<typename _ForwardIterator,
2909  typename _UniformRandomNumberGenerator>
2910  void
2912  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2913  _UniformRandomNumberGenerator& __urng,
2914  const param_type& __param)
2915  {
2916  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2917  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2918  __aurng(__urng);
2919 
2920  if (__param._M_cp.empty())
2921  {
2922  while (__f != __t)
2923  *__f++ = __aurng();
2924  return;
2925  }
2926 
2927  while (__f != __t)
2928  {
2929  const double __p = __aurng();
2930 
2931  auto __pos = std::lower_bound(__param._M_cp.begin(),
2932  __param._M_cp.end(), __p);
2933  const size_t __i = __pos - __param._M_cp.begin();
2934 
2935  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2936 
2937  *__f++ = (__param._M_int[__i]
2938  + (__p - __pref) / __param._M_den[__i]);
2939  }
2940  }
2941 
2942  template<typename _RealType, typename _CharT, typename _Traits>
2944  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2946  {
2947  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2948  typedef typename __ostream_type::ios_base __ios_base;
2949 
2950  const typename __ios_base::fmtflags __flags = __os.flags();
2951  const _CharT __fill = __os.fill();
2952  const std::streamsize __precision = __os.precision();
2953  const _CharT __space = __os.widen(' ');
2955  __os.fill(__space);
2957 
2958  std::vector<_RealType> __int = __x.intervals();
2959  __os << __int.size() - 1;
2960 
2961  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2962  __os << __space << *__xit;
2963 
2964  std::vector<double> __den = __x.densities();
2965  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2966  __os << __space << *__dit;
2967 
2968  __os.flags(__flags);
2969  __os.fill(__fill);
2970  __os.precision(__precision);
2971  return __os;
2972  }
2973 
2974  template<typename _RealType, typename _CharT, typename _Traits>
2978  {
2979  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2980  typedef typename __istream_type::ios_base __ios_base;
2981 
2982  const typename __ios_base::fmtflags __flags = __is.flags();
2984 
2985  size_t __n;
2986  __is >> __n;
2987 
2988  std::vector<_RealType> __int_vec;
2989  __int_vec.reserve(__n + 1);
2990  for (size_t __i = 0; __i <= __n; ++__i)
2991  {
2992  _RealType __int;
2993  __is >> __int;
2994  __int_vec.push_back(__int);
2995  }
2996 
2997  std::vector<double> __den_vec;
2998  __den_vec.reserve(__n);
2999  for (size_t __i = 0; __i < __n; ++__i)
3000  {
3001  double __den;
3002  __is >> __den;
3003  __den_vec.push_back(__den);
3004  }
3005 
3007  param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3008 
3009  __is.flags(__flags);
3010  return __is;
3011  }
3012 
3013 
3014  template<typename _RealType>
3015  void
3018  {
3019  if (_M_int.size() < 2
3020  || (_M_int.size() == 2
3021  && _M_int[0] == _RealType(0)
3022  && _M_int[1] == _RealType(1)
3023  && _M_den[0] == _M_den[1]))
3024  {
3025  _M_int.clear();
3026  _M_den.clear();
3027  return;
3028  }
3029 
3030  double __sum = 0.0;
3031  _M_cp.reserve(_M_int.size() - 1);
3032  _M_m.reserve(_M_int.size() - 1);
3033  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3034  {
3035  const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3036  __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3037  _M_cp.push_back(__sum);
3038  _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3039  }
3040 
3041  // Now normalize the densities...
3042  __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3043  __sum);
3044  // ... and partial sums...
3045  __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3046  // ... and slopes.
3047  __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3048 
3049  // Make sure the last cumulative probablility is one.
3050  _M_cp[_M_cp.size() - 1] = 1.0;
3051  }
3052 
3053  template<typename _RealType>
3054  template<typename _InputIteratorB, typename _InputIteratorW>
3056  param_type(_InputIteratorB __bbegin,
3057  _InputIteratorB __bend,
3058  _InputIteratorW __wbegin)
3059  : _M_int(), _M_den(), _M_cp(), _M_m()
3060  {
3061  for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
3062  {
3063  _M_int.push_back(*__bbegin);
3064  _M_den.push_back(*__wbegin);
3065  }
3066 
3067  _M_initialize();
3068  }
3069 
3070  template<typename _RealType>
3071  template<typename _Func>
3073  param_type(initializer_list<_RealType> __bl, _Func __fw)
3074  : _M_int(), _M_den(), _M_cp(), _M_m()
3075  {
3076  _M_int.reserve(__bl.size());
3077  _M_den.reserve(__bl.size());
3078  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3079  {
3080  _M_int.push_back(*__biter);
3081  _M_den.push_back(__fw(*__biter));
3082  }
3083 
3084  _M_initialize();
3085  }
3086 
3087  template<typename _RealType>
3088  template<typename _Func>
3090  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3091  : _M_int(), _M_den(), _M_cp(), _M_m()
3092  {
3093  const size_t __n = __nw == 0 ? 1 : __nw;
3094  const _RealType __delta = (__xmax - __xmin) / __n;
3095 
3096  _M_int.reserve(__n + 1);
3097  _M_den.reserve(__n + 1);
3098  for (size_t __k = 0; __k <= __nw; ++__k)
3099  {
3100  _M_int.push_back(__xmin + __k * __delta);
3101  _M_den.push_back(__fw(_M_int[__k] + __delta));
3102  }
3103 
3104  _M_initialize();
3105  }
3106 
3107  template<typename _RealType>
3108  template<typename _UniformRandomNumberGenerator>
3111  operator()(_UniformRandomNumberGenerator& __urng,
3112  const param_type& __param)
3113  {
3114  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3115  __aurng(__urng);
3116 
3117  const double __p = __aurng();
3118  if (__param._M_cp.empty())
3119  return __p;
3120 
3121  auto __pos = std::lower_bound(__param._M_cp.begin(),
3122  __param._M_cp.end(), __p);
3123  const size_t __i = __pos - __param._M_cp.begin();
3124 
3125  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3126 
3127  const double __a = 0.5 * __param._M_m[__i];
3128  const double __b = __param._M_den[__i];
3129  const double __cm = __p - __pref;
3130 
3131  _RealType __x = __param._M_int[__i];
3132  if (__a == 0)
3133  __x += __cm / __b;
3134  else
3135  {
3136  const double __d = __b * __b + 4.0 * __a * __cm;
3137  __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3138  }
3139 
3140  return __x;
3141  }
3142 
3143  template<typename _RealType>
3144  template<typename _ForwardIterator,
3145  typename _UniformRandomNumberGenerator>
3146  void
3148  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3149  _UniformRandomNumberGenerator& __urng,
3150  const param_type& __param)
3151  {
3152  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3153  // We could duplicate everything from operator()...
3154  while (__f != __t)
3155  *__f++ = this->operator()(__urng, __param);
3156  }
3157 
3158  template<typename _RealType, typename _CharT, typename _Traits>
3160  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3162  {
3163  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
3164  typedef typename __ostream_type::ios_base __ios_base;
3165 
3166  const typename __ios_base::fmtflags __flags = __os.flags();
3167  const _CharT __fill = __os.fill();
3168  const std::streamsize __precision = __os.precision();
3169  const _CharT __space = __os.widen(' ');
3171  __os.fill(__space);
3173 
3174  std::vector<_RealType> __int = __x.intervals();
3175  __os << __int.size() - 1;
3176 
3177  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3178  __os << __space << *__xit;
3179 
3180  std::vector<double> __den = __x.densities();
3181  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3182  __os << __space << *__dit;
3183 
3184  __os.flags(__flags);
3185  __os.fill(__fill);
3186  __os.precision(__precision);
3187  return __os;
3188  }
3189 
3190  template<typename _RealType, typename _CharT, typename _Traits>
3194  {
3195  typedef std::basic_istream<_CharT, _Traits> __istream_type;
3196  typedef typename __istream_type::ios_base __ios_base;
3197 
3198  const typename __ios_base::fmtflags __flags = __is.flags();
3200 
3201  size_t __n;
3202  __is >> __n;
3203 
3204  std::vector<_RealType> __int_vec;
3205  __int_vec.reserve(__n + 1);
3206  for (size_t __i = 0; __i <= __n; ++__i)
3207  {
3208  _RealType __int;
3209  __is >> __int;
3210  __int_vec.push_back(__int);
3211  }
3212 
3213  std::vector<double> __den_vec;
3214  __den_vec.reserve(__n + 1);
3215  for (size_t __i = 0; __i <= __n; ++__i)
3216  {
3217  double __den;
3218  __is >> __den;
3219  __den_vec.push_back(__den);
3220  }
3221 
3223  param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3224 
3225  __is.flags(__flags);
3226  return __is;
3227  }
3228 
3229 
3230  template<typename _IntType>
3232  {
3233  for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3234  _M_v.push_back(__detail::__mod<result_type,
3235  __detail::_Shift<result_type, 32>::__value>(*__iter));
3236  }
3237 
3238  template<typename _InputIterator>
3239  seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3240  {
3241  for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3242  _M_v.push_back(__detail::__mod<result_type,
3243  __detail::_Shift<result_type, 32>::__value>(*__iter));
3244  }
3245 
3246  template<typename _RandomAccessIterator>
3247  void
3248  seed_seq::generate(_RandomAccessIterator __begin,
3249  _RandomAccessIterator __end)
3250  {
3251  typedef typename iterator_traits<_RandomAccessIterator>::value_type
3252  _Type;
3253 
3254  if (__begin == __end)
3255  return;
3256 
3257  std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3258 
3259  const size_t __n = __end - __begin;
3260  const size_t __s = _M_v.size();
3261  const size_t __t = (__n >= 623) ? 11
3262  : (__n >= 68) ? 7
3263  : (__n >= 39) ? 5
3264  : (__n >= 7) ? 3
3265  : (__n - 1) / 2;
3266  const size_t __p = (__n - __t) / 2;
3267  const size_t __q = __p + __t;
3268  const size_t __m = std::max(size_t(__s + 1), __n);
3269 
3270  for (size_t __k = 0; __k < __m; ++__k)
3271  {
3272  _Type __arg = (__begin[__k % __n]
3273  ^ __begin[(__k + __p) % __n]
3274  ^ __begin[(__k - 1) % __n]);
3275  _Type __r1 = __arg ^ (__arg >> 27);
3276  __r1 = __detail::__mod<_Type,
3277  __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
3278  _Type __r2 = __r1;
3279  if (__k == 0)
3280  __r2 += __s;
3281  else if (__k <= __s)
3282  __r2 += __k % __n + _M_v[__k - 1];
3283  else
3284  __r2 += __k % __n;
3285  __r2 = __detail::__mod<_Type,
3286  __detail::_Shift<_Type, 32>::__value>(__r2);
3287  __begin[(__k + __p) % __n] += __r1;
3288  __begin[(__k + __q) % __n] += __r2;
3289  __begin[__k % __n] = __r2;
3290  }
3291 
3292  for (size_t __k = __m; __k < __m + __n; ++__k)
3293  {
3294  _Type __arg = (__begin[__k % __n]
3295  + __begin[(__k + __p) % __n]
3296  + __begin[(__k - 1) % __n]);
3297  _Type __r3 = __arg ^ (__arg >> 27);
3298  __r3 = __detail::__mod<_Type,
3299  __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
3300  _Type __r4 = __r3 - __k % __n;
3301  __r4 = __detail::__mod<_Type,
3302  __detail::_Shift<_Type, 32>::__value>(__r4);
3303  __begin[(__k + __p) % __n] ^= __r3;
3304  __begin[(__k + __q) % __n] ^= __r4;
3305  __begin[__k % __n] = __r4;
3306  }
3307  }
3308 
3309  template<typename _RealType, size_t __bits,
3310  typename _UniformRandomNumberGenerator>
3311  _RealType
3312  generate_canonical(_UniformRandomNumberGenerator& __urng)
3313  {
3315  "template argument must be a floating point type");
3316 
3317  const size_t __b
3318  = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3319  __bits);
3320  const long double __r = static_cast<long double>(__urng.max())
3321  - static_cast<long double>(__urng.min()) + 1.0L;
3322  const size_t __log2r = std::log(__r) / std::log(2.0L);
3323  const size_t __m = std::max<size_t>(1UL,
3324  (__b + __log2r - 1UL) / __log2r);
3325  _RealType __ret;
3326  _RealType __sum = _RealType(0);
3327  _RealType __tmp = _RealType(1);
3328  for (size_t __k = __m; __k != 0; --__k)
3329  {
3330  __sum += _RealType(__urng() - __urng.min()) * __tmp;
3331  __tmp *= __r;
3332  }
3333  __ret = __sum / __tmp;
3334  if (__builtin_expect(__ret >= _RealType(1), 0))
3335  {
3336 #if _GLIBCXX_USE_C99_MATH_TR1
3337  __ret = std::nextafter(_RealType(1), _RealType(0));
3338 #else
3339  __ret = _RealType(1)
3340  - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3341 #endif
3342  }
3343  return __ret;
3344  }
3345 
3346 _GLIBCXX_END_NAMESPACE_VERSION
3347 } // namespace
3348 
3349 #endif
static constexpr result_type modulus
Definition: random.h:252
A piecewise_linear_distribution random number distribution.
Definition: random.h:5673
ios_base & left(ios_base &__base)
Calls base.setf(ios_base::left, ios_base::adjustfield).
Definition: ios_base.h:999
A chi_squared_distribution random number distribution.
Definition: random.h:2574
_Tp abs(const complex< _Tp > &)
Return magnitude of z.
Definition: complex:601
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:3750
param_type param() const
Returns the parameter set of the distribution.
constexpr int __lg(int __n)
This is a helper function for the sort routines and for random.tcc.
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
Definition: complex:792
A Bernoulli random number distribution.
Definition: random.h:3447
char_type widen(char __c) const
Widens characters.
Definition: basic_ios.h:449
_ForwardIterator lower_bound(_ForwardIterator __first, _ForwardIterator __last, const _Tp &__val)
Finds the first position in which val could be inserted without changing the ordering.
Definition: stl_algobase.h:984
complex< _Tp > tan(const complex< _Tp > &)
Return complex tangent of z.
Definition: complex:928
Properties of fundamental types.
Definition: limits:312
result_type operator()()
Gets the next random number in the sequence.
_Tp accumulate(_InputIterator __first, _InputIterator __last, _Tp __init)
Accumulate values in a range.
Definition: stl_numeric.h:120
const _RandomNumberEngine & base() const noexcept
Gets a const reference to the underlying generator engine object.
Definition: random.h:951
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:1783
complex< _Tp > pow(const complex< _Tp > &, int)
Return x to the y&#39;th power.
Definition: complex:987
initializer_list
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:3290
ISO C++ entities toplevel namespace is std.
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2893
static constexpr _Tp epsilon() noexcept
Definition: limits:333
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4830
streamsize precision() const
Flags access.
Definition: ios_base.h:689
A discrete geometric random number distribution.
Definition: random.h:3893
static constexpr result_type multiplier
Definition: random.h:248
iterator end() noexcept
Definition: stl_vector.h:581
A model of a linear congruential random number generator.
Definition: random.h:236
char_type fill() const
Retrieves the empty character.
Definition: basic_ios.h:370
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:2863
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:2444
Template class basic_ostream.
Definition: iosfwd:86
back_insert_iterator< _Container > back_inserter(_Container &__x)
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:2010
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:2217
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4432
Uniform continuous distribution for random numbers.
Definition: random.h:1702
_RealType result_type
Definition: random.h:2355
Produces random numbers by combining random numbers from some base engine to produce random numbers w...
Definition: random.h:1284
_GLIBCXX14_CONSTEXPR const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:219
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:3998
A negative_binomial_distribution random number distribution.
Definition: random.h:4098
static constexpr _Tp min() noexcept
Definition: limits:317
_OutputIterator partial_sum(_InputIterator __first, _InputIterator __last, _OutputIterator __result)
Return list of partial sums.
Definition: stl_numeric.h:237
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2040
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4175
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2474
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4860
A extreme_value_distribution random number distribution.
Definition: random.h:4961
iterator begin() noexcept
Definition: stl_vector.h:563
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:2637
A normal continuous distribution for random numbers.
Definition: random.h:1925
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:3968
seed_seq() noexcept
Definition: random.h:5961
result_type operator()()
Gets the next value in the generated random number sequence.
void seed(result_type __sd=default_seed)
Seeds the initial state of the random number generator.
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:5068
void push_back(const value_type &__x)
Add data to the end of the vector.
Definition: stl_vector.h:939
A cauchy_distribution random number distribution.
Definition: random.h:2789
The Marsaglia-Zaman generator.
Definition: random.h:659
_RealType generate_canonical(_UniformRandomNumberGenerator &__g)
A function template for converting the output of a (integral) uniform random number generator to a fl...
_RandomNumberEngine::result_type result_type
Definition: random.h:1287
fmtflags flags() const
Access to format flags.
Definition: ios_base.h:619
ios_base & scientific(ios_base &__base)
Calls base.setf(ios_base::scientific, ios_base::floatfield).
Definition: ios_base.h:1049
A discrete_distribution random number distribution.
Definition: random.h:5166
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:5573
ios_base & dec(ios_base &__base)
Calls base.setf(ios_base::dec, ios_base::basefield).
Definition: ios_base.h:1016
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:3072
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:5537
ptrdiff_t streamsize
Integral type for I/O operation counts and buffer sizes.
Definition: postypes.h:98
A discrete Poisson random number distribution.
Definition: random.h:4325
A fisher_f_distribution random number distribution.
Definition: random.h:2995
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:5038
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:5269
void reserve(size_type __n)
Attempt to preallocate enough memory for specified number of elements.
Definition: vector.tcc:66
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:3780
Uniform discrete distribution for random numbers. A discrete random distribution on the range with e...
A discrete binomial random number distribution.
Definition: random.h:3657
result_type operator()()
Gets the next value in the generated random number sequence.
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4617
static constexpr result_type increment
Definition: random.h:250
std::vector< double > densities() const
Returns a vector of the probability densities.
Definition: random.h:5527
A lognormal_distribution random number distribution.
Definition: random.h:2143
An exponential continuous distribution for random numbers.
Definition: random.h:4546
static constexpr _Tp max() noexcept
Definition: limits:321
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:5847
size_type size() const noexcept
Definition: stl_vector.h:670
std::vector< _RealType > intervals() const
Returns a vector of the intervals.
Definition: random.h:5511
std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, const _Quoted_string< basic_string< _CharT, _Traits, _Alloc > &, _CharT > &__str)
Extractor for delimited strings. The left and right delimiters can be different.
A piecewise_constant_distribution random number distribution.
Definition: random.h:5401
A gamma continuous distribution for random numbers.
Definition: random.h:2352
ios_base & fixed(ios_base &__base)
Calls base.setf(ios_base::fixed, ios_base::floatfield).
Definition: ios_base.h:1041
Template class basic_istream.
Definition: iosfwd:83
A weibull_distribution random number distribution.
Definition: random.h:4753
_GLIBCXX14_CONSTEXPR const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:195
ios_base & skipws(ios_base &__base)
Calls base.setf(ios_base::skipws).
Definition: ios_base.h:942
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4402
complex< _Tp > sqrt(const complex< _Tp > &)
Return complex square root of z.
Definition: complex:901
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.
Definition: complex:765
is_floating_point
Definition: type_traits:356
A student_t_distribution random number distribution.
Definition: random.h:3224
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:5811
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:5302
void seed(result_type __s=default_seed)
Reseeds the linear_congruential_engine random number generator engine sequence to the seed __s...