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