random.tcc revision 1.4
1// Random number extensions -*- C++ -*- 2 3// Copyright (C) 2012-2016 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 ext/random.tcc 26 * This is an internal header file, included by other library headers. 27 * Do not attempt to use it directly. @headername{ext/random} 28 */ 29 30#ifndef _EXT_RANDOM_TCC 31#define _EXT_RANDOM_TCC 1 32 33#pragma GCC system_header 34 35namespace __gnu_cxx _GLIBCXX_VISIBILITY(default) 36{ 37_GLIBCXX_BEGIN_NAMESPACE_VERSION 38 39#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ 40 41 template<typename _UIntType, size_t __m, 42 size_t __pos1, size_t __sl1, size_t __sl2, 43 size_t __sr1, size_t __sr2, 44 uint32_t __msk1, uint32_t __msk2, 45 uint32_t __msk3, uint32_t __msk4, 46 uint32_t __parity1, uint32_t __parity2, 47 uint32_t __parity3, uint32_t __parity4> 48 void simd_fast_mersenne_twister_engine<_UIntType, __m, 49 __pos1, __sl1, __sl2, __sr1, __sr2, 50 __msk1, __msk2, __msk3, __msk4, 51 __parity1, __parity2, __parity3, 52 __parity4>:: 53 seed(_UIntType __seed) 54 { 55 _M_state32[0] = static_cast<uint32_t>(__seed); 56 for (size_t __i = 1; __i < _M_nstate32; ++__i) 57 _M_state32[__i] = (1812433253UL 58 * (_M_state32[__i - 1] ^ (_M_state32[__i - 1] >> 30)) 59 + __i); 60 _M_pos = state_size; 61 _M_period_certification(); 62 } 63 64 65 namespace { 66 67 inline uint32_t _Func1(uint32_t __x) 68 { 69 return (__x ^ (__x >> 27)) * UINT32_C(1664525); 70 } 71 72 inline uint32_t _Func2(uint32_t __x) 73 { 74 return (__x ^ (__x >> 27)) * UINT32_C(1566083941); 75 } 76 77 } 78 79 80 template<typename _UIntType, size_t __m, 81 size_t __pos1, size_t __sl1, size_t __sl2, 82 size_t __sr1, size_t __sr2, 83 uint32_t __msk1, uint32_t __msk2, 84 uint32_t __msk3, uint32_t __msk4, 85 uint32_t __parity1, uint32_t __parity2, 86 uint32_t __parity3, uint32_t __parity4> 87 template<typename _Sseq> 88 typename std::enable_if<std::is_class<_Sseq>::value>::type 89 simd_fast_mersenne_twister_engine<_UIntType, __m, 90 __pos1, __sl1, __sl2, __sr1, __sr2, 91 __msk1, __msk2, __msk3, __msk4, 92 __parity1, __parity2, __parity3, 93 __parity4>:: 94 seed(_Sseq& __q) 95 { 96 size_t __lag; 97 98 if (_M_nstate32 >= 623) 99 __lag = 11; 100 else if (_M_nstate32 >= 68) 101 __lag = 7; 102 else if (_M_nstate32 >= 39) 103 __lag = 5; 104 else 105 __lag = 3; 106 const size_t __mid = (_M_nstate32 - __lag) / 2; 107 108 std::fill(_M_state32, _M_state32 + _M_nstate32, UINT32_C(0x8b8b8b8b)); 109 uint32_t __arr[_M_nstate32]; 110 __q.generate(__arr + 0, __arr + _M_nstate32); 111 112 uint32_t __r = _Func1(_M_state32[0] ^ _M_state32[__mid] 113 ^ _M_state32[_M_nstate32 - 1]); 114 _M_state32[__mid] += __r; 115 __r += _M_nstate32; 116 _M_state32[__mid + __lag] += __r; 117 _M_state32[0] = __r; 118 119 for (size_t __i = 1, __j = 0; __j < _M_nstate32; ++__j) 120 { 121 __r = _Func1(_M_state32[__i] 122 ^ _M_state32[(__i + __mid) % _M_nstate32] 123 ^ _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]); 124 _M_state32[(__i + __mid) % _M_nstate32] += __r; 125 __r += __arr[__j] + __i; 126 _M_state32[(__i + __mid + __lag) % _M_nstate32] += __r; 127 _M_state32[__i] = __r; 128 __i = (__i + 1) % _M_nstate32; 129 } 130 for (size_t __j = 0; __j < _M_nstate32; ++__j) 131 { 132 const size_t __i = (__j + 1) % _M_nstate32; 133 __r = _Func2(_M_state32[__i] 134 + _M_state32[(__i + __mid) % _M_nstate32] 135 + _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]); 136 _M_state32[(__i + __mid) % _M_nstate32] ^= __r; 137 __r -= __i; 138 _M_state32[(__i + __mid + __lag) % _M_nstate32] ^= __r; 139 _M_state32[__i] = __r; 140 } 141 142 _M_pos = state_size; 143 _M_period_certification(); 144 } 145 146 147 template<typename _UIntType, size_t __m, 148 size_t __pos1, size_t __sl1, size_t __sl2, 149 size_t __sr1, size_t __sr2, 150 uint32_t __msk1, uint32_t __msk2, 151 uint32_t __msk3, uint32_t __msk4, 152 uint32_t __parity1, uint32_t __parity2, 153 uint32_t __parity3, uint32_t __parity4> 154 void simd_fast_mersenne_twister_engine<_UIntType, __m, 155 __pos1, __sl1, __sl2, __sr1, __sr2, 156 __msk1, __msk2, __msk3, __msk4, 157 __parity1, __parity2, __parity3, 158 __parity4>:: 159 _M_period_certification(void) 160 { 161 static const uint32_t __parity[4] = { __parity1, __parity2, 162 __parity3, __parity4 }; 163 uint32_t __inner = 0; 164 for (size_t __i = 0; __i < 4; ++__i) 165 if (__parity[__i] != 0) 166 __inner ^= _M_state32[__i] & __parity[__i]; 167 168 if (__builtin_parity(__inner) & 1) 169 return; 170 for (size_t __i = 0; __i < 4; ++__i) 171 if (__parity[__i] != 0) 172 { 173 _M_state32[__i] ^= 1 << (__builtin_ffs(__parity[__i]) - 1); 174 return; 175 } 176 __builtin_unreachable(); 177 } 178 179 180 template<typename _UIntType, size_t __m, 181 size_t __pos1, size_t __sl1, size_t __sl2, 182 size_t __sr1, size_t __sr2, 183 uint32_t __msk1, uint32_t __msk2, 184 uint32_t __msk3, uint32_t __msk4, 185 uint32_t __parity1, uint32_t __parity2, 186 uint32_t __parity3, uint32_t __parity4> 187 void simd_fast_mersenne_twister_engine<_UIntType, __m, 188 __pos1, __sl1, __sl2, __sr1, __sr2, 189 __msk1, __msk2, __msk3, __msk4, 190 __parity1, __parity2, __parity3, 191 __parity4>:: 192 discard(unsigned long long __z) 193 { 194 while (__z > state_size - _M_pos) 195 { 196 __z -= state_size - _M_pos; 197 198 _M_gen_rand(); 199 } 200 201 _M_pos += __z; 202 } 203 204 205#ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_GEN_READ 206 207 namespace { 208 209 template<size_t __shift> 210 inline void __rshift(uint32_t *__out, const uint32_t *__in) 211 { 212 uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32) 213 | static_cast<uint64_t>(__in[2])); 214 uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32) 215 | static_cast<uint64_t>(__in[0])); 216 217 uint64_t __oh = __th >> (__shift * 8); 218 uint64_t __ol = __tl >> (__shift * 8); 219 __ol |= __th << (64 - __shift * 8); 220 __out[1] = static_cast<uint32_t>(__ol >> 32); 221 __out[0] = static_cast<uint32_t>(__ol); 222 __out[3] = static_cast<uint32_t>(__oh >> 32); 223 __out[2] = static_cast<uint32_t>(__oh); 224 } 225 226 227 template<size_t __shift> 228 inline void __lshift(uint32_t *__out, const uint32_t *__in) 229 { 230 uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32) 231 | static_cast<uint64_t>(__in[2])); 232 uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32) 233 | static_cast<uint64_t>(__in[0])); 234 235 uint64_t __oh = __th << (__shift * 8); 236 uint64_t __ol = __tl << (__shift * 8); 237 __oh |= __tl >> (64 - __shift * 8); 238 __out[1] = static_cast<uint32_t>(__ol >> 32); 239 __out[0] = static_cast<uint32_t>(__ol); 240 __out[3] = static_cast<uint32_t>(__oh >> 32); 241 __out[2] = static_cast<uint32_t>(__oh); 242 } 243 244 245 template<size_t __sl1, size_t __sl2, size_t __sr1, size_t __sr2, 246 uint32_t __msk1, uint32_t __msk2, uint32_t __msk3, uint32_t __msk4> 247 inline void __recursion(uint32_t *__r, 248 const uint32_t *__a, const uint32_t *__b, 249 const uint32_t *__c, const uint32_t *__d) 250 { 251 uint32_t __x[4]; 252 uint32_t __y[4]; 253 254 __lshift<__sl2>(__x, __a); 255 __rshift<__sr2>(__y, __c); 256 __r[0] = (__a[0] ^ __x[0] ^ ((__b[0] >> __sr1) & __msk1) 257 ^ __y[0] ^ (__d[0] << __sl1)); 258 __r[1] = (__a[1] ^ __x[1] ^ ((__b[1] >> __sr1) & __msk2) 259 ^ __y[1] ^ (__d[1] << __sl1)); 260 __r[2] = (__a[2] ^ __x[2] ^ ((__b[2] >> __sr1) & __msk3) 261 ^ __y[2] ^ (__d[2] << __sl1)); 262 __r[3] = (__a[3] ^ __x[3] ^ ((__b[3] >> __sr1) & __msk4) 263 ^ __y[3] ^ (__d[3] << __sl1)); 264 } 265 266 } 267 268 269 template<typename _UIntType, size_t __m, 270 size_t __pos1, size_t __sl1, size_t __sl2, 271 size_t __sr1, size_t __sr2, 272 uint32_t __msk1, uint32_t __msk2, 273 uint32_t __msk3, uint32_t __msk4, 274 uint32_t __parity1, uint32_t __parity2, 275 uint32_t __parity3, uint32_t __parity4> 276 void simd_fast_mersenne_twister_engine<_UIntType, __m, 277 __pos1, __sl1, __sl2, __sr1, __sr2, 278 __msk1, __msk2, __msk3, __msk4, 279 __parity1, __parity2, __parity3, 280 __parity4>:: 281 _M_gen_rand(void) 282 { 283 const uint32_t *__r1 = &_M_state32[_M_nstate32 - 8]; 284 const uint32_t *__r2 = &_M_state32[_M_nstate32 - 4]; 285 static constexpr size_t __pos1_32 = __pos1 * 4; 286 287 size_t __i; 288 for (__i = 0; __i < _M_nstate32 - __pos1_32; __i += 4) 289 { 290 __recursion<__sl1, __sl2, __sr1, __sr2, 291 __msk1, __msk2, __msk3, __msk4> 292 (&_M_state32[__i], &_M_state32[__i], 293 &_M_state32[__i + __pos1_32], __r1, __r2); 294 __r1 = __r2; 295 __r2 = &_M_state32[__i]; 296 } 297 298 for (; __i < _M_nstate32; __i += 4) 299 { 300 __recursion<__sl1, __sl2, __sr1, __sr2, 301 __msk1, __msk2, __msk3, __msk4> 302 (&_M_state32[__i], &_M_state32[__i], 303 &_M_state32[__i + __pos1_32 - _M_nstate32], __r1, __r2); 304 __r1 = __r2; 305 __r2 = &_M_state32[__i]; 306 } 307 308 _M_pos = 0; 309 } 310 311#endif 312 313#ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_OPERATOREQUAL 314 template<typename _UIntType, size_t __m, 315 size_t __pos1, size_t __sl1, size_t __sl2, 316 size_t __sr1, size_t __sr2, 317 uint32_t __msk1, uint32_t __msk2, 318 uint32_t __msk3, uint32_t __msk4, 319 uint32_t __parity1, uint32_t __parity2, 320 uint32_t __parity3, uint32_t __parity4> 321 bool 322 operator==(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 323 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 324 __msk1, __msk2, __msk3, __msk4, 325 __parity1, __parity2, __parity3, __parity4>& __lhs, 326 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 327 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 328 __msk1, __msk2, __msk3, __msk4, 329 __parity1, __parity2, __parity3, __parity4>& __rhs) 330 { 331 typedef __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 332 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 333 __msk1, __msk2, __msk3, __msk4, 334 __parity1, __parity2, __parity3, __parity4> __engine; 335 return (std::equal(__lhs._M_stateT, 336 __lhs._M_stateT + __engine::state_size, 337 __rhs._M_stateT) 338 && __lhs._M_pos == __rhs._M_pos); 339 } 340#endif 341 342 template<typename _UIntType, size_t __m, 343 size_t __pos1, size_t __sl1, size_t __sl2, 344 size_t __sr1, size_t __sr2, 345 uint32_t __msk1, uint32_t __msk2, 346 uint32_t __msk3, uint32_t __msk4, 347 uint32_t __parity1, uint32_t __parity2, 348 uint32_t __parity3, uint32_t __parity4, 349 typename _CharT, typename _Traits> 350 std::basic_ostream<_CharT, _Traits>& 351 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 352 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 353 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 354 __msk1, __msk2, __msk3, __msk4, 355 __parity1, __parity2, __parity3, __parity4>& __x) 356 { 357 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 358 typedef typename __ostream_type::ios_base __ios_base; 359 360 const typename __ios_base::fmtflags __flags = __os.flags(); 361 const _CharT __fill = __os.fill(); 362 const _CharT __space = __os.widen(' '); 363 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 364 __os.fill(__space); 365 366 for (size_t __i = 0; __i < __x._M_nstate32; ++__i) 367 __os << __x._M_state32[__i] << __space; 368 __os << __x._M_pos; 369 370 __os.flags(__flags); 371 __os.fill(__fill); 372 return __os; 373 } 374 375 376 template<typename _UIntType, size_t __m, 377 size_t __pos1, size_t __sl1, size_t __sl2, 378 size_t __sr1, size_t __sr2, 379 uint32_t __msk1, uint32_t __msk2, 380 uint32_t __msk3, uint32_t __msk4, 381 uint32_t __parity1, uint32_t __parity2, 382 uint32_t __parity3, uint32_t __parity4, 383 typename _CharT, typename _Traits> 384 std::basic_istream<_CharT, _Traits>& 385 operator>>(std::basic_istream<_CharT, _Traits>& __is, 386 __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 387 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 388 __msk1, __msk2, __msk3, __msk4, 389 __parity1, __parity2, __parity3, __parity4>& __x) 390 { 391 typedef std::basic_istream<_CharT, _Traits> __istream_type; 392 typedef typename __istream_type::ios_base __ios_base; 393 394 const typename __ios_base::fmtflags __flags = __is.flags(); 395 __is.flags(__ios_base::dec | __ios_base::skipws); 396 397 for (size_t __i = 0; __i < __x._M_nstate32; ++__i) 398 __is >> __x._M_state32[__i]; 399 __is >> __x._M_pos; 400 401 __is.flags(__flags); 402 return __is; 403 } 404 405#endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ 406 407 /** 408 * Iteration method due to M.D. J<o:>hnk. 409 * 410 * M.D. J<o:>hnk, Erzeugung von betaverteilten und gammaverteilten 411 * Zufallszahlen, Metrika, Volume 8, 1964 412 */ 413 template<typename _RealType> 414 template<typename _UniformRandomNumberGenerator> 415 typename beta_distribution<_RealType>::result_type 416 beta_distribution<_RealType>:: 417 operator()(_UniformRandomNumberGenerator& __urng, 418 const param_type& __param) 419 { 420 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 421 __aurng(__urng); 422 423 result_type __x, __y; 424 do 425 { 426 __x = std::exp(std::log(__aurng()) / __param.alpha()); 427 __y = std::exp(std::log(__aurng()) / __param.beta()); 428 } 429 while (__x + __y > result_type(1)); 430 431 return __x / (__x + __y); 432 } 433 434 template<typename _RealType> 435 template<typename _OutputIterator, 436 typename _UniformRandomNumberGenerator> 437 void 438 beta_distribution<_RealType>:: 439 __generate_impl(_OutputIterator __f, _OutputIterator __t, 440 _UniformRandomNumberGenerator& __urng, 441 const param_type& __param) 442 { 443 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 444 445 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 446 __aurng(__urng); 447 448 while (__f != __t) 449 { 450 result_type __x, __y; 451 do 452 { 453 __x = std::exp(std::log(__aurng()) / __param.alpha()); 454 __y = std::exp(std::log(__aurng()) / __param.beta()); 455 } 456 while (__x + __y > result_type(1)); 457 458 *__f++ = __x / (__x + __y); 459 } 460 } 461 462 template<typename _RealType, typename _CharT, typename _Traits> 463 std::basic_ostream<_CharT, _Traits>& 464 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 465 const __gnu_cxx::beta_distribution<_RealType>& __x) 466 { 467 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 468 typedef typename __ostream_type::ios_base __ios_base; 469 470 const typename __ios_base::fmtflags __flags = __os.flags(); 471 const _CharT __fill = __os.fill(); 472 const std::streamsize __precision = __os.precision(); 473 const _CharT __space = __os.widen(' '); 474 __os.flags(__ios_base::scientific | __ios_base::left); 475 __os.fill(__space); 476 __os.precision(std::numeric_limits<_RealType>::max_digits10); 477 478 __os << __x.alpha() << __space << __x.beta(); 479 480 __os.flags(__flags); 481 __os.fill(__fill); 482 __os.precision(__precision); 483 return __os; 484 } 485 486 template<typename _RealType, typename _CharT, typename _Traits> 487 std::basic_istream<_CharT, _Traits>& 488 operator>>(std::basic_istream<_CharT, _Traits>& __is, 489 __gnu_cxx::beta_distribution<_RealType>& __x) 490 { 491 typedef std::basic_istream<_CharT, _Traits> __istream_type; 492 typedef typename __istream_type::ios_base __ios_base; 493 494 const typename __ios_base::fmtflags __flags = __is.flags(); 495 __is.flags(__ios_base::dec | __ios_base::skipws); 496 497 _RealType __alpha_val, __beta_val; 498 __is >> __alpha_val >> __beta_val; 499 __x.param(typename __gnu_cxx::beta_distribution<_RealType>:: 500 param_type(__alpha_val, __beta_val)); 501 502 __is.flags(__flags); 503 return __is; 504 } 505 506 507 template<std::size_t _Dimen, typename _RealType> 508 template<typename _InputIterator1, typename _InputIterator2> 509 void 510 normal_mv_distribution<_Dimen, _RealType>::param_type:: 511 _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 512 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend) 513 { 514 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 515 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 516 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 517 _M_mean.end(), _RealType(0)); 518 519 // Perform the Cholesky decomposition 520 auto __w = _M_t.begin(); 521 for (size_t __j = 0; __j < _Dimen; ++__j) 522 { 523 _RealType __sum = _RealType(0); 524 525 auto __slitbegin = __w; 526 auto __cit = _M_t.begin(); 527 for (size_t __i = 0; __i < __j; ++__i) 528 { 529 auto __slit = __slitbegin; 530 _RealType __s = *__varcovbegin++; 531 for (size_t __k = 0; __k < __i; ++__k) 532 __s -= *__slit++ * *__cit++; 533 534 *__w++ = __s /= *__cit++; 535 __sum += __s * __s; 536 } 537 538 __sum = *__varcovbegin - __sum; 539 if (__builtin_expect(__sum <= _RealType(0), 0)) 540 std::__throw_runtime_error(__N("normal_mv_distribution::" 541 "param_type::_M_init_full")); 542 *__w++ = std::sqrt(__sum); 543 544 std::advance(__varcovbegin, _Dimen - __j); 545 } 546 } 547 548 template<std::size_t _Dimen, typename _RealType> 549 template<typename _InputIterator1, typename _InputIterator2> 550 void 551 normal_mv_distribution<_Dimen, _RealType>::param_type:: 552 _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 553 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend) 554 { 555 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 556 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 557 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 558 _M_mean.end(), _RealType(0)); 559 560 // Perform the Cholesky decomposition 561 auto __w = _M_t.begin(); 562 for (size_t __j = 0; __j < _Dimen; ++__j) 563 { 564 _RealType __sum = _RealType(0); 565 566 auto __slitbegin = __w; 567 auto __cit = _M_t.begin(); 568 for (size_t __i = 0; __i < __j; ++__i) 569 { 570 auto __slit = __slitbegin; 571 _RealType __s = *__varcovbegin++; 572 for (size_t __k = 0; __k < __i; ++__k) 573 __s -= *__slit++ * *__cit++; 574 575 *__w++ = __s /= *__cit++; 576 __sum += __s * __s; 577 } 578 579 __sum = *__varcovbegin++ - __sum; 580 if (__builtin_expect(__sum <= _RealType(0), 0)) 581 std::__throw_runtime_error(__N("normal_mv_distribution::" 582 "param_type::_M_init_full")); 583 *__w++ = std::sqrt(__sum); 584 } 585 } 586 587 template<std::size_t _Dimen, typename _RealType> 588 template<typename _InputIterator1, typename _InputIterator2> 589 void 590 normal_mv_distribution<_Dimen, _RealType>::param_type:: 591 _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 592 _InputIterator2 __varbegin, _InputIterator2 __varend) 593 { 594 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 595 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 596 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 597 _M_mean.end(), _RealType(0)); 598 599 auto __w = _M_t.begin(); 600 size_t __step = 0; 601 while (__varbegin != __varend) 602 { 603 std::fill_n(__w, __step, _RealType(0)); 604 __w += __step++; 605 if (__builtin_expect(*__varbegin < _RealType(0), 0)) 606 std::__throw_runtime_error(__N("normal_mv_distribution::" 607 "param_type::_M_init_diagonal")); 608 *__w++ = std::sqrt(*__varbegin++); 609 } 610 } 611 612 template<std::size_t _Dimen, typename _RealType> 613 template<typename _UniformRandomNumberGenerator> 614 typename normal_mv_distribution<_Dimen, _RealType>::result_type 615 normal_mv_distribution<_Dimen, _RealType>:: 616 operator()(_UniformRandomNumberGenerator& __urng, 617 const param_type& __param) 618 { 619 result_type __ret; 620 621 _M_nd.__generate(__ret.begin(), __ret.end(), __urng); 622 623 auto __t_it = __param._M_t.crbegin(); 624 for (size_t __i = _Dimen; __i > 0; --__i) 625 { 626 _RealType __sum = _RealType(0); 627 for (size_t __j = __i; __j > 0; --__j) 628 __sum += __ret[__j - 1] * *__t_it++; 629 __ret[__i - 1] = __sum; 630 } 631 632 return __ret; 633 } 634 635 template<std::size_t _Dimen, typename _RealType> 636 template<typename _ForwardIterator, typename _UniformRandomNumberGenerator> 637 void 638 normal_mv_distribution<_Dimen, _RealType>:: 639 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 640 _UniformRandomNumberGenerator& __urng, 641 const param_type& __param) 642 { 643 __glibcxx_function_requires(_Mutable_ForwardIteratorConcept< 644 _ForwardIterator>) 645 while (__f != __t) 646 *__f++ = this->operator()(__urng, __param); 647 } 648 649 template<size_t _Dimen, typename _RealType> 650 bool 651 operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 652 __d1, 653 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 654 __d2) 655 { 656 return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd; 657 } 658 659 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits> 660 std::basic_ostream<_CharT, _Traits>& 661 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 662 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x) 663 { 664 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 665 typedef typename __ostream_type::ios_base __ios_base; 666 667 const typename __ios_base::fmtflags __flags = __os.flags(); 668 const _CharT __fill = __os.fill(); 669 const std::streamsize __precision = __os.precision(); 670 const _CharT __space = __os.widen(' '); 671 __os.flags(__ios_base::scientific | __ios_base::left); 672 __os.fill(__space); 673 __os.precision(std::numeric_limits<_RealType>::max_digits10); 674 675 auto __mean = __x._M_param.mean(); 676 for (auto __it : __mean) 677 __os << __it << __space; 678 auto __t = __x._M_param.varcov(); 679 for (auto __it : __t) 680 __os << __it << __space; 681 682 __os << __x._M_nd; 683 684 __os.flags(__flags); 685 __os.fill(__fill); 686 __os.precision(__precision); 687 return __os; 688 } 689 690 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits> 691 std::basic_istream<_CharT, _Traits>& 692 operator>>(std::basic_istream<_CharT, _Traits>& __is, 693 __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x) 694 { 695 typedef std::basic_istream<_CharT, _Traits> __istream_type; 696 typedef typename __istream_type::ios_base __ios_base; 697 698 const typename __ios_base::fmtflags __flags = __is.flags(); 699 __is.flags(__ios_base::dec | __ios_base::skipws); 700 701 std::array<_RealType, _Dimen> __mean; 702 for (auto& __it : __mean) 703 __is >> __it; 704 std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov; 705 for (auto& __it : __varcov) 706 __is >> __it; 707 708 __is >> __x._M_nd; 709 710 __x.param(typename normal_mv_distribution<_Dimen, _RealType>:: 711 param_type(__mean.begin(), __mean.end(), 712 __varcov.begin(), __varcov.end())); 713 714 __is.flags(__flags); 715 return __is; 716 } 717 718 719 template<typename _RealType> 720 template<typename _OutputIterator, 721 typename _UniformRandomNumberGenerator> 722 void 723 rice_distribution<_RealType>:: 724 __generate_impl(_OutputIterator __f, _OutputIterator __t, 725 _UniformRandomNumberGenerator& __urng, 726 const param_type& __p) 727 { 728 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 729 730 while (__f != __t) 731 { 732 typename std::normal_distribution<result_type>::param_type 733 __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma()); 734 result_type __x = this->_M_ndx(__px, __urng); 735 result_type __y = this->_M_ndy(__py, __urng); 736#if _GLIBCXX_USE_C99_MATH_TR1 737 *__f++ = std::hypot(__x, __y); 738#else 739 *__f++ = std::sqrt(__x * __x + __y * __y); 740#endif 741 } 742 } 743 744 template<typename _RealType, typename _CharT, typename _Traits> 745 std::basic_ostream<_CharT, _Traits>& 746 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 747 const rice_distribution<_RealType>& __x) 748 { 749 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 750 typedef typename __ostream_type::ios_base __ios_base; 751 752 const typename __ios_base::fmtflags __flags = __os.flags(); 753 const _CharT __fill = __os.fill(); 754 const std::streamsize __precision = __os.precision(); 755 const _CharT __space = __os.widen(' '); 756 __os.flags(__ios_base::scientific | __ios_base::left); 757 __os.fill(__space); 758 __os.precision(std::numeric_limits<_RealType>::max_digits10); 759 760 __os << __x.nu() << __space << __x.sigma(); 761 __os << __space << __x._M_ndx; 762 __os << __space << __x._M_ndy; 763 764 __os.flags(__flags); 765 __os.fill(__fill); 766 __os.precision(__precision); 767 return __os; 768 } 769 770 template<typename _RealType, typename _CharT, typename _Traits> 771 std::basic_istream<_CharT, _Traits>& 772 operator>>(std::basic_istream<_CharT, _Traits>& __is, 773 rice_distribution<_RealType>& __x) 774 { 775 typedef std::basic_istream<_CharT, _Traits> __istream_type; 776 typedef typename __istream_type::ios_base __ios_base; 777 778 const typename __ios_base::fmtflags __flags = __is.flags(); 779 __is.flags(__ios_base::dec | __ios_base::skipws); 780 781 _RealType __nu_val, __sigma_val; 782 __is >> __nu_val >> __sigma_val; 783 __is >> __x._M_ndx; 784 __is >> __x._M_ndy; 785 __x.param(typename rice_distribution<_RealType>:: 786 param_type(__nu_val, __sigma_val)); 787 788 __is.flags(__flags); 789 return __is; 790 } 791 792 793 template<typename _RealType> 794 template<typename _OutputIterator, 795 typename _UniformRandomNumberGenerator> 796 void 797 nakagami_distribution<_RealType>:: 798 __generate_impl(_OutputIterator __f, _OutputIterator __t, 799 _UniformRandomNumberGenerator& __urng, 800 const param_type& __p) 801 { 802 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 803 804 typename std::gamma_distribution<result_type>::param_type 805 __pg(__p.mu(), __p.omega() / __p.mu()); 806 while (__f != __t) 807 *__f++ = std::sqrt(this->_M_gd(__pg, __urng)); 808 } 809 810 template<typename _RealType, typename _CharT, typename _Traits> 811 std::basic_ostream<_CharT, _Traits>& 812 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 813 const nakagami_distribution<_RealType>& __x) 814 { 815 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 816 typedef typename __ostream_type::ios_base __ios_base; 817 818 const typename __ios_base::fmtflags __flags = __os.flags(); 819 const _CharT __fill = __os.fill(); 820 const std::streamsize __precision = __os.precision(); 821 const _CharT __space = __os.widen(' '); 822 __os.flags(__ios_base::scientific | __ios_base::left); 823 __os.fill(__space); 824 __os.precision(std::numeric_limits<_RealType>::max_digits10); 825 826 __os << __x.mu() << __space << __x.omega(); 827 __os << __space << __x._M_gd; 828 829 __os.flags(__flags); 830 __os.fill(__fill); 831 __os.precision(__precision); 832 return __os; 833 } 834 835 template<typename _RealType, typename _CharT, typename _Traits> 836 std::basic_istream<_CharT, _Traits>& 837 operator>>(std::basic_istream<_CharT, _Traits>& __is, 838 nakagami_distribution<_RealType>& __x) 839 { 840 typedef std::basic_istream<_CharT, _Traits> __istream_type; 841 typedef typename __istream_type::ios_base __ios_base; 842 843 const typename __ios_base::fmtflags __flags = __is.flags(); 844 __is.flags(__ios_base::dec | __ios_base::skipws); 845 846 _RealType __mu_val, __omega_val; 847 __is >> __mu_val >> __omega_val; 848 __is >> __x._M_gd; 849 __x.param(typename nakagami_distribution<_RealType>:: 850 param_type(__mu_val, __omega_val)); 851 852 __is.flags(__flags); 853 return __is; 854 } 855 856 857 template<typename _RealType> 858 template<typename _OutputIterator, 859 typename _UniformRandomNumberGenerator> 860 void 861 pareto_distribution<_RealType>:: 862 __generate_impl(_OutputIterator __f, _OutputIterator __t, 863 _UniformRandomNumberGenerator& __urng, 864 const param_type& __p) 865 { 866 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 867 868 result_type __mu_val = __p.mu(); 869 result_type __malphinv = -result_type(1) / __p.alpha(); 870 while (__f != __t) 871 *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv); 872 } 873 874 template<typename _RealType, typename _CharT, typename _Traits> 875 std::basic_ostream<_CharT, _Traits>& 876 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 877 const pareto_distribution<_RealType>& __x) 878 { 879 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 880 typedef typename __ostream_type::ios_base __ios_base; 881 882 const typename __ios_base::fmtflags __flags = __os.flags(); 883 const _CharT __fill = __os.fill(); 884 const std::streamsize __precision = __os.precision(); 885 const _CharT __space = __os.widen(' '); 886 __os.flags(__ios_base::scientific | __ios_base::left); 887 __os.fill(__space); 888 __os.precision(std::numeric_limits<_RealType>::max_digits10); 889 890 __os << __x.alpha() << __space << __x.mu(); 891 __os << __space << __x._M_ud; 892 893 __os.flags(__flags); 894 __os.fill(__fill); 895 __os.precision(__precision); 896 return __os; 897 } 898 899 template<typename _RealType, typename _CharT, typename _Traits> 900 std::basic_istream<_CharT, _Traits>& 901 operator>>(std::basic_istream<_CharT, _Traits>& __is, 902 pareto_distribution<_RealType>& __x) 903 { 904 typedef std::basic_istream<_CharT, _Traits> __istream_type; 905 typedef typename __istream_type::ios_base __ios_base; 906 907 const typename __ios_base::fmtflags __flags = __is.flags(); 908 __is.flags(__ios_base::dec | __ios_base::skipws); 909 910 _RealType __alpha_val, __mu_val; 911 __is >> __alpha_val >> __mu_val; 912 __is >> __x._M_ud; 913 __x.param(typename pareto_distribution<_RealType>:: 914 param_type(__alpha_val, __mu_val)); 915 916 __is.flags(__flags); 917 return __is; 918 } 919 920 921 template<typename _RealType> 922 template<typename _UniformRandomNumberGenerator> 923 typename k_distribution<_RealType>::result_type 924 k_distribution<_RealType>:: 925 operator()(_UniformRandomNumberGenerator& __urng) 926 { 927 result_type __x = this->_M_gd1(__urng); 928 result_type __y = this->_M_gd2(__urng); 929 return std::sqrt(__x * __y); 930 } 931 932 template<typename _RealType> 933 template<typename _UniformRandomNumberGenerator> 934 typename k_distribution<_RealType>::result_type 935 k_distribution<_RealType>:: 936 operator()(_UniformRandomNumberGenerator& __urng, 937 const param_type& __p) 938 { 939 typename std::gamma_distribution<result_type>::param_type 940 __p1(__p.lambda(), result_type(1) / __p.lambda()), 941 __p2(__p.nu(), __p.mu() / __p.nu()); 942 result_type __x = this->_M_gd1(__p1, __urng); 943 result_type __y = this->_M_gd2(__p2, __urng); 944 return std::sqrt(__x * __y); 945 } 946 947 template<typename _RealType> 948 template<typename _OutputIterator, 949 typename _UniformRandomNumberGenerator> 950 void 951 k_distribution<_RealType>:: 952 __generate_impl(_OutputIterator __f, _OutputIterator __t, 953 _UniformRandomNumberGenerator& __urng, 954 const param_type& __p) 955 { 956 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 957 958 typename std::gamma_distribution<result_type>::param_type 959 __p1(__p.lambda(), result_type(1) / __p.lambda()), 960 __p2(__p.nu(), __p.mu() / __p.nu()); 961 while (__f != __t) 962 { 963 result_type __x = this->_M_gd1(__p1, __urng); 964 result_type __y = this->_M_gd2(__p2, __urng); 965 *__f++ = std::sqrt(__x * __y); 966 } 967 } 968 969 template<typename _RealType, typename _CharT, typename _Traits> 970 std::basic_ostream<_CharT, _Traits>& 971 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 972 const k_distribution<_RealType>& __x) 973 { 974 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 975 typedef typename __ostream_type::ios_base __ios_base; 976 977 const typename __ios_base::fmtflags __flags = __os.flags(); 978 const _CharT __fill = __os.fill(); 979 const std::streamsize __precision = __os.precision(); 980 const _CharT __space = __os.widen(' '); 981 __os.flags(__ios_base::scientific | __ios_base::left); 982 __os.fill(__space); 983 __os.precision(std::numeric_limits<_RealType>::max_digits10); 984 985 __os << __x.lambda() << __space << __x.mu() << __space << __x.nu(); 986 __os << __space << __x._M_gd1; 987 __os << __space << __x._M_gd2; 988 989 __os.flags(__flags); 990 __os.fill(__fill); 991 __os.precision(__precision); 992 return __os; 993 } 994 995 template<typename _RealType, typename _CharT, typename _Traits> 996 std::basic_istream<_CharT, _Traits>& 997 operator>>(std::basic_istream<_CharT, _Traits>& __is, 998 k_distribution<_RealType>& __x) 999 { 1000 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1001 typedef typename __istream_type::ios_base __ios_base; 1002 1003 const typename __ios_base::fmtflags __flags = __is.flags(); 1004 __is.flags(__ios_base::dec | __ios_base::skipws); 1005 1006 _RealType __lambda_val, __mu_val, __nu_val; 1007 __is >> __lambda_val >> __mu_val >> __nu_val; 1008 __is >> __x._M_gd1; 1009 __is >> __x._M_gd2; 1010 __x.param(typename k_distribution<_RealType>:: 1011 param_type(__lambda_val, __mu_val, __nu_val)); 1012 1013 __is.flags(__flags); 1014 return __is; 1015 } 1016 1017 1018 template<typename _RealType> 1019 template<typename _OutputIterator, 1020 typename _UniformRandomNumberGenerator> 1021 void 1022 arcsine_distribution<_RealType>:: 1023 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1024 _UniformRandomNumberGenerator& __urng, 1025 const param_type& __p) 1026 { 1027 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 1028 1029 result_type __dif = __p.b() - __p.a(); 1030 result_type __sum = __p.a() + __p.b(); 1031 while (__f != __t) 1032 { 1033 result_type __x = std::sin(this->_M_ud(__urng)); 1034 *__f++ = (__x * __dif + __sum) / result_type(2); 1035 } 1036 } 1037 1038 template<typename _RealType, typename _CharT, typename _Traits> 1039 std::basic_ostream<_CharT, _Traits>& 1040 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1041 const arcsine_distribution<_RealType>& __x) 1042 { 1043 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1044 typedef typename __ostream_type::ios_base __ios_base; 1045 1046 const typename __ios_base::fmtflags __flags = __os.flags(); 1047 const _CharT __fill = __os.fill(); 1048 const std::streamsize __precision = __os.precision(); 1049 const _CharT __space = __os.widen(' '); 1050 __os.flags(__ios_base::scientific | __ios_base::left); 1051 __os.fill(__space); 1052 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1053 1054 __os << __x.a() << __space << __x.b(); 1055 __os << __space << __x._M_ud; 1056 1057 __os.flags(__flags); 1058 __os.fill(__fill); 1059 __os.precision(__precision); 1060 return __os; 1061 } 1062 1063 template<typename _RealType, typename _CharT, typename _Traits> 1064 std::basic_istream<_CharT, _Traits>& 1065 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1066 arcsine_distribution<_RealType>& __x) 1067 { 1068 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1069 typedef typename __istream_type::ios_base __ios_base; 1070 1071 const typename __ios_base::fmtflags __flags = __is.flags(); 1072 __is.flags(__ios_base::dec | __ios_base::skipws); 1073 1074 _RealType __a, __b; 1075 __is >> __a >> __b; 1076 __is >> __x._M_ud; 1077 __x.param(typename arcsine_distribution<_RealType>:: 1078 param_type(__a, __b)); 1079 1080 __is.flags(__flags); 1081 return __is; 1082 } 1083 1084 1085 template<typename _RealType> 1086 template<typename _UniformRandomNumberGenerator> 1087 typename hoyt_distribution<_RealType>::result_type 1088 hoyt_distribution<_RealType>:: 1089 operator()(_UniformRandomNumberGenerator& __urng) 1090 { 1091 result_type __x = this->_M_ad(__urng); 1092 result_type __y = this->_M_ed(__urng); 1093 return (result_type(2) * this->q() 1094 / (result_type(1) + this->q() * this->q())) 1095 * std::sqrt(this->omega() * __x * __y); 1096 } 1097 1098 template<typename _RealType> 1099 template<typename _UniformRandomNumberGenerator> 1100 typename hoyt_distribution<_RealType>::result_type 1101 hoyt_distribution<_RealType>:: 1102 operator()(_UniformRandomNumberGenerator& __urng, 1103 const param_type& __p) 1104 { 1105 result_type __q2 = __p.q() * __p.q(); 1106 result_type __num = result_type(0.5L) * (result_type(1) + __q2); 1107 typename __gnu_cxx::arcsine_distribution<result_type>::param_type 1108 __pa(__num, __num / __q2); 1109 result_type __x = this->_M_ad(__pa, __urng); 1110 result_type __y = this->_M_ed(__urng); 1111 return (result_type(2) * __p.q() / (result_type(1) + __q2)) 1112 * std::sqrt(__p.omega() * __x * __y); 1113 } 1114 1115 template<typename _RealType> 1116 template<typename _OutputIterator, 1117 typename _UniformRandomNumberGenerator> 1118 void 1119 hoyt_distribution<_RealType>:: 1120 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1121 _UniformRandomNumberGenerator& __urng, 1122 const param_type& __p) 1123 { 1124 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 1125 1126 result_type __2q = result_type(2) * __p.q(); 1127 result_type __q2 = __p.q() * __p.q(); 1128 result_type __q2p1 = result_type(1) + __q2; 1129 result_type __num = result_type(0.5L) * __q2p1; 1130 result_type __omega = __p.omega(); 1131 typename __gnu_cxx::arcsine_distribution<result_type>::param_type 1132 __pa(__num, __num / __q2); 1133 while (__f != __t) 1134 { 1135 result_type __x = this->_M_ad(__pa, __urng); 1136 result_type __y = this->_M_ed(__urng); 1137 *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y); 1138 } 1139 } 1140 1141 template<typename _RealType, typename _CharT, typename _Traits> 1142 std::basic_ostream<_CharT, _Traits>& 1143 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1144 const hoyt_distribution<_RealType>& __x) 1145 { 1146 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1147 typedef typename __ostream_type::ios_base __ios_base; 1148 1149 const typename __ios_base::fmtflags __flags = __os.flags(); 1150 const _CharT __fill = __os.fill(); 1151 const std::streamsize __precision = __os.precision(); 1152 const _CharT __space = __os.widen(' '); 1153 __os.flags(__ios_base::scientific | __ios_base::left); 1154 __os.fill(__space); 1155 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1156 1157 __os << __x.q() << __space << __x.omega(); 1158 __os << __space << __x._M_ad; 1159 __os << __space << __x._M_ed; 1160 1161 __os.flags(__flags); 1162 __os.fill(__fill); 1163 __os.precision(__precision); 1164 return __os; 1165 } 1166 1167 template<typename _RealType, typename _CharT, typename _Traits> 1168 std::basic_istream<_CharT, _Traits>& 1169 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1170 hoyt_distribution<_RealType>& __x) 1171 { 1172 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1173 typedef typename __istream_type::ios_base __ios_base; 1174 1175 const typename __ios_base::fmtflags __flags = __is.flags(); 1176 __is.flags(__ios_base::dec | __ios_base::skipws); 1177 1178 _RealType __q, __omega; 1179 __is >> __q >> __omega; 1180 __is >> __x._M_ad; 1181 __is >> __x._M_ed; 1182 __x.param(typename hoyt_distribution<_RealType>:: 1183 param_type(__q, __omega)); 1184 1185 __is.flags(__flags); 1186 return __is; 1187 } 1188 1189 1190 template<typename _RealType> 1191 template<typename _OutputIterator, 1192 typename _UniformRandomNumberGenerator> 1193 void 1194 triangular_distribution<_RealType>:: 1195 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1196 _UniformRandomNumberGenerator& __urng, 1197 const param_type& __param) 1198 { 1199 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 1200 1201 while (__f != __t) 1202 *__f++ = this->operator()(__urng, __param); 1203 } 1204 1205 template<typename _RealType, typename _CharT, typename _Traits> 1206 std::basic_ostream<_CharT, _Traits>& 1207 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1208 const __gnu_cxx::triangular_distribution<_RealType>& __x) 1209 { 1210 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1211 typedef typename __ostream_type::ios_base __ios_base; 1212 1213 const typename __ios_base::fmtflags __flags = __os.flags(); 1214 const _CharT __fill = __os.fill(); 1215 const std::streamsize __precision = __os.precision(); 1216 const _CharT __space = __os.widen(' '); 1217 __os.flags(__ios_base::scientific | __ios_base::left); 1218 __os.fill(__space); 1219 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1220 1221 __os << __x.a() << __space << __x.b() << __space << __x.c(); 1222 1223 __os.flags(__flags); 1224 __os.fill(__fill); 1225 __os.precision(__precision); 1226 return __os; 1227 } 1228 1229 template<typename _RealType, typename _CharT, typename _Traits> 1230 std::basic_istream<_CharT, _Traits>& 1231 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1232 __gnu_cxx::triangular_distribution<_RealType>& __x) 1233 { 1234 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1235 typedef typename __istream_type::ios_base __ios_base; 1236 1237 const typename __ios_base::fmtflags __flags = __is.flags(); 1238 __is.flags(__ios_base::dec | __ios_base::skipws); 1239 1240 _RealType __a, __b, __c; 1241 __is >> __a >> __b >> __c; 1242 __x.param(typename __gnu_cxx::triangular_distribution<_RealType>:: 1243 param_type(__a, __b, __c)); 1244 1245 __is.flags(__flags); 1246 return __is; 1247 } 1248 1249 1250 template<typename _RealType> 1251 template<typename _UniformRandomNumberGenerator> 1252 typename von_mises_distribution<_RealType>::result_type 1253 von_mises_distribution<_RealType>:: 1254 operator()(_UniformRandomNumberGenerator& __urng, 1255 const param_type& __p) 1256 { 1257 const result_type __pi 1258 = __gnu_cxx::__math_constants<result_type>::__pi; 1259 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1260 __aurng(__urng); 1261 1262 result_type __f; 1263 while (1) 1264 { 1265 result_type __rnd = std::cos(__pi * __aurng()); 1266 __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd); 1267 result_type __c = __p._M_kappa * (__p._M_r - __f); 1268 1269 result_type __rnd2 = __aurng(); 1270 if (__c * (result_type(2) - __c) > __rnd2) 1271 break; 1272 if (std::log(__c / __rnd2) >= __c - result_type(1)) 1273 break; 1274 } 1275 1276 result_type __res = std::acos(__f); 1277#if _GLIBCXX_USE_C99_MATH_TR1 1278 __res = std::copysign(__res, __aurng() - result_type(0.5)); 1279#else 1280 if (__aurng() < result_type(0.5)) 1281 __res = -__res; 1282#endif 1283 __res += __p._M_mu; 1284 if (__res > __pi) 1285 __res -= result_type(2) * __pi; 1286 else if (__res < -__pi) 1287 __res += result_type(2) * __pi; 1288 return __res; 1289 } 1290 1291 template<typename _RealType> 1292 template<typename _OutputIterator, 1293 typename _UniformRandomNumberGenerator> 1294 void 1295 von_mises_distribution<_RealType>:: 1296 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1297 _UniformRandomNumberGenerator& __urng, 1298 const param_type& __param) 1299 { 1300 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 1301 1302 while (__f != __t) 1303 *__f++ = this->operator()(__urng, __param); 1304 } 1305 1306 template<typename _RealType, typename _CharT, typename _Traits> 1307 std::basic_ostream<_CharT, _Traits>& 1308 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1309 const __gnu_cxx::von_mises_distribution<_RealType>& __x) 1310 { 1311 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1312 typedef typename __ostream_type::ios_base __ios_base; 1313 1314 const typename __ios_base::fmtflags __flags = __os.flags(); 1315 const _CharT __fill = __os.fill(); 1316 const std::streamsize __precision = __os.precision(); 1317 const _CharT __space = __os.widen(' '); 1318 __os.flags(__ios_base::scientific | __ios_base::left); 1319 __os.fill(__space); 1320 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1321 1322 __os << __x.mu() << __space << __x.kappa(); 1323 1324 __os.flags(__flags); 1325 __os.fill(__fill); 1326 __os.precision(__precision); 1327 return __os; 1328 } 1329 1330 template<typename _RealType, typename _CharT, typename _Traits> 1331 std::basic_istream<_CharT, _Traits>& 1332 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1333 __gnu_cxx::von_mises_distribution<_RealType>& __x) 1334 { 1335 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1336 typedef typename __istream_type::ios_base __ios_base; 1337 1338 const typename __ios_base::fmtflags __flags = __is.flags(); 1339 __is.flags(__ios_base::dec | __ios_base::skipws); 1340 1341 _RealType __mu, __kappa; 1342 __is >> __mu >> __kappa; 1343 __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>:: 1344 param_type(__mu, __kappa)); 1345 1346 __is.flags(__flags); 1347 return __is; 1348 } 1349 1350 1351 template<typename _UIntType> 1352 template<typename _UniformRandomNumberGenerator> 1353 typename hypergeometric_distribution<_UIntType>::result_type 1354 hypergeometric_distribution<_UIntType>:: 1355 operator()(_UniformRandomNumberGenerator& __urng, 1356 const param_type& __param) 1357 { 1358 std::__detail::_Adaptor<_UniformRandomNumberGenerator, double> 1359 __aurng(__urng); 1360 1361 result_type __a = __param.successful_size(); 1362 result_type __b = __param.total_size(); 1363 result_type __k = 0; 1364 1365 if (__param.total_draws() < __param.total_size() / 2) 1366 { 1367 for (result_type __i = 0; __i < __param.total_draws(); ++__i) 1368 { 1369 if (__b * __aurng() < __a) 1370 { 1371 ++__k; 1372 if (__k == __param.successful_size()) 1373 return __k; 1374 --__a; 1375 } 1376 --__b; 1377 } 1378 return __k; 1379 } 1380 else 1381 { 1382 for (result_type __i = 0; __i < __param.unsuccessful_size(); ++__i) 1383 { 1384 if (__b * __aurng() < __a) 1385 { 1386 ++__k; 1387 if (__k == __param.successful_size()) 1388 return __param.successful_size() - __k; 1389 --__a; 1390 } 1391 --__b; 1392 } 1393 return __param.successful_size() - __k; 1394 } 1395 } 1396 1397 template<typename _UIntType> 1398 template<typename _OutputIterator, 1399 typename _UniformRandomNumberGenerator> 1400 void 1401 hypergeometric_distribution<_UIntType>:: 1402 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1403 _UniformRandomNumberGenerator& __urng, 1404 const param_type& __param) 1405 { 1406 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 1407 1408 while (__f != __t) 1409 *__f++ = this->operator()(__urng); 1410 } 1411 1412 template<typename _UIntType, typename _CharT, typename _Traits> 1413 std::basic_ostream<_CharT, _Traits>& 1414 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1415 const __gnu_cxx::hypergeometric_distribution<_UIntType>& __x) 1416 { 1417 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1418 typedef typename __ostream_type::ios_base __ios_base; 1419 1420 const typename __ios_base::fmtflags __flags = __os.flags(); 1421 const _CharT __fill = __os.fill(); 1422 const std::streamsize __precision = __os.precision(); 1423 const _CharT __space = __os.widen(' '); 1424 __os.flags(__ios_base::scientific | __ios_base::left); 1425 __os.fill(__space); 1426 __os.precision(std::numeric_limits<_UIntType>::max_digits10); 1427 1428 __os << __x.total_size() << __space << __x.successful_size() << __space 1429 << __x.total_draws(); 1430 1431 __os.flags(__flags); 1432 __os.fill(__fill); 1433 __os.precision(__precision); 1434 return __os; 1435 } 1436 1437 template<typename _UIntType, typename _CharT, typename _Traits> 1438 std::basic_istream<_CharT, _Traits>& 1439 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1440 __gnu_cxx::hypergeometric_distribution<_UIntType>& __x) 1441 { 1442 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1443 typedef typename __istream_type::ios_base __ios_base; 1444 1445 const typename __ios_base::fmtflags __flags = __is.flags(); 1446 __is.flags(__ios_base::dec | __ios_base::skipws); 1447 1448 _UIntType __total_size, __successful_size, __total_draws; 1449 __is >> __total_size >> __successful_size >> __total_draws; 1450 __x.param(typename __gnu_cxx::hypergeometric_distribution<_UIntType>:: 1451 param_type(__total_size, __successful_size, __total_draws)); 1452 1453 __is.flags(__flags); 1454 return __is; 1455 } 1456 1457 1458 template<typename _RealType> 1459 template<typename _UniformRandomNumberGenerator> 1460 typename logistic_distribution<_RealType>::result_type 1461 logistic_distribution<_RealType>:: 1462 operator()(_UniformRandomNumberGenerator& __urng, 1463 const param_type& __p) 1464 { 1465 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1466 __aurng(__urng); 1467 1468 result_type __arg = result_type(1); 1469 while (__arg == result_type(1) || __arg == result_type(0)) 1470 __arg = __aurng(); 1471 return __p.a() 1472 + __p.b() * std::log(__arg / (result_type(1) - __arg)); 1473 } 1474 1475 template<typename _RealType> 1476 template<typename _OutputIterator, 1477 typename _UniformRandomNumberGenerator> 1478 void 1479 logistic_distribution<_RealType>:: 1480 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1481 _UniformRandomNumberGenerator& __urng, 1482 const param_type& __p) 1483 { 1484 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 1485 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1486 __aurng(__urng); 1487 1488 while (__f != __t) 1489 { 1490 result_type __arg = result_type(1); 1491 while (__arg == result_type(1) || __arg == result_type(0)) 1492 __arg = __aurng(); 1493 *__f++ = __p.a() 1494 + __p.b() * std::log(__arg / (result_type(1) - __arg)); 1495 } 1496 } 1497 1498 template<typename _RealType, typename _CharT, typename _Traits> 1499 std::basic_ostream<_CharT, _Traits>& 1500 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1501 const logistic_distribution<_RealType>& __x) 1502 { 1503 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1504 typedef typename __ostream_type::ios_base __ios_base; 1505 1506 const typename __ios_base::fmtflags __flags = __os.flags(); 1507 const _CharT __fill = __os.fill(); 1508 const std::streamsize __precision = __os.precision(); 1509 const _CharT __space = __os.widen(' '); 1510 __os.flags(__ios_base::scientific | __ios_base::left); 1511 __os.fill(__space); 1512 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1513 1514 __os << __x.a() << __space << __x.b(); 1515 1516 __os.flags(__flags); 1517 __os.fill(__fill); 1518 __os.precision(__precision); 1519 return __os; 1520 } 1521 1522 template<typename _RealType, typename _CharT, typename _Traits> 1523 std::basic_istream<_CharT, _Traits>& 1524 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1525 logistic_distribution<_RealType>& __x) 1526 { 1527 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1528 typedef typename __istream_type::ios_base __ios_base; 1529 1530 const typename __ios_base::fmtflags __flags = __is.flags(); 1531 __is.flags(__ios_base::dec | __ios_base::skipws); 1532 1533 _RealType __a, __b; 1534 __is >> __a >> __b; 1535 __x.param(typename logistic_distribution<_RealType>:: 1536 param_type(__a, __b)); 1537 1538 __is.flags(__flags); 1539 return __is; 1540 } 1541 1542 1543 namespace { 1544 1545 // Helper class for the uniform_on_sphere_distribution generation 1546 // function. 1547 template<std::size_t _Dimen, typename _RealType> 1548 class uniform_on_sphere_helper 1549 { 1550 typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>:: 1551 result_type result_type; 1552 1553 public: 1554 template<typename _NormalDistribution, 1555 typename _UniformRandomNumberGenerator> 1556 result_type operator()(_NormalDistribution& __nd, 1557 _UniformRandomNumberGenerator& __urng) 1558 { 1559 result_type __ret; 1560 typename result_type::value_type __norm; 1561 1562 do 1563 { 1564 auto __sum = _RealType(0); 1565 1566 std::generate(__ret.begin(), __ret.end(), 1567 [&__nd, &__urng, &__sum](){ 1568 _RealType __t = __nd(__urng); 1569 __sum += __t * __t; 1570 return __t; }); 1571 __norm = std::sqrt(__sum); 1572 } 1573 while (__norm == _RealType(0) || ! __builtin_isfinite(__norm)); 1574 1575 std::transform(__ret.begin(), __ret.end(), __ret.begin(), 1576 [__norm](_RealType __val){ return __val / __norm; }); 1577 1578 return __ret; 1579 } 1580 }; 1581 1582 1583 template<typename _RealType> 1584 class uniform_on_sphere_helper<2, _RealType> 1585 { 1586 typedef typename uniform_on_sphere_distribution<2, _RealType>:: 1587 result_type result_type; 1588 1589 public: 1590 template<typename _NormalDistribution, 1591 typename _UniformRandomNumberGenerator> 1592 result_type operator()(_NormalDistribution&, 1593 _UniformRandomNumberGenerator& __urng) 1594 { 1595 result_type __ret; 1596 _RealType __sq; 1597 std::__detail::_Adaptor<_UniformRandomNumberGenerator, 1598 _RealType> __aurng(__urng); 1599 1600 do 1601 { 1602 __ret[0] = _RealType(2) * __aurng() - _RealType(1); 1603 __ret[1] = _RealType(2) * __aurng() - _RealType(1); 1604 1605 __sq = __ret[0] * __ret[0] + __ret[1] * __ret[1]; 1606 } 1607 while (__sq == _RealType(0) || __sq > _RealType(1)); 1608 1609#if _GLIBCXX_USE_C99_MATH_TR1 1610 // Yes, we do not just use sqrt(__sq) because hypot() is more 1611 // accurate. 1612 auto __norm = std::hypot(__ret[0], __ret[1]); 1613#else 1614 auto __norm = std::sqrt(__sq); 1615#endif 1616 __ret[0] /= __norm; 1617 __ret[1] /= __norm; 1618 1619 return __ret; 1620 } 1621 }; 1622 1623 } 1624 1625 1626 template<std::size_t _Dimen, typename _RealType> 1627 template<typename _UniformRandomNumberGenerator> 1628 typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type 1629 uniform_on_sphere_distribution<_Dimen, _RealType>:: 1630 operator()(_UniformRandomNumberGenerator& __urng, 1631 const param_type& __p) 1632 { 1633 uniform_on_sphere_helper<_Dimen, _RealType> __helper; 1634 return __helper(_M_nd, __urng); 1635 } 1636 1637 template<std::size_t _Dimen, typename _RealType> 1638 template<typename _OutputIterator, 1639 typename _UniformRandomNumberGenerator> 1640 void 1641 uniform_on_sphere_distribution<_Dimen, _RealType>:: 1642 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1643 _UniformRandomNumberGenerator& __urng, 1644 const param_type& __param) 1645 { 1646 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 1647 1648 while (__f != __t) 1649 *__f++ = this->operator()(__urng, __param); 1650 } 1651 1652 template<std::size_t _Dimen, typename _RealType, typename _CharT, 1653 typename _Traits> 1654 std::basic_ostream<_CharT, _Traits>& 1655 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1656 const __gnu_cxx::uniform_on_sphere_distribution<_Dimen, 1657 _RealType>& __x) 1658 { 1659 return __os << __x._M_nd; 1660 } 1661 1662 template<std::size_t _Dimen, typename _RealType, typename _CharT, 1663 typename _Traits> 1664 std::basic_istream<_CharT, _Traits>& 1665 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1666 __gnu_cxx::uniform_on_sphere_distribution<_Dimen, 1667 _RealType>& __x) 1668 { 1669 return __is >> __x._M_nd; 1670 } 1671 1672_GLIBCXX_END_NAMESPACE_VERSION 1673} // namespace 1674 1675 1676#endif // _EXT_RANDOM_TCC 1677