1169691Skan// random number generation (out of line) -*- C++ -*- 2169691Skan 3169691Skan// Copyright (C) 2006, 2007 Free Software Foundation, Inc. 4169691Skan// 5169691Skan// This file is part of the GNU ISO C++ Library. This library is free 6169691Skan// software; you can redistribute it and/or modify it under the 7169691Skan// terms of the GNU General Public License as published by the 8169691Skan// Free Software Foundation; either version 2, or (at your option) 9169691Skan// any later version. 10169691Skan 11169691Skan// This library is distributed in the hope that it will be useful, 12169691Skan// but WITHOUT ANY WARRANTY; without even the implied warranty of 13169691Skan// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14169691Skan// GNU General Public License for more details. 15169691Skan 16169691Skan// You should have received a copy of the GNU General Public License along 17169691Skan// with this library; see the file COPYING. If not, write to the Free 18169691Skan// Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 19169691Skan// USA. 20169691Skan 21169691Skan// As a special exception, you may use this file as part of a free software 22169691Skan// library without restriction. Specifically, if other files instantiate 23169691Skan// templates or use macros or inline functions from this file, or you compile 24169691Skan// this file and link it with other files to produce an executable, this 25169691Skan// file does not by itself cause the resulting executable to be covered by 26169691Skan// the GNU General Public License. This exception does not however 27169691Skan// invalidate any other reasons why the executable file might be covered by 28169691Skan// the GNU General Public License. 29169691Skan 30169691Skan/** @file tr1/random.tcc 31169691Skan * This is a TR1 C++ Library header. 32169691Skan */ 33169691Skan 34169691Skannamespace std 35169691Skan{ 36169691Skan_GLIBCXX_BEGIN_NAMESPACE(tr1) 37169691Skan 38169691Skan /* 39169691Skan * (Further) implementation-space details. 40169691Skan */ 41169691Skan namespace __detail 42169691Skan { 43169691Skan // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid 44169691Skan // integer overflow. 45169691Skan // 46169691Skan // Because a and c are compile-time integral constants the compiler kindly 47169691Skan // elides any unreachable paths. 48169691Skan // 49169691Skan // Preconditions: a > 0, m > 0. 50169691Skan // 51169691Skan template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool> 52169691Skan struct _Mod 53169691Skan { 54169691Skan static _Tp 55169691Skan __calc(_Tp __x) 56169691Skan { 57169691Skan if (__a == 1) 58169691Skan __x %= __m; 59169691Skan else 60169691Skan { 61169691Skan static const _Tp __q = __m / __a; 62169691Skan static const _Tp __r = __m % __a; 63169691Skan 64169691Skan _Tp __t1 = __a * (__x % __q); 65169691Skan _Tp __t2 = __r * (__x / __q); 66169691Skan if (__t1 >= __t2) 67169691Skan __x = __t1 - __t2; 68169691Skan else 69169691Skan __x = __m - __t2 + __t1; 70169691Skan } 71169691Skan 72169691Skan if (__c != 0) 73169691Skan { 74169691Skan const _Tp __d = __m - __x; 75169691Skan if (__d > __c) 76169691Skan __x += __c; 77169691Skan else 78169691Skan __x = __c - __d; 79169691Skan } 80169691Skan return __x; 81169691Skan } 82169691Skan }; 83169691Skan 84169691Skan // Special case for m == 0 -- use unsigned integer overflow as modulo 85169691Skan // operator. 86169691Skan template<typename _Tp, _Tp __a, _Tp __c, _Tp __m> 87169691Skan struct _Mod<_Tp, __a, __c, __m, true> 88169691Skan { 89169691Skan static _Tp 90169691Skan __calc(_Tp __x) 91169691Skan { return __a * __x + __c; } 92169691Skan }; 93169691Skan } // namespace __detail 94169691Skan 95169691Skan /** 96169691Skan * Seeds the LCR with integral value @p __x0, adjusted so that the 97169691Skan * ring identity is never a member of the convergence set. 98169691Skan */ 99169691Skan template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 100169691Skan void 101169691Skan linear_congruential<_UIntType, __a, __c, __m>:: 102169691Skan seed(unsigned long __x0) 103169691Skan { 104169691Skan if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) 105169691Skan && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) 106169691Skan _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); 107169691Skan else 108169691Skan _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); 109169691Skan } 110169691Skan 111169691Skan /** 112169691Skan * Seeds the LCR engine with a value generated by @p __g. 113169691Skan */ 114169691Skan template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 115169691Skan template<class _Gen> 116169691Skan void 117169691Skan linear_congruential<_UIntType, __a, __c, __m>:: 118169691Skan seed(_Gen& __g, false_type) 119169691Skan { 120169691Skan _UIntType __x0 = __g(); 121169691Skan if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0) 122169691Skan && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0)) 123169691Skan _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1); 124169691Skan else 125169691Skan _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0); 126169691Skan } 127169691Skan 128169691Skan /** 129169691Skan * Gets the next generated value in sequence. 130169691Skan */ 131169691Skan template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 132169691Skan typename linear_congruential<_UIntType, __a, __c, __m>::result_type 133169691Skan linear_congruential<_UIntType, __a, __c, __m>:: 134169691Skan operator()() 135169691Skan { 136169691Skan _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x); 137169691Skan return _M_x; 138169691Skan } 139169691Skan 140169691Skan template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 141169691Skan typename _CharT, typename _Traits> 142169691Skan std::basic_ostream<_CharT, _Traits>& 143169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 144169691Skan const linear_congruential<_UIntType, __a, __c, __m>& __lcr) 145169691Skan { 146169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 147169691Skan typedef typename __ostream_type::ios_base __ios_base; 148169691Skan 149169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 150169691Skan const _CharT __fill = __os.fill(); 151169691Skan __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 152169691Skan __os.fill(__os.widen(' ')); 153169691Skan 154169691Skan __os << __lcr._M_x; 155169691Skan 156169691Skan __os.flags(__flags); 157169691Skan __os.fill(__fill); 158169691Skan return __os; 159169691Skan } 160169691Skan 161169691Skan template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 162169691Skan typename _CharT, typename _Traits> 163169691Skan std::basic_istream<_CharT, _Traits>& 164169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 165169691Skan linear_congruential<_UIntType, __a, __c, __m>& __lcr) 166169691Skan { 167169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 168169691Skan typedef typename __istream_type::ios_base __ios_base; 169169691Skan 170169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 171169691Skan __is.flags(__ios_base::dec); 172169691Skan 173169691Skan __is >> __lcr._M_x; 174169691Skan 175169691Skan __is.flags(__flags); 176169691Skan return __is; 177169691Skan } 178169691Skan 179169691Skan 180169691Skan template<class _UIntType, int __w, int __n, int __m, int __r, 181169691Skan _UIntType __a, int __u, int __s, 182169691Skan _UIntType __b, int __t, _UIntType __c, int __l> 183169691Skan void 184169691Skan mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 185169691Skan __b, __t, __c, __l>:: 186169691Skan seed(unsigned long __value) 187169691Skan { 188169691Skan _M_x[0] = __detail::__mod<_UIntType, 1, 0, 189169691Skan __detail::_Shift<_UIntType, __w>::__value>(__value); 190169691Skan 191169691Skan for (int __i = 1; __i < state_size; ++__i) 192169691Skan { 193169691Skan _UIntType __x = _M_x[__i - 1]; 194169691Skan __x ^= __x >> (__w - 2); 195169691Skan __x *= 1812433253ul; 196169691Skan __x += __i; 197169691Skan _M_x[__i] = __detail::__mod<_UIntType, 1, 0, 198169691Skan __detail::_Shift<_UIntType, __w>::__value>(__x); 199169691Skan } 200169691Skan _M_p = state_size; 201169691Skan } 202169691Skan 203169691Skan template<class _UIntType, int __w, int __n, int __m, int __r, 204169691Skan _UIntType __a, int __u, int __s, 205169691Skan _UIntType __b, int __t, _UIntType __c, int __l> 206169691Skan template<class _Gen> 207169691Skan void 208169691Skan mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 209169691Skan __b, __t, __c, __l>:: 210169691Skan seed(_Gen& __gen, false_type) 211169691Skan { 212169691Skan for (int __i = 0; __i < state_size; ++__i) 213169691Skan _M_x[__i] = __detail::__mod<_UIntType, 1, 0, 214169691Skan __detail::_Shift<_UIntType, __w>::__value>(__gen()); 215169691Skan _M_p = state_size; 216169691Skan } 217169691Skan 218169691Skan template<class _UIntType, int __w, int __n, int __m, int __r, 219169691Skan _UIntType __a, int __u, int __s, 220169691Skan _UIntType __b, int __t, _UIntType __c, int __l> 221169691Skan typename 222169691Skan mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 223169691Skan __b, __t, __c, __l>::result_type 224169691Skan mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s, 225169691Skan __b, __t, __c, __l>:: 226169691Skan operator()() 227169691Skan { 228169691Skan // Reload the vector - cost is O(n) amortized over n calls. 229169691Skan if (_M_p >= state_size) 230169691Skan { 231169691Skan const _UIntType __upper_mask = (~_UIntType()) << __r; 232169691Skan const _UIntType __lower_mask = ~__upper_mask; 233169691Skan 234169691Skan for (int __k = 0; __k < (__n - __m); ++__k) 235169691Skan { 236169691Skan _UIntType __y = ((_M_x[__k] & __upper_mask) 237169691Skan | (_M_x[__k + 1] & __lower_mask)); 238169691Skan _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) 239169691Skan ^ ((__y & 0x01) ? __a : 0)); 240169691Skan } 241169691Skan 242169691Skan for (int __k = (__n - __m); __k < (__n - 1); ++__k) 243169691Skan { 244169691Skan _UIntType __y = ((_M_x[__k] & __upper_mask) 245169691Skan | (_M_x[__k + 1] & __lower_mask)); 246169691Skan _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) 247169691Skan ^ ((__y & 0x01) ? __a : 0)); 248169691Skan } 249169691Skan 250169691Skan _UIntType __y = ((_M_x[__n - 1] & __upper_mask) 251169691Skan | (_M_x[0] & __lower_mask)); 252169691Skan _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) 253169691Skan ^ ((__y & 0x01) ? __a : 0)); 254169691Skan _M_p = 0; 255169691Skan } 256169691Skan 257169691Skan // Calculate o(x(i)). 258169691Skan result_type __z = _M_x[_M_p++]; 259169691Skan __z ^= (__z >> __u); 260169691Skan __z ^= (__z << __s) & __b; 261169691Skan __z ^= (__z << __t) & __c; 262169691Skan __z ^= (__z >> __l); 263169691Skan 264169691Skan return __z; 265169691Skan } 266169691Skan 267169691Skan template<class _UIntType, int __w, int __n, int __m, int __r, 268169691Skan _UIntType __a, int __u, int __s, _UIntType __b, int __t, 269169691Skan _UIntType __c, int __l, 270169691Skan typename _CharT, typename _Traits> 271169691Skan std::basic_ostream<_CharT, _Traits>& 272169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 273169691Skan const mersenne_twister<_UIntType, __w, __n, __m, 274169691Skan __r, __a, __u, __s, __b, __t, __c, __l>& __x) 275169691Skan { 276169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 277169691Skan typedef typename __ostream_type::ios_base __ios_base; 278169691Skan 279169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 280169691Skan const _CharT __fill = __os.fill(); 281169691Skan const _CharT __space = __os.widen(' '); 282169691Skan __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 283169691Skan __os.fill(__space); 284169691Skan 285169691Skan for (int __i = 0; __i < __n - 1; ++__i) 286169691Skan __os << __x._M_x[__i] << __space; 287169691Skan __os << __x._M_x[__n - 1]; 288169691Skan 289169691Skan __os.flags(__flags); 290169691Skan __os.fill(__fill); 291169691Skan return __os; 292169691Skan } 293169691Skan 294169691Skan template<class _UIntType, int __w, int __n, int __m, int __r, 295169691Skan _UIntType __a, int __u, int __s, _UIntType __b, int __t, 296169691Skan _UIntType __c, int __l, 297169691Skan typename _CharT, typename _Traits> 298169691Skan std::basic_istream<_CharT, _Traits>& 299169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 300169691Skan mersenne_twister<_UIntType, __w, __n, __m, 301169691Skan __r, __a, __u, __s, __b, __t, __c, __l>& __x) 302169691Skan { 303169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 304169691Skan typedef typename __istream_type::ios_base __ios_base; 305169691Skan 306169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 307169691Skan __is.flags(__ios_base::dec | __ios_base::skipws); 308169691Skan 309169691Skan for (int __i = 0; __i < __n; ++__i) 310169691Skan __is >> __x._M_x[__i]; 311169691Skan 312169691Skan __is.flags(__flags); 313169691Skan return __is; 314169691Skan } 315169691Skan 316169691Skan 317169691Skan template<typename _IntType, _IntType __m, int __s, int __r> 318169691Skan void 319169691Skan subtract_with_carry<_IntType, __m, __s, __r>:: 320169691Skan seed(unsigned long __value) 321169691Skan { 322169691Skan if (__value == 0) 323169691Skan __value = 19780503; 324169691Skan 325169691Skan std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563> 326169691Skan __lcg(__value); 327169691Skan 328169691Skan for (int __i = 0; __i < long_lag; ++__i) 329169691Skan _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg()); 330169691Skan 331169691Skan _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 332169691Skan _M_p = 0; 333169691Skan } 334169691Skan 335169691Skan template<typename _IntType, _IntType __m, int __s, int __r> 336169691Skan template<class _Gen> 337169691Skan void 338169691Skan subtract_with_carry<_IntType, __m, __s, __r>:: 339169691Skan seed(_Gen& __gen, false_type) 340169691Skan { 341169691Skan const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32; 342169691Skan 343169691Skan for (int __i = 0; __i < long_lag; ++__i) 344169691Skan { 345169691Skan _UIntType __tmp = 0; 346169691Skan _UIntType __factor = 1; 347169691Skan for (int __j = 0; __j < __n; ++__j) 348169691Skan { 349169691Skan __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0> 350169691Skan (__gen()) * __factor; 351169691Skan __factor *= __detail::_Shift<_UIntType, 32>::__value; 352169691Skan } 353169691Skan _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp); 354169691Skan } 355169691Skan _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 356169691Skan _M_p = 0; 357169691Skan } 358169691Skan 359169691Skan template<typename _IntType, _IntType __m, int __s, int __r> 360169691Skan typename subtract_with_carry<_IntType, __m, __s, __r>::result_type 361169691Skan subtract_with_carry<_IntType, __m, __s, __r>:: 362169691Skan operator()() 363169691Skan { 364169691Skan // Derive short lag index from current index. 365169691Skan int __ps = _M_p - short_lag; 366169691Skan if (__ps < 0) 367169691Skan __ps += long_lag; 368169691Skan 369169691Skan // Calculate new x(i) without overflow or division. 370169691Skan // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry 371169691Skan // cannot overflow. 372169691Skan _UIntType __xi; 373169691Skan if (_M_x[__ps] >= _M_x[_M_p] + _M_carry) 374169691Skan { 375169691Skan __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry; 376169691Skan _M_carry = 0; 377169691Skan } 378169691Skan else 379169691Skan { 380169691Skan __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps]; 381169691Skan _M_carry = 1; 382169691Skan } 383169691Skan _M_x[_M_p] = __xi; 384169691Skan 385169691Skan // Adjust current index to loop around in ring buffer. 386169691Skan if (++_M_p >= long_lag) 387169691Skan _M_p = 0; 388169691Skan 389169691Skan return __xi; 390169691Skan } 391169691Skan 392169691Skan template<typename _IntType, _IntType __m, int __s, int __r, 393169691Skan typename _CharT, typename _Traits> 394169691Skan std::basic_ostream<_CharT, _Traits>& 395169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 396169691Skan const subtract_with_carry<_IntType, __m, __s, __r>& __x) 397169691Skan { 398169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 399169691Skan typedef typename __ostream_type::ios_base __ios_base; 400169691Skan 401169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 402169691Skan const _CharT __fill = __os.fill(); 403169691Skan const _CharT __space = __os.widen(' '); 404169691Skan __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 405169691Skan __os.fill(__space); 406169691Skan 407169691Skan for (int __i = 0; __i < __r; ++__i) 408169691Skan __os << __x._M_x[__i] << __space; 409169691Skan __os << __x._M_carry; 410169691Skan 411169691Skan __os.flags(__flags); 412169691Skan __os.fill(__fill); 413169691Skan return __os; 414169691Skan } 415169691Skan 416169691Skan template<typename _IntType, _IntType __m, int __s, int __r, 417169691Skan typename _CharT, typename _Traits> 418169691Skan std::basic_istream<_CharT, _Traits>& 419169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 420169691Skan subtract_with_carry<_IntType, __m, __s, __r>& __x) 421169691Skan { 422169691Skan typedef std::basic_ostream<_CharT, _Traits> __istream_type; 423169691Skan typedef typename __istream_type::ios_base __ios_base; 424169691Skan 425169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 426169691Skan __is.flags(__ios_base::dec | __ios_base::skipws); 427169691Skan 428169691Skan for (int __i = 0; __i < __r; ++__i) 429169691Skan __is >> __x._M_x[__i]; 430169691Skan __is >> __x._M_carry; 431169691Skan 432169691Skan __is.flags(__flags); 433169691Skan return __is; 434169691Skan } 435169691Skan 436169691Skan 437169691Skan template<typename _RealType, int __w, int __s, int __r> 438169691Skan void 439169691Skan subtract_with_carry_01<_RealType, __w, __s, __r>:: 440169691Skan _M_initialize_npows() 441169691Skan { 442169691Skan for (int __j = 0; __j < __n; ++__j) 443169691Skan#if _GLIBCXX_USE_C99_MATH_TR1 444169691Skan _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32); 445169691Skan#else 446169691Skan _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32); 447169691Skan#endif 448169691Skan } 449169691Skan 450169691Skan template<typename _RealType, int __w, int __s, int __r> 451169691Skan void 452169691Skan subtract_with_carry_01<_RealType, __w, __s, __r>:: 453169691Skan seed(unsigned long __value) 454169691Skan { 455169691Skan if (__value == 0) 456169691Skan __value = 19780503; 457169691Skan 458169691Skan // _GLIBCXX_RESOLVE_LIB_DEFECTS 459169691Skan // 512. Seeding subtract_with_carry_01 from a single unsigned long. 460169691Skan std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563> 461169691Skan __lcg(__value); 462169691Skan 463169691Skan this->seed(__lcg); 464169691Skan } 465169691Skan 466169691Skan template<typename _RealType, int __w, int __s, int __r> 467169691Skan template<class _Gen> 468169691Skan void 469169691Skan subtract_with_carry_01<_RealType, __w, __s, __r>:: 470169691Skan seed(_Gen& __gen, false_type) 471169691Skan { 472169691Skan for (int __i = 0; __i < long_lag; ++__i) 473169691Skan { 474169691Skan for (int __j = 0; __j < __n - 1; ++__j) 475169691Skan _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen()); 476169691Skan _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, 477169691Skan __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen()); 478169691Skan } 479169691Skan 480169691Skan _M_carry = 1; 481169691Skan for (int __j = 0; __j < __n; ++__j) 482169691Skan if (_M_x[long_lag - 1][__j] != 0) 483169691Skan { 484169691Skan _M_carry = 0; 485169691Skan break; 486169691Skan } 487169691Skan 488169691Skan _M_p = 0; 489169691Skan } 490169691Skan 491169691Skan template<typename _RealType, int __w, int __s, int __r> 492169691Skan typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type 493169691Skan subtract_with_carry_01<_RealType, __w, __s, __r>:: 494169691Skan operator()() 495169691Skan { 496169691Skan // Derive short lag index from current index. 497169691Skan int __ps = _M_p - short_lag; 498169691Skan if (__ps < 0) 499169691Skan __ps += long_lag; 500169691Skan 501169691Skan _UInt32Type __new_carry; 502169691Skan for (int __j = 0; __j < __n - 1; ++__j) 503169691Skan { 504169691Skan if (_M_x[__ps][__j] > _M_x[_M_p][__j] 505169691Skan || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0)) 506169691Skan __new_carry = 0; 507169691Skan else 508169691Skan __new_carry = 1; 509169691Skan 510169691Skan _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry; 511169691Skan _M_carry = __new_carry; 512169691Skan } 513169691Skan 514169691Skan if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1] 515169691Skan || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0)) 516169691Skan __new_carry = 0; 517169691Skan else 518169691Skan __new_carry = 1; 519169691Skan 520169691Skan _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0, 521169691Skan __detail::_Shift<_UInt32Type, __w % 32>::__value> 522169691Skan (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry); 523169691Skan _M_carry = __new_carry; 524169691Skan 525169691Skan result_type __ret = 0.0; 526169691Skan for (int __j = 0; __j < __n; ++__j) 527169691Skan __ret += _M_x[_M_p][__j] * _M_npows[__j]; 528169691Skan 529169691Skan // Adjust current index to loop around in ring buffer. 530169691Skan if (++_M_p >= long_lag) 531169691Skan _M_p = 0; 532169691Skan 533169691Skan return __ret; 534169691Skan } 535169691Skan 536169691Skan template<typename _RealType, int __w, int __s, int __r, 537169691Skan typename _CharT, typename _Traits> 538169691Skan std::basic_ostream<_CharT, _Traits>& 539169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 540169691Skan const subtract_with_carry_01<_RealType, __w, __s, __r>& __x) 541169691Skan { 542169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 543169691Skan typedef typename __ostream_type::ios_base __ios_base; 544169691Skan 545169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 546169691Skan const _CharT __fill = __os.fill(); 547169691Skan const _CharT __space = __os.widen(' '); 548169691Skan __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 549169691Skan __os.fill(__space); 550169691Skan 551169691Skan for (int __i = 0; __i < __r; ++__i) 552169691Skan for (int __j = 0; __j < __x.__n; ++__j) 553169691Skan __os << __x._M_x[__i][__j] << __space; 554169691Skan __os << __x._M_carry; 555169691Skan 556169691Skan __os.flags(__flags); 557169691Skan __os.fill(__fill); 558169691Skan return __os; 559169691Skan } 560169691Skan 561169691Skan template<typename _RealType, int __w, int __s, int __r, 562169691Skan typename _CharT, typename _Traits> 563169691Skan std::basic_istream<_CharT, _Traits>& 564169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 565169691Skan subtract_with_carry_01<_RealType, __w, __s, __r>& __x) 566169691Skan { 567169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 568169691Skan typedef typename __istream_type::ios_base __ios_base; 569169691Skan 570169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 571169691Skan __is.flags(__ios_base::dec | __ios_base::skipws); 572169691Skan 573169691Skan for (int __i = 0; __i < __r; ++__i) 574169691Skan for (int __j = 0; __j < __x.__n; ++__j) 575169691Skan __is >> __x._M_x[__i][__j]; 576169691Skan __is >> __x._M_carry; 577169691Skan 578169691Skan __is.flags(__flags); 579169691Skan return __is; 580169691Skan } 581169691Skan 582169691Skan 583169691Skan template<class _UniformRandomNumberGenerator, int __p, int __r> 584169691Skan typename discard_block<_UniformRandomNumberGenerator, 585169691Skan __p, __r>::result_type 586169691Skan discard_block<_UniformRandomNumberGenerator, __p, __r>:: 587169691Skan operator()() 588169691Skan { 589169691Skan if (_M_n >= used_block) 590169691Skan { 591169691Skan while (_M_n < block_size) 592169691Skan { 593169691Skan _M_b(); 594169691Skan ++_M_n; 595169691Skan } 596169691Skan _M_n = 0; 597169691Skan } 598169691Skan ++_M_n; 599169691Skan return _M_b(); 600169691Skan } 601169691Skan 602169691Skan template<class _UniformRandomNumberGenerator, int __p, int __r, 603169691Skan typename _CharT, typename _Traits> 604169691Skan std::basic_ostream<_CharT, _Traits>& 605169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 606169691Skan const discard_block<_UniformRandomNumberGenerator, 607169691Skan __p, __r>& __x) 608169691Skan { 609169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 610169691Skan typedef typename __ostream_type::ios_base __ios_base; 611169691Skan 612169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 613169691Skan const _CharT __fill = __os.fill(); 614169691Skan const _CharT __space = __os.widen(' '); 615169691Skan __os.flags(__ios_base::dec | __ios_base::fixed 616169691Skan | __ios_base::left); 617169691Skan __os.fill(__space); 618169691Skan 619169691Skan __os << __x._M_b << __space << __x._M_n; 620169691Skan 621169691Skan __os.flags(__flags); 622169691Skan __os.fill(__fill); 623169691Skan return __os; 624169691Skan } 625169691Skan 626169691Skan template<class _UniformRandomNumberGenerator, int __p, int __r, 627169691Skan typename _CharT, typename _Traits> 628169691Skan std::basic_istream<_CharT, _Traits>& 629169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 630169691Skan discard_block<_UniformRandomNumberGenerator, __p, __r>& __x) 631169691Skan { 632169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 633169691Skan typedef typename __istream_type::ios_base __ios_base; 634169691Skan 635169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 636169691Skan __is.flags(__ios_base::dec | __ios_base::skipws); 637169691Skan 638169691Skan __is >> __x._M_b >> __x._M_n; 639169691Skan 640169691Skan __is.flags(__flags); 641169691Skan return __is; 642169691Skan } 643169691Skan 644169691Skan 645169691Skan template<class _UniformRandomNumberGenerator1, int __s1, 646169691Skan class _UniformRandomNumberGenerator2, int __s2> 647169691Skan void 648169691Skan xor_combine<_UniformRandomNumberGenerator1, __s1, 649169691Skan _UniformRandomNumberGenerator2, __s2>:: 650169691Skan _M_initialize_max() 651169691Skan { 652169691Skan const int __w = std::numeric_limits<result_type>::digits; 653169691Skan 654169691Skan const result_type __m1 = 655169691Skan std::min(result_type(_M_b1.max() - _M_b1.min()), 656169691Skan __detail::_Shift<result_type, __w - __s1>::__value - 1); 657169691Skan 658169691Skan const result_type __m2 = 659169691Skan std::min(result_type(_M_b2.max() - _M_b2.min()), 660169691Skan __detail::_Shift<result_type, __w - __s2>::__value - 1); 661169691Skan 662169691Skan // NB: In TR1 s1 is not required to be >= s2. 663169691Skan if (__s1 < __s2) 664169691Skan _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1; 665169691Skan else 666169691Skan _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2; 667169691Skan } 668169691Skan 669169691Skan template<class _UniformRandomNumberGenerator1, int __s1, 670169691Skan class _UniformRandomNumberGenerator2, int __s2> 671169691Skan typename xor_combine<_UniformRandomNumberGenerator1, __s1, 672169691Skan _UniformRandomNumberGenerator2, __s2>::result_type 673169691Skan xor_combine<_UniformRandomNumberGenerator1, __s1, 674169691Skan _UniformRandomNumberGenerator2, __s2>:: 675169691Skan _M_initialize_max_aux(result_type __a, result_type __b, int __d) 676169691Skan { 677169691Skan const result_type __two2d = result_type(1) << __d; 678169691Skan const result_type __c = __a * __two2d; 679169691Skan 680169691Skan if (__a == 0 || __b < __two2d) 681169691Skan return __c + __b; 682169691Skan 683169691Skan const result_type __t = std::max(__c, __b); 684169691Skan const result_type __u = std::min(__c, __b); 685169691Skan 686169691Skan result_type __ub = __u; 687169691Skan result_type __p; 688169691Skan for (__p = 0; __ub != 1; __ub >>= 1) 689169691Skan ++__p; 690169691Skan 691169691Skan const result_type __two2p = result_type(1) << __p; 692169691Skan const result_type __k = __t / __two2p; 693169691Skan 694169691Skan if (__k & 1) 695169691Skan return (__k + 1) * __two2p - 1; 696169691Skan 697169691Skan if (__c >= __b) 698169691Skan return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p) 699169691Skan / __two2d, 700169691Skan __u % __two2p, __d); 701169691Skan else 702169691Skan return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p) 703169691Skan / __two2d, 704169691Skan __t % __two2p, __d); 705169691Skan } 706169691Skan 707169691Skan template<class _UniformRandomNumberGenerator1, int __s1, 708169691Skan class _UniformRandomNumberGenerator2, int __s2, 709169691Skan typename _CharT, typename _Traits> 710169691Skan std::basic_ostream<_CharT, _Traits>& 711169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 712169691Skan const xor_combine<_UniformRandomNumberGenerator1, __s1, 713169691Skan _UniformRandomNumberGenerator2, __s2>& __x) 714169691Skan { 715169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 716169691Skan typedef typename __ostream_type::ios_base __ios_base; 717169691Skan 718169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 719169691Skan const _CharT __fill = __os.fill(); 720169691Skan const _CharT __space = __os.widen(' '); 721169691Skan __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 722169691Skan __os.fill(__space); 723169691Skan 724169691Skan __os << __x.base1() << __space << __x.base2(); 725169691Skan 726169691Skan __os.flags(__flags); 727169691Skan __os.fill(__fill); 728169691Skan return __os; 729169691Skan } 730169691Skan 731169691Skan template<class _UniformRandomNumberGenerator1, int __s1, 732169691Skan class _UniformRandomNumberGenerator2, int __s2, 733169691Skan typename _CharT, typename _Traits> 734169691Skan std::basic_istream<_CharT, _Traits>& 735169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 736169691Skan xor_combine<_UniformRandomNumberGenerator1, __s1, 737169691Skan _UniformRandomNumberGenerator2, __s2>& __x) 738169691Skan { 739169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 740169691Skan typedef typename __istream_type::ios_base __ios_base; 741169691Skan 742169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 743169691Skan __is.flags(__ios_base::skipws); 744169691Skan 745169691Skan __is >> __x._M_b1 >> __x._M_b2; 746169691Skan 747169691Skan __is.flags(__flags); 748169691Skan return __is; 749169691Skan } 750169691Skan 751169691Skan 752169691Skan template<typename _IntType, typename _CharT, typename _Traits> 753169691Skan std::basic_ostream<_CharT, _Traits>& 754169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 755169691Skan const uniform_int<_IntType>& __x) 756169691Skan { 757169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 758169691Skan typedef typename __ostream_type::ios_base __ios_base; 759169691Skan 760169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 761169691Skan const _CharT __fill = __os.fill(); 762169691Skan const _CharT __space = __os.widen(' '); 763169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 764169691Skan __os.fill(__space); 765169691Skan 766169691Skan __os << __x.min() << __space << __x.max(); 767169691Skan 768169691Skan __os.flags(__flags); 769169691Skan __os.fill(__fill); 770169691Skan return __os; 771169691Skan } 772169691Skan 773169691Skan template<typename _IntType, typename _CharT, typename _Traits> 774169691Skan std::basic_istream<_CharT, _Traits>& 775169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 776169691Skan uniform_int<_IntType>& __x) 777169691Skan { 778169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 779169691Skan typedef typename __istream_type::ios_base __ios_base; 780169691Skan 781169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 782169691Skan __is.flags(__ios_base::dec | __ios_base::skipws); 783169691Skan 784169691Skan __is >> __x._M_min >> __x._M_max; 785169691Skan 786169691Skan __is.flags(__flags); 787169691Skan return __is; 788169691Skan } 789169691Skan 790169691Skan 791169691Skan template<typename _CharT, typename _Traits> 792169691Skan std::basic_ostream<_CharT, _Traits>& 793169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 794169691Skan const bernoulli_distribution& __x) 795169691Skan { 796169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 797169691Skan typedef typename __ostream_type::ios_base __ios_base; 798169691Skan 799169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 800169691Skan const _CharT __fill = __os.fill(); 801169691Skan const std::streamsize __precision = __os.precision(); 802169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 803169691Skan __os.fill(__os.widen(' ')); 804169691Skan __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10); 805169691Skan 806169691Skan __os << __x.p(); 807169691Skan 808169691Skan __os.flags(__flags); 809169691Skan __os.fill(__fill); 810169691Skan __os.precision(__precision); 811169691Skan return __os; 812169691Skan } 813169691Skan 814169691Skan 815169691Skan template<typename _IntType, typename _RealType> 816169691Skan template<class _UniformRandomNumberGenerator> 817169691Skan typename geometric_distribution<_IntType, _RealType>::result_type 818169691Skan geometric_distribution<_IntType, _RealType>:: 819169691Skan operator()(_UniformRandomNumberGenerator& __urng) 820169691Skan { 821169691Skan // About the epsilon thing see this thread: 822169691Skan // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html 823169691Skan const _RealType __naf = 824169691Skan (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 825169691Skan // The largest _RealType convertible to _IntType. 826169691Skan const _RealType __thr = 827169691Skan std::numeric_limits<_IntType>::max() + __naf; 828169691Skan 829169691Skan _RealType __cand; 830169691Skan do 831169691Skan __cand = std::ceil(std::log(__urng()) / _M_log_p); 832169691Skan while (__cand >= __thr); 833169691Skan 834169691Skan return result_type(__cand + __naf); 835169691Skan } 836169691Skan 837169691Skan template<typename _IntType, typename _RealType, 838169691Skan typename _CharT, typename _Traits> 839169691Skan std::basic_ostream<_CharT, _Traits>& 840169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 841169691Skan const geometric_distribution<_IntType, _RealType>& __x) 842169691Skan { 843169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 844169691Skan typedef typename __ostream_type::ios_base __ios_base; 845169691Skan 846169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 847169691Skan const _CharT __fill = __os.fill(); 848169691Skan const std::streamsize __precision = __os.precision(); 849169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 850169691Skan __os.fill(__os.widen(' ')); 851169691Skan __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 852169691Skan 853169691Skan __os << __x.p(); 854169691Skan 855169691Skan __os.flags(__flags); 856169691Skan __os.fill(__fill); 857169691Skan __os.precision(__precision); 858169691Skan return __os; 859169691Skan } 860169691Skan 861169691Skan 862169691Skan template<typename _IntType, typename _RealType> 863169691Skan void 864169691Skan poisson_distribution<_IntType, _RealType>:: 865169691Skan _M_initialize() 866169691Skan { 867169691Skan#if _GLIBCXX_USE_C99_MATH_TR1 868169691Skan if (_M_mean >= 12) 869169691Skan { 870169691Skan const _RealType __m = std::floor(_M_mean); 871169691Skan _M_lm_thr = std::log(_M_mean); 872169691Skan _M_lfm = std::tr1::lgamma(__m + 1); 873169691Skan _M_sm = std::sqrt(__m); 874169691Skan 875169691Skan const _RealType __pi_4 = 0.7853981633974483096156608458198757L; 876169691Skan const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m 877169691Skan / __pi_4)); 878169691Skan _M_d = std::tr1::round(std::max(_RealType(6), 879169691Skan std::min(__m, __dx))); 880169691Skan const _RealType __cx = 2 * __m + _M_d; 881169691Skan _M_scx = std::sqrt(__cx / 2); 882169691Skan _M_1cx = 1 / __cx; 883169691Skan 884169691Skan _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); 885169691Skan _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d; 886169691Skan } 887169691Skan else 888169691Skan#endif 889169691Skan _M_lm_thr = std::exp(-_M_mean); 890169691Skan } 891169691Skan 892169691Skan /** 893169691Skan * A rejection algorithm when mean >= 12 and a simple method based 894169691Skan * upon the multiplication of uniform random variates otherwise. 895169691Skan * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 896169691Skan * is defined. 897169691Skan * 898169691Skan * Reference: 899169691Skan * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 900169691Skan * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). 901169691Skan */ 902169691Skan template<typename _IntType, typename _RealType> 903169691Skan template<class _UniformRandomNumberGenerator> 904169691Skan typename poisson_distribution<_IntType, _RealType>::result_type 905169691Skan poisson_distribution<_IntType, _RealType>:: 906169691Skan operator()(_UniformRandomNumberGenerator& __urng) 907169691Skan { 908169691Skan#if _GLIBCXX_USE_C99_MATH_TR1 909169691Skan if (_M_mean >= 12) 910169691Skan { 911169691Skan _RealType __x; 912169691Skan 913169691Skan // See comments above... 914169691Skan const _RealType __naf = 915169691Skan (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 916169691Skan const _RealType __thr = 917169691Skan std::numeric_limits<_IntType>::max() + __naf; 918169691Skan 919169691Skan const _RealType __m = std::floor(_M_mean); 920169691Skan // sqrt(pi / 2) 921169691Skan const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 922169691Skan const _RealType __c1 = _M_sm * __spi_2; 923169691Skan const _RealType __c2 = _M_c2b + __c1; 924169691Skan const _RealType __c3 = __c2 + 1; 925169691Skan const _RealType __c4 = __c3 + 1; 926169691Skan // e^(1 / 78) 927169691Skan const _RealType __e178 = 1.0129030479320018583185514777512983L; 928169691Skan const _RealType __c5 = __c4 + __e178; 929169691Skan const _RealType __c = _M_cb + __c5; 930169691Skan const _RealType __2cx = 2 * (2 * __m + _M_d); 931169691Skan 932169691Skan bool __reject = true; 933169691Skan do 934169691Skan { 935169691Skan const _RealType __u = __c * __urng(); 936169691Skan const _RealType __e = -std::log(__urng()); 937169691Skan 938169691Skan _RealType __w = 0.0; 939169691Skan 940169691Skan if (__u <= __c1) 941169691Skan { 942169691Skan const _RealType __n = _M_nd(__urng); 943169691Skan const _RealType __y = -std::abs(__n) * _M_sm - 1; 944169691Skan __x = std::floor(__y); 945169691Skan __w = -__n * __n / 2; 946169691Skan if (__x < -__m) 947169691Skan continue; 948169691Skan } 949169691Skan else if (__u <= __c2) 950169691Skan { 951169691Skan const _RealType __n = _M_nd(__urng); 952169691Skan const _RealType __y = 1 + std::abs(__n) * _M_scx; 953169691Skan __x = std::ceil(__y); 954169691Skan __w = __y * (2 - __y) * _M_1cx; 955169691Skan if (__x > _M_d) 956169691Skan continue; 957169691Skan } 958169691Skan else if (__u <= __c3) 959169691Skan // NB: This case not in the book, nor in the Errata, 960169691Skan // but should be ok... 961169691Skan __x = -1; 962169691Skan else if (__u <= __c4) 963169691Skan __x = 0; 964169691Skan else if (__u <= __c5) 965169691Skan __x = 1; 966169691Skan else 967169691Skan { 968169691Skan const _RealType __v = -std::log(__urng()); 969169691Skan const _RealType __y = _M_d + __v * __2cx / _M_d; 970169691Skan __x = std::ceil(__y); 971169691Skan __w = -_M_d * _M_1cx * (1 + __y / 2); 972169691Skan } 973169691Skan 974169691Skan __reject = (__w - __e - __x * _M_lm_thr 975169691Skan > _M_lfm - std::tr1::lgamma(__x + __m + 1)); 976169691Skan 977169691Skan __reject |= __x + __m >= __thr; 978169691Skan 979169691Skan } while (__reject); 980169691Skan 981169691Skan return result_type(__x + __m + __naf); 982169691Skan } 983169691Skan else 984169691Skan#endif 985169691Skan { 986169691Skan _IntType __x = 0; 987169691Skan _RealType __prod = 1.0; 988169691Skan 989169691Skan do 990169691Skan { 991169691Skan __prod *= __urng(); 992169691Skan __x += 1; 993169691Skan } 994169691Skan while (__prod > _M_lm_thr); 995169691Skan 996169691Skan return __x - 1; 997169691Skan } 998169691Skan } 999169691Skan 1000169691Skan template<typename _IntType, typename _RealType, 1001169691Skan typename _CharT, typename _Traits> 1002169691Skan std::basic_ostream<_CharT, _Traits>& 1003169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1004169691Skan const poisson_distribution<_IntType, _RealType>& __x) 1005169691Skan { 1006169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1007169691Skan typedef typename __ostream_type::ios_base __ios_base; 1008169691Skan 1009169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 1010169691Skan const _CharT __fill = __os.fill(); 1011169691Skan const std::streamsize __precision = __os.precision(); 1012169691Skan const _CharT __space = __os.widen(' '); 1013169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 1014169691Skan __os.fill(__space); 1015169691Skan __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1016169691Skan 1017169691Skan __os << __x.mean() << __space << __x._M_nd; 1018169691Skan 1019169691Skan __os.flags(__flags); 1020169691Skan __os.fill(__fill); 1021169691Skan __os.precision(__precision); 1022169691Skan return __os; 1023169691Skan } 1024169691Skan 1025169691Skan template<typename _IntType, typename _RealType, 1026169691Skan typename _CharT, typename _Traits> 1027169691Skan std::basic_istream<_CharT, _Traits>& 1028169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 1029169691Skan poisson_distribution<_IntType, _RealType>& __x) 1030169691Skan { 1031169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 1032169691Skan typedef typename __istream_type::ios_base __ios_base; 1033169691Skan 1034169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 1035169691Skan __is.flags(__ios_base::skipws); 1036169691Skan 1037169691Skan __is >> __x._M_mean >> __x._M_nd; 1038169691Skan __x._M_initialize(); 1039169691Skan 1040169691Skan __is.flags(__flags); 1041169691Skan return __is; 1042169691Skan } 1043169691Skan 1044169691Skan 1045169691Skan template<typename _IntType, typename _RealType> 1046169691Skan void 1047169691Skan binomial_distribution<_IntType, _RealType>:: 1048169691Skan _M_initialize() 1049169691Skan { 1050169691Skan const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 1051169691Skan 1052169691Skan _M_easy = true; 1053169691Skan 1054169691Skan#if _GLIBCXX_USE_C99_MATH_TR1 1055169691Skan if (_M_t * __p12 >= 8) 1056169691Skan { 1057169691Skan _M_easy = false; 1058169691Skan const _RealType __np = std::floor(_M_t * __p12); 1059169691Skan const _RealType __pa = __np / _M_t; 1060169691Skan const _RealType __1p = 1 - __pa; 1061169691Skan 1062169691Skan const _RealType __pi_4 = 0.7853981633974483096156608458198757L; 1063169691Skan const _RealType __d1x = 1064169691Skan std::sqrt(__np * __1p * std::log(32 * __np 1065169691Skan / (81 * __pi_4 * __1p))); 1066169691Skan _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x)); 1067169691Skan const _RealType __d2x = 1068169691Skan std::sqrt(__np * __1p * std::log(32 * _M_t * __1p 1069169691Skan / (__pi_4 * __pa))); 1070169691Skan _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x)); 1071169691Skan 1072169691Skan // sqrt(pi / 2) 1073169691Skan const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 1074169691Skan _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); 1075169691Skan _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); 1076169691Skan _M_c = 2 * _M_d1 / __np; 1077169691Skan _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; 1078169691Skan const _RealType __a12 = _M_a1 + _M_s2 * __spi_2; 1079169691Skan const _RealType __s1s = _M_s1 * _M_s1; 1080169691Skan _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) 1081169691Skan * 2 * __s1s / _M_d1 1082169691Skan * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); 1083169691Skan const _RealType __s2s = _M_s2 * _M_s2; 1084169691Skan _M_s = (_M_a123 + 2 * __s2s / _M_d2 1085169691Skan * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); 1086169691Skan _M_lf = (std::tr1::lgamma(__np + 1) 1087169691Skan + std::tr1::lgamma(_M_t - __np + 1)); 1088169691Skan _M_lp1p = std::log(__pa / __1p); 1089169691Skan 1090169691Skan _M_q = -std::log(1 - (__p12 - __pa) / __1p); 1091169691Skan } 1092169691Skan else 1093169691Skan#endif 1094169691Skan _M_q = -std::log(1 - __p12); 1095169691Skan } 1096169691Skan 1097169691Skan template<typename _IntType, typename _RealType> 1098169691Skan template<class _UniformRandomNumberGenerator> 1099169691Skan typename binomial_distribution<_IntType, _RealType>::result_type 1100169691Skan binomial_distribution<_IntType, _RealType>:: 1101169691Skan _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) 1102169691Skan { 1103169691Skan _IntType __x = 0; 1104169691Skan _RealType __sum = 0; 1105169691Skan 1106169691Skan do 1107169691Skan { 1108169691Skan const _RealType __e = -std::log(__urng()); 1109169691Skan __sum += __e / (__t - __x); 1110169691Skan __x += 1; 1111169691Skan } 1112169691Skan while (__sum <= _M_q); 1113169691Skan 1114169691Skan return __x - 1; 1115169691Skan } 1116169691Skan 1117169691Skan /** 1118169691Skan * A rejection algorithm when t * p >= 8 and a simple waiting time 1119169691Skan * method - the second in the referenced book - otherwise. 1120169691Skan * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 1121169691Skan * is defined. 1122169691Skan * 1123169691Skan * Reference: 1124169691Skan * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 1125169691Skan * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 1126169691Skan */ 1127169691Skan template<typename _IntType, typename _RealType> 1128169691Skan template<class _UniformRandomNumberGenerator> 1129169691Skan typename binomial_distribution<_IntType, _RealType>::result_type 1130169691Skan binomial_distribution<_IntType, _RealType>:: 1131169691Skan operator()(_UniformRandomNumberGenerator& __urng) 1132169691Skan { 1133169691Skan result_type __ret; 1134169691Skan const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 1135169691Skan 1136169691Skan#if _GLIBCXX_USE_C99_MATH_TR1 1137169691Skan if (!_M_easy) 1138169691Skan { 1139169691Skan _RealType __x; 1140169691Skan 1141169691Skan // See comments above... 1142169691Skan const _RealType __naf = 1143169691Skan (1 - std::numeric_limits<_RealType>::epsilon()) / 2; 1144169691Skan const _RealType __thr = 1145169691Skan std::numeric_limits<_IntType>::max() + __naf; 1146169691Skan 1147169691Skan const _RealType __np = std::floor(_M_t * __p12); 1148169691Skan const _RealType __pa = __np / _M_t; 1149169691Skan 1150169691Skan // sqrt(pi / 2) 1151169691Skan const _RealType __spi_2 = 1.2533141373155002512078826424055226L; 1152169691Skan const _RealType __a1 = _M_a1; 1153169691Skan const _RealType __a12 = __a1 + _M_s2 * __spi_2; 1154169691Skan const _RealType __a123 = _M_a123; 1155169691Skan const _RealType __s1s = _M_s1 * _M_s1; 1156169691Skan const _RealType __s2s = _M_s2 * _M_s2; 1157169691Skan 1158169691Skan bool __reject; 1159169691Skan do 1160169691Skan { 1161169691Skan const _RealType __u = _M_s * __urng(); 1162169691Skan 1163169691Skan _RealType __v; 1164169691Skan 1165169691Skan if (__u <= __a1) 1166169691Skan { 1167169691Skan const _RealType __n = _M_nd(__urng); 1168169691Skan const _RealType __y = _M_s1 * std::abs(__n); 1169169691Skan __reject = __y >= _M_d1; 1170169691Skan if (!__reject) 1171169691Skan { 1172169691Skan const _RealType __e = -std::log(__urng()); 1173169691Skan __x = std::floor(__y); 1174169691Skan __v = -__e - __n * __n / 2 + _M_c; 1175169691Skan } 1176169691Skan } 1177169691Skan else if (__u <= __a12) 1178169691Skan { 1179169691Skan const _RealType __n = _M_nd(__urng); 1180169691Skan const _RealType __y = _M_s2 * std::abs(__n); 1181169691Skan __reject = __y >= _M_d2; 1182169691Skan if (!__reject) 1183169691Skan { 1184169691Skan const _RealType __e = -std::log(__urng()); 1185169691Skan __x = std::floor(-__y); 1186169691Skan __v = -__e - __n * __n / 2; 1187169691Skan } 1188169691Skan } 1189169691Skan else if (__u <= __a123) 1190169691Skan { 1191169691Skan const _RealType __e1 = -std::log(__urng()); 1192169691Skan const _RealType __e2 = -std::log(__urng()); 1193169691Skan 1194169691Skan const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1; 1195169691Skan __x = std::floor(__y); 1196169691Skan __v = (-__e2 + _M_d1 * (1 / (_M_t - __np) 1197169691Skan -__y / (2 * __s1s))); 1198169691Skan __reject = false; 1199169691Skan } 1200169691Skan else 1201169691Skan { 1202169691Skan const _RealType __e1 = -std::log(__urng()); 1203169691Skan const _RealType __e2 = -std::log(__urng()); 1204169691Skan 1205169691Skan const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2; 1206169691Skan __x = std::floor(-__y); 1207169691Skan __v = -__e2 - _M_d2 * __y / (2 * __s2s); 1208169691Skan __reject = false; 1209169691Skan } 1210169691Skan 1211169691Skan __reject = __reject || __x < -__np || __x > _M_t - __np; 1212169691Skan if (!__reject) 1213169691Skan { 1214169691Skan const _RealType __lfx = 1215169691Skan std::tr1::lgamma(__np + __x + 1) 1216169691Skan + std::tr1::lgamma(_M_t - (__np + __x) + 1); 1217169691Skan __reject = __v > _M_lf - __lfx + __x * _M_lp1p; 1218169691Skan } 1219169691Skan 1220169691Skan __reject |= __x + __np >= __thr; 1221169691Skan } 1222169691Skan while (__reject); 1223169691Skan 1224169691Skan __x += __np + __naf; 1225169691Skan 1226169691Skan const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); 1227169691Skan __ret = _IntType(__x) + __z; 1228169691Skan } 1229169691Skan else 1230169691Skan#endif 1231169691Skan __ret = _M_waiting(__urng, _M_t); 1232169691Skan 1233169691Skan if (__p12 != _M_p) 1234169691Skan __ret = _M_t - __ret; 1235169691Skan return __ret; 1236169691Skan } 1237169691Skan 1238169691Skan template<typename _IntType, typename _RealType, 1239169691Skan typename _CharT, typename _Traits> 1240169691Skan std::basic_ostream<_CharT, _Traits>& 1241169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1242169691Skan const binomial_distribution<_IntType, _RealType>& __x) 1243169691Skan { 1244169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1245169691Skan typedef typename __ostream_type::ios_base __ios_base; 1246169691Skan 1247169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 1248169691Skan const _CharT __fill = __os.fill(); 1249169691Skan const std::streamsize __precision = __os.precision(); 1250169691Skan const _CharT __space = __os.widen(' '); 1251169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 1252169691Skan __os.fill(__space); 1253169691Skan __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1254169691Skan 1255169691Skan __os << __x.t() << __space << __x.p() 1256169691Skan << __space << __x._M_nd; 1257169691Skan 1258169691Skan __os.flags(__flags); 1259169691Skan __os.fill(__fill); 1260169691Skan __os.precision(__precision); 1261169691Skan return __os; 1262169691Skan } 1263169691Skan 1264169691Skan template<typename _IntType, typename _RealType, 1265169691Skan typename _CharT, typename _Traits> 1266169691Skan std::basic_istream<_CharT, _Traits>& 1267169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 1268169691Skan binomial_distribution<_IntType, _RealType>& __x) 1269169691Skan { 1270169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 1271169691Skan typedef typename __istream_type::ios_base __ios_base; 1272169691Skan 1273169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 1274169691Skan __is.flags(__ios_base::dec | __ios_base::skipws); 1275169691Skan 1276169691Skan __is >> __x._M_t >> __x._M_p >> __x._M_nd; 1277169691Skan __x._M_initialize(); 1278169691Skan 1279169691Skan __is.flags(__flags); 1280169691Skan return __is; 1281169691Skan } 1282169691Skan 1283169691Skan 1284169691Skan template<typename _RealType, typename _CharT, typename _Traits> 1285169691Skan std::basic_ostream<_CharT, _Traits>& 1286169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1287169691Skan const uniform_real<_RealType>& __x) 1288169691Skan { 1289169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1290169691Skan typedef typename __ostream_type::ios_base __ios_base; 1291169691Skan 1292169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 1293169691Skan const _CharT __fill = __os.fill(); 1294169691Skan const std::streamsize __precision = __os.precision(); 1295169691Skan const _CharT __space = __os.widen(' '); 1296169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 1297169691Skan __os.fill(__space); 1298169691Skan __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1299169691Skan 1300169691Skan __os << __x.min() << __space << __x.max(); 1301169691Skan 1302169691Skan __os.flags(__flags); 1303169691Skan __os.fill(__fill); 1304169691Skan __os.precision(__precision); 1305169691Skan return __os; 1306169691Skan } 1307169691Skan 1308169691Skan template<typename _RealType, typename _CharT, typename _Traits> 1309169691Skan std::basic_istream<_CharT, _Traits>& 1310169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 1311169691Skan uniform_real<_RealType>& __x) 1312169691Skan { 1313169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 1314169691Skan typedef typename __istream_type::ios_base __ios_base; 1315169691Skan 1316169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 1317169691Skan __is.flags(__ios_base::skipws); 1318169691Skan 1319169691Skan __is >> __x._M_min >> __x._M_max; 1320169691Skan 1321169691Skan __is.flags(__flags); 1322169691Skan return __is; 1323169691Skan } 1324169691Skan 1325169691Skan 1326169691Skan template<typename _RealType, typename _CharT, typename _Traits> 1327169691Skan std::basic_ostream<_CharT, _Traits>& 1328169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1329169691Skan const exponential_distribution<_RealType>& __x) 1330169691Skan { 1331169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1332169691Skan typedef typename __ostream_type::ios_base __ios_base; 1333169691Skan 1334169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 1335169691Skan const _CharT __fill = __os.fill(); 1336169691Skan const std::streamsize __precision = __os.precision(); 1337169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 1338169691Skan __os.fill(__os.widen(' ')); 1339169691Skan __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1340169691Skan 1341169691Skan __os << __x.lambda(); 1342169691Skan 1343169691Skan __os.flags(__flags); 1344169691Skan __os.fill(__fill); 1345169691Skan __os.precision(__precision); 1346169691Skan return __os; 1347169691Skan } 1348169691Skan 1349169691Skan 1350169691Skan /** 1351169691Skan * Polar method due to Marsaglia. 1352169691Skan * 1353169691Skan * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 1354169691Skan * New York, 1986, Ch. V, Sect. 4.4. 1355169691Skan */ 1356169691Skan template<typename _RealType> 1357169691Skan template<class _UniformRandomNumberGenerator> 1358169691Skan typename normal_distribution<_RealType>::result_type 1359169691Skan normal_distribution<_RealType>:: 1360169691Skan operator()(_UniformRandomNumberGenerator& __urng) 1361169691Skan { 1362169691Skan result_type __ret; 1363169691Skan 1364169691Skan if (_M_saved_available) 1365169691Skan { 1366169691Skan _M_saved_available = false; 1367169691Skan __ret = _M_saved; 1368169691Skan } 1369169691Skan else 1370169691Skan { 1371169691Skan result_type __x, __y, __r2; 1372169691Skan do 1373169691Skan { 1374169691Skan __x = result_type(2.0) * __urng() - 1.0; 1375169691Skan __y = result_type(2.0) * __urng() - 1.0; 1376169691Skan __r2 = __x * __x + __y * __y; 1377169691Skan } 1378169691Skan while (__r2 > 1.0 || __r2 == 0.0); 1379169691Skan 1380169691Skan const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 1381169691Skan _M_saved = __x * __mult; 1382169691Skan _M_saved_available = true; 1383169691Skan __ret = __y * __mult; 1384169691Skan } 1385169691Skan 1386169691Skan __ret = __ret * _M_sigma + _M_mean; 1387169691Skan return __ret; 1388169691Skan } 1389169691Skan 1390169691Skan template<typename _RealType, typename _CharT, typename _Traits> 1391169691Skan std::basic_ostream<_CharT, _Traits>& 1392169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1393169691Skan const normal_distribution<_RealType>& __x) 1394169691Skan { 1395169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1396169691Skan typedef typename __ostream_type::ios_base __ios_base; 1397169691Skan 1398169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 1399169691Skan const _CharT __fill = __os.fill(); 1400169691Skan const std::streamsize __precision = __os.precision(); 1401169691Skan const _CharT __space = __os.widen(' '); 1402169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 1403169691Skan __os.fill(__space); 1404169691Skan __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1405169691Skan 1406169691Skan __os << __x._M_saved_available << __space 1407169691Skan << __x.mean() << __space 1408169691Skan << __x.sigma(); 1409169691Skan if (__x._M_saved_available) 1410169691Skan __os << __space << __x._M_saved; 1411169691Skan 1412169691Skan __os.flags(__flags); 1413169691Skan __os.fill(__fill); 1414169691Skan __os.precision(__precision); 1415169691Skan return __os; 1416169691Skan } 1417169691Skan 1418169691Skan template<typename _RealType, typename _CharT, typename _Traits> 1419169691Skan std::basic_istream<_CharT, _Traits>& 1420169691Skan operator>>(std::basic_istream<_CharT, _Traits>& __is, 1421169691Skan normal_distribution<_RealType>& __x) 1422169691Skan { 1423169691Skan typedef std::basic_istream<_CharT, _Traits> __istream_type; 1424169691Skan typedef typename __istream_type::ios_base __ios_base; 1425169691Skan 1426169691Skan const typename __ios_base::fmtflags __flags = __is.flags(); 1427169691Skan __is.flags(__ios_base::dec | __ios_base::skipws); 1428169691Skan 1429169691Skan __is >> __x._M_saved_available >> __x._M_mean 1430169691Skan >> __x._M_sigma; 1431169691Skan if (__x._M_saved_available) 1432169691Skan __is >> __x._M_saved; 1433169691Skan 1434169691Skan __is.flags(__flags); 1435169691Skan return __is; 1436169691Skan } 1437169691Skan 1438169691Skan 1439169691Skan template<typename _RealType> 1440169691Skan void 1441169691Skan gamma_distribution<_RealType>:: 1442169691Skan _M_initialize() 1443169691Skan { 1444169691Skan if (_M_alpha >= 1) 1445169691Skan _M_l_d = std::sqrt(2 * _M_alpha - 1); 1446169691Skan else 1447169691Skan _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha)) 1448169691Skan * (1 - _M_alpha)); 1449169691Skan } 1450169691Skan 1451169691Skan /** 1452169691Skan * Cheng's rejection algorithm GB for alpha >= 1 and a modification 1453169691Skan * of Vaduva's rejection from Weibull algorithm due to Devroye for 1454169691Skan * alpha < 1. 1455169691Skan * 1456169691Skan * References: 1457169691Skan * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral 1458169691Skan * Shape Parameter." Applied Statistics, 26, 71-75, 1977. 1459169691Skan * 1460169691Skan * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection 1461169691Skan * and Composition Procedures." Math. Operationsforschung and Statistik, 1462169691Skan * Series in Statistics, 8, 545-576, 1977. 1463169691Skan * 1464169691Skan * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, 1465169691Skan * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!). 1466169691Skan */ 1467169691Skan template<typename _RealType> 1468169691Skan template<class _UniformRandomNumberGenerator> 1469169691Skan typename gamma_distribution<_RealType>::result_type 1470169691Skan gamma_distribution<_RealType>:: 1471169691Skan operator()(_UniformRandomNumberGenerator& __urng) 1472169691Skan { 1473169691Skan result_type __x; 1474169691Skan 1475169691Skan bool __reject; 1476169691Skan if (_M_alpha >= 1) 1477169691Skan { 1478169691Skan // alpha - log(4) 1479169691Skan const result_type __b = _M_alpha 1480169691Skan - result_type(1.3862943611198906188344642429163531L); 1481169691Skan const result_type __c = _M_alpha + _M_l_d; 1482169691Skan const result_type __1l = 1 / _M_l_d; 1483169691Skan 1484169691Skan // 1 + log(9 / 2) 1485169691Skan const result_type __k = 2.5040773967762740733732583523868748L; 1486169691Skan 1487169691Skan do 1488169691Skan { 1489169691Skan const result_type __u = __urng(); 1490169691Skan const result_type __v = __urng(); 1491169691Skan 1492169691Skan const result_type __y = __1l * std::log(__v / (1 - __v)); 1493169691Skan __x = _M_alpha * std::exp(__y); 1494169691Skan 1495169691Skan const result_type __z = __u * __v * __v; 1496169691Skan const result_type __r = __b + __c * __y - __x; 1497169691Skan 1498169691Skan __reject = __r < result_type(4.5) * __z - __k; 1499169691Skan if (__reject) 1500169691Skan __reject = __r < std::log(__z); 1501169691Skan } 1502169691Skan while (__reject); 1503169691Skan } 1504169691Skan else 1505169691Skan { 1506169691Skan const result_type __c = 1 / _M_alpha; 1507169691Skan 1508169691Skan do 1509169691Skan { 1510169691Skan const result_type __z = -std::log(__urng()); 1511169691Skan const result_type __e = -std::log(__urng()); 1512169691Skan 1513169691Skan __x = std::pow(__z, __c); 1514169691Skan 1515169691Skan __reject = __z + __e < _M_l_d + __x; 1516169691Skan } 1517169691Skan while (__reject); 1518169691Skan } 1519169691Skan 1520169691Skan return __x; 1521169691Skan } 1522169691Skan 1523169691Skan template<typename _RealType, typename _CharT, typename _Traits> 1524169691Skan std::basic_ostream<_CharT, _Traits>& 1525169691Skan operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1526169691Skan const gamma_distribution<_RealType>& __x) 1527169691Skan { 1528169691Skan typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1529169691Skan typedef typename __ostream_type::ios_base __ios_base; 1530169691Skan 1531169691Skan const typename __ios_base::fmtflags __flags = __os.flags(); 1532169691Skan const _CharT __fill = __os.fill(); 1533169691Skan const std::streamsize __precision = __os.precision(); 1534169691Skan __os.flags(__ios_base::scientific | __ios_base::left); 1535169691Skan __os.fill(__os.widen(' ')); 1536169691Skan __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10); 1537169691Skan 1538169691Skan __os << __x.alpha(); 1539169691Skan 1540169691Skan __os.flags(__flags); 1541169691Skan __os.fill(__fill); 1542169691Skan __os.precision(__precision); 1543169691Skan return __os; 1544169691Skan } 1545169691Skan 1546169691Skan_GLIBCXX_END_NAMESPACE 1547169691Skan} 1548