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