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