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