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

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