l4/pkg/libstdc++-v3/contrib/libstdc++-v3-4.3.3/include/tr1_impl/random.tcc

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