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