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, 01652 _IntType __t, double __q) 01653 { 01654 _IntType __x = 0; 01655 double __sum = 0.0; 01656 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01657 __aurng(__urng); 01658 01659 do 01660 { 01661 if (__t == __x) 01662 return __x; 01663 const double __e = -std::log(1.0 - __aurng()); 01664 __sum += __e / (__t - __x); 01665 __x += 1; 01666 } 01667 while (__sum <= __q); 01668 01669 return __x - 1; 01670 } 01671 01672 /** 01673 * A rejection algorithm when t * p >= 8 and a simple waiting time 01674 * method - the second in the referenced book - otherwise. 01675 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 01676 * is defined. 01677 * 01678 * Reference: 01679 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 01680 * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 01681 */ 01682 template<typename _IntType> 01683 template<typename _UniformRandomNumberGenerator> 01684 typename binomial_distribution<_IntType>::result_type 01685 binomial_distribution<_IntType>:: 01686 operator()(_UniformRandomNumberGenerator& __urng, 01687 const param_type& __param) 01688 { 01689 result_type __ret; 01690 const _IntType __t = __param.t(); 01691 const double __p = __param.p(); 01692 const double __p12 = __p <= 0.5 ? __p : 1.0 - __p; 01693 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 01694 __aurng(__urng); 01695 01696 #if _GLIBCXX_USE_C99_MATH_TR1 01697 if (!__param._M_easy) 01698 { 01699 double __x; 01700 01701 // See comments above... 01702 const double __naf = 01703 (1 - std::numeric_limits<double>::epsilon()) / 2; 01704 const double __thr = 01705 std::numeric_limits<_IntType>::max() + __naf; 01706 01707 const double __np = std::floor(__t * __p12); 01708 01709 // sqrt(pi / 2) 01710 const double __spi_2 = 1.2533141373155002512078826424055226L; 01711 const double __a1 = __param._M_a1; 01712 const double __a12 = __a1 + __param._M_s2 * __spi_2; 01713 const double __a123 = __param._M_a123; 01714 const double __s1s = __param._M_s1 * __param._M_s1; 01715 const double __s2s = __param._M_s2 * __param._M_s2; 01716 01717 bool __reject; 01718 do 01719 { 01720 const double __u = __param._M_s * __aurng(); 01721 01722 double __v; 01723 01724 if (__u <= __a1) 01725 { 01726 const double __n = _M_nd(__urng); 01727 const double __y = __param._M_s1 * std::abs(__n); 01728 __reject = __y >= __param._M_d1; 01729 if (!__reject) 01730 { 01731 const double __e = -std::log(1.0 - __aurng()); 01732 __x = std::floor(__y); 01733 __v = -__e - __n * __n / 2 + __param._M_c; 01734 } 01735 } 01736 else if (__u <= __a12) 01737 { 01738 const double __n = _M_nd(__urng); 01739 const double __y = __param._M_s2 * std::abs(__n); 01740 __reject = __y >= __param._M_d2; 01741 if (!__reject) 01742 { 01743 const double __e = -std::log(1.0 - __aurng()); 01744 __x = std::floor(-__y); 01745 __v = -__e - __n * __n / 2; 01746 } 01747 } 01748 else if (__u <= __a123) 01749 { 01750 const double __e1 = -std::log(1.0 - __aurng()); 01751 const double __e2 = -std::log(1.0 - __aurng()); 01752 01753 const double __y = __param._M_d1 01754 + 2 * __s1s * __e1 / __param._M_d1; 01755 __x = std::floor(__y); 01756 __v = (-__e2 + __param._M_d1 * (1 / (__t - __np) 01757 -__y / (2 * __s1s))); 01758 __reject = false; 01759 } 01760 else 01761 { 01762 const double __e1 = -std::log(1.0 - __aurng()); 01763 const double __e2 = -std::log(1.0 - __aurng()); 01764 01765 const double __y = __param._M_d2 01766 + 2 * __s2s * __e1 / __param._M_d2; 01767 __x = std::floor(-__y); 01768 __v = -__e2 - __param._M_d2 * __y / (2 * __s2s); 01769 __reject = false; 01770 } 01771 01772 __reject = __reject || __x < -__np || __x > __t - __np; 01773 if (!__reject) 01774 { 01775 const double __lfx = 01776 std::lgamma(__np + __x + 1) 01777 + std::lgamma(__t - (__np + __x) + 1); 01778 __reject = __v > __param._M_lf - __lfx 01779 + __x * __param._M_lp1p; 01780 } 01781 01782 __reject |= __x + __np >= __thr; 01783 } 01784 while (__reject); 01785 01786 __x += __np + __naf; 01787 01788 const _IntType __z = _M_waiting(__urng, __t - _IntType(__x), 01789 __param._M_q); 01790 __ret = _IntType(__x) + __z; 01791 } 01792 else 01793 #endif 01794 __ret = _M_waiting(__urng, __t, __param._M_q); 01795 01796 if (__p12 != __p) 01797 __ret = __t - __ret; 01798 return __ret; 01799 } 01800 01801 template<typename _IntType> 01802 template<typename _ForwardIterator, 01803 typename _UniformRandomNumberGenerator> 01804 void 01805 binomial_distribution<_IntType>:: 01806 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01807 _UniformRandomNumberGenerator& __urng, 01808 const param_type& __param) 01809 { 01810 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01811 // We could duplicate everything from operator()... 01812 while (__f != __t) 01813 *__f++ = this->operator()(__urng, __param); 01814 } 01815 01816 template<typename _IntType, 01817 typename _CharT, typename _Traits> 01818 std::basic_ostream<_CharT, _Traits>& 01819 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01820 const binomial_distribution<_IntType>& __x) 01821 { 01822 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01823 typedef typename __ostream_type::ios_base __ios_base; 01824 01825 const typename __ios_base::fmtflags __flags = __os.flags(); 01826 const _CharT __fill = __os.fill(); 01827 const std::streamsize __precision = __os.precision(); 01828 const _CharT __space = __os.widen(' '); 01829 __os.flags(__ios_base::scientific | __ios_base::left); 01830 __os.fill(__space); 01831 __os.precision(std::numeric_limits<double>::max_digits10); 01832 01833 __os << __x.t() << __space << __x.p() 01834 << __space << __x._M_nd; 01835 01836 __os.flags(__flags); 01837 __os.fill(__fill); 01838 __os.precision(__precision); 01839 return __os; 01840 } 01841 01842 template<typename _IntType, 01843 typename _CharT, typename _Traits> 01844 std::basic_istream<_CharT, _Traits>& 01845 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01846 binomial_distribution<_IntType>& __x) 01847 { 01848 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01849 typedef typename __istream_type::ios_base __ios_base; 01850 01851 const typename __ios_base::fmtflags __flags = __is.flags(); 01852 __is.flags(__ios_base::dec | __ios_base::skipws); 01853 01854 _IntType __t; 01855 double __p; 01856 __is >> __t >> __p >> __x._M_nd; 01857 __x.param(typename binomial_distribution<_IntType>:: 01858 param_type(__t, __p)); 01859 01860 __is.flags(__flags); 01861 return __is; 01862 } 01863 01864 01865 template<typename _RealType> 01866 template<typename _ForwardIterator, 01867 typename _UniformRandomNumberGenerator> 01868 void 01869 std::exponential_distribution<_RealType>:: 01870 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01871 _UniformRandomNumberGenerator& __urng, 01872 const param_type& __p) 01873 { 01874 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01875 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01876 __aurng(__urng); 01877 while (__f != __t) 01878 *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda(); 01879 } 01880 01881 template<typename _RealType, typename _CharT, typename _Traits> 01882 std::basic_ostream<_CharT, _Traits>& 01883 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 01884 const exponential_distribution<_RealType>& __x) 01885 { 01886 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 01887 typedef typename __ostream_type::ios_base __ios_base; 01888 01889 const typename __ios_base::fmtflags __flags = __os.flags(); 01890 const _CharT __fill = __os.fill(); 01891 const std::streamsize __precision = __os.precision(); 01892 __os.flags(__ios_base::scientific | __ios_base::left); 01893 __os.fill(__os.widen(' ')); 01894 __os.precision(std::numeric_limits<_RealType>::max_digits10); 01895 01896 __os << __x.lambda(); 01897 01898 __os.flags(__flags); 01899 __os.fill(__fill); 01900 __os.precision(__precision); 01901 return __os; 01902 } 01903 01904 template<typename _RealType, typename _CharT, typename _Traits> 01905 std::basic_istream<_CharT, _Traits>& 01906 operator>>(std::basic_istream<_CharT, _Traits>& __is, 01907 exponential_distribution<_RealType>& __x) 01908 { 01909 typedef std::basic_istream<_CharT, _Traits> __istream_type; 01910 typedef typename __istream_type::ios_base __ios_base; 01911 01912 const typename __ios_base::fmtflags __flags = __is.flags(); 01913 __is.flags(__ios_base::dec | __ios_base::skipws); 01914 01915 _RealType __lambda; 01916 __is >> __lambda; 01917 __x.param(typename exponential_distribution<_RealType>:: 01918 param_type(__lambda)); 01919 01920 __is.flags(__flags); 01921 return __is; 01922 } 01923 01924 01925 /** 01926 * Polar method due to Marsaglia. 01927 * 01928 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 01929 * New York, 1986, Ch. V, Sect. 4.4. 01930 */ 01931 template<typename _RealType> 01932 template<typename _UniformRandomNumberGenerator> 01933 typename normal_distribution<_RealType>::result_type 01934 normal_distribution<_RealType>:: 01935 operator()(_UniformRandomNumberGenerator& __urng, 01936 const param_type& __param) 01937 { 01938 result_type __ret; 01939 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01940 __aurng(__urng); 01941 01942 if (_M_saved_available) 01943 { 01944 _M_saved_available = false; 01945 __ret = _M_saved; 01946 } 01947 else 01948 { 01949 result_type __x, __y, __r2; 01950 do 01951 { 01952 __x = result_type(2.0) * __aurng() - 1.0; 01953 __y = result_type(2.0) * __aurng() - 1.0; 01954 __r2 = __x * __x + __y * __y; 01955 } 01956 while (__r2 > 1.0 || __r2 == 0.0); 01957 01958 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 01959 _M_saved = __x * __mult; 01960 _M_saved_available = true; 01961 __ret = __y * __mult; 01962 } 01963 01964 __ret = __ret * __param.stddev() + __param.mean(); 01965 return __ret; 01966 } 01967 01968 template<typename _RealType> 01969 template<typename _ForwardIterator, 01970 typename _UniformRandomNumberGenerator> 01971 void 01972 normal_distribution<_RealType>:: 01973 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 01974 _UniformRandomNumberGenerator& __urng, 01975 const param_type& __param) 01976 { 01977 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 01978 01979 if (__f == __t) 01980 return; 01981 01982 if (_M_saved_available) 01983 { 01984 _M_saved_available = false; 01985 *__f++ = _M_saved * __param.stddev() + __param.mean(); 01986 01987 if (__f == __t) 01988 return; 01989 } 01990 01991 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 01992 __aurng(__urng); 01993 01994 while (__f + 1 < __t) 01995 { 01996 result_type __x, __y, __r2; 01997 do 01998 { 01999 __x = result_type(2.0) * __aurng() - 1.0; 02000 __y = result_type(2.0) * __aurng() - 1.0; 02001 __r2 = __x * __x + __y * __y; 02002 } 02003 while (__r2 > 1.0 || __r2 == 0.0); 02004 02005 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 02006 *__f++ = __y * __mult * __param.stddev() + __param.mean(); 02007 *__f++ = __x * __mult * __param.stddev() + __param.mean(); 02008 } 02009 02010 if (__f != __t) 02011 { 02012 result_type __x, __y, __r2; 02013 do 02014 { 02015 __x = result_type(2.0) * __aurng() - 1.0; 02016 __y = result_type(2.0) * __aurng() - 1.0; 02017 __r2 = __x * __x + __y * __y; 02018 } 02019 while (__r2 > 1.0 || __r2 == 0.0); 02020 02021 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 02022 _M_saved = __x * __mult; 02023 _M_saved_available = true; 02024 *__f = __y * __mult * __param.stddev() + __param.mean(); 02025 } 02026 } 02027 02028 template<typename _RealType> 02029 bool 02030 operator==(const std::normal_distribution<_RealType>& __d1, 02031 const std::normal_distribution<_RealType>& __d2) 02032 { 02033 if (__d1._M_param == __d2._M_param 02034 && __d1._M_saved_available == __d2._M_saved_available) 02035 { 02036 if (__d1._M_saved_available 02037 && __d1._M_saved == __d2._M_saved) 02038 return true; 02039 else if(!__d1._M_saved_available) 02040 return true; 02041 else 02042 return false; 02043 } 02044 else 02045 return false; 02046 } 02047 02048 template<typename _RealType, typename _CharT, typename _Traits> 02049 std::basic_ostream<_CharT, _Traits>& 02050 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02051 const normal_distribution<_RealType>& __x) 02052 { 02053 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02054 typedef typename __ostream_type::ios_base __ios_base; 02055 02056 const typename __ios_base::fmtflags __flags = __os.flags(); 02057 const _CharT __fill = __os.fill(); 02058 const std::streamsize __precision = __os.precision(); 02059 const _CharT __space = __os.widen(' '); 02060 __os.flags(__ios_base::scientific | __ios_base::left); 02061 __os.fill(__space); 02062 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02063 02064 __os << __x.mean() << __space << __x.stddev() 02065 << __space << __x._M_saved_available; 02066 if (__x._M_saved_available) 02067 __os << __space << __x._M_saved; 02068 02069 __os.flags(__flags); 02070 __os.fill(__fill); 02071 __os.precision(__precision); 02072 return __os; 02073 } 02074 02075 template<typename _RealType, typename _CharT, typename _Traits> 02076 std::basic_istream<_CharT, _Traits>& 02077 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02078 normal_distribution<_RealType>& __x) 02079 { 02080 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02081 typedef typename __istream_type::ios_base __ios_base; 02082 02083 const typename __ios_base::fmtflags __flags = __is.flags(); 02084 __is.flags(__ios_base::dec | __ios_base::skipws); 02085 02086 double __mean, __stddev; 02087 __is >> __mean >> __stddev 02088 >> __x._M_saved_available; 02089 if (__x._M_saved_available) 02090 __is >> __x._M_saved; 02091 __x.param(typename normal_distribution<_RealType>:: 02092 param_type(__mean, __stddev)); 02093 02094 __is.flags(__flags); 02095 return __is; 02096 } 02097 02098 02099 template<typename _RealType> 02100 template<typename _ForwardIterator, 02101 typename _UniformRandomNumberGenerator> 02102 void 02103 lognormal_distribution<_RealType>:: 02104 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02105 _UniformRandomNumberGenerator& __urng, 02106 const param_type& __p) 02107 { 02108 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02109 while (__f != __t) 02110 *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m()); 02111 } 02112 02113 template<typename _RealType, typename _CharT, typename _Traits> 02114 std::basic_ostream<_CharT, _Traits>& 02115 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02116 const lognormal_distribution<_RealType>& __x) 02117 { 02118 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02119 typedef typename __ostream_type::ios_base __ios_base; 02120 02121 const typename __ios_base::fmtflags __flags = __os.flags(); 02122 const _CharT __fill = __os.fill(); 02123 const std::streamsize __precision = __os.precision(); 02124 const _CharT __space = __os.widen(' '); 02125 __os.flags(__ios_base::scientific | __ios_base::left); 02126 __os.fill(__space); 02127 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02128 02129 __os << __x.m() << __space << __x.s() 02130 << __space << __x._M_nd; 02131 02132 __os.flags(__flags); 02133 __os.fill(__fill); 02134 __os.precision(__precision); 02135 return __os; 02136 } 02137 02138 template<typename _RealType, typename _CharT, typename _Traits> 02139 std::basic_istream<_CharT, _Traits>& 02140 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02141 lognormal_distribution<_RealType>& __x) 02142 { 02143 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02144 typedef typename __istream_type::ios_base __ios_base; 02145 02146 const typename __ios_base::fmtflags __flags = __is.flags(); 02147 __is.flags(__ios_base::dec | __ios_base::skipws); 02148 02149 _RealType __m, __s; 02150 __is >> __m >> __s >> __x._M_nd; 02151 __x.param(typename lognormal_distribution<_RealType>:: 02152 param_type(__m, __s)); 02153 02154 __is.flags(__flags); 02155 return __is; 02156 } 02157 02158 template<typename _RealType> 02159 template<typename _ForwardIterator, 02160 typename _UniformRandomNumberGenerator> 02161 void 02162 std::chi_squared_distribution<_RealType>:: 02163 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02164 _UniformRandomNumberGenerator& __urng) 02165 { 02166 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02167 while (__f != __t) 02168 *__f++ = 2 * _M_gd(__urng); 02169 } 02170 02171 template<typename _RealType> 02172 template<typename _ForwardIterator, 02173 typename _UniformRandomNumberGenerator> 02174 void 02175 std::chi_squared_distribution<_RealType>:: 02176 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02177 _UniformRandomNumberGenerator& __urng, 02178 const typename 02179 std::gamma_distribution<result_type>::param_type& __p) 02180 { 02181 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02182 while (__f != __t) 02183 *__f++ = 2 * _M_gd(__urng, __p); 02184 } 02185 02186 template<typename _RealType, typename _CharT, typename _Traits> 02187 std::basic_ostream<_CharT, _Traits>& 02188 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02189 const chi_squared_distribution<_RealType>& __x) 02190 { 02191 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02192 typedef typename __ostream_type::ios_base __ios_base; 02193 02194 const typename __ios_base::fmtflags __flags = __os.flags(); 02195 const _CharT __fill = __os.fill(); 02196 const std::streamsize __precision = __os.precision(); 02197 const _CharT __space = __os.widen(' '); 02198 __os.flags(__ios_base::scientific | __ios_base::left); 02199 __os.fill(__space); 02200 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02201 02202 __os << __x.n() << __space << __x._M_gd; 02203 02204 __os.flags(__flags); 02205 __os.fill(__fill); 02206 __os.precision(__precision); 02207 return __os; 02208 } 02209 02210 template<typename _RealType, typename _CharT, typename _Traits> 02211 std::basic_istream<_CharT, _Traits>& 02212 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02213 chi_squared_distribution<_RealType>& __x) 02214 { 02215 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02216 typedef typename __istream_type::ios_base __ios_base; 02217 02218 const typename __ios_base::fmtflags __flags = __is.flags(); 02219 __is.flags(__ios_base::dec | __ios_base::skipws); 02220 02221 _RealType __n; 02222 __is >> __n >> __x._M_gd; 02223 __x.param(typename chi_squared_distribution<_RealType>:: 02224 param_type(__n)); 02225 02226 __is.flags(__flags); 02227 return __is; 02228 } 02229 02230 02231 template<typename _RealType> 02232 template<typename _UniformRandomNumberGenerator> 02233 typename cauchy_distribution<_RealType>::result_type 02234 cauchy_distribution<_RealType>:: 02235 operator()(_UniformRandomNumberGenerator& __urng, 02236 const param_type& __p) 02237 { 02238 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02239 __aurng(__urng); 02240 _RealType __u; 02241 do 02242 __u = __aurng(); 02243 while (__u == 0.5); 02244 02245 const _RealType __pi = 3.1415926535897932384626433832795029L; 02246 return __p.a() + __p.b() * std::tan(__pi * __u); 02247 } 02248 02249 template<typename _RealType> 02250 template<typename _ForwardIterator, 02251 typename _UniformRandomNumberGenerator> 02252 void 02253 cauchy_distribution<_RealType>:: 02254 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02255 _UniformRandomNumberGenerator& __urng, 02256 const param_type& __p) 02257 { 02258 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02259 const _RealType __pi = 3.1415926535897932384626433832795029L; 02260 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02261 __aurng(__urng); 02262 while (__f != __t) 02263 { 02264 _RealType __u; 02265 do 02266 __u = __aurng(); 02267 while (__u == 0.5); 02268 02269 *__f++ = __p.a() + __p.b() * std::tan(__pi * __u); 02270 } 02271 } 02272 02273 template<typename _RealType, typename _CharT, typename _Traits> 02274 std::basic_ostream<_CharT, _Traits>& 02275 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02276 const cauchy_distribution<_RealType>& __x) 02277 { 02278 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02279 typedef typename __ostream_type::ios_base __ios_base; 02280 02281 const typename __ios_base::fmtflags __flags = __os.flags(); 02282 const _CharT __fill = __os.fill(); 02283 const std::streamsize __precision = __os.precision(); 02284 const _CharT __space = __os.widen(' '); 02285 __os.flags(__ios_base::scientific | __ios_base::left); 02286 __os.fill(__space); 02287 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02288 02289 __os << __x.a() << __space << __x.b(); 02290 02291 __os.flags(__flags); 02292 __os.fill(__fill); 02293 __os.precision(__precision); 02294 return __os; 02295 } 02296 02297 template<typename _RealType, typename _CharT, typename _Traits> 02298 std::basic_istream<_CharT, _Traits>& 02299 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02300 cauchy_distribution<_RealType>& __x) 02301 { 02302 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02303 typedef typename __istream_type::ios_base __ios_base; 02304 02305 const typename __ios_base::fmtflags __flags = __is.flags(); 02306 __is.flags(__ios_base::dec | __ios_base::skipws); 02307 02308 _RealType __a, __b; 02309 __is >> __a >> __b; 02310 __x.param(typename cauchy_distribution<_RealType>:: 02311 param_type(__a, __b)); 02312 02313 __is.flags(__flags); 02314 return __is; 02315 } 02316 02317 02318 template<typename _RealType> 02319 template<typename _ForwardIterator, 02320 typename _UniformRandomNumberGenerator> 02321 void 02322 std::fisher_f_distribution<_RealType>:: 02323 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02324 _UniformRandomNumberGenerator& __urng) 02325 { 02326 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02327 while (__f != __t) 02328 *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m())); 02329 } 02330 02331 template<typename _RealType> 02332 template<typename _ForwardIterator, 02333 typename _UniformRandomNumberGenerator> 02334 void 02335 std::fisher_f_distribution<_RealType>:: 02336 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02337 _UniformRandomNumberGenerator& __urng, 02338 const param_type& __p) 02339 { 02340 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02341 typedef typename std::gamma_distribution<result_type>::param_type 02342 param_type; 02343 param_type __p1(__p.m() / 2); 02344 param_type __p2(__p.n() / 2); 02345 while (__f != __t) 02346 *__f++ = ((_M_gd_x(__urng, __p1) * n()) 02347 / (_M_gd_y(__urng, __p2) * m())); 02348 } 02349 02350 template<typename _RealType, typename _CharT, typename _Traits> 02351 std::basic_ostream<_CharT, _Traits>& 02352 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02353 const fisher_f_distribution<_RealType>& __x) 02354 { 02355 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02356 typedef typename __ostream_type::ios_base __ios_base; 02357 02358 const typename __ios_base::fmtflags __flags = __os.flags(); 02359 const _CharT __fill = __os.fill(); 02360 const std::streamsize __precision = __os.precision(); 02361 const _CharT __space = __os.widen(' '); 02362 __os.flags(__ios_base::scientific | __ios_base::left); 02363 __os.fill(__space); 02364 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02365 02366 __os << __x.m() << __space << __x.n() 02367 << __space << __x._M_gd_x << __space << __x._M_gd_y; 02368 02369 __os.flags(__flags); 02370 __os.fill(__fill); 02371 __os.precision(__precision); 02372 return __os; 02373 } 02374 02375 template<typename _RealType, typename _CharT, typename _Traits> 02376 std::basic_istream<_CharT, _Traits>& 02377 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02378 fisher_f_distribution<_RealType>& __x) 02379 { 02380 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02381 typedef typename __istream_type::ios_base __ios_base; 02382 02383 const typename __ios_base::fmtflags __flags = __is.flags(); 02384 __is.flags(__ios_base::dec | __ios_base::skipws); 02385 02386 _RealType __m, __n; 02387 __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y; 02388 __x.param(typename fisher_f_distribution<_RealType>:: 02389 param_type(__m, __n)); 02390 02391 __is.flags(__flags); 02392 return __is; 02393 } 02394 02395 02396 template<typename _RealType> 02397 template<typename _ForwardIterator, 02398 typename _UniformRandomNumberGenerator> 02399 void 02400 std::student_t_distribution<_RealType>:: 02401 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02402 _UniformRandomNumberGenerator& __urng) 02403 { 02404 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02405 while (__f != __t) 02406 *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng)); 02407 } 02408 02409 template<typename _RealType> 02410 template<typename _ForwardIterator, 02411 typename _UniformRandomNumberGenerator> 02412 void 02413 std::student_t_distribution<_RealType>:: 02414 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02415 _UniformRandomNumberGenerator& __urng, 02416 const param_type& __p) 02417 { 02418 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02419 typename std::gamma_distribution<result_type>::param_type 02420 __p2(__p.n() / 2, 2); 02421 while (__f != __t) 02422 *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2)); 02423 } 02424 02425 template<typename _RealType, typename _CharT, typename _Traits> 02426 std::basic_ostream<_CharT, _Traits>& 02427 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02428 const student_t_distribution<_RealType>& __x) 02429 { 02430 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02431 typedef typename __ostream_type::ios_base __ios_base; 02432 02433 const typename __ios_base::fmtflags __flags = __os.flags(); 02434 const _CharT __fill = __os.fill(); 02435 const std::streamsize __precision = __os.precision(); 02436 const _CharT __space = __os.widen(' '); 02437 __os.flags(__ios_base::scientific | __ios_base::left); 02438 __os.fill(__space); 02439 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02440 02441 __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd; 02442 02443 __os.flags(__flags); 02444 __os.fill(__fill); 02445 __os.precision(__precision); 02446 return __os; 02447 } 02448 02449 template<typename _RealType, typename _CharT, typename _Traits> 02450 std::basic_istream<_CharT, _Traits>& 02451 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02452 student_t_distribution<_RealType>& __x) 02453 { 02454 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02455 typedef typename __istream_type::ios_base __ios_base; 02456 02457 const typename __ios_base::fmtflags __flags = __is.flags(); 02458 __is.flags(__ios_base::dec | __ios_base::skipws); 02459 02460 _RealType __n; 02461 __is >> __n >> __x._M_nd >> __x._M_gd; 02462 __x.param(typename student_t_distribution<_RealType>::param_type(__n)); 02463 02464 __is.flags(__flags); 02465 return __is; 02466 } 02467 02468 02469 template<typename _RealType> 02470 void 02471 gamma_distribution<_RealType>::param_type:: 02472 _M_initialize() 02473 { 02474 _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha; 02475 02476 const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0); 02477 _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1); 02478 } 02479 02480 /** 02481 * Marsaglia, G. and Tsang, W. W. 02482 * "A Simple Method for Generating Gamma Variables" 02483 * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000. 02484 */ 02485 template<typename _RealType> 02486 template<typename _UniformRandomNumberGenerator> 02487 typename gamma_distribution<_RealType>::result_type 02488 gamma_distribution<_RealType>:: 02489 operator()(_UniformRandomNumberGenerator& __urng, 02490 const param_type& __param) 02491 { 02492 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02493 __aurng(__urng); 02494 02495 result_type __u, __v, __n; 02496 const result_type __a1 = (__param._M_malpha 02497 - _RealType(1.0) / _RealType(3.0)); 02498 02499 do 02500 { 02501 do 02502 { 02503 __n = _M_nd(__urng); 02504 __v = result_type(1.0) + __param._M_a2 * __n; 02505 } 02506 while (__v <= 0.0); 02507 02508 __v = __v * __v * __v; 02509 __u = __aurng(); 02510 } 02511 while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 02512 && (std::log(__u) > (0.5 * __n * __n + __a1 02513 * (1.0 - __v + std::log(__v))))); 02514 02515 if (__param.alpha() == __param._M_malpha) 02516 return __a1 * __v * __param.beta(); 02517 else 02518 { 02519 do 02520 __u = __aurng(); 02521 while (__u == 0.0); 02522 02523 return (std::pow(__u, result_type(1.0) / __param.alpha()) 02524 * __a1 * __v * __param.beta()); 02525 } 02526 } 02527 02528 template<typename _RealType> 02529 template<typename _ForwardIterator, 02530 typename _UniformRandomNumberGenerator> 02531 void 02532 gamma_distribution<_RealType>:: 02533 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02534 _UniformRandomNumberGenerator& __urng, 02535 const param_type& __param) 02536 { 02537 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02538 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02539 __aurng(__urng); 02540 02541 result_type __u, __v, __n; 02542 const result_type __a1 = (__param._M_malpha 02543 - _RealType(1.0) / _RealType(3.0)); 02544 02545 if (__param.alpha() == __param._M_malpha) 02546 while (__f != __t) 02547 { 02548 do 02549 { 02550 do 02551 { 02552 __n = _M_nd(__urng); 02553 __v = result_type(1.0) + __param._M_a2 * __n; 02554 } 02555 while (__v <= 0.0); 02556 02557 __v = __v * __v * __v; 02558 __u = __aurng(); 02559 } 02560 while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 02561 && (std::log(__u) > (0.5 * __n * __n + __a1 02562 * (1.0 - __v + std::log(__v))))); 02563 02564 *__f++ = __a1 * __v * __param.beta(); 02565 } 02566 else 02567 while (__f != __t) 02568 { 02569 do 02570 { 02571 do 02572 { 02573 __n = _M_nd(__urng); 02574 __v = result_type(1.0) + __param._M_a2 * __n; 02575 } 02576 while (__v <= 0.0); 02577 02578 __v = __v * __v * __v; 02579 __u = __aurng(); 02580 } 02581 while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 02582 && (std::log(__u) > (0.5 * __n * __n + __a1 02583 * (1.0 - __v + std::log(__v))))); 02584 02585 do 02586 __u = __aurng(); 02587 while (__u == 0.0); 02588 02589 *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha()) 02590 * __a1 * __v * __param.beta()); 02591 } 02592 } 02593 02594 template<typename _RealType, typename _CharT, typename _Traits> 02595 std::basic_ostream<_CharT, _Traits>& 02596 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02597 const gamma_distribution<_RealType>& __x) 02598 { 02599 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02600 typedef typename __ostream_type::ios_base __ios_base; 02601 02602 const typename __ios_base::fmtflags __flags = __os.flags(); 02603 const _CharT __fill = __os.fill(); 02604 const std::streamsize __precision = __os.precision(); 02605 const _CharT __space = __os.widen(' '); 02606 __os.flags(__ios_base::scientific | __ios_base::left); 02607 __os.fill(__space); 02608 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02609 02610 __os << __x.alpha() << __space << __x.beta() 02611 << __space << __x._M_nd; 02612 02613 __os.flags(__flags); 02614 __os.fill(__fill); 02615 __os.precision(__precision); 02616 return __os; 02617 } 02618 02619 template<typename _RealType, typename _CharT, typename _Traits> 02620 std::basic_istream<_CharT, _Traits>& 02621 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02622 gamma_distribution<_RealType>& __x) 02623 { 02624 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02625 typedef typename __istream_type::ios_base __ios_base; 02626 02627 const typename __ios_base::fmtflags __flags = __is.flags(); 02628 __is.flags(__ios_base::dec | __ios_base::skipws); 02629 02630 _RealType __alpha_val, __beta_val; 02631 __is >> __alpha_val >> __beta_val >> __x._M_nd; 02632 __x.param(typename gamma_distribution<_RealType>:: 02633 param_type(__alpha_val, __beta_val)); 02634 02635 __is.flags(__flags); 02636 return __is; 02637 } 02638 02639 02640 template<typename _RealType> 02641 template<typename _UniformRandomNumberGenerator> 02642 typename weibull_distribution<_RealType>::result_type 02643 weibull_distribution<_RealType>:: 02644 operator()(_UniformRandomNumberGenerator& __urng, 02645 const param_type& __p) 02646 { 02647 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02648 __aurng(__urng); 02649 return __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 02650 result_type(1) / __p.a()); 02651 } 02652 02653 template<typename _RealType> 02654 template<typename _ForwardIterator, 02655 typename _UniformRandomNumberGenerator> 02656 void 02657 weibull_distribution<_RealType>:: 02658 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02659 _UniformRandomNumberGenerator& __urng, 02660 const param_type& __p) 02661 { 02662 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02663 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02664 __aurng(__urng); 02665 auto __inv_a = result_type(1) / __p.a(); 02666 02667 while (__f != __t) 02668 *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 02669 __inv_a); 02670 } 02671 02672 template<typename _RealType, typename _CharT, typename _Traits> 02673 std::basic_ostream<_CharT, _Traits>& 02674 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02675 const weibull_distribution<_RealType>& __x) 02676 { 02677 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02678 typedef typename __ostream_type::ios_base __ios_base; 02679 02680 const typename __ios_base::fmtflags __flags = __os.flags(); 02681 const _CharT __fill = __os.fill(); 02682 const std::streamsize __precision = __os.precision(); 02683 const _CharT __space = __os.widen(' '); 02684 __os.flags(__ios_base::scientific | __ios_base::left); 02685 __os.fill(__space); 02686 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02687 02688 __os << __x.a() << __space << __x.b(); 02689 02690 __os.flags(__flags); 02691 __os.fill(__fill); 02692 __os.precision(__precision); 02693 return __os; 02694 } 02695 02696 template<typename _RealType, typename _CharT, typename _Traits> 02697 std::basic_istream<_CharT, _Traits>& 02698 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02699 weibull_distribution<_RealType>& __x) 02700 { 02701 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02702 typedef typename __istream_type::ios_base __ios_base; 02703 02704 const typename __ios_base::fmtflags __flags = __is.flags(); 02705 __is.flags(__ios_base::dec | __ios_base::skipws); 02706 02707 _RealType __a, __b; 02708 __is >> __a >> __b; 02709 __x.param(typename weibull_distribution<_RealType>:: 02710 param_type(__a, __b)); 02711 02712 __is.flags(__flags); 02713 return __is; 02714 } 02715 02716 02717 template<typename _RealType> 02718 template<typename _UniformRandomNumberGenerator> 02719 typename extreme_value_distribution<_RealType>::result_type 02720 extreme_value_distribution<_RealType>:: 02721 operator()(_UniformRandomNumberGenerator& __urng, 02722 const param_type& __p) 02723 { 02724 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02725 __aurng(__urng); 02726 return __p.a() - __p.b() * std::log(-std::log(result_type(1) 02727 - __aurng())); 02728 } 02729 02730 template<typename _RealType> 02731 template<typename _ForwardIterator, 02732 typename _UniformRandomNumberGenerator> 02733 void 02734 extreme_value_distribution<_RealType>:: 02735 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02736 _UniformRandomNumberGenerator& __urng, 02737 const param_type& __p) 02738 { 02739 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02740 __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 02741 __aurng(__urng); 02742 02743 while (__f != __t) 02744 *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1) 02745 - __aurng())); 02746 } 02747 02748 template<typename _RealType, typename _CharT, typename _Traits> 02749 std::basic_ostream<_CharT, _Traits>& 02750 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02751 const extreme_value_distribution<_RealType>& __x) 02752 { 02753 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02754 typedef typename __ostream_type::ios_base __ios_base; 02755 02756 const typename __ios_base::fmtflags __flags = __os.flags(); 02757 const _CharT __fill = __os.fill(); 02758 const std::streamsize __precision = __os.precision(); 02759 const _CharT __space = __os.widen(' '); 02760 __os.flags(__ios_base::scientific | __ios_base::left); 02761 __os.fill(__space); 02762 __os.precision(std::numeric_limits<_RealType>::max_digits10); 02763 02764 __os << __x.a() << __space << __x.b(); 02765 02766 __os.flags(__flags); 02767 __os.fill(__fill); 02768 __os.precision(__precision); 02769 return __os; 02770 } 02771 02772 template<typename _RealType, typename _CharT, typename _Traits> 02773 std::basic_istream<_CharT, _Traits>& 02774 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02775 extreme_value_distribution<_RealType>& __x) 02776 { 02777 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02778 typedef typename __istream_type::ios_base __ios_base; 02779 02780 const typename __ios_base::fmtflags __flags = __is.flags(); 02781 __is.flags(__ios_base::dec | __ios_base::skipws); 02782 02783 _RealType __a, __b; 02784 __is >> __a >> __b; 02785 __x.param(typename extreme_value_distribution<_RealType>:: 02786 param_type(__a, __b)); 02787 02788 __is.flags(__flags); 02789 return __is; 02790 } 02791 02792 02793 template<typename _IntType> 02794 void 02795 discrete_distribution<_IntType>::param_type:: 02796 _M_initialize() 02797 { 02798 if (_M_prob.size() < 2) 02799 { 02800 _M_prob.clear(); 02801 return; 02802 } 02803 02804 const double __sum = std::accumulate(_M_prob.begin(), 02805 _M_prob.end(), 0.0); 02806 // Now normalize the probabilites. 02807 __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), 02808 __sum); 02809 // Accumulate partial sums. 02810 _M_cp.reserve(_M_prob.size()); 02811 std::partial_sum(_M_prob.begin(), _M_prob.end(), 02812 std::back_inserter(_M_cp)); 02813 // Make sure the last cumulative probability is one. 02814 _M_cp[_M_cp.size() - 1] = 1.0; 02815 } 02816 02817 template<typename _IntType> 02818 template<typename _Func> 02819 discrete_distribution<_IntType>::param_type:: 02820 param_type(size_t __nw, double __xmin, double __xmax, _Func __fw) 02821 : _M_prob(), _M_cp() 02822 { 02823 const size_t __n = __nw == 0 ? 1 : __nw; 02824 const double __delta = (__xmax - __xmin) / __n; 02825 02826 _M_prob.reserve(__n); 02827 for (size_t __k = 0; __k < __nw; ++__k) 02828 _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta)); 02829 02830 _M_initialize(); 02831 } 02832 02833 template<typename _IntType> 02834 template<typename _UniformRandomNumberGenerator> 02835 typename discrete_distribution<_IntType>::result_type 02836 discrete_distribution<_IntType>:: 02837 operator()(_UniformRandomNumberGenerator& __urng, 02838 const param_type& __param) 02839 { 02840 if (__param._M_cp.empty()) 02841 return result_type(0); 02842 02843 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02844 __aurng(__urng); 02845 02846 const double __p = __aurng(); 02847 auto __pos = std::lower_bound(__param._M_cp.begin(), 02848 __param._M_cp.end(), __p); 02849 02850 return __pos - __param._M_cp.begin(); 02851 } 02852 02853 template<typename _IntType> 02854 template<typename _ForwardIterator, 02855 typename _UniformRandomNumberGenerator> 02856 void 02857 discrete_distribution<_IntType>:: 02858 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 02859 _UniformRandomNumberGenerator& __urng, 02860 const param_type& __param) 02861 { 02862 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 02863 02864 if (__param._M_cp.empty()) 02865 { 02866 while (__f != __t) 02867 *__f++ = result_type(0); 02868 return; 02869 } 02870 02871 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 02872 __aurng(__urng); 02873 02874 while (__f != __t) 02875 { 02876 const double __p = __aurng(); 02877 auto __pos = std::lower_bound(__param._M_cp.begin(), 02878 __param._M_cp.end(), __p); 02879 02880 *__f++ = __pos - __param._M_cp.begin(); 02881 } 02882 } 02883 02884 template<typename _IntType, typename _CharT, typename _Traits> 02885 std::basic_ostream<_CharT, _Traits>& 02886 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 02887 const discrete_distribution<_IntType>& __x) 02888 { 02889 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 02890 typedef typename __ostream_type::ios_base __ios_base; 02891 02892 const typename __ios_base::fmtflags __flags = __os.flags(); 02893 const _CharT __fill = __os.fill(); 02894 const std::streamsize __precision = __os.precision(); 02895 const _CharT __space = __os.widen(' '); 02896 __os.flags(__ios_base::scientific | __ios_base::left); 02897 __os.fill(__space); 02898 __os.precision(std::numeric_limits<double>::max_digits10); 02899 02900 std::vector<double> __prob = __x.probabilities(); 02901 __os << __prob.size(); 02902 for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit) 02903 __os << __space << *__dit; 02904 02905 __os.flags(__flags); 02906 __os.fill(__fill); 02907 __os.precision(__precision); 02908 return __os; 02909 } 02910 02911 template<typename _IntType, typename _CharT, typename _Traits> 02912 std::basic_istream<_CharT, _Traits>& 02913 operator>>(std::basic_istream<_CharT, _Traits>& __is, 02914 discrete_distribution<_IntType>& __x) 02915 { 02916 typedef std::basic_istream<_CharT, _Traits> __istream_type; 02917 typedef typename __istream_type::ios_base __ios_base; 02918 02919 const typename __ios_base::fmtflags __flags = __is.flags(); 02920 __is.flags(__ios_base::dec | __ios_base::skipws); 02921 02922 size_t __n; 02923 __is >> __n; 02924 02925 std::vector<double> __prob_vec; 02926 __prob_vec.reserve(__n); 02927 for (; __n != 0; --__n) 02928 { 02929 double __prob; 02930 __is >> __prob; 02931 __prob_vec.push_back(__prob); 02932 } 02933 02934 __x.param(typename discrete_distribution<_IntType>:: 02935 param_type(__prob_vec.begin(), __prob_vec.end())); 02936 02937 __is.flags(__flags); 02938 return __is; 02939 } 02940 02941 02942 template<typename _RealType> 02943 void 02944 piecewise_constant_distribution<_RealType>::param_type:: 02945 _M_initialize() 02946 { 02947 if (_M_int.size() < 2 02948 || (_M_int.size() == 2 02949 && _M_int[0] == _RealType(0) 02950 && _M_int[1] == _RealType(1))) 02951 { 02952 _M_int.clear(); 02953 _M_den.clear(); 02954 return; 02955 } 02956 02957 const double __sum = std::accumulate(_M_den.begin(), 02958 _M_den.end(), 0.0); 02959 02960 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 02961 __sum); 02962 02963 _M_cp.reserve(_M_den.size()); 02964 std::partial_sum(_M_den.begin(), _M_den.end(), 02965 std::back_inserter(_M_cp)); 02966 02967 // Make sure the last cumulative probability is one. 02968 _M_cp[_M_cp.size() - 1] = 1.0; 02969 02970 for (size_t __k = 0; __k < _M_den.size(); ++__k) 02971 _M_den[__k] /= _M_int[__k + 1] - _M_int[__k]; 02972 } 02973 02974 template<typename _RealType> 02975 template<typename _InputIteratorB, typename _InputIteratorW> 02976 piecewise_constant_distribution<_RealType>::param_type:: 02977 param_type(_InputIteratorB __bbegin, 02978 _InputIteratorB __bend, 02979 _InputIteratorW __wbegin) 02980 : _M_int(), _M_den(), _M_cp() 02981 { 02982 if (__bbegin != __bend) 02983 { 02984 for (;;) 02985 { 02986 _M_int.push_back(*__bbegin); 02987 ++__bbegin; 02988 if (__bbegin == __bend) 02989 break; 02990 02991 _M_den.push_back(*__wbegin); 02992 ++__wbegin; 02993 } 02994 } 02995 02996 _M_initialize(); 02997 } 02998 02999 template<typename _RealType> 03000 template<typename _Func> 03001 piecewise_constant_distribution<_RealType>::param_type:: 03002 param_type(initializer_list<_RealType> __bl, _Func __fw) 03003 : _M_int(), _M_den(), _M_cp() 03004 { 03005 _M_int.reserve(__bl.size()); 03006 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 03007 _M_int.push_back(*__biter); 03008 03009 _M_den.reserve(_M_int.size() - 1); 03010 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 03011 _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k]))); 03012 03013 _M_initialize(); 03014 } 03015 03016 template<typename _RealType> 03017 template<typename _Func> 03018 piecewise_constant_distribution<_RealType>::param_type:: 03019 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 03020 : _M_int(), _M_den(), _M_cp() 03021 { 03022 const size_t __n = __nw == 0 ? 1 : __nw; 03023 const _RealType __delta = (__xmax - __xmin) / __n; 03024 03025 _M_int.reserve(__n + 1); 03026 for (size_t __k = 0; __k <= __nw; ++__k) 03027 _M_int.push_back(__xmin + __k * __delta); 03028 03029 _M_den.reserve(__n); 03030 for (size_t __k = 0; __k < __nw; ++__k) 03031 _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta)); 03032 03033 _M_initialize(); 03034 } 03035 03036 template<typename _RealType> 03037 template<typename _UniformRandomNumberGenerator> 03038 typename piecewise_constant_distribution<_RealType>::result_type 03039 piecewise_constant_distribution<_RealType>:: 03040 operator()(_UniformRandomNumberGenerator& __urng, 03041 const param_type& __param) 03042 { 03043 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 03044 __aurng(__urng); 03045 03046 const double __p = __aurng(); 03047 if (__param._M_cp.empty()) 03048 return __p; 03049 03050 auto __pos = std::lower_bound(__param._M_cp.begin(), 03051 __param._M_cp.end(), __p); 03052 const size_t __i = __pos - __param._M_cp.begin(); 03053 03054 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 03055 03056 return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i]; 03057 } 03058 03059 template<typename _RealType> 03060 template<typename _ForwardIterator, 03061 typename _UniformRandomNumberGenerator> 03062 void 03063 piecewise_constant_distribution<_RealType>:: 03064 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 03065 _UniformRandomNumberGenerator& __urng, 03066 const param_type& __param) 03067 { 03068 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 03069 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 03070 __aurng(__urng); 03071 03072 if (__param._M_cp.empty()) 03073 { 03074 while (__f != __t) 03075 *__f++ = __aurng(); 03076 return; 03077 } 03078 03079 while (__f != __t) 03080 { 03081 const double __p = __aurng(); 03082 03083 auto __pos = std::lower_bound(__param._M_cp.begin(), 03084 __param._M_cp.end(), __p); 03085 const size_t __i = __pos - __param._M_cp.begin(); 03086 03087 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 03088 03089 *__f++ = (__param._M_int[__i] 03090 + (__p - __pref) / __param._M_den[__i]); 03091 } 03092 } 03093 03094 template<typename _RealType, typename _CharT, typename _Traits> 03095 std::basic_ostream<_CharT, _Traits>& 03096 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 03097 const piecewise_constant_distribution<_RealType>& __x) 03098 { 03099 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 03100 typedef typename __ostream_type::ios_base __ios_base; 03101 03102 const typename __ios_base::fmtflags __flags = __os.flags(); 03103 const _CharT __fill = __os.fill(); 03104 const std::streamsize __precision = __os.precision(); 03105 const _CharT __space = __os.widen(' '); 03106 __os.flags(__ios_base::scientific | __ios_base::left); 03107 __os.fill(__space); 03108 __os.precision(std::numeric_limits<_RealType>::max_digits10); 03109 03110 std::vector<_RealType> __int = __x.intervals(); 03111 __os << __int.size() - 1; 03112 03113 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 03114 __os << __space << *__xit; 03115 03116 std::vector<double> __den = __x.densities(); 03117 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 03118 __os << __space << *__dit; 03119 03120 __os.flags(__flags); 03121 __os.fill(__fill); 03122 __os.precision(__precision); 03123 return __os; 03124 } 03125 03126 template<typename _RealType, typename _CharT, typename _Traits> 03127 std::basic_istream<_CharT, _Traits>& 03128 operator>>(std::basic_istream<_CharT, _Traits>& __is, 03129 piecewise_constant_distribution<_RealType>& __x) 03130 { 03131 typedef std::basic_istream<_CharT, _Traits> __istream_type; 03132 typedef typename __istream_type::ios_base __ios_base; 03133 03134 const typename __ios_base::fmtflags __flags = __is.flags(); 03135 __is.flags(__ios_base::dec | __ios_base::skipws); 03136 03137 size_t __n; 03138 __is >> __n; 03139 03140 std::vector<_RealType> __int_vec; 03141 __int_vec.reserve(__n + 1); 03142 for (size_t __i = 0; __i <= __n; ++__i) 03143 { 03144 _RealType __int; 03145 __is >> __int; 03146 __int_vec.push_back(__int); 03147 } 03148 03149 std::vector<double> __den_vec; 03150 __den_vec.reserve(__n); 03151 for (size_t __i = 0; __i < __n; ++__i) 03152 { 03153 double __den; 03154 __is >> __den; 03155 __den_vec.push_back(__den); 03156 } 03157 03158 __x.param(typename piecewise_constant_distribution<_RealType>:: 03159 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 03160 03161 __is.flags(__flags); 03162 return __is; 03163 } 03164 03165 03166 template<typename _RealType> 03167 void 03168 piecewise_linear_distribution<_RealType>::param_type:: 03169 _M_initialize() 03170 { 03171 if (_M_int.size() < 2 03172 || (_M_int.size() == 2 03173 && _M_int[0] == _RealType(0) 03174 && _M_int[1] == _RealType(1) 03175 && _M_den[0] == _M_den[1])) 03176 { 03177 _M_int.clear(); 03178 _M_den.clear(); 03179 return; 03180 } 03181 03182 double __sum = 0.0; 03183 _M_cp.reserve(_M_int.size() - 1); 03184 _M_m.reserve(_M_int.size() - 1); 03185 for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 03186 { 03187 const _RealType __delta = _M_int[__k + 1] - _M_int[__k]; 03188 __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta; 03189 _M_cp.push_back(__sum); 03190 _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta); 03191 } 03192 03193 // Now normalize the densities... 03194 __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), 03195 __sum); 03196 // ... and partial sums... 03197 __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum); 03198 // ... and slopes. 03199 __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum); 03200 03201 // Make sure the last cumulative probablility is one. 03202 _M_cp[_M_cp.size() - 1] = 1.0; 03203 } 03204 03205 template<typename _RealType> 03206 template<typename _InputIteratorB, typename _InputIteratorW> 03207 piecewise_linear_distribution<_RealType>::param_type:: 03208 param_type(_InputIteratorB __bbegin, 03209 _InputIteratorB __bend, 03210 _InputIteratorW __wbegin) 03211 : _M_int(), _M_den(), _M_cp(), _M_m() 03212 { 03213 for (; __bbegin != __bend; ++__bbegin, ++__wbegin) 03214 { 03215 _M_int.push_back(*__bbegin); 03216 _M_den.push_back(*__wbegin); 03217 } 03218 03219 _M_initialize(); 03220 } 03221 03222 template<typename _RealType> 03223 template<typename _Func> 03224 piecewise_linear_distribution<_RealType>::param_type:: 03225 param_type(initializer_list<_RealType> __bl, _Func __fw) 03226 : _M_int(), _M_den(), _M_cp(), _M_m() 03227 { 03228 _M_int.reserve(__bl.size()); 03229 _M_den.reserve(__bl.size()); 03230 for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 03231 { 03232 _M_int.push_back(*__biter); 03233 _M_den.push_back(__fw(*__biter)); 03234 } 03235 03236 _M_initialize(); 03237 } 03238 03239 template<typename _RealType> 03240 template<typename _Func> 03241 piecewise_linear_distribution<_RealType>::param_type:: 03242 param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 03243 : _M_int(), _M_den(), _M_cp(), _M_m() 03244 { 03245 const size_t __n = __nw == 0 ? 1 : __nw; 03246 const _RealType __delta = (__xmax - __xmin) / __n; 03247 03248 _M_int.reserve(__n + 1); 03249 _M_den.reserve(__n + 1); 03250 for (size_t __k = 0; __k <= __nw; ++__k) 03251 { 03252 _M_int.push_back(__xmin + __k * __delta); 03253 _M_den.push_back(__fw(_M_int[__k] + __delta)); 03254 } 03255 03256 _M_initialize(); 03257 } 03258 03259 template<typename _RealType> 03260 template<typename _UniformRandomNumberGenerator> 03261 typename piecewise_linear_distribution<_RealType>::result_type 03262 piecewise_linear_distribution<_RealType>:: 03263 operator()(_UniformRandomNumberGenerator& __urng, 03264 const param_type& __param) 03265 { 03266 __detail::_Adaptor<_UniformRandomNumberGenerator, double> 03267 __aurng(__urng); 03268 03269 const double __p = __aurng(); 03270 if (__param._M_cp.empty()) 03271 return __p; 03272 03273 auto __pos = std::lower_bound(__param._M_cp.begin(), 03274 __param._M_cp.end(), __p); 03275 const size_t __i = __pos - __param._M_cp.begin(); 03276 03277 const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 03278 03279 const double __a = 0.5 * __param._M_m[__i]; 03280 const double __b = __param._M_den[__i]; 03281 const double __cm = __p - __pref; 03282 03283 _RealType __x = __param._M_int[__i]; 03284 if (__a == 0) 03285 __x += __cm / __b; 03286 else 03287 { 03288 const double __d = __b * __b + 4.0 * __a * __cm; 03289 __x += 0.5 * (std::sqrt(__d) - __b) / __a; 03290 } 03291 03292 return __x; 03293 } 03294 03295 template<typename _RealType> 03296 template<typename _ForwardIterator, 03297 typename _UniformRandomNumberGenerator> 03298 void 03299 piecewise_linear_distribution<_RealType>:: 03300 __generate_impl(_ForwardIterator __f, _ForwardIterator __t, 03301 _UniformRandomNumberGenerator& __urng, 03302 const param_type& __param) 03303 { 03304 __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) 03305 // We could duplicate everything from operator()... 03306 while (__f != __t) 03307 *__f++ = this->operator()(__urng, __param); 03308 } 03309 03310 template<typename _RealType, typename _CharT, typename _Traits> 03311 std::basic_ostream<_CharT, _Traits>& 03312 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 03313 const piecewise_linear_distribution<_RealType>& __x) 03314 { 03315 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 03316 typedef typename __ostream_type::ios_base __ios_base; 03317 03318 const typename __ios_base::fmtflags __flags = __os.flags(); 03319 const _CharT __fill = __os.fill(); 03320 const std::streamsize __precision = __os.precision(); 03321 const _CharT __space = __os.widen(' '); 03322 __os.flags(__ios_base::scientific | __ios_base::left); 03323 __os.fill(__space); 03324 __os.precision(std::numeric_limits<_RealType>::max_digits10); 03325 03326 std::vector<_RealType> __int = __x.intervals(); 03327 __os << __int.size() - 1; 03328 03329 for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 03330 __os << __space << *__xit; 03331 03332 std::vector<double> __den = __x.densities(); 03333 for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 03334 __os << __space << *__dit; 03335 03336 __os.flags(__flags); 03337 __os.fill(__fill); 03338 __os.precision(__precision); 03339 return __os; 03340 } 03341 03342 template<typename _RealType, typename _CharT, typename _Traits> 03343 std::basic_istream<_CharT, _Traits>& 03344 operator>>(std::basic_istream<_CharT, _Traits>& __is, 03345 piecewise_linear_distribution<_RealType>& __x) 03346 { 03347 typedef std::basic_istream<_CharT, _Traits> __istream_type; 03348 typedef typename __istream_type::ios_base __ios_base; 03349 03350 const typename __ios_base::fmtflags __flags = __is.flags(); 03351 __is.flags(__ios_base::dec | __ios_base::skipws); 03352 03353 size_t __n; 03354 __is >> __n; 03355 03356 std::vector<_RealType> __int_vec; 03357 __int_vec.reserve(__n + 1); 03358 for (size_t __i = 0; __i <= __n; ++__i) 03359 { 03360 _RealType __int; 03361 __is >> __int; 03362 __int_vec.push_back(__int); 03363 } 03364 03365 std::vector<double> __den_vec; 03366 __den_vec.reserve(__n + 1); 03367 for (size_t __i = 0; __i <= __n; ++__i) 03368 { 03369 double __den; 03370 __is >> __den; 03371 __den_vec.push_back(__den); 03372 } 03373 03374 __x.param(typename piecewise_linear_distribution<_RealType>:: 03375 param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 03376 03377 __is.flags(__flags); 03378 return __is; 03379 } 03380 03381 03382 template<typename _IntType> 03383 seed_seq::seed_seq(std::initializer_list<_IntType> __il) 03384 { 03385 for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter) 03386 _M_v.push_back(__detail::__mod<result_type, 03387 __detail::_Shift<result_type, 32>::__value>(*__iter)); 03388 } 03389 03390 template<typename _InputIterator> 03391 seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end) 03392 { 03393 for (_InputIterator __iter = __begin; __iter != __end; ++__iter) 03394 _M_v.push_back(__detail::__mod<result_type, 03395 __detail::_Shift<result_type, 32>::__value>(*__iter)); 03396 } 03397 03398 template<typename _RandomAccessIterator> 03399 void 03400 seed_seq::generate(_RandomAccessIterator __begin, 03401 _RandomAccessIterator __end) 03402 { 03403 typedef typename iterator_traits<_RandomAccessIterator>::value_type 03404 _Type; 03405 03406 if (__begin == __end) 03407 return; 03408 03409 std::fill(__begin, __end, _Type(0x8b8b8b8bu)); 03410 03411 const size_t __n = __end - __begin; 03412 const size_t __s = _M_v.size(); 03413 const size_t __t = (__n >= 623) ? 11 03414 : (__n >= 68) ? 7 03415 : (__n >= 39) ? 5 03416 : (__n >= 7) ? 3 03417 : (__n - 1) / 2; 03418 const size_t __p = (__n - __t) / 2; 03419 const size_t __q = __p + __t; 03420 const size_t __m = std::max(size_t(__s + 1), __n); 03421 03422 for (size_t __k = 0; __k < __m; ++__k) 03423 { 03424 _Type __arg = (__begin[__k % __n] 03425 ^ __begin[(__k + __p) % __n] 03426 ^ __begin[(__k - 1) % __n]); 03427 _Type __r1 = __arg ^ (__arg >> 27); 03428 __r1 = __detail::__mod<_Type, 03429 __detail::_Shift<_Type, 32>::__value>(1664525u * __r1); 03430 _Type __r2 = __r1; 03431 if (__k == 0) 03432 __r2 += __s; 03433 else if (__k <= __s) 03434 __r2 += __k % __n + _M_v[__k - 1]; 03435 else 03436 __r2 += __k % __n; 03437 __r2 = __detail::__mod<_Type, 03438 __detail::_Shift<_Type, 32>::__value>(__r2); 03439 __begin[(__k + __p) % __n] += __r1; 03440 __begin[(__k + __q) % __n] += __r2; 03441 __begin[__k % __n] = __r2; 03442 } 03443 03444 for (size_t __k = __m; __k < __m + __n; ++__k) 03445 { 03446 _Type __arg = (__begin[__k % __n] 03447 + __begin[(__k + __p) % __n] 03448 + __begin[(__k - 1) % __n]); 03449 _Type __r3 = __arg ^ (__arg >> 27); 03450 __r3 = __detail::__mod<_Type, 03451 __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3); 03452 _Type __r4 = __r3 - __k % __n; 03453 __r4 = __detail::__mod<_Type, 03454 __detail::_Shift<_Type, 32>::__value>(__r4); 03455 __begin[(__k + __p) % __n] ^= __r3; 03456 __begin[(__k + __q) % __n] ^= __r4; 03457 __begin[__k % __n] = __r4; 03458 } 03459 } 03460 03461 template<typename _RealType, size_t __bits, 03462 typename _UniformRandomNumberGenerator> 03463 _RealType 03464 generate_canonical(_UniformRandomNumberGenerator& __urng) 03465 { 03466 const size_t __b 03467 = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits), 03468 __bits); 03469 const long double __r = static_cast<long double>(__urng.max()) 03470 - static_cast<long double>(__urng.min()) + 1.0L; 03471 const size_t __log2r = std::log(__r) / std::log(2.0L); 03472 size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r); 03473 _RealType __sum = _RealType(0); 03474 _RealType __tmp = _RealType(1); 03475 for (; __k != 0; --__k) 03476 { 03477 __sum += _RealType(__urng() - __urng.min()) * __tmp; 03478 __tmp *= __r; 03479 } 03480 return __sum / __tmp; 03481 } 03482 03483 _GLIBCXX_END_NAMESPACE_VERSION 03484 } // namespace 03485 03486 #endif