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