libstdc++-v3/include/tr1/random.tcc

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