libstdc++
|
00001 // random number generation (out of line) -*- C++ -*- 00002 00003 // Copyright (C) 2009-2019 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 _GLIBCXX_BEGIN_NAMESPACE_VERSION 00038 00039 /* 00040 * (Further) implementation-space details. 00041 */ 00042 namespace __detail 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 } // namespace __detail 00093 00094 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00095 constexpr _UIntType 00096 linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier; 00097 00098 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00099 constexpr _UIntType 00100 linear_congruential_engine<_UIntType, __a, __c, __m>::increment; 00101 00102 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00103 constexpr _UIntType 00104 linear_congruential_engine<_UIntType, __a, __c, __m>::modulus; 00105 00106 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00107 constexpr _UIntType 00108 linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed; 00109 00110 /** 00111 * Seeds the LCR with integral value @p __s, adjusted so that the 00112 * ring identity is never a member of the convergence set. 00113 */ 00114 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00115 void 00116 linear_congruential_engine<_UIntType, __a, __c, __m>:: 00117 seed(result_type __s) 00118 { 00119 if ((__detail::__mod<_UIntType, __m>(__c) == 0) 00120 && (__detail::__mod<_UIntType, __m>(__s) == 0)) 00121 _M_x = 1; 00122 else 00123 _M_x = __detail::__mod<_UIntType, __m>(__s); 00124 } 00125 00126 /** 00127 * Seeds the LCR engine with a value generated by @p __q. 00128 */ 00129 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 00130 template<typename _Sseq> 00131 auto 00132 linear_congruential_engine<_UIntType, __a, __c, __m>:: 00133 seed(_Sseq& __q) 00134 -> _If_seed_seq<_Sseq> 00135 { 00136 const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits 00137 : std::__lg(__m); 00138 const _UIntType __k = (__k0 + 31) / 32; 00139 uint_least32_t __arr[__k + 3]; 00140 __q.generate(__arr + 0, __arr + __k + 3); 00141 _UIntType __factor = 1u; 00142 _UIntType __sum = 0u; 00143 for (size_t __j = 0; __j < __k; ++__j) 00144 { 00145 __sum += __arr[__j + 3] * __factor; 00146 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00147 } 00148 seed(__sum); 00149 } 00150 00151 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 00152 typename _CharT, typename _Traits> 00153 std::basic_ostream<_CharT, _Traits>& 00154 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00155 const linear_congruential_engine<_UIntType, 00156 __a, __c, __m>& __lcr) 00157 { 00158 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00159 typedef typename __ostream_type::ios_base __ios_base; 00160 00161 const typename __ios_base::fmtflags __flags = __os.flags(); 00162 const _CharT __fill = __os.fill(); 00163 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00164 __os.fill(__os.widen(' ')); 00165 00166 __os << __lcr._M_x; 00167 00168 __os.flags(__flags); 00169 __os.fill(__fill); 00170 return __os; 00171 } 00172 00173 template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 00174 typename _CharT, typename _Traits> 00175 std::basic_istream<_CharT, _Traits>& 00176 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00177 linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr) 00178 { 00179 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00180 typedef typename __istream_type::ios_base __ios_base; 00181 00182 const typename __ios_base::fmtflags __flags = __is.flags(); 00183 __is.flags(__ios_base::dec); 00184 00185 __is >> __lcr._M_x; 00186 00187 __is.flags(__flags); 00188 return __is; 00189 } 00190 00191 00192 template<typename _UIntType, 00193 size_t __w, size_t __n, size_t __m, size_t __r, 00194 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00195 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00196 _UIntType __f> 00197 constexpr size_t 00198 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00199 __s, __b, __t, __c, __l, __f>::word_size; 00200 00201 template<typename _UIntType, 00202 size_t __w, size_t __n, size_t __m, size_t __r, 00203 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00204 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00205 _UIntType __f> 00206 constexpr size_t 00207 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00208 __s, __b, __t, __c, __l, __f>::state_size; 00209 00210 template<typename _UIntType, 00211 size_t __w, size_t __n, size_t __m, size_t __r, 00212 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00213 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00214 _UIntType __f> 00215 constexpr size_t 00216 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00217 __s, __b, __t, __c, __l, __f>::shift_size; 00218 00219 template<typename _UIntType, 00220 size_t __w, size_t __n, size_t __m, size_t __r, 00221 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00222 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00223 _UIntType __f> 00224 constexpr size_t 00225 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00226 __s, __b, __t, __c, __l, __f>::mask_bits; 00227 00228 template<typename _UIntType, 00229 size_t __w, size_t __n, size_t __m, size_t __r, 00230 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00231 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00232 _UIntType __f> 00233 constexpr _UIntType 00234 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00235 __s, __b, __t, __c, __l, __f>::xor_mask; 00236 00237 template<typename _UIntType, 00238 size_t __w, size_t __n, size_t __m, size_t __r, 00239 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00240 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00241 _UIntType __f> 00242 constexpr size_t 00243 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00244 __s, __b, __t, __c, __l, __f>::tempering_u; 00245 00246 template<typename _UIntType, 00247 size_t __w, size_t __n, size_t __m, size_t __r, 00248 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00249 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00250 _UIntType __f> 00251 constexpr _UIntType 00252 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00253 __s, __b, __t, __c, __l, __f>::tempering_d; 00254 00255 template<typename _UIntType, 00256 size_t __w, size_t __n, size_t __m, size_t __r, 00257 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00258 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00259 _UIntType __f> 00260 constexpr size_t 00261 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00262 __s, __b, __t, __c, __l, __f>::tempering_s; 00263 00264 template<typename _UIntType, 00265 size_t __w, size_t __n, size_t __m, size_t __r, 00266 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00267 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00268 _UIntType __f> 00269 constexpr _UIntType 00270 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00271 __s, __b, __t, __c, __l, __f>::tempering_b; 00272 00273 template<typename _UIntType, 00274 size_t __w, size_t __n, size_t __m, size_t __r, 00275 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00276 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00277 _UIntType __f> 00278 constexpr size_t 00279 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00280 __s, __b, __t, __c, __l, __f>::tempering_t; 00281 00282 template<typename _UIntType, 00283 size_t __w, size_t __n, size_t __m, size_t __r, 00284 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00285 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00286 _UIntType __f> 00287 constexpr _UIntType 00288 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00289 __s, __b, __t, __c, __l, __f>::tempering_c; 00290 00291 template<typename _UIntType, 00292 size_t __w, size_t __n, size_t __m, size_t __r, 00293 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00294 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00295 _UIntType __f> 00296 constexpr size_t 00297 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00298 __s, __b, __t, __c, __l, __f>::tempering_l; 00299 00300 template<typename _UIntType, 00301 size_t __w, size_t __n, size_t __m, size_t __r, 00302 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00303 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00304 _UIntType __f> 00305 constexpr _UIntType 00306 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00307 __s, __b, __t, __c, __l, __f>:: 00308 initialization_multiplier; 00309 00310 template<typename _UIntType, 00311 size_t __w, size_t __n, size_t __m, size_t __r, 00312 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00313 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00314 _UIntType __f> 00315 constexpr _UIntType 00316 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00317 __s, __b, __t, __c, __l, __f>::default_seed; 00318 00319 template<typename _UIntType, 00320 size_t __w, size_t __n, size_t __m, size_t __r, 00321 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00322 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00323 _UIntType __f> 00324 void 00325 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00326 __s, __b, __t, __c, __l, __f>:: 00327 seed(result_type __sd) 00328 { 00329 _M_x[0] = __detail::__mod<_UIntType, 00330 __detail::_Shift<_UIntType, __w>::__value>(__sd); 00331 00332 for (size_t __i = 1; __i < state_size; ++__i) 00333 { 00334 _UIntType __x = _M_x[__i - 1]; 00335 __x ^= __x >> (__w - 2); 00336 __x *= __f; 00337 __x += __detail::__mod<_UIntType, __n>(__i); 00338 _M_x[__i] = __detail::__mod<_UIntType, 00339 __detail::_Shift<_UIntType, __w>::__value>(__x); 00340 } 00341 _M_p = state_size; 00342 } 00343 00344 template<typename _UIntType, 00345 size_t __w, size_t __n, size_t __m, size_t __r, 00346 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00347 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00348 _UIntType __f> 00349 template<typename _Sseq> 00350 auto 00351 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00352 __s, __b, __t, __c, __l, __f>:: 00353 seed(_Sseq& __q) 00354 -> _If_seed_seq<_Sseq> 00355 { 00356 const _UIntType __upper_mask = (~_UIntType()) << __r; 00357 const size_t __k = (__w + 31) / 32; 00358 uint_least32_t __arr[__n * __k]; 00359 __q.generate(__arr + 0, __arr + __n * __k); 00360 00361 bool __zero = true; 00362 for (size_t __i = 0; __i < state_size; ++__i) 00363 { 00364 _UIntType __factor = 1u; 00365 _UIntType __sum = 0u; 00366 for (size_t __j = 0; __j < __k; ++__j) 00367 { 00368 __sum += __arr[__k * __i + __j] * __factor; 00369 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00370 } 00371 _M_x[__i] = __detail::__mod<_UIntType, 00372 __detail::_Shift<_UIntType, __w>::__value>(__sum); 00373 00374 if (__zero) 00375 { 00376 if (__i == 0) 00377 { 00378 if ((_M_x[0] & __upper_mask) != 0u) 00379 __zero = false; 00380 } 00381 else if (_M_x[__i] != 0u) 00382 __zero = false; 00383 } 00384 } 00385 if (__zero) 00386 _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value; 00387 _M_p = state_size; 00388 } 00389 00390 template<typename _UIntType, size_t __w, 00391 size_t __n, size_t __m, size_t __r, 00392 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00393 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00394 _UIntType __f> 00395 void 00396 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00397 __s, __b, __t, __c, __l, __f>:: 00398 _M_gen_rand(void) 00399 { 00400 const _UIntType __upper_mask = (~_UIntType()) << __r; 00401 const _UIntType __lower_mask = ~__upper_mask; 00402 00403 for (size_t __k = 0; __k < (__n - __m); ++__k) 00404 { 00405 _UIntType __y = ((_M_x[__k] & __upper_mask) 00406 | (_M_x[__k + 1] & __lower_mask)); 00407 _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) 00408 ^ ((__y & 0x01) ? __a : 0)); 00409 } 00410 00411 for (size_t __k = (__n - __m); __k < (__n - 1); ++__k) 00412 { 00413 _UIntType __y = ((_M_x[__k] & __upper_mask) 00414 | (_M_x[__k + 1] & __lower_mask)); 00415 _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) 00416 ^ ((__y & 0x01) ? __a : 0)); 00417 } 00418 00419 _UIntType __y = ((_M_x[__n - 1] & __upper_mask) 00420 | (_M_x[0] & __lower_mask)); 00421 _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) 00422 ^ ((__y & 0x01) ? __a : 0)); 00423 _M_p = 0; 00424 } 00425 00426 template<typename _UIntType, size_t __w, 00427 size_t __n, size_t __m, size_t __r, 00428 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00429 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00430 _UIntType __f> 00431 void 00432 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00433 __s, __b, __t, __c, __l, __f>:: 00434 discard(unsigned long long __z) 00435 { 00436 while (__z > state_size - _M_p) 00437 { 00438 __z -= state_size - _M_p; 00439 _M_gen_rand(); 00440 } 00441 _M_p += __z; 00442 } 00443 00444 template<typename _UIntType, size_t __w, 00445 size_t __n, size_t __m, size_t __r, 00446 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00447 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00448 _UIntType __f> 00449 typename 00450 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00451 __s, __b, __t, __c, __l, __f>::result_type 00452 mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 00453 __s, __b, __t, __c, __l, __f>:: 00454 operator()() 00455 { 00456 // Reload the vector - cost is O(n) amortized over n calls. 00457 if (_M_p >= state_size) 00458 _M_gen_rand(); 00459 00460 // Calculate o(x(i)). 00461 result_type __z = _M_x[_M_p++]; 00462 __z ^= (__z >> __u) & __d; 00463 __z ^= (__z << __s) & __b; 00464 __z ^= (__z << __t) & __c; 00465 __z ^= (__z >> __l); 00466 00467 return __z; 00468 } 00469 00470 template<typename _UIntType, size_t __w, 00471 size_t __n, size_t __m, size_t __r, 00472 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00473 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00474 _UIntType __f, typename _CharT, typename _Traits> 00475 std::basic_ostream<_CharT, _Traits>& 00476 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 00477 const mersenne_twister_engine<_UIntType, __w, __n, __m, 00478 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 00479 { 00480 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 00481 typedef typename __ostream_type::ios_base __ios_base; 00482 00483 const typename __ios_base::fmtflags __flags = __os.flags(); 00484 const _CharT __fill = __os.fill(); 00485 const _CharT __space = __os.widen(' '); 00486 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 00487 __os.fill(__space); 00488 00489 for (size_t __i = 0; __i < __n; ++__i) 00490 __os << __x._M_x[__i] << __space; 00491 __os << __x._M_p; 00492 00493 __os.flags(__flags); 00494 __os.fill(__fill); 00495 return __os; 00496 } 00497 00498 template<typename _UIntType, size_t __w, 00499 size_t __n, size_t __m, size_t __r, 00500 _UIntType __a, size_t __u, _UIntType __d, size_t __s, 00501 _UIntType __b, size_t __t, _UIntType __c, size_t __l, 00502 _UIntType __f, typename _CharT, typename _Traits> 00503 std::basic_istream<_CharT, _Traits>& 00504 operator>>(std::basic_istream<_CharT, _Traits>& __is, 00505 mersenne_twister_engine<_UIntType, __w, __n, __m, 00506 __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 00507 { 00508 typedef std::basic_istream<_CharT, _Traits> __istream_type; 00509 typedef typename __istream_type::ios_base __ios_base; 00510 00511 const typename __ios_base::fmtflags __flags = __is.flags(); 00512 __is.flags(__ios_base::dec | __ios_base::skipws); 00513 00514 for (size_t __i = 0; __i < __n; ++__i) 00515 __is >> __x._M_x[__i]; 00516 __is >> __x._M_p; 00517 00518 __is.flags(__flags); 00519 return __is; 00520 } 00521 00522 00523 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00524 constexpr size_t 00525 subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size; 00526 00527 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00528 constexpr size_t 00529 subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag; 00530 00531 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00532 constexpr size_t 00533 subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag; 00534 00535 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00536 constexpr _UIntType 00537 subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed; 00538 00539 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00540 void 00541 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 00542 seed(result_type __value) 00543 { 00544 std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u> 00545 __lcg(__value == 0u ? default_seed : __value); 00546 00547 const size_t __n = (__w + 31) / 32; 00548 00549 for (size_t __i = 0; __i < long_lag; ++__i) 00550 { 00551 _UIntType __sum = 0u; 00552 _UIntType __factor = 1u; 00553 for (size_t __j = 0; __j < __n; ++__j) 00554 { 00555 __sum += __detail::__mod<uint_least32_t, 00556 __detail::_Shift<uint_least32_t, 32>::__value> 00557 (__lcg()) * __factor; 00558 __factor *= __detail::_Shift<_UIntType, 32>::__value; 00559 } 00560 _M_x[__i] = __detail::__mod<_UIntType, 00561 __detail::_Shift<_UIntType, __w>::__value>(__sum); 00562 } 00563 _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 00564 _M_p = 0; 00565 } 00566 00567 template<typename _UIntType, size_t __w, size_t __s, size_t __r> 00568 template<typename _Sseq> 00569 auto 00570 subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 00571 seed(_Sseq& __q) 00572 -> _If_seed_seq<_Sseq> 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<double>(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 // 1 / 78 01308 const double __178 = 0.0128205128205128205128205128205128L; 01309 // e^(1 / 78) 01310 const double __e178 = 1.0129030479320018583185514777512983L; 01311 const double __c5 = __c4 + __e178; 01312 const double __c = __param._M_cb + __c5; 01313 const double __2cx = 2 * (2 * __m + __param._M_d); 01314 01315 bool __reject = true; 01316 do 01317 { 01318 const double __u = __c * __aurng(); 01319 const double __e = -std::log(1.0 - __aurng()); 01320 01321 double __w = 0.0; 01322 01323 if (__u <= __c1) 01324 { 01325 const double __n = _M_nd(__urng); 01326 const double __y = -std::abs(__n) * __param._M_sm - 1; 01327 __x = std::floor(__y); 01328 __w = -__n * __n / 2; 01329 if (__x < -__m) 01330 continue; 01331 } 01332 else if (__u <= __c2) 01333 { 01334 const double __n = _M_nd(__urng); 01335 const double __y = 1 + std::abs(__n) * __param._M_scx; 01336 __x = std::ceil(__y); 01337 __w = __y * (2 - __y) * __param._M_1cx; 01338 if (__x > __param._M_d) 01339 continue; 01340 } 01341 else if (__u <= __c3) 01342 // NB: This case not in the book, nor in the Errata, 01343 // but should be ok... 01344 __x = -1; 01345 else if (__u <= __c4) 01346 __x = 0; 01347 else if (__u <= __c5) 01348 { 01349 __x = 1; 01350 // Only in the Errata, see libstdc++/83237. 01351 __w = __178; 01352 } 01353 else 01354 { 01355 const double __v = -std::log(1.0 - __aurng()); 01356 const double __y = __param._M_d 01357 + __v * __2cx / __param._M_d; 01358 __x = std::ceil(__y); 01359 __w = -__param._M_d * __param._M_1cx * (1 + __y / 2); 01360 } 01361 01362 __reject = (__w - __e - __x * __param._M_lm_thr 01363 > __param._M_lfm - std::lgamma(__x + __m + 1)); 01364 01365 __reject |= __x + __m >= __thr; 01366 01367 } while (__reject); 01368 01369 return result_type(__x + __m + __naf); 01370 } 01371 else 01372 #endif 01373 { 01374 _IntType __x = 0; 01375 double __prod = 1.0; 01376 01377 do 01378 { 01379 __prod *= __aurng(); 01380 __x += 1; 01381 } 01382 while (__prod > __param._M_lm_thr); 01383 01384 return __x - 1; 01385 } 01386 } 01387 01388 template<typename _IntType> 01389 template<typename _ForwardIterator, 01390 typename _UniformRandomNumberGenerator> 01391 void 01392 poisson_distribution<_IntType>:: 01393 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01394 _UniformRandomNumberGenerator& __urng, 01395 const param_type& __param) 01396 { 01397 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01398 // We could duplicate everything from operator()... 01399 while (__f != __t) 01400 *__f++ = this->operator()(__urng, __param); 01401 } 01402 01403 template<typename _IntType, 01404 typename _CharT, typename _Traits> 01405 std::basic_ostream<_CharT, _Traits>& 01406 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01407 const poisson_distribution<_IntType>& __x) 01408 { 01409 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01410 typedef typename __ostream_type::ios_base __ios_base; 01411 01412 const typename __ios_base::fmtflags __flags = __os.flags(); 01413 const _CharT __fill = __os.fill(); 01414 const std::streamsize __precision = __os.precision(); 01415 const _CharT __space = __os.widen(' '); 01416 __os.flags(__ios_base::scientific | __ios_base::left); 01417 __os.fill(__space); 01418 __os.precision(std::numeric_limits<double>::max_digits10); 01419 01420 __os << __x.mean() << __space << __x._M_nd; 01421 01422 __os.flags(__flags); 01423 __os.fill(__fill); 01424 __os.precision(__precision); 01425 return __os; 01426 } 01427 01428 template<typename _IntType, 01429 typename _CharT, typename _Traits> 01430 std::basic_istream<_CharT, _Traits>& 01431 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01432 poisson_distribution<_IntType>& __x) 01433 { 01434 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01435 typedef typename __istream_type::ios_base __ios_base; 01436 01437 const typename __ios_base::fmtflags __flags = __is.flags(); 01438 __is.flags(__ios_base::skipws); 01439 01440 double __mean; 01441 __is >> __mean >> __x._M_nd; 01442 __x.param(typename poisson_distribution<_IntType>::param_type(__mean)); 01443 01444 __is.flags(__flags); 01445 return __is; 01446 } 01447 01448 01449 template<typename _IntType> 01450 void 01451 binomial_distribution<_IntType>::param_type:: 01452 _M_initialize() 01453 { 01454 const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 01455 01456 _M_easy = true; 01457 01458 #if _GLIBCXX_USE_C99_MATH_TR1 01459 if (_M_t * __p12 >= 8) 01460 { 01461 _M_easy = false; 01462 const double __np = std::floor(_M_t * __p12); 01463 const double __pa = __np / _M_t; 01464 const double __1p = 1 - __pa; 01465 01466 const double __pi_4 = 0.7853981633974483096156608458198757L; 01467 const double __d1x = 01468 std::sqrt(__np * __1p * std::log(32 * __np 01469 / (81 * __pi_4 * __1p))); 01470 _M_d1 = std::round(std::max<double>(1.0, __d1x)); 01471 const double __d2x = 01472 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p 01473 / (__pi_4 * __pa))); 01474 _M_d2 = std::round(std::max<double>(1.0, __d2x)); 01475 01476 // sqrt(pi / 2) 01477 const double __spi_2 = 1.2533141373155002512078826424055226L; 01478 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); 01479 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); 01480 _M_c = 2 * _M_d1 / __np; 01481 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; 01482 const double __a12 = _M_a1 + _M_s2 * __spi_2; 01483 const double __s1s = _M_s1 * _M_s1; 01484 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) 01485 * 2 * __s1s / _M_d1 01486 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); 01487 const double __s2s = _M_s2 * _M_s2; 01488 _M_s = (_M_a123 + 2 * __s2s / _M_d2 01489 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); 01490 _M_lf = (std::lgamma(__np + 1) 01491 + std::lgamma(_M_t - __np + 1)); 01492 _M_lp1p = std::log(__pa / __1p); 01493 01494 _M_q = -std::log(1 - (__p12 - __pa) / __1p); 01495 } 01496 else 01497 #endif 01498 _M_q = -std::log(1 - __p12); 01499 } 01500 01501 template<typename _IntType> 01502 template<typename _UniformRandomNumberGenerator> 01503 typename binomial_distribution<_IntType>::result_type 01504 binomial_distribution<_IntType>:: 01505 _M_waiting(_UniformRandomNumberGenerator& __urng, 01506 _IntType __t, double __q) 01507 { 01508 _IntType __x = 0; 01509 double __sum = 0.0; 01510 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01511 __aurng(__urng); 01512 01513 do 01514 { 01515 if (__t == __x) 01516 return __x; 01517 const double __e = -std::log(1.0 - __aurng()); 01518 __sum += __e / (__t - __x); 01519 __x += 1; 01520 } 01521 while (__sum <= __q); 01522 01523 return __x - 1; 01524 } 01525 01526 /** 01527 * A rejection algorithm when t * p >= 8 and a simple waiting time 01528 * method - the second in the referenced book - otherwise. 01529 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 01530 * is defined. 01531 * 01532 * Reference: 01533 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 01534 * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 01535 */ 01536 template<typename _IntType> 01537 template<typename _UniformRandomNumberGenerator> 01538 typename binomial_distribution<_IntType>::result_type 01539 binomial_distribution<_IntType>:: 01540 operator()(_UniformRandomNumberGenerator& __urng, 01541 const param_type& __param) 01542 { 01543 result_type __ret; 01544 const _IntType __t = __param.t(); 01545 const double __p = __param.p(); 01546 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p; 01547 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01548 __aurng(__urng); 01549 01550 #if _GLIBCXX_USE_C99_MATH_TR1 01551 if (!__param._M_easy) 01552 { 01553 double __x; 01554 01555 // See comments above... 01556 const double __naf = 01557 (1 - std::numeric_limits<double>::epsilon()) / 2; 01558 const double __thr = 01559 std::numeric_limits<_IntType>::max() + __naf; 01560 01561 const double __np = std::floor(__t * __p12); 01562 01563 // sqrt(pi / 2) 01564 const double __spi_2 = 1.2533141373155002512078826424055226L; 01565 const double __a1 = __param._M_a1; 01566 const double __a12 = __a1 + __param._M_s2 * __spi_2; 01567 const double __a123 = __param._M_a123; 01568 const double __s1s = __param._M_s1 * __param._M_s1; 01569 const double __s2s = __param._M_s2 * __param._M_s2; 01570 01571 bool __reject; 01572 do 01573 { 01574 const double __u = __param._M_s * __aurng(); 01575 01576 double __v; 01577 01578 if (__u <= __a1) 01579 { 01580 const double __n = _M_nd(__urng); 01581 const double __y = __param._M_s1 * std::abs(__n); 01582 __reject = __y >= __param._M_d1; 01583 if (!__reject) 01584 { 01585 const double __e = -std::log(1.0 - __aurng()); 01586 __x = std::floor(__y); 01587 __v = -__e - __n * __n / 2 + __param._M_c; 01588 } 01589 } 01590 else if (__u <= __a12) 01591 { 01592 const double __n = _M_nd(__urng); 01593 const double __y = __param._M_s2 * std::abs(__n); 01594 __reject = __y >= __param._M_d2; 01595 if (!__reject) 01596 { 01597 const double __e = -std::log(1.0 - __aurng()); 01598 __x = std::floor(-__y); 01599 __v = -__e - __n * __n / 2; 01600 } 01601 } 01602 else if (__u <= __a123) 01603 { 01604 const double __e1 = -std::log(1.0 - __aurng()); 01605 const double __e2 = -std::log(1.0 - __aurng()); 01606 01607 const double __y = __param._M_d1 01608 + 2 * __s1s * __e1 / __param._M_d1; 01609 __x = std::floor(__y); 01610 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np) 01611 -__y / (2 * __s1s))); 01612 __reject = false; 01613 } 01614 else 01615 { 01616 const double __e1 = -std::log(1.0 - __aurng()); 01617 const double __e2 = -std::log(1.0 - __aurng()); 01618 01619 const double __y = __param._M_d2 01620 + 2 * __s2s * __e1 / __param._M_d2; 01621 __x = std::floor(-__y); 01622 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s); 01623 __reject = false; 01624 } 01625 01626 __reject = __reject || __x < -__np || __x > __t - __np; 01627 if (!__reject) 01628 { 01629 const double __lfx = 01630 std::lgamma(__np + __x + 1) 01631 + std::lgamma(__t - (__np + __x) + 1); 01632 __reject = __v > __param._M_lf - __lfx 01633 + __x * __param._M_lp1p; 01634 } 01635 01636 __reject |= __x + __np >= __thr; 01637 } 01638 while (__reject); 01639 01640 __x += __np + __naf; 01641 01642 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x), 01643 __param._M_q); 01644 __ret = _IntType(__x) + __z; 01645 } 01646 else 01647 #endif 01648 __ret = _M_waiting(__urng, __t, __param._M_q); 01649 01650 if (__p12 != __p) 01651 __ret = __t - __ret; 01652 return __ret; 01653 } 01654 01655 template<typename _IntType> 01656 template<typename _ForwardIterator, 01657 typename _UniformRandomNumberGenerator> 01658 void 01659 binomial_distribution<_IntType>:: 01660 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01661 _UniformRandomNumberGenerator& __urng, 01662 const param_type& __param) 01663 { 01664 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01665 // We could duplicate everything from operator()... 01666 while (__f != __t) 01667 *__f++ = this->operator()(__urng, __param); 01668 } 01669 01670 template<typename _IntType, 01671 typename _CharT, typename _Traits> 01672 std::basic_ostream<_CharT, _Traits>& 01673 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01674 const binomial_distribution<_IntType>& __x) 01675 { 01676 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01677 typedef typename __ostream_type::ios_base __ios_base; 01678 01679 const typename __ios_base::fmtflags __flags = __os.flags(); 01680 const _CharT __fill = __os.fill(); 01681 const std::streamsize __precision = __os.precision(); 01682 const _CharT __space = __os.widen(' '); 01683 __os.flags(__ios_base::scientific | __ios_base::left); 01684 __os.fill(__space); 01685 __os.precision(std::numeric_limits<double>::max_digits10); 01686 01687 __os << __x.t() << __space << __x.p() 01688 << __space << __x._M_nd; 01689 01690 __os.flags(__flags); 01691 __os.fill(__fill); 01692 __os.precision(__precision); 01693 return __os; 01694 } 01695 01696 template<typename _IntType, 01697 typename _CharT, typename _Traits> 01698 std::basic_istream<_CharT, _Traits>& 01699 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01700 binomial_distribution<_IntType>& __x) 01701 { 01702 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01703 typedef typename __istream_type::ios_base __ios_base; 01704 01705 const typename __ios_base::fmtflags __flags = __is.flags(); 01706 __is.flags(__ios_base::dec | __ios_base::skipws); 01707 01708 _IntType __t; 01709 double __p; 01710 __is >> __t >> __p >> __x._M_nd; 01711 __x.param(typename binomial_distribution<_IntType>:: 01712 param_type(__t, __p)); 01713 01714 __is.flags(__flags); 01715 return __is; 01716 } 01717 01718 01719 template<typename _RealType> 01720 template<typename _ForwardIterator, 01721 typename _UniformRandomNumberGenerator> 01722 void 01723 std::exponential_distribution<_RealType>:: 01724 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01725 _UniformRandomNumberGenerator& __urng, 01726 const param_type& __p) 01727 { 01728 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01729 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01730 __aurng(__urng); 01731 while (__f != __t) 01732 *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda(); 01733 } 01734 01735 template<typename _RealType, typename _CharT, typename _Traits> 01736 std::basic_ostream<_CharT, _Traits>& 01737 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01738 const exponential_distribution<_RealType>& __x) 01739 { 01740 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01741 typedef typename __ostream_type::ios_base __ios_base; 01742 01743 const typename __ios_base::fmtflags __flags = __os.flags(); 01744 const _CharT __fill = __os.fill(); 01745 const std::streamsize __precision = __os.precision(); 01746 __os.flags(__ios_base::scientific | __ios_base::left); 01747 __os.fill(__os.widen(' ')); 01748 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01749 01750 __os << __x.lambda(); 01751 01752 __os.flags(__flags); 01753 __os.fill(__fill); 01754 __os.precision(__precision); 01755 return __os; 01756 } 01757 01758 template<typename _RealType, typename _CharT, typename _Traits> 01759 std::basic_istream<_CharT, _Traits>& 01760 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01761 exponential_distribution<_RealType>& __x) 01762 { 01763 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01764 typedef typename __istream_type::ios_base __ios_base; 01765 01766 const typename __ios_base::fmtflags __flags = __is.flags(); 01767 __is.flags(__ios_base::dec | __ios_base::skipws); 01768 01769 _RealType __lambda; 01770 __is >> __lambda; 01771 __x.param(typename exponential_distribution<_RealType>:: 01772 param_type(__lambda)); 01773 01774 __is.flags(__flags); 01775 return __is; 01776 } 01777 01778 01779 /** 01780 * Polar method due to Marsaglia. 01781 * 01782 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 01783 * New York, 1986, Ch. V, Sect. 4.4. 01784 */ 01785 template<typename _RealType> 01786 template<typename _UniformRandomNumberGenerator> 01787 typename normal_distribution<_RealType>::result_type 01788 normal_distribution<_RealType>:: 01789 operator()(_UniformRandomNumberGenerator& __urng, 01790 const param_type& __param) 01791 { 01792 result_type __ret; 01793 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01794 __aurng(__urng); 01795 01796 if (_M_saved_available) 01797 { 01798 _M_saved_available = false; 01799 __ret = _M_saved; 01800 } 01801 else 01802 { 01803 result_type __x, __y, __r2; 01804 do 01805 { 01806 __x = result_type(2.0) * __aurng() - 1.0; 01807 __y = result_type(2.0) * __aurng() - 1.0; 01808 __r2 = __x * __x + __y * __y; 01809 } 01810 while (__r2 > 1.0 || __r2 == 0.0); 01811 01812 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01813 _M_saved = __x * __mult; 01814 _M_saved_available = true; 01815 __ret = __y * __mult; 01816 } 01817 01818 __ret = __ret * __param.stddev() + __param.mean(); 01819 return __ret; 01820 } 01821 01822 template<typename _RealType> 01823 template<typename _ForwardIterator, 01824 typename _UniformRandomNumberGenerator> 01825 void 01826 normal_distribution<_RealType>:: 01827 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01828 _UniformRandomNumberGenerator& __urng, 01829 const param_type& __param) 01830 { 01831 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01832 01833 if (__f == __t) 01834 return; 01835 01836 if (_M_saved_available) 01837 { 01838 _M_saved_available = false; 01839 *__f++ = _M_saved * __param.stddev() + __param.mean(); 01840 01841 if (__f == __t) 01842 return; 01843 } 01844 01845 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01846 __aurng(__urng); 01847 01848 while (__f + 1 < __t) 01849 { 01850 result_type __x, __y, __r2; 01851 do 01852 { 01853 __x = result_type(2.0) * __aurng() - 1.0; 01854 __y = result_type(2.0) * __aurng() - 1.0; 01855 __r2 = __x * __x + __y * __y; 01856 } 01857 while (__r2 > 1.0 || __r2 == 0.0); 01858 01859 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01860 *__f++ = __y * __mult * __param.stddev() + __param.mean(); 01861 *__f++ = __x * __mult * __param.stddev() + __param.mean(); 01862 } 01863 01864 if (__f != __t) 01865 { 01866 result_type __x, __y, __r2; 01867 do 01868 { 01869 __x = result_type(2.0) * __aurng() - 1.0; 01870 __y = result_type(2.0) * __aurng() - 1.0; 01871 __r2 = __x * __x + __y * __y; 01872 } 01873 while (__r2 > 1.0 || __r2 == 0.0); 01874 01875 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01876 _M_saved = __x * __mult; 01877 _M_saved_available = true; 01878 *__f = __y * __mult * __param.stddev() + __param.mean(); 01879 } 01880 } 01881 01882 template<typename _RealType> 01883 bool 01884 operator==(const std::normal_distribution<_RealType>& __d1, 01885 const std::normal_distribution<_RealType>& __d2) 01886 { 01887 if (__d1._M_param == __d2._M_param 01888 && __d1._M_saved_available == __d2._M_saved_available) 01889 { 01890 if (__d1._M_saved_available 01891 && __d1._M_saved == __d2._M_saved) 01892 return true; 01893 else if(!__d1._M_saved_available) 01894 return true; 01895 else 01896 return false; 01897 } 01898 else 01899 return false; 01900 } 01901 01902 template<typename _RealType, typename _CharT, typename _Traits> 01903 std::basic_ostream<_CharT, _Traits>& 01904 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01905 const normal_distribution<_RealType>& __x) 01906 { 01907 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01908 typedef typename __ostream_type::ios_base __ios_base; 01909 01910 const typename __ios_base::fmtflags __flags = __os.flags(); 01911 const _CharT __fill = __os.fill(); 01912 const std::streamsize __precision = __os.precision(); 01913 const _CharT __space = __os.widen(' '); 01914 __os.flags(__ios_base::scientific | __ios_base::left); 01915 __os.fill(__space); 01916 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01917 01918 __os << __x.mean() << __space << __x.stddev() 01919 << __space << __x._M_saved_available; 01920 if (__x._M_saved_available) 01921 __os << __space << __x._M_saved; 01922 01923 __os.flags(__flags); 01924 __os.fill(__fill); 01925 __os.precision(__precision); 01926 return __os; 01927 } 01928 01929 template<typename _RealType, typename _CharT, typename _Traits> 01930 std::basic_istream<_CharT, _Traits>& 01931 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01932 normal_distribution<_RealType>& __x) 01933 { 01934 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01935 typedef typename __istream_type::ios_base __ios_base; 01936 01937 const typename __ios_base::fmtflags __flags = __is.flags(); 01938 __is.flags(__ios_base::dec | __ios_base::skipws); 01939 01940 double __mean, __stddev; 01941 __is >> __mean >> __stddev 01942 >> __x._M_saved_available; 01943 if (__x._M_saved_available) 01944 __is >> __x._M_saved; 01945 __x.param(typename normal_distribution<_RealType>:: 01946 param_type(__mean, __stddev)); 01947 01948 __is.flags(__flags); 01949 return __is; 01950 } 01951 01952 01953 template<typename _RealType> 01954 template<typename _ForwardIterator, 01955 typename _UniformRandomNumberGenerator> 01956 void 01957 lognormal_distribution<_RealType>:: 01958 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01959 _UniformRandomNumberGenerator& __urng, 01960 const param_type& __p) 01961 { 01962 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01963 while (__f != __t) 01964 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m()); 01965 } 01966 01967 template<typename _RealType, typename _CharT, typename _Traits> 01968 std::basic_ostream<_CharT, _Traits>& 01969 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01970 const lognormal_distribution<_RealType>& __x) 01971 { 01972 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01973 typedef typename __ostream_type::ios_base __ios_base; 01974 01975 const typename __ios_base::fmtflags __flags = __os.flags(); 01976 const _CharT __fill = __os.fill(); 01977 const std::streamsize __precision = __os.precision(); 01978 const _CharT __space = __os.widen(' '); 01979 __os.flags(__ios_base::scientific | __ios_base::left); 01980 __os.fill(__space); 01981 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01982 01983 __os << __x.m() << __space << __x.s() 01984 << __space << __x._M_nd; 01985 01986 __os.flags(__flags); 01987 __os.fill(__fill); 01988 __os.precision(__precision); 01989 return __os; 01990 } 01991 01992 template<typename _RealType, typename _CharT, typename _Traits> 01993 std::basic_istream<_CharT, _Traits>& 01994 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01995 lognormal_distribution<_RealType>& __x) 01996 { 01997 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01998 typedef typename __istream_type::ios_base __ios_base; 01999 02000 const typename __ios_base::fmtflags __flags = __is.flags(); 02001 __is.flags(__ios_base::dec | __ios_base::skipws); 02002 02003 _RealType __m, __s; 02004 __is >> __m >> __s >> __x._M_nd; 02005 __x.param(typename lognormal_distribution<_RealType>:: 02006 param_type(__m, __s)); 02007 02008 __is.flags(__flags); 02009 return __is; 02010 } 02011 02012 template<typename _RealType> 02013 template<typename _ForwardIterator, 02014 typename _UniformRandomNumberGenerator> 02015 void 02016 std::chi_squared_distribution<_RealType>:: 02017 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02018 _UniformRandomNumberGenerator& __urng) 02019 { 02020 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02021 while (__f != __t) 02022 *__f++ = 2 * _M_gd(__urng); 02023 } 02024 02025 template<typename _RealType> 02026 template<typename _ForwardIterator, 02027 typename _UniformRandomNumberGenerator> 02028 void 02029 std::chi_squared_distribution<_RealType>:: 02030 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02031 _UniformRandomNumberGenerator& __urng, 02032 const typename 02033 std::gamma_distribution<result_type>::param_type& __p) 02034 { 02035 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02036 while (__f != __t) 02037 *__f++ = 2 * _M_gd(__urng, __p); 02038 } 02039 02040 template<typename _RealType, typename _CharT, typename _Traits> 02041 std::basic_ostream<_CharT, _Traits>& 02042 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02043 const chi_squared_distribution<_RealType>& __x) 02044 { 02045 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02046 typedef typename __ostream_type::ios_base __ios_base; 02047 02048 const typename __ios_base::fmtflags __flags = __os.flags(); 02049 const _CharT __fill = __os.fill(); 02050 const std::streamsize __precision = __os.precision(); 02051 const _CharT __space = __os.widen(' '); 02052 __os.flags(__ios_base::scientific | __ios_base::left); 02053 __os.fill(__space); 02054 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02055 02056 __os << __x.n() << __space << __x._M_gd; 02057 02058 __os.flags(__flags); 02059 __os.fill(__fill); 02060 __os.precision(__precision); 02061 return __os; 02062 } 02063 02064 template<typename _RealType, typename _CharT, typename _Traits> 02065 std::basic_istream<_CharT, _Traits>& 02066 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02067 chi_squared_distribution<_RealType>& __x) 02068 { 02069 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02070 typedef typename __istream_type::ios_base __ios_base; 02071 02072 const typename __ios_base::fmtflags __flags = __is.flags(); 02073 __is.flags(__ios_base::dec | __ios_base::skipws); 02074 02075 _RealType __n; 02076 __is >> __n >> __x._M_gd; 02077 __x.param(typename chi_squared_distribution<_RealType>:: 02078 param_type(__n)); 02079 02080 __is.flags(__flags); 02081 return __is; 02082 } 02083 02084 02085 template<typename _RealType> 02086 template<typename _UniformRandomNumberGenerator> 02087 typename cauchy_distribution<_RealType>::result_type 02088 cauchy_distribution<_RealType>:: 02089 operator()(_UniformRandomNumberGenerator& __urng, 02090 const param_type& __p) 02091 { 02092 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02093 __aurng(__urng); 02094 _RealType __u; 02095 do 02096 __u = __aurng(); 02097 while (__u == 0.5); 02098 02099 const _RealType __pi = 3.1415926535897932384626433832795029L; 02100 return __p.a() + __p.b() * std::tan(__pi * __u); 02101 } 02102 02103 template<typename _RealType> 02104 template<typename _ForwardIterator, 02105 typename _UniformRandomNumberGenerator> 02106 void 02107 cauchy_distribution<_RealType>:: 02108 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02109 _UniformRandomNumberGenerator& __urng, 02110 const param_type& __p) 02111 { 02112 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02113 const _RealType __pi = 3.1415926535897932384626433832795029L; 02114 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02115 __aurng(__urng); 02116 while (__f != __t) 02117 { 02118 _RealType __u; 02119 do 02120 __u = __aurng(); 02121 while (__u == 0.5); 02122 02123 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u); 02124 } 02125 } 02126 02127 template<typename _RealType, typename _CharT, typename _Traits> 02128 std::basic_ostream<_CharT, _Traits>& 02129 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02130 const cauchy_distribution<_RealType>& __x) 02131 { 02132 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02133 typedef typename __ostream_type::ios_base __ios_base; 02134 02135 const typename __ios_base::fmtflags __flags = __os.flags(); 02136 const _CharT __fill = __os.fill(); 02137 const std::streamsize __precision = __os.precision(); 02138 const _CharT __space = __os.widen(' '); 02139 __os.flags(__ios_base::scientific | __ios_base::left); 02140 __os.fill(__space); 02141 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02142 02143 __os << __x.a() << __space << __x.b(); 02144 02145 __os.flags(__flags); 02146 __os.fill(__fill); 02147 __os.precision(__precision); 02148 return __os; 02149 } 02150 02151 template<typename _RealType, typename _CharT, typename _Traits> 02152 std::basic_istream<_CharT, _Traits>& 02153 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02154 cauchy_distribution<_RealType>& __x) 02155 { 02156 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02157 typedef typename __istream_type::ios_base __ios_base; 02158 02159 const typename __ios_base::fmtflags __flags = __is.flags(); 02160 __is.flags(__ios_base::dec | __ios_base::skipws); 02161 02162 _RealType __a, __b; 02163 __is >> __a >> __b; 02164 __x.param(typename cauchy_distribution<_RealType>:: 02165 param_type(__a, __b)); 02166 02167 __is.flags(__flags); 02168 return __is; 02169 } 02170 02171 02172 template<typename _RealType> 02173 template<typename _ForwardIterator, 02174 typename _UniformRandomNumberGenerator> 02175 void 02176 std::fisher_f_distribution<_RealType>:: 02177 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02178 _UniformRandomNumberGenerator& __urng) 02179 { 02180 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02181 while (__f != __t) 02182 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m())); 02183 } 02184 02185 template<typename _RealType> 02186 template<typename _ForwardIterator, 02187 typename _UniformRandomNumberGenerator> 02188 void 02189 std::fisher_f_distribution<_RealType>:: 02190 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02191 _UniformRandomNumberGenerator& __urng, 02192 const param_type& __p) 02193 { 02194 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02195 typedef typename std::gamma_distribution<result_type>::param_type 02196 param_type; 02197 param_type __p1(__p.m() / 2); 02198 param_type __p2(__p.n() / 2); 02199 while (__f != __t) 02200 *__f++ = ((_M_gd_x(__urng, __p1) * n()) 02201 / (_M_gd_y(__urng, __p2) * m())); 02202 } 02203 02204 template<typename _RealType, typename _CharT, typename _Traits> 02205 std::basic_ostream<_CharT, _Traits>& 02206 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02207 const fisher_f_distribution<_RealType>& __x) 02208 { 02209 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02210 typedef typename __ostream_type::ios_base __ios_base; 02211 02212 const typename __ios_base::fmtflags __flags = __os.flags(); 02213 const _CharT __fill = __os.fill(); 02214 const std::streamsize __precision = __os.precision(); 02215 const _CharT __space = __os.widen(' '); 02216 __os.flags(__ios_base::scientific | __ios_base::left); 02217 __os.fill(__space); 02218 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02219 02220 __os << __x.m() << __space << __x.n() 02221 << __space << __x._M_gd_x << __space << __x._M_gd_y; 02222 02223 __os.flags(__flags); 02224 __os.fill(__fill); 02225 __os.precision(__precision); 02226 return __os; 02227 } 02228 02229 template<typename _RealType, typename _CharT, typename _Traits> 02230 std::basic_istream<_CharT, _Traits>& 02231 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02232 fisher_f_distribution<_RealType>& __x) 02233 { 02234 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02235 typedef typename __istream_type::ios_base __ios_base; 02236 02237 const typename __ios_base::fmtflags __flags = __is.flags(); 02238 __is.flags(__ios_base::dec | __ios_base::skipws); 02239 02240 _RealType __m, __n; 02241 __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y; 02242 __x.param(typename fisher_f_distribution<_RealType>:: 02243 param_type(__m, __n)); 02244 02245 __is.flags(__flags); 02246 return __is; 02247 } 02248 02249 02250 template<typename _RealType> 02251 template<typename _ForwardIterator, 02252 typename _UniformRandomNumberGenerator> 02253 void 02254 std::student_t_distribution<_RealType>:: 02255 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02256 _UniformRandomNumberGenerator& __urng) 02257 { 02258 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02259 while (__f != __t) 02260 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng)); 02261 } 02262 02263 template<typename _RealType> 02264 template<typename _ForwardIterator, 02265 typename _UniformRandomNumberGenerator> 02266 void 02267 std::student_t_distribution<_RealType>:: 02268 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02269 _UniformRandomNumberGenerator& __urng, 02270 const param_type& __p) 02271 { 02272 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02273 typename std::gamma_distribution<result_type>::param_type 02274 __p2(__p.n() / 2, 2); 02275 while (__f != __t) 02276 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2)); 02277 } 02278 02279 template<typename _RealType, typename _CharT, typename _Traits> 02280 std::basic_ostream<_CharT, _Traits>& 02281 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02282 const student_t_distribution<_RealType>& __x) 02283 { 02284 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02285 typedef typename __ostream_type::ios_base __ios_base; 02286 02287 const typename __ios_base::fmtflags __flags = __os.flags(); 02288 const _CharT __fill = __os.fill(); 02289 const std::streamsize __precision = __os.precision(); 02290 const _CharT __space = __os.widen(' '); 02291 __os.flags(__ios_base::scientific | __ios_base::left); 02292 __os.fill(__space); 02293 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02294 02295 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd; 02296 02297 __os.flags(__flags); 02298 __os.fill(__fill); 02299 __os.precision(__precision); 02300 return __os; 02301 } 02302 02303 template<typename _RealType, typename _CharT, typename _Traits> 02304 std::basic_istream<_CharT, _Traits>& 02305 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02306 student_t_distribution<_RealType>& __x) 02307 { 02308 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02309 typedef typename __istream_type::ios_base __ios_base; 02310 02311 const typename __ios_base::fmtflags __flags = __is.flags(); 02312 __is.flags(__ios_base::dec | __ios_base::skipws); 02313 02314 _RealType __n; 02315 __is >> __n >> __x._M_nd >> __x._M_gd; 02316 __x.param(typename student_t_distribution<_RealType>::param_type(__n)); 02317 02318 __is.flags(__flags); 02319 return __is; 02320 } 02321 02322 02323 template<typename _RealType> 02324 void 02325 gamma_distribution<_RealType>::param_type:: 02326 _M_initialize() 02327 { 02328 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha; 02329 02330 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0); 02331 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1); 02332 } 02333 02334 /** 02335 * Marsaglia, G. and Tsang, W. W. 02336 * "A Simple Method for Generating Gamma Variables" 02337 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000. 02338 */ 02339 template<typename _RealType> 02340 template<typename _UniformRandomNumberGenerator> 02341 typename gamma_distribution<_RealType>::result_type 02342 gamma_distribution<_RealType>:: 02343 operator()(_UniformRandomNumberGenerator& __urng, 02344 const param_type& __param) 02345 { 02346 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02347 __aurng(__urng); 02348 02349 result_type __u, __v, __n; 02350 const result_type __a1 = (__param._M_malpha 02351 - _RealType(1.0) / _RealType(3.0)); 02352 02353 do 02354 { 02355 do 02356 { 02357 __n = _M_nd(__urng); 02358 __v = result_type(1.0) + __param._M_a2 * __n; 02359 } 02360 while (__v <= 0.0); 02361 02362 __v = __v * __v * __v; 02363 __u = __aurng(); 02364 } 02365 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n 02366 && (std::log(__u) > (0.5 * __n * __n + __a1 02367 * (1.0 - __v + std::log(__v))))); 02368 02369 if (__param.alpha() == __param._M_malpha) 02370 return __a1 * __v * __param.beta(); 02371 else 02372 { 02373 do 02374 __u = __aurng(); 02375 while (__u == 0.0); 02376 02377 return (std::pow(__u, result_type(1.0) / __param.alpha()) 02378 * __a1 * __v * __param.beta()); 02379 } 02380 } 02381 02382 template<typename _RealType> 02383 template<typename _ForwardIterator, 02384 typename _UniformRandomNumberGenerator> 02385 void 02386 gamma_distribution<_RealType>:: 02387 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02388 _UniformRandomNumberGenerator& __urng, 02389 const param_type& __param) 02390 { 02391 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02392 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02393 __aurng(__urng); 02394 02395 result_type __u, __v, __n; 02396 const result_type __a1 = (__param._M_malpha 02397 - _RealType(1.0) / _RealType(3.0)); 02398 02399 if (__param.alpha() == __param._M_malpha) 02400 while (__f != __t) 02401 { 02402 do 02403 { 02404 do 02405 { 02406 __n = _M_nd(__urng); 02407 __v = result_type(1.0) + __param._M_a2 * __n; 02408 } 02409 while (__v <= 0.0); 02410 02411 __v = __v * __v * __v; 02412 __u = __aurng(); 02413 } 02414 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n 02415 && (std::log(__u) > (0.5 * __n * __n + __a1 02416 * (1.0 - __v + std::log(__v))))); 02417 02418 *__f++ = __a1 * __v * __param.beta(); 02419 } 02420 else 02421 while (__f != __t) 02422 { 02423 do 02424 { 02425 do 02426 { 02427 __n = _M_nd(__urng); 02428 __v = result_type(1.0) + __param._M_a2 * __n; 02429 } 02430 while (__v <= 0.0); 02431 02432 __v = __v * __v * __v; 02433 __u = __aurng(); 02434 } 02435 while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n 02436 && (std::log(__u) > (0.5 * __n * __n + __a1 02437 * (1.0 - __v + std::log(__v))))); 02438 02439 do 02440 __u = __aurng(); 02441 while (__u == 0.0); 02442 02443 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha()) 02444 * __a1 * __v * __param.beta()); 02445 } 02446 } 02447 02448 template<typename _RealType, typename _CharT, typename _Traits> 02449 std::basic_ostream<_CharT, _Traits>& 02450 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02451 const gamma_distribution<_RealType>& __x) 02452 { 02453 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02454 typedef typename __ostream_type::ios_base __ios_base; 02455 02456 const typename __ios_base::fmtflags __flags = __os.flags(); 02457 const _CharT __fill = __os.fill(); 02458 const std::streamsize __precision = __os.precision(); 02459 const _CharT __space = __os.widen(' '); 02460 __os.flags(__ios_base::scientific | __ios_base::left); 02461 __os.fill(__space); 02462 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02463 02464 __os << __x.alpha() << __space << __x.beta() 02465 << __space << __x._M_nd; 02466 02467 __os.flags(__flags); 02468 __os.fill(__fill); 02469 __os.precision(__precision); 02470 return __os; 02471 } 02472 02473 template<typename _RealType, typename _CharT, typename _Traits> 02474 std::basic_istream<_CharT, _Traits>& 02475 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02476 gamma_distribution<_RealType>& __x) 02477 { 02478 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02479 typedef typename __istream_type::ios_base __ios_base; 02480 02481 const typename __ios_base::fmtflags __flags = __is.flags(); 02482 __is.flags(__ios_base::dec | __ios_base::skipws); 02483 02484 _RealType __alpha_val, __beta_val; 02485 __is >> __alpha_val >> __beta_val >> __x._M_nd; 02486 __x.param(typename gamma_distribution<_RealType>:: 02487 param_type(__alpha_val, __beta_val)); 02488 02489 __is.flags(__flags); 02490 return __is; 02491 } 02492 02493 02494 template<typename _RealType> 02495 template<typename _UniformRandomNumberGenerator> 02496 typename weibull_distribution<_RealType>::result_type 02497 weibull_distribution<_RealType>:: 02498 operator()(_UniformRandomNumberGenerator& __urng, 02499 const param_type& __p) 02500 { 02501 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02502 __aurng(__urng); 02503 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 02504 result_type(1) / __p.a()); 02505 } 02506 02507 template<typename _RealType> 02508 template<typename _ForwardIterator, 02509 typename _UniformRandomNumberGenerator> 02510 void 02511 weibull_distribution<_RealType>:: 02512 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02513 _UniformRandomNumberGenerator& __urng, 02514 const param_type& __p) 02515 { 02516 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02517 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02518 __aurng(__urng); 02519 auto __inv_a = result_type(1) / __p.a(); 02520 02521 while (__f != __t) 02522 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 02523 __inv_a); 02524 } 02525 02526 template<typename _RealType, typename _CharT, typename _Traits> 02527 std::basic_ostream<_CharT, _Traits>& 02528 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02529 const weibull_distribution<_RealType>& __x) 02530 { 02531 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02532 typedef typename __ostream_type::ios_base __ios_base; 02533 02534 const typename __ios_base::fmtflags __flags = __os.flags(); 02535 const _CharT __fill = __os.fill(); 02536 const std::streamsize __precision = __os.precision(); 02537 const _CharT __space = __os.widen(' '); 02538 __os.flags(__ios_base::scientific | __ios_base::left); 02539 __os.fill(__space); 02540 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02541 02542 __os << __x.a() << __space << __x.b(); 02543 02544 __os.flags(__flags); 02545 __os.fill(__fill); 02546 __os.precision(__precision); 02547 return __os; 02548 } 02549 02550 template<typename _RealType, typename _CharT, typename _Traits> 02551 std::basic_istream<_CharT, _Traits>& 02552 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02553 weibull_distribution<_RealType>& __x) 02554 { 02555 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02556 typedef typename __istream_type::ios_base __ios_base; 02557 02558 const typename __ios_base::fmtflags __flags = __is.flags(); 02559 __is.flags(__ios_base::dec | __ios_base::skipws); 02560 02561 _RealType __a, __b; 02562 __is >> __a >> __b; 02563 __x.param(typename weibull_distribution<_RealType>:: 02564 param_type(__a, __b)); 02565 02566 __is.flags(__flags); 02567 return __is; 02568 } 02569 02570 02571 template<typename _RealType> 02572 template<typename _UniformRandomNumberGenerator> 02573 typename extreme_value_distribution<_RealType>::result_type 02574 extreme_value_distribution<_RealType>:: 02575 operator()(_UniformRandomNumberGenerator& __urng, 02576 const param_type& __p) 02577 { 02578 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02579 __aurng(__urng); 02580 return __p.a() - __p.b() * std::log(-std::log(result_type(1) 02581 - __aurng())); 02582 } 02583 02584 template<typename _RealType> 02585 template<typename _ForwardIterator, 02586 typename _UniformRandomNumberGenerator> 02587 void 02588 extreme_value_distribution<_RealType>:: 02589 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02590 _UniformRandomNumberGenerator& __urng, 02591 const param_type& __p) 02592 { 02593 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02594 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02595 __aurng(__urng); 02596 02597 while (__f != __t) 02598 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1) 02599 - __aurng())); 02600 } 02601 02602 template<typename _RealType, typename _CharT, typename _Traits> 02603 std::basic_ostream<_CharT, _Traits>& 02604 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02605 const extreme_value_distribution<_RealType>& __x) 02606 { 02607 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02608 typedef typename __ostream_type::ios_base __ios_base; 02609 02610 const typename __ios_base::fmtflags __flags = __os.flags(); 02611 const _CharT __fill = __os.fill(); 02612 const std::streamsize __precision = __os.precision(); 02613 const _CharT __space = __os.widen(' '); 02614 __os.flags(__ios_base::scientific | __ios_base::left); 02615 __os.fill(__space); 02616 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02617 02618 __os << __x.a() << __space << __x.b(); 02619 02620 __os.flags(__flags); 02621 __os.fill(__fill); 02622 __os.precision(__precision); 02623 return __os; 02624 } 02625 02626 template<typename _RealType, typename _CharT, typename _Traits> 02627 std::basic_istream<_CharT, _Traits>& 02628 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02629 extreme_value_distribution<_RealType>& __x) 02630 { 02631 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02632 typedef typename __istream_type::ios_base __ios_base; 02633 02634 const typename __ios_base::fmtflags __flags = __is.flags(); 02635 __is.flags(__ios_base::dec | __ios_base::skipws); 02636 02637 _RealType __a, __b; 02638 __is >> __a >> __b; 02639 __x.param(typename extreme_value_distribution<_RealType>:: 02640 param_type(__a, __b)); 02641 02642 __is.flags(__flags); 02643 return __is; 02644 } 02645 02646 02647 template<typename _IntType> 02648 void 02649 discrete_distribution<_IntType>::param_type:: 02650 _M_initialize() 02651 { 02652 if (_M_prob.size() < 2) 02653 { 02654 _M_prob.clear(); 02655 return; 02656 } 02657 02658 const double __sum = std::accumulate(_M_prob.begin(), 02659 _M_prob.end(), 0.0); 02660 // Now normalize the probabilites. 02661 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), 02662 __sum); 02663 // Accumulate partial sums. 02664 _M_cp.reserve(_M_prob.size()); 02665 std::partial_sum(_M_prob.begin(), _M_prob.end(), 02666 std::back_inserter(_M_cp)); 02667 // Make sure the last cumulative probability is one. 02668 _M_cp[_M_cp.size() - 1] = 1.0; 02669 } 02670 02671 template<typename _IntType> 02672 template<typename _Func> 02673 discrete_distribution<_IntType>::param_type:: 02674 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw) 02675 : _M_prob(), _M_cp() 02676 { 02677 const size_t __n = __nw == 0 ? 1 : __nw; 02678 const double __delta = (__xmax - __xmin) / __n; 02679 02680 _M_prob.reserve(__n); 02681 for (size_t __k = 0; __k < __nw; ++__k) 02682 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta)); 02683 02684 _M_initialize(); 02685 } 02686 02687 template<typename _IntType> 02688 template<typename _UniformRandomNumberGenerator> 02689 typename discrete_distribution<_IntType>::result_type 02690 discrete_distribution<_IntType>:: 02691 operator()(_UniformRandomNumberGenerator& __urng, 02692 const param_type& __param) 02693 { 02694 if (__param._M_cp.empty()) 02695 return result_type(0); 02696 02697 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02698 __aurng(__urng); 02699 02700 const double __p = __aurng(); 02701 auto __pos = std::lower_bound(__param._M_cp.begin(), 02702 __param._M_cp.end(), __p); 02703 02704 return __pos - __param._M_cp.begin(); 02705 } 02706 02707 template<typename _IntType> 02708 template<typename _ForwardIterator, 02709 typename _UniformRandomNumberGenerator> 02710 void 02711 discrete_distribution<_IntType>:: 02712 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02713 _UniformRandomNumberGenerator& __urng, 02714 const param_type& __param) 02715 { 02716 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02717 02718 if (__param._M_cp.empty()) 02719 { 02720 while (__f != __t) 02721 *__f++ = result_type(0); 02722 return; 02723 } 02724 02725 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02726 __aurng(__urng); 02727 02728 while (__f != __t) 02729 { 02730 const double __p = __aurng(); 02731 auto __pos = std::lower_bound(__param._M_cp.begin(), 02732 __param._M_cp.end(), __p); 02733 02734 *__f++ = __pos - __param._M_cp.begin(); 02735 } 02736 } 02737 02738 template<typename _IntType, typename _CharT, typename _Traits> 02739 std::basic_ostream<_CharT, _Traits>& 02740 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02741 const discrete_distribution<_IntType>& __x) 02742 { 02743 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02744 typedef typename __ostream_type::ios_base __ios_base; 02745 02746 const typename __ios_base::fmtflags __flags = __os.flags(); 02747 const _CharT __fill = __os.fill(); 02748 const std::streamsize __precision = __os.precision(); 02749 const _CharT __space = __os.widen(' '); 02750 __os.flags(__ios_base::scientific | __ios_base::left); 02751 __os.fill(__space); 02752 __os.precision(std::numeric_limits<double>::max_digits10); 02753 02754 std::vector<double> __prob = __x.probabilities(); 02755 __os << __prob.size(); 02756 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit) 02757 __os << __space << *__dit; 02758 02759 __os.flags(__flags); 02760 __os.fill(__fill); 02761 __os.precision(__precision); 02762 return __os; 02763 } 02764 02765 template<typename _IntType, typename _CharT, typename _Traits> 02766 std::basic_istream<_CharT, _Traits>& 02767 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02768 discrete_distribution<_IntType>& __x) 02769 { 02770 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02771 typedef typename __istream_type::ios_base __ios_base; 02772 02773 const typename __ios_base::fmtflags __flags = __is.flags(); 02774 __is.flags(__ios_base::dec | __ios_base::skipws); 02775 02776 size_t __n; 02777 __is >> __n; 02778 02779 std::vector<double> __prob_vec; 02780 __prob_vec.reserve(__n); 02781 for (; __n != 0; --__n) 02782 { 02783 double __prob; 02784 __is >> __prob; 02785 __prob_vec.push_back(__prob); 02786 } 02787 02788 __x.param(typename discrete_distribution<_IntType>:: 02789 param_type(__prob_vec.begin(), __prob_vec.end())); 02790 02791 __is.flags(__flags); 02792 return __is; 02793 } 02794 02795 02796 template<typename _RealType> 02797 void 02798 piecewise_constant_distribution<_RealType>::param_type:: 02799 _M_initialize() 02800 { 02801 if (_M_int.size() < 2 02802 || (_M_int.size() == 2 02803 && _M_int[0] == _RealType(0) 02804 && _M_int[1] == _RealType(1))) 02805 { 02806 _M_int.clear(); 02807 _M_den.clear(); 02808 return; 02809 } 02810 02811 const double __sum = std::accumulate(_M_den.begin(), 02812 _M_den.end(), 0.0); 02813 02814 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 02815 __sum); 02816 02817 _M_cp.reserve(_M_den.size()); 02818 std::partial_sum(_M_den.begin(), _M_den.end(), 02819 std::back_inserter(_M_cp)); 02820 02821 // Make sure the last cumulative probability is one. 02822 _M_cp[_M_cp.size() - 1] = 1.0; 02823 02824 for (size_t __k = 0; __k < _M_den.size(); ++__k) 02825 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k]; 02826 } 02827 02828 template<typename _RealType> 02829 template<typename _InputIteratorB, typename _InputIteratorW> 02830 piecewise_constant_distribution<_RealType>::param_type:: 02831 param_type(_InputIteratorB __bbegin, 02832 _InputIteratorB __bend, 02833 _InputIteratorW __wbegin) 02834 : _M_int(), _M_den(), _M_cp() 02835 { 02836 if (__bbegin != __bend) 02837 { 02838 for (;;) 02839 { 02840 _M_int.push_back(*__bbegin); 02841 ++__bbegin; 02842 if (__bbegin == __bend) 02843 break; 02844 02845 _M_den.push_back(*__wbegin); 02846 ++__wbegin; 02847 } 02848 } 02849 02850 _M_initialize(); 02851 } 02852 02853 template<typename _RealType> 02854 template<typename _Func> 02855 piecewise_constant_distribution<_RealType>::param_type:: 02856 param_type(initializer_list<_RealType> __bl, _Func __fw) 02857 : _M_int(), _M_den(), _M_cp() 02858 { 02859 _M_int.reserve(__bl.size()); 02860 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 02861 _M_int.push_back(*__biter); 02862 02863 _M_den.reserve(_M_int.size() - 1); 02864 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 02865 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k]))); 02866 02867 _M_initialize(); 02868 } 02869 02870 template<typename _RealType> 02871 template<typename _Func> 02872 piecewise_constant_distribution<_RealType>::param_type:: 02873 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 02874 : _M_int(), _M_den(), _M_cp() 02875 { 02876 const size_t __n = __nw == 0 ? 1 : __nw; 02877 const _RealType __delta = (__xmax - __xmin) / __n; 02878 02879 _M_int.reserve(__n + 1); 02880 for (size_t __k = 0; __k <= __nw; ++__k) 02881 _M_int.push_back(__xmin + __k * __delta); 02882 02883 _M_den.reserve(__n); 02884 for (size_t __k = 0; __k < __nw; ++__k) 02885 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta)); 02886 02887 _M_initialize(); 02888 } 02889 02890 template<typename _RealType> 02891 template<typename _UniformRandomNumberGenerator> 02892 typename piecewise_constant_distribution<_RealType>::result_type 02893 piecewise_constant_distribution<_RealType>:: 02894 operator()(_UniformRandomNumberGenerator& __urng, 02895 const param_type& __param) 02896 { 02897 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02898 __aurng(__urng); 02899 02900 const double __p = __aurng(); 02901 if (__param._M_cp.empty()) 02902 return __p; 02903 02904 auto __pos = std::lower_bound(__param._M_cp.begin(), 02905 __param._M_cp.end(), __p); 02906 const size_t __i = __pos - __param._M_cp.begin(); 02907 02908 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 02909 02910 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i]; 02911 } 02912 02913 template<typename _RealType> 02914 template<typename _ForwardIterator, 02915 typename _UniformRandomNumberGenerator> 02916 void 02917 piecewise_constant_distribution<_RealType>:: 02918 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02919 _UniformRandomNumberGenerator& __urng, 02920 const param_type& __param) 02921 { 02922 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02923 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02924 __aurng(__urng); 02925 02926 if (__param._M_cp.empty()) 02927 { 02928 while (__f != __t) 02929 *__f++ = __aurng(); 02930 return; 02931 } 02932 02933 while (__f != __t) 02934 { 02935 const double __p = __aurng(); 02936 02937 auto __pos = std::lower_bound(__param._M_cp.begin(), 02938 __param._M_cp.end(), __p); 02939 const size_t __i = __pos - __param._M_cp.begin(); 02940 02941 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 02942 02943 *__f++ = (__param._M_int[__i] 02944 + (__p - __pref) / __param._M_den[__i]); 02945 } 02946 } 02947 02948 template<typename _RealType, typename _CharT, typename _Traits> 02949 std::basic_ostream<_CharT, _Traits>& 02950 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02951 const piecewise_constant_distribution<_RealType>& __x) 02952 { 02953 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02954 typedef typename __ostream_type::ios_base __ios_base; 02955 02956 const typename __ios_base::fmtflags __flags = __os.flags(); 02957 const _CharT __fill = __os.fill(); 02958 const std::streamsize __precision = __os.precision(); 02959 const _CharT __space = __os.widen(' '); 02960 __os.flags(__ios_base::scientific | __ios_base::left); 02961 __os.fill(__space); 02962 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02963 02964 std::vector<_RealType> __int = __x.intervals(); 02965 __os << __int.size() - 1; 02966 02967 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 02968 __os << __space << *__xit; 02969 02970 std::vector<double> __den = __x.densities(); 02971 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 02972 __os << __space << *__dit; 02973 02974 __os.flags(__flags); 02975 __os.fill(__fill); 02976 __os.precision(__precision); 02977 return __os; 02978 } 02979 02980 template<typename _RealType, typename _CharT, typename _Traits> 02981 std::basic_istream<_CharT, _Traits>& 02982 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02983 piecewise_constant_distribution<_RealType>& __x) 02984 { 02985 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02986 typedef typename __istream_type::ios_base __ios_base; 02987 02988 const typename __ios_base::fmtflags __flags = __is.flags(); 02989 __is.flags(__ios_base::dec | __ios_base::skipws); 02990 02991 size_t __n; 02992 __is >> __n; 02993 02994 std::vector<_RealType> __int_vec; 02995 __int_vec.reserve(__n + 1); 02996 for (size_t __i = 0; __i <= __n; ++__i) 02997 { 02998 _RealType __int; 02999 __is >> __int; 03000 __int_vec.push_back(__int); 03001 } 03002 03003 std::vector<double> __den_vec; 03004 __den_vec.reserve(__n); 03005 for (size_t __i = 0; __i < __n; ++__i) 03006 { 03007 double __den; 03008 __is >> __den; 03009 __den_vec.push_back(__den); 03010 } 03011 03012 __x.param(typename piecewise_constant_distribution<_RealType>:: 03013 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 03014 03015 __is.flags(__flags); 03016 return __is; 03017 } 03018 03019 03020 template<typename _RealType> 03021 void 03022 piecewise_linear_distribution<_RealType>::param_type:: 03023 _M_initialize() 03024 { 03025 if (_M_int.size() < 2 03026 || (_M_int.size() == 2 03027 && _M_int[0] == _RealType(0) 03028 && _M_int[1] == _RealType(1) 03029 && _M_den[0] == _M_den[1])) 03030 { 03031 _M_int.clear(); 03032 _M_den.clear(); 03033 return; 03034 } 03035 03036 double __sum = 0.0; 03037 _M_cp.reserve(_M_int.size() - 1); 03038 _M_m.reserve(_M_int.size() - 1); 03039 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 03040 { 03041 const _RealType __delta = _M_int[__k + 1] - _M_int[__k]; 03042 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta; 03043 _M_cp.push_back(__sum); 03044 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta); 03045 } 03046 03047 // Now normalize the densities... 03048 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 03049 __sum); 03050 // ... and partial sums... 03051 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum); 03052 // ... and slopes. 03053 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum); 03054 03055 // Make sure the last cumulative probablility is one. 03056 _M_cp[_M_cp.size() - 1] = 1.0; 03057 } 03058 03059 template<typename _RealType> 03060 template<typename _InputIteratorB, typename _InputIteratorW> 03061 piecewise_linear_distribution<_RealType>::param_type:: 03062 param_type(_InputIteratorB __bbegin, 03063 _InputIteratorB __bend, 03064 _InputIteratorW __wbegin) 03065 : _M_int(), _M_den(), _M_cp(), _M_m() 03066 { 03067 for (; __bbegin != __bend; ++__bbegin, ++__wbegin) 03068 { 03069 _M_int.push_back(*__bbegin); 03070 _M_den.push_back(*__wbegin); 03071 } 03072 03073 _M_initialize(); 03074 } 03075 03076 template<typename _RealType> 03077 template<typename _Func> 03078 piecewise_linear_distribution<_RealType>::param_type:: 03079 param_type(initializer_list<_RealType> __bl, _Func __fw) 03080 : _M_int(), _M_den(), _M_cp(), _M_m() 03081 { 03082 _M_int.reserve(__bl.size()); 03083 _M_den.reserve(__bl.size()); 03084 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 03085 { 03086 _M_int.push_back(*__biter); 03087 _M_den.push_back(__fw(*__biter)); 03088 } 03089 03090 _M_initialize(); 03091 } 03092 03093 template<typename _RealType> 03094 template<typename _Func> 03095 piecewise_linear_distribution<_RealType>::param_type:: 03096 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 03097 : _M_int(), _M_den(), _M_cp(), _M_m() 03098 { 03099 const size_t __n = __nw == 0 ? 1 : __nw; 03100 const _RealType __delta = (__xmax - __xmin) / __n; 03101 03102 _M_int.reserve(__n + 1); 03103 _M_den.reserve(__n + 1); 03104 for (size_t __k = 0; __k <= __nw; ++__k) 03105 { 03106 _M_int.push_back(__xmin + __k * __delta); 03107 _M_den.push_back(__fw(_M_int[__k] + __delta)); 03108 } 03109 03110 _M_initialize(); 03111 } 03112 03113 template<typename _RealType> 03114 template<typename _UniformRandomNumberGenerator> 03115 typename piecewise_linear_distribution<_RealType>::result_type 03116 piecewise_linear_distribution<_RealType>:: 03117 operator()(_UniformRandomNumberGenerator& __urng, 03118 const param_type& __param) 03119 { 03120 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 03121 __aurng(__urng); 03122 03123 const double __p = __aurng(); 03124 if (__param._M_cp.empty()) 03125 return __p; 03126 03127 auto __pos = std::lower_bound(__param._M_cp.begin(), 03128 __param._M_cp.end(), __p); 03129 const size_t __i = __pos - __param._M_cp.begin(); 03130 03131 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 03132 03133 const double __a = 0.5 * __param._M_m[__i]; 03134 const double __b = __param._M_den[__i]; 03135 const double __cm = __p - __pref; 03136 03137 _RealType __x = __param._M_int[__i]; 03138 if (__a == 0) 03139 __x += __cm / __b; 03140 else 03141 { 03142 const double __d = __b * __b + 4.0 * __a * __cm; 03143 __x += 0.5 * (std::sqrt(__d) - __b) / __a; 03144 } 03145 03146 return __x; 03147 } 03148 03149 template<typename _RealType> 03150 template<typename _ForwardIterator, 03151 typename _UniformRandomNumberGenerator> 03152 void 03153 piecewise_linear_distribution<_RealType>:: 03154 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 03155 _UniformRandomNumberGenerator& __urng, 03156 const param_type& __param) 03157 { 03158 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 03159 // We could duplicate everything from operator()... 03160 while (__f != __t) 03161 *__f++ = this->operator()(__urng, __param); 03162 } 03163 03164 template<typename _RealType, typename _CharT, typename _Traits> 03165 std::basic_ostream<_CharT, _Traits>& 03166 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 03167 const piecewise_linear_distribution<_RealType>& __x) 03168 { 03169 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 03170 typedef typename __ostream_type::ios_base __ios_base; 03171 03172 const typename __ios_base::fmtflags __flags = __os.flags(); 03173 const _CharT __fill = __os.fill(); 03174 const std::streamsize __precision = __os.precision(); 03175 const _CharT __space = __os.widen(' '); 03176 __os.flags(__ios_base::scientific | __ios_base::left); 03177 __os.fill(__space); 03178 __os.precision(std::numeric_limits<_RealType>::max_digits10); 03179 03180 std::vector<_RealType> __int = __x.intervals(); 03181 __os << __int.size() - 1; 03182 03183 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 03184 __os << __space << *__xit; 03185 03186 std::vector<double> __den = __x.densities(); 03187 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 03188 __os << __space << *__dit; 03189 03190 __os.flags(__flags); 03191 __os.fill(__fill); 03192 __os.precision(__precision); 03193 return __os; 03194 } 03195 03196 template<typename _RealType, typename _CharT, typename _Traits> 03197 std::basic_istream<_CharT, _Traits>& 03198 operator>>(std::basic_istream<_CharT, _Traits>& __is, 03199 piecewise_linear_distribution<_RealType>& __x) 03200 { 03201 typedef std::basic_istream<_CharT, _Traits> __istream_type; 03202 typedef typename __istream_type::ios_base __ios_base; 03203 03204 const typename __ios_base::fmtflags __flags = __is.flags(); 03205 __is.flags(__ios_base::dec | __ios_base::skipws); 03206 03207 size_t __n; 03208 __is >> __n; 03209 03210 std::vector<_RealType> __int_vec; 03211 __int_vec.reserve(__n + 1); 03212 for (size_t __i = 0; __i <= __n; ++__i) 03213 { 03214 _RealType __int; 03215 __is >> __int; 03216 __int_vec.push_back(__int); 03217 } 03218 03219 std::vector<double> __den_vec; 03220 __den_vec.reserve(__n + 1); 03221 for (size_t __i = 0; __i <= __n; ++__i) 03222 { 03223 double __den; 03224 __is >> __den; 03225 __den_vec.push_back(__den); 03226 } 03227 03228 __x.param(typename piecewise_linear_distribution<_RealType>:: 03229 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 03230 03231 __is.flags(__flags); 03232 return __is; 03233 } 03234 03235 03236 template<typename _IntType> 03237 seed_seq::seed_seq(std::initializer_list<_IntType> __il) 03238 { 03239 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter) 03240 _M_v.push_back(__detail::__mod<result_type, 03241 __detail::_Shift<result_type, 32>::__value>(*__iter)); 03242 } 03243 03244 template<typename _InputIterator> 03245 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end) 03246 { 03247 for (_InputIterator __iter = __begin; __iter != __end; ++__iter) 03248 _M_v.push_back(__detail::__mod<result_type, 03249 __detail::_Shift<result_type, 32>::__value>(*__iter)); 03250 } 03251 03252 template<typename _RandomAccessIterator> 03253 void 03254 seed_seq::generate(_RandomAccessIterator __begin, 03255 _RandomAccessIterator __end) 03256 { 03257 typedef typename iterator_traits<_RandomAccessIterator>::value_type 03258 _Type; 03259 03260 if (__begin == __end) 03261 return; 03262 03263 std::fill(__begin, __end, _Type(0x8b8b8b8bu)); 03264 03265 const size_t __n = __end - __begin; 03266 const size_t __s = _M_v.size(); 03267 const size_t __t = (__n >= 623) ? 11 03268 : (__n >= 68) ? 7 03269 : (__n >= 39) ? 5 03270 : (__n >= 7) ? 3 03271 : (__n - 1) / 2; 03272 const size_t __p = (__n - __t) / 2; 03273 const size_t __q = __p + __t; 03274 const size_t __m = std::max(size_t(__s + 1), __n); 03275 03276 for (size_t __k = 0; __k < __m; ++__k) 03277 { 03278 _Type __arg = (__begin[__k % __n] 03279 ^ __begin[(__k + __p) % __n] 03280 ^ __begin[(__k - 1) % __n]); 03281 _Type __r1 = __arg ^ (__arg >> 27); 03282 __r1 = __detail::__mod<_Type, 03283 __detail::_Shift<_Type, 32>::__value>(1664525u * __r1); 03284 _Type __r2 = __r1; 03285 if (__k == 0) 03286 __r2 += __s; 03287 else if (__k <= __s) 03288 __r2 += __k % __n + _M_v[__k - 1]; 03289 else 03290 __r2 += __k % __n; 03291 __r2 = __detail::__mod<_Type, 03292 __detail::_Shift<_Type, 32>::__value>(__r2); 03293 __begin[(__k + __p) % __n] += __r1; 03294 __begin[(__k + __q) % __n] += __r2; 03295 __begin[__k % __n] = __r2; 03296 } 03297 03298 for (size_t __k = __m; __k < __m + __n; ++__k) 03299 { 03300 _Type __arg = (__begin[__k % __n] 03301 + __begin[(__k + __p) % __n] 03302 + __begin[(__k - 1) % __n]); 03303 _Type __r3 = __arg ^ (__arg >> 27); 03304 __r3 = __detail::__mod<_Type, 03305 __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3); 03306 _Type __r4 = __r3 - __k % __n; 03307 __r4 = __detail::__mod<_Type, 03308 __detail::_Shift<_Type, 32>::__value>(__r4); 03309 __begin[(__k + __p) % __n] ^= __r3; 03310 __begin[(__k + __q) % __n] ^= __r4; 03311 __begin[__k % __n] = __r4; 03312 } 03313 } 03314 03315 template<typename _RealType, size_t __bits, 03316 typename _UniformRandomNumberGenerator> 03317 _RealType 03318 generate_canonical(_UniformRandomNumberGenerator& __urng) 03319 { 03320 static_assert(std::is_floating_point<_RealType>::value, 03321 "template argument must be a floating point type"); 03322 03323 const size_t __b 03324 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits), 03325 __bits); 03326 const long double __r = static_cast<long double>(__urng.max()) 03327 - static_cast<long double>(__urng.min()) + 1.0L; 03328 const size_t __log2r = std::log(__r) / std::log(2.0L); 03329 const size_t __m = std::max<size_t>(1UL, 03330 (__b + __log2r - 1UL) / __log2r); 03331 _RealType __ret; 03332 _RealType __sum = _RealType(0); 03333 _RealType __tmp = _RealType(1); 03334 for (size_t __k = __m; __k != 0; --__k) 03335 { 03336 __sum += _RealType(__urng() - __urng.min()) * __tmp; 03337 __tmp *= __r; 03338 } 03339 __ret = __sum / __tmp; 03340 if (__builtin_expect(__ret >= _RealType(1), 0)) 03341 { 03342 #if _GLIBCXX_USE_C99_MATH_TR1 03343 __ret = std::nextafter(_RealType(1), _RealType(0)); 03344 #else 03345 __ret = _RealType(1) 03346 - std::numeric_limits<_RealType>::epsilon() / _RealType(2); 03347 #endif 03348 } 03349 return __ret; 03350 } 03351 03352 _GLIBCXX_END_NAMESPACE_VERSION 03353 } // namespace 03354 03355 #endif