# 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 }