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