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 namespace std
00031 {
00032 _GLIBCXX_BEGIN_NAMESPACE_TR1
00033
00034
00035
00036
00037 namespace __detail
00038 {
00039
00040
00041
00042
00043
00044
00045
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
00081
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 }
00090
00091
00092
00093
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
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
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
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
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
00361 int __ps = _M_p - short_lag;
00362 if (__ps < 0)
00363 __ps += long_lag;
00364
00365
00366
00367
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
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
00455
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
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
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
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
00756
00757
00758
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
00852
00853 const _RealType __naf =
00854 (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
00855
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
00924
00925
00926
00927
00928
00929
00930
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
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
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
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
00990
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
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
01149
01150
01151
01152
01153
01154
01155
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
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
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
01382
01383
01384
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
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
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
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
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 }