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

   1 // random number generation (out of line) -*- C++ -*-
   2
   3 // Copyright (C) 2009 Free Software Foundation, Inc.
   4 //
   5 // This file is part of the GNU ISO C++ Library.  This library is free
   6 // software; you can redistribute it and/or modify it under the
   7 // terms of the GNU General Public License as published by the
   8 // Free Software Foundation; either version 3, or (at your option)
   9 // any later version.
  10
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 // GNU General Public License for more details.
  15
  16 // Under Section 7 of GPL version 3, you are granted additional
  17 // permissions described in the GCC Runtime Library Exception, version
  18 // 3.1, as published by the Free Software Foundation.
  19
  20 // You should have received a copy of the GNU General Public License and
  21 // a copy of the GCC Runtime Library Exception along with this program;
  22 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  23 // <http://www.gnu.org/licenses/>.
  24
  25 /** @file bits/random.tcc
  26  *  This is an internal header file, included by other library headers.
  27  *  You should not attempt to use it directly.
  28  */
  29
  30 #include <numeric>
  31 #include <algorithm>
  32
  33 namespace std
  34 {
  35   /*
  36    * (Further) implementation-space details.
  37    */
  38   namespace __detail
  39   {
  40     // General case for x = (ax + c) mod m -- use Schrage's algorithm to
  41     // avoid integer overflow.
  42     //
  43     // Because a and c are compile-time integral constants the compiler
  44     // kindly elides any unreachable paths.
  45     //
  46     // Preconditions:  a > 0, m > 0.
  47     //
  48     template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool>
  49       struct _Mod
  50       {
  51         static _Tp
  52         __calc(_Tp __x)
  53         {
  54           if (__a == 1)
  55             __x %= __m;
  56           else
  57             {
  58               static const _Tp __q = __m / __a;
  59               static const _Tp __r = __m % __a;
  60
  61               _Tp __t1 = __a * (__x % __q);
  62               _Tp __t2 = __r * (__x / __q);
  63               if (__t1 >= __t2)
  64                 __x = __t1 - __t2;
  65               else
  66                 __x = __m - __t2 + __t1;
  67             }
  68
  69           if (__c != 0)
  70             {
  71               const _Tp __d = __m - __x;
  72               if (__d > __c)
  73                 __x += __c;
  74               else
  75                 __x = __c - __d;
  76             }
  77           return __x;
  78         }
  79       };
  80
  81     // Special case for m == 0 -- use unsigned integer overflow as modulo
  82     // operator.
  83     template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
  84       struct _Mod<_Tp, __m, __a, __c, true>
  85       {
  86         static _Tp
  87         __calc(_Tp __x)
  88         { return __a * __x + __c; }
  89       };
  90   } // namespace __detail
  91
  92   /**
  93    * Seeds the LCR with integral value @p __s, adjusted so that the
  94    * ring identity is never a member of the convergence set.
  95    */
  96   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  97     void
  98     linear_congruential_engine<_UIntType, __a, __c, __m>::
  99     seed(result_type __s)
 100     {
 101       if ((__detail::__mod<_UIntType, __m>(__c) == 0)
 102           && (__detail::__mod<_UIntType, __m>(__s) == 0))
 103         _M_x = 1;
 104       else
 105         _M_x = __detail::__mod<_UIntType, __m>(__s);
 106     }
 107
 108   /**
 109    * Seeds the LCR engine with a value generated by @p __q.
 110    */
 111   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 112     void
 113     linear_congruential_engine<_UIntType, __a, __c, __m>::
 114     seed(seed_seq& __q)
 115     {
 116       const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
 117                                       : std::__lg(__m);
 118       const _UIntType __k = (__k0 + 31) / 32;
 119       uint_least32_t __arr[__k + 3];
 120       __q.generate(__arr + 0, __arr + __k + 3);
 121       _UIntType __factor = 1u;
 122       _UIntType __sum = 0u;
 123       for (size_t __j = 0; __j < __k; ++__j)
 124         {
 125           __sum += __arr[__j + 3] * __factor;
 126           __factor *= __detail::_Shift<_UIntType, 32>::__value;
 127         }
 128       seed(__sum);
 129     }
 130
 131   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 132            typename _CharT, typename _Traits>
 133     std::basic_ostream<_CharT, _Traits>&
 134     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 135                const linear_congruential_engine<_UIntType,
 136                                                 __a, __c, __m>& __lcr)
 137     {
 138       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 139       typedef typename __ostream_type::ios_base    __ios_base;
 140
 141       const typename __ios_base::fmtflags __flags = __os.flags();
 142       const _CharT __fill = __os.fill();
 143       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 144       __os.fill(__os.widen(' '));
 145
 146       __os << __lcr._M_x;
 147
 148       __os.flags(__flags);
 149       __os.fill(__fill);
 150       return __os;
 151     }
 152
 153   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 154            typename _CharT, typename _Traits>
 155     std::basic_istream<_CharT, _Traits>&
 156     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 157                linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
 158     {
 159       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 160       typedef typename __istream_type::ios_base    __ios_base;
 161
 162       const typename __ios_base::fmtflags __flags = __is.flags();
 163       __is.flags(__ios_base::dec);
 164
 165       __is >> __lcr._M_x;
 166
 167       __is.flags(__flags);
 168       return __is;
 169     }
 170
 171
 172   template<typename _UIntType,
 173            size_t __w, size_t __n, size_t __m, size_t __r,
 174            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 175            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 176            _UIntType __f>
 177     void
 178     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 179                             __s, __b, __t, __c, __l, __f>::
 180     seed(result_type __sd)
 181     {
 182       _M_x[0] = __detail::__mod<_UIntType,
 183         __detail::_Shift<_UIntType, __w>::__value>(__sd);
 184
 185       for (size_t __i = 1; __i < state_size; ++__i)
 186         {
 187           _UIntType __x = _M_x[__i - 1];
 188           __x ^= __x >> (__w - 2);
 189           __x *= __f;
 190           __x += __detail::__mod<_UIntType, __n>(__i);
 191           _M_x[__i] = __detail::__mod<_UIntType,
 192             __detail::_Shift<_UIntType, __w>::__value>(__x);
 193         }
 194       _M_p = state_size;
 195     }
 196
 197   template<typename _UIntType,
 198            size_t __w, size_t __n, size_t __m, size_t __r,
 199            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 200            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 201            _UIntType __f>
 202     void
 203     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 204                               __s, __b, __t, __c, __l, __f>::
 205     seed(seed_seq& __q)
 206     {
 207       const _UIntType __upper_mask = (~_UIntType()) << __r;
 208       const size_t __k = (__w + 31) / 32;
 209       uint_least32_t __arr[__n * __k];
 210       __q.generate(__arr + 0, __arr + __n * __k);
 211
 212       bool __zero = true;
 213       for (size_t __i = 0; __i < state_size; ++__i)
 214         {
 215           _UIntType __factor = 1u;
 216           _UIntType __sum = 0u;
 217           for (size_t __j = 0; __j < __k; ++__j)
 218             {
 219               __sum += __arr[__k * __i + __j] * __factor;
 220               __factor *= __detail::_Shift<_UIntType, 32>::__value;
 221             }
 222           _M_x[__i] = __detail::__mod<_UIntType,
 223             __detail::_Shift<_UIntType, __w>::__value>(__sum);
 224
 225           if (__zero)
 226             {
 227               if (__i == 0)
 228                 {
 229                   if ((_M_x[0] & __upper_mask) != 0u)
 230                     __zero = false;
 231                 }
 232               else if (_M_x[__i] != 0u)
 233                 __zero = false;
 234             }
 235         }
 236         if (__zero)
 237           _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
 238     }
 239
 240   template<typename _UIntType, size_t __w,
 241            size_t __n, size_t __m, size_t __r,
 242            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 243            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 244            _UIntType __f>
 245     typename
 246     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 247                             __s, __b, __t, __c, __l, __f>::result_type
 248     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 249                             __s, __b, __t, __c, __l, __f>::
 250     operator()()
 251     {
 252       // Reload the vector - cost is O(n) amortized over n calls.
 253       if (_M_p >= state_size)
 254         {
 255           const _UIntType __upper_mask = (~_UIntType()) << __r;
 256           const _UIntType __lower_mask = ~__upper_mask;
 257
 258           for (size_t __k = 0; __k < (__n - __m); ++__k)
 259             {
 260               _UIntType __y = ((_M_x[__k] & __upper_mask)
 261                                | (_M_x[__k + 1] & __lower_mask));
 262               _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
 263                            ^ ((__y & 0x01) ? __a : 0));
 264             }
 265
 266           for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
 267             {
 268               _UIntType __y = ((_M_x[__k] & __upper_mask)
 269                                | (_M_x[__k + 1] & __lower_mask));
 270               _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
 271                            ^ ((__y & 0x01) ? __a : 0));
 272             }
 273
 274           _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
 275                            | (_M_x[0] & __lower_mask));
 276           _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
 277                            ^ ((__y & 0x01) ? __a : 0));
 278           _M_p = 0;
 279         }
 280
 281       // Calculate o(x(i)).
 282       result_type __z = _M_x[_M_p++];
 283       __z ^= (__z >> __u) & __d;
 284       __z ^= (__z << __s) & __b;
 285       __z ^= (__z << __t) & __c;
 286       __z ^= (__z >> __l);
 287
 288       return __z;
 289     }
 290
 291   template<typename _UIntType, size_t __w,
 292            size_t __n, size_t __m, size_t __r,
 293            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 294            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 295            _UIntType __f, typename _CharT, typename _Traits>
 296     std::basic_ostream<_CharT, _Traits>&
 297     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 298                const mersenne_twister_engine<_UIntType, __w, __n, __m,
 299                __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
 300     {
 301       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 302       typedef typename __ostream_type::ios_base    __ios_base;
 303
 304       const typename __ios_base::fmtflags __flags = __os.flags();
 305       const _CharT __fill = __os.fill();
 306       const _CharT __space = __os.widen(' ');
 307       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 308       __os.fill(__space);
 309
 310       for (size_t __i = 0; __i < __n - 1; ++__i)
 311         __os << __x._M_x[__i] << __space;
 312       __os << __x._M_x[__n - 1];
 313
 314       __os.flags(__flags);
 315       __os.fill(__fill);
 316       return __os;
 317     }
 318
 319   template<typename _UIntType, size_t __w,
 320            size_t __n, size_t __m, size_t __r,
 321            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 322            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 323            _UIntType __f, typename _CharT, typename _Traits>
 324     std::basic_istream<_CharT, _Traits>&
 325     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 326                mersenne_twister_engine<_UIntType, __w, __n, __m,
 327                __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
 328     {
 329       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 330       typedef typename __istream_type::ios_base    __ios_base;
 331
 332       const typename __ios_base::fmtflags __flags = __is.flags();
 333       __is.flags(__ios_base::dec | __ios_base::skipws);
 334
 335       for (size_t __i = 0; __i < __n; ++__i)
 336         __is >> __x._M_x[__i];
 337
 338       __is.flags(__flags);
 339       return __is;
 340     }
 341
 342
 343   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 344     void
 345     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 346     seed(result_type __value)
 347     {
 348       std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
 349         __lcg(__value == 0u ? default_seed : __value);
 350
 351       const size_t __n = (__w + 31) / 32;
 352
 353       for (size_t __i = 0; __i < long_lag; ++__i)
 354         {
 355           _UIntType __sum = 0u;
 356           _UIntType __factor = 1u;
 357           for (size_t __j = 0; __j < __n; ++__j)
 358             {
 359               __sum += __detail::__mod<uint_least32_t,
 360                        __detail::_Shift<uint_least32_t, 32>::__value>
 361                          (__lcg()) * __factor;
 362               __factor *= __detail::_Shift<_UIntType, 32>::__value;
 363             }
 364           _M_x[__i] = __detail::__mod<_UIntType,
 365             __detail::_Shift<_UIntType, __w>::__value>(__sum);
 366         }
 367       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 368       _M_p = 0;
 369     }
 370
 371   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 372     void
 373     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 374     seed(seed_seq& __q)
 375     {
 376       const size_t __k = (__w + 31) / 32;
 377       uint_least32_t __arr[__r * __k];
 378       __q.generate(__arr + 0, __arr + __r * __k);
 379
 380       for (size_t __i = 0; __i < long_lag; ++__i)
 381         {
 382           _UIntType __sum = 0u;
 383           _UIntType __factor = 1u;
 384           for (size_t __j = 0; __j < __k; ++__j)
 385             {
 386               __sum += __arr[__k * __i + __j] * __factor;
 387               __factor *= __detail::_Shift<_UIntType, 32>::__value;
 388             }
 389           _M_x[__i] = __detail::__mod<_UIntType,
 390             __detail::_Shift<_UIntType, __w>::__value>(__sum);
 391         }
 392       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 393       _M_p = 0;
 394     }
 395
 396   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 397     typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 398              result_type
 399     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 400     operator()()
 401     {
 402       // Derive short lag index from current index.
 403       long __ps = _M_p - short_lag;
 404       if (__ps < 0)
 405         __ps += long_lag;
 406
 407       // Calculate new x(i) without overflow or division.
 408       // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
 409       // cannot overflow.
 410       _UIntType __xi;
 411       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
 412         {
 413           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
 414           _M_carry = 0;
 415         }
 416       else
 417         {
 418           __xi = (__detail::_Shift<_UIntType, __w>::__value
 419                   - _M_x[_M_p] - _M_carry + _M_x[__ps]);
 420           _M_carry = 1;
 421         }
 422       _M_x[_M_p] = __xi;
 423
 424       // Adjust current index to loop around in ring buffer.
 425       if (++_M_p >= long_lag)
 426         _M_p = 0;
 427
 428       return __xi;
 429     }
 430
 431   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
 432            typename _CharT, typename _Traits>
 433     std::basic_ostream<_CharT, _Traits>&
 434     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 435                const subtract_with_carry_engine<_UIntType,
 436                                                 __w, __s, __r>& __x)
 437     {
 438       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 439       typedef typename __ostream_type::ios_base    __ios_base;
 440
 441       const typename __ios_base::fmtflags __flags = __os.flags();
 442       const _CharT __fill = __os.fill();
 443       const _CharT __space = __os.widen(' ');
 444       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 445       __os.fill(__space);
 446
 447       for (size_t __i = 0; __i < __r; ++__i)
 448         __os << __x._M_x[__i] << __space;
 449       __os << __x._M_carry;
 450
 451       __os.flags(__flags);
 452       __os.fill(__fill);
 453       return __os;
 454     }
 455
 456   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
 457            typename _CharT, typename _Traits>
 458     std::basic_istream<_CharT, _Traits>&
 459     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 460                subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
 461     {
 462       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
 463       typedef typename __istream_type::ios_base    __ios_base;
 464
 465       const typename __ios_base::fmtflags __flags = __is.flags();
 466       __is.flags(__ios_base::dec | __ios_base::skipws);
 467
 468       for (size_t __i = 0; __i < __r; ++__i)
 469         __is >> __x._M_x[__i];
 470       __is >> __x._M_carry;
 471
 472       __is.flags(__flags);
 473       return __is;
 474     }
 475
 476
 477   template<typename _RandomNumberEngine, size_t __p, size_t __r>
 478     typename discard_block_engine<_RandomNumberEngine,
 479                            __p, __r>::result_type
 480     discard_block_engine<_RandomNumberEngine, __p, __r>::
 481     operator()()
 482     {
 483       if (_M_n >= used_block)
 484         {
 485           _M_b.discard(block_size - _M_n);
 486           _M_n = 0;
 487         }
 488       ++_M_n;
 489       return _M_b();
 490     }
 491
 492   template<typename _RandomNumberEngine, size_t __p, size_t __r,
 493            typename _CharT, typename _Traits>
 494     std::basic_ostream<_CharT, _Traits>&
 495     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 496                const discard_block_engine<_RandomNumberEngine,
 497                __p, __r>& __x)
 498     {
 499       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 500       typedef typename __ostream_type::ios_base    __ios_base;
 501
 502       const typename __ios_base::fmtflags __flags = __os.flags();
 503       const _CharT __fill = __os.fill();
 504       const _CharT __space = __os.widen(' ');
 505       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 506       __os.fill(__space);
 507
 508       __os << __x.base() << __space << __x._M_n;
 509
 510       __os.flags(__flags);
 511       __os.fill(__fill);
 512       return __os;
 513     }
 514
 515   template<typename _RandomNumberEngine, size_t __p, size_t __r,
 516            typename _CharT, typename _Traits>
 517     std::basic_istream<_CharT, _Traits>&
 518     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 519                discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
 520     {
 521       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 522       typedef typename __istream_type::ios_base    __ios_base;
 523
 524       const typename __ios_base::fmtflags __flags = __is.flags();
 525       __is.flags(__ios_base::dec | __ios_base::skipws);
 526
 527       __is >> __x._M_b >> __x._M_n;
 528
 529       __is.flags(__flags);
 530       return __is;
 531     }
 532
 533
 534   template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
 535     typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
 536       result_type
 537     independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
 538     operator()()
 539     {
 540       const long double __r = static_cast<long double>(_M_b.max())
 541                             - static_cast<long double>(_M_b.min()) + 1.0L;
 542       const result_type __m = std::log(__r) / std::log(2.0L);
 543       result_type __n, __n0, __y0, __y1, __s0, __s1;
 544       for (size_t __i = 0; __i < 2; ++__i)
 545         {
 546           __n = (__w + __m - 1) / __m + __i;
 547           __n0 = __n - __w % __n;
 548           const result_type __w0 = __w / __n;
 549           const result_type __w1 = __w0 + 1;
 550           __s0 = result_type(1) << __w0;
 551           __s1 = result_type(1) << __w1;
 552           __y0 = __s0 * (__r / __s0);
 553           __y1 = __s1 * (__r / __s1);
 554           if (__r - __y0 <= __y0 / __n)
 555             break;
 556         }
 557
 558       result_type __sum = 0;
 559       for (size_t __k = 0; __k < __n0; ++__k)
 560         {
 561           result_type __u;
 562           do
 563             __u = _M_b() - _M_b.min();
 564           while (__u >= __y0);
 565           __sum = __s0 * __sum + __u % __s0;
 566         }
 567       for (size_t __k = __n0; __k < __n; ++__k)
 568         {
 569           result_type __u;
 570           do
 571             __u = _M_b() - _M_b.min();
 572           while (__u >= __y1);
 573           __sum = __s1 * __sum + __u % __s1;
 574         }
 575       return __sum;
 576     }
 577
 578
 579   template<typename _RandomNumberEngine, size_t __k>
 580     typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
 581     shuffle_order_engine<_RandomNumberEngine, __k>::
 582     operator()()
 583     {
 584       size_t __j = __k * ((_M_y - _M_b.min())
 585                           / (_M_b.max() - _M_b.min() + 1.0L));
 586       _M_y = _M_v[__j];
 587       _M_v[__j] = _M_b();
 588
 589       return _M_y;
 590     }
 591
 592   template<typename _RandomNumberEngine, size_t __k,
 593            typename _CharT, typename _Traits>
 594     std::basic_ostream<_CharT, _Traits>&
 595     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 596                const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
 597     {
 598       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 599       typedef typename __ostream_type::ios_base    __ios_base;
 600
 601       const typename __ios_base::fmtflags __flags = __os.flags();
 602       const _CharT __fill = __os.fill();
 603       const _CharT __space = __os.widen(' ');
 604       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 605       __os.fill(__space);
 606
 607       __os << __x.base();
 608       for (size_t __i = 0; __i < __k; ++__i)
 609         __os << __space << __x._M_v[__i];
 610       __os << __space << __x._M_y;
 611
 612       __os.flags(__flags);
 613       __os.fill(__fill);
 614       return __os;
 615     }
 616
 617   template<typename _RandomNumberEngine, size_t __k,
 618            typename _CharT, typename _Traits>
 619     std::basic_istream<_CharT, _Traits>&
 620     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 621                shuffle_order_engine<_RandomNumberEngine, __k>& __x)
 622     {
 623       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 624       typedef typename __istream_type::ios_base    __ios_base;
 625
 626       const typename __ios_base::fmtflags __flags = __is.flags();
 627       __is.flags(__ios_base::dec | __ios_base::skipws);
 628
 629       __is >> __x._M_b;
 630       for (size_t __i = 0; __i < __k; ++__i)
 631         __is >> __x._M_v[__i];
 632       __is >> __x._M_y;
 633
 634       __is.flags(__flags);
 635       return __is;
 636     }
 637
 638
 639   template<typename _IntType>
 640     template<typename _UniformRandomNumberGenerator>
 641       typename uniform_int_distribution<_IntType>::result_type
 642       uniform_int_distribution<_IntType>::
 643       operator()(_UniformRandomNumberGenerator& __urng,
 644                  const param_type& __param)
 645       {
 646         // XXX Must be fixed to work well for *arbitrary* __urng.max(),
 647         // __urng.min(), __param.b(), __param.a().  Currently works fine only
 648         // in the most common case __urng.max() - __urng.min() >=
 649         // __param.b() - __param.a(), with __urng.max() > __urng.min() >= 0.
 650         typedef typename __gnu_cxx::__add_unsigned<typename
 651           _UniformRandomNumberGenerator::result_type>::__type __urntype;
 652         typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
 653                                                               __utype;
 654         typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
 655                                                         > sizeof(__utype)),
 656           __urntype, __utype>::__type                         __uctype;
 657
 658         result_type __ret;
 659
 660         const __urntype __urnmin = __urng.min();
 661         const __urntype __urnmax = __urng.max();
 662         const __urntype __urnrange = __urnmax - __urnmin;
 663         const __uctype __urange = __param.b() - __param.a();
 664         const __uctype __udenom = (__urnrange <= __urange
 665                                    ? 1 : __urnrange / (__urange + 1));
 666         do
 667           __ret = (__urntype(__urng()) -  __urnmin) / __udenom;
 668         while (__ret > __param.b() - __param.a());
 669
 670         return __ret + __param.a();
 671       }
 672
 673   template<typename _IntType, typename _CharT, typename _Traits>
 674     std::basic_ostream<_CharT, _Traits>&
 675     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 676                const uniform_int_distribution<_IntType>& __x)
 677     {
 678       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 679       typedef typename __ostream_type::ios_base    __ios_base;
 680
 681       const typename __ios_base::fmtflags __flags = __os.flags();
 682       const _CharT __fill = __os.fill();
 683       const _CharT __space = __os.widen(' ');
 684       __os.flags(__ios_base::scientific | __ios_base::left);
 685       __os.fill(__space);
 686
 687       __os << __x.a() << __space << __x.b();
 688
 689       __os.flags(__flags);
 690       __os.fill(__fill);
 691       return __os;
 692     }
 693
 694   template<typename _IntType, typename _CharT, typename _Traits>
 695     std::basic_istream<_CharT, _Traits>&
 696     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 697                uniform_int_distribution<_IntType>& __x)
 698     {
 699       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 700       typedef typename __istream_type::ios_base    __ios_base;
 701
 702       const typename __ios_base::fmtflags __flags = __is.flags();
 703       __is.flags(__ios_base::dec | __ios_base::skipws);
 704
 705       _IntType __a, __b;
 706       __is >> __a >> __b;
 707       __x.param(typename uniform_int_distribution<_IntType>::
 708                 param_type(__a, __b));
 709
 710       __is.flags(__flags);
 711       return __is;
 712     }
 713
 714
 715   template<typename _RealType, typename _CharT, typename _Traits>
 716     std::basic_ostream<_CharT, _Traits>&
 717     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 718                const uniform_real_distribution<_RealType>& __x)
 719     {
 720       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 721       typedef typename __ostream_type::ios_base    __ios_base;
 722
 723       const typename __ios_base::fmtflags __flags = __os.flags();
 724       const _CharT __fill = __os.fill();
 725       const std::streamsize __precision = __os.precision();
 726       const _CharT __space = __os.widen(' ');
 727       __os.flags(__ios_base::scientific | __ios_base::left);
 728       __os.fill(__space);
 729       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
 730
 731       __os << __x.a() << __space << __x.b();
 732
 733       __os.flags(__flags);
 734       __os.fill(__fill);
 735       __os.precision(__precision);
 736       return __os;
 737     }
 738
 739   template<typename _RealType, typename _CharT, typename _Traits>
 740     std::basic_istream<_CharT, _Traits>&
 741     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 742                uniform_real_distribution<_RealType>& __x)
 743     {
 744       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 745       typedef typename __istream_type::ios_base    __ios_base;
 746
 747       const typename __ios_base::fmtflags __flags = __is.flags();
 748       __is.flags(__ios_base::skipws);
 749
 750       _RealType __a, __b;
 751       __is >> __a >> __b;
 752       __x.param(typename uniform_real_distribution<_RealType>::
 753                 param_type(__a, __b));
 754
 755       __is.flags(__flags);
 756       return __is;
 757     }
 758
 759
 760   template<typename _CharT, typename _Traits>
 761     std::basic_ostream<_CharT, _Traits>&
 762     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 763                const bernoulli_distribution& __x)
 764     {
 765       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 766       typedef typename __ostream_type::ios_base    __ios_base;
 767
 768       const typename __ios_base::fmtflags __flags = __os.flags();
 769       const _CharT __fill = __os.fill();
 770       const std::streamsize __precision = __os.precision();
 771       __os.flags(__ios_base::scientific | __ios_base::left);
 772       __os.fill(__os.widen(' '));
 773       __os.precision(std::numeric_limits<double>::digits10 + 1);
 774
 775       __os << __x.p();
 776
 777       __os.flags(__flags);
 778       __os.fill(__fill);
 779       __os.precision(__precision);
 780       return __os;
 781     }
 782
 783
 784   template<typename _IntType>
 785     template<typename _UniformRandomNumberGenerator>
 786       typename geometric_distribution<_IntType>::result_type
 787       geometric_distribution<_IntType>::
 788       operator()(_UniformRandomNumberGenerator& __urng,
 789                  const param_type& __param)
 790       {
 791         // About the epsilon thing see this thread:
 792         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
 793         const double __naf =
 794           (1 - std::numeric_limits<double>::epsilon()) / 2;
 795         // The largest _RealType convertible to _IntType.
 796         const double __thr =
 797           std::numeric_limits<_IntType>::max() + __naf;
 798         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
 799           __aurng(__urng);
 800
 801         double __cand;
 802         do
 803           __cand = std::ceil(std::log(__aurng()) / __param._M_log_p);
 804         while (__cand >= __thr);
 805
 806         return result_type(__cand + __naf);
 807       }
 808
 809   template<typename _IntType,
 810            typename _CharT, typename _Traits>
 811     std::basic_ostream<_CharT, _Traits>&
 812     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 813                const geometric_distribution<_IntType>& __x)
 814     {
 815       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 816       typedef typename __ostream_type::ios_base    __ios_base;
 817
 818       const typename __ios_base::fmtflags __flags = __os.flags();
 819       const _CharT __fill = __os.fill();
 820       const std::streamsize __precision = __os.precision();
 821       __os.flags(__ios_base::scientific | __ios_base::left);
 822       __os.fill(__os.widen(' '));
 823       __os.precision(std::numeric_limits<double>::digits10 + 1);
 824
 825       __os << __x.p();
 826
 827       __os.flags(__flags);
 828       __os.fill(__fill);
 829       __os.precision(__precision);
 830       return __os;
 831     }
 832
 833   template<typename _IntType,
 834            typename _CharT, typename _Traits>
 835     std::basic_istream<_CharT, _Traits>&
 836     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 837                geometric_distribution<_IntType>& __x)
 838     {
 839       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 840       typedef typename __istream_type::ios_base    __ios_base;
 841
 842       const typename __ios_base::fmtflags __flags = __is.flags();
 843       __is.flags(__ios_base::skipws);
 844
 845       double __p;
 846       __is >> __p;
 847       __x.param(typename geometric_distribution<_IntType>::param_type(__p));
 848
 849       __is.flags(__flags);
 850       return __is;
 851     }
 852
 853
 854   template<typename _IntType>
 855     template<typename _UniformRandomNumberGenerator>
 856       typename negative_binomial_distribution<_IntType>::result_type
 857       negative_binomial_distribution<_IntType>::
 858       operator()(_UniformRandomNumberGenerator& __urng)
 859       {
 860         const double __y = _M_gd(__urng);
 861
 862         // XXX Is the constructor too slow?
 863         std::poisson_distribution<result_type> __poisson(__y);
 864         return __poisson(__urng);
 865       }
 866
 867   template<typename _IntType>
 868     template<typename _UniformRandomNumberGenerator>
 869       typename negative_binomial_distribution<_IntType>::result_type
 870       negative_binomial_distribution<_IntType>::
 871       operator()(_UniformRandomNumberGenerator& __urng,
 872                  const param_type& __p)
 873       {
 874         typedef typename std::gamma_distribution<result_type>::param_type
 875           param_type;
 876
 877         const double __y =
 878           _M_gd(__urng, param_type(__p.k(), __p.p() / (1.0 - __p.p())));
 879
 880         std::poisson_distribution<result_type> __poisson(__y);
 881         return __poisson(__urng);
 882       }
 883
 884   template<typename _IntType, typename _CharT, typename _Traits>
 885     std::basic_ostream<_CharT, _Traits>&
 886     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 887                const negative_binomial_distribution<_IntType>& __x)
 888     {
 889       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 890       typedef typename __ostream_type::ios_base    __ios_base;
 891
 892       const typename __ios_base::fmtflags __flags = __os.flags();
 893       const _CharT __fill = __os.fill();
 894       const std::streamsize __precision = __os.precision();
 895       const _CharT __space = __os.widen(' ');
 896       __os.flags(__ios_base::scientific | __ios_base::left);
 897       __os.fill(__os.widen(' '));
 898       __os.precision(std::numeric_limits<double>::digits10 + 1);
 899
 900       __os << __x.k() << __space << __x.p()
 901            << __space << __x._M_gd;
 902
 903       __os.flags(__flags);
 904       __os.fill(__fill);
 905       __os.precision(__precision);
 906       return __os;
 907     }
 908
 909   template<typename _IntType, typename _CharT, typename _Traits>
 910     std::basic_istream<_CharT, _Traits>&
 911     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 912                negative_binomial_distribution<_IntType>& __x)
 913     {
 914       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 915       typedef typename __istream_type::ios_base    __ios_base;
 916
 917       const typename __ios_base::fmtflags __flags = __is.flags();
 918       __is.flags(__ios_base::skipws);
 919
 920       _IntType __k;
 921       double __p;
 922       __is >> __k >> __p >> __x._M_gd;
 923       __x.param(typename negative_binomial_distribution<_IntType>::
 924                 param_type(__k, __p));
 925
 926       __is.flags(__flags);
 927       return __is;
 928     }
 929
 930
 931   template<typename _IntType>
 932     void
 933     poisson_distribution<_IntType>::param_type::
 934     _M_initialize()
 935     {
 936 #if _GLIBCXX_USE_C99_MATH_TR1
 937       if (_M_mean >= 12)
 938         {
 939           const double __m = std::floor(_M_mean);
 940           _M_lm_thr = std::log(_M_mean);
 941           _M_lfm = std::lgamma(__m + 1);
 942           _M_sm = std::sqrt(__m);
 943
 944           const double __pi_4 = 0.7853981633974483096156608458198757L;
 945           const double __dx = std::sqrt(2 * __m * std::log(32 * __m
 946                                                               / __pi_4));
 947           _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
 948           const double __cx = 2 * __m + _M_d;
 949           _M_scx = std::sqrt(__cx / 2);
 950           _M_1cx = 1 / __cx;
 951
 952           _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
 953           _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
 954                 / _M_d;
 955         }
 956       else
 957 #endif
 958         _M_lm_thr = std::exp(-_M_mean);
 959       }
 960
 961   /**
 962    * A rejection algorithm when mean >= 12 and a simple method based
 963    * upon the multiplication of uniform random variates otherwise.
 964    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
 965    * is defined.
 966    *
 967    * Reference:
 968    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
 969    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
 970    */
 971   template<typename _IntType>
 972     template<typename _UniformRandomNumberGenerator>
 973       typename poisson_distribution<_IntType>::result_type
 974       poisson_distribution<_IntType>::
 975       operator()(_UniformRandomNumberGenerator& __urng,
 976                  const param_type& __param)
 977       {
 978         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
 979           __aurng(__urng);
 980 #if _GLIBCXX_USE_C99_MATH_TR1
 981         if (__param.mean() >= 12)
 982           {
 983             double __x;
 984
 985             // See comments above...
 986             const double __naf =
 987               (1 - std::numeric_limits<double>::epsilon()) / 2;
 988             const double __thr =
 989               std::numeric_limits<_IntType>::max() + __naf;
 990
 991             const double __m = std::floor(__param.mean());
 992             // sqrt(pi / 2)
 993             const double __spi_2 = 1.2533141373155002512078826424055226L;
 994             const double __c1 = __param._M_sm * __spi_2;
 995             const double __c2 = __param._M_c2b + __c1;
 996             const double __c3 = __c2 + 1;
 997             const double __c4 = __c3 + 1;
 998             // e^(1 / 78)
 999             const double __e178 = 1.0129030479320018583185514777512983L;
1000             const double __c5 = __c4 + __e178;
1001             const double __c = __param._M_cb + __c5;
1002             const double __2cx = 2 * (2 * __m + __param._M_d);
1003
1004             bool __reject = true;
1005             do
1006               {
1007                 const double __u = __c * __aurng();
1008                 const double __e = -std::log(__aurng());
1009
1010                 double __w = 0.0;
1011
1012                 if (__u <= __c1)
1013                   {
1014                     const double __n = _M_nd(__urng);
1015                     const double __y = -std::abs(__n) * __param._M_sm - 1;
1016                     __x = std::floor(__y);
1017                     __w = -__n * __n / 2;
1018                     if (__x < -__m)
1019                       continue;
1020                   }
1021                 else if (__u <= __c2)
1022                   {
1023                     const double __n = _M_nd(__urng);
1024                     const double __y = 1 + std::abs(__n) * __param._M_scx;
1025                     __x = std::ceil(__y);
1026                     __w = __y * (2 - __y) * __param._M_1cx;
1027                     if (__x > __param._M_d)
1028                       continue;
1029                   }
1030                 else if (__u <= __c3)
1031                   // NB: This case not in the book, nor in the Errata,
1032                   // but should be ok...
1033                   __x = -1;
1034                 else if (__u <= __c4)
1035                   __x = 0;
1036                 else if (__u <= __c5)
1037                   __x = 1;
1038                 else
1039                   {
1040                     const double __v = -std::log(__aurng());
1041                     const double __y = __param._M_d
1042                                      + __v * __2cx / __param._M_d;
1043                     __x = std::ceil(__y);
1044                     __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1045                   }
1046
1047                 __reject = (__w - __e - __x * __param._M_lm_thr
1048                             > __param._M_lfm - std::lgamma(__x + __m + 1));
1049
1050                 __reject |= __x + __m >= __thr;
1051
1052               } while (__reject);
1053
1054             return result_type(__x + __m + __naf);
1055           }
1056         else
1057 #endif
1058           {
1059             _IntType     __x = 0;
1060             double __prod = 1.0;
1061
1062             do
1063               {
1064                 __prod *= __aurng();
1065                 __x += 1;
1066               }
1067             while (__prod > __param._M_lm_thr);
1068
1069             return __x - 1;
1070           }
1071       }
1072
1073   template<typename _IntType,
1074            typename _CharT, typename _Traits>
1075     std::basic_ostream<_CharT, _Traits>&
1076     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1077                const poisson_distribution<_IntType>& __x)
1078     {
1079       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1080       typedef typename __ostream_type::ios_base    __ios_base;
1081
1082       const typename __ios_base::fmtflags __flags = __os.flags();
1083       const _CharT __fill = __os.fill();
1084       const std::streamsize __precision = __os.precision();
1085       const _CharT __space = __os.widen(' ');
1086       __os.flags(__ios_base::scientific | __ios_base::left);
1087       __os.fill(__space);
1088       __os.precision(std::numeric_limits<double>::digits10 + 1);
1089
1090       __os << __x.mean() << __space << __x._M_nd;
1091
1092       __os.flags(__flags);
1093       __os.fill(__fill);
1094       __os.precision(__precision);
1095       return __os;
1096     }
1097
1098   template<typename _IntType,
1099            typename _CharT, typename _Traits>
1100     std::basic_istream<_CharT, _Traits>&
1101     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1102                poisson_distribution<_IntType>& __x)
1103     {
1104       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1105       typedef typename __istream_type::ios_base    __ios_base;
1106
1107       const typename __ios_base::fmtflags __flags = __is.flags();
1108       __is.flags(__ios_base::skipws);
1109
1110       double __mean;
1111       __is >> __mean >> __x._M_nd;
1112       __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1113
1114       __is.flags(__flags);
1115       return __is;
1116     }
1117
1118
1119   template<typename _IntType>
1120     void
1121     binomial_distribution<_IntType>::param_type::
1122     _M_initialize()
1123     {
1124       const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1125
1126       _M_easy = true;
1127
1128 #if _GLIBCXX_USE_C99_MATH_TR1
1129       if (_M_t * __p12 >= 8)
1130         {
1131           _M_easy = false;
1132           const double __np = std::floor(_M_t * __p12);
1133           const double __pa = __np / _M_t;
1134           const double __1p = 1 - __pa;
1135
1136           const double __pi_4 = 0.7853981633974483096156608458198757L;
1137           const double __d1x =
1138             std::sqrt(__np * __1p * std::log(32 * __np
1139                                              / (81 * __pi_4 * __1p)));
1140           _M_d1 = std::round(std::max(1.0, __d1x));
1141           const double __d2x =
1142             std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1143                                              / (__pi_4 * __pa)));
1144           _M_d2 = std::round(std::max(1.0, __d2x));
1145
1146           // sqrt(pi / 2)
1147           const double __spi_2 = 1.2533141373155002512078826424055226L;
1148           _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1149           _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1150           _M_c = 2 * _M_d1 / __np;
1151           _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1152           const double __a12 = _M_a1 + _M_s2 * __spi_2;
1153           const double __s1s = _M_s1 * _M_s1;
1154           _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1155                              * 2 * __s1s / _M_d1
1156                              * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1157           const double __s2s = _M_s2 * _M_s2;
1158           _M_s = (_M_a123 + 2 * __s2s / _M_d2
1159                   * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1160           _M_lf = (std::lgamma(__np + 1)
1161                    + std::lgamma(_M_t - __np + 1));
1162           _M_lp1p = std::log(__pa / __1p);
1163
1164           _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1165         }
1166       else
1167 #endif
1168         _M_q = -std::log(1 - __p12);
1169     }
1170
1171   template<typename _IntType>
1172     template<typename _UniformRandomNumberGenerator>
1173       typename binomial_distribution<_IntType>::result_type
1174       binomial_distribution<_IntType>::
1175       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1176       {
1177         _IntType __x = 0;
1178         double __sum = 0.0;
1179         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1180           __aurng(__urng);
1181
1182         do
1183           {
1184             const double __e = -std::log(__aurng());
1185             __sum += __e / (__t - __x);
1186             __x += 1;
1187           }
1188         while (__sum <= _M_param._M_q);
1189
1190         return __x - 1;
1191       }
1192
1193   /**
1194    * A rejection algorithm when t * p >= 8 and a simple waiting time
1195    * method - the second in the referenced book - otherwise.
1196    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1197    * is defined.
1198    *
1199    * Reference:
1200    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1201    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1202    */
1203   template<typename _IntType>
1204     template<typename _UniformRandomNumberGenerator>
1205       typename binomial_distribution<_IntType>::result_type
1206       binomial_distribution<_IntType>::
1207       operator()(_UniformRandomNumberGenerator& __urng,
1208                  const param_type& __param)
1209       {
1210         result_type __ret;
1211         const _IntType __t = __param.t();
1212         const _IntType __p = __param.p();
1213         const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1214         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1215           __aurng(__urng);
1216
1217 #if _GLIBCXX_USE_C99_MATH_TR1
1218         if (!__param._M_easy)
1219           {
1220             double __x;
1221
1222             // See comments above...
1223             const double __naf =
1224               (1 - std::numeric_limits<double>::epsilon()) / 2;
1225             const double __thr =
1226               std::numeric_limits<_IntType>::max() + __naf;
1227
1228             const double __np = std::floor(__t * __p12);
1229
1230             // sqrt(pi / 2)
1231             const double __spi_2 = 1.2533141373155002512078826424055226L;
1232             const double __a1 = __param._M_a1;
1233             const double __a12 = __a1 + __param._M_s2 * __spi_2;
1234             const double __a123 = __param._M_a123;
1235             const double __s1s = __param._M_s1 * __param._M_s1;
1236             const double __s2s = __param._M_s2 * __param._M_s2;
1237
1238             bool __reject;
1239             do
1240               {
1241                 const double __u = __param._M_s * __aurng();
1242
1243                 double __v;
1244
1245                 if (__u <= __a1)
1246                   {
1247                     const double __n = _M_nd(__urng);
1248                     const double __y = __param._M_s1 * std::abs(__n);
1249                     __reject = __y >= __param._M_d1;
1250                     if (!__reject)
1251                       {
1252                         const double __e = -std::log(__aurng());
1253                         __x = std::floor(__y);
1254                         __v = -__e - __n * __n / 2 + __param._M_c;
1255                       }
1256                   }
1257                 else if (__u <= __a12)
1258                   {
1259                     const double __n = _M_nd(__urng);
1260                     const double __y = __param._M_s2 * std::abs(__n);
1261                     __reject = __y >= __param._M_d2;
1262                     if (!__reject)
1263                       {
1264                         const double __e = -std::log(__aurng());
1265                         __x = std::floor(-__y);
1266                         __v = -__e - __n * __n / 2;
1267                       }
1268                   }
1269                 else if (__u <= __a123)
1270                   {
1271                     const double __e1 = -std::log(__aurng());
1272                     const double __e2 = -std::log(__aurng());
1273
1274                     const double __y = __param._M_d1
1275                                      + 2 * __s1s * __e1 / __param._M_d1;
1276                     __x = std::floor(__y);
1277                     __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1278                                                     -__y / (2 * __s1s)));
1279                     __reject = false;
1280                   }
1281                 else
1282                   {
1283                     const double __e1 = -std::log(__aurng());
1284                     const double __e2 = -std::log(__aurng());
1285
1286                     const double __y = __param._M_d2
1287                                      + 2 * __s2s * __e1 / __param._M_d2;
1288                     __x = std::floor(-__y);
1289                     __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1290                     __reject = false;
1291                   }
1292
1293                 __reject = __reject || __x < -__np || __x > __t - __np;
1294                 if (!__reject)
1295                   {
1296                     const double __lfx =
1297                       std::lgamma(__np + __x + 1)
1298                       + std::lgamma(__t - (__np + __x) + 1);
1299                     __reject = __v > __param._M_lf - __lfx
1300                              + __x * __param._M_lp1p;
1301                   }
1302
1303                 __reject |= __x + __np >= __thr;
1304               }
1305             while (__reject);
1306
1307             __x += __np + __naf;
1308
1309             const _IntType __z = _M_waiting(__urng, __t - _IntType(__x));
1310             __ret = _IntType(__x) + __z;
1311           }
1312         else
1313 #endif
1314           __ret = _M_waiting(__urng, __t);
1315
1316         if (__p12 != __p)
1317           __ret = __t - __ret;
1318         return __ret;
1319       }
1320
1321   template<typename _IntType,
1322            typename _CharT, typename _Traits>
1323     std::basic_ostream<_CharT, _Traits>&
1324     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1325                const binomial_distribution<_IntType>& __x)
1326     {
1327       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1328       typedef typename __ostream_type::ios_base    __ios_base;
1329
1330       const typename __ios_base::fmtflags __flags = __os.flags();
1331       const _CharT __fill = __os.fill();
1332       const std::streamsize __precision = __os.precision();
1333       const _CharT __space = __os.widen(' ');
1334       __os.flags(__ios_base::scientific | __ios_base::left);
1335       __os.fill(__space);
1336       __os.precision(std::numeric_limits<double>::digits10 + 1);
1337
1338       __os << __x.t() << __space << __x.p()
1339            << __space << __x._M_nd;
1340
1341       __os.flags(__flags);
1342       __os.fill(__fill);
1343       __os.precision(__precision);
1344       return __os;
1345     }
1346
1347   template<typename _IntType,
1348            typename _CharT, typename _Traits>
1349     std::basic_istream<_CharT, _Traits>&
1350     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1351                binomial_distribution<_IntType>& __x)
1352     {
1353       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1354       typedef typename __istream_type::ios_base    __ios_base;
1355
1356       const typename __ios_base::fmtflags __flags = __is.flags();
1357       __is.flags(__ios_base::dec | __ios_base::skipws);
1358
1359       _IntType __t;
1360       double __p;
1361       __is >> __t >> __p >> __x._M_nd;
1362       __x.param(typename binomial_distribution<_IntType>::
1363                 param_type(__t, __p));
1364
1365       __is.flags(__flags);
1366       return __is;
1367     }
1368
1369
1370   template<typename _RealType, typename _CharT, typename _Traits>
1371     std::basic_ostream<_CharT, _Traits>&
1372     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1373                const exponential_distribution<_RealType>& __x)
1374     {
1375       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1376       typedef typename __ostream_type::ios_base    __ios_base;
1377
1378       const typename __ios_base::fmtflags __flags = __os.flags();
1379       const _CharT __fill = __os.fill();
1380       const std::streamsize __precision = __os.precision();
1381       __os.flags(__ios_base::scientific | __ios_base::left);
1382       __os.fill(__os.widen(' '));
1383       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1384
1385       __os << __x.lambda();
1386
1387       __os.flags(__flags);
1388       __os.fill(__fill);
1389       __os.precision(__precision);
1390       return __os;
1391     }
1392
1393   template<typename _RealType, typename _CharT, typename _Traits>
1394     std::basic_istream<_CharT, _Traits>&
1395     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1396                exponential_distribution<_RealType>& __x)
1397     {
1398       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1399       typedef typename __istream_type::ios_base    __ios_base;
1400
1401       const typename __ios_base::fmtflags __flags = __is.flags();
1402       __is.flags(__ios_base::dec | __ios_base::skipws);
1403
1404       _RealType __lambda;
1405       __is >> __lambda;
1406       __x.param(typename exponential_distribution<_RealType>::
1407                 param_type(__lambda));
1408
1409       __is.flags(__flags);
1410       return __is;
1411     }
1412
1413
1414   /**
1415    * Polar method due to Marsaglia.
1416    *
1417    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
1418    * New York, 1986, Ch. V, Sect. 4.4.
1419    */
1420   template<typename _RealType>
1421     template<typename _UniformRandomNumberGenerator>
1422       typename normal_distribution<_RealType>::result_type
1423       normal_distribution<_RealType>::
1424       operator()(_UniformRandomNumberGenerator& __urng,
1425                  const param_type& __param)
1426       {
1427         result_type __ret;
1428         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1429           __aurng(__urng);
1430
1431         if (_M_saved_available)
1432           {
1433             _M_saved_available = false;
1434             __ret = _M_saved;
1435           }
1436         else
1437           {
1438             result_type __x, __y, __r2;
1439             do
1440               {
1441                 __x = result_type(2.0) * __aurng() - 1.0;
1442                 __y = result_type(2.0) * __aurng() - 1.0;
1443                 __r2 = __x * __x + __y * __y;
1444               }
1445             while (__r2 > 1.0 || __r2 == 0.0);
1446
1447             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1448             _M_saved = __x * __mult;
1449             _M_saved_available = true;
1450             __ret = __y * __mult;
1451           }
1452
1453         __ret = __ret * __param.stddev() + __param.mean();
1454         return __ret;
1455       }
1456
1457   template<typename _RealType, typename _CharT, typename _Traits>
1458     std::basic_ostream<_CharT, _Traits>&
1459     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1460                const normal_distribution<_RealType>& __x)
1461     {
1462       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1463       typedef typename __ostream_type::ios_base    __ios_base;
1464
1465       const typename __ios_base::fmtflags __flags = __os.flags();
1466       const _CharT __fill = __os.fill();
1467       const std::streamsize __precision = __os.precision();
1468       const _CharT __space = __os.widen(' ');
1469       __os.flags(__ios_base::scientific | __ios_base::left);
1470       __os.fill(__space);
1471       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1472
1473       __os << __x.mean() << __space << __x.stddev()
1474            << __space << __x._M_saved_available;
1475       if (__x._M_saved_available)
1476         __os << __space << __x._M_saved;
1477
1478       __os.flags(__flags);
1479       __os.fill(__fill);
1480       __os.precision(__precision);
1481       return __os;
1482     }
1483
1484   template<typename _RealType, typename _CharT, typename _Traits>
1485     std::basic_istream<_CharT, _Traits>&
1486     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1487                normal_distribution<_RealType>& __x)
1488     {
1489       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1490       typedef typename __istream_type::ios_base    __ios_base;
1491
1492       const typename __ios_base::fmtflags __flags = __is.flags();
1493       __is.flags(__ios_base::dec | __ios_base::skipws);
1494
1495       double __mean, __stddev;
1496       __is >> __mean >> __stddev
1497            >> __x._M_saved_available;
1498       if (__x._M_saved_available)
1499         __is >> __x._M_saved;
1500       __x.param(typename normal_distribution<_RealType>::
1501                 param_type(__mean, __stddev));
1502
1503       __is.flags(__flags);
1504       return __is;
1505     }
1506
1507
1508   template<typename _RealType, typename _CharT, typename _Traits>
1509     std::basic_ostream<_CharT, _Traits>&
1510     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1511                const lognormal_distribution<_RealType>& __x)
1512     {
1513       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1514       typedef typename __ostream_type::ios_base    __ios_base;
1515
1516       const typename __ios_base::fmtflags __flags = __os.flags();
1517       const _CharT __fill = __os.fill();
1518       const std::streamsize __precision = __os.precision();
1519       const _CharT __space = __os.widen(' ');
1520       __os.flags(__ios_base::scientific | __ios_base::left);
1521       __os.fill(__space);
1522       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1523
1524       __os << __x.m() << __space << __x.s()
1525            << __space << __x._M_nd;
1526
1527       __os.flags(__flags);
1528       __os.fill(__fill);
1529       __os.precision(__precision);
1530       return __os;
1531     }
1532
1533   template<typename _RealType, typename _CharT, typename _Traits>
1534     std::basic_istream<_CharT, _Traits>&
1535     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1536                lognormal_distribution<_RealType>& __x)
1537     {
1538       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1539       typedef typename __istream_type::ios_base    __ios_base;
1540
1541       const typename __ios_base::fmtflags __flags = __is.flags();
1542       __is.flags(__ios_base::dec | __ios_base::skipws);
1543
1544       _RealType __m, __s;
1545       __is >> __m >> __s >> __x._M_nd;
1546       __x.param(typename lognormal_distribution<_RealType>::
1547                 param_type(__m, __s));
1548
1549       __is.flags(__flags);
1550       return __is;
1551     }
1552
1553
1554   template<typename _RealType, typename _CharT, typename _Traits>
1555     std::basic_ostream<_CharT, _Traits>&
1556     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1557                const chi_squared_distribution<_RealType>& __x)
1558     {
1559       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1560       typedef typename __ostream_type::ios_base    __ios_base;
1561
1562       const typename __ios_base::fmtflags __flags = __os.flags();
1563       const _CharT __fill = __os.fill();
1564       const std::streamsize __precision = __os.precision();
1565       const _CharT __space = __os.widen(' ');
1566       __os.flags(__ios_base::scientific | __ios_base::left);
1567       __os.fill(__space);
1568       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1569
1570       __os << __x.n() << __space << __x._M_gd;
1571
1572       __os.flags(__flags);
1573       __os.fill(__fill);
1574       __os.precision(__precision);
1575       return __os;
1576     }
1577
1578   template<typename _RealType, typename _CharT, typename _Traits>
1579     std::basic_istream<_CharT, _Traits>&
1580     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1581                chi_squared_distribution<_RealType>& __x)
1582     {
1583       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1584       typedef typename __istream_type::ios_base    __ios_base;
1585
1586       const typename __ios_base::fmtflags __flags = __is.flags();
1587       __is.flags(__ios_base::dec | __ios_base::skipws);
1588
1589       _RealType __n;
1590       __is >> __n >> __x._M_gd;
1591       __x.param(typename chi_squared_distribution<_RealType>::
1592                 param_type(__n));
1593
1594       __is.flags(__flags);
1595       return __is;
1596     }
1597
1598
1599   template<typename _RealType>
1600     template<typename _UniformRandomNumberGenerator>
1601       typename cauchy_distribution<_RealType>::result_type
1602       cauchy_distribution<_RealType>::
1603       operator()(_UniformRandomNumberGenerator& __urng,
1604                  const param_type& __p)
1605       {
1606         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1607           __aurng(__urng);
1608         _RealType __u;
1609         do
1610           __u = __aurng();
1611         while (__u == 0.5);
1612
1613         const _RealType __pi = 3.1415926535897932384626433832795029L;
1614         return __p.a() + __p.b() * std::tan(__pi * __u);
1615       }
1616
1617   template<typename _RealType, typename _CharT, typename _Traits>
1618     std::basic_ostream<_CharT, _Traits>&
1619     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1620                const cauchy_distribution<_RealType>& __x)
1621     {
1622       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1623       typedef typename __ostream_type::ios_base    __ios_base;
1624
1625       const typename __ios_base::fmtflags __flags = __os.flags();
1626       const _CharT __fill = __os.fill();
1627       const std::streamsize __precision = __os.precision();
1628       const _CharT __space = __os.widen(' ');
1629       __os.flags(__ios_base::scientific | __ios_base::left);
1630       __os.fill(__space);
1631       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1632
1633       __os << __x.a() << __space << __x.b();
1634
1635       __os.flags(__flags);
1636       __os.fill(__fill);
1637       __os.precision(__precision);
1638       return __os;
1639     }
1640
1641   template<typename _RealType, typename _CharT, typename _Traits>
1642     std::basic_istream<_CharT, _Traits>&
1643     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1644                cauchy_distribution<_RealType>& __x)
1645     {
1646       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1647       typedef typename __istream_type::ios_base    __ios_base;
1648
1649       const typename __ios_base::fmtflags __flags = __is.flags();
1650       __is.flags(__ios_base::dec | __ios_base::skipws);
1651
1652       _RealType __a, __b;
1653       __is >> __a >> __b;
1654       __x.param(typename cauchy_distribution<_RealType>::
1655                 param_type(__a, __b));
1656
1657       __is.flags(__flags);
1658       return __is;
1659     }
1660
1661
1662   template<typename _RealType, typename _CharT, typename _Traits>
1663     std::basic_ostream<_CharT, _Traits>&
1664     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1665                const fisher_f_distribution<_RealType>& __x)
1666     {
1667       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1668       typedef typename __ostream_type::ios_base    __ios_base;
1669
1670       const typename __ios_base::fmtflags __flags = __os.flags();
1671       const _CharT __fill = __os.fill();
1672       const std::streamsize __precision = __os.precision();
1673       const _CharT __space = __os.widen(' ');
1674       __os.flags(__ios_base::scientific | __ios_base::left);
1675       __os.fill(__space);
1676       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1677
1678       __os << __x.m() << __space << __x.n()
1679            << __space << __x._M_gd_x << __space << __x._M_gd_y;
1680
1681       __os.flags(__flags);
1682       __os.fill(__fill);
1683       __os.precision(__precision);
1684       return __os;
1685     }
1686
1687   template<typename _RealType, typename _CharT, typename _Traits>
1688     std::basic_istream<_CharT, _Traits>&
1689     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1690                fisher_f_distribution<_RealType>& __x)
1691     {
1692       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1693       typedef typename __istream_type::ios_base    __ios_base;
1694
1695       const typename __ios_base::fmtflags __flags = __is.flags();
1696       __is.flags(__ios_base::dec | __ios_base::skipws);
1697
1698       _RealType __m, __n;
1699       __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
1700       __x.param(typename fisher_f_distribution<_RealType>::
1701                 param_type(__m, __n));
1702
1703       __is.flags(__flags);
1704       return __is;
1705     }
1706
1707
1708   template<typename _RealType, typename _CharT, typename _Traits>
1709     std::basic_ostream<_CharT, _Traits>&
1710     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1711                const student_t_distribution<_RealType>& __x)
1712     {
1713       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1714       typedef typename __ostream_type::ios_base    __ios_base;
1715
1716       const typename __ios_base::fmtflags __flags = __os.flags();
1717       const _CharT __fill = __os.fill();
1718       const std::streamsize __precision = __os.precision();
1719       const _CharT __space = __os.widen(' ');
1720       __os.flags(__ios_base::scientific | __ios_base::left);
1721       __os.fill(__space);
1722       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1723
1724       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
1725
1726       __os.flags(__flags);
1727       __os.fill(__fill);
1728       __os.precision(__precision);
1729       return __os;
1730     }
1731
1732   template<typename _RealType, typename _CharT, typename _Traits>
1733     std::basic_istream<_CharT, _Traits>&
1734     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1735                student_t_distribution<_RealType>& __x)
1736     {
1737       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1738       typedef typename __istream_type::ios_base    __ios_base;
1739
1740       const typename __ios_base::fmtflags __flags = __is.flags();
1741       __is.flags(__ios_base::dec | __ios_base::skipws);
1742
1743       _RealType __n;
1744       __is >> __n >> __x._M_nd >> __x._M_gd;
1745       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
1746
1747       __is.flags(__flags);
1748       return __is;
1749     }
1750
1751
1752   template<typename _RealType>
1753     void
1754     gamma_distribution<_RealType>::param_type::
1755     _M_initialize()
1756     {
1757       _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
1758
1759       const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
1760       _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
1761     }
1762
1763   /**
1764    * Marsaglia, G. and Tsang, W. W.
1765    * "A Simple Method for Generating Gamma Variables"
1766    * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
1767    */
1768   template<typename _RealType>
1769     template<typename _UniformRandomNumberGenerator>
1770       typename gamma_distribution<_RealType>::result_type
1771       gamma_distribution<_RealType>::
1772       operator()(_UniformRandomNumberGenerator& __urng,
1773                  const param_type& __param)
1774       {
1775         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1776           __aurng(__urng);
1777
1778         result_type __u, __v, __n;
1779         const result_type __a1 = (__param._M_malpha
1780                                   - _RealType(1.0) / _RealType(3.0));
1781
1782         do
1783           {
1784             do
1785               {
1786                 __n = _M_nd(__urng);
1787                 __v = result_type(1.0) + __param._M_a2 * __n;
1788               }
1789             while (__v <= 0.0);
1790
1791             __v = __v * __v * __v;
1792             __u = __aurng();
1793           }
1794         while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
1795                && (std::log(__u) > (0.5 * __n * __n + __a1
1796                                     * (1.0 - __v + std::log(__v)))));
1797
1798         if (__param.alpha() == __param._M_malpha)
1799           return __a1 * __v * __param.beta();
1800         else
1801           {
1802             do
1803               __u = __aurng();
1804             while (__u == 0.0);
1805
1806             return (std::pow(__u, result_type(1.0) / __param.alpha())
1807                     * __a1 * __v * __param.beta());
1808           }
1809       }
1810
1811   template<typename _RealType, typename _CharT, typename _Traits>
1812     std::basic_ostream<_CharT, _Traits>&
1813     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1814                const gamma_distribution<_RealType>& __x)
1815     {
1816       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1817       typedef typename __ostream_type::ios_base    __ios_base;
1818
1819       const typename __ios_base::fmtflags __flags = __os.flags();
1820       const _CharT __fill = __os.fill();
1821       const std::streamsize __precision = __os.precision();
1822       const _CharT __space = __os.widen(' ');
1823       __os.flags(__ios_base::scientific | __ios_base::left);
1824       __os.fill(__space);
1825       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1826
1827       __os << __x.alpha() << __space << __x.beta()
1828            << __space << __x._M_nd;
1829
1830       __os.flags(__flags);
1831       __os.fill(__fill);
1832       __os.precision(__precision);
1833       return __os;
1834     }
1835
1836   template<typename _RealType, typename _CharT, typename _Traits>
1837     std::basic_istream<_CharT, _Traits>&
1838     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1839                gamma_distribution<_RealType>& __x)
1840     {
1841       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1842       typedef typename __istream_type::ios_base    __ios_base;
1843
1844       const typename __ios_base::fmtflags __flags = __is.flags();
1845       __is.flags(__ios_base::dec | __ios_base::skipws);
1846
1847       _RealType __alpha_val, __beta_val;
1848       __is >> __alpha_val >> __beta_val >> __x._M_nd;
1849       __x.param(typename gamma_distribution<_RealType>::
1850                 param_type(__alpha_val, __beta_val));
1851
1852       __is.flags(__flags);
1853       return __is;
1854     }
1855
1856
1857   template<typename _RealType>
1858     template<typename _UniformRandomNumberGenerator>
1859       typename weibull_distribution<_RealType>::result_type
1860       weibull_distribution<_RealType>::
1861       operator()(_UniformRandomNumberGenerator& __urng,
1862                  const param_type& __p)
1863       {
1864         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1865           __aurng(__urng);
1866         return __p.b() * std::pow(-std::log(__aurng()),
1867                                   result_type(1) / __p.a());
1868       }
1869
1870   template<typename _RealType, typename _CharT, typename _Traits>
1871     std::basic_ostream<_CharT, _Traits>&
1872     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1873                const weibull_distribution<_RealType>& __x)
1874     {
1875       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1876       typedef typename __ostream_type::ios_base    __ios_base;
1877
1878       const typename __ios_base::fmtflags __flags = __os.flags();
1879       const _CharT __fill = __os.fill();
1880       const std::streamsize __precision = __os.precision();
1881       const _CharT __space = __os.widen(' ');
1882       __os.flags(__ios_base::scientific | __ios_base::left);
1883       __os.fill(__space);
1884       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1885
1886       __os << __x.a() << __space << __x.b();
1887
1888       __os.flags(__flags);
1889       __os.fill(__fill);
1890       __os.precision(__precision);
1891       return __os;
1892     }
1893
1894   template<typename _RealType, typename _CharT, typename _Traits>
1895     std::basic_istream<_CharT, _Traits>&
1896     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1897                weibull_distribution<_RealType>& __x)
1898     {
1899       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1900       typedef typename __istream_type::ios_base    __ios_base;
1901
1902       const typename __ios_base::fmtflags __flags = __is.flags();
1903       __is.flags(__ios_base::dec | __ios_base::skipws);
1904
1905       _RealType __a, __b;
1906       __is >> __a >> __b;
1907       __x.param(typename weibull_distribution<_RealType>::
1908                 param_type(__a, __b));
1909
1910       __is.flags(__flags);
1911       return __is;
1912     }
1913
1914
1915   template<typename _RealType>
1916     template<typename _UniformRandomNumberGenerator>
1917       typename extreme_value_distribution<_RealType>::result_type
1918       extreme_value_distribution<_RealType>::
1919       operator()(_UniformRandomNumberGenerator& __urng,
1920                  const param_type& __p)
1921       {
1922         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1923           __aurng(__urng);
1924         return __p.a() - __p.b() * std::log(-std::log(__aurng()));
1925       }
1926
1927   template<typename _RealType, typename _CharT, typename _Traits>
1928     std::basic_ostream<_CharT, _Traits>&
1929     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1930                const extreme_value_distribution<_RealType>& __x)
1931     {
1932       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1933       typedef typename __ostream_type::ios_base    __ios_base;
1934
1935       const typename __ios_base::fmtflags __flags = __os.flags();
1936       const _CharT __fill = __os.fill();
1937       const std::streamsize __precision = __os.precision();
1938       const _CharT __space = __os.widen(' ');
1939       __os.flags(__ios_base::scientific | __ios_base::left);
1940       __os.fill(__space);
1941       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
1942
1943       __os << __x.a() << __space << __x.b();
1944
1945       __os.flags(__flags);
1946       __os.fill(__fill);
1947       __os.precision(__precision);
1948       return __os;
1949     }
1950
1951   template<typename _RealType, typename _CharT, typename _Traits>
1952     std::basic_istream<_CharT, _Traits>&
1953     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1954                extreme_value_distribution<_RealType>& __x)
1955     {
1956       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1957       typedef typename __istream_type::ios_base    __ios_base;
1958
1959       const typename __ios_base::fmtflags __flags = __is.flags();
1960       __is.flags(__ios_base::dec | __ios_base::skipws);
1961
1962       _RealType __a, __b;
1963       __is >> __a >> __b;
1964       __x.param(typename extreme_value_distribution<_RealType>::
1965                 param_type(__a, __b));
1966
1967       __is.flags(__flags);
1968       return __is;
1969     }
1970
1971
1972   template<typename _IntType>
1973     void
1974     discrete_distribution<_IntType>::param_type::
1975     _M_initialize()
1976     {
1977       if (_M_prob.size() < 2)
1978         {
1979           _M_prob.clear();
1980           _M_prob.push_back(1.0);
1981           return;
1982         }
1983
1984       const double __sum = std::accumulate(_M_prob.begin(),
1985                                            _M_prob.end(), 0.0);
1986       // Now normalize the probabilites.
1987       std::transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
1988                      std::bind2nd(std::divides<double>(), __sum));
1989       // Accumulate partial sums.
1990       _M_cp.reserve(_M_prob.size());
1991       std::partial_sum(_M_prob.begin(), _M_prob.end(),
1992                        std::back_inserter(_M_cp));
1993       // Make sure the last cumulative probability is one.
1994       _M_cp[_M_cp.size() - 1] = 1.0;
1995     }
1996
1997   template<typename _IntType>
1998     template<typename _Func>
1999       discrete_distribution<_IntType>::param_type::
2000       param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2001       : _M_prob(), _M_cp()
2002       {
2003         const size_t __n = __nw == 0 ? 1 : __nw;
2004         const double __delta = (__xmax - __xmin) / __n;
2005
2006         _M_prob.reserve(__n);
2007         for (size_t __k = 0; __k < __nw; ++__k)
2008           _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2009
2010         _M_initialize();
2011       }
2012
2013   template<typename _IntType>
2014     template<typename _UniformRandomNumberGenerator>
2015       typename discrete_distribution<_IntType>::result_type
2016       discrete_distribution<_IntType>::
2017       operator()(_UniformRandomNumberGenerator& __urng,
2018                  const param_type& __param)
2019       {
2020         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2021           __aurng(__urng);
2022
2023         const double __p = __aurng();
2024         auto __pos = std::lower_bound(__param._M_cp.begin(),
2025                                       __param._M_cp.end(), __p);
2026
2027         return __pos - __param._M_cp.begin();
2028       }
2029
2030   template<typename _IntType, typename _CharT, typename _Traits>
2031     std::basic_ostream<_CharT, _Traits>&
2032     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2033                const discrete_distribution<_IntType>& __x)
2034     {
2035       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2036       typedef typename __ostream_type::ios_base    __ios_base;
2037
2038       const typename __ios_base::fmtflags __flags = __os.flags();
2039       const _CharT __fill = __os.fill();
2040       const std::streamsize __precision = __os.precision();
2041       const _CharT __space = __os.widen(' ');
2042       __os.flags(__ios_base::scientific | __ios_base::left);
2043       __os.fill(__space);
2044       __os.precision(std::numeric_limits<double>::digits10 + 1);
2045
2046       std::vector<double> __prob = __x.probabilities();
2047       __os << __prob.size();
2048       for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2049         __os << __space << *__dit;
2050
2051       __os.flags(__flags);
2052       __os.fill(__fill);
2053       __os.precision(__precision);
2054       return __os;
2055     }
2056
2057   template<typename _IntType, typename _CharT, typename _Traits>
2058     std::basic_istream<_CharT, _Traits>&
2059     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2060                discrete_distribution<_IntType>& __x)
2061     {
2062       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2063       typedef typename __istream_type::ios_base    __ios_base;
2064
2065       const typename __ios_base::fmtflags __flags = __is.flags();
2066       __is.flags(__ios_base::dec | __ios_base::skipws);
2067
2068       size_t __n;
2069       __is >> __n;
2070
2071       std::vector<double> __prob_vec;
2072       __prob_vec.reserve(__n);
2073       for (; __n != 0; --__n)
2074         {
2075           double __prob;
2076           __is >> __prob;
2077           __prob_vec.push_back(__prob);
2078         }
2079
2080       __x.param(typename discrete_distribution<_IntType>::
2081                 param_type(__prob_vec.begin(), __prob_vec.end()));
2082
2083       __is.flags(__flags);
2084       return __is;
2085     }
2086
2087
2088   template<typename _RealType>
2089     void
2090     piecewise_constant_distribution<_RealType>::param_type::
2091     _M_initialize()
2092     {
2093       if (_M_int.size() < 2)
2094         {
2095           _M_int.clear();
2096           _M_int.reserve(2);
2097           _M_int.push_back(_RealType(0));
2098           _M_int.push_back(_RealType(1));
2099
2100           _M_den.clear();
2101           _M_den.push_back(1.0);
2102
2103           return;
2104         }
2105
2106       const double __sum = std::accumulate(_M_den.begin(),
2107                                            _M_den.end(), 0.0);
2108
2109       std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2110                      std::bind2nd(std::divides<double>(), __sum));
2111
2112       _M_cp.reserve(_M_den.size());
2113       std::partial_sum(_M_den.begin(), _M_den.end(),
2114                        std::back_inserter(_M_cp));
2115
2116       // Make sure the last cumulative probability is one.
2117       _M_cp[_M_cp.size() - 1] = 1.0;
2118
2119       for (size_t __k = 0; __k < _M_den.size(); ++__k)
2120         _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2121     }
2122
2123   template<typename _RealType>
2124     template<typename _InputIteratorB, typename _InputIteratorW>
2125       piecewise_constant_distribution<_RealType>::param_type::
2126       param_type(_InputIteratorB __bbegin,
2127                  _InputIteratorB __bend,
2128                  _InputIteratorW __wbegin)
2129       : _M_int(), _M_den(), _M_cp()
2130       {
2131         if (__bbegin != __bend)
2132           {
2133             for (;;)
2134               {
2135                 _M_int.push_back(*__bbegin);
2136                 ++__bbegin;
2137                 if (__bbegin == __bend)
2138                   break;
2139
2140                 _M_den.push_back(*__wbegin);
2141                 ++__wbegin;
2142               }
2143           }
2144
2145         _M_initialize();
2146       }
2147
2148   template<typename _RealType>
2149     template<typename _Func>
2150       piecewise_constant_distribution<_RealType>::param_type::
2151       param_type(initializer_list<_RealType> __bl, _Func __fw)
2152       : _M_int(), _M_den(), _M_cp()
2153       {
2154         _M_int.reserve(__bl.size());
2155         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2156           _M_int.push_back(*__biter);
2157
2158         _M_den.reserve(_M_int.size() - 1);
2159         for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2160           _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2161
2162         _M_initialize();
2163       }
2164
2165   template<typename _RealType>
2166     template<typename _Func>
2167       piecewise_constant_distribution<_RealType>::param_type::
2168       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2169       : _M_int(), _M_den(), _M_cp()
2170       {
2171         const size_t __n = __nw == 0 ? 1 : __nw;
2172         const _RealType __delta = (__xmax - __xmin) / __n;
2173
2174         _M_int.reserve(__n + 1);
2175         for (size_t __k = 0; __k <= __nw; ++__k)
2176           _M_int.push_back(__xmin + __k * __delta);
2177
2178         _M_den.reserve(__n);
2179         for (size_t __k = 0; __k < __nw; ++__k)
2180           _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2181
2182         _M_initialize();
2183       }
2184
2185   template<typename _RealType>
2186     template<typename _UniformRandomNumberGenerator>
2187       typename piecewise_constant_distribution<_RealType>::result_type
2188       piecewise_constant_distribution<_RealType>::
2189       operator()(_UniformRandomNumberGenerator& __urng,
2190                  const param_type& __param)
2191       {
2192         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2193           __aurng(__urng);
2194
2195         const double __p = __aurng();
2196         auto __pos = std::lower_bound(__param._M_cp.begin(),
2197                                       __param._M_cp.end(), __p);
2198         const size_t __i = __pos - __param._M_cp.begin();
2199
2200         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2201
2202         return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2203       }
2204
2205   template<typename _RealType, typename _CharT, typename _Traits>
2206     std::basic_ostream<_CharT, _Traits>&
2207     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2208                const piecewise_constant_distribution<_RealType>& __x)
2209     {
2210       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2211       typedef typename __ostream_type::ios_base    __ios_base;
2212
2213       const typename __ios_base::fmtflags __flags = __os.flags();
2214       const _CharT __fill = __os.fill();
2215       const std::streamsize __precision = __os.precision();
2216       const _CharT __space = __os.widen(' ');
2217       __os.flags(__ios_base::scientific | __ios_base::left);
2218       __os.fill(__space);
2219       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
2220
2221       std::vector<_RealType> __int = __x.intervals();
2222       __os << __int.size() - 1;
2223
2224       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2225         __os << __space << *__xit;
2226
2227       std::vector<double> __den = __x.densities();
2228       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2229         __os << __space << *__dit;
2230
2231       __os.flags(__flags);
2232       __os.fill(__fill);
2233       __os.precision(__precision);
2234       return __os;
2235     }
2236
2237   template<typename _RealType, typename _CharT, typename _Traits>
2238     std::basic_istream<_CharT, _Traits>&
2239     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2240                piecewise_constant_distribution<_RealType>& __x)
2241     {
2242       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2243       typedef typename __istream_type::ios_base    __ios_base;
2244
2245       const typename __ios_base::fmtflags __flags = __is.flags();
2246       __is.flags(__ios_base::dec | __ios_base::skipws);
2247
2248       size_t __n;
2249       __is >> __n;
2250
2251       std::vector<_RealType> __int_vec;
2252       __int_vec.reserve(__n + 1);
2253       for (size_t __i = 0; __i <= __n; ++__i)
2254         {
2255           _RealType __int;
2256           __is >> __int;
2257           __int_vec.push_back(__int);
2258         }
2259
2260       std::vector<double> __den_vec;
2261       __den_vec.reserve(__n);
2262       for (size_t __i = 0; __i < __n; ++__i)
2263         {
2264           double __den;
2265           __is >> __den;
2266           __den_vec.push_back(__den);
2267         }
2268
2269       __x.param(typename piecewise_constant_distribution<_RealType>::
2270           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2271
2272       __is.flags(__flags);
2273       return __is;
2274     }
2275
2276
2277   template<typename _RealType>
2278     void
2279     piecewise_linear_distribution<_RealType>::param_type::
2280     _M_initialize()
2281     {
2282       if (_M_int.size() < 2)
2283         {
2284           _M_int.clear();
2285           _M_int.reserve(2);
2286           _M_int.push_back(_RealType(0));
2287           _M_int.push_back(_RealType(1));
2288
2289           _M_den.clear();
2290           _M_den.reserve(2);
2291           _M_den.push_back(1.0);
2292           _M_den.push_back(1.0);
2293
2294           return;
2295         }
2296
2297       double __sum = 0.0;
2298       _M_cp.reserve(_M_int.size() - 1);
2299       _M_m.reserve(_M_int.size() - 1);
2300       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2301         {
2302           const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
2303           __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
2304           _M_cp.push_back(__sum);
2305           _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
2306         }
2307
2308       //  Now normalize the densities...
2309       std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2310                      std::bind2nd(std::divides<double>(), __sum));
2311       //  ... and partial sums...
2312       std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
2313                      std::bind2nd(std::divides<double>(), __sum));
2314       //  ... and slopes.
2315       std::transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
2316                      std::bind2nd(std::divides<double>(), __sum));
2317       //  Make sure the last cumulative probablility is one.
2318       _M_cp[_M_cp.size() - 1] = 1.0;
2319      }
2320
2321   template<typename _RealType>
2322     template<typename _InputIteratorB, typename _InputIteratorW>
2323       piecewise_linear_distribution<_RealType>::param_type::
2324       param_type(_InputIteratorB __bbegin,
2325                  _InputIteratorB __bend,
2326                  _InputIteratorW __wbegin)
2327       : _M_int(), _M_den(), _M_cp(), _M_m()
2328       {
2329         for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
2330           {
2331             _M_int.push_back(*__bbegin);
2332             _M_den.push_back(*__wbegin);
2333           }
2334
2335         _M_initialize();
2336       }
2337
2338   template<typename _RealType>
2339     template<typename _Func>
2340       piecewise_linear_distribution<_RealType>::param_type::
2341       param_type(initializer_list<_RealType> __bl, _Func __fw)
2342       : _M_int(), _M_den(), _M_cp(), _M_m()
2343       {
2344         _M_int.reserve(__bl.size());
2345         _M_den.reserve(__bl.size());
2346         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2347           {
2348             _M_int.push_back(*__biter);
2349             _M_den.push_back(__fw(*__biter));
2350           }
2351
2352         _M_initialize();
2353       }
2354
2355   template<typename _RealType>
2356     template<typename _Func>
2357       piecewise_linear_distribution<_RealType>::param_type::
2358       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2359       : _M_int(), _M_den(), _M_cp(), _M_m()
2360       {
2361         const size_t __n = __nw == 0 ? 1 : __nw;
2362         const _RealType __delta = (__xmax - __xmin) / __n;
2363
2364         _M_int.reserve(__n + 1);
2365         _M_den.reserve(__n + 1);
2366         for (size_t __k = 0; __k <= __nw; ++__k)
2367           {
2368             _M_int.push_back(__xmin + __k * __delta);
2369             _M_den.push_back(__fw(_M_int[__k] + __delta));
2370           }
2371
2372         _M_initialize();
2373       }
2374
2375   template<typename _RealType>
2376     template<typename _UniformRandomNumberGenerator>
2377       typename piecewise_linear_distribution<_RealType>::result_type
2378       piecewise_linear_distribution<_RealType>::
2379       operator()(_UniformRandomNumberGenerator& __urng,
2380                  const param_type& __param)
2381       {
2382         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2383           __aurng(__urng);
2384
2385         const double __p = __aurng();
2386         auto __pos = std::lower_bound(__param._M_cp.begin(),
2387                                       __param._M_cp.end(), __p);
2388         const size_t __i = __pos - __param._M_cp.begin();
2389
2390         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2391
2392         const double __a = 0.5 * __param._M_m[__i];
2393         const double __b = __param._M_den[__i];
2394         const double __cm = __p - __pref;
2395
2396         _RealType __x = __param._M_int[__i];
2397         if (__a == 0)
2398           __x += __cm / __b;
2399         else
2400           {
2401             const double __d = __b * __b + 4.0 * __a * __cm;
2402             __x += 0.5 * (std::sqrt(__d) - __b) / __a;
2403           }
2404
2405         return __x;
2406       }
2407
2408   template<typename _RealType, typename _CharT, typename _Traits>
2409     std::basic_ostream<_CharT, _Traits>&
2410     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2411                const piecewise_linear_distribution<_RealType>& __x)
2412     {
2413       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2414       typedef typename __ostream_type::ios_base    __ios_base;
2415
2416       const typename __ios_base::fmtflags __flags = __os.flags();
2417       const _CharT __fill = __os.fill();
2418       const std::streamsize __precision = __os.precision();
2419       const _CharT __space = __os.widen(' ');
2420       __os.flags(__ios_base::scientific | __ios_base::left);
2421       __os.fill(__space);
2422       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
2423
2424       std::vector<_RealType> __int = __x.intervals();
2425       __os << __int.size() - 1;
2426
2427       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2428         __os << __space << *__xit;
2429
2430       std::vector<double> __den = __x.densities();
2431       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2432         __os << __space << *__dit;
2433
2434       __os.flags(__flags);
2435       __os.fill(__fill);
2436       __os.precision(__precision);
2437       return __os;
2438     }
2439
2440   template<typename _RealType, typename _CharT, typename _Traits>
2441     std::basic_istream<_CharT, _Traits>&
2442     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2443                piecewise_linear_distribution<_RealType>& __x)
2444     {
2445       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2446       typedef typename __istream_type::ios_base    __ios_base;
2447
2448       const typename __ios_base::fmtflags __flags = __is.flags();
2449       __is.flags(__ios_base::dec | __ios_base::skipws);
2450
2451       size_t __n;
2452       __is >> __n;
2453
2454       std::vector<_RealType> __int_vec;
2455       __int_vec.reserve(__n + 1);
2456       for (size_t __i = 0; __i <= __n; ++__i)
2457         {
2458           _RealType __int;
2459           __is >> __int;
2460           __int_vec.push_back(__int);
2461         }
2462
2463       std::vector<double> __den_vec;
2464       __den_vec.reserve(__n + 1);
2465       for (size_t __i = 0; __i <= __n; ++__i)
2466         {
2467           double __den;
2468           __is >> __den;
2469           __den_vec.push_back(__den);
2470         }
2471
2472       __x.param(typename piecewise_linear_distribution<_RealType>::
2473           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2474
2475       __is.flags(__flags);
2476       return __is;
2477     }
2478
2479
2480   template<typename _IntType>
2481     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
2482     {
2483       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
2484         _M_v.push_back(__detail::__mod<result_type,
2485                        __detail::_Shift<result_type, 32>::__value>(*__iter));
2486     }
2487
2488   template<typename _InputIterator>
2489     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
2490     {
2491       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
2492         _M_v.push_back(__detail::__mod<result_type,
2493                        __detail::_Shift<result_type, 32>::__value>(*__iter));
2494     }
2495
2496   template<typename _RandomAccessIterator>
2497     void
2498     seed_seq::generate(_RandomAccessIterator __begin,
2499                        _RandomAccessIterator __end)
2500     {
2501       typedef typename iterator_traits<_RandomAccessIterator>::value_type
2502         _Type;
2503
2504       if (__begin == __end)
2505         return;
2506
2507       std::fill(__begin, __end, _Type(0x8b8b8b8bu));
2508
2509       const size_t __n = __end - __begin;
2510       const size_t __s = _M_v.size();
2511       const size_t __t = (__n >= 623) ? 11
2512                        : (__n >=  68) ? 7
2513                        : (__n >=  39) ? 5
2514                        : (__n >=   7) ? 3
2515                        : (__n - 1) / 2;
2516       const size_t __p = (__n - __t) / 2;
2517       const size_t __q = __p + __t;
2518       const size_t __m = std::max(__s + 1, __n);
2519
2520       for (size_t __k = 0; __k < __m; ++__k)
2521         {
2522           _Type __arg = (__begin[__k % __n]
2523                          ^ __begin[(__k + __p) % __n]
2524                          ^ __begin[(__k - 1) % __n]);
2525           _Type __r1 = __arg ^ (__arg << 27);
2526           __r1 = __detail::__mod<_Type, __detail::_Shift<_Type, 32>::__value,
2527                                  1664525u, 0u>(__r1);
2528           _Type __r2 = __r1;
2529           if (__k == 0)
2530             __r2 += __s;
2531           else if (__k <= __s)
2532             __r2 += __k % __n + _M_v[__k - 1];
2533           else
2534             __r2 += __k % __n;
2535           __r2 = __detail::__mod<_Type,
2536                    __detail::_Shift<_Type, 32>::__value>(__r2);
2537           __begin[(__k + __p) % __n] += __r1;
2538           __begin[(__k + __q) % __n] += __r2;
2539           __begin[__k % __n] = __r2;
2540         }
2541
2542       for (size_t __k = __m; __k < __m + __n; ++__k)
2543         {
2544           _Type __arg = (__begin[__k % __n]
2545                          + __begin[(__k + __p) % __n]
2546                          + __begin[(__k - 1) % __n]);
2547           _Type __r3 = __arg ^ (__arg << 27);
2548           __r3 = __detail::__mod<_Type, __detail::_Shift<_Type, 32>::__value,
2549                                  1566083941u, 0u>(__r3);
2550           _Type __r4 = __r3 - __k % __n;
2551           __r4 = __detail::__mod<_Type,
2552                    __detail::_Shift<_Type, 32>::__value>(__r4);
2553           __begin[(__k + __p) % __n] ^= __r4;
2554           __begin[(__k + __q) % __n] ^= __r3;
2555           __begin[__k % __n] = __r4;
2556         }
2557     }
2558
2559   template<typename _RealType, size_t __bits,
2560            typename _UniformRandomNumberGenerator>
2561     _RealType
2562     generate_canonical(_UniformRandomNumberGenerator& __urng)
2563     {
2564       const size_t __b
2565         = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
2566                    __bits);
2567       const long double __r = static_cast<long double>(__urng.max())
2568                             - static_cast<long double>(__urng.min()) + 1.0L;
2569       const size_t __log2r = std::log(__r) / std::log(2.0L);
2570       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
2571       _RealType __sum = _RealType(0);
2572       _RealType __tmp = _RealType(1);
2573       for (; __k != 0; --__k)
2574         {
2575           __sum += _RealType(__urng() - __urng.min()) * __tmp;
2576           __tmp *= __r;
2577         }
2578       return __sum / __tmp;
2579     }
2580 }