1233294Sstas// random number generation (out of line) -*- C++ -*- 2102644Snectar 357416Smarkm// Copyright (C) 2009-2022 Free Software Foundation, Inc. 4142403Snectar// 5233294Sstas// This file is part of the GNU ISO C++ Library. This library is free 6233294Sstas// software; you can redistribute it and/or modify it under the 757416Smarkm// terms of the GNU General Public License as published by the 857416Smarkm// Free Software Foundation; either version 3, or (at your option) 957416Smarkm// any later version. 1057416Smarkm 1157416Smarkm// This library is distributed in the hope that it will be useful, 1257416Smarkm// but WITHOUT ANY WARRANTY; without even the implied warranty of 1357416Smarkm// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1457416Smarkm// GNU General Public License for more details. 1557416Smarkm 1690926Snectar// Under Section 7 of GPL version 3, you are granted additional 1790926Snectar// permissions described in the GCC Runtime Library Exception, version 18233294Sstas// 3.1, as published by the Free Software Foundation. 1990926Snectar 20233294Sstas// You should have received a copy of the GNU General Public License and 2190926Snectar// a copy of the GCC Runtime Library Exception along with this program; 22233294Sstas// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 2357416Smarkm// <http://www.gnu.org/licenses/>. 2457416Smarkm 2557416Smarkm 26233294Sstas/** @file tr1/random.tcc 2757416Smarkm * This is an internal header file, included by other library headers. 28233294Sstas * Do not attempt to use it directly. @headername{tr1/random} 29102644Snectar */ 30102644Snectar 31102644Snectar#ifndef _GLIBCXX_TR1_RANDOM_TCC 32127808Snectar#define _GLIBCXX_TR1_RANDOM_TCC 1 3390926Snectar 34127808Snectarnamespace std _GLIBCXX_VISIBILITY(default) 3557416Smarkm{ 3657416Smarkm_GLIBCXX_BEGIN_NAMESPACE_VERSION 3757416Smarkm 3857416Smarkmnamespace tr1 3957416Smarkm{ 4057416Smarkm /* 41178825Sdfr * (Further) implementation-space details. 4257416Smarkm */ 43142403Snectar namespace __detail 44142403Snectar { 45142403Snectar // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid 46142403Snectar // integer overflow. 47142403Snectar // 48142403Snectar // Because a and c are compile-time integral constants the compiler kindly 49142403Snectar // elides any unreachable paths. 50233294Sstas // 51142403Snectar // Preconditions: a > 0, m > 0. 52142403Snectar // 53142403Snectar template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool> 54142403Snectar struct _Mod 55142403Snectar { 56142403Snectar static _Tp 57142403Snectar __calc(_Tp __x) 58142403Snectar { 59142403Snectar if (__a == 1) 60142403Snectar __x %= __m; 61142403Snectar else 62142403Snectar { 63142403Snectar static const _Tp __q = __m / __a; 64142403Snectar static const _Tp __r = __m % __a; 65233294Sstas 66142403Snectar _Tp __t1 = __a * (__x % __q); 67142403Snectar _Tp __t2 = __r * (__x / __q); 68142403Snectar if (__t1 >= __t2) 69142403Snectar __x = __t1 - __t2; 70178825Sdfr else 71142403Snectar __x = __m - __t2 + __t1; 72142403Snectar } 73142403Snectar 74142403Snectar if (__c != 0) 75142403Snectar { 76142403Snectar const _Tp __d = __m - __x; 77142403Snectar if (__d > __c) 78142403Snectar __x += __c; 79233294Sstas else 80233294Sstas __x = __c - __d; 81233294Sstas } 82233294Sstas return __x; 83233294Sstas } 84233294Sstas }; 85178825Sdfr 86178825Sdfr // Special case for m == 0 -- use unsigned integer overflow as modulo 87178825Sdfr // operator. 88178825Sdfr template<typename _Tp, _Tp __a, _Tp __c, _Tp __m> 89178825Sdfr struct _Mod<_Tp, __a, __c, __m, true> 90178825Sdfr { 91178825Sdfr static _Tp 92233294Sstas __calc(_Tp __x) 93142403Snectar { return __a * __x + __c; } 94142403Snectar }; 95178825Sdfr } // namespace __detail 96142403Snectar 97142403Snectar template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 98233294Sstas const _UIntType 99142403Snectar linear_congruential<_UIntType, __a, __c, __m>::multiplier; 100142403Snectar 101142403Snectar template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 102142403Snectar const _UIntType 103142403Snectar linear_congruential<_UIntType, __a, __c, __m>::increment; 104142403Snectar 105142403Snectar template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 106178825Sdfr const _UIntType 107178825Sdfr linear_congruential<_UIntType, __a, __c, __m>::modulus; 108178825Sdfr 109178825Sdfr /** 110233294Sstas * Seeds the LCR with integral value @p __x0, adjusted so that the 111233294Sstas * ring identity is never a member of the convergence set. 112233294Sstas */ 113233294Sstas template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 114142403Snectar void 115142403Snectar linear_congruential<_UIntType, __a, __c, __m>:: 116178825Sdfr seed(unsigned long __x0) 117178825Sdfr { 118178825Sdfr if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) 119142403Snectar && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) 120178825Sdfr _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); 121178825Sdfr else 122178825Sdfr _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); 123142403Snectar } 124142403Snectar 125233294Sstas /** 126233294Sstas * Seeds the LCR engine with a value generated by @p __g. 127233294Sstas */ 128233294Sstas template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 129233294Sstas template<class _Gen> 130233294Sstas void 131233294Sstas linear_congruential<_UIntType, __a, __c, __m>:: 132233294Sstas seed(_Gen& __g, false_type) 133233294Sstas { 134233294Sstas _UIntType __x0 = __g(); 135233294Sstas if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) 136233294Sstas && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) 137233294Sstas _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); 138233294Sstas else 139233294Sstas _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); 140233294Sstas } 141233294Sstas 142233294Sstas /** 143233294Sstas * Gets the next generated value in sequence. 144233294Sstas */ 145233294Sstas template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 146142403Snectar typename linear_congruential<_UIntType, __a, __c, __m>::result_type 147142403Snectar linear_congruential<_UIntType, __a, __c, __m>:: 148142403Snectar operator()() 149142403Snectar { 150142403Snectar _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x); 151127808Snectar return _M_x; 15257416Smarkm } 15372445Sassar 154127808Snectar template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 155233294Sstas typename _CharT, typename _Traits> 156233294Sstas std::basic_ostream<_CharT, _Traits>& 157127808Snectar operator<<(std::basic_ostream<_CharT, _Traits>& __os, 158127808Snectar const linear_congruential<_UIntType, __a, __c, __m>& __lcr) 159127808Snectar { 16057416Smarkm typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 16157416Smarkm typedef typename __ostream_type::ios_base __ios_base; 162233294Sstas 163233294Sstas const typename __ios_base::fmtflags __flags = __os.flags(); 16457416Smarkm const _CharT __fill = __os.fill(); 16557416Smarkm __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 16657416Smarkm __os.fill(__os.widen(' ')); 167233294Sstas 168127808Snectar __os << __lcr._M_x; 16990926Snectar 17072445Sassar __os.flags(__flags); 171127808Snectar __os.fill(__fill); 172127808Snectar return __os; 173233294Sstas } 17457416Smarkm 175127808Snectar template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 176233294Sstas typename _CharT, typename _Traits> 17790926Snectar std::basic_istream<_CharT, _Traits>& 178178825Sdfr operator>>(std::basic_istream<_CharT, _Traits>& __is, 179178825Sdfr linear_congruential<_UIntType, __a, __c, __m>& __lcr) 18072445Sassar { 181233294Sstas typedef std::basic_istream<_CharT, _Traits> __istream_type; 182233294Sstas typedef typename __istream_type::ios_base __ios_base; 183233294Sstas 184127808Snectar const typename __ios_base::fmtflags __flags = __is.flags(); 185127808Snectar __is.flags(__ios_base::dec); 186127808Snectar 187127808Snectar __is >> __lcr._M_x; 188127808Snectar 189233294Sstas __is.flags(__flags); 190178825Sdfr return __is; 19157416Smarkm } 19272445Sassar 193178825Sdfr 194127808Snectar template<class _UIntType, int __w, int __n, int __m, int __r, 195127808Snectar _UIntType __a, int __u, int __s, 196233294Sstas _UIntType __b, int __t, _UIntType __c, int __l> 197233294Sstas const int 198127808Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 199127808Snectar __b, __t, __c, __l>::word_size; 200233294Sstas 201178825Sdfr template<class _UIntType, int __w, int __n, int __m, int __r, 202127808Snectar _UIntType __a, int __u, int __s, 203127808Snectar _UIntType __b, int __t, _UIntType __c, int __l> 204127808Snectar const int 20590926Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 206233294Sstas __b, __t, __c, __l>::state_size; 207127808Snectar 208178825Sdfr template<class _UIntType, int __w, int __n, int __m, int __r, 20957416Smarkm _UIntType __a, int __u, int __s, 210102644Snectar _UIntType __b, int __t, _UIntType __c, int __l> 211102644Snectar const int 212178825Sdfr mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 213127808Snectar __b, __t, __c, __l>::shift_size; 214127808Snectar 21557416Smarkm template<class _UIntType, int __w, int __n, int __m, int __r, 21657416Smarkm _UIntType __a, int __u, int __s, 21790926Snectar _UIntType __b, int __t, _UIntType __c, int __l> 218127808Snectar const int 219127808Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 220127808Snectar __b, __t, __c, __l>::mask_bits; 221127808Snectar 222127808Snectar template<class _UIntType, int __w, int __n, int __m, int __r, 22390926Snectar _UIntType __a, int __u, int __s, 22490926Snectar _UIntType __b, int __t, _UIntType __c, int __l> 22590926Snectar const _UIntType 226127808Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 227127808Snectar __b, __t, __c, __l>::parameter_a; 228127808Snectar 229127808Snectar template<class _UIntType, int __w, int __n, int __m, int __r, 230233294Sstas _UIntType __a, int __u, int __s, 231127808Snectar _UIntType __b, int __t, _UIntType __c, int __l> 232127808Snectar const int 233233294Sstas mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 234178825Sdfr __b, __t, __c, __l>::output_u; 235127808Snectar 236127808Snectar template<class _UIntType, int __w, int __n, int __m, int __r, 237127808Snectar _UIntType __a, int __u, int __s, 238127808Snectar _UIntType __b, int __t, _UIntType __c, int __l> 239127808Snectar const int 240127808Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 241127808Snectar __b, __t, __c, __l>::output_s; 242127808Snectar 243178825Sdfr template<class _UIntType, int __w, int __n, int __m, int __r, 244178825Sdfr _UIntType __a, int __u, int __s, 245178825Sdfr _UIntType __b, int __t, _UIntType __c, int __l> 246178825Sdfr const _UIntType 247127808Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 248127808Snectar __b, __t, __c, __l>::output_b; 24957416Smarkm 250127808Snectar template<class _UIntType, int __w, int __n, int __m, int __r, 251233294Sstas _UIntType __a, int __u, int __s, 252233294Sstas _UIntType __b, int __t, _UIntType __c, int __l> 253127808Snectar const int 254127808Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 255127808Snectar __b, __t, __c, __l>::output_t; 256127808Snectar 257127808Snectar template<class _UIntType, int __w, int __n, int __m, int __r, 25857416Smarkm _UIntType __a, int __u, int __s, 259127808Snectar _UIntType __b, int __t, _UIntType __c, int __l> 260127808Snectar const _UIntType 261178825Sdfr mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 262127808Snectar __b, __t, __c, __l>::output_c; 263127808Snectar 26457416Smarkm template<class _UIntType, int __w, int __n, int __m, int __r, 26557416Smarkm _UIntType __a, int __u, int __s, 266127808Snectar _UIntType __b, int __t, _UIntType __c, int __l> 267127808Snectar const int 268233294Sstas mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 269127808Snectar __b, __t, __c, __l>::output_l; 270127808Snectar 271233294Sstas template<class _UIntType, int __w, int __n, int __m, int __r, 27257416Smarkm _UIntType __a, int __u, int __s, 27357416Smarkm _UIntType __b, int __t, _UIntType __c, int __l> 274120945Snectar void 275127808Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 276233294Sstas __b, __t, __c, __l>:: 277178825Sdfr seed(unsigned long __value) 278233294Sstas { 279233294Sstas _M_x[0] = __detail::__mod<_UIntType, 1, 0, 280233294Sstas __detail::_Shift<_UIntType, __w>::__value>(__value); 28157416Smarkm 282233294Sstas for (int __i = 1; __i < state_size; ++__i) 283127808Snectar { 284233294Sstas _UIntType __x = _M_x[__i - 1]; 285233294Sstas __x ^= __x >> (__w - 2); 28657416Smarkm __x *= 1812433253ul; 287127808Snectar __x += __i; 288127808Snectar _M_x[__i] = __detail::__mod<_UIntType, 1, 0, 289127808Snectar __detail::_Shift<_UIntType, __w>::__value>(__x); 290127808Snectar } 291233294Sstas _M_p = state_size; 292127808Snectar } 293127808Snectar 294233294Sstas template<class _UIntType, int __w, int __n, int __m, int __r, 295233294Sstas _UIntType __a, int __u, int __s, 296233294Sstas _UIntType __b, int __t, _UIntType __c, int __l> 297233294Sstas template<class _Gen> 29857416Smarkm void 299233294Sstas mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 300127808Snectar __b, __t, __c, __l>:: 301127808Snectar seed(_Gen& __gen, false_type) 302233294Sstas { 303233294Sstas for (int __i = 0; __i < state_size; ++__i) 304102644Snectar _M_x[__i] = __detail::__mod<_UIntType, 1, 0, 30557416Smarkm __detail::_Shift<_UIntType, __w>::__value>(__gen()); 306178825Sdfr _M_p = state_size; 30757416Smarkm } 30857416Smarkm 30957416Smarkm template<class _UIntType, int __w, int __n, int __m, int __r, 310178825Sdfr _UIntType __a, int __u, int __s, 31190926Snectar _UIntType __b, int __t, _UIntType __c, int __l> 31290926Snectar typename 31390926Snectar mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 31490926Snectar __b, __t, __c, __l>::result_type 31557416Smarkm mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 316178825Sdfr __b, __t, __c, __l>:: 317178825Sdfr operator()() 318178825Sdfr { 319178825Sdfr // Reload the vector - cost is O(n) amortized over n calls. 320178825Sdfr if (_M_p >= state_size) 321233294Sstas { 322127808Snectar const _UIntType __upper_mask = (~_UIntType()) << __r; 323233294Sstas const _UIntType __lower_mask = ~__upper_mask; 324233294Sstas 325127808Snectar for (int __k = 0; __k < (__n - __m); ++__k) 326233294Sstas { 327178825Sdfr _UIntType __y = ((_M_x[__k] & __upper_mask) 328178825Sdfr | (_M_x[__k + 1] & __lower_mask)); 329127808Snectar _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) 330127808Snectar ^ ((__y & 0x01) ? __a : 0)); 331127808Snectar } 332127808Snectar 333127808Snectar for (int __k = (__n - __m); __k < (__n - 1); ++__k) 334127808Snectar { 335178825Sdfr _UIntType __y = ((_M_x[__k] & __upper_mask) 336127808Snectar | (_M_x[__k + 1] & __lower_mask)); 337178825Sdfr _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) 338178825Sdfr ^ ((__y & 0x01) ? __a : 0)); 339102644Snectar } 340102644Snectar 341102644Snectar _UIntType __y = ((_M_x[__n - 1] & __upper_mask) 342178825Sdfr | (_M_x[0] & __lower_mask)); 343127808Snectar _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) 344127808Snectar ^ ((__y & 0x01) ? __a : 0)); 345127808Snectar _M_p = 0; 346127808Snectar } 347127808Snectar 348127808Snectar // Calculate o(x(i)). 349178825Sdfr result_type __z = _M_x[_M_p++]; 350127808Snectar __z ^= (__z >> __u); 351127808Snectar __z ^= (__z << __s) & __b; 35272445Sassar __z ^= (__z << __t) & __c; 353127808Snectar __z ^= (__z >> __l); 354127808Snectar 355178825Sdfr return __z; 356127808Snectar } 357127808Snectar 358142403Snectar template<class _UIntType, int __w, int __n, int __m, int __r, 359127808Snectar _UIntType __a, int __u, int __s, _UIntType __b, int __t, 360178825Sdfr _UIntType __c, int __l, 361127808Snectar typename _CharT, typename _Traits> 362127808Snectar std::basic_ostream<_CharT, _Traits>& 363178825Sdfr operator<<(std::basic_ostream<_CharT, _Traits>& __os, 364127808Snectar const mersenne_twister<_UIntType, __w, __n, __m, 365127808Snectar __r, __a, __u, __s, __b, __t, __c, __l>& __x) 366178825Sdfr { 367233294Sstas typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 368127808Snectar typedef typename __ostream_type::ios_base __ios_base; 369127808Snectar 370233294Sstas const typename __ios_base::fmtflags __flags = __os.flags(); 371178825Sdfr const _CharT __fill = __os.fill(); 372178825Sdfr const _CharT __space = __os.widen(' '); 373233294Sstas __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 374233294Sstas __os.fill(__space); 375233294Sstas 376102644Snectar for (int __i = 0; __i < __n - 1; ++__i) 37790926Snectar __os << __x._M_x[__i] << __space; 37872445Sassar __os << __x._M_x[__n - 1]; 37957416Smarkm 380233294Sstas __os.flags(__flags); 38157416Smarkm __os.fill(__fill); 38257416Smarkm return __os; 38357416Smarkm } 38457416Smarkm 38557416Smarkm template<class _UIntType, int __w, int __n, int __m, int __r, 38657416Smarkm _UIntType __a, int __u, int __s, _UIntType __b, int __t, 387233294Sstas _UIntType __c, int __l, 38857416Smarkm typename _CharT, typename _Traits> 389120945Snectar std::basic_istream<_CharT, _Traits>& 39090926Snectar operator>>(std::basic_istream<_CharT, _Traits>& __is, 39172445Sassar mersenne_twister<_UIntType, __w, __n, __m, 39257416Smarkm __r, __a, __u, __s, __b, __t, __c, __l>& __x) 39390926Snectar { 394233294Sstas typedef std::basic_istream<_CharT, _Traits> __istream_type; 39590926Snectar typedef typename __istream_type::ios_base __ios_base; 39657416Smarkm 39772445Sassar const typename __ios_base::fmtflags __flags = __is.flags(); 39872445Sassar __is.flags(__ios_base::dec | __ios_base::skipws); 39957416Smarkm 40057416Smarkm for (int __i = 0; __i < __n; ++__i) 40172445Sassar __is >> __x._M_x[__i]; 40272445Sassar 40372445Sassar __is.flags(__flags); 404178825Sdfr return __is; 40572445Sassar } 40672445Sassar 40790926Snectar 40890926Snectar template<typename _IntType, _IntType __m, int __s, int __r> 40978527Sassar const _IntType 41072445Sassar subtract_with_carry<_IntType, __m, __s, __r>::modulus; 41157416Smarkm 412233294Sstas template<typename _IntType, _IntType __m, int __s, int __r> 41390926Snectar const int 41457416Smarkm subtract_with_carry<_IntType, __m, __s, __r>::long_lag; 41557416Smarkm 416233294Sstas template<typename _IntType, _IntType __m, int __s, int __r> 417142403Snectar const int 418142403Snectar subtract_with_carry<_IntType, __m, __s, __r>::short_lag; 419142403Snectar 420142403Snectar template<typename _IntType, _IntType __m, int __s, int __r> 421233294Sstas void 422233294Sstas subtract_with_carry<_IntType, __m, __s, __r>:: 423142403Snectar seed(unsigned long __value) 424142403Snectar { 425142403Snectar if (__value == 0) 426233294Sstas __value = 19780503; 427233294Sstas 428233294Sstas std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563> 429142403Snectar __lcg(__value); 430142403Snectar 431142403Snectar for (int __i = 0; __i < long_lag; ++__i) 432142403Snectar _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg()); 433142403Snectar 434142403Snectar _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 435142403Snectar _M_p = 0; 436142403Snectar } 437142403Snectar 438142403Snectar template<typename _IntType, _IntType __m, int __s, int __r> 439142403Snectar template<class _Gen> 440142403Snectar void 441142403Snectar subtract_with_carry<_IntType, __m, __s, __r>:: 442142403Snectar seed(_Gen& __gen, false_type) 443142403Snectar { 444142403Snectar const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32; 445142403Snectar 446233294Sstas for (int __i = 0; __i < long_lag; ++__i) 44757416Smarkm { 44857416Smarkm _UIntType __tmp = 0; 449178825Sdfr _UIntType __factor = 1; 450233294Sstas for (int __j = 0; __j < __n; ++__j) 451233294Sstas { 452233294Sstas __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0> 453233294Sstas (__gen()) * __factor; 454233294Sstas __factor *= __detail::_Shift<_UIntType, 32>::__value; 455233294Sstas } 456233294Sstas _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp); 457233294Sstas } 458233294Sstas _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 459233294Sstas _M_p = 0; 460233294Sstas } 461233294Sstas 462233294Sstas template<typename _IntType, _IntType __m, int __s, int __r> 463233294Sstas typename subtract_with_carry<_IntType, __m, __s, __r>::result_type 464233294Sstas subtract_with_carry<_IntType, __m, __s, __r>:: 465233294Sstas operator()() 466233294Sstas { 467233294Sstas // Derive short lag index from current index. 468233294Sstas int __ps = _M_p - short_lag; 469233294Sstas if (__ps < 0) 470233294Sstas __ps += long_lag; 47157416Smarkm 47257416Smarkm // Calculate new x(i) without overflow or division. 47357416Smarkm // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry 474233294Sstas // cannot overflow. 475233294Sstas _UIntType __xi; 476233294Sstas if (_M_x[__ps] >= _M_x[_M_p] + _M_carry) 477233294Sstas { 478233294Sstas __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry; 479233294Sstas _M_carry = 0; 480233294Sstas } 48157416Smarkm else 48290926Snectar { 483233294Sstas __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps]; 484233294Sstas _M_carry = 1; 485233294Sstas } 486233294Sstas _M_x[_M_p] = __xi; 487233294Sstas 488233294Sstas // Adjust current index to loop around in ring buffer. 489233294Sstas if (++_M_p >= long_lag) 49090926Snectar _M_p = 0; 49190926Snectar 492178825Sdfr return __xi; 49390926Snectar } 49457416Smarkm 495142403Snectar template<typename _IntType, _IntType __m, int __s, int __r, 49657416Smarkm typename _CharT, typename _Traits> 49757416Smarkm std::basic_ostream<_CharT, _Traits>& 49857416Smarkm operator<<(std::basic_ostream<_CharT, _Traits>& __os, 49957416Smarkm const subtract_with_carry<_IntType, __m, __s, __r>& __x) 500233294Sstas { 501233294Sstas typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 502233294Sstas typedef typename __ostream_type::ios_base __ios_base; 503233294Sstas 504233294Sstas const typename __ios_base::fmtflags __flags = __os.flags(); 505233294Sstas const _CharT __fill = __os.fill(); 506233294Sstas const _CharT __space = __os.widen(' '); 507233294Sstas __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 508233294Sstas __os.fill(__space); 50990926Snectar 510233294Sstas for (int __i = 0; __i < __r; ++__i) 511233294Sstas __os << __x._M_x[__i] << __space; 512233294Sstas __os << __x._M_carry; 513233294Sstas 514233294Sstas __os.flags(__flags); 51557416Smarkm __os.fill(__fill); 51672445Sassar return __os; 517233294Sstas } 518233294Sstas 519233294Sstas template<typename _IntType, _IntType __m, int __s, int __r, 520233294Sstas typename _CharT, typename _Traits> 521233294Sstas std::basic_istream<_CharT, _Traits>& 52290926Snectar operator>>(std::basic_istream<_CharT, _Traits>& __is, 52372445Sassar subtract_with_carry<_IntType, __m, __s, __r>& __x) 524233294Sstas { 525233294Sstas typedef std::basic_ostream<_CharT, _Traits> __istream_type; 526233294Sstas typedef typename __istream_type::ios_base __ios_base; 527233294Sstas 528233294Sstas const typename __ios_base::fmtflags __flags = __is.flags(); 529102644Snectar __is.flags(__ios_base::dec | __ios_base::skipws); 530102644Snectar 531102644Snectar for (int __i = 0; __i < __r; ++__i) 532102644Snectar __is >> __x._M_x[__i]; 533102644Snectar __is >> __x._M_carry; 534102644Snectar 535233294Sstas __is.flags(__flags); 53690926Snectar return __is; 537178825Sdfr } 538233294Sstas 539233294Sstas 540233294Sstas template<typename _RealType, int __w, int __s, int __r> 541233294Sstas const int 542233294Sstas subtract_with_carry_01<_RealType, __w, __s, __r>::word_size; 543233294Sstas 544233294Sstas template<typename _RealType, int __w, int __s, int __r> 545233294Sstas const int 546233294Sstas subtract_with_carry_01<_RealType, __w, __s, __r>::long_lag; 547233294Sstas 548233294Sstas template<typename _RealType, int __w, int __s, int __r> 549233294Sstas const int 550233294Sstas subtract_with_carry_01<_RealType, __w, __s, __r>::short_lag; 551233294Sstas 552233294Sstas template<typename _RealType, int __w, int __s, int __r> 553233294Sstas void 55457416Smarkm subtract_with_carry_01<_RealType, __w, __s, __r>:: 555233294Sstas _M_initialize_npows() 556233294Sstas { 557233294Sstas for (int __j = 0; __j < __n; ++__j) 558233294Sstas#if _GLIBCXX_USE_C99_MATH_TR1 559233294Sstas _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32); 560233294Sstas#else 561233294Sstas _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32); 56257416Smarkm#endif 56390926Snectar } 564233294Sstas 565233294Sstas template<typename _RealType, int __w, int __s, int __r> 566233294Sstas void 567233294Sstas subtract_with_carry_01<_RealType, __w, __s, __r>:: 568233294Sstas seed(unsigned long __value) 569233294Sstas { 570233294Sstas if (__value == 0) 571233294Sstas __value = 19780503; 572233294Sstas 57357416Smarkm // _GLIBCXX_RESOLVE_LIB_DEFECTS 57472445Sassar // 512. Seeding subtract_with_carry_01 from a single unsigned long. 575102644Snectar std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563> 57672445Sassar __lcg(__value); 57772445Sassar 57872445Sassar this->seed(__lcg); 579233294Sstas } 580233294Sstas 581102644Snectar template<typename _RealType, int __w, int __s, int __r> 582142403Snectar template<class _Gen> 58357416Smarkm void 58472445Sassar subtract_with_carry_01<_RealType, __w, __s, __r>:: 58572445Sassar seed(_Gen& __gen, false_type) 586233294Sstas { 58757416Smarkm for (int __i = 0; __i < long_lag; ++__i) 588102644Snectar { 58972445Sassar for (int __j = 0; __j < __n - 1; ++__j) 59072445Sassar _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen()); 59172445Sassar _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, 592233294Sstas __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen()); 593233294Sstas } 594233294Sstas 595233294Sstas _M_carry = 1; 596178825Sdfr for (int __j = 0; __j < __n; ++__j) 597233294Sstas if (_M_x[long_lag - 1][__j] != 0) 598233294Sstas { 599233294Sstas _M_carry = 0; 600233294Sstas break; 601233294Sstas } 602233294Sstas 603233294Sstas _M_p = 0; 604178825Sdfr } 605127808Snectar 606127808Snectar template<typename _RealType, int __w, int __s, int __r> 607127808Snectar typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type 608127808Snectar subtract_with_carry_01<_RealType, __w, __s, __r>:: 609127808Snectar operator()() 610127808Snectar { 611127808Snectar // Derive short lag index from current index. 612233294Sstas int __ps = _M_p - short_lag; 613233294Sstas if (__ps < 0) 614233294Sstas __ps += long_lag; 615127808Snectar 616233294Sstas _UInt32Type __new_carry; 617127808Snectar for (int __j = 0; __j < __n - 1; ++__j) 61878527Sassar { 619102644Snectar if (_M_x[__ps][__j] > _M_x[_M_p][__j] 620233294Sstas || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0)) 621233294Sstas __new_carry = 0; 62278527Sassar else 62357416Smarkm __new_carry = 1; 624127808Snectar 62557416Smarkm _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry; 62657416Smarkm _M_carry = __new_carry; 627233294Sstas } 628233294Sstas 629233294Sstas if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1] 630233294Sstas || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0)) 631233294Sstas __new_carry = 0; 632233294Sstas else 633233294Sstas __new_carry = 1; 634233294Sstas 635233294Sstas _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, 636233294Sstas __detail::_Shift<_UInt32Type, __w % 32>::__value> 637233294Sstas (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry); 638233294Sstas _M_carry = __new_carry; 639233294Sstas 640178825Sdfr result_type __ret = 0.0; 641178825Sdfr for (int __j = 0; __j < __n; ++__j) 642178825Sdfr __ret += _M_x[_M_p][__j] * _M_npows[__j]; 643178825Sdfr 644178825Sdfr // Adjust current index to loop around in ring buffer. 645178825Sdfr if (++_M_p >= long_lag) 646178825Sdfr _M_p = 0; 647178825Sdfr 648178825Sdfr return __ret; 649178825Sdfr } 650178825Sdfr 651178825Sdfr template<typename _RealType, int __w, int __s, int __r, 652102644Snectar typename _CharT, typename _Traits> 65357416Smarkm std::basic_ostream<_CharT, _Traits>& 654178825Sdfr operator<<(std::basic_ostream<_CharT, _Traits>& __os, 655233294Sstas const subtract_with_carry_01<_RealType, __w, __s, __r>& __x) 656233294Sstas { 657233294Sstas typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 658102644Snectar typedef typename __ostream_type::ios_base __ios_base; 659233294Sstas 660233294Sstas const typename __ios_base::fmtflags __flags = __os.flags(); 661102644Snectar const _CharT __fill = __os.fill(); 662233294Sstas const _CharT __space = __os.widen(' '); 66357416Smarkm __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 664233294Sstas __os.fill(__space); 665233294Sstas 66672445Sassar for (int __i = 0; __i < __r; ++__i) 66757416Smarkm for (int __j = 0; __j < __x.__n; ++__j) 66857416Smarkm __os << __x._M_x[__i][__j] << __space; 66990926Snectar __os << __x._M_carry; 670127808Snectar 67190926Snectar __os.flags(__flags); 67257416Smarkm __os.fill(__fill); 67357416Smarkm return __os; 67457416Smarkm } 67590926Snectar 67690926Snectar template<typename _RealType, int __w, int __s, int __r, 677142403Snectar typename _CharT, typename _Traits> 678178825Sdfr std::basic_istream<_CharT, _Traits>& 679142403Snectar operator>>(std::basic_istream<_CharT, _Traits>& __is, 68090926Snectar subtract_with_carry_01<_RealType, __w, __s, __r>& __x) 68157416Smarkm { 68257416Smarkm typedef std::basic_istream<_CharT, _Traits> __istream_type; 68390926Snectar typedef typename __istream_type::ios_base __ios_base; 68457416Smarkm 68557416Smarkm const typename __ios_base::fmtflags __flags = __is.flags(); 68657416Smarkm __is.flags(__ios_base::dec | __ios_base::skipws); 68790926Snectar 68890926Snectar for (int __i = 0; __i < __r; ++__i) 68957416Smarkm for (int __j = 0; __j < __x.__n; ++__j) 69090926Snectar __is >> __x._M_x[__i][__j]; 691127808Snectar __is >> __x._M_carry; 69290926Snectar 69390926Snectar __is.flags(__flags); 69457416Smarkm return __is; 69557416Smarkm } 69657416Smarkm 69757416Smarkm template<class _UniformRandomNumberGenerator, int __p, int __r> 69857416Smarkm const int 699178825Sdfr discard_block<_UniformRandomNumberGenerator, __p, __r>::block_size; 700233294Sstas 70157416Smarkm template<class _UniformRandomNumberGenerator, int __p, int __r> 70257416Smarkm const int 70390926Snectar discard_block<_UniformRandomNumberGenerator, __p, __r>::used_block; 70490926Snectar 70590926Snectar template<class _UniformRandomNumberGenerator, int __p, int __r> 70657416Smarkm typename discard_block<_UniformRandomNumberGenerator, 70790926Snectar __p, __r>::result_type 70890926Snectar discard_block<_UniformRandomNumberGenerator, __p, __r>:: 70957416Smarkm operator()() 71090926Snectar { 711233294Sstas if (_M_n >= used_block) 712127808Snectar { 71390926Snectar while (_M_n < block_size) 714178825Sdfr { 71557416Smarkm _M_b(); 71690926Snectar ++_M_n; 71757416Smarkm } 71890926Snectar _M_n = 0; 71957416Smarkm } 720142403Snectar ++_M_n; 721142403Snectar return _M_b(); 722233294Sstas } 723233294Sstas 72490926Snectar template<class _UniformRandomNumberGenerator, int __p, int __r, 72557416Smarkm typename _CharT, typename _Traits> 72690926Snectar std::basic_ostream<_CharT, _Traits>& 72790926Snectar operator<<(std::basic_ostream<_CharT, _Traits>& __os, 728120945Snectar const discard_block<_UniformRandomNumberGenerator, 729120945Snectar __p, __r>& __x) 730120945Snectar { 731178825Sdfr typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 732178825Sdfr typedef typename __ostream_type::ios_base __ios_base; 733233294Sstas 734233294Sstas const typename __ios_base::fmtflags __flags = __os.flags(); 73590926Snectar const _CharT __fill = __os.fill(); 73690926Snectar const _CharT __space = __os.widen(' '); 73790926Snectar __os.flags(__ios_base::dec | __ios_base::fixed 738178825Sdfr | __ios_base::left); 739178825Sdfr __os.fill(__space); 740233294Sstas 741233294Sstas __os << __x._M_b << __space << __x._M_n; 74290926Snectar 74390926Snectar __os.flags(__flags); 744233294Sstas __os.fill(__fill); 745233294Sstas return __os; 74690926Snectar } 74790926Snectar 748178825Sdfr template<class _UniformRandomNumberGenerator, int __p, int __r, 749178825Sdfr typename _CharT, typename _Traits> 750233294Sstas std::basic_istream<_CharT, _Traits>& 751233294Sstas operator>>(std::basic_istream<_CharT, _Traits>& __is, 752178825Sdfr discard_block<_UniformRandomNumberGenerator, __p, __r>& __x) 753178825Sdfr { 754233294Sstas typedef std::basic_istream<_CharT, _Traits> __istream_type; 755233294Sstas typedef typename __istream_type::ios_base __ios_base; 75690926Snectar 75790926Snectar const typename __ios_base::fmtflags __flags = __is.flags(); 75857416Smarkm __is.flags(__ios_base::dec | __ios_base::skipws); 759233294Sstas 760127808Snectar __is >> __x._M_b >> __x._M_n; 76190926Snectar 76257416Smarkm __is.flags(__flags); 76390926Snectar return __is; 76457416Smarkm } 76590926Snectar 76690926Snectar 76790926Snectar template<class _UniformRandomNumberGenerator1, int __s1, 768127808Snectar class _UniformRandomNumberGenerator2, int __s2> 769127808Snectar const int 770127808Snectar xor_combine<_UniformRandomNumberGenerator1, __s1, 771127808Snectar _UniformRandomNumberGenerator2, __s2>::shift1; 772127808Snectar 773127808Snectar template<class _UniformRandomNumberGenerator1, int __s1, 774127808Snectar class _UniformRandomNumberGenerator2, int __s2> 775127808Snectar const int 776178825Sdfr xor_combine<_UniformRandomNumberGenerator1, __s1, 777178825Sdfr _UniformRandomNumberGenerator2, __s2>::shift2; 778178825Sdfr 77990926Snectar template<class _UniformRandomNumberGenerator1, int __s1, 78090926Snectar class _UniformRandomNumberGenerator2, int __s2> 781233294Sstas void 782233294Sstas xor_combine<_UniformRandomNumberGenerator1, __s1, 783178825Sdfr _UniformRandomNumberGenerator2, __s2>:: 784127808Snectar _M_initialize_max() 785127808Snectar { 786178825Sdfr const int __w = std::numeric_limits<result_type>::digits; 787142403Snectar 788142403Snectar const result_type __m1 = 789178825Sdfr std::min(result_type(_M_b1.max() - _M_b1.min()), 790178825Sdfr __detail::_Shift<result_type, __w - __s1>::__value - 1); 791178825Sdfr 792178825Sdfr const result_type __m2 = 793178825Sdfr std::min(result_type(_M_b2.max() - _M_b2.min()), 794178825Sdfr __detail::_Shift<result_type, __w - __s2>::__value - 1); 795178825Sdfr 796178825Sdfr // NB: In TR1 s1 is not required to be >= s2. 797178825Sdfr if (__s1 < __s2) 798178825Sdfr _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1; 79990926Snectar else 80090926Snectar _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2; 80157416Smarkm } 80257416Smarkm 80357416Smarkm template<class _UniformRandomNumberGenerator1, int __s1, 80457416Smarkm class _UniformRandomNumberGenerator2, int __s2> 80557416Smarkm typename xor_combine<_UniformRandomNumberGenerator1, __s1, 80672445Sassar _UniformRandomNumberGenerator2, __s2>::result_type 80772445Sassar xor_combine<_UniformRandomNumberGenerator1, __s1, 80872445Sassar _UniformRandomNumberGenerator2, __s2>:: 80972445Sassar _M_initialize_max_aux(result_type __a, result_type __b, int __d) 81057416Smarkm { 81157416Smarkm const result_type __two2d = result_type(1) << __d; 81257416Smarkm const result_type __c = __a * __two2d; 813178825Sdfr 814178825Sdfr if (__a == 0 || __b < __two2d) 81557416Smarkm return __c + __b; 81657416Smarkm 81757416Smarkm const result_type __t = std::max(__c, __b); 81857416Smarkm const result_type __u = std::min(__c, __b); 81957416Smarkm 82057416Smarkm result_type __ub = __u; 82172445Sassar result_type __p; 82272445Sassar for (__p = 0; __ub != 1; __ub >>= 1) 82357416Smarkm ++__p; 824178825Sdfr 825178825Sdfr const result_type __two2p = result_type(1) << __p; 826178825Sdfr const result_type __k = __t / __two2p; 827178825Sdfr 828178825Sdfr if (__k & 1) 829178825Sdfr return (__k + 1) * __two2p - 1; 830178825Sdfr 831178825Sdfr if (__c >= __b) 832178825Sdfr return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p) 833178825Sdfr / __two2d, 834178825Sdfr __u % __two2p, __d); 83557416Smarkm else 83657416Smarkm return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p) 83757416Smarkm / __two2d, 838102644Snectar __t % __two2p, __d); 839102644Snectar } 840178825Sdfr 841178825Sdfr template<class _UniformRandomNumberGenerator1, int __s1, 842102644Snectar class _UniformRandomNumberGenerator2, int __s2, 843102644Snectar typename _CharT, typename _Traits> 844102644Snectar std::basic_ostream<_CharT, _Traits>& 845102644Snectar operator<<(std::basic_ostream<_CharT, _Traits>& __os, 846102644Snectar const xor_combine<_UniformRandomNumberGenerator1, __s1, 847102644Snectar _UniformRandomNumberGenerator2, __s2>& __x) 848178825Sdfr { 849102644Snectar typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 850102644Snectar typedef typename __ostream_type::ios_base __ios_base; 851102644Snectar 852102644Snectar const typename __ios_base::fmtflags __flags = __os.flags(); 853102644Snectar const _CharT __fill = __os.fill(); 854102644Snectar const _CharT __space = __os.widen(' '); 855102644Snectar __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 856102644Snectar __os.fill(__space); 857102644Snectar 858102644Snectar __os << __x.base1() << __space << __x.base2(); 859102644Snectar 860102644Snectar __os.flags(__flags); 861102644Snectar __os.fill(__fill); 862102644Snectar return __os; 863102644Snectar } 864178825Sdfr 865102644Snectar template<class _UniformRandomNumberGenerator1, int __s1, 866102644Snectar class _UniformRandomNumberGenerator2, int __s2, 867102644Snectar typename _CharT, typename _Traits> 868102644Snectar std::basic_istream<_CharT, _Traits>& 869233294Sstas operator>>(std::basic_istream<_CharT, _Traits>& __is, 870233294Sstas xor_combine<_UniformRandomNumberGenerator1, __s1, 871233294Sstas _UniformRandomNumberGenerator2, __s2>& __x) 87257416Smarkm { 87357416Smarkm typedef std::basic_istream<_CharT, _Traits> __istream_type; 87457416Smarkm typedef typename __istream_type::ios_base __ios_base; 87557416Smarkm 87657416Smarkm const typename __ios_base::fmtflags __flags = __is.flags(); 87757416Smarkm __is.flags(__ios_base::skipws); 87857416Smarkm 87957416Smarkm __is >> __x._M_b1 >> __x._M_b2; 88057416Smarkm 88157416Smarkm __is.flags(__flags); 88257416Smarkm return __is; 88357416Smarkm } 88457416Smarkm 88557416Smarkm 88657416Smarkm template<typename _IntType> 88757416Smarkm template<typename _UniformRandomNumberGenerator> 88857416Smarkm typename uniform_int<_IntType>::result_type 88957416Smarkm uniform_int<_IntType>:: 89057416Smarkm _M_call(_UniformRandomNumberGenerator& __urng, 89157416Smarkm result_type __min, result_type __max, true_type) 89257416Smarkm { 89357416Smarkm // XXX Must be fixed to work well for *arbitrary* __urng.max(), 89457416Smarkm // __urng.min(), __max, __min. Currently works fine only in the 89557416Smarkm // most common case __urng.max() - __urng.min() >= __max - __min, 89657416Smarkm // with __urng.max() > __urng.min() >= 0. 89757416Smarkm typedef typename __gnu_cxx::__add_unsigned<typename 89857416Smarkm _UniformRandomNumberGenerator::result_type>::__type __urntype; 89957416Smarkm typedef typename __gnu_cxx::__add_unsigned<result_type>::__type 90057416Smarkm __utype; 90157416Smarkm typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype) 90257416Smarkm > sizeof(__utype)), 90357416Smarkm __urntype, __utype>::__type __uctype; 90457416Smarkm 90557416Smarkm result_type __ret; 90657416Smarkm 90757416Smarkm const __urntype __urnmin = __urng.min(); 90857416Smarkm const __urntype __urnmax = __urng.max(); 90957416Smarkm const __urntype __urnrange = __urnmax - __urnmin; 91057416Smarkm const __uctype __urange = __max - __min; 91157416Smarkm const __uctype __udenom = (__urnrange <= __urange 91257416Smarkm ? 1 : __urnrange / (__urange + 1)); 91357416Smarkm do 91457416Smarkm __ret = (__urntype(__urng()) - __urnmin) / __udenom; 91557416Smarkm while (__ret > __max - __min); 91657416Smarkm 91757416Smarkm return __ret + __min; 91857416Smarkm } 91957416Smarkm 92057416Smarkm template<typename _IntType, typename _CharT, typename _Traits> 92157416Smarkm std::basic_ostream<_CharT, _Traits>& 92257416Smarkm operator<<(std::basic_ostream<_CharT, _Traits>& __os, 92357416Smarkm const uniform_int<_IntType>& __x) 92457416Smarkm { 92557416Smarkm typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 92657416Smarkm typedef typename __ostream_type::ios_base __ios_base; 92757416Smarkm 92857416Smarkm const typename __ios_base::fmtflags __flags = __os.flags(); 92957416Smarkm const _CharT __fill = __os.fill(); 93057416Smarkm const _CharT __space = __os.widen(' '); 93157416Smarkm __os.flags(__ios_base::scientific | __ios_base::left); 93257416Smarkm __os.fill(__space); 93357416Smarkm 93457416Smarkm __os << __x.min() << __space << __x.max(); 93557416Smarkm 93672445Sassar __os.flags(__flags); 937178825Sdfr __os.fill(__fill); 93857416Smarkm return __os; 939178825Sdfr } 940178825Sdfr 941178825Sdfr template<typename _IntType, typename _CharT, typename _Traits> 942120945Snectar std::basic_istream<_CharT, _Traits>& 943178825Sdfr operator>>(std::basic_istream<_CharT, _Traits>& __is, 94457416Smarkm uniform_int<_IntType>& __x) 94557416Smarkm { 94657416Smarkm typedef std::basic_istream<_CharT, _Traits> __istream_type; 94757416Smarkm typedef typename __istream_type::ios_base __ios_base; 94857416Smarkm 949178825Sdfr const typename __ios_base::fmtflags __flags = __is.flags(); 950178825Sdfr __is.flags(__ios_base::dec | __ios_base::skipws); 951178825Sdfr 952178825Sdfr __is >> __x._M_min >> __x._M_max; 953178825Sdfr 954178825Sdfr __is.flags(__flags); 955178825Sdfr return __is; 956178825Sdfr } 957233294Sstas 958178825Sdfr 959178825Sdfr template<typename _CharT, typename _Traits> 960178825Sdfr std::basic_ostream<_CharT, _Traits>& 961178825Sdfr operator<<(std::basic_ostream<_CharT, _Traits>& __os, 962178825Sdfr const bernoulli_distribution& __x) 963178825Sdfr { 964178825Sdfr typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 965178825Sdfr typedef typename __ostream_type::ios_base __ios_base; 966178825Sdfr 967178825Sdfr const typename __ios_base::fmtflags __flags = __os.flags(); 968178825Sdfr const _CharT __fill = __os.fill(); 969178825Sdfr const std::streamsize __precision = __os.precision(); 970233294Sstas __os.flags(__ios_base::scientific | __ios_base::left); 97157416Smarkm __os.fill(__os.widen(' ')); 97257416Smarkm __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10); 97357416Smarkm 974 __os << __x.p(); 975 976 __os.flags(__flags); 977 __os.fill(__fill); 978 __os.precision(__precision); 979 return __os; 980 } 981 982 983 template<typename _IntType, typename _RealType> 984 template<class _UniformRandomNumberGenerator> 985 typename geometric_distribution<_IntType, _RealType>::result_type 986 geometric_distribution<_IntType, _RealType>:: 987 operator()(_UniformRandomNumberGenerator& __urng) 988 { 989 // About the epsilon thing see this thread: 990 // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html 991 const _RealType __naf = 992 (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 993 // The largest _RealType convertible to _IntType. 994 const _RealType __thr = 995 std::numeric_limits<_IntType>::max() + __naf; 996 997 _RealType __cand; 998 do 999 __cand = std::ceil(std::log(__urng()) / _M_log_p); 1000 while (__cand >= __thr); 1001 1002 return result_type(__cand + __naf); 1003 } 1004 1005 template<typename _IntType, typename _RealType, 1006 typename _CharT, typename _Traits> 1007 std::basic_ostream<_CharT, _Traits>& 1008 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1009 const geometric_distribution<_IntType, _RealType>& __x) 1010 { 1011 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1012 typedef typename __ostream_type::ios_base __ios_base; 1013 1014 const typename __ios_base::fmtflags __flags = __os.flags(); 1015 const _CharT __fill = __os.fill(); 1016 const std::streamsize __precision = __os.precision(); 1017 __os.flags(__ios_base::scientific | __ios_base::left); 1018 __os.fill(__os.widen(' ')); 1019 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1020 1021 __os << __x.p(); 1022 1023 __os.flags(__flags); 1024 __os.fill(__fill); 1025 __os.precision(__precision); 1026 return __os; 1027 } 1028 1029 1030 template<typename _IntType, typename _RealType> 1031 void 1032 poisson_distribution<_IntType, _RealType>:: 1033 _M_initialize() 1034 { 1035#if _GLIBCXX_USE_C99_MATH_TR1 1036 if (_M_mean >= 12) 1037 { 1038 const _RealType __m = std::floor(_M_mean); 1039 _M_lm_thr = std::log(_M_mean); 1040 _M_lfm = std::tr1::lgamma(__m + 1); 1041 _M_sm = std::sqrt(__m); 1042 1043 const _RealType __pi_4 = 0.7853981633974483096156608458198757L; 1044 const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m 1045 / __pi_4)); 1046 _M_d = std::tr1::round(std::max(_RealType(6), 1047 std::min(__m, __dx))); 1048 const _RealType __cx = 2 * __m + _M_d; 1049 _M_scx = std::sqrt(__cx / 2); 1050 _M_1cx = 1 / __cx; 1051 1052 _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); 1053 _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d; 1054 } 1055 else 1056#endif 1057 _M_lm_thr = std::exp(-_M_mean); 1058 } 1059 1060 /** 1061 * A rejection algorithm when mean >= 12 and a simple method based 1062 * upon the multiplication of uniform random variates otherwise. 1063 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 1064 * is defined. 1065 * 1066 * Reference: 1067 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1068 * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). 1069 */ 1070 template<typename _IntType, typename _RealType> 1071 template<class _UniformRandomNumberGenerator> 1072 typename poisson_distribution<_IntType, _RealType>::result_type 1073 poisson_distribution<_IntType, _RealType>:: 1074 operator()(_UniformRandomNumberGenerator& __urng) 1075 { 1076#if _GLIBCXX_USE_C99_MATH_TR1 1077 if (_M_mean >= 12) 1078 { 1079 _RealType __x; 1080 1081 // See comments above... 1082 const _RealType __naf = 1083 (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 1084 const _RealType __thr = 1085 std::numeric_limits<_IntType>::max() + __naf; 1086 1087 const _RealType __m = std::floor(_M_mean); 1088 // sqrt(pi / 2) 1089 const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 1090 const _RealType __c1 = _M_sm * __spi_2; 1091 const _RealType __c2 = _M_c2b + __c1; 1092 const _RealType __c3 = __c2 + 1; 1093 const _RealType __c4 = __c3 + 1; 1094 // e^(1 / 78) 1095 const _RealType __e178 = 1.0129030479320018583185514777512983L; 1096 const _RealType __c5 = __c4 + __e178; 1097 const _RealType __c = _M_cb + __c5; 1098 const _RealType __2cx = 2 * (2 * __m + _M_d); 1099 1100 bool __reject = true; 1101 do 1102 { 1103 const _RealType __u = __c * __urng(); 1104 const _RealType __e = -std::log(__urng()); 1105 1106 _RealType __w = 0.0; 1107 1108 if (__u <= __c1) 1109 { 1110 const _RealType __n = _M_nd(__urng); 1111 const _RealType __y = -std::abs(__n) * _M_sm - 1; 1112 __x = std::floor(__y); 1113 __w = -__n * __n / 2; 1114 if (__x < -__m) 1115 continue; 1116 } 1117 else if (__u <= __c2) 1118 { 1119 const _RealType __n = _M_nd(__urng); 1120 const _RealType __y = 1 + std::abs(__n) * _M_scx; 1121 __x = std::ceil(__y); 1122 __w = __y * (2 - __y) * _M_1cx; 1123 if (__x > _M_d) 1124 continue; 1125 } 1126 else if (__u <= __c3) 1127 // NB: This case not in the book, nor in the Errata, 1128 // but should be ok... 1129 __x = -1; 1130 else if (__u <= __c4) 1131 __x = 0; 1132 else if (__u <= __c5) 1133 __x = 1; 1134 else 1135 { 1136 const _RealType __v = -std::log(__urng()); 1137 const _RealType __y = _M_d + __v * __2cx / _M_d; 1138 __x = std::ceil(__y); 1139 __w = -_M_d * _M_1cx * (1 + __y / 2); 1140 } 1141 1142 __reject = (__w - __e - __x * _M_lm_thr 1143 > _M_lfm - std::tr1::lgamma(__x + __m + 1)); 1144 1145 __reject |= __x + __m >= __thr; 1146 1147 } while (__reject); 1148 1149 return result_type(__x + __m + __naf); 1150 } 1151 else 1152#endif 1153 { 1154 _IntType __x = 0; 1155 _RealType __prod = 1.0; 1156 1157 do 1158 { 1159 __prod *= __urng(); 1160 __x += 1; 1161 } 1162 while (__prod > _M_lm_thr); 1163 1164 return __x - 1; 1165 } 1166 } 1167 1168 template<typename _IntType, typename _RealType, 1169 typename _CharT, typename _Traits> 1170 std::basic_ostream<_CharT, _Traits>& 1171 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1172 const poisson_distribution<_IntType, _RealType>& __x) 1173 { 1174 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1175 typedef typename __ostream_type::ios_base __ios_base; 1176 1177 const typename __ios_base::fmtflags __flags = __os.flags(); 1178 const _CharT __fill = __os.fill(); 1179 const std::streamsize __precision = __os.precision(); 1180 const _CharT __space = __os.widen(' '); 1181 __os.flags(__ios_base::scientific | __ios_base::left); 1182 __os.fill(__space); 1183 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1184 1185 __os << __x.mean() << __space << __x._M_nd; 1186 1187 __os.flags(__flags); 1188 __os.fill(__fill); 1189 __os.precision(__precision); 1190 return __os; 1191 } 1192 1193 template<typename _IntType, typename _RealType, 1194 typename _CharT, typename _Traits> 1195 std::basic_istream<_CharT, _Traits>& 1196 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1197 poisson_distribution<_IntType, _RealType>& __x) 1198 { 1199 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1200 typedef typename __istream_type::ios_base __ios_base; 1201 1202 const typename __ios_base::fmtflags __flags = __is.flags(); 1203 __is.flags(__ios_base::skipws); 1204 1205 __is >> __x._M_mean >> __x._M_nd; 1206 __x._M_initialize(); 1207 1208 __is.flags(__flags); 1209 return __is; 1210 } 1211 1212 1213 template<typename _IntType, typename _RealType> 1214 void 1215 binomial_distribution<_IntType, _RealType>:: 1216 _M_initialize() 1217 { 1218 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 1219 1220 _M_easy = true; 1221 1222#if _GLIBCXX_USE_C99_MATH_TR1 1223 if (_M_t * __p12 >= 8) 1224 { 1225 _M_easy = false; 1226 const _RealType __np = std::floor(_M_t * __p12); 1227 const _RealType __pa = __np / _M_t; 1228 const _RealType __1p = 1 - __pa; 1229 1230 const _RealType __pi_4 = 0.7853981633974483096156608458198757L; 1231 const _RealType __d1x = 1232 std::sqrt(__np * __1p * std::log(32 * __np 1233 / (81 * __pi_4 * __1p))); 1234 _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x)); 1235 const _RealType __d2x = 1236 std::sqrt(__np * __1p * std::log(32 * _M_t * __1p 1237 / (__pi_4 * __pa))); 1238 _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x)); 1239 1240 // sqrt(pi / 2) 1241 const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 1242 _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); 1243 _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); 1244 _M_c = 2 * _M_d1 / __np; 1245 _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; 1246 const _RealType __a12 = _M_a1 + _M_s2 * __spi_2; 1247 const _RealType __s1s = _M_s1 * _M_s1; 1248 _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) 1249 * 2 * __s1s / _M_d1 1250 * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); 1251 const _RealType __s2s = _M_s2 * _M_s2; 1252 _M_s = (_M_a123 + 2 * __s2s / _M_d2 1253 * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); 1254 _M_lf = (std::tr1::lgamma(__np + 1) 1255 + std::tr1::lgamma(_M_t - __np + 1)); 1256 _M_lp1p = std::log(__pa / __1p); 1257 1258 _M_q = -std::log(1 - (__p12 - __pa) / __1p); 1259 } 1260 else 1261#endif 1262 _M_q = -std::log(1 - __p12); 1263 } 1264 1265 template<typename _IntType, typename _RealType> 1266 template<class _UniformRandomNumberGenerator> 1267 typename binomial_distribution<_IntType, _RealType>::result_type 1268 binomial_distribution<_IntType, _RealType>:: 1269 _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) 1270 { 1271 _IntType __x = 0; 1272 _RealType __sum = 0; 1273 1274 do 1275 { 1276 const _RealType __e = -std::log(__urng()); 1277 __sum += __e / (__t - __x); 1278 __x += 1; 1279 } 1280 while (__sum <= _M_q); 1281 1282 return __x - 1; 1283 } 1284 1285 /** 1286 * A rejection algorithm when t * p >= 8 and a simple waiting time 1287 * method - the second in the referenced book - otherwise. 1288 * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 1289 * is defined. 1290 * 1291 * Reference: 1292 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1293 * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 1294 */ 1295 template<typename _IntType, typename _RealType> 1296 template<class _UniformRandomNumberGenerator> 1297 typename binomial_distribution<_IntType, _RealType>::result_type 1298 binomial_distribution<_IntType, _RealType>:: 1299 operator()(_UniformRandomNumberGenerator& __urng) 1300 { 1301 result_type __ret; 1302 const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 1303 1304#if _GLIBCXX_USE_C99_MATH_TR1 1305 if (!_M_easy) 1306 { 1307 _RealType __x; 1308 1309 // See comments above... 1310 const _RealType __naf = 1311 (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 1312 const _RealType __thr = 1313 std::numeric_limits<_IntType>::max() + __naf; 1314 1315 const _RealType __np = std::floor(_M_t * __p12); 1316 const _RealType __pa = __np / _M_t; 1317 1318 // sqrt(pi / 2) 1319 const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 1320 const _RealType __a1 = _M_a1; 1321 const _RealType __a12 = __a1 + _M_s2 * __spi_2; 1322 const _RealType __a123 = _M_a123; 1323 const _RealType __s1s = _M_s1 * _M_s1; 1324 const _RealType __s2s = _M_s2 * _M_s2; 1325 1326 bool __reject; 1327 do 1328 { 1329 const _RealType __u = _M_s * __urng(); 1330 1331 _RealType __v; 1332 1333 if (__u <= __a1) 1334 { 1335 const _RealType __n = _M_nd(__urng); 1336 const _RealType __y = _M_s1 * std::abs(__n); 1337 __reject = __y >= _M_d1; 1338 if (!__reject) 1339 { 1340 const _RealType __e = -std::log(__urng()); 1341 __x = std::floor(__y); 1342 __v = -__e - __n * __n / 2 + _M_c; 1343 } 1344 } 1345 else if (__u <= __a12) 1346 { 1347 const _RealType __n = _M_nd(__urng); 1348 const _RealType __y = _M_s2 * std::abs(__n); 1349 __reject = __y >= _M_d2; 1350 if (!__reject) 1351 { 1352 const _RealType __e = -std::log(__urng()); 1353 __x = std::floor(-__y); 1354 __v = -__e - __n * __n / 2; 1355 } 1356 } 1357 else if (__u <= __a123) 1358 { 1359 const _RealType __e1 = -std::log(__urng()); 1360 const _RealType __e2 = -std::log(__urng()); 1361 1362 const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1; 1363 __x = std::floor(__y); 1364 __v = (-__e2 + _M_d1 * (1 / (_M_t - __np) 1365 -__y / (2 * __s1s))); 1366 __reject = false; 1367 } 1368 else 1369 { 1370 const _RealType __e1 = -std::log(__urng()); 1371 const _RealType __e2 = -std::log(__urng()); 1372 1373 const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2; 1374 __x = std::floor(-__y); 1375 __v = -__e2 - _M_d2 * __y / (2 * __s2s); 1376 __reject = false; 1377 } 1378 1379 __reject = __reject || __x < -__np || __x > _M_t - __np; 1380 if (!__reject) 1381 { 1382 const _RealType __lfx = 1383 std::tr1::lgamma(__np + __x + 1) 1384 + std::tr1::lgamma(_M_t - (__np + __x) + 1); 1385 __reject = __v > _M_lf - __lfx + __x * _M_lp1p; 1386 } 1387 1388 __reject |= __x + __np >= __thr; 1389 } 1390 while (__reject); 1391 1392 __x += __np + __naf; 1393 1394 const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); 1395 __ret = _IntType(__x) + __z; 1396 } 1397 else 1398#endif 1399 __ret = _M_waiting(__urng, _M_t); 1400 1401 if (__p12 != _M_p) 1402 __ret = _M_t - __ret; 1403 return __ret; 1404 } 1405 1406 template<typename _IntType, typename _RealType, 1407 typename _CharT, typename _Traits> 1408 std::basic_ostream<_CharT, _Traits>& 1409 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1410 const binomial_distribution<_IntType, _RealType>& __x) 1411 { 1412 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1413 typedef typename __ostream_type::ios_base __ios_base; 1414 1415 const typename __ios_base::fmtflags __flags = __os.flags(); 1416 const _CharT __fill = __os.fill(); 1417 const std::streamsize __precision = __os.precision(); 1418 const _CharT __space = __os.widen(' '); 1419 __os.flags(__ios_base::scientific | __ios_base::left); 1420 __os.fill(__space); 1421 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1422 1423 __os << __x.t() << __space << __x.p() 1424 << __space << __x._M_nd; 1425 1426 __os.flags(__flags); 1427 __os.fill(__fill); 1428 __os.precision(__precision); 1429 return __os; 1430 } 1431 1432 template<typename _IntType, typename _RealType, 1433 typename _CharT, typename _Traits> 1434 std::basic_istream<_CharT, _Traits>& 1435 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1436 binomial_distribution<_IntType, _RealType>& __x) 1437 { 1438 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1439 typedef typename __istream_type::ios_base __ios_base; 1440 1441 const typename __ios_base::fmtflags __flags = __is.flags(); 1442 __is.flags(__ios_base::dec | __ios_base::skipws); 1443 1444 __is >> __x._M_t >> __x._M_p >> __x._M_nd; 1445 __x._M_initialize(); 1446 1447 __is.flags(__flags); 1448 return __is; 1449 } 1450 1451 1452 template<typename _RealType, typename _CharT, typename _Traits> 1453 std::basic_ostream<_CharT, _Traits>& 1454 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1455 const uniform_real<_RealType>& __x) 1456 { 1457 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1458 typedef typename __ostream_type::ios_base __ios_base; 1459 1460 const typename __ios_base::fmtflags __flags = __os.flags(); 1461 const _CharT __fill = __os.fill(); 1462 const std::streamsize __precision = __os.precision(); 1463 const _CharT __space = __os.widen(' '); 1464 __os.flags(__ios_base::scientific | __ios_base::left); 1465 __os.fill(__space); 1466 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1467 1468 __os << __x.min() << __space << __x.max(); 1469 1470 __os.flags(__flags); 1471 __os.fill(__fill); 1472 __os.precision(__precision); 1473 return __os; 1474 } 1475 1476 template<typename _RealType, typename _CharT, typename _Traits> 1477 std::basic_istream<_CharT, _Traits>& 1478 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1479 uniform_real<_RealType>& __x) 1480 { 1481 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1482 typedef typename __istream_type::ios_base __ios_base; 1483 1484 const typename __ios_base::fmtflags __flags = __is.flags(); 1485 __is.flags(__ios_base::skipws); 1486 1487 __is >> __x._M_min >> __x._M_max; 1488 1489 __is.flags(__flags); 1490 return __is; 1491 } 1492 1493 1494 template<typename _RealType, typename _CharT, typename _Traits> 1495 std::basic_ostream<_CharT, _Traits>& 1496 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1497 const exponential_distribution<_RealType>& __x) 1498 { 1499 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1500 typedef typename __ostream_type::ios_base __ios_base; 1501 1502 const typename __ios_base::fmtflags __flags = __os.flags(); 1503 const _CharT __fill = __os.fill(); 1504 const std::streamsize __precision = __os.precision(); 1505 __os.flags(__ios_base::scientific | __ios_base::left); 1506 __os.fill(__os.widen(' ')); 1507 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1508 1509 __os << __x.lambda(); 1510 1511 __os.flags(__flags); 1512 __os.fill(__fill); 1513 __os.precision(__precision); 1514 return __os; 1515 } 1516 1517 1518 /** 1519 * Polar method due to Marsaglia. 1520 * 1521 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1522 * New York, 1986, Ch. V, Sect. 4.4. 1523 */ 1524 template<typename _RealType> 1525 template<class _UniformRandomNumberGenerator> 1526 typename normal_distribution<_RealType>::result_type 1527 normal_distribution<_RealType>:: 1528 operator()(_UniformRandomNumberGenerator& __urng) 1529 { 1530 result_type __ret; 1531 1532 if (_M_saved_available) 1533 { 1534 _M_saved_available = false; 1535 __ret = _M_saved; 1536 } 1537 else 1538 { 1539 result_type __x, __y, __r2; 1540 do 1541 { 1542 __x = result_type(2.0) * __urng() - 1.0; 1543 __y = result_type(2.0) * __urng() - 1.0; 1544 __r2 = __x * __x + __y * __y; 1545 } 1546 while (__r2 > 1.0 || __r2 == 0.0); 1547 1548 const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 1549 _M_saved = __x * __mult; 1550 _M_saved_available = true; 1551 __ret = __y * __mult; 1552 } 1553 1554 __ret = __ret * _M_sigma + _M_mean; 1555 return __ret; 1556 } 1557 1558 template<typename _RealType, typename _CharT, typename _Traits> 1559 std::basic_ostream<_CharT, _Traits>& 1560 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1561 const normal_distribution<_RealType>& __x) 1562 { 1563 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1564 typedef typename __ostream_type::ios_base __ios_base; 1565 1566 const typename __ios_base::fmtflags __flags = __os.flags(); 1567 const _CharT __fill = __os.fill(); 1568 const std::streamsize __precision = __os.precision(); 1569 const _CharT __space = __os.widen(' '); 1570 __os.flags(__ios_base::scientific | __ios_base::left); 1571 __os.fill(__space); 1572 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1573 1574 __os << __x._M_saved_available << __space 1575 << __x.mean() << __space 1576 << __x.sigma(); 1577 if (__x._M_saved_available) 1578 __os << __space << __x._M_saved; 1579 1580 __os.flags(__flags); 1581 __os.fill(__fill); 1582 __os.precision(__precision); 1583 return __os; 1584 } 1585 1586 template<typename _RealType, typename _CharT, typename _Traits> 1587 std::basic_istream<_CharT, _Traits>& 1588 operator>>(std::basic_istream<_CharT, _Traits>& __is, 1589 normal_distribution<_RealType>& __x) 1590 { 1591 typedef std::basic_istream<_CharT, _Traits> __istream_type; 1592 typedef typename __istream_type::ios_base __ios_base; 1593 1594 const typename __ios_base::fmtflags __flags = __is.flags(); 1595 __is.flags(__ios_base::dec | __ios_base::skipws); 1596 1597 __is >> __x._M_saved_available >> __x._M_mean 1598 >> __x._M_sigma; 1599 if (__x._M_saved_available) 1600 __is >> __x._M_saved; 1601 1602 __is.flags(__flags); 1603 return __is; 1604 } 1605 1606 1607 template<typename _RealType> 1608 void 1609 gamma_distribution<_RealType>:: 1610 _M_initialize() 1611 { 1612 if (_M_alpha >= 1) 1613 _M_l_d = std::sqrt(2 * _M_alpha - 1); 1614 else 1615 _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha)) 1616 * (1 - _M_alpha)); 1617 } 1618 1619 /** 1620 * Cheng's rejection algorithm GB for alpha >= 1 and a modification 1621 * of Vaduva's rejection from Weibull algorithm due to Devroye for 1622 * alpha < 1. 1623 * 1624 * References: 1625 * Cheng, R. C. The Generation of Gamma Random Variables with Non-integral 1626 * Shape Parameter. Applied Statistics, 26, 71-75, 1977. 1627 * 1628 * Vaduva, I. Computer Generation of Gamma Gandom Variables by Rejection 1629 * and Composition Procedures. Math. Operationsforschung and Statistik, 1630 * Series in Statistics, 8, 545-576, 1977. 1631 * 1632 * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1633 * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!). 1634 */ 1635 template<typename _RealType> 1636 template<class _UniformRandomNumberGenerator> 1637 typename gamma_distribution<_RealType>::result_type 1638 gamma_distribution<_RealType>:: 1639 operator()(_UniformRandomNumberGenerator& __urng) 1640 { 1641 result_type __x; 1642 1643 bool __reject; 1644 if (_M_alpha >= 1) 1645 { 1646 // alpha - log(4) 1647 const result_type __b = _M_alpha 1648 - result_type(1.3862943611198906188344642429163531L); 1649 const result_type __c = _M_alpha + _M_l_d; 1650 const result_type __1l = 1 / _M_l_d; 1651 1652 // 1 + log(9 / 2) 1653 const result_type __k = 2.5040773967762740733732583523868748L; 1654 1655 do 1656 { 1657 const result_type __u = __urng(); 1658 const result_type __v = __urng(); 1659 1660 const result_type __y = __1l * std::log(__v / (1 - __v)); 1661 __x = _M_alpha * std::exp(__y); 1662 1663 const result_type __z = __u * __v * __v; 1664 const result_type __r = __b + __c * __y - __x; 1665 1666 __reject = __r < result_type(4.5) * __z - __k; 1667 if (__reject) 1668 __reject = __r < std::log(__z); 1669 } 1670 while (__reject); 1671 } 1672 else 1673 { 1674 const result_type __c = 1 / _M_alpha; 1675 1676 do 1677 { 1678 const result_type __z = -std::log(__urng()); 1679 const result_type __e = -std::log(__urng()); 1680 1681 __x = std::pow(__z, __c); 1682 1683 __reject = __z + __e < _M_l_d + __x; 1684 } 1685 while (__reject); 1686 } 1687 1688 return __x; 1689 } 1690 1691 template<typename _RealType, typename _CharT, typename _Traits> 1692 std::basic_ostream<_CharT, _Traits>& 1693 operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1694 const gamma_distribution<_RealType>& __x) 1695 { 1696 typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1697 typedef typename __ostream_type::ios_base __ios_base; 1698 1699 const typename __ios_base::fmtflags __flags = __os.flags(); 1700 const _CharT __fill = __os.fill(); 1701 const std::streamsize __precision = __os.precision(); 1702 __os.flags(__ios_base::scientific | __ios_base::left); 1703 __os.fill(__os.widen(' ')); 1704 __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1705 1706 __os << __x.alpha(); 1707 1708 __os.flags(__flags); 1709 __os.fill(__fill); 1710 __os.precision(__precision); 1711 return __os; 1712 } 1713} 1714 1715_GLIBCXX_END_NAMESPACE_VERSION 1716} 1717 1718#endif 1719