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