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