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