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