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

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