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>
 130       typename std::enable_if<std::is_class<_Sseq>::value>::type
 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>
 348       typename std::enable_if<std::is_class<_Sseq>::value>::type
 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>
 535       typename std::enable_if<std::is_class<_Sseq>::value>::type
 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 std::make_unsigned<typename
 826           _UniformRandomNumberGenerator::result_type>::type __urntype;
 827         typedef typename std::make_unsigned<result_type>::type __utype;
 828         typedef typename std::conditional<(sizeof(__urntype) > sizeof(__utype)),
 829           __urntype, __utype>::type __uctype;
 830
 831         result_type __ret;
 832
 833         const __urntype __urnmin = __urng.min();
 834         const __urntype __urnmax = __urng.max();
 835         const __urntype __urnrange = __urnmax - __urnmin;
 836         const __uctype __urange = __param.b() - __param.a();
 837         const __uctype __udenom = (__urnrange <= __urange
 838                                    ? 1 : __urnrange / (__urange + 1));
 839         do
 840           __ret = (__urntype(__urng()) -  __urnmin) / __udenom;
 841         while (__ret > __param.b() - __param.a());
 842
 843         return __ret + __param.a();
 844       }
 845
 846   template<typename _IntType, typename _CharT, typename _Traits>
 847     std::basic_ostream<_CharT, _Traits>&
 848     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 849                const uniform_int_distribution<_IntType>& __x)
 850     {
 851       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 852       typedef typename __ostream_type::ios_base    __ios_base;
 853
 854       const typename __ios_base::fmtflags __flags = __os.flags();
 855       const _CharT __fill = __os.fill();
 856       const _CharT __space = __os.widen(' ');
 857       __os.flags(__ios_base::scientific | __ios_base::left);
 858       __os.fill(__space);
 859
 860       __os << __x.a() << __space << __x.b();
 861
 862       __os.flags(__flags);
 863       __os.fill(__fill);
 864       return __os;
 865     }
 866
 867   template<typename _IntType, typename _CharT, typename _Traits>
 868     std::basic_istream<_CharT, _Traits>&
 869     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 870                uniform_int_distribution<_IntType>& __x)
 871     {
 872       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 873       typedef typename __istream_type::ios_base    __ios_base;
 874
 875       const typename __ios_base::fmtflags __flags = __is.flags();
 876       __is.flags(__ios_base::dec | __ios_base::skipws);
 877
 878       _IntType __a, __b;
 879       __is >> __a >> __b;
 880       __x.param(typename uniform_int_distribution<_IntType>::
 881                 param_type(__a, __b));
 882
 883       __is.flags(__flags);
 884       return __is;
 885     }
 886
 887
 888   template<typename _RealType, typename _CharT, typename _Traits>
 889     std::basic_ostream<_CharT, _Traits>&
 890     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 891                const uniform_real_distribution<_RealType>& __x)
 892     {
 893       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 894       typedef typename __ostream_type::ios_base    __ios_base;
 895
 896       const typename __ios_base::fmtflags __flags = __os.flags();
 897       const _CharT __fill = __os.fill();
 898       const std::streamsize __precision = __os.precision();
 899       const _CharT __space = __os.widen(' ');
 900       __os.flags(__ios_base::scientific | __ios_base::left);
 901       __os.fill(__space);
 902       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 903
 904       __os << __x.a() << __space << __x.b();
 905
 906       __os.flags(__flags);
 907       __os.fill(__fill);
 908       __os.precision(__precision);
 909       return __os;
 910     }
 911
 912   template<typename _RealType, typename _CharT, typename _Traits>
 913     std::basic_istream<_CharT, _Traits>&
 914     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 915                uniform_real_distribution<_RealType>& __x)
 916     {
 917       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 918       typedef typename __istream_type::ios_base    __ios_base;
 919
 920       const typename __ios_base::fmtflags __flags = __is.flags();
 921       __is.flags(__ios_base::skipws);
 922
 923       _RealType __a, __b;
 924       __is >> __a >> __b;
 925       __x.param(typename uniform_real_distribution<_RealType>::
 926                 param_type(__a, __b));
 927
 928       __is.flags(__flags);
 929       return __is;
 930     }
 931
 932
 933   template<typename _CharT, typename _Traits>
 934     std::basic_ostream<_CharT, _Traits>&
 935     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 936                const bernoulli_distribution& __x)
 937     {
 938       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 939       typedef typename __ostream_type::ios_base    __ios_base;
 940
 941       const typename __ios_base::fmtflags __flags = __os.flags();
 942       const _CharT __fill = __os.fill();
 943       const std::streamsize __precision = __os.precision();
 944       __os.flags(__ios_base::scientific | __ios_base::left);
 945       __os.fill(__os.widen(' '));
 946       __os.precision(std::numeric_limits<double>::max_digits10);
 947
 948       __os << __x.p();
 949
 950       __os.flags(__flags);
 951       __os.fill(__fill);
 952       __os.precision(__precision);
 953       return __os;
 954     }
 955
 956
 957   template<typename _IntType>
 958     template<typename _UniformRandomNumberGenerator>
 959       typename geometric_distribution<_IntType>::result_type
 960       geometric_distribution<_IntType>::
 961       operator()(_UniformRandomNumberGenerator& __urng,
 962                  const param_type& __param)
 963       {
 964         // About the epsilon thing see this thread:
 965         // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
 966         const double __naf =
 967           (1 - std::numeric_limits<double>::epsilon()) / 2;
 968         // The largest _RealType convertible to _IntType.
 969         const double __thr =
 970           std::numeric_limits<_IntType>::max() + __naf;
 971         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
 972           __aurng(__urng);
 973
 974         double __cand;
 975         do
 976           __cand = std::ceil(std::log(__aurng()) / __param._M_log_p);
 977         while (__cand >= __thr);
 978
 979         return result_type(__cand + __naf);
 980       }
 981
 982   template<typename _IntType,
 983            typename _CharT, typename _Traits>
 984     std::basic_ostream<_CharT, _Traits>&
 985     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 986                const geometric_distribution<_IntType>& __x)
 987     {
 988       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 989       typedef typename __ostream_type::ios_base    __ios_base;
 990
 991       const typename __ios_base::fmtflags __flags = __os.flags();
 992       const _CharT __fill = __os.fill();
 993       const std::streamsize __precision = __os.precision();
 994       __os.flags(__ios_base::scientific | __ios_base::left);
 995       __os.fill(__os.widen(' '));
 996       __os.precision(std::numeric_limits<double>::max_digits10);
 997
 998       __os << __x.p();
 999
1000       __os.flags(__flags);
1001       __os.fill(__fill);
1002       __os.precision(__precision);
1003       return __os;
1004     }
1005
1006   template<typename _IntType,
1007            typename _CharT, typename _Traits>
1008     std::basic_istream<_CharT, _Traits>&
1009     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1010                geometric_distribution<_IntType>& __x)
1011     {
1012       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1013       typedef typename __istream_type::ios_base    __ios_base;
1014
1015       const typename __ios_base::fmtflags __flags = __is.flags();
1016       __is.flags(__ios_base::skipws);
1017
1018       double __p;
1019       __is >> __p;
1020       __x.param(typename geometric_distribution<_IntType>::param_type(__p));
1021
1022       __is.flags(__flags);
1023       return __is;
1024     }
1025
1026
1027   template<typename _IntType>
1028     template<typename _UniformRandomNumberGenerator>
1029       typename negative_binomial_distribution<_IntType>::result_type
1030       negative_binomial_distribution<_IntType>::
1031       operator()(_UniformRandomNumberGenerator& __urng)
1032       {
1033         const double __y = _M_gd(__urng);
1034
1035         // XXX Is the constructor too slow?
1036         std::poisson_distribution<result_type> __poisson(__y);
1037         return __poisson(__urng);
1038       }
1039
1040   template<typename _IntType>
1041     template<typename _UniformRandomNumberGenerator>
1042       typename negative_binomial_distribution<_IntType>::result_type
1043       negative_binomial_distribution<_IntType>::
1044       operator()(_UniformRandomNumberGenerator& __urng,
1045                  const param_type& __p)
1046       {
1047         typedef typename std::gamma_distribution<result_type>::param_type
1048           param_type;
1049
1050         const double __y =
1051           _M_gd(__urng, param_type(__p.k(), __p.p() / (1.0 - __p.p())));
1052
1053         std::poisson_distribution<result_type> __poisson(__y);
1054         return __poisson(__urng);
1055       }
1056
1057   template<typename _IntType, typename _CharT, typename _Traits>
1058     std::basic_ostream<_CharT, _Traits>&
1059     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1060                const negative_binomial_distribution<_IntType>& __x)
1061     {
1062       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1063       typedef typename __ostream_type::ios_base    __ios_base;
1064
1065       const typename __ios_base::fmtflags __flags = __os.flags();
1066       const _CharT __fill = __os.fill();
1067       const std::streamsize __precision = __os.precision();
1068       const _CharT __space = __os.widen(' ');
1069       __os.flags(__ios_base::scientific | __ios_base::left);
1070       __os.fill(__os.widen(' '));
1071       __os.precision(std::numeric_limits<double>::max_digits10);
1072
1073       __os << __x.k() << __space << __x.p()
1074            << __space << __x._M_gd;
1075
1076       __os.flags(__flags);
1077       __os.fill(__fill);
1078       __os.precision(__precision);
1079       return __os;
1080     }
1081
1082   template<typename _IntType, typename _CharT, typename _Traits>
1083     std::basic_istream<_CharT, _Traits>&
1084     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1085                negative_binomial_distribution<_IntType>& __x)
1086     {
1087       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1088       typedef typename __istream_type::ios_base    __ios_base;
1089
1090       const typename __ios_base::fmtflags __flags = __is.flags();
1091       __is.flags(__ios_base::skipws);
1092
1093       _IntType __k;
1094       double __p;
1095       __is >> __k >> __p >> __x._M_gd;
1096       __x.param(typename negative_binomial_distribution<_IntType>::
1097                 param_type(__k, __p));
1098
1099       __is.flags(__flags);
1100       return __is;
1101     }
1102
1103
1104   template<typename _IntType>
1105     void
1106     poisson_distribution<_IntType>::param_type::
1107     _M_initialize()
1108     {
1109 #if _GLIBCXX_USE_C99_MATH_TR1
1110       if (_M_mean >= 12)
1111         {
1112           const double __m = std::floor(_M_mean);
1113           _M_lm_thr = std::log(_M_mean);
1114           _M_lfm = std::lgamma(__m + 1);
1115           _M_sm = std::sqrt(__m);
1116
1117           const double __pi_4 = 0.7853981633974483096156608458198757L;
1118           const double __dx = std::sqrt(2 * __m * std::log(32 * __m
1119                                                               / __pi_4));
1120           _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
1121           const double __cx = 2 * __m + _M_d;
1122           _M_scx = std::sqrt(__cx / 2);
1123           _M_1cx = 1 / __cx;
1124
1125           _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1126           _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
1127                 / _M_d;
1128         }
1129       else
1130 #endif
1131         _M_lm_thr = std::exp(-_M_mean);
1132       }
1133
1134   /**
1135    * A rejection algorithm when mean >= 12 and a simple method based
1136    * upon the multiplication of uniform random variates otherwise.
1137    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1138    * is defined.
1139    *
1140    * Reference:
1141    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1142    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1143    */
1144   template<typename _IntType>
1145     template<typename _UniformRandomNumberGenerator>
1146       typename poisson_distribution<_IntType>::result_type
1147       poisson_distribution<_IntType>::
1148       operator()(_UniformRandomNumberGenerator& __urng,
1149                  const param_type& __param)
1150       {
1151         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1152           __aurng(__urng);
1153 #if _GLIBCXX_USE_C99_MATH_TR1
1154         if (__param.mean() >= 12)
1155           {
1156             double __x;
1157
1158             // See comments above...
1159             const double __naf =
1160               (1 - std::numeric_limits<double>::epsilon()) / 2;
1161             const double __thr =
1162               std::numeric_limits<_IntType>::max() + __naf;
1163
1164             const double __m = std::floor(__param.mean());
1165             // sqrt(pi / 2)
1166             const double __spi_2 = 1.2533141373155002512078826424055226L;
1167             const double __c1 = __param._M_sm * __spi_2;
1168             const double __c2 = __param._M_c2b + __c1;
1169             const double __c3 = __c2 + 1;
1170             const double __c4 = __c3 + 1;
1171             // e^(1 / 78)
1172             const double __e178 = 1.0129030479320018583185514777512983L;
1173             const double __c5 = __c4 + __e178;
1174             const double __c = __param._M_cb + __c5;
1175             const double __2cx = 2 * (2 * __m + __param._M_d);
1176
1177             bool __reject = true;
1178             do
1179               {
1180                 const double __u = __c * __aurng();
1181                 const double __e = -std::log(__aurng());
1182
1183                 double __w = 0.0;
1184
1185                 if (__u <= __c1)
1186                   {
1187                     const double __n = _M_nd(__urng);
1188                     const double __y = -std::abs(__n) * __param._M_sm - 1;
1189                     __x = std::floor(__y);
1190                     __w = -__n * __n / 2;
1191                     if (__x < -__m)
1192                       continue;
1193                   }
1194                 else if (__u <= __c2)
1195                   {
1196                     const double __n = _M_nd(__urng);
1197                     const double __y = 1 + std::abs(__n) * __param._M_scx;
1198                     __x = std::ceil(__y);
1199                     __w = __y * (2 - __y) * __param._M_1cx;
1200                     if (__x > __param._M_d)
1201                       continue;
1202                   }
1203                 else if (__u <= __c3)
1204                   // NB: This case not in the book, nor in the Errata,
1205                   // but should be ok...
1206                   __x = -1;
1207                 else if (__u <= __c4)
1208                   __x = 0;
1209                 else if (__u <= __c5)
1210                   __x = 1;
1211                 else
1212                   {
1213                     const double __v = -std::log(__aurng());
1214                     const double __y = __param._M_d
1215                                      + __v * __2cx / __param._M_d;
1216                     __x = std::ceil(__y);
1217                     __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1218                   }
1219
1220                 __reject = (__w - __e - __x * __param._M_lm_thr
1221                             > __param._M_lfm - std::lgamma(__x + __m + 1));
1222
1223                 __reject |= __x + __m >= __thr;
1224
1225               } while (__reject);
1226
1227             return result_type(__x + __m + __naf);
1228           }
1229         else
1230 #endif
1231           {
1232             _IntType     __x = 0;
1233             double __prod = 1.0;
1234
1235             do
1236               {
1237                 __prod *= __aurng();
1238                 __x += 1;
1239               }
1240             while (__prod > __param._M_lm_thr);
1241
1242             return __x - 1;
1243           }
1244       }
1245
1246   template<typename _IntType,
1247            typename _CharT, typename _Traits>
1248     std::basic_ostream<_CharT, _Traits>&
1249     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1250                const poisson_distribution<_IntType>& __x)
1251     {
1252       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1253       typedef typename __ostream_type::ios_base    __ios_base;
1254
1255       const typename __ios_base::fmtflags __flags = __os.flags();
1256       const _CharT __fill = __os.fill();
1257       const std::streamsize __precision = __os.precision();
1258       const _CharT __space = __os.widen(' ');
1259       __os.flags(__ios_base::scientific | __ios_base::left);
1260       __os.fill(__space);
1261       __os.precision(std::numeric_limits<double>::max_digits10);
1262
1263       __os << __x.mean() << __space << __x._M_nd;
1264
1265       __os.flags(__flags);
1266       __os.fill(__fill);
1267       __os.precision(__precision);
1268       return __os;
1269     }
1270
1271   template<typename _IntType,
1272            typename _CharT, typename _Traits>
1273     std::basic_istream<_CharT, _Traits>&
1274     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1275                poisson_distribution<_IntType>& __x)
1276     {
1277       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1278       typedef typename __istream_type::ios_base    __ios_base;
1279
1280       const typename __ios_base::fmtflags __flags = __is.flags();
1281       __is.flags(__ios_base::skipws);
1282
1283       double __mean;
1284       __is >> __mean >> __x._M_nd;
1285       __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
1286
1287       __is.flags(__flags);
1288       return __is;
1289     }
1290
1291
1292   template<typename _IntType>
1293     void
1294     binomial_distribution<_IntType>::param_type::
1295     _M_initialize()
1296     {
1297       const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1298
1299       _M_easy = true;
1300
1301 #if _GLIBCXX_USE_C99_MATH_TR1
1302       if (_M_t * __p12 >= 8)
1303         {
1304           _M_easy = false;
1305           const double __np = std::floor(_M_t * __p12);
1306           const double __pa = __np / _M_t;
1307           const double __1p = 1 - __pa;
1308
1309           const double __pi_4 = 0.7853981633974483096156608458198757L;
1310           const double __d1x =
1311             std::sqrt(__np * __1p * std::log(32 * __np
1312                                              / (81 * __pi_4 * __1p)));
1313           _M_d1 = std::round(std::max(1.0, __d1x));
1314           const double __d2x =
1315             std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1316                                              / (__pi_4 * __pa)));
1317           _M_d2 = std::round(std::max(1.0, __d2x));
1318
1319           // sqrt(pi / 2)
1320           const double __spi_2 = 1.2533141373155002512078826424055226L;
1321           _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1322           _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1323           _M_c = 2 * _M_d1 / __np;
1324           _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1325           const double __a12 = _M_a1 + _M_s2 * __spi_2;
1326           const double __s1s = _M_s1 * _M_s1;
1327           _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1328                              * 2 * __s1s / _M_d1
1329                              * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1330           const double __s2s = _M_s2 * _M_s2;
1331           _M_s = (_M_a123 + 2 * __s2s / _M_d2
1332                   * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1333           _M_lf = (std::lgamma(__np + 1)
1334                    + std::lgamma(_M_t - __np + 1));
1335           _M_lp1p = std::log(__pa / __1p);
1336
1337           _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1338         }
1339       else
1340 #endif
1341         _M_q = -std::log(1 - __p12);
1342     }
1343
1344   template<typename _IntType>
1345     template<typename _UniformRandomNumberGenerator>
1346       typename binomial_distribution<_IntType>::result_type
1347       binomial_distribution<_IntType>::
1348       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1349       {
1350         _IntType __x = 0;
1351         double __sum = 0.0;
1352         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1353           __aurng(__urng);
1354
1355         do
1356           {
1357             const double __e = -std::log(__aurng());
1358             __sum += __e / (__t - __x);
1359             __x += 1;
1360           }
1361         while (__sum <= _M_param._M_q);
1362
1363         return __x - 1;
1364       }
1365
1366   /**
1367    * A rejection algorithm when t * p >= 8 and a simple waiting time
1368    * method - the second in the referenced book - otherwise.
1369    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1370    * is defined.
1371    *
1372    * Reference:
1373    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1374    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1375    */
1376   template<typename _IntType>
1377     template<typename _UniformRandomNumberGenerator>
1378       typename binomial_distribution<_IntType>::result_type
1379       binomial_distribution<_IntType>::
1380       operator()(_UniformRandomNumberGenerator& __urng,
1381                  const param_type& __param)
1382       {
1383         result_type __ret;
1384         const _IntType __t = __param.t();
1385         const _IntType __p = __param.p();
1386         const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1387         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
1388           __aurng(__urng);
1389
1390 #if _GLIBCXX_USE_C99_MATH_TR1
1391         if (!__param._M_easy)
1392           {
1393             double __x;
1394
1395             // See comments above...
1396             const double __naf =
1397               (1 - std::numeric_limits<double>::epsilon()) / 2;
1398             const double __thr =
1399               std::numeric_limits<_IntType>::max() + __naf;
1400
1401             const double __np = std::floor(__t * __p12);
1402
1403             // sqrt(pi / 2)
1404             const double __spi_2 = 1.2533141373155002512078826424055226L;
1405             const double __a1 = __param._M_a1;
1406             const double __a12 = __a1 + __param._M_s2 * __spi_2;
1407             const double __a123 = __param._M_a123;
1408             const double __s1s = __param._M_s1 * __param._M_s1;
1409             const double __s2s = __param._M_s2 * __param._M_s2;
1410
1411             bool __reject;
1412             do
1413               {
1414                 const double __u = __param._M_s * __aurng();
1415
1416                 double __v;
1417
1418                 if (__u <= __a1)
1419                   {
1420                     const double __n = _M_nd(__urng);
1421                     const double __y = __param._M_s1 * std::abs(__n);
1422                     __reject = __y >= __param._M_d1;
1423                     if (!__reject)
1424                       {
1425                         const double __e = -std::log(__aurng());
1426                         __x = std::floor(__y);
1427                         __v = -__e - __n * __n / 2 + __param._M_c;
1428                       }
1429                   }
1430                 else if (__u <= __a12)
1431                   {
1432                     const double __n = _M_nd(__urng);
1433                     const double __y = __param._M_s2 * std::abs(__n);
1434                     __reject = __y >= __param._M_d2;
1435                     if (!__reject)
1436                       {
1437                         const double __e = -std::log(__aurng());
1438                         __x = std::floor(-__y);
1439                         __v = -__e - __n * __n / 2;
1440                       }
1441                   }
1442                 else if (__u <= __a123)
1443                   {
1444                     const double __e1 = -std::log(__aurng());
1445                     const double __e2 = -std::log(__aurng());
1446
1447                     const double __y = __param._M_d1
1448                                      + 2 * __s1s * __e1 / __param._M_d1;
1449                     __x = std::floor(__y);
1450                     __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1451                                                     -__y / (2 * __s1s)));
1452                     __reject = false;
1453                   }
1454                 else
1455                   {
1456                     const double __e1 = -std::log(__aurng());
1457                     const double __e2 = -std::log(__aurng());
1458
1459                     const double __y = __param._M_d2
1460                                      + 2 * __s2s * __e1 / __param._M_d2;
1461                     __x = std::floor(-__y);
1462                     __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1463                     __reject = false;
1464                   }
1465
1466                 __reject = __reject || __x < -__np || __x > __t - __np;
1467                 if (!__reject)
1468                   {
1469                     const double __lfx =
1470                       std::lgamma(__np + __x + 1)
1471                       + std::lgamma(__t - (__np + __x) + 1);
1472                     __reject = __v > __param._M_lf - __lfx
1473                              + __x * __param._M_lp1p;
1474                   }
1475
1476                 __reject |= __x + __np >= __thr;
1477               }
1478             while (__reject);
1479
1480             __x += __np + __naf;
1481
1482             const _IntType __z = _M_waiting(__urng, __t - _IntType(__x));
1483             __ret = _IntType(__x) + __z;
1484           }
1485         else
1486 #endif
1487           __ret = _M_waiting(__urng, __t);
1488
1489         if (__p12 != __p)
1490           __ret = __t - __ret;
1491         return __ret;
1492       }
1493
1494   template<typename _IntType,
1495            typename _CharT, typename _Traits>
1496     std::basic_ostream<_CharT, _Traits>&
1497     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1498                const binomial_distribution<_IntType>& __x)
1499     {
1500       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1501       typedef typename __ostream_type::ios_base    __ios_base;
1502
1503       const typename __ios_base::fmtflags __flags = __os.flags();
1504       const _CharT __fill = __os.fill();
1505       const std::streamsize __precision = __os.precision();
1506       const _CharT __space = __os.widen(' ');
1507       __os.flags(__ios_base::scientific | __ios_base::left);
1508       __os.fill(__space);
1509       __os.precision(std::numeric_limits<double>::max_digits10);
1510
1511       __os << __x.t() << __space << __x.p()
1512            << __space << __x._M_nd;
1513
1514       __os.flags(__flags);
1515       __os.fill(__fill);
1516       __os.precision(__precision);
1517       return __os;
1518     }
1519
1520   template<typename _IntType,
1521            typename _CharT, typename _Traits>
1522     std::basic_istream<_CharT, _Traits>&
1523     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1524                binomial_distribution<_IntType>& __x)
1525     {
1526       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1527       typedef typename __istream_type::ios_base    __ios_base;
1528
1529       const typename __ios_base::fmtflags __flags = __is.flags();
1530       __is.flags(__ios_base::dec | __ios_base::skipws);
1531
1532       _IntType __t;
1533       double __p;
1534       __is >> __t >> __p >> __x._M_nd;
1535       __x.param(typename binomial_distribution<_IntType>::
1536                 param_type(__t, __p));
1537
1538       __is.flags(__flags);
1539       return __is;
1540     }
1541
1542
1543   template<typename _RealType, typename _CharT, typename _Traits>
1544     std::basic_ostream<_CharT, _Traits>&
1545     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1546                const exponential_distribution<_RealType>& __x)
1547     {
1548       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1549       typedef typename __ostream_type::ios_base    __ios_base;
1550
1551       const typename __ios_base::fmtflags __flags = __os.flags();
1552       const _CharT __fill = __os.fill();
1553       const std::streamsize __precision = __os.precision();
1554       __os.flags(__ios_base::scientific | __ios_base::left);
1555       __os.fill(__os.widen(' '));
1556       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1557
1558       __os << __x.lambda();
1559
1560       __os.flags(__flags);
1561       __os.fill(__fill);
1562       __os.precision(__precision);
1563       return __os;
1564     }
1565
1566   template<typename _RealType, typename _CharT, typename _Traits>
1567     std::basic_istream<_CharT, _Traits>&
1568     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1569                exponential_distribution<_RealType>& __x)
1570     {
1571       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1572       typedef typename __istream_type::ios_base    __ios_base;
1573
1574       const typename __ios_base::fmtflags __flags = __is.flags();
1575       __is.flags(__ios_base::dec | __ios_base::skipws);
1576
1577       _RealType __lambda;
1578       __is >> __lambda;
1579       __x.param(typename exponential_distribution<_RealType>::
1580                 param_type(__lambda));
1581
1582       __is.flags(__flags);
1583       return __is;
1584     }
1585
1586
1587   /**
1588    * Polar method due to Marsaglia.
1589    *
1590    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1591    * New York, 1986, Ch. V, Sect. 4.4.
1592    */
1593   template<typename _RealType>
1594     template<typename _UniformRandomNumberGenerator>
1595       typename normal_distribution<_RealType>::result_type
1596       normal_distribution<_RealType>::
1597       operator()(_UniformRandomNumberGenerator& __urng,
1598                  const param_type& __param)
1599       {
1600         result_type __ret;
1601         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1602           __aurng(__urng);
1603
1604         if (_M_saved_available)
1605           {
1606             _M_saved_available = false;
1607             __ret = _M_saved;
1608           }
1609         else
1610           {
1611             result_type __x, __y, __r2;
1612             do
1613               {
1614                 __x = result_type(2.0) * __aurng() - 1.0;
1615                 __y = result_type(2.0) * __aurng() - 1.0;
1616                 __r2 = __x * __x + __y * __y;
1617               }
1618             while (__r2 > 1.0 || __r2 == 0.0);
1619
1620             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1621             _M_saved = __x * __mult;
1622             _M_saved_available = true;
1623             __ret = __y * __mult;
1624           }
1625
1626         __ret = __ret * __param.stddev() + __param.mean();
1627         return __ret;
1628       }
1629
1630   template<typename _RealType, typename _CharT, typename _Traits>
1631     std::basic_ostream<_CharT, _Traits>&
1632     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1633                const normal_distribution<_RealType>& __x)
1634     {
1635       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1636       typedef typename __ostream_type::ios_base    __ios_base;
1637
1638       const typename __ios_base::fmtflags __flags = __os.flags();
1639       const _CharT __fill = __os.fill();
1640       const std::streamsize __precision = __os.precision();
1641       const _CharT __space = __os.widen(' ');
1642       __os.flags(__ios_base::scientific | __ios_base::left);
1643       __os.fill(__space);
1644       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1645
1646       __os << __x.mean() << __space << __x.stddev()
1647            << __space << __x._M_saved_available;
1648       if (__x._M_saved_available)
1649         __os << __space << __x._M_saved;
1650
1651       __os.flags(__flags);
1652       __os.fill(__fill);
1653       __os.precision(__precision);
1654       return __os;
1655     }
1656
1657   template<typename _RealType, typename _CharT, typename _Traits>
1658     std::basic_istream<_CharT, _Traits>&
1659     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1660                normal_distribution<_RealType>& __x)
1661     {
1662       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1663       typedef typename __istream_type::ios_base    __ios_base;
1664
1665       const typename __ios_base::fmtflags __flags = __is.flags();
1666       __is.flags(__ios_base::dec | __ios_base::skipws);
1667
1668       double __mean, __stddev;
1669       __is >> __mean >> __stddev
1670            >> __x._M_saved_available;
1671       if (__x._M_saved_available)
1672         __is >> __x._M_saved;
1673       __x.param(typename normal_distribution<_RealType>::
1674                 param_type(__mean, __stddev));
1675
1676       __is.flags(__flags);
1677       return __is;
1678     }
1679
1680
1681   template<typename _RealType, typename _CharT, typename _Traits>
1682     std::basic_ostream<_CharT, _Traits>&
1683     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1684                const lognormal_distribution<_RealType>& __x)
1685     {
1686       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1687       typedef typename __ostream_type::ios_base    __ios_base;
1688
1689       const typename __ios_base::fmtflags __flags = __os.flags();
1690       const _CharT __fill = __os.fill();
1691       const std::streamsize __precision = __os.precision();
1692       const _CharT __space = __os.widen(' ');
1693       __os.flags(__ios_base::scientific | __ios_base::left);
1694       __os.fill(__space);
1695       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1696
1697       __os << __x.m() << __space << __x.s()
1698            << __space << __x._M_nd;
1699
1700       __os.flags(__flags);
1701       __os.fill(__fill);
1702       __os.precision(__precision);
1703       return __os;
1704     }
1705
1706   template<typename _RealType, typename _CharT, typename _Traits>
1707     std::basic_istream<_CharT, _Traits>&
1708     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1709                lognormal_distribution<_RealType>& __x)
1710     {
1711       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1712       typedef typename __istream_type::ios_base    __ios_base;
1713
1714       const typename __ios_base::fmtflags __flags = __is.flags();
1715       __is.flags(__ios_base::dec | __ios_base::skipws);
1716
1717       _RealType __m, __s;
1718       __is >> __m >> __s >> __x._M_nd;
1719       __x.param(typename lognormal_distribution<_RealType>::
1720                 param_type(__m, __s));
1721
1722       __is.flags(__flags);
1723       return __is;
1724     }
1725
1726
1727   template<typename _RealType, typename _CharT, typename _Traits>
1728     std::basic_ostream<_CharT, _Traits>&
1729     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1730                const chi_squared_distribution<_RealType>& __x)
1731     {
1732       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1733       typedef typename __ostream_type::ios_base    __ios_base;
1734
1735       const typename __ios_base::fmtflags __flags = __os.flags();
1736       const _CharT __fill = __os.fill();
1737       const std::streamsize __precision = __os.precision();
1738       const _CharT __space = __os.widen(' ');
1739       __os.flags(__ios_base::scientific | __ios_base::left);
1740       __os.fill(__space);
1741       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1742
1743       __os << __x.n() << __space << __x._M_gd;
1744
1745       __os.flags(__flags);
1746       __os.fill(__fill);
1747       __os.precision(__precision);
1748       return __os;
1749     }
1750
1751   template<typename _RealType, typename _CharT, typename _Traits>
1752     std::basic_istream<_CharT, _Traits>&
1753     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1754                chi_squared_distribution<_RealType>& __x)
1755     {
1756       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1757       typedef typename __istream_type::ios_base    __ios_base;
1758
1759       const typename __ios_base::fmtflags __flags = __is.flags();
1760       __is.flags(__ios_base::dec | __ios_base::skipws);
1761
1762       _RealType __n;
1763       __is >> __n >> __x._M_gd;
1764       __x.param(typename chi_squared_distribution<_RealType>::
1765                 param_type(__n));
1766
1767       __is.flags(__flags);
1768       return __is;
1769     }
1770
1771
1772   template<typename _RealType>
1773     template<typename _UniformRandomNumberGenerator>
1774       typename cauchy_distribution<_RealType>::result_type
1775       cauchy_distribution<_RealType>::
1776       operator()(_UniformRandomNumberGenerator& __urng,
1777                  const param_type& __p)
1778       {
1779         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1780           __aurng(__urng);
1781         _RealType __u;
1782         do
1783           __u = __aurng();
1784         while (__u == 0.5);
1785
1786         const _RealType __pi = 3.1415926535897932384626433832795029L;
1787         return __p.a() + __p.b() * std::tan(__pi * __u);
1788       }
1789
1790   template<typename _RealType, typename _CharT, typename _Traits>
1791     std::basic_ostream<_CharT, _Traits>&
1792     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1793                const cauchy_distribution<_RealType>& __x)
1794     {
1795       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1796       typedef typename __ostream_type::ios_base    __ios_base;
1797
1798       const typename __ios_base::fmtflags __flags = __os.flags();
1799       const _CharT __fill = __os.fill();
1800       const std::streamsize __precision = __os.precision();
1801       const _CharT __space = __os.widen(' ');
1802       __os.flags(__ios_base::scientific | __ios_base::left);
1803       __os.fill(__space);
1804       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1805
1806       __os << __x.a() << __space << __x.b();
1807
1808       __os.flags(__flags);
1809       __os.fill(__fill);
1810       __os.precision(__precision);
1811       return __os;
1812     }
1813
1814   template<typename _RealType, typename _CharT, typename _Traits>
1815     std::basic_istream<_CharT, _Traits>&
1816     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1817                cauchy_distribution<_RealType>& __x)
1818     {
1819       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1820       typedef typename __istream_type::ios_base    __ios_base;
1821
1822       const typename __ios_base::fmtflags __flags = __is.flags();
1823       __is.flags(__ios_base::dec | __ios_base::skipws);
1824
1825       _RealType __a, __b;
1826       __is >> __a >> __b;
1827       __x.param(typename cauchy_distribution<_RealType>::
1828                 param_type(__a, __b));
1829
1830       __is.flags(__flags);
1831       return __is;
1832     }
1833
1834
1835   template<typename _RealType, typename _CharT, typename _Traits>
1836     std::basic_ostream<_CharT, _Traits>&
1837     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1838                const fisher_f_distribution<_RealType>& __x)
1839     {
1840       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1841       typedef typename __ostream_type::ios_base    __ios_base;
1842
1843       const typename __ios_base::fmtflags __flags = __os.flags();
1844       const _CharT __fill = __os.fill();
1845       const std::streamsize __precision = __os.precision();
1846       const _CharT __space = __os.widen(' ');
1847       __os.flags(__ios_base::scientific | __ios_base::left);
1848       __os.fill(__space);
1849       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1850
1851       __os << __x.m() << __space << __x.n()
1852            << __space << __x._M_gd_x << __space << __x._M_gd_y;
1853
1854       __os.flags(__flags);
1855       __os.fill(__fill);
1856       __os.precision(__precision);
1857       return __os;
1858     }
1859
1860   template<typename _RealType, typename _CharT, typename _Traits>
1861     std::basic_istream<_CharT, _Traits>&
1862     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1863                fisher_f_distribution<_RealType>& __x)
1864     {
1865       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1866       typedef typename __istream_type::ios_base    __ios_base;
1867
1868       const typename __ios_base::fmtflags __flags = __is.flags();
1869       __is.flags(__ios_base::dec | __ios_base::skipws);
1870
1871       _RealType __m, __n;
1872       __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
1873       __x.param(typename fisher_f_distribution<_RealType>::
1874                 param_type(__m, __n));
1875
1876       __is.flags(__flags);
1877       return __is;
1878     }
1879
1880
1881   template<typename _RealType, typename _CharT, typename _Traits>
1882     std::basic_ostream<_CharT, _Traits>&
1883     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1884                const student_t_distribution<_RealType>& __x)
1885     {
1886       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1887       typedef typename __ostream_type::ios_base    __ios_base;
1888
1889       const typename __ios_base::fmtflags __flags = __os.flags();
1890       const _CharT __fill = __os.fill();
1891       const std::streamsize __precision = __os.precision();
1892       const _CharT __space = __os.widen(' ');
1893       __os.flags(__ios_base::scientific | __ios_base::left);
1894       __os.fill(__space);
1895       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1896
1897       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
1898
1899       __os.flags(__flags);
1900       __os.fill(__fill);
1901       __os.precision(__precision);
1902       return __os;
1903     }
1904
1905   template<typename _RealType, typename _CharT, typename _Traits>
1906     std::basic_istream<_CharT, _Traits>&
1907     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1908                student_t_distribution<_RealType>& __x)
1909     {
1910       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1911       typedef typename __istream_type::ios_base    __ios_base;
1912
1913       const typename __ios_base::fmtflags __flags = __is.flags();
1914       __is.flags(__ios_base::dec | __ios_base::skipws);
1915
1916       _RealType __n;
1917       __is >> __n >> __x._M_nd >> __x._M_gd;
1918       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
1919
1920       __is.flags(__flags);
1921       return __is;
1922     }
1923
1924
1925   template<typename _RealType>
1926     void
1927     gamma_distribution<_RealType>::param_type::
1928     _M_initialize()
1929     {
1930       _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
1931
1932       const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
1933       _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
1934     }
1935
1936   /**
1937    * Marsaglia, G. and Tsang, W. W.
1938    * "A Simple Method for Generating Gamma Variables"
1939    * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
1940    */
1941   template<typename _RealType>
1942     template<typename _UniformRandomNumberGenerator>
1943       typename gamma_distribution<_RealType>::result_type
1944       gamma_distribution<_RealType>::
1945       operator()(_UniformRandomNumberGenerator& __urng,
1946                  const param_type& __param)
1947       {
1948         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1949           __aurng(__urng);
1950
1951         result_type __u, __v, __n;
1952         const result_type __a1 = (__param._M_malpha
1953                                   - _RealType(1.0) / _RealType(3.0));
1954
1955         do
1956           {
1957             do
1958               {
1959                 __n = _M_nd(__urng);
1960                 __v = result_type(1.0) + __param._M_a2 * __n;
1961               }
1962             while (__v <= 0.0);
1963
1964             __v = __v * __v * __v;
1965             __u = __aurng();
1966           }
1967         while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
1968                && (std::log(__u) > (0.5 * __n * __n + __a1
1969                                     * (1.0 - __v + std::log(__v)))));
1970
1971         if (__param.alpha() == __param._M_malpha)
1972           return __a1 * __v * __param.beta();
1973         else
1974           {
1975             do
1976               __u = __aurng();
1977             while (__u == 0.0);
1978
1979             return (std::pow(__u, result_type(1.0) / __param.alpha())
1980                     * __a1 * __v * __param.beta());
1981           }
1982       }
1983
1984   template<typename _RealType, typename _CharT, typename _Traits>
1985     std::basic_ostream<_CharT, _Traits>&
1986     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1987                const gamma_distribution<_RealType>& __x)
1988     {
1989       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1990       typedef typename __ostream_type::ios_base    __ios_base;
1991
1992       const typename __ios_base::fmtflags __flags = __os.flags();
1993       const _CharT __fill = __os.fill();
1994       const std::streamsize __precision = __os.precision();
1995       const _CharT __space = __os.widen(' ');
1996       __os.flags(__ios_base::scientific | __ios_base::left);
1997       __os.fill(__space);
1998       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1999
2000       __os << __x.alpha() << __space << __x.beta()
2001            << __space << __x._M_nd;
2002
2003       __os.flags(__flags);
2004       __os.fill(__fill);
2005       __os.precision(__precision);
2006       return __os;
2007     }
2008
2009   template<typename _RealType, typename _CharT, typename _Traits>
2010     std::basic_istream<_CharT, _Traits>&
2011     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2012                gamma_distribution<_RealType>& __x)
2013     {
2014       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2015       typedef typename __istream_type::ios_base    __ios_base;
2016
2017       const typename __ios_base::fmtflags __flags = __is.flags();
2018       __is.flags(__ios_base::dec | __ios_base::skipws);
2019
2020       _RealType __alpha_val, __beta_val;
2021       __is >> __alpha_val >> __beta_val >> __x._M_nd;
2022       __x.param(typename gamma_distribution<_RealType>::
2023                 param_type(__alpha_val, __beta_val));
2024
2025       __is.flags(__flags);
2026       return __is;
2027     }
2028
2029
2030   template<typename _RealType>
2031     template<typename _UniformRandomNumberGenerator>
2032       typename weibull_distribution<_RealType>::result_type
2033       weibull_distribution<_RealType>::
2034       operator()(_UniformRandomNumberGenerator& __urng,
2035                  const param_type& __p)
2036       {
2037         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2038           __aurng(__urng);
2039         return __p.b() * std::pow(-std::log(__aurng()),
2040                                   result_type(1) / __p.a());
2041       }
2042
2043   template<typename _RealType, typename _CharT, typename _Traits>
2044     std::basic_ostream<_CharT, _Traits>&
2045     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2046                const weibull_distribution<_RealType>& __x)
2047     {
2048       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2049       typedef typename __ostream_type::ios_base    __ios_base;
2050
2051       const typename __ios_base::fmtflags __flags = __os.flags();
2052       const _CharT __fill = __os.fill();
2053       const std::streamsize __precision = __os.precision();
2054       const _CharT __space = __os.widen(' ');
2055       __os.flags(__ios_base::scientific | __ios_base::left);
2056       __os.fill(__space);
2057       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2058
2059       __os << __x.a() << __space << __x.b();
2060
2061       __os.flags(__flags);
2062       __os.fill(__fill);
2063       __os.precision(__precision);
2064       return __os;
2065     }
2066
2067   template<typename _RealType, typename _CharT, typename _Traits>
2068     std::basic_istream<_CharT, _Traits>&
2069     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2070                weibull_distribution<_RealType>& __x)
2071     {
2072       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2073       typedef typename __istream_type::ios_base    __ios_base;
2074
2075       const typename __ios_base::fmtflags __flags = __is.flags();
2076       __is.flags(__ios_base::dec | __ios_base::skipws);
2077
2078       _RealType __a, __b;
2079       __is >> __a >> __b;
2080       __x.param(typename weibull_distribution<_RealType>::
2081                 param_type(__a, __b));
2082
2083       __is.flags(__flags);
2084       return __is;
2085     }
2086
2087
2088   template<typename _RealType>
2089     template<typename _UniformRandomNumberGenerator>
2090       typename extreme_value_distribution<_RealType>::result_type
2091       extreme_value_distribution<_RealType>::
2092       operator()(_UniformRandomNumberGenerator& __urng,
2093                  const param_type& __p)
2094       {
2095         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2096           __aurng(__urng);
2097         return __p.a() - __p.b() * std::log(-std::log(__aurng()));
2098       }
2099
2100   template<typename _RealType, typename _CharT, typename _Traits>
2101     std::basic_ostream<_CharT, _Traits>&
2102     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2103                const extreme_value_distribution<_RealType>& __x)
2104     {
2105       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2106       typedef typename __ostream_type::ios_base    __ios_base;
2107
2108       const typename __ios_base::fmtflags __flags = __os.flags();
2109       const _CharT __fill = __os.fill();
2110       const std::streamsize __precision = __os.precision();
2111       const _CharT __space = __os.widen(' ');
2112       __os.flags(__ios_base::scientific | __ios_base::left);
2113       __os.fill(__space);
2114       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2115
2116       __os << __x.a() << __space << __x.b();
2117
2118       __os.flags(__flags);
2119       __os.fill(__fill);
2120       __os.precision(__precision);
2121       return __os;
2122     }
2123
2124   template<typename _RealType, typename _CharT, typename _Traits>
2125     std::basic_istream<_CharT, _Traits>&
2126     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2127                extreme_value_distribution<_RealType>& __x)
2128     {
2129       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2130       typedef typename __istream_type::ios_base    __ios_base;
2131
2132       const typename __ios_base::fmtflags __flags = __is.flags();
2133       __is.flags(__ios_base::dec | __ios_base::skipws);
2134
2135       _RealType __a, __b;
2136       __is >> __a >> __b;
2137       __x.param(typename extreme_value_distribution<_RealType>::
2138                 param_type(__a, __b));
2139
2140       __is.flags(__flags);
2141       return __is;
2142     }
2143
2144
2145   template<typename _IntType>
2146     void
2147     discrete_distribution<_IntType>::param_type::
2148     _M_initialize()
2149     {
2150       if (_M_prob.size() < 2)
2151         {
2152           _M_prob.clear();
2153           _M_prob.push_back(1.0);
2154           return;
2155         }
2156
2157       const double __sum = std::accumulate(_M_prob.begin(),
2158                                            _M_prob.end(), 0.0);
2159       // Now normalize the probabilites.
2160       std::transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
2161                      std::bind2nd(std::divides<double>(), __sum));
2162       // Accumulate partial sums.
2163       _M_cp.reserve(_M_prob.size());
2164       std::partial_sum(_M_prob.begin(), _M_prob.end(),
2165                        std::back_inserter(_M_cp));
2166       // Make sure the last cumulative probability is one.
2167       _M_cp[_M_cp.size() - 1] = 1.0;
2168     }
2169
2170   template<typename _IntType>
2171     template<typename _Func>
2172       discrete_distribution<_IntType>::param_type::
2173       param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2174       : _M_prob(), _M_cp()
2175       {
2176         const size_t __n = __nw == 0 ? 1 : __nw;
2177         const double __delta = (__xmax - __xmin) / __n;
2178
2179         _M_prob.reserve(__n);
2180         for (size_t __k = 0; __k < __nw; ++__k)
2181           _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2182
2183         _M_initialize();
2184       }
2185
2186   template<typename _IntType>
2187     template<typename _UniformRandomNumberGenerator>
2188       typename discrete_distribution<_IntType>::result_type
2189       discrete_distribution<_IntType>::
2190       operator()(_UniformRandomNumberGenerator& __urng,
2191                  const param_type& __param)
2192       {
2193         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2194           __aurng(__urng);
2195
2196         const double __p = __aurng();
2197         auto __pos = std::lower_bound(__param._M_cp.begin(),
2198                                       __param._M_cp.end(), __p);
2199
2200         return __pos - __param._M_cp.begin();
2201       }
2202
2203   template<typename _IntType, typename _CharT, typename _Traits>
2204     std::basic_ostream<_CharT, _Traits>&
2205     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2206                const discrete_distribution<_IntType>& __x)
2207     {
2208       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2209       typedef typename __ostream_type::ios_base    __ios_base;
2210
2211       const typename __ios_base::fmtflags __flags = __os.flags();
2212       const _CharT __fill = __os.fill();
2213       const std::streamsize __precision = __os.precision();
2214       const _CharT __space = __os.widen(' ');
2215       __os.flags(__ios_base::scientific | __ios_base::left);
2216       __os.fill(__space);
2217       __os.precision(std::numeric_limits<double>::max_digits10);
2218
2219       std::vector<double> __prob = __x.probabilities();
2220       __os << __prob.size();
2221       for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2222         __os << __space << *__dit;
2223
2224       __os.flags(__flags);
2225       __os.fill(__fill);
2226       __os.precision(__precision);
2227       return __os;
2228     }
2229
2230   template<typename _IntType, typename _CharT, typename _Traits>
2231     std::basic_istream<_CharT, _Traits>&
2232     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2233                discrete_distribution<_IntType>& __x)
2234     {
2235       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2236       typedef typename __istream_type::ios_base    __ios_base;
2237
2238       const typename __ios_base::fmtflags __flags = __is.flags();
2239       __is.flags(__ios_base::dec | __ios_base::skipws);
2240
2241       size_t __n;
2242       __is >> __n;
2243
2244       std::vector<double> __prob_vec;
2245       __prob_vec.reserve(__n);
2246       for (; __n != 0; --__n)
2247         {
2248           double __prob;
2249           __is >> __prob;
2250           __prob_vec.push_back(__prob);
2251         }
2252
2253       __x.param(typename discrete_distribution<_IntType>::
2254                 param_type(__prob_vec.begin(), __prob_vec.end()));
2255
2256       __is.flags(__flags);
2257       return __is;
2258     }
2259
2260
2261   template<typename _RealType>
2262     void
2263     piecewise_constant_distribution<_RealType>::param_type::
2264     _M_initialize()
2265     {
2266       if (_M_int.size() < 2)
2267         {
2268           _M_int.clear();
2269           _M_int.reserve(2);
2270           _M_int.push_back(_RealType(0));
2271           _M_int.push_back(_RealType(1));
2272
2273           _M_den.clear();
2274           _M_den.push_back(1.0);
2275
2276           return;
2277         }
2278
2279       const double __sum = std::accumulate(_M_den.begin(),
2280                                            _M_den.end(), 0.0);
2281
2282       std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2283                      std::bind2nd(std::divides<double>(), __sum));
2284
2285       _M_cp.reserve(_M_den.size());
2286       std::partial_sum(_M_den.begin(), _M_den.end(),
2287                        std::back_inserter(_M_cp));
2288
2289       // Make sure the last cumulative probability is one.
2290       _M_cp[_M_cp.size() - 1] = 1.0;
2291
2292       for (size_t __k = 0; __k < _M_den.size(); ++__k)
2293         _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2294     }
2295
2296   template<typename _RealType>
2297     template<typename _InputIteratorB, typename _InputIteratorW>
2298       piecewise_constant_distribution<_RealType>::param_type::
2299       param_type(_InputIteratorB __bbegin,
2300                  _InputIteratorB __bend,
2301                  _InputIteratorW __wbegin)
2302       : _M_int(), _M_den(), _M_cp()
2303       {
2304         if (__bbegin != __bend)
2305           {
2306             for (;;)
2307               {
2308                 _M_int.push_back(*__bbegin);
2309                 ++__bbegin;
2310                 if (__bbegin == __bend)
2311                   break;
2312
2313                 _M_den.push_back(*__wbegin);
2314                 ++__wbegin;
2315               }
2316           }
2317
2318         _M_initialize();
2319       }
2320
2321   template<typename _RealType>
2322     template<typename _Func>
2323       piecewise_constant_distribution<_RealType>::param_type::
2324       param_type(initializer_list<_RealType> __bl, _Func __fw)
2325       : _M_int(), _M_den(), _M_cp()
2326       {
2327         _M_int.reserve(__bl.size());
2328         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2329           _M_int.push_back(*__biter);
2330
2331         _M_den.reserve(_M_int.size() - 1);
2332         for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2333           _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2334
2335         _M_initialize();
2336       }
2337
2338   template<typename _RealType>
2339     template<typename _Func>
2340       piecewise_constant_distribution<_RealType>::param_type::
2341       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2342       : _M_int(), _M_den(), _M_cp()
2343       {
2344         const size_t __n = __nw == 0 ? 1 : __nw;
2345         const _RealType __delta = (__xmax - __xmin) / __n;
2346
2347         _M_int.reserve(__n + 1);
2348         for (size_t __k = 0; __k <= __nw; ++__k)
2349           _M_int.push_back(__xmin + __k * __delta);
2350
2351         _M_den.reserve(__n);
2352         for (size_t __k = 0; __k < __nw; ++__k)
2353           _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2354
2355         _M_initialize();
2356       }
2357
2358   template<typename _RealType>
2359     template<typename _UniformRandomNumberGenerator>
2360       typename piecewise_constant_distribution<_RealType>::result_type
2361       piecewise_constant_distribution<_RealType>::
2362       operator()(_UniformRandomNumberGenerator& __urng,
2363                  const param_type& __param)
2364       {
2365         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2366           __aurng(__urng);
2367
2368         const double __p = __aurng();
2369         auto __pos = std::lower_bound(__param._M_cp.begin(),
2370                                       __param._M_cp.end(), __p);
2371         const size_t __i = __pos - __param._M_cp.begin();
2372
2373         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2374
2375         return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2376       }
2377
2378   template<typename _RealType, typename _CharT, typename _Traits>
2379     std::basic_ostream<_CharT, _Traits>&
2380     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2381                const piecewise_constant_distribution<_RealType>& __x)
2382     {
2383       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2384       typedef typename __ostream_type::ios_base    __ios_base;
2385
2386       const typename __ios_base::fmtflags __flags = __os.flags();
2387       const _CharT __fill = __os.fill();
2388       const std::streamsize __precision = __os.precision();
2389       const _CharT __space = __os.widen(' ');
2390       __os.flags(__ios_base::scientific | __ios_base::left);
2391       __os.fill(__space);
2392       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2393
2394       std::vector<_RealType> __int = __x.intervals();
2395       __os << __int.size() - 1;
2396
2397       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2398         __os << __space << *__xit;
2399
2400       std::vector<double> __den = __x.densities();
2401       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2402         __os << __space << *__dit;
2403
2404       __os.flags(__flags);
2405       __os.fill(__fill);
2406       __os.precision(__precision);
2407       return __os;
2408     }
2409
2410   template<typename _RealType, typename _CharT, typename _Traits>
2411     std::basic_istream<_CharT, _Traits>&
2412     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2413                piecewise_constant_distribution<_RealType>& __x)
2414     {
2415       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2416       typedef typename __istream_type::ios_base    __ios_base;
2417
2418       const typename __ios_base::fmtflags __flags = __is.flags();
2419       __is.flags(__ios_base::dec | __ios_base::skipws);
2420
2421       size_t __n;
2422       __is >> __n;
2423
2424       std::vector<_RealType> __int_vec;
2425       __int_vec.reserve(__n + 1);
2426       for (size_t __i = 0; __i <= __n; ++__i)
2427         {
2428           _RealType __int;
2429           __is >> __int;
2430           __int_vec.push_back(__int);
2431         }
2432
2433       std::vector<double> __den_vec;
2434       __den_vec.reserve(__n);
2435       for (size_t __i = 0; __i < __n; ++__i)
2436         {
2437           double __den;
2438           __is >> __den;
2439           __den_vec.push_back(__den);
2440         }
2441
2442       __x.param(typename piecewise_constant_distribution<_RealType>::
2443           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2444
2445       __is.flags(__flags);
2446       return __is;
2447     }
2448
2449
2450   template<typename _RealType>
2451     void
2452     piecewise_linear_distribution<_RealType>::param_type::
2453     _M_initialize()
2454     {
2455       if (_M_int.size() < 2)
2456         {
2457           _M_int.clear();
2458           _M_int.reserve(2);
2459           _M_int.push_back(_RealType(0));
2460           _M_int.push_back(_RealType(1));
2461
2462           _M_den.clear();
2463           _M_den.reserve(2);
2464           _M_den.push_back(1.0);
2465           _M_den.push_back(1.0);
2466
2467           return;
2468         }
2469
2470       double __sum = 0.0;
2471       _M_cp.reserve(_M_int.size() - 1);
2472       _M_m.reserve(_M_int.size() - 1);
2473       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2474         {
2475           const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
2476           __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
2477           _M_cp.push_back(__sum);
2478           _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
2479         }
2480
2481       //  Now normalize the densities...
2482       std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
2483                      std::bind2nd(std::divides<double>(), __sum));
2484       //  ... and partial sums...
2485       std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
2486                      std::bind2nd(std::divides<double>(), __sum));
2487       //  ... and slopes.
2488       std::transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
2489                      std::bind2nd(std::divides<double>(), __sum));
2490       //  Make sure the last cumulative probablility is one.
2491       _M_cp[_M_cp.size() - 1] = 1.0;
2492      }
2493
2494   template<typename _RealType>
2495     template<typename _InputIteratorB, typename _InputIteratorW>
2496       piecewise_linear_distribution<_RealType>::param_type::
2497       param_type(_InputIteratorB __bbegin,
2498                  _InputIteratorB __bend,
2499                  _InputIteratorW __wbegin)
2500       : _M_int(), _M_den(), _M_cp(), _M_m()
2501       {
2502         for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
2503           {
2504             _M_int.push_back(*__bbegin);
2505             _M_den.push_back(*__wbegin);
2506           }
2507
2508         _M_initialize();
2509       }
2510
2511   template<typename _RealType>
2512     template<typename _Func>
2513       piecewise_linear_distribution<_RealType>::param_type::
2514       param_type(initializer_list<_RealType> __bl, _Func __fw)
2515       : _M_int(), _M_den(), _M_cp(), _M_m()
2516       {
2517         _M_int.reserve(__bl.size());
2518         _M_den.reserve(__bl.size());
2519         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2520           {
2521             _M_int.push_back(*__biter);
2522             _M_den.push_back(__fw(*__biter));
2523           }
2524
2525         _M_initialize();
2526       }
2527
2528   template<typename _RealType>
2529     template<typename _Func>
2530       piecewise_linear_distribution<_RealType>::param_type::
2531       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2532       : _M_int(), _M_den(), _M_cp(), _M_m()
2533       {
2534         const size_t __n = __nw == 0 ? 1 : __nw;
2535         const _RealType __delta = (__xmax - __xmin) / __n;
2536
2537         _M_int.reserve(__n + 1);
2538         _M_den.reserve(__n + 1);
2539         for (size_t __k = 0; __k <= __nw; ++__k)
2540           {
2541             _M_int.push_back(__xmin + __k * __delta);
2542             _M_den.push_back(__fw(_M_int[__k] + __delta));
2543           }
2544
2545         _M_initialize();
2546       }
2547
2548   template<typename _RealType>
2549     template<typename _UniformRandomNumberGenerator>
2550       typename piecewise_linear_distribution<_RealType>::result_type
2551       piecewise_linear_distribution<_RealType>::
2552       operator()(_UniformRandomNumberGenerator& __urng,
2553                  const param_type& __param)
2554       {
2555         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2556           __aurng(__urng);
2557
2558         const double __p = __aurng();
2559         auto __pos = std::lower_bound(__param._M_cp.begin(),
2560                                       __param._M_cp.end(), __p);
2561         const size_t __i = __pos - __param._M_cp.begin();
2562
2563         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2564
2565         const double __a = 0.5 * __param._M_m[__i];
2566         const double __b = __param._M_den[__i];
2567         const double __cm = __p - __pref;
2568
2569         _RealType __x = __param._M_int[__i];
2570         if (__a == 0)
2571           __x += __cm / __b;
2572         else
2573           {
2574             const double __d = __b * __b + 4.0 * __a * __cm;
2575             __x += 0.5 * (std::sqrt(__d) - __b) / __a;
2576           }
2577
2578         return __x;
2579       }
2580
2581   template<typename _RealType, typename _CharT, typename _Traits>
2582     std::basic_ostream<_CharT, _Traits>&
2583     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2584                const piecewise_linear_distribution<_RealType>& __x)
2585     {
2586       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
2587       typedef typename __ostream_type::ios_base    __ios_base;
2588
2589       const typename __ios_base::fmtflags __flags = __os.flags();
2590       const _CharT __fill = __os.fill();
2591       const std::streamsize __precision = __os.precision();
2592       const _CharT __space = __os.widen(' ');
2593       __os.flags(__ios_base::scientific | __ios_base::left);
2594       __os.fill(__space);
2595       __os.precision(std::numeric_limits<_RealType>::max_digits10);
2596
2597       std::vector<_RealType> __int = __x.intervals();
2598       __os << __int.size() - 1;
2599
2600       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2601         __os << __space << *__xit;
2602
2603       std::vector<double> __den = __x.densities();
2604       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2605         __os << __space << *__dit;
2606
2607       __os.flags(__flags);
2608       __os.fill(__fill);
2609       __os.precision(__precision);
2610       return __os;
2611     }
2612
2613   template<typename _RealType, typename _CharT, typename _Traits>
2614     std::basic_istream<_CharT, _Traits>&
2615     operator>>(std::basic_istream<_CharT, _Traits>& __is,
2616                piecewise_linear_distribution<_RealType>& __x)
2617     {
2618       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
2619       typedef typename __istream_type::ios_base    __ios_base;
2620
2621       const typename __ios_base::fmtflags __flags = __is.flags();
2622       __is.flags(__ios_base::dec | __ios_base::skipws);
2623
2624       size_t __n;
2625       __is >> __n;
2626
2627       std::vector<_RealType> __int_vec;
2628       __int_vec.reserve(__n + 1);
2629       for (size_t __i = 0; __i <= __n; ++__i)
2630         {
2631           _RealType __int;
2632           __is >> __int;
2633           __int_vec.push_back(__int);
2634         }
2635
2636       std::vector<double> __den_vec;
2637       __den_vec.reserve(__n + 1);
2638       for (size_t __i = 0; __i <= __n; ++__i)
2639         {
2640           double __den;
2641           __is >> __den;
2642           __den_vec.push_back(__den);
2643         }
2644
2645       __x.param(typename piecewise_linear_distribution<_RealType>::
2646           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
2647
2648       __is.flags(__flags);
2649       return __is;
2650     }
2651
2652
2653   template<typename _IntType>
2654     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
2655     {
2656       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
2657         _M_v.push_back(__detail::__mod<result_type,
2658                        __detail::_Shift<result_type, 32>::__value>(*__iter));
2659     }
2660
2661   template<typename _InputIterator>
2662     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
2663     {
2664       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
2665         _M_v.push_back(__detail::__mod<result_type,
2666                        __detail::_Shift<result_type, 32>::__value>(*__iter));
2667     }
2668
2669   template<typename _RandomAccessIterator>
2670     void
2671     seed_seq::generate(_RandomAccessIterator __begin,
2672                        _RandomAccessIterator __end)
2673     {
2674       typedef typename iterator_traits<_RandomAccessIterator>::value_type
2675         _Type;
2676
2677       if (__begin == __end)
2678         return;
2679
2680       std::fill(__begin, __end, _Type(0x8b8b8b8bu));
2681
2682       const size_t __n = __end - __begin;
2683       const size_t __s = _M_v.size();
2684       const size_t __t = (__n >= 623) ? 11
2685                        : (__n >=  68) ? 7
2686                        : (__n >=  39) ? 5
2687                        : (__n >=   7) ? 3
2688                        : (__n - 1) / 2;
2689       const size_t __p = (__n - __t) / 2;
2690       const size_t __q = __p + __t;
2691       const size_t __m = std::max(__s + 1, __n);
2692
2693       for (size_t __k = 0; __k < __m; ++__k)
2694         {
2695           _Type __arg = (__begin[__k % __n]
2696                          ^ __begin[(__k + __p) % __n]
2697                          ^ __begin[(__k - 1) % __n]);
2698           _Type __r1 = __arg ^ (__arg << 27);
2699           __r1 = __detail::__mod<_Type, __detail::_Shift<_Type, 32>::__value,
2700                                  1664525u, 0u>(__r1);
2701           _Type __r2 = __r1;
2702           if (__k == 0)
2703             __r2 += __s;
2704           else if (__k <= __s)
2705             __r2 += __k % __n + _M_v[__k - 1];
2706           else
2707             __r2 += __k % __n;
2708           __r2 = __detail::__mod<_Type,
2709                    __detail::_Shift<_Type, 32>::__value>(__r2);
2710           __begin[(__k + __p) % __n] += __r1;
2711           __begin[(__k + __q) % __n] += __r2;
2712           __begin[__k % __n] = __r2;
2713         }
2714
2715       for (size_t __k = __m; __k < __m + __n; ++__k)
2716         {
2717           _Type __arg = (__begin[__k % __n]
2718                          + __begin[(__k + __p) % __n]
2719                          + __begin[(__k - 1) % __n]);
2720           _Type __r3 = __arg ^ (__arg << 27);
2721           __r3 = __detail::__mod<_Type, __detail::_Shift<_Type, 32>::__value,
2722                                  1566083941u, 0u>(__r3);
2723           _Type __r4 = __r3 - __k % __n;
2724           __r4 = __detail::__mod<_Type,
2725                    __detail::_Shift<_Type, 32>::__value>(__r4);
2726           __begin[(__k + __p) % __n] ^= __r4;
2727           __begin[(__k + __q) % __n] ^= __r3;
2728           __begin[__k % __n] = __r4;
2729         }
2730     }
2731
2732   template<typename _RealType, size_t __bits,
2733            typename _UniformRandomNumberGenerator>
2734     _RealType
2735     generate_canonical(_UniformRandomNumberGenerator& __urng)
2736     {
2737       const size_t __b
2738         = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
2739                    __bits);
2740       const long double __r = static_cast<long double>(__urng.max())
2741                             - static_cast<long double>(__urng.min()) + 1.0L;
2742       const size_t __log2r = std::log(__r) / std::log(2.0L);
2743       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
2744       _RealType __sum = _RealType(0);
2745       _RealType __tmp = _RealType(1);
2746       for (; __k != 0; --__k)
2747         {
2748           __sum += _RealType(__urng() - __urng.min()) * __tmp;
2749           __tmp *= __r;
2750         }
2751       return __sum / __tmp;
2752     }
2753 }