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