1 // random number generation (out of line) -*- C++ -*- 2 3 // Copyright (C) 2009-2012 Free Software Foundation, Inc. 4 // 5 // This file is part of the GNU ISO C++ Library. This library is free 6 // software; you can redistribute it and/or modify it under the 7 // terms of the GNU General Public License as published by the 8 // Free Software Foundation; either version 3, or (at your option) 9 // any later version. 10 11 // This library is distributed in the hope that it will be useful, 12 // but WITHOUT ANY WARRANTY; without even the implied warranty of 13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 // GNU General Public License for more details. 15 16 // Under Section 7 of GPL version 3, you are granted additional 17 // permissions described in the GCC Runtime Library Exception, version 18 // 3.1, as published by the Free Software Foundation. 19 20 // You should have received a copy of the GNU General Public License and 21 // a copy of the GCC Runtime Library Exception along with this program; 22 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 23 // <http://www.gnu.org/licenses/>. 24 25 /** @file bits/random.tcc 26 * This is an internal header file, included by other library headers. 27 * Do not attempt to use it directly. @headername{random} 28 */ 29 30 #ifndef _RANDOM_TCC 31 #define _RANDOM_TCC 1 32 33 #include <numeric> // std::accumulate and std::partial_sum 34 35 namespace std _GLIBCXX_VISIBILITY(default) 36 { 37 /* 38 * (Further) implementation-space details. 39 */ 40 namespace __detail 41 { 42 _GLIBCXX_BEGIN_NAMESPACE_VERSION 43 44 // General case for x = (ax + c) mod m -- use Schrage's algorithm to 45 // avoid integer overflow. 46 // 47 // Because a and c are compile-time integral constants the compiler 48 // kindly elides any unreachable paths. 49 // 50 // Preconditions: a > 0, m > 0. 51 // 52 // XXX FIXME: as-is, only works correctly for __m % __a < __m / __a. 53 // 54 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool> 55 struct _Mod 56 { 57 static _Tp 58 __calc(_Tp __x) 59 { 60 if (__a == 1) 61 __x %= __m; 62 else 63 { 64 static const _Tp __q = __m / __a; 65 static const _Tp __r = __m % __a; 66 67 _Tp __t1 = __a * (__x % __q); 68 _Tp __t2 = __r * (__x / __q); 69 if (__t1 >= __t2) 70 __x = __t1 - __t2; 71 else 72 __x = __m - __t2 + __t1; 73 } 74 75 if (__c != 0) 76 { 77 const _Tp __d = __m - __x; 78 if (__d > __c) 79 __x += __c; 80 else 81 __x = __c - __d; 82 } 83 return __x; 84 } 85 }; 86 87 // Special case for m == 0 -- use unsigned integer overflow as modulo 88 // operator. 89 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c> 90 struct _Mod<_Tp, __m, __a, __c, true> 91 { 92 static _Tp 93 __calc(_Tp __x) 94 { return __a * __x + __c; } 95 }; 96 97 template<typename _InputIterator, typename _OutputIterator, 98 typename _UnaryOperation> 99 _OutputIterator 100 __transform(_InputIterator __first, _InputIterator __last, 101 _OutputIterator __result, _UnaryOperation __unary_op) 102 { 103 for (; __first != __last; ++__first, ++__result) 104 *__result = __unary_op(*__first); 105 return __result; 106 } 107 108 _GLIBCXX_END_NAMESPACE_VERSION 109 } // namespace __detail 110 111 _GLIBCXX_BEGIN_NAMESPACE_VERSION 112 113 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 114 constexpr _UIntType 115 linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier; 116 117 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 118 constexpr _UIntType 119 linear_congruential_engine<_UIntType, __a, __c, __m>::increment; 120 121 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 122 constexpr _UIntType 123 linear_congruential_engine<_UIntType, __a, __c, __m>::modulus; 124 125 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 126 constexpr _UIntType 127 linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed; 128 129 /** 130 * Seeds the LCR with integral value @p __s, adjusted so that the 131 * ring identity is never a member of the convergence set. 132 */ 133 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 134 void 135 linear_congruential_engine<_UIntType, __a, __c, __m>:: 136 seed(result_type __s) 137 { 138 if ((__detail::__mod<_UIntType, __m>(__c) == 0) 139 && (__detail::__mod<_UIntType, __m>(__s) == 0)) 140 _M_x = 1; 141 else 142 _M_x = __detail::__mod<_UIntType, __m>(__s); 143 } 144 145 /** 146 * Seeds the LCR engine with a value generated by @p __q. 147 */ 148 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 149 template<typename _Sseq> 150 typename std::enable_if<std::is_class<_Sseq>::value>::type 151 linear_congruential_engine<_UIntType, __a, __c, __m>:: 152 seed(_Sseq& __q) 153 { 154 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits 155 : std::__lg(__m); 156 const _UIntType __k = (__k0 + 31) / 32; 157 uint_least32_t __arr[__k + 3]; 158 __q.generate(__arr + 0, __arr + __k + 3); 159 _UIntType __factor = 1u; 160 _UIntType __sum = 0u; 161 for (size_t __j = 0; __j < __k; ++__j) 162 { 163 __sum += __arr[__j + 3] * __factor; 164 __factor *= __detail::_Shift<_UIntType, 32>::__value; 165 } 166 seed(__sum); 167 } 168 169 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 170 typename _CharT, typename _Traits> 171 std::basic_ostream<_CharT, _Traits>& 172 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 173 const linear_congruential_engine<_UIntType, 174 __a, __c, __m>& __lcr) 175 { 176 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 177 typedef typename __ostream_type::ios_base __ios_base; 178 179 const typename __ios_base::fmtflags __flags = __os.flags(); 180 const _CharT __fill = __os.fill(); 181 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 182 __os.fill(__os.widen(' ')); 183 184 __os << __lcr._M_x; 185 186 __os.flags(__flags); 187 __os.fill(__fill); 188 return __os; 189 } 190 191 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 192 typename _CharT, typename _Traits> 193 std::basic_istream<_CharT, _Traits>& 194 operator>>(std::basic_istream<_CharT, _Traits>& __is, 195 linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr) 196 { 197 typedef std::basic_istream<_CharT, _Traits> __istream_type; 198 typedef typename __istream_type::ios_base __ios_base; 199 200 const typename __ios_base::fmtflags __flags = __is.flags(); 201 __is.flags(__ios_base::dec); 202 203 __is >> __lcr._M_x; 204 205 __is.flags(__flags); 206 return __is; 207 } 208 209 210 template<typename _UIntType, 211 size_t __w, size_t __n, size_t __m, size_t __r, 212 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 213 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 214 _UIntType __f> 215 constexpr size_t 216 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 217 __s, __b, __t, __c, __l, __f>::word_size; 218 219 template<typename _UIntType, 220 size_t __w, size_t __n, size_t __m, size_t __r, 221 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 222 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 223 _UIntType __f> 224 constexpr size_t 225 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 226 __s, __b, __t, __c, __l, __f>::state_size; 227 228 template<typename _UIntType, 229 size_t __w, size_t __n, size_t __m, size_t __r, 230 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 231 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 232 _UIntType __f> 233 constexpr size_t 234 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 235 __s, __b, __t, __c, __l, __f>::shift_size; 236 237 template<typename _UIntType, 238 size_t __w, size_t __n, size_t __m, size_t __r, 239 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 240 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 241 _UIntType __f> 242 constexpr size_t 243 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 244 __s, __b, __t, __c, __l, __f>::mask_bits; 245 246 template<typename _UIntType, 247 size_t __w, size_t __n, size_t __m, size_t __r, 248 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 249 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 250 _UIntType __f> 251 constexpr _UIntType 252 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 253 __s, __b, __t, __c, __l, __f>::xor_mask; 254 255 template<typename _UIntType, 256 size_t __w, size_t __n, size_t __m, size_t __r, 257 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 258 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 259 _UIntType __f> 260 constexpr size_t 261 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 262 __s, __b, __t, __c, __l, __f>::tempering_u; 263 264 template<typename _UIntType, 265 size_t __w, size_t __n, size_t __m, size_t __r, 266 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 267 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 268 _UIntType __f> 269 constexpr _UIntType 270 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 271 __s, __b, __t, __c, __l, __f>::tempering_d; 272 273 template<typename _UIntType, 274 size_t __w, size_t __n, size_t __m, size_t __r, 275 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 276 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 277 _UIntType __f> 278 constexpr size_t 279 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 280 __s, __b, __t, __c, __l, __f>::tempering_s; 281 282 template<typename _UIntType, 283 size_t __w, size_t __n, size_t __m, size_t __r, 284 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 285 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 286 _UIntType __f> 287 constexpr _UIntType 288 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 289 __s, __b, __t, __c, __l, __f>::tempering_b; 290 291 template<typename _UIntType, 292 size_t __w, size_t __n, size_t __m, size_t __r, 293 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 294 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 295 _UIntType __f> 296 constexpr size_t 297 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 298 __s, __b, __t, __c, __l, __f>::tempering_t; 299 300 template<typename _UIntType, 301 size_t __w, size_t __n, size_t __m, size_t __r, 302 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 303 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 304 _UIntType __f> 305 constexpr _UIntType 306 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 307 __s, __b, __t, __c, __l, __f>::tempering_c; 308 309 template<typename _UIntType, 310 size_t __w, size_t __n, size_t __m, size_t __r, 311 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 312 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 313 _UIntType __f> 314 constexpr size_t 315 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 316 __s, __b, __t, __c, __l, __f>::tempering_l; 317 318 template<typename _UIntType, 319 size_t __w, size_t __n, size_t __m, size_t __r, 320 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 321 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 322 _UIntType __f> 323 constexpr _UIntType 324 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 325 __s, __b, __t, __c, __l, __f>:: 326 initialization_multiplier; 327 328 template<typename _UIntType, 329 size_t __w, size_t __n, size_t __m, size_t __r, 330 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 331 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 332 _UIntType __f> 333 constexpr _UIntType 334 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 335 __s, __b, __t, __c, __l, __f>::default_seed; 336 337 template<typename _UIntType, 338 size_t __w, size_t __n, size_t __m, size_t __r, 339 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 340 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 341 _UIntType __f> 342 void 343 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 344 __s, __b, __t, __c, __l, __f>:: 345 seed(result_type __sd) 346 { 347 _M_x[0] = __detail::__mod<_UIntType, 348 __detail::_Shift<_UIntType, __w>::__value>(__sd); 349 350 for (size_t __i = 1; __i < state_size; ++__i) 351 { 352 _UIntType __x = _M_x[__i - 1]; 353 __x ^= __x >> (__w - 2); 354 __x *= __f; 355 __x += __detail::__mod<_UIntType, __n>(__i); 356 _M_x[__i] = __detail::__mod<_UIntType, 357 __detail::_Shift<_UIntType, __w>::__value>(__x); 358 } 359 _M_p = state_size; 360 } 361 362 template<typename _UIntType, 363 size_t __w, size_t __n, size_t __m, size_t __r, 364 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 365 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 366 _UIntType __f> 367 template<typename _Sseq> 368 typename std::enable_if<std::is_class<_Sseq>::value>::type 369 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 370 __s, __b, __t, __c, __l, __f>:: 371 seed(_Sseq& __q) 372 { 373 const _UIntType __upper_mask = (~_UIntType()) << __r; 374 const size_t __k = (__w + 31) / 32; 375 uint_least32_t __arr[__n * __k]; 376 __q.generate(__arr + 0, __arr + __n * __k); 377 378 bool __zero = true; 379 for (size_t __i = 0; __i < state_size; ++__i) 380 { 381 _UIntType __factor = 1u; 382 _UIntType __sum = 0u; 383 for (size_t __j = 0; __j < __k; ++__j) 384 { 385 __sum += __arr[__k * __i + __j] * __factor; 386 __factor *= __detail::_Shift<_UIntType, 32>::__value; 387 } 388 _M_x[__i] = __detail::__mod<_UIntType, 389 __detail::_Shift<_UIntType, __w>::__value>(__sum); 390 391 if (__zero) 392 { 393 if (__i == 0) 394 { 395 if ((_M_x[0] & __upper_mask) != 0u) 396 __zero = false; 397 } 398 else if (_M_x[__i] != 0u) 399 __zero = false; 400 } 401 } 402 if (__zero) 403 _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value; 404 _M_p = state_size; 405 } 406 407 template<typename _UIntType, size_t __w, 408 size_t __n, size_t __m, size_t __r, 409 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 410 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 411 _UIntType __f> 412 typename 413 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 414 __s, __b, __t, __c, __l, __f>::result_type 415 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 416 __s, __b, __t, __c, __l, __f>:: 417 operator()() 418 { 419 // Reload the vector - cost is O(n) amortized over n calls. 420 if (_M_p >= state_size) 421 { 422 const _UIntType __upper_mask = (~_UIntType()) << __r; 423 const _UIntType __lower_mask = ~__upper_mask; 424 425 for (size_t __k = 0; __k < (__n - __m); ++__k) 426 { 427 _UIntType __y = ((_M_x[__k] & __upper_mask) 428 | (_M_x[__k + 1] & __lower_mask)); 429 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) 430 ^ ((__y & 0x01) ? __a : 0)); 431 } 432 433 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k) 434 { 435 _UIntType __y = ((_M_x[__k] & __upper_mask) 436 | (_M_x[__k + 1] & __lower_mask)); 437 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) 438 ^ ((__y & 0x01) ? __a : 0)); 439 } 440 441 _UIntType __y = ((_M_x[__n - 1] & __upper_mask) 442 | (_M_x[0] & __lower_mask)); 443 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) 444 ^ ((__y & 0x01) ? __a : 0)); 445 _M_p = 0; 446 } 447 448 // Calculate o(x(i)). 449 result_type __z = _M_x[_M_p++]; 450 __z ^= (__z >> __u) & __d; 451 __z ^= (__z << __s) & __b; 452 __z ^= (__z << __t) & __c; 453 __z ^= (__z >> __l); 454 455 return __z; 456 } 457 458 template<typename _UIntType, size_t __w, 459 size_t __n, size_t __m, size_t __r, 460 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 461 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 462 _UIntType __f, typename _CharT, typename _Traits> 463 std::basic_ostream<_CharT, _Traits>& 464 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 465 const mersenne_twister_engine<_UIntType, __w, __n, __m, 466 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 467 { 468 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 469 typedef typename __ostream_type::ios_base __ios_base; 470 471 const typename __ios_base::fmtflags __flags = __os.flags(); 472 const _CharT __fill = __os.fill(); 473 const _CharT __space = __os.widen(' '); 474 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 475 __os.fill(__space); 476 477 for (size_t __i = 0; __i < __n; ++__i) 478 __os << __x._M_x[__i] << __space; 479 __os << __x._M_p; 480 481 __os.flags(__flags); 482 __os.fill(__fill); 483 return __os; 484 } 485 486 template<typename _UIntType, size_t __w, 487 size_t __n, size_t __m, size_t __r, 488 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 489 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 490 _UIntType __f, typename _CharT, typename _Traits> 491 std::basic_istream<_CharT, _Traits>& 492 operator>>(std::basic_istream<_CharT, _Traits>& __is, 493 mersenne_twister_engine<_UIntType, __w, __n, __m, 494 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 495 { 496 typedef std::basic_istream<_CharT, _Traits> __istream_type; 497 typedef typename __istream_type::ios_base __ios_base; 498 499 const typename __ios_base::fmtflags __flags = __is.flags(); 500 __is.flags(__ios_base::dec | __ios_base::skipws); 501 502 for (size_t __i = 0; __i < __n; ++__i) 503 __is >> __x._M_x[__i]; 504 __is >> __x._M_p; 505 506 __is.flags(__flags); 507 return __is; 508 } 509 510 511 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 512 constexpr size_t 513 subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size; 514 515 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 516 constexpr size_t 517 subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag; 518 519 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 520 constexpr size_t 521 subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag; 522 523 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 524 constexpr _UIntType 525 subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed; 526 527 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 528 void 529 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 530 seed(result_type __value) 531 { 532 std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u> 533 __lcg(__value == 0u ? default_seed : __value); 534 535 const size_t __n = (__w + 31) / 32; 536 537 for (size_t __i = 0; __i < long_lag; ++__i) 538 { 539 _UIntType __sum = 0u; 540 _UIntType __factor = 1u; 541 for (size_t __j = 0; __j < __n; ++__j) 542 { 543 __sum += __detail::__mod<uint_least32_t, 544 __detail::_Shift<uint_least32_t, 32>::__value> 545 (__lcg()) * __factor; 546 __factor *= __detail::_Shift<_UIntType, 32>::__value; 547 } 548 _M_x[__i] = __detail::__mod<_UIntType, 549 __detail::_Shift<_UIntType, __w>::__value>(__sum); 550 } 551 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 552 _M_p = 0; 553 } 554 555 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 556 template<typename _Sseq> 557 typename std::enable_if<std::is_class<_Sseq>::value>::type 558 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 559 seed(_Sseq& __q) 560 { 561 const size_t __k = (__w + 31) / 32; 562 uint_least32_t __arr[__r * __k]; 563 __q.generate(__arr + 0, __arr + __r * __k); 564 565 for (size_t __i = 0; __i < long_lag; ++__i) 566 { 567 _UIntType __sum = 0u; 568 _UIntType __factor = 1u; 569 for (size_t __j = 0; __j < __k; ++__j) 570 { 571 __sum += __arr[__k * __i + __j] * __factor; 572 __factor *= __detail::_Shift<_UIntType, 32>::__value; 573 } 574 _M_x[__i] = __detail::__mod<_UIntType, 575 __detail::_Shift<_UIntType, __w>::__value>(__sum); 576 } 577 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 578 _M_p = 0; 579 } 580 581 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 582 typename subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 583 result_type 584 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 585 operator()() 586 { 587 // Derive short lag index from current index. 588 long __ps = _M_p - short_lag; 589 if (__ps < 0) 590 __ps += long_lag; 591 592 // Calculate new x(i) without overflow or division. 593 // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry 594 // cannot overflow. 595 _UIntType __xi; 596 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry) 597 { 598 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry; 599 _M_carry = 0; 600 } 601 else 602 { 603 __xi = (__detail::_Shift<_UIntType, __w>::__value 604 - _M_x[_M_p] - _M_carry + _M_x[__ps]); 605 _M_carry = 1; 606 } 607 _M_x[_M_p] = __xi; 608 609 // Adjust current index to loop around in ring buffer. 610 if (++_M_p >= long_lag) 611 _M_p = 0; 612 613 return __xi; 614 } 615 616 template<typename _UIntType, size_t __w, size_t __s, size_t __r, 617 typename _CharT, typename _Traits> 618 std::basic_ostream<_CharT, _Traits>& 619 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 620 const subtract_with_carry_engine<_UIntType, 621 __w, __s, __r>& __x) 622 { 623 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 624 typedef typename __ostream_type::ios_base __ios_base; 625 626 const typename __ios_base::fmtflags __flags = __os.flags(); 627 const _CharT __fill = __os.fill(); 628 const _CharT __space = __os.widen(' '); 629 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 630 __os.fill(__space); 631 632 for (size_t __i = 0; __i < __r; ++__i) 633 __os << __x._M_x[__i] << __space; 634 __os << __x._M_carry << __space << __x._M_p; 635 636 __os.flags(__flags); 637 __os.fill(__fill); 638 return __os; 639 } 640 641 template<typename _UIntType, size_t __w, size_t __s, size_t __r, 642 typename _CharT, typename _Traits> 643 std::basic_istream<_CharT, _Traits>& 644 operator>>(std::basic_istream<_CharT, _Traits>& __is, 645 subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x) 646 { 647 typedef std::basic_ostream<_CharT, _Traits> __istream_type; 648 typedef typename __istream_type::ios_base __ios_base; 649 650 const typename __ios_base::fmtflags __flags = __is.flags(); 651 __is.flags(__ios_base::dec | __ios_base::skipws); 652 653 for (size_t __i = 0; __i < __r; ++__i) 654 __is >> __x._M_x[__i]; 655 __is >> __x._M_carry; 656 __is >> __x._M_p; 657 658 __is.flags(__flags); 659 return __is; 660 } 661 662 663 template<typename _RandomNumberEngine, size_t __p, size_t __r> 664 constexpr size_t 665 discard_block_engine<_RandomNumberEngine, __p, __r>::block_size; 666 667 template<typename _RandomNumberEngine, size_t __p, size_t __r> 668 constexpr size_t 669 discard_block_engine<_RandomNumberEngine, __p, __r>::used_block; 670 671 template<typename _RandomNumberEngine, size_t __p, size_t __r> 672 typename discard_block_engine<_RandomNumberEngine, 673 __p, __r>::result_type 674 discard_block_engine<_RandomNumberEngine, __p, __r>:: 675 operator()() 676 { 677 if (_M_n >= used_block) 678 { 679 _M_b.discard(block_size - _M_n); 680 _M_n = 0; 681 } 682 ++_M_n; 683 return _M_b(); 684 } 685 686 template<typename _RandomNumberEngine, size_t __p, size_t __r, 687 typename _CharT, typename _Traits> 688 std::basic_ostream<_CharT, _Traits>& 689 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 690 const discard_block_engine<_RandomNumberEngine, 691 __p, __r>& __x) 692 { 693 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 694 typedef typename __ostream_type::ios_base __ios_base; 695 696 const typename __ios_base::fmtflags __flags = __os.flags(); 697 const _CharT __fill = __os.fill(); 698 const _CharT __space = __os.widen(' '); 699 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 700 __os.fill(__space); 701 702 __os << __x.base() << __space << __x._M_n; 703 704 __os.flags(__flags); 705 __os.fill(__fill); 706 return __os; 707 } 708 709 template<typename _RandomNumberEngine, size_t __p, size_t __r, 710 typename _CharT, typename _Traits> 711 std::basic_istream<_CharT, _Traits>& 712 operator>>(std::basic_istream<_CharT, _Traits>& __is, 713 discard_block_engine<_RandomNumberEngine, __p, __r>& __x) 714 { 715 typedef std::basic_istream<_CharT, _Traits> __istream_type; 716 typedef typename __istream_type::ios_base __ios_base; 717 718 const typename __ios_base::fmtflags __flags = __is.flags(); 719 __is.flags(__ios_base::dec | __ios_base::skipws); 720 721 __is >> __x._M_b >> __x._M_n; 722 723 __is.flags(__flags); 724 return __is; 725 } 726 727 728 template<typename _RandomNumberEngine, size_t __w, typename _UIntType> 729 typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>:: 730 result_type 731 independent_bits_engine<_RandomNumberEngine, __w, _UIntType>:: 732 operator()() 733 { 734 typedef typename _RandomNumberEngine::result_type _Eresult_type; 735 const _Eresult_type __r 736 = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max() 737 ? _M_b.max() - _M_b.min() + 1 : 0); 738 const unsigned __edig = std::numeric_limits<_Eresult_type>::digits; 739 const unsigned __m = __r ? std::__lg(__r) : __edig; 740 741 typedef typename std::common_type<_Eresult_type, result_type>::type 742 __ctype; 743 const unsigned __cdig = std::numeric_limits<__ctype>::digits; 744 745 unsigned __n, __n0; 746 __ctype __s0, __s1, __y0, __y1; 747 748 for (size_t __i = 0; __i < 2; ++__i) 749 { 750 __n = (__w + __m - 1) / __m + __i; 751 __n0 = __n - __w % __n; 752 const unsigned __w0 = __w / __n; // __w0 <= __m 753 754 __s0 = 0; 755 __s1 = 0; 756 if (__w0 < __cdig) 757 { 758 __s0 = __ctype(1) << __w0; 759 __s1 = __s0 << 1; 760 } 761 762 __y0 = 0; 763 __y1 = 0; 764 if (__r) 765 { 766 __y0 = __s0 * (__r / __s0); 767 if (__s1) 768 __y1 = __s1 * (__r / __s1); 769 770 if (__r - __y0 <= __y0 / __n) 771 break; 772 } 773 else 774 break; 775 } 776 777 result_type __sum = 0; 778 for (size_t __k = 0; __k < __n0; ++__k) 779 { 780 __ctype __u; 781 do 782 __u = _M_b() - _M_b.min(); 783 while (__y0 && __u >= __y0); 784 __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u); 785 } 786 for (size_t __k = __n0; __k < __n; ++__k) 787 { 788 __ctype __u; 789 do 790 __u = _M_b() - _M_b.min(); 791 while (__y1 && __u >= __y1); 792 __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u); 793 } 794 return __sum; 795 } 796 797 798 template<typename _RandomNumberEngine, size_t __k> 799 constexpr size_t 800 shuffle_order_engine<_RandomNumberEngine, __k>::table_size; 801 802 template<typename _RandomNumberEngine, size_t __k> 803 typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type 804 shuffle_order_engine<_RandomNumberEngine, __k>:: 805 operator()() 806 { 807 size_t __j = __k * ((_M_y - _M_b.min()) 808 / (_M_b.max() - _M_b.min() + 1.0L)); 809 _M_y = _M_v[__j]; 810 _M_v[__j] = _M_b(); 811 812 return _M_y; 813 } 814 815 template<typename _RandomNumberEngine, size_t __k, 816 typename _CharT, typename _Traits> 817 std::basic_ostream<_CharT, _Traits>& 818 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 819 const shuffle_order_engine<_RandomNumberEngine, __k>& __x) 820 { 821 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 822 typedef typename __ostream_type::ios_base __ios_base; 823 824 const typename __ios_base::fmtflags __flags = __os.flags(); 825 const _CharT __fill = __os.fill(); 826 const _CharT __space = __os.widen(' '); 827 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 828 __os.fill(__space); 829 830 __os << __x.base(); 831 for (size_t __i = 0; __i < __k; ++__i) 832 __os << __space << __x._M_v[__i]; 833 __os << __space << __x._M_y; 834 835 __os.flags(__flags); 836 __os.fill(__fill); 837 return __os; 838 } 839 840 template<typename _RandomNumberEngine, size_t __k, 841 typename _CharT, typename _Traits> 842 std::basic_istream<_CharT, _Traits>& 843 operator>>(std::basic_istream<_CharT, _Traits>& __is, 844 shuffle_order_engine<_RandomNumberEngine, __k>& __x) 845 { 846 typedef std::basic_istream<_CharT, _Traits> __istream_type; 847 typedef typename __istream_type::ios_base __ios_base; 848 849 const typename __ios_base::fmtflags __flags = __is.flags(); 850 __is.flags(__ios_base::dec | __ios_base::skipws); 851 852 __is >> __x._M_b; 853 for (size_t __i = 0; __i < __k; ++__i) 854 __is >> __x._M_v[__i]; 855 __is >> __x._M_y; 856 857 __is.flags(__flags); 858 return __is; 859 } 860 861 862 template<typename _IntType> 863 template<typename _UniformRandomNumberGenerator> 864 typename uniform_int_distribution<_IntType>::result_type 865 uniform_int_distribution<_IntType>:: 866 operator()(_UniformRandomNumberGenerator& __urng, 867 const param_type& __param) 868 { 869 typedef typename _UniformRandomNumberGenerator::result_type 870 _Gresult_type; 871 typedef typename std::make_unsigned<result_type>::type __utype; 872 typedef typename std::common_type<_Gresult_type, __utype>::type 873 __uctype; 874 875 const __uctype __urngmin = __urng.min(); 876 const __uctype __urngmax = __urng.max(); 877 const __uctype __urngrange = __urngmax - __urngmin; 878 const __uctype __urange 879 = __uctype(__param.b()) - __uctype(__param.a()); 880 881 __uctype __ret; 882 883 if (__urngrange > __urange) 884 { 885 // downscaling 886 const __uctype __uerange = __urange + 1; // __urange can be zero 887 const __uctype __scaling = __urngrange / __uerange; 888 const __uctype __past = __uerange * __scaling; 889 do 890 __ret = __uctype(__urng()) - __urngmin; 891 while (__ret >= __past); 892 __ret /= __scaling; 893 } 894 else if (__urngrange < __urange) 895 { 896 // upscaling 897 /* 898 Note that every value in [0, urange] 899 can be written uniquely as 900 901 (urngrange + 1) * high + low 902 903 where 904 905 high in [0, urange / (urngrange + 1)] 906 907 and 908 909 low in [0, urngrange]. 910 */ 911 __uctype __tmp; // wraparound control 912 do 913 { 914 const __uctype __uerngrange = __urngrange + 1; 915 __tmp = (__uerngrange * operator() 916 (__urng, param_type(0, __urange / __uerngrange))); 917 __ret = __tmp + (__uctype(__urng()) - __urngmin); 918 } 919 while (__ret > __urange || __ret < __tmp); 920 } 921 else 922 __ret = __uctype(__urng()) - __urngmin; 923 924 return __ret + __param.a(); 925 } 926 927 template<typename _IntType, typename _CharT, typename _Traits> 928 std::basic_ostream<_CharT, _Traits>& 929 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 930 const uniform_int_distribution<_IntType>& __x) 931 { 932 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 933 typedef typename __ostream_type::ios_base __ios_base; 934 935 const typename __ios_base::fmtflags __flags = __os.flags(); 936 const _CharT __fill = __os.fill(); 937 const _CharT __space = __os.widen(' '); 938 __os.flags(__ios_base::scientific | __ios_base::left); 939 __os.fill(__space); 940 941 __os << __x.a() << __space << __x.b(); 942 943 __os.flags(__flags); 944 __os.fill(__fill); 945 return __os; 946 } 947 948 template<typename _IntType, typename _CharT, typename _Traits> 949 std::basic_istream<_CharT, _Traits>& 950 operator>>(std::basic_istream<_CharT, _Traits>& __is, 951 uniform_int_distribution<_IntType>& __x) 952 { 953 typedef std::basic_istream<_CharT, _Traits> __istream_type; 954 typedef typename __istream_type::ios_base __ios_base; 955 956 const typename __ios_base::fmtflags __flags = __is.flags(); 957 __is.flags(__ios_base::dec | __ios_base::skipws); 958 959 _IntType __a, __b; 960 __is >> __a >> __b; 961 __x.param(typename uniform_int_distribution<_IntType>:: 962 param_type(__a, __b)); 963 964 __is.flags(__flags); 965 return __is; 966 } 967 968 969 template<typename _RealType, typename _CharT, typename _Traits> 970 std::basic_ostream<_CharT, _Traits>& 971 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 972 const uniform_real_distribution<_RealType>& __x) 973 { 974 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 975 typedef typename __ostream_type::ios_base __ios_base; 976 977 const typename __ios_base::fmtflags __flags = __os.flags(); 978 const _CharT __fill = __os.fill(); 979 const std::streamsize __precision = __os.precision(); 980 const _CharT __space = __os.widen(' '); 981 __os.flags(__ios_base::scientific | __ios_base::left); 982 __os.fill(__space); 983 __os.precision(std::numeric_limits<_RealType>::max_digits10); 984 985 __os << __x.a() << __space << __x.b(); 986 987 __os.flags(__flags); 988 __os.fill(__fill); 989 __os.precision(__precision); 990 return __os; 991 } 992 993 template<typename _RealType, typename _CharT, typename _Traits> 994 std::basic_istream<_CharT, _Traits>& 995 operator>>(std::basic_istream<_CharT, _Traits>& __is, 996 uniform_real_distribution<_RealType>& __x) 997 { 998 typedef std::basic_istream<_CharT, _Traits> __istream_type; 999 typedef typename __istream_type::ios_base __ios_base; 1000 1001 const typename __ios_base::fmtflags __flags = __is.flags(); 1002 __is.flags(__ios_base::skipws); 1003 1004 _RealType __a, __b; 1005 __is >> __a >> __b; 1006 __x.param(typename uniform_real_distribution<_RealType>:: 1007 param_type(__a, __b)); 1008 1009 __is.flags(__flags); 1010 return __is; 1011 } 1012 1013 1014 template<typename _CharT, typename _Traits> 1015 std::basic_ostream<_CharT, _Traits>& 1016 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1017 const bernoulli_distribution& __x) 1018 { 1019 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1020 typedef typename __ostream_type::ios_base __ios_base; 1021 1022 const typename __ios_base::fmtflags __flags = __os.flags(); 1023 const _CharT __fill = __os.fill(); 1024 const std::streamsize __precision = __os.precision(); 1025 __os.flags(__ios_base::scientific | __ios_base::left); 1026 __os.fill(__os.widen(' ')); 1027 __os.precision(std::numeric_limits<double>::max_digits10); 1028 1029 __os << __x.p(); 1030 1031 __os.flags(__flags); 1032 __os.fill(__fill); 1033 __os.precision(__precision); 1034 return __os; 1035 } 1036 1037 1038 template<typename _IntType> 1039 template<typename _UniformRandomNumberGenerator> 1040 typename geometric_distribution<_IntType>::result_type 1041 geometric_distribution<_IntType>:: 1042 operator()(_UniformRandomNumberGenerator& __urng, 1043 const param_type& __param) 1044 { 1045 // About the epsilon thing see this thread: 1046 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html 1047 const double __naf = 1048 (1 - std::numeric_limits<double>::epsilon()) / 2; 1049 // The largest _RealType convertible to _IntType. 1050 const double __thr = 1051 std::numeric_limits<_IntType>::max() + __naf; 1052 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1053 __aurng(__urng); 1054 1055 double __cand; 1056 do 1057 __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p); 1058 while (__cand >= __thr); 1059 1060 return result_type(__cand + __naf); 1061 } 1062 1063 template<typename _IntType, 1064 typename _CharT, typename _Traits> 1065 std::basic_ostream<_CharT, _Traits>& 1066 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1067 const geometric_distribution<_IntType>& __x) 1068 { 1069 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1070 typedef typename __ostream_type::ios_base __ios_base; 1071 1072 const typename __ios_base::fmtflags __flags = __os.flags(); 1073 const _CharT __fill = __os.fill(); 1074 const std::streamsize __precision = __os.precision(); 1075 __os.flags(__ios_base::scientific | __ios_base::left); 1076 __os.fill(__os.widen(' ')); 1077 __os.precision(std::numeric_limits<double>::max_digits10); 1078 1079 __os << __x.p(); 1080 1081 __os.flags(__flags); 1082 __os.fill(__fill); 1083 __os.precision(__precision); 1084 return __os; 1085 } 1086 1087 template<typename _IntType, 1088 typename _CharT, typename _Traits> 1089 std::basic_istream<_CharT, _Traits>& 1090 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1091 geometric_distribution<_IntType>& __x) 1092 { 1093 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1094 typedef typename __istream_type::ios_base __ios_base; 1095 1096 const typename __ios_base::fmtflags __flags = __is.flags(); 1097 __is.flags(__ios_base::skipws); 1098 1099 double __p; 1100 __is >> __p; 1101 __x.param(typename geometric_distribution<_IntType>::param_type(__p)); 1102 1103 __is.flags(__flags); 1104 return __is; 1105 } 1106 1107 // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5. 1108 template<typename _IntType> 1109 template<typename _UniformRandomNumberGenerator> 1110 typename negative_binomial_distribution<_IntType>::result_type 1111 negative_binomial_distribution<_IntType>:: 1112 operator()(_UniformRandomNumberGenerator& __urng) 1113 { 1114 const double __y = _M_gd(__urng); 1115 1116 // XXX Is the constructor too slow? 1117 std::poisson_distribution<result_type> __poisson(__y); 1118 return __poisson(__urng); 1119 } 1120 1121 template<typename _IntType> 1122 template<typename _UniformRandomNumberGenerator> 1123 typename negative_binomial_distribution<_IntType>::result_type 1124 negative_binomial_distribution<_IntType>:: 1125 operator()(_UniformRandomNumberGenerator& __urng, 1126 const param_type& __p) 1127 { 1128 typedef typename std::gamma_distribution<result_type>::param_type 1129 param_type; 1130 1131 const double __y = 1132 _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p())); 1133 1134 std::poisson_distribution<result_type> __poisson(__y); 1135 return __poisson(__urng); 1136 } 1137 1138 template<typename _IntType, typename _CharT, typename _Traits> 1139 std::basic_ostream<_CharT, _Traits>& 1140 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1141 const negative_binomial_distribution<_IntType>& __x) 1142 { 1143 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1144 typedef typename __ostream_type::ios_base __ios_base; 1145 1146 const typename __ios_base::fmtflags __flags = __os.flags(); 1147 const _CharT __fill = __os.fill(); 1148 const std::streamsize __precision = __os.precision(); 1149 const _CharT __space = __os.widen(' '); 1150 __os.flags(__ios_base::scientific | __ios_base::left); 1151 __os.fill(__os.widen(' ')); 1152 __os.precision(std::numeric_limits<double>::max_digits10); 1153 1154 __os << __x.k() << __space << __x.p() 1155 << __space << __x._M_gd; 1156 1157 __os.flags(__flags); 1158 __os.fill(__fill); 1159 __os.precision(__precision); 1160 return __os; 1161 } 1162 1163 template<typename _IntType, typename _CharT, typename _Traits> 1164 std::basic_istream<_CharT, _Traits>& 1165 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1166 negative_binomial_distribution<_IntType>& __x) 1167 { 1168 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1169 typedef typename __istream_type::ios_base __ios_base; 1170 1171 const typename __ios_base::fmtflags __flags = __is.flags(); 1172 __is.flags(__ios_base::skipws); 1173 1174 _IntType __k; 1175 double __p; 1176 __is >> __k >> __p >> __x._M_gd; 1177 __x.param(typename negative_binomial_distribution<_IntType>:: 1178 param_type(__k, __p)); 1179 1180 __is.flags(__flags); 1181 return __is; 1182 } 1183 1184 1185 template<typename _IntType> 1186 void 1187 poisson_distribution<_IntType>::param_type:: 1188 _M_initialize() 1189 { 1190 #if _GLIBCXX_USE_C99_MATH_TR1 1191 if (_M_mean >= 12) 1192 { 1193 const double __m = std::floor(_M_mean); 1194 _M_lm_thr = std::log(_M_mean); 1195 _M_lfm = std::lgamma(__m + 1); 1196 _M_sm = std::sqrt(__m); 1197 1198 const double __pi_4 = 0.7853981633974483096156608458198757L; 1199 const double __dx = std::sqrt(2 * __m * std::log(32 * __m 1200 / __pi_4)); 1201 _M_d = std::round(std::max(6.0, std::min(__m, __dx))); 1202 const double __cx = 2 * __m + _M_d; 1203 _M_scx = std::sqrt(__cx / 2); 1204 _M_1cx = 1 / __cx; 1205 1206 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); 1207 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) 1208 / _M_d; 1209 } 1210 else 1211 #endif 1212 _M_lm_thr = std::exp(-_M_mean); 1213 } 1214 1215 /** 1216 * A rejection algorithm when mean >= 12 and a simple method based 1217 * upon the multiplication of uniform random variates otherwise. 1218 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 1219 * is defined. 1220 * 1221 * Reference: 1222 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1223 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). 1224 */ 1225 template<typename _IntType> 1226 template<typename _UniformRandomNumberGenerator> 1227 typename poisson_distribution<_IntType>::result_type 1228 poisson_distribution<_IntType>:: 1229 operator()(_UniformRandomNumberGenerator& __urng, 1230 const param_type& __param) 1231 { 1232 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1233 __aurng(__urng); 1234 #if _GLIBCXX_USE_C99_MATH_TR1 1235 if (__param.mean() >= 12) 1236 { 1237 double __x; 1238 1239 // See comments above... 1240 const double __naf = 1241 (1 - std::numeric_limits<double>::epsilon()) / 2; 1242 const double __thr = 1243 std::numeric_limits<_IntType>::max() + __naf; 1244 1245 const double __m = std::floor(__param.mean()); 1246 // sqrt(pi / 2) 1247 const double __spi_2 = 1.2533141373155002512078826424055226L; 1248 const double __c1 = __param._M_sm * __spi_2; 1249 const double __c2 = __param._M_c2b + __c1; 1250 const double __c3 = __c2 + 1; 1251 const double __c4 = __c3 + 1; 1252 // e^(1 / 78) 1253 const double __e178 = 1.0129030479320018583185514777512983L; 1254 const double __c5 = __c4 + __e178; 1255 const double __c = __param._M_cb + __c5; 1256 const double __2cx = 2 * (2 * __m + __param._M_d); 1257 1258 bool __reject = true; 1259 do 1260 { 1261 const double __u = __c * __aurng(); 1262 const double __e = -std::log(1.0 - __aurng()); 1263 1264 double __w = 0.0; 1265 1266 if (__u <= __c1) 1267 { 1268 const double __n = _M_nd(__urng); 1269 const double __y = -std::abs(__n) * __param._M_sm - 1; 1270 __x = std::floor(__y); 1271 __w = -__n * __n / 2; 1272 if (__x < -__m) 1273 continue; 1274 } 1275 else if (__u <= __c2) 1276 { 1277 const double __n = _M_nd(__urng); 1278 const double __y = 1 + std::abs(__n) * __param._M_scx; 1279 __x = std::ceil(__y); 1280 __w = __y * (2 - __y) * __param._M_1cx; 1281 if (__x > __param._M_d) 1282 continue; 1283 } 1284 else if (__u <= __c3) 1285 // NB: This case not in the book, nor in the Errata, 1286 // but should be ok... 1287 __x = -1; 1288 else if (__u <= __c4) 1289 __x = 0; 1290 else if (__u <= __c5) 1291 __x = 1; 1292 else 1293 { 1294 const double __v = -std::log(1.0 - __aurng()); 1295 const double __y = __param._M_d 1296 + __v * __2cx / __param._M_d; 1297 __x = std::ceil(__y); 1298 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2); 1299 } 1300 1301 __reject = (__w - __e - __x * __param._M_lm_thr 1302 > __param._M_lfm - std::lgamma(__x + __m + 1)); 1303 1304 __reject |= __x + __m >= __thr; 1305 1306 } while (__reject); 1307 1308 return result_type(__x + __m + __naf); 1309 } 1310 else 1311 #endif 1312 { 1313 _IntType __x = 0; 1314 double __prod = 1.0; 1315 1316 do 1317 { 1318 __prod *= __aurng(); 1319 __x += 1; 1320 } 1321 while (__prod > __param._M_lm_thr); 1322 1323 return __x - 1; 1324 } 1325 } 1326 1327 template<typename _IntType, 1328 typename _CharT, typename _Traits> 1329 std::basic_ostream<_CharT, _Traits>& 1330 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1331 const poisson_distribution<_IntType>& __x) 1332 { 1333 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1334 typedef typename __ostream_type::ios_base __ios_base; 1335 1336 const typename __ios_base::fmtflags __flags = __os.flags(); 1337 const _CharT __fill = __os.fill(); 1338 const std::streamsize __precision = __os.precision(); 1339 const _CharT __space = __os.widen(' '); 1340 __os.flags(__ios_base::scientific | __ios_base::left); 1341 __os.fill(__space); 1342 __os.precision(std::numeric_limits<double>::max_digits10); 1343 1344 __os << __x.mean() << __space << __x._M_nd; 1345 1346 __os.flags(__flags); 1347 __os.fill(__fill); 1348 __os.precision(__precision); 1349 return __os; 1350 } 1351 1352 template<typename _IntType, 1353 typename _CharT, typename _Traits> 1354 std::basic_istream<_CharT, _Traits>& 1355 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1356 poisson_distribution<_IntType>& __x) 1357 { 1358 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1359 typedef typename __istream_type::ios_base __ios_base; 1360 1361 const typename __ios_base::fmtflags __flags = __is.flags(); 1362 __is.flags(__ios_base::skipws); 1363 1364 double __mean; 1365 __is >> __mean >> __x._M_nd; 1366 __x.param(typename poisson_distribution<_IntType>::param_type(__mean)); 1367 1368 __is.flags(__flags); 1369 return __is; 1370 } 1371 1372 1373 template<typename _IntType> 1374 void 1375 binomial_distribution<_IntType>::param_type:: 1376 _M_initialize() 1377 { 1378 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 1379 1380 _M_easy = true; 1381 1382 #if _GLIBCXX_USE_C99_MATH_TR1 1383 if (_M_t * __p12 >= 8) 1384 { 1385 _M_easy = false; 1386 const double __np = std::floor(_M_t * __p12); 1387 const double __pa = __np / _M_t; 1388 const double __1p = 1 - __pa; 1389 1390 const double __pi_4 = 0.7853981633974483096156608458198757L; 1391 const double __d1x = 1392 std::sqrt(__np * __1p * std::log(32 * __np 1393 / (81 * __pi_4 * __1p))); 1394 _M_d1 = std::round(std::max(1.0, __d1x)); 1395 const double __d2x = 1396 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p 1397 / (__pi_4 * __pa))); 1398 _M_d2 = std::round(std::max(1.0, __d2x)); 1399 1400 // sqrt(pi / 2) 1401 const double __spi_2 = 1.2533141373155002512078826424055226L; 1402 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); 1403 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); 1404 _M_c = 2 * _M_d1 / __np; 1405 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; 1406 const double __a12 = _M_a1 + _M_s2 * __spi_2; 1407 const double __s1s = _M_s1 * _M_s1; 1408 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) 1409 * 2 * __s1s / _M_d1 1410 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); 1411 const double __s2s = _M_s2 * _M_s2; 1412 _M_s = (_M_a123 + 2 * __s2s / _M_d2 1413 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); 1414 _M_lf = (std::lgamma(__np + 1) 1415 + std::lgamma(_M_t - __np + 1)); 1416 _M_lp1p = std::log(__pa / __1p); 1417 1418 _M_q = -std::log(1 - (__p12 - __pa) / __1p); 1419 } 1420 else 1421 #endif 1422 _M_q = -std::log(1 - __p12); 1423 } 1424 1425 template<typename _IntType> 1426 template<typename _UniformRandomNumberGenerator> 1427 typename binomial_distribution<_IntType>::result_type 1428 binomial_distribution<_IntType>:: 1429 _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) 1430 { 1431 _IntType __x = 0; 1432 double __sum = 0.0; 1433 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1434 __aurng(__urng); 1435 1436 do 1437 { 1438 const double __e = -std::log(1.0 - __aurng()); 1439 __sum += __e / (__t - __x); 1440 __x += 1; 1441 } 1442 while (__sum <= _M_param._M_q); 1443 1444 return __x - 1; 1445 } 1446 1447 /** 1448 * A rejection algorithm when t * p >= 8 and a simple waiting time 1449 * method - the second in the referenced book - otherwise. 1450 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 1451 * is defined. 1452 * 1453 * Reference: 1454 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1455 * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 1456 */ 1457 template<typename _IntType> 1458 template<typename _UniformRandomNumberGenerator> 1459 typename binomial_distribution<_IntType>::result_type 1460 binomial_distribution<_IntType>:: 1461 operator()(_UniformRandomNumberGenerator& __urng, 1462 const param_type& __param) 1463 { 1464 result_type __ret; 1465 const _IntType __t = __param.t(); 1466 const double __p = __param.p(); 1467 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p; 1468 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1469 __aurng(__urng); 1470 1471 #if _GLIBCXX_USE_C99_MATH_TR1 1472 if (!__param._M_easy) 1473 { 1474 double __x; 1475 1476 // See comments above... 1477 const double __naf = 1478 (1 - std::numeric_limits<double>::epsilon()) / 2; 1479 const double __thr = 1480 std::numeric_limits<_IntType>::max() + __naf; 1481 1482 const double __np = std::floor(__t * __p12); 1483 1484 // sqrt(pi / 2) 1485 const double __spi_2 = 1.2533141373155002512078826424055226L; 1486 const double __a1 = __param._M_a1; 1487 const double __a12 = __a1 + __param._M_s2 * __spi_2; 1488 const double __a123 = __param._M_a123; 1489 const double __s1s = __param._M_s1 * __param._M_s1; 1490 const double __s2s = __param._M_s2 * __param._M_s2; 1491 1492 bool __reject; 1493 do 1494 { 1495 const double __u = __param._M_s * __aurng(); 1496 1497 double __v; 1498 1499 if (__u <= __a1) 1500 { 1501 const double __n = _M_nd(__urng); 1502 const double __y = __param._M_s1 * std::abs(__n); 1503 __reject = __y >= __param._M_d1; 1504 if (!__reject) 1505 { 1506 const double __e = -std::log(1.0 - __aurng()); 1507 __x = std::floor(__y); 1508 __v = -__e - __n * __n / 2 + __param._M_c; 1509 } 1510 } 1511 else if (__u <= __a12) 1512 { 1513 const double __n = _M_nd(__urng); 1514 const double __y = __param._M_s2 * std::abs(__n); 1515 __reject = __y >= __param._M_d2; 1516 if (!__reject) 1517 { 1518 const double __e = -std::log(1.0 - __aurng()); 1519 __x = std::floor(-__y); 1520 __v = -__e - __n * __n / 2; 1521 } 1522 } 1523 else if (__u <= __a123) 1524 { 1525 const double __e1 = -std::log(1.0 - __aurng()); 1526 const double __e2 = -std::log(1.0 - __aurng()); 1527 1528 const double __y = __param._M_d1 1529 + 2 * __s1s * __e1 / __param._M_d1; 1530 __x = std::floor(__y); 1531 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np) 1532 -__y / (2 * __s1s))); 1533 __reject = false; 1534 } 1535 else 1536 { 1537 const double __e1 = -std::log(1.0 - __aurng()); 1538 const double __e2 = -std::log(1.0 - __aurng()); 1539 1540 const double __y = __param._M_d2 1541 + 2 * __s2s * __e1 / __param._M_d2; 1542 __x = std::floor(-__y); 1543 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s); 1544 __reject = false; 1545 } 1546 1547 __reject = __reject || __x < -__np || __x > __t - __np; 1548 if (!__reject) 1549 { 1550 const double __lfx = 1551 std::lgamma(__np + __x + 1) 1552 + std::lgamma(__t - (__np + __x) + 1); 1553 __reject = __v > __param._M_lf - __lfx 1554 + __x * __param._M_lp1p; 1555 } 1556 1557 __reject |= __x + __np >= __thr; 1558 } 1559 while (__reject); 1560 1561 __x += __np + __naf; 1562 1563 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x)); 1564 __ret = _IntType(__x) + __z; 1565 } 1566 else 1567 #endif 1568 __ret = _M_waiting(__urng, __t); 1569 1570 if (__p12 != __p) 1571 __ret = __t - __ret; 1572 return __ret; 1573 } 1574 1575 template<typename _IntType, 1576 typename _CharT, typename _Traits> 1577 std::basic_ostream<_CharT, _Traits>& 1578 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1579 const binomial_distribution<_IntType>& __x) 1580 { 1581 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1582 typedef typename __ostream_type::ios_base __ios_base; 1583 1584 const typename __ios_base::fmtflags __flags = __os.flags(); 1585 const _CharT __fill = __os.fill(); 1586 const std::streamsize __precision = __os.precision(); 1587 const _CharT __space = __os.widen(' '); 1588 __os.flags(__ios_base::scientific | __ios_base::left); 1589 __os.fill(__space); 1590 __os.precision(std::numeric_limits<double>::max_digits10); 1591 1592 __os << __x.t() << __space << __x.p() 1593 << __space << __x._M_nd; 1594 1595 __os.flags(__flags); 1596 __os.fill(__fill); 1597 __os.precision(__precision); 1598 return __os; 1599 } 1600 1601 template<typename _IntType, 1602 typename _CharT, typename _Traits> 1603 std::basic_istream<_CharT, _Traits>& 1604 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1605 binomial_distribution<_IntType>& __x) 1606 { 1607 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1608 typedef typename __istream_type::ios_base __ios_base; 1609 1610 const typename __ios_base::fmtflags __flags = __is.flags(); 1611 __is.flags(__ios_base::dec | __ios_base::skipws); 1612 1613 _IntType __t; 1614 double __p; 1615 __is >> __t >> __p >> __x._M_nd; 1616 __x.param(typename binomial_distribution<_IntType>:: 1617 param_type(__t, __p)); 1618 1619 __is.flags(__flags); 1620 return __is; 1621 } 1622 1623 1624 template<typename _RealType, typename _CharT, typename _Traits> 1625 std::basic_ostream<_CharT, _Traits>& 1626 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1627 const exponential_distribution<_RealType>& __x) 1628 { 1629 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1630 typedef typename __ostream_type::ios_base __ios_base; 1631 1632 const typename __ios_base::fmtflags __flags = __os.flags(); 1633 const _CharT __fill = __os.fill(); 1634 const std::streamsize __precision = __os.precision(); 1635 __os.flags(__ios_base::scientific | __ios_base::left); 1636 __os.fill(__os.widen(' ')); 1637 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1638 1639 __os << __x.lambda(); 1640 1641 __os.flags(__flags); 1642 __os.fill(__fill); 1643 __os.precision(__precision); 1644 return __os; 1645 } 1646 1647 template<typename _RealType, typename _CharT, typename _Traits> 1648 std::basic_istream<_CharT, _Traits>& 1649 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1650 exponential_distribution<_RealType>& __x) 1651 { 1652 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1653 typedef typename __istream_type::ios_base __ios_base; 1654 1655 const typename __ios_base::fmtflags __flags = __is.flags(); 1656 __is.flags(__ios_base::dec | __ios_base::skipws); 1657 1658 _RealType __lambda; 1659 __is >> __lambda; 1660 __x.param(typename exponential_distribution<_RealType>:: 1661 param_type(__lambda)); 1662 1663 __is.flags(__flags); 1664 return __is; 1665 } 1666 1667 1668 /** 1669 * Polar method due to Marsaglia. 1670 * 1671 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1672 * New York, 1986, Ch. V, Sect. 4.4. 1673 */ 1674 template<typename _RealType> 1675 template<typename _UniformRandomNumberGenerator> 1676 typename normal_distribution<_RealType>::result_type 1677 normal_distribution<_RealType>:: 1678 operator()(_UniformRandomNumberGenerator& __urng, 1679 const param_type& __param) 1680 { 1681 result_type __ret; 1682 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1683 __aurng(__urng); 1684 1685 if (_M_saved_available) 1686 { 1687 _M_saved_available = false; 1688 __ret = _M_saved; 1689 } 1690 else 1691 { 1692 result_type __x, __y, __r2; 1693 do 1694 { 1695 __x = result_type(2.0) * __aurng() - 1.0; 1696 __y = result_type(2.0) * __aurng() - 1.0; 1697 __r2 = __x * __x + __y * __y; 1698 } 1699 while (__r2 > 1.0 || __r2 == 0.0); 1700 1701 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 1702 _M_saved = __x * __mult; 1703 _M_saved_available = true; 1704 __ret = __y * __mult; 1705 } 1706 1707 __ret = __ret * __param.stddev() + __param.mean(); 1708 return __ret; 1709 } 1710 1711 template<typename _RealType> 1712 bool 1713 operator==(const std::normal_distribution<_RealType>& __d1, 1714 const std::normal_distribution<_RealType>& __d2) 1715 { 1716 if (__d1._M_param == __d2._M_param 1717 && __d1._M_saved_available == __d2._M_saved_available) 1718 { 1719 if (__d1._M_saved_available 1720 && __d1._M_saved == __d2._M_saved) 1721 return true; 1722 else if(!__d1._M_saved_available) 1723 return true; 1724 else 1725 return false; 1726 } 1727 else 1728 return false; 1729 } 1730 1731 template<typename _RealType, typename _CharT, typename _Traits> 1732 std::basic_ostream<_CharT, _Traits>& 1733 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1734 const normal_distribution<_RealType>& __x) 1735 { 1736 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1737 typedef typename __ostream_type::ios_base __ios_base; 1738 1739 const typename __ios_base::fmtflags __flags = __os.flags(); 1740 const _CharT __fill = __os.fill(); 1741 const std::streamsize __precision = __os.precision(); 1742 const _CharT __space = __os.widen(' '); 1743 __os.flags(__ios_base::scientific | __ios_base::left); 1744 __os.fill(__space); 1745 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1746 1747 __os << __x.mean() << __space << __x.stddev() 1748 << __space << __x._M_saved_available; 1749 if (__x._M_saved_available) 1750 __os << __space << __x._M_saved; 1751 1752 __os.flags(__flags); 1753 __os.fill(__fill); 1754 __os.precision(__precision); 1755 return __os; 1756 } 1757 1758 template<typename _RealType, typename _CharT, typename _Traits> 1759 std::basic_istream<_CharT, _Traits>& 1760 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1761 normal_distribution<_RealType>& __x) 1762 { 1763 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1764 typedef typename __istream_type::ios_base __ios_base; 1765 1766 const typename __ios_base::fmtflags __flags = __is.flags(); 1767 __is.flags(__ios_base::dec | __ios_base::skipws); 1768 1769 double __mean, __stddev; 1770 __is >> __mean >> __stddev 1771 >> __x._M_saved_available; 1772 if (__x._M_saved_available) 1773 __is >> __x._M_saved; 1774 __x.param(typename normal_distribution<_RealType>:: 1775 param_type(__mean, __stddev)); 1776 1777 __is.flags(__flags); 1778 return __is; 1779 } 1780 1781 1782 template<typename _RealType, typename _CharT, typename _Traits> 1783 std::basic_ostream<_CharT, _Traits>& 1784 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1785 const lognormal_distribution<_RealType>& __x) 1786 { 1787 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1788 typedef typename __ostream_type::ios_base __ios_base; 1789 1790 const typename __ios_base::fmtflags __flags = __os.flags(); 1791 const _CharT __fill = __os.fill(); 1792 const std::streamsize __precision = __os.precision(); 1793 const _CharT __space = __os.widen(' '); 1794 __os.flags(__ios_base::scientific | __ios_base::left); 1795 __os.fill(__space); 1796 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1797 1798 __os << __x.m() << __space << __x.s() 1799 << __space << __x._M_nd; 1800 1801 __os.flags(__flags); 1802 __os.fill(__fill); 1803 __os.precision(__precision); 1804 return __os; 1805 } 1806 1807 template<typename _RealType, typename _CharT, typename _Traits> 1808 std::basic_istream<_CharT, _Traits>& 1809 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1810 lognormal_distribution<_RealType>& __x) 1811 { 1812 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1813 typedef typename __istream_type::ios_base __ios_base; 1814 1815 const typename __ios_base::fmtflags __flags = __is.flags(); 1816 __is.flags(__ios_base::dec | __ios_base::skipws); 1817 1818 _RealType __m, __s; 1819 __is >> __m >> __s >> __x._M_nd; 1820 __x.param(typename lognormal_distribution<_RealType>:: 1821 param_type(__m, __s)); 1822 1823 __is.flags(__flags); 1824 return __is; 1825 } 1826 1827 1828 template<typename _RealType, typename _CharT, typename _Traits> 1829 std::basic_ostream<_CharT, _Traits>& 1830 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1831 const chi_squared_distribution<_RealType>& __x) 1832 { 1833 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1834 typedef typename __ostream_type::ios_base __ios_base; 1835 1836 const typename __ios_base::fmtflags __flags = __os.flags(); 1837 const _CharT __fill = __os.fill(); 1838 const std::streamsize __precision = __os.precision(); 1839 const _CharT __space = __os.widen(' '); 1840 __os.flags(__ios_base::scientific | __ios_base::left); 1841 __os.fill(__space); 1842 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1843 1844 __os << __x.n() << __space << __x._M_gd; 1845 1846 __os.flags(__flags); 1847 __os.fill(__fill); 1848 __os.precision(__precision); 1849 return __os; 1850 } 1851 1852 template<typename _RealType, typename _CharT, typename _Traits> 1853 std::basic_istream<_CharT, _Traits>& 1854 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1855 chi_squared_distribution<_RealType>& __x) 1856 { 1857 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1858 typedef typename __istream_type::ios_base __ios_base; 1859 1860 const typename __ios_base::fmtflags __flags = __is.flags(); 1861 __is.flags(__ios_base::dec | __ios_base::skipws); 1862 1863 _RealType __n; 1864 __is >> __n >> __x._M_gd; 1865 __x.param(typename chi_squared_distribution<_RealType>:: 1866 param_type(__n)); 1867 1868 __is.flags(__flags); 1869 return __is; 1870 } 1871 1872 1873 template<typename _RealType> 1874 template<typename _UniformRandomNumberGenerator> 1875 typename cauchy_distribution<_RealType>::result_type 1876 cauchy_distribution<_RealType>:: 1877 operator()(_UniformRandomNumberGenerator& __urng, 1878 const param_type& __p) 1879 { 1880 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1881 __aurng(__urng); 1882 _RealType __u; 1883 do 1884 __u = __aurng(); 1885 while (__u == 0.5); 1886 1887 const _RealType __pi = 3.1415926535897932384626433832795029L; 1888 return __p.a() + __p.b() * std::tan(__pi * __u); 1889 } 1890 1891 template<typename _RealType, typename _CharT, typename _Traits> 1892 std::basic_ostream<_CharT, _Traits>& 1893 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1894 const cauchy_distribution<_RealType>& __x) 1895 { 1896 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1897 typedef typename __ostream_type::ios_base __ios_base; 1898 1899 const typename __ios_base::fmtflags __flags = __os.flags(); 1900 const _CharT __fill = __os.fill(); 1901 const std::streamsize __precision = __os.precision(); 1902 const _CharT __space = __os.widen(' '); 1903 __os.flags(__ios_base::scientific | __ios_base::left); 1904 __os.fill(__space); 1905 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1906 1907 __os << __x.a() << __space << __x.b(); 1908 1909 __os.flags(__flags); 1910 __os.fill(__fill); 1911 __os.precision(__precision); 1912 return __os; 1913 } 1914 1915 template<typename _RealType, typename _CharT, typename _Traits> 1916 std::basic_istream<_CharT, _Traits>& 1917 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1918 cauchy_distribution<_RealType>& __x) 1919 { 1920 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1921 typedef typename __istream_type::ios_base __ios_base; 1922 1923 const typename __ios_base::fmtflags __flags = __is.flags(); 1924 __is.flags(__ios_base::dec | __ios_base::skipws); 1925 1926 _RealType __a, __b; 1927 __is >> __a >> __b; 1928 __x.param(typename cauchy_distribution<_RealType>:: 1929 param_type(__a, __b)); 1930 1931 __is.flags(__flags); 1932 return __is; 1933 } 1934 1935 1936 template<typename _RealType, typename _CharT, typename _Traits> 1937 std::basic_ostream<_CharT, _Traits>& 1938 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1939 const fisher_f_distribution<_RealType>& __x) 1940 { 1941 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1942 typedef typename __ostream_type::ios_base __ios_base; 1943 1944 const typename __ios_base::fmtflags __flags = __os.flags(); 1945 const _CharT __fill = __os.fill(); 1946 const std::streamsize __precision = __os.precision(); 1947 const _CharT __space = __os.widen(' '); 1948 __os.flags(__ios_base::scientific | __ios_base::left); 1949 __os.fill(__space); 1950 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1951 1952 __os << __x.m() << __space << __x.n() 1953 << __space << __x._M_gd_x << __space << __x._M_gd_y; 1954 1955 __os.flags(__flags); 1956 __os.fill(__fill); 1957 __os.precision(__precision); 1958 return __os; 1959 } 1960 1961 template<typename _RealType, typename _CharT, typename _Traits> 1962 std::basic_istream<_CharT, _Traits>& 1963 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1964 fisher_f_distribution<_RealType>& __x) 1965 { 1966 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1967 typedef typename __istream_type::ios_base __ios_base; 1968 1969 const typename __ios_base::fmtflags __flags = __is.flags(); 1970 __is.flags(__ios_base::dec | __ios_base::skipws); 1971 1972 _RealType __m, __n; 1973 __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y; 1974 __x.param(typename fisher_f_distribution<_RealType>:: 1975 param_type(__m, __n)); 1976 1977 __is.flags(__flags); 1978 return __is; 1979 } 1980 1981 1982 template<typename _RealType, typename _CharT, typename _Traits> 1983 std::basic_ostream<_CharT, _Traits>& 1984 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1985 const student_t_distribution<_RealType>& __x) 1986 { 1987 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1988 typedef typename __ostream_type::ios_base __ios_base; 1989 1990 const typename __ios_base::fmtflags __flags = __os.flags(); 1991 const _CharT __fill = __os.fill(); 1992 const std::streamsize __precision = __os.precision(); 1993 const _CharT __space = __os.widen(' '); 1994 __os.flags(__ios_base::scientific | __ios_base::left); 1995 __os.fill(__space); 1996 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1997 1998 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd; 1999 2000 __os.flags(__flags); 2001 __os.fill(__fill); 2002 __os.precision(__precision); 2003 return __os; 2004 } 2005 2006 template<typename _RealType, typename _CharT, typename _Traits> 2007 std::basic_istream<_CharT, _Traits>& 2008 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2009 student_t_distribution<_RealType>& __x) 2010 { 2011 typedef std::basic_istream<_CharT, _Traits> __istream_type; 2012 typedef typename __istream_type::ios_base __ios_base; 2013 2014 const typename __ios_base::fmtflags __flags = __is.flags(); 2015 __is.flags(__ios_base::dec | __ios_base::skipws); 2016 2017 _RealType __n; 2018 __is >> __n >> __x._M_nd >> __x._M_gd; 2019 __x.param(typename student_t_distribution<_RealType>::param_type(__n)); 2020 2021 __is.flags(__flags); 2022 return __is; 2023 } 2024 2025 2026 template<typename _RealType> 2027 void 2028 gamma_distribution<_RealType>::param_type:: 2029 _M_initialize() 2030 { 2031 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha; 2032 2033 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0); 2034 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1); 2035 } 2036 2037 /** 2038 * Marsaglia, G. and Tsang, W. W. 2039 * "A Simple Method for Generating Gamma Variables" 2040 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000. 2041 */ 2042 template<typename _RealType> 2043 template<typename _UniformRandomNumberGenerator> 2044 typename gamma_distribution<_RealType>::result_type 2045 gamma_distribution<_RealType>:: 2046 operator()(_UniformRandomNumberGenerator& __urng, 2047 const param_type& __param) 2048 { 2049 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2050 __aurng(__urng); 2051 2052 result_type __u, __v, __n; 2053 const result_type __a1 = (__param._M_malpha 2054 - _RealType(1.0) / _RealType(3.0)); 2055 2056 do 2057 { 2058 do 2059 { 2060 __n = _M_nd(__urng); 2061 __v = result_type(1.0) + __param._M_a2 * __n; 2062 } 2063 while (__v <= 0.0); 2064 2065 __v = __v * __v * __v; 2066 __u = __aurng(); 2067 } 2068 while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 2069 && (std::log(__u) > (0.5 * __n * __n + __a1 2070 * (1.0 - __v + std::log(__v))))); 2071 2072 if (__param.alpha() == __param._M_malpha) 2073 return __a1 * __v * __param.beta(); 2074 else 2075 { 2076 do 2077 __u = __aurng(); 2078 while (__u == 0.0); 2079 2080 return (std::pow(__u, result_type(1.0) / __param.alpha()) 2081 * __a1 * __v * __param.beta()); 2082 } 2083 } 2084 2085 template<typename _RealType, typename _CharT, typename _Traits> 2086 std::basic_ostream<_CharT, _Traits>& 2087 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2088 const gamma_distribution<_RealType>& __x) 2089 { 2090 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2091 typedef typename __ostream_type::ios_base __ios_base; 2092 2093 const typename __ios_base::fmtflags __flags = __os.flags(); 2094 const _CharT __fill = __os.fill(); 2095 const std::streamsize __precision = __os.precision(); 2096 const _CharT __space = __os.widen(' '); 2097 __os.flags(__ios_base::scientific | __ios_base::left); 2098 __os.fill(__space); 2099 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2100 2101 __os << __x.alpha() << __space << __x.beta() 2102 << __space << __x._M_nd; 2103 2104 __os.flags(__flags); 2105 __os.fill(__fill); 2106 __os.precision(__precision); 2107 return __os; 2108 } 2109 2110 template<typename _RealType, typename _CharT, typename _Traits> 2111 std::basic_istream<_CharT, _Traits>& 2112 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2113 gamma_distribution<_RealType>& __x) 2114 { 2115 typedef std::basic_istream<_CharT, _Traits> __istream_type; 2116 typedef typename __istream_type::ios_base __ios_base; 2117 2118 const typename __ios_base::fmtflags __flags = __is.flags(); 2119 __is.flags(__ios_base::dec | __ios_base::skipws); 2120 2121 _RealType __alpha_val, __beta_val; 2122 __is >> __alpha_val >> __beta_val >> __x._M_nd; 2123 __x.param(typename gamma_distribution<_RealType>:: 2124 param_type(__alpha_val, __beta_val)); 2125 2126 __is.flags(__flags); 2127 return __is; 2128 } 2129 2130 2131 template<typename _RealType> 2132 template<typename _UniformRandomNumberGenerator> 2133 typename weibull_distribution<_RealType>::result_type 2134 weibull_distribution<_RealType>:: 2135 operator()(_UniformRandomNumberGenerator& __urng, 2136 const param_type& __p) 2137 { 2138 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2139 __aurng(__urng); 2140 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 2141 result_type(1) / __p.a()); 2142 } 2143 2144 template<typename _RealType, typename _CharT, typename _Traits> 2145 std::basic_ostream<_CharT, _Traits>& 2146 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2147 const weibull_distribution<_RealType>& __x) 2148 { 2149 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2150 typedef typename __ostream_type::ios_base __ios_base; 2151 2152 const typename __ios_base::fmtflags __flags = __os.flags(); 2153 const _CharT __fill = __os.fill(); 2154 const std::streamsize __precision = __os.precision(); 2155 const _CharT __space = __os.widen(' '); 2156 __os.flags(__ios_base::scientific | __ios_base::left); 2157 __os.fill(__space); 2158 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2159 2160 __os << __x.a() << __space << __x.b(); 2161 2162 __os.flags(__flags); 2163 __os.fill(__fill); 2164 __os.precision(__precision); 2165 return __os; 2166 } 2167 2168 template<typename _RealType, typename _CharT, typename _Traits> 2169 std::basic_istream<_CharT, _Traits>& 2170 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2171 weibull_distribution<_RealType>& __x) 2172 { 2173 typedef std::basic_istream<_CharT, _Traits> __istream_type; 2174 typedef typename __istream_type::ios_base __ios_base; 2175 2176 const typename __ios_base::fmtflags __flags = __is.flags(); 2177 __is.flags(__ios_base::dec | __ios_base::skipws); 2178 2179 _RealType __a, __b; 2180 __is >> __a >> __b; 2181 __x.param(typename weibull_distribution<_RealType>:: 2182 param_type(__a, __b)); 2183 2184 __is.flags(__flags); 2185 return __is; 2186 } 2187 2188 2189 template<typename _RealType> 2190 template<typename _UniformRandomNumberGenerator> 2191 typename extreme_value_distribution<_RealType>::result_type 2192 extreme_value_distribution<_RealType>:: 2193 operator()(_UniformRandomNumberGenerator& __urng, 2194 const param_type& __p) 2195 { 2196 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2197 __aurng(__urng); 2198 return __p.a() - __p.b() * std::log(-std::log(result_type(1) 2199 - __aurng())); 2200 } 2201 2202 template<typename _RealType, typename _CharT, typename _Traits> 2203 std::basic_ostream<_CharT, _Traits>& 2204 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2205 const extreme_value_distribution<_RealType>& __x) 2206 { 2207 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2208 typedef typename __ostream_type::ios_base __ios_base; 2209 2210 const typename __ios_base::fmtflags __flags = __os.flags(); 2211 const _CharT __fill = __os.fill(); 2212 const std::streamsize __precision = __os.precision(); 2213 const _CharT __space = __os.widen(' '); 2214 __os.flags(__ios_base::scientific | __ios_base::left); 2215 __os.fill(__space); 2216 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2217 2218 __os << __x.a() << __space << __x.b(); 2219 2220 __os.flags(__flags); 2221 __os.fill(__fill); 2222 __os.precision(__precision); 2223 return __os; 2224 } 2225 2226 template<typename _RealType, typename _CharT, typename _Traits> 2227 std::basic_istream<_CharT, _Traits>& 2228 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2229 extreme_value_distribution<_RealType>& __x) 2230 { 2231 typedef std::basic_istream<_CharT, _Traits> __istream_type; 2232 typedef typename __istream_type::ios_base __ios_base; 2233 2234 const typename __ios_base::fmtflags __flags = __is.flags(); 2235 __is.flags(__ios_base::dec | __ios_base::skipws); 2236 2237 _RealType __a, __b; 2238 __is >> __a >> __b; 2239 __x.param(typename extreme_value_distribution<_RealType>:: 2240 param_type(__a, __b)); 2241 2242 __is.flags(__flags); 2243 return __is; 2244 } 2245 2246 2247 template<typename _IntType> 2248 void 2249 discrete_distribution<_IntType>::param_type:: 2250 _M_initialize() 2251 { 2252 if (_M_prob.size() < 2) 2253 { 2254 _M_prob.clear(); 2255 return; 2256 } 2257 2258 const double __sum = std::accumulate(_M_prob.begin(), 2259 _M_prob.end(), 0.0); 2260 // Now normalize the probabilites. 2261 __detail::__transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), 2262 std::bind2nd(std::divides<double>(), __sum)); 2263 // Accumulate partial sums. 2264 _M_cp.reserve(_M_prob.size()); 2265 std::partial_sum(_M_prob.begin(), _M_prob.end(), 2266 std::back_inserter(_M_cp)); 2267 // Make sure the last cumulative probability is one. 2268 _M_cp[_M_cp.size() - 1] = 1.0; 2269 } 2270 2271 template<typename _IntType> 2272 template<typename _Func> 2273 discrete_distribution<_IntType>::param_type:: 2274 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw) 2275 : _M_prob(), _M_cp() 2276 { 2277 const size_t __n = __nw == 0 ? 1 : __nw; 2278 const double __delta = (__xmax - __xmin) / __n; 2279 2280 _M_prob.reserve(__n); 2281 for (size_t __k = 0; __k < __nw; ++__k) 2282 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta)); 2283 2284 _M_initialize(); 2285 } 2286 2287 template<typename _IntType> 2288 template<typename _UniformRandomNumberGenerator> 2289 typename discrete_distribution<_IntType>::result_type 2290 discrete_distribution<_IntType>:: 2291 operator()(_UniformRandomNumberGenerator& __urng, 2292 const param_type& __param) 2293 { 2294 if (__param._M_cp.empty()) 2295 return result_type(0); 2296 2297 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2298 __aurng(__urng); 2299 2300 const double __p = __aurng(); 2301 auto __pos = std::lower_bound(__param._M_cp.begin(), 2302 __param._M_cp.end(), __p); 2303 2304 return __pos - __param._M_cp.begin(); 2305 } 2306 2307 template<typename _IntType, typename _CharT, typename _Traits> 2308 std::basic_ostream<_CharT, _Traits>& 2309 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2310 const discrete_distribution<_IntType>& __x) 2311 { 2312 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2313 typedef typename __ostream_type::ios_base __ios_base; 2314 2315 const typename __ios_base::fmtflags __flags = __os.flags(); 2316 const _CharT __fill = __os.fill(); 2317 const std::streamsize __precision = __os.precision(); 2318 const _CharT __space = __os.widen(' '); 2319 __os.flags(__ios_base::scientific | __ios_base::left); 2320 __os.fill(__space); 2321 __os.precision(std::numeric_limits<double>::max_digits10); 2322 2323 std::vector<double> __prob = __x.probabilities(); 2324 __os << __prob.size(); 2325 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit) 2326 __os << __space << *__dit; 2327 2328 __os.flags(__flags); 2329 __os.fill(__fill); 2330 __os.precision(__precision); 2331 return __os; 2332 } 2333 2334 template<typename _IntType, typename _CharT, typename _Traits> 2335 std::basic_istream<_CharT, _Traits>& 2336 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2337 discrete_distribution<_IntType>& __x) 2338 { 2339 typedef std::basic_istream<_CharT, _Traits> __istream_type; 2340 typedef typename __istream_type::ios_base __ios_base; 2341 2342 const typename __ios_base::fmtflags __flags = __is.flags(); 2343 __is.flags(__ios_base::dec | __ios_base::skipws); 2344 2345 size_t __n; 2346 __is >> __n; 2347 2348 std::vector<double> __prob_vec; 2349 __prob_vec.reserve(__n); 2350 for (; __n != 0; --__n) 2351 { 2352 double __prob; 2353 __is >> __prob; 2354 __prob_vec.push_back(__prob); 2355 } 2356 2357 __x.param(typename discrete_distribution<_IntType>:: 2358 param_type(__prob_vec.begin(), __prob_vec.end())); 2359 2360 __is.flags(__flags); 2361 return __is; 2362 } 2363 2364 2365 template<typename _RealType> 2366 void 2367 piecewise_constant_distribution<_RealType>::param_type:: 2368 _M_initialize() 2369 { 2370 if (_M_int.size() < 2 2371 || (_M_int.size() == 2 2372 && _M_int[0] == _RealType(0) 2373 && _M_int[1] == _RealType(1))) 2374 { 2375 _M_int.clear(); 2376 _M_den.clear(); 2377 return; 2378 } 2379 2380 const double __sum = std::accumulate(_M_den.begin(), 2381 _M_den.end(), 0.0); 2382 2383 __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(), 2384 std::bind2nd(std::divides<double>(), __sum)); 2385 2386 _M_cp.reserve(_M_den.size()); 2387 std::partial_sum(_M_den.begin(), _M_den.end(), 2388 std::back_inserter(_M_cp)); 2389 2390 // Make sure the last cumulative probability is one. 2391 _M_cp[_M_cp.size() - 1] = 1.0; 2392 2393 for (size_t __k = 0; __k < _M_den.size(); ++__k) 2394 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k]; 2395 } 2396 2397 template<typename _RealType> 2398 template<typename _InputIteratorB, typename _InputIteratorW> 2399 piecewise_constant_distribution<_RealType>::param_type:: 2400 param_type(_InputIteratorB __bbegin, 2401 _InputIteratorB __bend, 2402 _InputIteratorW __wbegin) 2403 : _M_int(), _M_den(), _M_cp() 2404 { 2405 if (__bbegin != __bend) 2406 { 2407 for (;;) 2408 { 2409 _M_int.push_back(*__bbegin); 2410 ++__bbegin; 2411 if (__bbegin == __bend) 2412 break; 2413 2414 _M_den.push_back(*__wbegin); 2415 ++__wbegin; 2416 } 2417 } 2418 2419 _M_initialize(); 2420 } 2421 2422 template<typename _RealType> 2423 template<typename _Func> 2424 piecewise_constant_distribution<_RealType>::param_type:: 2425 param_type(initializer_list<_RealType> __bl, _Func __fw) 2426 : _M_int(), _M_den(), _M_cp() 2427 { 2428 _M_int.reserve(__bl.size()); 2429 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 2430 _M_int.push_back(*__biter); 2431 2432 _M_den.reserve(_M_int.size() - 1); 2433 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 2434 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k]))); 2435 2436 _M_initialize(); 2437 } 2438 2439 template<typename _RealType> 2440 template<typename _Func> 2441 piecewise_constant_distribution<_RealType>::param_type:: 2442 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 2443 : _M_int(), _M_den(), _M_cp() 2444 { 2445 const size_t __n = __nw == 0 ? 1 : __nw; 2446 const _RealType __delta = (__xmax - __xmin) / __n; 2447 2448 _M_int.reserve(__n + 1); 2449 for (size_t __k = 0; __k <= __nw; ++__k) 2450 _M_int.push_back(__xmin + __k * __delta); 2451 2452 _M_den.reserve(__n); 2453 for (size_t __k = 0; __k < __nw; ++__k) 2454 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta)); 2455 2456 _M_initialize(); 2457 } 2458 2459 template<typename _RealType> 2460 template<typename _UniformRandomNumberGenerator> 2461 typename piecewise_constant_distribution<_RealType>::result_type 2462 piecewise_constant_distribution<_RealType>:: 2463 operator()(_UniformRandomNumberGenerator& __urng, 2464 const param_type& __param) 2465 { 2466 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2467 __aurng(__urng); 2468 2469 const double __p = __aurng(); 2470 if (__param._M_cp.empty()) 2471 return __p; 2472 2473 auto __pos = std::lower_bound(__param._M_cp.begin(), 2474 __param._M_cp.end(), __p); 2475 const size_t __i = __pos - __param._M_cp.begin(); 2476 2477 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 2478 2479 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i]; 2480 } 2481 2482 template<typename _RealType, typename _CharT, typename _Traits> 2483 std::basic_ostream<_CharT, _Traits>& 2484 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2485 const piecewise_constant_distribution<_RealType>& __x) 2486 { 2487 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2488 typedef typename __ostream_type::ios_base __ios_base; 2489 2490 const typename __ios_base::fmtflags __flags = __os.flags(); 2491 const _CharT __fill = __os.fill(); 2492 const std::streamsize __precision = __os.precision(); 2493 const _CharT __space = __os.widen(' '); 2494 __os.flags(__ios_base::scientific | __ios_base::left); 2495 __os.fill(__space); 2496 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2497 2498 std::vector<_RealType> __int = __x.intervals(); 2499 __os << __int.size() - 1; 2500 2501 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 2502 __os << __space << *__xit; 2503 2504 std::vector<double> __den = __x.densities(); 2505 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 2506 __os << __space << *__dit; 2507 2508 __os.flags(__flags); 2509 __os.fill(__fill); 2510 __os.precision(__precision); 2511 return __os; 2512 } 2513 2514 template<typename _RealType, typename _CharT, typename _Traits> 2515 std::basic_istream<_CharT, _Traits>& 2516 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2517 piecewise_constant_distribution<_RealType>& __x) 2518 { 2519 typedef std::basic_istream<_CharT, _Traits> __istream_type; 2520 typedef typename __istream_type::ios_base __ios_base; 2521 2522 const typename __ios_base::fmtflags __flags = __is.flags(); 2523 __is.flags(__ios_base::dec | __ios_base::skipws); 2524 2525 size_t __n; 2526 __is >> __n; 2527 2528 std::vector<_RealType> __int_vec; 2529 __int_vec.reserve(__n + 1); 2530 for (size_t __i = 0; __i <= __n; ++__i) 2531 { 2532 _RealType __int; 2533 __is >> __int; 2534 __int_vec.push_back(__int); 2535 } 2536 2537 std::vector<double> __den_vec; 2538 __den_vec.reserve(__n); 2539 for (size_t __i = 0; __i < __n; ++__i) 2540 { 2541 double __den; 2542 __is >> __den; 2543 __den_vec.push_back(__den); 2544 } 2545 2546 __x.param(typename piecewise_constant_distribution<_RealType>:: 2547 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 2548 2549 __is.flags(__flags); 2550 return __is; 2551 } 2552 2553 2554 template<typename _RealType> 2555 void 2556 piecewise_linear_distribution<_RealType>::param_type:: 2557 _M_initialize() 2558 { 2559 if (_M_int.size() < 2 2560 || (_M_int.size() == 2 2561 && _M_int[0] == _RealType(0) 2562 && _M_int[1] == _RealType(1) 2563 && _M_den[0] == _M_den[1])) 2564 { 2565 _M_int.clear(); 2566 _M_den.clear(); 2567 return; 2568 } 2569 2570 double __sum = 0.0; 2571 _M_cp.reserve(_M_int.size() - 1); 2572 _M_m.reserve(_M_int.size() - 1); 2573 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 2574 { 2575 const _RealType __delta = _M_int[__k + 1] - _M_int[__k]; 2576 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta; 2577 _M_cp.push_back(__sum); 2578 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta); 2579 } 2580 2581 // Now normalize the densities... 2582 __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(), 2583 std::bind2nd(std::divides<double>(), __sum)); 2584 // ... and partial sums... 2585 __detail::__transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), 2586 std::bind2nd(std::divides<double>(), __sum)); 2587 // ... and slopes. 2588 __detail::__transform(_M_m.begin(), _M_m.end(), _M_m.begin(), 2589 std::bind2nd(std::divides<double>(), __sum)); 2590 // Make sure the last cumulative probablility is one. 2591 _M_cp[_M_cp.size() - 1] = 1.0; 2592 } 2593 2594 template<typename _RealType> 2595 template<typename _InputIteratorB, typename _InputIteratorW> 2596 piecewise_linear_distribution<_RealType>::param_type:: 2597 param_type(_InputIteratorB __bbegin, 2598 _InputIteratorB __bend, 2599 _InputIteratorW __wbegin) 2600 : _M_int(), _M_den(), _M_cp(), _M_m() 2601 { 2602 for (; __bbegin != __bend; ++__bbegin, ++__wbegin) 2603 { 2604 _M_int.push_back(*__bbegin); 2605 _M_den.push_back(*__wbegin); 2606 } 2607 2608 _M_initialize(); 2609 } 2610 2611 template<typename _RealType> 2612 template<typename _Func> 2613 piecewise_linear_distribution<_RealType>::param_type:: 2614 param_type(initializer_list<_RealType> __bl, _Func __fw) 2615 : _M_int(), _M_den(), _M_cp(), _M_m() 2616 { 2617 _M_int.reserve(__bl.size()); 2618 _M_den.reserve(__bl.size()); 2619 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 2620 { 2621 _M_int.push_back(*__biter); 2622 _M_den.push_back(__fw(*__biter)); 2623 } 2624 2625 _M_initialize(); 2626 } 2627 2628 template<typename _RealType> 2629 template<typename _Func> 2630 piecewise_linear_distribution<_RealType>::param_type:: 2631 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 2632 : _M_int(), _M_den(), _M_cp(), _M_m() 2633 { 2634 const size_t __n = __nw == 0 ? 1 : __nw; 2635 const _RealType __delta = (__xmax - __xmin) / __n; 2636 2637 _M_int.reserve(__n + 1); 2638 _M_den.reserve(__n + 1); 2639 for (size_t __k = 0; __k <= __nw; ++__k) 2640 { 2641 _M_int.push_back(__xmin + __k * __delta); 2642 _M_den.push_back(__fw(_M_int[__k] + __delta)); 2643 } 2644 2645 _M_initialize(); 2646 } 2647 2648 template<typename _RealType> 2649 template<typename _UniformRandomNumberGenerator> 2650 typename piecewise_linear_distribution<_RealType>::result_type 2651 piecewise_linear_distribution<_RealType>:: 2652 operator()(_UniformRandomNumberGenerator& __urng, 2653 const param_type& __param) 2654 { 2655 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2656 __aurng(__urng); 2657 2658 const double __p = __aurng(); 2659 if (__param._M_cp.empty()) 2660 return __p; 2661 2662 auto __pos = std::lower_bound(__param._M_cp.begin(), 2663 __param._M_cp.end(), __p); 2664 const size_t __i = __pos - __param._M_cp.begin(); 2665 2666 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 2667 2668 const double __a = 0.5 * __param._M_m[__i]; 2669 const double __b = __param._M_den[__i]; 2670 const double __cm = __p - __pref; 2671 2672 _RealType __x = __param._M_int[__i]; 2673 if (__a == 0) 2674 __x += __cm / __b; 2675 else 2676 { 2677 const double __d = __b * __b + 4.0 * __a * __cm; 2678 __x += 0.5 * (std::sqrt(__d) - __b) / __a; 2679 } 2680 2681 return __x; 2682 } 2683 2684 template<typename _RealType, typename _CharT, typename _Traits> 2685 std::basic_ostream<_CharT, _Traits>& 2686 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2687 const piecewise_linear_distribution<_RealType>& __x) 2688 { 2689 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2690 typedef typename __ostream_type::ios_base __ios_base; 2691 2692 const typename __ios_base::fmtflags __flags = __os.flags(); 2693 const _CharT __fill = __os.fill(); 2694 const std::streamsize __precision = __os.precision(); 2695 const _CharT __space = __os.widen(' '); 2696 __os.flags(__ios_base::scientific | __ios_base::left); 2697 __os.fill(__space); 2698 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2699 2700 std::vector<_RealType> __int = __x.intervals(); 2701 __os << __int.size() - 1; 2702 2703 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 2704 __os << __space << *__xit; 2705 2706 std::vector<double> __den = __x.densities(); 2707 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 2708 __os << __space << *__dit; 2709 2710 __os.flags(__flags); 2711 __os.fill(__fill); 2712 __os.precision(__precision); 2713 return __os; 2714 } 2715 2716 template<typename _RealType, typename _CharT, typename _Traits> 2717 std::basic_istream<_CharT, _Traits>& 2718 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2719 piecewise_linear_distribution<_RealType>& __x) 2720 { 2721 typedef std::basic_istream<_CharT, _Traits> __istream_type; 2722 typedef typename __istream_type::ios_base __ios_base; 2723 2724 const typename __ios_base::fmtflags __flags = __is.flags(); 2725 __is.flags(__ios_base::dec | __ios_base::skipws); 2726 2727 size_t __n; 2728 __is >> __n; 2729 2730 std::vector<_RealType> __int_vec; 2731 __int_vec.reserve(__n + 1); 2732 for (size_t __i = 0; __i <= __n; ++__i) 2733 { 2734 _RealType __int; 2735 __is >> __int; 2736 __int_vec.push_back(__int); 2737 } 2738 2739 std::vector<double> __den_vec; 2740 __den_vec.reserve(__n + 1); 2741 for (size_t __i = 0; __i <= __n; ++__i) 2742 { 2743 double __den; 2744 __is >> __den; 2745 __den_vec.push_back(__den); 2746 } 2747 2748 __x.param(typename piecewise_linear_distribution<_RealType>:: 2749 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 2750 2751 __is.flags(__flags); 2752 return __is; 2753 } 2754 2755 2756 template<typename _IntType> 2757 seed_seq::seed_seq(std::initializer_list<_IntType> __il) 2758 { 2759 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter) 2760 _M_v.push_back(__detail::__mod<result_type, 2761 __detail::_Shift<result_type, 32>::__value>(*__iter)); 2762 } 2763 2764 template<typename _InputIterator> 2765 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end) 2766 { 2767 for (_InputIterator __iter = __begin; __iter != __end; ++__iter) 2768 _M_v.push_back(__detail::__mod<result_type, 2769 __detail::_Shift<result_type, 32>::__value>(*__iter)); 2770 } 2771 2772 template<typename _RandomAccessIterator> 2773 void 2774 seed_seq::generate(_RandomAccessIterator __begin, 2775 _RandomAccessIterator __end) 2776 { 2777 typedef typename iterator_traits<_RandomAccessIterator>::value_type 2778 _Type; 2779 2780 if (__begin == __end) 2781 return; 2782 2783 std::fill(__begin, __end, _Type(0x8b8b8b8bu)); 2784 2785 const size_t __n = __end - __begin; 2786 const size_t __s = _M_v.size(); 2787 const size_t __t = (__n >= 623) ? 11 2788 : (__n >= 68) ? 7 2789 : (__n >= 39) ? 5 2790 : (__n >= 7) ? 3 2791 : (__n - 1) / 2; 2792 const size_t __p = (__n - __t) / 2; 2793 const size_t __q = __p + __t; 2794 const size_t __m = std::max(size_t(__s + 1), __n); 2795 2796 for (size_t __k = 0; __k < __m; ++__k) 2797 { 2798 _Type __arg = (__begin[__k % __n] 2799 ^ __begin[(__k + __p) % __n] 2800 ^ __begin[(__k - 1) % __n]); 2801 _Type __r1 = __arg ^ (__arg >> 27); 2802 __r1 = __detail::__mod<_Type, 2803 __detail::_Shift<_Type, 32>::__value>(1664525u * __r1); 2804 _Type __r2 = __r1; 2805 if (__k == 0) 2806 __r2 += __s; 2807 else if (__k <= __s) 2808 __r2 += __k % __n + _M_v[__k - 1]; 2809 else 2810 __r2 += __k % __n; 2811 __r2 = __detail::__mod<_Type, 2812 __detail::_Shift<_Type, 32>::__value>(__r2); 2813 __begin[(__k + __p) % __n] += __r1; 2814 __begin[(__k + __q) % __n] += __r2; 2815 __begin[__k % __n] = __r2; 2816 } 2817 2818 for (size_t __k = __m; __k < __m + __n; ++__k) 2819 { 2820 _Type __arg = (__begin[__k % __n] 2821 + __begin[(__k + __p) % __n] 2822 + __begin[(__k - 1) % __n]); 2823 _Type __r3 = __arg ^ (__arg >> 27); 2824 __r3 = __detail::__mod<_Type, 2825 __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3); 2826 _Type __r4 = __r3 - __k % __n; 2827 __r4 = __detail::__mod<_Type, 2828 __detail::_Shift<_Type, 32>::__value>(__r4); 2829 __begin[(__k + __p) % __n] ^= __r3; 2830 __begin[(__k + __q) % __n] ^= __r4; 2831 __begin[__k % __n] = __r4; 2832 } 2833 } 2834 2835 template<typename _RealType, size_t __bits, 2836 typename _UniformRandomNumberGenerator> 2837 _RealType 2838 generate_canonical(_UniformRandomNumberGenerator& __urng) 2839 { 2840 const size_t __b 2841 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits), 2842 __bits); 2843 const long double __r = static_cast<long double>(__urng.max()) 2844 - static_cast<long double>(__urng.min()) + 1.0L; 2845 const size_t __log2r = std::log(__r) / std::log(2.0L); 2846 size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r); 2847 _RealType __sum = _RealType(0); 2848 _RealType __tmp = _RealType(1); 2849 for (; __k != 0; --__k) 2850 { 2851 __sum += _RealType(__urng() - __urng.min()) * __tmp; 2852 __tmp *= __r; 2853 } 2854 return __sum / __tmp; 2855 } 2856 2857 _GLIBCXX_END_NAMESPACE_VERSION 2858 } // namespace 2859 2860 #endif 2861