random.tcc

Go to the documentation of this file.
00001 // random number generation (out of line) -*- C++ -*-
00002 
00003 // Copyright (C) 2006, 2007 Free Software Foundation, Inc.
00004 //
00005 // This file is part of the GNU ISO C++ Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the
00007 // terms of the GNU General Public License as published by the
00008 // Free Software Foundation; either version 2, or (at your option)
00009 // any later version.
00010 
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 
00016 // You should have received a copy of the GNU General Public License along
00017 // with this library; see the file COPYING.  If not, write to the Free
00018 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
00019 // USA.
00020 
00021 // As a special exception, you may use this file as part of a free software
00022 // library without restriction.  Specifically, if other files instantiate
00023 // templates or use macros or inline functions from this file, or you compile
00024 // this file and link it with other files to produce an executable, this
00025 // file does not by itself cause the resulting executable to be covered by
00026 // the GNU General Public License.  This exception does not however
00027 // invalidate any other reasons why the executable file might be covered by
00028 // the GNU General Public License.
00029 
00030 /** @file tr1/random.tcc
00031  *  This is a TR1 C++ Library header. 
00032  */
00033 
00034 namespace std
00035 {
00036 _GLIBCXX_BEGIN_NAMESPACE(tr1)
00037 
00038   /*
00039    * (Further) implementation-space details.
00040    */
00041   namespace __detail
00042   {
00043     // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
00044     // integer overflow.
00045     //
00046     // Because a and c are compile-time integral constants the compiler kindly
00047     // elides any unreachable paths.
00048     //
00049     // Preconditions:  a > 0, m > 0.
00050     //
00051     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
00052       struct _Mod
00053       {
00054     static _Tp
00055     __calc(_Tp __x)
00056     {
00057       if (__a == 1)
00058         __x %= __m;
00059       else
00060         {
00061           static const _Tp __q = __m / __a;
00062           static const _Tp __r = __m % __a;
00063           
00064           _Tp __t1 = __a * (__x % __q);
00065           _Tp __t2 = __r * (__x / __q);
00066           if (__t1 >= __t2)
00067         __x = __t1 - __t2;
00068           else
00069         __x = __m - __t2 + __t1;
00070         }
00071 
00072       if (__c != 0)
00073         {
00074           const _Tp __d = __m - __x;
00075           if (__d > __c)
00076         __x += __c;
00077           else
00078         __x = __c - __d;
00079         }
00080       return __x;
00081     }
00082       };
00083 
00084     // Special case for m == 0 -- use unsigned integer overflow as modulo
00085     // operator.
00086     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
00087       struct _Mod<_Tp, __a, __c, __m, true>
00088       {
00089     static _Tp
00090     __calc(_Tp __x)
00091     { return __a * __x + __c; }
00092       };
00093   } // namespace __detail
00094 
00095   /**
00096    * Seeds the LCR with integral value @p __x0, adjusted so that the 
00097    * ring identity is never a member of the convergence set.
00098    */
00099   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00100     void
00101     linear_congruential<_UIntType, __a, __c, __m>::
00102     seed(unsigned long __x0)
00103     {
00104       if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
00105       && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
00106     _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
00107       else
00108     _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
00109     }
00110 
00111   /**
00112    * Seeds the LCR engine with a value generated by @p __g.
00113    */
00114   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00115     template<class _Gen>
00116       void
00117       linear_congruential<_UIntType, __a, __c, __m>::
00118       seed(_Gen& __g, false_type)
00119       {
00120     _UIntType __x0 = __g();
00121     if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
00122         && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
00123       _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
00124     else
00125       _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
00126       }
00127 
00128   /**
00129    * Gets the next generated value in sequence.
00130    */
00131   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00132     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
00133     linear_congruential<_UIntType, __a, __c, __m>::
00134     operator()()
00135     {
00136       _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
00137       return _M_x;
00138     }
00139 
00140   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
00141        typename _CharT, typename _Traits>
00142     std::basic_ostream<_CharT, _Traits>&
00143     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00144            const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
00145     {
00146       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00147       typedef typename __ostream_type::ios_base    __ios_base;
00148 
00149       const typename __ios_base::fmtflags __flags = __os.flags();
00150       const _CharT __fill = __os.fill();
00151       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00152       __os.fill(__os.widen(' '));
00153 
00154       __os << __lcr._M_x;
00155 
00156       __os.flags(__flags);
00157       __os.fill(__fill);
00158       return __os;
00159     }
00160 
00161   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
00162        typename _CharT, typename _Traits>
00163     std::basic_istream<_CharT, _Traits>&
00164     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00165            linear_congruential<_UIntType, __a, __c, __m>& __lcr)
00166     {
00167       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00168       typedef typename __istream_type::ios_base    __ios_base;
00169 
00170       const typename __ios_base::fmtflags __flags = __is.flags();
00171       __is.flags(__ios_base::dec);
00172 
00173       __is >> __lcr._M_x;
00174 
00175       __is.flags(__flags);
00176       return __is;
00177     } 
00178 
00179 
00180   template<class _UIntType, int __w, int __n, int __m, int __r,
00181        _UIntType __a, int __u, int __s,
00182        _UIntType __b, int __t, _UIntType __c, int __l>
00183     void
00184     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
00185              __b, __t, __c, __l>::
00186     seed(unsigned long __value)
00187     {
00188       _M_x[0] = __detail::__mod<_UIntType, 1, 0,
00189     __detail::_Shift<_UIntType, __w>::__value>(__value);
00190 
00191       for (int __i = 1; __i < state_size; ++__i)
00192     {
00193       _UIntType __x = _M_x[__i - 1];
00194       __x ^= __x >> (__w - 2);
00195       __x *= 1812433253ul;
00196       __x += __i;
00197       _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
00198         __detail::_Shift<_UIntType, __w>::__value>(__x);      
00199     }
00200       _M_p = state_size;
00201     }
00202 
00203   template<class _UIntType, int __w, int __n, int __m, int __r,
00204        _UIntType __a, int __u, int __s,
00205        _UIntType __b, int __t, _UIntType __c, int __l>
00206     template<class _Gen>
00207       void
00208       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
00209                __b, __t, __c, __l>::
00210       seed(_Gen& __gen, false_type)
00211       {
00212     for (int __i = 0; __i < state_size; ++__i)
00213       _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
00214         __detail::_Shift<_UIntType, __w>::__value>(__gen());
00215     _M_p = state_size;
00216       }
00217 
00218   template<class _UIntType, int __w, int __n, int __m, int __r,
00219        _UIntType __a, int __u, int __s,
00220        _UIntType __b, int __t, _UIntType __c, int __l>
00221     typename
00222     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
00223              __b, __t, __c, __l>::result_type
00224     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
00225              __b, __t, __c, __l>::
00226     operator()()
00227     {
00228       // Reload the vector - cost is O(n) amortized over n calls.
00229       if (_M_p >= state_size)
00230     {
00231       const _UIntType __upper_mask = (~_UIntType()) << __r;
00232       const _UIntType __lower_mask = ~__upper_mask;
00233 
00234       for (int __k = 0; __k < (__n - __m); ++__k)
00235         {
00236           _UIntType __y = ((_M_x[__k] & __upper_mask)
00237                    | (_M_x[__k + 1] & __lower_mask));
00238           _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
00239                ^ ((__y & 0x01) ? __a : 0));
00240         }
00241 
00242       for (int __k = (__n - __m); __k < (__n - 1); ++__k)
00243         {
00244           _UIntType __y = ((_M_x[__k] & __upper_mask)
00245                    | (_M_x[__k + 1] & __lower_mask));
00246           _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
00247                ^ ((__y & 0x01) ? __a : 0));
00248         }
00249 
00250       _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
00251                | (_M_x[0] & __lower_mask));
00252       _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
00253                ^ ((__y & 0x01) ? __a : 0));
00254       _M_p = 0;
00255     }
00256 
00257       // Calculate o(x(i)).
00258       result_type __z = _M_x[_M_p++];
00259       __z ^= (__z >> __u);
00260       __z ^= (__z << __s) & __b;
00261       __z ^= (__z << __t) & __c;
00262       __z ^= (__z >> __l);
00263 
00264       return __z;
00265     }
00266 
00267   template<class _UIntType, int __w, int __n, int __m, int __r,
00268        _UIntType __a, int __u, int __s, _UIntType __b, int __t,
00269        _UIntType __c, int __l,
00270        typename _CharT, typename _Traits>
00271     std::basic_ostream<_CharT, _Traits>&
00272     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00273            const mersenne_twister<_UIntType, __w, __n, __m,
00274            __r, __a, __u, __s, __b, __t, __c, __l>& __x)
00275     {
00276       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00277       typedef typename __ostream_type::ios_base    __ios_base;
00278 
00279       const typename __ios_base::fmtflags __flags = __os.flags();
00280       const _CharT __fill = __os.fill();
00281       const _CharT __space = __os.widen(' ');
00282       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00283       __os.fill(__space);
00284 
00285       for (int __i = 0; __i < __n - 1; ++__i)
00286     __os << __x._M_x[__i] << __space;
00287       __os << __x._M_x[__n - 1];
00288 
00289       __os.flags(__flags);
00290       __os.fill(__fill);
00291       return __os;
00292     }
00293 
00294   template<class _UIntType, int __w, int __n, int __m, int __r,
00295        _UIntType __a, int __u, int __s, _UIntType __b, int __t,
00296        _UIntType __c, int __l,
00297        typename _CharT, typename _Traits>
00298     std::basic_istream<_CharT, _Traits>&
00299     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00300            mersenne_twister<_UIntType, __w, __n, __m,
00301            __r, __a, __u, __s, __b, __t, __c, __l>& __x)
00302     {
00303       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00304       typedef typename __istream_type::ios_base    __ios_base;
00305 
00306       const typename __ios_base::fmtflags __flags = __is.flags();
00307       __is.flags(__ios_base::dec | __ios_base::skipws);
00308 
00309       for (int __i = 0; __i < __n; ++__i)
00310     __is >> __x._M_x[__i];
00311 
00312       __is.flags(__flags);
00313       return __is;
00314     }
00315 
00316 
00317   template<typename _IntType, _IntType __m, int __s, int __r>
00318     void
00319     subtract_with_carry<_IntType, __m, __s, __r>::
00320     seed(unsigned long __value)
00321     {
00322       if (__value == 0)
00323     __value = 19780503;
00324 
00325       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
00326     __lcg(__value);
00327 
00328       for (int __i = 0; __i < long_lag; ++__i)
00329     _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());
00330 
00331       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
00332       _M_p = 0;
00333     }
00334 
00335   template<typename _IntType, _IntType __m, int __s, int __r>
00336     template<class _Gen>
00337       void
00338       subtract_with_carry<_IntType, __m, __s, __r>::
00339       seed(_Gen& __gen, false_type)
00340       {
00341     const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;
00342 
00343     for (int __i = 0; __i < long_lag; ++__i)
00344       {
00345         _UIntType __tmp = 0;
00346         _UIntType __factor = 1;
00347         for (int __j = 0; __j < __n; ++__j)
00348           {
00349         __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
00350                  (__gen()) * __factor;
00351         __factor *= __detail::_Shift<_UIntType, 32>::__value;
00352           }
00353         _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
00354       }
00355     _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
00356     _M_p = 0;
00357       }
00358 
00359   template<typename _IntType, _IntType __m, int __s, int __r>
00360     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
00361     subtract_with_carry<_IntType, __m, __s, __r>::
00362     operator()()
00363     {
00364       // Derive short lag index from current index.
00365       int __ps = _M_p - short_lag;
00366       if (__ps < 0)
00367     __ps += long_lag;
00368 
00369       // Calculate new x(i) without overflow or division.
00370       // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
00371       // cannot overflow.
00372       _UIntType __xi;
00373       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
00374     {
00375       __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
00376       _M_carry = 0;
00377     }
00378       else
00379     {
00380       __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
00381       _M_carry = 1;
00382     }
00383       _M_x[_M_p] = __xi;
00384 
00385       // Adjust current index to loop around in ring buffer.
00386       if (++_M_p >= long_lag)
00387     _M_p = 0;
00388 
00389       return __xi;
00390     }
00391 
00392   template<typename _IntType, _IntType __m, int __s, int __r,
00393        typename _CharT, typename _Traits>
00394     std::basic_ostream<_CharT, _Traits>&
00395     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00396            const subtract_with_carry<_IntType, __m, __s, __r>& __x)
00397     {
00398       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00399       typedef typename __ostream_type::ios_base    __ios_base;
00400 
00401       const typename __ios_base::fmtflags __flags = __os.flags();
00402       const _CharT __fill = __os.fill();
00403       const _CharT __space = __os.widen(' ');
00404       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00405       __os.fill(__space);
00406 
00407       for (int __i = 0; __i < __r; ++__i)
00408     __os << __x._M_x[__i] << __space;
00409       __os << __x._M_carry;
00410 
00411       __os.flags(__flags);
00412       __os.fill(__fill);
00413       return __os;
00414     }
00415 
00416   template<typename _IntType, _IntType __m, int __s, int __r,
00417        typename _CharT, typename _Traits>
00418     std::basic_istream<_CharT, _Traits>&
00419     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00420            subtract_with_carry<_IntType, __m, __s, __r>& __x)
00421     {
00422       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
00423       typedef typename __istream_type::ios_base    __ios_base;
00424 
00425       const typename __ios_base::fmtflags __flags = __is.flags();
00426       __is.flags(__ios_base::dec | __ios_base::skipws);
00427 
00428       for (int __i = 0; __i < __r; ++__i)
00429     __is >> __x._M_x[__i];
00430       __is >> __x._M_carry;
00431 
00432       __is.flags(__flags);
00433       return __is;
00434     }
00435 
00436 
00437   template<typename _RealType, int __w, int __s, int __r>
00438     void
00439     subtract_with_carry_01<_RealType, __w, __s, __r>::
00440     _M_initialize_npows()
00441     {
00442       for (int __j = 0; __j < __n; ++__j)
00443 #if _GLIBCXX_USE_C99_MATH_TR1
00444     _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32);
00445 #else
00446         _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
00447 #endif
00448     }
00449 
00450   template<typename _RealType, int __w, int __s, int __r>
00451     void
00452     subtract_with_carry_01<_RealType, __w, __s, __r>::
00453     seed(unsigned long __value)
00454     {
00455       if (__value == 0)
00456     __value = 19780503;
00457 
00458       // _GLIBCXX_RESOLVE_LIB_DEFECTS
00459       // 512. Seeding subtract_with_carry_01 from a single unsigned long.
00460       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
00461     __lcg(__value);
00462 
00463       this->seed(__lcg);
00464     }
00465 
00466   template<typename _RealType, int __w, int __s, int __r>
00467     template<class _Gen>
00468       void
00469       subtract_with_carry_01<_RealType, __w, __s, __r>::
00470       seed(_Gen& __gen, false_type)
00471       {
00472     for (int __i = 0; __i < long_lag; ++__i)
00473       {
00474         for (int __j = 0; __j < __n - 1; ++__j)
00475           _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
00476         _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
00477           __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
00478       }
00479 
00480     _M_carry = 1;
00481     for (int __j = 0; __j < __n; ++__j)
00482       if (_M_x[long_lag - 1][__j] != 0)
00483         {
00484           _M_carry = 0;
00485           break;
00486         }
00487 
00488     _M_p = 0;
00489       }
00490 
00491   template<typename _RealType, int __w, int __s, int __r>
00492     typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
00493     subtract_with_carry_01<_RealType, __w, __s, __r>::
00494     operator()()
00495     {
00496       // Derive short lag index from current index.
00497       int __ps = _M_p - short_lag;
00498       if (__ps < 0)
00499     __ps += long_lag;
00500 
00501       _UInt32Type __new_carry;
00502       for (int __j = 0; __j < __n - 1; ++__j)
00503     {
00504       if (_M_x[__ps][__j] > _M_x[_M_p][__j]
00505           || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
00506         __new_carry = 0;
00507       else
00508         __new_carry = 1;
00509 
00510       _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
00511       _M_carry = __new_carry;
00512     }
00513 
00514       if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
00515       || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
00516     __new_carry = 0;
00517       else
00518     __new_carry = 1;
00519       
00520       _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
00521     __detail::_Shift<_UInt32Type, __w % 32>::__value>
00522     (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
00523       _M_carry = __new_carry;
00524 
00525       result_type __ret = 0.0;
00526       for (int __j = 0; __j < __n; ++__j)
00527     __ret += _M_x[_M_p][__j] * _M_npows[__j];
00528 
00529       // Adjust current index to loop around in ring buffer.
00530       if (++_M_p >= long_lag)
00531     _M_p = 0;
00532 
00533       return __ret;
00534     }
00535 
00536   template<typename _RealType, int __w, int __s, int __r,
00537        typename _CharT, typename _Traits>
00538     std::basic_ostream<_CharT, _Traits>&
00539     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00540            const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
00541     {
00542       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00543       typedef typename __ostream_type::ios_base    __ios_base;
00544 
00545       const typename __ios_base::fmtflags __flags = __os.flags();
00546       const _CharT __fill = __os.fill();
00547       const _CharT __space = __os.widen(' ');
00548       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00549       __os.fill(__space);
00550 
00551       for (int __i = 0; __i < __r; ++__i)
00552     for (int __j = 0; __j < __x.__n; ++__j)
00553       __os << __x._M_x[__i][__j] << __space;
00554       __os << __x._M_carry;
00555 
00556       __os.flags(__flags);
00557       __os.fill(__fill);
00558       return __os;
00559     }
00560 
00561   template<typename _RealType, int __w, int __s, int __r,
00562        typename _CharT, typename _Traits>
00563     std::basic_istream<_CharT, _Traits>&
00564     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00565            subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
00566     {
00567       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00568       typedef typename __istream_type::ios_base    __ios_base;
00569 
00570       const typename __ios_base::fmtflags __flags = __is.flags();
00571       __is.flags(__ios_base::dec | __ios_base::skipws);
00572 
00573       for (int __i = 0; __i < __r; ++__i)
00574     for (int __j = 0; __j < __x.__n; ++__j)
00575       __is >> __x._M_x[__i][__j];
00576       __is >> __x._M_carry;
00577 
00578       __is.flags(__flags);
00579       return __is;
00580     }
00581 
00582 
00583   template<class _UniformRandomNumberGenerator, int __p, int __r>
00584     typename discard_block<_UniformRandomNumberGenerator,
00585                __p, __r>::result_type
00586     discard_block<_UniformRandomNumberGenerator, __p, __r>::
00587     operator()()
00588     {
00589       if (_M_n >= used_block)
00590     {
00591       while (_M_n < block_size)
00592         {
00593           _M_b();
00594           ++_M_n;
00595         }
00596       _M_n = 0;
00597     }
00598       ++_M_n;
00599       return _M_b();
00600     }
00601 
00602   template<class _UniformRandomNumberGenerator, int __p, int __r,
00603        typename _CharT, typename _Traits>
00604     std::basic_ostream<_CharT, _Traits>&
00605     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00606            const discard_block<_UniformRandomNumberGenerator,
00607            __p, __r>& __x)
00608     {
00609       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00610       typedef typename __ostream_type::ios_base    __ios_base;
00611 
00612       const typename __ios_base::fmtflags __flags = __os.flags();
00613       const _CharT __fill = __os.fill();
00614       const _CharT __space = __os.widen(' ');
00615       __os.flags(__ios_base::dec | __ios_base::fixed
00616          | __ios_base::left);
00617       __os.fill(__space);
00618 
00619       __os << __x._M_b << __space << __x._M_n;
00620 
00621       __os.flags(__flags);
00622       __os.fill(__fill);
00623       return __os;
00624     }
00625 
00626   template<class _UniformRandomNumberGenerator, int __p, int __r,
00627        typename _CharT, typename _Traits>
00628     std::basic_istream<_CharT, _Traits>&
00629     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00630            discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
00631     {
00632       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00633       typedef typename __istream_type::ios_base    __ios_base;
00634 
00635       const typename __ios_base::fmtflags __flags = __is.flags();
00636       __is.flags(__ios_base::dec | __ios_base::skipws);
00637 
00638       __is >> __x._M_b >> __x._M_n;
00639 
00640       __is.flags(__flags);
00641       return __is;
00642     }
00643 
00644 
00645   template<class _UniformRandomNumberGenerator1, int __s1,
00646        class _UniformRandomNumberGenerator2, int __s2>
00647     void
00648     xor_combine<_UniformRandomNumberGenerator1, __s1,
00649         _UniformRandomNumberGenerator2, __s2>::
00650     _M_initialize_max()
00651     {
00652       const int __w = std::numeric_limits<result_type>::digits;
00653 
00654       const result_type __m1 =
00655     std::min(result_type(_M_b1.max() - _M_b1.min()),
00656          __detail::_Shift<result_type, __w - __s1>::__value - 1);
00657 
00658       const result_type __m2 =
00659     std::min(result_type(_M_b2.max() - _M_b2.min()),
00660          __detail::_Shift<result_type, __w - __s2>::__value - 1);
00661 
00662       // NB: In TR1 s1 is not required to be >= s2.
00663       if (__s1 < __s2)
00664     _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
00665       else
00666     _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
00667     }
00668 
00669   template<class _UniformRandomNumberGenerator1, int __s1,
00670        class _UniformRandomNumberGenerator2, int __s2>
00671     typename xor_combine<_UniformRandomNumberGenerator1, __s1,
00672              _UniformRandomNumberGenerator2, __s2>::result_type
00673     xor_combine<_UniformRandomNumberGenerator1, __s1,
00674         _UniformRandomNumberGenerator2, __s2>::
00675     _M_initialize_max_aux(result_type __a, result_type __b, int __d)
00676     {
00677       const result_type __two2d = result_type(1) << __d;
00678       const result_type __c = __a * __two2d;
00679 
00680       if (__a == 0 || __b < __two2d)
00681     return __c + __b;
00682 
00683       const result_type __t = std::max(__c, __b);
00684       const result_type __u = std::min(__c, __b);
00685 
00686       result_type __ub = __u;
00687       result_type __p;
00688       for (__p = 0; __ub != 1; __ub >>= 1)
00689     ++__p;
00690 
00691       const result_type __two2p = result_type(1) << __p;
00692       const result_type __k = __t / __two2p;
00693 
00694       if (__k & 1)
00695     return (__k + 1) * __two2p - 1;
00696 
00697       if (__c >= __b)
00698     return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
00699                                / __two2d,
00700                                __u % __two2p, __d);
00701       else
00702     return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
00703                                / __two2d,
00704                                __t % __two2p, __d);
00705     }
00706 
00707   template<class _UniformRandomNumberGenerator1, int __s1,
00708        class _UniformRandomNumberGenerator2, int __s2,
00709        typename _CharT, typename _Traits>
00710     std::basic_ostream<_CharT, _Traits>&
00711     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00712            const xor_combine<_UniformRandomNumberGenerator1, __s1,
00713            _UniformRandomNumberGenerator2, __s2>& __x)
00714     {
00715       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00716       typedef typename __ostream_type::ios_base    __ios_base;
00717 
00718       const typename __ios_base::fmtflags __flags = __os.flags();
00719       const _CharT __fill = __os.fill();
00720       const _CharT __space = __os.widen(' ');
00721       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00722       __os.fill(__space);
00723 
00724       __os << __x.base1() << __space << __x.base2();
00725 
00726       __os.flags(__flags);
00727       __os.fill(__fill);
00728       return __os; 
00729     }
00730 
00731   template<class _UniformRandomNumberGenerator1, int __s1,
00732        class _UniformRandomNumberGenerator2, int __s2,
00733        typename _CharT, typename _Traits>
00734     std::basic_istream<_CharT, _Traits>&
00735     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00736            xor_combine<_UniformRandomNumberGenerator1, __s1,
00737            _UniformRandomNumberGenerator2, __s2>& __x)
00738     {
00739       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00740       typedef typename __istream_type::ios_base    __ios_base;
00741 
00742       const typename __ios_base::fmtflags __flags = __is.flags();
00743       __is.flags(__ios_base::skipws);
00744 
00745       __is >> __x._M_b1 >> __x._M_b2;
00746 
00747       __is.flags(__flags);
00748       return __is;
00749     }
00750 
00751 
00752   template<typename _IntType, typename _CharT, typename _Traits>
00753     std::basic_ostream<_CharT, _Traits>&
00754     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00755            const uniform_int<_IntType>& __x)
00756     {
00757       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00758       typedef typename __ostream_type::ios_base    __ios_base;
00759 
00760       const typename __ios_base::fmtflags __flags = __os.flags();
00761       const _CharT __fill = __os.fill();
00762       const _CharT __space = __os.widen(' ');
00763       __os.flags(__ios_base::scientific | __ios_base::left);
00764       __os.fill(__space);
00765 
00766       __os << __x.min() << __space << __x.max();
00767 
00768       __os.flags(__flags);
00769       __os.fill(__fill);
00770       return __os;
00771     }
00772 
00773   template<typename _IntType, typename _CharT, typename _Traits>
00774     std::basic_istream<_CharT, _Traits>&
00775     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00776            uniform_int<_IntType>& __x)
00777     {
00778       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00779       typedef typename __istream_type::ios_base    __ios_base;
00780 
00781       const typename __ios_base::fmtflags __flags = __is.flags();
00782       __is.flags(__ios_base::dec | __ios_base::skipws);
00783 
00784       __is >> __x._M_min >> __x._M_max;
00785 
00786       __is.flags(__flags);
00787       return __is;
00788     }
00789 
00790   
00791   template<typename _CharT, typename _Traits>
00792     std::basic_ostream<_CharT, _Traits>&
00793     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00794            const bernoulli_distribution& __x)
00795     {
00796       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00797       typedef typename __ostream_type::ios_base    __ios_base;
00798 
00799       const typename __ios_base::fmtflags __flags = __os.flags();
00800       const _CharT __fill = __os.fill();
00801       const std::streamsize __precision = __os.precision();
00802       __os.flags(__ios_base::scientific | __ios_base::left);
00803       __os.fill(__os.widen(' '));
00804       __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);
00805 
00806       __os << __x.p();
00807 
00808       __os.flags(__flags);
00809       __os.fill(__fill);
00810       __os.precision(__precision);
00811       return __os;
00812     }
00813 
00814 
00815   template<typename _IntType, typename _RealType>
00816     template<class _UniformRandomNumberGenerator>
00817       typename geometric_distribution<_IntType, _RealType>::result_type
00818       geometric_distribution<_IntType, _RealType>::
00819       operator()(_UniformRandomNumberGenerator& __urng)
00820       {
00821     // About the epsilon thing see this thread:
00822         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
00823     const _RealType __naf =
00824       (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
00825     // The largest _RealType convertible to _IntType.
00826     const _RealType __thr =
00827       std::numeric_limits<_IntType>::max() + __naf;
00828 
00829     _RealType __cand;
00830     do
00831       __cand = std::ceil(std::log(__urng()) / _M_log_p);
00832     while (__cand >= __thr);
00833 
00834     return result_type(__cand + __naf);
00835       }
00836 
00837   template<typename _IntType, typename _RealType,
00838        typename _CharT, typename _Traits>
00839     std::basic_ostream<_CharT, _Traits>&
00840     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00841            const geometric_distribution<_IntType, _RealType>& __x)
00842     {
00843       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00844       typedef typename __ostream_type::ios_base    __ios_base;
00845 
00846       const typename __ios_base::fmtflags __flags = __os.flags();
00847       const _CharT __fill = __os.fill();
00848       const std::streamsize __precision = __os.precision();
00849       __os.flags(__ios_base::scientific | __ios_base::left);
00850       __os.fill(__os.widen(' '));
00851       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
00852 
00853       __os << __x.p();
00854 
00855       __os.flags(__flags);
00856       __os.fill(__fill);
00857       __os.precision(__precision);
00858       return __os;
00859     }
00860 
00861 
00862   template<typename _IntType, typename _RealType>
00863     void
00864     poisson_distribution<_IntType, _RealType>::
00865     _M_initialize()
00866     {
00867 #if _GLIBCXX_USE_C99_MATH_TR1
00868       if (_M_mean >= 12)
00869     {
00870       const _RealType __m = std::floor(_M_mean);
00871       _M_lm_thr = std::log(_M_mean);
00872       _M_lfm = std::tr1::lgamma(__m + 1);
00873       _M_sm = std::sqrt(__m);
00874 
00875       const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
00876       const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
00877                                   / __pi_4));
00878       _M_d = std::tr1::round(std::max(_RealType(6),
00879                       std::min(__m, __dx)));
00880       const _RealType __cx = 2 * __m + _M_d;
00881       _M_scx = std::sqrt(__cx / 2);
00882       _M_1cx = 1 / __cx;
00883 
00884       _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
00885       _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
00886     }
00887       else
00888 #endif
00889     _M_lm_thr = std::exp(-_M_mean);
00890       }
00891 
00892   /**
00893    * A rejection algorithm when mean >= 12 and a simple method based
00894    * upon the multiplication of uniform random variates otherwise.
00895    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
00896    * is defined.
00897    *
00898    * Reference:
00899    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
00900    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
00901    */
00902   template<typename _IntType, typename _RealType>
00903     template<class _UniformRandomNumberGenerator>
00904       typename poisson_distribution<_IntType, _RealType>::result_type
00905       poisson_distribution<_IntType, _RealType>::
00906       operator()(_UniformRandomNumberGenerator& __urng)
00907       {
00908 #if _GLIBCXX_USE_C99_MATH_TR1
00909     if (_M_mean >= 12)
00910       {
00911         _RealType __x;
00912 
00913         // See comments above...
00914         const _RealType __naf =
00915           (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
00916         const _RealType __thr =
00917           std::numeric_limits<_IntType>::max() + __naf;
00918 
00919         const _RealType __m = std::floor(_M_mean);
00920         // sqrt(pi / 2)
00921         const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
00922         const _RealType __c1 = _M_sm * __spi_2;
00923         const _RealType __c2 = _M_c2b + __c1; 
00924         const _RealType __c3 = __c2 + 1;
00925         const _RealType __c4 = __c3 + 1;
00926         // e^(1 / 78)
00927         const _RealType __e178 = 1.0129030479320018583185514777512983L;
00928         const _RealType __c5 = __c4 + __e178;
00929         const _RealType __c = _M_cb + __c5;
00930         const _RealType __2cx = 2 * (2 * __m + _M_d);
00931 
00932         bool __reject = true;
00933         do
00934           {
00935         const _RealType __u = __c * __urng();
00936         const _RealType __e = -std::log(__urng());
00937 
00938         _RealType __w = 0.0;
00939         
00940         if (__u <= __c1)
00941           {
00942             const _RealType __n = _M_nd(__urng);
00943             const _RealType __y = -std::abs(__n) * _M_sm - 1;
00944             __x = std::floor(__y);
00945             __w = -__n * __n / 2;
00946             if (__x < -__m)
00947               continue;
00948           }
00949         else if (__u <= __c2)
00950           {
00951             const _RealType __n = _M_nd(__urng);
00952             const _RealType __y = 1 + std::abs(__n) * _M_scx;
00953             __x = std::ceil(__y);
00954             __w = __y * (2 - __y) * _M_1cx;
00955             if (__x > _M_d)
00956               continue;
00957           }
00958         else if (__u <= __c3)
00959           // NB: This case not in the book, nor in the Errata,
00960           // but should be ok...
00961           __x = -1;
00962         else if (__u <= __c4)
00963           __x = 0;
00964         else if (__u <= __c5)
00965           __x = 1;
00966         else
00967           {
00968             const _RealType __v = -std::log(__urng());
00969             const _RealType __y = _M_d + __v * __2cx / _M_d;
00970             __x = std::ceil(__y);
00971             __w = -_M_d * _M_1cx * (1 + __y / 2);
00972           }
00973 
00974         __reject = (__w - __e - __x * _M_lm_thr
00975                 > _M_lfm - std::tr1::lgamma(__x + __m + 1));
00976 
00977         __reject |= __x + __m >= __thr;
00978 
00979           } while (__reject);
00980 
00981         return result_type(__x + __m + __naf);
00982       }
00983     else
00984 #endif
00985       {
00986         _IntType     __x = 0;
00987         _RealType __prod = 1.0;
00988 
00989         do
00990           {
00991         __prod *= __urng();
00992         __x += 1;
00993           }
00994         while (__prod > _M_lm_thr);
00995 
00996         return __x - 1;
00997       }
00998       }
00999 
01000   template<typename _IntType, typename _RealType,
01001        typename _CharT, typename _Traits>
01002     std::basic_ostream<_CharT, _Traits>&
01003     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01004            const poisson_distribution<_IntType, _RealType>& __x)
01005     {
01006       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01007       typedef typename __ostream_type::ios_base    __ios_base;
01008 
01009       const typename __ios_base::fmtflags __flags = __os.flags();
01010       const _CharT __fill = __os.fill();
01011       const std::streamsize __precision = __os.precision();
01012       const _CharT __space = __os.widen(' ');
01013       __os.flags(__ios_base::scientific | __ios_base::left);
01014       __os.fill(__space);
01015       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
01016 
01017       __os << __x.mean() << __space << __x._M_nd;
01018 
01019       __os.flags(__flags);
01020       __os.fill(__fill);
01021       __os.precision(__precision);
01022       return __os;
01023     }
01024 
01025   template<typename _IntType, typename _RealType,
01026        typename _CharT, typename _Traits>
01027     std::basic_istream<_CharT, _Traits>&
01028     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01029            poisson_distribution<_IntType, _RealType>& __x)
01030     {
01031       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01032       typedef typename __istream_type::ios_base    __ios_base;
01033 
01034       const typename __ios_base::fmtflags __flags = __is.flags();
01035       __is.flags(__ios_base::skipws);
01036 
01037       __is >> __x._M_mean >> __x._M_nd;
01038       __x._M_initialize();
01039 
01040       __is.flags(__flags);
01041       return __is;
01042     }
01043 
01044 
01045   template<typename _IntType, typename _RealType>
01046     void
01047     binomial_distribution<_IntType, _RealType>::
01048     _M_initialize()
01049     {
01050       const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
01051 
01052       _M_easy = true;
01053 
01054 #if _GLIBCXX_USE_C99_MATH_TR1
01055       if (_M_t * __p12 >= 8)
01056     {
01057       _M_easy = false;
01058       const _RealType __np = std::floor(_M_t * __p12);
01059       const _RealType __pa = __np / _M_t;
01060       const _RealType __1p = 1 - __pa;
01061       
01062       const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
01063       const _RealType __d1x =
01064         std::sqrt(__np * __1p * std::log(32 * __np
01065                          / (81 * __pi_4 * __1p)));
01066       _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));
01067       const _RealType __d2x =
01068         std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
01069                          / (__pi_4 * __pa)));
01070       _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));
01071       
01072       // sqrt(pi / 2)
01073       const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
01074       _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
01075       _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
01076       _M_c = 2 * _M_d1 / __np;
01077       _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
01078       const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
01079       const _RealType __s1s = _M_s1 * _M_s1;
01080       _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
01081                  * 2 * __s1s / _M_d1
01082                  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
01083       const _RealType __s2s = _M_s2 * _M_s2;
01084       _M_s = (_M_a123 + 2 * __s2s / _M_d2
01085           * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
01086       _M_lf = (std::tr1::lgamma(__np + 1)
01087            + std::tr1::lgamma(_M_t - __np + 1));
01088       _M_lp1p = std::log(__pa / __1p);
01089 
01090       _M_q = -std::log(1 - (__p12 - __pa) / __1p);
01091     }
01092       else
01093 #endif
01094     _M_q = -std::log(1 - __p12);
01095     }
01096 
01097   template<typename _IntType, typename _RealType>
01098     template<class _UniformRandomNumberGenerator>
01099       typename binomial_distribution<_IntType, _RealType>::result_type
01100       binomial_distribution<_IntType, _RealType>::
01101       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
01102       {
01103     _IntType    __x = 0;
01104     _RealType __sum = 0;
01105 
01106     do
01107       {
01108         const _RealType __e = -std::log(__urng());
01109         __sum += __e / (__t - __x);
01110         __x += 1;
01111       }
01112     while (__sum <= _M_q);
01113 
01114     return __x - 1;
01115       }
01116 
01117   /**
01118    * A rejection algorithm when t * p >= 8 and a simple waiting time
01119    * method - the second in the referenced book - otherwise.
01120    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
01121    * is defined.
01122    *
01123    * Reference:
01124    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
01125    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
01126    */
01127   template<typename _IntType, typename _RealType>
01128     template<class _UniformRandomNumberGenerator>
01129       typename binomial_distribution<_IntType, _RealType>::result_type
01130       binomial_distribution<_IntType, _RealType>::
01131       operator()(_UniformRandomNumberGenerator& __urng)
01132       {
01133     result_type __ret;
01134     const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
01135 
01136 #if _GLIBCXX_USE_C99_MATH_TR1
01137     if (!_M_easy)
01138       {
01139         _RealType __x;
01140 
01141         // See comments above...
01142         const _RealType __naf =
01143           (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
01144         const _RealType __thr =
01145           std::numeric_limits<_IntType>::max() + __naf;
01146 
01147         const _RealType __np = std::floor(_M_t * __p12);
01148         const _RealType __pa = __np / _M_t;
01149 
01150         // sqrt(pi / 2)
01151         const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
01152         const _RealType __a1 = _M_a1;
01153         const _RealType __a12 = __a1 + _M_s2 * __spi_2;
01154         const _RealType __a123 = _M_a123;
01155         const _RealType __s1s = _M_s1 * _M_s1;
01156         const _RealType __s2s = _M_s2 * _M_s2;
01157 
01158         bool __reject;
01159         do
01160           {
01161         const _RealType __u = _M_s * __urng();
01162 
01163         _RealType __v;
01164 
01165         if (__u <= __a1)
01166           {
01167             const _RealType __n = _M_nd(__urng);
01168             const _RealType __y = _M_s1 * std::abs(__n);
01169             __reject = __y >= _M_d1;
01170             if (!__reject)
01171               {
01172             const _RealType __e = -std::log(__urng());
01173             __x = std::floor(__y);
01174             __v = -__e - __n * __n / 2 + _M_c;
01175               }
01176           }
01177         else if (__u <= __a12)
01178           {
01179             const _RealType __n = _M_nd(__urng);
01180             const _RealType __y = _M_s2 * std::abs(__n);
01181             __reject = __y >= _M_d2;
01182             if (!__reject)
01183               {
01184             const _RealType __e = -std::log(__urng());
01185             __x = std::floor(-__y);
01186             __v = -__e - __n * __n / 2;
01187               }
01188           }
01189         else if (__u <= __a123)
01190           {
01191             const _RealType __e1 = -std::log(__urng());         
01192             const _RealType __e2 = -std::log(__urng());
01193 
01194             const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
01195             __x = std::floor(__y);
01196             __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
01197                         -__y / (2 * __s1s)));
01198             __reject = false;
01199           }
01200         else
01201           {
01202             const _RealType __e1 = -std::log(__urng());         
01203             const _RealType __e2 = -std::log(__urng());
01204 
01205             const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
01206             __x = std::floor(-__y);
01207             __v = -__e2 - _M_d2 * __y / (2 * __s2s);
01208             __reject = false;
01209           }
01210 
01211         __reject = __reject || __x < -__np || __x > _M_t - __np;
01212         if (!__reject)
01213           {
01214             const _RealType __lfx =
01215               std::tr1::lgamma(__np + __x + 1)
01216               + std::tr1::lgamma(_M_t - (__np + __x) + 1);
01217             __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
01218           }
01219 
01220         __reject |= __x + __np >= __thr;
01221           }
01222         while (__reject);
01223 
01224         __x += __np + __naf;
01225 
01226         const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); 
01227         __ret = _IntType(__x) + __z;
01228       }
01229     else
01230 #endif
01231       __ret = _M_waiting(__urng, _M_t);
01232 
01233     if (__p12 != _M_p)
01234       __ret = _M_t - __ret;
01235     return __ret;
01236       }
01237 
01238   template<typename _IntType, typename _RealType,
01239        typename _CharT, typename _Traits>
01240     std::basic_ostream<_CharT, _Traits>&
01241     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01242            const binomial_distribution<_IntType, _RealType>& __x)
01243     {
01244       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01245       typedef typename __ostream_type::ios_base    __ios_base;
01246 
01247       const typename __ios_base::fmtflags __flags = __os.flags();
01248       const _CharT __fill = __os.fill();
01249       const std::streamsize __precision = __os.precision();
01250       const _CharT __space = __os.widen(' ');
01251       __os.flags(__ios_base::scientific | __ios_base::left);
01252       __os.fill(__space);
01253       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
01254 
01255       __os << __x.t() << __space << __x.p() 
01256        << __space << __x._M_nd;
01257 
01258       __os.flags(__flags);
01259       __os.fill(__fill);
01260       __os.precision(__precision);
01261       return __os;
01262     }
01263 
01264   template<typename _IntType, typename _RealType,
01265        typename _CharT, typename _Traits>
01266     std::basic_istream<_CharT, _Traits>&
01267     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01268            binomial_distribution<_IntType, _RealType>& __x)
01269     {
01270       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01271       typedef typename __istream_type::ios_base    __ios_base;
01272 
01273       const typename __ios_base::fmtflags __flags = __is.flags();
01274       __is.flags(__ios_base::dec | __ios_base::skipws);
01275 
01276       __is >> __x._M_t >> __x._M_p >> __x._M_nd;
01277       __x._M_initialize();
01278 
01279       __is.flags(__flags);
01280       return __is;
01281     }
01282 
01283 
01284   template<typename _RealType, typename _CharT, typename _Traits>
01285     std::basic_ostream<_CharT, _Traits>&
01286     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01287            const uniform_real<_RealType>& __x)
01288     {
01289       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01290       typedef typename __ostream_type::ios_base    __ios_base;
01291 
01292       const typename __ios_base::fmtflags __flags = __os.flags();
01293       const _CharT __fill = __os.fill();
01294       const std::streamsize __precision = __os.precision();
01295       const _CharT __space = __os.widen(' ');
01296       __os.flags(__ios_base::scientific | __ios_base::left);
01297       __os.fill(__space);
01298       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
01299 
01300       __os << __x.min() << __space << __x.max();
01301 
01302       __os.flags(__flags);
01303       __os.fill(__fill);
01304       __os.precision(__precision);
01305       return __os;
01306     }
01307 
01308   template<typename _RealType, typename _CharT, typename _Traits>
01309     std::basic_istream<_CharT, _Traits>&
01310     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01311            uniform_real<_RealType>& __x)
01312     {
01313       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01314       typedef typename __istream_type::ios_base    __ios_base;
01315 
01316       const typename __ios_base::fmtflags __flags = __is.flags();
01317       __is.flags(__ios_base::skipws);
01318 
01319       __is >> __x._M_min >> __x._M_max;
01320 
01321       __is.flags(__flags);
01322       return __is;
01323     }
01324 
01325 
01326   template<typename _RealType, typename _CharT, typename _Traits>
01327     std::basic_ostream<_CharT, _Traits>&
01328     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01329            const exponential_distribution<_RealType>& __x)
01330     {
01331       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01332       typedef typename __ostream_type::ios_base    __ios_base;
01333 
01334       const typename __ios_base::fmtflags __flags = __os.flags();
01335       const _CharT __fill = __os.fill();
01336       const std::streamsize __precision = __os.precision();
01337       __os.flags(__ios_base::scientific | __ios_base::left);
01338       __os.fill(__os.widen(' '));
01339       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
01340 
01341       __os << __x.lambda();
01342 
01343       __os.flags(__flags);
01344       __os.fill(__fill);
01345       __os.precision(__precision);
01346       return __os;
01347     }
01348 
01349 
01350   /**
01351    * Polar method due to Marsaglia.
01352    *
01353    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
01354    * New York, 1986, Ch. V, Sect. 4.4.
01355    */
01356   template<typename _RealType>
01357     template<class _UniformRandomNumberGenerator>
01358       typename normal_distribution<_RealType>::result_type
01359       normal_distribution<_RealType>::
01360       operator()(_UniformRandomNumberGenerator& __urng)
01361       {
01362     result_type __ret;
01363 
01364     if (_M_saved_available)
01365       {
01366         _M_saved_available = false;
01367         __ret = _M_saved;
01368       }
01369     else
01370       {
01371         result_type __x, __y, __r2;
01372         do
01373           {
01374         __x = result_type(2.0) * __urng() - 1.0;
01375         __y = result_type(2.0) * __urng() - 1.0;
01376         __r2 = __x * __x + __y * __y;
01377           }
01378         while (__r2 > 1.0 || __r2 == 0.0);
01379 
01380         const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
01381         _M_saved = __x * __mult;
01382         _M_saved_available = true;
01383         __ret = __y * __mult;
01384       }
01385     
01386     __ret = __ret * _M_sigma + _M_mean;
01387     return __ret;
01388       }
01389 
01390   template<typename _RealType, typename _CharT, typename _Traits>
01391     std::basic_ostream<_CharT, _Traits>&
01392     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01393            const normal_distribution<_RealType>& __x)
01394     {
01395       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01396       typedef typename __ostream_type::ios_base    __ios_base;
01397 
01398       const typename __ios_base::fmtflags __flags = __os.flags();
01399       const _CharT __fill = __os.fill();
01400       const std::streamsize __precision = __os.precision();
01401       const _CharT __space = __os.widen(' ');
01402       __os.flags(__ios_base::scientific | __ios_base::left);
01403       __os.fill(__space);
01404       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
01405 
01406       __os << __x._M_saved_available << __space
01407        << __x.mean() << __space
01408        << __x.sigma();
01409       if (__x._M_saved_available)
01410     __os << __space << __x._M_saved;
01411 
01412       __os.flags(__flags);
01413       __os.fill(__fill);
01414       __os.precision(__precision);
01415       return __os;
01416     }
01417 
01418   template<typename _RealType, typename _CharT, typename _Traits>
01419     std::basic_istream<_CharT, _Traits>&
01420     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01421            normal_distribution<_RealType>& __x)
01422     {
01423       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01424       typedef typename __istream_type::ios_base    __ios_base;
01425 
01426       const typename __ios_base::fmtflags __flags = __is.flags();
01427       __is.flags(__ios_base::dec | __ios_base::skipws);
01428 
01429       __is >> __x._M_saved_available >> __x._M_mean
01430        >> __x._M_sigma;
01431       if (__x._M_saved_available)
01432     __is >> __x._M_saved;
01433 
01434       __is.flags(__flags);
01435       return __is;
01436     }
01437 
01438 
01439   template<typename _RealType>
01440     void
01441     gamma_distribution<_RealType>::
01442     _M_initialize()
01443     {
01444       if (_M_alpha >= 1)
01445     _M_l_d = std::sqrt(2 * _M_alpha - 1);
01446       else
01447     _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
01448           * (1 - _M_alpha));
01449     }
01450 
01451   /**
01452    * Cheng's rejection algorithm GB for alpha >= 1 and a modification
01453    * of Vaduva's rejection from Weibull algorithm due to Devroye for
01454    * alpha < 1.
01455    *
01456    * References:
01457    * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
01458    * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
01459    *
01460    * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
01461    * and Composition Procedures." Math. Operationsforschung and Statistik,
01462    * Series in Statistics, 8, 545-576, 1977.
01463    *
01464    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
01465    * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
01466    */
01467   template<typename _RealType>
01468     template<class _UniformRandomNumberGenerator>
01469       typename gamma_distribution<_RealType>::result_type
01470       gamma_distribution<_RealType>::
01471       operator()(_UniformRandomNumberGenerator& __urng)
01472       {
01473     result_type __x;
01474 
01475     bool __reject;
01476     if (_M_alpha >= 1)
01477       {
01478         // alpha - log(4)
01479         const result_type __b = _M_alpha
01480           - result_type(1.3862943611198906188344642429163531L);
01481         const result_type __c = _M_alpha + _M_l_d;
01482         const result_type __1l = 1 / _M_l_d;
01483 
01484         // 1 + log(9 / 2)
01485         const result_type __k = 2.5040773967762740733732583523868748L;
01486 
01487         do
01488           {
01489         const result_type __u = __urng();
01490         const result_type __v = __urng();
01491 
01492         const result_type __y = __1l * std::log(__v / (1 - __v));
01493         __x = _M_alpha * std::exp(__y);
01494 
01495         const result_type __z = __u * __v * __v;
01496         const result_type __r = __b + __c * __y - __x;
01497 
01498         __reject = __r < result_type(4.5) * __z - __k;
01499         if (__reject)
01500           __reject = __r < std::log(__z);
01501           }
01502         while (__reject);
01503       }
01504     else
01505       {
01506         const result_type __c = 1 / _M_alpha;
01507 
01508         do
01509           {
01510         const result_type __z = -std::log(__urng());
01511         const result_type __e = -std::log(__urng());
01512 
01513         __x = std::pow(__z, __c);
01514 
01515         __reject = __z + __e < _M_l_d + __x;
01516           }
01517         while (__reject);
01518       }
01519 
01520     return __x;
01521       }
01522 
01523   template<typename _RealType, typename _CharT, typename _Traits>
01524     std::basic_ostream<_CharT, _Traits>&
01525     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01526            const gamma_distribution<_RealType>& __x)
01527     {
01528       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01529       typedef typename __ostream_type::ios_base    __ios_base;
01530 
01531       const typename __ios_base::fmtflags __flags = __os.flags();
01532       const _CharT __fill = __os.fill();
01533       const std::streamsize __precision = __os.precision();
01534       __os.flags(__ios_base::scientific | __ios_base::left);
01535       __os.fill(__os.widen(' '));
01536       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
01537 
01538       __os << __x.alpha();
01539 
01540       __os.flags(__flags);
01541       __os.fill(__fill);
01542       __os.precision(__precision);
01543       return __os;
01544     }
01545 
01546 _GLIBCXX_END_NAMESPACE
01547 }

Generated on Thu Nov 1 13:12:21 2007 for libstdc++ by  doxygen 1.5.1