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