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