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

   1 // random number generation (out of line) -*- C++ -*-
   2
   3 // Copyright (C) 2006 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 namespace std
  31 {
  32 _GLIBCXX_BEGIN_NAMESPACE(tr1)
  33
  34   /*
  35    * Implementation-space details.
  36    */
  37   namespace
  38   {
  39     // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
  40     // integer overflow.
  41     //
  42     // Because a and c are compile-time integral constants the compiler kindly
  43     // elides any unreachable paths.
  44     //
  45     // Preconditions:  a > 0, m > 0.
  46     //
  47     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
  48       struct _Mod
  49       {
  50         static _Tp
  51         __calc(_Tp __x)
  52         {
  53           if (__a == 1)
  54             __x %= __m;
  55           else
  56             {
  57               static const _Tp __q = __m / __a;
  58               static const _Tp __r = __m % __a;
  59
  60               _Tp __t1 = __a * (__x % __q);
  61               _Tp __t2 = __r * (__x / __q);
  62               if (__t1 >= __t2)
  63                 __x = __t1 - __t2;
  64               else
  65                 __x = __m - __t2 + __t1;
  66             }
  67
  68           if (__c != 0)
  69             {
  70               const _Tp __d = __m - __x;
  71               if (__d > __c)
  72                 __x += __c;
  73               else
  74                 __x = __c - __d;
  75             }
  76           return __x;
  77         }
  78       };
  79
  80     // Special case for m == 0 -- use unsigned integer overflow as modulo
  81     // operator.
  82     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
  83       struct _Mod<_Tp, __a, __c, __m, true>
  84       {
  85         static _Tp
  86         __calc(_Tp __x)
  87         { return __a * __x + __c; }
  88       };
  89
  90     // Dispatch based on modulus value to prevent divide-by-zero compile-time
  91     // errors when m == 0.
  92     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
  93       inline _Tp
  94       __mod(_Tp __x)
  95       { return _Mod<_Tp, __a, __c, __m, __m == 0>::__calc(__x); }
  96
  97     // See N1822.
  98     template<typename _RealType>
  99       struct _Max_digits10
 100       {
 101         static const std::streamsize __value =
 102           2 + std::numeric_limits<_RealType>::digits * 3010/10000;
 103       };
 104
 105     template<typename _ValueT>
 106       struct _To_Unsigned_Type
 107       { typedef _ValueT _Type; };
 108
 109     template<>
 110       struct _To_Unsigned_Type<short>
 111       { typedef unsigned short _Type; };
 112
 113     template<>
 114       struct _To_Unsigned_Type<int>
 115       { typedef unsigned int _Type; };
 116
 117     template<>
 118       struct _To_Unsigned_Type<long>
 119       { typedef unsigned long _Type; };
 120
 121 #ifdef _GLIBCXX_USE_LONG_LONG
 122     template<>
 123       struct _To_Unsigned_Type<long long>
 124       { typedef unsigned long long _Type; };
 125 #endif
 126
 127   } // anonymous namespace
 128
 129
 130   /**
 131    * Seeds the LCR with integral value @p __x0, adjusted so that the
 132    * ring identity is never a member of the convergence set.
 133    */
 134   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 135     void
 136     linear_congruential<_UIntType, __a, __c, __m>::
 137     seed(unsigned long __x0)
 138     {
 139       if ((__mod<_UIntType, 1, 0, __m>(__c) == 0)
 140           && (__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 141         _M_x = __mod<_UIntType, 1, 0, __m>(1);
 142       else
 143         _M_x = __mod<_UIntType, 1, 0, __m>(__x0);
 144     }
 145
 146   /**
 147    * Seeds the LCR engine with a value generated by @p __g.
 148    */
 149   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 150     template<class _Gen>
 151       void
 152       linear_congruential<_UIntType, __a, __c, __m>::
 153       seed(_Gen& __g, false_type)
 154       {
 155         _UIntType __x0 = __g();
 156         if ((__mod<_UIntType, 1, 0, __m>(__c) == 0)
 157             && (__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 158           _M_x = __mod<_UIntType, 1, 0, __m>(1);
 159         else
 160           _M_x = __mod<_UIntType, 1, 0, __m>(__x0);
 161       }
 162
 163   /**
 164    * Returns a value that is less than or equal to all values potentially
 165    * returned by operator(). The return value of this function does not
 166    * change during the lifetime of the object..
 167    *
 168    * The minumum depends on the @p __c parameter: if it is zero, the
 169    * minimum generated must be > 0, otherwise 0 is allowed.
 170    */
 171   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 172     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 173     linear_congruential<_UIntType, __a, __c, __m>::
 174     min() const
 175     { return (__mod<_UIntType, 1, 0, __m>(__c) == 0) ? 1 : 0; }
 176
 177   /**
 178    * Gets the maximum possible value of the generated range.
 179    *
 180    * For a linear congruential generator, the maximum is always @p __m - 1.
 181    */
 182   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 183     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 184     linear_congruential<_UIntType, __a, __c, __m>::
 185     max() const
 186     { return (__m == 0) ? std::numeric_limits<_UIntType>::max() : (__m - 1); }
 187
 188   /**
 189    * Gets the next generated value in sequence.
 190    */
 191   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 192     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 193     linear_congruential<_UIntType, __a, __c, __m>::
 194     operator()()
 195     {
 196       _M_x = __mod<_UIntType, __a, __c, __m>(_M_x);
 197       return _M_x;
 198     }
 199
 200   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 201            typename _CharT, typename _Traits>
 202     std::basic_ostream<_CharT, _Traits>&
 203     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 204                const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
 205     {
 206       const std::ios_base::fmtflags __flags = __os.flags();
 207       const _CharT __fill = __os.fill();
 208       __os.flags(std::ios_base::dec | std::ios_base::fixed
 209                  | std::ios_base::left);
 210       __os.fill(__os.widen(' '));
 211
 212       __os << __lcr._M_x;
 213
 214       __os.flags(__flags);
 215       __os.fill(__fill);
 216       return __os;
 217     }
 218
 219   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 220            typename _CharT, typename _Traits>
 221     std::basic_istream<_CharT, _Traits>&
 222     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 223                linear_congruential<_UIntType, __a, __c, __m>& __lcr)
 224     {
 225       const std::ios_base::fmtflags __flags = __is.flags();
 226       __is.flags(std::ios_base::dec);
 227
 228       __is >> __lcr._M_x;
 229
 230       __is.flags(__flags);
 231       return __is;
 232     }
 233
 234
 235   template<class _UIntType, int __w, int __n, int __m, int __r,
 236            _UIntType __a, int __u, int __s,
 237            _UIntType __b, int __t, _UIntType __c, int __l>
 238     void
 239     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 240                      __b, __t, __c, __l>::
 241     seed(unsigned long __value)
 242     {
 243       _M_x[0] = __mod<_UIntType, 1, 0,
 244         _Shift<_UIntType, __w>::__value>(__value);
 245
 246       for (int __i = 1; __i < state_size; ++__i)
 247         {
 248           _UIntType __x = _M_x[__i - 1];
 249           __x ^= __x >> (__w - 2);
 250           __x *= 1812433253ul;
 251           __x += __i;
 252           _M_x[__i] = __mod<_UIntType, 1, 0,
 253             _Shift<_UIntType, __w>::__value>(__x);
 254         }
 255       _M_p = state_size;
 256     }
 257
 258   template<class _UIntType, int __w, int __n, int __m, int __r,
 259            _UIntType __a, int __u, int __s,
 260            _UIntType __b, int __t, _UIntType __c, int __l>
 261     template<class _Gen>
 262       void
 263       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 264                        __b, __t, __c, __l>::
 265       seed(_Gen& __gen, false_type)
 266       {
 267         for (int __i = 0; __i < state_size; ++__i)
 268           _M_x[__i] = __mod<_UIntType, 1, 0,
 269             _Shift<_UIntType, __w>::__value>(__gen());
 270         _M_p = state_size;
 271       }
 272
 273   template<class _UIntType, int __w, int __n, int __m, int __r,
 274            _UIntType __a, int __u, int __s,
 275            _UIntType __b, int __t, _UIntType __c, int __l>
 276     typename
 277     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 278                      __b, __t, __c, __l>::result_type
 279     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 280                      __b, __t, __c, __l>::
 281     operator()()
 282     {
 283       // Reload the vector - cost is O(n) amortized over n calls.
 284       if (_M_p >= state_size)
 285         {
 286           const _UIntType __upper_mask = (~_UIntType()) << __r;
 287           const _UIntType __lower_mask = ~__upper_mask;
 288           const _UIntType __fx[2] = { 0, __a };
 289
 290           for (int __k = 0; __k < (__n - __m); ++__k)
 291             {
 292               _UIntType __y = ((_M_x[__k] & __upper_mask)
 293                                | (_M_x[__k + 1] & __lower_mask));
 294               _M_x[__k] = _M_x[__k + __m] ^ (__y >> 1) ^ __fx[__y & 0x01];
 295             }
 296
 297           for (int __k = (__n - __m); __k < (__n - 1); ++__k)
 298             {
 299               _UIntType __y = ((_M_x[__k] & __upper_mask)
 300                                | (_M_x[__k + 1] & __lower_mask));
 301               _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
 302                            ^ __fx[__y & 0x01]);
 303             }
 304
 305           _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
 306                            | (_M_x[0] & __lower_mask));
 307           _M_x[__n - 1] = _M_x[__m - 1] ^ (__y >> 1) ^ __fx[__y & 0x01];
 308           _M_p = 0;
 309         }
 310
 311       // Calculate o(x(i)).
 312       result_type __z = _M_x[_M_p++];
 313       __z ^= (__z >> __u);
 314       __z ^= (__z << __s) & __b;
 315       __z ^= (__z << __t) & __c;
 316       __z ^= (__z >> __l);
 317
 318       return __z;
 319     }
 320
 321   template<class _UIntType, int __w, int __n, int __m, int __r,
 322            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
 323            _UIntType __c, int __l,
 324            typename _CharT, typename _Traits>
 325     std::basic_ostream<_CharT, _Traits>&
 326     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 327                const mersenne_twister<_UIntType, __w, __n, __m,
 328                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
 329     {
 330       const std::ios_base::fmtflags __flags = __os.flags();
 331       const _CharT __fill = __os.fill();
 332       const _CharT __space = __os.widen(' ');
 333       __os.flags(std::ios_base::dec | std::ios_base::fixed
 334                  | std::ios_base::left);
 335       __os.fill(__space);
 336
 337       for (int __i = 0; __i < __n - 1; ++__i)
 338         __os << __x._M_x[__i] << __space;
 339       __os << __x._M_x[__n - 1];
 340
 341       __os.flags(__flags);
 342       __os.fill(__fill);
 343       return __os;
 344     }
 345
 346   template<class _UIntType, int __w, int __n, int __m, int __r,
 347            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
 348            _UIntType __c, int __l,
 349            typename _CharT, typename _Traits>
 350     std::basic_istream<_CharT, _Traits>&
 351     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 352                mersenne_twister<_UIntType, __w, __n, __m,
 353                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
 354     {
 355       const std::ios_base::fmtflags __flags = __is.flags();
 356       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 357
 358       for (int __i = 0; __i < __n; ++__i)
 359         __is >> __x._M_x[__i];
 360
 361       __is.flags(__flags);
 362       return __is;
 363     }
 364
 365
 366   template<typename _IntType, _IntType __m, int __s, int __r>
 367     void
 368     subtract_with_carry<_IntType, __m, __s, __r>::
 369     seed(unsigned long __value)
 370     {
 371       if (__value == 0)
 372         __value = 19780503;
 373
 374       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
 375         __lcg(__value);
 376
 377       for (int __i = 0; __i < long_lag; ++__i)
 378         _M_x[__i] = __mod<_IntType, 1, 0, modulus>(__lcg());
 379
 380       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 381       _M_p = 0;
 382     }
 383
 384   template<typename _IntType, _IntType __m, int __s, int __r>
 385     template<class _Gen>
 386       void
 387       subtract_with_carry<_IntType, __m, __s, __r>::
 388       seed(_Gen& __gen, false_type)
 389       {
 390         const int __n = (std::numeric_limits<_IntType>::digits + 31) / 32;
 391
 392         typedef typename _Select<(sizeof(unsigned) == 4),
 393           unsigned, unsigned long>::_Type _UInt32Type;
 394
 395         typedef typename _To_Unsigned_Type<_IntType>::_Type
 396           _UIntType;
 397
 398         for (int __i = 0; __i < long_lag; ++__i)
 399           {
 400             _UIntType __tmp = 0;
 401             _UIntType __factor = 1;
 402             for (int __j = 0; __j < __n; ++__j)
 403               {
 404                 __tmp += (__mod<_UInt32Type, 1, 0, 0>(__gen())
 405                           * __factor);
 406                 __factor *= _Shift<_UIntType, 32>::__value;
 407               }
 408             _M_x[__i] = __mod<_UIntType, 1, 0, modulus>(__tmp);
 409           }
 410         _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 411         _M_p = 0;
 412       }
 413
 414   template<typename _IntType, _IntType __m, int __s, int __r>
 415     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
 416     subtract_with_carry<_IntType, __m, __s, __r>::
 417     operator()()
 418     {
 419       // Derive short lag index from current index.
 420       int __ps = _M_p - short_lag;
 421       if (__ps < 0)
 422         __ps += long_lag;
 423
 424       // Calculate new x(i) without overflow or division.
 425       _IntType __xi;
 426       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
 427         {
 428           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
 429           _M_carry = 0;
 430         }
 431       else
 432         {
 433           __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
 434           _M_carry = 1;
 435         }
 436       _M_x[_M_p++] = __xi;
 437
 438       // Adjust current index to loop around in ring buffer.
 439       if (_M_p >= long_lag)
 440         _M_p = 0;
 441
 442       return __xi;
 443     }
 444
 445   template<typename _IntType, _IntType __m, int __s, int __r,
 446            typename _CharT, typename _Traits>
 447     std::basic_ostream<_CharT, _Traits>&
 448     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 449                const subtract_with_carry<_IntType, __m, __s, __r>& __x)
 450     {
 451       const std::ios_base::fmtflags __flags = __os.flags();
 452       const _CharT __fill = __os.fill();
 453       const _CharT __space = __os.widen(' ');
 454       __os.flags(std::ios_base::dec | std::ios_base::fixed
 455                  | std::ios_base::left);
 456       __os.fill(__space);
 457
 458       for (int __i = 0; __i < __r; ++__i)
 459         __os << __x._M_x[__i] << __space;
 460       __os << __x._M_carry;
 461
 462       __os.flags(__flags);
 463       __os.fill(__fill);
 464       return __os;
 465     }
 466
 467   template<typename _IntType, _IntType __m, int __s, int __r,
 468            typename _CharT, typename _Traits>
 469     std::basic_istream<_CharT, _Traits>&
 470     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 471                subtract_with_carry<_IntType, __m, __s, __r>& __x)
 472     {
 473       const std::ios_base::fmtflags __flags = __is.flags();
 474       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 475
 476       for (int __i = 0; __i < __r; ++__i)
 477         __is >> __x._M_x[__i];
 478       __is >> __x._M_carry;
 479
 480       __is.flags(__flags);
 481       return __is;
 482     }
 483
 484
 485   template<class _UniformRandomNumberGenerator, int __p, int __r>
 486     typename discard_block<_UniformRandomNumberGenerator,
 487                            __p, __r>::result_type
 488     discard_block<_UniformRandomNumberGenerator, __p, __r>::
 489     operator()()
 490     {
 491       if (_M_n >= used_block)
 492         {
 493           while (_M_n < block_size)
 494             {
 495               _M_b();
 496               ++_M_n;
 497             }
 498           _M_n = 0;
 499         }
 500       ++_M_n;
 501       return _M_b();
 502     }
 503
 504   template<class _UniformRandomNumberGenerator, int __p, int __r,
 505            typename _CharT, typename _Traits>
 506     std::basic_ostream<_CharT, _Traits>&
 507     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 508                const discard_block<_UniformRandomNumberGenerator,
 509                __p, __r>& __x)
 510     {
 511       const std::ios_base::fmtflags __flags = __os.flags();
 512       const _CharT __fill = __os.fill();
 513       const _CharT __space = __os.widen(' ');
 514       __os.flags(std::ios_base::dec | std::ios_base::fixed
 515                  | std::ios_base::left);
 516       __os.fill(__space);
 517
 518       __os << __x._M_b << __space << __x._M_n;
 519
 520       __os.flags(__flags);
 521       __os.fill(__fill);
 522       return __os;
 523     }
 524
 525   template<class _UniformRandomNumberGenerator, int __p, int __r,
 526            typename _CharT, typename _Traits>
 527     std::basic_istream<_CharT, _Traits>&
 528     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 529                discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
 530     {
 531       const std::ios_base::fmtflags __flags = __is.flags();
 532       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 533
 534       __is >> __x._M_b >> __x._M_n;
 535
 536       __is.flags(__flags);
 537       return __is;
 538     }
 539
 540
 541   template<class _UniformRandomNumberGenerator1, int __s1,
 542            class _UniformRandomNumberGenerator2, int __s2,
 543            typename _CharT, typename _Traits>
 544     std::basic_ostream<_CharT, _Traits>&
 545     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 546                const xor_combine<_UniformRandomNumberGenerator1, __s1,
 547                _UniformRandomNumberGenerator2, __s2>& __x)
 548     {
 549       const std::ios_base::fmtflags __flags = __os.flags();
 550       const _CharT __fill = __os.fill();
 551       const _CharT __space = __os.widen(' ');
 552       __os.flags(std::ios_base::dec | std::ios_base::fixed
 553                  | std::ios_base::left);
 554       __os.fill(__space);
 555
 556       __os << __x.base1() << __space << __x.base2();
 557
 558       __os.flags(__flags);
 559       __os.fill(__fill);
 560       return __os;
 561     }
 562
 563   template<class _UniformRandomNumberGenerator1, int __s1,
 564            class _UniformRandomNumberGenerator2, int __s2,
 565            typename _CharT, typename _Traits>
 566     std::basic_istream<_CharT, _Traits>&
 567     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 568                xor_combine<_UniformRandomNumberGenerator1, __s1,
 569                _UniformRandomNumberGenerator2, __s2>& __x)
 570     {
 571       const std::ios_base::fmtflags __flags = __is.flags();
 572       __is.flags(std::ios_base::skipws);
 573
 574       __is >> __x._M_b1 >> __x._M_b2;
 575
 576       __is.flags(__flags);
 577       return __is;
 578     }
 579
 580
 581   template<typename _IntType, typename _CharT, typename _Traits>
 582     std::basic_ostream<_CharT, _Traits>&
 583     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 584                const uniform_int<_IntType>& __x)
 585     {
 586       const std::ios_base::fmtflags __flags = __os.flags();
 587       const _CharT __fill = __os.fill();
 588       const _CharT __space = __os.widen(' ');
 589       __os.flags(std::ios_base::scientific | std::ios_base::left);
 590       __os.fill(__space);
 591
 592       __os << __x.min() << __space << __x.max();
 593
 594       __os.flags(__flags);
 595       __os.fill(__fill);
 596       return __os;
 597     }
 598
 599   template<typename _IntType, typename _CharT, typename _Traits>
 600     std::basic_istream<_CharT, _Traits>&
 601     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 602                uniform_int<_IntType>& __x)
 603     {
 604       const std::ios_base::fmtflags __flags = __is.flags();
 605       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 606
 607       __is >> __x._M_min >> __x._M_max;
 608
 609       __is.flags(__flags);
 610       return __is;
 611     }
 612
 613
 614   template<typename _CharT, typename _Traits>
 615     std::basic_ostream<_CharT, _Traits>&
 616     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 617                const bernoulli_distribution& __x)
 618     {
 619       const std::ios_base::fmtflags __flags = __os.flags();
 620       const _CharT __fill = __os.fill();
 621       const std::streamsize __precision = __os.precision();
 622       __os.flags(std::ios_base::scientific | std::ios_base::left);
 623       __os.fill(__os.widen(' '));
 624       __os.precision(_Max_digits10<double>::__value);
 625
 626       __os << __x.p();
 627
 628       __os.flags(__flags);
 629       __os.fill(__fill);
 630       __os.precision(__precision);
 631       return __os;
 632     }
 633
 634
 635   template<typename _IntType, typename _RealType,
 636            typename _CharT, typename _Traits>
 637     std::basic_ostream<_CharT, _Traits>&
 638     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 639                const geometric_distribution<_IntType, _RealType>& __x)
 640     {
 641       const std::ios_base::fmtflags __flags = __os.flags();
 642       const _CharT __fill = __os.fill();
 643       const std::streamsize __precision = __os.precision();
 644       __os.flags(std::ios_base::scientific | std::ios_base::left);
 645       __os.fill(__os.widen(' '));
 646       __os.precision(_Max_digits10<_RealType>::__value);
 647
 648       __os << __x.p();
 649
 650       __os.flags(__flags);
 651       __os.fill(__fill);
 652       __os.precision(__precision);
 653       return __os;
 654     }
 655
 656
 657   template<typename _IntType, typename _RealType>
 658     poisson_distribution<_IntType, _RealType>::
 659     poisson_distribution(const _RealType& __mean)
 660     : _M_mean(__mean), _M_large(false)
 661     {
 662       _GLIBCXX_DEBUG_ASSERT(_M_mean > 0.0);
 663
 664 #if _GLIBCXX_USE_C99_MATH_TR1
 665       if (_M_mean >= 12)
 666         {
 667           _M_large = true;
 668           const _RealType __m = std::floor(_M_mean);
 669           _M_lm_thr = std::log(_M_mean);
 670           _M_lfm = std::tr1::lgamma(__m + 1);
 671           _M_sm = std::sqrt(__m);
 672
 673           const _RealType __dx =
 674             std::sqrt(2 * __m
 675             * std::log(_RealType(40.743665431525205956834243423363677L)
 676                        * __m));
 677           _M_d = std::tr1::round(std::max(_RealType(6),
 678                                           std::min(__m, __dx)));
 679           const _RealType __cx = 2 * (2 * __m + _M_d);
 680           const _RealType __cx4 = __cx / 4;
 681           _M_scx4 = std::sqrt(__cx4);
 682           _M_2cx = 2 / __cx;
 683
 684           const _RealType __pi_2 = 1.5707963267948966192313216916397514L;
 685           _M_c2b = std::sqrt(__pi_2 * __cx4) * std::exp(_M_2cx);
 686           _M_cb = __cx * std::exp(-_M_d * _M_2cx * (1 + _M_d / 2)) / _M_d;
 687         }
 688       else
 689 #endif
 690         _M_lm_thr = std::exp(-_M_mean);
 691       }
 692
 693   /**
 694    * A rejection algorithm when mean >= 12 and a simple method based
 695    * upon the multiplication of uniform random variates otherwise.
 696    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
 697    * is defined.
 698    *
 699    * Reference:
 700    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
 701    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
 702    */
 703   template<typename _IntType, typename _RealType>
 704     template<class _UniformRandomNumberGenerator>
 705       typename poisson_distribution<_IntType, _RealType>::result_type
 706       poisson_distribution<_IntType, _RealType>::
 707       operator()(_UniformRandomNumberGenerator& __urng)
 708       {
 709 #if _GLIBCXX_USE_C99_MATH_TR1
 710         if (_M_large)
 711           {
 712             _RealType __x;
 713
 714             const _RealType __m = std::floor(_M_mean);
 715             // sqrt(mu * pi / 2)
 716             const _RealType __c1 = (_M_sm
 717                                     * 1.2533141373155002512078826424055226L);
 718             const _RealType __c2 = _M_c2b + __c1;
 719             const _RealType __c3 = __c2 + 1;
 720             const _RealType __c4 = __c3 + 1;
 721             // c4 + e^(1 / 78)
 722             const _RealType __c5 = (__c4
 723                                     + 1.0129030479320018583185514777512983L);
 724             const _RealType __c = _M_cb + __c5;
 725             const _RealType __cx = 2 * (2 * __m + _M_d);
 726
 727             normal_distribution<_RealType> __nd;
 728
 729             bool __keepgoing = true;
 730             do
 731               {
 732                 const _RealType __u = __c * __urng();
 733                 const _RealType __e = -std::log(__urng());
 734
 735                 _RealType __w = 0.0;
 736
 737                 if (__u <= __c1)
 738                   {
 739                     const _RealType __n = __nd(__urng);
 740                     const _RealType __y = -std::abs(__n) * _M_sm - 1;
 741                     __x = std::floor(__y);
 742                     __w = -__n * __n / 2;
 743                     if (__x < -__m)
 744                       continue;
 745                   }
 746                 else if (__u <= __c2)
 747                   {
 748                     const _RealType __n = __nd(__urng);
 749                     const _RealType __y = 1 + std::abs(__n) * _M_scx4;
 750                     __x = std::ceil(__y);
 751                     __w = __y * (2 - __y) * _M_2cx;
 752                     if (__x > _M_d)
 753                       continue;
 754                   }
 755                 else if (__u <= __c3)
 756                   // XXX This case not in the book, nor in the Errata...
 757                   __x = -1;
 758                 else if (__u <= __c4)
 759                   __x = 0;
 760                 else if (__u <= __c5)
 761                   __x = 1;
 762                 else
 763                   {
 764                     const _RealType __v = -std::log(__urng());
 765                     const _RealType __y = _M_d + __v * __cx / _M_d;
 766                     __x = std::ceil(__y);
 767                     __w = -_M_d * _M_2cx * (1 + __y / 2);
 768                   }
 769
 770                 __keepgoing = (__w - __e - __x * _M_lm_thr
 771                                > _M_lfm - std::tr1::lgamma(__x + __m + 1));
 772
 773               } while (__keepgoing);
 774
 775             return _IntType(std::tr1::round(__x + __m));
 776           }
 777         else
 778 #endif
 779           {
 780             _IntType __x = -1;
 781             _RealType __prod = 1.0;
 782
 783             do
 784               {
 785                 __prod *= __urng();
 786                 __x += 1;
 787               }
 788             while (__prod > _M_lm_thr);
 789
 790             return __x;
 791           }
 792       }
 793
 794   template<typename _IntType, typename _RealType,
 795            typename _CharT, typename _Traits>
 796     std::basic_ostream<_CharT, _Traits>&
 797     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 798                const poisson_distribution<_IntType, _RealType>& __x)
 799     {
 800       const std::ios_base::fmtflags __flags = __os.flags();
 801       const _CharT __fill = __os.fill();
 802       const std::streamsize __precision = __os.precision();
 803       const _CharT __space = __os.widen(' ');
 804       __os.flags(std::ios_base::scientific | std::ios_base::left);
 805       __os.fill(__space);
 806       __os.precision(_Max_digits10<_RealType>::__value);
 807
 808       __os << __x._M_large << __space << __x.mean()
 809            << __space << __x._M_lm_thr;
 810 #if _GLIBCXX_USE_C99_MATH_TR1
 811       if (__x._M_large)
 812         __os << __space << __x._M_lfm << __space << __x._M_sm
 813              << __space << __x._M_d << __space << __x._M_scx4
 814              << __space << __x._M_2cx << __space << __x._M_c2b
 815              << __space << __x._M_cb;
 816 #endif
 817
 818       __os.flags(__flags);
 819       __os.fill(__fill);
 820       __os.precision(__precision);
 821       return __os;
 822     }
 823
 824   template<typename _IntType, typename _RealType,
 825            typename _CharT, typename _Traits>
 826     std::basic_istream<_CharT, _Traits>&
 827     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 828                poisson_distribution<_IntType, _RealType>& __x)
 829     {
 830       const std::ios_base::fmtflags __flags = __is.flags();
 831       __is.flags(std::ios_base::skipws);
 832
 833       __is >> __x._M_large >> __x._M_mean >> __x._M_lm_thr;
 834 #if _GLIBCXX_USE_C99_MATH_TR1
 835       if (__x._M_large)
 836         __is >> __x._M_lfm >> __x._M_sm >> __x._M_d >> __x._M_scx4
 837              >> __x._M_2cx >> __x._M_c2b >> __x._M_cb;
 838 #endif
 839
 840       __is.flags(__flags);
 841       return __is;
 842     }
 843
 844
 845   template<typename _RealType, typename _CharT, typename _Traits>
 846     std::basic_ostream<_CharT, _Traits>&
 847     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 848                const uniform_real<_RealType>& __x)
 849     {
 850       const std::ios_base::fmtflags __flags = __os.flags();
 851       const _CharT __fill = __os.fill();
 852       const std::streamsize __precision = __os.precision();
 853       const _CharT __space = __os.widen(' ');
 854       __os.flags(std::ios_base::scientific | std::ios_base::left);
 855       __os.fill(__space);
 856       __os.precision(_Max_digits10<_RealType>::__value);
 857
 858       __os << __x.min() << __space << __x.max();
 859
 860       __os.flags(__flags);
 861       __os.fill(__fill);
 862       __os.precision(__precision);
 863       return __os;
 864     }
 865
 866   template<typename _RealType, typename _CharT, typename _Traits>
 867     std::basic_istream<_CharT, _Traits>&
 868     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 869                uniform_real<_RealType>& __x)
 870     {
 871       const std::ios_base::fmtflags __flags = __is.flags();
 872       __is.flags(std::ios_base::skipws);
 873
 874       __is >> __x._M_min >> __x._M_max;
 875
 876       __is.flags(__flags);
 877       return __is;
 878     }
 879
 880
 881   template<typename _RealType, typename _CharT, typename _Traits>
 882     std::basic_ostream<_CharT, _Traits>&
 883     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 884                const exponential_distribution<_RealType>& __x)
 885     {
 886       const std::ios_base::fmtflags __flags = __os.flags();
 887       const _CharT __fill = __os.fill();
 888       const std::streamsize __precision = __os.precision();
 889       __os.flags(std::ios_base::scientific | std::ios_base::left);
 890       __os.fill(__os.widen(' '));
 891       __os.precision(_Max_digits10<_RealType>::__value);
 892
 893       __os << __x.lambda();
 894
 895       __os.flags(__flags);
 896       __os.fill(__fill);
 897       __os.precision(__precision);
 898       return __os;
 899     }
 900
 901
 902   /**
 903    * Classic Box-Muller method.
 904    *
 905    * Reference:
 906    * Box, G. E. P. and Muller, M. E. "A Note on the Generation of
 907    * Random Normal Deviates." Ann. Math. Stat. 29, 610-611, 1958.
 908    */
 909   template<typename _RealType>
 910     template<class _UniformRandomNumberGenerator>
 911       typename normal_distribution<_RealType>::result_type
 912       normal_distribution<_RealType>::
 913       operator()(_UniformRandomNumberGenerator& __urng)
 914       {
 915         result_type __ret;
 916
 917         if (_M_saved_available)
 918           {
 919             _M_saved_available = false;
 920             __ret = _M_saved;
 921           }
 922         else
 923           {
 924             result_type __x, __y, __r2;
 925             do
 926               {
 927                 __x = result_type(2.0) * __urng() - result_type(1.0);
 928                 __y = result_type(2.0) * __urng() - result_type(1.0);
 929                 __r2 = __x * __x + __y * __y;
 930               }
 931             while (__r2 > result_type(1.0) || __r2 == result_type(0));
 932
 933             const result_type __mult = std::sqrt(-result_type(2.0)
 934                                                  * std::log(__r2) / __r2);
 935             _M_saved = __x * __mult;
 936             _M_saved_available = true;
 937             __ret = __y * __mult;
 938           }
 939
 940         return __ret * _M_sigma + _M_mean;
 941       }
 942
 943   template<typename _RealType, typename _CharT, typename _Traits>
 944     std::basic_ostream<_CharT, _Traits>&
 945     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 946                const normal_distribution<_RealType>& __x)
 947     {
 948       const std::ios_base::fmtflags __flags = __os.flags();
 949       const _CharT __fill = __os.fill();
 950       const std::streamsize __precision = __os.precision();
 951       const _CharT __space = __os.widen(' ');
 952       __os.flags(std::ios_base::scientific | std::ios_base::left);
 953       __os.fill(__space);
 954       __os.precision(_Max_digits10<_RealType>::__value);
 955
 956       __os << __x._M_saved_available << __space
 957            << __x.mean() << __space
 958            << __x.sigma();
 959       if (__x._M_saved_available)
 960         __os << __space << __x._M_saved;
 961
 962       __os.flags(__flags);
 963       __os.fill(__fill);
 964       __os.precision(__precision);
 965       return __os;
 966     }
 967
 968   template<typename _RealType, typename _CharT, typename _Traits>
 969     std::basic_istream<_CharT, _Traits>&
 970     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 971                normal_distribution<_RealType>& __x)
 972     {
 973       const std::ios_base::fmtflags __flags = __is.flags();
 974       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 975
 976       __is >> __x._M_saved_available >> __x._M_mean
 977            >> __x._M_sigma;
 978       if (__x._M_saved_available)
 979         __is >> __x._M_saved;
 980
 981       __is.flags(__flags);
 982       return __is;
 983     }
 984
 985
 986   /**
 987    * Cheng's rejection algorithm GB for alpha >= 1 and a modification
 988    * of Vaduva's rejection from Weibull algorithm due to Devroye for
 989    * alpha < 1.
 990    *
 991    * References:
 992    * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
 993    * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
 994    *
 995    * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
 996    * and Composition Procedures." Math. Operationsforschung and Statistik,
 997    * Series in Statistics, 8, 545-576, 1977.
 998    *
 999    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1000    * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
1001    */
1002   template<typename _RealType>
1003     template<class _UniformRandomNumberGenerator>
1004       typename gamma_distribution<_RealType>::result_type
1005       gamma_distribution<_RealType>::
1006       operator()(_UniformRandomNumberGenerator& __urng)
1007       {
1008         result_type __x;
1009
1010         if (_M_alpha >= 1)
1011           {
1012             // alpha - log(4)
1013             const result_type __b = _M_alpha
1014               - result_type(1.3862943611198906188344642429163531L);
1015             const result_type __l = std::sqrt(2 * _M_alpha - 1);
1016             const result_type __c = _M_alpha + __l;
1017             const result_type __1l = 1 / __l;
1018
1019             // 1 + log(9 / 2)
1020             const result_type __k = 2.5040773967762740733732583523868748L;
1021
1022             result_type __z, __r;
1023             do
1024               {
1025                 const result_type __u = __urng();
1026                 const result_type __v = __urng();
1027
1028                 const result_type __y = __1l * std::log(__v / (1 - __v));
1029                 __x = _M_alpha * std::exp(__y);
1030
1031                 __z = __u * __v * __v;
1032                 __r = __b + __c * __y - __x;
1033               }
1034             while (__r < result_type(4.5) * __z - __k
1035                    && __r < std::log(__z));
1036           }
1037         else
1038           {
1039             const result_type __c = 1 / _M_alpha;
1040             const result_type __d =
1041               std::pow(_M_alpha, _M_alpha / (1 - _M_alpha)) * (1 - _M_alpha);
1042
1043             result_type __z, __e;
1044             do
1045               {
1046                 __z = -std::log(__urng());
1047                 __e = -std::log(__urng());
1048
1049                 __x = std::pow(__z, __c);
1050               }
1051             while (__z + __e < __d + __x);
1052           }
1053
1054         return __x;
1055       }
1056
1057   template<typename _RealType, typename _CharT, typename _Traits>
1058     std::basic_ostream<_CharT, _Traits>&
1059     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1060                const gamma_distribution<_RealType>& __x)
1061     {
1062       const std::ios_base::fmtflags __flags = __os.flags();
1063       const _CharT __fill = __os.fill();
1064       const std::streamsize __precision = __os.precision();
1065       __os.flags(std::ios_base::scientific | std::ios_base::left);
1066       __os.fill(__os.widen(' '));
1067       __os.precision(_Max_digits10<_RealType>::__value);
1068
1069       __os << __x.alpha();
1070
1071       __os.flags(__flags);
1072       __os.fill(__fill);
1073       __os.precision(__precision);
1074       return __os;
1075     }
1076
1077 _GLIBCXX_END_NAMESPACE
1078 }