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