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

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