libstdc++
|
00001 // random number generation (out of line) -*- C++ -*- 00002 00003 // Copyright (C) 2007, 2009 Free Software Foundation, Inc. 00004 // 00005 // This file is part of the GNU ISO C++ Library. This library is free 00006 // software; you can redistribute it and/or modify it under the 00007 // terms of the GNU General Public License as published by the 00008 // Free Software Foundation; either version 3, or (at your option) 00009 // any later version. 00010 00011 // This library is distributed in the hope that it will be useful, 00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 // GNU General Public License for more details. 00015 00016 // Under Section 7 of GPL version 3, you are granted additional 00017 // permissions described in the GCC Runtime Library Exception, version 00018 // 3.1, as published by the Free Software Foundation. 00019 00020 // You should have received a copy of the GNU General Public License and 00021 // a copy of the GCC Runtime Library Exception along with this program; 00022 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 00023 // <http://www.gnu.org/licenses/>. 00024 00025 /** @file tr1_impl/random.tcc 00026 * This is an internal header file, included by other library headers. 00027 * You should not attempt to use it directly. 00028 */ 00029 00030 namespace std 00031 { 00032 _GLIBCXX_BEGIN_NAMESPACE_TR1 00033 00034 /* 00035 * (Further) implementation-space details. 00036 */ 00037 namespace __detail 00038 { 00039 // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid 00040 // integer overflow. 00041 // 00042 // Because a and c are compile-time integral constants the compiler kindly 00043 // elides any unreachable paths. 00044 // 00045 // Preconditions: a > 0, m > 0. 00046 // 00047 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool> 00048 struct _Mod 00049 { 00050 static _Tp 00051 __calc(_Tp __x) 00052 { 00053 if (__a == 1) 00054 __x %= __m; 00055 else 00056 { 00057 static const _Tp __q = __m / __a; 00058 static const _Tp __r = __m % __a; 00059 00060 _Tp __t1 = __a * (__x % __q); 00061 _Tp __t2 = __r * (__x / __q); 00062 if (__t1 >= __t2) 00063 __x = __t1 - __t2; 00064 else 00065 __x = __m - __t2 + __t1; 00066 } 00067 00068 if (__c != 0) 00069 { 00070 const _Tp __d = __m - __x; 00071 if (__d > __c) 00072 __x += __c; 00073 else 00074 __x = __c - __d; 00075 } 00076 return __x; 00077 } 00078 }; 00079 00080 // Special case for m == 0 -- use unsigned integer overflow as modulo 00081 // operator. 00082 template<typename _Tp, _Tp __a, _Tp __c, _Tp __m> 00083 struct _Mod<_Tp, __a, __c, __m, true> 00084 { 00085 static _Tp 00086 __calc(_Tp __x) 00087 { return __a * __x + __c; } 00088 }; 00089 } // namespace __detail 00090 00091 /** 00092 * Seeds the LCR with integral value @p __x0, adjusted so that the 00093 * ring identity is never a member of the convergence set. 00094 */ 00095 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00096 void 00097 linear_congruential<_UIntType, __a, __c, __m>:: 00098 seed(unsigned long __x0) 00099 { 00100 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) 00101 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) 00102 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); 00103 else 00104 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); 00105 } 00106 00107 /** 00108 * Seeds the LCR engine with a value generated by @p __g. 00109 */ 00110 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00111 template<class _Gen> 00112 void 00113 linear_congruential<_UIntType, __a, __c, __m>:: 00114 seed(_Gen& __g, false_type) 00115 { 00116 _UIntType __x0 = __g(); 00117 if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) 00118 && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) 00119 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); 00120 else 00121 _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); 00122 } 00123 00124 /** 00125 * Gets the next generated value in sequence. 00126 */ 00127 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00128 typename linear_congruential<_UIntType, __a, __c, __m>::result_type 00129 linear_congruential<_UIntType, __a, __c, __m>:: 00130 operator()() 00131 { 00132 _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x); 00133 return _M_x; 00134 } 00135 00136 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 00137 typename _CharT, typename _Traits> 00138 std::basic_ostream<_CharT, _Traits>& 00139 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00140 const linear_congruential<_UIntType, __a, __c, __m>& __lcr) 00141 { 00142 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00143 typedef typename __ostream_type::ios_base __ios_base; 00144 00145 const typename __ios_base::fmtflags __flags = __os.flags(); 00146 const _CharT __fill = __os.fill(); 00147 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00148 __os.fill(__os.widen(' ')); 00149 00150 __os << __lcr._M_x; 00151 00152 __os.flags(__flags); 00153 __os.fill(__fill); 00154 return __os; 00155 } 00156 00157 template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 00158 typename _CharT, typename _Traits> 00159 std::basic_istream<_CharT, _Traits>& 00160 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00161 linear_congruential<_UIntType, __a, __c, __m>& __lcr) 00162 { 00163 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00164 typedef typename __istream_type::ios_base __ios_base; 00165 00166 const typename __ios_base::fmtflags __flags = __is.flags(); 00167 __is.flags(__ios_base::dec); 00168 00169 __is >> __lcr._M_x; 00170 00171 __is.flags(__flags); 00172 return __is; 00173 } 00174 00175 00176 template<class _UIntType, int __w, int __n, int __m, int __r, 00177 _UIntType __a, int __u, int __s, 00178 _UIntType __b, int __t, _UIntType __c, int __l> 00179 void 00180 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 00181 __b, __t, __c, __l>:: 00182 seed(unsigned long __value) 00183 { 00184 _M_x[0] = __detail::__mod<_UIntType, 1, 0, 00185 __detail::_Shift<_UIntType, __w>::__value>(__value); 00186 00187 for (int __i = 1; __i < state_size; ++__i) 00188 { 00189 _UIntType __x = _M_x[__i - 1]; 00190 __x ^= __x >> (__w - 2); 00191 __x *= 1812433253ul; 00192 __x += __i; 00193 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, 00194 __detail::_Shift<_UIntType, __w>::__value>(__x); 00195 } 00196 _M_p = state_size; 00197 } 00198 00199 template<class _UIntType, int __w, int __n, int __m, int __r, 00200 _UIntType __a, int __u, int __s, 00201 _UIntType __b, int __t, _UIntType __c, int __l> 00202 template<class _Gen> 00203 void 00204 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 00205 __b, __t, __c, __l>:: 00206 seed(_Gen& __gen, false_type) 00207 { 00208 for (int __i = 0; __i < state_size; ++__i) 00209 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, 00210 __detail::_Shift<_UIntType, __w>::__value>(__gen()); 00211 _M_p = state_size; 00212 } 00213 00214 template<class _UIntType, int __w, int __n, int __m, int __r, 00215 _UIntType __a, int __u, int __s, 00216 _UIntType __b, int __t, _UIntType __c, int __l> 00217 typename 00218 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 00219 __b, __t, __c, __l>::result_type 00220 mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 00221 __b, __t, __c, __l>:: 00222 operator()() 00223 { 00224 // Reload the vector - cost is O(n) amortized over n calls. 00225 if (_M_p >= state_size) 00226 { 00227 const _UIntType __upper_mask = (~_UIntType()) << __r; 00228 const _UIntType __lower_mask = ~__upper_mask; 00229 00230 for (int __k = 0; __k < (__n - __m); ++__k) 00231 { 00232 _UIntType __y = ((_M_x[__k] & __upper_mask) 00233 | (_M_x[__k + 1] & __lower_mask)); 00234 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) 00235 ^ ((__y & 0x01) ? __a : 0)); 00236 } 00237 00238 for (int __k = (__n - __m); __k < (__n - 1); ++__k) 00239 { 00240 _UIntType __y = ((_M_x[__k] & __upper_mask) 00241 | (_M_x[__k + 1] & __lower_mask)); 00242 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) 00243 ^ ((__y & 0x01) ? __a : 0)); 00244 } 00245 00246 _UIntType __y = ((_M_x[__n - 1] & __upper_mask) 00247 | (_M_x[0] & __lower_mask)); 00248 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) 00249 ^ ((__y & 0x01) ? __a : 0)); 00250 _M_p = 0; 00251 } 00252 00253 // Calculate o(x(i)). 00254 result_type __z = _M_x[_M_p++]; 00255 __z ^= (__z >> __u); 00256 __z ^= (__z << __s) & __b; 00257 __z ^= (__z << __t) & __c; 00258 __z ^= (__z >> __l); 00259 00260 return __z; 00261 } 00262 00263 template<class _UIntType, int __w, int __n, int __m, int __r, 00264 _UIntType __a, int __u, int __s, _UIntType __b, int __t, 00265 _UIntType __c, int __l, 00266 typename _CharT, typename _Traits> 00267 std::basic_ostream<_CharT, _Traits>& 00268 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00269 const mersenne_twister<_UIntType, __w, __n, __m, 00270 __r, __a, __u, __s, __b, __t, __c, __l>& __x) 00271 { 00272 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00273 typedef typename __ostream_type::ios_base __ios_base; 00274 00275 const typename __ios_base::fmtflags __flags = __os.flags(); 00276 const _CharT __fill = __os.fill(); 00277 const _CharT __space = __os.widen(' '); 00278 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00279 __os.fill(__space); 00280 00281 for (int __i = 0; __i < __n - 1; ++__i) 00282 __os << __x._M_x[__i] << __space; 00283 __os << __x._M_x[__n - 1]; 00284 00285 __os.flags(__flags); 00286 __os.fill(__fill); 00287 return __os; 00288 } 00289 00290 template<class _UIntType, int __w, int __n, int __m, int __r, 00291 _UIntType __a, int __u, int __s, _UIntType __b, int __t, 00292 _UIntType __c, int __l, 00293 typename _CharT, typename _Traits> 00294 std::basic_istream<_CharT, _Traits>& 00295 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00296 mersenne_twister<_UIntType, __w, __n, __m, 00297 __r, __a, __u, __s, __b, __t, __c, __l>& __x) 00298 { 00299 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00300 typedef typename __istream_type::ios_base __ios_base; 00301 00302 const typename __ios_base::fmtflags __flags = __is.flags(); 00303 __is.flags(__ios_base::dec | __ios_base::skipws); 00304 00305 for (int __i = 0; __i < __n; ++__i) 00306 __is >> __x._M_x[__i]; 00307 00308 __is.flags(__flags); 00309 return __is; 00310 } 00311 00312 00313 template<typename _IntType, _IntType __m, int __s, int __r> 00314 void 00315 subtract_with_carry<_IntType, __m, __s, __r>:: 00316 seed(unsigned long __value) 00317 { 00318 if (__value == 0) 00319 __value = 19780503; 00320 00321 std::_GLIBCXX_TR1 linear_congruential<unsigned long, 40014, 0, 2147483563> 00322 __lcg(__value); 00323 00324 for (int __i = 0; __i < long_lag; ++__i) 00325 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg()); 00326 00327 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 00328 _M_p = 0; 00329 } 00330 00331 template<typename _IntType, _IntType __m, int __s, int __r> 00332 template<class _Gen> 00333 void 00334 subtract_with_carry<_IntType, __m, __s, __r>:: 00335 seed(_Gen& __gen, false_type) 00336 { 00337 const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32; 00338 00339 for (int __i = 0; __i < long_lag; ++__i) 00340 { 00341 _UIntType __tmp = 0; 00342 _UIntType __factor = 1; 00343 for (int __j = 0; __j < __n; ++__j) 00344 { 00345 __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0> 00346 (__gen()) * __factor; 00347 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00348 } 00349 _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp); 00350 } 00351 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 00352 _M_p = 0; 00353 } 00354 00355 template<typename _IntType, _IntType __m, int __s, int __r> 00356 typename subtract_with_carry<_IntType, __m, __s, __r>::result_type 00357 subtract_with_carry<_IntType, __m, __s, __r>:: 00358 operator()() 00359 { 00360 // Derive short lag index from current index. 00361 int __ps = _M_p - short_lag; 00362 if (__ps < 0) 00363 __ps += long_lag; 00364 00365 // Calculate new x(i) without overflow or division. 00366 // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry 00367 // cannot overflow. 00368 _UIntType __xi; 00369 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry) 00370 { 00371 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry; 00372 _M_carry = 0; 00373 } 00374 else 00375 { 00376 __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps]; 00377 _M_carry = 1; 00378 } 00379 _M_x[_M_p] = __xi; 00380 00381 // Adjust current index to loop around in ring buffer. 00382 if (++_M_p >= long_lag) 00383 _M_p = 0; 00384 00385 return __xi; 00386 } 00387 00388 template<typename _IntType, _IntType __m, int __s, int __r, 00389 typename _CharT, typename _Traits> 00390 std::basic_ostream<_CharT, _Traits>& 00391 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00392 const subtract_with_carry<_IntType, __m, __s, __r>& __x) 00393 { 00394 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00395 typedef typename __ostream_type::ios_base __ios_base; 00396 00397 const typename __ios_base::fmtflags __flags = __os.flags(); 00398 const _CharT __fill = __os.fill(); 00399 const _CharT __space = __os.widen(' '); 00400 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00401 __os.fill(__space); 00402 00403 for (int __i = 0; __i < __r; ++__i) 00404 __os << __x._M_x[__i] << __space; 00405 __os << __x._M_carry; 00406 00407 __os.flags(__flags); 00408 __os.fill(__fill); 00409 return __os; 00410 } 00411 00412 template<typename _IntType, _IntType __m, int __s, int __r, 00413 typename _CharT, typename _Traits> 00414 std::basic_istream<_CharT, _Traits>& 00415 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00416 subtract_with_carry<_IntType, __m, __s, __r>& __x) 00417 { 00418 typedef std::basic_ostream<_CharT, _Traits> __istream_type; 00419 typedef typename __istream_type::ios_base __ios_base; 00420 00421 const typename __ios_base::fmtflags __flags = __is.flags(); 00422 __is.flags(__ios_base::dec | __ios_base::skipws); 00423 00424 for (int __i = 0; __i < __r; ++__i) 00425 __is >> __x._M_x[__i]; 00426 __is >> __x._M_carry; 00427 00428 __is.flags(__flags); 00429 return __is; 00430 } 00431 00432 00433 template<typename _RealType, int __w, int __s, int __r> 00434 void 00435 subtract_with_carry_01<_RealType, __w, __s, __r>:: 00436 _M_initialize_npows() 00437 { 00438 for (int __j = 0; __j < __n; ++__j) 00439 #if _GLIBCXX_USE_C99_MATH_TR1 00440 _M_npows[__j] = std::_GLIBCXX_TR1 ldexp(_RealType(1), -__w + __j * 32); 00441 #else 00442 _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32); 00443 #endif 00444 } 00445 00446 template<typename _RealType, int __w, int __s, int __r> 00447 void 00448 subtract_with_carry_01<_RealType, __w, __s, __r>:: 00449 seed(unsigned long __value) 00450 { 00451 if (__value == 0) 00452 __value = 19780503; 00453 00454 // _GLIBCXX_RESOLVE_LIB_DEFECTS 00455 // 512. Seeding subtract_with_carry_01 from a single unsigned long. 00456 std::_GLIBCXX_TR1 linear_congruential<unsigned long, 40014, 0, 2147483563> 00457 __lcg(__value); 00458 00459 this->seed(__lcg); 00460 } 00461 00462 template<typename _RealType, int __w, int __s, int __r> 00463 template<class _Gen> 00464 void 00465 subtract_with_carry_01<_RealType, __w, __s, __r>:: 00466 seed(_Gen& __gen, false_type) 00467 { 00468 for (int __i = 0; __i < long_lag; ++__i) 00469 { 00470 for (int __j = 0; __j < __n - 1; ++__j) 00471 _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen()); 00472 _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, 00473 __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen()); 00474 } 00475 00476 _M_carry = 1; 00477 for (int __j = 0; __j < __n; ++__j) 00478 if (_M_x[long_lag - 1][__j] != 0) 00479 { 00480 _M_carry = 0; 00481 break; 00482 } 00483 00484 _M_p = 0; 00485 } 00486 00487 template<typename _RealType, int __w, int __s, int __r> 00488 typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type 00489 subtract_with_carry_01<_RealType, __w, __s, __r>:: 00490 operator()() 00491 { 00492 // Derive short lag index from current index. 00493 int __ps = _M_p - short_lag; 00494 if (__ps < 0) 00495 __ps += long_lag; 00496 00497 _UInt32Type __new_carry; 00498 for (int __j = 0; __j < __n - 1; ++__j) 00499 { 00500 if (_M_x[__ps][__j] > _M_x[_M_p][__j] 00501 || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0)) 00502 __new_carry = 0; 00503 else 00504 __new_carry = 1; 00505 00506 _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry; 00507 _M_carry = __new_carry; 00508 } 00509 00510 if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1] 00511 || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0)) 00512 __new_carry = 0; 00513 else 00514 __new_carry = 1; 00515 00516 _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, 00517 __detail::_Shift<_UInt32Type, __w % 32>::__value> 00518 (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry); 00519 _M_carry = __new_carry; 00520 00521 result_type __ret = 0.0; 00522 for (int __j = 0; __j < __n; ++__j) 00523 __ret += _M_x[_M_p][__j] * _M_npows[__j]; 00524 00525 // Adjust current index to loop around in ring buffer. 00526 if (++_M_p >= long_lag) 00527 _M_p = 0; 00528 00529 return __ret; 00530 } 00531 00532 template<typename _RealType, int __w, int __s, int __r, 00533 typename _CharT, typename _Traits> 00534 std::basic_ostream<_CharT, _Traits>& 00535 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00536 const subtract_with_carry_01<_RealType, __w, __s, __r>& __x) 00537 { 00538 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00539 typedef typename __ostream_type::ios_base __ios_base; 00540 00541 const typename __ios_base::fmtflags __flags = __os.flags(); 00542 const _CharT __fill = __os.fill(); 00543 const _CharT __space = __os.widen(' '); 00544 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00545 __os.fill(__space); 00546 00547 for (int __i = 0; __i < __r; ++__i) 00548 for (int __j = 0; __j < __x.__n; ++__j) 00549 __os << __x._M_x[__i][__j] << __space; 00550 __os << __x._M_carry; 00551 00552 __os.flags(__flags); 00553 __os.fill(__fill); 00554 return __os; 00555 } 00556 00557 template<typename _RealType, int __w, int __s, int __r, 00558 typename _CharT, typename _Traits> 00559 std::basic_istream<_CharT, _Traits>& 00560 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00561 subtract_with_carry_01<_RealType, __w, __s, __r>& __x) 00562 { 00563 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00564 typedef typename __istream_type::ios_base __ios_base; 00565 00566 const typename __ios_base::fmtflags __flags = __is.flags(); 00567 __is.flags(__ios_base::dec | __ios_base::skipws); 00568 00569 for (int __i = 0; __i < __r; ++__i) 00570 for (int __j = 0; __j < __x.__n; ++__j) 00571 __is >> __x._M_x[__i][__j]; 00572 __is >> __x._M_carry; 00573 00574 __is.flags(__flags); 00575 return __is; 00576 } 00577 00578 00579 template<class _UniformRandomNumberGenerator, int __p, int __r> 00580 typename discard_block<_UniformRandomNumberGenerator, 00581 __p, __r>::result_type 00582 discard_block<_UniformRandomNumberGenerator, __p, __r>:: 00583 operator()() 00584 { 00585 if (_M_n >= used_block) 00586 { 00587 while (_M_n < block_size) 00588 { 00589 _M_b(); 00590 ++_M_n; 00591 } 00592 _M_n = 0; 00593 } 00594 ++_M_n; 00595 return _M_b(); 00596 } 00597 00598 template<class _UniformRandomNumberGenerator, int __p, int __r, 00599 typename _CharT, typename _Traits> 00600 std::basic_ostream<_CharT, _Traits>& 00601 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00602 const discard_block<_UniformRandomNumberGenerator, 00603 __p, __r>& __x) 00604 { 00605 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00606 typedef typename __ostream_type::ios_base __ios_base; 00607 00608 const typename __ios_base::fmtflags __flags = __os.flags(); 00609 const _CharT __fill = __os.fill(); 00610 const _CharT __space = __os.widen(' '); 00611 __os.flags(__ios_base::dec | __ios_base::fixed 00612 | __ios_base::left); 00613 __os.fill(__space); 00614 00615 __os << __x._M_b << __space << __x._M_n; 00616 00617 __os.flags(__flags); 00618 __os.fill(__fill); 00619 return __os; 00620 } 00621 00622 template<class _UniformRandomNumberGenerator, int __p, int __r, 00623 typename _CharT, typename _Traits> 00624 std::basic_istream<_CharT, _Traits>& 00625 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00626 discard_block<_UniformRandomNumberGenerator, __p, __r>& __x) 00627 { 00628 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00629 typedef typename __istream_type::ios_base __ios_base; 00630 00631 const typename __ios_base::fmtflags __flags = __is.flags(); 00632 __is.flags(__ios_base::dec | __ios_base::skipws); 00633 00634 __is >> __x._M_b >> __x._M_n; 00635 00636 __is.flags(__flags); 00637 return __is; 00638 } 00639 00640 00641 template<class _UniformRandomNumberGenerator1, int __s1, 00642 class _UniformRandomNumberGenerator2, int __s2> 00643 void 00644 xor_combine<_UniformRandomNumberGenerator1, __s1, 00645 _UniformRandomNumberGenerator2, __s2>:: 00646 _M_initialize_max() 00647 { 00648 const int __w = std::numeric_limits<result_type>::digits; 00649 00650 const result_type __m1 = 00651 std::min(result_type(_M_b1.max() - _M_b1.min()), 00652 __detail::_Shift<result_type, __w - __s1>::__value - 1); 00653 00654 const result_type __m2 = 00655 std::min(result_type(_M_b2.max() - _M_b2.min()), 00656 __detail::_Shift<result_type, __w - __s2>::__value - 1); 00657 00658 // NB: In TR1 s1 is not required to be >= s2. 00659 if (__s1 < __s2) 00660 _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1; 00661 else 00662 _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2; 00663 } 00664 00665 template<class _UniformRandomNumberGenerator1, int __s1, 00666 class _UniformRandomNumberGenerator2, int __s2> 00667 typename xor_combine<_UniformRandomNumberGenerator1, __s1, 00668 _UniformRandomNumberGenerator2, __s2>::result_type 00669 xor_combine<_UniformRandomNumberGenerator1, __s1, 00670 _UniformRandomNumberGenerator2, __s2>:: 00671 _M_initialize_max_aux(result_type __a, result_type __b, int __d) 00672 { 00673 const result_type __two2d = result_type(1) << __d; 00674 const result_type __c = __a * __two2d; 00675 00676 if (__a == 0 || __b < __two2d) 00677 return __c + __b; 00678 00679 const result_type __t = std::max(__c, __b); 00680 const result_type __u = std::min(__c, __b); 00681 00682 result_type __ub = __u; 00683 result_type __p; 00684 for (__p = 0; __ub != 1; __ub >>= 1) 00685 ++__p; 00686 00687 const result_type __two2p = result_type(1) << __p; 00688 const result_type __k = __t / __two2p; 00689 00690 if (__k & 1) 00691 return (__k + 1) * __two2p - 1; 00692 00693 if (__c >= __b) 00694 return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p) 00695 / __two2d, 00696 __u % __two2p, __d); 00697 else 00698 return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p) 00699 / __two2d, 00700 __t % __two2p, __d); 00701 } 00702 00703 template<class _UniformRandomNumberGenerator1, int __s1, 00704 class _UniformRandomNumberGenerator2, int __s2, 00705 typename _CharT, typename _Traits> 00706 std::basic_ostream<_CharT, _Traits>& 00707 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00708 const xor_combine<_UniformRandomNumberGenerator1, __s1, 00709 _UniformRandomNumberGenerator2, __s2>& __x) 00710 { 00711 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00712 typedef typename __ostream_type::ios_base __ios_base; 00713 00714 const typename __ios_base::fmtflags __flags = __os.flags(); 00715 const _CharT __fill = __os.fill(); 00716 const _CharT __space = __os.widen(' '); 00717 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00718 __os.fill(__space); 00719 00720 __os << __x.base1() << __space << __x.base2(); 00721 00722 __os.flags(__flags); 00723 __os.fill(__fill); 00724 return __os; 00725 } 00726 00727 template<class _UniformRandomNumberGenerator1, int __s1, 00728 class _UniformRandomNumberGenerator2, int __s2, 00729 typename _CharT, typename _Traits> 00730 std::basic_istream<_CharT, _Traits>& 00731 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00732 xor_combine<_UniformRandomNumberGenerator1, __s1, 00733 _UniformRandomNumberGenerator2, __s2>& __x) 00734 { 00735 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00736 typedef typename __istream_type::ios_base __ios_base; 00737 00738 const typename __ios_base::fmtflags __flags = __is.flags(); 00739 __is.flags(__ios_base::skipws); 00740 00741 __is >> __x._M_b1 >> __x._M_b2; 00742 00743 __is.flags(__flags); 00744 return __is; 00745 } 00746 00747 00748 template<typename _IntType> 00749 template<typename _UniformRandomNumberGenerator> 00750 typename uniform_int<_IntType>::result_type 00751 uniform_int<_IntType>:: 00752 _M_call(_UniformRandomNumberGenerator& __urng, 00753 result_type __min, result_type __max, true_type) 00754 { 00755 // XXX Must be fixed to work well for *arbitrary* __urng.max(), 00756 // __urng.min(), __max, __min. Currently works fine only in the 00757 // most common case __urng.max() - __urng.min() >= __max - __min, 00758 // with __urng.max() > __urng.min() >= 0. 00759 typedef typename __gnu_cxx::__add_unsigned<typename 00760 _UniformRandomNumberGenerator::result_type>::__type __urntype; 00761 typedef typename __gnu_cxx::__add_unsigned<result_type>::__type 00762 __utype; 00763 typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype) 00764 > sizeof(__utype)), 00765 __urntype, __utype>::__type __uctype; 00766 00767 result_type __ret; 00768 00769 const __urntype __urnmin = __urng.min(); 00770 const __urntype __urnmax = __urng.max(); 00771 const __urntype __urnrange = __urnmax - __urnmin; 00772 const __uctype __urange = __max - __min; 00773 const __uctype __udenom = (__urnrange <= __urange 00774 ? 1 : __urnrange / (__urange + 1)); 00775 do 00776 __ret = (__urntype(__urng()) - __urnmin) / __udenom; 00777 while (__ret > __max - __min); 00778 00779 return __ret + __min; 00780 } 00781 00782 template<typename _IntType, typename _CharT, typename _Traits> 00783 std::basic_ostream<_CharT, _Traits>& 00784 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00785 const uniform_int<_IntType>& __x) 00786 { 00787 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00788 typedef typename __ostream_type::ios_base __ios_base; 00789 00790 const typename __ios_base::fmtflags __flags = __os.flags(); 00791 const _CharT __fill = __os.fill(); 00792 const _CharT __space = __os.widen(' '); 00793 __os.flags(__ios_base::scientific | __ios_base::left); 00794 __os.fill(__space); 00795 00796 __os << __x.min() << __space << __x.max(); 00797 00798 __os.flags(__flags); 00799 __os.fill(__fill); 00800 return __os; 00801 } 00802 00803 template<typename _IntType, typename _CharT, typename _Traits> 00804 std::basic_istream<_CharT, _Traits>& 00805 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00806 uniform_int<_IntType>& __x) 00807 { 00808 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00809 typedef typename __istream_type::ios_base __ios_base; 00810 00811 const typename __ios_base::fmtflags __flags = __is.flags(); 00812 __is.flags(__ios_base::dec | __ios_base::skipws); 00813 00814 __is >> __x._M_min >> __x._M_max; 00815 00816 __is.flags(__flags); 00817 return __is; 00818 } 00819 00820 00821 template<typename _CharT, typename _Traits> 00822 std::basic_ostream<_CharT, _Traits>& 00823 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00824 const bernoulli_distribution& __x) 00825 { 00826 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00827 typedef typename __ostream_type::ios_base __ios_base; 00828 00829 const typename __ios_base::fmtflags __flags = __os.flags(); 00830 const _CharT __fill = __os.fill(); 00831 const std::streamsize __precision = __os.precision(); 00832 __os.flags(__ios_base::scientific | __ios_base::left); 00833 __os.fill(__os.widen(' ')); 00834 __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10); 00835 00836 __os << __x.p(); 00837 00838 __os.flags(__flags); 00839 __os.fill(__fill); 00840 __os.precision(__precision); 00841 return __os; 00842 } 00843 00844 00845 template<typename _IntType, typename _RealType> 00846 template<class _UniformRandomNumberGenerator> 00847 typename geometric_distribution<_IntType, _RealType>::result_type 00848 geometric_distribution<_IntType, _RealType>:: 00849 operator()(_UniformRandomNumberGenerator& __urng) 00850 { 00851 // About the epsilon thing see this thread: 00852 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html 00853 const _RealType __naf = 00854 (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 00855 // The largest _RealType convertible to _IntType. 00856 const _RealType __thr = 00857 std::numeric_limits<_IntType>::max() + __naf; 00858 00859 _RealType __cand; 00860 do 00861 __cand = std::ceil(std::log(__urng()) / _M_log_p); 00862 while (__cand >= __thr); 00863 00864 return result_type(__cand + __naf); 00865 } 00866 00867 template<typename _IntType, typename _RealType, 00868 typename _CharT, typename _Traits> 00869 std::basic_ostream<_CharT, _Traits>& 00870 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00871 const geometric_distribution<_IntType, _RealType>& __x) 00872 { 00873 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00874 typedef typename __ostream_type::ios_base __ios_base; 00875 00876 const typename __ios_base::fmtflags __flags = __os.flags(); 00877 const _CharT __fill = __os.fill(); 00878 const std::streamsize __precision = __os.precision(); 00879 __os.flags(__ios_base::scientific | __ios_base::left); 00880 __os.fill(__os.widen(' ')); 00881 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 00882 00883 __os << __x.p(); 00884 00885 __os.flags(__flags); 00886 __os.fill(__fill); 00887 __os.precision(__precision); 00888 return __os; 00889 } 00890 00891 00892 template<typename _IntType, typename _RealType> 00893 void 00894 poisson_distribution<_IntType, _RealType>:: 00895 _M_initialize() 00896 { 00897 #if _GLIBCXX_USE_C99_MATH_TR1 00898 if (_M_mean >= 12) 00899 { 00900 const _RealType __m = std::floor(_M_mean); 00901 _M_lm_thr = std::log(_M_mean); 00902 _M_lfm = std::_GLIBCXX_TR1 lgamma(__m + 1); 00903 _M_sm = std::sqrt(__m); 00904 00905 const _RealType __pi_4 = 0.7853981633974483096156608458198757L; 00906 const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m 00907 / __pi_4)); 00908 _M_d = std::_GLIBCXX_TR1 round(std::max(_RealType(6), 00909 std::min(__m, __dx))); 00910 const _RealType __cx = 2 * __m + _M_d; 00911 _M_scx = std::sqrt(__cx / 2); 00912 _M_1cx = 1 / __cx; 00913 00914 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); 00915 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d; 00916 } 00917 else 00918 #endif 00919 _M_lm_thr = std::exp(-_M_mean); 00920 } 00921 00922 /** 00923 * A rejection algorithm when mean >= 12 and a simple method based 00924 * upon the multiplication of uniform random variates otherwise. 00925 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 00926 * is defined. 00927 * 00928 * Reference: 00929 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 00930 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). 00931 */ 00932 template<typename _IntType, typename _RealType> 00933 template<class _UniformRandomNumberGenerator> 00934 typename poisson_distribution<_IntType, _RealType>::result_type 00935 poisson_distribution<_IntType, _RealType>:: 00936 operator()(_UniformRandomNumberGenerator& __urng) 00937 { 00938 #if _GLIBCXX_USE_C99_MATH_TR1 00939 if (_M_mean >= 12) 00940 { 00941 _RealType __x; 00942 00943 // See comments above... 00944 const _RealType __naf = 00945 (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 00946 const _RealType __thr = 00947 std::numeric_limits<_IntType>::max() + __naf; 00948 00949 const _RealType __m = std::floor(_M_mean); 00950 // sqrt(pi / 2) 00951 const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 00952 const _RealType __c1 = _M_sm * __spi_2; 00953 const _RealType __c2 = _M_c2b + __c1; 00954 const _RealType __c3 = __c2 + 1; 00955 const _RealType __c4 = __c3 + 1; 00956 // e^(1 / 78) 00957 const _RealType __e178 = 1.0129030479320018583185514777512983L; 00958 const _RealType __c5 = __c4 + __e178; 00959 const _RealType __c = _M_cb + __c5; 00960 const _RealType __2cx = 2 * (2 * __m + _M_d); 00961 00962 bool __reject = true; 00963 do 00964 { 00965 const _RealType __u = __c * __urng(); 00966 const _RealType __e = -std::log(__urng()); 00967 00968 _RealType __w = 0.0; 00969 00970 if (__u <= __c1) 00971 { 00972 const _RealType __n = _M_nd(__urng); 00973 const _RealType __y = -std::abs(__n) * _M_sm - 1; 00974 __x = std::floor(__y); 00975 __w = -__n * __n / 2; 00976 if (__x < -__m) 00977 continue; 00978 } 00979 else if (__u <= __c2) 00980 { 00981 const _RealType __n = _M_nd(__urng); 00982 const _RealType __y = 1 + std::abs(__n) * _M_scx; 00983 __x = std::ceil(__y); 00984 __w = __y * (2 - __y) * _M_1cx; 00985 if (__x > _M_d) 00986 continue; 00987 } 00988 else if (__u <= __c3) 00989 // NB: This case not in the book, nor in the Errata, 00990 // but should be ok... 00991 __x = -1; 00992 else if (__u <= __c4) 00993 __x = 0; 00994 else if (__u <= __c5) 00995 __x = 1; 00996 else 00997 { 00998 const _RealType __v = -std::log(__urng()); 00999 const _RealType __y = _M_d + __v * __2cx / _M_d; 01000 __x = std::ceil(__y); 01001 __w = -_M_d * _M_1cx * (1 + __y / 2); 01002 } 01003 01004 __reject = (__w - __e - __x * _M_lm_thr 01005 > _M_lfm - std::_GLIBCXX_TR1 lgamma(__x + __m + 1)); 01006 01007 __reject |= __x + __m >= __thr; 01008 01009 } while (__reject); 01010 01011 return result_type(__x + __m + __naf); 01012 } 01013 else 01014 #endif 01015 { 01016 _IntType __x = 0; 01017 _RealType __prod = 1.0; 01018 01019 do 01020 { 01021 __prod *= __urng(); 01022 __x += 1; 01023 } 01024 while (__prod > _M_lm_thr); 01025 01026 return __x - 1; 01027 } 01028 } 01029 01030 template<typename _IntType, typename _RealType, 01031 typename _CharT, typename _Traits> 01032 std::basic_ostream<_CharT, _Traits>& 01033 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01034 const poisson_distribution<_IntType, _RealType>& __x) 01035 { 01036 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01037 typedef typename __ostream_type::ios_base __ios_base; 01038 01039 const typename __ios_base::fmtflags __flags = __os.flags(); 01040 const _CharT __fill = __os.fill(); 01041 const std::streamsize __precision = __os.precision(); 01042 const _CharT __space = __os.widen(' '); 01043 __os.flags(__ios_base::scientific | __ios_base::left); 01044 __os.fill(__space); 01045 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 01046 01047 __os << __x.mean() << __space << __x._M_nd; 01048 01049 __os.flags(__flags); 01050 __os.fill(__fill); 01051 __os.precision(__precision); 01052 return __os; 01053 } 01054 01055 template<typename _IntType, typename _RealType, 01056 typename _CharT, typename _Traits> 01057 std::basic_istream<_CharT, _Traits>& 01058 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01059 poisson_distribution<_IntType, _RealType>& __x) 01060 { 01061 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01062 typedef typename __istream_type::ios_base __ios_base; 01063 01064 const typename __ios_base::fmtflags __flags = __is.flags(); 01065 __is.flags(__ios_base::skipws); 01066 01067 __is >> __x._M_mean >> __x._M_nd; 01068 __x._M_initialize(); 01069 01070 __is.flags(__flags); 01071 return __is; 01072 } 01073 01074 01075 template<typename _IntType, typename _RealType> 01076 void 01077 binomial_distribution<_IntType, _RealType>:: 01078 _M_initialize() 01079 { 01080 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 01081 01082 _M_easy = true; 01083 01084 #if _GLIBCXX_USE_C99_MATH_TR1 01085 if (_M_t * __p12 >= 8) 01086 { 01087 _M_easy = false; 01088 const _RealType __np = std::floor(_M_t * __p12); 01089 const _RealType __pa = __np / _M_t; 01090 const _RealType __1p = 1 - __pa; 01091 01092 const _RealType __pi_4 = 0.7853981633974483096156608458198757L; 01093 const _RealType __d1x = 01094 std::sqrt(__np * __1p * std::log(32 * __np 01095 / (81 * __pi_4 * __1p))); 01096 _M_d1 = std::_GLIBCXX_TR1 round(std::max(_RealType(1), __d1x)); 01097 const _RealType __d2x = 01098 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p 01099 / (__pi_4 * __pa))); 01100 _M_d2 = std::_GLIBCXX_TR1 round(std::max(_RealType(1), __d2x)); 01101 01102 // sqrt(pi / 2) 01103 const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 01104 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); 01105 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); 01106 _M_c = 2 * _M_d1 / __np; 01107 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; 01108 const _RealType __a12 = _M_a1 + _M_s2 * __spi_2; 01109 const _RealType __s1s = _M_s1 * _M_s1; 01110 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) 01111 * 2 * __s1s / _M_d1 01112 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); 01113 const _RealType __s2s = _M_s2 * _M_s2; 01114 _M_s = (_M_a123 + 2 * __s2s / _M_d2 01115 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); 01116 _M_lf = (std::_GLIBCXX_TR1 lgamma(__np + 1) 01117 + std::_GLIBCXX_TR1 lgamma(_M_t - __np + 1)); 01118 _M_lp1p = std::log(__pa / __1p); 01119 01120 _M_q = -std::log(1 - (__p12 - __pa) / __1p); 01121 } 01122 else 01123 #endif 01124 _M_q = -std::log(1 - __p12); 01125 } 01126 01127 template<typename _IntType, typename _RealType> 01128 template<class _UniformRandomNumberGenerator> 01129 typename binomial_distribution<_IntType, _RealType>::result_type 01130 binomial_distribution<_IntType, _RealType>:: 01131 _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) 01132 { 01133 _IntType __x = 0; 01134 _RealType __sum = 0; 01135 01136 do 01137 { 01138 const _RealType __e = -std::log(__urng()); 01139 __sum += __e / (__t - __x); 01140 __x += 1; 01141 } 01142 while (__sum <= _M_q); 01143 01144 return __x - 1; 01145 } 01146 01147 /** 01148 * A rejection algorithm when t * p >= 8 and a simple waiting time 01149 * method - the second in the referenced book - otherwise. 01150 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 01151 * is defined. 01152 * 01153 * Reference: 01154 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 01155 * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 01156 */ 01157 template<typename _IntType, typename _RealType> 01158 template<class _UniformRandomNumberGenerator> 01159 typename binomial_distribution<_IntType, _RealType>::result_type 01160 binomial_distribution<_IntType, _RealType>:: 01161 operator()(_UniformRandomNumberGenerator& __urng) 01162 { 01163 result_type __ret; 01164 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 01165 01166 #if _GLIBCXX_USE_C99_MATH_TR1 01167 if (!_M_easy) 01168 { 01169 _RealType __x; 01170 01171 // See comments above... 01172 const _RealType __naf = 01173 (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 01174 const _RealType __thr = 01175 std::numeric_limits<_IntType>::max() + __naf; 01176 01177 const _RealType __np = std::floor(_M_t * __p12); 01178 const _RealType __pa = __np / _M_t; 01179 01180 // sqrt(pi / 2) 01181 const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 01182 const _RealType __a1 = _M_a1; 01183 const _RealType __a12 = __a1 + _M_s2 * __spi_2; 01184 const _RealType __a123 = _M_a123; 01185 const _RealType __s1s = _M_s1 * _M_s1; 01186 const _RealType __s2s = _M_s2 * _M_s2; 01187 01188 bool __reject; 01189 do 01190 { 01191 const _RealType __u = _M_s * __urng(); 01192 01193 _RealType __v; 01194 01195 if (__u <= __a1) 01196 { 01197 const _RealType __n = _M_nd(__urng); 01198 const _RealType __y = _M_s1 * std::abs(__n); 01199 __reject = __y >= _M_d1; 01200 if (!__reject) 01201 { 01202 const _RealType __e = -std::log(__urng()); 01203 __x = std::floor(__y); 01204 __v = -__e - __n * __n / 2 + _M_c; 01205 } 01206 } 01207 else if (__u <= __a12) 01208 { 01209 const _RealType __n = _M_nd(__urng); 01210 const _RealType __y = _M_s2 * std::abs(__n); 01211 __reject = __y >= _M_d2; 01212 if (!__reject) 01213 { 01214 const _RealType __e = -std::log(__urng()); 01215 __x = std::floor(-__y); 01216 __v = -__e - __n * __n / 2; 01217 } 01218 } 01219 else if (__u <= __a123) 01220 { 01221 const _RealType __e1 = -std::log(__urng()); 01222 const _RealType __e2 = -std::log(__urng()); 01223 01224 const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1; 01225 __x = std::floor(__y); 01226 __v = (-__e2 + _M_d1 * (1 / (_M_t - __np) 01227 -__y / (2 * __s1s))); 01228 __reject = false; 01229 } 01230 else 01231 { 01232 const _RealType __e1 = -std::log(__urng()); 01233 const _RealType __e2 = -std::log(__urng()); 01234 01235 const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2; 01236 __x = std::floor(-__y); 01237 __v = -__e2 - _M_d2 * __y / (2 * __s2s); 01238 __reject = false; 01239 } 01240 01241 __reject = __reject || __x < -__np || __x > _M_t - __np; 01242 if (!__reject) 01243 { 01244 const _RealType __lfx = 01245 std::_GLIBCXX_TR1 lgamma(__np + __x + 1) 01246 + std::_GLIBCXX_TR1 lgamma(_M_t - (__np + __x) + 1); 01247 __reject = __v > _M_lf - __lfx + __x * _M_lp1p; 01248 } 01249 01250 __reject |= __x + __np >= __thr; 01251 } 01252 while (__reject); 01253 01254 __x += __np + __naf; 01255 01256 const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); 01257 __ret = _IntType(__x) + __z; 01258 } 01259 else 01260 #endif 01261 __ret = _M_waiting(__urng, _M_t); 01262 01263 if (__p12 != _M_p) 01264 __ret = _M_t - __ret; 01265 return __ret; 01266 } 01267 01268 template<typename _IntType, typename _RealType, 01269 typename _CharT, typename _Traits> 01270 std::basic_ostream<_CharT, _Traits>& 01271 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01272 const binomial_distribution<_IntType, _RealType>& __x) 01273 { 01274 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01275 typedef typename __ostream_type::ios_base __ios_base; 01276 01277 const typename __ios_base::fmtflags __flags = __os.flags(); 01278 const _CharT __fill = __os.fill(); 01279 const std::streamsize __precision = __os.precision(); 01280 const _CharT __space = __os.widen(' '); 01281 __os.flags(__ios_base::scientific | __ios_base::left); 01282 __os.fill(__space); 01283 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 01284 01285 __os << __x.t() << __space << __x.p() 01286 << __space << __x._M_nd; 01287 01288 __os.flags(__flags); 01289 __os.fill(__fill); 01290 __os.precision(__precision); 01291 return __os; 01292 } 01293 01294 template<typename _IntType, typename _RealType, 01295 typename _CharT, typename _Traits> 01296 std::basic_istream<_CharT, _Traits>& 01297 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01298 binomial_distribution<_IntType, _RealType>& __x) 01299 { 01300 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01301 typedef typename __istream_type::ios_base __ios_base; 01302 01303 const typename __ios_base::fmtflags __flags = __is.flags(); 01304 __is.flags(__ios_base::dec | __ios_base::skipws); 01305 01306 __is >> __x._M_t >> __x._M_p >> __x._M_nd; 01307 __x._M_initialize(); 01308 01309 __is.flags(__flags); 01310 return __is; 01311 } 01312 01313 01314 template<typename _RealType, typename _CharT, typename _Traits> 01315 std::basic_ostream<_CharT, _Traits>& 01316 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01317 const uniform_real<_RealType>& __x) 01318 { 01319 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01320 typedef typename __ostream_type::ios_base __ios_base; 01321 01322 const typename __ios_base::fmtflags __flags = __os.flags(); 01323 const _CharT __fill = __os.fill(); 01324 const std::streamsize __precision = __os.precision(); 01325 const _CharT __space = __os.widen(' '); 01326 __os.flags(__ios_base::scientific | __ios_base::left); 01327 __os.fill(__space); 01328 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 01329 01330 __os << __x.min() << __space << __x.max(); 01331 01332 __os.flags(__flags); 01333 __os.fill(__fill); 01334 __os.precision(__precision); 01335 return __os; 01336 } 01337 01338 template<typename _RealType, typename _CharT, typename _Traits> 01339 std::basic_istream<_CharT, _Traits>& 01340 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01341 uniform_real<_RealType>& __x) 01342 { 01343 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01344 typedef typename __istream_type::ios_base __ios_base; 01345 01346 const typename __ios_base::fmtflags __flags = __is.flags(); 01347 __is.flags(__ios_base::skipws); 01348 01349 __is >> __x._M_min >> __x._M_max; 01350 01351 __is.flags(__flags); 01352 return __is; 01353 } 01354 01355 01356 template<typename _RealType, typename _CharT, typename _Traits> 01357 std::basic_ostream<_CharT, _Traits>& 01358 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01359 const exponential_distribution<_RealType>& __x) 01360 { 01361 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01362 typedef typename __ostream_type::ios_base __ios_base; 01363 01364 const typename __ios_base::fmtflags __flags = __os.flags(); 01365 const _CharT __fill = __os.fill(); 01366 const std::streamsize __precision = __os.precision(); 01367 __os.flags(__ios_base::scientific | __ios_base::left); 01368 __os.fill(__os.widen(' ')); 01369 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 01370 01371 __os << __x.lambda(); 01372 01373 __os.flags(__flags); 01374 __os.fill(__fill); 01375 __os.precision(__precision); 01376 return __os; 01377 } 01378 01379 01380 /** 01381 * Polar method due to Marsaglia. 01382 * 01383 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 01384 * New York, 1986, Ch. V, Sect. 4.4. 01385 */ 01386 template<typename _RealType> 01387 template<class _UniformRandomNumberGenerator> 01388 typename normal_distribution<_RealType>::result_type 01389 normal_distribution<_RealType>:: 01390 operator()(_UniformRandomNumberGenerator& __urng) 01391 { 01392 result_type __ret; 01393 01394 if (_M_saved_available) 01395 { 01396 _M_saved_available = false; 01397 __ret = _M_saved; 01398 } 01399 else 01400 { 01401 result_type __x, __y, __r2; 01402 do 01403 { 01404 __x = result_type(2.0) * __urng() - 1.0; 01405 __y = result_type(2.0) * __urng() - 1.0; 01406 __r2 = __x * __x + __y * __y; 01407 } 01408 while (__r2 > 1.0 || __r2 == 0.0); 01409 01410 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01411 _M_saved = __x * __mult; 01412 _M_saved_available = true; 01413 __ret = __y * __mult; 01414 } 01415 01416 __ret = __ret * _M_sigma + _M_mean; 01417 return __ret; 01418 } 01419 01420 template<typename _RealType, typename _CharT, typename _Traits> 01421 std::basic_ostream<_CharT, _Traits>& 01422 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01423 const normal_distribution<_RealType>& __x) 01424 { 01425 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01426 typedef typename __ostream_type::ios_base __ios_base; 01427 01428 const typename __ios_base::fmtflags __flags = __os.flags(); 01429 const _CharT __fill = __os.fill(); 01430 const std::streamsize __precision = __os.precision(); 01431 const _CharT __space = __os.widen(' '); 01432 __os.flags(__ios_base::scientific | __ios_base::left); 01433 __os.fill(__space); 01434 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 01435 01436 __os << __x._M_saved_available << __space 01437 << __x.mean() << __space 01438 << __x.sigma(); 01439 if (__x._M_saved_available) 01440 __os << __space << __x._M_saved; 01441 01442 __os.flags(__flags); 01443 __os.fill(__fill); 01444 __os.precision(__precision); 01445 return __os; 01446 } 01447 01448 template<typename _RealType, typename _CharT, typename _Traits> 01449 std::basic_istream<_CharT, _Traits>& 01450 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01451 normal_distribution<_RealType>& __x) 01452 { 01453 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01454 typedef typename __istream_type::ios_base __ios_base; 01455 01456 const typename __ios_base::fmtflags __flags = __is.flags(); 01457 __is.flags(__ios_base::dec | __ios_base::skipws); 01458 01459 __is >> __x._M_saved_available >> __x._M_mean 01460 >> __x._M_sigma; 01461 if (__x._M_saved_available) 01462 __is >> __x._M_saved; 01463 01464 __is.flags(__flags); 01465 return __is; 01466 } 01467 01468 01469 template<typename _RealType> 01470 void 01471 gamma_distribution<_RealType>:: 01472 _M_initialize() 01473 { 01474 if (_M_alpha >= 1) 01475 _M_l_d = std::sqrt(2 * _M_alpha - 1); 01476 else 01477 _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha)) 01478 * (1 - _M_alpha)); 01479 } 01480 01481 /** 01482 * Cheng's rejection algorithm GB for alpha >= 1 and a modification 01483 * of Vaduva's rejection from Weibull algorithm due to Devroye for 01484 * alpha < 1. 01485 * 01486 * References: 01487 * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral 01488 * Shape Parameter." Applied Statistics, 26, 71-75, 1977. 01489 * 01490 * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection 01491 * and Composition Procedures." Math. Operationsforschung and Statistik, 01492 * Series in Statistics, 8, 545-576, 1977. 01493 * 01494 * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 01495 * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!). 01496 */ 01497 template<typename _RealType> 01498 template<class _UniformRandomNumberGenerator> 01499 typename gamma_distribution<_RealType>::result_type 01500 gamma_distribution<_RealType>:: 01501 operator()(_UniformRandomNumberGenerator& __urng) 01502 { 01503 result_type __x; 01504 01505 bool __reject; 01506 if (_M_alpha >= 1) 01507 { 01508 // alpha - log(4) 01509 const result_type __b = _M_alpha 01510 - result_type(1.3862943611198906188344642429163531L); 01511 const result_type __c = _M_alpha + _M_l_d; 01512 const result_type __1l = 1 / _M_l_d; 01513 01514 // 1 + log(9 / 2) 01515 const result_type __k = 2.5040773967762740733732583523868748L; 01516 01517 do 01518 { 01519 const result_type __u = __urng(); 01520 const result_type __v = __urng(); 01521 01522 const result_type __y = __1l * std::log(__v / (1 - __v)); 01523 __x = _M_alpha * std::exp(__y); 01524 01525 const result_type __z = __u * __v * __v; 01526 const result_type __r = __b + __c * __y - __x; 01527 01528 __reject = __r < result_type(4.5) * __z - __k; 01529 if (__reject) 01530 __reject = __r < std::log(__z); 01531 } 01532 while (__reject); 01533 } 01534 else 01535 { 01536 const result_type __c = 1 / _M_alpha; 01537 01538 do 01539 { 01540 const result_type __z = -std::log(__urng()); 01541 const result_type __e = -std::log(__urng()); 01542 01543 __x = std::pow(__z, __c); 01544 01545 __reject = __z + __e < _M_l_d + __x; 01546 } 01547 while (__reject); 01548 } 01549 01550 return __x; 01551 } 01552 01553 template<typename _RealType, typename _CharT, typename _Traits> 01554 std::basic_ostream<_CharT, _Traits>& 01555 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01556 const gamma_distribution<_RealType>& __x) 01557 { 01558 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01559 typedef typename __ostream_type::ios_base __ios_base; 01560 01561 const typename __ios_base::fmtflags __flags = __os.flags(); 01562 const _CharT __fill = __os.fill(); 01563 const std::streamsize __precision = __os.precision(); 01564 __os.flags(__ios_base::scientific | __ios_base::left); 01565 __os.fill(__os.widen(' ')); 01566 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 01567 01568 __os << __x.alpha(); 01569 01570 __os.flags(__flags); 01571 __os.fill(__fill); 01572 __os.precision(__precision); 01573 return __os; 01574 } 01575 01576 _GLIBCXX_END_NAMESPACE_TR1 01577 }