libstdc++
random.tcc
Go to the documentation of this file.
00001 // random number generation (out of line) -*- C++ -*-
00002 
00003 // Copyright (C) 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
00004 //
00005 // This file is part of the GNU ISO C++ Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the
00007 // terms of the GNU General Public License as published by the
00008 // Free Software Foundation; either version 3, or (at your option)
00009 // any later version.
00010 
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 
00016 // Under Section 7 of GPL version 3, you are granted additional
00017 // permissions described in the GCC Runtime Library Exception, version
00018 // 3.1, as published by the Free Software Foundation.
00019 
00020 // You should have received a copy of the GNU General Public License and
00021 // a copy of the GCC Runtime Library Exception along with this program;
00022 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
00023 // <http://www.gnu.org/licenses/>.
00024 
00025 /** @file bits/random.tcc
00026  *  This is an internal header file, included by other library headers.
00027  *  Do not attempt to use it directly. @headername{random}
00028  */
00029 
00030 #ifndef _RANDOM_TCC
00031 #define _RANDOM_TCC 1
00032 
00033 #include <numeric> // std::accumulate and std::partial_sum
00034 
00035 namespace std _GLIBCXX_VISIBILITY(default)
00036 {
00037   /*
00038    * (Further) implementation-space details.
00039    */
00040   namespace __detail
00041   {
00042   _GLIBCXX_BEGIN_NAMESPACE_VERSION
00043 
00044     // General case for x = (ax + c) mod m -- use Schrage's algorithm to
00045     // avoid integer overflow.
00046     //
00047     // Because a and c are compile-time integral constants the compiler
00048     // kindly elides any unreachable paths.
00049     //
00050     // Preconditions:  a > 0, m > 0.
00051     //
00052     // XXX FIXME: as-is, only works correctly for __m % __a < __m / __a. 
00053     //
00054     template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool>
00055       struct _Mod
00056       {
00057     static _Tp
00058     __calc(_Tp __x)
00059     {
00060       if (__a == 1)
00061         __x %= __m;
00062       else
00063         {
00064           static const _Tp __q = __m / __a;
00065           static const _Tp __r = __m % __a;
00066 
00067           _Tp __t1 = __a * (__x % __q);
00068           _Tp __t2 = __r * (__x / __q);
00069           if (__t1 >= __t2)
00070         __x = __t1 - __t2;
00071           else
00072         __x = __m - __t2 + __t1;
00073         }
00074 
00075       if (__c != 0)
00076         {
00077           const _Tp __d = __m - __x;
00078           if (__d > __c)
00079         __x += __c;
00080           else
00081         __x = __c - __d;
00082         }
00083       return __x;
00084     }
00085       };
00086 
00087     // Special case for m == 0 -- use unsigned integer overflow as modulo
00088     // operator.
00089     template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
00090       struct _Mod<_Tp, __m, __a, __c, true>
00091       {
00092     static _Tp
00093     __calc(_Tp __x)
00094     { return __a * __x + __c; }
00095       };
00096 
00097     template<typename _InputIterator, typename _OutputIterator,
00098          typename _UnaryOperation>
00099       _OutputIterator
00100       __transform(_InputIterator __first, _InputIterator __last,
00101           _OutputIterator __result, _UnaryOperation __unary_op)
00102       {
00103     for (; __first != __last; ++__first, ++__result)
00104       *__result = __unary_op(*__first);
00105     return __result;
00106       }
00107 
00108   _GLIBCXX_END_NAMESPACE_VERSION
00109   } // namespace __detail
00110 
00111 _GLIBCXX_BEGIN_NAMESPACE_VERSION
00112 
00113   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00114     constexpr _UIntType
00115     linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
00116 
00117   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00118     constexpr _UIntType
00119     linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
00120 
00121   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00122     constexpr _UIntType
00123     linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
00124 
00125   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00126     constexpr _UIntType
00127     linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
00128 
00129   /**
00130    * Seeds the LCR with integral value @p __s, adjusted so that the
00131    * ring identity is never a member of the convergence set.
00132    */
00133   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00134     void
00135     linear_congruential_engine<_UIntType, __a, __c, __m>::
00136     seed(result_type __s)
00137     {
00138       if ((__detail::__mod<_UIntType, __m>(__c) == 0)
00139       && (__detail::__mod<_UIntType, __m>(__s) == 0))
00140     _M_x = 1;
00141       else
00142     _M_x = __detail::__mod<_UIntType, __m>(__s);
00143     }
00144 
00145   /**
00146    * Seeds the LCR engine with a value generated by @p __q.
00147    */
00148   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00149     template<typename _Sseq>
00150       typename std::enable_if<std::is_class<_Sseq>::value>::type
00151       linear_congruential_engine<_UIntType, __a, __c, __m>::
00152       seed(_Sseq& __q)
00153       {
00154     const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
00155                                     : std::__lg(__m);
00156     const _UIntType __k = (__k0 + 31) / 32;
00157     uint_least32_t __arr[__k + 3];
00158     __q.generate(__arr + 0, __arr + __k + 3);
00159     _UIntType __factor = 1u;
00160     _UIntType __sum = 0u;
00161     for (size_t __j = 0; __j < __k; ++__j)
00162       {
00163         __sum += __arr[__j + 3] * __factor;
00164         __factor *= __detail::_Shift<_UIntType, 32>::__value;
00165       }
00166     seed(__sum);
00167       }
00168 
00169   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
00170        typename _CharT, typename _Traits>
00171     std::basic_ostream<_CharT, _Traits>&
00172     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00173            const linear_congruential_engine<_UIntType,
00174                         __a, __c, __m>& __lcr)
00175     {
00176       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00177       typedef typename __ostream_type::ios_base    __ios_base;
00178 
00179       const typename __ios_base::fmtflags __flags = __os.flags();
00180       const _CharT __fill = __os.fill();
00181       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00182       __os.fill(__os.widen(' '));
00183 
00184       __os << __lcr._M_x;
00185 
00186       __os.flags(__flags);
00187       __os.fill(__fill);
00188       return __os;
00189     }
00190 
00191   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
00192        typename _CharT, typename _Traits>
00193     std::basic_istream<_CharT, _Traits>&
00194     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00195            linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
00196     {
00197       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00198       typedef typename __istream_type::ios_base    __ios_base;
00199 
00200       const typename __ios_base::fmtflags __flags = __is.flags();
00201       __is.flags(__ios_base::dec);
00202 
00203       __is >> __lcr._M_x;
00204 
00205       __is.flags(__flags);
00206       return __is;
00207     }
00208 
00209 
00210   template<typename _UIntType,
00211        size_t __w, size_t __n, size_t __m, size_t __r,
00212        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00213        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00214        _UIntType __f>
00215     constexpr size_t
00216     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00217                 __s, __b, __t, __c, __l, __f>::word_size;
00218 
00219   template<typename _UIntType,
00220        size_t __w, size_t __n, size_t __m, size_t __r,
00221        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00222        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00223        _UIntType __f>
00224     constexpr size_t
00225     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00226                 __s, __b, __t, __c, __l, __f>::state_size;
00227 
00228   template<typename _UIntType,
00229        size_t __w, size_t __n, size_t __m, size_t __r,
00230        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00231        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00232        _UIntType __f>
00233     constexpr size_t
00234     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00235                 __s, __b, __t, __c, __l, __f>::shift_size;
00236 
00237   template<typename _UIntType,
00238        size_t __w, size_t __n, size_t __m, size_t __r,
00239        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00240        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00241        _UIntType __f>
00242     constexpr size_t
00243     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00244                 __s, __b, __t, __c, __l, __f>::mask_bits;
00245 
00246   template<typename _UIntType,
00247        size_t __w, size_t __n, size_t __m, size_t __r,
00248        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00249        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00250        _UIntType __f>
00251     constexpr _UIntType
00252     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00253                 __s, __b, __t, __c, __l, __f>::xor_mask;
00254 
00255   template<typename _UIntType,
00256        size_t __w, size_t __n, size_t __m, size_t __r,
00257        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00258        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00259        _UIntType __f>
00260     constexpr size_t
00261     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00262                 __s, __b, __t, __c, __l, __f>::tempering_u;
00263    
00264   template<typename _UIntType,
00265        size_t __w, size_t __n, size_t __m, size_t __r,
00266        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00267        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00268        _UIntType __f>
00269     constexpr _UIntType
00270     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00271                 __s, __b, __t, __c, __l, __f>::tempering_d;
00272 
00273   template<typename _UIntType,
00274        size_t __w, size_t __n, size_t __m, size_t __r,
00275        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00276        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00277        _UIntType __f>
00278     constexpr size_t
00279     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00280                 __s, __b, __t, __c, __l, __f>::tempering_s;
00281 
00282   template<typename _UIntType,
00283        size_t __w, size_t __n, size_t __m, size_t __r,
00284        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00285        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00286        _UIntType __f>
00287     constexpr _UIntType
00288     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00289                 __s, __b, __t, __c, __l, __f>::tempering_b;
00290 
00291   template<typename _UIntType,
00292        size_t __w, size_t __n, size_t __m, size_t __r,
00293        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00294        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00295        _UIntType __f>
00296     constexpr size_t
00297     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00298                 __s, __b, __t, __c, __l, __f>::tempering_t;
00299 
00300   template<typename _UIntType,
00301        size_t __w, size_t __n, size_t __m, size_t __r,
00302        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00303        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00304        _UIntType __f>
00305     constexpr _UIntType
00306     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00307                 __s, __b, __t, __c, __l, __f>::tempering_c;
00308 
00309   template<typename _UIntType,
00310        size_t __w, size_t __n, size_t __m, size_t __r,
00311        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00312        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00313        _UIntType __f>
00314     constexpr size_t
00315     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00316                 __s, __b, __t, __c, __l, __f>::tempering_l;
00317 
00318   template<typename _UIntType,
00319        size_t __w, size_t __n, size_t __m, size_t __r,
00320        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00321        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00322        _UIntType __f>
00323     constexpr _UIntType
00324     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00325                 __s, __b, __t, __c, __l, __f>::
00326                                               initialization_multiplier;
00327 
00328   template<typename _UIntType,
00329        size_t __w, size_t __n, size_t __m, size_t __r,
00330        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00331        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00332        _UIntType __f>
00333     constexpr _UIntType
00334     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00335                 __s, __b, __t, __c, __l, __f>::default_seed;
00336 
00337   template<typename _UIntType,
00338        size_t __w, size_t __n, size_t __m, size_t __r,
00339        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00340        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00341        _UIntType __f>
00342     void
00343     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00344                 __s, __b, __t, __c, __l, __f>::
00345     seed(result_type __sd)
00346     {
00347       _M_x[0] = __detail::__mod<_UIntType,
00348     __detail::_Shift<_UIntType, __w>::__value>(__sd);
00349 
00350       for (size_t __i = 1; __i < state_size; ++__i)
00351     {
00352       _UIntType __x = _M_x[__i - 1];
00353       __x ^= __x >> (__w - 2);
00354       __x *= __f;
00355       __x += __detail::__mod<_UIntType, __n>(__i);
00356       _M_x[__i] = __detail::__mod<_UIntType,
00357         __detail::_Shift<_UIntType, __w>::__value>(__x);
00358     }
00359       _M_p = state_size;
00360     }
00361 
00362   template<typename _UIntType,
00363        size_t __w, size_t __n, size_t __m, size_t __r,
00364        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00365        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00366        _UIntType __f>
00367     template<typename _Sseq>
00368       typename std::enable_if<std::is_class<_Sseq>::value>::type
00369       mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00370                   __s, __b, __t, __c, __l, __f>::
00371       seed(_Sseq& __q)
00372       {
00373     const _UIntType __upper_mask = (~_UIntType()) << __r;
00374     const size_t __k = (__w + 31) / 32;
00375     uint_least32_t __arr[__n * __k];
00376     __q.generate(__arr + 0, __arr + __n * __k);
00377 
00378     bool __zero = true;
00379     for (size_t __i = 0; __i < state_size; ++__i)
00380       {
00381         _UIntType __factor = 1u;
00382         _UIntType __sum = 0u;
00383         for (size_t __j = 0; __j < __k; ++__j)
00384           {
00385         __sum += __arr[__k * __i + __j] * __factor;
00386         __factor *= __detail::_Shift<_UIntType, 32>::__value;
00387           }
00388         _M_x[__i] = __detail::__mod<_UIntType,
00389           __detail::_Shift<_UIntType, __w>::__value>(__sum);
00390 
00391         if (__zero)
00392           {
00393         if (__i == 0)
00394           {
00395             if ((_M_x[0] & __upper_mask) != 0u)
00396               __zero = false;
00397           }
00398         else if (_M_x[__i] != 0u)
00399           __zero = false;
00400           }
00401       }
00402         if (__zero)
00403           _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
00404       }
00405 
00406   template<typename _UIntType, size_t __w,
00407        size_t __n, size_t __m, size_t __r,
00408        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00409        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00410        _UIntType __f>
00411     typename
00412     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00413                 __s, __b, __t, __c, __l, __f>::result_type
00414     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00415                 __s, __b, __t, __c, __l, __f>::
00416     operator()()
00417     {
00418       // Reload the vector - cost is O(n) amortized over n calls.
00419       if (_M_p >= state_size)
00420     {
00421       const _UIntType __upper_mask = (~_UIntType()) << __r;
00422       const _UIntType __lower_mask = ~__upper_mask;
00423 
00424       for (size_t __k = 0; __k < (__n - __m); ++__k)
00425         {
00426           _UIntType __y = ((_M_x[__k] & __upper_mask)
00427                    | (_M_x[__k + 1] & __lower_mask));
00428           _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
00429                ^ ((__y & 0x01) ? __a : 0));
00430         }
00431 
00432       for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
00433         {
00434           _UIntType __y = ((_M_x[__k] & __upper_mask)
00435                    | (_M_x[__k + 1] & __lower_mask));
00436           _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
00437                ^ ((__y & 0x01) ? __a : 0));
00438         }
00439 
00440       _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
00441                | (_M_x[0] & __lower_mask));
00442       _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
00443                ^ ((__y & 0x01) ? __a : 0));
00444       _M_p = 0;
00445     }
00446 
00447       // Calculate o(x(i)).
00448       result_type __z = _M_x[_M_p++];
00449       __z ^= (__z >> __u) & __d;
00450       __z ^= (__z << __s) & __b;
00451       __z ^= (__z << __t) & __c;
00452       __z ^= (__z >> __l);
00453 
00454       return __z;
00455     }
00456 
00457   template<typename _UIntType, size_t __w,
00458        size_t __n, size_t __m, size_t __r,
00459        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00460        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00461        _UIntType __f, typename _CharT, typename _Traits>
00462     std::basic_ostream<_CharT, _Traits>&
00463     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00464            const mersenne_twister_engine<_UIntType, __w, __n, __m,
00465            __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
00466     {
00467       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00468       typedef typename __ostream_type::ios_base    __ios_base;
00469 
00470       const typename __ios_base::fmtflags __flags = __os.flags();
00471       const _CharT __fill = __os.fill();
00472       const _CharT __space = __os.widen(' ');
00473       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00474       __os.fill(__space);
00475 
00476       for (size_t __i = 0; __i < __n; ++__i)
00477     __os << __x._M_x[__i] << __space;
00478       __os << __x._M_p;
00479 
00480       __os.flags(__flags);
00481       __os.fill(__fill);
00482       return __os;
00483     }
00484 
00485   template<typename _UIntType, size_t __w,
00486        size_t __n, size_t __m, size_t __r,
00487        _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00488        _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00489        _UIntType __f, typename _CharT, typename _Traits>
00490     std::basic_istream<_CharT, _Traits>&
00491     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00492            mersenne_twister_engine<_UIntType, __w, __n, __m,
00493            __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
00494     {
00495       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00496       typedef typename __istream_type::ios_base    __ios_base;
00497 
00498       const typename __ios_base::fmtflags __flags = __is.flags();
00499       __is.flags(__ios_base::dec | __ios_base::skipws);
00500 
00501       for (size_t __i = 0; __i < __n; ++__i)
00502     __is >> __x._M_x[__i];
00503       __is >> __x._M_p;
00504 
00505       __is.flags(__flags);
00506       return __is;
00507     }
00508 
00509 
00510   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00511     constexpr size_t
00512     subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
00513 
00514   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00515     constexpr size_t
00516     subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
00517 
00518   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00519     constexpr size_t
00520     subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
00521 
00522   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00523     constexpr _UIntType
00524     subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
00525 
00526   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00527     void
00528     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00529     seed(result_type __value)
00530     {
00531       std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
00532     __lcg(__value == 0u ? default_seed : __value);
00533 
00534       const size_t __n = (__w + 31) / 32;
00535 
00536       for (size_t __i = 0; __i < long_lag; ++__i)
00537     {
00538       _UIntType __sum = 0u;
00539       _UIntType __factor = 1u;
00540       for (size_t __j = 0; __j < __n; ++__j)
00541         {
00542           __sum += __detail::__mod<uint_least32_t,
00543                __detail::_Shift<uint_least32_t, 32>::__value>
00544              (__lcg()) * __factor;
00545           __factor *= __detail::_Shift<_UIntType, 32>::__value;
00546         }
00547       _M_x[__i] = __detail::__mod<_UIntType,
00548         __detail::_Shift<_UIntType, __w>::__value>(__sum);
00549     }
00550       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
00551       _M_p = 0;
00552     }
00553 
00554   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00555     template<typename _Sseq>
00556       typename std::enable_if<std::is_class<_Sseq>::value>::type
00557       subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00558       seed(_Sseq& __q)
00559       {
00560     const size_t __k = (__w + 31) / 32;
00561     uint_least32_t __arr[__r * __k];
00562     __q.generate(__arr + 0, __arr + __r * __k);
00563 
00564     for (size_t __i = 0; __i < long_lag; ++__i)
00565       {
00566         _UIntType __sum = 0u;
00567         _UIntType __factor = 1u;
00568         for (size_t __j = 0; __j < __k; ++__j)
00569           {
00570         __sum += __arr[__k * __i + __j] * __factor;
00571         __factor *= __detail::_Shift<_UIntType, 32>::__value;
00572           }
00573         _M_x[__i] = __detail::__mod<_UIntType,
00574           __detail::_Shift<_UIntType, __w>::__value>(__sum);
00575       }
00576     _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
00577     _M_p = 0;
00578       }
00579 
00580   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00581     typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00582          result_type
00583     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00584     operator()()
00585     {
00586       // Derive short lag index from current index.
00587       long __ps = _M_p - short_lag;
00588       if (__ps < 0)
00589     __ps += long_lag;
00590 
00591       // Calculate new x(i) without overflow or division.
00592       // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
00593       // cannot overflow.
00594       _UIntType __xi;
00595       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
00596     {
00597       __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
00598       _M_carry = 0;
00599     }
00600       else
00601     {
00602       __xi = (__detail::_Shift<_UIntType, __w>::__value
00603           - _M_x[_M_p] - _M_carry + _M_x[__ps]);
00604       _M_carry = 1;
00605     }
00606       _M_x[_M_p] = __xi;
00607 
00608       // Adjust current index to loop around in ring buffer.
00609       if (++_M_p >= long_lag)
00610     _M_p = 0;
00611 
00612       return __xi;
00613     }
00614 
00615   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
00616        typename _CharT, typename _Traits>
00617     std::basic_ostream<_CharT, _Traits>&
00618     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00619            const subtract_with_carry_engine<_UIntType,
00620                         __w, __s, __r>& __x)
00621     {
00622       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00623       typedef typename __ostream_type::ios_base    __ios_base;
00624 
00625       const typename __ios_base::fmtflags __flags = __os.flags();
00626       const _CharT __fill = __os.fill();
00627       const _CharT __space = __os.widen(' ');
00628       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00629       __os.fill(__space);
00630 
00631       for (size_t __i = 0; __i < __r; ++__i)
00632     __os << __x._M_x[__i] << __space;
00633       __os << __x._M_carry << __space << __x._M_p;
00634 
00635       __os.flags(__flags);
00636       __os.fill(__fill);
00637       return __os;
00638     }
00639 
00640   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
00641        typename _CharT, typename _Traits>
00642     std::basic_istream<_CharT, _Traits>&
00643     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00644            subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
00645     {
00646       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
00647       typedef typename __istream_type::ios_base    __ios_base;
00648 
00649       const typename __ios_base::fmtflags __flags = __is.flags();
00650       __is.flags(__ios_base::dec | __ios_base::skipws);
00651 
00652       for (size_t __i = 0; __i < __r; ++__i)
00653     __is >> __x._M_x[__i];
00654       __is >> __x._M_carry;
00655       __is >> __x._M_p;
00656 
00657       __is.flags(__flags);
00658       return __is;
00659     }
00660 
00661 
00662   template<typename _RandomNumberEngine, size_t __p, size_t __r>
00663     constexpr size_t
00664     discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
00665 
00666   template<typename _RandomNumberEngine, size_t __p, size_t __r>
00667     constexpr size_t
00668     discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
00669 
00670   template<typename _RandomNumberEngine, size_t __p, size_t __r>
00671     typename discard_block_engine<_RandomNumberEngine,
00672                __p, __r>::result_type
00673     discard_block_engine<_RandomNumberEngine, __p, __r>::
00674     operator()()
00675     {
00676       if (_M_n >= used_block)
00677     {
00678       _M_b.discard(block_size - _M_n);
00679       _M_n = 0;
00680     }
00681       ++_M_n;
00682       return _M_b();
00683     }
00684 
00685   template<typename _RandomNumberEngine, size_t __p, size_t __r,
00686        typename _CharT, typename _Traits>
00687     std::basic_ostream<_CharT, _Traits>&
00688     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00689            const discard_block_engine<_RandomNumberEngine,
00690            __p, __r>& __x)
00691     {
00692       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00693       typedef typename __ostream_type::ios_base    __ios_base;
00694 
00695       const typename __ios_base::fmtflags __flags = __os.flags();
00696       const _CharT __fill = __os.fill();
00697       const _CharT __space = __os.widen(' ');
00698       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00699       __os.fill(__space);
00700 
00701       __os << __x.base() << __space << __x._M_n;
00702 
00703       __os.flags(__flags);
00704       __os.fill(__fill);
00705       return __os;
00706     }
00707 
00708   template<typename _RandomNumberEngine, size_t __p, size_t __r,
00709        typename _CharT, typename _Traits>
00710     std::basic_istream<_CharT, _Traits>&
00711     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00712            discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
00713     {
00714       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00715       typedef typename __istream_type::ios_base    __ios_base;
00716 
00717       const typename __ios_base::fmtflags __flags = __is.flags();
00718       __is.flags(__ios_base::dec | __ios_base::skipws);
00719 
00720       __is >> __x._M_b >> __x._M_n;
00721 
00722       __is.flags(__flags);
00723       return __is;
00724     }
00725 
00726 
00727   template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
00728     typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
00729       result_type
00730     independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
00731     operator()()
00732     {
00733       typedef typename _RandomNumberEngine::result_type _Eresult_type;
00734       const _Eresult_type __r
00735     = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
00736        ? _M_b.max() - _M_b.min() + 1 : 0);
00737       const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
00738       const unsigned __m = __r ? std::__lg(__r) : __edig;
00739 
00740       typedef typename std::common_type<_Eresult_type, result_type>::type
00741     __ctype;
00742       const unsigned __cdig = std::numeric_limits<__ctype>::digits;
00743 
00744       unsigned __n, __n0;
00745       __ctype __s0, __s1, __y0, __y1;
00746 
00747       for (size_t __i = 0; __i < 2; ++__i)
00748     {
00749       __n = (__w + __m - 1) / __m + __i;
00750       __n0 = __n - __w % __n;
00751       const unsigned __w0 = __w / __n;  // __w0 <= __m
00752 
00753       __s0 = 0;
00754       __s1 = 0;
00755       if (__w0 < __cdig)
00756         {
00757           __s0 = __ctype(1) << __w0;
00758           __s1 = __s0 << 1;
00759         }
00760 
00761       __y0 = 0;
00762       __y1 = 0;
00763       if (__r)
00764         {
00765           __y0 = __s0 * (__r / __s0);
00766           if (__s1)
00767         __y1 = __s1 * (__r / __s1);
00768 
00769           if (__r - __y0 <= __y0 / __n)
00770         break;
00771         }
00772       else
00773         break;
00774     }
00775 
00776       result_type __sum = 0;
00777       for (size_t __k = 0; __k < __n0; ++__k)
00778     {
00779       __ctype __u;
00780       do
00781         __u = _M_b() - _M_b.min();
00782       while (__y0 && __u >= __y0);
00783       __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
00784     }
00785       for (size_t __k = __n0; __k < __n; ++__k)
00786     {
00787       __ctype __u;
00788       do
00789         __u = _M_b() - _M_b.min();
00790       while (__y1 && __u >= __y1);
00791       __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
00792     }
00793       return __sum;
00794     }
00795 
00796 
00797   template<typename _RandomNumberEngine, size_t __k>
00798     constexpr size_t
00799     shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
00800 
00801   template<typename _RandomNumberEngine, size_t __k>
00802     typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
00803     shuffle_order_engine<_RandomNumberEngine, __k>::
00804     operator()()
00805     {
00806       size_t __j = __k * ((_M_y - _M_b.min())
00807               / (_M_b.max() - _M_b.min() + 1.0L));
00808       _M_y = _M_v[__j];
00809       _M_v[__j] = _M_b();
00810 
00811       return _M_y;
00812     }
00813 
00814   template<typename _RandomNumberEngine, size_t __k,
00815        typename _CharT, typename _Traits>
00816     std::basic_ostream<_CharT, _Traits>&
00817     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00818            const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
00819     {
00820       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00821       typedef typename __ostream_type::ios_base    __ios_base;
00822 
00823       const typename __ios_base::fmtflags __flags = __os.flags();
00824       const _CharT __fill = __os.fill();
00825       const _CharT __space = __os.widen(' ');
00826       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00827       __os.fill(__space);
00828 
00829       __os << __x.base();
00830       for (size_t __i = 0; __i < __k; ++__i)
00831     __os << __space << __x._M_v[__i];
00832       __os << __space << __x._M_y;
00833 
00834       __os.flags(__flags);
00835       __os.fill(__fill);
00836       return __os;
00837     }
00838 
00839   template<typename _RandomNumberEngine, size_t __k,
00840        typename _CharT, typename _Traits>
00841     std::basic_istream<_CharT, _Traits>&
00842     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00843            shuffle_order_engine<_RandomNumberEngine, __k>& __x)
00844     {
00845       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00846       typedef typename __istream_type::ios_base    __ios_base;
00847 
00848       const typename __ios_base::fmtflags __flags = __is.flags();
00849       __is.flags(__ios_base::dec | __ios_base::skipws);
00850 
00851       __is >> __x._M_b;
00852       for (size_t __i = 0; __i < __k; ++__i)
00853     __is >> __x._M_v[__i];
00854       __is >> __x._M_y;
00855 
00856       __is.flags(__flags);
00857       return __is;
00858     }
00859 
00860 
00861   template<typename _IntType>
00862     template<typename _UniformRandomNumberGenerator>
00863       typename uniform_int_distribution<_IntType>::result_type
00864       uniform_int_distribution<_IntType>::
00865       operator()(_UniformRandomNumberGenerator& __urng,
00866          const param_type& __param)
00867       {
00868     typedef typename _UniformRandomNumberGenerator::result_type
00869       _Gresult_type;
00870     typedef typename std::make_unsigned<result_type>::type __utype;
00871     typedef typename std::common_type<_Gresult_type, __utype>::type
00872       __uctype;
00873 
00874     const __uctype __urngmin = __urng.min();
00875     const __uctype __urngmax = __urng.max();
00876     const __uctype __urngrange = __urngmax - __urngmin;
00877     const __uctype __urange
00878       = __uctype(__param.b()) - __uctype(__param.a());
00879 
00880     __uctype __ret;
00881 
00882     if (__urngrange > __urange)
00883       {
00884         // downscaling
00885         const __uctype __uerange = __urange + 1; // __urange can be zero
00886         const __uctype __scaling = __urngrange / __uerange;
00887         const __uctype __past = __uerange * __scaling;
00888         do
00889           __ret = __uctype(__urng()) - __urngmin;
00890         while (__ret >= __past);
00891         __ret /= __scaling;
00892       }
00893     else if (__urngrange < __urange)
00894       {
00895         // upscaling
00896         /*
00897           Note that every value in [0, urange]
00898           can be written uniquely as
00899 
00900           (urngrange + 1) * high + low
00901 
00902           where
00903 
00904           high in [0, urange / (urngrange + 1)]
00905 
00906           and
00907     
00908           low in [0, urngrange].
00909         */
00910         __uctype __tmp; // wraparound control
00911         do
00912           {
00913         const __uctype __uerngrange = __urngrange + 1;
00914         __tmp = (__uerngrange * operator()
00915              (__urng, param_type(0, __urange / __uerngrange)));
00916         __ret = __tmp + (__uctype(__urng()) - __urngmin);
00917           }
00918         while (__ret > __urange || __ret < __tmp);
00919       }
00920     else
00921       __ret = __uctype(__urng()) - __urngmin;
00922 
00923     return __ret + __param.a();
00924       }
00925 
00926   template<typename _IntType, typename _CharT, typename _Traits>
00927     std::basic_ostream<_CharT, _Traits>&
00928     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00929            const uniform_int_distribution<_IntType>& __x)
00930     {
00931       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00932       typedef typename __ostream_type::ios_base    __ios_base;
00933 
00934       const typename __ios_base::fmtflags __flags = __os.flags();
00935       const _CharT __fill = __os.fill();
00936       const _CharT __space = __os.widen(' ');
00937       __os.flags(__ios_base::scientific | __ios_base::left);
00938       __os.fill(__space);
00939 
00940       __os << __x.a() << __space << __x.b();
00941 
00942       __os.flags(__flags);
00943       __os.fill(__fill);
00944       return __os;
00945     }
00946 
00947   template<typename _IntType, typename _CharT, typename _Traits>
00948     std::basic_istream<_CharT, _Traits>&
00949     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00950            uniform_int_distribution<_IntType>& __x)
00951     {
00952       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00953       typedef typename __istream_type::ios_base    __ios_base;
00954 
00955       const typename __ios_base::fmtflags __flags = __is.flags();
00956       __is.flags(__ios_base::dec | __ios_base::skipws);
00957 
00958       _IntType __a, __b;
00959       __is >> __a >> __b;
00960       __x.param(typename uniform_int_distribution<_IntType>::
00961         param_type(__a, __b));
00962 
00963       __is.flags(__flags);
00964       return __is;
00965     }
00966 
00967 
00968   template<typename _RealType, typename _CharT, typename _Traits>
00969     std::basic_ostream<_CharT, _Traits>&
00970     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00971            const uniform_real_distribution<_RealType>& __x)
00972     {
00973       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00974       typedef typename __ostream_type::ios_base    __ios_base;
00975 
00976       const typename __ios_base::fmtflags __flags = __os.flags();
00977       const _CharT __fill = __os.fill();
00978       const std::streamsize __precision = __os.precision();
00979       const _CharT __space = __os.widen(' ');
00980       __os.flags(__ios_base::scientific | __ios_base::left);
00981       __os.fill(__space);
00982       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00983 
00984       __os << __x.a() << __space << __x.b();
00985 
00986       __os.flags(__flags);
00987       __os.fill(__fill);
00988       __os.precision(__precision);
00989       return __os;
00990     }
00991 
00992   template<typename _RealType, typename _CharT, typename _Traits>
00993     std::basic_istream<_CharT, _Traits>&
00994     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00995            uniform_real_distribution<_RealType>& __x)
00996     {
00997       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00998       typedef typename __istream_type::ios_base    __ios_base;
00999 
01000       const typename __ios_base::fmtflags __flags = __is.flags();
01001       __is.flags(__ios_base::skipws);
01002 
01003       _RealType __a, __b;
01004       __is >> __a >> __b;
01005       __x.param(typename uniform_real_distribution<_RealType>::
01006         param_type(__a, __b));
01007 
01008       __is.flags(__flags);
01009       return __is;
01010     }
01011 
01012 
01013   template<typename _CharT, typename _Traits>
01014     std::basic_ostream<_CharT, _Traits>&
01015     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01016            const bernoulli_distribution& __x)
01017     {
01018       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01019       typedef typename __ostream_type::ios_base    __ios_base;
01020 
01021       const typename __ios_base::fmtflags __flags = __os.flags();
01022       const _CharT __fill = __os.fill();
01023       const std::streamsize __precision = __os.precision();
01024       __os.flags(__ios_base::scientific | __ios_base::left);
01025       __os.fill(__os.widen(' '));
01026       __os.precision(std::numeric_limits<double>::max_digits10);
01027 
01028       __os << __x.p();
01029 
01030       __os.flags(__flags);
01031       __os.fill(__fill);
01032       __os.precision(__precision);
01033       return __os;
01034     }
01035 
01036 
01037   template<typename _IntType>
01038     template<typename _UniformRandomNumberGenerator>
01039       typename geometric_distribution<_IntType>::result_type
01040       geometric_distribution<_IntType>::
01041       operator()(_UniformRandomNumberGenerator& __urng,
01042          const param_type& __param)
01043       {
01044     // About the epsilon thing see this thread:
01045     // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
01046     const double __naf =
01047       (1 - std::numeric_limits<double>::epsilon()) / 2;
01048     // The largest _RealType convertible to _IntType.
01049     const double __thr =
01050       std::numeric_limits<_IntType>::max() + __naf;
01051     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01052       __aurng(__urng);
01053 
01054     double __cand;
01055     do
01056       __cand = std::floor(std::log(__aurng()) / __param._M_log_1_p);
01057     while (__cand >= __thr);
01058 
01059     return result_type(__cand + __naf);
01060       }
01061 
01062   template<typename _IntType,
01063        typename _CharT, typename _Traits>
01064     std::basic_ostream<_CharT, _Traits>&
01065     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01066            const geometric_distribution<_IntType>& __x)
01067     {
01068       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01069       typedef typename __ostream_type::ios_base    __ios_base;
01070 
01071       const typename __ios_base::fmtflags __flags = __os.flags();
01072       const _CharT __fill = __os.fill();
01073       const std::streamsize __precision = __os.precision();
01074       __os.flags(__ios_base::scientific | __ios_base::left);
01075       __os.fill(__os.widen(' '));
01076       __os.precision(std::numeric_limits<double>::max_digits10);
01077 
01078       __os << __x.p();
01079 
01080       __os.flags(__flags);
01081       __os.fill(__fill);
01082       __os.precision(__precision);
01083       return __os;
01084     }
01085 
01086   template<typename _IntType,
01087        typename _CharT, typename _Traits>
01088     std::basic_istream<_CharT, _Traits>&
01089     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01090            geometric_distribution<_IntType>& __x)
01091     {
01092       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01093       typedef typename __istream_type::ios_base    __ios_base;
01094 
01095       const typename __ios_base::fmtflags __flags = __is.flags();
01096       __is.flags(__ios_base::skipws);
01097 
01098       double __p;
01099       __is >> __p;
01100       __x.param(typename geometric_distribution<_IntType>::param_type(__p));
01101 
01102       __is.flags(__flags);
01103       return __is;
01104     }
01105 
01106   // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
01107   template<typename _IntType>
01108     template<typename _UniformRandomNumberGenerator>
01109       typename negative_binomial_distribution<_IntType>::result_type
01110       negative_binomial_distribution<_IntType>::
01111       operator()(_UniformRandomNumberGenerator& __urng)
01112       {
01113     const double __y = _M_gd(__urng);
01114 
01115     // XXX Is the constructor too slow?
01116     std::poisson_distribution<result_type> __poisson(__y);
01117     return __poisson(__urng);
01118       }
01119 
01120   template<typename _IntType>
01121     template<typename _UniformRandomNumberGenerator>
01122       typename negative_binomial_distribution<_IntType>::result_type
01123       negative_binomial_distribution<_IntType>::
01124       operator()(_UniformRandomNumberGenerator& __urng,
01125          const param_type& __p)
01126       {
01127     typedef typename std::gamma_distribution<result_type>::param_type
01128       param_type;
01129     
01130     const double __y =
01131       _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
01132 
01133     std::poisson_distribution<result_type> __poisson(__y);
01134     return __poisson(__urng);
01135       }
01136 
01137   template<typename _IntType, typename _CharT, typename _Traits>
01138     std::basic_ostream<_CharT, _Traits>&
01139     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01140            const negative_binomial_distribution<_IntType>& __x)
01141     {
01142       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01143       typedef typename __ostream_type::ios_base    __ios_base;
01144 
01145       const typename __ios_base::fmtflags __flags = __os.flags();
01146       const _CharT __fill = __os.fill();
01147       const std::streamsize __precision = __os.precision();
01148       const _CharT __space = __os.widen(' ');
01149       __os.flags(__ios_base::scientific | __ios_base::left);
01150       __os.fill(__os.widen(' '));
01151       __os.precision(std::numeric_limits<double>::max_digits10);
01152 
01153       __os << __x.k() << __space << __x.p()
01154        << __space << __x._M_gd;
01155 
01156       __os.flags(__flags);
01157       __os.fill(__fill);
01158       __os.precision(__precision);
01159       return __os;
01160     }
01161 
01162   template<typename _IntType, typename _CharT, typename _Traits>
01163     std::basic_istream<_CharT, _Traits>&
01164     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01165            negative_binomial_distribution<_IntType>& __x)
01166     {
01167       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01168       typedef typename __istream_type::ios_base    __ios_base;
01169 
01170       const typename __ios_base::fmtflags __flags = __is.flags();
01171       __is.flags(__ios_base::skipws);
01172 
01173       _IntType __k;
01174       double __p;
01175       __is >> __k >> __p >> __x._M_gd;
01176       __x.param(typename negative_binomial_distribution<_IntType>::
01177         param_type(__k, __p));
01178 
01179       __is.flags(__flags);
01180       return __is;
01181     }
01182 
01183 
01184   template<typename _IntType>
01185     void
01186     poisson_distribution<_IntType>::param_type::
01187     _M_initialize()
01188     {
01189 #if _GLIBCXX_USE_C99_MATH_TR1
01190       if (_M_mean >= 12)
01191     {
01192       const double __m = std::floor(_M_mean);
01193       _M_lm_thr = std::log(_M_mean);
01194       _M_lfm = std::lgamma(__m + 1);
01195       _M_sm = std::sqrt(__m);
01196 
01197       const double __pi_4 = 0.7853981633974483096156608458198757L;
01198       const double __dx = std::sqrt(2 * __m * std::log(32 * __m
01199                                   / __pi_4));
01200       _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
01201       const double __cx = 2 * __m + _M_d;
01202       _M_scx = std::sqrt(__cx / 2);
01203       _M_1cx = 1 / __cx;
01204 
01205       _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
01206       _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
01207         / _M_d;
01208     }
01209       else
01210 #endif
01211     _M_lm_thr = std::exp(-_M_mean);
01212       }
01213 
01214   /**
01215    * A rejection algorithm when mean >= 12 and a simple method based
01216    * upon the multiplication of uniform random variates otherwise.
01217    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
01218    * is defined.
01219    *
01220    * Reference:
01221    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01222    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
01223    */
01224   template<typename _IntType>
01225     template<typename _UniformRandomNumberGenerator>
01226       typename poisson_distribution<_IntType>::result_type
01227       poisson_distribution<_IntType>::
01228       operator()(_UniformRandomNumberGenerator& __urng,
01229          const param_type& __param)
01230       {
01231     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01232       __aurng(__urng);
01233 #if _GLIBCXX_USE_C99_MATH_TR1
01234     if (__param.mean() >= 12)
01235       {
01236         double __x;
01237 
01238         // See comments above...
01239         const double __naf =
01240           (1 - std::numeric_limits<double>::epsilon()) / 2;
01241         const double __thr =
01242           std::numeric_limits<_IntType>::max() + __naf;
01243 
01244         const double __m = std::floor(__param.mean());
01245         // sqrt(pi / 2)
01246         const double __spi_2 = 1.2533141373155002512078826424055226L;
01247         const double __c1 = __param._M_sm * __spi_2;
01248         const double __c2 = __param._M_c2b + __c1;
01249         const double __c3 = __c2 + 1;
01250         const double __c4 = __c3 + 1;
01251         // e^(1 / 78)
01252         const double __e178 = 1.0129030479320018583185514777512983L;
01253         const double __c5 = __c4 + __e178;
01254         const double __c = __param._M_cb + __c5;
01255         const double __2cx = 2 * (2 * __m + __param._M_d);
01256 
01257         bool __reject = true;
01258         do
01259           {
01260         const double __u = __c * __aurng();
01261         const double __e = -std::log(__aurng());
01262 
01263         double __w = 0.0;
01264 
01265         if (__u <= __c1)
01266           {
01267             const double __n = _M_nd(__urng);
01268             const double __y = -std::abs(__n) * __param._M_sm - 1;
01269             __x = std::floor(__y);
01270             __w = -__n * __n / 2;
01271             if (__x < -__m)
01272               continue;
01273           }
01274         else if (__u <= __c2)
01275           {
01276             const double __n = _M_nd(__urng);
01277             const double __y = 1 + std::abs(__n) * __param._M_scx;
01278             __x = std::ceil(__y);
01279             __w = __y * (2 - __y) * __param._M_1cx;
01280             if (__x > __param._M_d)
01281               continue;
01282           }
01283         else if (__u <= __c3)
01284           // NB: This case not in the book, nor in the Errata,
01285           // but should be ok...
01286           __x = -1;
01287         else if (__u <= __c4)
01288           __x = 0;
01289         else if (__u <= __c5)
01290           __x = 1;
01291         else
01292           {
01293             const double __v = -std::log(__aurng());
01294             const double __y = __param._M_d
01295                      + __v * __2cx / __param._M_d;
01296             __x = std::ceil(__y);
01297             __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
01298           }
01299 
01300         __reject = (__w - __e - __x * __param._M_lm_thr
01301                 > __param._M_lfm - std::lgamma(__x + __m + 1));
01302 
01303         __reject |= __x + __m >= __thr;
01304 
01305           } while (__reject);
01306 
01307         return result_type(__x + __m + __naf);
01308       }
01309     else
01310 #endif
01311       {
01312         _IntType     __x = 0;
01313         double __prod = 1.0;
01314 
01315         do
01316           {
01317         __prod *= __aurng();
01318         __x += 1;
01319           }
01320         while (__prod > __param._M_lm_thr);
01321 
01322         return __x - 1;
01323       }
01324       }
01325 
01326   template<typename _IntType,
01327        typename _CharT, typename _Traits>
01328     std::basic_ostream<_CharT, _Traits>&
01329     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01330            const poisson_distribution<_IntType>& __x)
01331     {
01332       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01333       typedef typename __ostream_type::ios_base    __ios_base;
01334 
01335       const typename __ios_base::fmtflags __flags = __os.flags();
01336       const _CharT __fill = __os.fill();
01337       const std::streamsize __precision = __os.precision();
01338       const _CharT __space = __os.widen(' ');
01339       __os.flags(__ios_base::scientific | __ios_base::left);
01340       __os.fill(__space);
01341       __os.precision(std::numeric_limits<double>::max_digits10);
01342 
01343       __os << __x.mean() << __space << __x._M_nd;
01344 
01345       __os.flags(__flags);
01346       __os.fill(__fill);
01347       __os.precision(__precision);
01348       return __os;
01349     }
01350 
01351   template<typename _IntType,
01352        typename _CharT, typename _Traits>
01353     std::basic_istream<_CharT, _Traits>&
01354     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01355            poisson_distribution<_IntType>& __x)
01356     {
01357       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01358       typedef typename __istream_type::ios_base    __ios_base;
01359 
01360       const typename __ios_base::fmtflags __flags = __is.flags();
01361       __is.flags(__ios_base::skipws);
01362 
01363       double __mean;
01364       __is >> __mean >> __x._M_nd;
01365       __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
01366 
01367       __is.flags(__flags);
01368       return __is;
01369     }
01370 
01371 
01372   template<typename _IntType>
01373     void
01374     binomial_distribution<_IntType>::param_type::
01375     _M_initialize()
01376     {
01377       const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
01378 
01379       _M_easy = true;
01380 
01381 #if _GLIBCXX_USE_C99_MATH_TR1
01382       if (_M_t * __p12 >= 8)
01383     {
01384       _M_easy = false;
01385       const double __np = std::floor(_M_t * __p12);
01386       const double __pa = __np / _M_t;
01387       const double __1p = 1 - __pa;
01388 
01389       const double __pi_4 = 0.7853981633974483096156608458198757L;
01390       const double __d1x =
01391         std::sqrt(__np * __1p * std::log(32 * __np
01392                          / (81 * __pi_4 * __1p)));
01393       _M_d1 = std::round(std::max(1.0, __d1x));
01394       const double __d2x =
01395         std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
01396                          / (__pi_4 * __pa)));
01397       _M_d2 = std::round(std::max(1.0, __d2x));
01398 
01399       // sqrt(pi / 2)
01400       const double __spi_2 = 1.2533141373155002512078826424055226L;
01401       _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
01402       _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
01403       _M_c = 2 * _M_d1 / __np;
01404       _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
01405       const double __a12 = _M_a1 + _M_s2 * __spi_2;
01406       const double __s1s = _M_s1 * _M_s1;
01407       _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
01408                  * 2 * __s1s / _M_d1
01409                  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
01410       const double __s2s = _M_s2 * _M_s2;
01411       _M_s = (_M_a123 + 2 * __s2s / _M_d2
01412           * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
01413       _M_lf = (std::lgamma(__np + 1)
01414            + std::lgamma(_M_t - __np + 1));
01415       _M_lp1p = std::log(__pa / __1p);
01416 
01417       _M_q = -std::log(1 - (__p12 - __pa) / __1p);
01418     }
01419       else
01420 #endif
01421     _M_q = -std::log(1 - __p12);
01422     }
01423 
01424   template<typename _IntType>
01425     template<typename _UniformRandomNumberGenerator>
01426       typename binomial_distribution<_IntType>::result_type
01427       binomial_distribution<_IntType>::
01428       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
01429       {
01430     _IntType __x = 0;
01431     double __sum = 0.0;
01432     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01433       __aurng(__urng);
01434 
01435     do
01436       {
01437         const double __e = -std::log(__aurng());
01438         __sum += __e / (__t - __x);
01439         __x += 1;
01440       }
01441     while (__sum <= _M_param._M_q);
01442 
01443     return __x - 1;
01444       }
01445 
01446   /**
01447    * A rejection algorithm when t * p >= 8 and a simple waiting time
01448    * method - the second in the referenced book - otherwise.
01449    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
01450    * is defined.
01451    *
01452    * Reference:
01453    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01454    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
01455    */
01456   template<typename _IntType>
01457     template<typename _UniformRandomNumberGenerator>
01458       typename binomial_distribution<_IntType>::result_type
01459       binomial_distribution<_IntType>::
01460       operator()(_UniformRandomNumberGenerator& __urng,
01461          const param_type& __param)
01462       {
01463     result_type __ret;
01464     const _IntType __t = __param.t();
01465     const double __p = __param.p();
01466     const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
01467     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01468       __aurng(__urng);
01469 
01470 #if _GLIBCXX_USE_C99_MATH_TR1
01471     if (!__param._M_easy)
01472       {
01473         double __x;
01474 
01475         // See comments above...
01476         const double __naf =
01477           (1 - std::numeric_limits<double>::epsilon()) / 2;
01478         const double __thr =
01479           std::numeric_limits<_IntType>::max() + __naf;
01480 
01481         const double __np = std::floor(__t * __p12);
01482 
01483         // sqrt(pi / 2)
01484         const double __spi_2 = 1.2533141373155002512078826424055226L;
01485         const double __a1 = __param._M_a1;
01486         const double __a12 = __a1 + __param._M_s2 * __spi_2;
01487         const double __a123 = __param._M_a123;
01488         const double __s1s = __param._M_s1 * __param._M_s1;
01489         const double __s2s = __param._M_s2 * __param._M_s2;
01490 
01491         bool __reject;
01492         do
01493           {
01494         const double __u = __param._M_s * __aurng();
01495 
01496         double __v;
01497 
01498         if (__u <= __a1)
01499           {
01500             const double __n = _M_nd(__urng);
01501             const double __y = __param._M_s1 * std::abs(__n);
01502             __reject = __y >= __param._M_d1;
01503             if (!__reject)
01504               {
01505             const double __e = -std::log(__aurng());
01506             __x = std::floor(__y);
01507             __v = -__e - __n * __n / 2 + __param._M_c;
01508               }
01509           }
01510         else if (__u <= __a12)
01511           {
01512             const double __n = _M_nd(__urng);
01513             const double __y = __param._M_s2 * std::abs(__n);
01514             __reject = __y >= __param._M_d2;
01515             if (!__reject)
01516               {
01517             const double __e = -std::log(__aurng());
01518             __x = std::floor(-__y);
01519             __v = -__e - __n * __n / 2;
01520               }
01521           }
01522         else if (__u <= __a123)
01523           {
01524             const double __e1 = -std::log(__aurng());
01525             const double __e2 = -std::log(__aurng());
01526 
01527             const double __y = __param._M_d1
01528                      + 2 * __s1s * __e1 / __param._M_d1;
01529             __x = std::floor(__y);
01530             __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
01531                             -__y / (2 * __s1s)));
01532             __reject = false;
01533           }
01534         else
01535           {
01536             const double __e1 = -std::log(__aurng());
01537             const double __e2 = -std::log(__aurng());
01538 
01539             const double __y = __param._M_d2
01540                      + 2 * __s2s * __e1 / __param._M_d2;
01541             __x = std::floor(-__y);
01542             __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
01543             __reject = false;
01544           }
01545 
01546         __reject = __reject || __x < -__np || __x > __t - __np;
01547         if (!__reject)
01548           {
01549             const double __lfx =
01550               std::lgamma(__np + __x + 1)
01551               + std::lgamma(__t - (__np + __x) + 1);
01552             __reject = __v > __param._M_lf - __lfx
01553                  + __x * __param._M_lp1p;
01554           }
01555 
01556         __reject |= __x + __np >= __thr;
01557           }
01558         while (__reject);
01559 
01560         __x += __np + __naf;
01561 
01562         const _IntType __z = _M_waiting(__urng, __t - _IntType(__x));
01563         __ret = _IntType(__x) + __z;
01564       }
01565     else
01566 #endif
01567       __ret = _M_waiting(__urng, __t);
01568 
01569     if (__p12 != __p)
01570       __ret = __t - __ret;
01571     return __ret;
01572       }
01573 
01574   template<typename _IntType,
01575        typename _CharT, typename _Traits>
01576     std::basic_ostream<_CharT, _Traits>&
01577     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01578            const binomial_distribution<_IntType>& __x)
01579     {
01580       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01581       typedef typename __ostream_type::ios_base    __ios_base;
01582 
01583       const typename __ios_base::fmtflags __flags = __os.flags();
01584       const _CharT __fill = __os.fill();
01585       const std::streamsize __precision = __os.precision();
01586       const _CharT __space = __os.widen(' ');
01587       __os.flags(__ios_base::scientific | __ios_base::left);
01588       __os.fill(__space);
01589       __os.precision(std::numeric_limits<double>::max_digits10);
01590 
01591       __os << __x.t() << __space << __x.p()
01592        << __space << __x._M_nd;
01593 
01594       __os.flags(__flags);
01595       __os.fill(__fill);
01596       __os.precision(__precision);
01597       return __os;
01598     }
01599 
01600   template<typename _IntType,
01601        typename _CharT, typename _Traits>
01602     std::basic_istream<_CharT, _Traits>&
01603     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01604            binomial_distribution<_IntType>& __x)
01605     {
01606       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01607       typedef typename __istream_type::ios_base    __ios_base;
01608 
01609       const typename __ios_base::fmtflags __flags = __is.flags();
01610       __is.flags(__ios_base::dec | __ios_base::skipws);
01611 
01612       _IntType __t;
01613       double __p;
01614       __is >> __t >> __p >> __x._M_nd;
01615       __x.param(typename binomial_distribution<_IntType>::
01616         param_type(__t, __p));
01617 
01618       __is.flags(__flags);
01619       return __is;
01620     }
01621 
01622 
01623   template<typename _RealType, typename _CharT, typename _Traits>
01624     std::basic_ostream<_CharT, _Traits>&
01625     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01626            const exponential_distribution<_RealType>& __x)
01627     {
01628       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01629       typedef typename __ostream_type::ios_base    __ios_base;
01630 
01631       const typename __ios_base::fmtflags __flags = __os.flags();
01632       const _CharT __fill = __os.fill();
01633       const std::streamsize __precision = __os.precision();
01634       __os.flags(__ios_base::scientific | __ios_base::left);
01635       __os.fill(__os.widen(' '));
01636       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01637 
01638       __os << __x.lambda();
01639 
01640       __os.flags(__flags);
01641       __os.fill(__fill);
01642       __os.precision(__precision);
01643       return __os;
01644     }
01645 
01646   template<typename _RealType, typename _CharT, typename _Traits>
01647     std::basic_istream<_CharT, _Traits>&
01648     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01649            exponential_distribution<_RealType>& __x)
01650     {
01651       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01652       typedef typename __istream_type::ios_base    __ios_base;
01653 
01654       const typename __ios_base::fmtflags __flags = __is.flags();
01655       __is.flags(__ios_base::dec | __ios_base::skipws);
01656 
01657       _RealType __lambda;
01658       __is >> __lambda;
01659       __x.param(typename exponential_distribution<_RealType>::
01660         param_type(__lambda));
01661 
01662       __is.flags(__flags);
01663       return __is;
01664     }
01665 
01666 
01667   /**
01668    * Polar method due to Marsaglia.
01669    *
01670    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01671    * New York, 1986, Ch. V, Sect. 4.4.
01672    */
01673   template<typename _RealType>
01674     template<typename _UniformRandomNumberGenerator>
01675       typename normal_distribution<_RealType>::result_type
01676       normal_distribution<_RealType>::
01677       operator()(_UniformRandomNumberGenerator& __urng,
01678          const param_type& __param)
01679       {
01680     result_type __ret;
01681     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01682       __aurng(__urng);
01683 
01684     if (_M_saved_available)
01685       {
01686         _M_saved_available = false;
01687         __ret = _M_saved;
01688       }
01689     else
01690       {
01691         result_type __x, __y, __r2;
01692         do
01693           {
01694         __x = result_type(2.0) * __aurng() - 1.0;
01695         __y = result_type(2.0) * __aurng() - 1.0;
01696         __r2 = __x * __x + __y * __y;
01697           }
01698         while (__r2 > 1.0 || __r2 == 0.0);
01699 
01700         const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
01701         _M_saved = __x * __mult;
01702         _M_saved_available = true;
01703         __ret = __y * __mult;
01704       }
01705 
01706     __ret = __ret * __param.stddev() + __param.mean();
01707     return __ret;
01708       }
01709 
01710   template<typename _RealType>
01711     bool
01712     operator==(const std::normal_distribution<_RealType>& __d1,
01713            const std::normal_distribution<_RealType>& __d2)
01714     {
01715       if (__d1._M_param == __d2._M_param
01716       && __d1._M_saved_available == __d2._M_saved_available)
01717     {
01718       if (__d1._M_saved_available
01719           && __d1._M_saved == __d2._M_saved)
01720         return true;
01721       else if(!__d1._M_saved_available)
01722         return true;
01723       else
01724         return false;
01725     }
01726       else
01727     return false;
01728     }
01729 
01730   template<typename _RealType, typename _CharT, typename _Traits>
01731     std::basic_ostream<_CharT, _Traits>&
01732     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01733            const normal_distribution<_RealType>& __x)
01734     {
01735       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01736       typedef typename __ostream_type::ios_base    __ios_base;
01737 
01738       const typename __ios_base::fmtflags __flags = __os.flags();
01739       const _CharT __fill = __os.fill();
01740       const std::streamsize __precision = __os.precision();
01741       const _CharT __space = __os.widen(' ');
01742       __os.flags(__ios_base::scientific | __ios_base::left);
01743       __os.fill(__space);
01744       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01745 
01746       __os << __x.mean() << __space << __x.stddev()
01747        << __space << __x._M_saved_available;
01748       if (__x._M_saved_available)
01749     __os << __space << __x._M_saved;
01750 
01751       __os.flags(__flags);
01752       __os.fill(__fill);
01753       __os.precision(__precision);
01754       return __os;
01755     }
01756 
01757   template<typename _RealType, typename _CharT, typename _Traits>
01758     std::basic_istream<_CharT, _Traits>&
01759     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01760            normal_distribution<_RealType>& __x)
01761     {
01762       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01763       typedef typename __istream_type::ios_base    __ios_base;
01764 
01765       const typename __ios_base::fmtflags __flags = __is.flags();
01766       __is.flags(__ios_base::dec | __ios_base::skipws);
01767 
01768       double __mean, __stddev;
01769       __is >> __mean >> __stddev
01770        >> __x._M_saved_available;
01771       if (__x._M_saved_available)
01772     __is >> __x._M_saved;
01773       __x.param(typename normal_distribution<_RealType>::
01774         param_type(__mean, __stddev));
01775 
01776       __is.flags(__flags);
01777       return __is;
01778     }
01779 
01780 
01781   template<typename _RealType, typename _CharT, typename _Traits>
01782     std::basic_ostream<_CharT, _Traits>&
01783     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01784            const lognormal_distribution<_RealType>& __x)
01785     {
01786       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01787       typedef typename __ostream_type::ios_base    __ios_base;
01788 
01789       const typename __ios_base::fmtflags __flags = __os.flags();
01790       const _CharT __fill = __os.fill();
01791       const std::streamsize __precision = __os.precision();
01792       const _CharT __space = __os.widen(' ');
01793       __os.flags(__ios_base::scientific | __ios_base::left);
01794       __os.fill(__space);
01795       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01796 
01797       __os << __x.m() << __space << __x.s()
01798        << __space << __x._M_nd;
01799 
01800       __os.flags(__flags);
01801       __os.fill(__fill);
01802       __os.precision(__precision);
01803       return __os;
01804     }
01805 
01806   template<typename _RealType, typename _CharT, typename _Traits>
01807     std::basic_istream<_CharT, _Traits>&
01808     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01809            lognormal_distribution<_RealType>& __x)
01810     {
01811       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01812       typedef typename __istream_type::ios_base    __ios_base;
01813 
01814       const typename __ios_base::fmtflags __flags = __is.flags();
01815       __is.flags(__ios_base::dec | __ios_base::skipws);
01816 
01817       _RealType __m, __s;
01818       __is >> __m >> __s >> __x._M_nd;
01819       __x.param(typename lognormal_distribution<_RealType>::
01820         param_type(__m, __s));
01821 
01822       __is.flags(__flags);
01823       return __is;
01824     }
01825 
01826 
01827   template<typename _RealType, typename _CharT, typename _Traits>
01828     std::basic_ostream<_CharT, _Traits>&
01829     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01830            const chi_squared_distribution<_RealType>& __x)
01831     {
01832       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01833       typedef typename __ostream_type::ios_base    __ios_base;
01834 
01835       const typename __ios_base::fmtflags __flags = __os.flags();
01836       const _CharT __fill = __os.fill();
01837       const std::streamsize __precision = __os.precision();
01838       const _CharT __space = __os.widen(' ');
01839       __os.flags(__ios_base::scientific | __ios_base::left);
01840       __os.fill(__space);
01841       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01842 
01843       __os << __x.n() << __space << __x._M_gd;
01844 
01845       __os.flags(__flags);
01846       __os.fill(__fill);
01847       __os.precision(__precision);
01848       return __os;
01849     }
01850 
01851   template<typename _RealType, typename _CharT, typename _Traits>
01852     std::basic_istream<_CharT, _Traits>&
01853     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01854            chi_squared_distribution<_RealType>& __x)
01855     {
01856       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01857       typedef typename __istream_type::ios_base    __ios_base;
01858 
01859       const typename __ios_base::fmtflags __flags = __is.flags();
01860       __is.flags(__ios_base::dec | __ios_base::skipws);
01861 
01862       _RealType __n;
01863       __is >> __n >> __x._M_gd;
01864       __x.param(typename chi_squared_distribution<_RealType>::
01865         param_type(__n));
01866 
01867       __is.flags(__flags);
01868       return __is;
01869     }
01870 
01871 
01872   template<typename _RealType>
01873     template<typename _UniformRandomNumberGenerator>
01874       typename cauchy_distribution<_RealType>::result_type
01875       cauchy_distribution<_RealType>::
01876       operator()(_UniformRandomNumberGenerator& __urng,
01877          const param_type& __p)
01878       {
01879     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01880       __aurng(__urng);
01881     _RealType __u;
01882     do
01883       __u = __aurng();
01884     while (__u == 0.5);
01885 
01886     const _RealType __pi = 3.1415926535897932384626433832795029L;
01887     return __p.a() + __p.b() * std::tan(__pi * __u);
01888       }
01889 
01890   template<typename _RealType, typename _CharT, typename _Traits>
01891     std::basic_ostream<_CharT, _Traits>&
01892     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01893            const cauchy_distribution<_RealType>& __x)
01894     {
01895       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01896       typedef typename __ostream_type::ios_base    __ios_base;
01897 
01898       const typename __ios_base::fmtflags __flags = __os.flags();
01899       const _CharT __fill = __os.fill();
01900       const std::streamsize __precision = __os.precision();
01901       const _CharT __space = __os.widen(' ');
01902       __os.flags(__ios_base::scientific | __ios_base::left);
01903       __os.fill(__space);
01904       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01905 
01906       __os << __x.a() << __space << __x.b();
01907 
01908       __os.flags(__flags);
01909       __os.fill(__fill);
01910       __os.precision(__precision);
01911       return __os;
01912     }
01913 
01914   template<typename _RealType, typename _CharT, typename _Traits>
01915     std::basic_istream<_CharT, _Traits>&
01916     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01917            cauchy_distribution<_RealType>& __x)
01918     {
01919       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01920       typedef typename __istream_type::ios_base    __ios_base;
01921 
01922       const typename __ios_base::fmtflags __flags = __is.flags();
01923       __is.flags(__ios_base::dec | __ios_base::skipws);
01924 
01925       _RealType __a, __b;
01926       __is >> __a >> __b;
01927       __x.param(typename cauchy_distribution<_RealType>::
01928         param_type(__a, __b));
01929 
01930       __is.flags(__flags);
01931       return __is;
01932     }
01933 
01934 
01935   template<typename _RealType, typename _CharT, typename _Traits>
01936     std::basic_ostream<_CharT, _Traits>&
01937     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01938            const fisher_f_distribution<_RealType>& __x)
01939     {
01940       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01941       typedef typename __ostream_type::ios_base    __ios_base;
01942 
01943       const typename __ios_base::fmtflags __flags = __os.flags();
01944       const _CharT __fill = __os.fill();
01945       const std::streamsize __precision = __os.precision();
01946       const _CharT __space = __os.widen(' ');
01947       __os.flags(__ios_base::scientific | __ios_base::left);
01948       __os.fill(__space);
01949       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01950 
01951       __os << __x.m() << __space << __x.n()
01952        << __space << __x._M_gd_x << __space << __x._M_gd_y;
01953 
01954       __os.flags(__flags);
01955       __os.fill(__fill);
01956       __os.precision(__precision);
01957       return __os;
01958     }
01959 
01960   template<typename _RealType, typename _CharT, typename _Traits>
01961     std::basic_istream<_CharT, _Traits>&
01962     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01963            fisher_f_distribution<_RealType>& __x)
01964     {
01965       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01966       typedef typename __istream_type::ios_base    __ios_base;
01967 
01968       const typename __ios_base::fmtflags __flags = __is.flags();
01969       __is.flags(__ios_base::dec | __ios_base::skipws);
01970 
01971       _RealType __m, __n;
01972       __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
01973       __x.param(typename fisher_f_distribution<_RealType>::
01974         param_type(__m, __n));
01975 
01976       __is.flags(__flags);
01977       return __is;
01978     }
01979 
01980 
01981   template<typename _RealType, typename _CharT, typename _Traits>
01982     std::basic_ostream<_CharT, _Traits>&
01983     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01984            const student_t_distribution<_RealType>& __x)
01985     {
01986       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01987       typedef typename __ostream_type::ios_base    __ios_base;
01988 
01989       const typename __ios_base::fmtflags __flags = __os.flags();
01990       const _CharT __fill = __os.fill();
01991       const std::streamsize __precision = __os.precision();
01992       const _CharT __space = __os.widen(' ');
01993       __os.flags(__ios_base::scientific | __ios_base::left);
01994       __os.fill(__space);
01995       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01996 
01997       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
01998 
01999       __os.flags(__flags);
02000       __os.fill(__fill);
02001       __os.precision(__precision);
02002       return __os;
02003     }
02004 
02005   template<typename _RealType, typename _CharT, typename _Traits>
02006     std::basic_istream<_CharT, _Traits>&
02007     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02008            student_t_distribution<_RealType>& __x)
02009     {
02010       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02011       typedef typename __istream_type::ios_base    __ios_base;
02012 
02013       const typename __ios_base::fmtflags __flags = __is.flags();
02014       __is.flags(__ios_base::dec | __ios_base::skipws);
02015 
02016       _RealType __n;
02017       __is >> __n >> __x._M_nd >> __x._M_gd;
02018       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
02019 
02020       __is.flags(__flags);
02021       return __is;
02022     }
02023 
02024 
02025   template<typename _RealType>
02026     void
02027     gamma_distribution<_RealType>::param_type::
02028     _M_initialize()
02029     {
02030       _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
02031 
02032       const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
02033       _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
02034     }
02035 
02036   /**
02037    * Marsaglia, G. and Tsang, W. W.
02038    * "A Simple Method for Generating Gamma Variables"
02039    * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
02040    */
02041   template<typename _RealType>
02042     template<typename _UniformRandomNumberGenerator>
02043       typename gamma_distribution<_RealType>::result_type
02044       gamma_distribution<_RealType>::
02045       operator()(_UniformRandomNumberGenerator& __urng,
02046          const param_type& __param)
02047       {
02048     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02049       __aurng(__urng);
02050 
02051     result_type __u, __v, __n;
02052     const result_type __a1 = (__param._M_malpha
02053                   - _RealType(1.0) / _RealType(3.0));
02054 
02055     do
02056       {
02057         do
02058           {
02059         __n = _M_nd(__urng);
02060         __v = result_type(1.0) + __param._M_a2 * __n; 
02061           }
02062         while (__v <= 0.0);
02063 
02064         __v = __v * __v * __v;
02065         __u = __aurng();
02066       }
02067     while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
02068            && (std::log(__u) > (0.5 * __n * __n + __a1
02069                     * (1.0 - __v + std::log(__v)))));
02070 
02071     if (__param.alpha() == __param._M_malpha)
02072       return __a1 * __v * __param.beta();
02073     else
02074       {
02075         do
02076           __u = __aurng();
02077         while (__u == 0.0);
02078         
02079         return (std::pow(__u, result_type(1.0) / __param.alpha())
02080             * __a1 * __v * __param.beta());
02081       }
02082       }
02083 
02084   template<typename _RealType, typename _CharT, typename _Traits>
02085     std::basic_ostream<_CharT, _Traits>&
02086     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02087            const gamma_distribution<_RealType>& __x)
02088     {
02089       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02090       typedef typename __ostream_type::ios_base    __ios_base;
02091 
02092       const typename __ios_base::fmtflags __flags = __os.flags();
02093       const _CharT __fill = __os.fill();
02094       const std::streamsize __precision = __os.precision();
02095       const _CharT __space = __os.widen(' ');
02096       __os.flags(__ios_base::scientific | __ios_base::left);
02097       __os.fill(__space);
02098       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02099 
02100       __os << __x.alpha() << __space << __x.beta()
02101        << __space << __x._M_nd;
02102 
02103       __os.flags(__flags);
02104       __os.fill(__fill);
02105       __os.precision(__precision);
02106       return __os;
02107     }
02108 
02109   template<typename _RealType, typename _CharT, typename _Traits>
02110     std::basic_istream<_CharT, _Traits>&
02111     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02112            gamma_distribution<_RealType>& __x)
02113     {
02114       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02115       typedef typename __istream_type::ios_base    __ios_base;
02116 
02117       const typename __ios_base::fmtflags __flags = __is.flags();
02118       __is.flags(__ios_base::dec | __ios_base::skipws);
02119 
02120       _RealType __alpha_val, __beta_val;
02121       __is >> __alpha_val >> __beta_val >> __x._M_nd;
02122       __x.param(typename gamma_distribution<_RealType>::
02123         param_type(__alpha_val, __beta_val));
02124 
02125       __is.flags(__flags);
02126       return __is;
02127     }
02128 
02129 
02130   template<typename _RealType>
02131     template<typename _UniformRandomNumberGenerator>
02132       typename weibull_distribution<_RealType>::result_type
02133       weibull_distribution<_RealType>::
02134       operator()(_UniformRandomNumberGenerator& __urng,
02135          const param_type& __p)
02136       {
02137     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02138       __aurng(__urng);
02139     return __p.b() * std::pow(-std::log(__aurng()),
02140                   result_type(1) / __p.a());
02141       }
02142 
02143   template<typename _RealType, typename _CharT, typename _Traits>
02144     std::basic_ostream<_CharT, _Traits>&
02145     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02146            const weibull_distribution<_RealType>& __x)
02147     {
02148       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02149       typedef typename __ostream_type::ios_base    __ios_base;
02150 
02151       const typename __ios_base::fmtflags __flags = __os.flags();
02152       const _CharT __fill = __os.fill();
02153       const std::streamsize __precision = __os.precision();
02154       const _CharT __space = __os.widen(' ');
02155       __os.flags(__ios_base::scientific | __ios_base::left);
02156       __os.fill(__space);
02157       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02158 
02159       __os << __x.a() << __space << __x.b();
02160 
02161       __os.flags(__flags);
02162       __os.fill(__fill);
02163       __os.precision(__precision);
02164       return __os;
02165     }
02166 
02167   template<typename _RealType, typename _CharT, typename _Traits>
02168     std::basic_istream<_CharT, _Traits>&
02169     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02170            weibull_distribution<_RealType>& __x)
02171     {
02172       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02173       typedef typename __istream_type::ios_base    __ios_base;
02174 
02175       const typename __ios_base::fmtflags __flags = __is.flags();
02176       __is.flags(__ios_base::dec | __ios_base::skipws);
02177 
02178       _RealType __a, __b;
02179       __is >> __a >> __b;
02180       __x.param(typename weibull_distribution<_RealType>::
02181         param_type(__a, __b));
02182 
02183       __is.flags(__flags);
02184       return __is;
02185     }
02186 
02187 
02188   template<typename _RealType>
02189     template<typename _UniformRandomNumberGenerator>
02190       typename extreme_value_distribution<_RealType>::result_type
02191       extreme_value_distribution<_RealType>::
02192       operator()(_UniformRandomNumberGenerator& __urng,
02193          const param_type& __p)
02194       {
02195     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02196       __aurng(__urng);
02197     return __p.a() - __p.b() * std::log(-std::log(__aurng()));
02198       }
02199 
02200   template<typename _RealType, typename _CharT, typename _Traits>
02201     std::basic_ostream<_CharT, _Traits>&
02202     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02203            const extreme_value_distribution<_RealType>& __x)
02204     {
02205       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02206       typedef typename __ostream_type::ios_base    __ios_base;
02207 
02208       const typename __ios_base::fmtflags __flags = __os.flags();
02209       const _CharT __fill = __os.fill();
02210       const std::streamsize __precision = __os.precision();
02211       const _CharT __space = __os.widen(' ');
02212       __os.flags(__ios_base::scientific | __ios_base::left);
02213       __os.fill(__space);
02214       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02215 
02216       __os << __x.a() << __space << __x.b();
02217 
02218       __os.flags(__flags);
02219       __os.fill(__fill);
02220       __os.precision(__precision);
02221       return __os;
02222     }
02223 
02224   template<typename _RealType, typename _CharT, typename _Traits>
02225     std::basic_istream<_CharT, _Traits>&
02226     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02227            extreme_value_distribution<_RealType>& __x)
02228     {
02229       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02230       typedef typename __istream_type::ios_base    __ios_base;
02231 
02232       const typename __ios_base::fmtflags __flags = __is.flags();
02233       __is.flags(__ios_base::dec | __ios_base::skipws);
02234 
02235       _RealType __a, __b;
02236       __is >> __a >> __b;
02237       __x.param(typename extreme_value_distribution<_RealType>::
02238         param_type(__a, __b));
02239 
02240       __is.flags(__flags);
02241       return __is;
02242     }
02243 
02244 
02245   template<typename _IntType>
02246     void
02247     discrete_distribution<_IntType>::param_type::
02248     _M_initialize()
02249     {
02250       if (_M_prob.size() < 2)
02251     {
02252       _M_prob.clear();
02253       return;
02254     }
02255 
02256       const double __sum = std::accumulate(_M_prob.begin(),
02257                        _M_prob.end(), 0.0);
02258       // Now normalize the probabilites.
02259       __detail::__transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
02260               std::bind2nd(std::divides<double>(), __sum));
02261       // Accumulate partial sums.
02262       _M_cp.reserve(_M_prob.size());
02263       std::partial_sum(_M_prob.begin(), _M_prob.end(),
02264                std::back_inserter(_M_cp));
02265       // Make sure the last cumulative probability is one.
02266       _M_cp[_M_cp.size() - 1] = 1.0;
02267     }
02268 
02269   template<typename _IntType>
02270     template<typename _Func>
02271       discrete_distribution<_IntType>::param_type::
02272       param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
02273       : _M_prob(), _M_cp()
02274       {
02275     const size_t __n = __nw == 0 ? 1 : __nw;
02276     const double __delta = (__xmax - __xmin) / __n;
02277 
02278     _M_prob.reserve(__n);
02279     for (size_t __k = 0; __k < __nw; ++__k)
02280       _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
02281 
02282     _M_initialize();
02283       }
02284 
02285   template<typename _IntType>
02286     template<typename _UniformRandomNumberGenerator>
02287       typename discrete_distribution<_IntType>::result_type
02288       discrete_distribution<_IntType>::
02289       operator()(_UniformRandomNumberGenerator& __urng,
02290          const param_type& __param)
02291       {
02292     if (__param._M_cp.empty())
02293       return result_type(0);
02294 
02295     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02296       __aurng(__urng);
02297 
02298     const double __p = __aurng();
02299     auto __pos = std::lower_bound(__param._M_cp.begin(),
02300                       __param._M_cp.end(), __p);
02301 
02302     return __pos - __param._M_cp.begin();
02303       }
02304 
02305   template<typename _IntType, typename _CharT, typename _Traits>
02306     std::basic_ostream<_CharT, _Traits>&
02307     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02308            const discrete_distribution<_IntType>& __x)
02309     {
02310       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02311       typedef typename __ostream_type::ios_base    __ios_base;
02312 
02313       const typename __ios_base::fmtflags __flags = __os.flags();
02314       const _CharT __fill = __os.fill();
02315       const std::streamsize __precision = __os.precision();
02316       const _CharT __space = __os.widen(' ');
02317       __os.flags(__ios_base::scientific | __ios_base::left);
02318       __os.fill(__space);
02319       __os.precision(std::numeric_limits<double>::max_digits10);
02320 
02321       std::vector<double> __prob = __x.probabilities();
02322       __os << __prob.size();
02323       for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
02324     __os << __space << *__dit;
02325 
02326       __os.flags(__flags);
02327       __os.fill(__fill);
02328       __os.precision(__precision);
02329       return __os;
02330     }
02331 
02332   template<typename _IntType, typename _CharT, typename _Traits>
02333     std::basic_istream<_CharT, _Traits>&
02334     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02335            discrete_distribution<_IntType>& __x)
02336     {
02337       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02338       typedef typename __istream_type::ios_base    __ios_base;
02339 
02340       const typename __ios_base::fmtflags __flags = __is.flags();
02341       __is.flags(__ios_base::dec | __ios_base::skipws);
02342 
02343       size_t __n;
02344       __is >> __n;
02345 
02346       std::vector<double> __prob_vec;
02347       __prob_vec.reserve(__n);
02348       for (; __n != 0; --__n)
02349     {
02350       double __prob;
02351       __is >> __prob;
02352       __prob_vec.push_back(__prob);
02353     }
02354 
02355       __x.param(typename discrete_distribution<_IntType>::
02356         param_type(__prob_vec.begin(), __prob_vec.end()));
02357 
02358       __is.flags(__flags);
02359       return __is;
02360     }
02361 
02362 
02363   template<typename _RealType>
02364     void
02365     piecewise_constant_distribution<_RealType>::param_type::
02366     _M_initialize()
02367     {
02368       if (_M_int.size() < 2
02369       || (_M_int.size() == 2
02370           && _M_int[0] == _RealType(0)
02371           && _M_int[1] == _RealType(1)))
02372     {
02373       _M_int.clear();
02374       _M_den.clear();
02375       return;
02376     }
02377 
02378       const double __sum = std::accumulate(_M_den.begin(),
02379                        _M_den.end(), 0.0);
02380 
02381       __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
02382                 std::bind2nd(std::divides<double>(), __sum));
02383 
02384       _M_cp.reserve(_M_den.size());
02385       std::partial_sum(_M_den.begin(), _M_den.end(),
02386                std::back_inserter(_M_cp));
02387 
02388       // Make sure the last cumulative probability is one.
02389       _M_cp[_M_cp.size() - 1] = 1.0;
02390 
02391       for (size_t __k = 0; __k < _M_den.size(); ++__k)
02392     _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
02393     }
02394 
02395   template<typename _RealType>
02396     template<typename _InputIteratorB, typename _InputIteratorW>
02397       piecewise_constant_distribution<_RealType>::param_type::
02398       param_type(_InputIteratorB __bbegin,
02399          _InputIteratorB __bend,
02400          _InputIteratorW __wbegin)
02401       : _M_int(), _M_den(), _M_cp()
02402       {
02403     if (__bbegin != __bend)
02404       {
02405         for (;;)
02406           {
02407         _M_int.push_back(*__bbegin);
02408         ++__bbegin;
02409         if (__bbegin == __bend)
02410           break;
02411 
02412         _M_den.push_back(*__wbegin);
02413         ++__wbegin;
02414           }
02415       }
02416 
02417     _M_initialize();
02418       }
02419 
02420   template<typename _RealType>
02421     template<typename _Func>
02422       piecewise_constant_distribution<_RealType>::param_type::
02423       param_type(initializer_list<_RealType> __bl, _Func __fw)
02424       : _M_int(), _M_den(), _M_cp()
02425       {
02426     _M_int.reserve(__bl.size());
02427     for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
02428       _M_int.push_back(*__biter);
02429 
02430     _M_den.reserve(_M_int.size() - 1);
02431     for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
02432       _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
02433 
02434     _M_initialize();
02435       }
02436 
02437   template<typename _RealType>
02438     template<typename _Func>
02439       piecewise_constant_distribution<_RealType>::param_type::
02440       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
02441       : _M_int(), _M_den(), _M_cp()
02442       {
02443     const size_t __n = __nw == 0 ? 1 : __nw;
02444     const _RealType __delta = (__xmax - __xmin) / __n;
02445 
02446     _M_int.reserve(__n + 1);
02447     for (size_t __k = 0; __k <= __nw; ++__k)
02448       _M_int.push_back(__xmin + __k * __delta);
02449 
02450     _M_den.reserve(__n);
02451     for (size_t __k = 0; __k < __nw; ++__k)
02452       _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
02453 
02454     _M_initialize();
02455       }
02456 
02457   template<typename _RealType>
02458     template<typename _UniformRandomNumberGenerator>
02459       typename piecewise_constant_distribution<_RealType>::result_type
02460       piecewise_constant_distribution<_RealType>::
02461       operator()(_UniformRandomNumberGenerator& __urng,
02462          const param_type& __param)
02463       {
02464     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02465       __aurng(__urng);
02466 
02467     const double __p = __aurng();
02468     if (__param._M_cp.empty())
02469       return __p;
02470 
02471     auto __pos = std::lower_bound(__param._M_cp.begin(),
02472                       __param._M_cp.end(), __p);
02473     const size_t __i = __pos - __param._M_cp.begin();
02474 
02475     const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
02476 
02477     return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
02478       }
02479 
02480   template<typename _RealType, typename _CharT, typename _Traits>
02481     std::basic_ostream<_CharT, _Traits>&
02482     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02483            const piecewise_constant_distribution<_RealType>& __x)
02484     {
02485       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02486       typedef typename __ostream_type::ios_base    __ios_base;
02487 
02488       const typename __ios_base::fmtflags __flags = __os.flags();
02489       const _CharT __fill = __os.fill();
02490       const std::streamsize __precision = __os.precision();
02491       const _CharT __space = __os.widen(' ');
02492       __os.flags(__ios_base::scientific | __ios_base::left);
02493       __os.fill(__space);
02494       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02495 
02496       std::vector<_RealType> __int = __x.intervals();
02497       __os << __int.size() - 1;
02498 
02499       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
02500     __os << __space << *__xit;
02501 
02502       std::vector<double> __den = __x.densities();
02503       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
02504     __os << __space << *__dit;
02505 
02506       __os.flags(__flags);
02507       __os.fill(__fill);
02508       __os.precision(__precision);
02509       return __os;
02510     }
02511 
02512   template<typename _RealType, typename _CharT, typename _Traits>
02513     std::basic_istream<_CharT, _Traits>&
02514     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02515            piecewise_constant_distribution<_RealType>& __x)
02516     {
02517       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02518       typedef typename __istream_type::ios_base    __ios_base;
02519 
02520       const typename __ios_base::fmtflags __flags = __is.flags();
02521       __is.flags(__ios_base::dec | __ios_base::skipws);
02522 
02523       size_t __n;
02524       __is >> __n;
02525 
02526       std::vector<_RealType> __int_vec;
02527       __int_vec.reserve(__n + 1);
02528       for (size_t __i = 0; __i <= __n; ++__i)
02529     {
02530       _RealType __int;
02531       __is >> __int;
02532       __int_vec.push_back(__int);
02533     }
02534 
02535       std::vector<double> __den_vec;
02536       __den_vec.reserve(__n);
02537       for (size_t __i = 0; __i < __n; ++__i)
02538     {
02539       double __den;
02540       __is >> __den;
02541       __den_vec.push_back(__den);
02542     }
02543 
02544       __x.param(typename piecewise_constant_distribution<_RealType>::
02545       param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
02546 
02547       __is.flags(__flags);
02548       return __is;
02549     }
02550 
02551 
02552   template<typename _RealType>
02553     void
02554     piecewise_linear_distribution<_RealType>::param_type::
02555     _M_initialize()
02556     {
02557       if (_M_int.size() < 2
02558       || (_M_int.size() == 2
02559           && _M_int[0] == _RealType(0)
02560           && _M_int[1] == _RealType(1)
02561           && _M_den[0] == _M_den[1]))
02562     {
02563       _M_int.clear();
02564       _M_den.clear();
02565       return;
02566     }
02567 
02568       double __sum = 0.0;
02569       _M_cp.reserve(_M_int.size() - 1);
02570       _M_m.reserve(_M_int.size() - 1);
02571       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
02572     {
02573       const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
02574       __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
02575       _M_cp.push_back(__sum);
02576       _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
02577     }
02578 
02579       //  Now normalize the densities...
02580       __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
02581               std::bind2nd(std::divides<double>(), __sum));
02582       //  ... and partial sums... 
02583       __detail::__transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
02584                 std::bind2nd(std::divides<double>(), __sum));
02585       //  ... and slopes.
02586       __detail::__transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
02587                 std::bind2nd(std::divides<double>(), __sum));
02588       //  Make sure the last cumulative probablility is one.
02589       _M_cp[_M_cp.size() - 1] = 1.0;
02590      }
02591 
02592   template<typename _RealType>
02593     template<typename _InputIteratorB, typename _InputIteratorW>
02594       piecewise_linear_distribution<_RealType>::param_type::
02595       param_type(_InputIteratorB __bbegin,
02596          _InputIteratorB __bend,
02597          _InputIteratorW __wbegin)
02598       : _M_int(), _M_den(), _M_cp(), _M_m()
02599       {
02600     for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
02601       {
02602         _M_int.push_back(*__bbegin);
02603         _M_den.push_back(*__wbegin);
02604       }
02605 
02606     _M_initialize();
02607       }
02608 
02609   template<typename _RealType>
02610     template<typename _Func>
02611       piecewise_linear_distribution<_RealType>::param_type::
02612       param_type(initializer_list<_RealType> __bl, _Func __fw)
02613       : _M_int(), _M_den(), _M_cp(), _M_m()
02614       {
02615     _M_int.reserve(__bl.size());
02616     _M_den.reserve(__bl.size());
02617     for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
02618       {
02619         _M_int.push_back(*__biter);
02620         _M_den.push_back(__fw(*__biter));
02621       }
02622 
02623     _M_initialize();
02624       }
02625 
02626   template<typename _RealType>
02627     template<typename _Func>
02628       piecewise_linear_distribution<_RealType>::param_type::
02629       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
02630       : _M_int(), _M_den(), _M_cp(), _M_m()
02631       {
02632     const size_t __n = __nw == 0 ? 1 : __nw;
02633     const _RealType __delta = (__xmax - __xmin) / __n;
02634 
02635     _M_int.reserve(__n + 1);
02636     _M_den.reserve(__n + 1);
02637     for (size_t __k = 0; __k <= __nw; ++__k)
02638       {
02639         _M_int.push_back(__xmin + __k * __delta);
02640         _M_den.push_back(__fw(_M_int[__k] + __delta));
02641       }
02642 
02643     _M_initialize();
02644       }
02645 
02646   template<typename _RealType>
02647     template<typename _UniformRandomNumberGenerator>
02648       typename piecewise_linear_distribution<_RealType>::result_type
02649       piecewise_linear_distribution<_RealType>::
02650       operator()(_UniformRandomNumberGenerator& __urng,
02651          const param_type& __param)
02652       {
02653     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02654       __aurng(__urng);
02655 
02656     const double __p = __aurng();
02657     if (__param._M_cp.empty())
02658       return __p;
02659 
02660     auto __pos = std::lower_bound(__param._M_cp.begin(),
02661                       __param._M_cp.end(), __p);
02662     const size_t __i = __pos - __param._M_cp.begin();
02663 
02664     const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
02665 
02666     const double __a = 0.5 * __param._M_m[__i];
02667     const double __b = __param._M_den[__i];
02668     const double __cm = __p - __pref;
02669 
02670     _RealType __x = __param._M_int[__i];
02671     if (__a == 0)
02672       __x += __cm / __b;
02673     else
02674       {
02675         const double __d = __b * __b + 4.0 * __a * __cm;
02676         __x += 0.5 * (std::sqrt(__d) - __b) / __a;
02677           }
02678 
02679         return __x;
02680       }
02681 
02682   template<typename _RealType, typename _CharT, typename _Traits>
02683     std::basic_ostream<_CharT, _Traits>&
02684     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02685            const piecewise_linear_distribution<_RealType>& __x)
02686     {
02687       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02688       typedef typename __ostream_type::ios_base    __ios_base;
02689 
02690       const typename __ios_base::fmtflags __flags = __os.flags();
02691       const _CharT __fill = __os.fill();
02692       const std::streamsize __precision = __os.precision();
02693       const _CharT __space = __os.widen(' ');
02694       __os.flags(__ios_base::scientific | __ios_base::left);
02695       __os.fill(__space);
02696       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02697 
02698       std::vector<_RealType> __int = __x.intervals();
02699       __os << __int.size() - 1;
02700 
02701       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
02702     __os << __space << *__xit;
02703 
02704       std::vector<double> __den = __x.densities();
02705       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
02706     __os << __space << *__dit;
02707 
02708       __os.flags(__flags);
02709       __os.fill(__fill);
02710       __os.precision(__precision);
02711       return __os;
02712     }
02713 
02714   template<typename _RealType, typename _CharT, typename _Traits>
02715     std::basic_istream<_CharT, _Traits>&
02716     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02717            piecewise_linear_distribution<_RealType>& __x)
02718     {
02719       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02720       typedef typename __istream_type::ios_base    __ios_base;
02721 
02722       const typename __ios_base::fmtflags __flags = __is.flags();
02723       __is.flags(__ios_base::dec | __ios_base::skipws);
02724 
02725       size_t __n;
02726       __is >> __n;
02727 
02728       std::vector<_RealType> __int_vec;
02729       __int_vec.reserve(__n + 1);
02730       for (size_t __i = 0; __i <= __n; ++__i)
02731     {
02732       _RealType __int;
02733       __is >> __int;
02734       __int_vec.push_back(__int);
02735     }
02736 
02737       std::vector<double> __den_vec;
02738       __den_vec.reserve(__n + 1);
02739       for (size_t __i = 0; __i <= __n; ++__i)
02740     {
02741       double __den;
02742       __is >> __den;
02743       __den_vec.push_back(__den);
02744     }
02745 
02746       __x.param(typename piecewise_linear_distribution<_RealType>::
02747       param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
02748 
02749       __is.flags(__flags);
02750       return __is;
02751     }
02752 
02753 
02754   template<typename _IntType>
02755     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
02756     {
02757       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
02758     _M_v.push_back(__detail::__mod<result_type,
02759                __detail::_Shift<result_type, 32>::__value>(*__iter));
02760     }
02761 
02762   template<typename _InputIterator>
02763     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
02764     {
02765       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
02766     _M_v.push_back(__detail::__mod<result_type,
02767                __detail::_Shift<result_type, 32>::__value>(*__iter));
02768     }
02769 
02770   template<typename _RandomAccessIterator>
02771     void
02772     seed_seq::generate(_RandomAccessIterator __begin,
02773                _RandomAccessIterator __end)
02774     {
02775       typedef typename iterator_traits<_RandomAccessIterator>::value_type
02776         _Type;
02777 
02778       if (__begin == __end)
02779     return;
02780 
02781       std::fill(__begin, __end, _Type(0x8b8b8b8bu));
02782 
02783       const size_t __n = __end - __begin;
02784       const size_t __s = _M_v.size();
02785       const size_t __t = (__n >= 623) ? 11
02786                : (__n >=  68) ? 7
02787                : (__n >=  39) ? 5
02788                : (__n >=   7) ? 3
02789                : (__n - 1) / 2;
02790       const size_t __p = (__n - __t) / 2;
02791       const size_t __q = __p + __t;
02792       const size_t __m = std::max(__s + 1, __n);
02793 
02794       for (size_t __k = 0; __k < __m; ++__k)
02795     {
02796       _Type __arg = (__begin[__k % __n]
02797              ^ __begin[(__k + __p) % __n]
02798              ^ __begin[(__k - 1) % __n]);
02799       _Type __r1 = __arg ^ (__arg >> 27);
02800       __r1 = __detail::__mod<_Type,
02801             __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
02802       _Type __r2 = __r1;
02803       if (__k == 0)
02804         __r2 += __s;
02805       else if (__k <= __s)
02806         __r2 += __k % __n + _M_v[__k - 1];
02807       else
02808         __r2 += __k % __n;
02809       __r2 = __detail::__mod<_Type,
02810                __detail::_Shift<_Type, 32>::__value>(__r2);
02811       __begin[(__k + __p) % __n] += __r1;
02812       __begin[(__k + __q) % __n] += __r2;
02813       __begin[__k % __n] = __r2;
02814     }
02815 
02816       for (size_t __k = __m; __k < __m + __n; ++__k)
02817     {
02818       _Type __arg = (__begin[__k % __n]
02819              + __begin[(__k + __p) % __n]
02820              + __begin[(__k - 1) % __n]);
02821       _Type __r3 = __arg ^ (__arg >> 27);
02822       __r3 = __detail::__mod<_Type,
02823            __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
02824       _Type __r4 = __r3 - __k % __n;
02825       __r4 = __detail::__mod<_Type,
02826                __detail::_Shift<_Type, 32>::__value>(__r4);
02827       __begin[(__k + __p) % __n] ^= __r3;
02828       __begin[(__k + __q) % __n] ^= __r4;
02829       __begin[__k % __n] = __r4;
02830     }
02831     }
02832 
02833   template<typename _RealType, size_t __bits,
02834        typename _UniformRandomNumberGenerator>
02835     _RealType
02836     generate_canonical(_UniformRandomNumberGenerator& __urng)
02837     {
02838       const size_t __b
02839     = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
02840                    __bits);
02841       const long double __r = static_cast<long double>(__urng.max())
02842                 - static_cast<long double>(__urng.min()) + 1.0L;
02843       const size_t __log2r = std::log(__r) / std::log(2.0L);
02844       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
02845       _RealType __sum = _RealType(0);
02846       _RealType __tmp = _RealType(1);
02847       for (; __k != 0; --__k)
02848     {
02849       __sum += _RealType(__urng() - __urng.min()) * __tmp;
02850       __tmp *= __r;
02851     }
02852       return __sum / __tmp;
02853     }
02854 
02855 _GLIBCXX_END_NAMESPACE_VERSION
02856 } // namespace
02857 
02858 #endif