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