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

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