1 // random number generation (out of line) -*- C++ -*- 2 3 // Copyright (C) 2009-2022 Free Software Foundation, Inc. 4 // 5 // This file is part of the GNU ISO C++ Library. This library is free 6 // software; you can redistribute it and/or modify it under the 7 // terms of the GNU General Public License as published by the 8 // Free Software Foundation; either version 3, or (at your option) 9 // any later version. 10 11 // This library is distributed in the hope that it will be useful, 12 // but WITHOUT ANY WARRANTY; without even the implied warranty of 13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 // GNU General Public License for more details. 15 16 // Under Section 7 of GPL version 3, you are granted additional 17 // permissions described in the GCC Runtime Library Exception, version 18 // 3.1, as published by the Free Software Foundation. 19 20 // You should have received a copy of the GNU General Public License and 21 // a copy of the GCC Runtime Library Exception along with this program; 22 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 23 // <http://www.gnu.org/licenses/>. 24 25 /** @file bits/random.tcc 26 * This is an internal header file, included by other library headers. 27 * Do not attempt to use it directly. @headername{random} 28 */ 29 30 #ifndef _RANDOM_TCC 31 #define _RANDOM_TCC 1 32 33 #include <numeric> // std::accumulate and std::partial_sum 34 35 namespace std _GLIBCXX_VISIBILITY(default) 36 { 37 _GLIBCXX_BEGIN_NAMESPACE_VERSION 38 39 /// @cond undocumented 40 // (Further) implementation-space details. 41 namespace __detail 42 { 43 // General case for x = (ax + c) mod m -- use Schrage's algorithm 44 // to avoid integer overflow. 45 // 46 // Preconditions: a > 0, m > 0. 47 // 48 // Note: only works correctly for __m % __a < __m / __a. 49 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c> 50 _Tp 51 _Mod<_Tp, __m, __a, __c, false, true>:: __calc(_Tp __x)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 __normalize(_InputIterator __first,_InputIterator __last,_OutputIterator __result,const _Tp & __factor)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>:: seed(result_type __s)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>:: seed(_Sseq & __q)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const linear_congruential_engine<_UIntType,__a,__c,__m> & __lcr)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,linear_congruential_engine<_UIntType,__a,__c,__m> & __lcr)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>:: seed(result_type __sd)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>:: seed(_Sseq & __q)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>:: _M_gen_rand(void)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>:: discard(unsigned long long __z)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>:: operator ()()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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const mersenne_twister_engine<_UIntType,__w,__n,__m,__r,__a,__u,__d,__s,__b,__t,__c,__l,__f> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,mersenne_twister_engine<_UIntType,__w,__n,__m,__r,__a,__u,__d,__s,__b,__t,__c,__l,__f> & __x)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>:: seed(result_type __value)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>:: seed(_Sseq & __q)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>:: operator ()()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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const subtract_with_carry_engine<_UIntType,__w,__s,__r> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,subtract_with_carry_engine<_UIntType,__w,__s,__r> & __x)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>:: operator ()()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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const discard_block_engine<_RandomNumberEngine,__p,__r> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,discard_block_engine<_RandomNumberEngine,__p,__r> & __x)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>:: operator ()()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 __representable_as_double(_Tp __x)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 __p1_representable_as_double(_Tp __x)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>:: operator ()()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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const shuffle_order_engine<_RandomNumberEngine,__k> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,shuffle_order_engine<_RandomNumberEngine,__k> & __x)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const uniform_int_distribution<_IntType> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,uniform_int_distribution<_IntType> & __x)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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const uniform_real_distribution<_RealType> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,uniform_real_distribution<_RealType> & __x)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:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const bernoulli_distribution & __x)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>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const geometric_distribution<_IntType> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,geometric_distribution<_IntType> & __x)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>:: operator ()(_UniformRandomNumberGenerator & __urng)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>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng)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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const negative_binomial_distribution<_IntType> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,negative_binomial_distribution<_IntType> & __x)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:: _M_initialize()1271 _M_initialize() 1272 { 1273 #if _GLIBCXX_USE_C99_MATH_TR1 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_TR1 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>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)1312 operator()(_UniformRandomNumberGenerator& __urng, 1313 const param_type& __param) 1314 { 1315 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1316 __aurng(__urng); 1317 #if _GLIBCXX_USE_C99_MATH_TR1 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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const poisson_distribution<_IntType> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,poisson_distribution<_IntType> & __x)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:: _M_initialize()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_TR1 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>:: _M_waiting(_UniformRandomNumberGenerator & __urng,_IntType __t,double __q)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_TR1 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>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)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_TR1 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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const binomial_distribution<_IntType> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,binomial_distribution<_IntType> & __x)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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)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>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const exponential_distribution<_RealType> & __x)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>& operator >>(std::basic_istream<_CharT,_Traits> & __is,exponential_distribution<_RealType> & __x)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>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)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>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)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 operator ==(const std::normal_distribution<_RealType> & __d1,const std::normal_distribution<_RealType> & __d2)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 { 1914 if (__d1._M_saved_available 1915 && __d1._M_saved == __d2._M_saved) 1916 return true; 1917 else if(!__d1._M_saved_available) 1918 return true; 1919 else 1920 return false; 1921 } 1922 else 1923 return false; 1924 } 1925 1926 template<typename _RealType, typename _CharT, typename _Traits> 1927 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const normal_distribution<_RealType> & __x)1928 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1929 const normal_distribution<_RealType>& __x) 1930 { 1931 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 1932 1933 const typename __ios_base::fmtflags __flags = __os.flags(); 1934 const _CharT __fill = __os.fill(); 1935 const std::streamsize __precision = __os.precision(); 1936 const _CharT __space = __os.widen(' '); 1937 __os.flags(__ios_base::scientific | __ios_base::left); 1938 __os.fill(__space); 1939 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1940 1941 __os << __x.mean() << __space << __x.stddev() 1942 << __space << __x._M_saved_available; 1943 if (__x._M_saved_available) 1944 __os << __space << __x._M_saved; 1945 1946 __os.flags(__flags); 1947 __os.fill(__fill); 1948 __os.precision(__precision); 1949 return __os; 1950 } 1951 1952 template<typename _RealType, typename _CharT, typename _Traits> 1953 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,normal_distribution<_RealType> & __x)1954 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1955 normal_distribution<_RealType>& __x) 1956 { 1957 using param_type = typename normal_distribution<_RealType>::param_type; 1958 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 1959 1960 const typename __ios_base::fmtflags __flags = __is.flags(); 1961 __is.flags(__ios_base::dec | __ios_base::skipws); 1962 1963 double __mean, __stddev; 1964 bool __saved_avail; 1965 if (__is >> __mean >> __stddev >> __saved_avail) 1966 { 1967 if (!__saved_avail || (__is >> __x._M_saved)) 1968 { 1969 __x._M_saved_available = __saved_avail; 1970 __x.param(param_type(__mean, __stddev)); 1971 } 1972 } 1973 1974 __is.flags(__flags); 1975 return __is; 1976 } 1977 1978 1979 template<typename _RealType> 1980 template<typename _ForwardIterator, 1981 typename _UniformRandomNumberGenerator> 1982 void 1983 lognormal_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)1984 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 1985 _UniformRandomNumberGenerator& __urng, 1986 const param_type& __p) 1987 { 1988 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 1989 while (__f != __t) 1990 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m()); 1991 } 1992 1993 template<typename _RealType, typename _CharT, typename _Traits> 1994 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const lognormal_distribution<_RealType> & __x)1995 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1996 const lognormal_distribution<_RealType>& __x) 1997 { 1998 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 1999 2000 const typename __ios_base::fmtflags __flags = __os.flags(); 2001 const _CharT __fill = __os.fill(); 2002 const std::streamsize __precision = __os.precision(); 2003 const _CharT __space = __os.widen(' '); 2004 __os.flags(__ios_base::scientific | __ios_base::left); 2005 __os.fill(__space); 2006 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2007 2008 __os << __x.m() << __space << __x.s() 2009 << __space << __x._M_nd; 2010 2011 __os.flags(__flags); 2012 __os.fill(__fill); 2013 __os.precision(__precision); 2014 return __os; 2015 } 2016 2017 template<typename _RealType, typename _CharT, typename _Traits> 2018 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,lognormal_distribution<_RealType> & __x)2019 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2020 lognormal_distribution<_RealType>& __x) 2021 { 2022 using param_type 2023 = typename lognormal_distribution<_RealType>::param_type; 2024 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2025 2026 const typename __ios_base::fmtflags __flags = __is.flags(); 2027 __is.flags(__ios_base::dec | __ios_base::skipws); 2028 2029 _RealType __m, __s; 2030 if (__is >> __m >> __s >> __x._M_nd) 2031 __x.param(param_type(__m, __s)); 2032 2033 __is.flags(__flags); 2034 return __is; 2035 } 2036 2037 template<typename _RealType> 2038 template<typename _ForwardIterator, 2039 typename _UniformRandomNumberGenerator> 2040 void 2041 std::chi_squared_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng)2042 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2043 _UniformRandomNumberGenerator& __urng) 2044 { 2045 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2046 while (__f != __t) 2047 *__f++ = 2 * _M_gd(__urng); 2048 } 2049 2050 template<typename _RealType> 2051 template<typename _ForwardIterator, 2052 typename _UniformRandomNumberGenerator> 2053 void 2054 std::chi_squared_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const typename std::gamma_distribution<result_type>::param_type & __p)2055 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2056 _UniformRandomNumberGenerator& __urng, 2057 const typename 2058 std::gamma_distribution<result_type>::param_type& __p) 2059 { 2060 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2061 while (__f != __t) 2062 *__f++ = 2 * _M_gd(__urng, __p); 2063 } 2064 2065 template<typename _RealType, typename _CharT, typename _Traits> 2066 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const chi_squared_distribution<_RealType> & __x)2067 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2068 const chi_squared_distribution<_RealType>& __x) 2069 { 2070 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2071 2072 const typename __ios_base::fmtflags __flags = __os.flags(); 2073 const _CharT __fill = __os.fill(); 2074 const std::streamsize __precision = __os.precision(); 2075 const _CharT __space = __os.widen(' '); 2076 __os.flags(__ios_base::scientific | __ios_base::left); 2077 __os.fill(__space); 2078 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2079 2080 __os << __x.n() << __space << __x._M_gd; 2081 2082 __os.flags(__flags); 2083 __os.fill(__fill); 2084 __os.precision(__precision); 2085 return __os; 2086 } 2087 2088 template<typename _RealType, typename _CharT, typename _Traits> 2089 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,chi_squared_distribution<_RealType> & __x)2090 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2091 chi_squared_distribution<_RealType>& __x) 2092 { 2093 using param_type 2094 = typename chi_squared_distribution<_RealType>::param_type; 2095 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2096 2097 const typename __ios_base::fmtflags __flags = __is.flags(); 2098 __is.flags(__ios_base::dec | __ios_base::skipws); 2099 2100 _RealType __n; 2101 if (__is >> __n >> __x._M_gd) 2102 __x.param(param_type(__n)); 2103 2104 __is.flags(__flags); 2105 return __is; 2106 } 2107 2108 2109 template<typename _RealType> 2110 template<typename _UniformRandomNumberGenerator> 2111 typename cauchy_distribution<_RealType>::result_type 2112 cauchy_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)2113 operator()(_UniformRandomNumberGenerator& __urng, 2114 const param_type& __p) 2115 { 2116 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2117 __aurng(__urng); 2118 _RealType __u; 2119 do 2120 __u = __aurng(); 2121 while (__u == 0.5); 2122 2123 const _RealType __pi = 3.1415926535897932384626433832795029L; 2124 return __p.a() + __p.b() * std::tan(__pi * __u); 2125 } 2126 2127 template<typename _RealType> 2128 template<typename _ForwardIterator, 2129 typename _UniformRandomNumberGenerator> 2130 void 2131 cauchy_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)2132 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2133 _UniformRandomNumberGenerator& __urng, 2134 const param_type& __p) 2135 { 2136 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2137 const _RealType __pi = 3.1415926535897932384626433832795029L; 2138 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2139 __aurng(__urng); 2140 while (__f != __t) 2141 { 2142 _RealType __u; 2143 do 2144 __u = __aurng(); 2145 while (__u == 0.5); 2146 2147 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u); 2148 } 2149 } 2150 2151 template<typename _RealType, typename _CharT, typename _Traits> 2152 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const cauchy_distribution<_RealType> & __x)2153 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2154 const cauchy_distribution<_RealType>& __x) 2155 { 2156 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2157 2158 const typename __ios_base::fmtflags __flags = __os.flags(); 2159 const _CharT __fill = __os.fill(); 2160 const std::streamsize __precision = __os.precision(); 2161 const _CharT __space = __os.widen(' '); 2162 __os.flags(__ios_base::scientific | __ios_base::left); 2163 __os.fill(__space); 2164 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2165 2166 __os << __x.a() << __space << __x.b(); 2167 2168 __os.flags(__flags); 2169 __os.fill(__fill); 2170 __os.precision(__precision); 2171 return __os; 2172 } 2173 2174 template<typename _RealType, typename _CharT, typename _Traits> 2175 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,cauchy_distribution<_RealType> & __x)2176 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2177 cauchy_distribution<_RealType>& __x) 2178 { 2179 using param_type = typename cauchy_distribution<_RealType>::param_type; 2180 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2181 2182 const typename __ios_base::fmtflags __flags = __is.flags(); 2183 __is.flags(__ios_base::dec | __ios_base::skipws); 2184 2185 _RealType __a, __b; 2186 if (__is >> __a >> __b) 2187 __x.param(param_type(__a, __b)); 2188 2189 __is.flags(__flags); 2190 return __is; 2191 } 2192 2193 2194 template<typename _RealType> 2195 template<typename _ForwardIterator, 2196 typename _UniformRandomNumberGenerator> 2197 void 2198 std::fisher_f_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng)2199 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2200 _UniformRandomNumberGenerator& __urng) 2201 { 2202 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2203 while (__f != __t) 2204 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m())); 2205 } 2206 2207 template<typename _RealType> 2208 template<typename _ForwardIterator, 2209 typename _UniformRandomNumberGenerator> 2210 void 2211 std::fisher_f_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)2212 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2213 _UniformRandomNumberGenerator& __urng, 2214 const param_type& __p) 2215 { 2216 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2217 typedef typename std::gamma_distribution<result_type>::param_type 2218 param_type; 2219 param_type __p1(__p.m() / 2); 2220 param_type __p2(__p.n() / 2); 2221 while (__f != __t) 2222 *__f++ = ((_M_gd_x(__urng, __p1) * n()) 2223 / (_M_gd_y(__urng, __p2) * m())); 2224 } 2225 2226 template<typename _RealType, typename _CharT, typename _Traits> 2227 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const fisher_f_distribution<_RealType> & __x)2228 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2229 const fisher_f_distribution<_RealType>& __x) 2230 { 2231 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2232 2233 const typename __ios_base::fmtflags __flags = __os.flags(); 2234 const _CharT __fill = __os.fill(); 2235 const std::streamsize __precision = __os.precision(); 2236 const _CharT __space = __os.widen(' '); 2237 __os.flags(__ios_base::scientific | __ios_base::left); 2238 __os.fill(__space); 2239 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2240 2241 __os << __x.m() << __space << __x.n() 2242 << __space << __x._M_gd_x << __space << __x._M_gd_y; 2243 2244 __os.flags(__flags); 2245 __os.fill(__fill); 2246 __os.precision(__precision); 2247 return __os; 2248 } 2249 2250 template<typename _RealType, typename _CharT, typename _Traits> 2251 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,fisher_f_distribution<_RealType> & __x)2252 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2253 fisher_f_distribution<_RealType>& __x) 2254 { 2255 using param_type 2256 = typename fisher_f_distribution<_RealType>::param_type; 2257 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2258 2259 const typename __ios_base::fmtflags __flags = __is.flags(); 2260 __is.flags(__ios_base::dec | __ios_base::skipws); 2261 2262 _RealType __m, __n; 2263 if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y) 2264 __x.param(param_type(__m, __n)); 2265 2266 __is.flags(__flags); 2267 return __is; 2268 } 2269 2270 2271 template<typename _RealType> 2272 template<typename _ForwardIterator, 2273 typename _UniformRandomNumberGenerator> 2274 void 2275 std::student_t_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng)2276 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2277 _UniformRandomNumberGenerator& __urng) 2278 { 2279 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2280 while (__f != __t) 2281 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng)); 2282 } 2283 2284 template<typename _RealType> 2285 template<typename _ForwardIterator, 2286 typename _UniformRandomNumberGenerator> 2287 void 2288 std::student_t_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)2289 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2290 _UniformRandomNumberGenerator& __urng, 2291 const param_type& __p) 2292 { 2293 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2294 typename std::gamma_distribution<result_type>::param_type 2295 __p2(__p.n() / 2, 2); 2296 while (__f != __t) 2297 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2)); 2298 } 2299 2300 template<typename _RealType, typename _CharT, typename _Traits> 2301 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const student_t_distribution<_RealType> & __x)2302 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2303 const student_t_distribution<_RealType>& __x) 2304 { 2305 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2306 2307 const typename __ios_base::fmtflags __flags = __os.flags(); 2308 const _CharT __fill = __os.fill(); 2309 const std::streamsize __precision = __os.precision(); 2310 const _CharT __space = __os.widen(' '); 2311 __os.flags(__ios_base::scientific | __ios_base::left); 2312 __os.fill(__space); 2313 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2314 2315 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd; 2316 2317 __os.flags(__flags); 2318 __os.fill(__fill); 2319 __os.precision(__precision); 2320 return __os; 2321 } 2322 2323 template<typename _RealType, typename _CharT, typename _Traits> 2324 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,student_t_distribution<_RealType> & __x)2325 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2326 student_t_distribution<_RealType>& __x) 2327 { 2328 using param_type 2329 = typename student_t_distribution<_RealType>::param_type; 2330 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2331 2332 const typename __ios_base::fmtflags __flags = __is.flags(); 2333 __is.flags(__ios_base::dec | __ios_base::skipws); 2334 2335 _RealType __n; 2336 if (__is >> __n >> __x._M_nd >> __x._M_gd) 2337 __x.param(param_type(__n)); 2338 2339 __is.flags(__flags); 2340 return __is; 2341 } 2342 2343 2344 template<typename _RealType> 2345 void 2346 gamma_distribution<_RealType>::param_type:: _M_initialize()2347 _M_initialize() 2348 { 2349 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha; 2350 2351 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0); 2352 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1); 2353 } 2354 2355 /** 2356 * Marsaglia, G. and Tsang, W. W. 2357 * "A Simple Method for Generating Gamma Variables" 2358 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000. 2359 */ 2360 template<typename _RealType> 2361 template<typename _UniformRandomNumberGenerator> 2362 typename gamma_distribution<_RealType>::result_type 2363 gamma_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)2364 operator()(_UniformRandomNumberGenerator& __urng, 2365 const param_type& __param) 2366 { 2367 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2368 __aurng(__urng); 2369 2370 result_type __u, __v, __n; 2371 const result_type __a1 = (__param._M_malpha 2372 - _RealType(1.0) / _RealType(3.0)); 2373 2374 do 2375 { 2376 do 2377 { 2378 __n = _M_nd(__urng); 2379 __v = result_type(1.0) + __param._M_a2 * __n; 2380 } 2381 while (__v <= 0.0); 2382 2383 __v = __v * __v * __v; 2384 __u = __aurng(); 2385 } 2386 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n 2387 && (std::log(__u) > (0.5 * __n * __n + __a1 2388 * (1.0 - __v + std::log(__v))))); 2389 2390 if (__param.alpha() == __param._M_malpha) 2391 return __a1 * __v * __param.beta(); 2392 else 2393 { 2394 do 2395 __u = __aurng(); 2396 while (__u == 0.0); 2397 2398 return (std::pow(__u, result_type(1.0) / __param.alpha()) 2399 * __a1 * __v * __param.beta()); 2400 } 2401 } 2402 2403 template<typename _RealType> 2404 template<typename _ForwardIterator, 2405 typename _UniformRandomNumberGenerator> 2406 void 2407 gamma_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)2408 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2409 _UniformRandomNumberGenerator& __urng, 2410 const param_type& __param) 2411 { 2412 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2413 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2414 __aurng(__urng); 2415 2416 result_type __u, __v, __n; 2417 const result_type __a1 = (__param._M_malpha 2418 - _RealType(1.0) / _RealType(3.0)); 2419 2420 if (__param.alpha() == __param._M_malpha) 2421 while (__f != __t) 2422 { 2423 do 2424 { 2425 do 2426 { 2427 __n = _M_nd(__urng); 2428 __v = result_type(1.0) + __param._M_a2 * __n; 2429 } 2430 while (__v <= 0.0); 2431 2432 __v = __v * __v * __v; 2433 __u = __aurng(); 2434 } 2435 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n 2436 && (std::log(__u) > (0.5 * __n * __n + __a1 2437 * (1.0 - __v + std::log(__v))))); 2438 2439 *__f++ = __a1 * __v * __param.beta(); 2440 } 2441 else 2442 while (__f != __t) 2443 { 2444 do 2445 { 2446 do 2447 { 2448 __n = _M_nd(__urng); 2449 __v = result_type(1.0) + __param._M_a2 * __n; 2450 } 2451 while (__v <= 0.0); 2452 2453 __v = __v * __v * __v; 2454 __u = __aurng(); 2455 } 2456 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n 2457 && (std::log(__u) > (0.5 * __n * __n + __a1 2458 * (1.0 - __v + std::log(__v))))); 2459 2460 do 2461 __u = __aurng(); 2462 while (__u == 0.0); 2463 2464 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha()) 2465 * __a1 * __v * __param.beta()); 2466 } 2467 } 2468 2469 template<typename _RealType, typename _CharT, typename _Traits> 2470 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const gamma_distribution<_RealType> & __x)2471 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2472 const gamma_distribution<_RealType>& __x) 2473 { 2474 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2475 2476 const typename __ios_base::fmtflags __flags = __os.flags(); 2477 const _CharT __fill = __os.fill(); 2478 const std::streamsize __precision = __os.precision(); 2479 const _CharT __space = __os.widen(' '); 2480 __os.flags(__ios_base::scientific | __ios_base::left); 2481 __os.fill(__space); 2482 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2483 2484 __os << __x.alpha() << __space << __x.beta() 2485 << __space << __x._M_nd; 2486 2487 __os.flags(__flags); 2488 __os.fill(__fill); 2489 __os.precision(__precision); 2490 return __os; 2491 } 2492 2493 template<typename _RealType, typename _CharT, typename _Traits> 2494 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,gamma_distribution<_RealType> & __x)2495 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2496 gamma_distribution<_RealType>& __x) 2497 { 2498 using param_type = typename gamma_distribution<_RealType>::param_type; 2499 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2500 2501 const typename __ios_base::fmtflags __flags = __is.flags(); 2502 __is.flags(__ios_base::dec | __ios_base::skipws); 2503 2504 _RealType __alpha_val, __beta_val; 2505 if (__is >> __alpha_val >> __beta_val >> __x._M_nd) 2506 __x.param(param_type(__alpha_val, __beta_val)); 2507 2508 __is.flags(__flags); 2509 return __is; 2510 } 2511 2512 2513 template<typename _RealType> 2514 template<typename _UniformRandomNumberGenerator> 2515 typename weibull_distribution<_RealType>::result_type 2516 weibull_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)2517 operator()(_UniformRandomNumberGenerator& __urng, 2518 const param_type& __p) 2519 { 2520 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2521 __aurng(__urng); 2522 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 2523 result_type(1) / __p.a()); 2524 } 2525 2526 template<typename _RealType> 2527 template<typename _ForwardIterator, 2528 typename _UniformRandomNumberGenerator> 2529 void 2530 weibull_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)2531 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2532 _UniformRandomNumberGenerator& __urng, 2533 const param_type& __p) 2534 { 2535 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2536 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2537 __aurng(__urng); 2538 auto __inv_a = result_type(1) / __p.a(); 2539 2540 while (__f != __t) 2541 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 2542 __inv_a); 2543 } 2544 2545 template<typename _RealType, typename _CharT, typename _Traits> 2546 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const weibull_distribution<_RealType> & __x)2547 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2548 const weibull_distribution<_RealType>& __x) 2549 { 2550 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2551 2552 const typename __ios_base::fmtflags __flags = __os.flags(); 2553 const _CharT __fill = __os.fill(); 2554 const std::streamsize __precision = __os.precision(); 2555 const _CharT __space = __os.widen(' '); 2556 __os.flags(__ios_base::scientific | __ios_base::left); 2557 __os.fill(__space); 2558 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2559 2560 __os << __x.a() << __space << __x.b(); 2561 2562 __os.flags(__flags); 2563 __os.fill(__fill); 2564 __os.precision(__precision); 2565 return __os; 2566 } 2567 2568 template<typename _RealType, typename _CharT, typename _Traits> 2569 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,weibull_distribution<_RealType> & __x)2570 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2571 weibull_distribution<_RealType>& __x) 2572 { 2573 using param_type = typename weibull_distribution<_RealType>::param_type; 2574 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2575 2576 const typename __ios_base::fmtflags __flags = __is.flags(); 2577 __is.flags(__ios_base::dec | __ios_base::skipws); 2578 2579 _RealType __a, __b; 2580 if (__is >> __a >> __b) 2581 __x.param(param_type(__a, __b)); 2582 2583 __is.flags(__flags); 2584 return __is; 2585 } 2586 2587 2588 template<typename _RealType> 2589 template<typename _UniformRandomNumberGenerator> 2590 typename extreme_value_distribution<_RealType>::result_type 2591 extreme_value_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)2592 operator()(_UniformRandomNumberGenerator& __urng, 2593 const param_type& __p) 2594 { 2595 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2596 __aurng(__urng); 2597 return __p.a() - __p.b() * std::log(-std::log(result_type(1) 2598 - __aurng())); 2599 } 2600 2601 template<typename _RealType> 2602 template<typename _ForwardIterator, 2603 typename _UniformRandomNumberGenerator> 2604 void 2605 extreme_value_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)2606 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2607 _UniformRandomNumberGenerator& __urng, 2608 const param_type& __p) 2609 { 2610 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2611 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2612 __aurng(__urng); 2613 2614 while (__f != __t) 2615 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1) 2616 - __aurng())); 2617 } 2618 2619 template<typename _RealType, typename _CharT, typename _Traits> 2620 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const extreme_value_distribution<_RealType> & __x)2621 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2622 const extreme_value_distribution<_RealType>& __x) 2623 { 2624 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2625 2626 const typename __ios_base::fmtflags __flags = __os.flags(); 2627 const _CharT __fill = __os.fill(); 2628 const std::streamsize __precision = __os.precision(); 2629 const _CharT __space = __os.widen(' '); 2630 __os.flags(__ios_base::scientific | __ios_base::left); 2631 __os.fill(__space); 2632 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2633 2634 __os << __x.a() << __space << __x.b(); 2635 2636 __os.flags(__flags); 2637 __os.fill(__fill); 2638 __os.precision(__precision); 2639 return __os; 2640 } 2641 2642 template<typename _RealType, typename _CharT, typename _Traits> 2643 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,extreme_value_distribution<_RealType> & __x)2644 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2645 extreme_value_distribution<_RealType>& __x) 2646 { 2647 using param_type 2648 = typename extreme_value_distribution<_RealType>::param_type; 2649 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2650 2651 const typename __ios_base::fmtflags __flags = __is.flags(); 2652 __is.flags(__ios_base::dec | __ios_base::skipws); 2653 2654 _RealType __a, __b; 2655 if (__is >> __a >> __b) 2656 __x.param(param_type(__a, __b)); 2657 2658 __is.flags(__flags); 2659 return __is; 2660 } 2661 2662 2663 template<typename _IntType> 2664 void 2665 discrete_distribution<_IntType>::param_type:: _M_initialize()2666 _M_initialize() 2667 { 2668 if (_M_prob.size() < 2) 2669 { 2670 _M_prob.clear(); 2671 return; 2672 } 2673 2674 const double __sum = std::accumulate(_M_prob.begin(), 2675 _M_prob.end(), 0.0); 2676 __glibcxx_assert(__sum > 0); 2677 // Now normalize the probabilites. 2678 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), 2679 __sum); 2680 // Accumulate partial sums. 2681 _M_cp.reserve(_M_prob.size()); 2682 std::partial_sum(_M_prob.begin(), _M_prob.end(), 2683 std::back_inserter(_M_cp)); 2684 // Make sure the last cumulative probability is one. 2685 _M_cp[_M_cp.size() - 1] = 1.0; 2686 } 2687 2688 template<typename _IntType> 2689 template<typename _Func> 2690 discrete_distribution<_IntType>::param_type:: param_type(size_t __nw,double __xmin,double __xmax,_Func __fw)2691 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw) 2692 : _M_prob(), _M_cp() 2693 { 2694 const size_t __n = __nw == 0 ? 1 : __nw; 2695 const double __delta = (__xmax - __xmin) / __n; 2696 2697 _M_prob.reserve(__n); 2698 for (size_t __k = 0; __k < __nw; ++__k) 2699 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta)); 2700 2701 _M_initialize(); 2702 } 2703 2704 template<typename _IntType> 2705 template<typename _UniformRandomNumberGenerator> 2706 typename discrete_distribution<_IntType>::result_type 2707 discrete_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)2708 operator()(_UniformRandomNumberGenerator& __urng, 2709 const param_type& __param) 2710 { 2711 if (__param._M_cp.empty()) 2712 return result_type(0); 2713 2714 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2715 __aurng(__urng); 2716 2717 const double __p = __aurng(); 2718 auto __pos = std::lower_bound(__param._M_cp.begin(), 2719 __param._M_cp.end(), __p); 2720 2721 return __pos - __param._M_cp.begin(); 2722 } 2723 2724 template<typename _IntType> 2725 template<typename _ForwardIterator, 2726 typename _UniformRandomNumberGenerator> 2727 void 2728 discrete_distribution<_IntType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)2729 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2730 _UniformRandomNumberGenerator& __urng, 2731 const param_type& __param) 2732 { 2733 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2734 2735 if (__param._M_cp.empty()) 2736 { 2737 while (__f != __t) 2738 *__f++ = result_type(0); 2739 return; 2740 } 2741 2742 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2743 __aurng(__urng); 2744 2745 while (__f != __t) 2746 { 2747 const double __p = __aurng(); 2748 auto __pos = std::lower_bound(__param._M_cp.begin(), 2749 __param._M_cp.end(), __p); 2750 2751 *__f++ = __pos - __param._M_cp.begin(); 2752 } 2753 } 2754 2755 template<typename _IntType, typename _CharT, typename _Traits> 2756 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const discrete_distribution<_IntType> & __x)2757 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2758 const discrete_distribution<_IntType>& __x) 2759 { 2760 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2761 2762 const typename __ios_base::fmtflags __flags = __os.flags(); 2763 const _CharT __fill = __os.fill(); 2764 const std::streamsize __precision = __os.precision(); 2765 const _CharT __space = __os.widen(' '); 2766 __os.flags(__ios_base::scientific | __ios_base::left); 2767 __os.fill(__space); 2768 __os.precision(std::numeric_limits<double>::max_digits10); 2769 2770 std::vector<double> __prob = __x.probabilities(); 2771 __os << __prob.size(); 2772 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit) 2773 __os << __space << *__dit; 2774 2775 __os.flags(__flags); 2776 __os.fill(__fill); 2777 __os.precision(__precision); 2778 return __os; 2779 } 2780 2781 namespace __detail 2782 { 2783 template<typename _ValT, typename _CharT, typename _Traits> 2784 basic_istream<_CharT, _Traits>& __extract_params(basic_istream<_CharT,_Traits> & __is,vector<_ValT> & __vals,size_t __n)2785 __extract_params(basic_istream<_CharT, _Traits>& __is, 2786 vector<_ValT>& __vals, size_t __n) 2787 { 2788 __vals.reserve(__n); 2789 while (__n--) 2790 { 2791 _ValT __val; 2792 if (__is >> __val) 2793 __vals.push_back(__val); 2794 else 2795 break; 2796 } 2797 return __is; 2798 } 2799 } // namespace __detail 2800 2801 template<typename _IntType, typename _CharT, typename _Traits> 2802 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,discrete_distribution<_IntType> & __x)2803 operator>>(std::basic_istream<_CharT, _Traits>& __is, 2804 discrete_distribution<_IntType>& __x) 2805 { 2806 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 2807 2808 const typename __ios_base::fmtflags __flags = __is.flags(); 2809 __is.flags(__ios_base::dec | __ios_base::skipws); 2810 2811 size_t __n; 2812 if (__is >> __n) 2813 { 2814 std::vector<double> __prob_vec; 2815 if (__detail::__extract_params(__is, __prob_vec, __n)) 2816 __x.param({__prob_vec.begin(), __prob_vec.end()}); 2817 } 2818 2819 __is.flags(__flags); 2820 return __is; 2821 } 2822 2823 2824 template<typename _RealType> 2825 void 2826 piecewise_constant_distribution<_RealType>::param_type:: _M_initialize()2827 _M_initialize() 2828 { 2829 if (_M_int.size() < 2 2830 || (_M_int.size() == 2 2831 && _M_int[0] == _RealType(0) 2832 && _M_int[1] == _RealType(1))) 2833 { 2834 _M_int.clear(); 2835 _M_den.clear(); 2836 return; 2837 } 2838 2839 const double __sum = std::accumulate(_M_den.begin(), 2840 _M_den.end(), 0.0); 2841 __glibcxx_assert(__sum > 0); 2842 2843 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 2844 __sum); 2845 2846 _M_cp.reserve(_M_den.size()); 2847 std::partial_sum(_M_den.begin(), _M_den.end(), 2848 std::back_inserter(_M_cp)); 2849 2850 // Make sure the last cumulative probability is one. 2851 _M_cp[_M_cp.size() - 1] = 1.0; 2852 2853 for (size_t __k = 0; __k < _M_den.size(); ++__k) 2854 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k]; 2855 } 2856 2857 template<typename _RealType> 2858 template<typename _InputIteratorB, typename _InputIteratorW> 2859 piecewise_constant_distribution<_RealType>::param_type:: param_type(_InputIteratorB __bbegin,_InputIteratorB __bend,_InputIteratorW __wbegin)2860 param_type(_InputIteratorB __bbegin, 2861 _InputIteratorB __bend, 2862 _InputIteratorW __wbegin) 2863 : _M_int(), _M_den(), _M_cp() 2864 { 2865 if (__bbegin != __bend) 2866 { 2867 for (;;) 2868 { 2869 _M_int.push_back(*__bbegin); 2870 ++__bbegin; 2871 if (__bbegin == __bend) 2872 break; 2873 2874 _M_den.push_back(*__wbegin); 2875 ++__wbegin; 2876 } 2877 } 2878 2879 _M_initialize(); 2880 } 2881 2882 template<typename _RealType> 2883 template<typename _Func> 2884 piecewise_constant_distribution<_RealType>::param_type:: param_type(initializer_list<_RealType> __bl,_Func __fw)2885 param_type(initializer_list<_RealType> __bl, _Func __fw) 2886 : _M_int(), _M_den(), _M_cp() 2887 { 2888 _M_int.reserve(__bl.size()); 2889 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 2890 _M_int.push_back(*__biter); 2891 2892 _M_den.reserve(_M_int.size() - 1); 2893 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 2894 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k]))); 2895 2896 _M_initialize(); 2897 } 2898 2899 template<typename _RealType> 2900 template<typename _Func> 2901 piecewise_constant_distribution<_RealType>::param_type:: param_type(size_t __nw,_RealType __xmin,_RealType __xmax,_Func __fw)2902 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 2903 : _M_int(), _M_den(), _M_cp() 2904 { 2905 const size_t __n = __nw == 0 ? 1 : __nw; 2906 const _RealType __delta = (__xmax - __xmin) / __n; 2907 2908 _M_int.reserve(__n + 1); 2909 for (size_t __k = 0; __k <= __nw; ++__k) 2910 _M_int.push_back(__xmin + __k * __delta); 2911 2912 _M_den.reserve(__n); 2913 for (size_t __k = 0; __k < __nw; ++__k) 2914 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta)); 2915 2916 _M_initialize(); 2917 } 2918 2919 template<typename _RealType> 2920 template<typename _UniformRandomNumberGenerator> 2921 typename piecewise_constant_distribution<_RealType>::result_type 2922 piecewise_constant_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)2923 operator()(_UniformRandomNumberGenerator& __urng, 2924 const param_type& __param) 2925 { 2926 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2927 __aurng(__urng); 2928 2929 const double __p = __aurng(); 2930 if (__param._M_cp.empty()) 2931 return __p; 2932 2933 auto __pos = std::lower_bound(__param._M_cp.begin(), 2934 __param._M_cp.end(), __p); 2935 const size_t __i = __pos - __param._M_cp.begin(); 2936 2937 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 2938 2939 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i]; 2940 } 2941 2942 template<typename _RealType> 2943 template<typename _ForwardIterator, 2944 typename _UniformRandomNumberGenerator> 2945 void 2946 piecewise_constant_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)2947 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 2948 _UniformRandomNumberGenerator& __urng, 2949 const param_type& __param) 2950 { 2951 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 2952 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2953 __aurng(__urng); 2954 2955 if (__param._M_cp.empty()) 2956 { 2957 while (__f != __t) 2958 *__f++ = __aurng(); 2959 return; 2960 } 2961 2962 while (__f != __t) 2963 { 2964 const double __p = __aurng(); 2965 2966 auto __pos = std::lower_bound(__param._M_cp.begin(), 2967 __param._M_cp.end(), __p); 2968 const size_t __i = __pos - __param._M_cp.begin(); 2969 2970 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 2971 2972 *__f++ = (__param._M_int[__i] 2973 + (__p - __pref) / __param._M_den[__i]); 2974 } 2975 } 2976 2977 template<typename _RealType, typename _CharT, typename _Traits> 2978 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const piecewise_constant_distribution<_RealType> & __x)2979 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2980 const piecewise_constant_distribution<_RealType>& __x) 2981 { 2982 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 2983 2984 const typename __ios_base::fmtflags __flags = __os.flags(); 2985 const _CharT __fill = __os.fill(); 2986 const std::streamsize __precision = __os.precision(); 2987 const _CharT __space = __os.widen(' '); 2988 __os.flags(__ios_base::scientific | __ios_base::left); 2989 __os.fill(__space); 2990 __os.precision(std::numeric_limits<_RealType>::max_digits10); 2991 2992 std::vector<_RealType> __int = __x.intervals(); 2993 __os << __int.size() - 1; 2994 2995 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 2996 __os << __space << *__xit; 2997 2998 std::vector<double> __den = __x.densities(); 2999 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 3000 __os << __space << *__dit; 3001 3002 __os.flags(__flags); 3003 __os.fill(__fill); 3004 __os.precision(__precision); 3005 return __os; 3006 } 3007 3008 template<typename _RealType, typename _CharT, typename _Traits> 3009 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,piecewise_constant_distribution<_RealType> & __x)3010 operator>>(std::basic_istream<_CharT, _Traits>& __is, 3011 piecewise_constant_distribution<_RealType>& __x) 3012 { 3013 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 3014 3015 const typename __ios_base::fmtflags __flags = __is.flags(); 3016 __is.flags(__ios_base::dec | __ios_base::skipws); 3017 3018 size_t __n; 3019 if (__is >> __n) 3020 { 3021 std::vector<_RealType> __int_vec; 3022 if (__detail::__extract_params(__is, __int_vec, __n + 1)) 3023 { 3024 std::vector<double> __den_vec; 3025 if (__detail::__extract_params(__is, __den_vec, __n)) 3026 { 3027 __x.param({ __int_vec.begin(), __int_vec.end(), 3028 __den_vec.begin() }); 3029 } 3030 } 3031 } 3032 3033 __is.flags(__flags); 3034 return __is; 3035 } 3036 3037 3038 template<typename _RealType> 3039 void 3040 piecewise_linear_distribution<_RealType>::param_type:: _M_initialize()3041 _M_initialize() 3042 { 3043 if (_M_int.size() < 2 3044 || (_M_int.size() == 2 3045 && _M_int[0] == _RealType(0) 3046 && _M_int[1] == _RealType(1) 3047 && _M_den[0] == _M_den[1])) 3048 { 3049 _M_int.clear(); 3050 _M_den.clear(); 3051 return; 3052 } 3053 3054 double __sum = 0.0; 3055 _M_cp.reserve(_M_int.size() - 1); 3056 _M_m.reserve(_M_int.size() - 1); 3057 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 3058 { 3059 const _RealType __delta = _M_int[__k + 1] - _M_int[__k]; 3060 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta; 3061 _M_cp.push_back(__sum); 3062 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta); 3063 } 3064 __glibcxx_assert(__sum > 0); 3065 3066 // Now normalize the densities... 3067 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 3068 __sum); 3069 // ... and partial sums... 3070 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum); 3071 // ... and slopes. 3072 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum); 3073 3074 // Make sure the last cumulative probablility is one. 3075 _M_cp[_M_cp.size() - 1] = 1.0; 3076 } 3077 3078 template<typename _RealType> 3079 template<typename _InputIteratorB, typename _InputIteratorW> 3080 piecewise_linear_distribution<_RealType>::param_type:: param_type(_InputIteratorB __bbegin,_InputIteratorB __bend,_InputIteratorW __wbegin)3081 param_type(_InputIteratorB __bbegin, 3082 _InputIteratorB __bend, 3083 _InputIteratorW __wbegin) 3084 : _M_int(), _M_den(), _M_cp(), _M_m() 3085 { 3086 for (; __bbegin != __bend; ++__bbegin, ++__wbegin) 3087 { 3088 _M_int.push_back(*__bbegin); 3089 _M_den.push_back(*__wbegin); 3090 } 3091 3092 _M_initialize(); 3093 } 3094 3095 template<typename _RealType> 3096 template<typename _Func> 3097 piecewise_linear_distribution<_RealType>::param_type:: param_type(initializer_list<_RealType> __bl,_Func __fw)3098 param_type(initializer_list<_RealType> __bl, _Func __fw) 3099 : _M_int(), _M_den(), _M_cp(), _M_m() 3100 { 3101 _M_int.reserve(__bl.size()); 3102 _M_den.reserve(__bl.size()); 3103 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 3104 { 3105 _M_int.push_back(*__biter); 3106 _M_den.push_back(__fw(*__biter)); 3107 } 3108 3109 _M_initialize(); 3110 } 3111 3112 template<typename _RealType> 3113 template<typename _Func> 3114 piecewise_linear_distribution<_RealType>::param_type:: param_type(size_t __nw,_RealType __xmin,_RealType __xmax,_Func __fw)3115 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 3116 : _M_int(), _M_den(), _M_cp(), _M_m() 3117 { 3118 const size_t __n = __nw == 0 ? 1 : __nw; 3119 const _RealType __delta = (__xmax - __xmin) / __n; 3120 3121 _M_int.reserve(__n + 1); 3122 _M_den.reserve(__n + 1); 3123 for (size_t __k = 0; __k <= __nw; ++__k) 3124 { 3125 _M_int.push_back(__xmin + __k * __delta); 3126 _M_den.push_back(__fw(_M_int[__k] + __delta)); 3127 } 3128 3129 _M_initialize(); 3130 } 3131 3132 template<typename _RealType> 3133 template<typename _UniformRandomNumberGenerator> 3134 typename piecewise_linear_distribution<_RealType>::result_type 3135 piecewise_linear_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)3136 operator()(_UniformRandomNumberGenerator& __urng, 3137 const param_type& __param) 3138 { 3139 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 3140 __aurng(__urng); 3141 3142 const double __p = __aurng(); 3143 if (__param._M_cp.empty()) 3144 return __p; 3145 3146 auto __pos = std::lower_bound(__param._M_cp.begin(), 3147 __param._M_cp.end(), __p); 3148 const size_t __i = __pos - __param._M_cp.begin(); 3149 3150 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 3151 3152 const double __a = 0.5 * __param._M_m[__i]; 3153 const double __b = __param._M_den[__i]; 3154 const double __cm = __p - __pref; 3155 3156 _RealType __x = __param._M_int[__i]; 3157 if (__a == 0) 3158 __x += __cm / __b; 3159 else 3160 { 3161 const double __d = __b * __b + 4.0 * __a * __cm; 3162 __x += 0.5 * (std::sqrt(__d) - __b) / __a; 3163 } 3164 3165 return __x; 3166 } 3167 3168 template<typename _RealType> 3169 template<typename _ForwardIterator, 3170 typename _UniformRandomNumberGenerator> 3171 void 3172 piecewise_linear_distribution<_RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)3173 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 3174 _UniformRandomNumberGenerator& __urng, 3175 const param_type& __param) 3176 { 3177 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 3178 // We could duplicate everything from operator()... 3179 while (__f != __t) 3180 *__f++ = this->operator()(__urng, __param); 3181 } 3182 3183 template<typename _RealType, typename _CharT, typename _Traits> 3184 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const piecewise_linear_distribution<_RealType> & __x)3185 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 3186 const piecewise_linear_distribution<_RealType>& __x) 3187 { 3188 using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; 3189 3190 const typename __ios_base::fmtflags __flags = __os.flags(); 3191 const _CharT __fill = __os.fill(); 3192 const std::streamsize __precision = __os.precision(); 3193 const _CharT __space = __os.widen(' '); 3194 __os.flags(__ios_base::scientific | __ios_base::left); 3195 __os.fill(__space); 3196 __os.precision(std::numeric_limits<_RealType>::max_digits10); 3197 3198 std::vector<_RealType> __int = __x.intervals(); 3199 __os << __int.size() - 1; 3200 3201 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 3202 __os << __space << *__xit; 3203 3204 std::vector<double> __den = __x.densities(); 3205 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 3206 __os << __space << *__dit; 3207 3208 __os.flags(__flags); 3209 __os.fill(__fill); 3210 __os.precision(__precision); 3211 return __os; 3212 } 3213 3214 template<typename _RealType, typename _CharT, typename _Traits> 3215 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,piecewise_linear_distribution<_RealType> & __x)3216 operator>>(std::basic_istream<_CharT, _Traits>& __is, 3217 piecewise_linear_distribution<_RealType>& __x) 3218 { 3219 using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; 3220 3221 const typename __ios_base::fmtflags __flags = __is.flags(); 3222 __is.flags(__ios_base::dec | __ios_base::skipws); 3223 3224 size_t __n; 3225 if (__is >> __n) 3226 { 3227 vector<_RealType> __int_vec; 3228 if (__detail::__extract_params(__is, __int_vec, __n + 1)) 3229 { 3230 vector<double> __den_vec; 3231 if (__detail::__extract_params(__is, __den_vec, __n + 1)) 3232 { 3233 __x.param({ __int_vec.begin(), __int_vec.end(), 3234 __den_vec.begin() }); 3235 } 3236 } 3237 } 3238 __is.flags(__flags); 3239 return __is; 3240 } 3241 3242 3243 template<typename _IntType, typename> seed_seq(std::initializer_list<_IntType> __il)3244 seed_seq::seed_seq(std::initializer_list<_IntType> __il) 3245 { 3246 _M_v.reserve(__il.size()); 3247 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter) 3248 _M_v.push_back(__detail::__mod<result_type, 3249 __detail::_Shift<result_type, 32>::__value>(*__iter)); 3250 } 3251 3252 template<typename _InputIterator> seed_seq(_InputIterator __begin,_InputIterator __end)3253 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end) 3254 { 3255 if _GLIBCXX17_CONSTEXPR (__is_random_access_iter<_InputIterator>::value) 3256 _M_v.reserve(std::distance(__begin, __end)); 3257 3258 for (_InputIterator __iter = __begin; __iter != __end; ++__iter) 3259 _M_v.push_back(__detail::__mod<result_type, 3260 __detail::_Shift<result_type, 32>::__value>(*__iter)); 3261 } 3262 3263 template<typename _RandomAccessIterator> 3264 void generate(_RandomAccessIterator __begin,_RandomAccessIterator __end)3265 seed_seq::generate(_RandomAccessIterator __begin, 3266 _RandomAccessIterator __end) 3267 { 3268 typedef typename iterator_traits<_RandomAccessIterator>::value_type 3269 _Type; 3270 3271 if (__begin == __end) 3272 return; 3273 3274 std::fill(__begin, __end, _Type(0x8b8b8b8bu)); 3275 3276 const size_t __n = __end - __begin; 3277 const size_t __s = _M_v.size(); 3278 const size_t __t = (__n >= 623) ? 11 3279 : (__n >= 68) ? 7 3280 : (__n >= 39) ? 5 3281 : (__n >= 7) ? 3 3282 : (__n - 1) / 2; 3283 const size_t __p = (__n - __t) / 2; 3284 const size_t __q = __p + __t; 3285 const size_t __m = std::max(size_t(__s + 1), __n); 3286 3287 #ifndef __UINT32_TYPE__ 3288 struct _Up 3289 { 3290 _Up(uint_least32_t v) : _M_v(v & 0xffffffffu) { } 3291 3292 operator uint_least32_t() const { return _M_v; } 3293 3294 uint_least32_t _M_v; 3295 }; 3296 using uint32_t = _Up; 3297 #endif 3298 3299 // k == 0, every element in [begin,end) equals 0x8b8b8b8bu 3300 { 3301 uint32_t __r1 = 1371501266u; 3302 uint32_t __r2 = __r1 + __s; 3303 __begin[__p] += __r1; 3304 __begin[__q] = (uint32_t)__begin[__q] + __r2; 3305 __begin[0] = __r2; 3306 } 3307 3308 for (size_t __k = 1; __k <= __s; ++__k) 3309 { 3310 const size_t __kn = __k % __n; 3311 const size_t __kpn = (__k + __p) % __n; 3312 const size_t __kqn = (__k + __q) % __n; 3313 uint32_t __arg = (__begin[__kn] 3314 ^ __begin[__kpn] 3315 ^ __begin[(__k - 1) % __n]); 3316 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27)); 3317 uint32_t __r2 = __r1 + (uint32_t)__kn + _M_v[__k - 1]; 3318 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1; 3319 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2; 3320 __begin[__kn] = __r2; 3321 } 3322 3323 for (size_t __k = __s + 1; __k < __m; ++__k) 3324 { 3325 const size_t __kn = __k % __n; 3326 const size_t __kpn = (__k + __p) % __n; 3327 const size_t __kqn = (__k + __q) % __n; 3328 uint32_t __arg = (__begin[__kn] 3329 ^ __begin[__kpn] 3330 ^ __begin[(__k - 1) % __n]); 3331 uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27)); 3332 uint32_t __r2 = __r1 + (uint32_t)__kn; 3333 __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1; 3334 __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2; 3335 __begin[__kn] = __r2; 3336 } 3337 3338 for (size_t __k = __m; __k < __m + __n; ++__k) 3339 { 3340 const size_t __kn = __k % __n; 3341 const size_t __kpn = (__k + __p) % __n; 3342 const size_t __kqn = (__k + __q) % __n; 3343 uint32_t __arg = (__begin[__kn] 3344 + __begin[__kpn] 3345 + __begin[(__k - 1) % __n]); 3346 uint32_t __r3 = 1566083941u * (__arg ^ (__arg >> 27)); 3347 uint32_t __r4 = __r3 - __kn; 3348 __begin[__kpn] ^= __r3; 3349 __begin[__kqn] ^= __r4; 3350 __begin[__kn] = __r4; 3351 } 3352 } 3353 3354 template<typename _RealType, size_t __bits, 3355 typename _UniformRandomNumberGenerator> 3356 _RealType generate_canonical(_UniformRandomNumberGenerator & __urng)3357 generate_canonical(_UniformRandomNumberGenerator& __urng) 3358 { 3359 static_assert(std::is_floating_point<_RealType>::value, 3360 "template argument must be a floating point type"); 3361 3362 const size_t __b 3363 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits), 3364 __bits); 3365 const long double __r = static_cast<long double>(__urng.max()) 3366 - static_cast<long double>(__urng.min()) + 1.0L; 3367 const size_t __log2r = std::log(__r) / std::log(2.0L); 3368 const size_t __m = std::max<size_t>(1UL, 3369 (__b + __log2r - 1UL) / __log2r); 3370 _RealType __ret; 3371 _RealType __sum = _RealType(0); 3372 _RealType __tmp = _RealType(1); 3373 for (size_t __k = __m; __k != 0; --__k) 3374 { 3375 __sum += _RealType(__urng() - __urng.min()) * __tmp; 3376 __tmp *= __r; 3377 } 3378 __ret = __sum / __tmp; 3379 if (__builtin_expect(__ret >= _RealType(1), 0)) 3380 { 3381 #if _GLIBCXX_USE_C99_MATH_TR1 3382 __ret = std::nextafter(_RealType(1), _RealType(0)); 3383 #else 3384 __ret = _RealType(1) 3385 - std::numeric_limits<_RealType>::epsilon() / _RealType(2); 3386 #endif 3387 } 3388 return __ret; 3389 } 3390 3391 _GLIBCXX_END_NAMESPACE_VERSION 3392 } // namespace 3393 3394 #endif 3395