contrib/gcc-5.0/libstdc++-v3/include/bits/random.tcc

   1 // random number generation (out of line) -*- C++ -*-
   2
   3 // Copyright (C) 2009-2015 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  *  Do not attempt to use it directly. @headername{random}
  28  */
  29
  30 #ifndef _RANDOM_TCC
  31 #define _RANDOM_TCC 1
  32
  33 #include <numeric> // std::accumulate and std::partial_sum
  34
  35 namespace std _GLIBCXX_VISIBILITY(default)
  36 {
  37   /*
  38    * (Further) implementation-space details.
  39    */
  40   namespace __detail
  41   {
  42   _GLIBCXX_BEGIN_NAMESPACE_VERSION
  43
  44     // General case for x = (ax + c) mod m -- use Schrage's algorithm
  45     // to avoid integer overflow.
  46     //
  47     // Preconditions:  a > 0, m > 0.
  48     //
  49     // Note: only works correctly for __m % __a < __m / __a.
  50     template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
  51       _Tp
  52       _Mod<_Tp, __m, __a, __c, false, true>::
  53       __calc(_Tp __x)
  54       {
  55         if (__a == 1)
  56           __x %= __m;
  57         else
  58           {
  59             static const _Tp __q = __m / __a;
  60             static const _Tp __r = __m % __a;
  61
  62             _Tp __t1 = __a * (__x % __q);
  63             _Tp __t2 = __r * (__x / __q);
  64             if (__t1 >= __t2)
  65               __x = __t1 - __t2;
  66             else
  67               __x = __m - __t2 + __t1;
  68           }
  69
  70         if (__c != 0)
  71           {
  72             const _Tp __d = __m - __x;
  73             if (__d > __c)
  74               __x += __c;
  75             else
  76               __x = __c - __d;
  77           }
  78         return __x;
  79       }
  80
  81     template<typename _InputIterator, typename _OutputIterator,
  82              typename _Tp>
  83       _OutputIterator
  84       __normalize(_InputIterator __first, _InputIterator __last,
  85                   _OutputIterator __result, const _Tp& __factor)
  86       {
  87         for (; __first != __last; ++__first, ++__result)
  88           *__result = *__first / __factor;
  89         return __result;
  90       }
  91
  92   _GLIBCXX_END_NAMESPACE_VERSION
  93   } // namespace __detail
  94
  95 _GLIBCXX_BEGIN_NAMESPACE_VERSION
  96
  97   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
  98     constexpr _UIntType
  99     linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
 100
 101   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 102     constexpr _UIntType
 103     linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
 104
 105   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 106     constexpr _UIntType
 107     linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
 108
 109   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 110     constexpr _UIntType
 111     linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
 112
 113   /**
 114    * Seeds the LCR with integral value @p __s, adjusted so that the
 115    * ring identity is never a member of the convergence set.
 116    */
 117   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 118     void
 119     linear_congruential_engine<_UIntType, __a, __c, __m>::
 120     seed(result_type __s)
 121     {
 122       if ((__detail::__mod<_UIntType, __m>(__c) == 0)
 123           && (__detail::__mod<_UIntType, __m>(__s) == 0))
 124         _M_x = 1;
 125       else
 126         _M_x = __detail::__mod<_UIntType, __m>(__s);
 127     }
 128
 129   /**
 130    * Seeds the LCR engine with a value generated by @p __q.
 131    */
 132   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 133     template<typename _Sseq>
 134       typename std::enable_if<std::is_class<_Sseq>::value>::type
 135       linear_congruential_engine<_UIntType, __a, __c, __m>::
 136       seed(_Sseq& __q)
 137       {
 138         const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
 139                                         : std::__lg(__m);
 140         const _UIntType __k = (__k0 + 31) / 32;
 141         uint_least32_t __arr[__k + 3];
 142         __q.generate(__arr + 0, __arr + __k + 3);
 143         _UIntType __factor = 1u;
 144         _UIntType __sum = 0u;
 145         for (size_t __j = 0; __j < __k; ++__j)
 146           {
 147             __sum += __arr[__j + 3] * __factor;
 148             __factor *= __detail::_Shift<_UIntType, 32>::__value;
 149           }
 150         seed(__sum);
 151       }
 152
 153   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 154            typename _CharT, typename _Traits>
 155     std::basic_ostream<_CharT, _Traits>&
 156     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 157                const linear_congruential_engine<_UIntType,
 158                                                 __a, __c, __m>& __lcr)
 159     {
 160       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 161       typedef typename __ostream_type::ios_base    __ios_base;
 162
 163       const typename __ios_base::fmtflags __flags = __os.flags();
 164       const _CharT __fill = __os.fill();
 165       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 166       __os.fill(__os.widen(' '));
 167
 168       __os << __lcr._M_x;
 169
 170       __os.flags(__flags);
 171       __os.fill(__fill);
 172       return __os;
 173     }
 174
 175   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
 176            typename _CharT, typename _Traits>
 177     std::basic_istream<_CharT, _Traits>&
 178     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 179                linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
 180     {
 181       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 182       typedef typename __istream_type::ios_base    __ios_base;
 183
 184       const typename __ios_base::fmtflags __flags = __is.flags();
 185       __is.flags(__ios_base::dec);
 186
 187       __is >> __lcr._M_x;
 188
 189       __is.flags(__flags);
 190       return __is;
 191     }
 192
 193
 194   template<typename _UIntType,
 195            size_t __w, size_t __n, size_t __m, size_t __r,
 196            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 197            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 198            _UIntType __f>
 199     constexpr size_t
 200     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 201                             __s, __b, __t, __c, __l, __f>::word_size;
 202
 203   template<typename _UIntType,
 204            size_t __w, size_t __n, size_t __m, size_t __r,
 205            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 206            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 207            _UIntType __f>
 208     constexpr size_t
 209     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 210                             __s, __b, __t, __c, __l, __f>::state_size;
 211
 212   template<typename _UIntType,
 213            size_t __w, size_t __n, size_t __m, size_t __r,
 214            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 215            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 216            _UIntType __f>
 217     constexpr size_t
 218     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 219                             __s, __b, __t, __c, __l, __f>::shift_size;
 220
 221   template<typename _UIntType,
 222            size_t __w, size_t __n, size_t __m, size_t __r,
 223            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 224            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 225            _UIntType __f>
 226     constexpr size_t
 227     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 228                             __s, __b, __t, __c, __l, __f>::mask_bits;
 229
 230   template<typename _UIntType,
 231            size_t __w, size_t __n, size_t __m, size_t __r,
 232            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 233            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 234            _UIntType __f>
 235     constexpr _UIntType
 236     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 237                             __s, __b, __t, __c, __l, __f>::xor_mask;
 238
 239   template<typename _UIntType,
 240            size_t __w, size_t __n, size_t __m, size_t __r,
 241            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 242            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 243            _UIntType __f>
 244     constexpr size_t
 245     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 246                             __s, __b, __t, __c, __l, __f>::tempering_u;
 247
 248   template<typename _UIntType,
 249            size_t __w, size_t __n, size_t __m, size_t __r,
 250            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 251            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 252            _UIntType __f>
 253     constexpr _UIntType
 254     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 255                             __s, __b, __t, __c, __l, __f>::tempering_d;
 256
 257   template<typename _UIntType,
 258            size_t __w, size_t __n, size_t __m, size_t __r,
 259            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 260            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 261            _UIntType __f>
 262     constexpr size_t
 263     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 264                             __s, __b, __t, __c, __l, __f>::tempering_s;
 265
 266   template<typename _UIntType,
 267            size_t __w, size_t __n, size_t __m, size_t __r,
 268            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 269            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 270            _UIntType __f>
 271     constexpr _UIntType
 272     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 273                             __s, __b, __t, __c, __l, __f>::tempering_b;
 274
 275   template<typename _UIntType,
 276            size_t __w, size_t __n, size_t __m, size_t __r,
 277            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 278            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 279            _UIntType __f>
 280     constexpr size_t
 281     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 282                             __s, __b, __t, __c, __l, __f>::tempering_t;
 283
 284   template<typename _UIntType,
 285            size_t __w, size_t __n, size_t __m, size_t __r,
 286            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 287            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 288            _UIntType __f>
 289     constexpr _UIntType
 290     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 291                             __s, __b, __t, __c, __l, __f>::tempering_c;
 292
 293   template<typename _UIntType,
 294            size_t __w, size_t __n, size_t __m, size_t __r,
 295            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 296            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 297            _UIntType __f>
 298     constexpr size_t
 299     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 300                             __s, __b, __t, __c, __l, __f>::tempering_l;
 301
 302   template<typename _UIntType,
 303            size_t __w, size_t __n, size_t __m, size_t __r,
 304            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 305            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 306            _UIntType __f>
 307     constexpr _UIntType
 308     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 309                             __s, __b, __t, __c, __l, __f>::
 310                                               initialization_multiplier;
 311
 312   template<typename _UIntType,
 313            size_t __w, size_t __n, size_t __m, size_t __r,
 314            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 315            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 316            _UIntType __f>
 317     constexpr _UIntType
 318     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 319                             __s, __b, __t, __c, __l, __f>::default_seed;
 320
 321   template<typename _UIntType,
 322            size_t __w, size_t __n, size_t __m, size_t __r,
 323            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 324            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 325            _UIntType __f>
 326     void
 327     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 328                             __s, __b, __t, __c, __l, __f>::
 329     seed(result_type __sd)
 330     {
 331       _M_x[0] = __detail::__mod<_UIntType,
 332         __detail::_Shift<_UIntType, __w>::__value>(__sd);
 333
 334       for (size_t __i = 1; __i < state_size; ++__i)
 335         {
 336           _UIntType __x = _M_x[__i - 1];
 337           __x ^= __x >> (__w - 2);
 338           __x *= __f;
 339           __x += __detail::__mod<_UIntType, __n>(__i);
 340           _M_x[__i] = __detail::__mod<_UIntType,
 341             __detail::_Shift<_UIntType, __w>::__value>(__x);
 342         }
 343       _M_p = state_size;
 344     }
 345
 346   template<typename _UIntType,
 347            size_t __w, size_t __n, size_t __m, size_t __r,
 348            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 349            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 350            _UIntType __f>
 351     template<typename _Sseq>
 352       typename std::enable_if<std::is_class<_Sseq>::value>::type
 353       mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 354                               __s, __b, __t, __c, __l, __f>::
 355       seed(_Sseq& __q)
 356       {
 357         const _UIntType __upper_mask = (~_UIntType()) << __r;
 358         const size_t __k = (__w + 31) / 32;
 359         uint_least32_t __arr[__n * __k];
 360         __q.generate(__arr + 0, __arr + __n * __k);
 361
 362         bool __zero = true;
 363         for (size_t __i = 0; __i < state_size; ++__i)
 364           {
 365             _UIntType __factor = 1u;
 366             _UIntType __sum = 0u;
 367             for (size_t __j = 0; __j < __k; ++__j)
 368               {
 369                 __sum += __arr[__k * __i + __j] * __factor;
 370                 __factor *= __detail::_Shift<_UIntType, 32>::__value;
 371               }
 372             _M_x[__i] = __detail::__mod<_UIntType,
 373               __detail::_Shift<_UIntType, __w>::__value>(__sum);
 374
 375             if (__zero)
 376               {
 377                 if (__i == 0)
 378                   {
 379                     if ((_M_x[0] & __upper_mask) != 0u)
 380                       __zero = false;
 381                   }
 382                 else if (_M_x[__i] != 0u)
 383                   __zero = false;
 384               }
 385           }
 386         if (__zero)
 387           _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
 388         _M_p = state_size;
 389       }
 390
 391   template<typename _UIntType, size_t __w,
 392            size_t __n, size_t __m, size_t __r,
 393            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 394            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 395            _UIntType __f>
 396     void
 397     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 398                             __s, __b, __t, __c, __l, __f>::
 399     _M_gen_rand(void)
 400     {
 401       const _UIntType __upper_mask = (~_UIntType()) << __r;
 402       const _UIntType __lower_mask = ~__upper_mask;
 403
 404       for (size_t __k = 0; __k < (__n - __m); ++__k)
 405         {
 406           _UIntType __y = ((_M_x[__k] & __upper_mask)
 407                            | (_M_x[__k + 1] & __lower_mask));
 408           _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
 409                        ^ ((__y & 0x01) ? __a : 0));
 410         }
 411
 412       for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
 413         {
 414           _UIntType __y = ((_M_x[__k] & __upper_mask)
 415                            | (_M_x[__k + 1] & __lower_mask));
 416           _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
 417                        ^ ((__y & 0x01) ? __a : 0));
 418         }
 419
 420       _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
 421                        | (_M_x[0] & __lower_mask));
 422       _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
 423                        ^ ((__y & 0x01) ? __a : 0));
 424       _M_p = 0;
 425     }
 426
 427   template<typename _UIntType, size_t __w,
 428            size_t __n, size_t __m, size_t __r,
 429            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 430            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 431            _UIntType __f>
 432     void
 433     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 434                             __s, __b, __t, __c, __l, __f>::
 435     discard(unsigned long long __z)
 436     {
 437       while (__z > state_size - _M_p)
 438         {
 439           __z -= state_size - _M_p;
 440           _M_gen_rand();
 441         }
 442       _M_p += __z;
 443     }
 444
 445   template<typename _UIntType, size_t __w,
 446            size_t __n, size_t __m, size_t __r,
 447            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 448            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 449            _UIntType __f>
 450     typename
 451     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 452                             __s, __b, __t, __c, __l, __f>::result_type
 453     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
 454                             __s, __b, __t, __c, __l, __f>::
 455     operator()()
 456     {
 457       // Reload the vector - cost is O(n) amortized over n calls.
 458       if (_M_p >= state_size)
 459         _M_gen_rand();
 460
 461       // Calculate o(x(i)).
 462       result_type __z = _M_x[_M_p++];
 463       __z ^= (__z >> __u) & __d;
 464       __z ^= (__z << __s) & __b;
 465       __z ^= (__z << __t) & __c;
 466       __z ^= (__z >> __l);
 467
 468       return __z;
 469     }
 470
 471   template<typename _UIntType, size_t __w,
 472            size_t __n, size_t __m, size_t __r,
 473            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 474            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 475            _UIntType __f, typename _CharT, typename _Traits>
 476     std::basic_ostream<_CharT, _Traits>&
 477     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 478                const mersenne_twister_engine<_UIntType, __w, __n, __m,
 479                __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
 480     {
 481       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 482       typedef typename __ostream_type::ios_base    __ios_base;
 483
 484       const typename __ios_base::fmtflags __flags = __os.flags();
 485       const _CharT __fill = __os.fill();
 486       const _CharT __space = __os.widen(' ');
 487       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 488       __os.fill(__space);
 489
 490       for (size_t __i = 0; __i < __n; ++__i)
 491         __os << __x._M_x[__i] << __space;
 492       __os << __x._M_p;
 493
 494       __os.flags(__flags);
 495       __os.fill(__fill);
 496       return __os;
 497     }
 498
 499   template<typename _UIntType, size_t __w,
 500            size_t __n, size_t __m, size_t __r,
 501            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
 502            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
 503            _UIntType __f, typename _CharT, typename _Traits>
 504     std::basic_istream<_CharT, _Traits>&
 505     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 506                mersenne_twister_engine<_UIntType, __w, __n, __m,
 507                __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
 508     {
 509       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 510       typedef typename __istream_type::ios_base    __ios_base;
 511
 512       const typename __ios_base::fmtflags __flags = __is.flags();
 513       __is.flags(__ios_base::dec | __ios_base::skipws);
 514
 515       for (size_t __i = 0; __i < __n; ++__i)
 516         __is >> __x._M_x[__i];
 517       __is >> __x._M_p;
 518
 519       __is.flags(__flags);
 520       return __is;
 521     }
 522
 523
 524   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 525     constexpr size_t
 526     subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
 527
 528   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 529     constexpr size_t
 530     subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
 531
 532   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 533     constexpr size_t
 534     subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
 535
 536   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 537     constexpr _UIntType
 538     subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
 539
 540   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 541     void
 542     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 543     seed(result_type __value)
 544     {
 545       std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
 546         __lcg(__value == 0u ? default_seed : __value);
 547
 548       const size_t __n = (__w + 31) / 32;
 549
 550       for (size_t __i = 0; __i < long_lag; ++__i)
 551         {
 552           _UIntType __sum = 0u;
 553           _UIntType __factor = 1u;
 554           for (size_t __j = 0; __j < __n; ++__j)
 555             {
 556               __sum += __detail::__mod<uint_least32_t,
 557                        __detail::_Shift<uint_least32_t, 32>::__value>
 558                          (__lcg()) * __factor;
 559               __factor *= __detail::_Shift<_UIntType, 32>::__value;
 560             }
 561           _M_x[__i] = __detail::__mod<_UIntType,
 562             __detail::_Shift<_UIntType, __w>::__value>(__sum);
 563         }
 564       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 565       _M_p = 0;
 566     }
 567
 568   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 569     template<typename _Sseq>
 570       typename std::enable_if<std::is_class<_Sseq>::value>::type
 571       subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 572       seed(_Sseq& __q)
 573       {
 574         const size_t __k = (__w + 31) / 32;
 575         uint_least32_t __arr[__r * __k];
 576         __q.generate(__arr + 0, __arr + __r * __k);
 577
 578         for (size_t __i = 0; __i < long_lag; ++__i)
 579           {
 580             _UIntType __sum = 0u;
 581             _UIntType __factor = 1u;
 582             for (size_t __j = 0; __j < __k; ++__j)
 583               {
 584                 __sum += __arr[__k * __i + __j] * __factor;
 585                 __factor *= __detail::_Shift<_UIntType, 32>::__value;
 586               }
 587             _M_x[__i] = __detail::__mod<_UIntType,
 588               __detail::_Shift<_UIntType, __w>::__value>(__sum);
 589           }
 590         _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 591         _M_p = 0;
 592       }
 593
 594   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
 595     typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 596              result_type
 597     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
 598     operator()()
 599     {
 600       // Derive short lag index from current index.
 601       long __ps = _M_p - short_lag;
 602       if (__ps < 0)
 603         __ps += long_lag;
 604
 605       // Calculate new x(i) without overflow or division.
 606       // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
 607       // cannot overflow.
 608       _UIntType __xi;
 609       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
 610         {
 611           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
 612           _M_carry = 0;
 613         }
 614       else
 615         {
 616           __xi = (__detail::_Shift<_UIntType, __w>::__value
 617                   - _M_x[_M_p] - _M_carry + _M_x[__ps]);
 618           _M_carry = 1;
 619         }
 620       _M_x[_M_p] = __xi;
 621
 622       // Adjust current index to loop around in ring buffer.
 623       if (++_M_p >= long_lag)
 624         _M_p = 0;
 625
 626       return __xi;
 627     }
 628
 629   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
 630            typename _CharT, typename _Traits>
 631     std::basic_ostream<_CharT, _Traits>&
 632     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 633                const subtract_with_carry_engine<_UIntType,
 634                                                 __w, __s, __r>& __x)
 635     {
 636       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 637       typedef typename __ostream_type::ios_base    __ios_base;
 638
 639       const typename __ios_base::fmtflags __flags = __os.flags();
 640       const _CharT __fill = __os.fill();
 641       const _CharT __space = __os.widen(' ');
 642       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 643       __os.fill(__space);
 644
 645       for (size_t __i = 0; __i < __r; ++__i)
 646         __os << __x._M_x[__i] << __space;
 647       __os << __x._M_carry << __space << __x._M_p;
 648
 649       __os.flags(__flags);
 650       __os.fill(__fill);
 651       return __os;
 652     }
 653
 654   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
 655            typename _CharT, typename _Traits>
 656     std::basic_istream<_CharT, _Traits>&
 657     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 658                subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
 659     {
 660       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
 661       typedef typename __istream_type::ios_base    __ios_base;
 662
 663       const typename __ios_base::fmtflags __flags = __is.flags();
 664       __is.flags(__ios_base::dec | __ios_base::skipws);
 665
 666       for (size_t __i = 0; __i < __r; ++__i)
 667         __is >> __x._M_x[__i];
 668       __is >> __x._M_carry;
 669       __is >> __x._M_p;
 670
 671       __is.flags(__flags);
 672       return __is;
 673     }
 674
 675
 676   template<typename _RandomNumberEngine, size_t __p, size_t __r>
 677     constexpr size_t
 678     discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
 679
 680   template<typename _RandomNumberEngine, size_t __p, size_t __r>
 681     constexpr size_t
 682     discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
 683
 684   template<typename _RandomNumberEngine, size_t __p, size_t __r>
 685     typename discard_block_engine<_RandomNumberEngine,
 686                            __p, __r>::result_type
 687     discard_block_engine<_RandomNumberEngine, __p, __r>::
 688     operator()()
 689     {
 690       if (_M_n >= used_block)
 691         {
 692           _M_b.discard(block_size - _M_n);
 693           _M_n = 0;
 694         }
 695       ++_M_n;
 696       return _M_b();
 697     }
 698
 699   template<typename _RandomNumberEngine, size_t __p, size_t __r,
 700            typename _CharT, typename _Traits>
 701     std::basic_ostream<_CharT, _Traits>&
 702     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 703                const discard_block_engine<_RandomNumberEngine,
 704                __p, __r>& __x)
 705     {
 706       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 707       typedef typename __ostream_type::ios_base    __ios_base;
 708
 709       const typename __ios_base::fmtflags __flags = __os.flags();
 710       const _CharT __fill = __os.fill();
 711       const _CharT __space = __os.widen(' ');
 712       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 713       __os.fill(__space);
 714
 715       __os << __x.base() << __space << __x._M_n;
 716
 717       __os.flags(__flags);
 718       __os.fill(__fill);
 719       return __os;
 720     }
 721
 722   template<typename _RandomNumberEngine, size_t __p, size_t __r,
 723            typename _CharT, typename _Traits>
 724     std::basic_istream<_CharT, _Traits>&
 725     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 726                discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
 727     {
 728       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 729       typedef typename __istream_type::ios_base    __ios_base;
 730
 731       const typename __ios_base::fmtflags __flags = __is.flags();
 732       __is.flags(__ios_base::dec | __ios_base::skipws);
 733
 734       __is >> __x._M_b >> __x._M_n;
 735
 736       __is.flags(__flags);
 737       return __is;
 738     }
 739
 740
 741   template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
 742     typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
 743       result_type
 744     independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
 745     operator()()
 746     {
 747       typedef typename _RandomNumberEngine::result_type _Eresult_type;
 748       const _Eresult_type __r
 749         = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
 750            ? _M_b.max() - _M_b.min() + 1 : 0);
 751       const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
 752       const unsigned __m = __r ? std::__lg(__r) : __edig;
 753
 754       typedef typename std::common_type<_Eresult_type, result_type>::type
 755         __ctype;
 756       const unsigned __cdig = std::numeric_limits<__ctype>::digits;
 757
 758       unsigned __n, __n0;
 759       __ctype __s0, __s1, __y0, __y1;
 760
 761       for (size_t __i = 0; __i < 2; ++__i)
 762         {
 763           __n = (__w + __m - 1) / __m + __i;
 764           __n0 = __n - __w % __n;
 765           const unsigned __w0 = __w / __n;  // __w0 <= __m
 766
 767           __s0 = 0;
 768           __s1 = 0;
 769           if (__w0 < __cdig)
 770             {
 771               __s0 = __ctype(1) << __w0;
 772               __s1 = __s0 << 1;
 773             }
 774
 775           __y0 = 0;
 776           __y1 = 0;
 777           if (__r)
 778             {
 779               __y0 = __s0 * (__r / __s0);
 780               if (__s1)
 781                 __y1 = __s1 * (__r / __s1);
 782
 783               if (__r - __y0 <= __y0 / __n)
 784                 break;
 785             }
 786           else
 787             break;
 788         }
 789
 790       result_type __sum = 0;
 791       for (size_t __k = 0; __k < __n0; ++__k)
 792         {
 793           __ctype __u;
 794           do
 795             __u = _M_b() - _M_b.min();
 796           while (__y0 && __u >= __y0);
 797           __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
 798         }
 799       for (size_t __k = __n0; __k < __n; ++__k)
 800         {
 801           __ctype __u;
 802           do
 803             __u = _M_b() - _M_b.min();
 804           while (__y1 && __u >= __y1);
 805           __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
 806         }
 807       return __sum;
 808     }
 809
 810
 811   template<typename _RandomNumberEngine, size_t __k>
 812     constexpr size_t
 813     shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
 814
 815   template<typename _RandomNumberEngine, size_t __k>
 816     typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
 817     shuffle_order_engine<_RandomNumberEngine, __k>::
 818     operator()()
 819     {
 820       size_t __j = __k * ((_M_y - _M_b.min())
 821                           / (_M_b.max() - _M_b.min() + 1.0L));
 822       _M_y = _M_v[__j];
 823       _M_v[__j] = _M_b();
 824
 825       return _M_y;
 826     }
 827
 828   template<typename _RandomNumberEngine, size_t __k,
 829            typename _CharT, typename _Traits>
 830     std::basic_ostream<_CharT, _Traits>&
 831     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 832                const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
 833     {
 834       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 835       typedef typename __ostream_type::ios_base    __ios_base;
 836
 837       const typename __ios_base::fmtflags __flags = __os.flags();
 838       const _CharT __fill = __os.fill();
 839       const _CharT __space = __os.widen(' ');
 840       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 841       __os.fill(__space);
 842
 843       __os << __x.base();
 844       for (size_t __i = 0; __i < __k; ++__i)
 845         __os << __space << __x._M_v[__i];
 846       __os << __space << __x._M_y;
 847
 848       __os.flags(__flags);
 849       __os.fill(__fill);
 850       return __os;
 851     }
 852
 853   template<typename _RandomNumberEngine, size_t __k,
 854            typename _CharT, typename _Traits>
 855     std::basic_istream<_CharT, _Traits>&
 856     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 857                shuffle_order_engine<_RandomNumberEngine, __k>& __x)
 858     {
 859       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 860       typedef typename __istream_type::ios_base    __ios_base;
 861
 862       const typename __ios_base::fmtflags __flags = __is.flags();
 863       __is.flags(__ios_base::dec | __ios_base::skipws);
 864
 865       __is >> __x._M_b;
 866       for (size_t __i = 0; __i < __k; ++__i)
 867         __is >> __x._M_v[__i];
 868       __is >> __x._M_y;
 869
 870       __is.flags(__flags);
 871       return __is;
 872     }
 873
 874
 875   template<typename _IntType>
 876     template<typename _UniformRandomNumberGenerator>
 877       typename uniform_int_distribution<_IntType>::result_type
 878       uniform_int_distribution<_IntType>::
 879       operator()(_UniformRandomNumberGenerator& __urng,
 880                  const param_type& __param)
 881       {
 882         typedef typename _UniformRandomNumberGenerator::result_type
 883           _Gresult_type;
 884         typedef typename std::make_unsigned<result_type>::type __utype;
 885         typedef typename std::common_type<_Gresult_type, __utype>::type
 886           __uctype;
 887
 888         const __uctype __urngmin = __urng.min();
 889         const __uctype __urngmax = __urng.max();
 890         const __uctype __urngrange = __urngmax - __urngmin;
 891         const __uctype __urange
 892           = __uctype(__param.b()) - __uctype(__param.a());
 893
 894         __uctype __ret;
 895
 896         if (__urngrange > __urange)
 897           {
 898             // downscaling
 899             const __uctype __uerange = __urange + 1; // __urange can be zero
 900             const __uctype __scaling = __urngrange / __uerange;
 901             const __uctype __past = __uerange * __scaling;
 902             do
 903               __ret = __uctype(__urng()) - __urngmin;
 904             while (__ret >= __past);
 905             __ret /= __scaling;
 906           }
 907         else if (__urngrange < __urange)
 908           {
 909             // upscaling
 910             /*
 911               Note that every value in [0, urange]
 912               can be written uniquely as
 913
 914               (urngrange + 1) * high + low
 915
 916               where
 917
 918               high in [0, urange / (urngrange + 1)]
 919
 920               and
 921
 922               low in [0, urngrange].
 923             */
 924             __uctype __tmp; // wraparound control
 925             do
 926               {
 927                 const __uctype __uerngrange = __urngrange + 1;
 928                 __tmp = (__uerngrange * operator()
 929                          (__urng, param_type(0, __urange / __uerngrange)));
 930                 __ret = __tmp + (__uctype(__urng()) - __urngmin);
 931               }
 932             while (__ret > __urange || __ret < __tmp);
 933           }
 934         else
 935           __ret = __uctype(__urng()) - __urngmin;
 936
 937         return __ret + __param.a();
 938       }
 939
 940
 941   template<typename _IntType>
 942     template<typename _ForwardIterator,
 943              typename _UniformRandomNumberGenerator>
 944       void
 945       uniform_int_distribution<_IntType>::
 946       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
 947                       _UniformRandomNumberGenerator& __urng,
 948                       const param_type& __param)
 949       {
 950         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
 951         typedef typename _UniformRandomNumberGenerator::result_type
 952           _Gresult_type;
 953         typedef typename std::make_unsigned<result_type>::type __utype;
 954         typedef typename std::common_type<_Gresult_type, __utype>::type
 955           __uctype;
 956
 957         const __uctype __urngmin = __urng.min();
 958         const __uctype __urngmax = __urng.max();
 959         const __uctype __urngrange = __urngmax - __urngmin;
 960         const __uctype __urange
 961           = __uctype(__param.b()) - __uctype(__param.a());
 962
 963         __uctype __ret;
 964
 965         if (__urngrange > __urange)
 966           {
 967             if (__detail::_Power_of_2(__urngrange + 1)
 968                 && __detail::_Power_of_2(__urange + 1))
 969               {
 970                 while (__f != __t)
 971                   {
 972                     __ret = __uctype(__urng()) - __urngmin;
 973                     *__f++ = (__ret & __urange) + __param.a();
 974                   }
 975               }
 976             else
 977               {
 978                 // downscaling
 979                 const __uctype __uerange = __urange + 1; // __urange can be zero
 980                 const __uctype __scaling = __urngrange / __uerange;
 981                 const __uctype __past = __uerange * __scaling;
 982                 while (__f != __t)
 983                   {
 984                     do
 985                       __ret = __uctype(__urng()) - __urngmin;
 986                     while (__ret >= __past);
 987                     *__f++ = __ret / __scaling + __param.a();
 988                   }
 989               }
 990           }
 991         else if (__urngrange < __urange)
 992           {
 993             // upscaling
 994             /*
 995               Note that every value in [0, urange]
 996               can be written uniquely as
 997
 998               (urngrange + 1) * high + low
 999
1000               where
1001
1002               high in [0, urange / (urngrange + 1)]
1003
1004               and
1005
1006               low in [0, urngrange].
1007             */
1008             __uctype __tmp; // wraparound control
1009             while (__f != __t)
1010               {
1011                 do
1012                   {
1013                     const __uctype __uerngrange = __urngrange + 1;
1014                     __tmp = (__uerngrange * operator()
1015                              (__urng, param_type(0, __urange / __uerngrange)));
1016                     __ret = __tmp + (__uctype(__urng()) - __urngmin);
1017                   }
1018                 while (__ret > __urange || __ret < __tmp);
1019                 *__f++ = __ret;
1020               }
1021           }
1022         else
1023           while (__f != __t)
1024             *__f++ = __uctype(__urng()) - __urngmin + __param.a();
1025       }
1026
1027   template<typename _IntType, typename _CharT, typename _Traits>
1028     std::basic_ostream<_CharT, _Traits>&
1029     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1030                const uniform_int_distribution<_IntType>& __x)
1031     {
1032       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1033       typedef typename __ostream_type::ios_base    __ios_base;
1034
1035       const typename __ios_base::fmtflags __flags = __os.flags();
1036       const _CharT __fill = __os.fill();
1037       const _CharT __space = __os.widen(' ');
1038       __os.flags(__ios_base::scientific | __ios_base::left);
1039       __os.fill(__space);
1040
1041       __os << __x.a() << __space << __x.b();
1042
1043       __os.flags(__flags);
1044       __os.fill(__fill);
1045       return __os;
1046     }
1047
1048   template<typename _IntType, typename _CharT, typename _Traits>
1049     std::basic_istream<_CharT, _Traits>&
1050     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1051                uniform_int_distribution<_IntType>& __x)
1052     {
1053       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1054       typedef typename __istream_type::ios_base    __ios_base;
1055
1056       const typename __ios_base::fmtflags __flags = __is.flags();
1057       __is.flags(__ios_base::dec | __ios_base::skipws);
1058
1059       _IntType __a, __b;
1060       __is >> __a >> __b;
1061       __x.param(typename uniform_int_distribution<_IntType>::
1062                 param_type(__a, __b));
1063
1064       __is.flags(__flags);
1065       return __is;
1066     }
1067
1068
1069   template<typename _RealType>
1070     template<typename _ForwardIterator,
1071              typename _UniformRandomNumberGenerator>
1072       void
1073       uniform_real_distribution<_RealType>::
1074       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1075                       _UniformRandomNumberGenerator& __urng,
1076                       const param_type& __p)
1077       {
1078         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1079         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1080           __aurng(__urng);
1081         auto __range = __p.b() - __p.a();
1082         while (__f != __t)
1083           *__f++ = __aurng() * __range + __p.a();
1084       }
1085
1086   template<typename _RealType, typename _CharT, typename _Traits>
1087     std::basic_ostream<_CharT, _Traits>&
1088     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1089                const uniform_real_distribution<_RealType>& __x)
1090     {
1091       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1092       typedef typename __ostream_type::ios_base    __ios_base;
1093
1094       const typename __ios_base::fmtflags __flags = __os.flags();
1095       const _CharT __fill = __os.fill();
1096       const std::streamsize __precision = __os.precision();
1097       const _CharT __space = __os.widen(' ');
1098       __os.flags(__ios_base::scientific | __ios_base::left);
1099       __os.fill(__space);
1100       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1101
1102       __os << __x.a() << __space << __x.b();
1103
1104       __os.flags(__flags);
1105       __os.fill(__fill);
1106       __os.precision(__precision);
1107       return __os;
1108     }
1109
1110   template<typename _RealType, typename _CharT, typename _Traits>
1111     std::basic_istream<_CharT, _Traits>&
1112     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1113                uniform_real_distribution<_RealType>& __x)
1114     {
1115       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1116       typedef typename __istream_type::ios_base    __ios_base;
1117
1118       const typename __ios_base::fmtflags __flags = __is.flags();
1119       __is.flags(__ios_base::skipws);
1120
1121       _RealType __a, __b;
1122       __is >> __a >> __b;
1123       __x.param(typename uniform_real_distribution<_RealType>::
1124                 param_type(__a, __b));
1125
1126       __is.flags(__flags);
1127       return __is;
1128     }
1129
1130
1131   template<typename _ForwardIterator,
1132            typename _UniformRandomNumberGenerator>
1133     void
1134     std::bernoulli_distribution::
1135     __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1136                     _UniformRandomNumberGenerator& __urng,
1137                     const param_type& __p)
1138     {
1139       __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1140       __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1141         __aurng(__urng);
1142       auto __limit = __p.p() * (__aurng.max() - __aurng.min());
1143
1144       while (__f != __t)
1145         *__f++ = (__aurng() - __aurng.min()) < __limit;
1146     }
1147
1148   template<typename _CharT, typename _Traits>
1149     std::basic_ostream<_CharT, _Traits>&
1150     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1151                const bernoulli_distribution& __x)
1152     {
1153       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1154       typedef typename __ostream_type::ios_base    __ios_base;
1155
1156       const typename __ios_base::fmtflags __flags = __os.flags();
1157       const _CharT __fill = __os.fill();
1158       const std::streamsize __precision = __os.precision();
1159       __os.flags(__ios_base::scientific | __ios_base::left);
1160       __os.fill(__os.widen(' '));
1161       __os.precision(std::numeric_limits<double>::max_digits10);
1162
1163       __os << __x.p();
1164
1165       __os.flags(__flags);
1166       __os.fill(__fill);
1167       __os.precision(__precision);
1168       return __os;
1169     }
1170
1171
1172   template<typename _IntType>
1173     template<typename _UniformRandomNumberGenerator>
1174       typename geometric_distribution<_IntType>::result_type
1175       geometric_distribution<_IntType>::
1176       operator()(_UniformRandomNumberGenerator& __urng,
1177                  const param_type& __param)
1178       {
1179         // About the epsilon thing see this thread:
1180         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1181         const double __naf =
1182           (1 - std::numeric_limits<double>::epsilon()) / 2;
1183         // The largest _RealType convertible to _IntType.
1184         const double __thr =
1185           std::numeric_limits<_IntType>::max() + __naf;
1186         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1187           __aurng(__urng);
1188
1189         double __cand;
1190         do
1191           __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1192         while (__cand >= __thr);
1193
1194         return result_type(__cand + __naf);
1195       }
1196
1197   template<typename _IntType>
1198     template<typename _ForwardIterator,
1199              typename _UniformRandomNumberGenerator>
1200       void
1201       geometric_distribution<_IntType>::
1202       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1203                       _UniformRandomNumberGenerator& __urng,
1204                       const param_type& __param)
1205       {
1206         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1207         // About the epsilon thing see this thread:
1208         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1209         const double __naf =
1210           (1 - std::numeric_limits<double>::epsilon()) / 2;
1211         // The largest _RealType convertible to _IntType.
1212         const double __thr =
1213           std::numeric_limits<_IntType>::max() + __naf;
1214         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1215           __aurng(__urng);
1216
1217         while (__f != __t)
1218           {
1219             double __cand;
1220             do
1221               __cand = std::floor(std::log(1.0 - __aurng())
1222                                   / __param._M_log_1_p);
1223             while (__cand >= __thr);
1224
1225             *__f++ = __cand + __naf;
1226           }
1227       }
1228
1229   template<typename _IntType,
1230            typename _CharT, typename _Traits>
1231     std::basic_ostream<_CharT, _Traits>&
1232     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1233                const geometric_distribution<_IntType>& __x)
1234     {
1235       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1236       typedef typename __ostream_type::ios_base    __ios_base;
1237
1238       const typename __ios_base::fmtflags __flags = __os.flags();
1239       const _CharT __fill = __os.fill();
1240       const std::streamsize __precision = __os.precision();
1241       __os.flags(__ios_base::scientific | __ios_base::left);
1242       __os.fill(__os.widen(' '));
1243       __os.precision(std::numeric_limits<double>::max_digits10);
1244
1245       __os << __x.p();
1246
1247       __os.flags(__flags);
1248       __os.fill(__fill);
1249       __os.precision(__precision);
1250       return __os;
1251     }
1252
1253   template<typename _IntType,
1254            typename _CharT, typename _Traits>
1255     std::basic_istream<_CharT, _Traits>&
1256     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1257                geometric_distribution<_IntType>& __x)
1258     {
1259       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1260       typedef typename __istream_type::ios_base    __ios_base;
1261
1262       const typename __ios_base::fmtflags __flags = __is.flags();
1263       __is.flags(__ios_base::skipws);
1264
1265       double __p;
1266       __is >> __p;
1267       __x.param(typename geometric_distribution<_IntType>::param_type(__p));
1268
1269       __is.flags(__flags);
1270       return __is;
1271     }
1272
1273   // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1274   template<typename _IntType>
1275     template<typename _UniformRandomNumberGenerator>
1276       typename negative_binomial_distribution<_IntType>::result_type
1277       negative_binomial_distribution<_IntType>::
1278       operator()(_UniformRandomNumberGenerator& __urng)
1279       {
1280         const double __y = _M_gd(__urng);
1281
1282         // XXX Is the constructor too slow?
1283         std::poisson_distribution<result_type> __poisson(__y);
1284         return __poisson(__urng);
1285       }
1286
1287   template<typename _IntType>
1288     template<typename _UniformRandomNumberGenerator>
1289       typename negative_binomial_distribution<_IntType>::result_type
1290       negative_binomial_distribution<_IntType>::
1291       operator()(_UniformRandomNumberGenerator& __urng,
1292                  const param_type& __p)
1293       {
1294         typedef typename std::gamma_distribution<double>::param_type
1295           param_type;
1296
1297         const double __y =
1298           _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1299
1300         std::poisson_distribution<result_type> __poisson(__y);
1301         return __poisson(__urng);
1302       }
1303
1304   template<typename _IntType>
1305     template<typename _ForwardIterator,
1306              typename _UniformRandomNumberGenerator>
1307       void
1308       negative_binomial_distribution<_IntType>::
1309       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1310                       _UniformRandomNumberGenerator& __urng)
1311       {
1312         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1313         while (__f != __t)
1314           {
1315             const double __y = _M_gd(__urng);
1316
1317             // XXX Is the constructor too slow?
1318             std::poisson_distribution<result_type> __poisson(__y);
1319             *__f++ = __poisson(__urng);
1320           }
1321       }
1322
1323   template<typename _IntType>
1324     template<typename _ForwardIterator,
1325              typename _UniformRandomNumberGenerator>
1326       void
1327       negative_binomial_distribution<_IntType>::
1328       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1329                       _UniformRandomNumberGenerator& __urng,
1330                       const param_type& __p)
1331       {
1332         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1333         typename std::gamma_distribution<result_type>::param_type
1334           __p2(__p.k(), (1.0 - __p.p()) / __p.p());
1335
1336         while (__f != __t)
1337           {
1338             const double __y = _M_gd(__urng, __p2);
1339
1340             std::poisson_distribution<result_type> __poisson(__y);
1341             *__f++ = __poisson(__urng);
1342           }
1343       }
1344
1345   template<typename _IntType, typename _CharT, typename _Traits>
1346     std::basic_ostream<_CharT, _Traits>&
1347     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1348                const negative_binomial_distribution<_IntType>& __x)
1349     {
1350       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1351       typedef typename __ostream_type::ios_base    __ios_base;
1352
1353       const typename __ios_base::fmtflags __flags = __os.flags();
1354       const _CharT __fill = __os.fill();
1355       const std::streamsize __precision = __os.precision();
1356       const _CharT __space = __os.widen(' ');
1357       __os.flags(__ios_base::scientific | __ios_base::left);
1358       __os.fill(__os.widen(' '));
1359       __os.precision(std::numeric_limits<double>::max_digits10);
1360
1361       __os << __x.k() << __space << __x.p()
1362            << __space << __x._M_gd;
1363
1364       __os.flags(__flags);
1365       __os.fill(__fill);
1366       __os.precision(__precision);
1367       return __os;
1368     }
1369
1370   template<typename _IntType, typename _CharT, typename _Traits>
1371     std::basic_istream<_CharT, _Traits>&
1372     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1373                negative_binomial_distribution<_IntType>& __x)
1374     {
1375       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1376       typedef typename __istream_type::ios_base    __ios_base;
1377
1378       const typename __ios_base::fmtflags __flags = __is.flags();
1379       __is.flags(__ios_base::skipws);
1380
1381       _IntType __k;
1382       double __p;
1383       __is >> __k >> __p >> __x._M_gd;
1384       __x.param(typename negative_binomial_distribution<_IntType>::
1385                 param_type(__k, __p));
1386
1387       __is.flags(__flags);
1388       return __is;
1389     }
1390
1391
1392   template<typename _IntType>
1393     void
1394     poisson_distribution<_IntType>::param_type::
1395     _M_initialize()
1396     {
1397 #if _GLIBCXX_USE_C99_MATH_TR1
1398       if (_M_mean >= 12)
1399         {
1400           const double __m = std::floor(_M_mean);
1401           _M_lm_thr = std::log(_M_mean);
1402           _M_lfm = std::lgamma(__m + 1);
1403           _M_sm = std::sqrt(__m);
1404
1405           const double __pi_4 = 0.7853981633974483096156608458198757L;
1406           const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1407                                                               / __pi_4));
1408           _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
1409           const double __cx = 2 * __m + _M_d;
1410           _M_scx = std::sqrt(__cx / 2);
1411           _M_1cx = 1 / __cx;
1412
1413           _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1414           _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1415                 / _M_d;
1416         }
1417       else
1418 #endif
1419         _M_lm_thr = std::exp(-_M_mean);
1420       }
1421
1422   /**
1423    * A rejection algorithm when mean >= 12 and a simple method based
1424    * upon the multiplication of uniform random variates otherwise.
1425    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1426    * is defined.
1427    *
1428    * Reference:
1429    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1430    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1431    */
1432   template<typename _IntType>
1433     template<typename _UniformRandomNumberGenerator>
1434       typename poisson_distribution<_IntType>::result_type
1435       poisson_distribution<_IntType>::
1436       operator()(_UniformRandomNumberGenerator& __urng,
1437                  const param_type& __param)
1438       {
1439         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1440           __aurng(__urng);
1441 #if _GLIBCXX_USE_C99_MATH_TR1
1442         if (__param.mean() >= 12)
1443           {
1444             double __x;
1445
1446             // See comments above...
1447             const double __naf =
1448               (1 - std::numeric_limits<double>::epsilon()) / 2;
1449             const double __thr =
1450               std::numeric_limits<_IntType>::max() + __naf;
1451
1452             const double __m = std::floor(__param.mean());
1453             // sqrt(pi / 2)
1454             const double __spi_2 = 1.2533141373155002512078826424055226L;
1455             const double __c1 = __param._M_sm * __spi_2;
1456             const double __c2 = __param._M_c2b + __c1;
1457             const double __c3 = __c2 + 1;
1458             const double __c4 = __c3 + 1;
1459             // e^(1 / 78)
1460             const double __e178 = 1.0129030479320018583185514777512983L;
1461             const double __c5 = __c4 + __e178;
1462             const double __c = __param._M_cb + __c5;
1463             const double __2cx = 2 * (2 * __m + __param._M_d);
1464
1465             bool __reject = true;
1466             do
1467               {
1468                 const double __u = __c * __aurng();
1469                 const double __e = -std::log(1.0 - __aurng());
1470
1471                 double __w = 0.0;
1472
1473                 if (__u <= __c1)
1474                   {
1475                     const double __n = _M_nd(__urng);
1476                     const double __y = -std::abs(__n) * __param._M_sm - 1;
1477                     __x = std::floor(__y);
1478                     __w = -__n * __n / 2;
1479                     if (__x < -__m)
1480                       continue;
1481                   }
1482                 else if (__u <= __c2)
1483                   {
1484                     const double __n = _M_nd(__urng);
1485                     const double __y = 1 + std::abs(__n) * __param._M_scx;
1486                     __x = std::ceil(__y);
1487                     __w = __y * (2 - __y) * __param._M_1cx;
1488                     if (__x > __param._M_d)
1489                       continue;
1490                   }
1491                 else if (__u <= __c3)
1492                   // NB: This case not in the book, nor in the Errata,
1493                   // but should be ok...
1494                   __x = -1;
1495                 else if (__u <= __c4)
1496                   __x = 0;
1497                 else if (__u <= __c5)
1498                   __x = 1;
1499                 else
1500                   {
1501                     const double __v = -std::log(1.0 - __aurng());
1502                     const double __y = __param._M_d
1503                                      + __v * __2cx / __param._M_d;
1504                     __x = std::ceil(__y);
1505                     __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1506                   }
1507
1508                 __reject = (__w - __e - __x * __param._M_lm_thr
1509                             > __param._M_lfm - std::lgamma(__x + __m + 1));
1510
1511                 __reject |= __x + __m >= __thr;
1512
1513               } while (__reject);
1514
1515             return result_type(__x + __m + __naf);
1516           }
1517         else
1518 #endif
1519           {
1520             _IntType     __x = 0;
1521             double __prod = 1.0;
1522
1523             do
1524               {
1525                 __prod *= __aurng();
1526                 __x += 1;
1527               }
1528             while (__prod > __param._M_lm_thr);
1529
1530             return __x - 1;
1531           }
1532       }
1533
1534   template<typename _IntType>
1535     template<typename _ForwardIterator,
1536              typename _UniformRandomNumberGenerator>
1537       void
1538       poisson_distribution<_IntType>::
1539       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1540                       _UniformRandomNumberGenerator& __urng,
1541                       const param_type& __param)
1542       {
1543         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1544         // We could duplicate everything from operator()...
1545         while (__f != __t)
1546           *__f++ = this->operator()(__urng, __param);
1547       }
1548
1549   template<typename _IntType,
1550            typename _CharT, typename _Traits>
1551     std::basic_ostream<_CharT, _Traits>&
1552     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1553                const poisson_distribution<_IntType>& __x)
1554     {
1555       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1556       typedef typename __ostream_type::ios_base    __ios_base;
1557
1558       const typename __ios_base::fmtflags __flags = __os.flags();
1559       const _CharT __fill = __os.fill();
1560       const std::streamsize __precision = __os.precision();
1561       const _CharT __space = __os.widen(' ');
1562       __os.flags(__ios_base::scientific | __ios_base::left);
1563       __os.fill(__space);
1564       __os.precision(std::numeric_limits<double>::max_digits10);
1565
1566       __os << __x.mean() << __space << __x._M_nd;
1567
1568       __os.flags(__flags);
1569       __os.fill(__fill);
1570       __os.precision(__precision);
1571       return __os;
1572     }
1573
1574   template<typename _IntType,
1575            typename _CharT, typename _Traits>
1576     std::basic_istream<_CharT, _Traits>&
1577     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1578                poisson_distribution<_IntType>& __x)
1579     {
1580       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1581       typedef typename __istream_type::ios_base    __ios_base;
1582
1583       const typename __ios_base::fmtflags __flags = __is.flags();
1584       __is.flags(__ios_base::skipws);
1585
1586       double __mean;
1587       __is >> __mean >> __x._M_nd;
1588       __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1589
1590       __is.flags(__flags);
1591       return __is;
1592     }
1593
1594
1595   template<typename _IntType>
1596     void
1597     binomial_distribution<_IntType>::param_type::
1598     _M_initialize()
1599     {
1600       const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1601
1602       _M_easy = true;
1603
1604 #if _GLIBCXX_USE_C99_MATH_TR1
1605       if (_M_t * __p12 >= 8)
1606         {
1607           _M_easy = false;
1608           const double __np = std::floor(_M_t * __p12);
1609           const double __pa = __np / _M_t;
1610           const double __1p = 1 - __pa;
1611
1612           const double __pi_4 = 0.7853981633974483096156608458198757L;
1613           const double __d1x =
1614             std::sqrt(__np * __1p * std::log(32 * __np
1615                                              / (81 * __pi_4 * __1p)));
1616           _M_d1 = std::round(std::max(1.0, __d1x));
1617           const double __d2x =
1618             std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1619                                              / (__pi_4 * __pa)));
1620           _M_d2 = std::round(std::max(1.0, __d2x));
1621
1622           // sqrt(pi / 2)
1623           const double __spi_2 = 1.2533141373155002512078826424055226L;
1624           _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1625           _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1626           _M_c = 2 * _M_d1 / __np;
1627           _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1628           const double __a12 = _M_a1 + _M_s2 * __spi_2;
1629           const double __s1s = _M_s1 * _M_s1;
1630           _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1631                              * 2 * __s1s / _M_d1
1632                              * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1633           const double __s2s = _M_s2 * _M_s2;
1634           _M_s = (_M_a123 + 2 * __s2s / _M_d2
1635                   * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1636           _M_lf = (std::lgamma(__np + 1)
1637                    + std::lgamma(_M_t - __np + 1));
1638           _M_lp1p = std::log(__pa / __1p);
1639
1640           _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1641         }
1642       else
1643 #endif
1644         _M_q = -std::log(1 - __p12);
1645     }
1646
1647   template<typename _IntType>
1648     template<typename _UniformRandomNumberGenerator>
1649       typename binomial_distribution<_IntType>::result_type
1650       binomial_distribution<_IntType>::
1651       _M_waiting(_UniformRandomNumberGenerator& __urng,
1652                  _IntType __t, double __q)
1653       {
1654         _IntType __x = 0;
1655         double __sum = 0.0;
1656         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1657           __aurng(__urng);
1658
1659         do
1660           {
1661             if (__t == __x)
1662               return __x;
1663             const double __e = -std::log(1.0 - __aurng());
1664             __sum += __e / (__t - __x);
1665             __x += 1;
1666           }
1667         while (__sum <= __q);
1668
1669         return __x - 1;
1670       }
1671
1672   /**
1673    * A rejection algorithm when t * p >= 8 and a simple waiting time
1674    * method - the second in the referenced book - otherwise.
1675    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1676    * is defined.
1677    *
1678    * Reference:
1679    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1680    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1681    */
1682   template<typename _IntType>
1683     template<typename _UniformRandomNumberGenerator>
1684       typename binomial_distribution<_IntType>::result_type
1685       binomial_distribution<_IntType>::
1686       operator()(_UniformRandomNumberGenerator& __urng,
1687                  const param_type& __param)
1688       {
1689         result_type __ret;
1690         const _IntType __t = __param.t();
1691         const double __p = __param.p();
1692         const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1693         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1694           __aurng(__urng);
1695
1696 #if _GLIBCXX_USE_C99_MATH_TR1
1697         if (!__param._M_easy)
1698           {
1699             double __x;
1700
1701             // See comments above...
1702             const double __naf =
1703               (1 - std::numeric_limits<double>::epsilon()) / 2;
1704             const double __thr =
1705               std::numeric_limits<_IntType>::max() + __naf;
1706
1707             const double __np = std::floor(__t * __p12);
1708
1709             // sqrt(pi / 2)
1710             const double __spi_2 = 1.2533141373155002512078826424055226L;
1711             const double __a1 = __param._M_a1;
1712             const double __a12 = __a1 + __param._M_s2 * __spi_2;
1713             const double __a123 = __param._M_a123;
1714             const double __s1s = __param._M_s1 * __param._M_s1;
1715             const double __s2s = __param._M_s2 * __param._M_s2;
1716
1717             bool __reject;
1718             do
1719               {
1720                 const double __u = __param._M_s * __aurng();
1721
1722                 double __v;
1723
1724                 if (__u <= __a1)
1725                   {
1726                     const double __n = _M_nd(__urng);
1727                     const double __y = __param._M_s1 * std::abs(__n);
1728                     __reject = __y >= __param._M_d1;
1729                     if (!__reject)
1730                       {
1731                         const double __e = -std::log(1.0 - __aurng());
1732                         __x = std::floor(__y);
1733                         __v = -__e - __n * __n / 2 + __param._M_c;
1734                       }
1735                   }
1736                 else if (__u <= __a12)
1737                   {
1738                     const double __n = _M_nd(__urng);
1739                     const double __y = __param._M_s2 * std::abs(__n);
1740                     __reject = __y >= __param._M_d2;
1741                     if (!__reject)
1742                       {
1743                         const double __e = -std::log(1.0 - __aurng());
1744                         __x = std::floor(-__y);
1745                         __v = -__e - __n * __n / 2;
1746                       }
1747                   }
1748                 else if (__u <= __a123)
1749                   {
1750                     const double __e1 = -std::log(1.0 - __aurng());
1751                     const double __e2 = -std::log(1.0 - __aurng());
1752
1753                     const double __y = __param._M_d1
1754                                      + 2 * __s1s * __e1 / __param._M_d1;
1755                     __x = std::floor(__y);
1756                     __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1757                                                     -__y / (2 * __s1s)));
1758                     __reject = false;
1759                   }
1760                 else
1761                   {
1762                     const double __e1 = -std::log(1.0 - __aurng());
1763                     const double __e2 = -std::log(1.0 - __aurng());
1764
1765                     const double __y = __param._M_d2
1766                                      + 2 * __s2s * __e1 / __param._M_d2;
1767                     __x = std::floor(-__y);
1768                     __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1769                     __reject = false;
1770                   }
1771
1772                 __reject = __reject || __x < -__np || __x > __t - __np;
1773                 if (!__reject)
1774                   {
1775                     const double __lfx =
1776                       std::lgamma(__np + __x + 1)
1777                       + std::lgamma(__t - (__np + __x) + 1);
1778                     __reject = __v > __param._M_lf - __lfx
1779                              + __x * __param._M_lp1p;
1780                   }
1781
1782                 __reject |= __x + __np >= __thr;
1783               }
1784             while (__reject);
1785
1786             __x += __np + __naf;
1787
1788             const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1789                                             __param._M_q);
1790             __ret = _IntType(__x) + __z;
1791           }
1792         else
1793 #endif
1794           __ret = _M_waiting(__urng, __t, __param._M_q);
1795
1796         if (__p12 != __p)
1797           __ret = __t - __ret;
1798         return __ret;
1799       }
1800
1801   template<typename _IntType>
1802     template<typename _ForwardIterator,
1803              typename _UniformRandomNumberGenerator>
1804       void
1805       binomial_distribution<_IntType>::
1806       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1807                       _UniformRandomNumberGenerator& __urng,
1808                       const param_type& __param)
1809       {
1810         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1811         // We could duplicate everything from operator()...
1812         while (__f != __t)
1813           *__f++ = this->operator()(__urng, __param);
1814       }
1815
1816   template<typename _IntType,
1817            typename _CharT, typename _Traits>
1818     std::basic_ostream<_CharT, _Traits>&
1819     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1820                const binomial_distribution<_IntType>& __x)
1821     {
1822       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1823       typedef typename __ostream_type::ios_base    __ios_base;
1824
1825       const typename __ios_base::fmtflags __flags = __os.flags();
1826       const _CharT __fill = __os.fill();
1827       const std::streamsize __precision = __os.precision();
1828       const _CharT __space = __os.widen(' ');
1829       __os.flags(__ios_base::scientific | __ios_base::left);
1830       __os.fill(__space);
1831       __os.precision(std::numeric_limits<double>::max_digits10);
1832
1833       __os << __x.t() << __space << __x.p()
1834            << __space << __x._M_nd;
1835
1836       __os.flags(__flags);
1837       __os.fill(__fill);
1838       __os.precision(__precision);
1839       return __os;
1840     }
1841
1842   template<typename _IntType,
1843            typename _CharT, typename _Traits>
1844     std::basic_istream<_CharT, _Traits>&
1845     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1846                binomial_distribution<_IntType>& __x)
1847     {
1848       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1849       typedef typename __istream_type::ios_base    __ios_base;
1850
1851       const typename __ios_base::fmtflags __flags = __is.flags();
1852       __is.flags(__ios_base::dec | __ios_base::skipws);
1853
1854       _IntType __t;
1855       double __p;
1856       __is >> __t >> __p >> __x._M_nd;
1857       __x.param(typename binomial_distribution<_IntType>::
1858                 param_type(__t, __p));
1859
1860       __is.flags(__flags);
1861       return __is;
1862     }
1863
1864
1865   template<typename _RealType>
1866     template<typename _ForwardIterator,
1867              typename _UniformRandomNumberGenerator>
1868       void
1869       std::exponential_distribution<_RealType>::
1870       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1871                       _UniformRandomNumberGenerator& __urng,
1872                       const param_type& __p)
1873       {
1874         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1875         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1876           __aurng(__urng);
1877         while (__f != __t)
1878           *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1879       }
1880
1881   template<typename _RealType, typename _CharT, typename _Traits>
1882     std::basic_ostream<_CharT, _Traits>&
1883     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1884                const exponential_distribution<_RealType>& __x)
1885     {
1886       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1887       typedef typename __ostream_type::ios_base    __ios_base;
1888
1889       const typename __ios_base::fmtflags __flags = __os.flags();
1890       const _CharT __fill = __os.fill();
1891       const std::streamsize __precision = __os.precision();
1892       __os.flags(__ios_base::scientific | __ios_base::left);
1893       __os.fill(__os.widen(' '));
1894       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1895
1896       __os << __x.lambda();
1897
1898       __os.flags(__flags);
1899       __os.fill(__fill);
1900       __os.precision(__precision);
1901       return __os;
1902     }
1903
1904   template<typename _RealType, typename _CharT, typename _Traits>
1905     std::basic_istream<_CharT, _Traits>&
1906     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1907                exponential_distribution<_RealType>& __x)
1908     {
1909       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1910       typedef typename __istream_type::ios_base    __ios_base;
1911
1912       const typename __ios_base::fmtflags __flags = __is.flags();
1913       __is.flags(__ios_base::dec | __ios_base::skipws);
1914
1915       _RealType __lambda;
1916       __is >> __lambda;
1917       __x.param(typename exponential_distribution<_RealType>::
1918                 param_type(__lambda));
1919
1920       __is.flags(__flags);
1921       return __is;
1922     }
1923
1924
1925   /**
1926    * Polar method due to Marsaglia.
1927    *
1928    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1929    * New York, 1986, Ch. V, Sect. 4.4.
1930    */
1931   template<typename _RealType>
1932     template<typename _UniformRandomNumberGenerator>
1933       typename normal_distribution<_RealType>::result_type
1934       normal_distribution<_RealType>::
1935       operator()(_UniformRandomNumberGenerator& __urng,
1936                  const param_type& __param)
1937       {
1938         result_type __ret;
1939         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1940           __aurng(__urng);
1941
1942         if (_M_saved_available)
1943           {
1944             _M_saved_available = false;
1945             __ret = _M_saved;
1946           }
1947         else
1948           {
1949             result_type __x, __y, __r2;
1950             do
1951               {
1952                 __x = result_type(2.0) * __aurng() - 1.0;
1953                 __y = result_type(2.0) * __aurng() - 1.0;
1954                 __r2 = __x * __x + __y * __y;
1955               }
1956             while (__r2 > 1.0 || __r2 == 0.0);
1957
1958             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1959             _M_saved = __x * __mult;
1960             _M_saved_available = true;
1961             __ret = __y * __mult;
1962           }
1963
1964         __ret = __ret * __param.stddev() + __param.mean();
1965         return __ret;
1966       }
1967
1968   template<typename _RealType>
1969     template<typename _ForwardIterator,
1970              typename _UniformRandomNumberGenerator>
1971       void
1972       normal_distribution<_RealType>::
1973       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1974                       _UniformRandomNumberGenerator& __urng,
1975                       const param_type& __param)
1976       {
1977         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1978
1979         if (__f == __t)
1980           return;
1981
1982         if (_M_saved_available)
1983           {
1984             _M_saved_available = false;
1985             *__f++ = _M_saved * __param.stddev() + __param.mean();
1986
1987             if (__f == __t)
1988               return;
1989           }
1990
1991         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1992           __aurng(__urng);
1993
1994         while (__f + 1 < __t)
1995           {
1996             result_type __x, __y, __r2;
1997             do
1998               {
1999                 __x = result_type(2.0) * __aurng() - 1.0;
2000                 __y = result_type(2.0) * __aurng() - 1.0;
2001                 __r2 = __x * __x + __y * __y;
2002               }
2003             while (__r2 > 1.0 || __r2 == 0.0);
2004
2005             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2006             *__f++ = __y * __mult * __param.stddev() + __param.mean();
2007             *__f++ = __x * __mult * __param.stddev() + __param.mean();
2008           }
2009
2010         if (__f != __t)
2011           {
2012             result_type __x, __y, __r2;
2013             do
2014               {
2015                 __x = result_type(2.0) * __aurng() - 1.0;
2016                 __y = result_type(2.0) * __aurng() - 1.0;
2017                 __r2 = __x * __x + __y * __y;
2018               }
2019             while (__r2 > 1.0 || __r2 == 0.0);
2020
2021             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
2022             _M_saved = __x * __mult;
2023             _M_saved_available = true;
2024             *__f = __y * __mult * __param.stddev() + __param.mean();
2025           }
2026       }
2027
2028   template<typename _RealType>
2029     bool
2030     operator==(const std::normal_distribution<_RealType>& __d1,
2031                const std::normal_distribution<_RealType>& __d2)
2032     {
2033       if (__d1._M_param == __d2._M_param
2034           && __d1._M_saved_available == __d2._M_saved_available)
2035         {
2036           if (__d1._M_saved_available
2037               && __d1._M_saved == __d2._M_saved)
2038             return true;
2039           else if(!__d1._M_saved_available)
2040             return true;
2041           else
2042             return false;
2043         }
2044       else
2045         return false;
2046     }
2047
2048   template<typename _RealType, typename _CharT, typename _Traits>
2049     std::basic_ostream<_CharT, _Traits>&
2050     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2051                const normal_distribution<_RealType>& __x)
2052     {
2053       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2054       typedef typename __ostream_type::ios_base    __ios_base;
2055
2056       const typename __ios_base::fmtflags __flags = __os.flags();
2057       const _CharT __fill = __os.fill();
2058       const std::streamsize __precision = __os.precision();
2059       const _CharT __space = __os.widen(' ');
2060       __os.flags(__ios_base::scientific | __ios_base::left);
2061       __os.fill(__space);
2062       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2063
2064       __os << __x.mean() << __space << __x.stddev()
2065            << __space << __x._M_saved_available;
2066       if (__x._M_saved_available)
2067         __os << __space << __x._M_saved;
2068
2069       __os.flags(__flags);
2070       __os.fill(__fill);
2071       __os.precision(__precision);
2072       return __os;
2073     }
2074
2075   template<typename _RealType, typename _CharT, typename _Traits>
2076     std::basic_istream<_CharT, _Traits>&
2077     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2078                normal_distribution<_RealType>& __x)
2079     {
2080       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2081       typedef typename __istream_type::ios_base    __ios_base;
2082
2083       const typename __ios_base::fmtflags __flags = __is.flags();
2084       __is.flags(__ios_base::dec | __ios_base::skipws);
2085
2086       double __mean, __stddev;
2087       __is >> __mean >> __stddev
2088            >> __x._M_saved_available;
2089       if (__x._M_saved_available)
2090         __is >> __x._M_saved;
2091       __x.param(typename normal_distribution<_RealType>::
2092                 param_type(__mean, __stddev));
2093
2094       __is.flags(__flags);
2095       return __is;
2096     }
2097
2098
2099   template<typename _RealType>
2100     template<typename _ForwardIterator,
2101              typename _UniformRandomNumberGenerator>
2102       void
2103       lognormal_distribution<_RealType>::
2104       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2105                       _UniformRandomNumberGenerator& __urng,
2106                       const param_type& __p)
2107       {
2108         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2109           while (__f != __t)
2110             *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
2111       }
2112
2113   template<typename _RealType, typename _CharT, typename _Traits>
2114     std::basic_ostream<_CharT, _Traits>&
2115     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2116                const lognormal_distribution<_RealType>& __x)
2117     {
2118       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2119       typedef typename __ostream_type::ios_base    __ios_base;
2120
2121       const typename __ios_base::fmtflags __flags = __os.flags();
2122       const _CharT __fill = __os.fill();
2123       const std::streamsize __precision = __os.precision();
2124       const _CharT __space = __os.widen(' ');
2125       __os.flags(__ios_base::scientific | __ios_base::left);
2126       __os.fill(__space);
2127       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2128
2129       __os << __x.m() << __space << __x.s()
2130            << __space << __x._M_nd;
2131
2132       __os.flags(__flags);
2133       __os.fill(__fill);
2134       __os.precision(__precision);
2135       return __os;
2136     }
2137
2138   template<typename _RealType, typename _CharT, typename _Traits>
2139     std::basic_istream<_CharT, _Traits>&
2140     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2141                lognormal_distribution<_RealType>& __x)
2142     {
2143       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2144       typedef typename __istream_type::ios_base    __ios_base;
2145
2146       const typename __ios_base::fmtflags __flags = __is.flags();
2147       __is.flags(__ios_base::dec | __ios_base::skipws);
2148
2149       _RealType __m, __s;
2150       __is >> __m >> __s >> __x._M_nd;
2151       __x.param(typename lognormal_distribution<_RealType>::
2152                 param_type(__m, __s));
2153
2154       __is.flags(__flags);
2155       return __is;
2156     }
2157
2158   template<typename _RealType>
2159     template<typename _ForwardIterator,
2160              typename _UniformRandomNumberGenerator>
2161       void
2162       std::chi_squared_distribution<_RealType>::
2163       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2164                       _UniformRandomNumberGenerator& __urng)
2165       {
2166         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2167         while (__f != __t)
2168           *__f++ = 2 * _M_gd(__urng);
2169       }
2170
2171   template<typename _RealType>
2172     template<typename _ForwardIterator,
2173              typename _UniformRandomNumberGenerator>
2174       void
2175       std::chi_squared_distribution<_RealType>::
2176       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2177                       _UniformRandomNumberGenerator& __urng,
2178                       const typename
2179                       std::gamma_distribution<result_type>::param_type& __p)
2180       {
2181         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2182         while (__f != __t)
2183           *__f++ = 2 * _M_gd(__urng, __p);
2184       }
2185
2186   template<typename _RealType, typename _CharT, typename _Traits>
2187     std::basic_ostream<_CharT, _Traits>&
2188     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2189                const chi_squared_distribution<_RealType>& __x)
2190     {
2191       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2192       typedef typename __ostream_type::ios_base    __ios_base;
2193
2194       const typename __ios_base::fmtflags __flags = __os.flags();
2195       const _CharT __fill = __os.fill();
2196       const std::streamsize __precision = __os.precision();
2197       const _CharT __space = __os.widen(' ');
2198       __os.flags(__ios_base::scientific | __ios_base::left);
2199       __os.fill(__space);
2200       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2201
2202       __os << __x.n() << __space << __x._M_gd;
2203
2204       __os.flags(__flags);
2205       __os.fill(__fill);
2206       __os.precision(__precision);
2207       return __os;
2208     }
2209
2210   template<typename _RealType, typename _CharT, typename _Traits>
2211     std::basic_istream<_CharT, _Traits>&
2212     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2213                chi_squared_distribution<_RealType>& __x)
2214     {
2215       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2216       typedef typename __istream_type::ios_base    __ios_base;
2217
2218       const typename __ios_base::fmtflags __flags = __is.flags();
2219       __is.flags(__ios_base::dec | __ios_base::skipws);
2220
2221       _RealType __n;
2222       __is >> __n >> __x._M_gd;
2223       __x.param(typename chi_squared_distribution<_RealType>::
2224                 param_type(__n));
2225
2226       __is.flags(__flags);
2227       return __is;
2228     }
2229
2230
2231   template<typename _RealType>
2232     template<typename _UniformRandomNumberGenerator>
2233       typename cauchy_distribution<_RealType>::result_type
2234       cauchy_distribution<_RealType>::
2235       operator()(_UniformRandomNumberGenerator& __urng,
2236                  const param_type& __p)
2237       {
2238         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2239           __aurng(__urng);
2240         _RealType __u;
2241         do
2242           __u = __aurng();
2243         while (__u == 0.5);
2244
2245         const _RealType __pi = 3.1415926535897932384626433832795029L;
2246         return __p.a() + __p.b() * std::tan(__pi * __u);
2247       }
2248
2249   template<typename _RealType>
2250     template<typename _ForwardIterator,
2251              typename _UniformRandomNumberGenerator>
2252       void
2253       cauchy_distribution<_RealType>::
2254       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2255                       _UniformRandomNumberGenerator& __urng,
2256                       const param_type& __p)
2257       {
2258         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2259         const _RealType __pi = 3.1415926535897932384626433832795029L;
2260         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2261           __aurng(__urng);
2262         while (__f != __t)
2263           {
2264             _RealType __u;
2265             do
2266               __u = __aurng();
2267             while (__u == 0.5);
2268
2269             *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2270           }
2271       }
2272
2273   template<typename _RealType, typename _CharT, typename _Traits>
2274     std::basic_ostream<_CharT, _Traits>&
2275     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2276                const cauchy_distribution<_RealType>& __x)
2277     {
2278       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2279       typedef typename __ostream_type::ios_base    __ios_base;
2280
2281       const typename __ios_base::fmtflags __flags = __os.flags();
2282       const _CharT __fill = __os.fill();
2283       const std::streamsize __precision = __os.precision();
2284       const _CharT __space = __os.widen(' ');
2285       __os.flags(__ios_base::scientific | __ios_base::left);
2286       __os.fill(__space);
2287       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2288
2289       __os << __x.a() << __space << __x.b();
2290
2291       __os.flags(__flags);
2292       __os.fill(__fill);
2293       __os.precision(__precision);
2294       return __os;
2295     }
2296
2297   template<typename _RealType, typename _CharT, typename _Traits>
2298     std::basic_istream<_CharT, _Traits>&
2299     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2300                cauchy_distribution<_RealType>& __x)
2301     {
2302       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2303       typedef typename __istream_type::ios_base    __ios_base;
2304
2305       const typename __ios_base::fmtflags __flags = __is.flags();
2306       __is.flags(__ios_base::dec | __ios_base::skipws);
2307
2308       _RealType __a, __b;
2309       __is >> __a >> __b;
2310       __x.param(typename cauchy_distribution<_RealType>::
2311                 param_type(__a, __b));
2312
2313       __is.flags(__flags);
2314       return __is;
2315     }
2316
2317
2318   template<typename _RealType>
2319     template<typename _ForwardIterator,
2320              typename _UniformRandomNumberGenerator>
2321       void
2322       std::fisher_f_distribution<_RealType>::
2323       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2324                       _UniformRandomNumberGenerator& __urng)
2325       {
2326         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2327         while (__f != __t)
2328           *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2329       }
2330
2331   template<typename _RealType>
2332     template<typename _ForwardIterator,
2333              typename _UniformRandomNumberGenerator>
2334       void
2335       std::fisher_f_distribution<_RealType>::
2336       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2337                       _UniformRandomNumberGenerator& __urng,
2338                       const param_type& __p)
2339       {
2340         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2341         typedef typename std::gamma_distribution<result_type>::param_type
2342           param_type;
2343         param_type __p1(__p.m() / 2);
2344         param_type __p2(__p.n() / 2);
2345         while (__f != __t)
2346           *__f++ = ((_M_gd_x(__urng, __p1) * n())
2347                     / (_M_gd_y(__urng, __p2) * m()));
2348       }
2349
2350   template<typename _RealType, typename _CharT, typename _Traits>
2351     std::basic_ostream<_CharT, _Traits>&
2352     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2353                const fisher_f_distribution<_RealType>& __x)
2354     {
2355       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2356       typedef typename __ostream_type::ios_base    __ios_base;
2357
2358       const typename __ios_base::fmtflags __flags = __os.flags();
2359       const _CharT __fill = __os.fill();
2360       const std::streamsize __precision = __os.precision();
2361       const _CharT __space = __os.widen(' ');
2362       __os.flags(__ios_base::scientific | __ios_base::left);
2363       __os.fill(__space);
2364       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2365
2366       __os << __x.m() << __space << __x.n()
2367            << __space << __x._M_gd_x << __space << __x._M_gd_y;
2368
2369       __os.flags(__flags);
2370       __os.fill(__fill);
2371       __os.precision(__precision);
2372       return __os;
2373     }
2374
2375   template<typename _RealType, typename _CharT, typename _Traits>
2376     std::basic_istream<_CharT, _Traits>&
2377     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2378                fisher_f_distribution<_RealType>& __x)
2379     {
2380       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2381       typedef typename __istream_type::ios_base    __ios_base;
2382
2383       const typename __ios_base::fmtflags __flags = __is.flags();
2384       __is.flags(__ios_base::dec | __ios_base::skipws);
2385
2386       _RealType __m, __n;
2387       __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
2388       __x.param(typename fisher_f_distribution<_RealType>::
2389                 param_type(__m, __n));
2390
2391       __is.flags(__flags);
2392       return __is;
2393     }
2394
2395
2396   template<typename _RealType>
2397     template<typename _ForwardIterator,
2398              typename _UniformRandomNumberGenerator>
2399       void
2400       std::student_t_distribution<_RealType>::
2401       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2402                       _UniformRandomNumberGenerator& __urng)
2403       {
2404         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2405         while (__f != __t)
2406           *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2407       }
2408
2409   template<typename _RealType>
2410     template<typename _ForwardIterator,
2411              typename _UniformRandomNumberGenerator>
2412       void
2413       std::student_t_distribution<_RealType>::
2414       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2415                       _UniformRandomNumberGenerator& __urng,
2416                       const param_type& __p)
2417       {
2418         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2419         typename std::gamma_distribution<result_type>::param_type
2420           __p2(__p.n() / 2, 2);
2421         while (__f != __t)
2422           *__f++ =  _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2423       }
2424
2425   template<typename _RealType, typename _CharT, typename _Traits>
2426     std::basic_ostream<_CharT, _Traits>&
2427     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2428                const student_t_distribution<_RealType>& __x)
2429     {
2430       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2431       typedef typename __ostream_type::ios_base    __ios_base;
2432
2433       const typename __ios_base::fmtflags __flags = __os.flags();
2434       const _CharT __fill = __os.fill();
2435       const std::streamsize __precision = __os.precision();
2436       const _CharT __space = __os.widen(' ');
2437       __os.flags(__ios_base::scientific | __ios_base::left);
2438       __os.fill(__space);
2439       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2440
2441       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2442
2443       __os.flags(__flags);
2444       __os.fill(__fill);
2445       __os.precision(__precision);
2446       return __os;
2447     }
2448
2449   template<typename _RealType, typename _CharT, typename _Traits>
2450     std::basic_istream<_CharT, _Traits>&
2451     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2452                student_t_distribution<_RealType>& __x)
2453     {
2454       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2455       typedef typename __istream_type::ios_base    __ios_base;
2456
2457       const typename __ios_base::fmtflags __flags = __is.flags();
2458       __is.flags(__ios_base::dec | __ios_base::skipws);
2459
2460       _RealType __n;
2461       __is >> __n >> __x._M_nd >> __x._M_gd;
2462       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
2463
2464       __is.flags(__flags);
2465       return __is;
2466     }
2467
2468
2469   template<typename _RealType>
2470     void
2471     gamma_distribution<_RealType>::param_type::
2472     _M_initialize()
2473     {
2474       _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2475
2476       const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2477       _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2478     }
2479
2480   /**
2481    * Marsaglia, G. and Tsang, W. W.
2482    * "A Simple Method for Generating Gamma Variables"
2483    * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2484    */
2485   template<typename _RealType>
2486     template<typename _UniformRandomNumberGenerator>
2487       typename gamma_distribution<_RealType>::result_type
2488       gamma_distribution<_RealType>::
2489       operator()(_UniformRandomNumberGenerator& __urng,
2490                  const param_type& __param)
2491       {
2492         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2493           __aurng(__urng);
2494
2495         result_type __u, __v, __n;
2496         const result_type __a1 = (__param._M_malpha
2497                                   - _RealType(1.0) / _RealType(3.0));
2498
2499         do
2500           {
2501             do
2502               {
2503                 __n = _M_nd(__urng);
2504                 __v = result_type(1.0) + __param._M_a2 * __n;
2505               }
2506             while (__v <= 0.0);
2507
2508             __v = __v * __v * __v;
2509             __u = __aurng();
2510           }
2511         while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2512                && (std::log(__u) > (0.5 * __n * __n + __a1
2513                                     * (1.0 - __v + std::log(__v)))));
2514
2515         if (__param.alpha() == __param._M_malpha)
2516           return __a1 * __v * __param.beta();
2517         else
2518           {
2519             do
2520               __u = __aurng();
2521             while (__u == 0.0);
2522
2523             return (std::pow(__u, result_type(1.0) / __param.alpha())
2524                     * __a1 * __v * __param.beta());
2525           }
2526       }
2527
2528   template<typename _RealType>
2529     template<typename _ForwardIterator,
2530              typename _UniformRandomNumberGenerator>
2531       void
2532       gamma_distribution<_RealType>::
2533       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2534                       _UniformRandomNumberGenerator& __urng,
2535                       const param_type& __param)
2536       {
2537         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2538         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2539           __aurng(__urng);
2540
2541         result_type __u, __v, __n;
2542         const result_type __a1 = (__param._M_malpha
2543                                   - _RealType(1.0) / _RealType(3.0));
2544
2545         if (__param.alpha() == __param._M_malpha)
2546           while (__f != __t)
2547             {
2548               do
2549                 {
2550                   do
2551                     {
2552                       __n = _M_nd(__urng);
2553                       __v = result_type(1.0) + __param._M_a2 * __n;
2554                     }
2555                   while (__v <= 0.0);
2556
2557                   __v = __v * __v * __v;
2558                   __u = __aurng();
2559                 }
2560               while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2561                      && (std::log(__u) > (0.5 * __n * __n + __a1
2562                                           * (1.0 - __v + std::log(__v)))));
2563
2564               *__f++ = __a1 * __v * __param.beta();
2565             }
2566         else
2567           while (__f != __t)
2568             {
2569               do
2570                 {
2571                   do
2572                     {
2573                       __n = _M_nd(__urng);
2574                       __v = result_type(1.0) + __param._M_a2 * __n;
2575                     }
2576                   while (__v <= 0.0);
2577
2578                   __v = __v * __v * __v;
2579                   __u = __aurng();
2580                 }
2581               while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
2582                      && (std::log(__u) > (0.5 * __n * __n + __a1
2583                                           * (1.0 - __v + std::log(__v)))));
2584
2585               do
2586                 __u = __aurng();
2587               while (__u == 0.0);
2588
2589               *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2590                         * __a1 * __v * __param.beta());
2591             }
2592       }
2593
2594   template<typename _RealType, typename _CharT, typename _Traits>
2595     std::basic_ostream<_CharT, _Traits>&
2596     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2597                const gamma_distribution<_RealType>& __x)
2598     {
2599       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2600       typedef typename __ostream_type::ios_base    __ios_base;
2601
2602       const typename __ios_base::fmtflags __flags = __os.flags();
2603       const _CharT __fill = __os.fill();
2604       const std::streamsize __precision = __os.precision();
2605       const _CharT __space = __os.widen(' ');
2606       __os.flags(__ios_base::scientific | __ios_base::left);
2607       __os.fill(__space);
2608       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2609
2610       __os << __x.alpha() << __space << __x.beta()
2611            << __space << __x._M_nd;
2612
2613       __os.flags(__flags);
2614       __os.fill(__fill);
2615       __os.precision(__precision);
2616       return __os;
2617     }
2618
2619   template<typename _RealType, typename _CharT, typename _Traits>
2620     std::basic_istream<_CharT, _Traits>&
2621     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2622                gamma_distribution<_RealType>& __x)
2623     {
2624       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2625       typedef typename __istream_type::ios_base    __ios_base;
2626
2627       const typename __ios_base::fmtflags __flags = __is.flags();
2628       __is.flags(__ios_base::dec | __ios_base::skipws);
2629
2630       _RealType __alpha_val, __beta_val;
2631       __is >> __alpha_val >> __beta_val >> __x._M_nd;
2632       __x.param(typename gamma_distribution<_RealType>::
2633                 param_type(__alpha_val, __beta_val));
2634
2635       __is.flags(__flags);
2636       return __is;
2637     }
2638
2639
2640   template<typename _RealType>
2641     template<typename _UniformRandomNumberGenerator>
2642       typename weibull_distribution<_RealType>::result_type
2643       weibull_distribution<_RealType>::
2644       operator()(_UniformRandomNumberGenerator& __urng,
2645                  const param_type& __p)
2646       {
2647         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2648           __aurng(__urng);
2649         return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2650                                   result_type(1) / __p.a());
2651       }
2652
2653   template<typename _RealType>
2654     template<typename _ForwardIterator,
2655              typename _UniformRandomNumberGenerator>
2656       void
2657       weibull_distribution<_RealType>::
2658       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2659                       _UniformRandomNumberGenerator& __urng,
2660                       const param_type& __p)
2661       {
2662         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2663         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2664           __aurng(__urng);
2665         auto __inv_a = result_type(1) / __p.a();
2666
2667         while (__f != __t)
2668           *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2669                                       __inv_a);
2670       }
2671
2672   template<typename _RealType, typename _CharT, typename _Traits>
2673     std::basic_ostream<_CharT, _Traits>&
2674     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2675                const weibull_distribution<_RealType>& __x)
2676     {
2677       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2678       typedef typename __ostream_type::ios_base    __ios_base;
2679
2680       const typename __ios_base::fmtflags __flags = __os.flags();
2681       const _CharT __fill = __os.fill();
2682       const std::streamsize __precision = __os.precision();
2683       const _CharT __space = __os.widen(' ');
2684       __os.flags(__ios_base::scientific | __ios_base::left);
2685       __os.fill(__space);
2686       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2687
2688       __os << __x.a() << __space << __x.b();
2689
2690       __os.flags(__flags);
2691       __os.fill(__fill);
2692       __os.precision(__precision);
2693       return __os;
2694     }
2695
2696   template<typename _RealType, typename _CharT, typename _Traits>
2697     std::basic_istream<_CharT, _Traits>&
2698     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2699                weibull_distribution<_RealType>& __x)
2700     {
2701       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2702       typedef typename __istream_type::ios_base    __ios_base;
2703
2704       const typename __ios_base::fmtflags __flags = __is.flags();
2705       __is.flags(__ios_base::dec | __ios_base::skipws);
2706
2707       _RealType __a, __b;
2708       __is >> __a >> __b;
2709       __x.param(typename weibull_distribution<_RealType>::
2710                 param_type(__a, __b));
2711
2712       __is.flags(__flags);
2713       return __is;
2714     }
2715
2716
2717   template<typename _RealType>
2718     template<typename _UniformRandomNumberGenerator>
2719       typename extreme_value_distribution<_RealType>::result_type
2720       extreme_value_distribution<_RealType>::
2721       operator()(_UniformRandomNumberGenerator& __urng,
2722                  const param_type& __p)
2723       {
2724         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2725           __aurng(__urng);
2726         return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2727                                                       - __aurng()));
2728       }
2729
2730   template<typename _RealType>
2731     template<typename _ForwardIterator,
2732              typename _UniformRandomNumberGenerator>
2733       void
2734       extreme_value_distribution<_RealType>::
2735       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2736                       _UniformRandomNumberGenerator& __urng,
2737                       const param_type& __p)
2738       {
2739         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2740         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2741           __aurng(__urng);
2742
2743         while (__f != __t)
2744           *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2745                                                           - __aurng()));
2746       }
2747
2748   template<typename _RealType, typename _CharT, typename _Traits>
2749     std::basic_ostream<_CharT, _Traits>&
2750     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2751                const extreme_value_distribution<_RealType>& __x)
2752     {
2753       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2754       typedef typename __ostream_type::ios_base    __ios_base;
2755
2756       const typename __ios_base::fmtflags __flags = __os.flags();
2757       const _CharT __fill = __os.fill();
2758       const std::streamsize __precision = __os.precision();
2759       const _CharT __space = __os.widen(' ');
2760       __os.flags(__ios_base::scientific | __ios_base::left);
2761       __os.fill(__space);
2762       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2763
2764       __os << __x.a() << __space << __x.b();
2765
2766       __os.flags(__flags);
2767       __os.fill(__fill);
2768       __os.precision(__precision);
2769       return __os;
2770     }
2771
2772   template<typename _RealType, typename _CharT, typename _Traits>
2773     std::basic_istream<_CharT, _Traits>&
2774     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2775                extreme_value_distribution<_RealType>& __x)
2776     {
2777       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2778       typedef typename __istream_type::ios_base    __ios_base;
2779
2780       const typename __ios_base::fmtflags __flags = __is.flags();
2781       __is.flags(__ios_base::dec | __ios_base::skipws);
2782
2783       _RealType __a, __b;
2784       __is >> __a >> __b;
2785       __x.param(typename extreme_value_distribution<_RealType>::
2786                 param_type(__a, __b));
2787
2788       __is.flags(__flags);
2789       return __is;
2790     }
2791
2792
2793   template<typename _IntType>
2794     void
2795     discrete_distribution<_IntType>::param_type::
2796     _M_initialize()
2797     {
2798       if (_M_prob.size() < 2)
2799         {
2800           _M_prob.clear();
2801           return;
2802         }
2803
2804       const double __sum = std::accumulate(_M_prob.begin(),
2805                                            _M_prob.end(), 0.0);
2806       // Now normalize the probabilites.
2807       __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2808                             __sum);
2809       // Accumulate partial sums.
2810       _M_cp.reserve(_M_prob.size());
2811       std::partial_sum(_M_prob.begin(), _M_prob.end(),
2812                        std::back_inserter(_M_cp));
2813       // Make sure the last cumulative probability is one.
2814       _M_cp[_M_cp.size() - 1] = 1.0;
2815     }
2816
2817   template<typename _IntType>
2818     template<typename _Func>
2819       discrete_distribution<_IntType>::param_type::
2820       param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2821       : _M_prob(), _M_cp()
2822       {
2823         const size_t __n = __nw == 0 ? 1 : __nw;
2824         const double __delta = (__xmax - __xmin) / __n;
2825
2826         _M_prob.reserve(__n);
2827         for (size_t __k = 0; __k < __nw; ++__k)
2828           _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2829
2830         _M_initialize();
2831       }
2832
2833   template<typename _IntType>
2834     template<typename _UniformRandomNumberGenerator>
2835       typename discrete_distribution<_IntType>::result_type
2836       discrete_distribution<_IntType>::
2837       operator()(_UniformRandomNumberGenerator& __urng,
2838                  const param_type& __param)
2839       {
2840         if (__param._M_cp.empty())
2841           return result_type(0);
2842
2843         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2844           __aurng(__urng);
2845
2846         const double __p = __aurng();
2847         auto __pos = std::lower_bound(__param._M_cp.begin(),
2848                                       __param._M_cp.end(), __p);
2849
2850         return __pos - __param._M_cp.begin();
2851       }
2852
2853   template<typename _IntType>
2854     template<typename _ForwardIterator,
2855              typename _UniformRandomNumberGenerator>
2856       void
2857       discrete_distribution<_IntType>::
2858       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2859                       _UniformRandomNumberGenerator& __urng,
2860                       const param_type& __param)
2861       {
2862         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2863
2864         if (__param._M_cp.empty())
2865           {
2866             while (__f != __t)
2867               *__f++ = result_type(0);
2868             return;
2869           }
2870
2871         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2872           __aurng(__urng);
2873
2874         while (__f != __t)
2875           {
2876             const double __p = __aurng();
2877             auto __pos = std::lower_bound(__param._M_cp.begin(),
2878                                           __param._M_cp.end(), __p);
2879
2880             *__f++ = __pos - __param._M_cp.begin();
2881           }
2882       }
2883
2884   template<typename _IntType, typename _CharT, typename _Traits>
2885     std::basic_ostream<_CharT, _Traits>&
2886     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2887                const discrete_distribution<_IntType>& __x)
2888     {
2889       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2890       typedef typename __ostream_type::ios_base    __ios_base;
2891
2892       const typename __ios_base::fmtflags __flags = __os.flags();
2893       const _CharT __fill = __os.fill();
2894       const std::streamsize __precision = __os.precision();
2895       const _CharT __space = __os.widen(' ');
2896       __os.flags(__ios_base::scientific | __ios_base::left);
2897       __os.fill(__space);
2898       __os.precision(std::numeric_limits<double>::max_digits10);
2899
2900       std::vector<double> __prob = __x.probabilities();
2901       __os << __prob.size();
2902       for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2903         __os << __space << *__dit;
2904
2905       __os.flags(__flags);
2906       __os.fill(__fill);
2907       __os.precision(__precision);
2908       return __os;
2909     }
2910
2911   template<typename _IntType, typename _CharT, typename _Traits>
2912     std::basic_istream<_CharT, _Traits>&
2913     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2914                discrete_distribution<_IntType>& __x)
2915     {
2916       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2917       typedef typename __istream_type::ios_base    __ios_base;
2918
2919       const typename __ios_base::fmtflags __flags = __is.flags();
2920       __is.flags(__ios_base::dec | __ios_base::skipws);
2921
2922       size_t __n;
2923       __is >> __n;
2924
2925       std::vector<double> __prob_vec;
2926       __prob_vec.reserve(__n);
2927       for (; __n != 0; --__n)
2928         {
2929           double __prob;
2930           __is >> __prob;
2931           __prob_vec.push_back(__prob);
2932         }
2933
2934       __x.param(typename discrete_distribution<_IntType>::
2935                 param_type(__prob_vec.begin(), __prob_vec.end()));
2936
2937       __is.flags(__flags);
2938       return __is;
2939     }
2940
2941
2942   template<typename _RealType>
2943     void
2944     piecewise_constant_distribution<_RealType>::param_type::
2945     _M_initialize()
2946     {
2947       if (_M_int.size() < 2
2948           || (_M_int.size() == 2
2949               && _M_int[0] == _RealType(0)
2950               && _M_int[1] == _RealType(1)))
2951         {
2952           _M_int.clear();
2953           _M_den.clear();
2954           return;
2955         }
2956
2957       const double __sum = std::accumulate(_M_den.begin(),
2958                                            _M_den.end(), 0.0);
2959
2960       __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
2961                             __sum);
2962
2963       _M_cp.reserve(_M_den.size());
2964       std::partial_sum(_M_den.begin(), _M_den.end(),
2965                        std::back_inserter(_M_cp));
2966
2967       // Make sure the last cumulative probability is one.
2968       _M_cp[_M_cp.size() - 1] = 1.0;
2969
2970       for (size_t __k = 0; __k < _M_den.size(); ++__k)
2971         _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2972     }
2973
2974   template<typename _RealType>
2975     template<typename _InputIteratorB, typename _InputIteratorW>
2976       piecewise_constant_distribution<_RealType>::param_type::
2977       param_type(_InputIteratorB __bbegin,
2978                  _InputIteratorB __bend,
2979                  _InputIteratorW __wbegin)
2980       : _M_int(), _M_den(), _M_cp()
2981       {
2982         if (__bbegin != __bend)
2983           {
2984             for (;;)
2985               {
2986                 _M_int.push_back(*__bbegin);
2987                 ++__bbegin;
2988                 if (__bbegin == __bend)
2989                   break;
2990
2991                 _M_den.push_back(*__wbegin);
2992                 ++__wbegin;
2993               }
2994           }
2995
2996         _M_initialize();
2997       }
2998
2999   template<typename _RealType>
3000     template<typename _Func>
3001       piecewise_constant_distribution<_RealType>::param_type::
3002       param_type(initializer_list<_RealType> __bl, _Func __fw)
3003       : _M_int(), _M_den(), _M_cp()
3004       {
3005         _M_int.reserve(__bl.size());
3006         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3007           _M_int.push_back(*__biter);
3008
3009         _M_den.reserve(_M_int.size() - 1);
3010         for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3011           _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
3012
3013         _M_initialize();
3014       }
3015
3016   template<typename _RealType>
3017     template<typename _Func>
3018       piecewise_constant_distribution<_RealType>::param_type::
3019       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3020       : _M_int(), _M_den(), _M_cp()
3021       {
3022         const size_t __n = __nw == 0 ? 1 : __nw;
3023         const _RealType __delta = (__xmax - __xmin) / __n;
3024
3025         _M_int.reserve(__n + 1);
3026         for (size_t __k = 0; __k <= __nw; ++__k)
3027           _M_int.push_back(__xmin + __k * __delta);
3028
3029         _M_den.reserve(__n);
3030         for (size_t __k = 0; __k < __nw; ++__k)
3031           _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
3032
3033         _M_initialize();
3034       }
3035
3036   template<typename _RealType>
3037     template<typename _UniformRandomNumberGenerator>
3038       typename piecewise_constant_distribution<_RealType>::result_type
3039       piecewise_constant_distribution<_RealType>::
3040       operator()(_UniformRandomNumberGenerator& __urng,
3041                  const param_type& __param)
3042       {
3043         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3044           __aurng(__urng);
3045
3046         const double __p = __aurng();
3047         if (__param._M_cp.empty())
3048           return __p;
3049
3050         auto __pos = std::lower_bound(__param._M_cp.begin(),
3051                                       __param._M_cp.end(), __p);
3052         const size_t __i = __pos - __param._M_cp.begin();
3053
3054         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3055
3056         return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
3057       }
3058
3059   template<typename _RealType>
3060     template<typename _ForwardIterator,
3061              typename _UniformRandomNumberGenerator>
3062       void
3063       piecewise_constant_distribution<_RealType>::
3064       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3065                       _UniformRandomNumberGenerator& __urng,
3066                       const param_type& __param)
3067       {
3068         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3069         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3070           __aurng(__urng);
3071
3072         if (__param._M_cp.empty())
3073           {
3074             while (__f != __t)
3075               *__f++ = __aurng();
3076             return;
3077           }
3078
3079         while (__f != __t)
3080           {
3081             const double __p = __aurng();
3082
3083             auto __pos = std::lower_bound(__param._M_cp.begin(),
3084                                           __param._M_cp.end(), __p);
3085             const size_t __i = __pos - __param._M_cp.begin();
3086
3087             const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3088
3089             *__f++ = (__param._M_int[__i]
3090                       + (__p - __pref) / __param._M_den[__i]);
3091           }
3092       }
3093
3094   template<typename _RealType, typename _CharT, typename _Traits>
3095     std::basic_ostream<_CharT, _Traits>&
3096     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3097                const piecewise_constant_distribution<_RealType>& __x)
3098     {
3099       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
3100       typedef typename __ostream_type::ios_base    __ios_base;
3101
3102       const typename __ios_base::fmtflags __flags = __os.flags();
3103       const _CharT __fill = __os.fill();
3104       const std::streamsize __precision = __os.precision();
3105       const _CharT __space = __os.widen(' ');
3106       __os.flags(__ios_base::scientific | __ios_base::left);
3107       __os.fill(__space);
3108       __os.precision(std::numeric_limits<_RealType>::max_digits10);
3109
3110       std::vector<_RealType> __int = __x.intervals();
3111       __os << __int.size() - 1;
3112
3113       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3114         __os << __space << *__xit;
3115
3116       std::vector<double> __den = __x.densities();
3117       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3118         __os << __space << *__dit;
3119
3120       __os.flags(__flags);
3121       __os.fill(__fill);
3122       __os.precision(__precision);
3123       return __os;
3124     }
3125
3126   template<typename _RealType, typename _CharT, typename _Traits>
3127     std::basic_istream<_CharT, _Traits>&
3128     operator>>(std::basic_istream<_CharT, _Traits>& __is,
3129                piecewise_constant_distribution<_RealType>& __x)
3130     {
3131       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
3132       typedef typename __istream_type::ios_base    __ios_base;
3133
3134       const typename __ios_base::fmtflags __flags = __is.flags();
3135       __is.flags(__ios_base::dec | __ios_base::skipws);
3136
3137       size_t __n;
3138       __is >> __n;
3139
3140       std::vector<_RealType> __int_vec;
3141       __int_vec.reserve(__n + 1);
3142       for (size_t __i = 0; __i <= __n; ++__i)
3143         {
3144           _RealType __int;
3145           __is >> __int;
3146           __int_vec.push_back(__int);
3147         }
3148
3149       std::vector<double> __den_vec;
3150       __den_vec.reserve(__n);
3151       for (size_t __i = 0; __i < __n; ++__i)
3152         {
3153           double __den;
3154           __is >> __den;
3155           __den_vec.push_back(__den);
3156         }
3157
3158       __x.param(typename piecewise_constant_distribution<_RealType>::
3159           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3160
3161       __is.flags(__flags);
3162       return __is;
3163     }
3164
3165
3166   template<typename _RealType>
3167     void
3168     piecewise_linear_distribution<_RealType>::param_type::
3169     _M_initialize()
3170     {
3171       if (_M_int.size() < 2
3172           || (_M_int.size() == 2
3173               && _M_int[0] == _RealType(0)
3174               && _M_int[1] == _RealType(1)
3175               && _M_den[0] == _M_den[1]))
3176         {
3177           _M_int.clear();
3178           _M_den.clear();
3179           return;
3180         }
3181
3182       double __sum = 0.0;
3183       _M_cp.reserve(_M_int.size() - 1);
3184       _M_m.reserve(_M_int.size() - 1);
3185       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3186         {
3187           const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3188           __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3189           _M_cp.push_back(__sum);
3190           _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3191         }
3192
3193       //  Now normalize the densities...
3194       __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
3195                             __sum);
3196       //  ... and partial sums...
3197       __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
3198       //  ... and slopes.
3199       __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
3200
3201       //  Make sure the last cumulative probablility is one.
3202       _M_cp[_M_cp.size() - 1] = 1.0;
3203      }
3204
3205   template<typename _RealType>
3206     template<typename _InputIteratorB, typename _InputIteratorW>
3207       piecewise_linear_distribution<_RealType>::param_type::
3208       param_type(_InputIteratorB __bbegin,
3209                  _InputIteratorB __bend,
3210                  _InputIteratorW __wbegin)
3211       : _M_int(), _M_den(), _M_cp(), _M_m()
3212       {
3213         for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
3214           {
3215             _M_int.push_back(*__bbegin);
3216             _M_den.push_back(*__wbegin);
3217           }
3218
3219         _M_initialize();
3220       }
3221
3222   template<typename _RealType>
3223     template<typename _Func>
3224       piecewise_linear_distribution<_RealType>::param_type::
3225       param_type(initializer_list<_RealType> __bl, _Func __fw)
3226       : _M_int(), _M_den(), _M_cp(), _M_m()
3227       {
3228         _M_int.reserve(__bl.size());
3229         _M_den.reserve(__bl.size());
3230         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3231           {
3232             _M_int.push_back(*__biter);
3233             _M_den.push_back(__fw(*__biter));
3234           }
3235
3236         _M_initialize();
3237       }
3238
3239   template<typename _RealType>
3240     template<typename _Func>
3241       piecewise_linear_distribution<_RealType>::param_type::
3242       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3243       : _M_int(), _M_den(), _M_cp(), _M_m()
3244       {
3245         const size_t __n = __nw == 0 ? 1 : __nw;
3246         const _RealType __delta = (__xmax - __xmin) / __n;
3247
3248         _M_int.reserve(__n + 1);
3249         _M_den.reserve(__n + 1);
3250         for (size_t __k = 0; __k <= __nw; ++__k)
3251           {
3252             _M_int.push_back(__xmin + __k * __delta);
3253             _M_den.push_back(__fw(_M_int[__k] + __delta));
3254           }
3255
3256         _M_initialize();
3257       }
3258
3259   template<typename _RealType>
3260     template<typename _UniformRandomNumberGenerator>
3261       typename piecewise_linear_distribution<_RealType>::result_type
3262       piecewise_linear_distribution<_RealType>::
3263       operator()(_UniformRandomNumberGenerator& __urng,
3264                  const param_type& __param)
3265       {
3266         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
3267           __aurng(__urng);
3268
3269         const double __p = __aurng();
3270         if (__param._M_cp.empty())
3271           return __p;
3272
3273         auto __pos = std::lower_bound(__param._M_cp.begin(),
3274                                       __param._M_cp.end(), __p);
3275         const size_t __i = __pos - __param._M_cp.begin();
3276
3277         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3278
3279         const double __a = 0.5 * __param._M_m[__i];
3280         const double __b = __param._M_den[__i];
3281         const double __cm = __p - __pref;
3282
3283         _RealType __x = __param._M_int[__i];
3284         if (__a == 0)
3285           __x += __cm / __b;
3286         else
3287           {
3288             const double __d = __b * __b + 4.0 * __a * __cm;
3289             __x += 0.5 * (std::sqrt(__d) - __b) / __a;
3290           }
3291
3292         return __x;
3293       }
3294
3295   template<typename _RealType>
3296     template<typename _ForwardIterator,
3297              typename _UniformRandomNumberGenerator>
3298       void
3299       piecewise_linear_distribution<_RealType>::
3300       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3301                       _UniformRandomNumberGenerator& __urng,
3302                       const param_type& __param)
3303       {
3304         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3305         // We could duplicate everything from operator()...
3306         while (__f != __t)
3307           *__f++ = this->operator()(__urng, __param);
3308       }
3309
3310   template<typename _RealType, typename _CharT, typename _Traits>
3311     std::basic_ostream<_CharT, _Traits>&
3312     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3313                const piecewise_linear_distribution<_RealType>& __x)
3314     {
3315       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
3316       typedef typename __ostream_type::ios_base    __ios_base;
3317
3318       const typename __ios_base::fmtflags __flags = __os.flags();
3319       const _CharT __fill = __os.fill();
3320       const std::streamsize __precision = __os.precision();
3321       const _CharT __space = __os.widen(' ');
3322       __os.flags(__ios_base::scientific | __ios_base::left);
3323       __os.fill(__space);
3324       __os.precision(std::numeric_limits<_RealType>::max_digits10);
3325
3326       std::vector<_RealType> __int = __x.intervals();
3327       __os << __int.size() - 1;
3328
3329       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3330         __os << __space << *__xit;
3331
3332       std::vector<double> __den = __x.densities();
3333       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3334         __os << __space << *__dit;
3335
3336       __os.flags(__flags);
3337       __os.fill(__fill);
3338       __os.precision(__precision);
3339       return __os;
3340     }
3341
3342   template<typename _RealType, typename _CharT, typename _Traits>
3343     std::basic_istream<_CharT, _Traits>&
3344     operator>>(std::basic_istream<_CharT, _Traits>& __is,
3345                piecewise_linear_distribution<_RealType>& __x)
3346     {
3347       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
3348       typedef typename __istream_type::ios_base    __ios_base;
3349
3350       const typename __ios_base::fmtflags __flags = __is.flags();
3351       __is.flags(__ios_base::dec | __ios_base::skipws);
3352
3353       size_t __n;
3354       __is >> __n;
3355
3356       std::vector<_RealType> __int_vec;
3357       __int_vec.reserve(__n + 1);
3358       for (size_t __i = 0; __i <= __n; ++__i)
3359         {
3360           _RealType __int;
3361           __is >> __int;
3362           __int_vec.push_back(__int);
3363         }
3364
3365       std::vector<double> __den_vec;
3366       __den_vec.reserve(__n + 1);
3367       for (size_t __i = 0; __i <= __n; ++__i)
3368         {
3369           double __den;
3370           __is >> __den;
3371           __den_vec.push_back(__den);
3372         }
3373
3374       __x.param(typename piecewise_linear_distribution<_RealType>::
3375           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
3376
3377       __is.flags(__flags);
3378       return __is;
3379     }
3380
3381
3382   template<typename _IntType>
3383     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3384     {
3385       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3386         _M_v.push_back(__detail::__mod<result_type,
3387                        __detail::_Shift<result_type, 32>::__value>(*__iter));
3388     }
3389
3390   template<typename _InputIterator>
3391     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3392     {
3393       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3394         _M_v.push_back(__detail::__mod<result_type,
3395                        __detail::_Shift<result_type, 32>::__value>(*__iter));
3396     }
3397
3398   template<typename _RandomAccessIterator>
3399     void
3400     seed_seq::generate(_RandomAccessIterator __begin,
3401                        _RandomAccessIterator __end)
3402     {
3403       typedef typename iterator_traits<_RandomAccessIterator>::value_type
3404         _Type;
3405
3406       if (__begin == __end)
3407         return;
3408
3409       std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3410
3411       const size_t __n = __end - __begin;
3412       const size_t __s = _M_v.size();
3413       const size_t __t = (__n >= 623) ? 11
3414                        : (__n >=  68) ? 7
3415                        : (__n >=  39) ? 5
3416                        : (__n >=   7) ? 3
3417                        : (__n - 1) / 2;
3418       const size_t __p = (__n - __t) / 2;
3419       const size_t __q = __p + __t;
3420       const size_t __m = std::max(size_t(__s + 1), __n);
3421
3422       for (size_t __k = 0; __k < __m; ++__k)
3423         {
3424           _Type __arg = (__begin[__k % __n]
3425                          ^ __begin[(__k + __p) % __n]
3426                          ^ __begin[(__k - 1) % __n]);
3427           _Type __r1 = __arg ^ (__arg >> 27);
3428           __r1 = __detail::__mod<_Type,
3429                     __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
3430           _Type __r2 = __r1;
3431           if (__k == 0)
3432             __r2 += __s;
3433           else if (__k <= __s)
3434             __r2 += __k % __n + _M_v[__k - 1];
3435           else
3436             __r2 += __k % __n;
3437           __r2 = __detail::__mod<_Type,
3438                    __detail::_Shift<_Type, 32>::__value>(__r2);
3439           __begin[(__k + __p) % __n] += __r1;
3440           __begin[(__k + __q) % __n] += __r2;
3441           __begin[__k % __n] = __r2;
3442         }
3443
3444       for (size_t __k = __m; __k < __m + __n; ++__k)
3445         {
3446           _Type __arg = (__begin[__k % __n]
3447                          + __begin[(__k + __p) % __n]
3448                          + __begin[(__k - 1) % __n]);
3449           _Type __r3 = __arg ^ (__arg >> 27);
3450           __r3 = __detail::__mod<_Type,
3451                    __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
3452           _Type __r4 = __r3 - __k % __n;
3453           __r4 = __detail::__mod<_Type,
3454                    __detail::_Shift<_Type, 32>::__value>(__r4);
3455           __begin[(__k + __p) % __n] ^= __r3;
3456           __begin[(__k + __q) % __n] ^= __r4;
3457           __begin[__k % __n] = __r4;
3458         }
3459     }
3460
3461   template<typename _RealType, size_t __bits,
3462            typename _UniformRandomNumberGenerator>
3463     _RealType
3464     generate_canonical(_UniformRandomNumberGenerator& __urng)
3465     {
3466       static_assert(std::is_floating_point<_RealType>::value,
3467                     "template argument not a floating point type");
3468
3469       const size_t __b
3470         = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3471                    __bits);
3472       const long double __r = static_cast<long double>(__urng.max())
3473                             - static_cast<long double>(__urng.min()) + 1.0L;
3474       const size_t __log2r = std::log(__r) / std::log(2.0L);
3475       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
3476       _RealType __sum = _RealType(0);
3477       _RealType __tmp = _RealType(1);
3478       for (; __k != 0; --__k)
3479         {
3480           __sum += _RealType(__urng() - __urng.min()) * __tmp;
3481           __tmp *= __r;
3482         }
3483       return __sum / __tmp;
3484     }
3485
3486 _GLIBCXX_END_NAMESPACE_VERSION
3487 } // namespace
3488
3489 #endif