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