random.tcc

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

Generated on Thu Jul 23 21:16:13 2009 for libstdc++ by  doxygen 1.5.8