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