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