]>
Commit | Line | Data |
---|---|---|
8e79468d BK |
1 | // random number generation (out of line) -*- C++ -*- |
2 | ||
99dee823 | 3 | // Copyright (C) 2009-2021 Free Software Foundation, Inc. |
8e79468d BK |
4 | // |
5 | // This file is part of the GNU ISO C++ Library. This library is free | |
6 | // software; you can redistribute it and/or modify it under the | |
7 | // terms of the GNU General Public License as published by the | |
748086b7 | 8 | // Free Software Foundation; either version 3, or (at your option) |
8e79468d BK |
9 | // any later version. |
10 | ||
11 | // This library is distributed in the hope that it will be useful, | |
12 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
14 | // GNU General Public License for more details. | |
15 | ||
748086b7 JJ |
16 | // Under Section 7 of GPL version 3, you are granted additional |
17 | // permissions described in the GCC Runtime Library Exception, version | |
18 | // 3.1, as published by the Free Software Foundation. | |
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/>. | |
8e79468d | 24 | |
8e79468d BK |
25 | |
26 | /** @file tr1/random.tcc | |
27 | * This is an internal header file, included by other library headers. | |
f910786b | 28 | * Do not attempt to use it directly. @headername{tr1/random} |
8e79468d BK |
29 | */ |
30 | ||
06f29237 PC |
31 | #ifndef _GLIBCXX_TR1_RANDOM_TCC |
32 | #define _GLIBCXX_TR1_RANDOM_TCC 1 | |
33 | ||
12ffa228 | 34 | namespace std _GLIBCXX_VISIBILITY(default) |
8e79468d | 35 | { |
4a15d842 FD |
36 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
37 | ||
8e79468d BK |
38 | namespace tr1 |
39 | { | |
8e79468d BK |
40 | /* |
41 | * (Further) implementation-space details. | |
42 | */ | |
43 | namespace __detail | |
44 | { | |
45 | // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid | |
46 | // integer overflow. | |
47 | // | |
48 | // Because a and c are compile-time integral constants the compiler kindly | |
49 | // elides any unreachable paths. | |
50 | // | |
51 | // Preconditions: a > 0, m > 0. | |
52 | // | |
53 | template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool> | |
54 | struct _Mod | |
55 | { | |
56 | static _Tp | |
57 | __calc(_Tp __x) | |
58 | { | |
59 | if (__a == 1) | |
60 | __x %= __m; | |
61 | else | |
62 | { | |
63 | static const _Tp __q = __m / __a; | |
64 | static const _Tp __r = __m % __a; | |
65 | ||
66 | _Tp __t1 = __a * (__x % __q); | |
67 | _Tp __t2 = __r * (__x / __q); | |
68 | if (__t1 >= __t2) | |
69 | __x = __t1 - __t2; | |
70 | else | |
71 | __x = __m - __t2 + __t1; | |
72 | } | |
73 | ||
74 | if (__c != 0) | |
75 | { | |
76 | const _Tp __d = __m - __x; | |
77 | if (__d > __c) | |
78 | __x += __c; | |
79 | else | |
80 | __x = __c - __d; | |
81 | } | |
82 | return __x; | |
83 | } | |
84 | }; | |
85 | ||
86 | // Special case for m == 0 -- use unsigned integer overflow as modulo | |
87 | // operator. | |
88 | template<typename _Tp, _Tp __a, _Tp __c, _Tp __m> | |
89 | struct _Mod<_Tp, __a, __c, __m, true> | |
90 | { | |
91 | static _Tp | |
92 | __calc(_Tp __x) | |
93 | { return __a * __x + __c; } | |
94 | }; | |
95 | } // namespace __detail | |
96 | ||
300ea283 PC |
97 | template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> |
98 | const _UIntType | |
99 | linear_congruential<_UIntType, __a, __c, __m>::multiplier; | |
100 | ||
101 | template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> | |
102 | const _UIntType | |
103 | linear_congruential<_UIntType, __a, __c, __m>::increment; | |
104 | ||
105 | template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> | |
106 | const _UIntType | |
107 | linear_congruential<_UIntType, __a, __c, __m>::modulus; | |
108 | ||
8e79468d BK |
109 | /** |
110 | * Seeds the LCR with integral value @p __x0, adjusted so that the | |
111 | * ring identity is never a member of the convergence set. | |
112 | */ | |
113 | template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> | |
114 | void | |
115 | linear_congruential<_UIntType, __a, __c, __m>:: | |
116 | seed(unsigned long __x0) | |
117 | { | |
118 | if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) | |
119 | && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) | |
120 | _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); | |
121 | else | |
122 | _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); | |
123 | } | |
124 | ||
125 | /** | |
126 | * Seeds the LCR engine with a value generated by @p __g. | |
127 | */ | |
128 | template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> | |
129 | template<class _Gen> | |
130 | void | |
131 | linear_congruential<_UIntType, __a, __c, __m>:: | |
132 | seed(_Gen& __g, false_type) | |
133 | { | |
134 | _UIntType __x0 = __g(); | |
135 | if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) | |
136 | && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) | |
137 | _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); | |
138 | else | |
139 | _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); | |
140 | } | |
141 | ||
142 | /** | |
143 | * Gets the next generated value in sequence. | |
144 | */ | |
145 | template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> | |
146 | typename linear_congruential<_UIntType, __a, __c, __m>::result_type | |
147 | linear_congruential<_UIntType, __a, __c, __m>:: | |
148 | operator()() | |
149 | { | |
150 | _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x); | |
151 | return _M_x; | |
152 | } | |
153 | ||
154 | template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, | |
155 | typename _CharT, typename _Traits> | |
156 | std::basic_ostream<_CharT, _Traits>& | |
157 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
158 | const linear_congruential<_UIntType, __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<class _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<_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 | ||
300ea283 PC |
194 | template<class _UIntType, int __w, int __n, int __m, int __r, |
195 | _UIntType __a, int __u, int __s, | |
196 | _UIntType __b, int __t, _UIntType __c, int __l> | |
197 | const int | |
198 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
199 | __b, __t, __c, __l>::word_size; | |
200 | ||
201 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
202 | _UIntType __a, int __u, int __s, | |
203 | _UIntType __b, int __t, _UIntType __c, int __l> | |
204 | const int | |
205 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
206 | __b, __t, __c, __l>::state_size; | |
207 | ||
208 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
209 | _UIntType __a, int __u, int __s, | |
210 | _UIntType __b, int __t, _UIntType __c, int __l> | |
211 | const int | |
212 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
213 | __b, __t, __c, __l>::shift_size; | |
214 | ||
215 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
216 | _UIntType __a, int __u, int __s, | |
217 | _UIntType __b, int __t, _UIntType __c, int __l> | |
218 | const int | |
219 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
220 | __b, __t, __c, __l>::mask_bits; | |
221 | ||
222 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
223 | _UIntType __a, int __u, int __s, | |
224 | _UIntType __b, int __t, _UIntType __c, int __l> | |
225 | const _UIntType | |
226 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
227 | __b, __t, __c, __l>::parameter_a; | |
228 | ||
229 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
230 | _UIntType __a, int __u, int __s, | |
231 | _UIntType __b, int __t, _UIntType __c, int __l> | |
232 | const int | |
233 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
234 | __b, __t, __c, __l>::output_u; | |
235 | ||
236 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
237 | _UIntType __a, int __u, int __s, | |
238 | _UIntType __b, int __t, _UIntType __c, int __l> | |
239 | const int | |
240 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
241 | __b, __t, __c, __l>::output_s; | |
242 | ||
243 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
244 | _UIntType __a, int __u, int __s, | |
245 | _UIntType __b, int __t, _UIntType __c, int __l> | |
246 | const _UIntType | |
247 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
248 | __b, __t, __c, __l>::output_b; | |
249 | ||
250 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
251 | _UIntType __a, int __u, int __s, | |
252 | _UIntType __b, int __t, _UIntType __c, int __l> | |
253 | const int | |
254 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
255 | __b, __t, __c, __l>::output_t; | |
256 | ||
257 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
258 | _UIntType __a, int __u, int __s, | |
259 | _UIntType __b, int __t, _UIntType __c, int __l> | |
260 | const _UIntType | |
261 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
262 | __b, __t, __c, __l>::output_c; | |
263 | ||
264 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
265 | _UIntType __a, int __u, int __s, | |
266 | _UIntType __b, int __t, _UIntType __c, int __l> | |
267 | const int | |
268 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
269 | __b, __t, __c, __l>::output_l; | |
270 | ||
8e79468d BK |
271 | template<class _UIntType, int __w, int __n, int __m, int __r, |
272 | _UIntType __a, int __u, int __s, | |
273 | _UIntType __b, int __t, _UIntType __c, int __l> | |
274 | void | |
275 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
276 | __b, __t, __c, __l>:: | |
277 | seed(unsigned long __value) | |
278 | { | |
279 | _M_x[0] = __detail::__mod<_UIntType, 1, 0, | |
280 | __detail::_Shift<_UIntType, __w>::__value>(__value); | |
281 | ||
282 | for (int __i = 1; __i < state_size; ++__i) | |
283 | { | |
284 | _UIntType __x = _M_x[__i - 1]; | |
285 | __x ^= __x >> (__w - 2); | |
286 | __x *= 1812433253ul; | |
287 | __x += __i; | |
288 | _M_x[__i] = __detail::__mod<_UIntType, 1, 0, | |
289 | __detail::_Shift<_UIntType, __w>::__value>(__x); | |
290 | } | |
291 | _M_p = state_size; | |
292 | } | |
293 | ||
294 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
295 | _UIntType __a, int __u, int __s, | |
296 | _UIntType __b, int __t, _UIntType __c, int __l> | |
297 | template<class _Gen> | |
298 | void | |
299 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
300 | __b, __t, __c, __l>:: | |
301 | seed(_Gen& __gen, false_type) | |
302 | { | |
303 | for (int __i = 0; __i < state_size; ++__i) | |
304 | _M_x[__i] = __detail::__mod<_UIntType, 1, 0, | |
305 | __detail::_Shift<_UIntType, __w>::__value>(__gen()); | |
306 | _M_p = state_size; | |
307 | } | |
308 | ||
309 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
310 | _UIntType __a, int __u, int __s, | |
311 | _UIntType __b, int __t, _UIntType __c, int __l> | |
312 | typename | |
313 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
314 | __b, __t, __c, __l>::result_type | |
315 | mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, | |
316 | __b, __t, __c, __l>:: | |
317 | operator()() | |
318 | { | |
319 | // Reload the vector - cost is O(n) amortized over n calls. | |
320 | if (_M_p >= state_size) | |
321 | { | |
322 | const _UIntType __upper_mask = (~_UIntType()) << __r; | |
323 | const _UIntType __lower_mask = ~__upper_mask; | |
324 | ||
325 | for (int __k = 0; __k < (__n - __m); ++__k) | |
326 | { | |
327 | _UIntType __y = ((_M_x[__k] & __upper_mask) | |
328 | | (_M_x[__k + 1] & __lower_mask)); | |
329 | _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) | |
330 | ^ ((__y & 0x01) ? __a : 0)); | |
331 | } | |
332 | ||
333 | for (int __k = (__n - __m); __k < (__n - 1); ++__k) | |
334 | { | |
335 | _UIntType __y = ((_M_x[__k] & __upper_mask) | |
336 | | (_M_x[__k + 1] & __lower_mask)); | |
337 | _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) | |
338 | ^ ((__y & 0x01) ? __a : 0)); | |
339 | } | |
340 | ||
341 | _UIntType __y = ((_M_x[__n - 1] & __upper_mask) | |
342 | | (_M_x[0] & __lower_mask)); | |
343 | _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) | |
344 | ^ ((__y & 0x01) ? __a : 0)); | |
345 | _M_p = 0; | |
346 | } | |
347 | ||
348 | // Calculate o(x(i)). | |
349 | result_type __z = _M_x[_M_p++]; | |
350 | __z ^= (__z >> __u); | |
351 | __z ^= (__z << __s) & __b; | |
352 | __z ^= (__z << __t) & __c; | |
353 | __z ^= (__z >> __l); | |
354 | ||
355 | return __z; | |
356 | } | |
357 | ||
358 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
359 | _UIntType __a, int __u, int __s, _UIntType __b, int __t, | |
360 | _UIntType __c, int __l, | |
361 | typename _CharT, typename _Traits> | |
362 | std::basic_ostream<_CharT, _Traits>& | |
363 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
364 | const mersenne_twister<_UIntType, __w, __n, __m, | |
365 | __r, __a, __u, __s, __b, __t, __c, __l>& __x) | |
366 | { | |
367 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
368 | typedef typename __ostream_type::ios_base __ios_base; | |
369 | ||
370 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
371 | const _CharT __fill = __os.fill(); | |
372 | const _CharT __space = __os.widen(' '); | |
373 | __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); | |
374 | __os.fill(__space); | |
375 | ||
376 | for (int __i = 0; __i < __n - 1; ++__i) | |
377 | __os << __x._M_x[__i] << __space; | |
378 | __os << __x._M_x[__n - 1]; | |
379 | ||
380 | __os.flags(__flags); | |
381 | __os.fill(__fill); | |
382 | return __os; | |
383 | } | |
384 | ||
385 | template<class _UIntType, int __w, int __n, int __m, int __r, | |
386 | _UIntType __a, int __u, int __s, _UIntType __b, int __t, | |
387 | _UIntType __c, int __l, | |
388 | typename _CharT, typename _Traits> | |
389 | std::basic_istream<_CharT, _Traits>& | |
390 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
391 | mersenne_twister<_UIntType, __w, __n, __m, | |
392 | __r, __a, __u, __s, __b, __t, __c, __l>& __x) | |
393 | { | |
394 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
395 | typedef typename __istream_type::ios_base __ios_base; | |
396 | ||
397 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
398 | __is.flags(__ios_base::dec | __ios_base::skipws); | |
399 | ||
400 | for (int __i = 0; __i < __n; ++__i) | |
401 | __is >> __x._M_x[__i]; | |
402 | ||
403 | __is.flags(__flags); | |
404 | return __is; | |
405 | } | |
406 | ||
407 | ||
300ea283 PC |
408 | template<typename _IntType, _IntType __m, int __s, int __r> |
409 | const _IntType | |
410 | subtract_with_carry<_IntType, __m, __s, __r>::modulus; | |
411 | ||
412 | template<typename _IntType, _IntType __m, int __s, int __r> | |
413 | const int | |
414 | subtract_with_carry<_IntType, __m, __s, __r>::long_lag; | |
415 | ||
416 | template<typename _IntType, _IntType __m, int __s, int __r> | |
417 | const int | |
418 | subtract_with_carry<_IntType, __m, __s, __r>::short_lag; | |
419 | ||
8e79468d BK |
420 | template<typename _IntType, _IntType __m, int __s, int __r> |
421 | void | |
422 | subtract_with_carry<_IntType, __m, __s, __r>:: | |
423 | seed(unsigned long __value) | |
424 | { | |
425 | if (__value == 0) | |
426 | __value = 19780503; | |
427 | ||
95fe602e | 428 | std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563> |
8e79468d BK |
429 | __lcg(__value); |
430 | ||
431 | for (int __i = 0; __i < long_lag; ++__i) | |
432 | _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg()); | |
433 | ||
434 | _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; | |
435 | _M_p = 0; | |
436 | } | |
437 | ||
438 | template<typename _IntType, _IntType __m, int __s, int __r> | |
439 | template<class _Gen> | |
440 | void | |
441 | subtract_with_carry<_IntType, __m, __s, __r>:: | |
442 | seed(_Gen& __gen, false_type) | |
443 | { | |
444 | const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32; | |
445 | ||
446 | for (int __i = 0; __i < long_lag; ++__i) | |
447 | { | |
448 | _UIntType __tmp = 0; | |
449 | _UIntType __factor = 1; | |
450 | for (int __j = 0; __j < __n; ++__j) | |
451 | { | |
452 | __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0> | |
453 | (__gen()) * __factor; | |
454 | __factor *= __detail::_Shift<_UIntType, 32>::__value; | |
455 | } | |
456 | _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp); | |
457 | } | |
458 | _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; | |
459 | _M_p = 0; | |
460 | } | |
461 | ||
462 | template<typename _IntType, _IntType __m, int __s, int __r> | |
463 | typename subtract_with_carry<_IntType, __m, __s, __r>::result_type | |
464 | subtract_with_carry<_IntType, __m, __s, __r>:: | |
465 | operator()() | |
466 | { | |
467 | // Derive short lag index from current index. | |
468 | int __ps = _M_p - short_lag; | |
469 | if (__ps < 0) | |
470 | __ps += long_lag; | |
471 | ||
472 | // Calculate new x(i) without overflow or division. | |
473 | // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry | |
474 | // cannot overflow. | |
475 | _UIntType __xi; | |
476 | if (_M_x[__ps] >= _M_x[_M_p] + _M_carry) | |
477 | { | |
478 | __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry; | |
479 | _M_carry = 0; | |
480 | } | |
481 | else | |
482 | { | |
483 | __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps]; | |
484 | _M_carry = 1; | |
485 | } | |
486 | _M_x[_M_p] = __xi; | |
487 | ||
488 | // Adjust current index to loop around in ring buffer. | |
489 | if (++_M_p >= long_lag) | |
490 | _M_p = 0; | |
491 | ||
492 | return __xi; | |
493 | } | |
494 | ||
495 | template<typename _IntType, _IntType __m, int __s, int __r, | |
496 | typename _CharT, typename _Traits> | |
497 | std::basic_ostream<_CharT, _Traits>& | |
498 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
499 | const subtract_with_carry<_IntType, __m, __s, __r>& __x) | |
500 | { | |
501 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
502 | typedef typename __ostream_type::ios_base __ios_base; | |
503 | ||
504 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
505 | const _CharT __fill = __os.fill(); | |
506 | const _CharT __space = __os.widen(' '); | |
507 | __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); | |
508 | __os.fill(__space); | |
509 | ||
510 | for (int __i = 0; __i < __r; ++__i) | |
511 | __os << __x._M_x[__i] << __space; | |
512 | __os << __x._M_carry; | |
513 | ||
514 | __os.flags(__flags); | |
515 | __os.fill(__fill); | |
516 | return __os; | |
517 | } | |
518 | ||
519 | template<typename _IntType, _IntType __m, int __s, int __r, | |
520 | typename _CharT, typename _Traits> | |
521 | std::basic_istream<_CharT, _Traits>& | |
522 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
523 | subtract_with_carry<_IntType, __m, __s, __r>& __x) | |
524 | { | |
525 | typedef std::basic_ostream<_CharT, _Traits> __istream_type; | |
526 | typedef typename __istream_type::ios_base __ios_base; | |
527 | ||
528 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
529 | __is.flags(__ios_base::dec | __ios_base::skipws); | |
530 | ||
531 | for (int __i = 0; __i < __r; ++__i) | |
532 | __is >> __x._M_x[__i]; | |
533 | __is >> __x._M_carry; | |
534 | ||
535 | __is.flags(__flags); | |
536 | return __is; | |
537 | } | |
538 | ||
539 | ||
300ea283 PC |
540 | template<typename _RealType, int __w, int __s, int __r> |
541 | const int | |
542 | subtract_with_carry_01<_RealType, __w, __s, __r>::word_size; | |
543 | ||
544 | template<typename _RealType, int __w, int __s, int __r> | |
545 | const int | |
546 | subtract_with_carry_01<_RealType, __w, __s, __r>::long_lag; | |
547 | ||
548 | template<typename _RealType, int __w, int __s, int __r> | |
549 | const int | |
550 | subtract_with_carry_01<_RealType, __w, __s, __r>::short_lag; | |
551 | ||
8e79468d BK |
552 | template<typename _RealType, int __w, int __s, int __r> |
553 | void | |
554 | subtract_with_carry_01<_RealType, __w, __s, __r>:: | |
555 | _M_initialize_npows() | |
556 | { | |
557 | for (int __j = 0; __j < __n; ++__j) | |
558 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
95fe602e | 559 | _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32); |
8e79468d BK |
560 | #else |
561 | _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32); | |
562 | #endif | |
563 | } | |
564 | ||
565 | template<typename _RealType, int __w, int __s, int __r> | |
566 | void | |
567 | subtract_with_carry_01<_RealType, __w, __s, __r>:: | |
568 | seed(unsigned long __value) | |
569 | { | |
570 | if (__value == 0) | |
571 | __value = 19780503; | |
572 | ||
573 | // _GLIBCXX_RESOLVE_LIB_DEFECTS | |
574 | // 512. Seeding subtract_with_carry_01 from a single unsigned long. | |
95fe602e | 575 | std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563> |
8e79468d BK |
576 | __lcg(__value); |
577 | ||
578 | this->seed(__lcg); | |
579 | } | |
580 | ||
581 | template<typename _RealType, int __w, int __s, int __r> | |
582 | template<class _Gen> | |
583 | void | |
584 | subtract_with_carry_01<_RealType, __w, __s, __r>:: | |
585 | seed(_Gen& __gen, false_type) | |
586 | { | |
587 | for (int __i = 0; __i < long_lag; ++__i) | |
588 | { | |
589 | for (int __j = 0; __j < __n - 1; ++__j) | |
590 | _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen()); | |
591 | _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, | |
592 | __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen()); | |
593 | } | |
594 | ||
595 | _M_carry = 1; | |
596 | for (int __j = 0; __j < __n; ++__j) | |
597 | if (_M_x[long_lag - 1][__j] != 0) | |
598 | { | |
599 | _M_carry = 0; | |
600 | break; | |
601 | } | |
602 | ||
603 | _M_p = 0; | |
604 | } | |
605 | ||
606 | template<typename _RealType, int __w, int __s, int __r> | |
607 | typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type | |
608 | subtract_with_carry_01<_RealType, __w, __s, __r>:: | |
609 | operator()() | |
610 | { | |
611 | // Derive short lag index from current index. | |
612 | int __ps = _M_p - short_lag; | |
613 | if (__ps < 0) | |
614 | __ps += long_lag; | |
615 | ||
616 | _UInt32Type __new_carry; | |
617 | for (int __j = 0; __j < __n - 1; ++__j) | |
618 | { | |
619 | if (_M_x[__ps][__j] > _M_x[_M_p][__j] | |
620 | || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0)) | |
621 | __new_carry = 0; | |
622 | else | |
623 | __new_carry = 1; | |
624 | ||
625 | _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry; | |
626 | _M_carry = __new_carry; | |
627 | } | |
628 | ||
629 | if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1] | |
630 | || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0)) | |
631 | __new_carry = 0; | |
632 | else | |
633 | __new_carry = 1; | |
634 | ||
635 | _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, | |
636 | __detail::_Shift<_UInt32Type, __w % 32>::__value> | |
637 | (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry); | |
638 | _M_carry = __new_carry; | |
639 | ||
640 | result_type __ret = 0.0; | |
641 | for (int __j = 0; __j < __n; ++__j) | |
642 | __ret += _M_x[_M_p][__j] * _M_npows[__j]; | |
643 | ||
644 | // Adjust current index to loop around in ring buffer. | |
645 | if (++_M_p >= long_lag) | |
646 | _M_p = 0; | |
647 | ||
648 | return __ret; | |
649 | } | |
650 | ||
651 | template<typename _RealType, int __w, int __s, int __r, | |
652 | typename _CharT, typename _Traits> | |
653 | std::basic_ostream<_CharT, _Traits>& | |
654 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
655 | const subtract_with_carry_01<_RealType, __w, __s, __r>& __x) | |
656 | { | |
657 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
658 | typedef typename __ostream_type::ios_base __ios_base; | |
659 | ||
660 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
661 | const _CharT __fill = __os.fill(); | |
662 | const _CharT __space = __os.widen(' '); | |
663 | __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); | |
664 | __os.fill(__space); | |
665 | ||
666 | for (int __i = 0; __i < __r; ++__i) | |
667 | for (int __j = 0; __j < __x.__n; ++__j) | |
668 | __os << __x._M_x[__i][__j] << __space; | |
669 | __os << __x._M_carry; | |
670 | ||
671 | __os.flags(__flags); | |
672 | __os.fill(__fill); | |
673 | return __os; | |
674 | } | |
675 | ||
676 | template<typename _RealType, int __w, int __s, int __r, | |
677 | typename _CharT, typename _Traits> | |
678 | std::basic_istream<_CharT, _Traits>& | |
679 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
680 | subtract_with_carry_01<_RealType, __w, __s, __r>& __x) | |
681 | { | |
682 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
683 | typedef typename __istream_type::ios_base __ios_base; | |
684 | ||
685 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
686 | __is.flags(__ios_base::dec | __ios_base::skipws); | |
687 | ||
688 | for (int __i = 0; __i < __r; ++__i) | |
689 | for (int __j = 0; __j < __x.__n; ++__j) | |
690 | __is >> __x._M_x[__i][__j]; | |
691 | __is >> __x._M_carry; | |
692 | ||
693 | __is.flags(__flags); | |
694 | return __is; | |
695 | } | |
696 | ||
300ea283 PC |
697 | template<class _UniformRandomNumberGenerator, int __p, int __r> |
698 | const int | |
699 | discard_block<_UniformRandomNumberGenerator, __p, __r>::block_size; | |
700 | ||
701 | template<class _UniformRandomNumberGenerator, int __p, int __r> | |
702 | const int | |
703 | discard_block<_UniformRandomNumberGenerator, __p, __r>::used_block; | |
8e79468d BK |
704 | |
705 | template<class _UniformRandomNumberGenerator, int __p, int __r> | |
706 | typename discard_block<_UniformRandomNumberGenerator, | |
707 | __p, __r>::result_type | |
708 | discard_block<_UniformRandomNumberGenerator, __p, __r>:: | |
709 | operator()() | |
710 | { | |
711 | if (_M_n >= used_block) | |
712 | { | |
713 | while (_M_n < block_size) | |
714 | { | |
715 | _M_b(); | |
716 | ++_M_n; | |
717 | } | |
718 | _M_n = 0; | |
719 | } | |
720 | ++_M_n; | |
721 | return _M_b(); | |
722 | } | |
723 | ||
724 | template<class _UniformRandomNumberGenerator, int __p, int __r, | |
725 | typename _CharT, typename _Traits> | |
726 | std::basic_ostream<_CharT, _Traits>& | |
727 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
728 | const discard_block<_UniformRandomNumberGenerator, | |
729 | __p, __r>& __x) | |
730 | { | |
731 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
732 | typedef typename __ostream_type::ios_base __ios_base; | |
733 | ||
734 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
735 | const _CharT __fill = __os.fill(); | |
736 | const _CharT __space = __os.widen(' '); | |
737 | __os.flags(__ios_base::dec | __ios_base::fixed | |
738 | | __ios_base::left); | |
739 | __os.fill(__space); | |
740 | ||
741 | __os << __x._M_b << __space << __x._M_n; | |
742 | ||
743 | __os.flags(__flags); | |
744 | __os.fill(__fill); | |
745 | return __os; | |
746 | } | |
747 | ||
748 | template<class _UniformRandomNumberGenerator, int __p, int __r, | |
749 | typename _CharT, typename _Traits> | |
750 | std::basic_istream<_CharT, _Traits>& | |
751 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
752 | discard_block<_UniformRandomNumberGenerator, __p, __r>& __x) | |
753 | { | |
754 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
755 | typedef typename __istream_type::ios_base __ios_base; | |
756 | ||
757 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
758 | __is.flags(__ios_base::dec | __ios_base::skipws); | |
759 | ||
760 | __is >> __x._M_b >> __x._M_n; | |
761 | ||
762 | __is.flags(__flags); | |
763 | return __is; | |
764 | } | |
765 | ||
766 | ||
300ea283 PC |
767 | template<class _UniformRandomNumberGenerator1, int __s1, |
768 | class _UniformRandomNumberGenerator2, int __s2> | |
769 | const int | |
770 | xor_combine<_UniformRandomNumberGenerator1, __s1, | |
771 | _UniformRandomNumberGenerator2, __s2>::shift1; | |
772 | ||
773 | template<class _UniformRandomNumberGenerator1, int __s1, | |
774 | class _UniformRandomNumberGenerator2, int __s2> | |
775 | const int | |
776 | xor_combine<_UniformRandomNumberGenerator1, __s1, | |
777 | _UniformRandomNumberGenerator2, __s2>::shift2; | |
778 | ||
8e79468d BK |
779 | template<class _UniformRandomNumberGenerator1, int __s1, |
780 | class _UniformRandomNumberGenerator2, int __s2> | |
781 | void | |
782 | xor_combine<_UniformRandomNumberGenerator1, __s1, | |
783 | _UniformRandomNumberGenerator2, __s2>:: | |
784 | _M_initialize_max() | |
785 | { | |
786 | const int __w = std::numeric_limits<result_type>::digits; | |
787 | ||
788 | const result_type __m1 = | |
789 | std::min(result_type(_M_b1.max() - _M_b1.min()), | |
790 | __detail::_Shift<result_type, __w - __s1>::__value - 1); | |
791 | ||
792 | const result_type __m2 = | |
793 | std::min(result_type(_M_b2.max() - _M_b2.min()), | |
794 | __detail::_Shift<result_type, __w - __s2>::__value - 1); | |
795 | ||
796 | // NB: In TR1 s1 is not required to be >= s2. | |
797 | if (__s1 < __s2) | |
798 | _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1; | |
799 | else | |
800 | _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2; | |
801 | } | |
802 | ||
803 | template<class _UniformRandomNumberGenerator1, int __s1, | |
804 | class _UniformRandomNumberGenerator2, int __s2> | |
805 | typename xor_combine<_UniformRandomNumberGenerator1, __s1, | |
806 | _UniformRandomNumberGenerator2, __s2>::result_type | |
807 | xor_combine<_UniformRandomNumberGenerator1, __s1, | |
808 | _UniformRandomNumberGenerator2, __s2>:: | |
809 | _M_initialize_max_aux(result_type __a, result_type __b, int __d) | |
810 | { | |
811 | const result_type __two2d = result_type(1) << __d; | |
812 | const result_type __c = __a * __two2d; | |
813 | ||
814 | if (__a == 0 || __b < __two2d) | |
815 | return __c + __b; | |
816 | ||
817 | const result_type __t = std::max(__c, __b); | |
818 | const result_type __u = std::min(__c, __b); | |
819 | ||
820 | result_type __ub = __u; | |
821 | result_type __p; | |
822 | for (__p = 0; __ub != 1; __ub >>= 1) | |
823 | ++__p; | |
824 | ||
825 | const result_type __two2p = result_type(1) << __p; | |
826 | const result_type __k = __t / __two2p; | |
827 | ||
828 | if (__k & 1) | |
829 | return (__k + 1) * __two2p - 1; | |
830 | ||
831 | if (__c >= __b) | |
832 | return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p) | |
833 | / __two2d, | |
834 | __u % __two2p, __d); | |
835 | else | |
836 | return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p) | |
837 | / __two2d, | |
838 | __t % __two2p, __d); | |
839 | } | |
840 | ||
841 | template<class _UniformRandomNumberGenerator1, int __s1, | |
842 | class _UniformRandomNumberGenerator2, int __s2, | |
843 | typename _CharT, typename _Traits> | |
844 | std::basic_ostream<_CharT, _Traits>& | |
845 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
846 | const xor_combine<_UniformRandomNumberGenerator1, __s1, | |
847 | _UniformRandomNumberGenerator2, __s2>& __x) | |
848 | { | |
849 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
850 | typedef typename __ostream_type::ios_base __ios_base; | |
851 | ||
852 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
853 | const _CharT __fill = __os.fill(); | |
854 | const _CharT __space = __os.widen(' '); | |
855 | __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); | |
856 | __os.fill(__space); | |
857 | ||
858 | __os << __x.base1() << __space << __x.base2(); | |
859 | ||
860 | __os.flags(__flags); | |
861 | __os.fill(__fill); | |
862 | return __os; | |
863 | } | |
864 | ||
865 | template<class _UniformRandomNumberGenerator1, int __s1, | |
866 | class _UniformRandomNumberGenerator2, int __s2, | |
867 | typename _CharT, typename _Traits> | |
868 | std::basic_istream<_CharT, _Traits>& | |
869 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
870 | xor_combine<_UniformRandomNumberGenerator1, __s1, | |
871 | _UniformRandomNumberGenerator2, __s2>& __x) | |
872 | { | |
873 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
874 | typedef typename __istream_type::ios_base __ios_base; | |
875 | ||
876 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
877 | __is.flags(__ios_base::skipws); | |
878 | ||
879 | __is >> __x._M_b1 >> __x._M_b2; | |
880 | ||
881 | __is.flags(__flags); | |
882 | return __is; | |
883 | } | |
884 | ||
885 | ||
886 | template<typename _IntType> | |
887 | template<typename _UniformRandomNumberGenerator> | |
888 | typename uniform_int<_IntType>::result_type | |
889 | uniform_int<_IntType>:: | |
890 | _M_call(_UniformRandomNumberGenerator& __urng, | |
891 | result_type __min, result_type __max, true_type) | |
892 | { | |
893 | // XXX Must be fixed to work well for *arbitrary* __urng.max(), | |
894 | // __urng.min(), __max, __min. Currently works fine only in the | |
895 | // most common case __urng.max() - __urng.min() >= __max - __min, | |
896 | // with __urng.max() > __urng.min() >= 0. | |
897 | typedef typename __gnu_cxx::__add_unsigned<typename | |
898 | _UniformRandomNumberGenerator::result_type>::__type __urntype; | |
899 | typedef typename __gnu_cxx::__add_unsigned<result_type>::__type | |
900 | __utype; | |
901 | typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype) | |
902 | > sizeof(__utype)), | |
903 | __urntype, __utype>::__type __uctype; | |
904 | ||
905 | result_type __ret; | |
906 | ||
907 | const __urntype __urnmin = __urng.min(); | |
908 | const __urntype __urnmax = __urng.max(); | |
909 | const __urntype __urnrange = __urnmax - __urnmin; | |
910 | const __uctype __urange = __max - __min; | |
911 | const __uctype __udenom = (__urnrange <= __urange | |
912 | ? 1 : __urnrange / (__urange + 1)); | |
913 | do | |
914 | __ret = (__urntype(__urng()) - __urnmin) / __udenom; | |
915 | while (__ret > __max - __min); | |
916 | ||
917 | return __ret + __min; | |
918 | } | |
919 | ||
920 | template<typename _IntType, typename _CharT, typename _Traits> | |
921 | std::basic_ostream<_CharT, _Traits>& | |
922 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
923 | const uniform_int<_IntType>& __x) | |
924 | { | |
925 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
926 | typedef typename __ostream_type::ios_base __ios_base; | |
927 | ||
928 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
929 | const _CharT __fill = __os.fill(); | |
930 | const _CharT __space = __os.widen(' '); | |
931 | __os.flags(__ios_base::scientific | __ios_base::left); | |
932 | __os.fill(__space); | |
933 | ||
934 | __os << __x.min() << __space << __x.max(); | |
935 | ||
936 | __os.flags(__flags); | |
937 | __os.fill(__fill); | |
938 | return __os; | |
939 | } | |
940 | ||
941 | template<typename _IntType, typename _CharT, typename _Traits> | |
942 | std::basic_istream<_CharT, _Traits>& | |
943 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
944 | uniform_int<_IntType>& __x) | |
945 | { | |
946 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
947 | typedef typename __istream_type::ios_base __ios_base; | |
948 | ||
949 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
950 | __is.flags(__ios_base::dec | __ios_base::skipws); | |
951 | ||
952 | __is >> __x._M_min >> __x._M_max; | |
953 | ||
954 | __is.flags(__flags); | |
955 | return __is; | |
956 | } | |
957 | ||
958 | ||
959 | template<typename _CharT, typename _Traits> | |
960 | std::basic_ostream<_CharT, _Traits>& | |
961 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
962 | const bernoulli_distribution& __x) | |
963 | { | |
964 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
965 | typedef typename __ostream_type::ios_base __ios_base; | |
966 | ||
967 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
968 | const _CharT __fill = __os.fill(); | |
969 | const std::streamsize __precision = __os.precision(); | |
970 | __os.flags(__ios_base::scientific | __ios_base::left); | |
971 | __os.fill(__os.widen(' ')); | |
972 | __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10); | |
973 | ||
974 | __os << __x.p(); | |
975 | ||
976 | __os.flags(__flags); | |
977 | __os.fill(__fill); | |
978 | __os.precision(__precision); | |
979 | return __os; | |
980 | } | |
981 | ||
982 | ||
983 | template<typename _IntType, typename _RealType> | |
984 | template<class _UniformRandomNumberGenerator> | |
985 | typename geometric_distribution<_IntType, _RealType>::result_type | |
986 | geometric_distribution<_IntType, _RealType>:: | |
987 | operator()(_UniformRandomNumberGenerator& __urng) | |
988 | { | |
989 | // About the epsilon thing see this thread: | |
990 | // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html | |
991 | const _RealType __naf = | |
992 | (1 - std::numeric_limits<_RealType>::epsilon()) / 2; | |
993 | // The largest _RealType convertible to _IntType. | |
994 | const _RealType __thr = | |
995 | std::numeric_limits<_IntType>::max() + __naf; | |
996 | ||
997 | _RealType __cand; | |
998 | do | |
999 | __cand = std::ceil(std::log(__urng()) / _M_log_p); | |
1000 | while (__cand >= __thr); | |
1001 | ||
1002 | return result_type(__cand + __naf); | |
1003 | } | |
1004 | ||
1005 | template<typename _IntType, typename _RealType, | |
1006 | typename _CharT, typename _Traits> | |
1007 | std::basic_ostream<_CharT, _Traits>& | |
1008 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
1009 | const geometric_distribution<_IntType, _RealType>& __x) | |
1010 | { | |
1011 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
1012 | typedef typename __ostream_type::ios_base __ios_base; | |
1013 | ||
1014 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
1015 | const _CharT __fill = __os.fill(); | |
1016 | const std::streamsize __precision = __os.precision(); | |
1017 | __os.flags(__ios_base::scientific | __ios_base::left); | |
1018 | __os.fill(__os.widen(' ')); | |
1019 | __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); | |
1020 | ||
1021 | __os << __x.p(); | |
1022 | ||
1023 | __os.flags(__flags); | |
1024 | __os.fill(__fill); | |
1025 | __os.precision(__precision); | |
1026 | return __os; | |
1027 | } | |
1028 | ||
1029 | ||
1030 | template<typename _IntType, typename _RealType> | |
1031 | void | |
1032 | poisson_distribution<_IntType, _RealType>:: | |
1033 | _M_initialize() | |
1034 | { | |
1035 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
1036 | if (_M_mean >= 12) | |
1037 | { | |
1038 | const _RealType __m = std::floor(_M_mean); | |
1039 | _M_lm_thr = std::log(_M_mean); | |
95fe602e | 1040 | _M_lfm = std::tr1::lgamma(__m + 1); |
8e79468d BK |
1041 | _M_sm = std::sqrt(__m); |
1042 | ||
1043 | const _RealType __pi_4 = 0.7853981633974483096156608458198757L; | |
1044 | const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m | |
1045 | / __pi_4)); | |
95fe602e PC |
1046 | _M_d = std::tr1::round(std::max(_RealType(6), |
1047 | std::min(__m, __dx))); | |
8e79468d BK |
1048 | const _RealType __cx = 2 * __m + _M_d; |
1049 | _M_scx = std::sqrt(__cx / 2); | |
1050 | _M_1cx = 1 / __cx; | |
1051 | ||
1052 | _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); | |
1053 | _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d; | |
1054 | } | |
1055 | else | |
1056 | #endif | |
1057 | _M_lm_thr = std::exp(-_M_mean); | |
1058 | } | |
1059 | ||
1060 | /** | |
1061 | * A rejection algorithm when mean >= 12 and a simple method based | |
1062 | * upon the multiplication of uniform random variates otherwise. | |
1063 | * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 | |
1064 | * is defined. | |
1065 | * | |
1066 | * Reference: | |
2a60a9f6 | 1067 | * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, |
8e79468d BK |
1068 | * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). |
1069 | */ | |
1070 | template<typename _IntType, typename _RealType> | |
1071 | template<class _UniformRandomNumberGenerator> | |
1072 | typename poisson_distribution<_IntType, _RealType>::result_type | |
1073 | poisson_distribution<_IntType, _RealType>:: | |
1074 | operator()(_UniformRandomNumberGenerator& __urng) | |
1075 | { | |
1076 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
1077 | if (_M_mean >= 12) | |
1078 | { | |
1079 | _RealType __x; | |
1080 | ||
1081 | // See comments above... | |
1082 | const _RealType __naf = | |
1083 | (1 - std::numeric_limits<_RealType>::epsilon()) / 2; | |
1084 | const _RealType __thr = | |
1085 | std::numeric_limits<_IntType>::max() + __naf; | |
1086 | ||
1087 | const _RealType __m = std::floor(_M_mean); | |
1088 | // sqrt(pi / 2) | |
1089 | const _RealType __spi_2 = 1.2533141373155002512078826424055226L; | |
1090 | const _RealType __c1 = _M_sm * __spi_2; | |
1091 | const _RealType __c2 = _M_c2b + __c1; | |
1092 | const _RealType __c3 = __c2 + 1; | |
1093 | const _RealType __c4 = __c3 + 1; | |
1094 | // e^(1 / 78) | |
1095 | const _RealType __e178 = 1.0129030479320018583185514777512983L; | |
1096 | const _RealType __c5 = __c4 + __e178; | |
1097 | const _RealType __c = _M_cb + __c5; | |
1098 | const _RealType __2cx = 2 * (2 * __m + _M_d); | |
1099 | ||
1100 | bool __reject = true; | |
1101 | do | |
1102 | { | |
1103 | const _RealType __u = __c * __urng(); | |
1104 | const _RealType __e = -std::log(__urng()); | |
1105 | ||
1106 | _RealType __w = 0.0; | |
1107 | ||
1108 | if (__u <= __c1) | |
1109 | { | |
1110 | const _RealType __n = _M_nd(__urng); | |
1111 | const _RealType __y = -std::abs(__n) * _M_sm - 1; | |
1112 | __x = std::floor(__y); | |
1113 | __w = -__n * __n / 2; | |
1114 | if (__x < -__m) | |
1115 | continue; | |
1116 | } | |
1117 | else if (__u <= __c2) | |
1118 | { | |
1119 | const _RealType __n = _M_nd(__urng); | |
1120 | const _RealType __y = 1 + std::abs(__n) * _M_scx; | |
1121 | __x = std::ceil(__y); | |
1122 | __w = __y * (2 - __y) * _M_1cx; | |
1123 | if (__x > _M_d) | |
1124 | continue; | |
1125 | } | |
1126 | else if (__u <= __c3) | |
1127 | // NB: This case not in the book, nor in the Errata, | |
1128 | // but should be ok... | |
1129 | __x = -1; | |
1130 | else if (__u <= __c4) | |
1131 | __x = 0; | |
1132 | else if (__u <= __c5) | |
1133 | __x = 1; | |
1134 | else | |
1135 | { | |
1136 | const _RealType __v = -std::log(__urng()); | |
1137 | const _RealType __y = _M_d + __v * __2cx / _M_d; | |
1138 | __x = std::ceil(__y); | |
1139 | __w = -_M_d * _M_1cx * (1 + __y / 2); | |
1140 | } | |
1141 | ||
1142 | __reject = (__w - __e - __x * _M_lm_thr | |
95fe602e | 1143 | > _M_lfm - std::tr1::lgamma(__x + __m + 1)); |
8e79468d BK |
1144 | |
1145 | __reject |= __x + __m >= __thr; | |
1146 | ||
1147 | } while (__reject); | |
1148 | ||
1149 | return result_type(__x + __m + __naf); | |
1150 | } | |
1151 | else | |
1152 | #endif | |
1153 | { | |
1154 | _IntType __x = 0; | |
1155 | _RealType __prod = 1.0; | |
1156 | ||
1157 | do | |
1158 | { | |
1159 | __prod *= __urng(); | |
1160 | __x += 1; | |
1161 | } | |
1162 | while (__prod > _M_lm_thr); | |
1163 | ||
1164 | return __x - 1; | |
1165 | } | |
1166 | } | |
1167 | ||
1168 | template<typename _IntType, typename _RealType, | |
1169 | typename _CharT, typename _Traits> | |
1170 | std::basic_ostream<_CharT, _Traits>& | |
1171 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
1172 | const poisson_distribution<_IntType, _RealType>& __x) | |
1173 | { | |
1174 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
1175 | typedef typename __ostream_type::ios_base __ios_base; | |
1176 | ||
1177 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
1178 | const _CharT __fill = __os.fill(); | |
1179 | const std::streamsize __precision = __os.precision(); | |
1180 | const _CharT __space = __os.widen(' '); | |
1181 | __os.flags(__ios_base::scientific | __ios_base::left); | |
1182 | __os.fill(__space); | |
1183 | __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); | |
1184 | ||
1185 | __os << __x.mean() << __space << __x._M_nd; | |
1186 | ||
1187 | __os.flags(__flags); | |
1188 | __os.fill(__fill); | |
1189 | __os.precision(__precision); | |
1190 | return __os; | |
1191 | } | |
1192 | ||
1193 | template<typename _IntType, typename _RealType, | |
1194 | typename _CharT, typename _Traits> | |
1195 | std::basic_istream<_CharT, _Traits>& | |
1196 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
1197 | poisson_distribution<_IntType, _RealType>& __x) | |
1198 | { | |
1199 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
1200 | typedef typename __istream_type::ios_base __ios_base; | |
1201 | ||
1202 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
1203 | __is.flags(__ios_base::skipws); | |
1204 | ||
1205 | __is >> __x._M_mean >> __x._M_nd; | |
1206 | __x._M_initialize(); | |
1207 | ||
1208 | __is.flags(__flags); | |
1209 | return __is; | |
1210 | } | |
1211 | ||
1212 | ||
1213 | template<typename _IntType, typename _RealType> | |
1214 | void | |
1215 | binomial_distribution<_IntType, _RealType>:: | |
1216 | _M_initialize() | |
1217 | { | |
1218 | const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; | |
1219 | ||
1220 | _M_easy = true; | |
1221 | ||
1222 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
1223 | if (_M_t * __p12 >= 8) | |
1224 | { | |
1225 | _M_easy = false; | |
1226 | const _RealType __np = std::floor(_M_t * __p12); | |
1227 | const _RealType __pa = __np / _M_t; | |
1228 | const _RealType __1p = 1 - __pa; | |
1229 | ||
1230 | const _RealType __pi_4 = 0.7853981633974483096156608458198757L; | |
1231 | const _RealType __d1x = | |
1232 | std::sqrt(__np * __1p * std::log(32 * __np | |
1233 | / (81 * __pi_4 * __1p))); | |
95fe602e | 1234 | _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x)); |
8e79468d BK |
1235 | const _RealType __d2x = |
1236 | std::sqrt(__np * __1p * std::log(32 * _M_t * __1p | |
1237 | / (__pi_4 * __pa))); | |
95fe602e | 1238 | _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x)); |
8e79468d BK |
1239 | |
1240 | // sqrt(pi / 2) | |
1241 | const _RealType __spi_2 = 1.2533141373155002512078826424055226L; | |
1242 | _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); | |
1243 | _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); | |
1244 | _M_c = 2 * _M_d1 / __np; | |
1245 | _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; | |
1246 | const _RealType __a12 = _M_a1 + _M_s2 * __spi_2; | |
1247 | const _RealType __s1s = _M_s1 * _M_s1; | |
1248 | _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) | |
1249 | * 2 * __s1s / _M_d1 | |
1250 | * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); | |
1251 | const _RealType __s2s = _M_s2 * _M_s2; | |
1252 | _M_s = (_M_a123 + 2 * __s2s / _M_d2 | |
1253 | * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); | |
95fe602e PC |
1254 | _M_lf = (std::tr1::lgamma(__np + 1) |
1255 | + std::tr1::lgamma(_M_t - __np + 1)); | |
8e79468d BK |
1256 | _M_lp1p = std::log(__pa / __1p); |
1257 | ||
1258 | _M_q = -std::log(1 - (__p12 - __pa) / __1p); | |
1259 | } | |
1260 | else | |
1261 | #endif | |
1262 | _M_q = -std::log(1 - __p12); | |
1263 | } | |
1264 | ||
1265 | template<typename _IntType, typename _RealType> | |
1266 | template<class _UniformRandomNumberGenerator> | |
1267 | typename binomial_distribution<_IntType, _RealType>::result_type | |
1268 | binomial_distribution<_IntType, _RealType>:: | |
1269 | _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) | |
1270 | { | |
1271 | _IntType __x = 0; | |
1272 | _RealType __sum = 0; | |
1273 | ||
1274 | do | |
1275 | { | |
1276 | const _RealType __e = -std::log(__urng()); | |
1277 | __sum += __e / (__t - __x); | |
1278 | __x += 1; | |
1279 | } | |
1280 | while (__sum <= _M_q); | |
1281 | ||
1282 | return __x - 1; | |
1283 | } | |
1284 | ||
1285 | /** | |
1286 | * A rejection algorithm when t * p >= 8 and a simple waiting time | |
1287 | * method - the second in the referenced book - otherwise. | |
1288 | * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 | |
1289 | * is defined. | |
1290 | * | |
1291 | * Reference: | |
2a60a9f6 | 1292 | * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, |
8e79468d BK |
1293 | * New York, 1986, Ch. X, Sect. 4 (+ Errata!). |
1294 | */ | |
1295 | template<typename _IntType, typename _RealType> | |
1296 | template<class _UniformRandomNumberGenerator> | |
1297 | typename binomial_distribution<_IntType, _RealType>::result_type | |
1298 | binomial_distribution<_IntType, _RealType>:: | |
1299 | operator()(_UniformRandomNumberGenerator& __urng) | |
1300 | { | |
1301 | result_type __ret; | |
1302 | const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; | |
1303 | ||
1304 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
1305 | if (!_M_easy) | |
1306 | { | |
1307 | _RealType __x; | |
1308 | ||
1309 | // See comments above... | |
1310 | const _RealType __naf = | |
1311 | (1 - std::numeric_limits<_RealType>::epsilon()) / 2; | |
1312 | const _RealType __thr = | |
1313 | std::numeric_limits<_IntType>::max() + __naf; | |
1314 | ||
1315 | const _RealType __np = std::floor(_M_t * __p12); | |
1316 | const _RealType __pa = __np / _M_t; | |
1317 | ||
1318 | // sqrt(pi / 2) | |
1319 | const _RealType __spi_2 = 1.2533141373155002512078826424055226L; | |
1320 | const _RealType __a1 = _M_a1; | |
1321 | const _RealType __a12 = __a1 + _M_s2 * __spi_2; | |
1322 | const _RealType __a123 = _M_a123; | |
1323 | const _RealType __s1s = _M_s1 * _M_s1; | |
1324 | const _RealType __s2s = _M_s2 * _M_s2; | |
1325 | ||
1326 | bool __reject; | |
1327 | do | |
1328 | { | |
1329 | const _RealType __u = _M_s * __urng(); | |
1330 | ||
1331 | _RealType __v; | |
1332 | ||
1333 | if (__u <= __a1) | |
1334 | { | |
1335 | const _RealType __n = _M_nd(__urng); | |
1336 | const _RealType __y = _M_s1 * std::abs(__n); | |
1337 | __reject = __y >= _M_d1; | |
1338 | if (!__reject) | |
1339 | { | |
1340 | const _RealType __e = -std::log(__urng()); | |
1341 | __x = std::floor(__y); | |
1342 | __v = -__e - __n * __n / 2 + _M_c; | |
1343 | } | |
1344 | } | |
1345 | else if (__u <= __a12) | |
1346 | { | |
1347 | const _RealType __n = _M_nd(__urng); | |
1348 | const _RealType __y = _M_s2 * std::abs(__n); | |
1349 | __reject = __y >= _M_d2; | |
1350 | if (!__reject) | |
1351 | { | |
1352 | const _RealType __e = -std::log(__urng()); | |
1353 | __x = std::floor(-__y); | |
1354 | __v = -__e - __n * __n / 2; | |
1355 | } | |
1356 | } | |
1357 | else if (__u <= __a123) | |
1358 | { | |
1359 | const _RealType __e1 = -std::log(__urng()); | |
1360 | const _RealType __e2 = -std::log(__urng()); | |
1361 | ||
1362 | const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1; | |
1363 | __x = std::floor(__y); | |
1364 | __v = (-__e2 + _M_d1 * (1 / (_M_t - __np) | |
1365 | -__y / (2 * __s1s))); | |
1366 | __reject = false; | |
1367 | } | |
1368 | else | |
1369 | { | |
1370 | const _RealType __e1 = -std::log(__urng()); | |
1371 | const _RealType __e2 = -std::log(__urng()); | |
1372 | ||
1373 | const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2; | |
1374 | __x = std::floor(-__y); | |
1375 | __v = -__e2 - _M_d2 * __y / (2 * __s2s); | |
1376 | __reject = false; | |
1377 | } | |
1378 | ||
1379 | __reject = __reject || __x < -__np || __x > _M_t - __np; | |
1380 | if (!__reject) | |
1381 | { | |
1382 | const _RealType __lfx = | |
95fe602e PC |
1383 | std::tr1::lgamma(__np + __x + 1) |
1384 | + std::tr1::lgamma(_M_t - (__np + __x) + 1); | |
8e79468d BK |
1385 | __reject = __v > _M_lf - __lfx + __x * _M_lp1p; |
1386 | } | |
1387 | ||
1388 | __reject |= __x + __np >= __thr; | |
1389 | } | |
1390 | while (__reject); | |
1391 | ||
1392 | __x += __np + __naf; | |
1393 | ||
1394 | const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); | |
1395 | __ret = _IntType(__x) + __z; | |
1396 | } | |
1397 | else | |
1398 | #endif | |
1399 | __ret = _M_waiting(__urng, _M_t); | |
1400 | ||
1401 | if (__p12 != _M_p) | |
1402 | __ret = _M_t - __ret; | |
1403 | return __ret; | |
1404 | } | |
1405 | ||
1406 | template<typename _IntType, typename _RealType, | |
1407 | typename _CharT, typename _Traits> | |
1408 | std::basic_ostream<_CharT, _Traits>& | |
1409 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
1410 | const binomial_distribution<_IntType, _RealType>& __x) | |
1411 | { | |
1412 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
1413 | typedef typename __ostream_type::ios_base __ios_base; | |
1414 | ||
1415 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
1416 | const _CharT __fill = __os.fill(); | |
1417 | const std::streamsize __precision = __os.precision(); | |
1418 | const _CharT __space = __os.widen(' '); | |
1419 | __os.flags(__ios_base::scientific | __ios_base::left); | |
1420 | __os.fill(__space); | |
1421 | __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); | |
1422 | ||
1423 | __os << __x.t() << __space << __x.p() | |
1424 | << __space << __x._M_nd; | |
1425 | ||
1426 | __os.flags(__flags); | |
1427 | __os.fill(__fill); | |
1428 | __os.precision(__precision); | |
1429 | return __os; | |
1430 | } | |
1431 | ||
1432 | template<typename _IntType, typename _RealType, | |
1433 | typename _CharT, typename _Traits> | |
1434 | std::basic_istream<_CharT, _Traits>& | |
1435 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
1436 | binomial_distribution<_IntType, _RealType>& __x) | |
1437 | { | |
1438 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
1439 | typedef typename __istream_type::ios_base __ios_base; | |
1440 | ||
1441 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
1442 | __is.flags(__ios_base::dec | __ios_base::skipws); | |
1443 | ||
1444 | __is >> __x._M_t >> __x._M_p >> __x._M_nd; | |
1445 | __x._M_initialize(); | |
1446 | ||
1447 | __is.flags(__flags); | |
1448 | return __is; | |
1449 | } | |
1450 | ||
1451 | ||
1452 | template<typename _RealType, typename _CharT, typename _Traits> | |
1453 | std::basic_ostream<_CharT, _Traits>& | |
1454 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
1455 | const uniform_real<_RealType>& __x) | |
1456 | { | |
1457 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
1458 | typedef typename __ostream_type::ios_base __ios_base; | |
1459 | ||
1460 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
1461 | const _CharT __fill = __os.fill(); | |
1462 | const std::streamsize __precision = __os.precision(); | |
1463 | const _CharT __space = __os.widen(' '); | |
1464 | __os.flags(__ios_base::scientific | __ios_base::left); | |
1465 | __os.fill(__space); | |
1466 | __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); | |
1467 | ||
1468 | __os << __x.min() << __space << __x.max(); | |
1469 | ||
1470 | __os.flags(__flags); | |
1471 | __os.fill(__fill); | |
1472 | __os.precision(__precision); | |
1473 | return __os; | |
1474 | } | |
1475 | ||
1476 | template<typename _RealType, typename _CharT, typename _Traits> | |
1477 | std::basic_istream<_CharT, _Traits>& | |
1478 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
1479 | uniform_real<_RealType>& __x) | |
1480 | { | |
1481 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
1482 | typedef typename __istream_type::ios_base __ios_base; | |
1483 | ||
1484 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
1485 | __is.flags(__ios_base::skipws); | |
1486 | ||
1487 | __is >> __x._M_min >> __x._M_max; | |
1488 | ||
1489 | __is.flags(__flags); | |
1490 | return __is; | |
1491 | } | |
1492 | ||
1493 | ||
1494 | template<typename _RealType, typename _CharT, typename _Traits> | |
1495 | std::basic_ostream<_CharT, _Traits>& | |
1496 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
1497 | const exponential_distribution<_RealType>& __x) | |
1498 | { | |
1499 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
1500 | typedef typename __ostream_type::ios_base __ios_base; | |
1501 | ||
1502 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
1503 | const _CharT __fill = __os.fill(); | |
1504 | const std::streamsize __precision = __os.precision(); | |
1505 | __os.flags(__ios_base::scientific | __ios_base::left); | |
1506 | __os.fill(__os.widen(' ')); | |
1507 | __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); | |
1508 | ||
1509 | __os << __x.lambda(); | |
1510 | ||
1511 | __os.flags(__flags); | |
1512 | __os.fill(__fill); | |
1513 | __os.precision(__precision); | |
1514 | return __os; | |
1515 | } | |
1516 | ||
1517 | ||
1518 | /** | |
1519 | * Polar method due to Marsaglia. | |
1520 | * | |
2a60a9f6 | 1521 | * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, |
8e79468d BK |
1522 | * New York, 1986, Ch. V, Sect. 4.4. |
1523 | */ | |
1524 | template<typename _RealType> | |
1525 | template<class _UniformRandomNumberGenerator> | |
1526 | typename normal_distribution<_RealType>::result_type | |
1527 | normal_distribution<_RealType>:: | |
1528 | operator()(_UniformRandomNumberGenerator& __urng) | |
1529 | { | |
1530 | result_type __ret; | |
1531 | ||
1532 | if (_M_saved_available) | |
1533 | { | |
1534 | _M_saved_available = false; | |
1535 | __ret = _M_saved; | |
1536 | } | |
1537 | else | |
1538 | { | |
1539 | result_type __x, __y, __r2; | |
1540 | do | |
1541 | { | |
1542 | __x = result_type(2.0) * __urng() - 1.0; | |
1543 | __y = result_type(2.0) * __urng() - 1.0; | |
1544 | __r2 = __x * __x + __y * __y; | |
1545 | } | |
1546 | while (__r2 > 1.0 || __r2 == 0.0); | |
1547 | ||
1548 | const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); | |
1549 | _M_saved = __x * __mult; | |
1550 | _M_saved_available = true; | |
1551 | __ret = __y * __mult; | |
1552 | } | |
1553 | ||
1554 | __ret = __ret * _M_sigma + _M_mean; | |
1555 | return __ret; | |
1556 | } | |
1557 | ||
1558 | template<typename _RealType, typename _CharT, typename _Traits> | |
1559 | std::basic_ostream<_CharT, _Traits>& | |
1560 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
1561 | const normal_distribution<_RealType>& __x) | |
1562 | { | |
1563 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
1564 | typedef typename __ostream_type::ios_base __ios_base; | |
1565 | ||
1566 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
1567 | const _CharT __fill = __os.fill(); | |
1568 | const std::streamsize __precision = __os.precision(); | |
1569 | const _CharT __space = __os.widen(' '); | |
1570 | __os.flags(__ios_base::scientific | __ios_base::left); | |
1571 | __os.fill(__space); | |
1572 | __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); | |
1573 | ||
1574 | __os << __x._M_saved_available << __space | |
1575 | << __x.mean() << __space | |
1576 | << __x.sigma(); | |
1577 | if (__x._M_saved_available) | |
1578 | __os << __space << __x._M_saved; | |
1579 | ||
1580 | __os.flags(__flags); | |
1581 | __os.fill(__fill); | |
1582 | __os.precision(__precision); | |
1583 | return __os; | |
1584 | } | |
1585 | ||
1586 | template<typename _RealType, typename _CharT, typename _Traits> | |
1587 | std::basic_istream<_CharT, _Traits>& | |
1588 | operator>>(std::basic_istream<_CharT, _Traits>& __is, | |
1589 | normal_distribution<_RealType>& __x) | |
1590 | { | |
1591 | typedef std::basic_istream<_CharT, _Traits> __istream_type; | |
1592 | typedef typename __istream_type::ios_base __ios_base; | |
1593 | ||
1594 | const typename __ios_base::fmtflags __flags = __is.flags(); | |
1595 | __is.flags(__ios_base::dec | __ios_base::skipws); | |
1596 | ||
1597 | __is >> __x._M_saved_available >> __x._M_mean | |
1598 | >> __x._M_sigma; | |
1599 | if (__x._M_saved_available) | |
1600 | __is >> __x._M_saved; | |
1601 | ||
1602 | __is.flags(__flags); | |
1603 | return __is; | |
1604 | } | |
1605 | ||
1606 | ||
1607 | template<typename _RealType> | |
1608 | void | |
1609 | gamma_distribution<_RealType>:: | |
1610 | _M_initialize() | |
1611 | { | |
1612 | if (_M_alpha >= 1) | |
1613 | _M_l_d = std::sqrt(2 * _M_alpha - 1); | |
1614 | else | |
1615 | _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha)) | |
1616 | * (1 - _M_alpha)); | |
1617 | } | |
1618 | ||
1619 | /** | |
1620 | * Cheng's rejection algorithm GB for alpha >= 1 and a modification | |
1621 | * of Vaduva's rejection from Weibull algorithm due to Devroye for | |
1622 | * alpha < 1. | |
1623 | * | |
1624 | * References: | |
2a60a9f6 BK |
1625 | * Cheng, R. C. The Generation of Gamma Random Variables with Non-integral |
1626 | * Shape Parameter. Applied Statistics, 26, 71-75, 1977. | |
8e79468d | 1627 | * |
2a60a9f6 BK |
1628 | * Vaduva, I. Computer Generation of Gamma Gandom Variables by Rejection |
1629 | * and Composition Procedures. Math. Operationsforschung and Statistik, | |
8e79468d BK |
1630 | * Series in Statistics, 8, 545-576, 1977. |
1631 | * | |
2a60a9f6 | 1632 | * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, |
8e79468d BK |
1633 | * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!). |
1634 | */ | |
1635 | template<typename _RealType> | |
1636 | template<class _UniformRandomNumberGenerator> | |
1637 | typename gamma_distribution<_RealType>::result_type | |
1638 | gamma_distribution<_RealType>:: | |
1639 | operator()(_UniformRandomNumberGenerator& __urng) | |
1640 | { | |
1641 | result_type __x; | |
1642 | ||
1643 | bool __reject; | |
1644 | if (_M_alpha >= 1) | |
1645 | { | |
1646 | // alpha - log(4) | |
1647 | const result_type __b = _M_alpha | |
1648 | - result_type(1.3862943611198906188344642429163531L); | |
1649 | const result_type __c = _M_alpha + _M_l_d; | |
1650 | const result_type __1l = 1 / _M_l_d; | |
1651 | ||
1652 | // 1 + log(9 / 2) | |
1653 | const result_type __k = 2.5040773967762740733732583523868748L; | |
1654 | ||
1655 | do | |
1656 | { | |
1657 | const result_type __u = __urng(); | |
1658 | const result_type __v = __urng(); | |
1659 | ||
1660 | const result_type __y = __1l * std::log(__v / (1 - __v)); | |
1661 | __x = _M_alpha * std::exp(__y); | |
1662 | ||
1663 | const result_type __z = __u * __v * __v; | |
1664 | const result_type __r = __b + __c * __y - __x; | |
1665 | ||
1666 | __reject = __r < result_type(4.5) * __z - __k; | |
1667 | if (__reject) | |
1668 | __reject = __r < std::log(__z); | |
1669 | } | |
1670 | while (__reject); | |
1671 | } | |
1672 | else | |
1673 | { | |
1674 | const result_type __c = 1 / _M_alpha; | |
1675 | ||
1676 | do | |
1677 | { | |
1678 | const result_type __z = -std::log(__urng()); | |
1679 | const result_type __e = -std::log(__urng()); | |
1680 | ||
1681 | __x = std::pow(__z, __c); | |
1682 | ||
1683 | __reject = __z + __e < _M_l_d + __x; | |
1684 | } | |
1685 | while (__reject); | |
1686 | } | |
1687 | ||
1688 | return __x; | |
1689 | } | |
1690 | ||
1691 | template<typename _RealType, typename _CharT, typename _Traits> | |
1692 | std::basic_ostream<_CharT, _Traits>& | |
1693 | operator<<(std::basic_ostream<_CharT, _Traits>& __os, | |
1694 | const gamma_distribution<_RealType>& __x) | |
1695 | { | |
1696 | typedef std::basic_ostream<_CharT, _Traits> __ostream_type; | |
1697 | typedef typename __ostream_type::ios_base __ios_base; | |
1698 | ||
1699 | const typename __ios_base::fmtflags __flags = __os.flags(); | |
1700 | const _CharT __fill = __os.fill(); | |
1701 | const std::streamsize __precision = __os.precision(); | |
1702 | __os.flags(__ios_base::scientific | __ios_base::left); | |
1703 | __os.fill(__os.widen(' ')); | |
1704 | __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); | |
1705 | ||
1706 | __os << __x.alpha(); | |
1707 | ||
1708 | __os.flags(__flags); | |
1709 | __os.fill(__fill); | |
1710 | __os.precision(__precision); | |
1711 | return __os; | |
1712 | } | |
4a15d842 | 1713 | } |
12ffa228 BK |
1714 | |
1715 | _GLIBCXX_END_NAMESPACE_VERSION | |
8e79468d | 1716 | } |
06f29237 PC |
1717 | |
1718 | #endif |