libstdc++
random.tcc
Go to the documentation of this file.
1 // random number generation (out of line) -*- C++ -*-
2 
3 // Copyright (C) 2009-2012 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 to
45  // avoid integer overflow.
46  //
47  // Because a and c are compile-time integral constants the compiler
48  // kindly elides any unreachable paths.
49  //
50  // Preconditions: a > 0, m > 0.
51  //
52  // XXX FIXME: as-is, only works correctly for __m % __a < __m / __a.
53  //
54  template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool>
55  struct _Mod
56  {
57  static _Tp
58  __calc(_Tp __x)
59  {
60  if (__a == 1)
61  __x %= __m;
62  else
63  {
64  static const _Tp __q = __m / __a;
65  static const _Tp __r = __m % __a;
66 
67  _Tp __t1 = __a * (__x % __q);
68  _Tp __t2 = __r * (__x / __q);
69  if (__t1 >= __t2)
70  __x = __t1 - __t2;
71  else
72  __x = __m - __t2 + __t1;
73  }
74 
75  if (__c != 0)
76  {
77  const _Tp __d = __m - __x;
78  if (__d > __c)
79  __x += __c;
80  else
81  __x = __c - __d;
82  }
83  return __x;
84  }
85  };
86 
87  // Special case for m == 0 -- use unsigned integer overflow as modulo
88  // operator.
89  template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
90  struct _Mod<_Tp, __m, __a, __c, true>
91  {
92  static _Tp
93  __calc(_Tp __x)
94  { return __a * __x + __c; }
95  };
96 
97  template<typename _InputIterator, typename _OutputIterator,
98  typename _UnaryOperation>
99  _OutputIterator
100  __transform(_InputIterator __first, _InputIterator __last,
101  _OutputIterator __result, _UnaryOperation __unary_op)
102  {
103  for (; __first != __last; ++__first, ++__result)
104  *__result = __unary_op(*__first);
105  return __result;
106  }
107 
108  _GLIBCXX_END_NAMESPACE_VERSION
109  } // namespace __detail
110 
111 _GLIBCXX_BEGIN_NAMESPACE_VERSION
112 
113  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
114  constexpr _UIntType
115  linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
116 
117  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
118  constexpr _UIntType
119  linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
120 
121  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
122  constexpr _UIntType
123  linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
124 
125  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
126  constexpr _UIntType
127  linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
128 
129  /**
130  * Seeds the LCR with integral value @p __s, adjusted so that the
131  * ring identity is never a member of the convergence set.
132  */
133  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
134  void
135  linear_congruential_engine<_UIntType, __a, __c, __m>::
136  seed(result_type __s)
137  {
138  if ((__detail::__mod<_UIntType, __m>(__c) == 0)
139  && (__detail::__mod<_UIntType, __m>(__s) == 0))
140  _M_x = 1;
141  else
142  _M_x = __detail::__mod<_UIntType, __m>(__s);
143  }
144 
145  /**
146  * Seeds the LCR engine with a value generated by @p __q.
147  */
148  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
149  template<typename _Sseq>
150  typename std::enable_if<std::is_class<_Sseq>::value>::type
151  linear_congruential_engine<_UIntType, __a, __c, __m>::
152  seed(_Sseq& __q)
153  {
154  const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
155  : std::__lg(__m);
156  const _UIntType __k = (__k0 + 31) / 32;
157  uint_least32_t __arr[__k + 3];
158  __q.generate(__arr + 0, __arr + __k + 3);
159  _UIntType __factor = 1u;
160  _UIntType __sum = 0u;
161  for (size_t __j = 0; __j < __k; ++__j)
162  {
163  __sum += __arr[__j + 3] * __factor;
164  __factor *= __detail::_Shift<_UIntType, 32>::__value;
165  }
166  seed(__sum);
167  }
168 
169  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
170  typename _CharT, typename _Traits>
171  std::basic_ostream<_CharT, _Traits>&
172  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
173  const linear_congruential_engine<_UIntType,
174  __a, __c, __m>& __lcr)
175  {
176  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
177  typedef typename __ostream_type::ios_base __ios_base;
178 
179  const typename __ios_base::fmtflags __flags = __os.flags();
180  const _CharT __fill = __os.fill();
181  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
182  __os.fill(__os.widen(' '));
183 
184  __os << __lcr._M_x;
185 
186  __os.flags(__flags);
187  __os.fill(__fill);
188  return __os;
189  }
190 
191  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
192  typename _CharT, typename _Traits>
193  std::basic_istream<_CharT, _Traits>&
194  operator>>(std::basic_istream<_CharT, _Traits>& __is,
195  linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
196  {
197  typedef std::basic_istream<_CharT, _Traits> __istream_type;
198  typedef typename __istream_type::ios_base __ios_base;
199 
200  const typename __ios_base::fmtflags __flags = __is.flags();
201  __is.flags(__ios_base::dec);
202 
203  __is >> __lcr._M_x;
204 
205  __is.flags(__flags);
206  return __is;
207  }
208 
209 
210  template<typename _UIntType,
211  size_t __w, size_t __n, size_t __m, size_t __r,
212  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
213  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
214  _UIntType __f>
215  constexpr size_t
216  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
217  __s, __b, __t, __c, __l, __f>::word_size;
218 
219  template<typename _UIntType,
220  size_t __w, size_t __n, size_t __m, size_t __r,
221  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
222  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
223  _UIntType __f>
224  constexpr size_t
225  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
226  __s, __b, __t, __c, __l, __f>::state_size;
227 
228  template<typename _UIntType,
229  size_t __w, size_t __n, size_t __m, size_t __r,
230  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
231  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
232  _UIntType __f>
233  constexpr size_t
234  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
235  __s, __b, __t, __c, __l, __f>::shift_size;
236 
237  template<typename _UIntType,
238  size_t __w, size_t __n, size_t __m, size_t __r,
239  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
240  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
241  _UIntType __f>
242  constexpr size_t
243  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
244  __s, __b, __t, __c, __l, __f>::mask_bits;
245 
246  template<typename _UIntType,
247  size_t __w, size_t __n, size_t __m, size_t __r,
248  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
249  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
250  _UIntType __f>
251  constexpr _UIntType
252  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
253  __s, __b, __t, __c, __l, __f>::xor_mask;
254 
255  template<typename _UIntType,
256  size_t __w, size_t __n, size_t __m, size_t __r,
257  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
258  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
259  _UIntType __f>
260  constexpr size_t
261  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
262  __s, __b, __t, __c, __l, __f>::tempering_u;
263 
264  template<typename _UIntType,
265  size_t __w, size_t __n, size_t __m, size_t __r,
266  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
267  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
268  _UIntType __f>
269  constexpr _UIntType
270  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
271  __s, __b, __t, __c, __l, __f>::tempering_d;
272 
273  template<typename _UIntType,
274  size_t __w, size_t __n, size_t __m, size_t __r,
275  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
276  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
277  _UIntType __f>
278  constexpr size_t
279  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
280  __s, __b, __t, __c, __l, __f>::tempering_s;
281 
282  template<typename _UIntType,
283  size_t __w, size_t __n, size_t __m, size_t __r,
284  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
285  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
286  _UIntType __f>
287  constexpr _UIntType
288  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
289  __s, __b, __t, __c, __l, __f>::tempering_b;
290 
291  template<typename _UIntType,
292  size_t __w, size_t __n, size_t __m, size_t __r,
293  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
294  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
295  _UIntType __f>
296  constexpr size_t
297  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
298  __s, __b, __t, __c, __l, __f>::tempering_t;
299 
300  template<typename _UIntType,
301  size_t __w, size_t __n, size_t __m, size_t __r,
302  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
303  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
304  _UIntType __f>
305  constexpr _UIntType
306  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
307  __s, __b, __t, __c, __l, __f>::tempering_c;
308 
309  template<typename _UIntType,
310  size_t __w, size_t __n, size_t __m, size_t __r,
311  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
312  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
313  _UIntType __f>
314  constexpr size_t
315  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
316  __s, __b, __t, __c, __l, __f>::tempering_l;
317 
318  template<typename _UIntType,
319  size_t __w, size_t __n, size_t __m, size_t __r,
320  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
321  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
322  _UIntType __f>
323  constexpr _UIntType
324  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
325  __s, __b, __t, __c, __l, __f>::
326  initialization_multiplier;
327 
328  template<typename _UIntType,
329  size_t __w, size_t __n, size_t __m, size_t __r,
330  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
331  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
332  _UIntType __f>
333  constexpr _UIntType
334  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
335  __s, __b, __t, __c, __l, __f>::default_seed;
336 
337  template<typename _UIntType,
338  size_t __w, size_t __n, size_t __m, size_t __r,
339  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
340  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
341  _UIntType __f>
342  void
343  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
344  __s, __b, __t, __c, __l, __f>::
345  seed(result_type __sd)
346  {
347  _M_x[0] = __detail::__mod<_UIntType,
348  __detail::_Shift<_UIntType, __w>::__value>(__sd);
349 
350  for (size_t __i = 1; __i < state_size; ++__i)
351  {
352  _UIntType __x = _M_x[__i - 1];
353  __x ^= __x >> (__w - 2);
354  __x *= __f;
355  __x += __detail::__mod<_UIntType, __n>(__i);
356  _M_x[__i] = __detail::__mod<_UIntType,
357  __detail::_Shift<_UIntType, __w>::__value>(__x);
358  }
359  _M_p = state_size;
360  }
361 
362  template<typename _UIntType,
363  size_t __w, size_t __n, size_t __m, size_t __r,
364  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
365  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
366  _UIntType __f>
367  template<typename _Sseq>
368  typename std::enable_if<std::is_class<_Sseq>::value>::type
369  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
370  __s, __b, __t, __c, __l, __f>::
371  seed(_Sseq& __q)
372  {
373  const _UIntType __upper_mask = (~_UIntType()) << __r;
374  const size_t __k = (__w + 31) / 32;
375  uint_least32_t __arr[__n * __k];
376  __q.generate(__arr + 0, __arr + __n * __k);
377 
378  bool __zero = true;
379  for (size_t __i = 0; __i < state_size; ++__i)
380  {
381  _UIntType __factor = 1u;
382  _UIntType __sum = 0u;
383  for (size_t __j = 0; __j < __k; ++__j)
384  {
385  __sum += __arr[__k * __i + __j] * __factor;
386  __factor *= __detail::_Shift<_UIntType, 32>::__value;
387  }
388  _M_x[__i] = __detail::__mod<_UIntType,
389  __detail::_Shift<_UIntType, __w>::__value>(__sum);
390 
391  if (__zero)
392  {
393  if (__i == 0)
394  {
395  if ((_M_x[0] & __upper_mask) != 0u)
396  __zero = false;
397  }
398  else if (_M_x[__i] != 0u)
399  __zero = false;
400  }
401  }
402  if (__zero)
403  _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
404  _M_p = state_size;
405  }
406 
407  template<typename _UIntType, size_t __w,
408  size_t __n, size_t __m, size_t __r,
409  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
410  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
411  _UIntType __f>
412  typename
413  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
414  __s, __b, __t, __c, __l, __f>::result_type
415  mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
416  __s, __b, __t, __c, __l, __f>::
417  operator()()
418  {
419  // Reload the vector - cost is O(n) amortized over n calls.
420  if (_M_p >= state_size)
421  {
422  const _UIntType __upper_mask = (~_UIntType()) << __r;
423  const _UIntType __lower_mask = ~__upper_mask;
424 
425  for (size_t __k = 0; __k < (__n - __m); ++__k)
426  {
427  _UIntType __y = ((_M_x[__k] & __upper_mask)
428  | (_M_x[__k + 1] & __lower_mask));
429  _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
430  ^ ((__y & 0x01) ? __a : 0));
431  }
432 
433  for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
434  {
435  _UIntType __y = ((_M_x[__k] & __upper_mask)
436  | (_M_x[__k + 1] & __lower_mask));
437  _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
438  ^ ((__y & 0x01) ? __a : 0));
439  }
440 
441  _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
442  | (_M_x[0] & __lower_mask));
443  _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
444  ^ ((__y & 0x01) ? __a : 0));
445  _M_p = 0;
446  }
447 
448  // Calculate o(x(i)).
449  result_type __z = _M_x[_M_p++];
450  __z ^= (__z >> __u) & __d;
451  __z ^= (__z << __s) & __b;
452  __z ^= (__z << __t) & __c;
453  __z ^= (__z >> __l);
454 
455  return __z;
456  }
457 
458  template<typename _UIntType, size_t __w,
459  size_t __n, size_t __m, size_t __r,
460  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
461  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
462  _UIntType __f, typename _CharT, typename _Traits>
463  std::basic_ostream<_CharT, _Traits>&
464  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
465  const mersenne_twister_engine<_UIntType, __w, __n, __m,
466  __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
467  {
468  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
469  typedef typename __ostream_type::ios_base __ios_base;
470 
471  const typename __ios_base::fmtflags __flags = __os.flags();
472  const _CharT __fill = __os.fill();
473  const _CharT __space = __os.widen(' ');
474  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
475  __os.fill(__space);
476 
477  for (size_t __i = 0; __i < __n; ++__i)
478  __os << __x._M_x[__i] << __space;
479  __os << __x._M_p;
480 
481  __os.flags(__flags);
482  __os.fill(__fill);
483  return __os;
484  }
485 
486  template<typename _UIntType, size_t __w,
487  size_t __n, size_t __m, size_t __r,
488  _UIntType __a, size_t __u, _UIntType __d, size_t __s,
489  _UIntType __b, size_t __t, _UIntType __c, size_t __l,
490  _UIntType __f, typename _CharT, typename _Traits>
491  std::basic_istream<_CharT, _Traits>&
492  operator>>(std::basic_istream<_CharT, _Traits>& __is,
493  mersenne_twister_engine<_UIntType, __w, __n, __m,
494  __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
495  {
496  typedef std::basic_istream<_CharT, _Traits> __istream_type;
497  typedef typename __istream_type::ios_base __ios_base;
498 
499  const typename __ios_base::fmtflags __flags = __is.flags();
500  __is.flags(__ios_base::dec | __ios_base::skipws);
501 
502  for (size_t __i = 0; __i < __n; ++__i)
503  __is >> __x._M_x[__i];
504  __is >> __x._M_p;
505 
506  __is.flags(__flags);
507  return __is;
508  }
509 
510 
511  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
512  constexpr size_t
513  subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
514 
515  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
516  constexpr size_t
517  subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
518 
519  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
520  constexpr size_t
521  subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
522 
523  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
524  constexpr _UIntType
525  subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
526 
527  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
528  void
529  subtract_with_carry_engine<_UIntType, __w, __s, __r>::
530  seed(result_type __value)
531  {
532  std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
533  __lcg(__value == 0u ? default_seed : __value);
534 
535  const size_t __n = (__w + 31) / 32;
536 
537  for (size_t __i = 0; __i < long_lag; ++__i)
538  {
539  _UIntType __sum = 0u;
540  _UIntType __factor = 1u;
541  for (size_t __j = 0; __j < __n; ++__j)
542  {
543  __sum += __detail::__mod<uint_least32_t,
544  __detail::_Shift<uint_least32_t, 32>::__value>
545  (__lcg()) * __factor;
546  __factor *= __detail::_Shift<_UIntType, 32>::__value;
547  }
548  _M_x[__i] = __detail::__mod<_UIntType,
549  __detail::_Shift<_UIntType, __w>::__value>(__sum);
550  }
551  _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
552  _M_p = 0;
553  }
554 
555  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
556  template<typename _Sseq>
557  typename std::enable_if<std::is_class<_Sseq>::value>::type
558  subtract_with_carry_engine<_UIntType, __w, __s, __r>::
559  seed(_Sseq& __q)
560  {
561  const size_t __k = (__w + 31) / 32;
562  uint_least32_t __arr[__r * __k];
563  __q.generate(__arr + 0, __arr + __r * __k);
564 
565  for (size_t __i = 0; __i < long_lag; ++__i)
566  {
567  _UIntType __sum = 0u;
568  _UIntType __factor = 1u;
569  for (size_t __j = 0; __j < __k; ++__j)
570  {
571  __sum += __arr[__k * __i + __j] * __factor;
572  __factor *= __detail::_Shift<_UIntType, 32>::__value;
573  }
574  _M_x[__i] = __detail::__mod<_UIntType,
575  __detail::_Shift<_UIntType, __w>::__value>(__sum);
576  }
577  _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
578  _M_p = 0;
579  }
580 
581  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
582  typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
583  result_type
584  subtract_with_carry_engine<_UIntType, __w, __s, __r>::
585  operator()()
586  {
587  // Derive short lag index from current index.
588  long __ps = _M_p - short_lag;
589  if (__ps < 0)
590  __ps += long_lag;
591 
592  // Calculate new x(i) without overflow or division.
593  // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
594  // cannot overflow.
595  _UIntType __xi;
596  if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
597  {
598  __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
599  _M_carry = 0;
600  }
601  else
602  {
603  __xi = (__detail::_Shift<_UIntType, __w>::__value
604  - _M_x[_M_p] - _M_carry + _M_x[__ps]);
605  _M_carry = 1;
606  }
607  _M_x[_M_p] = __xi;
608 
609  // Adjust current index to loop around in ring buffer.
610  if (++_M_p >= long_lag)
611  _M_p = 0;
612 
613  return __xi;
614  }
615 
616  template<typename _UIntType, size_t __w, size_t __s, size_t __r,
617  typename _CharT, typename _Traits>
618  std::basic_ostream<_CharT, _Traits>&
619  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
620  const subtract_with_carry_engine<_UIntType,
621  __w, __s, __r>& __x)
622  {
623  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
624  typedef typename __ostream_type::ios_base __ios_base;
625 
626  const typename __ios_base::fmtflags __flags = __os.flags();
627  const _CharT __fill = __os.fill();
628  const _CharT __space = __os.widen(' ');
629  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
630  __os.fill(__space);
631 
632  for (size_t __i = 0; __i < __r; ++__i)
633  __os << __x._M_x[__i] << __space;
634  __os << __x._M_carry << __space << __x._M_p;
635 
636  __os.flags(__flags);
637  __os.fill(__fill);
638  return __os;
639  }
640 
641  template<typename _UIntType, size_t __w, size_t __s, size_t __r,
642  typename _CharT, typename _Traits>
643  std::basic_istream<_CharT, _Traits>&
644  operator>>(std::basic_istream<_CharT, _Traits>& __is,
645  subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
646  {
647  typedef std::basic_ostream<_CharT, _Traits> __istream_type;
648  typedef typename __istream_type::ios_base __ios_base;
649 
650  const typename __ios_base::fmtflags __flags = __is.flags();
651  __is.flags(__ios_base::dec | __ios_base::skipws);
652 
653  for (size_t __i = 0; __i < __r; ++__i)
654  __is >> __x._M_x[__i];
655  __is >> __x._M_carry;
656  __is >> __x._M_p;
657 
658  __is.flags(__flags);
659  return __is;
660  }
661 
662 
663  template<typename _RandomNumberEngine, size_t __p, size_t __r>
664  constexpr size_t
665  discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
666 
667  template<typename _RandomNumberEngine, size_t __p, size_t __r>
668  constexpr size_t
669  discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
670 
671  template<typename _RandomNumberEngine, size_t __p, size_t __r>
672  typename discard_block_engine<_RandomNumberEngine,
673  __p, __r>::result_type
674  discard_block_engine<_RandomNumberEngine, __p, __r>::
675  operator()()
676  {
677  if (_M_n >= used_block)
678  {
679  _M_b.discard(block_size - _M_n);
680  _M_n = 0;
681  }
682  ++_M_n;
683  return _M_b();
684  }
685 
686  template<typename _RandomNumberEngine, size_t __p, size_t __r,
687  typename _CharT, typename _Traits>
688  std::basic_ostream<_CharT, _Traits>&
689  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
690  const discard_block_engine<_RandomNumberEngine,
691  __p, __r>& __x)
692  {
693  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
694  typedef typename __ostream_type::ios_base __ios_base;
695 
696  const typename __ios_base::fmtflags __flags = __os.flags();
697  const _CharT __fill = __os.fill();
698  const _CharT __space = __os.widen(' ');
699  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
700  __os.fill(__space);
701 
702  __os << __x.base() << __space << __x._M_n;
703 
704  __os.flags(__flags);
705  __os.fill(__fill);
706  return __os;
707  }
708 
709  template<typename _RandomNumberEngine, size_t __p, size_t __r,
710  typename _CharT, typename _Traits>
711  std::basic_istream<_CharT, _Traits>&
712  operator>>(std::basic_istream<_CharT, _Traits>& __is,
713  discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
714  {
715  typedef std::basic_istream<_CharT, _Traits> __istream_type;
716  typedef typename __istream_type::ios_base __ios_base;
717 
718  const typename __ios_base::fmtflags __flags = __is.flags();
719  __is.flags(__ios_base::dec | __ios_base::skipws);
720 
721  __is >> __x._M_b >> __x._M_n;
722 
723  __is.flags(__flags);
724  return __is;
725  }
726 
727 
728  template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
729  typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
730  result_type
731  independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
732  operator()()
733  {
734  typedef typename _RandomNumberEngine::result_type _Eresult_type;
735  const _Eresult_type __r
736  = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
737  ? _M_b.max() - _M_b.min() + 1 : 0);
738  const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
739  const unsigned __m = __r ? std::__lg(__r) : __edig;
740 
741  typedef typename std::common_type<_Eresult_type, result_type>::type
742  __ctype;
743  const unsigned __cdig = std::numeric_limits<__ctype>::digits;
744 
745  unsigned __n, __n0;
746  __ctype __s0, __s1, __y0, __y1;
747 
748  for (size_t __i = 0; __i < 2; ++__i)
749  {
750  __n = (__w + __m - 1) / __m + __i;
751  __n0 = __n - __w % __n;
752  const unsigned __w0 = __w / __n; // __w0 <= __m
753 
754  __s0 = 0;
755  __s1 = 0;
756  if (__w0 < __cdig)
757  {
758  __s0 = __ctype(1) << __w0;
759  __s1 = __s0 << 1;
760  }
761 
762  __y0 = 0;
763  __y1 = 0;
764  if (__r)
765  {
766  __y0 = __s0 * (__r / __s0);
767  if (__s1)
768  __y1 = __s1 * (__r / __s1);
769 
770  if (__r - __y0 <= __y0 / __n)
771  break;
772  }
773  else
774  break;
775  }
776 
777  result_type __sum = 0;
778  for (size_t __k = 0; __k < __n0; ++__k)
779  {
780  __ctype __u;
781  do
782  __u = _M_b() - _M_b.min();
783  while (__y0 && __u >= __y0);
784  __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
785  }
786  for (size_t __k = __n0; __k < __n; ++__k)
787  {
788  __ctype __u;
789  do
790  __u = _M_b() - _M_b.min();
791  while (__y1 && __u >= __y1);
792  __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
793  }
794  return __sum;
795  }
796 
797 
798  template<typename _RandomNumberEngine, size_t __k>
799  constexpr size_t
800  shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
801 
802  template<typename _RandomNumberEngine, size_t __k>
803  typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
804  shuffle_order_engine<_RandomNumberEngine, __k>::
805  operator()()
806  {
807  size_t __j = __k * ((_M_y - _M_b.min())
808  / (_M_b.max() - _M_b.min() + 1.0L));
809  _M_y = _M_v[__j];
810  _M_v[__j] = _M_b();
811 
812  return _M_y;
813  }
814 
815  template<typename _RandomNumberEngine, size_t __k,
816  typename _CharT, typename _Traits>
817  std::basic_ostream<_CharT, _Traits>&
818  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
819  const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
820  {
821  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
822  typedef typename __ostream_type::ios_base __ios_base;
823 
824  const typename __ios_base::fmtflags __flags = __os.flags();
825  const _CharT __fill = __os.fill();
826  const _CharT __space = __os.widen(' ');
827  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
828  __os.fill(__space);
829 
830  __os << __x.base();
831  for (size_t __i = 0; __i < __k; ++__i)
832  __os << __space << __x._M_v[__i];
833  __os << __space << __x._M_y;
834 
835  __os.flags(__flags);
836  __os.fill(__fill);
837  return __os;
838  }
839 
840  template<typename _RandomNumberEngine, size_t __k,
841  typename _CharT, typename _Traits>
842  std::basic_istream<_CharT, _Traits>&
843  operator>>(std::basic_istream<_CharT, _Traits>& __is,
844  shuffle_order_engine<_RandomNumberEngine, __k>& __x)
845  {
846  typedef std::basic_istream<_CharT, _Traits> __istream_type;
847  typedef typename __istream_type::ios_base __ios_base;
848 
849  const typename __ios_base::fmtflags __flags = __is.flags();
850  __is.flags(__ios_base::dec | __ios_base::skipws);
851 
852  __is >> __x._M_b;
853  for (size_t __i = 0; __i < __k; ++__i)
854  __is >> __x._M_v[__i];
855  __is >> __x._M_y;
856 
857  __is.flags(__flags);
858  return __is;
859  }
860 
861 
862  template<typename _IntType>
863  template<typename _UniformRandomNumberGenerator>
864  typename uniform_int_distribution<_IntType>::result_type
865  uniform_int_distribution<_IntType>::
866  operator()(_UniformRandomNumberGenerator& __urng,
867  const param_type& __param)
868  {
869  typedef typename _UniformRandomNumberGenerator::result_type
870  _Gresult_type;
871  typedef typename std::make_unsigned<result_type>::type __utype;
872  typedef typename std::common_type<_Gresult_type, __utype>::type
873  __uctype;
874 
875  const __uctype __urngmin = __urng.min();
876  const __uctype __urngmax = __urng.max();
877  const __uctype __urngrange = __urngmax - __urngmin;
878  const __uctype __urange
879  = __uctype(__param.b()) - __uctype(__param.a());
880 
881  __uctype __ret;
882 
883  if (__urngrange > __urange)
884  {
885  // downscaling
886  const __uctype __uerange = __urange + 1; // __urange can be zero
887  const __uctype __scaling = __urngrange / __uerange;
888  const __uctype __past = __uerange * __scaling;
889  do
890  __ret = __uctype(__urng()) - __urngmin;
891  while (__ret >= __past);
892  __ret /= __scaling;
893  }
894  else if (__urngrange < __urange)
895  {
896  // upscaling
897  /*
898  Note that every value in [0, urange]
899  can be written uniquely as
900 
901  (urngrange + 1) * high + low
902 
903  where
904 
905  high in [0, urange / (urngrange + 1)]
906 
907  and
908 
909  low in [0, urngrange].
910  */
911  __uctype __tmp; // wraparound control
912  do
913  {
914  const __uctype __uerngrange = __urngrange + 1;
915  __tmp = (__uerngrange * operator()
916  (__urng, param_type(0, __urange / __uerngrange)));
917  __ret = __tmp + (__uctype(__urng()) - __urngmin);
918  }
919  while (__ret > __urange || __ret < __tmp);
920  }
921  else
922  __ret = __uctype(__urng()) - __urngmin;
923 
924  return __ret + __param.a();
925  }
926 
927  template<typename _IntType, typename _CharT, typename _Traits>
928  std::basic_ostream<_CharT, _Traits>&
929  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
930  const uniform_int_distribution<_IntType>& __x)
931  {
932  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
933  typedef typename __ostream_type::ios_base __ios_base;
934 
935  const typename __ios_base::fmtflags __flags = __os.flags();
936  const _CharT __fill = __os.fill();
937  const _CharT __space = __os.widen(' ');
938  __os.flags(__ios_base::scientific | __ios_base::left);
939  __os.fill(__space);
940 
941  __os << __x.a() << __space << __x.b();
942 
943  __os.flags(__flags);
944  __os.fill(__fill);
945  return __os;
946  }
947 
948  template<typename _IntType, typename _CharT, typename _Traits>
949  std::basic_istream<_CharT, _Traits>&
950  operator>>(std::basic_istream<_CharT, _Traits>& __is,
951  uniform_int_distribution<_IntType>& __x)
952  {
953  typedef std::basic_istream<_CharT, _Traits> __istream_type;
954  typedef typename __istream_type::ios_base __ios_base;
955 
956  const typename __ios_base::fmtflags __flags = __is.flags();
957  __is.flags(__ios_base::dec | __ios_base::skipws);
958 
959  _IntType __a, __b;
960  __is >> __a >> __b;
961  __x.param(typename uniform_int_distribution<_IntType>::
962  param_type(__a, __b));
963 
964  __is.flags(__flags);
965  return __is;
966  }
967 
968 
969  template<typename _RealType, typename _CharT, typename _Traits>
970  std::basic_ostream<_CharT, _Traits>&
971  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
972  const uniform_real_distribution<_RealType>& __x)
973  {
974  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
975  typedef typename __ostream_type::ios_base __ios_base;
976 
977  const typename __ios_base::fmtflags __flags = __os.flags();
978  const _CharT __fill = __os.fill();
979  const std::streamsize __precision = __os.precision();
980  const _CharT __space = __os.widen(' ');
981  __os.flags(__ios_base::scientific | __ios_base::left);
982  __os.fill(__space);
983  __os.precision(std::numeric_limits<_RealType>::max_digits10);
984 
985  __os << __x.a() << __space << __x.b();
986 
987  __os.flags(__flags);
988  __os.fill(__fill);
989  __os.precision(__precision);
990  return __os;
991  }
992 
993  template<typename _RealType, typename _CharT, typename _Traits>
994  std::basic_istream<_CharT, _Traits>&
995  operator>>(std::basic_istream<_CharT, _Traits>& __is,
996  uniform_real_distribution<_RealType>& __x)
997  {
998  typedef std::basic_istream<_CharT, _Traits> __istream_type;
999  typedef typename __istream_type::ios_base __ios_base;
1000 
1001  const typename __ios_base::fmtflags __flags = __is.flags();
1002  __is.flags(__ios_base::skipws);
1003 
1004  _RealType __a, __b;
1005  __is >> __a >> __b;
1006  __x.param(typename uniform_real_distribution<_RealType>::
1007  param_type(__a, __b));
1008 
1009  __is.flags(__flags);
1010  return __is;
1011  }
1012 
1013 
1014  template<typename _CharT, typename _Traits>
1015  std::basic_ostream<_CharT, _Traits>&
1016  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1017  const bernoulli_distribution& __x)
1018  {
1019  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1020  typedef typename __ostream_type::ios_base __ios_base;
1021 
1022  const typename __ios_base::fmtflags __flags = __os.flags();
1023  const _CharT __fill = __os.fill();
1024  const std::streamsize __precision = __os.precision();
1025  __os.flags(__ios_base::scientific | __ios_base::left);
1026  __os.fill(__os.widen(' '));
1027  __os.precision(std::numeric_limits<double>::max_digits10);
1028 
1029  __os << __x.p();
1030 
1031  __os.flags(__flags);
1032  __os.fill(__fill);
1033  __os.precision(__precision);
1034  return __os;
1035  }
1036 
1037 
1038  template<typename _IntType>
1039  template<typename _UniformRandomNumberGenerator>
1040  typename geometric_distribution<_IntType>::result_type
1041  geometric_distribution<_IntType>::
1042  operator()(_UniformRandomNumberGenerator& __urng,
1043  const param_type& __param)
1044  {
1045  // About the epsilon thing see this thread:
1046  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1047  const double __naf =
1048  (1 - std::numeric_limits<double>::epsilon()) / 2;
1049  // The largest _RealType convertible to _IntType.
1050  const double __thr =
1051  std::numeric_limits<_IntType>::max() + __naf;
1052  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1053  __aurng(__urng);
1054 
1055  double __cand;
1056  do
1057  __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1058  while (__cand >= __thr);
1059 
1060  return result_type(__cand + __naf);
1061  }
1062 
1063  template<typename _IntType,
1064  typename _CharT, typename _Traits>
1065  std::basic_ostream<_CharT, _Traits>&
1066  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1067  const geometric_distribution<_IntType>& __x)
1068  {
1069  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1070  typedef typename __ostream_type::ios_base __ios_base;
1071 
1072  const typename __ios_base::fmtflags __flags = __os.flags();
1073  const _CharT __fill = __os.fill();
1074  const std::streamsize __precision = __os.precision();
1075  __os.flags(__ios_base::scientific | __ios_base::left);
1076  __os.fill(__os.widen(' '));
1077  __os.precision(std::numeric_limits<double>::max_digits10);
1078 
1079  __os << __x.p();
1080 
1081  __os.flags(__flags);
1082  __os.fill(__fill);
1083  __os.precision(__precision);
1084  return __os;
1085  }
1086 
1087  template<typename _IntType,
1088  typename _CharT, typename _Traits>
1089  std::basic_istream<_CharT, _Traits>&
1090  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1091  geometric_distribution<_IntType>& __x)
1092  {
1093  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1094  typedef typename __istream_type::ios_base __ios_base;
1095 
1096  const typename __ios_base::fmtflags __flags = __is.flags();
1097  __is.flags(__ios_base::skipws);
1098 
1099  double __p;
1100  __is >> __p;
1101  __x.param(typename geometric_distribution<_IntType>::param_type(__p));
1102 
1103  __is.flags(__flags);
1104  return __is;
1105  }
1106 
1107  // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1108  template<typename _IntType>
1109  template<typename _UniformRandomNumberGenerator>
1110  typename negative_binomial_distribution<_IntType>::result_type
1111  negative_binomial_distribution<_IntType>::
1112  operator()(_UniformRandomNumberGenerator& __urng)
1113  {
1114  const double __y = _M_gd(__urng);
1115 
1116  // XXX Is the constructor too slow?
1117  std::poisson_distribution<result_type> __poisson(__y);
1118  return __poisson(__urng);
1119  }
1120 
1121  template<typename _IntType>
1122  template<typename _UniformRandomNumberGenerator>
1123  typename negative_binomial_distribution<_IntType>::result_type
1124  negative_binomial_distribution<_IntType>::
1125  operator()(_UniformRandomNumberGenerator& __urng,
1126  const param_type& __p)
1127  {
1128  typedef typename std::gamma_distribution<double>::param_type
1129  param_type;
1130 
1131  const double __y =
1132  _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1133 
1134  std::poisson_distribution<result_type> __poisson(__y);
1135  return __poisson(__urng);
1136  }
1137 
1138  template<typename _IntType, typename _CharT, typename _Traits>
1139  std::basic_ostream<_CharT, _Traits>&
1140  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1141  const negative_binomial_distribution<_IntType>& __x)
1142  {
1143  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1144  typedef typename __ostream_type::ios_base __ios_base;
1145 
1146  const typename __ios_base::fmtflags __flags = __os.flags();
1147  const _CharT __fill = __os.fill();
1148  const std::streamsize __precision = __os.precision();
1149  const _CharT __space = __os.widen(' ');
1150  __os.flags(__ios_base::scientific | __ios_base::left);
1151  __os.fill(__os.widen(' '));
1152  __os.precision(std::numeric_limits<double>::max_digits10);
1153 
1154  __os << __x.k() << __space << __x.p()
1155  << __space << __x._M_gd;
1156 
1157  __os.flags(__flags);
1158  __os.fill(__fill);
1159  __os.precision(__precision);
1160  return __os;
1161  }
1162 
1163  template<typename _IntType, typename _CharT, typename _Traits>
1164  std::basic_istream<_CharT, _Traits>&
1165  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1166  negative_binomial_distribution<_IntType>& __x)
1167  {
1168  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1169  typedef typename __istream_type::ios_base __ios_base;
1170 
1171  const typename __ios_base::fmtflags __flags = __is.flags();
1172  __is.flags(__ios_base::skipws);
1173 
1174  _IntType __k;
1175  double __p;
1176  __is >> __k >> __p >> __x._M_gd;
1177  __x.param(typename negative_binomial_distribution<_IntType>::
1178  param_type(__k, __p));
1179 
1180  __is.flags(__flags);
1181  return __is;
1182  }
1183 
1184 
1185  template<typename _IntType>
1186  void
1187  poisson_distribution<_IntType>::param_type::
1188  _M_initialize()
1189  {
1190 #if _GLIBCXX_USE_C99_MATH_TR1
1191  if (_M_mean >= 12)
1192  {
1193  const double __m = std::floor(_M_mean);
1194  _M_lm_thr = std::log(_M_mean);
1195  _M_lfm = std::lgamma(__m + 1);
1196  _M_sm = std::sqrt(__m);
1197 
1198  const double __pi_4 = 0.7853981633974483096156608458198757L;
1199  const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1200  / __pi_4));
1201  _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
1202  const double __cx = 2 * __m + _M_d;
1203  _M_scx = std::sqrt(__cx / 2);
1204  _M_1cx = 1 / __cx;
1205 
1206  _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1207  _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1208  / _M_d;
1209  }
1210  else
1211 #endif
1212  _M_lm_thr = std::exp(-_M_mean);
1213  }
1214 
1215  /**
1216  * A rejection algorithm when mean >= 12 and a simple method based
1217  * upon the multiplication of uniform random variates otherwise.
1218  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1219  * is defined.
1220  *
1221  * Reference:
1222  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1223  * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1224  */
1225  template<typename _IntType>
1226  template<typename _UniformRandomNumberGenerator>
1227  typename poisson_distribution<_IntType>::result_type
1228  poisson_distribution<_IntType>::
1229  operator()(_UniformRandomNumberGenerator& __urng,
1230  const param_type& __param)
1231  {
1232  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1233  __aurng(__urng);
1234 #if _GLIBCXX_USE_C99_MATH_TR1
1235  if (__param.mean() >= 12)
1236  {
1237  double __x;
1238 
1239  // See comments above...
1240  const double __naf =
1241  (1 - std::numeric_limits<double>::epsilon()) / 2;
1242  const double __thr =
1243  std::numeric_limits<_IntType>::max() + __naf;
1244 
1245  const double __m = std::floor(__param.mean());
1246  // sqrt(pi / 2)
1247  const double __spi_2 = 1.2533141373155002512078826424055226L;
1248  const double __c1 = __param._M_sm * __spi_2;
1249  const double __c2 = __param._M_c2b + __c1;
1250  const double __c3 = __c2 + 1;
1251  const double __c4 = __c3 + 1;
1252  // e^(1 / 78)
1253  const double __e178 = 1.0129030479320018583185514777512983L;
1254  const double __c5 = __c4 + __e178;
1255  const double __c = __param._M_cb + __c5;
1256  const double __2cx = 2 * (2 * __m + __param._M_d);
1257 
1258  bool __reject = true;
1259  do
1260  {
1261  const double __u = __c * __aurng();
1262  const double __e = -std::log(1.0 - __aurng());
1263 
1264  double __w = 0.0;
1265 
1266  if (__u <= __c1)
1267  {
1268  const double __n = _M_nd(__urng);
1269  const double __y = -std::abs(__n) * __param._M_sm - 1;
1270  __x = std::floor(__y);
1271  __w = -__n * __n / 2;
1272  if (__x < -__m)
1273  continue;
1274  }
1275  else if (__u <= __c2)
1276  {
1277  const double __n = _M_nd(__urng);
1278  const double __y = 1 + std::abs(__n) * __param._M_scx;
1279  __x = std::ceil(__y);
1280  __w = __y * (2 - __y) * __param._M_1cx;
1281  if (__x > __param._M_d)
1282  continue;
1283  }
1284  else if (__u <= __c3)
1285  // NB: This case not in the book, nor in the Errata,
1286  // but should be ok...
1287  __x = -1;
1288  else if (__u <= __c4)
1289  __x = 0;
1290  else if (__u <= __c5)
1291  __x = 1;
1292  else
1293  {
1294  const double __v = -std::log(1.0 - __aurng());
1295  const double __y = __param._M_d
1296  + __v * __2cx / __param._M_d;
1297  __x = std::ceil(__y);
1298  __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1299  }
1300 
1301  __reject = (__w - __e - __x * __param._M_lm_thr
1302  > __param._M_lfm - std::lgamma(__x + __m + 1));
1303 
1304  __reject |= __x + __m >= __thr;
1305 
1306  } while (__reject);
1307 
1308  return result_type(__x + __m + __naf);
1309  }
1310  else
1311 #endif
1312  {
1313  _IntType __x = 0;
1314  double __prod = 1.0;
1315 
1316  do
1317  {
1318  __prod *= __aurng();
1319  __x += 1;
1320  }
1321  while (__prod > __param._M_lm_thr);
1322 
1323  return __x - 1;
1324  }
1325  }
1326 
1327  template<typename _IntType,
1328  typename _CharT, typename _Traits>
1329  std::basic_ostream<_CharT, _Traits>&
1330  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1331  const poisson_distribution<_IntType>& __x)
1332  {
1333  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1334  typedef typename __ostream_type::ios_base __ios_base;
1335 
1336  const typename __ios_base::fmtflags __flags = __os.flags();
1337  const _CharT __fill = __os.fill();
1338  const std::streamsize __precision = __os.precision();
1339  const _CharT __space = __os.widen(' ');
1340  __os.flags(__ios_base::scientific | __ios_base::left);
1341  __os.fill(__space);
1342  __os.precision(std::numeric_limits<double>::max_digits10);
1343 
1344  __os << __x.mean() << __space << __x._M_nd;
1345 
1346  __os.flags(__flags);
1347  __os.fill(__fill);
1348  __os.precision(__precision);
1349  return __os;
1350  }
1351 
1352  template<typename _IntType,
1353  typename _CharT, typename _Traits>
1354  std::basic_istream<_CharT, _Traits>&
1355  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1356  poisson_distribution<_IntType>& __x)
1357  {
1358  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1359  typedef typename __istream_type::ios_base __ios_base;
1360 
1361  const typename __ios_base::fmtflags __flags = __is.flags();
1362  __is.flags(__ios_base::skipws);
1363 
1364  double __mean;
1365  __is >> __mean >> __x._M_nd;
1366  __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1367 
1368  __is.flags(__flags);
1369  return __is;
1370  }
1371 
1372 
1373  template<typename _IntType>
1374  void
1375  binomial_distribution<_IntType>::param_type::
1376  _M_initialize()
1377  {
1378  const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1379 
1380  _M_easy = true;
1381 
1382 #if _GLIBCXX_USE_C99_MATH_TR1
1383  if (_M_t * __p12 >= 8)
1384  {
1385  _M_easy = false;
1386  const double __np = std::floor(_M_t * __p12);
1387  const double __pa = __np / _M_t;
1388  const double __1p = 1 - __pa;
1389 
1390  const double __pi_4 = 0.7853981633974483096156608458198757L;
1391  const double __d1x =
1392  std::sqrt(__np * __1p * std::log(32 * __np
1393  / (81 * __pi_4 * __1p)));
1394  _M_d1 = std::round(std::max(1.0, __d1x));
1395  const double __d2x =
1396  std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1397  / (__pi_4 * __pa)));
1398  _M_d2 = std::round(std::max(1.0, __d2x));
1399 
1400  // sqrt(pi / 2)
1401  const double __spi_2 = 1.2533141373155002512078826424055226L;
1402  _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1403  _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1404  _M_c = 2 * _M_d1 / __np;
1405  _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1406  const double __a12 = _M_a1 + _M_s2 * __spi_2;
1407  const double __s1s = _M_s1 * _M_s1;
1408  _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1409  * 2 * __s1s / _M_d1
1410  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1411  const double __s2s = _M_s2 * _M_s2;
1412  _M_s = (_M_a123 + 2 * __s2s / _M_d2
1413  * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1414  _M_lf = (std::lgamma(__np + 1)
1415  + std::lgamma(_M_t - __np + 1));
1416  _M_lp1p = std::log(__pa / __1p);
1417 
1418  _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1419  }
1420  else
1421 #endif
1422  _M_q = -std::log(1 - __p12);
1423  }
1424 
1425  template<typename _IntType>
1426  template<typename _UniformRandomNumberGenerator>
1427  typename binomial_distribution<_IntType>::result_type
1428  binomial_distribution<_IntType>::
1429  _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1430  {
1431  _IntType __x = 0;
1432  double __sum = 0.0;
1433  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1434  __aurng(__urng);
1435 
1436  do
1437  {
1438  const double __e = -std::log(1.0 - __aurng());
1439  __sum += __e / (__t - __x);
1440  __x += 1;
1441  }
1442  while (__sum <= _M_param._M_q);
1443 
1444  return __x - 1;
1445  }
1446 
1447  /**
1448  * A rejection algorithm when t * p >= 8 and a simple waiting time
1449  * method - the second in the referenced book - otherwise.
1450  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1451  * is defined.
1452  *
1453  * Reference:
1454  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1455  * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1456  */
1457  template<typename _IntType>
1458  template<typename _UniformRandomNumberGenerator>
1459  typename binomial_distribution<_IntType>::result_type
1460  binomial_distribution<_IntType>::
1461  operator()(_UniformRandomNumberGenerator& __urng,
1462  const param_type& __param)
1463  {
1464  result_type __ret;
1465  const _IntType __t = __param.t();
1466  const double __p = __param.p();
1467  const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1468  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1469  __aurng(__urng);
1470 
1471 #if _GLIBCXX_USE_C99_MATH_TR1
1472  if (!__param._M_easy)
1473  {
1474  double __x;
1475 
1476  // See comments above...
1477  const double __naf =
1478  (1 - std::numeric_limits<double>::epsilon()) / 2;
1479  const double __thr =
1480  std::numeric_limits<_IntType>::max() + __naf;
1481 
1482  const double __np = std::floor(__t * __p12);
1483 
1484  // sqrt(pi / 2)
1485  const double __spi_2 = 1.2533141373155002512078826424055226L;
1486  const double __a1 = __param._M_a1;
1487  const double __a12 = __a1 + __param._M_s2 * __spi_2;
1488  const double __a123 = __param._M_a123;
1489  const double __s1s = __param._M_s1 * __param._M_s1;
1490  const double __s2s = __param._M_s2 * __param._M_s2;
1491 
1492  bool __reject;
1493  do
1494  {
1495  const double __u = __param._M_s * __aurng();
1496 
1497  double __v;
1498 
1499  if (__u <= __a1)
1500  {
1501  const double __n = _M_nd(__urng);
1502  const double __y = __param._M_s1 * std::abs(__n);
1503  __reject = __y >= __param._M_d1;
1504  if (!__reject)
1505  {
1506  const double __e = -std::log(1.0 - __aurng());
1507  __x = std::floor(__y);
1508  __v = -__e - __n * __n / 2 + __param._M_c;
1509  }
1510  }
1511  else if (__u <= __a12)
1512  {
1513  const double __n = _M_nd(__urng);
1514  const double __y = __param._M_s2 * std::abs(__n);
1515  __reject = __y >= __param._M_d2;
1516  if (!__reject)
1517  {
1518  const double __e = -std::log(1.0 - __aurng());
1519  __x = std::floor(-__y);
1520  __v = -__e - __n * __n / 2;
1521  }
1522  }
1523  else if (__u <= __a123)
1524  {
1525  const double __e1 = -std::log(1.0 - __aurng());
1526  const double __e2 = -std::log(1.0 - __aurng());
1527 
1528  const double __y = __param._M_d1
1529  + 2 * __s1s * __e1 / __param._M_d1;
1530  __x = std::floor(__y);
1531  __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1532  -__y / (2 * __s1s)));
1533  __reject = false;
1534  }
1535  else
1536  {
1537  const double __e1 = -std::log(1.0 - __aurng());
1538  const double __e2 = -std::log(1.0 - __aurng());
1539 
1540  const double __y = __param._M_d2
1541  + 2 * __s2s * __e1 / __param._M_d2;
1542  __x = std::floor(-__y);
1543  __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1544  __reject = false;
1545  }
1546 
1547  __reject = __reject || __x < -__np || __x > __t - __np;
1548  if (!__reject)
1549  {
1550  const double __lfx =
1551  std::lgamma(__np + __x + 1)
1552  + std::lgamma(__t - (__np + __x) + 1);
1553  __reject = __v > __param._M_lf - __lfx
1554  + __x * __param._M_lp1p;
1555  }
1556 
1557  __reject |= __x + __np >= __thr;
1558  }
1559  while (__reject);
1560 
1561  __x += __np + __naf;
1562 
1563  const _IntType __z = _M_waiting(__urng, __t - _IntType(__x));
1564  __ret = _IntType(__x) + __z;
1565  }
1566  else
1567 #endif
1568  __ret = _M_waiting(__urng, __t);
1569 
1570  if (__p12 != __p)
1571  __ret = __t - __ret;
1572  return __ret;
1573  }
1574 
1575  template<typename _IntType,
1576  typename _CharT, typename _Traits>
1577  std::basic_ostream<_CharT, _Traits>&
1578  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1579  const binomial_distribution<_IntType>& __x)
1580  {
1581  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1582  typedef typename __ostream_type::ios_base __ios_base;
1583 
1584  const typename __ios_base::fmtflags __flags = __os.flags();
1585  const _CharT __fill = __os.fill();
1586  const std::streamsize __precision = __os.precision();
1587  const _CharT __space = __os.widen(' ');
1588  __os.flags(__ios_base::scientific | __ios_base::left);
1589  __os.fill(__space);
1590  __os.precision(std::numeric_limits<double>::max_digits10);
1591 
1592  __os << __x.t() << __space << __x.p()
1593  << __space << __x._M_nd;
1594 
1595  __os.flags(__flags);
1596  __os.fill(__fill);
1597  __os.precision(__precision);
1598  return __os;
1599  }
1600 
1601  template<typename _IntType,
1602  typename _CharT, typename _Traits>
1603  std::basic_istream<_CharT, _Traits>&
1604  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1605  binomial_distribution<_IntType>& __x)
1606  {
1607  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1608  typedef typename __istream_type::ios_base __ios_base;
1609 
1610  const typename __ios_base::fmtflags __flags = __is.flags();
1611  __is.flags(__ios_base::dec | __ios_base::skipws);
1612 
1613  _IntType __t;
1614  double __p;
1615  __is >> __t >> __p >> __x._M_nd;
1616  __x.param(typename binomial_distribution<_IntType>::
1617  param_type(__t, __p));
1618 
1619  __is.flags(__flags);
1620  return __is;
1621  }
1622 
1623 
1624  template<typename _RealType, typename _CharT, typename _Traits>
1625  std::basic_ostream<_CharT, _Traits>&
1626  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1627  const exponential_distribution<_RealType>& __x)
1628  {
1629  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1630  typedef typename __ostream_type::ios_base __ios_base;
1631 
1632  const typename __ios_base::fmtflags __flags = __os.flags();
1633  const _CharT __fill = __os.fill();
1634  const std::streamsize __precision = __os.precision();
1635  __os.flags(__ios_base::scientific | __ios_base::left);
1636  __os.fill(__os.widen(' '));
1637  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1638 
1639  __os << __x.lambda();
1640 
1641  __os.flags(__flags);
1642  __os.fill(__fill);
1643  __os.precision(__precision);
1644  return __os;
1645  }
1646 
1647  template<typename _RealType, typename _CharT, typename _Traits>
1648  std::basic_istream<_CharT, _Traits>&
1649  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1650  exponential_distribution<_RealType>& __x)
1651  {
1652  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1653  typedef typename __istream_type::ios_base __ios_base;
1654 
1655  const typename __ios_base::fmtflags __flags = __is.flags();
1656  __is.flags(__ios_base::dec | __ios_base::skipws);
1657 
1658  _RealType __lambda;
1659  __is >> __lambda;
1660  __x.param(typename exponential_distribution<_RealType>::
1661  param_type(__lambda));
1662 
1663  __is.flags(__flags);
1664  return __is;
1665  }
1666 
1667 
1668  /**
1669  * Polar method due to Marsaglia.
1670  *
1671  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1672  * New York, 1986, Ch. V, Sect. 4.4.
1673  */
1674  template<typename _RealType>
1675  template<typename _UniformRandomNumberGenerator>
1676  typename normal_distribution<_RealType>::result_type
1677  normal_distribution<_RealType>::
1678  operator()(_UniformRandomNumberGenerator& __urng,
1679  const param_type& __param)
1680  {
1681  result_type __ret;
1682  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1683  __aurng(__urng);
1684 
1685  if (_M_saved_available)
1686  {
1687  _M_saved_available = false;
1688  __ret = _M_saved;
1689  }
1690  else
1691  {
1692  result_type __x, __y, __r2;
1693  do
1694  {
1695  __x = result_type(2.0) * __aurng() - 1.0;
1696  __y = result_type(2.0) * __aurng() - 1.0;
1697  __r2 = __x * __x + __y * __y;
1698  }
1699  while (__r2 > 1.0 || __r2 == 0.0);
1700 
1701  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1702  _M_saved = __x * __mult;
1703  _M_saved_available = true;
1704  __ret = __y * __mult;
1705  }
1706 
1707  __ret = __ret * __param.stddev() + __param.mean();
1708  return __ret;
1709  }
1710 
1711  template<typename _RealType>
1712  bool
1713  operator==(const std::normal_distribution<_RealType>& __d1,
1714  const std::normal_distribution<_RealType>& __d2)
1715  {
1716  if (__d1._M_param == __d2._M_param
1717  && __d1._M_saved_available == __d2._M_saved_available)
1718  {
1719  if (__d1._M_saved_available
1720  && __d1._M_saved == __d2._M_saved)
1721  return true;
1722  else if(!__d1._M_saved_available)
1723  return true;
1724  else
1725  return false;
1726  }
1727  else
1728  return false;
1729  }
1730 
1731  template<typename _RealType, typename _CharT, typename _Traits>
1732  std::basic_ostream<_CharT, _Traits>&
1733  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1734  const normal_distribution<_RealType>& __x)
1735  {
1736  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1737  typedef typename __ostream_type::ios_base __ios_base;
1738 
1739  const typename __ios_base::fmtflags __flags = __os.flags();
1740  const _CharT __fill = __os.fill();
1741  const std::streamsize __precision = __os.precision();
1742  const _CharT __space = __os.widen(' ');
1743  __os.flags(__ios_base::scientific | __ios_base::left);
1744  __os.fill(__space);
1745  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1746 
1747  __os << __x.mean() << __space << __x.stddev()
1748  << __space << __x._M_saved_available;
1749  if (__x._M_saved_available)
1750  __os << __space << __x._M_saved;
1751 
1752  __os.flags(__flags);
1753  __os.fill(__fill);
1754  __os.precision(__precision);
1755  return __os;
1756  }
1757 
1758  template<typename _RealType, typename _CharT, typename _Traits>
1759  std::basic_istream<_CharT, _Traits>&
1760  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1761  normal_distribution<_RealType>& __x)
1762  {
1763  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1764  typedef typename __istream_type::ios_base __ios_base;
1765 
1766  const typename __ios_base::fmtflags __flags = __is.flags();
1767  __is.flags(__ios_base::dec | __ios_base::skipws);
1768 
1769  double __mean, __stddev;
1770  __is >> __mean >> __stddev
1771  >> __x._M_saved_available;
1772  if (__x._M_saved_available)
1773  __is >> __x._M_saved;
1774  __x.param(typename normal_distribution<_RealType>::
1775  param_type(__mean, __stddev));
1776 
1777  __is.flags(__flags);
1778  return __is;
1779  }
1780 
1781 
1782  template<typename _RealType, typename _CharT, typename _Traits>
1783  std::basic_ostream<_CharT, _Traits>&
1784  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1785  const lognormal_distribution<_RealType>& __x)
1786  {
1787  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1788  typedef typename __ostream_type::ios_base __ios_base;
1789 
1790  const typename __ios_base::fmtflags __flags = __os.flags();
1791  const _CharT __fill = __os.fill();
1792  const std::streamsize __precision = __os.precision();
1793  const _CharT __space = __os.widen(' ');
1794  __os.flags(__ios_base::scientific | __ios_base::left);
1795  __os.fill(__space);
1796  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1797 
1798  __os << __x.m() << __space << __x.s()
1799  << __space << __x._M_nd;
1800 
1801  __os.flags(__flags);
1802  __os.fill(__fill);
1803  __os.precision(__precision);
1804  return __os;
1805  }
1806 
1807  template<typename _RealType, typename _CharT, typename _Traits>
1808  std::basic_istream<_CharT, _Traits>&
1809  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1810  lognormal_distribution<_RealType>& __x)
1811  {
1812  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1813  typedef typename __istream_type::ios_base __ios_base;
1814 
1815  const typename __ios_base::fmtflags __flags = __is.flags();
1816  __is.flags(__ios_base::dec | __ios_base::skipws);
1817 
1818  _RealType __m, __s;
1819  __is >> __m >> __s >> __x._M_nd;
1820  __x.param(typename lognormal_distribution<_RealType>::
1821  param_type(__m, __s));
1822 
1823  __is.flags(__flags);
1824  return __is;
1825  }
1826 
1827 
1828  template<typename _RealType, typename _CharT, typename _Traits>
1829  std::basic_ostream<_CharT, _Traits>&
1830  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1831  const chi_squared_distribution<_RealType>& __x)
1832  {
1833  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1834  typedef typename __ostream_type::ios_base __ios_base;
1835 
1836  const typename __ios_base::fmtflags __flags = __os.flags();
1837  const _CharT __fill = __os.fill();
1838  const std::streamsize __precision = __os.precision();
1839  const _CharT __space = __os.widen(' ');
1840  __os.flags(__ios_base::scientific | __ios_base::left);
1841  __os.fill(__space);
1842  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1843 
1844  __os << __x.n() << __space << __x._M_gd;
1845 
1846  __os.flags(__flags);
1847  __os.fill(__fill);
1848  __os.precision(__precision);
1849  return __os;
1850  }
1851 
1852  template<typename _RealType, typename _CharT, typename _Traits>
1853  std::basic_istream<_CharT, _Traits>&
1854  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1855  chi_squared_distribution<_RealType>& __x)
1856  {
1857  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1858  typedef typename __istream_type::ios_base __ios_base;
1859 
1860  const typename __ios_base::fmtflags __flags = __is.flags();
1861  __is.flags(__ios_base::dec | __ios_base::skipws);
1862 
1863  _RealType __n;
1864  __is >> __n >> __x._M_gd;
1865  __x.param(typename chi_squared_distribution<_RealType>::
1866  param_type(__n));
1867 
1868  __is.flags(__flags);
1869  return __is;
1870  }
1871 
1872 
1873  template<typename _RealType>
1874  template<typename _UniformRandomNumberGenerator>
1875  typename cauchy_distribution<_RealType>::result_type
1876  cauchy_distribution<_RealType>::
1877  operator()(_UniformRandomNumberGenerator& __urng,
1878  const param_type& __p)
1879  {
1880  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1881  __aurng(__urng);
1882  _RealType __u;
1883  do
1884  __u = __aurng();
1885  while (__u == 0.5);
1886 
1887  const _RealType __pi = 3.1415926535897932384626433832795029L;
1888  return __p.a() + __p.b() * std::tan(__pi * __u);
1889  }
1890 
1891  template<typename _RealType, typename _CharT, typename _Traits>
1892  std::basic_ostream<_CharT, _Traits>&
1893  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1894  const cauchy_distribution<_RealType>& __x)
1895  {
1896  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1897  typedef typename __ostream_type::ios_base __ios_base;
1898 
1899  const typename __ios_base::fmtflags __flags = __os.flags();
1900  const _CharT __fill = __os.fill();
1901  const std::streamsize __precision = __os.precision();
1902  const _CharT __space = __os.widen(' ');
1903  __os.flags(__ios_base::scientific | __ios_base::left);
1904  __os.fill(__space);
1905  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1906 
1907  __os << __x.a() << __space << __x.b();
1908 
1909  __os.flags(__flags);
1910  __os.fill(__fill);
1911  __os.precision(__precision);
1912  return __os;
1913  }
1914 
1915  template<typename _RealType, typename _CharT, typename _Traits>
1916  std::basic_istream<_CharT, _Traits>&
1917  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1918  cauchy_distribution<_RealType>& __x)
1919  {
1920  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1921  typedef typename __istream_type::ios_base __ios_base;
1922 
1923  const typename __ios_base::fmtflags __flags = __is.flags();
1924  __is.flags(__ios_base::dec | __ios_base::skipws);
1925 
1926  _RealType __a, __b;
1927  __is >> __a >> __b;
1928  __x.param(typename cauchy_distribution<_RealType>::
1929  param_type(__a, __b));
1930 
1931  __is.flags(__flags);
1932  return __is;
1933  }
1934 
1935 
1936  template<typename _RealType, typename _CharT, typename _Traits>
1937  std::basic_ostream<_CharT, _Traits>&
1938  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1939  const fisher_f_distribution<_RealType>& __x)
1940  {
1941  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1942  typedef typename __ostream_type::ios_base __ios_base;
1943 
1944  const typename __ios_base::fmtflags __flags = __os.flags();
1945  const _CharT __fill = __os.fill();
1946  const std::streamsize __precision = __os.precision();
1947  const _CharT __space = __os.widen(' ');
1948  __os.flags(__ios_base::scientific | __ios_base::left);
1949  __os.fill(__space);
1950  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1951 
1952  __os << __x.m() << __space << __x.n()
1953  << __space << __x._M_gd_x << __space << __x._M_gd_y;
1954 
1955  __os.flags(__flags);
1956  __os.fill(__fill);
1957  __os.precision(__precision);
1958  return __os;
1959  }
1960 
1961  template<typename _RealType, typename _CharT, typename _Traits>
1962  std::basic_istream<_CharT, _Traits>&
1963  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1964  fisher_f_distribution<_RealType>& __x)
1965  {
1966  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1967  typedef typename __istream_type::ios_base __ios_base;
1968 
1969  const typename __ios_base::fmtflags __flags = __is.flags();
1970  __is.flags(__ios_base::dec | __ios_base::skipws);
1971 
1972  _RealType __m, __n;
1973  __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
1974  __x.param(typename fisher_f_distribution<_RealType>::
1975  param_type(__m, __n));
1976 
1977  __is.flags(__flags);
1978  return __is;
1979  }
1980 
1981 
1982  template<typename _RealType, typename _CharT, typename _Traits>
1983  std::basic_ostream<_CharT, _Traits>&
1984  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1985  const student_t_distribution<_RealType>& __x)
1986  {
1987  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1988  typedef typename __ostream_type::ios_base __ios_base;
1989 
1990  const typename __ios_base::fmtflags __flags = __os.flags();
1991  const _CharT __fill = __os.fill();
1992  const std::streamsize __precision = __os.precision();
1993  const _CharT __space = __os.widen(' ');
1994  __os.flags(__ios_base::scientific | __ios_base::left);
1995  __os.fill(__space);
1996  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1997 
1998  __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
1999 
2000  __os.flags(__flags);
2001  __os.fill(__fill);
2002  __os.precision(__precision);
2003  return __os;
2004  }
2005 
2006  template<typename _RealType, typename _CharT, typename _Traits>
2007  std::basic_istream<_CharT, _Traits>&
2008  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2009  student_t_distribution<_RealType>& __x)
2010  {
2011  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2012  typedef typename __istream_type::ios_base __ios_base;
2013 
2014  const typename __ios_base::fmtflags __flags = __is.flags();
2015  __is.flags(__ios_base::dec | __ios_base::skipws);
2016 
2017  _RealType __n;
2018  __is >> __n >> __x._M_nd >> __x._M_gd;
2019  __x.param(typename student_t_distribution<_RealType>::param_type(__n));
2020 
2021  __is.flags(__flags);
2022  return __is;
2023  }
2024 
2025 
2026  template<typename _RealType>
2027  void
2028  gamma_distribution<_RealType>::param_type::
2029  _M_initialize()
2030  {
2031  _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2032 
2033  const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2034  _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2035  }
2036 
2037  /**
2038  * Marsaglia, G. and Tsang, W. W.
2039  * "A Simple Method for Generating Gamma Variables"
2040  * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2041  */
2042  template<typename _RealType>
2043  template<typename _UniformRandomNumberGenerator>
2044  typename gamma_distribution<_RealType>::result_type
2045  gamma_distribution<_RealType>::
2046  operator()(_UniformRandomNumberGenerator& __urng,
2047  const param_type& __param)
2048  {
2049  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2050  __aurng(__urng);
2051 
2052  result_type __u, __v, __n;
2053  const result_type __a1 = (__param._M_malpha
2054  - _RealType(1.0) / _RealType(3.0));
2055 
2056  do
2057  {
2058  do
2059  {
2060  __n = _M_nd(__urng);
2061  __v = result_type(1.0) + __param._M_a2 * __n;
2062  }
2063  while (__v <= 0.0);
2064 
2065  __v = __v * __v * __v;
2066  __u = __aurng();
2067  }
2068  while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2069  && (std::log(__u) > (0.5 * __n * __n + __a1
2070  * (1.0 - __v + std::log(__v)))));
2071 
2072  if (__param.alpha() == __param._M_malpha)
2073  return __a1 * __v * __param.beta();
2074  else
2075  {
2076  do
2077  __u = __aurng();
2078  while (__u == 0.0);
2079 
2080  return (std::pow(__u, result_type(1.0) / __param.alpha())
2081  * __a1 * __v * __param.beta());
2082  }
2083  }
2084 
2085  template<typename _RealType, typename _CharT, typename _Traits>
2086  std::basic_ostream<_CharT, _Traits>&
2087  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2088  const gamma_distribution<_RealType>& __x)
2089  {
2090  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2091  typedef typename __ostream_type::ios_base __ios_base;
2092 
2093  const typename __ios_base::fmtflags __flags = __os.flags();
2094  const _CharT __fill = __os.fill();
2095  const std::streamsize __precision = __os.precision();
2096  const _CharT __space = __os.widen(' ');
2097  __os.flags(__ios_base::scientific | __ios_base::left);
2098  __os.fill(__space);
2099  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2100 
2101  __os << __x.alpha() << __space << __x.beta()
2102  << __space << __x._M_nd;
2103 
2104  __os.flags(__flags);
2105  __os.fill(__fill);
2106  __os.precision(__precision);
2107  return __os;
2108  }
2109 
2110  template<typename _RealType, typename _CharT, typename _Traits>
2111  std::basic_istream<_CharT, _Traits>&
2112  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2113  gamma_distribution<_RealType>& __x)
2114  {
2115  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2116  typedef typename __istream_type::ios_base __ios_base;
2117 
2118  const typename __ios_base::fmtflags __flags = __is.flags();
2119  __is.flags(__ios_base::dec | __ios_base::skipws);
2120 
2121  _RealType __alpha_val, __beta_val;
2122  __is >> __alpha_val >> __beta_val >> __x._M_nd;
2123  __x.param(typename gamma_distribution<_RealType>::
2124  param_type(__alpha_val, __beta_val));
2125 
2126  __is.flags(__flags);
2127  return __is;
2128  }
2129 
2130 
2131  template<typename _RealType>
2132  template<typename _UniformRandomNumberGenerator>
2133  typename weibull_distribution<_RealType>::result_type
2134  weibull_distribution<_RealType>::
2135  operator()(_UniformRandomNumberGenerator& __urng,
2136  const param_type& __p)
2137  {
2138  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2139  __aurng(__urng);
2140  return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2141  result_type(1) / __p.a());
2142  }
2143 
2144  template<typename _RealType, typename _CharT, typename _Traits>
2145  std::basic_ostream<_CharT, _Traits>&
2146  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2147  const weibull_distribution<_RealType>& __x)
2148  {
2149  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2150  typedef typename __ostream_type::ios_base __ios_base;
2151 
2152  const typename __ios_base::fmtflags __flags = __os.flags();
2153  const _CharT __fill = __os.fill();
2154  const std::streamsize __precision = __os.precision();
2155  const _CharT __space = __os.widen(' ');
2156  __os.flags(__ios_base::scientific | __ios_base::left);
2157  __os.fill(__space);
2158  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2159 
2160  __os << __x.a() << __space << __x.b();
2161 
2162  __os.flags(__flags);
2163  __os.fill(__fill);
2164  __os.precision(__precision);
2165  return __os;
2166  }
2167 
2168  template<typename _RealType, typename _CharT, typename _Traits>
2169  std::basic_istream<_CharT, _Traits>&
2170  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2171  weibull_distribution<_RealType>& __x)
2172  {
2173  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2174  typedef typename __istream_type::ios_base __ios_base;
2175 
2176  const typename __ios_base::fmtflags __flags = __is.flags();
2177  __is.flags(__ios_base::dec | __ios_base::skipws);
2178 
2179  _RealType __a, __b;
2180  __is >> __a >> __b;
2181  __x.param(typename weibull_distribution<_RealType>::
2182  param_type(__a, __b));
2183 
2184  __is.flags(__flags);
2185  return __is;
2186  }
2187 
2188 
2189  template<typename _RealType>
2190  template<typename _UniformRandomNumberGenerator>
2191  typename extreme_value_distribution<_RealType>::result_type
2192  extreme_value_distribution<_RealType>::
2193  operator()(_UniformRandomNumberGenerator& __urng,
2194  const param_type& __p)
2195  {
2196  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2197  __aurng(__urng);
2198  return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2199  - __aurng()));
2200  }
2201 
2202  template<typename _RealType, typename _CharT, typename _Traits>
2203  std::basic_ostream<_CharT, _Traits>&
2204  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2205  const extreme_value_distribution<_RealType>& __x)
2206  {
2207  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2208  typedef typename __ostream_type::ios_base __ios_base;
2209 
2210  const typename __ios_base::fmtflags __flags = __os.flags();
2211  const _CharT __fill = __os.fill();
2212  const std::streamsize __precision = __os.precision();
2213  const _CharT __space = __os.widen(' ');
2214  __os.flags(__ios_base::scientific | __ios_base::left);
2215  __os.fill(__space);
2216  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2217 
2218  __os << __x.a() << __space << __x.b();
2219 
2220  __os.flags(__flags);
2221  __os.fill(__fill);
2222  __os.precision(__precision);
2223  return __os;
2224  }
2225 
2226  template<typename _RealType, typename _CharT, typename _Traits>
2227  std::basic_istream<_CharT, _Traits>&
2228  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2229  extreme_value_distribution<_RealType>& __x)
2230  {
2231  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2232  typedef typename __istream_type::ios_base __ios_base;
2233 
2234  const typename __ios_base::fmtflags __flags = __is.flags();
2235  __is.flags(__ios_base::dec | __ios_base::skipws);
2236 
2237  _RealType __a, __b;
2238  __is >> __a >> __b;
2239  __x.param(typename extreme_value_distribution<_RealType>::
2240  param_type(__a, __b));
2241 
2242  __is.flags(__flags);
2243  return __is;
2244  }
2245 
2246 
2247  template<typename _IntType>
2248  void
2249  discrete_distribution<_IntType>::param_type::
2250  _M_initialize()
2251  {
2252  if (_M_prob.size() < 2)
2253  {
2254  _M_prob.clear();
2255  return;
2256  }
2257 
2258  const double __sum = std::accumulate(_M_prob.begin(),
2259  _M_prob.end(), 0.0);
2260  // Now normalize the probabilites.
2261  __detail::__transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2262  std::bind2nd(std::divides<double>(), __sum));
2263  // Accumulate partial sums.
2264  _M_cp.reserve(_M_prob.size());
2265  std::partial_sum(_M_prob.begin(), _M_prob.end(),
2266  std::back_inserter(_M_cp));
2267  // Make sure the last cumulative probability is one.
2268  _M_cp[_M_cp.size() - 1] = 1.0;
2269  }
2270 
2271  template<typename _IntType>
2272  template<typename _Func>
2273  discrete_distribution<_IntType>::param_type::
2274  param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2275  : _M_prob(), _M_cp()
2276  {
2277  const size_t __n = __nw == 0 ? 1 : __nw;
2278  const double __delta = (__xmax - __xmin) / __n;
2279 
2280  _M_prob.reserve(__n);
2281  for (size_t __k = 0; __k < __nw; ++__k)
2282  _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2283 
2284  _M_initialize();
2285  }
2286 
2287  template<typename _IntType>
2288  template<typename _UniformRandomNumberGenerator>
2289  typename discrete_distribution<_IntType>::result_type
2290  discrete_distribution<_IntType>::
2291  operator()(_UniformRandomNumberGenerator& __urng,
2292  const param_type& __param)
2293  {
2294  if (__param._M_cp.empty())
2295  return result_type(0);
2296 
2297  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2298  __aurng(__urng);
2299 
2300  const double __p = __aurng();
2301  auto __pos = std::lower_bound(__param._M_cp.begin(),
2302  __param._M_cp.end(), __p);
2303 
2304  return __pos - __param._M_cp.begin();
2305  }
2306 
2307  template<typename _IntType, typename _CharT, typename _Traits>
2308  std::basic_ostream<_CharT, _Traits>&
2309  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2310  const discrete_distribution<_IntType>& __x)
2311  {
2312  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2313  typedef typename __ostream_type::ios_base __ios_base;
2314 
2315  const typename __ios_base::fmtflags __flags = __os.flags();
2316  const _CharT __fill = __os.fill();
2317  const std::streamsize __precision = __os.precision();
2318  const _CharT __space = __os.widen(' ');
2319  __os.flags(__ios_base::scientific | __ios_base::left);
2320  __os.fill(__space);
2321  __os.precision(std::numeric_limits<double>::max_digits10);
2322 
2323  std::vector<double> __prob = __x.probabilities();
2324  __os << __prob.size();
2325  for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2326  __os << __space << *__dit;
2327 
2328  __os.flags(__flags);
2329  __os.fill(__fill);
2330  __os.precision(__precision);
2331  return __os;
2332  }
2333 
2334  template<typename _IntType, typename _CharT, typename _Traits>
2335  std::basic_istream<_CharT, _Traits>&
2336  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2337  discrete_distribution<_IntType>& __x)
2338  {
2339  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2340  typedef typename __istream_type::ios_base __ios_base;
2341 
2342  const typename __ios_base::fmtflags __flags = __is.flags();
2343  __is.flags(__ios_base::dec | __ios_base::skipws);
2344 
2345  size_t __n;
2346  __is >> __n;
2347 
2348  std::vector<double> __prob_vec;
2349  __prob_vec.reserve(__n);
2350  for (; __n != 0; --__n)
2351  {
2352  double __prob;
2353  __is >> __prob;
2354  __prob_vec.push_back(__prob);
2355  }
2356 
2357  __x.param(typename discrete_distribution<_IntType>::
2358  param_type(__prob_vec.begin(), __prob_vec.end()));
2359 
2360  __is.flags(__flags);
2361  return __is;
2362  }
2363 
2364 
2365  template<typename _RealType>
2366  void
2367  piecewise_constant_distribution<_RealType>::param_type::
2368  _M_initialize()
2369  {
2370  if (_M_int.size() < 2
2371  || (_M_int.size() == 2
2372  && _M_int[0] == _RealType(0)
2373  && _M_int[1] == _RealType(1)))
2374  {
2375  _M_int.clear();
2376  _M_den.clear();
2377  return;
2378  }
2379 
2380  const double __sum = std::accumulate(_M_den.begin(),
2381  _M_den.end(), 0.0);
2382 
2383  __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2384  std::bind2nd(std::divides<double>(), __sum));
2385 
2386  _M_cp.reserve(_M_den.size());
2387  std::partial_sum(_M_den.begin(), _M_den.end(),
2388  std::back_inserter(_M_cp));
2389 
2390  // Make sure the last cumulative probability is one.
2391  _M_cp[_M_cp.size() - 1] = 1.0;
2392 
2393  for (size_t __k = 0; __k < _M_den.size(); ++__k)
2394  _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2395  }
2396 
2397  template<typename _RealType>
2398  template<typename _InputIteratorB, typename _InputIteratorW>
2399  piecewise_constant_distribution<_RealType>::param_type::
2400  param_type(_InputIteratorB __bbegin,
2401  _InputIteratorB __bend,
2402  _InputIteratorW __wbegin)
2403  : _M_int(), _M_den(), _M_cp()
2404  {
2405  if (__bbegin != __bend)
2406  {
2407  for (;;)
2408  {
2409  _M_int.push_back(*__bbegin);
2410  ++__bbegin;
2411  if (__bbegin == __bend)
2412  break;
2413 
2414  _M_den.push_back(*__wbegin);
2415  ++__wbegin;
2416  }
2417  }
2418 
2419  _M_initialize();
2420  }
2421 
2422  template<typename _RealType>
2423  template<typename _Func>
2424  piecewise_constant_distribution<_RealType>::param_type::
2425  param_type(initializer_list<_RealType> __bl, _Func __fw)
2426  : _M_int(), _M_den(), _M_cp()
2427  {
2428  _M_int.reserve(__bl.size());
2429  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2430  _M_int.push_back(*__biter);
2431 
2432  _M_den.reserve(_M_int.size() - 1);
2433  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2434  _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2435 
2436  _M_initialize();
2437  }
2438 
2439  template<typename _RealType>
2440  template<typename _Func>
2441  piecewise_constant_distribution<_RealType>::param_type::
2442  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2443  : _M_int(), _M_den(), _M_cp()
2444  {
2445  const size_t __n = __nw == 0 ? 1 : __nw;
2446  const _RealType __delta = (__xmax - __xmin) / __n;
2447 
2448  _M_int.reserve(__n + 1);
2449  for (size_t __k = 0; __k <= __nw; ++__k)
2450  _M_int.push_back(__xmin + __k * __delta);
2451 
2452  _M_den.reserve(__n);
2453  for (size_t __k = 0; __k < __nw; ++__k)
2454  _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2455 
2456  _M_initialize();
2457  }
2458 
2459  template<typename _RealType>
2460  template<typename _UniformRandomNumberGenerator>
2461  typename piecewise_constant_distribution<_RealType>::result_type
2462  piecewise_constant_distribution<_RealType>::
2463  operator()(_UniformRandomNumberGenerator& __urng,
2464  const param_type& __param)
2465  {
2466  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2467  __aurng(__urng);
2468 
2469  const double __p = __aurng();
2470  if (__param._M_cp.empty())
2471  return __p;
2472 
2473  auto __pos = std::lower_bound(__param._M_cp.begin(),
2474  __param._M_cp.end(), __p);
2475  const size_t __i = __pos - __param._M_cp.begin();
2476 
2477  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2478 
2479  return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2480  }
2481 
2482  template<typename _RealType, typename _CharT, typename _Traits>
2483  std::basic_ostream<_CharT, _Traits>&
2484  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2485  const piecewise_constant_distribution<_RealType>& __x)
2486  {
2487  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2488  typedef typename __ostream_type::ios_base __ios_base;
2489 
2490  const typename __ios_base::fmtflags __flags = __os.flags();
2491  const _CharT __fill = __os.fill();
2492  const std::streamsize __precision = __os.precision();
2493  const _CharT __space = __os.widen(' ');
2494  __os.flags(__ios_base::scientific | __ios_base::left);
2495  __os.fill(__space);
2496  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2497 
2498  std::vector<_RealType> __int = __x.intervals();
2499  __os << __int.size() - 1;
2500 
2501  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2502  __os << __space << *__xit;
2503 
2504  std::vector<double> __den = __x.densities();
2505  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2506  __os << __space << *__dit;
2507 
2508  __os.flags(__flags);
2509  __os.fill(__fill);
2510  __os.precision(__precision);
2511  return __os;
2512  }
2513 
2514  template<typename _RealType, typename _CharT, typename _Traits>
2515  std::basic_istream<_CharT, _Traits>&
2516  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2517  piecewise_constant_distribution<_RealType>& __x)
2518  {
2519  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2520  typedef typename __istream_type::ios_base __ios_base;
2521 
2522  const typename __ios_base::fmtflags __flags = __is.flags();
2523  __is.flags(__ios_base::dec | __ios_base::skipws);
2524 
2525  size_t __n;
2526  __is >> __n;
2527 
2528  std::vector<_RealType> __int_vec;
2529  __int_vec.reserve(__n + 1);
2530  for (size_t __i = 0; __i <= __n; ++__i)
2531  {
2532  _RealType __int;
2533  __is >> __int;
2534  __int_vec.push_back(__int);
2535  }
2536 
2537  std::vector<double> __den_vec;
2538  __den_vec.reserve(__n);
2539  for (size_t __i = 0; __i < __n; ++__i)
2540  {
2541  double __den;
2542  __is >> __den;
2543  __den_vec.push_back(__den);
2544  }
2545 
2546  __x.param(typename piecewise_constant_distribution<_RealType>::
2547  param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2548 
2549  __is.flags(__flags);
2550  return __is;
2551  }
2552 
2553 
2554  template<typename _RealType>
2555  void
2556  piecewise_linear_distribution<_RealType>::param_type::
2557  _M_initialize()
2558  {
2559  if (_M_int.size() < 2
2560  || (_M_int.size() == 2
2561  && _M_int[0] == _RealType(0)
2562  && _M_int[1] == _RealType(1)
2563  && _M_den[0] == _M_den[1]))
2564  {
2565  _M_int.clear();
2566  _M_den.clear();
2567  return;
2568  }
2569 
2570  double __sum = 0.0;
2571  _M_cp.reserve(_M_int.size() - 1);
2572  _M_m.reserve(_M_int.size() - 1);
2573  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2574  {
2575  const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
2576  __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
2577  _M_cp.push_back(__sum);
2578  _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
2579  }
2580 
2581  // Now normalize the densities...
2582  __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2583  std::bind2nd(std::divides<double>(), __sum));
2584  // ... and partial sums...
2585  __detail::__transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
2586  std::bind2nd(std::divides<double>(), __sum));
2587  // ... and slopes.
2588  __detail::__transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
2589  std::bind2nd(std::divides<double>(), __sum));
2590  // Make sure the last cumulative probablility is one.
2591  _M_cp[_M_cp.size() - 1] = 1.0;
2592  }
2593 
2594  template<typename _RealType>
2595  template<typename _InputIteratorB, typename _InputIteratorW>
2596  piecewise_linear_distribution<_RealType>::param_type::
2597  param_type(_InputIteratorB __bbegin,
2598  _InputIteratorB __bend,
2599  _InputIteratorW __wbegin)
2600  : _M_int(), _M_den(), _M_cp(), _M_m()
2601  {
2602  for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
2603  {
2604  _M_int.push_back(*__bbegin);
2605  _M_den.push_back(*__wbegin);
2606  }
2607 
2608  _M_initialize();
2609  }
2610 
2611  template<typename _RealType>
2612  template<typename _Func>
2613  piecewise_linear_distribution<_RealType>::param_type::
2614  param_type(initializer_list<_RealType> __bl, _Func __fw)
2615  : _M_int(), _M_den(), _M_cp(), _M_m()
2616  {
2617  _M_int.reserve(__bl.size());
2618  _M_den.reserve(__bl.size());
2619  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2620  {
2621  _M_int.push_back(*__biter);
2622  _M_den.push_back(__fw(*__biter));
2623  }
2624 
2625  _M_initialize();
2626  }
2627 
2628  template<typename _RealType>
2629  template<typename _Func>
2630  piecewise_linear_distribution<_RealType>::param_type::
2631  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2632  : _M_int(), _M_den(), _M_cp(), _M_m()
2633  {
2634  const size_t __n = __nw == 0 ? 1 : __nw;
2635  const _RealType __delta = (__xmax - __xmin) / __n;
2636 
2637  _M_int.reserve(__n + 1);
2638  _M_den.reserve(__n + 1);
2639  for (size_t __k = 0; __k <= __nw; ++__k)
2640  {
2641  _M_int.push_back(__xmin + __k * __delta);
2642  _M_den.push_back(__fw(_M_int[__k] + __delta));
2643  }
2644 
2645  _M_initialize();
2646  }
2647 
2648  template<typename _RealType>
2649  template<typename _UniformRandomNumberGenerator>
2650  typename piecewise_linear_distribution<_RealType>::result_type
2651  piecewise_linear_distribution<_RealType>::
2652  operator()(_UniformRandomNumberGenerator& __urng,
2653  const param_type& __param)
2654  {
2655  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2656  __aurng(__urng);
2657 
2658  const double __p = __aurng();
2659  if (__param._M_cp.empty())
2660  return __p;
2661 
2662  auto __pos = std::lower_bound(__param._M_cp.begin(),
2663  __param._M_cp.end(), __p);
2664  const size_t __i = __pos - __param._M_cp.begin();
2665 
2666  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2667 
2668  const double __a = 0.5 * __param._M_m[__i];
2669  const double __b = __param._M_den[__i];
2670  const double __cm = __p - __pref;
2671 
2672  _RealType __x = __param._M_int[__i];
2673  if (__a == 0)
2674  __x += __cm / __b;
2675  else
2676  {
2677  const double __d = __b * __b + 4.0 * __a * __cm;
2678  __x += 0.5 * (std::sqrt(__d) - __b) / __a;
2679  }
2680 
2681  return __x;
2682  }
2683 
2684  template<typename _RealType, typename _CharT, typename _Traits>
2685  std::basic_ostream<_CharT, _Traits>&
2686  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2687  const piecewise_linear_distribution<_RealType>& __x)
2688  {
2689  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
2690  typedef typename __ostream_type::ios_base __ios_base;
2691 
2692  const typename __ios_base::fmtflags __flags = __os.flags();
2693  const _CharT __fill = __os.fill();
2694  const std::streamsize __precision = __os.precision();
2695  const _CharT __space = __os.widen(' ');
2696  __os.flags(__ios_base::scientific | __ios_base::left);
2697  __os.fill(__space);
2698  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2699 
2700  std::vector<_RealType> __int = __x.intervals();
2701  __os << __int.size() - 1;
2702 
2703  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2704  __os << __space << *__xit;
2705 
2706  std::vector<double> __den = __x.densities();
2707  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2708  __os << __space << *__dit;
2709 
2710  __os.flags(__flags);
2711  __os.fill(__fill);
2712  __os.precision(__precision);
2713  return __os;
2714  }
2715 
2716  template<typename _RealType, typename _CharT, typename _Traits>
2717  std::basic_istream<_CharT, _Traits>&
2718  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2719  piecewise_linear_distribution<_RealType>& __x)
2720  {
2721  typedef std::basic_istream<_CharT, _Traits> __istream_type;
2722  typedef typename __istream_type::ios_base __ios_base;
2723 
2724  const typename __ios_base::fmtflags __flags = __is.flags();
2725  __is.flags(__ios_base::dec | __ios_base::skipws);
2726 
2727  size_t __n;
2728  __is >> __n;
2729 
2730  std::vector<_RealType> __int_vec;
2731  __int_vec.reserve(__n + 1);
2732  for (size_t __i = 0; __i <= __n; ++__i)
2733  {
2734  _RealType __int;
2735  __is >> __int;
2736  __int_vec.push_back(__int);
2737  }
2738 
2739  std::vector<double> __den_vec;
2740  __den_vec.reserve(__n + 1);
2741  for (size_t __i = 0; __i <= __n; ++__i)
2742  {
2743  double __den;
2744  __is >> __den;
2745  __den_vec.push_back(__den);
2746  }
2747 
2748  __x.param(typename piecewise_linear_distribution<_RealType>::
2749  param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2750 
2751  __is.flags(__flags);
2752  return __is;
2753  }
2754 
2755 
2756  template<typename _IntType>
2757  seed_seq::seed_seq(std::initializer_list<_IntType> __il)
2758  {
2759  for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
2760  _M_v.push_back(__detail::__mod<result_type,
2761  __detail::_Shift<result_type, 32>::__value>(*__iter));
2762  }
2763 
2764  template<typename _InputIterator>
2765  seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
2766  {
2767  for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
2768  _M_v.push_back(__detail::__mod<result_type,
2769  __detail::_Shift<result_type, 32>::__value>(*__iter));
2770  }
2771 
2772  template<typename _RandomAccessIterator>
2773  void
2774  seed_seq::generate(_RandomAccessIterator __begin,
2775  _RandomAccessIterator __end)
2776  {
2777  typedef typename iterator_traits<_RandomAccessIterator>::value_type
2778  _Type;
2779 
2780  if (__begin == __end)
2781  return;
2782 
2783  std::fill(__begin, __end, _Type(0x8b8b8b8bu));
2784 
2785  const size_t __n = __end - __begin;
2786  const size_t __s = _M_v.size();
2787  const size_t __t = (__n >= 623) ? 11
2788  : (__n >= 68) ? 7
2789  : (__n >= 39) ? 5
2790  : (__n >= 7) ? 3
2791  : (__n - 1) / 2;
2792  const size_t __p = (__n - __t) / 2;
2793  const size_t __q = __p + __t;
2794  const size_t __m = std::max(size_t(__s + 1), __n);
2795 
2796  for (size_t __k = 0; __k < __m; ++__k)
2797  {
2798  _Type __arg = (__begin[__k % __n]
2799  ^ __begin[(__k + __p) % __n]
2800  ^ __begin[(__k - 1) % __n]);
2801  _Type __r1 = __arg ^ (__arg >> 27);
2802  __r1 = __detail::__mod<_Type,
2803  __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
2804  _Type __r2 = __r1;
2805  if (__k == 0)
2806  __r2 += __s;
2807  else if (__k <= __s)
2808  __r2 += __k % __n + _M_v[__k - 1];
2809  else
2810  __r2 += __k % __n;
2811  __r2 = __detail::__mod<_Type,
2812  __detail::_Shift<_Type, 32>::__value>(__r2);
2813  __begin[(__k + __p) % __n] += __r1;
2814  __begin[(__k + __q) % __n] += __r2;
2815  __begin[__k % __n] = __r2;
2816  }
2817 
2818  for (size_t __k = __m; __k < __m + __n; ++__k)
2819  {
2820  _Type __arg = (__begin[__k % __n]
2821  + __begin[(__k + __p) % __n]
2822  + __begin[(__k - 1) % __n]);
2823  _Type __r3 = __arg ^ (__arg >> 27);
2824  __r3 = __detail::__mod<_Type,
2825  __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
2826  _Type __r4 = __r3 - __k % __n;
2827  __r4 = __detail::__mod<_Type,
2828  __detail::_Shift<_Type, 32>::__value>(__r4);
2829  __begin[(__k + __p) % __n] ^= __r3;
2830  __begin[(__k + __q) % __n] ^= __r4;
2831  __begin[__k % __n] = __r4;
2832  }
2833  }
2834 
2835  template<typename _RealType, size_t __bits,
2836  typename _UniformRandomNumberGenerator>
2837  _RealType
2838  generate_canonical(_UniformRandomNumberGenerator& __urng)
2839  {
2840  const size_t __b
2841  = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
2842  __bits);
2843  const long double __r = static_cast<long double>(__urng.max())
2844  - static_cast<long double>(__urng.min()) + 1.0L;
2845  const size_t __log2r = std::log(__r) / std::log(2.0L);
2846  size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
2847  _RealType __sum = _RealType(0);
2848  _RealType __tmp = _RealType(1);
2849  for (; __k != 0; --__k)
2850  {
2851  __sum += _RealType(__urng() - __urng.min()) * __tmp;
2852  __tmp *= __r;
2853  }
2854  return __sum / __tmp;
2855  }
2856 
2857 _GLIBCXX_END_NAMESPACE_VERSION
2858 } // namespace
2859 
2860 #endif
std::initializer_list
initializer_list
Definition: initializer_list:47
std::binomial_distribution
A discrete binomial random number distribution.
Definition: random.h:3400
std::uniform_int_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:1721
std::discard_block_engine::operator()
result_type operator()()
Gets the next value in the generated random number sequence.
Definition: random.tcc:675
std::normal_distribution
A normal continuous distribution for random numbers.
Definition: random.h:1979
std::scientific
ios_base & scientific(ios_base &__base)
Calls base.setf(ios_base::scientific, ios_base::floatfield).
Definition: ios_base.h:968
std::gamma_distribution
A gamma continuous distribution for random numbers.
Definition: random.h:2338
std::shuffle_order_engine::operator()
result_type operator()()
Definition: random.tcc:805
std::skipws
ios_base & skipws(ios_base &__base)
Calls base.setf(ios_base::skipws).
Definition: ios_base.h:861
std::numeric_limits::max
static constexpr _Tp max() _GLIBCXX_USE_NOEXCEPT
Definition: limits:313
std::common_type
common_type
Definition: type_traits:1741
std::bernoulli_distribution
A Bernoulli random number distribution.
Definition: random.h:3223
std::uniform_real_distribution
Uniform continuous distribution for random numbers.
Definition: random.h:1790
std::divides
One of the math functors.
Definition: stl_function.h:169
std::bind2nd
binder2nd< _Operation > bind2nd(const _Operation &__fn, const _Tp &__x)
One of the binder functors.
Definition: binders.h:167
std::weibull_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4416
std::extreme_value_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4560
std::exponential_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4207
std::basic_ostream
Template class basic_ostream.This is the base class for all output streams. It provides text formatti...
Definition: iosfwd:88
std::numeric_limits
Properties of fundamental types.
Definition: limits:304
std::normal_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2089
std::vector::end
iterator end() _GLIBCXX_NOEXCEPT
Definition: stl_vector.h:536
std::accumulate
_Tp accumulate(_InputIterator __first, _InputIterator __last, _Tp __init)
Accumulate values in a range.
Definition: stl_numeric.h:121
std::exponential_distribution
An exponential continuous distribution for random numbers.
Definition: random.h:4141
std::enable_if
Define a member typedef type only if a boolean constant is true.
Definition: type_traits:1718
std::poisson_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4056
std::tan
complex< _Tp > tan(const complex< _Tp > &)
Return complex tangent of z.
Definition: complex:918
std::abs
_Tp abs(const complex< _Tp > &)
Return magnitude of z.
Definition: complex:594
std::poisson_distribution::operator==
friend bool operator==(const poisson_distribution &__d1, const poisson_distribution &__d2)
Return true if two Poisson distributions have the same parameters and the sequences that would be gen...
Definition: random.h:4070
std::extreme_value_distribution
A extreme_value_distribution random number distribution.
Definition: random.h:4488
std::gamma_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2455
std::dec
ios_base & dec(ios_base &__base)
Calls base.setf(ios_base::dec, ios_base::basefield).
Definition: ios_base.h:935
std::__lg
_Size __lg(_Size __n)
This is a helper function for the sort routines and for random.tcc.
Definition: stl_algobase.h:973
std::vector::reserve
void reserve(size_type __n)
Attempt to preallocate enough memory for specified number of elements.
Definition: vector.tcc:67
std::uniform_int_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:1691
std::linear_congruential_engine::seed
void seed(result_type __s=default_seed)
Reseeds the linear_congruential_engine random number generator engine sequence to the seed __s...
Definition: random.tcc:136
std::geometric_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:3671
std::numeric_limits::epsilon
static constexpr _Tp epsilon() _GLIBCXX_USE_NOEXCEPT
Definition: limits:325
std::cauchy_distribution
A cauchy_distribution random number distribution.
Definition: random.h:2690
std::poisson_distribution
A discrete Poisson random number distribution.
Definition: random.h:3954
std::gamma_distribution::result_type
_RealType result_type
Definition: random.h:2341
std::fixed
ios_base & fixed(ios_base &__base)
Calls base.setf(ios_base::fixed, ios_base::floatfield).
Definition: ios_base.h:960
std::exp
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.
Definition: complex:755
std::geometric_distribution
A discrete geometric random number distribution.
Definition: random.h:3601
std::vector::size
size_type size() const _GLIBCXX_NOEXCEPT
Definition: stl_vector.h:625
std::uniform_real_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:1866
std::log
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
Definition: complex:782
std::generate_canonical
_RealType generate_canonical(_UniformRandomNumberGenerator &__g)
A function template for converting the output of a (integral) uniform random number generator to a fl...
Definition: random.tcc:2838
std::min
const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:187
std::shuffle_order_engine
Produces random numbers by combining random numbers from some base engine to produce random numbers w...
Definition: random.h:1228
std::weibull_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:4386
std::streamsize
ptrdiff_t streamsize
Integral type for I/O operation counts and buffer sizes.
Definition: postypes.h:100
std::basic_istream
Template class basic_istream.This is the base class for all input streams. It provides text formattin...
Definition: iosfwd:85
std::extreme_value_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:4590
std::ios_base::flags
fmtflags flags() const
Access to format flags.
Definition: ios_base.h:553
std::sqrt
complex< _Tp > sqrt(const complex< _Tp > &)
Return complex square root of z.
Definition: complex:891
std::binomial_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:3518
std::partial_sum
_OutputIterator partial_sum(_InputIterator __first, _InputIterator __last, _OutputIterator __result)
Return list of partial sums.
Definition: stl_numeric.h:238
std::negative_binomial_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
std::basic_ios::clear
void clear(iostate __state=goodbit)
[Re]sets the error state.
Definition: basic_ios.tcc:42
std::vector::begin
iterator begin() _GLIBCXX_NOEXCEPT
Definition: stl_vector.h:518
std::uniform_int_distribution
Uniform discrete distribution for random numbers. A discrete random distribution on the range with e...
Definition: random.h:1618
std::pow
complex< _Tp > pow(const complex< _Tp > &, const _Tp &)
Return x to the y'th power.
Definition: complex:984
std::poisson_distribution::operator>>
friend std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, std::poisson_distribution< _IntType1 > &__x)
Extracts a poisson_distribution random number distribution __x from the input stream __is...
std::weibull_distribution
A weibull_distribution random number distribution.
Definition: random.h:4314
std::operator>>
bitset< _Nb > operator>>(size_t __position) const _GLIBCXX_NOEXCEPT
Self-explanatory.
Definition: bitset:1345
std::linear_congruential_engine
A model of a linear congruential random number generator.
Definition: random.h:170
std::max
const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:210
std::left
ios_base & left(ios_base &__base)
Calls base.setf(ios_base::left, ios_base::adjustfield).
Definition: ios_base.h:918
std::vector::push_back
void push_back(const value_type &__x)
Add data to the end of the vector.
Definition: stl_vector.h:881
std::back_inserter
back_insert_iterator< _Container > back_inserter(_Container &__x)
Definition: stl_iterator.h:481
std::lower_bound
_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:937
std::independent_bits_engine::operator()
result_type operator()()
Gets the next value in the generated random number sequence.
Definition: random.tcc:732
std::numeric_limits::min
static constexpr _Tp min() _GLIBCXX_USE_NOEXCEPT
Definition: limits:309
std::geometric_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:3701
std::cauchy_distribution::param
param_type param() const
Returns the parameter set of the distribution.
Definition: random.h:2759
std::shuffle_order_engine::result_type
_RandomNumberEngine::result_type result_type
Definition: random.h:1231
std::cauchy_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: random.h:2789