contrib/libstdc++/include/tr1/random.tcc

   1 // random number generation (out of line) -*- C++ -*-
   2
   3 // Copyright (C) 2006, 2007 Free Software Foundation, Inc.
   4 //
   5 // This file is part of the GNU ISO C++ Library.  This library is free
   6 // software; you can redistribute it and/or modify it under the
   7 // terms of the GNU General Public License as published by the
   8 // Free Software Foundation; either version 2, or (at your option)
   9 // any later version.
  10
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 // GNU General Public License for more details.
  15
  16 // You should have received a copy of the GNU General Public License along
  17 // with this library; see the file COPYING.  If not, write to the Free
  18 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
  19 // USA.
  20
  21 // As a special exception, you may use this file as part of a free software
  22 // library without restriction.  Specifically, if other files instantiate
  23 // templates or use macros or inline functions from this file, or you compile
  24 // this file and link it with other files to produce an executable, this
  25 // file does not by itself cause the resulting executable to be covered by
  26 // the GNU General Public License.  This exception does not however
  27 // invalidate any other reasons why the executable file might be covered by
  28 // the GNU General Public License.
  29
  30 /** @file tr1/random.tcc
  31  *  This is a TR1 C++ Library header.
  32  */
  33
  34 namespace std
  35 {
  36 _GLIBCXX_BEGIN_NAMESPACE(tr1)
  37
  38   /*
  39    * (Further) implementation-space details.
  40    */
  41   namespace __detail
  42   {
  43     // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
  44     // integer overflow.
  45     //
  46     // Because a and c are compile-time integral constants the compiler kindly
  47     // elides any unreachable paths.
  48     //
  49     // Preconditions:  a > 0, m > 0.
  50     //
  51     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
  52       struct _Mod
  53       {
  54         static _Tp
  55         __calc(_Tp __x)
  56         {
  57           if (__a == 1)
  58             __x %= __m;
  59           else
  60             {
  61               static const _Tp __q = __m / __a;
  62               static const _Tp __r = __m % __a;
  63
  64               _Tp __t1 = __a * (__x % __q);
  65               _Tp __t2 = __r * (__x / __q);
  66               if (__t1 >= __t2)
  67                 __x = __t1 - __t2;
  68               else
  69                 __x = __m - __t2 + __t1;
  70             }
  71
  72           if (__c != 0)
  73             {
  74               const _Tp __d = __m - __x;
  75               if (__d > __c)
  76                 __x += __c;
  77               else
  78                 __x = __c - __d;
  79             }
  80           return __x;
  81         }
  82       };
  83
  84     // Special case for m == 0 -- use unsigned integer overflow as modulo
  85     // operator.
  86     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
  87       struct _Mod<_Tp, __a, __c, __m, true>
  88       {
  89         static _Tp
  90         __calc(_Tp __x)
  91         { return __a * __x + __c; }
  92       };
  93   } // namespace __detail
  94
  95   /**
  96    * Seeds the LCR with integral value @p __x0, adjusted so that the
  97    * ring identity is never a member of the convergence set.
  98    */
  99   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 100     void
 101     linear_congruential<_UIntType, __a, __c, __m>::
 102     seed(unsigned long __x0)
 103     {
 104       if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
 105           && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 106         _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
 107       else
 108         _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
 109     }
 110
 111   /**
 112    * Seeds the LCR engine with a value generated by @p __g.
 113    */
 114   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 115     template<class _Gen>
 116       void
 117       linear_congruential<_UIntType, __a, __c, __m>::
 118       seed(_Gen& __g, false_type)
 119       {
 120         _UIntType __x0 = __g();
 121         if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
 122             && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 123           _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
 124         else
 125           _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
 126       }
 127
 128   /**
 129    * Gets the next generated value in sequence.
 130    */
 131   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 132     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 133     linear_congruential<_UIntType, __a, __c, __m>::
 134     operator()()
 135     {
 136       _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
 137       return _M_x;
 138     }
 139
 140   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 141            typename _CharT, typename _Traits>
 142     std::basic_ostream<_CharT, _Traits>&
 143     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 144                const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
 145     {
 146       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 147       typedef typename __ostream_type::ios_base    __ios_base;
 148
 149       const typename __ios_base::fmtflags __flags = __os.flags();
 150       const _CharT __fill = __os.fill();
 151       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 152       __os.fill(__os.widen(' '));
 153
 154       __os << __lcr._M_x;
 155
 156       __os.flags(__flags);
 157       __os.fill(__fill);
 158       return __os;
 159     }
 160
 161   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 162            typename _CharT, typename _Traits>
 163     std::basic_istream<_CharT, _Traits>&
 164     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 165                linear_congruential<_UIntType, __a, __c, __m>& __lcr)
 166     {
 167       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 168       typedef typename __istream_type::ios_base    __ios_base;
 169
 170       const typename __ios_base::fmtflags __flags = __is.flags();
 171       __is.flags(__ios_base::dec);
 172
 173       __is >> __lcr._M_x;
 174
 175       __is.flags(__flags);
 176       return __is;
 177     }
 178
 179
 180   template<class _UIntType, int __w, int __n, int __m, int __r,
 181            _UIntType __a, int __u, int __s,
 182            _UIntType __b, int __t, _UIntType __c, int __l>
 183     void
 184     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 185                      __b, __t, __c, __l>::
 186     seed(unsigned long __value)
 187     {
 188       _M_x[0] = __detail::__mod<_UIntType, 1, 0,
 189         __detail::_Shift<_UIntType, __w>::__value>(__value);
 190
 191       for (int __i = 1; __i < state_size; ++__i)
 192         {
 193           _UIntType __x = _M_x[__i - 1];
 194           __x ^= __x >> (__w - 2);
 195           __x *= 1812433253ul;
 196           __x += __i;
 197           _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
 198             __detail::_Shift<_UIntType, __w>::__value>(__x);
 199         }
 200       _M_p = state_size;
 201     }
 202
 203   template<class _UIntType, int __w, int __n, int __m, int __r,
 204            _UIntType __a, int __u, int __s,
 205            _UIntType __b, int __t, _UIntType __c, int __l>
 206     template<class _Gen>
 207       void
 208       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 209                        __b, __t, __c, __l>::
 210       seed(_Gen& __gen, false_type)
 211       {
 212         for (int __i = 0; __i < state_size; ++__i)
 213           _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
 214             __detail::_Shift<_UIntType, __w>::__value>(__gen());
 215         _M_p = state_size;
 216       }
 217
 218   template<class _UIntType, int __w, int __n, int __m, int __r,
 219            _UIntType __a, int __u, int __s,
 220            _UIntType __b, int __t, _UIntType __c, int __l>
 221     typename
 222     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 223                      __b, __t, __c, __l>::result_type
 224     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 225                      __b, __t, __c, __l>::
 226     operator()()
 227     {
 228       // Reload the vector - cost is O(n) amortized over n calls.
 229       if (_M_p >= state_size)
 230         {
 231           const _UIntType __upper_mask = (~_UIntType()) << __r;
 232           const _UIntType __lower_mask = ~__upper_mask;
 233
 234           for (int __k = 0; __k < (__n - __m); ++__k)
 235             {
 236               _UIntType __y = ((_M_x[__k] & __upper_mask)
 237                                | (_M_x[__k + 1] & __lower_mask));
 238               _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
 239                            ^ ((__y & 0x01) ? __a : 0));
 240             }
 241
 242           for (int __k = (__n - __m); __k < (__n - 1); ++__k)
 243             {
 244               _UIntType __y = ((_M_x[__k] & __upper_mask)
 245                                | (_M_x[__k + 1] & __lower_mask));
 246               _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
 247                            ^ ((__y & 0x01) ? __a : 0));
 248             }
 249
 250           _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
 251                            | (_M_x[0] & __lower_mask));
 252           _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
 253                            ^ ((__y & 0x01) ? __a : 0));
 254           _M_p = 0;
 255         }
 256
 257       // Calculate o(x(i)).
 258       result_type __z = _M_x[_M_p++];
 259       __z ^= (__z >> __u);
 260       __z ^= (__z << __s) & __b;
 261       __z ^= (__z << __t) & __c;
 262       __z ^= (__z >> __l);
 263
 264       return __z;
 265     }
 266
 267   template<class _UIntType, int __w, int __n, int __m, int __r,
 268            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
 269            _UIntType __c, int __l,
 270            typename _CharT, typename _Traits>
 271     std::basic_ostream<_CharT, _Traits>&
 272     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 273                const mersenne_twister<_UIntType, __w, __n, __m,
 274                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
 275     {
 276       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 277       typedef typename __ostream_type::ios_base    __ios_base;
 278
 279       const typename __ios_base::fmtflags __flags = __os.flags();
 280       const _CharT __fill = __os.fill();
 281       const _CharT __space = __os.widen(' ');
 282       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 283       __os.fill(__space);
 284
 285       for (int __i = 0; __i < __n - 1; ++__i)
 286         __os << __x._M_x[__i] << __space;
 287       __os << __x._M_x[__n - 1];
 288
 289       __os.flags(__flags);
 290       __os.fill(__fill);
 291       return __os;
 292     }
 293
 294   template<class _UIntType, int __w, int __n, int __m, int __r,
 295            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
 296            _UIntType __c, int __l,
 297            typename _CharT, typename _Traits>
 298     std::basic_istream<_CharT, _Traits>&
 299     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 300                mersenne_twister<_UIntType, __w, __n, __m,
 301                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
 302     {
 303       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 304       typedef typename __istream_type::ios_base    __ios_base;
 305
 306       const typename __ios_base::fmtflags __flags = __is.flags();
 307       __is.flags(__ios_base::dec | __ios_base::skipws);
 308
 309       for (int __i = 0; __i < __n; ++__i)
 310         __is >> __x._M_x[__i];
 311
 312       __is.flags(__flags);
 313       return __is;
 314     }
 315
 316
 317   template<typename _IntType, _IntType __m, int __s, int __r>
 318     void
 319     subtract_with_carry<_IntType, __m, __s, __r>::
 320     seed(unsigned long __value)
 321     {
 322       if (__value == 0)
 323         __value = 19780503;
 324
 325       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
 326         __lcg(__value);
 327
 328       for (int __i = 0; __i < long_lag; ++__i)
 329         _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());
 330
 331       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 332       _M_p = 0;
 333     }
 334
 335   template<typename _IntType, _IntType __m, int __s, int __r>
 336     template<class _Gen>
 337       void
 338       subtract_with_carry<_IntType, __m, __s, __r>::
 339       seed(_Gen& __gen, false_type)
 340       {
 341         const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;
 342
 343         for (int __i = 0; __i < long_lag; ++__i)
 344           {
 345             _UIntType __tmp = 0;
 346             _UIntType __factor = 1;
 347             for (int __j = 0; __j < __n; ++__j)
 348               {
 349                 __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
 350                          (__gen()) * __factor;
 351                 __factor *= __detail::_Shift<_UIntType, 32>::__value;
 352               }
 353             _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
 354           }
 355         _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 356         _M_p = 0;
 357       }
 358
 359   template<typename _IntType, _IntType __m, int __s, int __r>
 360     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
 361     subtract_with_carry<_IntType, __m, __s, __r>::
 362     operator()()
 363     {
 364       // Derive short lag index from current index.
 365       int __ps = _M_p - short_lag;
 366       if (__ps < 0)
 367         __ps += long_lag;
 368
 369       // Calculate new x(i) without overflow or division.
 370       // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
 371       // cannot overflow.
 372       _UIntType __xi;
 373       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
 374         {
 375           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
 376           _M_carry = 0;
 377         }
 378       else
 379         {
 380           __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
 381           _M_carry = 1;
 382         }
 383       _M_x[_M_p] = __xi;
 384
 385       // Adjust current index to loop around in ring buffer.
 386       if (++_M_p >= long_lag)
 387         _M_p = 0;
 388
 389       return __xi;
 390     }
 391
 392   template<typename _IntType, _IntType __m, int __s, int __r,
 393            typename _CharT, typename _Traits>
 394     std::basic_ostream<_CharT, _Traits>&
 395     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 396                const subtract_with_carry<_IntType, __m, __s, __r>& __x)
 397     {
 398       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 399       typedef typename __ostream_type::ios_base    __ios_base;
 400
 401       const typename __ios_base::fmtflags __flags = __os.flags();
 402       const _CharT __fill = __os.fill();
 403       const _CharT __space = __os.widen(' ');
 404       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 405       __os.fill(__space);
 406
 407       for (int __i = 0; __i < __r; ++__i)
 408         __os << __x._M_x[__i] << __space;
 409       __os << __x._M_carry;
 410
 411       __os.flags(__flags);
 412       __os.fill(__fill);
 413       return __os;
 414     }
 415
 416   template<typename _IntType, _IntType __m, int __s, int __r,
 417            typename _CharT, typename _Traits>
 418     std::basic_istream<_CharT, _Traits>&
 419     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 420                subtract_with_carry<_IntType, __m, __s, __r>& __x)
 421     {
 422       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
 423       typedef typename __istream_type::ios_base    __ios_base;
 424
 425       const typename __ios_base::fmtflags __flags = __is.flags();
 426       __is.flags(__ios_base::dec | __ios_base::skipws);
 427
 428       for (int __i = 0; __i < __r; ++__i)
 429         __is >> __x._M_x[__i];
 430       __is >> __x._M_carry;
 431
 432       __is.flags(__flags);
 433       return __is;
 434     }
 435
 436
 437   template<typename _RealType, int __w, int __s, int __r>
 438     void
 439     subtract_with_carry_01<_RealType, __w, __s, __r>::
 440     _M_initialize_npows()
 441     {
 442       for (int __j = 0; __j < __n; ++__j)
 443 #if _GLIBCXX_USE_C99_MATH_TR1
 444         _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32);
 445 #else
 446         _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
 447 #endif
 448     }
 449
 450   template<typename _RealType, int __w, int __s, int __r>
 451     void
 452     subtract_with_carry_01<_RealType, __w, __s, __r>::
 453     seed(unsigned long __value)
 454     {
 455       if (__value == 0)
 456         __value = 19780503;
 457
 458       // _GLIBCXX_RESOLVE_LIB_DEFECTS
 459       // 512. Seeding subtract_with_carry_01 from a single unsigned long.
 460       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
 461         __lcg(__value);
 462
 463       this->seed(__lcg);
 464     }
 465
 466   template<typename _RealType, int __w, int __s, int __r>
 467     template<class _Gen>
 468       void
 469       subtract_with_carry_01<_RealType, __w, __s, __r>::
 470       seed(_Gen& __gen, false_type)
 471       {
 472         for (int __i = 0; __i < long_lag; ++__i)
 473           {
 474             for (int __j = 0; __j < __n - 1; ++__j)
 475               _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
 476             _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
 477               __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
 478           }
 479
 480         _M_carry = 1;
 481         for (int __j = 0; __j < __n; ++__j)
 482           if (_M_x[long_lag - 1][__j] != 0)
 483             {
 484               _M_carry = 0;
 485               break;
 486             }
 487
 488         _M_p = 0;
 489       }
 490
 491   template<typename _RealType, int __w, int __s, int __r>
 492     typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
 493     subtract_with_carry_01<_RealType, __w, __s, __r>::
 494     operator()()
 495     {
 496       // Derive short lag index from current index.
 497       int __ps = _M_p - short_lag;
 498       if (__ps < 0)
 499         __ps += long_lag;
 500
 501       _UInt32Type __new_carry;
 502       for (int __j = 0; __j < __n - 1; ++__j)
 503         {
 504           if (_M_x[__ps][__j] > _M_x[_M_p][__j]
 505               || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
 506             __new_carry = 0;
 507           else
 508             __new_carry = 1;
 509
 510           _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
 511           _M_carry = __new_carry;
 512         }
 513
 514       if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
 515           || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
 516         __new_carry = 0;
 517       else
 518         __new_carry = 1;
 519
 520       _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
 521         __detail::_Shift<_UInt32Type, __w % 32>::__value>
 522         (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
 523       _M_carry = __new_carry;
 524
 525       result_type __ret = 0.0;
 526       for (int __j = 0; __j < __n; ++__j)
 527         __ret += _M_x[_M_p][__j] * _M_npows[__j];
 528
 529       // Adjust current index to loop around in ring buffer.
 530       if (++_M_p >= long_lag)
 531         _M_p = 0;
 532
 533       return __ret;
 534     }
 535
 536   template<typename _RealType, int __w, int __s, int __r,
 537            typename _CharT, typename _Traits>
 538     std::basic_ostream<_CharT, _Traits>&
 539     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 540                const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
 541     {
 542       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 543       typedef typename __ostream_type::ios_base    __ios_base;
 544
 545       const typename __ios_base::fmtflags __flags = __os.flags();
 546       const _CharT __fill = __os.fill();
 547       const _CharT __space = __os.widen(' ');
 548       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 549       __os.fill(__space);
 550
 551       for (int __i = 0; __i < __r; ++__i)
 552         for (int __j = 0; __j < __x.__n; ++__j)
 553           __os << __x._M_x[__i][__j] << __space;
 554       __os << __x._M_carry;
 555
 556       __os.flags(__flags);
 557       __os.fill(__fill);
 558       return __os;
 559     }
 560
 561   template<typename _RealType, int __w, int __s, int __r,
 562            typename _CharT, typename _Traits>
 563     std::basic_istream<_CharT, _Traits>&
 564     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 565                subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
 566     {
 567       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 568       typedef typename __istream_type::ios_base    __ios_base;
 569
 570       const typename __ios_base::fmtflags __flags = __is.flags();
 571       __is.flags(__ios_base::dec | __ios_base::skipws);
 572
 573       for (int __i = 0; __i < __r; ++__i)
 574         for (int __j = 0; __j < __x.__n; ++__j)
 575           __is >> __x._M_x[__i][__j];
 576       __is >> __x._M_carry;
 577
 578       __is.flags(__flags);
 579       return __is;
 580     }
 581
 582
 583   template<class _UniformRandomNumberGenerator, int __p, int __r>
 584     typename discard_block<_UniformRandomNumberGenerator,
 585                            __p, __r>::result_type
 586     discard_block<_UniformRandomNumberGenerator, __p, __r>::
 587     operator()()
 588     {
 589       if (_M_n >= used_block)
 590         {
 591           while (_M_n < block_size)
 592             {
 593               _M_b();
 594               ++_M_n;
 595             }
 596           _M_n = 0;
 597         }
 598       ++_M_n;
 599       return _M_b();
 600     }
 601
 602   template<class _UniformRandomNumberGenerator, int __p, int __r,
 603            typename _CharT, typename _Traits>
 604     std::basic_ostream<_CharT, _Traits>&
 605     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 606                const discard_block<_UniformRandomNumberGenerator,
 607                __p, __r>& __x)
 608     {
 609       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 610       typedef typename __ostream_type::ios_base    __ios_base;
 611
 612       const typename __ios_base::fmtflags __flags = __os.flags();
 613       const _CharT __fill = __os.fill();
 614       const _CharT __space = __os.widen(' ');
 615       __os.flags(__ios_base::dec | __ios_base::fixed
 616                  | __ios_base::left);
 617       __os.fill(__space);
 618
 619       __os << __x._M_b << __space << __x._M_n;
 620
 621       __os.flags(__flags);
 622       __os.fill(__fill);
 623       return __os;
 624     }
 625
 626   template<class _UniformRandomNumberGenerator, int __p, int __r,
 627            typename _CharT, typename _Traits>
 628     std::basic_istream<_CharT, _Traits>&
 629     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 630                discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
 631     {
 632       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 633       typedef typename __istream_type::ios_base    __ios_base;
 634
 635       const typename __ios_base::fmtflags __flags = __is.flags();
 636       __is.flags(__ios_base::dec | __ios_base::skipws);
 637
 638       __is >> __x._M_b >> __x._M_n;
 639
 640       __is.flags(__flags);
 641       return __is;
 642     }
 643
 644
 645   template<class _UniformRandomNumberGenerator1, int __s1,
 646            class _UniformRandomNumberGenerator2, int __s2>
 647     void
 648     xor_combine<_UniformRandomNumberGenerator1, __s1,
 649                 _UniformRandomNumberGenerator2, __s2>::
 650     _M_initialize_max()
 651     {
 652       const int __w = std::numeric_limits<result_type>::digits;
 653
 654       const result_type __m1 =
 655         std::min(result_type(_M_b1.max() - _M_b1.min()),
 656                  __detail::_Shift<result_type, __w - __s1>::__value - 1);
 657
 658       const result_type __m2 =
 659         std::min(result_type(_M_b2.max() - _M_b2.min()),
 660                  __detail::_Shift<result_type, __w - __s2>::__value - 1);
 661
 662       // NB: In TR1 s1 is not required to be >= s2.
 663       if (__s1 < __s2)
 664         _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
 665       else
 666         _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
 667     }
 668
 669   template<class _UniformRandomNumberGenerator1, int __s1,
 670            class _UniformRandomNumberGenerator2, int __s2>
 671     typename xor_combine<_UniformRandomNumberGenerator1, __s1,
 672                          _UniformRandomNumberGenerator2, __s2>::result_type
 673     xor_combine<_UniformRandomNumberGenerator1, __s1,
 674                 _UniformRandomNumberGenerator2, __s2>::
 675     _M_initialize_max_aux(result_type __a, result_type __b, int __d)
 676     {
 677       const result_type __two2d = result_type(1) << __d;
 678       const result_type __c = __a * __two2d;
 679
 680       if (__a == 0 || __b < __two2d)
 681         return __c + __b;
 682
 683       const result_type __t = std::max(__c, __b);
 684       const result_type __u = std::min(__c, __b);
 685
 686       result_type __ub = __u;
 687       result_type __p;
 688       for (__p = 0; __ub != 1; __ub >>= 1)
 689         ++__p;
 690
 691       const result_type __two2p = result_type(1) << __p;
 692       const result_type __k = __t / __two2p;
 693
 694       if (__k & 1)
 695         return (__k + 1) * __two2p - 1;
 696
 697       if (__c >= __b)
 698         return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
 699                                                            / __two2d,
 700                                                            __u % __two2p, __d);
 701       else
 702         return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
 703                                                            / __two2d,
 704                                                            __t % __two2p, __d);
 705     }
 706
 707   template<class _UniformRandomNumberGenerator1, int __s1,
 708            class _UniformRandomNumberGenerator2, int __s2,
 709            typename _CharT, typename _Traits>
 710     std::basic_ostream<_CharT, _Traits>&
 711     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 712                const xor_combine<_UniformRandomNumberGenerator1, __s1,
 713                _UniformRandomNumberGenerator2, __s2>& __x)
 714     {
 715       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 716       typedef typename __ostream_type::ios_base    __ios_base;
 717
 718       const typename __ios_base::fmtflags __flags = __os.flags();
 719       const _CharT __fill = __os.fill();
 720       const _CharT __space = __os.widen(' ');
 721       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 722       __os.fill(__space);
 723
 724       __os << __x.base1() << __space << __x.base2();
 725
 726       __os.flags(__flags);
 727       __os.fill(__fill);
 728       return __os;
 729     }
 730
 731   template<class _UniformRandomNumberGenerator1, int __s1,
 732            class _UniformRandomNumberGenerator2, int __s2,
 733            typename _CharT, typename _Traits>
 734     std::basic_istream<_CharT, _Traits>&
 735     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 736                xor_combine<_UniformRandomNumberGenerator1, __s1,
 737                _UniformRandomNumberGenerator2, __s2>& __x)
 738     {
 739       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 740       typedef typename __istream_type::ios_base    __ios_base;
 741
 742       const typename __ios_base::fmtflags __flags = __is.flags();
 743       __is.flags(__ios_base::skipws);
 744
 745       __is >> __x._M_b1 >> __x._M_b2;
 746
 747       __is.flags(__flags);
 748       return __is;
 749     }
 750
 751
 752   template<typename _IntType, typename _CharT, typename _Traits>
 753     std::basic_ostream<_CharT, _Traits>&
 754     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 755                const uniform_int<_IntType>& __x)
 756     {
 757       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 758       typedef typename __ostream_type::ios_base    __ios_base;
 759
 760       const typename __ios_base::fmtflags __flags = __os.flags();
 761       const _CharT __fill = __os.fill();
 762       const _CharT __space = __os.widen(' ');
 763       __os.flags(__ios_base::scientific | __ios_base::left);
 764       __os.fill(__space);
 765
 766       __os << __x.min() << __space << __x.max();
 767
 768       __os.flags(__flags);
 769       __os.fill(__fill);
 770       return __os;
 771     }
 772
 773   template<typename _IntType, typename _CharT, typename _Traits>
 774     std::basic_istream<_CharT, _Traits>&
 775     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 776                uniform_int<_IntType>& __x)
 777     {
 778       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 779       typedef typename __istream_type::ios_base    __ios_base;
 780
 781       const typename __ios_base::fmtflags __flags = __is.flags();
 782       __is.flags(__ios_base::dec | __ios_base::skipws);
 783
 784       __is >> __x._M_min >> __x._M_max;
 785
 786       __is.flags(__flags);
 787       return __is;
 788     }
 789
 790
 791   template<typename _CharT, typename _Traits>
 792     std::basic_ostream<_CharT, _Traits>&
 793     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 794                const bernoulli_distribution& __x)
 795     {
 796       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 797       typedef typename __ostream_type::ios_base    __ios_base;
 798
 799       const typename __ios_base::fmtflags __flags = __os.flags();
 800       const _CharT __fill = __os.fill();
 801       const std::streamsize __precision = __os.precision();
 802       __os.flags(__ios_base::scientific | __ios_base::left);
 803       __os.fill(__os.widen(' '));
 804       __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);
 805
 806       __os << __x.p();
 807
 808       __os.flags(__flags);
 809       __os.fill(__fill);
 810       __os.precision(__precision);
 811       return __os;
 812     }
 813
 814
 815   template<typename _IntType, typename _RealType>
 816     template<class _UniformRandomNumberGenerator>
 817       typename geometric_distribution<_IntType, _RealType>::result_type
 818       geometric_distribution<_IntType, _RealType>::
 819       operator()(_UniformRandomNumberGenerator& __urng)
 820       {
 821         // About the epsilon thing see this thread:
 822         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
 823         const _RealType __naf =
 824           (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
 825         // The largest _RealType convertible to _IntType.
 826         const _RealType __thr =
 827           std::numeric_limits<_IntType>::max() + __naf;
 828
 829         _RealType __cand;
 830         do
 831           __cand = std::ceil(std::log(__urng()) / _M_log_p);
 832         while (__cand >= __thr);
 833
 834         return result_type(__cand + __naf);
 835       }
 836
 837   template<typename _IntType, typename _RealType,
 838            typename _CharT, typename _Traits>
 839     std::basic_ostream<_CharT, _Traits>&
 840     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 841                const geometric_distribution<_IntType, _RealType>& __x)
 842     {
 843       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 844       typedef typename __ostream_type::ios_base    __ios_base;
 845
 846       const typename __ios_base::fmtflags __flags = __os.flags();
 847       const _CharT __fill = __os.fill();
 848       const std::streamsize __precision = __os.precision();
 849       __os.flags(__ios_base::scientific | __ios_base::left);
 850       __os.fill(__os.widen(' '));
 851       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
 852
 853       __os << __x.p();
 854
 855       __os.flags(__flags);
 856       __os.fill(__fill);
 857       __os.precision(__precision);
 858       return __os;
 859     }
 860
 861
 862   template<typename _IntType, typename _RealType>
 863     void
 864     poisson_distribution<_IntType, _RealType>::
 865     _M_initialize()
 866     {
 867 #if _GLIBCXX_USE_C99_MATH_TR1
 868       if (_M_mean >= 12)
 869         {
 870           const _RealType __m = std::floor(_M_mean);
 871           _M_lm_thr = std::log(_M_mean);
 872           _M_lfm = std::tr1::lgamma(__m + 1);
 873           _M_sm = std::sqrt(__m);
 874
 875           const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
 876           const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
 877                                                               / __pi_4));
 878           _M_d = std::tr1::round(std::max(_RealType(6),
 879                                           std::min(__m, __dx)));
 880           const _RealType __cx = 2 * __m + _M_d;
 881           _M_scx = std::sqrt(__cx / 2);
 882           _M_1cx = 1 / __cx;
 883
 884           _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
 885           _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
 886         }
 887       else
 888 #endif
 889         _M_lm_thr = std::exp(-_M_mean);
 890       }
 891
 892   /**
 893    * A rejection algorithm when mean >= 12 and a simple method based
 894    * upon the multiplication of uniform random variates otherwise.
 895    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
 896    * is defined.
 897    *
 898    * Reference:
 899    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
 900    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
 901    */
 902   template<typename _IntType, typename _RealType>
 903     template<class _UniformRandomNumberGenerator>
 904       typename poisson_distribution<_IntType, _RealType>::result_type
 905       poisson_distribution<_IntType, _RealType>::
 906       operator()(_UniformRandomNumberGenerator& __urng)
 907       {
 908 #if _GLIBCXX_USE_C99_MATH_TR1
 909         if (_M_mean >= 12)
 910           {
 911             _RealType __x;
 912
 913             // See comments above...
 914             const _RealType __naf =
 915               (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
 916             const _RealType __thr =
 917               std::numeric_limits<_IntType>::max() + __naf;
 918
 919             const _RealType __m = std::floor(_M_mean);
 920             // sqrt(pi / 2)
 921             const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
 922             const _RealType __c1 = _M_sm * __spi_2;
 923             const _RealType __c2 = _M_c2b + __c1;
 924             const _RealType __c3 = __c2 + 1;
 925             const _RealType __c4 = __c3 + 1;
 926             // e^(1 / 78)
 927             const _RealType __e178 = 1.0129030479320018583185514777512983L;
 928             const _RealType __c5 = __c4 + __e178;
 929             const _RealType __c = _M_cb + __c5;
 930             const _RealType __2cx = 2 * (2 * __m + _M_d);
 931
 932             bool __reject = true;
 933             do
 934               {
 935                 const _RealType __u = __c * __urng();
 936                 const _RealType __e = -std::log(__urng());
 937
 938                 _RealType __w = 0.0;
 939
 940                 if (__u <= __c1)
 941                   {
 942                     const _RealType __n = _M_nd(__urng);
 943                     const _RealType __y = -std::abs(__n) * _M_sm - 1;
 944                     __x = std::floor(__y);
 945                     __w = -__n * __n / 2;
 946                     if (__x < -__m)
 947                       continue;
 948                   }
 949                 else if (__u <= __c2)
 950                   {
 951                     const _RealType __n = _M_nd(__urng);
 952                     const _RealType __y = 1 + std::abs(__n) * _M_scx;
 953                     __x = std::ceil(__y);
 954                     __w = __y * (2 - __y) * _M_1cx;
 955                     if (__x > _M_d)
 956                       continue;
 957                   }
 958                 else if (__u <= __c3)
 959                   // NB: This case not in the book, nor in the Errata,
 960                   // but should be ok...
 961                   __x = -1;
 962                 else if (__u <= __c4)
 963                   __x = 0;
 964                 else if (__u <= __c5)
 965                   __x = 1;
 966                 else
 967                   {
 968                     const _RealType __v = -std::log(__urng());
 969                     const _RealType __y = _M_d + __v * __2cx / _M_d;
 970                     __x = std::ceil(__y);
 971                     __w = -_M_d * _M_1cx * (1 + __y / 2);
 972                   }
 973
 974                 __reject = (__w - __e - __x * _M_lm_thr
 975                             > _M_lfm - std::tr1::lgamma(__x + __m + 1));
 976
 977                 __reject |= __x + __m >= __thr;
 978
 979               } while (__reject);
 980
 981             return result_type(__x + __m + __naf);
 982           }
 983         else
 984 #endif
 985           {
 986             _IntType     __x = 0;
 987             _RealType __prod = 1.0;
 988
 989             do
 990               {
 991                 __prod *= __urng();
 992                 __x += 1;
 993               }
 994             while (__prod > _M_lm_thr);
 995
 996             return __x - 1;
 997           }
 998       }
 999
1000   template<typename _IntType, typename _RealType,
1001            typename _CharT, typename _Traits>
1002     std::basic_ostream<_CharT, _Traits>&
1003     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1004                const poisson_distribution<_IntType, _RealType>& __x)
1005     {
1006       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1007       typedef typename __ostream_type::ios_base    __ios_base;
1008
1009       const typename __ios_base::fmtflags __flags = __os.flags();
1010       const _CharT __fill = __os.fill();
1011       const std::streamsize __precision = __os.precision();
1012       const _CharT __space = __os.widen(' ');
1013       __os.flags(__ios_base::scientific | __ios_base::left);
1014       __os.fill(__space);
1015       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1016
1017       __os << __x.mean() << __space << __x._M_nd;
1018
1019       __os.flags(__flags);
1020       __os.fill(__fill);
1021       __os.precision(__precision);
1022       return __os;
1023     }
1024
1025   template<typename _IntType, typename _RealType,
1026            typename _CharT, typename _Traits>
1027     std::basic_istream<_CharT, _Traits>&
1028     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1029                poisson_distribution<_IntType, _RealType>& __x)
1030     {
1031       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1032       typedef typename __istream_type::ios_base    __ios_base;
1033
1034       const typename __ios_base::fmtflags __flags = __is.flags();
1035       __is.flags(__ios_base::skipws);
1036
1037       __is >> __x._M_mean >> __x._M_nd;
1038       __x._M_initialize();
1039
1040       __is.flags(__flags);
1041       return __is;
1042     }
1043
1044
1045   template<typename _IntType, typename _RealType>
1046     void
1047     binomial_distribution<_IntType, _RealType>::
1048     _M_initialize()
1049     {
1050       const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1051
1052       _M_easy = true;
1053
1054 #if _GLIBCXX_USE_C99_MATH_TR1
1055       if (_M_t * __p12 >= 8)
1056         {
1057           _M_easy = false;
1058           const _RealType __np = std::floor(_M_t * __p12);
1059           const _RealType __pa = __np / _M_t;
1060           const _RealType __1p = 1 - __pa;
1061
1062           const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
1063           const _RealType __d1x =
1064             std::sqrt(__np * __1p * std::log(32 * __np
1065                                              / (81 * __pi_4 * __1p)));
1066           _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));
1067           const _RealType __d2x =
1068             std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1069                                              / (__pi_4 * __pa)));
1070           _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));
1071
1072           // sqrt(pi / 2)
1073           const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1074           _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1075           _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1076           _M_c = 2 * _M_d1 / __np;
1077           _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1078           const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
1079           const _RealType __s1s = _M_s1 * _M_s1;
1080           _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1081                              * 2 * __s1s / _M_d1
1082                              * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1083           const _RealType __s2s = _M_s2 * _M_s2;
1084           _M_s = (_M_a123 + 2 * __s2s / _M_d2
1085                   * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1086           _M_lf = (std::tr1::lgamma(__np + 1)
1087                    + std::tr1::lgamma(_M_t - __np + 1));
1088           _M_lp1p = std::log(__pa / __1p);
1089
1090           _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1091         }
1092       else
1093 #endif
1094         _M_q = -std::log(1 - __p12);
1095     }
1096
1097   template<typename _IntType, typename _RealType>
1098     template<class _UniformRandomNumberGenerator>
1099       typename binomial_distribution<_IntType, _RealType>::result_type
1100       binomial_distribution<_IntType, _RealType>::
1101       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1102       {
1103         _IntType    __x = 0;
1104         _RealType __sum = 0;
1105
1106         do
1107           {
1108             const _RealType __e = -std::log(__urng());
1109             __sum += __e / (__t - __x);
1110             __x += 1;
1111           }
1112         while (__sum <= _M_q);
1113
1114         return __x - 1;
1115       }
1116
1117   /**
1118    * A rejection algorithm when t * p >= 8 and a simple waiting time
1119    * method - the second in the referenced book - otherwise.
1120    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1121    * is defined.
1122    *
1123    * Reference:
1124    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1125    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1126    */
1127   template<typename _IntType, typename _RealType>
1128     template<class _UniformRandomNumberGenerator>
1129       typename binomial_distribution<_IntType, _RealType>::result_type
1130       binomial_distribution<_IntType, _RealType>::
1131       operator()(_UniformRandomNumberGenerator& __urng)
1132       {
1133         result_type __ret;
1134         const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1135
1136 #if _GLIBCXX_USE_C99_MATH_TR1
1137         if (!_M_easy)
1138           {
1139             _RealType __x;
1140
1141             // See comments above...
1142             const _RealType __naf =
1143               (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
1144             const _RealType __thr =
1145               std::numeric_limits<_IntType>::max() + __naf;
1146
1147             const _RealType __np = std::floor(_M_t * __p12);
1148             const _RealType __pa = __np / _M_t;
1149
1150             // sqrt(pi / 2)
1151             const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1152             const _RealType __a1 = _M_a1;
1153             const _RealType __a12 = __a1 + _M_s2 * __spi_2;
1154             const _RealType __a123 = _M_a123;
1155             const _RealType __s1s = _M_s1 * _M_s1;
1156             const _RealType __s2s = _M_s2 * _M_s2;
1157
1158             bool __reject;
1159             do
1160               {
1161                 const _RealType __u = _M_s * __urng();
1162
1163                 _RealType __v;
1164
1165                 if (__u <= __a1)
1166                   {
1167                     const _RealType __n = _M_nd(__urng);
1168                     const _RealType __y = _M_s1 * std::abs(__n);
1169                     __reject = __y >= _M_d1;
1170                     if (!__reject)
1171                       {
1172                         const _RealType __e = -std::log(__urng());
1173                         __x = std::floor(__y);
1174                         __v = -__e - __n * __n / 2 + _M_c;
1175                       }
1176                   }
1177                 else if (__u <= __a12)
1178                   {
1179                     const _RealType __n = _M_nd(__urng);
1180                     const _RealType __y = _M_s2 * std::abs(__n);
1181                     __reject = __y >= _M_d2;
1182                     if (!__reject)
1183                       {
1184                         const _RealType __e = -std::log(__urng());
1185                         __x = std::floor(-__y);
1186                         __v = -__e - __n * __n / 2;
1187                       }
1188                   }
1189                 else if (__u <= __a123)
1190                   {
1191                     const _RealType __e1 = -std::log(__urng());
1192                     const _RealType __e2 = -std::log(__urng());
1193
1194                     const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
1195                     __x = std::floor(__y);
1196                     __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
1197                                             -__y / (2 * __s1s)));
1198                     __reject = false;
1199                   }
1200                 else
1201                   {
1202                     const _RealType __e1 = -std::log(__urng());
1203                     const _RealType __e2 = -std::log(__urng());
1204
1205                     const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
1206                     __x = std::floor(-__y);
1207                     __v = -__e2 - _M_d2 * __y / (2 * __s2s);
1208                     __reject = false;
1209                   }
1210
1211                 __reject = __reject || __x < -__np || __x > _M_t - __np;
1212                 if (!__reject)
1213                   {
1214                     const _RealType __lfx =
1215                       std::tr1::lgamma(__np + __x + 1)
1216                       + std::tr1::lgamma(_M_t - (__np + __x) + 1);
1217                     __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
1218                   }
1219
1220                 __reject |= __x + __np >= __thr;
1221               }
1222             while (__reject);
1223
1224             __x += __np + __naf;
1225
1226             const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x));
1227             __ret = _IntType(__x) + __z;
1228           }
1229         else
1230 #endif
1231           __ret = _M_waiting(__urng, _M_t);
1232
1233         if (__p12 != _M_p)
1234           __ret = _M_t - __ret;
1235         return __ret;
1236       }
1237
1238   template<typename _IntType, typename _RealType,
1239            typename _CharT, typename _Traits>
1240     std::basic_ostream<_CharT, _Traits>&
1241     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1242                const binomial_distribution<_IntType, _RealType>& __x)
1243     {
1244       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1245       typedef typename __ostream_type::ios_base    __ios_base;
1246
1247       const typename __ios_base::fmtflags __flags = __os.flags();
1248       const _CharT __fill = __os.fill();
1249       const std::streamsize __precision = __os.precision();
1250       const _CharT __space = __os.widen(' ');
1251       __os.flags(__ios_base::scientific | __ios_base::left);
1252       __os.fill(__space);
1253       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1254
1255       __os << __x.t() << __space << __x.p()
1256            << __space << __x._M_nd;
1257
1258       __os.flags(__flags);
1259       __os.fill(__fill);
1260       __os.precision(__precision);
1261       return __os;
1262     }
1263
1264   template<typename _IntType, typename _RealType,
1265            typename _CharT, typename _Traits>
1266     std::basic_istream<_CharT, _Traits>&
1267     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1268                binomial_distribution<_IntType, _RealType>& __x)
1269     {
1270       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1271       typedef typename __istream_type::ios_base    __ios_base;
1272
1273       const typename __ios_base::fmtflags __flags = __is.flags();
1274       __is.flags(__ios_base::dec | __ios_base::skipws);
1275
1276       __is >> __x._M_t >> __x._M_p >> __x._M_nd;
1277       __x._M_initialize();
1278
1279       __is.flags(__flags);
1280       return __is;
1281     }
1282
1283
1284   template<typename _RealType, typename _CharT, typename _Traits>
1285     std::basic_ostream<_CharT, _Traits>&
1286     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1287                const uniform_real<_RealType>& __x)
1288     {
1289       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1290       typedef typename __ostream_type::ios_base    __ios_base;
1291
1292       const typename __ios_base::fmtflags __flags = __os.flags();
1293       const _CharT __fill = __os.fill();
1294       const std::streamsize __precision = __os.precision();
1295       const _CharT __space = __os.widen(' ');
1296       __os.flags(__ios_base::scientific | __ios_base::left);
1297       __os.fill(__space);
1298       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1299
1300       __os << __x.min() << __space << __x.max();
1301
1302       __os.flags(__flags);
1303       __os.fill(__fill);
1304       __os.precision(__precision);
1305       return __os;
1306     }
1307
1308   template<typename _RealType, typename _CharT, typename _Traits>
1309     std::basic_istream<_CharT, _Traits>&
1310     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1311                uniform_real<_RealType>& __x)
1312     {
1313       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1314       typedef typename __istream_type::ios_base    __ios_base;
1315
1316       const typename __ios_base::fmtflags __flags = __is.flags();
1317       __is.flags(__ios_base::skipws);
1318
1319       __is >> __x._M_min >> __x._M_max;
1320
1321       __is.flags(__flags);
1322       return __is;
1323     }
1324
1325
1326   template<typename _RealType, typename _CharT, typename _Traits>
1327     std::basic_ostream<_CharT, _Traits>&
1328     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1329                const exponential_distribution<_RealType>& __x)
1330     {
1331       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1332       typedef typename __ostream_type::ios_base    __ios_base;
1333
1334       const typename __ios_base::fmtflags __flags = __os.flags();
1335       const _CharT __fill = __os.fill();
1336       const std::streamsize __precision = __os.precision();
1337       __os.flags(__ios_base::scientific | __ios_base::left);
1338       __os.fill(__os.widen(' '));
1339       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1340
1341       __os << __x.lambda();
1342
1343       __os.flags(__flags);
1344       __os.fill(__fill);
1345       __os.precision(__precision);
1346       return __os;
1347     }
1348
1349
1350   /**
1351    * Polar method due to Marsaglia.
1352    *
1353    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1354    * New York, 1986, Ch. V, Sect. 4.4.
1355    */
1356   template<typename _RealType>
1357     template<class _UniformRandomNumberGenerator>
1358       typename normal_distribution<_RealType>::result_type
1359       normal_distribution<_RealType>::
1360       operator()(_UniformRandomNumberGenerator& __urng)
1361       {
1362         result_type __ret;
1363
1364         if (_M_saved_available)
1365           {
1366             _M_saved_available = false;
1367             __ret = _M_saved;
1368           }
1369         else
1370           {
1371             result_type __x, __y, __r2;
1372             do
1373               {
1374                 __x = result_type(2.0) * __urng() - 1.0;
1375                 __y = result_type(2.0) * __urng() - 1.0;
1376                 __r2 = __x * __x + __y * __y;
1377               }
1378             while (__r2 > 1.0 || __r2 == 0.0);
1379
1380             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1381             _M_saved = __x * __mult;
1382             _M_saved_available = true;
1383             __ret = __y * __mult;
1384           }
1385
1386         __ret = __ret * _M_sigma + _M_mean;
1387         return __ret;
1388       }
1389
1390   template<typename _RealType, typename _CharT, typename _Traits>
1391     std::basic_ostream<_CharT, _Traits>&
1392     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1393                const normal_distribution<_RealType>& __x)
1394     {
1395       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1396       typedef typename __ostream_type::ios_base    __ios_base;
1397
1398       const typename __ios_base::fmtflags __flags = __os.flags();
1399       const _CharT __fill = __os.fill();
1400       const std::streamsize __precision = __os.precision();
1401       const _CharT __space = __os.widen(' ');
1402       __os.flags(__ios_base::scientific | __ios_base::left);
1403       __os.fill(__space);
1404       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1405
1406       __os << __x._M_saved_available << __space
1407            << __x.mean() << __space
1408            << __x.sigma();
1409       if (__x._M_saved_available)
1410         __os << __space << __x._M_saved;
1411
1412       __os.flags(__flags);
1413       __os.fill(__fill);
1414       __os.precision(__precision);
1415       return __os;
1416     }
1417
1418   template<typename _RealType, typename _CharT, typename _Traits>
1419     std::basic_istream<_CharT, _Traits>&
1420     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1421                normal_distribution<_RealType>& __x)
1422     {
1423       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1424       typedef typename __istream_type::ios_base    __ios_base;
1425
1426       const typename __ios_base::fmtflags __flags = __is.flags();
1427       __is.flags(__ios_base::dec | __ios_base::skipws);
1428
1429       __is >> __x._M_saved_available >> __x._M_mean
1430            >> __x._M_sigma;
1431       if (__x._M_saved_available)
1432         __is >> __x._M_saved;
1433
1434       __is.flags(__flags);
1435       return __is;
1436     }
1437
1438
1439   template<typename _RealType>
1440     void
1441     gamma_distribution<_RealType>::
1442     _M_initialize()
1443     {
1444       if (_M_alpha >= 1)
1445         _M_l_d = std::sqrt(2 * _M_alpha - 1);
1446       else
1447         _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
1448                   * (1 - _M_alpha));
1449     }
1450
1451   /**
1452    * Cheng's rejection algorithm GB for alpha >= 1 and a modification
1453    * of Vaduva's rejection from Weibull algorithm due to Devroye for
1454    * alpha < 1.
1455    *
1456    * References:
1457    * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
1458    * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
1459    *
1460    * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
1461    * and Composition Procedures." Math. Operationsforschung and Statistik,
1462    * Series in Statistics, 8, 545-576, 1977.
1463    *
1464    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1465    * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
1466    */
1467   template<typename _RealType>
1468     template<class _UniformRandomNumberGenerator>
1469       typename gamma_distribution<_RealType>::result_type
1470       gamma_distribution<_RealType>::
1471       operator()(_UniformRandomNumberGenerator& __urng)
1472       {
1473         result_type __x;
1474
1475         bool __reject;
1476         if (_M_alpha >= 1)
1477           {
1478             // alpha - log(4)
1479             const result_type __b = _M_alpha
1480               - result_type(1.3862943611198906188344642429163531L);
1481             const result_type __c = _M_alpha + _M_l_d;
1482             const result_type __1l = 1 / _M_l_d;
1483
1484             // 1 + log(9 / 2)
1485             const result_type __k = 2.5040773967762740733732583523868748L;
1486
1487             do
1488               {
1489                 const result_type __u = __urng();
1490                 const result_type __v = __urng();
1491
1492                 const result_type __y = __1l * std::log(__v / (1 - __v));
1493                 __x = _M_alpha * std::exp(__y);
1494
1495                 const result_type __z = __u * __v * __v;
1496                 const result_type __r = __b + __c * __y - __x;
1497
1498                 __reject = __r < result_type(4.5) * __z - __k;
1499                 if (__reject)
1500                   __reject = __r < std::log(__z);
1501               }
1502             while (__reject);
1503           }
1504         else
1505           {
1506             const result_type __c = 1 / _M_alpha;
1507
1508             do
1509               {
1510                 const result_type __z = -std::log(__urng());
1511                 const result_type __e = -std::log(__urng());
1512
1513                 __x = std::pow(__z, __c);
1514
1515                 __reject = __z + __e < _M_l_d + __x;
1516               }
1517             while (__reject);
1518           }
1519
1520         return __x;
1521       }
1522
1523   template<typename _RealType, typename _CharT, typename _Traits>
1524     std::basic_ostream<_CharT, _Traits>&
1525     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1526                const gamma_distribution<_RealType>& __x)
1527     {
1528       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1529       typedef typename __ostream_type::ios_base    __ios_base;
1530
1531       const typename __ios_base::fmtflags __flags = __os.flags();
1532       const _CharT __fill = __os.fill();
1533       const std::streamsize __precision = __os.precision();
1534       __os.flags(__ios_base::scientific | __ios_base::left);
1535       __os.fill(__os.widen(' '));
1536       __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1537
1538       __os << __x.alpha();
1539
1540       __os.flags(__flags);
1541       __os.fill(__fill);
1542       __os.precision(__precision);
1543       return __os;
1544     }
1545
1546 _GLIBCXX_END_NAMESPACE
1547 }