contrib/gcc-4.4/libstdc++-v3/include/tr1_impl/random.tcc

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