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 _Private
  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   } // namespace _Private
 106
 107
 108   /**
 109    * Seeds the LCR with integral value @p __x0, adjusted so that the
 110    * ring identity is never a member of the convergence set.
 111    */
 112   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 113     void
 114     linear_congruential<_UIntType, __a, __c, __m>::
 115     seed(unsigned long __x0)
 116     {
 117       if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
 118           && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 119         _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
 120       else
 121         _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
 122     }
 123
 124   /**
 125    * Seeds the LCR engine with a value generated by @p __g.
 126    */
 127   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 128     template<class _Gen>
 129       void
 130       linear_congruential<_UIntType, __a, __c, __m>::
 131       seed(_Gen& __g, false_type)
 132       {
 133         _UIntType __x0 = __g();
 134         if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
 135             && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 136           _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
 137         else
 138           _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
 139       }
 140
 141   /**
 142    * Returns a value that is less than or equal to all values potentially
 143    * returned by operator(). The return value of this function does not
 144    * change during the lifetime of the object..
 145    *
 146    * The minumum depends on the @p __c parameter: if it is zero, the
 147    * minimum generated must be > 0, otherwise 0 is allowed.
 148    */
 149   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 150     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 151     linear_congruential<_UIntType, __a, __c, __m>::
 152     min() const
 153     { return (_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0) ? 1 : 0; }
 154
 155   /**
 156    * Gets the maximum possible value of the generated range.
 157    *
 158    * For a linear congruential generator, the maximum is always @p __m - 1.
 159    */
 160   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 161     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 162     linear_congruential<_UIntType, __a, __c, __m>::
 163     max() const
 164     { return (__m == 0) ? std::numeric_limits<_UIntType>::max() : (__m - 1); }
 165
 166   /**
 167    * Gets the next generated value in sequence.
 168    */
 169   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 170     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 171     linear_congruential<_UIntType, __a, __c, __m>::
 172     operator()()
 173     {
 174       _M_x = _Private::__mod<_UIntType, __a, __c, __m>(_M_x);
 175       return _M_x;
 176     }
 177
 178   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 179            typename _CharT, typename _Traits>
 180     std::basic_ostream<_CharT, _Traits>&
 181     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 182                const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
 183     {
 184       const std::ios_base::fmtflags __flags = __os.flags();
 185       const _CharT __fill = __os.fill();
 186       __os.flags(std::ios_base::dec | std::ios_base::fixed
 187                  | std::ios_base::left);
 188       __os.fill(__os.widen(' '));
 189
 190       __os << __lcr._M_x;
 191
 192       __os.flags(__flags);
 193       __os.fill(__fill);
 194       return __os;
 195     }
 196
 197   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 198            typename _CharT, typename _Traits>
 199     std::basic_istream<_CharT, _Traits>&
 200     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 201                linear_congruential<_UIntType, __a, __c, __m>& __lcr)
 202     {
 203       const std::ios_base::fmtflags __flags = __is.flags();
 204       __is.flags(std::ios_base::dec);
 205
 206       __is >> __lcr._M_x;
 207
 208       __is.flags(__flags);
 209       return __is;
 210     }
 211
 212
 213   template<class _UIntType, int __w, int __n, int __m, int __r,
 214            _UIntType __a, int __u, int __s,
 215            _UIntType __b, int __t, _UIntType __c, int __l>
 216     void
 217     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 218                      __b, __t, __c, __l>::
 219     seed(unsigned long __value)
 220     {
 221       _M_x[0] = _Private::__mod<_UIntType, 1, 0,
 222         _Private::_Shift<_UIntType, __w>::__value>(__value);
 223
 224       for (int __i = 1; __i < state_size; ++__i)
 225         {
 226           _UIntType __x = _M_x[__i - 1];
 227           __x ^= __x >> (__w - 2);
 228           __x *= 1812433253ul;
 229           __x += __i;
 230           _M_x[__i] = _Private::__mod<_UIntType, 1, 0,
 231             _Private::_Shift<_UIntType, __w>::__value>(__x);
 232         }
 233       _M_p = state_size;
 234     }
 235
 236   template<class _UIntType, int __w, int __n, int __m, int __r,
 237            _UIntType __a, int __u, int __s,
 238            _UIntType __b, int __t, _UIntType __c, int __l>
 239     template<class _Gen>
 240       void
 241       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 242                        __b, __t, __c, __l>::
 243       seed(_Gen& __gen, false_type)
 244       {
 245         for (int __i = 0; __i < state_size; ++__i)
 246           _M_x[__i] = _Private::__mod<_UIntType, 1, 0,
 247             _Private::_Shift<_UIntType, __w>::__value>(__gen());
 248         _M_p = state_size;
 249       }
 250
 251   template<class _UIntType, int __w, int __n, int __m, int __r,
 252            _UIntType __a, int __u, int __s,
 253            _UIntType __b, int __t, _UIntType __c, int __l>
 254     typename
 255     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 256                      __b, __t, __c, __l>::result_type
 257     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 258                      __b, __t, __c, __l>::
 259     operator()()
 260     {
 261       // Reload the vector - cost is O(n) amortized over n calls.
 262       if (_M_p >= state_size)
 263         {
 264           const _UIntType __upper_mask = (~_UIntType()) << __r;
 265           const _UIntType __lower_mask = ~__upper_mask;
 266
 267           for (int __k = 0; __k < (__n - __m); ++__k)
 268             {
 269               _UIntType __y = ((_M_x[__k] & __upper_mask)
 270                                | (_M_x[__k + 1] & __lower_mask));
 271               _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
 272                            ^ ((__y & 0x01) ? __a : 0));
 273             }
 274
 275           for (int __k = (__n - __m); __k < (__n - 1); ++__k)
 276             {
 277               _UIntType __y = ((_M_x[__k] & __upper_mask)
 278                                | (_M_x[__k + 1] & __lower_mask));
 279               _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
 280                            ^ ((__y & 0x01) ? __a : 0));
 281             }
 282
 283           _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
 284                            | (_M_x[0] & __lower_mask));
 285           _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
 286                            ^ ((__y & 0x01) ? __a : 0));
 287           _M_p = 0;
 288         }
 289
 290       // Calculate o(x(i)).
 291       result_type __z = _M_x[_M_p++];
 292       __z ^= (__z >> __u);
 293       __z ^= (__z << __s) & __b;
 294       __z ^= (__z << __t) & __c;
 295       __z ^= (__z >> __l);
 296
 297       return __z;
 298     }
 299
 300   template<class _UIntType, int __w, int __n, int __m, int __r,
 301            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
 302            _UIntType __c, int __l,
 303            typename _CharT, typename _Traits>
 304     std::basic_ostream<_CharT, _Traits>&
 305     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 306                const mersenne_twister<_UIntType, __w, __n, __m,
 307                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
 308     {
 309       const std::ios_base::fmtflags __flags = __os.flags();
 310       const _CharT __fill = __os.fill();
 311       const _CharT __space = __os.widen(' ');
 312       __os.flags(std::ios_base::dec | std::ios_base::fixed
 313                  | std::ios_base::left);
 314       __os.fill(__space);
 315
 316       for (int __i = 0; __i < __n - 1; ++__i)
 317         __os << __x._M_x[__i] << __space;
 318       __os << __x._M_x[__n - 1];
 319
 320       __os.flags(__flags);
 321       __os.fill(__fill);
 322       return __os;
 323     }
 324
 325   template<class _UIntType, int __w, int __n, int __m, int __r,
 326            _UIntType __a, int __u, int __s, _UIntType __b, int __t,
 327            _UIntType __c, int __l,
 328            typename _CharT, typename _Traits>
 329     std::basic_istream<_CharT, _Traits>&
 330     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 331                mersenne_twister<_UIntType, __w, __n, __m,
 332                __r, __a, __u, __s, __b, __t, __c, __l>& __x)
 333     {
 334       const std::ios_base::fmtflags __flags = __is.flags();
 335       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 336
 337       for (int __i = 0; __i < __n; ++__i)
 338         __is >> __x._M_x[__i];
 339
 340       __is.flags(__flags);
 341       return __is;
 342     }
 343
 344
 345   template<typename _IntType, _IntType __m, int __s, int __r>
 346     void
 347     subtract_with_carry<_IntType, __m, __s, __r>::
 348     seed(unsigned long __value)
 349     {
 350       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
 351         __lcg(__value);
 352
 353       for (int __i = 0; __i < long_lag; ++__i)
 354         _M_x[__i] = _Private::__mod<_IntType, 1, 0, modulus>(__lcg());
 355
 356       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 357       _M_p = 0;
 358     }
 359
 360   //
 361   // This implementation differs from the tr1 spec because the tr1 spec refused
 362   // to make any sense to me:  the exponent of the factor in the spec goes from
 363   // 1 to (n-1), but it would only make sense to me if it went from 0 to (n-1).
 364   //
 365   // This algorithm is still problematic because it can overflow left right and
 366   // center.
 367   //
 368   template<typename _IntType, _IntType __m, int __s, int __r>
 369     template<class _Gen>
 370     void
 371     subtract_with_carry<_IntType, __m, __s, __r>::
 372     seed(_Gen& __gen, false_type)
 373     {
 374       const int __n = (std::numeric_limits<_IntType>::digits + 31) / 32;
 375       for (int __i = 0; __i < long_lag; ++__i)
 376         {
 377           _M_x[__i] = 0;
 378           unsigned long __factor = 1;
 379           for (int __j = 0; __j < __n; ++__j)
 380             {
 381               _M_x[__i] += __gen() * __factor;
 382               __factor *= 0x80000000;
 383             }
 384           _M_x[__i] = _Private::__mod<_IntType, 1, 0, modulus>(_M_x[__i]);
 385         }
 386       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 387       _M_p = 0;
 388     }
 389
 390   template<typename _IntType, _IntType __m, int __s, int __r>
 391     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
 392     subtract_with_carry<_IntType, __m, __s, __r>::
 393     operator()()
 394     {
 395       // Derive short lag index from current index.
 396       int __ps = _M_p - short_lag;
 397       if (__ps < 0)
 398         __ps += long_lag;
 399
 400       // Calculate new x(i) without overflow or division.
 401       _IntType __xi;
 402       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
 403         {
 404           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
 405           _M_carry = 0;
 406         }
 407       else
 408         {
 409           __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
 410           _M_carry = 1;
 411         }
 412       _M_x[_M_p++] = __xi;
 413
 414       // Adjust current index to loop around in ring buffer.
 415       if (_M_p >= long_lag)
 416         _M_p = 0;
 417
 418       return __xi;
 419     }
 420
 421   template<typename _IntType, _IntType __m, int __s, int __r,
 422            typename _CharT, typename _Traits>
 423     std::basic_ostream<_CharT, _Traits>&
 424     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 425                const subtract_with_carry<_IntType, __m, __s, __r>& __x)
 426     {
 427       const std::ios_base::fmtflags __flags = __os.flags();
 428       const _CharT __fill = __os.fill();
 429       const _CharT __space = __os.widen(' ');
 430       __os.flags(std::ios_base::dec | std::ios_base::fixed
 431                  | std::ios_base::left);
 432       __os.fill(__space);
 433
 434       for (int __i = 0; __i < __r; ++__i)
 435         __os << __x._M_x[__i] << __space;
 436       __os << __x._M_carry;
 437
 438       __os.flags(__flags);
 439       __os.fill(__fill);
 440       return __os;
 441     }
 442
 443   template<typename _IntType, _IntType __m, int __s, int __r,
 444            typename _CharT, typename _Traits>
 445     std::basic_istream<_CharT, _Traits>&
 446     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 447                subtract_with_carry<_IntType, __m, __s, __r>& __x)
 448     {
 449       const std::ios_base::fmtflags __flags = __is.flags();
 450       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 451
 452       for (int __i = 0; __i < __r; ++__i)
 453         __is >> __x._M_x[__i];
 454       __is >> __x._M_carry;
 455
 456       __is.flags(__flags);
 457       return __is;
 458     }
 459
 460
 461   template<class _UniformRandomNumberGenerator, int __p, int __r>
 462     typename discard_block<_UniformRandomNumberGenerator,
 463                            __p, __r>::result_type
 464     discard_block<_UniformRandomNumberGenerator, __p, __r>::
 465     operator()()
 466     {
 467       if (_M_n >= used_block)
 468         {
 469           while (_M_n < block_size)
 470             {
 471               _M_b();
 472               ++_M_n;
 473             }
 474           _M_n = 0;
 475         }
 476       ++_M_n;
 477       return _M_b();
 478     }
 479
 480   template<class _UniformRandomNumberGenerator, int __p, int __r,
 481            typename _CharT, typename _Traits>
 482     std::basic_ostream<_CharT, _Traits>&
 483     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 484                const discard_block<_UniformRandomNumberGenerator,
 485                __p, __r>& __x)
 486     {
 487       const std::ios_base::fmtflags __flags = __os.flags();
 488       const _CharT __fill = __os.fill();
 489       const _CharT __space = __os.widen(' ');
 490       __os.flags(std::ios_base::dec | std::ios_base::fixed
 491                  | std::ios_base::left);
 492       __os.fill(__space);
 493
 494       __os << __x._M_b << __space << __x._M_n;
 495
 496       __os.flags(__flags);
 497       __os.fill(__fill);
 498       return __os;
 499     }
 500
 501   template<class _UniformRandomNumberGenerator, int __p, int __r,
 502            typename _CharT, typename _Traits>
 503     std::basic_istream<_CharT, _Traits>&
 504     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 505                discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
 506     {
 507       const std::ios_base::fmtflags __flags = __is.flags();
 508       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 509
 510       __is >> __x._M_b >> __x._M_n;
 511
 512       __is.flags(__flags);
 513       return __is;
 514     }
 515
 516
 517   template<class _UniformRandomNumberGenerator1, int __s1,
 518            class _UniformRandomNumberGenerator2, int __s2,
 519            typename _CharT, typename _Traits>
 520     std::basic_ostream<_CharT, _Traits>&
 521     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 522                const xor_combine<_UniformRandomNumberGenerator1, __s1,
 523                _UniformRandomNumberGenerator2, __s2>& __x)
 524     {
 525       const std::ios_base::fmtflags __flags = __os.flags();
 526       const _CharT __fill = __os.fill();
 527       const _CharT __space = __os.widen(' ');
 528       __os.flags(std::ios_base::dec | std::ios_base::fixed
 529                  | std::ios_base::left);
 530       __os.fill(__space);
 531
 532       __os << __x.base1() << __space << __x.base2();
 533
 534       __os.flags(__flags);
 535       __os.fill(__fill);
 536       return __os;
 537     }
 538
 539   template<class _UniformRandomNumberGenerator1, int __s1,
 540            class _UniformRandomNumberGenerator2, int __s2,
 541            typename _CharT, typename _Traits>
 542     std::basic_istream<_CharT, _Traits>&
 543     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 544                xor_combine<_UniformRandomNumberGenerator1, __s1,
 545                _UniformRandomNumberGenerator2, __s2>& __x)
 546     {
 547       const std::ios_base::fmtflags __flags = __is.flags();
 548       __is.flags(std::ios_base::skipws);
 549
 550       __is >> __x._M_b1 >> __x._M_b2;
 551
 552       __is.flags(__flags);
 553       return __is;
 554     }
 555
 556
 557   template<typename _IntType, typename _CharT, typename _Traits>
 558     std::basic_ostream<_CharT, _Traits>&
 559     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 560                const uniform_int<_IntType>& __x)
 561     {
 562       const std::ios_base::fmtflags __flags = __os.flags();
 563       const _CharT __fill = __os.fill();
 564       const _CharT __space = __os.widen(' ');
 565       __os.flags(std::ios_base::scientific | std::ios_base::left);
 566       __os.fill(__space);
 567
 568       __os << __x.min() << __space << __x.max();
 569
 570       __os.flags(__flags);
 571       __os.fill(__fill);
 572       return __os;
 573     }
 574
 575   template<typename _IntType, typename _CharT, typename _Traits>
 576     std::basic_istream<_CharT, _Traits>&
 577     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 578                uniform_int<_IntType>& __x)
 579     {
 580       const std::ios_base::fmtflags __flags = __is.flags();
 581       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 582
 583       __is >> __x._M_min >> __x._M_max;
 584
 585       __is.flags(__flags);
 586       return __is;
 587     }
 588
 589
 590   template<typename _CharT, typename _Traits>
 591     std::basic_ostream<_CharT, _Traits>&
 592     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 593                const bernoulli_distribution& __x)
 594     {
 595       const std::ios_base::fmtflags __flags = __os.flags();
 596       const _CharT __fill = __os.fill();
 597       const std::streamsize __precision = __os.precision();
 598       __os.flags(std::ios_base::scientific | std::ios_base::left);
 599       __os.fill(__os.widen(' '));
 600       __os.precision(_Private::_Max_digits10<double>::__value);
 601
 602       __os << __x.p();
 603
 604       __os.flags(__flags);
 605       __os.fill(__fill);
 606       __os.precision(__precision);
 607       return __os;
 608     }
 609
 610
 611   template<typename _IntType, typename _RealType,
 612            typename _CharT, typename _Traits>
 613     std::basic_ostream<_CharT, _Traits>&
 614     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 615                const geometric_distribution<_IntType, _RealType>& __x)
 616     {
 617       const std::ios_base::fmtflags __flags = __os.flags();
 618       const _CharT __fill = __os.fill();
 619       const std::streamsize __precision = __os.precision();
 620       __os.flags(std::ios_base::scientific | std::ios_base::left);
 621       __os.fill(__os.widen(' '));
 622       __os.precision(_Private::_Max_digits10<_RealType>::__value);
 623
 624       __os << __x.p();
 625
 626       __os.flags(__flags);
 627       __os.fill(__fill);
 628       __os.precision(__precision);
 629       return __os;
 630     }
 631
 632
 633   template<typename _RealType, typename _CharT, typename _Traits>
 634     std::basic_ostream<_CharT, _Traits>&
 635     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 636                const uniform_real<_RealType>& __x)
 637     {
 638       const std::ios_base::fmtflags __flags = __os.flags();
 639       const _CharT __fill = __os.fill();
 640       const std::streamsize __precision = __os.precision();
 641       const _CharT __space = __os.widen(' ');
 642       __os.flags(std::ios_base::scientific | std::ios_base::left);
 643       __os.fill(__space);
 644       __os.precision(_Private::_Max_digits10<_RealType>::__value);
 645
 646       __os << __x.min() << __space << __x.max();
 647
 648       __os.flags(__flags);
 649       __os.fill(__fill);
 650       __os.precision(__precision);
 651       return __os;
 652     }
 653
 654   template<typename _RealType, typename _CharT, typename _Traits>
 655     std::basic_istream<_CharT, _Traits>&
 656     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 657                uniform_real<_RealType>& __x)
 658     {
 659       const std::ios_base::fmtflags __flags = __is.flags();
 660       __is.flags(std::ios_base::skipws);
 661
 662       __is >> __x._M_min >> __x._M_max;
 663
 664       __is.flags(__flags);
 665       return __is;
 666     }
 667
 668
 669   template<typename _RealType, typename _CharT, typename _Traits>
 670     std::basic_ostream<_CharT, _Traits>&
 671     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 672                const exponential_distribution<_RealType>& __x)
 673     {
 674       const std::ios_base::fmtflags __flags = __os.flags();
 675       const _CharT __fill = __os.fill();
 676       const std::streamsize __precision = __os.precision();
 677       __os.flags(std::ios_base::scientific | std::ios_base::left);
 678       __os.fill(__os.widen(' '));
 679       __os.precision(_Private::_Max_digits10<_RealType>::__value);
 680
 681       __os << __x.lambda();
 682
 683       __os.flags(__flags);
 684       __os.fill(__fill);
 685       __os.precision(__precision);
 686       return __os;
 687     }
 688
 689
 690   /**
 691    * Classic Box-Muller method.
 692    *
 693    * Reference:
 694    * Box, G. E. P. and Muller, M. E. "A Note on the Generation of
 695    * Random Normal Deviates." Ann. Math. Stat. 29, 610-611, 1958.
 696    */
 697   template<typename _RealType>
 698     template<class _UniformRandomNumberGenerator>
 699       typename normal_distribution<_RealType>::result_type
 700       normal_distribution<_RealType>::
 701       operator()(_UniformRandomNumberGenerator& __urng)
 702       {
 703         result_type __ret;
 704
 705         if (_M_saved_available)
 706           {
 707             _M_saved_available = false;
 708             __ret = _M_saved;
 709           }
 710         else
 711           {
 712             result_type __x, __y, __r2;
 713             do
 714               {
 715                 __x = result_type(2.0) * __urng() - result_type(1.0);
 716                 __y = result_type(2.0) * __urng() - result_type(1.0);
 717                 __r2 = __x * __x + __y * __y;
 718               }
 719             while (__r2 > result_type(1.0) || __r2 == result_type(0));
 720
 721             const result_type __mult = std::sqrt(-result_type(2.0)
 722                                                  * std::log(__r2) / __r2);
 723             _M_saved = __x * __mult;
 724             _M_saved_available = true;
 725             __ret = __y * __mult;
 726           }
 727
 728         return __ret * _M_sigma + _M_mean;
 729       }
 730
 731   template<typename _RealType, typename _CharT, typename _Traits>
 732     std::basic_ostream<_CharT, _Traits>&
 733     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 734                const normal_distribution<_RealType>& __x)
 735     {
 736       const std::ios_base::fmtflags __flags = __os.flags();
 737       const _CharT __fill = __os.fill();
 738       const std::streamsize __precision = __os.precision();
 739       const _CharT __space = __os.widen(' ');
 740       __os.flags(std::ios_base::scientific | std::ios_base::left);
 741       __os.fill(__space);
 742       __os.precision(_Private::_Max_digits10<_RealType>::__value);
 743
 744       __os << __x.mean() << __space
 745            << __x.sigma() << __space
 746            << __x._M_saved << __space
 747            << __x._M_saved_available;
 748
 749       __os.flags(__flags);
 750       __os.fill(__fill);
 751       __os.precision(__precision);
 752       return __os;
 753     }
 754
 755   template<typename _RealType, typename _CharT, typename _Traits>
 756     std::basic_istream<_CharT, _Traits>&
 757     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 758                normal_distribution<_RealType>& __x)
 759     {
 760       const std::ios_base::fmtflags __flags = __is.flags();
 761       __is.flags(std::ios_base::dec | std::ios_base::skipws);
 762
 763       __is >> __x._M_mean >> __x._M_sigma
 764            >> __x._M_saved >> __x._M_saved_available;
 765
 766       __is.flags(__flags);
 767       return __is;
 768     }
 769
 770 _GLIBCXX_END_NAMESPACE
 771 }