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