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