libstdc++
|
00001 // Random number extensions -*- C++ -*- 00002 00003 // Copyright (C) 2012-2015 Free Software Foundation, Inc. 00004 // 00005 // This file is part of the GNU ISO C++ Library. This library is free 00006 // software; you can redistribute it and/or modify it under the 00007 // terms of the GNU General Public License as published by the 00008 // Free Software Foundation; either version 3, or (at your option) 00009 // any later version. 00010 00011 // This library is distributed in the hope that it will be useful, 00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 // GNU General Public License for more details. 00015 00016 // Under Section 7 of GPL version 3, you are granted additional 00017 // permissions described in the GCC Runtime Library Exception, version 00018 // 3.1, as published by the Free Software Foundation. 00019 00020 // You should have received a copy of the GNU General Public License and 00021 // a copy of the GCC Runtime Library Exception along with this program; 00022 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 00023 // <http://www.gnu.org/licenses/>. 00024 00025 /** @file ext/random.tcc 00026 * This is an internal header file, included by other library headers. 00027 * Do not attempt to use it directly. @headername{ext/random} 00028 */ 00029 00030 #ifndef _EXT_RANDOM_TCC 00031 #define _EXT_RANDOM_TCC 1 00032 00033 #pragma GCC system_header 00034 00035 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default) 00036 { 00037 _GLIBCXX_BEGIN_NAMESPACE_VERSION 00038 00039 #if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ 00040 00041 template<typename _UIntType, size_t __m, 00042 size_t __pos1, size_t __sl1, size_t __sl2, 00043 size_t __sr1, size_t __sr2, 00044 uint32_t __msk1, uint32_t __msk2, 00045 uint32_t __msk3, uint32_t __msk4, 00046 uint32_t __parity1, uint32_t __parity2, 00047 uint32_t __parity3, uint32_t __parity4> 00048 void simd_fast_mersenne_twister_engine<_UIntType, __m, 00049 __pos1, __sl1, __sl2, __sr1, __sr2, 00050 __msk1, __msk2, __msk3, __msk4, 00051 __parity1, __parity2, __parity3, 00052 __parity4>:: 00053 seed(_UIntType __seed) 00054 { 00055 _M_state32[0] = static_cast<uint32_t>(__seed); 00056 for (size_t __i = 1; __i < _M_nstate32; ++__i) 00057 _M_state32[__i] = (1812433253UL 00058 * (_M_state32[__i - 1] ^ (_M_state32[__i - 1] >> 30)) 00059 + __i); 00060 _M_pos = state_size; 00061 _M_period_certification(); 00062 } 00063 00064 00065 namespace { 00066 00067 inline uint32_t _Func1(uint32_t __x) 00068 { 00069 return (__x ^ (__x >> 27)) * UINT32_C(1664525); 00070 } 00071 00072 inline uint32_t _Func2(uint32_t __x) 00073 { 00074 return (__x ^ (__x >> 27)) * UINT32_C(1566083941); 00075 } 00076 00077 } 00078 00079 00080 template<typename _UIntType, size_t __m, 00081 size_t __pos1, size_t __sl1, size_t __sl2, 00082 size_t __sr1, size_t __sr2, 00083 uint32_t __msk1, uint32_t __msk2, 00084 uint32_t __msk3, uint32_t __msk4, 00085 uint32_t __parity1, uint32_t __parity2, 00086 uint32_t __parity3, uint32_t __parity4> 00087 template<typename _Sseq> 00088 typename std::enable_if<std::is_class<_Sseq>::value>::type 00089 simd_fast_mersenne_twister_engine<_UIntType, __m, 00090 __pos1, __sl1, __sl2, __sr1, __sr2, 00091 __msk1, __msk2, __msk3, __msk4, 00092 __parity1, __parity2, __parity3, 00093 __parity4>:: 00094 seed(_Sseq& __q) 00095 { 00096 size_t __lag; 00097 00098 if (_M_nstate32 >= 623) 00099 __lag = 11; 00100 else if (_M_nstate32 >= 68) 00101 __lag = 7; 00102 else if (_M_nstate32 >= 39) 00103 __lag = 5; 00104 else 00105 __lag = 3; 00106 const size_t __mid = (_M_nstate32 - __lag) / 2; 00107 00108 std::fill(_M_state32, _M_state32 + _M_nstate32, UINT32_C(0x8b8b8b8b)); 00109 uint32_t __arr[_M_nstate32]; 00110 __q.generate(__arr + 0, __arr + _M_nstate32); 00111 00112 uint32_t __r = _Func1(_M_state32[0] ^ _M_state32[__mid] 00113 ^ _M_state32[_M_nstate32 - 1]); 00114 _M_state32[__mid] += __r; 00115 __r += _M_nstate32; 00116 _M_state32[__mid + __lag] += __r; 00117 _M_state32[0] = __r; 00118 00119 for (size_t __i = 1, __j = 0; __j < _M_nstate32; ++__j) 00120 { 00121 __r = _Func1(_M_state32[__i] 00122 ^ _M_state32[(__i + __mid) % _M_nstate32] 00123 ^ _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]); 00124 _M_state32[(__i + __mid) % _M_nstate32] += __r; 00125 __r += __arr[__j] + __i; 00126 _M_state32[(__i + __mid + __lag) % _M_nstate32] += __r; 00127 _M_state32[__i] = __r; 00128 __i = (__i + 1) % _M_nstate32; 00129 } 00130 for (size_t __j = 0; __j < _M_nstate32; ++__j) 00131 { 00132 const size_t __i = (__j + 1) % _M_nstate32; 00133 __r = _Func2(_M_state32[__i] 00134 + _M_state32[(__i + __mid) % _M_nstate32] 00135 + _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]); 00136 _M_state32[(__i + __mid) % _M_nstate32] ^= __r; 00137 __r -= __i; 00138 _M_state32[(__i + __mid + __lag) % _M_nstate32] ^= __r; 00139 _M_state32[__i] = __r; 00140 } 00141 00142 _M_pos = state_size; 00143 _M_period_certification(); 00144 } 00145 00146 00147 template<typename _UIntType, size_t __m, 00148 size_t __pos1, size_t __sl1, size_t __sl2, 00149 size_t __sr1, size_t __sr2, 00150 uint32_t __msk1, uint32_t __msk2, 00151 uint32_t __msk3, uint32_t __msk4, 00152 uint32_t __parity1, uint32_t __parity2, 00153 uint32_t __parity3, uint32_t __parity4> 00154 void simd_fast_mersenne_twister_engine<_UIntType, __m, 00155 __pos1, __sl1, __sl2, __sr1, __sr2, 00156 __msk1, __msk2, __msk3, __msk4, 00157 __parity1, __parity2, __parity3, 00158 __parity4>:: 00159 _M_period_certification(void) 00160 { 00161 static const uint32_t __parity[4] = { __parity1, __parity2, 00162 __parity3, __parity4 }; 00163 uint32_t __inner = 0; 00164 for (size_t __i = 0; __i < 4; ++__i) 00165 if (__parity[__i] != 0) 00166 __inner ^= _M_state32[__i] & __parity[__i]; 00167 00168 if (__builtin_parity(__inner) & 1) 00169 return; 00170 for (size_t __i = 0; __i < 4; ++__i) 00171 if (__parity[__i] != 0) 00172 { 00173 _M_state32[__i] ^= 1 << (__builtin_ffs(__parity[__i]) - 1); 00174 return; 00175 } 00176 __builtin_unreachable(); 00177 } 00178 00179 00180 template<typename _UIntType, size_t __m, 00181 size_t __pos1, size_t __sl1, size_t __sl2, 00182 size_t __sr1, size_t __sr2, 00183 uint32_t __msk1, uint32_t __msk2, 00184 uint32_t __msk3, uint32_t __msk4, 00185 uint32_t __parity1, uint32_t __parity2, 00186 uint32_t __parity3, uint32_t __parity4> 00187 void simd_fast_mersenne_twister_engine<_UIntType, __m, 00188 __pos1, __sl1, __sl2, __sr1, __sr2, 00189 __msk1, __msk2, __msk3, __msk4, 00190 __parity1, __parity2, __parity3, 00191 __parity4>:: 00192 discard(unsigned long long __z) 00193 { 00194 while (__z > state_size - _M_pos) 00195 { 00196 __z -= state_size - _M_pos; 00197 00198 _M_gen_rand(); 00199 } 00200 00201 _M_pos += __z; 00202 } 00203 00204 00205 #ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_GEN_READ 00206 00207 namespace { 00208 00209 template<size_t __shift> 00210 inline void __rshift(uint32_t *__out, const uint32_t *__in) 00211 { 00212 uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32) 00213 | static_cast<uint64_t>(__in[2])); 00214 uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32) 00215 | static_cast<uint64_t>(__in[0])); 00216 00217 uint64_t __oh = __th >> (__shift * 8); 00218 uint64_t __ol = __tl >> (__shift * 8); 00219 __ol |= __th << (64 - __shift * 8); 00220 __out[1] = static_cast<uint32_t>(__ol >> 32); 00221 __out[0] = static_cast<uint32_t>(__ol); 00222 __out[3] = static_cast<uint32_t>(__oh >> 32); 00223 __out[2] = static_cast<uint32_t>(__oh); 00224 } 00225 00226 00227 template<size_t __shift> 00228 inline void __lshift(uint32_t *__out, const uint32_t *__in) 00229 { 00230 uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32) 00231 | static_cast<uint64_t>(__in[2])); 00232 uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32) 00233 | static_cast<uint64_t>(__in[0])); 00234 00235 uint64_t __oh = __th << (__shift * 8); 00236 uint64_t __ol = __tl << (__shift * 8); 00237 __oh |= __tl >> (64 - __shift * 8); 00238 __out[1] = static_cast<uint32_t>(__ol >> 32); 00239 __out[0] = static_cast<uint32_t>(__ol); 00240 __out[3] = static_cast<uint32_t>(__oh >> 32); 00241 __out[2] = static_cast<uint32_t>(__oh); 00242 } 00243 00244 00245 template<size_t __sl1, size_t __sl2, size_t __sr1, size_t __sr2, 00246 uint32_t __msk1, uint32_t __msk2, uint32_t __msk3, uint32_t __msk4> 00247 inline void __recursion(uint32_t *__r, 00248 const uint32_t *__a, const uint32_t *__b, 00249 const uint32_t *__c, const uint32_t *__d) 00250 { 00251 uint32_t __x[4]; 00252 uint32_t __y[4]; 00253 00254 __lshift<__sl2>(__x, __a); 00255 __rshift<__sr2>(__y, __c); 00256 __r[0] = (__a[0] ^ __x[0] ^ ((__b[0] >> __sr1) & __msk1) 00257 ^ __y[0] ^ (__d[0] << __sl1)); 00258 __r[1] = (__a[1] ^ __x[1] ^ ((__b[1] >> __sr1) & __msk2) 00259 ^ __y[1] ^ (__d[1] << __sl1)); 00260 __r[2] = (__a[2] ^ __x[2] ^ ((__b[2] >> __sr1) & __msk3) 00261 ^ __y[2] ^ (__d[2] << __sl1)); 00262 __r[3] = (__a[3] ^ __x[3] ^ ((__b[3] >> __sr1) & __msk4) 00263 ^ __y[3] ^ (__d[3] << __sl1)); 00264 } 00265 00266 } 00267 00268 00269 template<typename _UIntType, size_t __m, 00270 size_t __pos1, size_t __sl1, size_t __sl2, 00271 size_t __sr1, size_t __sr2, 00272 uint32_t __msk1, uint32_t __msk2, 00273 uint32_t __msk3, uint32_t __msk4, 00274 uint32_t __parity1, uint32_t __parity2, 00275 uint32_t __parity3, uint32_t __parity4> 00276 void simd_fast_mersenne_twister_engine<_UIntType, __m, 00277 __pos1, __sl1, __sl2, __sr1, __sr2, 00278 __msk1, __msk2, __msk3, __msk4, 00279 __parity1, __parity2, __parity3, 00280 __parity4>:: 00281 _M_gen_rand(void) 00282 { 00283 const uint32_t *__r1 = &_M_state32[_M_nstate32 - 8]; 00284 const uint32_t *__r2 = &_M_state32[_M_nstate32 - 4]; 00285 static constexpr size_t __pos1_32 = __pos1 * 4; 00286 00287 size_t __i; 00288 for (__i = 0; __i < _M_nstate32 - __pos1_32; __i += 4) 00289 { 00290 __recursion<__sl1, __sl2, __sr1, __sr2, 00291 __msk1, __msk2, __msk3, __msk4> 00292 (&_M_state32[__i], &_M_state32[__i], 00293 &_M_state32[__i + __pos1_32], __r1, __r2); 00294 __r1 = __r2; 00295 __r2 = &_M_state32[__i]; 00296 } 00297 00298 for (; __i < _M_nstate32; __i += 4) 00299 { 00300 __recursion<__sl1, __sl2, __sr1, __sr2, 00301 __msk1, __msk2, __msk3, __msk4> 00302 (&_M_state32[__i], &_M_state32[__i], 00303 &_M_state32[__i + __pos1_32 - _M_nstate32], __r1, __r2); 00304 __r1 = __r2; 00305 __r2 = &_M_state32[__i]; 00306 } 00307 00308 _M_pos = 0; 00309 } 00310 00311 #endif 00312 00313 #ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_OPERATOREQUAL 00314 template<typename _UIntType, size_t __m, 00315 size_t __pos1, size_t __sl1, size_t __sl2, 00316 size_t __sr1, size_t __sr2, 00317 uint32_t __msk1, uint32_t __msk2, 00318 uint32_t __msk3, uint32_t __msk4, 00319 uint32_t __parity1, uint32_t __parity2, 00320 uint32_t __parity3, uint32_t __parity4> 00321 bool 00322 operator==(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 00323 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 00324 __msk1, __msk2, __msk3, __msk4, 00325 __parity1, __parity2, __parity3, __parity4>& __lhs, 00326 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 00327 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 00328 __msk1, __msk2, __msk3, __msk4, 00329 __parity1, __parity2, __parity3, __parity4>& __rhs) 00330 { 00331 typedef __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 00332 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 00333 __msk1, __msk2, __msk3, __msk4, 00334 __parity1, __parity2, __parity3, __parity4> __engine; 00335 return (std::equal(__lhs._M_stateT, 00336 __lhs._M_stateT + __engine::state_size, 00337 __rhs._M_stateT) 00338 && __lhs._M_pos == __rhs._M_pos); 00339 } 00340 #endif 00341 00342 template<typename _UIntType, size_t __m, 00343 size_t __pos1, size_t __sl1, size_t __sl2, 00344 size_t __sr1, size_t __sr2, 00345 uint32_t __msk1, uint32_t __msk2, 00346 uint32_t __msk3, uint32_t __msk4, 00347 uint32_t __parity1, uint32_t __parity2, 00348 uint32_t __parity3, uint32_t __parity4, 00349 typename _CharT, typename _Traits> 00350 std::basic_ostream<_CharT, _Traits>& 00351 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00352 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 00353 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 00354 __msk1, __msk2, __msk3, __msk4, 00355 __parity1, __parity2, __parity3, __parity4>& __x) 00356 { 00357 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00358 typedef typename __ostream_type::ios_base __ios_base; 00359 00360 const typename __ios_base::fmtflags __flags = __os.flags(); 00361 const _CharT __fill = __os.fill(); 00362 const _CharT __space = __os.widen(' '); 00363 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00364 __os.fill(__space); 00365 00366 for (size_t __i = 0; __i < __x._M_nstate32; ++__i) 00367 __os << __x._M_state32[__i] << __space; 00368 __os << __x._M_pos; 00369 00370 __os.flags(__flags); 00371 __os.fill(__fill); 00372 return __os; 00373 } 00374 00375 00376 template<typename _UIntType, size_t __m, 00377 size_t __pos1, size_t __sl1, size_t __sl2, 00378 size_t __sr1, size_t __sr2, 00379 uint32_t __msk1, uint32_t __msk2, 00380 uint32_t __msk3, uint32_t __msk4, 00381 uint32_t __parity1, uint32_t __parity2, 00382 uint32_t __parity3, uint32_t __parity4, 00383 typename _CharT, typename _Traits> 00384 std::basic_istream<_CharT, _Traits>& 00385 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00386 __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType, 00387 __m, __pos1, __sl1, __sl2, __sr1, __sr2, 00388 __msk1, __msk2, __msk3, __msk4, 00389 __parity1, __parity2, __parity3, __parity4>& __x) 00390 { 00391 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00392 typedef typename __istream_type::ios_base __ios_base; 00393 00394 const typename __ios_base::fmtflags __flags = __is.flags(); 00395 __is.flags(__ios_base::dec | __ios_base::skipws); 00396 00397 for (size_t __i = 0; __i < __x._M_nstate32; ++__i) 00398 __is >> __x._M_state32[__i]; 00399 __is >> __x._M_pos; 00400 00401 __is.flags(__flags); 00402 return __is; 00403 } 00404 00405 #endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ 00406 00407 /** 00408 * Iteration method due to M.D. J<o:>hnk. 00409 * 00410 * M.D. J<o:>hnk, Erzeugung von betaverteilten und gammaverteilten 00411 * Zufallszahlen, Metrika, Volume 8, 1964 00412 */ 00413 template<typename _RealType> 00414 template<typename _UniformRandomNumberGenerator> 00415 typename beta_distribution<_RealType>::result_type 00416 beta_distribution<_RealType>:: 00417 operator()(_UniformRandomNumberGenerator& __urng, 00418 const param_type& __param) 00419 { 00420 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 00421 __aurng(__urng); 00422 00423 result_type __x, __y; 00424 do 00425 { 00426 __x = std::exp(std::log(__aurng()) / __param.alpha()); 00427 __y = std::exp(std::log(__aurng()) / __param.beta()); 00428 } 00429 while (__x + __y > result_type(1)); 00430 00431 return __x / (__x + __y); 00432 } 00433 00434 template<typename _RealType> 00435 template<typename _OutputIterator, 00436 typename _UniformRandomNumberGenerator> 00437 void 00438 beta_distribution<_RealType>:: 00439 __generate_impl(_OutputIterator __f, _OutputIterator __t, 00440 _UniformRandomNumberGenerator& __urng, 00441 const param_type& __param) 00442 { 00443 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 00444 00445 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 00446 __aurng(__urng); 00447 00448 while (__f != __t) 00449 { 00450 result_type __x, __y; 00451 do 00452 { 00453 __x = std::exp(std::log(__aurng()) / __param.alpha()); 00454 __y = std::exp(std::log(__aurng()) / __param.beta()); 00455 } 00456 while (__x + __y > result_type(1)); 00457 00458 *__f++ = __x / (__x + __y); 00459 } 00460 } 00461 00462 template<typename _RealType, typename _CharT, typename _Traits> 00463 std::basic_ostream<_CharT, _Traits>& 00464 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00465 const __gnu_cxx::beta_distribution<_RealType>& __x) 00466 { 00467 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00468 typedef typename __ostream_type::ios_base __ios_base; 00469 00470 const typename __ios_base::fmtflags __flags = __os.flags(); 00471 const _CharT __fill = __os.fill(); 00472 const std::streamsize __precision = __os.precision(); 00473 const _CharT __space = __os.widen(' '); 00474 __os.flags(__ios_base::scientific | __ios_base::left); 00475 __os.fill(__space); 00476 __os.precision(std::numeric_limits<_RealType>::max_digits10); 00477 00478 __os << __x.alpha() << __space << __x.beta(); 00479 00480 __os.flags(__flags); 00481 __os.fill(__fill); 00482 __os.precision(__precision); 00483 return __os; 00484 } 00485 00486 template<typename _RealType, typename _CharT, typename _Traits> 00487 std::basic_istream<_CharT, _Traits>& 00488 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00489 __gnu_cxx::beta_distribution<_RealType>& __x) 00490 { 00491 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00492 typedef typename __istream_type::ios_base __ios_base; 00493 00494 const typename __ios_base::fmtflags __flags = __is.flags(); 00495 __is.flags(__ios_base::dec | __ios_base::skipws); 00496 00497 _RealType __alpha_val, __beta_val; 00498 __is >> __alpha_val >> __beta_val; 00499 __x.param(typename __gnu_cxx::beta_distribution<_RealType>:: 00500 param_type(__alpha_val, __beta_val)); 00501 00502 __is.flags(__flags); 00503 return __is; 00504 } 00505 00506 00507 template<std::size_t _Dimen, typename _RealType> 00508 template<typename _InputIterator1, typename _InputIterator2> 00509 void 00510 normal_mv_distribution<_Dimen, _RealType>::param_type:: 00511 _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 00512 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend) 00513 { 00514 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 00515 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 00516 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 00517 _M_mean.end(), _RealType(0)); 00518 00519 // Perform the Cholesky decomposition 00520 auto __w = _M_t.begin(); 00521 for (size_t __j = 0; __j < _Dimen; ++__j) 00522 { 00523 _RealType __sum = _RealType(0); 00524 00525 auto __slitbegin = __w; 00526 auto __cit = _M_t.begin(); 00527 for (size_t __i = 0; __i < __j; ++__i) 00528 { 00529 auto __slit = __slitbegin; 00530 _RealType __s = *__varcovbegin++; 00531 for (size_t __k = 0; __k < __i; ++__k) 00532 __s -= *__slit++ * *__cit++; 00533 00534 *__w++ = __s /= *__cit++; 00535 __sum += __s * __s; 00536 } 00537 00538 __sum = *__varcovbegin - __sum; 00539 if (__builtin_expect(__sum <= _RealType(0), 0)) 00540 std::__throw_runtime_error(__N("normal_mv_distribution::" 00541 "param_type::_M_init_full")); 00542 *__w++ = std::sqrt(__sum); 00543 00544 std::advance(__varcovbegin, _Dimen - __j); 00545 } 00546 } 00547 00548 template<std::size_t _Dimen, typename _RealType> 00549 template<typename _InputIterator1, typename _InputIterator2> 00550 void 00551 normal_mv_distribution<_Dimen, _RealType>::param_type:: 00552 _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 00553 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend) 00554 { 00555 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 00556 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 00557 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 00558 _M_mean.end(), _RealType(0)); 00559 00560 // Perform the Cholesky decomposition 00561 auto __w = _M_t.begin(); 00562 for (size_t __j = 0; __j < _Dimen; ++__j) 00563 { 00564 _RealType __sum = _RealType(0); 00565 00566 auto __slitbegin = __w; 00567 auto __cit = _M_t.begin(); 00568 for (size_t __i = 0; __i < __j; ++__i) 00569 { 00570 auto __slit = __slitbegin; 00571 _RealType __s = *__varcovbegin++; 00572 for (size_t __k = 0; __k < __i; ++__k) 00573 __s -= *__slit++ * *__cit++; 00574 00575 *__w++ = __s /= *__cit++; 00576 __sum += __s * __s; 00577 } 00578 00579 __sum = *__varcovbegin++ - __sum; 00580 if (__builtin_expect(__sum <= _RealType(0), 0)) 00581 std::__throw_runtime_error(__N("normal_mv_distribution::" 00582 "param_type::_M_init_full")); 00583 *__w++ = std::sqrt(__sum); 00584 } 00585 } 00586 00587 template<std::size_t _Dimen, typename _RealType> 00588 template<typename _InputIterator1, typename _InputIterator2> 00589 void 00590 normal_mv_distribution<_Dimen, _RealType>::param_type:: 00591 _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend, 00592 _InputIterator2 __varbegin, _InputIterator2 __varend) 00593 { 00594 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>) 00595 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>) 00596 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()), 00597 _M_mean.end(), _RealType(0)); 00598 00599 auto __w = _M_t.begin(); 00600 size_t __step = 0; 00601 while (__varbegin != __varend) 00602 { 00603 std::fill_n(__w, __step, _RealType(0)); 00604 __w += __step++; 00605 if (__builtin_expect(*__varbegin < _RealType(0), 0)) 00606 std::__throw_runtime_error(__N("normal_mv_distribution::" 00607 "param_type::_M_init_diagonal")); 00608 *__w++ = std::sqrt(*__varbegin++); 00609 } 00610 } 00611 00612 template<std::size_t _Dimen, typename _RealType> 00613 template<typename _UniformRandomNumberGenerator> 00614 typename normal_mv_distribution<_Dimen, _RealType>::result_type 00615 normal_mv_distribution<_Dimen, _RealType>:: 00616 operator()(_UniformRandomNumberGenerator& __urng, 00617 const param_type& __param) 00618 { 00619 result_type __ret; 00620 00621 _M_nd.__generate(__ret.begin(), __ret.end(), __urng); 00622 00623 auto __t_it = __param._M_t.crbegin(); 00624 for (size_t __i = _Dimen; __i > 0; --__i) 00625 { 00626 _RealType __sum = _RealType(0); 00627 for (size_t __j = __i; __j > 0; --__j) 00628 __sum += __ret[__j - 1] * *__t_it++; 00629 __ret[__i - 1] = __sum; 00630 } 00631 00632 return __ret; 00633 } 00634 00635 template<std::size_t _Dimen, typename _RealType> 00636 template<typename _ForwardIterator, typename _UniformRandomNumberGenerator> 00637 void 00638 normal_mv_distribution<_Dimen, _RealType>:: 00639 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 00640 _UniformRandomNumberGenerator& __urng, 00641 const param_type& __param) 00642 { 00643 __glibcxx_function_requires(_Mutable_ForwardIteratorConcept< 00644 _ForwardIterator>) 00645 while (__f != __t) 00646 *__f++ = this->operator()(__urng, __param); 00647 } 00648 00649 template<size_t _Dimen, typename _RealType> 00650 bool 00651 operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 00652 __d1, 00653 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& 00654 __d2) 00655 { 00656 return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd; 00657 } 00658 00659 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits> 00660 std::basic_ostream<_CharT, _Traits>& 00661 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00662 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x) 00663 { 00664 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00665 typedef typename __ostream_type::ios_base __ios_base; 00666 00667 const typename __ios_base::fmtflags __flags = __os.flags(); 00668 const _CharT __fill = __os.fill(); 00669 const std::streamsize __precision = __os.precision(); 00670 const _CharT __space = __os.widen(' '); 00671 __os.flags(__ios_base::scientific | __ios_base::left); 00672 __os.fill(__space); 00673 __os.precision(std::numeric_limits<_RealType>::max_digits10); 00674 00675 auto __mean = __x._M_param.mean(); 00676 for (auto __it : __mean) 00677 __os << __it << __space; 00678 auto __t = __x._M_param.varcov(); 00679 for (auto __it : __t) 00680 __os << __it << __space; 00681 00682 __os << __x._M_nd; 00683 00684 __os.flags(__flags); 00685 __os.fill(__fill); 00686 __os.precision(__precision); 00687 return __os; 00688 } 00689 00690 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits> 00691 std::basic_istream<_CharT, _Traits>& 00692 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00693 __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x) 00694 { 00695 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00696 typedef typename __istream_type::ios_base __ios_base; 00697 00698 const typename __ios_base::fmtflags __flags = __is.flags(); 00699 __is.flags(__ios_base::dec | __ios_base::skipws); 00700 00701 std::array<_RealType, _Dimen> __mean; 00702 for (auto& __it : __mean) 00703 __is >> __it; 00704 std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov; 00705 for (auto& __it : __varcov) 00706 __is >> __it; 00707 00708 __is >> __x._M_nd; 00709 00710 __x.param(typename normal_mv_distribution<_Dimen, _RealType>:: 00711 param_type(__mean.begin(), __mean.end(), 00712 __varcov.begin(), __varcov.end())); 00713 00714 __is.flags(__flags); 00715 return __is; 00716 } 00717 00718 00719 template<typename _RealType> 00720 template<typename _OutputIterator, 00721 typename _UniformRandomNumberGenerator> 00722 void 00723 rice_distribution<_RealType>:: 00724 __generate_impl(_OutputIterator __f, _OutputIterator __t, 00725 _UniformRandomNumberGenerator& __urng, 00726 const param_type& __p) 00727 { 00728 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 00729 00730 while (__f != __t) 00731 { 00732 typename std::normal_distribution<result_type>::param_type 00733 __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma()); 00734 result_type __x = this->_M_ndx(__px, __urng); 00735 result_type __y = this->_M_ndy(__py, __urng); 00736 #if _GLIBCXX_USE_C99_MATH_TR1 00737 *__f++ = std::hypot(__x, __y); 00738 #else 00739 *__f++ = std::sqrt(__x * __x + __y * __y); 00740 #endif 00741 } 00742 } 00743 00744 template<typename _RealType, typename _CharT, typename _Traits> 00745 std::basic_ostream<_CharT, _Traits>& 00746 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00747 const rice_distribution<_RealType>& __x) 00748 { 00749 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00750 typedef typename __ostream_type::ios_base __ios_base; 00751 00752 const typename __ios_base::fmtflags __flags = __os.flags(); 00753 const _CharT __fill = __os.fill(); 00754 const std::streamsize __precision = __os.precision(); 00755 const _CharT __space = __os.widen(' '); 00756 __os.flags(__ios_base::scientific | __ios_base::left); 00757 __os.fill(__space); 00758 __os.precision(std::numeric_limits<_RealType>::max_digits10); 00759 00760 __os << __x.nu() << __space << __x.sigma(); 00761 __os << __space << __x._M_ndx; 00762 __os << __space << __x._M_ndy; 00763 00764 __os.flags(__flags); 00765 __os.fill(__fill); 00766 __os.precision(__precision); 00767 return __os; 00768 } 00769 00770 template<typename _RealType, typename _CharT, typename _Traits> 00771 std::basic_istream<_CharT, _Traits>& 00772 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00773 rice_distribution<_RealType>& __x) 00774 { 00775 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00776 typedef typename __istream_type::ios_base __ios_base; 00777 00778 const typename __ios_base::fmtflags __flags = __is.flags(); 00779 __is.flags(__ios_base::dec | __ios_base::skipws); 00780 00781 _RealType __nu_val, __sigma_val; 00782 __is >> __nu_val >> __sigma_val; 00783 __is >> __x._M_ndx; 00784 __is >> __x._M_ndy; 00785 __x.param(typename rice_distribution<_RealType>:: 00786 param_type(__nu_val, __sigma_val)); 00787 00788 __is.flags(__flags); 00789 return __is; 00790 } 00791 00792 00793 template<typename _RealType> 00794 template<typename _OutputIterator, 00795 typename _UniformRandomNumberGenerator> 00796 void 00797 nakagami_distribution<_RealType>:: 00798 __generate_impl(_OutputIterator __f, _OutputIterator __t, 00799 _UniformRandomNumberGenerator& __urng, 00800 const param_type& __p) 00801 { 00802 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 00803 00804 typename std::gamma_distribution<result_type>::param_type 00805 __pg(__p.mu(), __p.omega() / __p.mu()); 00806 while (__f != __t) 00807 *__f++ = std::sqrt(this->_M_gd(__pg, __urng)); 00808 } 00809 00810 template<typename _RealType, typename _CharT, typename _Traits> 00811 std::basic_ostream<_CharT, _Traits>& 00812 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00813 const nakagami_distribution<_RealType>& __x) 00814 { 00815 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00816 typedef typename __ostream_type::ios_base __ios_base; 00817 00818 const typename __ios_base::fmtflags __flags = __os.flags(); 00819 const _CharT __fill = __os.fill(); 00820 const std::streamsize __precision = __os.precision(); 00821 const _CharT __space = __os.widen(' '); 00822 __os.flags(__ios_base::scientific | __ios_base::left); 00823 __os.fill(__space); 00824 __os.precision(std::numeric_limits<_RealType>::max_digits10); 00825 00826 __os << __x.mu() << __space << __x.omega(); 00827 __os << __space << __x._M_gd; 00828 00829 __os.flags(__flags); 00830 __os.fill(__fill); 00831 __os.precision(__precision); 00832 return __os; 00833 } 00834 00835 template<typename _RealType, typename _CharT, typename _Traits> 00836 std::basic_istream<_CharT, _Traits>& 00837 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00838 nakagami_distribution<_RealType>& __x) 00839 { 00840 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00841 typedef typename __istream_type::ios_base __ios_base; 00842 00843 const typename __ios_base::fmtflags __flags = __is.flags(); 00844 __is.flags(__ios_base::dec | __ios_base::skipws); 00845 00846 _RealType __mu_val, __omega_val; 00847 __is >> __mu_val >> __omega_val; 00848 __is >> __x._M_gd; 00849 __x.param(typename nakagami_distribution<_RealType>:: 00850 param_type(__mu_val, __omega_val)); 00851 00852 __is.flags(__flags); 00853 return __is; 00854 } 00855 00856 00857 template<typename _RealType> 00858 template<typename _OutputIterator, 00859 typename _UniformRandomNumberGenerator> 00860 void 00861 pareto_distribution<_RealType>:: 00862 __generate_impl(_OutputIterator __f, _OutputIterator __t, 00863 _UniformRandomNumberGenerator& __urng, 00864 const param_type& __p) 00865 { 00866 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 00867 00868 result_type __mu_val = __p.mu(); 00869 result_type __malphinv = -result_type(1) / __p.alpha(); 00870 while (__f != __t) 00871 *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv); 00872 } 00873 00874 template<typename _RealType, typename _CharT, typename _Traits> 00875 std::basic_ostream<_CharT, _Traits>& 00876 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00877 const pareto_distribution<_RealType>& __x) 00878 { 00879 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00880 typedef typename __ostream_type::ios_base __ios_base; 00881 00882 const typename __ios_base::fmtflags __flags = __os.flags(); 00883 const _CharT __fill = __os.fill(); 00884 const std::streamsize __precision = __os.precision(); 00885 const _CharT __space = __os.widen(' '); 00886 __os.flags(__ios_base::scientific | __ios_base::left); 00887 __os.fill(__space); 00888 __os.precision(std::numeric_limits<_RealType>::max_digits10); 00889 00890 __os << __x.alpha() << __space << __x.mu(); 00891 __os << __space << __x._M_ud; 00892 00893 __os.flags(__flags); 00894 __os.fill(__fill); 00895 __os.precision(__precision); 00896 return __os; 00897 } 00898 00899 template<typename _RealType, typename _CharT, typename _Traits> 00900 std::basic_istream<_CharT, _Traits>& 00901 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00902 pareto_distribution<_RealType>& __x) 00903 { 00904 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00905 typedef typename __istream_type::ios_base __ios_base; 00906 00907 const typename __ios_base::fmtflags __flags = __is.flags(); 00908 __is.flags(__ios_base::dec | __ios_base::skipws); 00909 00910 _RealType __alpha_val, __mu_val; 00911 __is >> __alpha_val >> __mu_val; 00912 __is >> __x._M_ud; 00913 __x.param(typename pareto_distribution<_RealType>:: 00914 param_type(__alpha_val, __mu_val)); 00915 00916 __is.flags(__flags); 00917 return __is; 00918 } 00919 00920 00921 template<typename _RealType> 00922 template<typename _UniformRandomNumberGenerator> 00923 typename k_distribution<_RealType>::result_type 00924 k_distribution<_RealType>:: 00925 operator()(_UniformRandomNumberGenerator& __urng) 00926 { 00927 result_type __x = this->_M_gd1(__urng); 00928 result_type __y = this->_M_gd2(__urng); 00929 return std::sqrt(__x * __y); 00930 } 00931 00932 template<typename _RealType> 00933 template<typename _UniformRandomNumberGenerator> 00934 typename k_distribution<_RealType>::result_type 00935 k_distribution<_RealType>:: 00936 operator()(_UniformRandomNumberGenerator& __urng, 00937 const param_type& __p) 00938 { 00939 typename std::gamma_distribution<result_type>::param_type 00940 __p1(__p.lambda(), result_type(1) / __p.lambda()), 00941 __p2(__p.nu(), __p.mu() / __p.nu()); 00942 result_type __x = this->_M_gd1(__p1, __urng); 00943 result_type __y = this->_M_gd2(__p2, __urng); 00944 return std::sqrt(__x * __y); 00945 } 00946 00947 template<typename _RealType> 00948 template<typename _OutputIterator, 00949 typename _UniformRandomNumberGenerator> 00950 void 00951 k_distribution<_RealType>:: 00952 __generate_impl(_OutputIterator __f, _OutputIterator __t, 00953 _UniformRandomNumberGenerator& __urng, 00954 const param_type& __p) 00955 { 00956 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 00957 00958 typename std::gamma_distribution<result_type>::param_type 00959 __p1(__p.lambda(), result_type(1) / __p.lambda()), 00960 __p2(__p.nu(), __p.mu() / __p.nu()); 00961 while (__f != __t) 00962 { 00963 result_type __x = this->_M_gd1(__p1, __urng); 00964 result_type __y = this->_M_gd2(__p2, __urng); 00965 *__f++ = std::sqrt(__x * __y); 00966 } 00967 } 00968 00969 template<typename _RealType, typename _CharT, typename _Traits> 00970 std::basic_ostream<_CharT, _Traits>& 00971 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00972 const k_distribution<_RealType>& __x) 00973 { 00974 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00975 typedef typename __ostream_type::ios_base __ios_base; 00976 00977 const typename __ios_base::fmtflags __flags = __os.flags(); 00978 const _CharT __fill = __os.fill(); 00979 const std::streamsize __precision = __os.precision(); 00980 const _CharT __space = __os.widen(' '); 00981 __os.flags(__ios_base::scientific | __ios_base::left); 00982 __os.fill(__space); 00983 __os.precision(std::numeric_limits<_RealType>::max_digits10); 00984 00985 __os << __x.lambda() << __space << __x.mu() << __space << __x.nu(); 00986 __os << __space << __x._M_gd1; 00987 __os << __space << __x._M_gd2; 00988 00989 __os.flags(__flags); 00990 __os.fill(__fill); 00991 __os.precision(__precision); 00992 return __os; 00993 } 00994 00995 template<typename _RealType, typename _CharT, typename _Traits> 00996 std::basic_istream<_CharT, _Traits>& 00997 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00998 k_distribution<_RealType>& __x) 00999 { 01000 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01001 typedef typename __istream_type::ios_base __ios_base; 01002 01003 const typename __ios_base::fmtflags __flags = __is.flags(); 01004 __is.flags(__ios_base::dec | __ios_base::skipws); 01005 01006 _RealType __lambda_val, __mu_val, __nu_val; 01007 __is >> __lambda_val >> __mu_val >> __nu_val; 01008 __is >> __x._M_gd1; 01009 __is >> __x._M_gd2; 01010 __x.param(typename k_distribution<_RealType>:: 01011 param_type(__lambda_val, __mu_val, __nu_val)); 01012 01013 __is.flags(__flags); 01014 return __is; 01015 } 01016 01017 01018 template<typename _RealType> 01019 template<typename _OutputIterator, 01020 typename _UniformRandomNumberGenerator> 01021 void 01022 arcsine_distribution<_RealType>:: 01023 __generate_impl(_OutputIterator __f, _OutputIterator __t, 01024 _UniformRandomNumberGenerator& __urng, 01025 const param_type& __p) 01026 { 01027 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 01028 01029 result_type __dif = __p.b() - __p.a(); 01030 result_type __sum = __p.a() + __p.b(); 01031 while (__f != __t) 01032 { 01033 result_type __x = std::sin(this->_M_ud(__urng)); 01034 *__f++ = (__x * __dif + __sum) / result_type(2); 01035 } 01036 } 01037 01038 template<typename _RealType, typename _CharT, typename _Traits> 01039 std::basic_ostream<_CharT, _Traits>& 01040 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01041 const arcsine_distribution<_RealType>& __x) 01042 { 01043 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01044 typedef typename __ostream_type::ios_base __ios_base; 01045 01046 const typename __ios_base::fmtflags __flags = __os.flags(); 01047 const _CharT __fill = __os.fill(); 01048 const std::streamsize __precision = __os.precision(); 01049 const _CharT __space = __os.widen(' '); 01050 __os.flags(__ios_base::scientific | __ios_base::left); 01051 __os.fill(__space); 01052 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01053 01054 __os << __x.a() << __space << __x.b(); 01055 __os << __space << __x._M_ud; 01056 01057 __os.flags(__flags); 01058 __os.fill(__fill); 01059 __os.precision(__precision); 01060 return __os; 01061 } 01062 01063 template<typename _RealType, typename _CharT, typename _Traits> 01064 std::basic_istream<_CharT, _Traits>& 01065 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01066 arcsine_distribution<_RealType>& __x) 01067 { 01068 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01069 typedef typename __istream_type::ios_base __ios_base; 01070 01071 const typename __ios_base::fmtflags __flags = __is.flags(); 01072 __is.flags(__ios_base::dec | __ios_base::skipws); 01073 01074 _RealType __a, __b; 01075 __is >> __a >> __b; 01076 __is >> __x._M_ud; 01077 __x.param(typename arcsine_distribution<_RealType>:: 01078 param_type(__a, __b)); 01079 01080 __is.flags(__flags); 01081 return __is; 01082 } 01083 01084 01085 template<typename _RealType> 01086 template<typename _UniformRandomNumberGenerator> 01087 typename hoyt_distribution<_RealType>::result_type 01088 hoyt_distribution<_RealType>:: 01089 operator()(_UniformRandomNumberGenerator& __urng) 01090 { 01091 result_type __x = this->_M_ad(__urng); 01092 result_type __y = this->_M_ed(__urng); 01093 return (result_type(2) * this->q() 01094 / (result_type(1) + this->q() * this->q())) 01095 * std::sqrt(this->omega() * __x * __y); 01096 } 01097 01098 template<typename _RealType> 01099 template<typename _UniformRandomNumberGenerator> 01100 typename hoyt_distribution<_RealType>::result_type 01101 hoyt_distribution<_RealType>:: 01102 operator()(_UniformRandomNumberGenerator& __urng, 01103 const param_type& __p) 01104 { 01105 result_type __q2 = __p.q() * __p.q(); 01106 result_type __num = result_type(0.5L) * (result_type(1) + __q2); 01107 typename __gnu_cxx::arcsine_distribution<result_type>::param_type 01108 __pa(__num, __num / __q2); 01109 result_type __x = this->_M_ad(__pa, __urng); 01110 result_type __y = this->_M_ed(__urng); 01111 return (result_type(2) * __p.q() / (result_type(1) + __q2)) 01112 * std::sqrt(__p.omega() * __x * __y); 01113 } 01114 01115 template<typename _RealType> 01116 template<typename _OutputIterator, 01117 typename _UniformRandomNumberGenerator> 01118 void 01119 hoyt_distribution<_RealType>:: 01120 __generate_impl(_OutputIterator __f, _OutputIterator __t, 01121 _UniformRandomNumberGenerator& __urng, 01122 const param_type& __p) 01123 { 01124 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 01125 01126 result_type __2q = result_type(2) * __p.q(); 01127 result_type __q2 = __p.q() * __p.q(); 01128 result_type __q2p1 = result_type(1) + __q2; 01129 result_type __num = result_type(0.5L) * __q2p1; 01130 result_type __omega = __p.omega(); 01131 typename __gnu_cxx::arcsine_distribution<result_type>::param_type 01132 __pa(__num, __num / __q2); 01133 while (__f != __t) 01134 { 01135 result_type __x = this->_M_ad(__pa, __urng); 01136 result_type __y = this->_M_ed(__urng); 01137 *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y); 01138 } 01139 } 01140 01141 template<typename _RealType, typename _CharT, typename _Traits> 01142 std::basic_ostream<_CharT, _Traits>& 01143 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01144 const hoyt_distribution<_RealType>& __x) 01145 { 01146 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01147 typedef typename __ostream_type::ios_base __ios_base; 01148 01149 const typename __ios_base::fmtflags __flags = __os.flags(); 01150 const _CharT __fill = __os.fill(); 01151 const std::streamsize __precision = __os.precision(); 01152 const _CharT __space = __os.widen(' '); 01153 __os.flags(__ios_base::scientific | __ios_base::left); 01154 __os.fill(__space); 01155 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01156 01157 __os << __x.q() << __space << __x.omega(); 01158 __os << __space << __x._M_ad; 01159 __os << __space << __x._M_ed; 01160 01161 __os.flags(__flags); 01162 __os.fill(__fill); 01163 __os.precision(__precision); 01164 return __os; 01165 } 01166 01167 template<typename _RealType, typename _CharT, typename _Traits> 01168 std::basic_istream<_CharT, _Traits>& 01169 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01170 hoyt_distribution<_RealType>& __x) 01171 { 01172 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01173 typedef typename __istream_type::ios_base __ios_base; 01174 01175 const typename __ios_base::fmtflags __flags = __is.flags(); 01176 __is.flags(__ios_base::dec | __ios_base::skipws); 01177 01178 _RealType __q, __omega; 01179 __is >> __q >> __omega; 01180 __is >> __x._M_ad; 01181 __is >> __x._M_ed; 01182 __x.param(typename hoyt_distribution<_RealType>:: 01183 param_type(__q, __omega)); 01184 01185 __is.flags(__flags); 01186 return __is; 01187 } 01188 01189 01190 template<typename _RealType> 01191 template<typename _OutputIterator, 01192 typename _UniformRandomNumberGenerator> 01193 void 01194 triangular_distribution<_RealType>:: 01195 __generate_impl(_OutputIterator __f, _OutputIterator __t, 01196 _UniformRandomNumberGenerator& __urng, 01197 const param_type& __param) 01198 { 01199 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 01200 01201 while (__f != __t) 01202 *__f++ = this->operator()(__urng, __param); 01203 } 01204 01205 template<typename _RealType, typename _CharT, typename _Traits> 01206 std::basic_ostream<_CharT, _Traits>& 01207 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01208 const __gnu_cxx::triangular_distribution<_RealType>& __x) 01209 { 01210 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01211 typedef typename __ostream_type::ios_base __ios_base; 01212 01213 const typename __ios_base::fmtflags __flags = __os.flags(); 01214 const _CharT __fill = __os.fill(); 01215 const std::streamsize __precision = __os.precision(); 01216 const _CharT __space = __os.widen(' '); 01217 __os.flags(__ios_base::scientific | __ios_base::left); 01218 __os.fill(__space); 01219 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01220 01221 __os << __x.a() << __space << __x.b() << __space << __x.c(); 01222 01223 __os.flags(__flags); 01224 __os.fill(__fill); 01225 __os.precision(__precision); 01226 return __os; 01227 } 01228 01229 template<typename _RealType, typename _CharT, typename _Traits> 01230 std::basic_istream<_CharT, _Traits>& 01231 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01232 __gnu_cxx::triangular_distribution<_RealType>& __x) 01233 { 01234 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01235 typedef typename __istream_type::ios_base __ios_base; 01236 01237 const typename __ios_base::fmtflags __flags = __is.flags(); 01238 __is.flags(__ios_base::dec | __ios_base::skipws); 01239 01240 _RealType __a, __b, __c; 01241 __is >> __a >> __b >> __c; 01242 __x.param(typename __gnu_cxx::triangular_distribution<_RealType>:: 01243 param_type(__a, __b, __c)); 01244 01245 __is.flags(__flags); 01246 return __is; 01247 } 01248 01249 01250 template<typename _RealType> 01251 template<typename _UniformRandomNumberGenerator> 01252 typename von_mises_distribution<_RealType>::result_type 01253 von_mises_distribution<_RealType>:: 01254 operator()(_UniformRandomNumberGenerator& __urng, 01255 const param_type& __p) 01256 { 01257 const result_type __pi 01258 = __gnu_cxx::__math_constants<result_type>::__pi; 01259 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01260 __aurng(__urng); 01261 01262 result_type __f; 01263 while (1) 01264 { 01265 result_type __rnd = std::cos(__pi * __aurng()); 01266 __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd); 01267 result_type __c = __p._M_kappa * (__p._M_r - __f); 01268 01269 result_type __rnd2 = __aurng(); 01270 if (__c * (result_type(2) - __c) > __rnd2) 01271 break; 01272 if (std::log(__c / __rnd2) >= __c - result_type(1)) 01273 break; 01274 } 01275 01276 result_type __res = std::acos(__f); 01277 #if _GLIBCXX_USE_C99_MATH_TR1 01278 __res = std::copysign(__res, __aurng() - result_type(0.5)); 01279 #else 01280 if (__aurng() < result_type(0.5)) 01281 __res = -__res; 01282 #endif 01283 __res += __p._M_mu; 01284 if (__res > __pi) 01285 __res -= result_type(2) * __pi; 01286 else if (__res < -__pi) 01287 __res += result_type(2) * __pi; 01288 return __res; 01289 } 01290 01291 template<typename _RealType> 01292 template<typename _OutputIterator, 01293 typename _UniformRandomNumberGenerator> 01294 void 01295 von_mises_distribution<_RealType>:: 01296 __generate_impl(_OutputIterator __f, _OutputIterator __t, 01297 _UniformRandomNumberGenerator& __urng, 01298 const param_type& __param) 01299 { 01300 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 01301 01302 while (__f != __t) 01303 *__f++ = this->operator()(__urng, __param); 01304 } 01305 01306 template<typename _RealType, typename _CharT, typename _Traits> 01307 std::basic_ostream<_CharT, _Traits>& 01308 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01309 const __gnu_cxx::von_mises_distribution<_RealType>& __x) 01310 { 01311 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01312 typedef typename __ostream_type::ios_base __ios_base; 01313 01314 const typename __ios_base::fmtflags __flags = __os.flags(); 01315 const _CharT __fill = __os.fill(); 01316 const std::streamsize __precision = __os.precision(); 01317 const _CharT __space = __os.widen(' '); 01318 __os.flags(__ios_base::scientific | __ios_base::left); 01319 __os.fill(__space); 01320 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01321 01322 __os << __x.mu() << __space << __x.kappa(); 01323 01324 __os.flags(__flags); 01325 __os.fill(__fill); 01326 __os.precision(__precision); 01327 return __os; 01328 } 01329 01330 template<typename _RealType, typename _CharT, typename _Traits> 01331 std::basic_istream<_CharT, _Traits>& 01332 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01333 __gnu_cxx::von_mises_distribution<_RealType>& __x) 01334 { 01335 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01336 typedef typename __istream_type::ios_base __ios_base; 01337 01338 const typename __ios_base::fmtflags __flags = __is.flags(); 01339 __is.flags(__ios_base::dec | __ios_base::skipws); 01340 01341 _RealType __mu, __kappa; 01342 __is >> __mu >> __kappa; 01343 __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>:: 01344 param_type(__mu, __kappa)); 01345 01346 __is.flags(__flags); 01347 return __is; 01348 } 01349 01350 01351 template<typename _UIntType> 01352 template<typename _UniformRandomNumberGenerator> 01353 typename hypergeometric_distribution<_UIntType>::result_type 01354 hypergeometric_distribution<_UIntType>:: 01355 operator()(_UniformRandomNumberGenerator& __urng, 01356 const param_type& __param) 01357 { 01358 std::__detail::_Adaptor<_UniformRandomNumberGenerator, double> 01359 __aurng(__urng); 01360 01361 result_type __a = __param.successful_size(); 01362 result_type __b = __param.total_size(); 01363 result_type __k = 0; 01364 01365 if (__param.total_draws() < __param.total_size() / 2) 01366 { 01367 for (result_type __i = 0; __i < __param.total_draws(); ++__i) 01368 { 01369 if (__b * __aurng() < __a) 01370 { 01371 ++__k; 01372 if (__k == __param.successful_size()) 01373 return __k; 01374 --__a; 01375 } 01376 --__b; 01377 } 01378 return __k; 01379 } 01380 else 01381 { 01382 for (result_type __i = 0; __i < __param.unsuccessful_size(); ++__i) 01383 { 01384 if (__b * __aurng() < __a) 01385 { 01386 ++__k; 01387 if (__k == __param.successful_size()) 01388 return __param.successful_size() - __k; 01389 --__a; 01390 } 01391 --__b; 01392 } 01393 return __param.successful_size() - __k; 01394 } 01395 } 01396 01397 template<typename _UIntType> 01398 template<typename _OutputIterator, 01399 typename _UniformRandomNumberGenerator> 01400 void 01401 hypergeometric_distribution<_UIntType>:: 01402 __generate_impl(_OutputIterator __f, _OutputIterator __t, 01403 _UniformRandomNumberGenerator& __urng, 01404 const param_type& __param) 01405 { 01406 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 01407 01408 while (__f != __t) 01409 *__f++ = this->operator()(__urng); 01410 } 01411 01412 template<typename _UIntType, typename _CharT, typename _Traits> 01413 std::basic_ostream<_CharT, _Traits>& 01414 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01415 const __gnu_cxx::hypergeometric_distribution<_UIntType>& __x) 01416 { 01417 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01418 typedef typename __ostream_type::ios_base __ios_base; 01419 01420 const typename __ios_base::fmtflags __flags = __os.flags(); 01421 const _CharT __fill = __os.fill(); 01422 const std::streamsize __precision = __os.precision(); 01423 const _CharT __space = __os.widen(' '); 01424 __os.flags(__ios_base::scientific | __ios_base::left); 01425 __os.fill(__space); 01426 __os.precision(std::numeric_limits<_UIntType>::max_digits10); 01427 01428 __os << __x.total_size() << __space << __x.successful_size() << __space 01429 << __x.total_draws(); 01430 01431 __os.flags(__flags); 01432 __os.fill(__fill); 01433 __os.precision(__precision); 01434 return __os; 01435 } 01436 01437 template<typename _UIntType, typename _CharT, typename _Traits> 01438 std::basic_istream<_CharT, _Traits>& 01439 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01440 __gnu_cxx::hypergeometric_distribution<_UIntType>& __x) 01441 { 01442 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01443 typedef typename __istream_type::ios_base __ios_base; 01444 01445 const typename __ios_base::fmtflags __flags = __is.flags(); 01446 __is.flags(__ios_base::dec | __ios_base::skipws); 01447 01448 _UIntType __total_size, __successful_size, __total_draws; 01449 __is >> __total_size >> __successful_size >> __total_draws; 01450 __x.param(typename __gnu_cxx::hypergeometric_distribution<_UIntType>:: 01451 param_type(__total_size, __successful_size, __total_draws)); 01452 01453 __is.flags(__flags); 01454 return __is; 01455 } 01456 01457 01458 template<typename _RealType> 01459 template<typename _UniformRandomNumberGenerator> 01460 typename logistic_distribution<_RealType>::result_type 01461 logistic_distribution<_RealType>:: 01462 operator()(_UniformRandomNumberGenerator& __urng, 01463 const param_type& __p) 01464 { 01465 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01466 __aurng(__urng); 01467 01468 result_type __arg = result_type(1); 01469 while (__arg == result_type(1) || __arg == result_type(0)) 01470 __arg = __aurng(); 01471 return __p.a() 01472 + __p.b() * std::log(__arg / (result_type(1) - __arg)); 01473 } 01474 01475 template<typename _RealType> 01476 template<typename _OutputIterator, 01477 typename _UniformRandomNumberGenerator> 01478 void 01479 logistic_distribution<_RealType>:: 01480 __generate_impl(_OutputIterator __f, _OutputIterator __t, 01481 _UniformRandomNumberGenerator& __urng, 01482 const param_type& __p) 01483 { 01484 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 01485 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01486 __aurng(__urng); 01487 01488 while (__f != __t) 01489 { 01490 result_type __arg = result_type(1); 01491 while (__arg == result_type(1) || __arg == result_type(0)) 01492 __arg = __aurng(); 01493 *__f++ = __p.a() 01494 + __p.b() * std::log(__arg / (result_type(1) - __arg)); 01495 } 01496 } 01497 01498 template<typename _RealType, typename _CharT, typename _Traits> 01499 std::basic_ostream<_CharT, _Traits>& 01500 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01501 const logistic_distribution<_RealType>& __x) 01502 { 01503 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01504 typedef typename __ostream_type::ios_base __ios_base; 01505 01506 const typename __ios_base::fmtflags __flags = __os.flags(); 01507 const _CharT __fill = __os.fill(); 01508 const std::streamsize __precision = __os.precision(); 01509 const _CharT __space = __os.widen(' '); 01510 __os.flags(__ios_base::scientific | __ios_base::left); 01511 __os.fill(__space); 01512 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01513 01514 __os << __x.a() << __space << __x.b(); 01515 01516 __os.flags(__flags); 01517 __os.fill(__fill); 01518 __os.precision(__precision); 01519 return __os; 01520 } 01521 01522 template<typename _RealType, typename _CharT, typename _Traits> 01523 std::basic_istream<_CharT, _Traits>& 01524 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01525 logistic_distribution<_RealType>& __x) 01526 { 01527 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01528 typedef typename __istream_type::ios_base __ios_base; 01529 01530 const typename __ios_base::fmtflags __flags = __is.flags(); 01531 __is.flags(__ios_base::dec | __ios_base::skipws); 01532 01533 _RealType __a, __b; 01534 __is >> __a >> __b; 01535 __x.param(typename logistic_distribution<_RealType>:: 01536 param_type(__a, __b)); 01537 01538 __is.flags(__flags); 01539 return __is; 01540 } 01541 01542 01543 namespace { 01544 01545 // Helper class for the uniform_on_sphere_distribution generation 01546 // function. 01547 template<std::size_t _Dimen, typename _RealType> 01548 class uniform_on_sphere_helper 01549 { 01550 typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>:: 01551 result_type result_type; 01552 01553 public: 01554 template<typename _NormalDistribution, 01555 typename _UniformRandomNumberGenerator> 01556 result_type operator()(_NormalDistribution& __nd, 01557 _UniformRandomNumberGenerator& __urng) 01558 { 01559 result_type __ret; 01560 typename result_type::value_type __norm; 01561 01562 do 01563 { 01564 auto __sum = _RealType(0); 01565 01566 std::generate(__ret.begin(), __ret.end(), 01567 [&__nd, &__urng, &__sum](){ 01568 _RealType __t = __nd(__urng); 01569 __sum += __t * __t; 01570 return __t; }); 01571 __norm = std::sqrt(__sum); 01572 } 01573 while (__norm == _RealType(0) || ! std::isfinite(__norm)); 01574 01575 std::transform(__ret.begin(), __ret.end(), __ret.begin(), 01576 [__norm](_RealType __val){ return __val / __norm; }); 01577 01578 return __ret; 01579 } 01580 }; 01581 01582 01583 template<typename _RealType> 01584 class uniform_on_sphere_helper<2, _RealType> 01585 { 01586 typedef typename uniform_on_sphere_distribution<2, _RealType>:: 01587 result_type result_type; 01588 01589 public: 01590 template<typename _NormalDistribution, 01591 typename _UniformRandomNumberGenerator> 01592 result_type operator()(_NormalDistribution&, 01593 _UniformRandomNumberGenerator& __urng) 01594 { 01595 result_type __ret; 01596 _RealType __sq; 01597 std::__detail::_Adaptor<_UniformRandomNumberGenerator, 01598 _RealType> __aurng(__urng); 01599 01600 do 01601 { 01602 __ret[0] = _RealType(2) * __aurng() - _RealType(1); 01603 __ret[1] = _RealType(2) * __aurng() - _RealType(1); 01604 01605 __sq = __ret[0] * __ret[0] + __ret[1] * __ret[1]; 01606 } 01607 while (__sq == _RealType(0) || __sq > _RealType(1)); 01608 01609 #if _GLIBCXX_USE_C99_MATH_TR1 01610 // Yes, we do not just use sqrt(__sq) because hypot() is more 01611 // accurate. 01612 auto __norm = std::hypot(__ret[0], __ret[1]); 01613 #else 01614 auto __norm = std::sqrt(__sq); 01615 #endif 01616 __ret[0] /= __norm; 01617 __ret[1] /= __norm; 01618 01619 return __ret; 01620 } 01621 }; 01622 01623 } 01624 01625 01626 template<std::size_t _Dimen, typename _RealType> 01627 template<typename _UniformRandomNumberGenerator> 01628 typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type 01629 uniform_on_sphere_distribution<_Dimen, _RealType>:: 01630 operator()(_UniformRandomNumberGenerator& __urng, 01631 const param_type& __p) 01632 { 01633 uniform_on_sphere_helper<_Dimen, _RealType> __helper; 01634 return __helper(_M_nd, __urng); 01635 } 01636 01637 template<std::size_t _Dimen, typename _RealType> 01638 template<typename _OutputIterator, 01639 typename _UniformRandomNumberGenerator> 01640 void 01641 uniform_on_sphere_distribution<_Dimen, _RealType>:: 01642 __generate_impl(_OutputIterator __f, _OutputIterator __t, 01643 _UniformRandomNumberGenerator& __urng, 01644 const param_type& __param) 01645 { 01646 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>) 01647 01648 while (__f != __t) 01649 *__f++ = this->operator()(__urng, __param); 01650 } 01651 01652 template<std::size_t _Dimen, typename _RealType, typename _CharT, 01653 typename _Traits> 01654 std::basic_ostream<_CharT, _Traits>& 01655 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01656 const __gnu_cxx::uniform_on_sphere_distribution<_Dimen, 01657 _RealType>& __x) 01658 { 01659 return __os << __x._M_nd; 01660 } 01661 01662 template<std::size_t _Dimen, typename _RealType, typename _CharT, 01663 typename _Traits> 01664 std::basic_istream<_CharT, _Traits>& 01665 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01666 __gnu_cxx::uniform_on_sphere_distribution<_Dimen, 01667 _RealType>& __x) 01668 { 01669 return __is >> __x._M_nd; 01670 } 01671 01672 _GLIBCXX_END_NAMESPACE_VERSION 01673 } // namespace 01674 01675 01676 #endif // _EXT_RANDOM_TCC