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