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