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