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