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