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