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