libstdc++
|
00001 // random number generation (out of line) -*- C++ -*- 00002 00003 // Copyright (C) 2009-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 bits/random.tcc 00026 * This is an internal header file, included by other library headers. 00027 * Do not attempt to use it directly. @headername{random} 00028 */ 00029 00030 #ifndef _RANDOM_TCC 00031 #define _RANDOM_TCC 1 00032 00033 #include <numeric> // std::accumulate and std::partial_sum 00034 00035 namespace std _GLIBCXX_VISIBILITY(default) 00036 { 00037 /* 00038 * (Further) implementation-space details. 00039 */ 00040 namespace __detail 00041 { 00042 _GLIBCXX_BEGIN_NAMESPACE_VERSION 00043 00044 // General case for x = (ax + c) mod m -- use Schrage's algorithm 00045 // to avoid integer overflow. 00046 // 00047 // Preconditions: a > 0, m > 0. 00048 // 00049 // Note: only works correctly for __m % __a < __m / __a. 00050 template<typename _Tp, _Tp __m, _Tp __a, _Tp __c> 00051 _Tp 00052 _Mod<_Tp, __m, __a, __c, false, true>:: 00053 __calc(_Tp __x) 00054 { 00055 if (__a == 1) 00056 __x %= __m; 00057 else 00058 { 00059 static const _Tp __q = __m / __a; 00060 static const _Tp __r = __m % __a; 00061 00062 _Tp __t1 = __a * (__x % __q); 00063 _Tp __t2 = __r * (__x / __q); 00064 if (__t1 >= __t2) 00065 __x = __t1 - __t2; 00066 else 00067 __x = __m - __t2 + __t1; 00068 } 00069 00070 if (__c != 0) 00071 { 00072 const _Tp __d = __m - __x; 00073 if (__d > __c) 00074 __x += __c; 00075 else 00076 __x = __c - __d; 00077 } 00078 return __x; 00079 } 00080 00081 template<typename _InputIterator, typename _OutputIterator, 00082 typename _Tp> 00083 _OutputIterator 00084 __normalize(_InputIterator __first, _InputIterator __last, 00085 _OutputIterator __result, const _Tp& __factor) 00086 { 00087 for (; __first != __last; ++__first, ++__result) 00088 *__result = *__first / __factor; 00089 return __result; 00090 } 00091 00092 _GLIBCXX_END_NAMESPACE_VERSION 00093 } // namespace __detail 00094 00095 _GLIBCXX_BEGIN_NAMESPACE_VERSION 00096 00097 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00098 constexpr _UIntType 00099 linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier; 00100 00101 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00102 constexpr _UIntType 00103 linear_congruential_engine<_UIntType, __a, __c, __m>::increment; 00104 00105 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00106 constexpr _UIntType 00107 linear_congruential_engine<_UIntType, __a, __c, __m>::modulus; 00108 00109 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00110 constexpr _UIntType 00111 linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed; 00112 00113 /** 00114 * Seeds the LCR with integral value @p __s, adjusted so that the 00115 * ring identity is never a member of the convergence set. 00116 */ 00117 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00118 void 00119 linear_congruential_engine<_UIntType, __a, __c, __m>:: 00120 seed(result_type __s) 00121 { 00122 if ((__detail::__mod<_UIntType, __m>(__c) == 0) 00123 && (__detail::__mod<_UIntType, __m>(__s) == 0)) 00124 _M_x = 1; 00125 else 00126 _M_x = __detail::__mod<_UIntType, __m>(__s); 00127 } 00128 00129 /** 00130 * Seeds the LCR engine with a value generated by @p __q. 00131 */ 00132 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00133 template<typename _Sseq> 00134 typename std::enable_if<std::is_class<_Sseq>::value>::type 00135 linear_congruential_engine<_UIntType, __a, __c, __m>:: 00136 seed(_Sseq& __q) 00137 { 00138 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits 00139 : std::__lg(__m); 00140 const _UIntType __k = (__k0 + 31) / 32; 00141 uint_least32_t __arr[__k + 3]; 00142 __q.generate(__arr + 0, __arr + __k + 3); 00143 _UIntType __factor = 1u; 00144 _UIntType __sum = 0u; 00145 for (size_t __j = 0; __j < __k; ++__j) 00146 { 00147 __sum += __arr[__j + 3] * __factor; 00148 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00149 } 00150 seed(__sum); 00151 } 00152 00153 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 00154 typename _CharT, typename _Traits> 00155 std::basic_ostream<_CharT, _Traits>& 00156 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00157 const linear_congruential_engine<_UIntType, 00158 __a, __c, __m>& __lcr) 00159 { 00160 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00161 typedef typename __ostream_type::ios_base __ios_base; 00162 00163 const typename __ios_base::fmtflags __flags = __os.flags(); 00164 const _CharT __fill = __os.fill(); 00165 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00166 __os.fill(__os.widen(' ')); 00167 00168 __os << __lcr._M_x; 00169 00170 __os.flags(__flags); 00171 __os.fill(__fill); 00172 return __os; 00173 } 00174 00175 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 00176 typename _CharT, typename _Traits> 00177 std::basic_istream<_CharT, _Traits>& 00178 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00179 linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr) 00180 { 00181 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00182 typedef typename __istream_type::ios_base __ios_base; 00183 00184 const typename __ios_base::fmtflags __flags = __is.flags(); 00185 __is.flags(__ios_base::dec); 00186 00187 __is >> __lcr._M_x; 00188 00189 __is.flags(__flags); 00190 return __is; 00191 } 00192 00193 00194 template<typename _UIntType, 00195 size_t __w, size_t __n, size_t __m, size_t __r, 00196 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00197 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00198 _UIntType __f> 00199 constexpr size_t 00200 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00201 __s, __b, __t, __c, __l, __f>::word_size; 00202 00203 template<typename _UIntType, 00204 size_t __w, size_t __n, size_t __m, size_t __r, 00205 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00206 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00207 _UIntType __f> 00208 constexpr size_t 00209 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00210 __s, __b, __t, __c, __l, __f>::state_size; 00211 00212 template<typename _UIntType, 00213 size_t __w, size_t __n, size_t __m, size_t __r, 00214 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00215 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00216 _UIntType __f> 00217 constexpr size_t 00218 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00219 __s, __b, __t, __c, __l, __f>::shift_size; 00220 00221 template<typename _UIntType, 00222 size_t __w, size_t __n, size_t __m, size_t __r, 00223 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00224 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00225 _UIntType __f> 00226 constexpr size_t 00227 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00228 __s, __b, __t, __c, __l, __f>::mask_bits; 00229 00230 template<typename _UIntType, 00231 size_t __w, size_t __n, size_t __m, size_t __r, 00232 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00233 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00234 _UIntType __f> 00235 constexpr _UIntType 00236 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00237 __s, __b, __t, __c, __l, __f>::xor_mask; 00238 00239 template<typename _UIntType, 00240 size_t __w, size_t __n, size_t __m, size_t __r, 00241 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00242 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00243 _UIntType __f> 00244 constexpr size_t 00245 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00246 __s, __b, __t, __c, __l, __f>::tempering_u; 00247 00248 template<typename _UIntType, 00249 size_t __w, size_t __n, size_t __m, size_t __r, 00250 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00251 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00252 _UIntType __f> 00253 constexpr _UIntType 00254 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00255 __s, __b, __t, __c, __l, __f>::tempering_d; 00256 00257 template<typename _UIntType, 00258 size_t __w, size_t __n, size_t __m, size_t __r, 00259 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00260 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00261 _UIntType __f> 00262 constexpr size_t 00263 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00264 __s, __b, __t, __c, __l, __f>::tempering_s; 00265 00266 template<typename _UIntType, 00267 size_t __w, size_t __n, size_t __m, size_t __r, 00268 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00269 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00270 _UIntType __f> 00271 constexpr _UIntType 00272 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00273 __s, __b, __t, __c, __l, __f>::tempering_b; 00274 00275 template<typename _UIntType, 00276 size_t __w, size_t __n, size_t __m, size_t __r, 00277 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00278 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00279 _UIntType __f> 00280 constexpr size_t 00281 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00282 __s, __b, __t, __c, __l, __f>::tempering_t; 00283 00284 template<typename _UIntType, 00285 size_t __w, size_t __n, size_t __m, size_t __r, 00286 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00287 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00288 _UIntType __f> 00289 constexpr _UIntType 00290 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00291 __s, __b, __t, __c, __l, __f>::tempering_c; 00292 00293 template<typename _UIntType, 00294 size_t __w, size_t __n, size_t __m, size_t __r, 00295 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00296 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00297 _UIntType __f> 00298 constexpr size_t 00299 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00300 __s, __b, __t, __c, __l, __f>::tempering_l; 00301 00302 template<typename _UIntType, 00303 size_t __w, size_t __n, size_t __m, size_t __r, 00304 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00305 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00306 _UIntType __f> 00307 constexpr _UIntType 00308 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00309 __s, __b, __t, __c, __l, __f>:: 00310 initialization_multiplier; 00311 00312 template<typename _UIntType, 00313 size_t __w, size_t __n, size_t __m, size_t __r, 00314 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00315 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00316 _UIntType __f> 00317 constexpr _UIntType 00318 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00319 __s, __b, __t, __c, __l, __f>::default_seed; 00320 00321 template<typename _UIntType, 00322 size_t __w, size_t __n, size_t __m, size_t __r, 00323 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00324 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00325 _UIntType __f> 00326 void 00327 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00328 __s, __b, __t, __c, __l, __f>:: 00329 seed(result_type __sd) 00330 { 00331 _M_x[0] = __detail::__mod<_UIntType, 00332 __detail::_Shift<_UIntType, __w>::__value>(__sd); 00333 00334 for (size_t __i = 1; __i < state_size; ++__i) 00335 { 00336 _UIntType __x = _M_x[__i - 1]; 00337 __x ^= __x >> (__w - 2); 00338 __x *= __f; 00339 __x += __detail::__mod<_UIntType, __n>(__i); 00340 _M_x[__i] = __detail::__mod<_UIntType, 00341 __detail::_Shift<_UIntType, __w>::__value>(__x); 00342 } 00343 _M_p = state_size; 00344 } 00345 00346 template<typename _UIntType, 00347 size_t __w, size_t __n, size_t __m, size_t __r, 00348 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00349 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00350 _UIntType __f> 00351 template<typename _Sseq> 00352 typename std::enable_if<std::is_class<_Sseq>::value>::type 00353 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00354 __s, __b, __t, __c, __l, __f>:: 00355 seed(_Sseq& __q) 00356 { 00357 const _UIntType __upper_mask = (~_UIntType()) << __r; 00358 const size_t __k = (__w + 31) / 32; 00359 uint_least32_t __arr[__n * __k]; 00360 __q.generate(__arr + 0, __arr + __n * __k); 00361 00362 bool __zero = true; 00363 for (size_t __i = 0; __i < state_size; ++__i) 00364 { 00365 _UIntType __factor = 1u; 00366 _UIntType __sum = 0u; 00367 for (size_t __j = 0; __j < __k; ++__j) 00368 { 00369 __sum += __arr[__k * __i + __j] * __factor; 00370 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00371 } 00372 _M_x[__i] = __detail::__mod<_UIntType, 00373 __detail::_Shift<_UIntType, __w>::__value>(__sum); 00374 00375 if (__zero) 00376 { 00377 if (__i == 0) 00378 { 00379 if ((_M_x[0] & __upper_mask) != 0u) 00380 __zero = false; 00381 } 00382 else if (_M_x[__i] != 0u) 00383 __zero = false; 00384 } 00385 } 00386 if (__zero) 00387 _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value; 00388 _M_p = state_size; 00389 } 00390 00391 template<typename _UIntType, size_t __w, 00392 size_t __n, size_t __m, size_t __r, 00393 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00394 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00395 _UIntType __f> 00396 void 00397 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00398 __s, __b, __t, __c, __l, __f>:: 00399 _M_gen_rand(void) 00400 { 00401 const _UIntType __upper_mask = (~_UIntType()) << __r; 00402 const _UIntType __lower_mask = ~__upper_mask; 00403 00404 for (size_t __k = 0; __k < (__n - __m); ++__k) 00405 { 00406 _UIntType __y = ((_M_x[__k] & __upper_mask) 00407 | (_M_x[__k + 1] & __lower_mask)); 00408 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) 00409 ^ ((__y & 0x01) ? __a : 0)); 00410 } 00411 00412 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k) 00413 { 00414 _UIntType __y = ((_M_x[__k] & __upper_mask) 00415 | (_M_x[__k + 1] & __lower_mask)); 00416 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) 00417 ^ ((__y & 0x01) ? __a : 0)); 00418 } 00419 00420 _UIntType __y = ((_M_x[__n - 1] & __upper_mask) 00421 | (_M_x[0] & __lower_mask)); 00422 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) 00423 ^ ((__y & 0x01) ? __a : 0)); 00424 _M_p = 0; 00425 } 00426 00427 template<typename _UIntType, size_t __w, 00428 size_t __n, size_t __m, size_t __r, 00429 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00430 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00431 _UIntType __f> 00432 void 00433 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00434 __s, __b, __t, __c, __l, __f>:: 00435 discard(unsigned long long __z) 00436 { 00437 while (__z > state_size - _M_p) 00438 { 00439 __z -= state_size - _M_p; 00440 _M_gen_rand(); 00441 } 00442 _M_p += __z; 00443 } 00444 00445 template<typename _UIntType, size_t __w, 00446 size_t __n, size_t __m, size_t __r, 00447 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00448 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00449 _UIntType __f> 00450 typename 00451 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00452 __s, __b, __t, __c, __l, __f>::result_type 00453 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00454 __s, __b, __t, __c, __l, __f>:: 00455 operator()() 00456 { 00457 // Reload the vector - cost is O(n) amortized over n calls. 00458 if (_M_p >= state_size) 00459 _M_gen_rand(); 00460 00461 // Calculate o(x(i)). 00462 result_type __z = _M_x[_M_p++]; 00463 __z ^= (__z >> __u) & __d; 00464 __z ^= (__z << __s) & __b; 00465 __z ^= (__z << __t) & __c; 00466 __z ^= (__z >> __l); 00467 00468 return __z; 00469 } 00470 00471 template<typename _UIntType, size_t __w, 00472 size_t __n, size_t __m, size_t __r, 00473 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00474 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00475 _UIntType __f, typename _CharT, typename _Traits> 00476 std::basic_ostream<_CharT, _Traits>& 00477 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00478 const mersenne_twister_engine<_UIntType, __w, __n, __m, 00479 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 00480 { 00481 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00482 typedef typename __ostream_type::ios_base __ios_base; 00483 00484 const typename __ios_base::fmtflags __flags = __os.flags(); 00485 const _CharT __fill = __os.fill(); 00486 const _CharT __space = __os.widen(' '); 00487 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00488 __os.fill(__space); 00489 00490 for (size_t __i = 0; __i < __n; ++__i) 00491 __os << __x._M_x[__i] << __space; 00492 __os << __x._M_p; 00493 00494 __os.flags(__flags); 00495 __os.fill(__fill); 00496 return __os; 00497 } 00498 00499 template<typename _UIntType, size_t __w, 00500 size_t __n, size_t __m, size_t __r, 00501 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00502 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00503 _UIntType __f, typename _CharT, typename _Traits> 00504 std::basic_istream<_CharT, _Traits>& 00505 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00506 mersenne_twister_engine<_UIntType, __w, __n, __m, 00507 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 00508 { 00509 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00510 typedef typename __istream_type::ios_base __ios_base; 00511 00512 const typename __ios_base::fmtflags __flags = __is.flags(); 00513 __is.flags(__ios_base::dec | __ios_base::skipws); 00514 00515 for (size_t __i = 0; __i < __n; ++__i) 00516 __is >> __x._M_x[__i]; 00517 __is >> __x._M_p; 00518 00519 __is.flags(__flags); 00520 return __is; 00521 } 00522 00523 00524 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00525 constexpr size_t 00526 subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size; 00527 00528 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00529 constexpr size_t 00530 subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag; 00531 00532 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00533 constexpr size_t 00534 subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag; 00535 00536 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00537 constexpr _UIntType 00538 subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed; 00539 00540 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00541 void 00542 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 00543 seed(result_type __value) 00544 { 00545 std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u> 00546 __lcg(__value == 0u ? default_seed : __value); 00547 00548 const size_t __n = (__w + 31) / 32; 00549 00550 for (size_t __i = 0; __i < long_lag; ++__i) 00551 { 00552 _UIntType __sum = 0u; 00553 _UIntType __factor = 1u; 00554 for (size_t __j = 0; __j < __n; ++__j) 00555 { 00556 __sum += __detail::__mod<uint_least32_t, 00557 __detail::_Shift<uint_least32_t, 32>::__value> 00558 (__lcg()) * __factor; 00559 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00560 } 00561 _M_x[__i] = __detail::__mod<_UIntType, 00562 __detail::_Shift<_UIntType, __w>::__value>(__sum); 00563 } 00564 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 00565 _M_p = 0; 00566 } 00567 00568 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00569 template<typename _Sseq> 00570 typename std::enable_if<std::is_class<_Sseq>::value>::type 00571 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 00572 seed(_Sseq& __q) 00573 { 00574 const size_t __k = (__w + 31) / 32; 00575 uint_least32_t __arr[__r * __k]; 00576 __q.generate(__arr + 0, __arr + __r * __k); 00577 00578 for (size_t __i = 0; __i < long_lag; ++__i) 00579 { 00580 _UIntType __sum = 0u; 00581 _UIntType __factor = 1u; 00582 for (size_t __j = 0; __j < __k; ++__j) 00583 { 00584 __sum += __arr[__k * __i + __j] * __factor; 00585 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00586 } 00587 _M_x[__i] = __detail::__mod<_UIntType, 00588 __detail::_Shift<_UIntType, __w>::__value>(__sum); 00589 } 00590 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 00591 _M_p = 0; 00592 } 00593 00594 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00595 typename subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 00596 result_type 00597 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 00598 operator()() 00599 { 00600 // Derive short lag index from current index. 00601 long __ps = _M_p - short_lag; 00602 if (__ps < 0) 00603 __ps += long_lag; 00604 00605 // Calculate new x(i) without overflow or division. 00606 // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry 00607 // cannot overflow. 00608 _UIntType __xi; 00609 if (_M_x[__ps] >= _M_x[_M_p] + _M_carry) 00610 { 00611 __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry; 00612 _M_carry = 0; 00613 } 00614 else 00615 { 00616 __xi = (__detail::_Shift<_UIntType, __w>::__value 00617 - _M_x[_M_p] - _M_carry + _M_x[__ps]); 00618 _M_carry = 1; 00619 } 00620 _M_x[_M_p] = __xi; 00621 00622 // Adjust current index to loop around in ring buffer. 00623 if (++_M_p >= long_lag) 00624 _M_p = 0; 00625 00626 return __xi; 00627 } 00628 00629 template<typename _UIntType, size_t __w, size_t __s, size_t __r, 00630 typename _CharT, typename _Traits> 00631 std::basic_ostream<_CharT, _Traits>& 00632 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00633 const subtract_with_carry_engine<_UIntType, 00634 __w, __s, __r>& __x) 00635 { 00636 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00637 typedef typename __ostream_type::ios_base __ios_base; 00638 00639 const typename __ios_base::fmtflags __flags = __os.flags(); 00640 const _CharT __fill = __os.fill(); 00641 const _CharT __space = __os.widen(' '); 00642 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00643 __os.fill(__space); 00644 00645 for (size_t __i = 0; __i < __r; ++__i) 00646 __os << __x._M_x[__i] << __space; 00647 __os << __x._M_carry << __space << __x._M_p; 00648 00649 __os.flags(__flags); 00650 __os.fill(__fill); 00651 return __os; 00652 } 00653 00654 template<typename _UIntType, size_t __w, size_t __s, size_t __r, 00655 typename _CharT, typename _Traits> 00656 std::basic_istream<_CharT, _Traits>& 00657 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00658 subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x) 00659 { 00660 typedef std::basic_ostream<_CharT, _Traits> __istream_type; 00661 typedef typename __istream_type::ios_base __ios_base; 00662 00663 const typename __ios_base::fmtflags __flags = __is.flags(); 00664 __is.flags(__ios_base::dec | __ios_base::skipws); 00665 00666 for (size_t __i = 0; __i < __r; ++__i) 00667 __is >> __x._M_x[__i]; 00668 __is >> __x._M_carry; 00669 __is >> __x._M_p; 00670 00671 __is.flags(__flags); 00672 return __is; 00673 } 00674 00675 00676 template<typename _RandomNumberEngine, size_t __p, size_t __r> 00677 constexpr size_t 00678 discard_block_engine<_RandomNumberEngine, __p, __r>::block_size; 00679 00680 template<typename _RandomNumberEngine, size_t __p, size_t __r> 00681 constexpr size_t 00682 discard_block_engine<_RandomNumberEngine, __p, __r>::used_block; 00683 00684 template<typename _RandomNumberEngine, size_t __p, size_t __r> 00685 typename discard_block_engine<_RandomNumberEngine, 00686 __p, __r>::result_type 00687 discard_block_engine<_RandomNumberEngine, __p, __r>:: 00688 operator()() 00689 { 00690 if (_M_n >= used_block) 00691 { 00692 _M_b.discard(block_size - _M_n); 00693 _M_n = 0; 00694 } 00695 ++_M_n; 00696 return _M_b(); 00697 } 00698 00699 template<typename _RandomNumberEngine, size_t __p, size_t __r, 00700 typename _CharT, typename _Traits> 00701 std::basic_ostream<_CharT, _Traits>& 00702 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00703 const discard_block_engine<_RandomNumberEngine, 00704 __p, __r>& __x) 00705 { 00706 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00707 typedef typename __ostream_type::ios_base __ios_base; 00708 00709 const typename __ios_base::fmtflags __flags = __os.flags(); 00710 const _CharT __fill = __os.fill(); 00711 const _CharT __space = __os.widen(' '); 00712 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00713 __os.fill(__space); 00714 00715 __os << __x.base() << __space << __x._M_n; 00716 00717 __os.flags(__flags); 00718 __os.fill(__fill); 00719 return __os; 00720 } 00721 00722 template<typename _RandomNumberEngine, size_t __p, size_t __r, 00723 typename _CharT, typename _Traits> 00724 std::basic_istream<_CharT, _Traits>& 00725 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00726 discard_block_engine<_RandomNumberEngine, __p, __r>& __x) 00727 { 00728 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00729 typedef typename __istream_type::ios_base __ios_base; 00730 00731 const typename __ios_base::fmtflags __flags = __is.flags(); 00732 __is.flags(__ios_base::dec | __ios_base::skipws); 00733 00734 __is >> __x._M_b >> __x._M_n; 00735 00736 __is.flags(__flags); 00737 return __is; 00738 } 00739 00740 00741 template<typename _RandomNumberEngine, size_t __w, typename _UIntType> 00742 typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>:: 00743 result_type 00744 independent_bits_engine<_RandomNumberEngine, __w, _UIntType>:: 00745 operator()() 00746 { 00747 typedef typename _RandomNumberEngine::result_type _Eresult_type; 00748 const _Eresult_type __r 00749 = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max() 00750 ? _M_b.max() - _M_b.min() + 1 : 0); 00751 const unsigned __edig = std::numeric_limits<_Eresult_type>::digits; 00752 const unsigned __m = __r ? std::__lg(__r) : __edig; 00753 00754 typedef typename std::common_type<_Eresult_type, result_type>::type 00755 __ctype; 00756 const unsigned __cdig = std::numeric_limits<__ctype>::digits; 00757 00758 unsigned __n, __n0; 00759 __ctype __s0, __s1, __y0, __y1; 00760 00761 for (size_t __i = 0; __i < 2; ++__i) 00762 { 00763 __n = (__w + __m - 1) / __m + __i; 00764 __n0 = __n - __w % __n; 00765 const unsigned __w0 = __w / __n; // __w0 <= __m 00766 00767 __s0 = 0; 00768 __s1 = 0; 00769 if (__w0 < __cdig) 00770 { 00771 __s0 = __ctype(1) << __w0; 00772 __s1 = __s0 << 1; 00773 } 00774 00775 __y0 = 0; 00776 __y1 = 0; 00777 if (__r) 00778 { 00779 __y0 = __s0 * (__r / __s0); 00780 if (__s1) 00781 __y1 = __s1 * (__r / __s1); 00782 00783 if (__r - __y0 <= __y0 / __n) 00784 break; 00785 } 00786 else 00787 break; 00788 } 00789 00790 result_type __sum = 0; 00791 for (size_t __k = 0; __k < __n0; ++__k) 00792 { 00793 __ctype __u; 00794 do 00795 __u = _M_b() - _M_b.min(); 00796 while (__y0 && __u >= __y0); 00797 __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u); 00798 } 00799 for (size_t __k = __n0; __k < __n; ++__k) 00800 { 00801 __ctype __u; 00802 do 00803 __u = _M_b() - _M_b.min(); 00804 while (__y1 && __u >= __y1); 00805 __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u); 00806 } 00807 return __sum; 00808 } 00809 00810 00811 template<typename _RandomNumberEngine, size_t __k> 00812 constexpr size_t 00813 shuffle_order_engine<_RandomNumberEngine, __k>::table_size; 00814 00815 template<typename _RandomNumberEngine, size_t __k> 00816 typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type 00817 shuffle_order_engine<_RandomNumberEngine, __k>:: 00818 operator()() 00819 { 00820 size_t __j = __k * ((_M_y - _M_b.min()) 00821 / (_M_b.max() - _M_b.min() + 1.0L)); 00822 _M_y = _M_v[__j]; 00823 _M_v[__j] = _M_b(); 00824 00825 return _M_y; 00826 } 00827 00828 template<typename _RandomNumberEngine, size_t __k, 00829 typename _CharT, typename _Traits> 00830 std::basic_ostream<_CharT, _Traits>& 00831 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00832 const shuffle_order_engine<_RandomNumberEngine, __k>& __x) 00833 { 00834 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00835 typedef typename __ostream_type::ios_base __ios_base; 00836 00837 const typename __ios_base::fmtflags __flags = __os.flags(); 00838 const _CharT __fill = __os.fill(); 00839 const _CharT __space = __os.widen(' '); 00840 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00841 __os.fill(__space); 00842 00843 __os << __x.base(); 00844 for (size_t __i = 0; __i < __k; ++__i) 00845 __os << __space << __x._M_v[__i]; 00846 __os << __space << __x._M_y; 00847 00848 __os.flags(__flags); 00849 __os.fill(__fill); 00850 return __os; 00851 } 00852 00853 template<typename _RandomNumberEngine, size_t __k, 00854 typename _CharT, typename _Traits> 00855 std::basic_istream<_CharT, _Traits>& 00856 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00857 shuffle_order_engine<_RandomNumberEngine, __k>& __x) 00858 { 00859 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00860 typedef typename __istream_type::ios_base __ios_base; 00861 00862 const typename __ios_base::fmtflags __flags = __is.flags(); 00863 __is.flags(__ios_base::dec | __ios_base::skipws); 00864 00865 __is >> __x._M_b; 00866 for (size_t __i = 0; __i < __k; ++__i) 00867 __is >> __x._M_v[__i]; 00868 __is >> __x._M_y; 00869 00870 __is.flags(__flags); 00871 return __is; 00872 } 00873 00874 00875 template<typename _IntType, typename _CharT, typename _Traits> 00876 std::basic_ostream<_CharT, _Traits>& 00877 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00878 const uniform_int_distribution<_IntType>& __x) 00879 { 00880 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00881 typedef typename __ostream_type::ios_base __ios_base; 00882 00883 const typename __ios_base::fmtflags __flags = __os.flags(); 00884 const _CharT __fill = __os.fill(); 00885 const _CharT __space = __os.widen(' '); 00886 __os.flags(__ios_base::scientific | __ios_base::left); 00887 __os.fill(__space); 00888 00889 __os << __x.a() << __space << __x.b(); 00890 00891 __os.flags(__flags); 00892 __os.fill(__fill); 00893 return __os; 00894 } 00895 00896 template<typename _IntType, typename _CharT, typename _Traits> 00897 std::basic_istream<_CharT, _Traits>& 00898 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00899 uniform_int_distribution<_IntType>& __x) 00900 { 00901 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00902 typedef typename __istream_type::ios_base __ios_base; 00903 00904 const typename __ios_base::fmtflags __flags = __is.flags(); 00905 __is.flags(__ios_base::dec | __ios_base::skipws); 00906 00907 _IntType __a, __b; 00908 __is >> __a >> __b; 00909 __x.param(typename uniform_int_distribution<_IntType>:: 00910 param_type(__a, __b)); 00911 00912 __is.flags(__flags); 00913 return __is; 00914 } 00915 00916 00917 template<typename _RealType> 00918 template<typename _ForwardIterator, 00919 typename _UniformRandomNumberGenerator> 00920 void 00921 uniform_real_distribution<_RealType>:: 00922 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 00923 _UniformRandomNumberGenerator& __urng, 00924 const param_type& __p) 00925 { 00926 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 00927 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 00928 __aurng(__urng); 00929 auto __range = __p.b() - __p.a(); 00930 while (__f != __t) 00931 *__f++ = __aurng() * __range + __p.a(); 00932 } 00933 00934 template<typename _RealType, typename _CharT, typename _Traits> 00935 std::basic_ostream<_CharT, _Traits>& 00936 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00937 const uniform_real_distribution<_RealType>& __x) 00938 { 00939 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00940 typedef typename __ostream_type::ios_base __ios_base; 00941 00942 const typename __ios_base::fmtflags __flags = __os.flags(); 00943 const _CharT __fill = __os.fill(); 00944 const std::streamsize __precision = __os.precision(); 00945 const _CharT __space = __os.widen(' '); 00946 __os.flags(__ios_base::scientific | __ios_base::left); 00947 __os.fill(__space); 00948 __os.precision(std::numeric_limits<_RealType>::max_digits10); 00949 00950 __os << __x.a() << __space << __x.b(); 00951 00952 __os.flags(__flags); 00953 __os.fill(__fill); 00954 __os.precision(__precision); 00955 return __os; 00956 } 00957 00958 template<typename _RealType, typename _CharT, typename _Traits> 00959 std::basic_istream<_CharT, _Traits>& 00960 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00961 uniform_real_distribution<_RealType>& __x) 00962 { 00963 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00964 typedef typename __istream_type::ios_base __ios_base; 00965 00966 const typename __ios_base::fmtflags __flags = __is.flags(); 00967 __is.flags(__ios_base::skipws); 00968 00969 _RealType __a, __b; 00970 __is >> __a >> __b; 00971 __x.param(typename uniform_real_distribution<_RealType>:: 00972 param_type(__a, __b)); 00973 00974 __is.flags(__flags); 00975 return __is; 00976 } 00977 00978 00979 template<typename _ForwardIterator, 00980 typename _UniformRandomNumberGenerator> 00981 void 00982 std::bernoulli_distribution:: 00983 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 00984 _UniformRandomNumberGenerator& __urng, 00985 const param_type& __p) 00986 { 00987 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 00988 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 00989 __aurng(__urng); 00990 auto __limit = __p.p() * (__aurng.max() - __aurng.min()); 00991 00992 while (__f != __t) 00993 *__f++ = (__aurng() - __aurng.min()) < __limit; 00994 } 00995 00996 template<typename _CharT, typename _Traits> 00997 std::basic_ostream<_CharT, _Traits>& 00998 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00999 const bernoulli_distribution& __x) 01000 { 01001 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01002 typedef typename __ostream_type::ios_base __ios_base; 01003 01004 const typename __ios_base::fmtflags __flags = __os.flags(); 01005 const _CharT __fill = __os.fill(); 01006 const std::streamsize __precision = __os.precision(); 01007 __os.flags(__ios_base::scientific | __ios_base::left); 01008 __os.fill(__os.widen(' ')); 01009 __os.precision(std::numeric_limits<double>::max_digits10); 01010 01011 __os << __x.p(); 01012 01013 __os.flags(__flags); 01014 __os.fill(__fill); 01015 __os.precision(__precision); 01016 return __os; 01017 } 01018 01019 01020 template<typename _IntType> 01021 template<typename _UniformRandomNumberGenerator> 01022 typename geometric_distribution<_IntType>::result_type 01023 geometric_distribution<_IntType>:: 01024 operator()(_UniformRandomNumberGenerator& __urng, 01025 const param_type& __param) 01026 { 01027 // About the epsilon thing see this thread: 01028 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html 01029 const double __naf = 01030 (1 - std::numeric_limits<double>::epsilon()) / 2; 01031 // The largest _RealType convertible to _IntType. 01032 const double __thr = 01033 std::numeric_limits<_IntType>::max() + __naf; 01034 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01035 __aurng(__urng); 01036 01037 double __cand; 01038 do 01039 __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p); 01040 while (__cand >= __thr); 01041 01042 return result_type(__cand + __naf); 01043 } 01044 01045 template<typename _IntType> 01046 template<typename _ForwardIterator, 01047 typename _UniformRandomNumberGenerator> 01048 void 01049 geometric_distribution<_IntType>:: 01050 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01051 _UniformRandomNumberGenerator& __urng, 01052 const param_type& __param) 01053 { 01054 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01055 // About the epsilon thing see this thread: 01056 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html 01057 const double __naf = 01058 (1 - std::numeric_limits<double>::epsilon()) / 2; 01059 // The largest _RealType convertible to _IntType. 01060 const double __thr = 01061 std::numeric_limits<_IntType>::max() + __naf; 01062 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01063 __aurng(__urng); 01064 01065 while (__f != __t) 01066 { 01067 double __cand; 01068 do 01069 __cand = std::floor(std::log(1.0 - __aurng()) 01070 / __param._M_log_1_p); 01071 while (__cand >= __thr); 01072 01073 *__f++ = __cand + __naf; 01074 } 01075 } 01076 01077 template<typename _IntType, 01078 typename _CharT, typename _Traits> 01079 std::basic_ostream<_CharT, _Traits>& 01080 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01081 const geometric_distribution<_IntType>& __x) 01082 { 01083 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01084 typedef typename __ostream_type::ios_base __ios_base; 01085 01086 const typename __ios_base::fmtflags __flags = __os.flags(); 01087 const _CharT __fill = __os.fill(); 01088 const std::streamsize __precision = __os.precision(); 01089 __os.flags(__ios_base::scientific | __ios_base::left); 01090 __os.fill(__os.widen(' ')); 01091 __os.precision(std::numeric_limits<double>::max_digits10); 01092 01093 __os << __x.p(); 01094 01095 __os.flags(__flags); 01096 __os.fill(__fill); 01097 __os.precision(__precision); 01098 return __os; 01099 } 01100 01101 template<typename _IntType, 01102 typename _CharT, typename _Traits> 01103 std::basic_istream<_CharT, _Traits>& 01104 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01105 geometric_distribution<_IntType>& __x) 01106 { 01107 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01108 typedef typename __istream_type::ios_base __ios_base; 01109 01110 const typename __ios_base::fmtflags __flags = __is.flags(); 01111 __is.flags(__ios_base::skipws); 01112 01113 double __p; 01114 __is >> __p; 01115 __x.param(typename geometric_distribution<_IntType>::param_type(__p)); 01116 01117 __is.flags(__flags); 01118 return __is; 01119 } 01120 01121 // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5. 01122 template<typename _IntType> 01123 template<typename _UniformRandomNumberGenerator> 01124 typename negative_binomial_distribution<_IntType>::result_type 01125 negative_binomial_distribution<_IntType>:: 01126 operator()(_UniformRandomNumberGenerator& __urng) 01127 { 01128 const double __y = _M_gd(__urng); 01129 01130 // XXX Is the constructor too slow? 01131 std::poisson_distribution<result_type> __poisson(__y); 01132 return __poisson(__urng); 01133 } 01134 01135 template<typename _IntType> 01136 template<typename _UniformRandomNumberGenerator> 01137 typename negative_binomial_distribution<_IntType>::result_type 01138 negative_binomial_distribution<_IntType>:: 01139 operator()(_UniformRandomNumberGenerator& __urng, 01140 const param_type& __p) 01141 { 01142 typedef typename std::gamma_distribution<double>::param_type 01143 param_type; 01144 01145 const double __y = 01146 _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p())); 01147 01148 std::poisson_distribution<result_type> __poisson(__y); 01149 return __poisson(__urng); 01150 } 01151 01152 template<typename _IntType> 01153 template<typename _ForwardIterator, 01154 typename _UniformRandomNumberGenerator> 01155 void 01156 negative_binomial_distribution<_IntType>:: 01157 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01158 _UniformRandomNumberGenerator& __urng) 01159 { 01160 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01161 while (__f != __t) 01162 { 01163 const double __y = _M_gd(__urng); 01164 01165 // XXX Is the constructor too slow? 01166 std::poisson_distribution<result_type> __poisson(__y); 01167 *__f++ = __poisson(__urng); 01168 } 01169 } 01170 01171 template<typename _IntType> 01172 template<typename _ForwardIterator, 01173 typename _UniformRandomNumberGenerator> 01174 void 01175 negative_binomial_distribution<_IntType>:: 01176 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01177 _UniformRandomNumberGenerator& __urng, 01178 const param_type& __p) 01179 { 01180 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01181 typename std::gamma_distribution<result_type>::param_type 01182 __p2(__p.k(), (1.0 - __p.p()) / __p.p()); 01183 01184 while (__f != __t) 01185 { 01186 const double __y = _M_gd(__urng, __p2); 01187 01188 std::poisson_distribution<result_type> __poisson(__y); 01189 *__f++ = __poisson(__urng); 01190 } 01191 } 01192 01193 template<typename _IntType, typename _CharT, typename _Traits> 01194 std::basic_ostream<_CharT, _Traits>& 01195 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01196 const negative_binomial_distribution<_IntType>& __x) 01197 { 01198 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01199 typedef typename __ostream_type::ios_base __ios_base; 01200 01201 const typename __ios_base::fmtflags __flags = __os.flags(); 01202 const _CharT __fill = __os.fill(); 01203 const std::streamsize __precision = __os.precision(); 01204 const _CharT __space = __os.widen(' '); 01205 __os.flags(__ios_base::scientific | __ios_base::left); 01206 __os.fill(__os.widen(' ')); 01207 __os.precision(std::numeric_limits<double>::max_digits10); 01208 01209 __os << __x.k() << __space << __x.p() 01210 << __space << __x._M_gd; 01211 01212 __os.flags(__flags); 01213 __os.fill(__fill); 01214 __os.precision(__precision); 01215 return __os; 01216 } 01217 01218 template<typename _IntType, typename _CharT, typename _Traits> 01219 std::basic_istream<_CharT, _Traits>& 01220 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01221 negative_binomial_distribution<_IntType>& __x) 01222 { 01223 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01224 typedef typename __istream_type::ios_base __ios_base; 01225 01226 const typename __ios_base::fmtflags __flags = __is.flags(); 01227 __is.flags(__ios_base::skipws); 01228 01229 _IntType __k; 01230 double __p; 01231 __is >> __k >> __p >> __x._M_gd; 01232 __x.param(typename negative_binomial_distribution<_IntType>:: 01233 param_type(__k, __p)); 01234 01235 __is.flags(__flags); 01236 return __is; 01237 } 01238 01239 01240 template<typename _IntType> 01241 void 01242 poisson_distribution<_IntType>::param_type:: 01243 _M_initialize() 01244 { 01245 #if _GLIBCXX_USE_C99_MATH_TR1 01246 if (_M_mean >= 12) 01247 { 01248 const double __m = std::floor(_M_mean); 01249 _M_lm_thr = std::log(_M_mean); 01250 _M_lfm = std::lgamma(__m + 1); 01251 _M_sm = std::sqrt(__m); 01252 01253 const double __pi_4 = 0.7853981633974483096156608458198757L; 01254 const double __dx = std::sqrt(2 * __m * std::log(32 * __m 01255 / __pi_4)); 01256 _M_d = std::round(std::max(6.0, std::min(__m, __dx))); 01257 const double __cx = 2 * __m + _M_d; 01258 _M_scx = std::sqrt(__cx / 2); 01259 _M_1cx = 1 / __cx; 01260 01261 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); 01262 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) 01263 / _M_d; 01264 } 01265 else 01266 #endif 01267 _M_lm_thr = std::exp(-_M_mean); 01268 } 01269 01270 /** 01271 * A rejection algorithm when mean >= 12 and a simple method based 01272 * upon the multiplication of uniform random variates otherwise. 01273 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 01274 * is defined. 01275 * 01276 * Reference: 01277 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 01278 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). 01279 */ 01280 template<typename _IntType> 01281 template<typename _UniformRandomNumberGenerator> 01282 typename poisson_distribution<_IntType>::result_type 01283 poisson_distribution<_IntType>:: 01284 operator()(_UniformRandomNumberGenerator& __urng, 01285 const param_type& __param) 01286 { 01287 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01288 __aurng(__urng); 01289 #if _GLIBCXX_USE_C99_MATH_TR1 01290 if (__param.mean() >= 12) 01291 { 01292 double __x; 01293 01294 // See comments above... 01295 const double __naf = 01296 (1 - std::numeric_limits<double>::epsilon()) / 2; 01297 const double __thr = 01298 std::numeric_limits<_IntType>::max() + __naf; 01299 01300 const double __m = std::floor(__param.mean()); 01301 // sqrt(pi / 2) 01302 const double __spi_2 = 1.2533141373155002512078826424055226L; 01303 const double __c1 = __param._M_sm * __spi_2; 01304 const double __c2 = __param._M_c2b + __c1; 01305 const double __c3 = __c2 + 1; 01306 const double __c4 = __c3 + 1; 01307 // e^(1 / 78) 01308 const double __e178 = 1.0129030479320018583185514777512983L; 01309 const double __c5 = __c4 + __e178; 01310 const double __c = __param._M_cb + __c5; 01311 const double __2cx = 2 * (2 * __m + __param._M_d); 01312 01313 bool __reject = true; 01314 do 01315 { 01316 const double __u = __c * __aurng(); 01317 const double __e = -std::log(1.0 - __aurng()); 01318 01319 double __w = 0.0; 01320 01321 if (__u <= __c1) 01322 { 01323 const double __n = _M_nd(__urng); 01324 const double __y = -std::abs(__n) * __param._M_sm - 1; 01325 __x = std::floor(__y); 01326 __w = -__n * __n / 2; 01327 if (__x < -__m) 01328 continue; 01329 } 01330 else if (__u <= __c2) 01331 { 01332 const double __n = _M_nd(__urng); 01333 const double __y = 1 + std::abs(__n) * __param._M_scx; 01334 __x = std::ceil(__y); 01335 __w = __y * (2 - __y) * __param._M_1cx; 01336 if (__x > __param._M_d) 01337 continue; 01338 } 01339 else if (__u <= __c3) 01340 // NB: This case not in the book, nor in the Errata, 01341 // but should be ok... 01342 __x = -1; 01343 else if (__u <= __c4) 01344 __x = 0; 01345 else if (__u <= __c5) 01346 __x = 1; 01347 else 01348 { 01349 const double __v = -std::log(1.0 - __aurng()); 01350 const double __y = __param._M_d 01351 + __v * __2cx / __param._M_d; 01352 __x = std::ceil(__y); 01353 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2); 01354 } 01355 01356 __reject = (__w - __e - __x * __param._M_lm_thr 01357 > __param._M_lfm - std::lgamma(__x + __m + 1)); 01358 01359 __reject |= __x + __m >= __thr; 01360 01361 } while (__reject); 01362 01363 return result_type(__x + __m + __naf); 01364 } 01365 else 01366 #endif 01367 { 01368 _IntType __x = 0; 01369 double __prod = 1.0; 01370 01371 do 01372 { 01373 __prod *= __aurng(); 01374 __x += 1; 01375 } 01376 while (__prod > __param._M_lm_thr); 01377 01378 return __x - 1; 01379 } 01380 } 01381 01382 template<typename _IntType> 01383 template<typename _ForwardIterator, 01384 typename _UniformRandomNumberGenerator> 01385 void 01386 poisson_distribution<_IntType>:: 01387 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01388 _UniformRandomNumberGenerator& __urng, 01389 const param_type& __param) 01390 { 01391 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01392 // We could duplicate everything from operator()... 01393 while (__f != __t) 01394 *__f++ = this->operator()(__urng, __param); 01395 } 01396 01397 template<typename _IntType, 01398 typename _CharT, typename _Traits> 01399 std::basic_ostream<_CharT, _Traits>& 01400 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01401 const poisson_distribution<_IntType>& __x) 01402 { 01403 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01404 typedef typename __ostream_type::ios_base __ios_base; 01405 01406 const typename __ios_base::fmtflags __flags = __os.flags(); 01407 const _CharT __fill = __os.fill(); 01408 const std::streamsize __precision = __os.precision(); 01409 const _CharT __space = __os.widen(' '); 01410 __os.flags(__ios_base::scientific | __ios_base::left); 01411 __os.fill(__space); 01412 __os.precision(std::numeric_limits<double>::max_digits10); 01413 01414 __os << __x.mean() << __space << __x._M_nd; 01415 01416 __os.flags(__flags); 01417 __os.fill(__fill); 01418 __os.precision(__precision); 01419 return __os; 01420 } 01421 01422 template<typename _IntType, 01423 typename _CharT, typename _Traits> 01424 std::basic_istream<_CharT, _Traits>& 01425 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01426 poisson_distribution<_IntType>& __x) 01427 { 01428 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01429 typedef typename __istream_type::ios_base __ios_base; 01430 01431 const typename __ios_base::fmtflags __flags = __is.flags(); 01432 __is.flags(__ios_base::skipws); 01433 01434 double __mean; 01435 __is >> __mean >> __x._M_nd; 01436 __x.param(typename poisson_distribution<_IntType>::param_type(__mean)); 01437 01438 __is.flags(__flags); 01439 return __is; 01440 } 01441 01442 01443 template<typename _IntType> 01444 void 01445 binomial_distribution<_IntType>::param_type:: 01446 _M_initialize() 01447 { 01448 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 01449 01450 _M_easy = true; 01451 01452 #if _GLIBCXX_USE_C99_MATH_TR1 01453 if (_M_t * __p12 >= 8) 01454 { 01455 _M_easy = false; 01456 const double __np = std::floor(_M_t * __p12); 01457 const double __pa = __np / _M_t; 01458 const double __1p = 1 - __pa; 01459 01460 const double __pi_4 = 0.7853981633974483096156608458198757L; 01461 const double __d1x = 01462 std::sqrt(__np * __1p * std::log(32 * __np 01463 / (81 * __pi_4 * __1p))); 01464 _M_d1 = std::round(std::max(1.0, __d1x)); 01465 const double __d2x = 01466 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p 01467 / (__pi_4 * __pa))); 01468 _M_d2 = std::round(std::max(1.0, __d2x)); 01469 01470 // sqrt(pi / 2) 01471 const double __spi_2 = 1.2533141373155002512078826424055226L; 01472 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); 01473 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); 01474 _M_c = 2 * _M_d1 / __np; 01475 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; 01476 const double __a12 = _M_a1 + _M_s2 * __spi_2; 01477 const double __s1s = _M_s1 * _M_s1; 01478 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) 01479 * 2 * __s1s / _M_d1 01480 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); 01481 const double __s2s = _M_s2 * _M_s2; 01482 _M_s = (_M_a123 + 2 * __s2s / _M_d2 01483 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); 01484 _M_lf = (std::lgamma(__np + 1) 01485 + std::lgamma(_M_t - __np + 1)); 01486 _M_lp1p = std::log(__pa / __1p); 01487 01488 _M_q = -std::log(1 - (__p12 - __pa) / __1p); 01489 } 01490 else 01491 #endif 01492 _M_q = -std::log(1 - __p12); 01493 } 01494 01495 template<typename _IntType> 01496 template<typename _UniformRandomNumberGenerator> 01497 typename binomial_distribution<_IntType>::result_type 01498 binomial_distribution<_IntType>:: 01499 _M_waiting(_UniformRandomNumberGenerator& __urng, 01500 _IntType __t, double __q) 01501 { 01502 _IntType __x = 0; 01503 double __sum = 0.0; 01504 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01505 __aurng(__urng); 01506 01507 do 01508 { 01509 if (__t == __x) 01510 return __x; 01511 const double __e = -std::log(1.0 - __aurng()); 01512 __sum += __e / (__t - __x); 01513 __x += 1; 01514 } 01515 while (__sum <= __q); 01516 01517 return __x - 1; 01518 } 01519 01520 /** 01521 * A rejection algorithm when t * p >= 8 and a simple waiting time 01522 * method - the second in the referenced book - otherwise. 01523 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 01524 * is defined. 01525 * 01526 * Reference: 01527 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 01528 * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 01529 */ 01530 template<typename _IntType> 01531 template<typename _UniformRandomNumberGenerator> 01532 typename binomial_distribution<_IntType>::result_type 01533 binomial_distribution<_IntType>:: 01534 operator()(_UniformRandomNumberGenerator& __urng, 01535 const param_type& __param) 01536 { 01537 result_type __ret; 01538 const _IntType __t = __param.t(); 01539 const double __p = __param.p(); 01540 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p; 01541 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01542 __aurng(__urng); 01543 01544 #if _GLIBCXX_USE_C99_MATH_TR1 01545 if (!__param._M_easy) 01546 { 01547 double __x; 01548 01549 // See comments above... 01550 const double __naf = 01551 (1 - std::numeric_limits<double>::epsilon()) / 2; 01552 const double __thr = 01553 std::numeric_limits<_IntType>::max() + __naf; 01554 01555 const double __np = std::floor(__t * __p12); 01556 01557 // sqrt(pi / 2) 01558 const double __spi_2 = 1.2533141373155002512078826424055226L; 01559 const double __a1 = __param._M_a1; 01560 const double __a12 = __a1 + __param._M_s2 * __spi_2; 01561 const double __a123 = __param._M_a123; 01562 const double __s1s = __param._M_s1 * __param._M_s1; 01563 const double __s2s = __param._M_s2 * __param._M_s2; 01564 01565 bool __reject; 01566 do 01567 { 01568 const double __u = __param._M_s * __aurng(); 01569 01570 double __v; 01571 01572 if (__u <= __a1) 01573 { 01574 const double __n = _M_nd(__urng); 01575 const double __y = __param._M_s1 * std::abs(__n); 01576 __reject = __y >= __param._M_d1; 01577 if (!__reject) 01578 { 01579 const double __e = -std::log(1.0 - __aurng()); 01580 __x = std::floor(__y); 01581 __v = -__e - __n * __n / 2 + __param._M_c; 01582 } 01583 } 01584 else if (__u <= __a12) 01585 { 01586 const double __n = _M_nd(__urng); 01587 const double __y = __param._M_s2 * std::abs(__n); 01588 __reject = __y >= __param._M_d2; 01589 if (!__reject) 01590 { 01591 const double __e = -std::log(1.0 - __aurng()); 01592 __x = std::floor(-__y); 01593 __v = -__e - __n * __n / 2; 01594 } 01595 } 01596 else if (__u <= __a123) 01597 { 01598 const double __e1 = -std::log(1.0 - __aurng()); 01599 const double __e2 = -std::log(1.0 - __aurng()); 01600 01601 const double __y = __param._M_d1 01602 + 2 * __s1s * __e1 / __param._M_d1; 01603 __x = std::floor(__y); 01604 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np) 01605 -__y / (2 * __s1s))); 01606 __reject = false; 01607 } 01608 else 01609 { 01610 const double __e1 = -std::log(1.0 - __aurng()); 01611 const double __e2 = -std::log(1.0 - __aurng()); 01612 01613 const double __y = __param._M_d2 01614 + 2 * __s2s * __e1 / __param._M_d2; 01615 __x = std::floor(-__y); 01616 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s); 01617 __reject = false; 01618 } 01619 01620 __reject = __reject || __x < -__np || __x > __t - __np; 01621 if (!__reject) 01622 { 01623 const double __lfx = 01624 std::lgamma(__np + __x + 1) 01625 + std::lgamma(__t - (__np + __x) + 1); 01626 __reject = __v > __param._M_lf - __lfx 01627 + __x * __param._M_lp1p; 01628 } 01629 01630 __reject |= __x + __np >= __thr; 01631 } 01632 while (__reject); 01633 01634 __x += __np + __naf; 01635 01636 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x), 01637 __param._M_q); 01638 __ret = _IntType(__x) + __z; 01639 } 01640 else 01641 #endif 01642 __ret = _M_waiting(__urng, __t, __param._M_q); 01643 01644 if (__p12 != __p) 01645 __ret = __t - __ret; 01646 return __ret; 01647 } 01648 01649 template<typename _IntType> 01650 template<typename _ForwardIterator, 01651 typename _UniformRandomNumberGenerator> 01652 void 01653 binomial_distribution<_IntType>:: 01654 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01655 _UniformRandomNumberGenerator& __urng, 01656 const param_type& __param) 01657 { 01658 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01659 // We could duplicate everything from operator()... 01660 while (__f != __t) 01661 *__f++ = this->operator()(__urng, __param); 01662 } 01663 01664 template<typename _IntType, 01665 typename _CharT, typename _Traits> 01666 std::basic_ostream<_CharT, _Traits>& 01667 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01668 const binomial_distribution<_IntType>& __x) 01669 { 01670 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01671 typedef typename __ostream_type::ios_base __ios_base; 01672 01673 const typename __ios_base::fmtflags __flags = __os.flags(); 01674 const _CharT __fill = __os.fill(); 01675 const std::streamsize __precision = __os.precision(); 01676 const _CharT __space = __os.widen(' '); 01677 __os.flags(__ios_base::scientific | __ios_base::left); 01678 __os.fill(__space); 01679 __os.precision(std::numeric_limits<double>::max_digits10); 01680 01681 __os << __x.t() << __space << __x.p() 01682 << __space << __x._M_nd; 01683 01684 __os.flags(__flags); 01685 __os.fill(__fill); 01686 __os.precision(__precision); 01687 return __os; 01688 } 01689 01690 template<typename _IntType, 01691 typename _CharT, typename _Traits> 01692 std::basic_istream<_CharT, _Traits>& 01693 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01694 binomial_distribution<_IntType>& __x) 01695 { 01696 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01697 typedef typename __istream_type::ios_base __ios_base; 01698 01699 const typename __ios_base::fmtflags __flags = __is.flags(); 01700 __is.flags(__ios_base::dec | __ios_base::skipws); 01701 01702 _IntType __t; 01703 double __p; 01704 __is >> __t >> __p >> __x._M_nd; 01705 __x.param(typename binomial_distribution<_IntType>:: 01706 param_type(__t, __p)); 01707 01708 __is.flags(__flags); 01709 return __is; 01710 } 01711 01712 01713 template<typename _RealType> 01714 template<typename _ForwardIterator, 01715 typename _UniformRandomNumberGenerator> 01716 void 01717 std::exponential_distribution<_RealType>:: 01718 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01719 _UniformRandomNumberGenerator& __urng, 01720 const param_type& __p) 01721 { 01722 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01723 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01724 __aurng(__urng); 01725 while (__f != __t) 01726 *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda(); 01727 } 01728 01729 template<typename _RealType, typename _CharT, typename _Traits> 01730 std::basic_ostream<_CharT, _Traits>& 01731 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01732 const exponential_distribution<_RealType>& __x) 01733 { 01734 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01735 typedef typename __ostream_type::ios_base __ios_base; 01736 01737 const typename __ios_base::fmtflags __flags = __os.flags(); 01738 const _CharT __fill = __os.fill(); 01739 const std::streamsize __precision = __os.precision(); 01740 __os.flags(__ios_base::scientific | __ios_base::left); 01741 __os.fill(__os.widen(' ')); 01742 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01743 01744 __os << __x.lambda(); 01745 01746 __os.flags(__flags); 01747 __os.fill(__fill); 01748 __os.precision(__precision); 01749 return __os; 01750 } 01751 01752 template<typename _RealType, typename _CharT, typename _Traits> 01753 std::basic_istream<_CharT, _Traits>& 01754 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01755 exponential_distribution<_RealType>& __x) 01756 { 01757 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01758 typedef typename __istream_type::ios_base __ios_base; 01759 01760 const typename __ios_base::fmtflags __flags = __is.flags(); 01761 __is.flags(__ios_base::dec | __ios_base::skipws); 01762 01763 _RealType __lambda; 01764 __is >> __lambda; 01765 __x.param(typename exponential_distribution<_RealType>:: 01766 param_type(__lambda)); 01767 01768 __is.flags(__flags); 01769 return __is; 01770 } 01771 01772 01773 /** 01774 * Polar method due to Marsaglia. 01775 * 01776 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 01777 * New York, 1986, Ch. V, Sect. 4.4. 01778 */ 01779 template<typename _RealType> 01780 template<typename _UniformRandomNumberGenerator> 01781 typename normal_distribution<_RealType>::result_type 01782 normal_distribution<_RealType>:: 01783 operator()(_UniformRandomNumberGenerator& __urng, 01784 const param_type& __param) 01785 { 01786 result_type __ret; 01787 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01788 __aurng(__urng); 01789 01790 if (_M_saved_available) 01791 { 01792 _M_saved_available = false; 01793 __ret = _M_saved; 01794 } 01795 else 01796 { 01797 result_type __x, __y, __r2; 01798 do 01799 { 01800 __x = result_type(2.0) * __aurng() - 1.0; 01801 __y = result_type(2.0) * __aurng() - 1.0; 01802 __r2 = __x * __x + __y * __y; 01803 } 01804 while (__r2 > 1.0 || __r2 == 0.0); 01805 01806 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01807 _M_saved = __x * __mult; 01808 _M_saved_available = true; 01809 __ret = __y * __mult; 01810 } 01811 01812 __ret = __ret * __param.stddev() + __param.mean(); 01813 return __ret; 01814 } 01815 01816 template<typename _RealType> 01817 template<typename _ForwardIterator, 01818 typename _UniformRandomNumberGenerator> 01819 void 01820 normal_distribution<_RealType>:: 01821 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01822 _UniformRandomNumberGenerator& __urng, 01823 const param_type& __param) 01824 { 01825 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01826 01827 if (__f == __t) 01828 return; 01829 01830 if (_M_saved_available) 01831 { 01832 _M_saved_available = false; 01833 *__f++ = _M_saved * __param.stddev() + __param.mean(); 01834 01835 if (__f == __t) 01836 return; 01837 } 01838 01839 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01840 __aurng(__urng); 01841 01842 while (__f + 1 < __t) 01843 { 01844 result_type __x, __y, __r2; 01845 do 01846 { 01847 __x = result_type(2.0) * __aurng() - 1.0; 01848 __y = result_type(2.0) * __aurng() - 1.0; 01849 __r2 = __x * __x + __y * __y; 01850 } 01851 while (__r2 > 1.0 || __r2 == 0.0); 01852 01853 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01854 *__f++ = __y * __mult * __param.stddev() + __param.mean(); 01855 *__f++ = __x * __mult * __param.stddev() + __param.mean(); 01856 } 01857 01858 if (__f != __t) 01859 { 01860 result_type __x, __y, __r2; 01861 do 01862 { 01863 __x = result_type(2.0) * __aurng() - 1.0; 01864 __y = result_type(2.0) * __aurng() - 1.0; 01865 __r2 = __x * __x + __y * __y; 01866 } 01867 while (__r2 > 1.0 || __r2 == 0.0); 01868 01869 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01870 _M_saved = __x * __mult; 01871 _M_saved_available = true; 01872 *__f = __y * __mult * __param.stddev() + __param.mean(); 01873 } 01874 } 01875 01876 template<typename _RealType> 01877 bool 01878 operator==(const std::normal_distribution<_RealType>& __d1, 01879 const std::normal_distribution<_RealType>& __d2) 01880 { 01881 if (__d1._M_param == __d2._M_param 01882 && __d1._M_saved_available == __d2._M_saved_available) 01883 { 01884 if (__d1._M_saved_available 01885 && __d1._M_saved == __d2._M_saved) 01886 return true; 01887 else if(!__d1._M_saved_available) 01888 return true; 01889 else 01890 return false; 01891 } 01892 else 01893 return false; 01894 } 01895 01896 template<typename _RealType, typename _CharT, typename _Traits> 01897 std::basic_ostream<_CharT, _Traits>& 01898 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01899 const normal_distribution<_RealType>& __x) 01900 { 01901 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01902 typedef typename __ostream_type::ios_base __ios_base; 01903 01904 const typename __ios_base::fmtflags __flags = __os.flags(); 01905 const _CharT __fill = __os.fill(); 01906 const std::streamsize __precision = __os.precision(); 01907 const _CharT __space = __os.widen(' '); 01908 __os.flags(__ios_base::scientific | __ios_base::left); 01909 __os.fill(__space); 01910 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01911 01912 __os << __x.mean() << __space << __x.stddev() 01913 << __space << __x._M_saved_available; 01914 if (__x._M_saved_available) 01915 __os << __space << __x._M_saved; 01916 01917 __os.flags(__flags); 01918 __os.fill(__fill); 01919 __os.precision(__precision); 01920 return __os; 01921 } 01922 01923 template<typename _RealType, typename _CharT, typename _Traits> 01924 std::basic_istream<_CharT, _Traits>& 01925 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01926 normal_distribution<_RealType>& __x) 01927 { 01928 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01929 typedef typename __istream_type::ios_base __ios_base; 01930 01931 const typename __ios_base::fmtflags __flags = __is.flags(); 01932 __is.flags(__ios_base::dec | __ios_base::skipws); 01933 01934 double __mean, __stddev; 01935 __is >> __mean >> __stddev 01936 >> __x._M_saved_available; 01937 if (__x._M_saved_available) 01938 __is >> __x._M_saved; 01939 __x.param(typename normal_distribution<_RealType>:: 01940 param_type(__mean, __stddev)); 01941 01942 __is.flags(__flags); 01943 return __is; 01944 } 01945 01946 01947 template<typename _RealType> 01948 template<typename _ForwardIterator, 01949 typename _UniformRandomNumberGenerator> 01950 void 01951 lognormal_distribution<_RealType>:: 01952 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01953 _UniformRandomNumberGenerator& __urng, 01954 const param_type& __p) 01955 { 01956 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01957 while (__f != __t) 01958 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m()); 01959 } 01960 01961 template<typename _RealType, typename _CharT, typename _Traits> 01962 std::basic_ostream<_CharT, _Traits>& 01963 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01964 const lognormal_distribution<_RealType>& __x) 01965 { 01966 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01967 typedef typename __ostream_type::ios_base __ios_base; 01968 01969 const typename __ios_base::fmtflags __flags = __os.flags(); 01970 const _CharT __fill = __os.fill(); 01971 const std::streamsize __precision = __os.precision(); 01972 const _CharT __space = __os.widen(' '); 01973 __os.flags(__ios_base::scientific | __ios_base::left); 01974 __os.fill(__space); 01975 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01976 01977 __os << __x.m() << __space << __x.s() 01978 << __space << __x._M_nd; 01979 01980 __os.flags(__flags); 01981 __os.fill(__fill); 01982 __os.precision(__precision); 01983 return __os; 01984 } 01985 01986 template<typename _RealType, typename _CharT, typename _Traits> 01987 std::basic_istream<_CharT, _Traits>& 01988 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01989 lognormal_distribution<_RealType>& __x) 01990 { 01991 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01992 typedef typename __istream_type::ios_base __ios_base; 01993 01994 const typename __ios_base::fmtflags __flags = __is.flags(); 01995 __is.flags(__ios_base::dec | __ios_base::skipws); 01996 01997 _RealType __m, __s; 01998 __is >> __m >> __s >> __x._M_nd; 01999 __x.param(typename lognormal_distribution<_RealType>:: 02000 param_type(__m, __s)); 02001 02002 __is.flags(__flags); 02003 return __is; 02004 } 02005 02006 template<typename _RealType> 02007 template<typename _ForwardIterator, 02008 typename _UniformRandomNumberGenerator> 02009 void 02010 std::chi_squared_distribution<_RealType>:: 02011 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02012 _UniformRandomNumberGenerator& __urng) 02013 { 02014 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02015 while (__f != __t) 02016 *__f++ = 2 * _M_gd(__urng); 02017 } 02018 02019 template<typename _RealType> 02020 template<typename _ForwardIterator, 02021 typename _UniformRandomNumberGenerator> 02022 void 02023 std::chi_squared_distribution<_RealType>:: 02024 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02025 _UniformRandomNumberGenerator& __urng, 02026 const typename 02027 std::gamma_distribution<result_type>::param_type& __p) 02028 { 02029 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02030 while (__f != __t) 02031 *__f++ = 2 * _M_gd(__urng, __p); 02032 } 02033 02034 template<typename _RealType, typename _CharT, typename _Traits> 02035 std::basic_ostream<_CharT, _Traits>& 02036 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02037 const chi_squared_distribution<_RealType>& __x) 02038 { 02039 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02040 typedef typename __ostream_type::ios_base __ios_base; 02041 02042 const typename __ios_base::fmtflags __flags = __os.flags(); 02043 const _CharT __fill = __os.fill(); 02044 const std::streamsize __precision = __os.precision(); 02045 const _CharT __space = __os.widen(' '); 02046 __os.flags(__ios_base::scientific | __ios_base::left); 02047 __os.fill(__space); 02048 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02049 02050 __os << __x.n() << __space << __x._M_gd; 02051 02052 __os.flags(__flags); 02053 __os.fill(__fill); 02054 __os.precision(__precision); 02055 return __os; 02056 } 02057 02058 template<typename _RealType, typename _CharT, typename _Traits> 02059 std::basic_istream<_CharT, _Traits>& 02060 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02061 chi_squared_distribution<_RealType>& __x) 02062 { 02063 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02064 typedef typename __istream_type::ios_base __ios_base; 02065 02066 const typename __ios_base::fmtflags __flags = __is.flags(); 02067 __is.flags(__ios_base::dec | __ios_base::skipws); 02068 02069 _RealType __n; 02070 __is >> __n >> __x._M_gd; 02071 __x.param(typename chi_squared_distribution<_RealType>:: 02072 param_type(__n)); 02073 02074 __is.flags(__flags); 02075 return __is; 02076 } 02077 02078 02079 template<typename _RealType> 02080 template<typename _UniformRandomNumberGenerator> 02081 typename cauchy_distribution<_RealType>::result_type 02082 cauchy_distribution<_RealType>:: 02083 operator()(_UniformRandomNumberGenerator& __urng, 02084 const param_type& __p) 02085 { 02086 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02087 __aurng(__urng); 02088 _RealType __u; 02089 do 02090 __u = __aurng(); 02091 while (__u == 0.5); 02092 02093 const _RealType __pi = 3.1415926535897932384626433832795029L; 02094 return __p.a() + __p.b() * std::tan(__pi * __u); 02095 } 02096 02097 template<typename _RealType> 02098 template<typename _ForwardIterator, 02099 typename _UniformRandomNumberGenerator> 02100 void 02101 cauchy_distribution<_RealType>:: 02102 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02103 _UniformRandomNumberGenerator& __urng, 02104 const param_type& __p) 02105 { 02106 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02107 const _RealType __pi = 3.1415926535897932384626433832795029L; 02108 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02109 __aurng(__urng); 02110 while (__f != __t) 02111 { 02112 _RealType __u; 02113 do 02114 __u = __aurng(); 02115 while (__u == 0.5); 02116 02117 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u); 02118 } 02119 } 02120 02121 template<typename _RealType, typename _CharT, typename _Traits> 02122 std::basic_ostream<_CharT, _Traits>& 02123 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02124 const cauchy_distribution<_RealType>& __x) 02125 { 02126 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02127 typedef typename __ostream_type::ios_base __ios_base; 02128 02129 const typename __ios_base::fmtflags __flags = __os.flags(); 02130 const _CharT __fill = __os.fill(); 02131 const std::streamsize __precision = __os.precision(); 02132 const _CharT __space = __os.widen(' '); 02133 __os.flags(__ios_base::scientific | __ios_base::left); 02134 __os.fill(__space); 02135 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02136 02137 __os << __x.a() << __space << __x.b(); 02138 02139 __os.flags(__flags); 02140 __os.fill(__fill); 02141 __os.precision(__precision); 02142 return __os; 02143 } 02144 02145 template<typename _RealType, typename _CharT, typename _Traits> 02146 std::basic_istream<_CharT, _Traits>& 02147 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02148 cauchy_distribution<_RealType>& __x) 02149 { 02150 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02151 typedef typename __istream_type::ios_base __ios_base; 02152 02153 const typename __ios_base::fmtflags __flags = __is.flags(); 02154 __is.flags(__ios_base::dec | __ios_base::skipws); 02155 02156 _RealType __a, __b; 02157 __is >> __a >> __b; 02158 __x.param(typename cauchy_distribution<_RealType>:: 02159 param_type(__a, __b)); 02160 02161 __is.flags(__flags); 02162 return __is; 02163 } 02164 02165 02166 template<typename _RealType> 02167 template<typename _ForwardIterator, 02168 typename _UniformRandomNumberGenerator> 02169 void 02170 std::fisher_f_distribution<_RealType>:: 02171 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02172 _UniformRandomNumberGenerator& __urng) 02173 { 02174 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02175 while (__f != __t) 02176 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m())); 02177 } 02178 02179 template<typename _RealType> 02180 template<typename _ForwardIterator, 02181 typename _UniformRandomNumberGenerator> 02182 void 02183 std::fisher_f_distribution<_RealType>:: 02184 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02185 _UniformRandomNumberGenerator& __urng, 02186 const param_type& __p) 02187 { 02188 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02189 typedef typename std::gamma_distribution<result_type>::param_type 02190 param_type; 02191 param_type __p1(__p.m() / 2); 02192 param_type __p2(__p.n() / 2); 02193 while (__f != __t) 02194 *__f++ = ((_M_gd_x(__urng, __p1) * n()) 02195 / (_M_gd_y(__urng, __p2) * m())); 02196 } 02197 02198 template<typename _RealType, typename _CharT, typename _Traits> 02199 std::basic_ostream<_CharT, _Traits>& 02200 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02201 const fisher_f_distribution<_RealType>& __x) 02202 { 02203 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02204 typedef typename __ostream_type::ios_base __ios_base; 02205 02206 const typename __ios_base::fmtflags __flags = __os.flags(); 02207 const _CharT __fill = __os.fill(); 02208 const std::streamsize __precision = __os.precision(); 02209 const _CharT __space = __os.widen(' '); 02210 __os.flags(__ios_base::scientific | __ios_base::left); 02211 __os.fill(__space); 02212 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02213 02214 __os << __x.m() << __space << __x.n() 02215 << __space << __x._M_gd_x << __space << __x._M_gd_y; 02216 02217 __os.flags(__flags); 02218 __os.fill(__fill); 02219 __os.precision(__precision); 02220 return __os; 02221 } 02222 02223 template<typename _RealType, typename _CharT, typename _Traits> 02224 std::basic_istream<_CharT, _Traits>& 02225 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02226 fisher_f_distribution<_RealType>& __x) 02227 { 02228 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02229 typedef typename __istream_type::ios_base __ios_base; 02230 02231 const typename __ios_base::fmtflags __flags = __is.flags(); 02232 __is.flags(__ios_base::dec | __ios_base::skipws); 02233 02234 _RealType __m, __n; 02235 __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y; 02236 __x.param(typename fisher_f_distribution<_RealType>:: 02237 param_type(__m, __n)); 02238 02239 __is.flags(__flags); 02240 return __is; 02241 } 02242 02243 02244 template<typename _RealType> 02245 template<typename _ForwardIterator, 02246 typename _UniformRandomNumberGenerator> 02247 void 02248 std::student_t_distribution<_RealType>:: 02249 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02250 _UniformRandomNumberGenerator& __urng) 02251 { 02252 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02253 while (__f != __t) 02254 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng)); 02255 } 02256 02257 template<typename _RealType> 02258 template<typename _ForwardIterator, 02259 typename _UniformRandomNumberGenerator> 02260 void 02261 std::student_t_distribution<_RealType>:: 02262 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02263 _UniformRandomNumberGenerator& __urng, 02264 const param_type& __p) 02265 { 02266 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02267 typename std::gamma_distribution<result_type>::param_type 02268 __p2(__p.n() / 2, 2); 02269 while (__f != __t) 02270 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2)); 02271 } 02272 02273 template<typename _RealType, typename _CharT, typename _Traits> 02274 std::basic_ostream<_CharT, _Traits>& 02275 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02276 const student_t_distribution<_RealType>& __x) 02277 { 02278 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02279 typedef typename __ostream_type::ios_base __ios_base; 02280 02281 const typename __ios_base::fmtflags __flags = __os.flags(); 02282 const _CharT __fill = __os.fill(); 02283 const std::streamsize __precision = __os.precision(); 02284 const _CharT __space = __os.widen(' '); 02285 __os.flags(__ios_base::scientific | __ios_base::left); 02286 __os.fill(__space); 02287 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02288 02289 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd; 02290 02291 __os.flags(__flags); 02292 __os.fill(__fill); 02293 __os.precision(__precision); 02294 return __os; 02295 } 02296 02297 template<typename _RealType, typename _CharT, typename _Traits> 02298 std::basic_istream<_CharT, _Traits>& 02299 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02300 student_t_distribution<_RealType>& __x) 02301 { 02302 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02303 typedef typename __istream_type::ios_base __ios_base; 02304 02305 const typename __ios_base::fmtflags __flags = __is.flags(); 02306 __is.flags(__ios_base::dec | __ios_base::skipws); 02307 02308 _RealType __n; 02309 __is >> __n >> __x._M_nd >> __x._M_gd; 02310 __x.param(typename student_t_distribution<_RealType>::param_type(__n)); 02311 02312 __is.flags(__flags); 02313 return __is; 02314 } 02315 02316 02317 template<typename _RealType> 02318 void 02319 gamma_distribution<_RealType>::param_type:: 02320 _M_initialize() 02321 { 02322 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha; 02323 02324 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0); 02325 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1); 02326 } 02327 02328 /** 02329 * Marsaglia, G. and Tsang, W. W. 02330 * "A Simple Method for Generating Gamma Variables" 02331 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000. 02332 */ 02333 template<typename _RealType> 02334 template<typename _UniformRandomNumberGenerator> 02335 typename gamma_distribution<_RealType>::result_type 02336 gamma_distribution<_RealType>:: 02337 operator()(_UniformRandomNumberGenerator& __urng, 02338 const param_type& __param) 02339 { 02340 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02341 __aurng(__urng); 02342 02343 result_type __u, __v, __n; 02344 const result_type __a1 = (__param._M_malpha 02345 - _RealType(1.0) / _RealType(3.0)); 02346 02347 do 02348 { 02349 do 02350 { 02351 __n = _M_nd(__urng); 02352 __v = result_type(1.0) + __param._M_a2 * __n; 02353 } 02354 while (__v <= 0.0); 02355 02356 __v = __v * __v * __v; 02357 __u = __aurng(); 02358 } 02359 while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 02360 && (std::log(__u) > (0.5 * __n * __n + __a1 02361 * (1.0 - __v + std::log(__v))))); 02362 02363 if (__param.alpha() == __param._M_malpha) 02364 return __a1 * __v * __param.beta(); 02365 else 02366 { 02367 do 02368 __u = __aurng(); 02369 while (__u == 0.0); 02370 02371 return (std::pow(__u, result_type(1.0) / __param.alpha()) 02372 * __a1 * __v * __param.beta()); 02373 } 02374 } 02375 02376 template<typename _RealType> 02377 template<typename _ForwardIterator, 02378 typename _UniformRandomNumberGenerator> 02379 void 02380 gamma_distribution<_RealType>:: 02381 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02382 _UniformRandomNumberGenerator& __urng, 02383 const param_type& __param) 02384 { 02385 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02386 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02387 __aurng(__urng); 02388 02389 result_type __u, __v, __n; 02390 const result_type __a1 = (__param._M_malpha 02391 - _RealType(1.0) / _RealType(3.0)); 02392 02393 if (__param.alpha() == __param._M_malpha) 02394 while (__f != __t) 02395 { 02396 do 02397 { 02398 do 02399 { 02400 __n = _M_nd(__urng); 02401 __v = result_type(1.0) + __param._M_a2 * __n; 02402 } 02403 while (__v <= 0.0); 02404 02405 __v = __v * __v * __v; 02406 __u = __aurng(); 02407 } 02408 while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 02409 && (std::log(__u) > (0.5 * __n * __n + __a1 02410 * (1.0 - __v + std::log(__v))))); 02411 02412 *__f++ = __a1 * __v * __param.beta(); 02413 } 02414 else 02415 while (__f != __t) 02416 { 02417 do 02418 { 02419 do 02420 { 02421 __n = _M_nd(__urng); 02422 __v = result_type(1.0) + __param._M_a2 * __n; 02423 } 02424 while (__v <= 0.0); 02425 02426 __v = __v * __v * __v; 02427 __u = __aurng(); 02428 } 02429 while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 02430 && (std::log(__u) > (0.5 * __n * __n + __a1 02431 * (1.0 - __v + std::log(__v))))); 02432 02433 do 02434 __u = __aurng(); 02435 while (__u == 0.0); 02436 02437 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha()) 02438 * __a1 * __v * __param.beta()); 02439 } 02440 } 02441 02442 template<typename _RealType, typename _CharT, typename _Traits> 02443 std::basic_ostream<_CharT, _Traits>& 02444 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02445 const gamma_distribution<_RealType>& __x) 02446 { 02447 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02448 typedef typename __ostream_type::ios_base __ios_base; 02449 02450 const typename __ios_base::fmtflags __flags = __os.flags(); 02451 const _CharT __fill = __os.fill(); 02452 const std::streamsize __precision = __os.precision(); 02453 const _CharT __space = __os.widen(' '); 02454 __os.flags(__ios_base::scientific | __ios_base::left); 02455 __os.fill(__space); 02456 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02457 02458 __os << __x.alpha() << __space << __x.beta() 02459 << __space << __x._M_nd; 02460 02461 __os.flags(__flags); 02462 __os.fill(__fill); 02463 __os.precision(__precision); 02464 return __os; 02465 } 02466 02467 template<typename _RealType, typename _CharT, typename _Traits> 02468 std::basic_istream<_CharT, _Traits>& 02469 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02470 gamma_distribution<_RealType>& __x) 02471 { 02472 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02473 typedef typename __istream_type::ios_base __ios_base; 02474 02475 const typename __ios_base::fmtflags __flags = __is.flags(); 02476 __is.flags(__ios_base::dec | __ios_base::skipws); 02477 02478 _RealType __alpha_val, __beta_val; 02479 __is >> __alpha_val >> __beta_val >> __x._M_nd; 02480 __x.param(typename gamma_distribution<_RealType>:: 02481 param_type(__alpha_val, __beta_val)); 02482 02483 __is.flags(__flags); 02484 return __is; 02485 } 02486 02487 02488 template<typename _RealType> 02489 template<typename _UniformRandomNumberGenerator> 02490 typename weibull_distribution<_RealType>::result_type 02491 weibull_distribution<_RealType>:: 02492 operator()(_UniformRandomNumberGenerator& __urng, 02493 const param_type& __p) 02494 { 02495 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02496 __aurng(__urng); 02497 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 02498 result_type(1) / __p.a()); 02499 } 02500 02501 template<typename _RealType> 02502 template<typename _ForwardIterator, 02503 typename _UniformRandomNumberGenerator> 02504 void 02505 weibull_distribution<_RealType>:: 02506 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02507 _UniformRandomNumberGenerator& __urng, 02508 const param_type& __p) 02509 { 02510 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02511 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02512 __aurng(__urng); 02513 auto __inv_a = result_type(1) / __p.a(); 02514 02515 while (__f != __t) 02516 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 02517 __inv_a); 02518 } 02519 02520 template<typename _RealType, typename _CharT, typename _Traits> 02521 std::basic_ostream<_CharT, _Traits>& 02522 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02523 const weibull_distribution<_RealType>& __x) 02524 { 02525 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02526 typedef typename __ostream_type::ios_base __ios_base; 02527 02528 const typename __ios_base::fmtflags __flags = __os.flags(); 02529 const _CharT __fill = __os.fill(); 02530 const std::streamsize __precision = __os.precision(); 02531 const _CharT __space = __os.widen(' '); 02532 __os.flags(__ios_base::scientific | __ios_base::left); 02533 __os.fill(__space); 02534 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02535 02536 __os << __x.a() << __space << __x.b(); 02537 02538 __os.flags(__flags); 02539 __os.fill(__fill); 02540 __os.precision(__precision); 02541 return __os; 02542 } 02543 02544 template<typename _RealType, typename _CharT, typename _Traits> 02545 std::basic_istream<_CharT, _Traits>& 02546 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02547 weibull_distribution<_RealType>& __x) 02548 { 02549 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02550 typedef typename __istream_type::ios_base __ios_base; 02551 02552 const typename __ios_base::fmtflags __flags = __is.flags(); 02553 __is.flags(__ios_base::dec | __ios_base::skipws); 02554 02555 _RealType __a, __b; 02556 __is >> __a >> __b; 02557 __x.param(typename weibull_distribution<_RealType>:: 02558 param_type(__a, __b)); 02559 02560 __is.flags(__flags); 02561 return __is; 02562 } 02563 02564 02565 template<typename _RealType> 02566 template<typename _UniformRandomNumberGenerator> 02567 typename extreme_value_distribution<_RealType>::result_type 02568 extreme_value_distribution<_RealType>:: 02569 operator()(_UniformRandomNumberGenerator& __urng, 02570 const param_type& __p) 02571 { 02572 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02573 __aurng(__urng); 02574 return __p.a() - __p.b() * std::log(-std::log(result_type(1) 02575 - __aurng())); 02576 } 02577 02578 template<typename _RealType> 02579 template<typename _ForwardIterator, 02580 typename _UniformRandomNumberGenerator> 02581 void 02582 extreme_value_distribution<_RealType>:: 02583 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02584 _UniformRandomNumberGenerator& __urng, 02585 const param_type& __p) 02586 { 02587 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02588 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02589 __aurng(__urng); 02590 02591 while (__f != __t) 02592 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1) 02593 - __aurng())); 02594 } 02595 02596 template<typename _RealType, typename _CharT, typename _Traits> 02597 std::basic_ostream<_CharT, _Traits>& 02598 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02599 const extreme_value_distribution<_RealType>& __x) 02600 { 02601 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02602 typedef typename __ostream_type::ios_base __ios_base; 02603 02604 const typename __ios_base::fmtflags __flags = __os.flags(); 02605 const _CharT __fill = __os.fill(); 02606 const std::streamsize __precision = __os.precision(); 02607 const _CharT __space = __os.widen(' '); 02608 __os.flags(__ios_base::scientific | __ios_base::left); 02609 __os.fill(__space); 02610 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02611 02612 __os << __x.a() << __space << __x.b(); 02613 02614 __os.flags(__flags); 02615 __os.fill(__fill); 02616 __os.precision(__precision); 02617 return __os; 02618 } 02619 02620 template<typename _RealType, typename _CharT, typename _Traits> 02621 std::basic_istream<_CharT, _Traits>& 02622 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02623 extreme_value_distribution<_RealType>& __x) 02624 { 02625 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02626 typedef typename __istream_type::ios_base __ios_base; 02627 02628 const typename __ios_base::fmtflags __flags = __is.flags(); 02629 __is.flags(__ios_base::dec | __ios_base::skipws); 02630 02631 _RealType __a, __b; 02632 __is >> __a >> __b; 02633 __x.param(typename extreme_value_distribution<_RealType>:: 02634 param_type(__a, __b)); 02635 02636 __is.flags(__flags); 02637 return __is; 02638 } 02639 02640 02641 template<typename _IntType> 02642 void 02643 discrete_distribution<_IntType>::param_type:: 02644 _M_initialize() 02645 { 02646 if (_M_prob.size() < 2) 02647 { 02648 _M_prob.clear(); 02649 return; 02650 } 02651 02652 const double __sum = std::accumulate(_M_prob.begin(), 02653 _M_prob.end(), 0.0); 02654 // Now normalize the probabilites. 02655 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), 02656 __sum); 02657 // Accumulate partial sums. 02658 _M_cp.reserve(_M_prob.size()); 02659 std::partial_sum(_M_prob.begin(), _M_prob.end(), 02660 std::back_inserter(_M_cp)); 02661 // Make sure the last cumulative probability is one. 02662 _M_cp[_M_cp.size() - 1] = 1.0; 02663 } 02664 02665 template<typename _IntType> 02666 template<typename _Func> 02667 discrete_distribution<_IntType>::param_type:: 02668 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw) 02669 : _M_prob(), _M_cp() 02670 { 02671 const size_t __n = __nw == 0 ? 1 : __nw; 02672 const double __delta = (__xmax - __xmin) / __n; 02673 02674 _M_prob.reserve(__n); 02675 for (size_t __k = 0; __k < __nw; ++__k) 02676 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta)); 02677 02678 _M_initialize(); 02679 } 02680 02681 template<typename _IntType> 02682 template<typename _UniformRandomNumberGenerator> 02683 typename discrete_distribution<_IntType>::result_type 02684 discrete_distribution<_IntType>:: 02685 operator()(_UniformRandomNumberGenerator& __urng, 02686 const param_type& __param) 02687 { 02688 if (__param._M_cp.empty()) 02689 return result_type(0); 02690 02691 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02692 __aurng(__urng); 02693 02694 const double __p = __aurng(); 02695 auto __pos = std::lower_bound(__param._M_cp.begin(), 02696 __param._M_cp.end(), __p); 02697 02698 return __pos - __param._M_cp.begin(); 02699 } 02700 02701 template<typename _IntType> 02702 template<typename _ForwardIterator, 02703 typename _UniformRandomNumberGenerator> 02704 void 02705 discrete_distribution<_IntType>:: 02706 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02707 _UniformRandomNumberGenerator& __urng, 02708 const param_type& __param) 02709 { 02710 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02711 02712 if (__param._M_cp.empty()) 02713 { 02714 while (__f != __t) 02715 *__f++ = result_type(0); 02716 return; 02717 } 02718 02719 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02720 __aurng(__urng); 02721 02722 while (__f != __t) 02723 { 02724 const double __p = __aurng(); 02725 auto __pos = std::lower_bound(__param._M_cp.begin(), 02726 __param._M_cp.end(), __p); 02727 02728 *__f++ = __pos - __param._M_cp.begin(); 02729 } 02730 } 02731 02732 template<typename _IntType, typename _CharT, typename _Traits> 02733 std::basic_ostream<_CharT, _Traits>& 02734 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02735 const discrete_distribution<_IntType>& __x) 02736 { 02737 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02738 typedef typename __ostream_type::ios_base __ios_base; 02739 02740 const typename __ios_base::fmtflags __flags = __os.flags(); 02741 const _CharT __fill = __os.fill(); 02742 const std::streamsize __precision = __os.precision(); 02743 const _CharT __space = __os.widen(' '); 02744 __os.flags(__ios_base::scientific | __ios_base::left); 02745 __os.fill(__space); 02746 __os.precision(std::numeric_limits<double>::max_digits10); 02747 02748 std::vector<double> __prob = __x.probabilities(); 02749 __os << __prob.size(); 02750 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit) 02751 __os << __space << *__dit; 02752 02753 __os.flags(__flags); 02754 __os.fill(__fill); 02755 __os.precision(__precision); 02756 return __os; 02757 } 02758 02759 template<typename _IntType, typename _CharT, typename _Traits> 02760 std::basic_istream<_CharT, _Traits>& 02761 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02762 discrete_distribution<_IntType>& __x) 02763 { 02764 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02765 typedef typename __istream_type::ios_base __ios_base; 02766 02767 const typename __ios_base::fmtflags __flags = __is.flags(); 02768 __is.flags(__ios_base::dec | __ios_base::skipws); 02769 02770 size_t __n; 02771 __is >> __n; 02772 02773 std::vector<double> __prob_vec; 02774 __prob_vec.reserve(__n); 02775 for (; __n != 0; --__n) 02776 { 02777 double __prob; 02778 __is >> __prob; 02779 __prob_vec.push_back(__prob); 02780 } 02781 02782 __x.param(typename discrete_distribution<_IntType>:: 02783 param_type(__prob_vec.begin(), __prob_vec.end())); 02784 02785 __is.flags(__flags); 02786 return __is; 02787 } 02788 02789 02790 template<typename _RealType> 02791 void 02792 piecewise_constant_distribution<_RealType>::param_type:: 02793 _M_initialize() 02794 { 02795 if (_M_int.size() < 2 02796 || (_M_int.size() == 2 02797 && _M_int[0] == _RealType(0) 02798 && _M_int[1] == _RealType(1))) 02799 { 02800 _M_int.clear(); 02801 _M_den.clear(); 02802 return; 02803 } 02804 02805 const double __sum = std::accumulate(_M_den.begin(), 02806 _M_den.end(), 0.0); 02807 02808 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 02809 __sum); 02810 02811 _M_cp.reserve(_M_den.size()); 02812 std::partial_sum(_M_den.begin(), _M_den.end(), 02813 std::back_inserter(_M_cp)); 02814 02815 // Make sure the last cumulative probability is one. 02816 _M_cp[_M_cp.size() - 1] = 1.0; 02817 02818 for (size_t __k = 0; __k < _M_den.size(); ++__k) 02819 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k]; 02820 } 02821 02822 template<typename _RealType> 02823 template<typename _InputIteratorB, typename _InputIteratorW> 02824 piecewise_constant_distribution<_RealType>::param_type:: 02825 param_type(_InputIteratorB __bbegin, 02826 _InputIteratorB __bend, 02827 _InputIteratorW __wbegin) 02828 : _M_int(), _M_den(), _M_cp() 02829 { 02830 if (__bbegin != __bend) 02831 { 02832 for (;;) 02833 { 02834 _M_int.push_back(*__bbegin); 02835 ++__bbegin; 02836 if (__bbegin == __bend) 02837 break; 02838 02839 _M_den.push_back(*__wbegin); 02840 ++__wbegin; 02841 } 02842 } 02843 02844 _M_initialize(); 02845 } 02846 02847 template<typename _RealType> 02848 template<typename _Func> 02849 piecewise_constant_distribution<_RealType>::param_type:: 02850 param_type(initializer_list<_RealType> __bl, _Func __fw) 02851 : _M_int(), _M_den(), _M_cp() 02852 { 02853 _M_int.reserve(__bl.size()); 02854 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 02855 _M_int.push_back(*__biter); 02856 02857 _M_den.reserve(_M_int.size() - 1); 02858 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 02859 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k]))); 02860 02861 _M_initialize(); 02862 } 02863 02864 template<typename _RealType> 02865 template<typename _Func> 02866 piecewise_constant_distribution<_RealType>::param_type:: 02867 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 02868 : _M_int(), _M_den(), _M_cp() 02869 { 02870 const size_t __n = __nw == 0 ? 1 : __nw; 02871 const _RealType __delta = (__xmax - __xmin) / __n; 02872 02873 _M_int.reserve(__n + 1); 02874 for (size_t __k = 0; __k <= __nw; ++__k) 02875 _M_int.push_back(__xmin + __k * __delta); 02876 02877 _M_den.reserve(__n); 02878 for (size_t __k = 0; __k < __nw; ++__k) 02879 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta)); 02880 02881 _M_initialize(); 02882 } 02883 02884 template<typename _RealType> 02885 template<typename _UniformRandomNumberGenerator> 02886 typename piecewise_constant_distribution<_RealType>::result_type 02887 piecewise_constant_distribution<_RealType>:: 02888 operator()(_UniformRandomNumberGenerator& __urng, 02889 const param_type& __param) 02890 { 02891 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02892 __aurng(__urng); 02893 02894 const double __p = __aurng(); 02895 if (__param._M_cp.empty()) 02896 return __p; 02897 02898 auto __pos = std::lower_bound(__param._M_cp.begin(), 02899 __param._M_cp.end(), __p); 02900 const size_t __i = __pos - __param._M_cp.begin(); 02901 02902 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 02903 02904 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i]; 02905 } 02906 02907 template<typename _RealType> 02908 template<typename _ForwardIterator, 02909 typename _UniformRandomNumberGenerator> 02910 void 02911 piecewise_constant_distribution<_RealType>:: 02912 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02913 _UniformRandomNumberGenerator& __urng, 02914 const param_type& __param) 02915 { 02916 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02917 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02918 __aurng(__urng); 02919 02920 if (__param._M_cp.empty()) 02921 { 02922 while (__f != __t) 02923 *__f++ = __aurng(); 02924 return; 02925 } 02926 02927 while (__f != __t) 02928 { 02929 const double __p = __aurng(); 02930 02931 auto __pos = std::lower_bound(__param._M_cp.begin(), 02932 __param._M_cp.end(), __p); 02933 const size_t __i = __pos - __param._M_cp.begin(); 02934 02935 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 02936 02937 *__f++ = (__param._M_int[__i] 02938 + (__p - __pref) / __param._M_den[__i]); 02939 } 02940 } 02941 02942 template<typename _RealType, typename _CharT, typename _Traits> 02943 std::basic_ostream<_CharT, _Traits>& 02944 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02945 const piecewise_constant_distribution<_RealType>& __x) 02946 { 02947 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02948 typedef typename __ostream_type::ios_base __ios_base; 02949 02950 const typename __ios_base::fmtflags __flags = __os.flags(); 02951 const _CharT __fill = __os.fill(); 02952 const std::streamsize __precision = __os.precision(); 02953 const _CharT __space = __os.widen(' '); 02954 __os.flags(__ios_base::scientific | __ios_base::left); 02955 __os.fill(__space); 02956 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02957 02958 std::vector<_RealType> __int = __x.intervals(); 02959 __os << __int.size() - 1; 02960 02961 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 02962 __os << __space << *__xit; 02963 02964 std::vector<double> __den = __x.densities(); 02965 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 02966 __os << __space << *__dit; 02967 02968 __os.flags(__flags); 02969 __os.fill(__fill); 02970 __os.precision(__precision); 02971 return __os; 02972 } 02973 02974 template<typename _RealType, typename _CharT, typename _Traits> 02975 std::basic_istream<_CharT, _Traits>& 02976 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02977 piecewise_constant_distribution<_RealType>& __x) 02978 { 02979 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02980 typedef typename __istream_type::ios_base __ios_base; 02981 02982 const typename __ios_base::fmtflags __flags = __is.flags(); 02983 __is.flags(__ios_base::dec | __ios_base::skipws); 02984 02985 size_t __n; 02986 __is >> __n; 02987 02988 std::vector<_RealType> __int_vec; 02989 __int_vec.reserve(__n + 1); 02990 for (size_t __i = 0; __i <= __n; ++__i) 02991 { 02992 _RealType __int; 02993 __is >> __int; 02994 __int_vec.push_back(__int); 02995 } 02996 02997 std::vector<double> __den_vec; 02998 __den_vec.reserve(__n); 02999 for (size_t __i = 0; __i < __n; ++__i) 03000 { 03001 double __den; 03002 __is >> __den; 03003 __den_vec.push_back(__den); 03004 } 03005 03006 __x.param(typename piecewise_constant_distribution<_RealType>:: 03007 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 03008 03009 __is.flags(__flags); 03010 return __is; 03011 } 03012 03013 03014 template<typename _RealType> 03015 void 03016 piecewise_linear_distribution<_RealType>::param_type:: 03017 _M_initialize() 03018 { 03019 if (_M_int.size() < 2 03020 || (_M_int.size() == 2 03021 && _M_int[0] == _RealType(0) 03022 && _M_int[1] == _RealType(1) 03023 && _M_den[0] == _M_den[1])) 03024 { 03025 _M_int.clear(); 03026 _M_den.clear(); 03027 return; 03028 } 03029 03030 double __sum = 0.0; 03031 _M_cp.reserve(_M_int.size() - 1); 03032 _M_m.reserve(_M_int.size() - 1); 03033 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 03034 { 03035 const _RealType __delta = _M_int[__k + 1] - _M_int[__k]; 03036 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta; 03037 _M_cp.push_back(__sum); 03038 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta); 03039 } 03040 03041 // Now normalize the densities... 03042 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 03043 __sum); 03044 // ... and partial sums... 03045 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum); 03046 // ... and slopes. 03047 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum); 03048 03049 // Make sure the last cumulative probablility is one. 03050 _M_cp[_M_cp.size() - 1] = 1.0; 03051 } 03052 03053 template<typename _RealType> 03054 template<typename _InputIteratorB, typename _InputIteratorW> 03055 piecewise_linear_distribution<_RealType>::param_type:: 03056 param_type(_InputIteratorB __bbegin, 03057 _InputIteratorB __bend, 03058 _InputIteratorW __wbegin) 03059 : _M_int(), _M_den(), _M_cp(), _M_m() 03060 { 03061 for (; __bbegin != __bend; ++__bbegin, ++__wbegin) 03062 { 03063 _M_int.push_back(*__bbegin); 03064 _M_den.push_back(*__wbegin); 03065 } 03066 03067 _M_initialize(); 03068 } 03069 03070 template<typename _RealType> 03071 template<typename _Func> 03072 piecewise_linear_distribution<_RealType>::param_type:: 03073 param_type(initializer_list<_RealType> __bl, _Func __fw) 03074 : _M_int(), _M_den(), _M_cp(), _M_m() 03075 { 03076 _M_int.reserve(__bl.size()); 03077 _M_den.reserve(__bl.size()); 03078 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 03079 { 03080 _M_int.push_back(*__biter); 03081 _M_den.push_back(__fw(*__biter)); 03082 } 03083 03084 _M_initialize(); 03085 } 03086 03087 template<typename _RealType> 03088 template<typename _Func> 03089 piecewise_linear_distribution<_RealType>::param_type:: 03090 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 03091 : _M_int(), _M_den(), _M_cp(), _M_m() 03092 { 03093 const size_t __n = __nw == 0 ? 1 : __nw; 03094 const _RealType __delta = (__xmax - __xmin) / __n; 03095 03096 _M_int.reserve(__n + 1); 03097 _M_den.reserve(__n + 1); 03098 for (size_t __k = 0; __k <= __nw; ++__k) 03099 { 03100 _M_int.push_back(__xmin + __k * __delta); 03101 _M_den.push_back(__fw(_M_int[__k] + __delta)); 03102 } 03103 03104 _M_initialize(); 03105 } 03106 03107 template<typename _RealType> 03108 template<typename _UniformRandomNumberGenerator> 03109 typename piecewise_linear_distribution<_RealType>::result_type 03110 piecewise_linear_distribution<_RealType>:: 03111 operator()(_UniformRandomNumberGenerator& __urng, 03112 const param_type& __param) 03113 { 03114 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 03115 __aurng(__urng); 03116 03117 const double __p = __aurng(); 03118 if (__param._M_cp.empty()) 03119 return __p; 03120 03121 auto __pos = std::lower_bound(__param._M_cp.begin(), 03122 __param._M_cp.end(), __p); 03123 const size_t __i = __pos - __param._M_cp.begin(); 03124 03125 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 03126 03127 const double __a = 0.5 * __param._M_m[__i]; 03128 const double __b = __param._M_den[__i]; 03129 const double __cm = __p - __pref; 03130 03131 _RealType __x = __param._M_int[__i]; 03132 if (__a == 0) 03133 __x += __cm / __b; 03134 else 03135 { 03136 const double __d = __b * __b + 4.0 * __a * __cm; 03137 __x += 0.5 * (std::sqrt(__d) - __b) / __a; 03138 } 03139 03140 return __x; 03141 } 03142 03143 template<typename _RealType> 03144 template<typename _ForwardIterator, 03145 typename _UniformRandomNumberGenerator> 03146 void 03147 piecewise_linear_distribution<_RealType>:: 03148 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 03149 _UniformRandomNumberGenerator& __urng, 03150 const param_type& __param) 03151 { 03152 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 03153 // We could duplicate everything from operator()... 03154 while (__f != __t) 03155 *__f++ = this->operator()(__urng, __param); 03156 } 03157 03158 template<typename _RealType, typename _CharT, typename _Traits> 03159 std::basic_ostream<_CharT, _Traits>& 03160 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 03161 const piecewise_linear_distribution<_RealType>& __x) 03162 { 03163 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 03164 typedef typename __ostream_type::ios_base __ios_base; 03165 03166 const typename __ios_base::fmtflags __flags = __os.flags(); 03167 const _CharT __fill = __os.fill(); 03168 const std::streamsize __precision = __os.precision(); 03169 const _CharT __space = __os.widen(' '); 03170 __os.flags(__ios_base::scientific | __ios_base::left); 03171 __os.fill(__space); 03172 __os.precision(std::numeric_limits<_RealType>::max_digits10); 03173 03174 std::vector<_RealType> __int = __x.intervals(); 03175 __os << __int.size() - 1; 03176 03177 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 03178 __os << __space << *__xit; 03179 03180 std::vector<double> __den = __x.densities(); 03181 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 03182 __os << __space << *__dit; 03183 03184 __os.flags(__flags); 03185 __os.fill(__fill); 03186 __os.precision(__precision); 03187 return __os; 03188 } 03189 03190 template<typename _RealType, typename _CharT, typename _Traits> 03191 std::basic_istream<_CharT, _Traits>& 03192 operator>>(std::basic_istream<_CharT, _Traits>& __is, 03193 piecewise_linear_distribution<_RealType>& __x) 03194 { 03195 typedef std::basic_istream<_CharT, _Traits> __istream_type; 03196 typedef typename __istream_type::ios_base __ios_base; 03197 03198 const typename __ios_base::fmtflags __flags = __is.flags(); 03199 __is.flags(__ios_base::dec | __ios_base::skipws); 03200 03201 size_t __n; 03202 __is >> __n; 03203 03204 std::vector<_RealType> __int_vec; 03205 __int_vec.reserve(__n + 1); 03206 for (size_t __i = 0; __i <= __n; ++__i) 03207 { 03208 _RealType __int; 03209 __is >> __int; 03210 __int_vec.push_back(__int); 03211 } 03212 03213 std::vector<double> __den_vec; 03214 __den_vec.reserve(__n + 1); 03215 for (size_t __i = 0; __i <= __n; ++__i) 03216 { 03217 double __den; 03218 __is >> __den; 03219 __den_vec.push_back(__den); 03220 } 03221 03222 __x.param(typename piecewise_linear_distribution<_RealType>:: 03223 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 03224 03225 __is.flags(__flags); 03226 return __is; 03227 } 03228 03229 03230 template<typename _IntType> 03231 seed_seq::seed_seq(std::initializer_list<_IntType> __il) 03232 { 03233 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter) 03234 _M_v.push_back(__detail::__mod<result_type, 03235 __detail::_Shift<result_type, 32>::__value>(*__iter)); 03236 } 03237 03238 template<typename _InputIterator> 03239 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end) 03240 { 03241 for (_InputIterator __iter = __begin; __iter != __end; ++__iter) 03242 _M_v.push_back(__detail::__mod<result_type, 03243 __detail::_Shift<result_type, 32>::__value>(*__iter)); 03244 } 03245 03246 template<typename _RandomAccessIterator> 03247 void 03248 seed_seq::generate(_RandomAccessIterator __begin, 03249 _RandomAccessIterator __end) 03250 { 03251 typedef typename iterator_traits<_RandomAccessIterator>::value_type 03252 _Type; 03253 03254 if (__begin == __end) 03255 return; 03256 03257 std::fill(__begin, __end, _Type(0x8b8b8b8bu)); 03258 03259 const size_t __n = __end - __begin; 03260 const size_t __s = _M_v.size(); 03261 const size_t __t = (__n >= 623) ? 11 03262 : (__n >= 68) ? 7 03263 : (__n >= 39) ? 5 03264 : (__n >= 7) ? 3 03265 : (__n - 1) / 2; 03266 const size_t __p = (__n - __t) / 2; 03267 const size_t __q = __p + __t; 03268 const size_t __m = std::max(size_t(__s + 1), __n); 03269 03270 for (size_t __k = 0; __k < __m; ++__k) 03271 { 03272 _Type __arg = (__begin[__k % __n] 03273 ^ __begin[(__k + __p) % __n] 03274 ^ __begin[(__k - 1) % __n]); 03275 _Type __r1 = __arg ^ (__arg >> 27); 03276 __r1 = __detail::__mod<_Type, 03277 __detail::_Shift<_Type, 32>::__value>(1664525u * __r1); 03278 _Type __r2 = __r1; 03279 if (__k == 0) 03280 __r2 += __s; 03281 else if (__k <= __s) 03282 __r2 += __k % __n + _M_v[__k - 1]; 03283 else 03284 __r2 += __k % __n; 03285 __r2 = __detail::__mod<_Type, 03286 __detail::_Shift<_Type, 32>::__value>(__r2); 03287 __begin[(__k + __p) % __n] += __r1; 03288 __begin[(__k + __q) % __n] += __r2; 03289 __begin[__k % __n] = __r2; 03290 } 03291 03292 for (size_t __k = __m; __k < __m + __n; ++__k) 03293 { 03294 _Type __arg = (__begin[__k % __n] 03295 + __begin[(__k + __p) % __n] 03296 + __begin[(__k - 1) % __n]); 03297 _Type __r3 = __arg ^ (__arg >> 27); 03298 __r3 = __detail::__mod<_Type, 03299 __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3); 03300 _Type __r4 = __r3 - __k % __n; 03301 __r4 = __detail::__mod<_Type, 03302 __detail::_Shift<_Type, 32>::__value>(__r4); 03303 __begin[(__k + __p) % __n] ^= __r3; 03304 __begin[(__k + __q) % __n] ^= __r4; 03305 __begin[__k % __n] = __r4; 03306 } 03307 } 03308 03309 template<typename _RealType, size_t __bits, 03310 typename _UniformRandomNumberGenerator> 03311 _RealType 03312 generate_canonical(_UniformRandomNumberGenerator& __urng) 03313 { 03314 static_assert(std::is_floating_point<_RealType>::value, 03315 "template argument not a floating point type"); 03316 03317 const size_t __b 03318 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits), 03319 __bits); 03320 const long double __r = static_cast<long double>(__urng.max()) 03321 - static_cast<long double>(__urng.min()) + 1.0L; 03322 const size_t __log2r = std::log(__r) / std::log(2.0L); 03323 size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r); 03324 _RealType __sum = _RealType(0); 03325 _RealType __tmp = _RealType(1); 03326 for (; __k != 0; --__k) 03327 { 03328 __sum += _RealType(__urng() - __urng.min()) * __tmp; 03329 __tmp *= __r; 03330 } 03331 return __sum / __tmp; 03332 } 03333 03334 _GLIBCXX_END_NAMESPACE_VERSION 03335 } // namespace 03336 03337 #endif