1 // Random number extensions -*- C++ -*- 2 3 // Copyright (C) 2012-2018 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 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>:: seed(_Sseq & __q)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>:: _M_period_certification(void)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>:: discard(unsigned long long __z)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> __rshift(uint32_t * __out,const uint32_t * __in)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> __lshift(uint32_t * __out,const uint32_t * __in)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> __recursion(uint32_t * __r,const uint32_t * __a,const uint32_t * __b,const uint32_t * __c,const uint32_t * __d)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>:: _M_gen_rand(void)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 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)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>& 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)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>& 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)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>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)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>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)439 __generate_impl(_OutputIterator __f, _OutputIterator __t, 440 _UniformRandomNumberGenerator& __urng, 441 const param_type& __param) 442 { 443 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 444 result_type>) 445 446 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 447 __aurng(__urng); 448 449 while (__f != __t) 450 { 451 result_type __x, __y; 452 do 453 { 454 __x = std::exp(std::log(__aurng()) / __param.alpha()); 455 __y = std::exp(std::log(__aurng()) / __param.beta()); 456 } 457 while (__x + __y > result_type(1)); 458 459 *__f++ = __x / (__x + __y); 460 } 461 } 462 463 template<typename _RealType, typename _CharT, typename _Traits> 464 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const __gnu_cxx::beta_distribution<_RealType> & __x)465 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 466 const __gnu_cxx::beta_distribution<_RealType>& __x) 467 { 468 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 469 typedef typename __ostream_type::ios_base __ios_base; 470 471 const typename __ios_base::fmtflags __flags = __os.flags(); 472 const _CharT __fill = __os.fill(); 473 const std::streamsize __precision = __os.precision(); 474 const _CharT __space = __os.widen(' '); 475 __os.flags(__ios_base::scientific | __ios_base::left); 476 __os.fill(__space); 477 __os.precision(std::numeric_limits<_RealType>::max_digits10); 478 479 __os << __x.alpha() << __space << __x.beta(); 480 481 __os.flags(__flags); 482 __os.fill(__fill); 483 __os.precision(__precision); 484 return __os; 485 } 486 487 template<typename _RealType, typename _CharT, typename _Traits> 488 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,__gnu_cxx::beta_distribution<_RealType> & __x)489 operator>>(std::basic_istream<_CharT, _Traits>& __is, 490 __gnu_cxx::beta_distribution<_RealType>& __x) 491 { 492 typedef std::basic_istream<_CharT, _Traits> __istream_type; 493 typedef typename __istream_type::ios_base __ios_base; 494 495 const typename __ios_base::fmtflags __flags = __is.flags(); 496 __is.flags(__ios_base::dec | __ios_base::skipws); 497 498 _RealType __alpha_val, __beta_val; 499 __is >> __alpha_val >> __beta_val; 500 __x.param(typename __gnu_cxx::beta_distribution<_RealType>:: 501 param_type(__alpha_val, __beta_val)); 502 503 __is.flags(__flags); 504 return __is; 505 } 506 507 508 template<std::size_t _Dimen, typename _RealType> 509 template<typename _InputIterator1, typename _InputIterator2> 510 void 511 normal_mv_distribution<_Dimen, _RealType>::param_type:: _M_init_full(_InputIterator1 __meanbegin,_InputIterator1 __meanend,_InputIterator2 __varcovbegin,_InputIterator2 __varcovend)512 _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 513 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend) 514 { 515 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 516 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 517 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 518 _M_mean.end(), _RealType(0)); 519 520 // Perform the Cholesky decomposition 521 auto __w = _M_t.begin(); 522 for (size_t __j = 0; __j < _Dimen; ++__j) 523 { 524 _RealType __sum = _RealType(0); 525 526 auto __slitbegin = __w; 527 auto __cit = _M_t.begin(); 528 for (size_t __i = 0; __i < __j; ++__i) 529 { 530 auto __slit = __slitbegin; 531 _RealType __s = *__varcovbegin++; 532 for (size_t __k = 0; __k < __i; ++__k) 533 __s -= *__slit++ * *__cit++; 534 535 *__w++ = __s /= *__cit++; 536 __sum += __s * __s; 537 } 538 539 __sum = *__varcovbegin - __sum; 540 if (__builtin_expect(__sum <= _RealType(0), 0)) 541 std::__throw_runtime_error(__N("normal_mv_distribution::" 542 "param_type::_M_init_full")); 543 *__w++ = std::sqrt(__sum); 544 545 std::advance(__varcovbegin, _Dimen - __j); 546 } 547 } 548 549 template<std::size_t _Dimen, typename _RealType> 550 template<typename _InputIterator1, typename _InputIterator2> 551 void 552 normal_mv_distribution<_Dimen, _RealType>::param_type:: _M_init_lower(_InputIterator1 __meanbegin,_InputIterator1 __meanend,_InputIterator2 __varcovbegin,_InputIterator2 __varcovend)553 _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 554 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend) 555 { 556 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 557 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 558 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 559 _M_mean.end(), _RealType(0)); 560 561 // Perform the Cholesky decomposition 562 auto __w = _M_t.begin(); 563 for (size_t __j = 0; __j < _Dimen; ++__j) 564 { 565 _RealType __sum = _RealType(0); 566 567 auto __slitbegin = __w; 568 auto __cit = _M_t.begin(); 569 for (size_t __i = 0; __i < __j; ++__i) 570 { 571 auto __slit = __slitbegin; 572 _RealType __s = *__varcovbegin++; 573 for (size_t __k = 0; __k < __i; ++__k) 574 __s -= *__slit++ * *__cit++; 575 576 *__w++ = __s /= *__cit++; 577 __sum += __s * __s; 578 } 579 580 __sum = *__varcovbegin++ - __sum; 581 if (__builtin_expect(__sum <= _RealType(0), 0)) 582 std::__throw_runtime_error(__N("normal_mv_distribution::" 583 "param_type::_M_init_full")); 584 *__w++ = std::sqrt(__sum); 585 } 586 } 587 588 template<std::size_t _Dimen, typename _RealType> 589 template<typename _InputIterator1, typename _InputIterator2> 590 void 591 normal_mv_distribution<_Dimen, _RealType>::param_type:: _M_init_diagonal(_InputIterator1 __meanbegin,_InputIterator1 __meanend,_InputIterator2 __varbegin,_InputIterator2 __varend)592 _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 593 _InputIterator2 __varbegin, _InputIterator2 __varend) 594 { 595 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 596 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 597 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 598 _M_mean.end(), _RealType(0)); 599 600 auto __w = _M_t.begin(); 601 size_t __step = 0; 602 while (__varbegin != __varend) 603 { 604 std::fill_n(__w, __step, _RealType(0)); 605 __w += __step++; 606 if (__builtin_expect(*__varbegin < _RealType(0), 0)) 607 std::__throw_runtime_error(__N("normal_mv_distribution::" 608 "param_type::_M_init_diagonal")); 609 *__w++ = std::sqrt(*__varbegin++); 610 } 611 } 612 613 template<std::size_t _Dimen, typename _RealType> 614 template<typename _UniformRandomNumberGenerator> 615 typename normal_mv_distribution<_Dimen, _RealType>::result_type 616 normal_mv_distribution<_Dimen, _RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)617 operator()(_UniformRandomNumberGenerator& __urng, 618 const param_type& __param) 619 { 620 result_type __ret; 621 622 _M_nd.__generate(__ret.begin(), __ret.end(), __urng); 623 624 auto __t_it = __param._M_t.crbegin(); 625 for (size_t __i = _Dimen; __i > 0; --__i) 626 { 627 _RealType __sum = _RealType(0); 628 for (size_t __j = __i; __j > 0; --__j) 629 __sum += __ret[__j - 1] * *__t_it++; 630 __ret[__i - 1] = __sum; 631 } 632 633 return __ret; 634 } 635 636 template<std::size_t _Dimen, typename _RealType> 637 template<typename _ForwardIterator, typename _UniformRandomNumberGenerator> 638 void 639 normal_mv_distribution<_Dimen, _RealType>:: __generate_impl(_ForwardIterator __f,_ForwardIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)640 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 641 _UniformRandomNumberGenerator& __urng, 642 const param_type& __param) 643 { 644 __glibcxx_function_requires(_Mutable_ForwardIteratorConcept< 645 _ForwardIterator>) 646 while (__f != __t) 647 *__f++ = this->operator()(__urng, __param); 648 } 649 650 template<size_t _Dimen, typename _RealType> 651 bool operator ==(const __gnu_cxx::normal_mv_distribution<_Dimen,_RealType> & __d1,const __gnu_cxx::normal_mv_distribution<_Dimen,_RealType> & __d2)652 operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 653 __d1, 654 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 655 __d2) 656 { 657 return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd; 658 } 659 660 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits> 661 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const __gnu_cxx::normal_mv_distribution<_Dimen,_RealType> & __x)662 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 663 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x) 664 { 665 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 666 typedef typename __ostream_type::ios_base __ios_base; 667 668 const typename __ios_base::fmtflags __flags = __os.flags(); 669 const _CharT __fill = __os.fill(); 670 const std::streamsize __precision = __os.precision(); 671 const _CharT __space = __os.widen(' '); 672 __os.flags(__ios_base::scientific | __ios_base::left); 673 __os.fill(__space); 674 __os.precision(std::numeric_limits<_RealType>::max_digits10); 675 676 auto __mean = __x._M_param.mean(); 677 for (auto __it : __mean) 678 __os << __it << __space; 679 auto __t = __x._M_param.varcov(); 680 for (auto __it : __t) 681 __os << __it << __space; 682 683 __os << __x._M_nd; 684 685 __os.flags(__flags); 686 __os.fill(__fill); 687 __os.precision(__precision); 688 return __os; 689 } 690 691 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits> 692 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,__gnu_cxx::normal_mv_distribution<_Dimen,_RealType> & __x)693 operator>>(std::basic_istream<_CharT, _Traits>& __is, 694 __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x) 695 { 696 typedef std::basic_istream<_CharT, _Traits> __istream_type; 697 typedef typename __istream_type::ios_base __ios_base; 698 699 const typename __ios_base::fmtflags __flags = __is.flags(); 700 __is.flags(__ios_base::dec | __ios_base::skipws); 701 702 std::array<_RealType, _Dimen> __mean; 703 for (auto& __it : __mean) 704 __is >> __it; 705 std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov; 706 for (auto& __it : __varcov) 707 __is >> __it; 708 709 __is >> __x._M_nd; 710 711 __x.param(typename normal_mv_distribution<_Dimen, _RealType>:: 712 param_type(__mean.begin(), __mean.end(), 713 __varcov.begin(), __varcov.end())); 714 715 __is.flags(__flags); 716 return __is; 717 } 718 719 720 template<typename _RealType> 721 template<typename _OutputIterator, 722 typename _UniformRandomNumberGenerator> 723 void 724 rice_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)725 __generate_impl(_OutputIterator __f, _OutputIterator __t, 726 _UniformRandomNumberGenerator& __urng, 727 const param_type& __p) 728 { 729 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 730 result_type>) 731 732 while (__f != __t) 733 { 734 typename std::normal_distribution<result_type>::param_type 735 __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma()); 736 result_type __x = this->_M_ndx(__px, __urng); 737 result_type __y = this->_M_ndy(__py, __urng); 738 #if _GLIBCXX_USE_C99_MATH_TR1 739 *__f++ = std::hypot(__x, __y); 740 #else 741 *__f++ = std::sqrt(__x * __x + __y * __y); 742 #endif 743 } 744 } 745 746 template<typename _RealType, typename _CharT, typename _Traits> 747 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const rice_distribution<_RealType> & __x)748 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 749 const rice_distribution<_RealType>& __x) 750 { 751 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 752 typedef typename __ostream_type::ios_base __ios_base; 753 754 const typename __ios_base::fmtflags __flags = __os.flags(); 755 const _CharT __fill = __os.fill(); 756 const std::streamsize __precision = __os.precision(); 757 const _CharT __space = __os.widen(' '); 758 __os.flags(__ios_base::scientific | __ios_base::left); 759 __os.fill(__space); 760 __os.precision(std::numeric_limits<_RealType>::max_digits10); 761 762 __os << __x.nu() << __space << __x.sigma(); 763 __os << __space << __x._M_ndx; 764 __os << __space << __x._M_ndy; 765 766 __os.flags(__flags); 767 __os.fill(__fill); 768 __os.precision(__precision); 769 return __os; 770 } 771 772 template<typename _RealType, typename _CharT, typename _Traits> 773 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,rice_distribution<_RealType> & __x)774 operator>>(std::basic_istream<_CharT, _Traits>& __is, 775 rice_distribution<_RealType>& __x) 776 { 777 typedef std::basic_istream<_CharT, _Traits> __istream_type; 778 typedef typename __istream_type::ios_base __ios_base; 779 780 const typename __ios_base::fmtflags __flags = __is.flags(); 781 __is.flags(__ios_base::dec | __ios_base::skipws); 782 783 _RealType __nu_val, __sigma_val; 784 __is >> __nu_val >> __sigma_val; 785 __is >> __x._M_ndx; 786 __is >> __x._M_ndy; 787 __x.param(typename rice_distribution<_RealType>:: 788 param_type(__nu_val, __sigma_val)); 789 790 __is.flags(__flags); 791 return __is; 792 } 793 794 795 template<typename _RealType> 796 template<typename _OutputIterator, 797 typename _UniformRandomNumberGenerator> 798 void 799 nakagami_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)800 __generate_impl(_OutputIterator __f, _OutputIterator __t, 801 _UniformRandomNumberGenerator& __urng, 802 const param_type& __p) 803 { 804 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 805 result_type>) 806 807 typename std::gamma_distribution<result_type>::param_type 808 __pg(__p.mu(), __p.omega() / __p.mu()); 809 while (__f != __t) 810 *__f++ = std::sqrt(this->_M_gd(__pg, __urng)); 811 } 812 813 template<typename _RealType, typename _CharT, typename _Traits> 814 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const nakagami_distribution<_RealType> & __x)815 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 816 const nakagami_distribution<_RealType>& __x) 817 { 818 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 819 typedef typename __ostream_type::ios_base __ios_base; 820 821 const typename __ios_base::fmtflags __flags = __os.flags(); 822 const _CharT __fill = __os.fill(); 823 const std::streamsize __precision = __os.precision(); 824 const _CharT __space = __os.widen(' '); 825 __os.flags(__ios_base::scientific | __ios_base::left); 826 __os.fill(__space); 827 __os.precision(std::numeric_limits<_RealType>::max_digits10); 828 829 __os << __x.mu() << __space << __x.omega(); 830 __os << __space << __x._M_gd; 831 832 __os.flags(__flags); 833 __os.fill(__fill); 834 __os.precision(__precision); 835 return __os; 836 } 837 838 template<typename _RealType, typename _CharT, typename _Traits> 839 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,nakagami_distribution<_RealType> & __x)840 operator>>(std::basic_istream<_CharT, _Traits>& __is, 841 nakagami_distribution<_RealType>& __x) 842 { 843 typedef std::basic_istream<_CharT, _Traits> __istream_type; 844 typedef typename __istream_type::ios_base __ios_base; 845 846 const typename __ios_base::fmtflags __flags = __is.flags(); 847 __is.flags(__ios_base::dec | __ios_base::skipws); 848 849 _RealType __mu_val, __omega_val; 850 __is >> __mu_val >> __omega_val; 851 __is >> __x._M_gd; 852 __x.param(typename nakagami_distribution<_RealType>:: 853 param_type(__mu_val, __omega_val)); 854 855 __is.flags(__flags); 856 return __is; 857 } 858 859 860 template<typename _RealType> 861 template<typename _OutputIterator, 862 typename _UniformRandomNumberGenerator> 863 void 864 pareto_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)865 __generate_impl(_OutputIterator __f, _OutputIterator __t, 866 _UniformRandomNumberGenerator& __urng, 867 const param_type& __p) 868 { 869 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 870 result_type>) 871 872 result_type __mu_val = __p.mu(); 873 result_type __malphinv = -result_type(1) / __p.alpha(); 874 while (__f != __t) 875 *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv); 876 } 877 878 template<typename _RealType, typename _CharT, typename _Traits> 879 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const pareto_distribution<_RealType> & __x)880 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 881 const pareto_distribution<_RealType>& __x) 882 { 883 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 884 typedef typename __ostream_type::ios_base __ios_base; 885 886 const typename __ios_base::fmtflags __flags = __os.flags(); 887 const _CharT __fill = __os.fill(); 888 const std::streamsize __precision = __os.precision(); 889 const _CharT __space = __os.widen(' '); 890 __os.flags(__ios_base::scientific | __ios_base::left); 891 __os.fill(__space); 892 __os.precision(std::numeric_limits<_RealType>::max_digits10); 893 894 __os << __x.alpha() << __space << __x.mu(); 895 __os << __space << __x._M_ud; 896 897 __os.flags(__flags); 898 __os.fill(__fill); 899 __os.precision(__precision); 900 return __os; 901 } 902 903 template<typename _RealType, typename _CharT, typename _Traits> 904 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,pareto_distribution<_RealType> & __x)905 operator>>(std::basic_istream<_CharT, _Traits>& __is, 906 pareto_distribution<_RealType>& __x) 907 { 908 typedef std::basic_istream<_CharT, _Traits> __istream_type; 909 typedef typename __istream_type::ios_base __ios_base; 910 911 const typename __ios_base::fmtflags __flags = __is.flags(); 912 __is.flags(__ios_base::dec | __ios_base::skipws); 913 914 _RealType __alpha_val, __mu_val; 915 __is >> __alpha_val >> __mu_val; 916 __is >> __x._M_ud; 917 __x.param(typename pareto_distribution<_RealType>:: 918 param_type(__alpha_val, __mu_val)); 919 920 __is.flags(__flags); 921 return __is; 922 } 923 924 925 template<typename _RealType> 926 template<typename _UniformRandomNumberGenerator> 927 typename k_distribution<_RealType>::result_type 928 k_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng)929 operator()(_UniformRandomNumberGenerator& __urng) 930 { 931 result_type __x = this->_M_gd1(__urng); 932 result_type __y = this->_M_gd2(__urng); 933 return std::sqrt(__x * __y); 934 } 935 936 template<typename _RealType> 937 template<typename _UniformRandomNumberGenerator> 938 typename k_distribution<_RealType>::result_type 939 k_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)940 operator()(_UniformRandomNumberGenerator& __urng, 941 const param_type& __p) 942 { 943 typename std::gamma_distribution<result_type>::param_type 944 __p1(__p.lambda(), result_type(1) / __p.lambda()), 945 __p2(__p.nu(), __p.mu() / __p.nu()); 946 result_type __x = this->_M_gd1(__p1, __urng); 947 result_type __y = this->_M_gd2(__p2, __urng); 948 return std::sqrt(__x * __y); 949 } 950 951 template<typename _RealType> 952 template<typename _OutputIterator, 953 typename _UniformRandomNumberGenerator> 954 void 955 k_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)956 __generate_impl(_OutputIterator __f, _OutputIterator __t, 957 _UniformRandomNumberGenerator& __urng, 958 const param_type& __p) 959 { 960 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 961 result_type>) 962 963 typename std::gamma_distribution<result_type>::param_type 964 __p1(__p.lambda(), result_type(1) / __p.lambda()), 965 __p2(__p.nu(), __p.mu() / __p.nu()); 966 while (__f != __t) 967 { 968 result_type __x = this->_M_gd1(__p1, __urng); 969 result_type __y = this->_M_gd2(__p2, __urng); 970 *__f++ = std::sqrt(__x * __y); 971 } 972 } 973 974 template<typename _RealType, typename _CharT, typename _Traits> 975 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const k_distribution<_RealType> & __x)976 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 977 const k_distribution<_RealType>& __x) 978 { 979 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 980 typedef typename __ostream_type::ios_base __ios_base; 981 982 const typename __ios_base::fmtflags __flags = __os.flags(); 983 const _CharT __fill = __os.fill(); 984 const std::streamsize __precision = __os.precision(); 985 const _CharT __space = __os.widen(' '); 986 __os.flags(__ios_base::scientific | __ios_base::left); 987 __os.fill(__space); 988 __os.precision(std::numeric_limits<_RealType>::max_digits10); 989 990 __os << __x.lambda() << __space << __x.mu() << __space << __x.nu(); 991 __os << __space << __x._M_gd1; 992 __os << __space << __x._M_gd2; 993 994 __os.flags(__flags); 995 __os.fill(__fill); 996 __os.precision(__precision); 997 return __os; 998 } 999 1000 template<typename _RealType, typename _CharT, typename _Traits> 1001 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,k_distribution<_RealType> & __x)1002 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1003 k_distribution<_RealType>& __x) 1004 { 1005 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1006 typedef typename __istream_type::ios_base __ios_base; 1007 1008 const typename __ios_base::fmtflags __flags = __is.flags(); 1009 __is.flags(__ios_base::dec | __ios_base::skipws); 1010 1011 _RealType __lambda_val, __mu_val, __nu_val; 1012 __is >> __lambda_val >> __mu_val >> __nu_val; 1013 __is >> __x._M_gd1; 1014 __is >> __x._M_gd2; 1015 __x.param(typename k_distribution<_RealType>:: 1016 param_type(__lambda_val, __mu_val, __nu_val)); 1017 1018 __is.flags(__flags); 1019 return __is; 1020 } 1021 1022 1023 template<typename _RealType> 1024 template<typename _OutputIterator, 1025 typename _UniformRandomNumberGenerator> 1026 void 1027 arcsine_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)1028 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1029 _UniformRandomNumberGenerator& __urng, 1030 const param_type& __p) 1031 { 1032 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1033 result_type>) 1034 1035 result_type __dif = __p.b() - __p.a(); 1036 result_type __sum = __p.a() + __p.b(); 1037 while (__f != __t) 1038 { 1039 result_type __x = std::sin(this->_M_ud(__urng)); 1040 *__f++ = (__x * __dif + __sum) / result_type(2); 1041 } 1042 } 1043 1044 template<typename _RealType, typename _CharT, typename _Traits> 1045 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const arcsine_distribution<_RealType> & __x)1046 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1047 const arcsine_distribution<_RealType>& __x) 1048 { 1049 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1050 typedef typename __ostream_type::ios_base __ios_base; 1051 1052 const typename __ios_base::fmtflags __flags = __os.flags(); 1053 const _CharT __fill = __os.fill(); 1054 const std::streamsize __precision = __os.precision(); 1055 const _CharT __space = __os.widen(' '); 1056 __os.flags(__ios_base::scientific | __ios_base::left); 1057 __os.fill(__space); 1058 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1059 1060 __os << __x.a() << __space << __x.b(); 1061 __os << __space << __x._M_ud; 1062 1063 __os.flags(__flags); 1064 __os.fill(__fill); 1065 __os.precision(__precision); 1066 return __os; 1067 } 1068 1069 template<typename _RealType, typename _CharT, typename _Traits> 1070 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,arcsine_distribution<_RealType> & __x)1071 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1072 arcsine_distribution<_RealType>& __x) 1073 { 1074 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1075 typedef typename __istream_type::ios_base __ios_base; 1076 1077 const typename __ios_base::fmtflags __flags = __is.flags(); 1078 __is.flags(__ios_base::dec | __ios_base::skipws); 1079 1080 _RealType __a, __b; 1081 __is >> __a >> __b; 1082 __is >> __x._M_ud; 1083 __x.param(typename arcsine_distribution<_RealType>:: 1084 param_type(__a, __b)); 1085 1086 __is.flags(__flags); 1087 return __is; 1088 } 1089 1090 1091 template<typename _RealType> 1092 template<typename _UniformRandomNumberGenerator> 1093 typename hoyt_distribution<_RealType>::result_type 1094 hoyt_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng)1095 operator()(_UniformRandomNumberGenerator& __urng) 1096 { 1097 result_type __x = this->_M_ad(__urng); 1098 result_type __y = this->_M_ed(__urng); 1099 return (result_type(2) * this->q() 1100 / (result_type(1) + this->q() * this->q())) 1101 * std::sqrt(this->omega() * __x * __y); 1102 } 1103 1104 template<typename _RealType> 1105 template<typename _UniformRandomNumberGenerator> 1106 typename hoyt_distribution<_RealType>::result_type 1107 hoyt_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)1108 operator()(_UniformRandomNumberGenerator& __urng, 1109 const param_type& __p) 1110 { 1111 result_type __q2 = __p.q() * __p.q(); 1112 result_type __num = result_type(0.5L) * (result_type(1) + __q2); 1113 typename __gnu_cxx::arcsine_distribution<result_type>::param_type 1114 __pa(__num, __num / __q2); 1115 result_type __x = this->_M_ad(__pa, __urng); 1116 result_type __y = this->_M_ed(__urng); 1117 return (result_type(2) * __p.q() / (result_type(1) + __q2)) 1118 * std::sqrt(__p.omega() * __x * __y); 1119 } 1120 1121 template<typename _RealType> 1122 template<typename _OutputIterator, 1123 typename _UniformRandomNumberGenerator> 1124 void 1125 hoyt_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)1126 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1127 _UniformRandomNumberGenerator& __urng, 1128 const param_type& __p) 1129 { 1130 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1131 result_type>) 1132 1133 result_type __2q = result_type(2) * __p.q(); 1134 result_type __q2 = __p.q() * __p.q(); 1135 result_type __q2p1 = result_type(1) + __q2; 1136 result_type __num = result_type(0.5L) * __q2p1; 1137 result_type __omega = __p.omega(); 1138 typename __gnu_cxx::arcsine_distribution<result_type>::param_type 1139 __pa(__num, __num / __q2); 1140 while (__f != __t) 1141 { 1142 result_type __x = this->_M_ad(__pa, __urng); 1143 result_type __y = this->_M_ed(__urng); 1144 *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y); 1145 } 1146 } 1147 1148 template<typename _RealType, typename _CharT, typename _Traits> 1149 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const hoyt_distribution<_RealType> & __x)1150 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1151 const hoyt_distribution<_RealType>& __x) 1152 { 1153 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1154 typedef typename __ostream_type::ios_base __ios_base; 1155 1156 const typename __ios_base::fmtflags __flags = __os.flags(); 1157 const _CharT __fill = __os.fill(); 1158 const std::streamsize __precision = __os.precision(); 1159 const _CharT __space = __os.widen(' '); 1160 __os.flags(__ios_base::scientific | __ios_base::left); 1161 __os.fill(__space); 1162 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1163 1164 __os << __x.q() << __space << __x.omega(); 1165 __os << __space << __x._M_ad; 1166 __os << __space << __x._M_ed; 1167 1168 __os.flags(__flags); 1169 __os.fill(__fill); 1170 __os.precision(__precision); 1171 return __os; 1172 } 1173 1174 template<typename _RealType, typename _CharT, typename _Traits> 1175 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,hoyt_distribution<_RealType> & __x)1176 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1177 hoyt_distribution<_RealType>& __x) 1178 { 1179 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1180 typedef typename __istream_type::ios_base __ios_base; 1181 1182 const typename __ios_base::fmtflags __flags = __is.flags(); 1183 __is.flags(__ios_base::dec | __ios_base::skipws); 1184 1185 _RealType __q, __omega; 1186 __is >> __q >> __omega; 1187 __is >> __x._M_ad; 1188 __is >> __x._M_ed; 1189 __x.param(typename hoyt_distribution<_RealType>:: 1190 param_type(__q, __omega)); 1191 1192 __is.flags(__flags); 1193 return __is; 1194 } 1195 1196 1197 template<typename _RealType> 1198 template<typename _OutputIterator, 1199 typename _UniformRandomNumberGenerator> 1200 void 1201 triangular_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)1202 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1203 _UniformRandomNumberGenerator& __urng, 1204 const param_type& __param) 1205 { 1206 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1207 result_type>) 1208 1209 while (__f != __t) 1210 *__f++ = this->operator()(__urng, __param); 1211 } 1212 1213 template<typename _RealType, typename _CharT, typename _Traits> 1214 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const __gnu_cxx::triangular_distribution<_RealType> & __x)1215 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1216 const __gnu_cxx::triangular_distribution<_RealType>& __x) 1217 { 1218 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1219 typedef typename __ostream_type::ios_base __ios_base; 1220 1221 const typename __ios_base::fmtflags __flags = __os.flags(); 1222 const _CharT __fill = __os.fill(); 1223 const std::streamsize __precision = __os.precision(); 1224 const _CharT __space = __os.widen(' '); 1225 __os.flags(__ios_base::scientific | __ios_base::left); 1226 __os.fill(__space); 1227 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1228 1229 __os << __x.a() << __space << __x.b() << __space << __x.c(); 1230 1231 __os.flags(__flags); 1232 __os.fill(__fill); 1233 __os.precision(__precision); 1234 return __os; 1235 } 1236 1237 template<typename _RealType, typename _CharT, typename _Traits> 1238 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,__gnu_cxx::triangular_distribution<_RealType> & __x)1239 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1240 __gnu_cxx::triangular_distribution<_RealType>& __x) 1241 { 1242 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1243 typedef typename __istream_type::ios_base __ios_base; 1244 1245 const typename __ios_base::fmtflags __flags = __is.flags(); 1246 __is.flags(__ios_base::dec | __ios_base::skipws); 1247 1248 _RealType __a, __b, __c; 1249 __is >> __a >> __b >> __c; 1250 __x.param(typename __gnu_cxx::triangular_distribution<_RealType>:: 1251 param_type(__a, __b, __c)); 1252 1253 __is.flags(__flags); 1254 return __is; 1255 } 1256 1257 1258 template<typename _RealType> 1259 template<typename _UniformRandomNumberGenerator> 1260 typename von_mises_distribution<_RealType>::result_type 1261 von_mises_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)1262 operator()(_UniformRandomNumberGenerator& __urng, 1263 const param_type& __p) 1264 { 1265 const result_type __pi 1266 = __gnu_cxx::__math_constants<result_type>::__pi; 1267 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1268 __aurng(__urng); 1269 1270 result_type __f; 1271 while (1) 1272 { 1273 result_type __rnd = std::cos(__pi * __aurng()); 1274 __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd); 1275 result_type __c = __p._M_kappa * (__p._M_r - __f); 1276 1277 result_type __rnd2 = __aurng(); 1278 if (__c * (result_type(2) - __c) > __rnd2) 1279 break; 1280 if (std::log(__c / __rnd2) >= __c - result_type(1)) 1281 break; 1282 } 1283 1284 result_type __res = std::acos(__f); 1285 #if _GLIBCXX_USE_C99_MATH_TR1 1286 __res = std::copysign(__res, __aurng() - result_type(0.5)); 1287 #else 1288 if (__aurng() < result_type(0.5)) 1289 __res = -__res; 1290 #endif 1291 __res += __p._M_mu; 1292 if (__res > __pi) 1293 __res -= result_type(2) * __pi; 1294 else if (__res < -__pi) 1295 __res += result_type(2) * __pi; 1296 return __res; 1297 } 1298 1299 template<typename _RealType> 1300 template<typename _OutputIterator, 1301 typename _UniformRandomNumberGenerator> 1302 void 1303 von_mises_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)1304 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1305 _UniformRandomNumberGenerator& __urng, 1306 const param_type& __param) 1307 { 1308 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1309 result_type>) 1310 1311 while (__f != __t) 1312 *__f++ = this->operator()(__urng, __param); 1313 } 1314 1315 template<typename _RealType, typename _CharT, typename _Traits> 1316 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const __gnu_cxx::von_mises_distribution<_RealType> & __x)1317 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1318 const __gnu_cxx::von_mises_distribution<_RealType>& __x) 1319 { 1320 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1321 typedef typename __ostream_type::ios_base __ios_base; 1322 1323 const typename __ios_base::fmtflags __flags = __os.flags(); 1324 const _CharT __fill = __os.fill(); 1325 const std::streamsize __precision = __os.precision(); 1326 const _CharT __space = __os.widen(' '); 1327 __os.flags(__ios_base::scientific | __ios_base::left); 1328 __os.fill(__space); 1329 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1330 1331 __os << __x.mu() << __space << __x.kappa(); 1332 1333 __os.flags(__flags); 1334 __os.fill(__fill); 1335 __os.precision(__precision); 1336 return __os; 1337 } 1338 1339 template<typename _RealType, typename _CharT, typename _Traits> 1340 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,__gnu_cxx::von_mises_distribution<_RealType> & __x)1341 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1342 __gnu_cxx::von_mises_distribution<_RealType>& __x) 1343 { 1344 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1345 typedef typename __istream_type::ios_base __ios_base; 1346 1347 const typename __ios_base::fmtflags __flags = __is.flags(); 1348 __is.flags(__ios_base::dec | __ios_base::skipws); 1349 1350 _RealType __mu, __kappa; 1351 __is >> __mu >> __kappa; 1352 __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>:: 1353 param_type(__mu, __kappa)); 1354 1355 __is.flags(__flags); 1356 return __is; 1357 } 1358 1359 1360 template<typename _UIntType> 1361 template<typename _UniformRandomNumberGenerator> 1362 typename hypergeometric_distribution<_UIntType>::result_type 1363 hypergeometric_distribution<_UIntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)1364 operator()(_UniformRandomNumberGenerator& __urng, 1365 const param_type& __param) 1366 { 1367 std::__detail::_Adaptor<_UniformRandomNumberGenerator, double> 1368 __aurng(__urng); 1369 1370 result_type __a = __param.successful_size(); 1371 result_type __b = __param.total_size(); 1372 result_type __k = 0; 1373 1374 if (__param.total_draws() < __param.total_size() / 2) 1375 { 1376 for (result_type __i = 0; __i < __param.total_draws(); ++__i) 1377 { 1378 if (__b * __aurng() < __a) 1379 { 1380 ++__k; 1381 if (__k == __param.successful_size()) 1382 return __k; 1383 --__a; 1384 } 1385 --__b; 1386 } 1387 return __k; 1388 } 1389 else 1390 { 1391 for (result_type __i = 0; __i < __param.unsuccessful_size(); ++__i) 1392 { 1393 if (__b * __aurng() < __a) 1394 { 1395 ++__k; 1396 if (__k == __param.successful_size()) 1397 return __param.successful_size() - __k; 1398 --__a; 1399 } 1400 --__b; 1401 } 1402 return __param.successful_size() - __k; 1403 } 1404 } 1405 1406 template<typename _UIntType> 1407 template<typename _OutputIterator, 1408 typename _UniformRandomNumberGenerator> 1409 void 1410 hypergeometric_distribution<_UIntType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)1411 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1412 _UniformRandomNumberGenerator& __urng, 1413 const param_type& __param) 1414 { 1415 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1416 result_type>) 1417 1418 while (__f != __t) 1419 *__f++ = this->operator()(__urng); 1420 } 1421 1422 template<typename _UIntType, typename _CharT, typename _Traits> 1423 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const __gnu_cxx::hypergeometric_distribution<_UIntType> & __x)1424 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1425 const __gnu_cxx::hypergeometric_distribution<_UIntType>& __x) 1426 { 1427 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1428 typedef typename __ostream_type::ios_base __ios_base; 1429 1430 const typename __ios_base::fmtflags __flags = __os.flags(); 1431 const _CharT __fill = __os.fill(); 1432 const std::streamsize __precision = __os.precision(); 1433 const _CharT __space = __os.widen(' '); 1434 __os.flags(__ios_base::scientific | __ios_base::left); 1435 __os.fill(__space); 1436 __os.precision(std::numeric_limits<_UIntType>::max_digits10); 1437 1438 __os << __x.total_size() << __space << __x.successful_size() << __space 1439 << __x.total_draws(); 1440 1441 __os.flags(__flags); 1442 __os.fill(__fill); 1443 __os.precision(__precision); 1444 return __os; 1445 } 1446 1447 template<typename _UIntType, typename _CharT, typename _Traits> 1448 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,__gnu_cxx::hypergeometric_distribution<_UIntType> & __x)1449 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1450 __gnu_cxx::hypergeometric_distribution<_UIntType>& __x) 1451 { 1452 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1453 typedef typename __istream_type::ios_base __ios_base; 1454 1455 const typename __ios_base::fmtflags __flags = __is.flags(); 1456 __is.flags(__ios_base::dec | __ios_base::skipws); 1457 1458 _UIntType __total_size, __successful_size, __total_draws; 1459 __is >> __total_size >> __successful_size >> __total_draws; 1460 __x.param(typename __gnu_cxx::hypergeometric_distribution<_UIntType>:: 1461 param_type(__total_size, __successful_size, __total_draws)); 1462 1463 __is.flags(__flags); 1464 return __is; 1465 } 1466 1467 1468 template<typename _RealType> 1469 template<typename _UniformRandomNumberGenerator> 1470 typename logistic_distribution<_RealType>::result_type 1471 logistic_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)1472 operator()(_UniformRandomNumberGenerator& __urng, 1473 const param_type& __p) 1474 { 1475 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1476 __aurng(__urng); 1477 1478 result_type __arg = result_type(1); 1479 while (__arg == result_type(1) || __arg == result_type(0)) 1480 __arg = __aurng(); 1481 return __p.a() 1482 + __p.b() * std::log(__arg / (result_type(1) - __arg)); 1483 } 1484 1485 template<typename _RealType> 1486 template<typename _OutputIterator, 1487 typename _UniformRandomNumberGenerator> 1488 void 1489 logistic_distribution<_RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __p)1490 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1491 _UniformRandomNumberGenerator& __urng, 1492 const param_type& __p) 1493 { 1494 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1495 result_type>) 1496 1497 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1498 __aurng(__urng); 1499 1500 while (__f != __t) 1501 { 1502 result_type __arg = result_type(1); 1503 while (__arg == result_type(1) || __arg == result_type(0)) 1504 __arg = __aurng(); 1505 *__f++ = __p.a() 1506 + __p.b() * std::log(__arg / (result_type(1) - __arg)); 1507 } 1508 } 1509 1510 template<typename _RealType, typename _CharT, typename _Traits> 1511 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const logistic_distribution<_RealType> & __x)1512 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1513 const logistic_distribution<_RealType>& __x) 1514 { 1515 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1516 typedef typename __ostream_type::ios_base __ios_base; 1517 1518 const typename __ios_base::fmtflags __flags = __os.flags(); 1519 const _CharT __fill = __os.fill(); 1520 const std::streamsize __precision = __os.precision(); 1521 const _CharT __space = __os.widen(' '); 1522 __os.flags(__ios_base::scientific | __ios_base::left); 1523 __os.fill(__space); 1524 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1525 1526 __os << __x.a() << __space << __x.b(); 1527 1528 __os.flags(__flags); 1529 __os.fill(__fill); 1530 __os.precision(__precision); 1531 return __os; 1532 } 1533 1534 template<typename _RealType, typename _CharT, typename _Traits> 1535 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,logistic_distribution<_RealType> & __x)1536 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1537 logistic_distribution<_RealType>& __x) 1538 { 1539 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1540 typedef typename __istream_type::ios_base __ios_base; 1541 1542 const typename __ios_base::fmtflags __flags = __is.flags(); 1543 __is.flags(__ios_base::dec | __ios_base::skipws); 1544 1545 _RealType __a, __b; 1546 __is >> __a >> __b; 1547 __x.param(typename logistic_distribution<_RealType>:: 1548 param_type(__a, __b)); 1549 1550 __is.flags(__flags); 1551 return __is; 1552 } 1553 1554 1555 namespace { 1556 1557 // Helper class for the uniform_on_sphere_distribution generation 1558 // function. 1559 template<std::size_t _Dimen, typename _RealType> 1560 class uniform_on_sphere_helper 1561 { 1562 typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>:: 1563 result_type result_type; 1564 1565 public: 1566 template<typename _NormalDistribution, 1567 typename _UniformRandomNumberGenerator> operator ()(_NormalDistribution & __nd,_UniformRandomNumberGenerator & __urng)1568 result_type operator()(_NormalDistribution& __nd, 1569 _UniformRandomNumberGenerator& __urng) 1570 { 1571 result_type __ret; 1572 typename result_type::value_type __norm; 1573 1574 do 1575 { 1576 auto __sum = _RealType(0); 1577 1578 std::generate(__ret.begin(), __ret.end(), 1579 [&__nd, &__urng, &__sum](){ 1580 _RealType __t = __nd(__urng); 1581 __sum += __t * __t; 1582 return __t; }); 1583 __norm = std::sqrt(__sum); 1584 } 1585 while (__norm == _RealType(0) || ! __builtin_isfinite(__norm)); 1586 1587 std::transform(__ret.begin(), __ret.end(), __ret.begin(), 1588 [__norm](_RealType __val){ return __val / __norm; }); 1589 1590 return __ret; 1591 } 1592 }; 1593 1594 1595 template<typename _RealType> 1596 class uniform_on_sphere_helper<2, _RealType> 1597 { 1598 typedef typename uniform_on_sphere_distribution<2, _RealType>:: 1599 result_type result_type; 1600 1601 public: 1602 template<typename _NormalDistribution, 1603 typename _UniformRandomNumberGenerator> operator ()(_NormalDistribution &,_UniformRandomNumberGenerator & __urng)1604 result_type operator()(_NormalDistribution&, 1605 _UniformRandomNumberGenerator& __urng) 1606 { 1607 result_type __ret; 1608 _RealType __sq; 1609 std::__detail::_Adaptor<_UniformRandomNumberGenerator, 1610 _RealType> __aurng(__urng); 1611 1612 do 1613 { 1614 __ret[0] = _RealType(2) * __aurng() - _RealType(1); 1615 __ret[1] = _RealType(2) * __aurng() - _RealType(1); 1616 1617 __sq = __ret[0] * __ret[0] + __ret[1] * __ret[1]; 1618 } 1619 while (__sq == _RealType(0) || __sq > _RealType(1)); 1620 1621 #if _GLIBCXX_USE_C99_MATH_TR1 1622 // Yes, we do not just use sqrt(__sq) because hypot() is more 1623 // accurate. 1624 auto __norm = std::hypot(__ret[0], __ret[1]); 1625 #else 1626 auto __norm = std::sqrt(__sq); 1627 #endif 1628 __ret[0] /= __norm; 1629 __ret[1] /= __norm; 1630 1631 return __ret; 1632 } 1633 }; 1634 1635 } 1636 1637 1638 template<std::size_t _Dimen, typename _RealType> 1639 template<typename _UniformRandomNumberGenerator> 1640 typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type 1641 uniform_on_sphere_distribution<_Dimen, _RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)1642 operator()(_UniformRandomNumberGenerator& __urng, 1643 const param_type& __p) 1644 { 1645 uniform_on_sphere_helper<_Dimen, _RealType> __helper; 1646 return __helper(_M_nd, __urng); 1647 } 1648 1649 template<std::size_t _Dimen, typename _RealType> 1650 template<typename _OutputIterator, 1651 typename _UniformRandomNumberGenerator> 1652 void 1653 uniform_on_sphere_distribution<_Dimen, _RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)1654 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1655 _UniformRandomNumberGenerator& __urng, 1656 const param_type& __param) 1657 { 1658 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1659 result_type>) 1660 1661 while (__f != __t) 1662 *__f++ = this->operator()(__urng, __param); 1663 } 1664 1665 template<std::size_t _Dimen, typename _RealType, typename _CharT, 1666 typename _Traits> 1667 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const __gnu_cxx::uniform_on_sphere_distribution<_Dimen,_RealType> & __x)1668 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1669 const __gnu_cxx::uniform_on_sphere_distribution<_Dimen, 1670 _RealType>& __x) 1671 { 1672 return __os << __x._M_nd; 1673 } 1674 1675 template<std::size_t _Dimen, typename _RealType, typename _CharT, 1676 typename _Traits> 1677 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,__gnu_cxx::uniform_on_sphere_distribution<_Dimen,_RealType> & __x)1678 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1679 __gnu_cxx::uniform_on_sphere_distribution<_Dimen, 1680 _RealType>& __x) 1681 { 1682 return __is >> __x._M_nd; 1683 } 1684 1685 1686 namespace { 1687 1688 // Helper class for the uniform_inside_sphere_distribution generation 1689 // function. 1690 template<std::size_t _Dimen, bool _SmallDimen, typename _RealType> 1691 class uniform_inside_sphere_helper; 1692 1693 template<std::size_t _Dimen, typename _RealType> 1694 class uniform_inside_sphere_helper<_Dimen, false, _RealType> 1695 { 1696 using result_type 1697 = typename uniform_inside_sphere_distribution<_Dimen, _RealType>:: 1698 result_type; 1699 1700 public: 1701 template<typename _UniformOnSphereDistribution, 1702 typename _UniformRandomNumberGenerator> 1703 result_type operator ()(_UniformOnSphereDistribution & __uosd,_UniformRandomNumberGenerator & __urng,_RealType __radius)1704 operator()(_UniformOnSphereDistribution& __uosd, 1705 _UniformRandomNumberGenerator& __urng, 1706 _RealType __radius) 1707 { 1708 std::__detail::_Adaptor<_UniformRandomNumberGenerator, 1709 _RealType> __aurng(__urng); 1710 1711 _RealType __pow = 1 / _RealType(_Dimen); 1712 _RealType __urt = __radius * std::pow(__aurng(), __pow); 1713 result_type __ret = __uosd(__aurng); 1714 1715 std::transform(__ret.begin(), __ret.end(), __ret.begin(), 1716 [__urt](_RealType __val) 1717 { return __val * __urt; }); 1718 1719 return __ret; 1720 } 1721 }; 1722 1723 // Helper class for the uniform_inside_sphere_distribution generation 1724 // function specialized for small dimensions. 1725 template<std::size_t _Dimen, typename _RealType> 1726 class uniform_inside_sphere_helper<_Dimen, true, _RealType> 1727 { 1728 using result_type 1729 = typename uniform_inside_sphere_distribution<_Dimen, _RealType>:: 1730 result_type; 1731 1732 public: 1733 template<typename _UniformOnSphereDistribution, 1734 typename _UniformRandomNumberGenerator> 1735 result_type operator ()(_UniformOnSphereDistribution &,_UniformRandomNumberGenerator & __urng,_RealType __radius)1736 operator()(_UniformOnSphereDistribution&, 1737 _UniformRandomNumberGenerator& __urng, 1738 _RealType __radius) 1739 { 1740 result_type __ret; 1741 _RealType __sq; 1742 _RealType __radsq = __radius * __radius; 1743 std::__detail::_Adaptor<_UniformRandomNumberGenerator, 1744 _RealType> __aurng(__urng); 1745 1746 do 1747 { 1748 __sq = _RealType(0); 1749 for (int i = 0; i < _Dimen; ++i) 1750 { 1751 __ret[i] = _RealType(2) * __aurng() - _RealType(1); 1752 __sq += __ret[i] * __ret[i]; 1753 } 1754 } 1755 while (__sq > _RealType(1)); 1756 1757 for (int i = 0; i < _Dimen; ++i) 1758 __ret[i] *= __radius; 1759 1760 return __ret; 1761 } 1762 }; 1763 } // namespace 1764 1765 // 1766 // Experiments have shown that rejection is more efficient than transform 1767 // for dimensions less than 8. 1768 // 1769 template<std::size_t _Dimen, typename _RealType> 1770 template<typename _UniformRandomNumberGenerator> 1771 typename uniform_inside_sphere_distribution<_Dimen, _RealType>::result_type 1772 uniform_inside_sphere_distribution<_Dimen, _RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)1773 operator()(_UniformRandomNumberGenerator& __urng, 1774 const param_type& __p) 1775 { 1776 uniform_inside_sphere_helper<_Dimen, _Dimen < 8, _RealType> __helper; 1777 return __helper(_M_uosd, __urng, __p.radius()); 1778 } 1779 1780 template<std::size_t _Dimen, typename _RealType> 1781 template<typename _OutputIterator, 1782 typename _UniformRandomNumberGenerator> 1783 void 1784 uniform_inside_sphere_distribution<_Dimen, _RealType>:: __generate_impl(_OutputIterator __f,_OutputIterator __t,_UniformRandomNumberGenerator & __urng,const param_type & __param)1785 __generate_impl(_OutputIterator __f, _OutputIterator __t, 1786 _UniformRandomNumberGenerator& __urng, 1787 const param_type& __param) 1788 { 1789 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator, 1790 result_type>) 1791 1792 while (__f != __t) 1793 *__f++ = this->operator()(__urng, __param); 1794 } 1795 1796 template<std::size_t _Dimen, typename _RealType, typename _CharT, 1797 typename _Traits> 1798 std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const __gnu_cxx::uniform_inside_sphere_distribution<_Dimen,_RealType> & __x)1799 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1800 const __gnu_cxx::uniform_inside_sphere_distribution<_Dimen, 1801 _RealType>& __x) 1802 { 1803 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1804 typedef typename __ostream_type::ios_base __ios_base; 1805 1806 const typename __ios_base::fmtflags __flags = __os.flags(); 1807 const _CharT __fill = __os.fill(); 1808 const std::streamsize __precision = __os.precision(); 1809 const _CharT __space = __os.widen(' '); 1810 __os.flags(__ios_base::scientific | __ios_base::left); 1811 __os.fill(__space); 1812 __os.precision(std::numeric_limits<_RealType>::max_digits10); 1813 1814 __os << __x.radius() << __space << __x._M_uosd; 1815 1816 __os.flags(__flags); 1817 __os.fill(__fill); 1818 __os.precision(__precision); 1819 1820 return __os; 1821 } 1822 1823 template<std::size_t _Dimen, typename _RealType, typename _CharT, 1824 typename _Traits> 1825 std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,__gnu_cxx::uniform_inside_sphere_distribution<_Dimen,_RealType> & __x)1826 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1827 __gnu_cxx::uniform_inside_sphere_distribution<_Dimen, 1828 _RealType>& __x) 1829 { 1830 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1831 typedef typename __istream_type::ios_base __ios_base; 1832 1833 const typename __ios_base::fmtflags __flags = __is.flags(); 1834 __is.flags(__ios_base::dec | __ios_base::skipws); 1835 1836 _RealType __radius_val; 1837 __is >> __radius_val >> __x._M_uosd; 1838 __x.param(typename uniform_inside_sphere_distribution<_Dimen, _RealType>:: 1839 param_type(__radius_val)); 1840 1841 __is.flags(__flags); 1842 1843 return __is; 1844 } 1845 1846 _GLIBCXX_END_NAMESPACE_VERSION 1847 } // namespace __gnu_cxx 1848 1849 1850 #endif // _EXT_RANDOM_TCC 1851