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