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

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