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