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