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