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