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