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