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