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