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]);
949  ++__uh;
950  __type __lh
951  = (static_cast<__type>(_M_x[3]) << __w)
952  | static_cast<__type>(_M_x[2]);
953  __type __bigMask
954  = ~__type(0) >> ((sizeof(__type) * __CHAR_BIT__) - (__w * 2));
955  if ((__uh & __bigMask) == 0)
956  {
957  ++__lh;
958  __uh = 0;
959  }
960  _M_x[0] = static_cast<_UIntType>(__uh & max());
961  _M_x[1] = static_cast<_UIntType>((__uh >> (__w)) & max());
962  _M_x[2] = static_cast<_UIntType>(__lh & max());
963  _M_x[3] = static_cast<_UIntType>((__lh >> (__w)) & max());
964  }
965  else
966  {
967  __type __num
968  = (static_cast<__type>(_M_x[1]) << __w)
969  | static_cast<__type>(_M_x[0]);
970  ++__num;
971  _M_x[0] = static_cast<_UIntType>(__num & max());
972  _M_x[1] = static_cast<_UIntType>((__num >> __w) & max());
973  }
974  _M_i = 0;
975  }
976 
977  template<typename _UIntType, size_t __w, size_t __n, size_t __r,
978  _UIntType... __consts>
979  void
980  philox_engine<_UIntType, __w, __n, __r, __consts...>::_M_philox()
981  {
982  array<_UIntType, __n> __outputSeq = _M_x;
983  for (size_t __j = 0; __j < __r; ++__j)
984  {
985  array<_UIntType, __n> __intermedSeq{};
986  if constexpr (__n == 4)
987  {
988  __intermedSeq[0] = __outputSeq[2];
989  __intermedSeq[1] = __outputSeq[1];
990  __intermedSeq[2] = __outputSeq[0];
991  __intermedSeq[3] = __outputSeq[3];
992  }
993  else
994  {
995  __intermedSeq[0] = __outputSeq[0];
996  __intermedSeq[1] = __outputSeq[1];
997  }
998  for (unsigned long __k = 0; __k < (__n/2); ++__k)
999  {
1000  __outputSeq[2*__k]
1001  = _S_mulhi(__intermedSeq[2*__k], multipliers[__k])
1002  ^ (((_M_k[__k] + (__j * round_consts[__k])) & max()))
1003  ^ __intermedSeq[2*__k+1];
1004 
1005  __outputSeq[(2*__k)+1]
1006  = _S_mullo(__intermedSeq[2*__k], multipliers[__k]);
1007  }
1008  }
1009  _M_y = __outputSeq;
1010  }
1011 
1012  template<typename _UIntType, size_t __w, size_t __n, size_t __r,
1013  _UIntType... __consts>
1014  template<typename _Sseq>
1015  void
1016  philox_engine<_UIntType, __w, __n, __r, __consts...>::seed(_Sseq& __q)
1017  requires __is_seed_seq<_Sseq>
1018  {
1019  seed(0);
1020 
1021  const unsigned __p = 1 + ((__w - 1) / 32);
1022  uint_least32_t __tmpArr[(__n/2) * __p];
1023  __q.generate(__tmpArr + 0, __tmpArr + ((__n/2) * __p));
1024  for (unsigned __k = 0; __k < (__n/2); ++__k)
1025  {
1026  unsigned long long __precalc = 0;
1027  for (unsigned __j = 0; __j < __p; ++__j)
1028  {
1029  unsigned long long __multiplicand = (1ull << (32 * __j));
1030  __precalc += (__tmpArr[__k * __p + __j] * __multiplicand) & max();
1031  }
1032  _M_k[__k] = __precalc;
1033  }
1034  }
1035 #endif // philox_engine
1036 
1037  template<typename _IntType, typename _CharT, typename _Traits>
1038  std::basic_ostream<_CharT, _Traits>&
1039  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1040  const uniform_int_distribution<_IntType>& __x)
1041  {
1042  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1043 
1044  const typename __ios_base::fmtflags __flags = __os.flags();
1045  const _CharT __fill = __os.fill();
1046  const _CharT __space = __os.widen(' ');
1047  __os.flags(__ios_base::scientific | __ios_base::left);
1048  __os.fill(__space);
1049 
1050  __os << __x.a() << __space << __x.b();
1051 
1052  __os.flags(__flags);
1053  __os.fill(__fill);
1054  return __os;
1055  }
1056 
1057  template<typename _IntType, typename _CharT, typename _Traits>
1058  std::basic_istream<_CharT, _Traits>&
1059  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1060  uniform_int_distribution<_IntType>& __x)
1061  {
1062  using param_type
1063  = typename uniform_int_distribution<_IntType>::param_type;
1064  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1065 
1066  const typename __ios_base::fmtflags __flags = __is.flags();
1067  __is.flags(__ios_base::dec | __ios_base::skipws);
1068 
1069  _IntType __a, __b;
1070  if (__is >> __a >> __b)
1071  __x.param(param_type(__a, __b));
1072 
1073  __is.flags(__flags);
1074  return __is;
1075  }
1076 
1077 
1078  template<typename _RealType>
1079  template<typename _ForwardIterator,
1080  typename _UniformRandomNumberGenerator>
1081  void
1082  uniform_real_distribution<_RealType>::
1083  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1084  _UniformRandomNumberGenerator& __urng,
1085  const param_type& __p)
1086  {
1087  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1088  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1089  __aurng(__urng);
1090  auto __range = __p.b() - __p.a();
1091  while (__f != __t)
1092  *__f++ = __aurng() * __range + __p.a();
1093  }
1094 
1095  template<typename _RealType, typename _CharT, typename _Traits>
1096  std::basic_ostream<_CharT, _Traits>&
1097  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1098  const uniform_real_distribution<_RealType>& __x)
1099  {
1100  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1101 
1102  const typename __ios_base::fmtflags __flags = __os.flags();
1103  const _CharT __fill = __os.fill();
1104  const std::streamsize __precision = __os.precision();
1105  const _CharT __space = __os.widen(' ');
1106  __os.flags(__ios_base::scientific | __ios_base::left);
1107  __os.fill(__space);
1108  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1109 
1110  __os << __x.a() << __space << __x.b();
1111 
1112  __os.flags(__flags);
1113  __os.fill(__fill);
1114  __os.precision(__precision);
1115  return __os;
1116  }
1117 
1118  template<typename _RealType, typename _CharT, typename _Traits>
1119  std::basic_istream<_CharT, _Traits>&
1120  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1121  uniform_real_distribution<_RealType>& __x)
1122  {
1123  using param_type
1124  = typename uniform_real_distribution<_RealType>::param_type;
1125  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1126 
1127  const typename __ios_base::fmtflags __flags = __is.flags();
1128  __is.flags(__ios_base::skipws);
1129 
1130  _RealType __a, __b;
1131  if (__is >> __a >> __b)
1132  __x.param(param_type(__a, __b));
1133 
1134  __is.flags(__flags);
1135  return __is;
1136  }
1137 
1138 
1139  template<typename _ForwardIterator,
1140  typename _UniformRandomNumberGenerator>
1141  void
1142  std::bernoulli_distribution::
1143  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1144  _UniformRandomNumberGenerator& __urng,
1145  const param_type& __p)
1146  {
1147  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1148  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1149  __aurng(__urng);
1150  auto __limit = __p.p() * (__aurng.max() - __aurng.min());
1151 
1152  while (__f != __t)
1153  *__f++ = (__aurng() - __aurng.min()) < __limit;
1154  }
1155 
1156  template<typename _CharT, typename _Traits>
1157  std::basic_ostream<_CharT, _Traits>&
1158  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1159  const bernoulli_distribution& __x)
1160  {
1161  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1162 
1163  const typename __ios_base::fmtflags __flags = __os.flags();
1164  const _CharT __fill = __os.fill();
1165  const std::streamsize __precision = __os.precision();
1166  __os.flags(__ios_base::scientific | __ios_base::left);
1167  __os.fill(__os.widen(' '));
1168  __os.precision(std::numeric_limits<double>::max_digits10);
1169 
1170  __os << __x.p();
1171 
1172  __os.flags(__flags);
1173  __os.fill(__fill);
1174  __os.precision(__precision);
1175  return __os;
1176  }
1177 
1178 
1179  template<typename _IntType>
1180  template<typename _UniformRandomNumberGenerator>
1181  typename geometric_distribution<_IntType>::result_type
1182  geometric_distribution<_IntType>::
1183  operator()(_UniformRandomNumberGenerator& __urng,
1184  const param_type& __param)
1185  {
1186  // About the epsilon thing see this thread:
1187  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1188  const double __naf =
1189  (1 - std::numeric_limits<double>::epsilon()) / 2;
1190  // The largest _RealType convertible to _IntType.
1191  const double __thr =
1192  std::numeric_limits<_IntType>::max() + __naf;
1193  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1194  __aurng(__urng);
1195 
1196  double __cand;
1197  do
1198  __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1199  while (__cand >= __thr);
1200 
1201  return result_type(__cand + __naf);
1202  }
1203 
1204  template<typename _IntType>
1205  template<typename _ForwardIterator,
1206  typename _UniformRandomNumberGenerator>
1207  void
1208  geometric_distribution<_IntType>::
1209  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1210  _UniformRandomNumberGenerator& __urng,
1211  const param_type& __param)
1212  {
1213  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1214  // About the epsilon thing see this thread:
1215  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1216  const double __naf =
1217  (1 - std::numeric_limits<double>::epsilon()) / 2;
1218  // The largest _RealType convertible to _IntType.
1219  const double __thr =
1220  std::numeric_limits<_IntType>::max() + __naf;
1221  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1222  __aurng(__urng);
1223 
1224  while (__f != __t)
1225  {
1226  double __cand;
1227  do
1228  __cand = std::floor(std::log(1.0 - __aurng())
1229  / __param._M_log_1_p);
1230  while (__cand >= __thr);
1231 
1232  *__f++ = __cand + __naf;
1233  }
1234  }
1235 
1236  template<typename _IntType,
1237  typename _CharT, typename _Traits>
1238  std::basic_ostream<_CharT, _Traits>&
1239  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1240  const geometric_distribution<_IntType>& __x)
1241  {
1242  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1243 
1244  const typename __ios_base::fmtflags __flags = __os.flags();
1245  const _CharT __fill = __os.fill();
1246  const std::streamsize __precision = __os.precision();
1247  __os.flags(__ios_base::scientific | __ios_base::left);
1248  __os.fill(__os.widen(' '));
1249  __os.precision(std::numeric_limits<double>::max_digits10);
1250 
1251  __os << __x.p();
1252 
1253  __os.flags(__flags);
1254  __os.fill(__fill);
1255  __os.precision(__precision);
1256  return __os;
1257  }
1258 
1259  template<typename _IntType,
1260  typename _CharT, typename _Traits>
1261  std::basic_istream<_CharT, _Traits>&
1262  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1263  geometric_distribution<_IntType>& __x)
1264  {
1265  using param_type = typename geometric_distribution<_IntType>::param_type;
1266  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1267 
1268  const typename __ios_base::fmtflags __flags = __is.flags();
1269  __is.flags(__ios_base::skipws);
1270 
1271  double __p;
1272  if (__is >> __p)
1273  __x.param(param_type(__p));
1274 
1275  __is.flags(__flags);
1276  return __is;
1277  }
1278 
1279  // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1280  template<typename _IntType>
1281  template<typename _UniformRandomNumberGenerator>
1282  typename negative_binomial_distribution<_IntType>::result_type
1283  negative_binomial_distribution<_IntType>::
1284  operator()(_UniformRandomNumberGenerator& __urng)
1285  {
1286  const double __y = _M_gd(__urng);
1287 
1288  // XXX Is the constructor too slow?
1289  std::poisson_distribution<result_type> __poisson(__y);
1290  return __poisson(__urng);
1291  }
1292 
1293  template<typename _IntType>
1294  template<typename _UniformRandomNumberGenerator>
1295  typename negative_binomial_distribution<_IntType>::result_type
1296  negative_binomial_distribution<_IntType>::
1297  operator()(_UniformRandomNumberGenerator& __urng,
1298  const param_type& __p)
1299  {
1300  typedef typename std::gamma_distribution<double>::param_type
1301  param_type;
1302 
1303  const double __y =
1304  _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1305 
1306  std::poisson_distribution<result_type> __poisson(__y);
1307  return __poisson(__urng);
1308  }
1309 
1310  template<typename _IntType>
1311  template<typename _ForwardIterator,
1312  typename _UniformRandomNumberGenerator>
1313  void
1314  negative_binomial_distribution<_IntType>::
1315  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1316  _UniformRandomNumberGenerator& __urng)
1317  {
1318  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1319  while (__f != __t)
1320  {
1321  const double __y = _M_gd(__urng);
1322 
1323  // XXX Is the constructor too slow?
1324  std::poisson_distribution<result_type> __poisson(__y);
1325  *__f++ = __poisson(__urng);
1326  }
1327  }
1328 
1329  template<typename _IntType>
1330  template<typename _ForwardIterator,
1331  typename _UniformRandomNumberGenerator>
1332  void
1333  negative_binomial_distribution<_IntType>::
1334  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1335  _UniformRandomNumberGenerator& __urng,
1336  const param_type& __p)
1337  {
1338  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1339  typename std::gamma_distribution<result_type>::param_type
1340  __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1341 
1342  while (__f != __t)
1343  {
1344  const double __y = _M_gd(__urng, __p2);
1345 
1346  std::poisson_distribution<result_type> __poisson(__y);
1347  *__f++ = __poisson(__urng);
1348  }
1349  }
1350 
1351  template<typename _IntType, typename _CharT, typename _Traits>
1352  std::basic_ostream<_CharT, _Traits>&
1353  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1354  const negative_binomial_distribution<_IntType>& __x)
1355  {
1356  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1357 
1358  const typename __ios_base::fmtflags __flags = __os.flags();
1359  const _CharT __fill = __os.fill();
1360  const std::streamsize __precision = __os.precision();
1361  const _CharT __space = __os.widen(' ');
1362  __os.flags(__ios_base::scientific | __ios_base::left);
1363  __os.fill(__os.widen(' '));
1364  __os.precision(std::numeric_limits<double>::max_digits10);
1365 
1366  __os << __x.k() << __space << __x.p()
1367  << __space << __x._M_gd;
1368 
1369  __os.flags(__flags);
1370  __os.fill(__fill);
1371  __os.precision(__precision);
1372  return __os;
1373  }
1374 
1375  template<typename _IntType, typename _CharT, typename _Traits>
1376  std::basic_istream<_CharT, _Traits>&
1377  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1378  negative_binomial_distribution<_IntType>& __x)
1379  {
1380  using param_type
1381  = typename negative_binomial_distribution<_IntType>::param_type;
1382  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1383 
1384  const typename __ios_base::fmtflags __flags = __is.flags();
1385  __is.flags(__ios_base::skipws);
1386 
1387  _IntType __k;
1388  double __p;
1389  if (__is >> __k >> __p >> __x._M_gd)
1390  __x.param(param_type(__k, __p));
1391 
1392  __is.flags(__flags);
1393  return __is;
1394  }
1395 
1396 
1397  template<typename _IntType>
1398  void
1399  poisson_distribution<_IntType>::param_type::
1400  _M_initialize()
1401  {
1402 #if _GLIBCXX_USE_C99_MATH_FUNCS
1403  if (_M_mean >= 12)
1404  {
1405  const double __m = std::floor(_M_mean);
1406  _M_lm_thr = std::log(_M_mean);
1407  _M_lfm = std::lgamma(__m + 1);
1408  _M_sm = std::sqrt(__m);
1409 
1410  const double __pi_4 = 0.7853981633974483096156608458198757L;
1411  const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1412  / __pi_4));
1413  _M_d = std::round(std::max<double>(6.0, std::min(__m, __dx)));
1414  const double __cx = 2 * __m + _M_d;
1415  _M_scx = std::sqrt(__cx / 2);
1416  _M_1cx = 1 / __cx;
1417 
1418  _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1419  _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1420  / _M_d;
1421  }
1422  else
1423 #endif
1424  _M_lm_thr = std::exp(-_M_mean);
1425  }
1426 
1427  /**
1428  * A rejection algorithm when mean >= 12 and a simple method based
1429  * upon the multiplication of uniform random variates otherwise.
1430  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_FUNCS
1431  * is defined.
1432  *
1433  * Reference:
1434  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1435  * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1436  */
1437  template<typename _IntType>
1438  template<typename _UniformRandomNumberGenerator>
1439  typename poisson_distribution<_IntType>::result_type
1440  poisson_distribution<_IntType>::
1441  operator()(_UniformRandomNumberGenerator& __urng,
1442  const param_type& __param)
1443  {
1444  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1445  __aurng(__urng);
1446 #if _GLIBCXX_USE_C99_MATH_FUNCS
1447  if (__param.mean() >= 12)
1448  {
1449  double __x;
1450 
1451  // See comments above...
1452  const double __naf =
1453  (1 - std::numeric_limits<double>::epsilon()) / 2;
1454  const double __thr =
1455  std::numeric_limits<_IntType>::max() + __naf;
1456 
1457  const double __m = std::floor(__param.mean());
1458  // sqrt(pi / 2)
1459  const double __spi_2 = 1.2533141373155002512078826424055226L;
1460  const double __c1 = __param._M_sm * __spi_2;
1461  const double __c2 = __param._M_c2b + __c1;
1462  const double __c3 = __c2 + 1;
1463  const double __c4 = __c3 + 1;
1464  // 1 / 78
1465  const double __178 = 0.0128205128205128205128205128205128L;
1466  // e^(1 / 78)
1467  const double __e178 = 1.0129030479320018583185514777512983L;
1468  const double __c5 = __c4 + __e178;
1469  const double __c = __param._M_cb + __c5;
1470  const double __2cx = 2 * (2 * __m + __param._M_d);
1471 
1472  bool __reject = true;
1473  do
1474  {
1475  const double __u = __c * __aurng();
1476  const double __e = -std::log(1.0 - __aurng());
1477 
1478  double __w = 0.0;
1479 
1480  if (__u <= __c1)
1481  {
1482  const double __n = _M_nd(__urng);
1483  const double __y = -std::abs(__n) * __param._M_sm - 1;
1484  __x = std::floor(__y);
1485  __w = -__n * __n / 2;
1486  if (__x < -__m)
1487  continue;
1488  }
1489  else if (__u <= __c2)
1490  {
1491  const double __n = _M_nd(__urng);
1492  const double __y = 1 + std::abs(__n) * __param._M_scx;
1493  __x = std::ceil(__y);
1494  __w = __y * (2 - __y) * __param._M_1cx;
1495  if (__x > __param._M_d)
1496  continue;
1497  }
1498  else if (__u <= __c3)
1499  // NB: This case not in the book, nor in the Errata,
1500  // but should be ok...
1501  __x = -1;
1502  else if (__u <= __c4)
1503  __x = 0;
1504  else if (__u <= __c5)
1505  {
1506  __x = 1;
1507  // Only in the Errata, see libstdc++/83237.
1508  __w = __178;
1509  }
1510  else
1511  {
1512  const double __v = -std::log(1.0 - __aurng());
1513  const double __y = __param._M_d
1514  + __v * __2cx / __param._M_d;
1515  __x = std::ceil(__y);
1516  __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1517  }
1518 
1519  __reject = (__w - __e - __x * __param._M_lm_thr
1520  > __param._M_lfm - std::lgamma(__x + __m + 1));
1521 
1522  __reject |= __x + __m >= __thr;
1523 
1524  } while (__reject);
1525 
1526  return result_type(__x + __m + __naf);
1527  }
1528  else
1529 #endif
1530  {
1531  _IntType __x = 0;
1532  double __prod = 1.0;
1533 
1534  do
1535  {
1536  __prod *= __aurng();
1537  __x += 1;
1538  }
1539  while (__prod > __param._M_lm_thr);
1540 
1541  return __x - 1;
1542  }
1543  }
1544 
1545  template<typename _IntType>
1546  template<typename _ForwardIterator,
1547  typename _UniformRandomNumberGenerator>
1548  void
1549  poisson_distribution<_IntType>::
1550  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1551  _UniformRandomNumberGenerator& __urng,
1552  const param_type& __param)
1553  {
1554  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1555  // We could duplicate everything from operator()...
1556  while (__f != __t)
1557  *__f++ = this->operator()(__urng, __param);
1558  }
1559 
1560  template<typename _IntType,
1561  typename _CharT, typename _Traits>
1562  std::basic_ostream<_CharT, _Traits>&
1563  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1564  const poisson_distribution<_IntType>& __x)
1565  {
1566  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1567 
1568  const typename __ios_base::fmtflags __flags = __os.flags();
1569  const _CharT __fill = __os.fill();
1570  const std::streamsize __precision = __os.precision();
1571  const _CharT __space = __os.widen(' ');
1572  __os.flags(__ios_base::scientific | __ios_base::left);
1573  __os.fill(__space);
1574  __os.precision(std::numeric_limits<double>::max_digits10);
1575 
1576  __os << __x.mean() << __space << __x._M_nd;
1577 
1578  __os.flags(__flags);
1579  __os.fill(__fill);
1580  __os.precision(__precision);
1581  return __os;
1582  }
1583 
1584  template<typename _IntType,
1585  typename _CharT, typename _Traits>
1586  std::basic_istream<_CharT, _Traits>&
1587  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1588  poisson_distribution<_IntType>& __x)
1589  {
1590  using param_type = typename poisson_distribution<_IntType>::param_type;
1591  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1592 
1593  const typename __ios_base::fmtflags __flags = __is.flags();
1594  __is.flags(__ios_base::skipws);
1595 
1596  double __mean;
1597  if (__is >> __mean >> __x._M_nd)
1598  __x.param(param_type(__mean));
1599 
1600  __is.flags(__flags);
1601  return __is;
1602  }
1603 
1604 
1605  template<typename _IntType>
1606  void
1607  binomial_distribution<_IntType>::param_type::
1608  _M_initialize()
1609  {
1610  const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1611 
1612  _M_easy = true;
1613 
1614 #if _GLIBCXX_USE_C99_MATH_FUNCS
1615  if (_M_t * __p12 >= 8)
1616  {
1617  _M_easy = false;
1618  const double __np = std::floor(_M_t * __p12);
1619  const double __pa = __np / _M_t;
1620  const double __1p = 1 - __pa;
1621 
1622  const double __pi_4 = 0.7853981633974483096156608458198757L;
1623  const double __d1x =
1624  std::sqrt(__np * __1p * std::log(32 * __np
1625  / (81 * __pi_4 * __1p)));
1626  _M_d1 = std::round(std::max<double>(1.0, __d1x));
1627  const double __d2x =
1628  std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1629  / (__pi_4 * __pa)));
1630  _M_d2 = std::round(std::max<double>(1.0, __d2x));
1631 
1632  // sqrt(pi / 2)
1633  const double __spi_2 = 1.2533141373155002512078826424055226L;
1634  _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1635  _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * (_M_t * __1p)));
1636  _M_c = 2 * _M_d1 / __np;
1637  _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1638  const double __a12 = _M_a1 + _M_s2 * __spi_2;
1639  const double __s1s = _M_s1 * _M_s1;
1640  _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1641  * 2 * __s1s / _M_d1
1642  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1643  const double __s2s = _M_s2 * _M_s2;
1644  _M_s = (_M_a123 + 2 * __s2s / _M_d2
1645  * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1646  _M_lf = (std::lgamma(__np + 1)
1647  + std::lgamma(_M_t - __np + 1));
1648  _M_lp1p = std::log(__pa / __1p);
1649 
1650  _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1651  }
1652  else
1653 #endif
1654  _M_q = -std::log(1 - __p12);
1655  }
1656 
1657  template<typename _IntType>
1658  template<typename _UniformRandomNumberGenerator>
1659  typename binomial_distribution<_IntType>::result_type
1660  binomial_distribution<_IntType>::
1661  _M_waiting(_UniformRandomNumberGenerator& __urng,
1662  _IntType __t, double __q)
1663  {
1664  _IntType __x = 0;
1665  double __sum = 0.0;
1666  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1667  __aurng(__urng);
1668 
1669  do
1670  {
1671  if (__t == __x)
1672  return __x;
1673  const double __e = -std::log(1.0 - __aurng());
1674  __sum += __e / (__t - __x);
1675  __x += 1;
1676  }
1677  while (__sum <= __q);
1678 
1679  return __x - 1;
1680  }
1681 
1682  /**
1683  * A rejection algorithm when t * p >= 8 and a simple waiting time
1684  * method - the second in the referenced book - otherwise.
1685  * NB: The former is available only if _GLIBCXX_USE_C99_MATH_FUNCS
1686  * is defined.
1687  *
1688  * Reference:
1689  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1690  * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1691  */
1692  template<typename _IntType>
1693  template<typename _UniformRandomNumberGenerator>
1694  typename binomial_distribution<_IntType>::result_type
1695  binomial_distribution<_IntType>::
1696  operator()(_UniformRandomNumberGenerator& __urng,
1697  const param_type& __param)
1698  {
1699  result_type __ret;
1700  const _IntType __t = __param.t();
1701  const double __p = __param.p();
1702  const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1703  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1704  __aurng(__urng);
1705 
1706 #if _GLIBCXX_USE_C99_MATH_FUNCS
1707  if (!__param._M_easy)
1708  {
1709  double __x;
1710 
1711  // See comments above...
1712  const double __naf =
1713  (1 - std::numeric_limits<double>::epsilon()) / 2;
1714  const double __thr =
1715  std::numeric_limits<_IntType>::max() + __naf;
1716 
1717  const double __np = std::floor(__t * __p12);
1718 
1719  // sqrt(pi / 2)
1720  const double __spi_2 = 1.2533141373155002512078826424055226L;
1721  const double __a1 = __param._M_a1;
1722  const double __a12 = __a1 + __param._M_s2 * __spi_2;
1723  const double __a123 = __param._M_a123;
1724  const double __s1s = __param._M_s1 * __param._M_s1;
1725  const double __s2s = __param._M_s2 * __param._M_s2;
1726 
1727  bool __reject;
1728  do
1729  {
1730  const double __u = __param._M_s * __aurng();
1731 
1732  double __v;
1733 
1734  if (__u <= __a1)
1735  {
1736  const double __n = _M_nd(__urng);
1737  const double __y = __param._M_s1 * std::abs(__n);
1738  __reject = __y >= __param._M_d1;
1739  if (!__reject)
1740  {
1741  const double __e = -std::log(1.0 - __aurng());
1742  __x = std::floor(__y);
1743  __v = -__e - __n * __n / 2 + __param._M_c;
1744  }
1745  }
1746  else if (__u <= __a12)
1747  {
1748  const double __n = _M_nd(__urng);
1749  const double __y = __param._M_s2 * std::abs(__n);
1750  __reject = __y >= __param._M_d2;
1751  if (!__reject)
1752  {
1753  const double __e = -std::log(1.0 - __aurng());
1754  __x = std::floor(-__y);
1755  __v = -__e - __n * __n / 2;
1756  }
1757  }
1758  else if (__u <= __a123)
1759  {
1760  const double __e1 = -std::log(1.0 - __aurng());
1761  const double __e2 = -std::log(1.0 - __aurng());
1762 
1763  const double __y = __param._M_d1
1764  + 2 * __s1s * __e1 / __param._M_d1;
1765  __x = std::floor(__y);
1766  __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1767  -__y / (2 * __s1s)));
1768  __reject = false;
1769  }
1770  else
1771  {
1772  const double __e1 = -std::log(1.0 - __aurng());
1773  const double __e2 = -std::log(1.0 - __aurng());
1774 
1775  const double __y = __param._M_d2
1776  + 2 * __s2s * __e1 / __param._M_d2;
1777  __x = std::floor(-__y);
1778  __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1779  __reject = false;
1780  }
1781 
1782  __reject = __reject || __x < -__np || __x > __t - __np;
1783  if (!__reject)
1784  {
1785  const double __lfx =
1786  std::lgamma(__np + __x + 1)
1787  + std::lgamma(__t - (__np + __x) + 1);
1788  __reject = __v > __param._M_lf - __lfx
1789  + __x * __param._M_lp1p;
1790  }
1791 
1792  __reject |= __x + __np >= __thr;
1793  }
1794  while (__reject);
1795 
1796  __x += __np + __naf;
1797 
1798  const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1799  __param._M_q);
1800  __ret = _IntType(__x) + __z;
1801  }
1802  else
1803 #endif
1804  __ret = _M_waiting(__urng, __t, __param._M_q);
1805 
1806  if (__p12 != __p)
1807  __ret = __t - __ret;
1808  return __ret;
1809  }
1810 
1811  template<typename _IntType>
1812  template<typename _ForwardIterator,
1813  typename _UniformRandomNumberGenerator>
1814  void
1815  binomial_distribution<_IntType>::
1816  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1817  _UniformRandomNumberGenerator& __urng,
1818  const param_type& __param)
1819  {
1820  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1821  // We could duplicate everything from operator()...
1822  while (__f != __t)
1823  *__f++ = this->operator()(__urng, __param);
1824  }
1825 
1826  template<typename _IntType,
1827  typename _CharT, typename _Traits>
1828  std::basic_ostream<_CharT, _Traits>&
1829  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1830  const binomial_distribution<_IntType>& __x)
1831  {
1832  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1833 
1834  const typename __ios_base::fmtflags __flags = __os.flags();
1835  const _CharT __fill = __os.fill();
1836  const std::streamsize __precision = __os.precision();
1837  const _CharT __space = __os.widen(' ');
1838  __os.flags(__ios_base::scientific | __ios_base::left);
1839  __os.fill(__space);
1840  __os.precision(std::numeric_limits<double>::max_digits10);
1841 
1842  __os << __x.t() << __space << __x.p()
1843  << __space << __x._M_nd;
1844 
1845  __os.flags(__flags);
1846  __os.fill(__fill);
1847  __os.precision(__precision);
1848  return __os;
1849  }
1850 
1851  template<typename _IntType,
1852  typename _CharT, typename _Traits>
1853  std::basic_istream<_CharT, _Traits>&
1854  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1855  binomial_distribution<_IntType>& __x)
1856  {
1857  using param_type = typename binomial_distribution<_IntType>::param_type;
1858  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1859 
1860  const typename __ios_base::fmtflags __flags = __is.flags();
1861  __is.flags(__ios_base::dec | __ios_base::skipws);
1862 
1863  _IntType __t;
1864  double __p;
1865  if (__is >> __t >> __p >> __x._M_nd)
1866  __x.param(param_type(__t, __p));
1867 
1868  __is.flags(__flags);
1869  return __is;
1870  }
1871 
1872 
1873  template<typename _RealType>
1874  template<typename _ForwardIterator,
1875  typename _UniformRandomNumberGenerator>
1876  void
1877  std::exponential_distribution<_RealType>::
1878  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1879  _UniformRandomNumberGenerator& __urng,
1880  const param_type& __p)
1881  {
1882  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1883  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1884  __aurng(__urng);
1885  while (__f != __t)
1886  *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1887  }
1888 
1889  template<typename _RealType, typename _CharT, typename _Traits>
1890  std::basic_ostream<_CharT, _Traits>&
1891  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1892  const exponential_distribution<_RealType>& __x)
1893  {
1894  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1895 
1896  const typename __ios_base::fmtflags __flags = __os.flags();
1897  const _CharT __fill = __os.fill();
1898  const std::streamsize __precision = __os.precision();
1899  __os.flags(__ios_base::scientific | __ios_base::left);
1900  __os.fill(__os.widen(' '));
1901  __os.precision(std::numeric_limits<_RealType>::max_digits10);
1902 
1903  __os << __x.lambda();
1904 
1905  __os.flags(__flags);
1906  __os.fill(__fill);
1907  __os.precision(__precision);
1908  return __os;
1909  }
1910 
1911  template<typename _RealType, typename _CharT, typename _Traits>
1912  std::basic_istream<_CharT, _Traits>&
1913  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1914  exponential_distribution<_RealType>& __x)
1915  {
1916  using param_type
1917  = typename exponential_distribution<_RealType>::param_type;
1918  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1919 
1920  const typename __ios_base::fmtflags __flags = __is.flags();
1921  __is.flags(__ios_base::dec | __ios_base::skipws);
1922 
1923  _RealType __lambda;
1924  if (__is >> __lambda)
1925  __x.param(param_type(__lambda));
1926 
1927  __is.flags(__flags);
1928  return __is;
1929  }
1930 
1931 
1932  /**
1933  * Polar method due to Marsaglia.
1934  *
1935  * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1936  * New York, 1986, Ch. V, Sect. 4.4.
1937  */
1938  template<typename _RealType>
1939  template<typename _UniformRandomNumberGenerator>
1940  typename normal_distribution<_RealType>::result_type
1941  normal_distribution<_RealType>::
1942  operator()(_UniformRandomNumberGenerator& __urng,
1943  const param_type& __param)
1944  {
1945  result_type __ret;
1946  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1947  __aurng(__urng);
1948 
1949  if (_M_saved_available)
1950  {
1951  _M_saved_available = false;
1952  __ret = _M_saved;
1953  }
1954  else
1955  {
1956  result_type __x, __y, __r2;
1957  do
1958  {
1959  __x = result_type(2.0) * __aurng() - 1.0;
1960  __y = result_type(2.0) * __aurng() - 1.0;
1961  __r2 = __x * __x + __y * __y;
1962  }
1963  while (__r2 > 1.0 || __r2 == 0.0);
1964 
1965  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1966  _M_saved = __x * __mult;
1967  _M_saved_available = true;
1968  __ret = __y * __mult;
1969  }
1970 
1971  __ret = __ret * __param.stddev() + __param.mean();
1972  return __ret;
1973  }
1974 
1975  template<typename _RealType>
1976  template<typename _ForwardIterator,
1977  typename _UniformRandomNumberGenerator>
1978  void
1979  normal_distribution<_RealType>::
1980  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1981  _UniformRandomNumberGenerator& __urng,
1982  const param_type& __param)
1983  {
1984  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1985 
1986  if (__f == __t)
1987  return;
1988 
1989  if (_M_saved_available)
1990  {
1991  _M_saved_available = false;
1992  *__f++ = _M_saved * __param.stddev() + __param.mean();
1993 
1994  if (__f == __t)
1995  return;
1996  }
1997 
1998  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1999  __aurng(__urng);
2000 
2001  while (__f + 1 < __t)
2002  {
2003  result_type __x, __y, __r2;
2004  do
2005  {
2006  __x = result_type(2.0) * __aurng() - 1.0;
2007  __y = result_type(2.0) * __aurng() - 1.0;
2008  __r2 = __x * __x + __y * __y;
2009  }
2010  while (__r2 > 1.0 || __r2 == 0.0);
2011 
2012  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2013  *__f++ = __y * __mult * __param.stddev() + __param.mean();
2014  *__f++ = __x * __mult * __param.stddev() + __param.mean();
2015  }
2016 
2017  if (__f != __t)
2018  {
2019  result_type __x, __y, __r2;
2020  do
2021  {
2022  __x = result_type(2.0) * __aurng() - 1.0;
2023  __y = result_type(2.0) * __aurng() - 1.0;
2024  __r2 = __x * __x + __y * __y;
2025  }
2026  while (__r2 > 1.0 || __r2 == 0.0);
2027 
2028  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2029  _M_saved = __x * __mult;
2030  _M_saved_available = true;
2031  *__f = __y * __mult * __param.stddev() + __param.mean();
2032  }
2033  }
2034 
2035  template<typename _RealType>
2036  bool
2037  operator==(const std::normal_distribution<_RealType>& __d1,
2038  const std::normal_distribution<_RealType>& __d2)
2039  {
2040  if (__d1._M_param == __d2._M_param
2041  && __d1._M_saved_available == __d2._M_saved_available)
2042  return __d1._M_saved_available ? __d1._M_saved == __d2._M_saved : true;
2043  else
2044  return false;
2045  }
2046 
2047  template<typename _RealType, typename _CharT, typename _Traits>
2048  std::basic_ostream<_CharT, _Traits>&
2049  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2050  const normal_distribution<_RealType>& __x)
2051  {
2052  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2053 
2054  const typename __ios_base::fmtflags __flags = __os.flags();
2055  const _CharT __fill = __os.fill();
2056  const std::streamsize __precision = __os.precision();
2057  const _CharT __space = __os.widen(' ');
2058  __os.flags(__ios_base::scientific | __ios_base::left);
2059  __os.fill(__space);
2060  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2061 
2062  __os << __x.mean() << __space << __x.stddev()
2063  << __space << __x._M_saved_available;
2064  if (__x._M_saved_available)
2065  __os << __space << __x._M_saved;
2066 
2067  __os.flags(__flags);
2068  __os.fill(__fill);
2069  __os.precision(__precision);
2070  return __os;
2071  }
2072 
2073  template<typename _RealType, typename _CharT, typename _Traits>
2074  std::basic_istream<_CharT, _Traits>&
2075  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2076  normal_distribution<_RealType>& __x)
2077  {
2078  using param_type = typename normal_distribution<_RealType>::param_type;
2079  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2080 
2081  const typename __ios_base::fmtflags __flags = __is.flags();
2082  __is.flags(__ios_base::dec | __ios_base::skipws);
2083 
2084  double __mean, __stddev;
2085  bool __saved_avail;
2086  if (__is >> __mean >> __stddev >> __saved_avail)
2087  {
2088  if (!__saved_avail || (__is >> __x._M_saved))
2089  {
2090  __x._M_saved_available = __saved_avail;
2091  __x.param(param_type(__mean, __stddev));
2092  }
2093  }
2094 
2095  __is.flags(__flags);
2096  return __is;
2097  }
2098 
2099 
2100  template<typename _RealType>
2101  template<typename _ForwardIterator,
2102  typename _UniformRandomNumberGenerator>
2103  void
2104  lognormal_distribution<_RealType>::
2105  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2106  _UniformRandomNumberGenerator& __urng,
2107  const param_type& __p)
2108  {
2109  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2110  while (__f != __t)
2111  *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
2112  }
2113 
2114  template<typename _RealType, typename _CharT, typename _Traits>
2115  std::basic_ostream<_CharT, _Traits>&
2116  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2117  const lognormal_distribution<_RealType>& __x)
2118  {
2119  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2120 
2121  const typename __ios_base::fmtflags __flags = __os.flags();
2122  const _CharT __fill = __os.fill();
2123  const std::streamsize __precision = __os.precision();
2124  const _CharT __space = __os.widen(' ');
2125  __os.flags(__ios_base::scientific | __ios_base::left);
2126  __os.fill(__space);
2127  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2128 
2129  __os << __x.m() << __space << __x.s()
2130  << __space << __x._M_nd;
2131 
2132  __os.flags(__flags);
2133  __os.fill(__fill);
2134  __os.precision(__precision);
2135  return __os;
2136  }
2137 
2138  template<typename _RealType, typename _CharT, typename _Traits>
2139  std::basic_istream<_CharT, _Traits>&
2140  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2141  lognormal_distribution<_RealType>& __x)
2142  {
2143  using param_type
2144  = typename lognormal_distribution<_RealType>::param_type;
2145  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2146 
2147  const typename __ios_base::fmtflags __flags = __is.flags();
2148  __is.flags(__ios_base::dec | __ios_base::skipws);
2149 
2150  _RealType __m, __s;
2151  if (__is >> __m >> __s >> __x._M_nd)
2152  __x.param(param_type(__m, __s));
2153 
2154  __is.flags(__flags);
2155  return __is;
2156  }
2157 
2158  template<typename _RealType>
2159  template<typename _ForwardIterator,
2160  typename _UniformRandomNumberGenerator>
2161  void
2162  std::chi_squared_distribution<_RealType>::
2163  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2164  _UniformRandomNumberGenerator& __urng)
2165  {
2166  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2167  while (__f != __t)
2168  *__f++ = 2 * _M_gd(__urng);
2169  }
2170 
2171  template<typename _RealType>
2172  template<typename _ForwardIterator,
2173  typename _UniformRandomNumberGenerator>
2174  void
2175  std::chi_squared_distribution<_RealType>::
2176  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2177  _UniformRandomNumberGenerator& __urng,
2178  const typename
2179  std::gamma_distribution<result_type>::param_type& __p)
2180  {
2181  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2182  while (__f != __t)
2183  *__f++ = 2 * _M_gd(__urng, __p);
2184  }
2185 
2186  template<typename _RealType, typename _CharT, typename _Traits>
2187  std::basic_ostream<_CharT, _Traits>&
2188  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2189  const chi_squared_distribution<_RealType>& __x)
2190  {
2191  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2192 
2193  const typename __ios_base::fmtflags __flags = __os.flags();
2194  const _CharT __fill = __os.fill();
2195  const std::streamsize __precision = __os.precision();
2196  const _CharT __space = __os.widen(' ');
2197  __os.flags(__ios_base::scientific | __ios_base::left);
2198  __os.fill(__space);
2199  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2200 
2201  __os << __x.n() << __space << __x._M_gd;
2202 
2203  __os.flags(__flags);
2204  __os.fill(__fill);
2205  __os.precision(__precision);
2206  return __os;
2207  }
2208 
2209  template<typename _RealType, typename _CharT, typename _Traits>
2210  std::basic_istream<_CharT, _Traits>&
2211  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2212  chi_squared_distribution<_RealType>& __x)
2213  {
2214  using param_type
2215  = typename chi_squared_distribution<_RealType>::param_type;
2216  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2217 
2218  const typename __ios_base::fmtflags __flags = __is.flags();
2219  __is.flags(__ios_base::dec | __ios_base::skipws);
2220 
2221  _RealType __n;
2222  if (__is >> __n >> __x._M_gd)
2223  __x.param(param_type(__n));
2224 
2225  __is.flags(__flags);
2226  return __is;
2227  }
2228 
2229 
2230  template<typename _RealType>
2231  template<typename _UniformRandomNumberGenerator>
2232  typename cauchy_distribution<_RealType>::result_type
2233  cauchy_distribution<_RealType>::
2234  operator()(_UniformRandomNumberGenerator& __urng,
2235  const param_type& __p)
2236  {
2237  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2238  __aurng(__urng);
2239  _RealType __u;
2240  do
2241  __u = __aurng();
2242  while (__u == 0.5);
2243 
2244  const _RealType __pi = 3.1415926535897932384626433832795029L;
2245  return __p.a() + __p.b() * std::tan(__pi * __u);
2246  }
2247 
2248  template<typename _RealType>
2249  template<typename _ForwardIterator,
2250  typename _UniformRandomNumberGenerator>
2251  void
2252  cauchy_distribution<_RealType>::
2253  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2254  _UniformRandomNumberGenerator& __urng,
2255  const param_type& __p)
2256  {
2257  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2258  const _RealType __pi = 3.1415926535897932384626433832795029L;
2259  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2260  __aurng(__urng);
2261  while (__f != __t)
2262  {
2263  _RealType __u;
2264  do
2265  __u = __aurng();
2266  while (__u == 0.5);
2267 
2268  *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2269  }
2270  }
2271 
2272  template<typename _RealType, typename _CharT, typename _Traits>
2273  std::basic_ostream<_CharT, _Traits>&
2274  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2275  const cauchy_distribution<_RealType>& __x)
2276  {
2277  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2278 
2279  const typename __ios_base::fmtflags __flags = __os.flags();
2280  const _CharT __fill = __os.fill();
2281  const std::streamsize __precision = __os.precision();
2282  const _CharT __space = __os.widen(' ');
2283  __os.flags(__ios_base::scientific | __ios_base::left);
2284  __os.fill(__space);
2285  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2286 
2287  __os << __x.a() << __space << __x.b();
2288 
2289  __os.flags(__flags);
2290  __os.fill(__fill);
2291  __os.precision(__precision);
2292  return __os;
2293  }
2294 
2295  template<typename _RealType, typename _CharT, typename _Traits>
2296  std::basic_istream<_CharT, _Traits>&
2297  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2298  cauchy_distribution<_RealType>& __x)
2299  {
2300  using param_type = typename cauchy_distribution<_RealType>::param_type;
2301  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2302 
2303  const typename __ios_base::fmtflags __flags = __is.flags();
2304  __is.flags(__ios_base::dec | __ios_base::skipws);
2305 
2306  _RealType __a, __b;
2307  if (__is >> __a >> __b)
2308  __x.param(param_type(__a, __b));
2309 
2310  __is.flags(__flags);
2311  return __is;
2312  }
2313 
2314 
2315  template<typename _RealType>
2316  template<typename _ForwardIterator,
2317  typename _UniformRandomNumberGenerator>
2318  void
2319  std::fisher_f_distribution<_RealType>::
2320  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2321  _UniformRandomNumberGenerator& __urng)
2322  {
2323  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2324  while (__f != __t)
2325  *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2326  }
2327 
2328  template<typename _RealType>
2329  template<typename _ForwardIterator,
2330  typename _UniformRandomNumberGenerator>
2331  void
2332  std::fisher_f_distribution<_RealType>::
2333  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2334  _UniformRandomNumberGenerator& __urng,
2335  const param_type& __p)
2336  {
2337  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2338  typedef typename std::gamma_distribution<result_type>::param_type
2339  param_type;
2340  param_type __p1(__p.m() / 2);
2341  param_type __p2(__p.n() / 2);
2342  while (__f != __t)
2343  *__f++ = ((_M_gd_x(__urng, __p1) * n())
2344  / (_M_gd_y(__urng, __p2) * m()));
2345  }
2346 
2347  template<typename _RealType, typename _CharT, typename _Traits>
2348  std::basic_ostream<_CharT, _Traits>&
2349  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2350  const fisher_f_distribution<_RealType>& __x)
2351  {
2352  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2353 
2354  const typename __ios_base::fmtflags __flags = __os.flags();
2355  const _CharT __fill = __os.fill();
2356  const std::streamsize __precision = __os.precision();
2357  const _CharT __space = __os.widen(' ');
2358  __os.flags(__ios_base::scientific | __ios_base::left);
2359  __os.fill(__space);
2360  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2361 
2362  __os << __x.m() << __space << __x.n()
2363  << __space << __x._M_gd_x << __space << __x._M_gd_y;
2364 
2365  __os.flags(__flags);
2366  __os.fill(__fill);
2367  __os.precision(__precision);
2368  return __os;
2369  }
2370 
2371  template<typename _RealType, typename _CharT, typename _Traits>
2372  std::basic_istream<_CharT, _Traits>&
2373  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2374  fisher_f_distribution<_RealType>& __x)
2375  {
2376  using param_type
2377  = typename fisher_f_distribution<_RealType>::param_type;
2378  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2379 
2380  const typename __ios_base::fmtflags __flags = __is.flags();
2381  __is.flags(__ios_base::dec | __ios_base::skipws);
2382 
2383  _RealType __m, __n;
2384  if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y)
2385  __x.param(param_type(__m, __n));
2386 
2387  __is.flags(__flags);
2388  return __is;
2389  }
2390 
2391 
2392  template<typename _RealType>
2393  template<typename _ForwardIterator,
2394  typename _UniformRandomNumberGenerator>
2395  void
2396  std::student_t_distribution<_RealType>::
2397  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2398  _UniformRandomNumberGenerator& __urng)
2399  {
2400  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2401  while (__f != __t)
2402  *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2403  }
2404 
2405  template<typename _RealType>
2406  template<typename _ForwardIterator,
2407  typename _UniformRandomNumberGenerator>
2408  void
2409  std::student_t_distribution<_RealType>::
2410  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2411  _UniformRandomNumberGenerator& __urng,
2412  const param_type& __p)
2413  {
2414  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2415  typename std::gamma_distribution<result_type>::param_type
2416  __p2(__p.n() / 2, 2);
2417  while (__f != __t)
2418  *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2419  }
2420 
2421  template<typename _RealType, typename _CharT, typename _Traits>
2422  std::basic_ostream<_CharT, _Traits>&
2423  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2424  const student_t_distribution<_RealType>& __x)
2425  {
2426  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2427 
2428  const typename __ios_base::fmtflags __flags = __os.flags();
2429  const _CharT __fill = __os.fill();
2430  const std::streamsize __precision = __os.precision();
2431  const _CharT __space = __os.widen(' ');
2432  __os.flags(__ios_base::scientific | __ios_base::left);
2433  __os.fill(__space);
2434  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2435 
2436  __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2437 
2438  __os.flags(__flags);
2439  __os.fill(__fill);
2440  __os.precision(__precision);
2441  return __os;
2442  }
2443 
2444  template<typename _RealType, typename _CharT, typename _Traits>
2445  std::basic_istream<_CharT, _Traits>&
2446  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2447  student_t_distribution<_RealType>& __x)
2448  {
2449  using param_type
2450  = typename student_t_distribution<_RealType>::param_type;
2451  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2452 
2453  const typename __ios_base::fmtflags __flags = __is.flags();
2454  __is.flags(__ios_base::dec | __ios_base::skipws);
2455 
2456  _RealType __n;
2457  if (__is >> __n >> __x._M_nd >> __x._M_gd)
2458  __x.param(param_type(__n));
2459 
2460  __is.flags(__flags);
2461  return __is;
2462  }
2463 
2464 
2465  template<typename _RealType>
2466  void
2467  gamma_distribution<_RealType>::param_type::
2468  _M_initialize()
2469  {
2470  _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2471 
2472  const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2473  _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2474  }
2475 
2476  /**
2477  * Marsaglia, G. and Tsang, W. W.
2478  * "A Simple Method for Generating Gamma Variables"
2479  * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2480  */
2481  template<typename _RealType>
2482  template<typename _UniformRandomNumberGenerator>
2483  typename gamma_distribution<_RealType>::result_type
2484  gamma_distribution<_RealType>::
2485  operator()(_UniformRandomNumberGenerator& __urng,
2486  const param_type& __param)
2487  {
2488  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2489  __aurng(__urng);
2490 
2491  result_type __u, __v, __n;
2492  const result_type __a1 = (__param._M_malpha
2493  - _RealType(1.0) / _RealType(3.0));
2494 
2495  do
2496  {
2497  do
2498  {
2499  __n = _M_nd(__urng);
2500  __v = result_type(1.0) + __param._M_a2 * __n;
2501  }
2502  while (__v <= 0.0);
2503 
2504  __v = __v * __v * __v;
2505  __u = __aurng();
2506  }
2507  while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2508  && (std::log(__u) > (0.5 * __n * __n + __a1
2509  * (1.0 - __v + std::log(__v)))));
2510 
2511  if (__param.alpha() == __param._M_malpha)
2512  return __a1 * __v * __param.beta();
2513  else
2514  {
2515  do
2516  __u = __aurng();
2517  while (__u == 0.0);
2518 
2519  return (std::pow(__u, result_type(1.0) / __param.alpha())
2520  * __a1 * __v * __param.beta());
2521  }
2522  }
2523 
2524  template<typename _RealType>
2525  template<typename _ForwardIterator,
2526  typename _UniformRandomNumberGenerator>
2527  void
2528  gamma_distribution<_RealType>::
2529  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2530  _UniformRandomNumberGenerator& __urng,
2531  const param_type& __param)
2532  {
2533  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2534  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2535  __aurng(__urng);
2536 
2537  result_type __u, __v, __n;
2538  const result_type __a1 = (__param._M_malpha
2539  - _RealType(1.0) / _RealType(3.0));
2540 
2541  if (__param.alpha() == __param._M_malpha)
2542  while (__f != __t)
2543  {
2544  do
2545  {
2546  do
2547  {
2548  __n = _M_nd(__urng);
2549  __v = result_type(1.0) + __param._M_a2 * __n;
2550  }
2551  while (__v <= 0.0);
2552 
2553  __v = __v * __v * __v;
2554  __u = __aurng();
2555  }
2556  while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2557  && (std::log(__u) > (0.5 * __n * __n + __a1
2558  * (1.0 - __v + std::log(__v)))));
2559 
2560  *__f++ = __a1 * __v * __param.beta();
2561  }
2562  else
2563  while (__f != __t)
2564  {
2565  do
2566  {
2567  do
2568  {
2569  __n = _M_nd(__urng);
2570  __v = result_type(1.0) + __param._M_a2 * __n;
2571  }
2572  while (__v <= 0.0);
2573 
2574  __v = __v * __v * __v;
2575  __u = __aurng();
2576  }
2577  while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2578  && (std::log(__u) > (0.5 * __n * __n + __a1
2579  * (1.0 - __v + std::log(__v)))));
2580 
2581  do
2582  __u = __aurng();
2583  while (__u == 0.0);
2584 
2585  *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2586  * __a1 * __v * __param.beta());
2587  }
2588  }
2589 
2590  template<typename _RealType, typename _CharT, typename _Traits>
2591  std::basic_ostream<_CharT, _Traits>&
2592  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2593  const gamma_distribution<_RealType>& __x)
2594  {
2595  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2596 
2597  const typename __ios_base::fmtflags __flags = __os.flags();
2598  const _CharT __fill = __os.fill();
2599  const std::streamsize __precision = __os.precision();
2600  const _CharT __space = __os.widen(' ');
2601  __os.flags(__ios_base::scientific | __ios_base::left);
2602  __os.fill(__space);
2603  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2604 
2605  __os << __x.alpha() << __space << __x.beta()
2606  << __space << __x._M_nd;
2607 
2608  __os.flags(__flags);
2609  __os.fill(__fill);
2610  __os.precision(__precision);
2611  return __os;
2612  }
2613 
2614  template<typename _RealType, typename _CharT, typename _Traits>
2615  std::basic_istream<_CharT, _Traits>&
2616  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2617  gamma_distribution<_RealType>& __x)
2618  {
2619  using param_type = typename gamma_distribution<_RealType>::param_type;
2620  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2621 
2622  const typename __ios_base::fmtflags __flags = __is.flags();
2623  __is.flags(__ios_base::dec | __ios_base::skipws);
2624 
2625  _RealType __alpha_val, __beta_val;
2626  if (__is >> __alpha_val >> __beta_val >> __x._M_nd)
2627  __x.param(param_type(__alpha_val, __beta_val));
2628 
2629  __is.flags(__flags);
2630  return __is;
2631  }
2632 
2633 
2634  template<typename _RealType>
2635  template<typename _UniformRandomNumberGenerator>
2636  typename weibull_distribution<_RealType>::result_type
2637  weibull_distribution<_RealType>::
2638  operator()(_UniformRandomNumberGenerator& __urng,
2639  const param_type& __p)
2640  {
2641  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2642  __aurng(__urng);
2643  return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2644  result_type(1) / __p.a());
2645  }
2646 
2647  template<typename _RealType>
2648  template<typename _ForwardIterator,
2649  typename _UniformRandomNumberGenerator>
2650  void
2651  weibull_distribution<_RealType>::
2652  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2653  _UniformRandomNumberGenerator& __urng,
2654  const param_type& __p)
2655  {
2656  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2657  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2658  __aurng(__urng);
2659  auto __inv_a = result_type(1) / __p.a();
2660 
2661  while (__f != __t)
2662  *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2663  __inv_a);
2664  }
2665 
2666  template<typename _RealType, typename _CharT, typename _Traits>
2667  std::basic_ostream<_CharT, _Traits>&
2668  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2669  const weibull_distribution<_RealType>& __x)
2670  {
2671  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2672 
2673  const typename __ios_base::fmtflags __flags = __os.flags();
2674  const _CharT __fill = __os.fill();
2675  const std::streamsize __precision = __os.precision();
2676  const _CharT __space = __os.widen(' ');
2677  __os.flags(__ios_base::scientific | __ios_base::left);
2678  __os.fill(__space);
2679  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2680 
2681  __os << __x.a() << __space << __x.b();
2682 
2683  __os.flags(__flags);
2684  __os.fill(__fill);
2685  __os.precision(__precision);
2686  return __os;
2687  }
2688 
2689  template<typename _RealType, typename _CharT, typename _Traits>
2690  std::basic_istream<_CharT, _Traits>&
2691  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2692  weibull_distribution<_RealType>& __x)
2693  {
2694  using param_type = typename weibull_distribution<_RealType>::param_type;
2695  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2696 
2697  const typename __ios_base::fmtflags __flags = __is.flags();
2698  __is.flags(__ios_base::dec | __ios_base::skipws);
2699 
2700  _RealType __a, __b;
2701  if (__is >> __a >> __b)
2702  __x.param(param_type(__a, __b));
2703 
2704  __is.flags(__flags);
2705  return __is;
2706  }
2707 
2708 
2709  template<typename _RealType>
2710  template<typename _UniformRandomNumberGenerator>
2711  typename extreme_value_distribution<_RealType>::result_type
2712  extreme_value_distribution<_RealType>::
2713  operator()(_UniformRandomNumberGenerator& __urng,
2714  const param_type& __p)
2715  {
2716  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2717  __aurng(__urng);
2718  return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2719  - __aurng()));
2720  }
2721 
2722  template<typename _RealType>
2723  template<typename _ForwardIterator,
2724  typename _UniformRandomNumberGenerator>
2725  void
2726  extreme_value_distribution<_RealType>::
2727  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2728  _UniformRandomNumberGenerator& __urng,
2729  const param_type& __p)
2730  {
2731  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2732  __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2733  __aurng(__urng);
2734 
2735  while (__f != __t)
2736  *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2737  - __aurng()));
2738  }
2739 
2740  template<typename _RealType, typename _CharT, typename _Traits>
2741  std::basic_ostream<_CharT, _Traits>&
2742  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2743  const extreme_value_distribution<_RealType>& __x)
2744  {
2745  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2746 
2747  const typename __ios_base::fmtflags __flags = __os.flags();
2748  const _CharT __fill = __os.fill();
2749  const std::streamsize __precision = __os.precision();
2750  const _CharT __space = __os.widen(' ');
2751  __os.flags(__ios_base::scientific | __ios_base::left);
2752  __os.fill(__space);
2753  __os.precision(std::numeric_limits<_RealType>::max_digits10);
2754 
2755  __os << __x.a() << __space << __x.b();
2756 
2757  __os.flags(__flags);
2758  __os.fill(__fill);
2759  __os.precision(__precision);
2760  return __os;
2761  }
2762 
2763  template<typename _RealType, typename _CharT, typename _Traits>
2764  std::basic_istream<_CharT, _Traits>&
2765  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2766  extreme_value_distribution<_RealType>& __x)
2767  {
2768  using param_type
2769  = typename extreme_value_distribution<_RealType>::param_type;
2770  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2771 
2772  const typename __ios_base::fmtflags __flags = __is.flags();
2773  __is.flags(__ios_base::dec | __ios_base::skipws);
2774 
2775  _RealType __a, __b;
2776  if (__is >> __a >> __b)
2777  __x.param(param_type(__a, __b));
2778 
2779  __is.flags(__flags);
2780  return __is;
2781  }
2782 
2783 
2784  template<typename _IntType>
2785  void
2786  discrete_distribution<_IntType>::param_type::
2787  _M_initialize()
2788  {
2789  if (_M_prob.size() < 2)
2790  {
2791  _M_prob.clear();
2792  return;
2793  }
2794 
2795  const double __sum = std::accumulate(_M_prob.begin(),
2796  _M_prob.end(), 0.0);
2797  __glibcxx_assert(__sum > 0);
2798  // Now normalize the probabilites.
2799  __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2800  __sum);
2801  // Accumulate partial sums.
2802  _M_cp.reserve(_M_prob.size());
2803  std::partial_sum(_M_prob.begin(), _M_prob.end(),
2804  std::back_inserter(_M_cp));
2805  // Make sure the last cumulative probability is one.
2806  _M_cp[_M_cp.size() - 1] = 1.0;
2807  }
2808 
2809  template<typename _IntType>
2810  template<typename _Func>
2811  discrete_distribution<_IntType>::param_type::
2812  param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2813  : _M_prob(), _M_cp()
2814  {
2815  const size_t __n = __nw == 0 ? 1 : __nw;
2816  const double __delta = (__xmax - __xmin) / __n;
2817 
2818  _M_prob.reserve(__n);
2819  for (size_t __k = 0; __k < __nw; ++__k)
2820  _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2821 
2822  _M_initialize();
2823  }
2824 
2825  template<typename _IntType>
2826  template<typename _UniformRandomNumberGenerator>
2827  typename discrete_distribution<_IntType>::result_type
2828  discrete_distribution<_IntType>::
2829  operator()(_UniformRandomNumberGenerator& __urng,
2830  const param_type& __param)
2831  {
2832  if (__param._M_cp.empty())
2833  return result_type(0);
2834 
2835  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2836  __aurng(__urng);
2837 
2838  const double __p = __aurng();
2839  auto __pos = std::lower_bound(__param._M_cp.begin(),
2840  __param._M_cp.end(), __p);
2841 
2842  return __pos - __param._M_cp.begin();
2843  }
2844 
2845  template<typename _IntType>
2846  template<typename _ForwardIterator,
2847  typename _UniformRandomNumberGenerator>
2848  void
2849  discrete_distribution<_IntType>::
2850  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2851  _UniformRandomNumberGenerator& __urng,
2852  const param_type& __param)
2853  {
2854  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2855 
2856  if (__param._M_cp.empty())
2857  {
2858  while (__f != __t)
2859  *__f++ = result_type(0);
2860  return;
2861  }
2862 
2863  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2864  __aurng(__urng);
2865 
2866  while (__f != __t)
2867  {
2868  const double __p = __aurng();
2869  auto __pos = std::lower_bound(__param._M_cp.begin(),
2870  __param._M_cp.end(), __p);
2871 
2872  *__f++ = __pos - __param._M_cp.begin();
2873  }
2874  }
2875 
2876  template<typename _IntType, typename _CharT, typename _Traits>
2877  std::basic_ostream<_CharT, _Traits>&
2878  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2879  const discrete_distribution<_IntType>& __x)
2880  {
2881  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2882 
2883  const typename __ios_base::fmtflags __flags = __os.flags();
2884  const _CharT __fill = __os.fill();
2885  const std::streamsize __precision = __os.precision();
2886  const _CharT __space = __os.widen(' ');
2887  __os.flags(__ios_base::scientific | __ios_base::left);
2888  __os.fill(__space);
2889  __os.precision(std::numeric_limits<double>::max_digits10);
2890 
2891  std::vector<double> __prob = __x.probabilities();
2892  __os << __prob.size();
2893  for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2894  __os << __space << *__dit;
2895 
2896  __os.flags(__flags);
2897  __os.fill(__fill);
2898  __os.precision(__precision);
2899  return __os;
2900  }
2901 
2902 namespace __detail
2903 {
2904  template<typename _ValT, typename _CharT, typename _Traits>
2905  basic_istream<_CharT, _Traits>&
2906  __extract_params(basic_istream<_CharT, _Traits>& __is,
2907  vector<_ValT>& __vals, size_t __n)
2908  {
2909  __vals.reserve(__n);
2910  while (__n--)
2911  {
2912  _ValT __val;
2913  if (__is >> __val)
2914  __vals.push_back(__val);
2915  else
2916  break;
2917  }
2918  return __is;
2919  }
2920 } // namespace __detail
2921 
2922  template<typename _IntType, typename _CharT, typename _Traits>
2923  std::basic_istream<_CharT, _Traits>&
2924  operator>>(std::basic_istream<_CharT, _Traits>& __is,
2925  discrete_distribution<_IntType>& __x)
2926  {
2927  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2928 
2929  const typename __ios_base::fmtflags __flags = __is.flags();
2930  __is.flags(__ios_base::dec | __ios_base::skipws);
2931 
2932  size_t __n;
2933  if (__is >> __n)
2934  {
2935  std::vector<double> __prob_vec;
2936  if (__detail::__extract_params(__is, __prob_vec, __n))
2937  __x.param({__prob_vec.begin(), __prob_vec.end()});
2938  }
2939 
2940  __is.flags(__flags);
2941  return __is;
2942  }
2943 
2944 
2945  template<typename _RealType>
2946  void
2947  piecewise_constant_distribution<_RealType>::param_type::
2948  _M_initialize()
2949  {
2950  if (_M_int.size() < 2
2951  || (_M_int.size() == 2
2952  && _M_int[0] == _RealType(0)
2953  && _M_int[1] == _RealType(1)))
2954  {
2955  _M_int.clear();
2956  _M_den.clear();
2957  return;
2958  }
2959 
2960  const double __sum = std::accumulate(_M_den.begin(),
2961  _M_den.end(), 0.0);
2962  __glibcxx_assert(__sum > 0);
2963 
2964  __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2965  __sum);
2966 
2967  _M_cp.reserve(_M_den.size());
2968  std::partial_sum(_M_den.begin(), _M_den.end(),
2969  std::back_inserter(_M_cp));
2970 
2971  // Make sure the last cumulative probability is one.
2972  _M_cp[_M_cp.size() - 1] = 1.0;
2973 
2974  for (size_t __k = 0; __k < _M_den.size(); ++__k)
2975  _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2976  }
2977 
2978  template<typename _RealType>
2979  template<typename _InputIteratorB, typename _InputIteratorW>
2980  piecewise_constant_distribution<_RealType>::param_type::
2981  param_type(_InputIteratorB __bbegin,
2982  _InputIteratorB __bend,
2983  _InputIteratorW __wbegin)
2984  : _M_int(), _M_den(), _M_cp()
2985  {
2986  if (__bbegin != __bend)
2987  {
2988  for (;;)
2989  {
2990  _M_int.push_back(*__bbegin);
2991  ++__bbegin;
2992  if (__bbegin == __bend)
2993  break;
2994 
2995  _M_den.push_back(*__wbegin);
2996  ++__wbegin;
2997  }
2998  }
2999 
3000  _M_initialize();
3001  }
3002 
3003  template<typename _RealType>
3004  template<typename _Func>
3005  piecewise_constant_distribution<_RealType>::param_type::
3006  param_type(initializer_list<_RealType> __bl, _Func __fw)
3007  : _M_int(), _M_den(), _M_cp()
3008  {
3009  _M_int.reserve(__bl.size());
3010  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3011  _M_int.push_back(*__biter);
3012 
3013  _M_den.reserve(_M_int.size() - 1);
3014  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3015  _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
3016 
3017  _M_initialize();
3018  }
3019 
3020  template<typename _RealType>
3021  template<typename _Func>
3022  piecewise_constant_distribution<_RealType>::param_type::
3023  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3024  : _M_int(), _M_den(), _M_cp()
3025  {
3026  const size_t __n = __nw == 0 ? 1 : __nw;
3027  const _RealType __delta = (__xmax - __xmin) / __n;
3028 
3029  _M_int.reserve(__n + 1);
3030  for (size_t __k = 0; __k <= __nw; ++__k)
3031  _M_int.push_back(__xmin + __k * __delta);
3032 
3033  _M_den.reserve(__n);
3034  for (size_t __k = 0; __k < __nw; ++__k)
3035  _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
3036 
3037  _M_initialize();
3038  }
3039 
3040  template<typename _RealType>
3041  template<typename _UniformRandomNumberGenerator>
3042  typename piecewise_constant_distribution<_RealType>::result_type
3043  piecewise_constant_distribution<_RealType>::
3044  operator()(_UniformRandomNumberGenerator& __urng,
3045  const param_type& __param)
3046  {
3047  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3048  __aurng(__urng);
3049 
3050  const double __p = __aurng();
3051  if (__param._M_cp.empty())
3052  return __p;
3053 
3054  auto __pos = std::lower_bound(__param._M_cp.begin(),
3055  __param._M_cp.end(), __p);
3056  const size_t __i = __pos - __param._M_cp.begin();
3057 
3058  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3059 
3060  return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
3061  }
3062 
3063  template<typename _RealType>
3064  template<typename _ForwardIterator,
3065  typename _UniformRandomNumberGenerator>
3066  void
3067  piecewise_constant_distribution<_RealType>::
3068  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3069  _UniformRandomNumberGenerator& __urng,
3070  const param_type& __param)
3071  {
3072  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3073  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3074  __aurng(__urng);
3075 
3076  if (__param._M_cp.empty())
3077  {
3078  while (__f != __t)
3079  *__f++ = __aurng();
3080  return;
3081  }
3082 
3083  while (__f != __t)
3084  {
3085  const double __p = __aurng();
3086 
3087  auto __pos = std::lower_bound(__param._M_cp.begin(),
3088  __param._M_cp.end(), __p);
3089  const size_t __i = __pos - __param._M_cp.begin();
3090 
3091  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3092 
3093  *__f++ = (__param._M_int[__i]
3094  + (__p - __pref) / __param._M_den[__i]);
3095  }
3096  }
3097 
3098  template<typename _RealType, typename _CharT, typename _Traits>
3099  std::basic_ostream<_CharT, _Traits>&
3100  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3101  const piecewise_constant_distribution<_RealType>& __x)
3102  {
3103  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3104 
3105  const typename __ios_base::fmtflags __flags = __os.flags();
3106  const _CharT __fill = __os.fill();
3107  const std::streamsize __precision = __os.precision();
3108  const _CharT __space = __os.widen(' ');
3109  __os.flags(__ios_base::scientific | __ios_base::left);
3110  __os.fill(__space);
3111  __os.precision(std::numeric_limits<_RealType>::max_digits10);
3112 
3113  std::vector<_RealType> __int = __x.intervals();
3114  __os << __int.size() - 1;
3115 
3116  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3117  __os << __space << *__xit;
3118 
3119  std::vector<double> __den = __x.densities();
3120  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3121  __os << __space << *__dit;
3122 
3123  __os.flags(__flags);
3124  __os.fill(__fill);
3125  __os.precision(__precision);
3126  return __os;
3127  }
3128 
3129  template<typename _RealType, typename _CharT, typename _Traits>
3130  std::basic_istream<_CharT, _Traits>&
3131  operator>>(std::basic_istream<_CharT, _Traits>& __is,
3132  piecewise_constant_distribution<_RealType>& __x)
3133  {
3134  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3135 
3136  const typename __ios_base::fmtflags __flags = __is.flags();
3137  __is.flags(__ios_base::dec | __ios_base::skipws);
3138 
3139  size_t __n;
3140  if (__is >> __n)
3141  {
3142  std::vector<_RealType> __int_vec;
3143  if (__detail::__extract_params(__is, __int_vec, __n + 1))
3144  {
3145  std::vector<double> __den_vec;
3146  if (__detail::__extract_params(__is, __den_vec, __n))
3147  {
3148  __x.param({ __int_vec.begin(), __int_vec.end(),
3149  __den_vec.begin() });
3150  }
3151  }
3152  }
3153 
3154  __is.flags(__flags);
3155  return __is;
3156  }
3157 
3158 
3159  template<typename _RealType>
3160  void
3161  piecewise_linear_distribution<_RealType>::param_type::
3162  _M_initialize()
3163  {
3164  if (_M_int.size() < 2
3165  || (_M_int.size() == 2
3166  && _M_int[0] == _RealType(0)
3167  && _M_int[1] == _RealType(1)
3168  && _M_den[0] == _M_den[1]))
3169  {
3170  _M_int.clear();
3171  _M_den.clear();
3172  return;
3173  }
3174 
3175  double __sum = 0.0;
3176  _M_cp.reserve(_M_int.size() - 1);
3177  _M_m.reserve(_M_int.size() - 1);
3178  for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3179  {
3180  const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3181  __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3182  _M_cp.push_back(__sum);
3183  _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3184  }
3185  __glibcxx_assert(__sum > 0);
3186 
3187  // Now normalize the densities...
3188  __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3189  __sum);
3190  // ... and partial sums...
3191  __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3192  // ... and slopes.
3193  __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3194 
3195  // Make sure the last cumulative probablility is one.
3196  _M_cp[_M_cp.size() - 1] = 1.0;
3197  }
3198 
3199  template<typename _RealType>
3200  template<typename _InputIteratorB, typename _InputIteratorW>
3201  piecewise_linear_distribution<_RealType>::param_type::
3202  param_type(_InputIteratorB __bbegin,
3203  _InputIteratorB __bend,
3204  _InputIteratorW __wbegin)
3205  : _M_int(), _M_den(), _M_cp(), _M_m()
3206  {
3207  for (; __bbegin != __bend; ++__bbegin, (void) ++__wbegin)
3208  {
3209  _M_int.push_back(*__bbegin);
3210  _M_den.push_back(*__wbegin);
3211  }
3212 
3213  _M_initialize();
3214  }
3215 
3216  template<typename _RealType>
3217  template<typename _Func>
3218  piecewise_linear_distribution<_RealType>::param_type::
3219  param_type(initializer_list<_RealType> __bl, _Func __fw)
3220  : _M_int(), _M_den(), _M_cp(), _M_m()
3221  {
3222  _M_int.reserve(__bl.size());
3223  _M_den.reserve(__bl.size());
3224  for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3225  {
3226  _M_int.push_back(*__biter);
3227  _M_den.push_back(__fw(*__biter));
3228  }
3229 
3230  _M_initialize();
3231  }
3232 
3233  template<typename _RealType>
3234  template<typename _Func>
3235  piecewise_linear_distribution<_RealType>::param_type::
3236  param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3237  : _M_int(), _M_den(), _M_cp(), _M_m()
3238  {
3239  const size_t __n = __nw == 0 ? 1 : __nw;
3240  const _RealType __delta = (__xmax - __xmin) / __n;
3241 
3242  _M_int.reserve(__n + 1);
3243  _M_den.reserve(__n + 1);
3244  for (size_t __k = 0; __k <= __nw; ++__k)
3245  {
3246  _M_int.push_back(__xmin + __k * __delta);
3247  _M_den.push_back(__fw(_M_int[__k] + __delta));
3248  }
3249 
3250  _M_initialize();
3251  }
3252 
3253  template<typename _RealType>
3254  template<typename _UniformRandomNumberGenerator>
3255  typename piecewise_linear_distribution<_RealType>::result_type
3256  piecewise_linear_distribution<_RealType>::
3257  operator()(_UniformRandomNumberGenerator& __urng,
3258  const param_type& __param)
3259  {
3260  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3261  __aurng(__urng);
3262 
3263  const double __p = __aurng();
3264  if (__param._M_cp.empty())
3265  return __p;
3266 
3267  auto __pos = std::lower_bound(__param._M_cp.begin(),
3268  __param._M_cp.end(), __p);
3269  const size_t __i = __pos - __param._M_cp.begin();
3270 
3271  const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3272 
3273  const double __a = 0.5 * __param._M_m[__i];
3274  const double __b = __param._M_den[__i];
3275  const double __cm = __p - __pref;
3276 
3277  _RealType __x = __param._M_int[__i];
3278  if (__a == 0)
3279  __x += __cm / __b;
3280  else
3281  {
3282  const double __d = __b * __b + 4.0 * __a * __cm;
3283  __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3284  }
3285 
3286  return __x;
3287  }
3288 
3289  template<typename _RealType>
3290  template<typename _ForwardIterator,
3291  typename _UniformRandomNumberGenerator>
3292  void
3293  piecewise_linear_distribution<_RealType>::
3294  __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3295  _UniformRandomNumberGenerator& __urng,
3296  const param_type& __param)
3297  {
3298  __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3299  // We could duplicate everything from operator()...
3300  while (__f != __t)
3301  *__f++ = this->operator()(__urng, __param);
3302  }
3303 
3304  template<typename _RealType, typename _CharT, typename _Traits>
3305  std::basic_ostream<_CharT, _Traits>&
3306  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3307  const piecewise_linear_distribution<_RealType>& __x)
3308  {
3309  using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3310 
3311  const typename __ios_base::fmtflags __flags = __os.flags();
3312  const _CharT __fill = __os.fill();
3313  const std::streamsize __precision = __os.precision();
3314  const _CharT __space = __os.widen(' ');
3315  __os.flags(__ios_base::scientific | __ios_base::left);
3316  __os.fill(__space);
3317  __os.precision(std::numeric_limits<_RealType>::max_digits10);
3318 
3319  std::vector<_RealType> __int = __x.intervals();
3320  __os << __int.size() - 1;
3321 
3322  for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3323  __os << __space << *__xit;
3324 
3325  std::vector<double> __den = __x.densities();
3326  for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3327  __os << __space << *__dit;
3328 
3329  __os.flags(__flags);
3330  __os.fill(__fill);
3331  __os.precision(__precision);
3332  return __os;
3333  }
3334 
3335  template<typename _RealType, typename _CharT, typename _Traits>
3336  std::basic_istream<_CharT, _Traits>&
3337  operator>>(std::basic_istream<_CharT, _Traits>& __is,
3338  piecewise_linear_distribution<_RealType>& __x)
3339  {
3340  using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3341 
3342  const typename __ios_base::fmtflags __flags = __is.flags();
3343  __is.flags(__ios_base::dec | __ios_base::skipws);
3344 
3345  size_t __n;
3346  if (__is >> __n)
3347  {
3348  vector<_RealType> __int_vec;
3349  if (__detail::__extract_params(__is, __int_vec, __n + 1))
3350  {
3351  vector<double> __den_vec;
3352  if (__detail::__extract_params(__is, __den_vec, __n + 1))
3353  {
3354  __x.param({ __int_vec.begin(), __int_vec.end(),
3355  __den_vec.begin() });
3356  }
3357  }
3358  }
3359  __is.flags(__flags);
3360  return __is;
3361  }
3362 
3363 
3364  template<typename _IntType, typename>
3365  seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3366  {
3367  _M_v.reserve(__il.size());
3368  for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3369  _M_v.push_back(__detail::__mod<result_type,
3370  __detail::_Shift<result_type, 32>::__value>(*__iter));
3371  }
3372 
3373  template<typename _InputIterator>
3374  seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3375  {
3376 #pragma GCC diagnostic push
3377 #pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3378  if constexpr (__is_random_access_iter<_InputIterator>::value)
3379  _M_v.reserve(std::distance(__begin, __end));
3380 #pragma GCC diagnostic pop
3381 
3382  for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3383  _M_v.push_back(__detail::__mod<result_type,
3384  __detail::_Shift<result_type, 32>::__value>(*__iter));
3385  }
3386 
3387  template<typename _RandomAccessIterator>
3388  void
3389  seed_seq::generate(_RandomAccessIterator __begin,
3390  _RandomAccessIterator __end)
3391  {
3392  typedef typename iterator_traits<_RandomAccessIterator>::value_type
3393  _Type;
3394 
3395  if (__begin == __end)
3396  return;
3397 
3398  std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3399 
3400  const size_t __n = __end - __begin;
3401  const size_t __s = _M_v.size();
3402  const size_t __t = (__n >= 623) ? 11
3403  : (__n >= 68) ? 7
3404  : (__n >= 39) ? 5
3405  : (__n >= 7) ? 3
3406  : (__n - 1) / 2;
3407  const size_t __p = (__n - __t) / 2;
3408  const size_t __q = __p + __t;
3409  const size_t __m = std::max(size_t(__s + 1), __n);
3410 
3411 #ifndef __UINT32_TYPE__
3412  struct _Up
3413  {
3414  _Up(uint_least32_t v) : _M_v(v & 0xffffffffu) { }
3415 
3416  operator uint_least32_t() const { return _M_v; }
3417 
3418  uint_least32_t _M_v;
3419  };
3420  using uint32_t = _Up;
3421 #endif
3422 
3423  // k == 0, every element in [begin,end) equals 0x8b8b8b8bu
3424  {
3425  uint32_t __r1 = 1371501266u;
3426  uint32_t __r2 = __r1 + __s;
3427  __begin[__p] += __r1;
3428  __begin[__q] = (uint32_t)__begin[__q] + __r2;
3429  __begin[0] = __r2;
3430  }
3431 
3432  for (size_t __k = 1; __k <= __s; ++__k)
3433  {
3434  const size_t __kn = __k % __n;
3435  const size_t __kpn = (__k + __p) % __n;
3436  const size_t __kqn = (__k + __q) % __n;
3437  uint32_t __arg = (__begin[__kn]
3438  ^ __begin[__kpn]
3439  ^ __begin[(__k - 1) % __n]);
3440  uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3441  uint32_t __r2 = __r1 + (uint32_t)__kn + _M_v[__k - 1];
3442  __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3443  __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3444  __begin[__kn] = __r2;
3445  }
3446 
3447  for (size_t __k = __s + 1; __k < __m; ++__k)
3448  {
3449  const size_t __kn = __k % __n;
3450  const size_t __kpn = (__k + __p) % __n;
3451  const size_t __kqn = (__k + __q) % __n;
3452  uint32_t __arg = (__begin[__kn]
3453  ^ __begin[__kpn]
3454  ^ __begin[(__k - 1) % __n]);
3455  uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3456  uint32_t __r2 = __r1 + (uint32_t)__kn;
3457  __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3458  __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3459  __begin[__kn] = __r2;
3460  }
3461 
3462  for (size_t __k = __m; __k < __m + __n; ++__k)
3463  {
3464  const size_t __kn = __k % __n;
3465  const size_t __kpn = (__k + __p) % __n;
3466  const size_t __kqn = (__k + __q) % __n;
3467  uint32_t __arg = (__begin[__kn]
3468  + __begin[__kpn]
3469  + __begin[(__k - 1) % __n]);
3470  uint32_t __r3 = 1566083941u * (__arg ^ (__arg >> 27));
3471  uint32_t __r4 = __r3 - __kn;
3472  __begin[__kpn] ^= __r3;
3473  __begin[__kqn] ^= __r4;
3474  __begin[__kn] = __r4;
3475  }
3476  }
3477 
3478 // [rand.util.canonical]
3479 // generate_canonical(RNG&)
3480 
3481 #ifndef _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3482 
3483 #pragma GCC diagnostic push
3484 #pragma GCC diagnostic ignored "-Wc++14-extensions" // for variable templates
3485 #pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
3486 
3487  // __generate_canonical_pow2 is used when Urbg::max()-Urbg::min() is
3488  // a power of two less 1. It works by calling urng() as many times as
3489  // needed to fill the target mantissa, accumulating entropy into an
3490  // integer value, converting that to the float type, and then dividing
3491  // by the range of the integer value (a constexpr power of 2,
3492  // so only adjusts the exponent) to produce a result in [0..1].
3493  // In case of an exact 1.0 result, we re-try.
3494  //
3495  // It needs to work even when the integer type used is only as big
3496  // as the float mantissa, such as uint64_t for long double. So,
3497  // commented-out assignments represent computations the Standard
3498  // prescribes but cannot be performed, or are not used. Names are
3499  // chosen to match the description in the Standard.
3500  //
3501  // When the result is close to zero, the strict Standard-prescribed
3502  // calculation may leave more low-order zeros in the mantissa than
3503  // is usually necessary. When spare entropy has been extracted, as
3504  // is usual for float and double, some or all of the spare entropy
3505  // can commonly be pulled into the result for better randomness.
3506  // Defining _GLIBCXX_GENERATE_CANONICAL_STRICT discards it instead.
3507  //
3508  // When k calls to urng() yield more bits of entropy, log2_Rk_max,
3509  // than fit into UInt, we discard some of it by overflowing, which
3510  // is OK. On converting the integer representation of the sample
3511  // to the float value, we must divide out the (possibly-truncated)
3512  // size log2_Rk.
3513  //
3514  // This implementation works with std::bfloat16, which can exactly
3515  // represent 2^32, but not with std::float16_t, limited to 2^15.
3516 
3517  template<typename _RealT, size_t __d, typename _Urbg>
3518  _RealT
3519  __generate_canonical_pow2(_Urbg& __urng)
3520  {
3521  using _UInt = typename __detail::_Select_uint_least_t<__d>::type;
3522 
3523  // Parameter __d is the actual target number of bits.
3524  // Commented-out assignments below are of values specified in
3525  // the Standard, but not used here for reasons noted.
3526  // r = 2; // Redundant, we only support radix 2.
3527  using _Rng = decltype(_Urbg::max());
3528  const _Rng __rng_range_less_1 = _Urbg::max() - _Urbg::min();
3529  // R = _UInt(__rng_range_less_1) + 1; // May wrap to 0.
3530  const auto __log2_R = __builtin_popcountg(__rng_range_less_1);
3531  const auto __log2_uint_max = sizeof(_UInt) * __CHAR_BIT__;
3532  // rd = _UInt(1) << __d; // Could overflow, UB.
3533  const unsigned __k = (__d + __log2_R - 1) / __log2_R;
3534  const unsigned __log2_Rk_max = __k * __log2_R;
3535  const unsigned __log2_Rk = // Bits of entropy actually obtained:
3536  __log2_uint_max < __log2_Rk_max ? __log2_uint_max : __log2_Rk_max;
3537  // Rk = _UInt(1) << __log2_Rk; // Likely overflows, UB.
3538  _GLIBCXX14_CONSTEXPR const _RealT __Rk
3539  = _RealT(_UInt(1) << (__log2_Rk - 1)) * _RealT(2.0);
3540 #if defined(_GLIBCXX_GENERATE_CANONICAL_STRICT)
3541  const unsigned __log2_x = __log2_Rk - __d; // # of spare entropy bits.
3542 #else
3543  const unsigned __log2_x = 0;
3544 #endif
3545  _GLIBCXX14_CONSTEXPR const _UInt __x = _UInt(1) << __log2_x;
3546  _GLIBCXX14_CONSTEXPR const _RealT __rd = __Rk / _RealT(__x);
3547  // xrd = __x << __d; // Could overflow.
3548 
3549  while (true)
3550  {
3551  _UInt __sum = _UInt(__urng() - _Urbg::min());
3552  for (unsigned __i = __k - 1, __shift = 0; __i > 0; --__i)
3553  {
3554  __shift += __log2_R;
3555  __sum |= _UInt(__urng() - _Urbg::min()) << __shift;
3556  }
3557  const _RealT __ret = _RealT(__sum >> __log2_x) / _RealT(__rd);
3558  if (__ret < _RealT(1.0))
3559  return __ret;
3560  }
3561  }
3562 
3563 
3564  template<typename _UInt>
3565  struct __gen_canon_log_res
3566  {
3567  unsigned __floor_log;
3568  _UInt __floor_pow;
3569 
3570  constexpr __gen_canon_log_res
3571  update(_UInt __base) const
3572  { return {__floor_log + 1, __floor_pow * __base}; }
3573  };
3574 
3575 
3576  template <typename _UInt1, typename _UInt2,
3577  typename _UComm = __conditional_t<(sizeof(_UInt2) > sizeof(_UInt1)),
3578  _UInt2, _UInt1>>
3579  constexpr __gen_canon_log_res<_UInt1>
3580  __gen_canon_log(_UInt1 __val, _UInt2 __base)
3581  {
3582 #if __cplusplus >= 201402L
3583  __gen_canon_log_res<_UInt1> __res{0, _UInt1(1)};
3584  if (_UComm(__base) > _UComm(__val))
3585  return __res;
3586 
3587  const _UInt1 __base1(__base);
3588  do
3589  {
3590  __val /= __base1;
3591  __res = __res.update(__base1);
3592  }
3593  while (__val >= __base1);
3594  return __res;
3595 #else
3596  return (_UComm(__val) >= _UComm(__base))
3597  ? __gen_canon_log(__val / _UInt1(__base), _UInt1(__base))
3598  .update(_UInt1(__base))
3599  : __gen_canon_log_res<_UInt1>{0, _UInt1(1)};
3600 #endif
3601  }
3602 
3603  // This version must be used when the range of possible RNG results,
3604  // Urbg::max()-Urbg::min(), is not a power of two less one. The UInt
3605  // type passed must be big enough to represent Rk, R^k, a power of R
3606  // (the range of values produced by the rng) up to twice the length
3607  // of the mantissa.
3608 
3609  template<typename _RealT, size_t __d, typename _Urbg>
3610  _RealT
3611  __generate_canonical_any(_Urbg& __urng)
3612  {
3613  // Names below are chosen to match the description in the Standard.
3614  // Parameter d is the actual target number of bits.
3615 #if (__cplusplus >= 201402L) || defined(__SIZEOF_INT128__)
3616 # define _GLIBCXX_GEN_CANON_CONST constexpr
3617 #else
3618 # define _GLIBCXX_GEN_CANON_CONST const
3619 #endif
3620 
3621  using _UIntR = typename make_unsigned<decltype(_Urbg::max())>::type;
3622  // Cannot overflow, as _Urbg::max() - _Urbg::min() is not power of
3623  // two minus one
3624  constexpr _UIntR __R = _UIntR(_Urbg::max() - _Urbg::min()) + 1;
3625  constexpr unsigned __log2R
3626  = sizeof(_UIntR) * __CHAR_BIT__ - __builtin_clzg(__R) - 1;
3627  // We overstimate number of required bits, by computing
3628  // r such that l * log2(R) >= d, so:
3629  // R^l >= (2 ^ log2(R)) ^ l == 2 ^ (log2(r) * l) >= 2^d
3630  // And then requiring l * bit_width(R) bits.
3631  constexpr unsigned __l = (__d + __log2R - 1) / __log2R;
3632  constexpr unsigned __bits = (__log2R + 1) * __l;
3633  using _UInt = typename __detail::_Select_uint_least_t<__bits>::type;
3634 
3635  _GLIBCXX_GEN_CANON_CONST _UInt __rd = _UInt(1) << __d;
3636  _GLIBCXX_GEN_CANON_CONST auto __logRrd = __gen_canon_log(__rd, __R);
3637  _GLIBCXX_GEN_CANON_CONST unsigned __k
3638  = __logRrd.__floor_log + (__rd > __logRrd.__floor_pow);
3639 
3640  _GLIBCXX_GEN_CANON_CONST _UInt __Rk
3641  = (__k > __logRrd.__floor_log)
3642  ? _UInt(__logRrd.__floor_pow) * _UInt(__R)
3643  : _UInt(__logRrd.__floor_pow);
3644  _GLIBCXX_GEN_CANON_CONST _UInt __x = __Rk / __rd;
3645 
3646  while (true)
3647  {
3648  _UInt __Ri{1};
3649  _UInt __sum(__urng() - _Urbg::min());
3650  for (int __i = __k - 1; __i > 0; --__i)
3651  {
3652  __Ri *= _UInt(__R);
3653  __sum += _UInt(__urng() - _Urbg::min()) * __Ri;
3654  }
3655  const _RealT __ret = _RealT(__sum / __x) / _RealT(__rd);
3656  if (__ret < _RealT(1.0))
3657  return __ret;
3658  }
3659 #undef _GLIBCXX_GEN_CANON_CONST
3660  }
3661 
3662 #if !defined(_GLIBCXX_GENERATE_CANONICAL_STRICT)
3663  template <typename _Tp>
3664  const bool __is_rand_dist_float_v = is_floating_point<_Tp>::value;
3665 #else
3666  template <typename _Tp> const bool __is_rand_dist_float_v = false;
3667  template <> const bool __is_rand_dist_float_v<float> = true;
3668  template <> const bool __is_rand_dist_float_v<double> = true;
3669  template <> const bool __is_rand_dist_float_v<long double> = true;
3670 #endif
3671 
3672  // Note, this works even when (__range + 1) overflows:
3673  template <typename _Rng>
3674  constexpr bool __is_power_of_2_less_1(_Rng __range)
3675  { return ((__range + 1) & __range) == 0; };
3676 
3677 _GLIBCXX_BEGIN_INLINE_ABI_NAMESPACE(_V2)
3678  /** Produce a random floating-point value in the range [0..1)
3679  *
3680  * The result of `std::generate_canonical<RealT,digits>(urng)` is a
3681  * random floating-point value of type `RealT` in the range [0..1),
3682  * using entropy provided by the uniform random bit generator `urng`.
3683  * A value for `digits` may be passed to limit the precision of the
3684  * result to so many bits, but normally `-1u` is passed to get the
3685  * native precision of `RealT`. As many `urng()` calls are made as
3686  * needed to obtain the required entropy. On rare occasions, more
3687  * `urng()` calls are used. It is fastest when the value of
3688  * `Urbg::max()` is a power of two less one, such as from a
3689  * `std::philox4x32` (for `float`) or `philox4x64` (for `double`).
3690  *
3691  * @since C++11
3692  */
3693  template<typename _RealT, size_t __digits,
3694  typename _Urbg>
3695  _RealT
3696  generate_canonical(_Urbg& __urng)
3697  {
3698 #ifdef __glibcxx_concepts
3699  static_assert(uniform_random_bit_generator<_Urbg>);
3700 #endif
3701  static_assert(__is_rand_dist_float_v<_RealT>,
3702  "template argument must be a floating point type");
3703  static_assert(__digits != 0 && _Urbg::max() > _Urbg::min(),
3704  "random samples with 0 bits are not meaningful");
3705  static_assert(std::numeric_limits<_RealT>::radix == 2,
3706  "only base-2 float types are supported");
3707 #if defined(__STDCPP_FLOAT16_T__)
3708  static_assert(! is_same_v<_RealT, _Float16>,
3709  "float16_t type is not supported, consider using bfloat16_t");
3710 #endif
3711 
3712  const unsigned __d_max = std::numeric_limits<_RealT>::digits;
3713  const unsigned __d = __digits > __d_max ? __d_max : __digits;
3714 
3715  // If the RNG range is a power of 2 less 1, the float type mantissa
3716  // is enough bits. If not, we need more.
3717  if constexpr (__is_power_of_2_less_1(_Urbg::max() - _Urbg::min()))
3718  return __generate_canonical_pow2<_RealT, __d>(__urng);
3719  else // Need up to 2x bits.
3720  return __generate_canonical_any<_RealT, __d>(__urng);
3721  }
3722 _GLIBCXX_END_INLINE_ABI_NAMESPACE(_V2)
3723 
3724 #pragma GCC diagnostic pop
3725 
3726 #else // _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3727 
3728  // This is the pre-P0952 definition, to reproduce old results.
3729 
3730  template<typename _RealType, size_t __bits,
3731  typename _UniformRandomNumberGenerator>
3732  _RealType
3733  generate_canonical(_UniformRandomNumberGenerator& __urng)
3734  {
3735  static_assert(std::is_floating_point<_RealType>::value,
3736  "template argument must be a floating point type");
3737 
3738  const size_t __b
3739  = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3740  __bits);
3741  const long double __r = static_cast<long double>(__urng.max())
3742  - static_cast<long double>(__urng.min()) + 1.0L;
3743  const size_t __log2r = std::log(__r) / std::log(2.0L);
3744  const size_t __m = std::max<size_t>(1UL,
3745  (__b + __log2r - 1UL) / __log2r);
3746  _RealType __ret;
3747  _RealType __sum = _RealType(0);
3748  _RealType __tmp = _RealType(1);
3749  for (size_t __k = __m; __k != 0; --__k)
3750  {
3751  __sum += _RealType(__urng() - __urng.min()) * __tmp;
3752  __tmp *= __r;
3753  }
3754  __ret = __sum / __tmp;
3755  if (__builtin_expect(__ret >= _RealType(1), 0))
3756  {
3757 # if _GLIBCXX_USE_C99_MATH_FUNCS
3758  __ret = std::nextafter(_RealType(1), _RealType(0));
3759 # else
3760  __ret = _RealType(1)
3761  - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3762 # endif
3763  }
3764  return __ret;
3765  }
3766 
3767 #endif // _GLIBCXX_USE_OLD_GENERATE_CANONICAL
3768 
3769 _GLIBCXX_END_NAMESPACE_VERSION
3770 } // namespace
3771 
3772 #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::requires
requires requires
Definition: complex:1948
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:256
std::min
constexpr const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:232
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:133
std::partial_sum
constexpr _OutputIterator partial_sum(_InputIterator __first, _InputIterator __last, _OutputIterator __result)
Return list of partial sums.
Definition: stl_numeric.h:255
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:1558
std::generate_canonical
_RealT generate_canonical(_Urbg &__urng)
Definition: bits/random.tcc:3696
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:2579
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::philox_engine
A discrete pseudorandom number generator with weak cryptographic properties.
Definition: random.h:2060
std::uniform_real_distribution
Uniform continuous distribution for random numbers.
Definition: random.h:2493
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::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::negative_binomial_distribution::operator()
result_type operator()(_UniformRandomNumberGenerator &__urng)
Generating functions.
Definition: bits/random.tcc:1284
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::weibull_distribution
A weibull_distribution random number distribution.
Definition: random.h:5719
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::piecewise_linear_distribution
A piecewise_linear_distribution random number distribution.
Definition: random.h:6686
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
Generated by doxygen 1.9.1