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 - 1; ++__i)
00477     __os << __x._M_x[__i] << __space;
00478       __os << __x._M_x[__n - 1];
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 
00504       __is.flags(__flags);
00505       return __is;
00506     }
00507 
00508 
00509   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00510     constexpr size_t
00511     subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
00512 
00513   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00514     constexpr size_t
00515     subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
00516 
00517   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00518     constexpr size_t
00519     subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
00520 
00521   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00522     constexpr _UIntType
00523     subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
00524 
00525   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00526     void
00527     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00528     seed(result_type __value)
00529     {
00530       std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
00531     __lcg(__value == 0u ? default_seed : __value);
00532 
00533       const size_t __n = (__w + 31) / 32;
00534 
00535       for (size_t __i = 0; __i < long_lag; ++__i)
00536     {
00537       _UIntType __sum = 0u;
00538       _UIntType __factor = 1u;
00539       for (size_t __j = 0; __j < __n; ++__j)
00540         {
00541           __sum += __detail::__mod<uint_least32_t,
00542                __detail::_Shift<uint_least32_t, 32>::__value>
00543              (__lcg()) * __factor;
00544           __factor *= __detail::_Shift<_UIntType, 32>::__value;
00545         }
00546       _M_x[__i] = __detail::__mod<_UIntType,
00547         __detail::_Shift<_UIntType, __w>::__value>(__sum);
00548     }
00549       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
00550       _M_p = 0;
00551     }
00552 
00553   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00554     template<typename _Sseq>
00555       typename std::enable_if<std::is_class<_Sseq>::value>::type
00556       subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00557       seed(_Sseq& __q)
00558       {
00559     const size_t __k = (__w + 31) / 32;
00560     uint_least32_t __arr[__r * __k];
00561     __q.generate(__arr + 0, __arr + __r * __k);
00562 
00563     for (size_t __i = 0; __i < long_lag; ++__i)
00564       {
00565         _UIntType __sum = 0u;
00566         _UIntType __factor = 1u;
00567         for (size_t __j = 0; __j < __k; ++__j)
00568           {
00569         __sum += __arr[__k * __i + __j] * __factor;
00570         __factor *= __detail::_Shift<_UIntType, 32>::__value;
00571           }
00572         _M_x[__i] = __detail::__mod<_UIntType,
00573           __detail::_Shift<_UIntType, __w>::__value>(__sum);
00574       }
00575     _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
00576     _M_p = 0;
00577       }
00578 
00579   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00580     typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00581          result_type
00582     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00583     operator()()
00584     {
00585       // Derive short lag index from current index.
00586       long __ps = _M_p - short_lag;
00587       if (__ps < 0)
00588     __ps += long_lag;
00589 
00590       // Calculate new x(i) without overflow or division.
00591       // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
00592       // cannot overflow.
00593       _UIntType __xi;
00594       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
00595     {
00596       __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
00597       _M_carry = 0;
00598     }
00599       else
00600     {
00601       __xi = (__detail::_Shift<_UIntType, __w>::__value
00602           - _M_x[_M_p] - _M_carry + _M_x[__ps]);
00603       _M_carry = 1;
00604     }
00605       _M_x[_M_p] = __xi;
00606 
00607       // Adjust current index to loop around in ring buffer.
00608       if (++_M_p >= long_lag)
00609     _M_p = 0;
00610 
00611       return __xi;
00612     }
00613 
00614   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
00615        typename _CharT, typename _Traits>
00616     std::basic_ostream<_CharT, _Traits>&
00617     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00618            const subtract_with_carry_engine<_UIntType,
00619                         __w, __s, __r>& __x)
00620     {
00621       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00622       typedef typename __ostream_type::ios_base    __ios_base;
00623 
00624       const typename __ios_base::fmtflags __flags = __os.flags();
00625       const _CharT __fill = __os.fill();
00626       const _CharT __space = __os.widen(' ');
00627       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00628       __os.fill(__space);
00629 
00630       for (size_t __i = 0; __i < __r; ++__i)
00631     __os << __x._M_x[__i] << __space;
00632       __os << __x._M_carry;
00633 
00634       __os.flags(__flags);
00635       __os.fill(__fill);
00636       return __os;
00637     }
00638 
00639   template<typename _UIntType, size_t __w, size_t __s, size_t __r,
00640        typename _CharT, typename _Traits>
00641     std::basic_istream<_CharT, _Traits>&
00642     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00643            subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
00644     {
00645       typedef std::basic_ostream<_CharT, _Traits>  __istream_type;
00646       typedef typename __istream_type::ios_base    __ios_base;
00647 
00648       const typename __ios_base::fmtflags __flags = __is.flags();
00649       __is.flags(__ios_base::dec | __ios_base::skipws);
00650 
00651       for (size_t __i = 0; __i < __r; ++__i)
00652     __is >> __x._M_x[__i];
00653       __is >> __x._M_carry;
00654 
00655       __is.flags(__flags);
00656       return __is;
00657     }
00658 
00659 
00660   template<typename _RandomNumberEngine, size_t __p, size_t __r>
00661     constexpr size_t
00662     discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
00663 
00664   template<typename _RandomNumberEngine, size_t __p, size_t __r>
00665     constexpr size_t
00666     discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
00667 
00668   template<typename _RandomNumberEngine, size_t __p, size_t __r>
00669     typename discard_block_engine<_RandomNumberEngine,
00670                __p, __r>::result_type
00671     discard_block_engine<_RandomNumberEngine, __p, __r>::
00672     operator()()
00673     {
00674       if (_M_n >= used_block)
00675     {
00676       _M_b.discard(block_size - _M_n);
00677       _M_n = 0;
00678     }
00679       ++_M_n;
00680       return _M_b();
00681     }
00682 
00683   template<typename _RandomNumberEngine, size_t __p, size_t __r,
00684        typename _CharT, typename _Traits>
00685     std::basic_ostream<_CharT, _Traits>&
00686     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00687            const discard_block_engine<_RandomNumberEngine,
00688            __p, __r>& __x)
00689     {
00690       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00691       typedef typename __ostream_type::ios_base    __ios_base;
00692 
00693       const typename __ios_base::fmtflags __flags = __os.flags();
00694       const _CharT __fill = __os.fill();
00695       const _CharT __space = __os.widen(' ');
00696       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00697       __os.fill(__space);
00698 
00699       __os << __x.base() << __space << __x._M_n;
00700 
00701       __os.flags(__flags);
00702       __os.fill(__fill);
00703       return __os;
00704     }
00705 
00706   template<typename _RandomNumberEngine, size_t __p, size_t __r,
00707        typename _CharT, typename _Traits>
00708     std::basic_istream<_CharT, _Traits>&
00709     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00710            discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
00711     {
00712       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00713       typedef typename __istream_type::ios_base    __ios_base;
00714 
00715       const typename __ios_base::fmtflags __flags = __is.flags();
00716       __is.flags(__ios_base::dec | __ios_base::skipws);
00717 
00718       __is >> __x._M_b >> __x._M_n;
00719 
00720       __is.flags(__flags);
00721       return __is;
00722     }
00723 
00724 
00725   template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
00726     typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
00727       result_type
00728     independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
00729     operator()()
00730     {
00731       const long double __r = static_cast<long double>(_M_b.max())
00732                 - static_cast<long double>(_M_b.min()) + 1.0L;
00733       const result_type __m = std::log(__r) / std::log(2.0L);
00734       result_type __n, __n0, __y0, __y1, __s0, __s1;
00735       for (size_t __i = 0; __i < 2; ++__i)
00736     {
00737       __n = (__w + __m - 1) / __m + __i;
00738       __n0 = __n - __w % __n;
00739       const result_type __w0 = __w / __n;
00740       const result_type __w1 = __w0 + 1;
00741       __s0 = result_type(1) << __w0;
00742       __s1 = result_type(1) << __w1;
00743       __y0 = __s0 * (__r / __s0);
00744       __y1 = __s1 * (__r / __s1);
00745       if (__r - __y0 <= __y0 / __n)
00746         break;
00747     }
00748 
00749       result_type __sum = 0;
00750       for (size_t __k = 0; __k < __n0; ++__k)
00751     {
00752       result_type __u;
00753       do
00754         __u = _M_b() - _M_b.min();
00755       while (__u >= __y0);
00756       __sum = __s0 * __sum + __u % __s0;
00757     }
00758       for (size_t __k = __n0; __k < __n; ++__k)
00759     {
00760       result_type __u;
00761       do
00762         __u = _M_b() - _M_b.min();
00763       while (__u >= __y1);
00764       __sum = __s1 * __sum + __u % __s1;
00765     }
00766       return __sum;
00767     }
00768 
00769 
00770   template<typename _RandomNumberEngine, size_t __k>
00771     constexpr size_t
00772     shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
00773 
00774   template<typename _RandomNumberEngine, size_t __k>
00775     typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
00776     shuffle_order_engine<_RandomNumberEngine, __k>::
00777     operator()()
00778     {
00779       size_t __j = __k * ((_M_y - _M_b.min())
00780               / (_M_b.max() - _M_b.min() + 1.0L));
00781       _M_y = _M_v[__j];
00782       _M_v[__j] = _M_b();
00783 
00784       return _M_y;
00785     }
00786 
00787   template<typename _RandomNumberEngine, size_t __k,
00788        typename _CharT, typename _Traits>
00789     std::basic_ostream<_CharT, _Traits>&
00790     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00791            const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
00792     {
00793       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00794       typedef typename __ostream_type::ios_base    __ios_base;
00795 
00796       const typename __ios_base::fmtflags __flags = __os.flags();
00797       const _CharT __fill = __os.fill();
00798       const _CharT __space = __os.widen(' ');
00799       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00800       __os.fill(__space);
00801 
00802       __os << __x.base();
00803       for (size_t __i = 0; __i < __k; ++__i)
00804     __os << __space << __x._M_v[__i];
00805       __os << __space << __x._M_y;
00806 
00807       __os.flags(__flags);
00808       __os.fill(__fill);
00809       return __os;
00810     }
00811 
00812   template<typename _RandomNumberEngine, size_t __k,
00813        typename _CharT, typename _Traits>
00814     std::basic_istream<_CharT, _Traits>&
00815     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00816            shuffle_order_engine<_RandomNumberEngine, __k>& __x)
00817     {
00818       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00819       typedef typename __istream_type::ios_base    __ios_base;
00820 
00821       const typename __ios_base::fmtflags __flags = __is.flags();
00822       __is.flags(__ios_base::dec | __ios_base::skipws);
00823 
00824       __is >> __x._M_b;
00825       for (size_t __i = 0; __i < __k; ++__i)
00826     __is >> __x._M_v[__i];
00827       __is >> __x._M_y;
00828 
00829       __is.flags(__flags);
00830       return __is;
00831     }
00832 
00833 
00834   template<typename _IntType>
00835     template<typename _UniformRandomNumberGenerator>
00836       typename uniform_int_distribution<_IntType>::result_type
00837       uniform_int_distribution<_IntType>::
00838       operator()(_UniformRandomNumberGenerator& __urng,
00839          const param_type& __param)
00840       {
00841     typedef typename std::make_unsigned<typename
00842       _UniformRandomNumberGenerator::result_type>::type __urngtype;
00843     typedef typename std::make_unsigned<result_type>::type __utype;
00844     typedef typename std::conditional<(sizeof(__urngtype)
00845                        > sizeof(__utype)),
00846       __urngtype, __utype>::type __uctype;
00847 
00848     const __uctype __urngmin = __urng.min();
00849     const __uctype __urngmax = __urng.max();
00850     const __uctype __urngrange = __urngmax - __urngmin;
00851     const __uctype __urange
00852       = __uctype(__param.b()) - __uctype(__param.a());
00853 
00854     __uctype __ret;
00855 
00856     if (__urngrange > __urange)
00857       {
00858         // downscaling
00859         const __uctype __uerange = __urange + 1; // __urange can be zero
00860         const __uctype __scaling = __urngrange / __uerange;
00861         const __uctype __past = __uerange * __scaling;
00862         do
00863           __ret = __uctype(__urng()) - __urngmin;
00864         while (__ret >= __past);
00865         __ret /= __scaling;
00866       }
00867     else if (__urngrange < __urange)
00868       {
00869         // upscaling
00870         /*
00871           Note that every value in [0, urange]
00872           can be written uniquely as
00873 
00874           (urngrange + 1) * high + low
00875 
00876           where
00877 
00878           high in [0, urange / (urngrange + 1)]
00879 
00880           and
00881     
00882           low in [0, urngrange].
00883         */
00884         __uctype __tmp; // wraparound control
00885         do
00886           {
00887         const __uctype __uerngrange = __urngrange + 1;
00888         __tmp = (__uerngrange * operator()
00889              (__urng, param_type(0, __urange / __uerngrange)));
00890         __ret = __tmp + (__uctype(__urng()) - __urngmin);
00891           }
00892         while (__ret > __urange || __ret < __tmp);
00893       }
00894     else
00895       __ret = __uctype(__urng()) - __urngmin;
00896 
00897     return __ret + __param.a();
00898       }
00899 
00900   template<typename _IntType, typename _CharT, typename _Traits>
00901     std::basic_ostream<_CharT, _Traits>&
00902     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00903            const uniform_int_distribution<_IntType>& __x)
00904     {
00905       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00906       typedef typename __ostream_type::ios_base    __ios_base;
00907 
00908       const typename __ios_base::fmtflags __flags = __os.flags();
00909       const _CharT __fill = __os.fill();
00910       const _CharT __space = __os.widen(' ');
00911       __os.flags(__ios_base::scientific | __ios_base::left);
00912       __os.fill(__space);
00913 
00914       __os << __x.a() << __space << __x.b();
00915 
00916       __os.flags(__flags);
00917       __os.fill(__fill);
00918       return __os;
00919     }
00920 
00921   template<typename _IntType, typename _CharT, typename _Traits>
00922     std::basic_istream<_CharT, _Traits>&
00923     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00924            uniform_int_distribution<_IntType>& __x)
00925     {
00926       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00927       typedef typename __istream_type::ios_base    __ios_base;
00928 
00929       const typename __ios_base::fmtflags __flags = __is.flags();
00930       __is.flags(__ios_base::dec | __ios_base::skipws);
00931 
00932       _IntType __a, __b;
00933       __is >> __a >> __b;
00934       __x.param(typename uniform_int_distribution<_IntType>::
00935         param_type(__a, __b));
00936 
00937       __is.flags(__flags);
00938       return __is;
00939     }
00940 
00941 
00942   template<typename _RealType, typename _CharT, typename _Traits>
00943     std::basic_ostream<_CharT, _Traits>&
00944     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00945            const uniform_real_distribution<_RealType>& __x)
00946     {
00947       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00948       typedef typename __ostream_type::ios_base    __ios_base;
00949 
00950       const typename __ios_base::fmtflags __flags = __os.flags();
00951       const _CharT __fill = __os.fill();
00952       const std::streamsize __precision = __os.precision();
00953       const _CharT __space = __os.widen(' ');
00954       __os.flags(__ios_base::scientific | __ios_base::left);
00955       __os.fill(__space);
00956       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00957 
00958       __os << __x.a() << __space << __x.b();
00959 
00960       __os.flags(__flags);
00961       __os.fill(__fill);
00962       __os.precision(__precision);
00963       return __os;
00964     }
00965 
00966   template<typename _RealType, typename _CharT, typename _Traits>
00967     std::basic_istream<_CharT, _Traits>&
00968     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00969            uniform_real_distribution<_RealType>& __x)
00970     {
00971       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00972       typedef typename __istream_type::ios_base    __ios_base;
00973 
00974       const typename __ios_base::fmtflags __flags = __is.flags();
00975       __is.flags(__ios_base::skipws);
00976 
00977       _RealType __a, __b;
00978       __is >> __a >> __b;
00979       __x.param(typename uniform_real_distribution<_RealType>::
00980         param_type(__a, __b));
00981 
00982       __is.flags(__flags);
00983       return __is;
00984     }
00985 
00986 
00987   template<typename _CharT, typename _Traits>
00988     std::basic_ostream<_CharT, _Traits>&
00989     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00990            const bernoulli_distribution& __x)
00991     {
00992       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00993       typedef typename __ostream_type::ios_base    __ios_base;
00994 
00995       const typename __ios_base::fmtflags __flags = __os.flags();
00996       const _CharT __fill = __os.fill();
00997       const std::streamsize __precision = __os.precision();
00998       __os.flags(__ios_base::scientific | __ios_base::left);
00999       __os.fill(__os.widen(' '));
01000       __os.precision(std::numeric_limits<double>::max_digits10);
01001 
01002       __os << __x.p();
01003 
01004       __os.flags(__flags);
01005       __os.fill(__fill);
01006       __os.precision(__precision);
01007       return __os;
01008     }
01009 
01010 
01011   template<typename _IntType>
01012     template<typename _UniformRandomNumberGenerator>
01013       typename geometric_distribution<_IntType>::result_type
01014       geometric_distribution<_IntType>::
01015       operator()(_UniformRandomNumberGenerator& __urng,
01016          const param_type& __param)
01017       {
01018     // About the epsilon thing see this thread:
01019     // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
01020     const double __naf =
01021       (1 - std::numeric_limits<double>::epsilon()) / 2;
01022     // The largest _RealType convertible to _IntType.
01023     const double __thr =
01024       std::numeric_limits<_IntType>::max() + __naf;
01025     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01026       __aurng(__urng);
01027 
01028     double __cand;
01029     do
01030       __cand = std::floor(std::log(__aurng()) / __param._M_log_1_p);
01031     while (__cand >= __thr);
01032 
01033     return result_type(__cand + __naf);
01034       }
01035 
01036   template<typename _IntType,
01037        typename _CharT, typename _Traits>
01038     std::basic_ostream<_CharT, _Traits>&
01039     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01040            const geometric_distribution<_IntType>& __x)
01041     {
01042       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01043       typedef typename __ostream_type::ios_base    __ios_base;
01044 
01045       const typename __ios_base::fmtflags __flags = __os.flags();
01046       const _CharT __fill = __os.fill();
01047       const std::streamsize __precision = __os.precision();
01048       __os.flags(__ios_base::scientific | __ios_base::left);
01049       __os.fill(__os.widen(' '));
01050       __os.precision(std::numeric_limits<double>::max_digits10);
01051 
01052       __os << __x.p();
01053 
01054       __os.flags(__flags);
01055       __os.fill(__fill);
01056       __os.precision(__precision);
01057       return __os;
01058     }
01059 
01060   template<typename _IntType,
01061        typename _CharT, typename _Traits>
01062     std::basic_istream<_CharT, _Traits>&
01063     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01064            geometric_distribution<_IntType>& __x)
01065     {
01066       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01067       typedef typename __istream_type::ios_base    __ios_base;
01068 
01069       const typename __ios_base::fmtflags __flags = __is.flags();
01070       __is.flags(__ios_base::skipws);
01071 
01072       double __p;
01073       __is >> __p;
01074       __x.param(typename geometric_distribution<_IntType>::param_type(__p));
01075 
01076       __is.flags(__flags);
01077       return __is;
01078     }
01079 
01080 
01081   template<typename _IntType>
01082     template<typename _UniformRandomNumberGenerator>
01083       typename negative_binomial_distribution<_IntType>::result_type
01084       negative_binomial_distribution<_IntType>::
01085       operator()(_UniformRandomNumberGenerator& __urng)
01086       {
01087     const double __y = _M_gd(__urng);
01088 
01089     // XXX Is the constructor too slow?
01090     std::poisson_distribution<result_type> __poisson(__y);
01091     return __poisson(__urng);
01092       }
01093 
01094   template<typename _IntType>
01095     template<typename _UniformRandomNumberGenerator>
01096       typename negative_binomial_distribution<_IntType>::result_type
01097       negative_binomial_distribution<_IntType>::
01098       operator()(_UniformRandomNumberGenerator& __urng,
01099          const param_type& __p)
01100       {
01101     typedef typename std::gamma_distribution<result_type>::param_type
01102       param_type;
01103     
01104     const double __y =
01105       _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
01106 
01107     std::poisson_distribution<result_type> __poisson(__y);
01108     return __poisson(__urng);
01109       }
01110 
01111   template<typename _IntType, typename _CharT, typename _Traits>
01112     std::basic_ostream<_CharT, _Traits>&
01113     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01114            const negative_binomial_distribution<_IntType>& __x)
01115     {
01116       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01117       typedef typename __ostream_type::ios_base    __ios_base;
01118 
01119       const typename __ios_base::fmtflags __flags = __os.flags();
01120       const _CharT __fill = __os.fill();
01121       const std::streamsize __precision = __os.precision();
01122       const _CharT __space = __os.widen(' ');
01123       __os.flags(__ios_base::scientific | __ios_base::left);
01124       __os.fill(__os.widen(' '));
01125       __os.precision(std::numeric_limits<double>::max_digits10);
01126 
01127       __os << __x.k() << __space << __x.p()
01128        << __space << __x._M_gd;
01129 
01130       __os.flags(__flags);
01131       __os.fill(__fill);
01132       __os.precision(__precision);
01133       return __os;
01134     }
01135 
01136   template<typename _IntType, typename _CharT, typename _Traits>
01137     std::basic_istream<_CharT, _Traits>&
01138     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01139            negative_binomial_distribution<_IntType>& __x)
01140     {
01141       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01142       typedef typename __istream_type::ios_base    __ios_base;
01143 
01144       const typename __ios_base::fmtflags __flags = __is.flags();
01145       __is.flags(__ios_base::skipws);
01146 
01147       _IntType __k;
01148       double __p;
01149       __is >> __k >> __p >> __x._M_gd;
01150       __x.param(typename negative_binomial_distribution<_IntType>::
01151         param_type(__k, __p));
01152 
01153       __is.flags(__flags);
01154       return __is;
01155     }
01156 
01157 
01158   template<typename _IntType>
01159     void
01160     poisson_distribution<_IntType>::param_type::
01161     _M_initialize()
01162     {
01163 #if _GLIBCXX_USE_C99_MATH_TR1
01164       if (_M_mean >= 12)
01165     {
01166       const double __m = std::floor(_M_mean);
01167       _M_lm_thr = std::log(_M_mean);
01168       _M_lfm = std::lgamma(__m + 1);
01169       _M_sm = std::sqrt(__m);
01170 
01171       const double __pi_4 = 0.7853981633974483096156608458198757L;
01172       const double __dx = std::sqrt(2 * __m * std::log(32 * __m
01173                                   / __pi_4));
01174       _M_d = std::round(std::max(6.0, std::min(__m, __dx)));
01175       const double __cx = 2 * __m + _M_d;
01176       _M_scx = std::sqrt(__cx / 2);
01177       _M_1cx = 1 / __cx;
01178 
01179       _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
01180       _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
01181         / _M_d;
01182     }
01183       else
01184 #endif
01185     _M_lm_thr = std::exp(-_M_mean);
01186       }
01187 
01188   /**
01189    * A rejection algorithm when mean >= 12 and a simple method based
01190    * upon the multiplication of uniform random variates otherwise.
01191    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
01192    * is defined.
01193    *
01194    * Reference:
01195    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01196    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
01197    */
01198   template<typename _IntType>
01199     template<typename _UniformRandomNumberGenerator>
01200       typename poisson_distribution<_IntType>::result_type
01201       poisson_distribution<_IntType>::
01202       operator()(_UniformRandomNumberGenerator& __urng,
01203          const param_type& __param)
01204       {
01205     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01206       __aurng(__urng);
01207 #if _GLIBCXX_USE_C99_MATH_TR1
01208     if (__param.mean() >= 12)
01209       {
01210         double __x;
01211 
01212         // See comments above...
01213         const double __naf =
01214           (1 - std::numeric_limits<double>::epsilon()) / 2;
01215         const double __thr =
01216           std::numeric_limits<_IntType>::max() + __naf;
01217 
01218         const double __m = std::floor(__param.mean());
01219         // sqrt(pi / 2)
01220         const double __spi_2 = 1.2533141373155002512078826424055226L;
01221         const double __c1 = __param._M_sm * __spi_2;
01222         const double __c2 = __param._M_c2b + __c1;
01223         const double __c3 = __c2 + 1;
01224         const double __c4 = __c3 + 1;
01225         // e^(1 / 78)
01226         const double __e178 = 1.0129030479320018583185514777512983L;
01227         const double __c5 = __c4 + __e178;
01228         const double __c = __param._M_cb + __c5;
01229         const double __2cx = 2 * (2 * __m + __param._M_d);
01230 
01231         bool __reject = true;
01232         do
01233           {
01234         const double __u = __c * __aurng();
01235         const double __e = -std::log(__aurng());
01236 
01237         double __w = 0.0;
01238 
01239         if (__u <= __c1)
01240           {
01241             const double __n = _M_nd(__urng);
01242             const double __y = -std::abs(__n) * __param._M_sm - 1;
01243             __x = std::floor(__y);
01244             __w = -__n * __n / 2;
01245             if (__x < -__m)
01246               continue;
01247           }
01248         else if (__u <= __c2)
01249           {
01250             const double __n = _M_nd(__urng);
01251             const double __y = 1 + std::abs(__n) * __param._M_scx;
01252             __x = std::ceil(__y);
01253             __w = __y * (2 - __y) * __param._M_1cx;
01254             if (__x > __param._M_d)
01255               continue;
01256           }
01257         else if (__u <= __c3)
01258           // NB: This case not in the book, nor in the Errata,
01259           // but should be ok...
01260           __x = -1;
01261         else if (__u <= __c4)
01262           __x = 0;
01263         else if (__u <= __c5)
01264           __x = 1;
01265         else
01266           {
01267             const double __v = -std::log(__aurng());
01268             const double __y = __param._M_d
01269                      + __v * __2cx / __param._M_d;
01270             __x = std::ceil(__y);
01271             __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
01272           }
01273 
01274         __reject = (__w - __e - __x * __param._M_lm_thr
01275                 > __param._M_lfm - std::lgamma(__x + __m + 1));
01276 
01277         __reject |= __x + __m >= __thr;
01278 
01279           } while (__reject);
01280 
01281         return result_type(__x + __m + __naf);
01282       }
01283     else
01284 #endif
01285       {
01286         _IntType     __x = 0;
01287         double __prod = 1.0;
01288 
01289         do
01290           {
01291         __prod *= __aurng();
01292         __x += 1;
01293           }
01294         while (__prod > __param._M_lm_thr);
01295 
01296         return __x - 1;
01297       }
01298       }
01299 
01300   template<typename _IntType,
01301        typename _CharT, typename _Traits>
01302     std::basic_ostream<_CharT, _Traits>&
01303     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01304            const poisson_distribution<_IntType>& __x)
01305     {
01306       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01307       typedef typename __ostream_type::ios_base    __ios_base;
01308 
01309       const typename __ios_base::fmtflags __flags = __os.flags();
01310       const _CharT __fill = __os.fill();
01311       const std::streamsize __precision = __os.precision();
01312       const _CharT __space = __os.widen(' ');
01313       __os.flags(__ios_base::scientific | __ios_base::left);
01314       __os.fill(__space);
01315       __os.precision(std::numeric_limits<double>::max_digits10);
01316 
01317       __os << __x.mean() << __space << __x._M_nd;
01318 
01319       __os.flags(__flags);
01320       __os.fill(__fill);
01321       __os.precision(__precision);
01322       return __os;
01323     }
01324 
01325   template<typename _IntType,
01326        typename _CharT, typename _Traits>
01327     std::basic_istream<_CharT, _Traits>&
01328     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01329            poisson_distribution<_IntType>& __x)
01330     {
01331       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01332       typedef typename __istream_type::ios_base    __ios_base;
01333 
01334       const typename __ios_base::fmtflags __flags = __is.flags();
01335       __is.flags(__ios_base::skipws);
01336 
01337       double __mean;
01338       __is >> __mean >> __x._M_nd;
01339       __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
01340 
01341       __is.flags(__flags);
01342       return __is;
01343     }
01344 
01345 
01346   template<typename _IntType>
01347     void
01348     binomial_distribution<_IntType>::param_type::
01349     _M_initialize()
01350     {
01351       const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
01352 
01353       _M_easy = true;
01354 
01355 #if _GLIBCXX_USE_C99_MATH_TR1
01356       if (_M_t * __p12 >= 8)
01357     {
01358       _M_easy = false;
01359       const double __np = std::floor(_M_t * __p12);
01360       const double __pa = __np / _M_t;
01361       const double __1p = 1 - __pa;
01362 
01363       const double __pi_4 = 0.7853981633974483096156608458198757L;
01364       const double __d1x =
01365         std::sqrt(__np * __1p * std::log(32 * __np
01366                          / (81 * __pi_4 * __1p)));
01367       _M_d1 = std::round(std::max(1.0, __d1x));
01368       const double __d2x =
01369         std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
01370                          / (__pi_4 * __pa)));
01371       _M_d2 = std::round(std::max(1.0, __d2x));
01372 
01373       // sqrt(pi / 2)
01374       const double __spi_2 = 1.2533141373155002512078826424055226L;
01375       _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
01376       _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
01377       _M_c = 2 * _M_d1 / __np;
01378       _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
01379       const double __a12 = _M_a1 + _M_s2 * __spi_2;
01380       const double __s1s = _M_s1 * _M_s1;
01381       _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
01382                  * 2 * __s1s / _M_d1
01383                  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
01384       const double __s2s = _M_s2 * _M_s2;
01385       _M_s = (_M_a123 + 2 * __s2s / _M_d2
01386           * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
01387       _M_lf = (std::lgamma(__np + 1)
01388            + std::lgamma(_M_t - __np + 1));
01389       _M_lp1p = std::log(__pa / __1p);
01390 
01391       _M_q = -std::log(1 - (__p12 - __pa) / __1p);
01392     }
01393       else
01394 #endif
01395     _M_q = -std::log(1 - __p12);
01396     }
01397 
01398   template<typename _IntType>
01399     template<typename _UniformRandomNumberGenerator>
01400       typename binomial_distribution<_IntType>::result_type
01401       binomial_distribution<_IntType>::
01402       _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
01403       {
01404     _IntType __x = 0;
01405     double __sum = 0.0;
01406     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01407       __aurng(__urng);
01408 
01409     do
01410       {
01411         const double __e = -std::log(__aurng());
01412         __sum += __e / (__t - __x);
01413         __x += 1;
01414       }
01415     while (__sum <= _M_param._M_q);
01416 
01417     return __x - 1;
01418       }
01419 
01420   /**
01421    * A rejection algorithm when t * p >= 8 and a simple waiting time
01422    * method - the second in the referenced book - otherwise.
01423    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
01424    * is defined.
01425    *
01426    * Reference:
01427    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01428    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
01429    */
01430   template<typename _IntType>
01431     template<typename _UniformRandomNumberGenerator>
01432       typename binomial_distribution<_IntType>::result_type
01433       binomial_distribution<_IntType>::
01434       operator()(_UniformRandomNumberGenerator& __urng,
01435          const param_type& __param)
01436       {
01437     result_type __ret;
01438     const _IntType __t = __param.t();
01439     const double __p = __param.p();
01440     const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
01441     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01442       __aurng(__urng);
01443 
01444 #if _GLIBCXX_USE_C99_MATH_TR1
01445     if (!__param._M_easy)
01446       {
01447         double __x;
01448 
01449         // See comments above...
01450         const double __naf =
01451           (1 - std::numeric_limits<double>::epsilon()) / 2;
01452         const double __thr =
01453           std::numeric_limits<_IntType>::max() + __naf;
01454 
01455         const double __np = std::floor(__t * __p12);
01456 
01457         // sqrt(pi / 2)
01458         const double __spi_2 = 1.2533141373155002512078826424055226L;
01459         const double __a1 = __param._M_a1;
01460         const double __a12 = __a1 + __param._M_s2 * __spi_2;
01461         const double __a123 = __param._M_a123;
01462         const double __s1s = __param._M_s1 * __param._M_s1;
01463         const double __s2s = __param._M_s2 * __param._M_s2;
01464 
01465         bool __reject;
01466         do
01467           {
01468         const double __u = __param._M_s * __aurng();
01469 
01470         double __v;
01471 
01472         if (__u <= __a1)
01473           {
01474             const double __n = _M_nd(__urng);
01475             const double __y = __param._M_s1 * std::abs(__n);
01476             __reject = __y >= __param._M_d1;
01477             if (!__reject)
01478               {
01479             const double __e = -std::log(__aurng());
01480             __x = std::floor(__y);
01481             __v = -__e - __n * __n / 2 + __param._M_c;
01482               }
01483           }
01484         else if (__u <= __a12)
01485           {
01486             const double __n = _M_nd(__urng);
01487             const double __y = __param._M_s2 * std::abs(__n);
01488             __reject = __y >= __param._M_d2;
01489             if (!__reject)
01490               {
01491             const double __e = -std::log(__aurng());
01492             __x = std::floor(-__y);
01493             __v = -__e - __n * __n / 2;
01494               }
01495           }
01496         else if (__u <= __a123)
01497           {
01498             const double __e1 = -std::log(__aurng());
01499             const double __e2 = -std::log(__aurng());
01500 
01501             const double __y = __param._M_d1
01502                      + 2 * __s1s * __e1 / __param._M_d1;
01503             __x = std::floor(__y);
01504             __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
01505                             -__y / (2 * __s1s)));
01506             __reject = false;
01507           }
01508         else
01509           {
01510             const double __e1 = -std::log(__aurng());
01511             const double __e2 = -std::log(__aurng());
01512 
01513             const double __y = __param._M_d2
01514                      + 2 * __s2s * __e1 / __param._M_d2;
01515             __x = std::floor(-__y);
01516             __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
01517             __reject = false;
01518           }
01519 
01520         __reject = __reject || __x < -__np || __x > __t - __np;
01521         if (!__reject)
01522           {
01523             const double __lfx =
01524               std::lgamma(__np + __x + 1)
01525               + std::lgamma(__t - (__np + __x) + 1);
01526             __reject = __v > __param._M_lf - __lfx
01527                  + __x * __param._M_lp1p;
01528           }
01529 
01530         __reject |= __x + __np >= __thr;
01531           }
01532         while (__reject);
01533 
01534         __x += __np + __naf;
01535 
01536         const _IntType __z = _M_waiting(__urng, __t - _IntType(__x));
01537         __ret = _IntType(__x) + __z;
01538       }
01539     else
01540 #endif
01541       __ret = _M_waiting(__urng, __t);
01542 
01543     if (__p12 != __p)
01544       __ret = __t - __ret;
01545     return __ret;
01546       }
01547 
01548   template<typename _IntType,
01549        typename _CharT, typename _Traits>
01550     std::basic_ostream<_CharT, _Traits>&
01551     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01552            const binomial_distribution<_IntType>& __x)
01553     {
01554       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01555       typedef typename __ostream_type::ios_base    __ios_base;
01556 
01557       const typename __ios_base::fmtflags __flags = __os.flags();
01558       const _CharT __fill = __os.fill();
01559       const std::streamsize __precision = __os.precision();
01560       const _CharT __space = __os.widen(' ');
01561       __os.flags(__ios_base::scientific | __ios_base::left);
01562       __os.fill(__space);
01563       __os.precision(std::numeric_limits<double>::max_digits10);
01564 
01565       __os << __x.t() << __space << __x.p()
01566        << __space << __x._M_nd;
01567 
01568       __os.flags(__flags);
01569       __os.fill(__fill);
01570       __os.precision(__precision);
01571       return __os;
01572     }
01573 
01574   template<typename _IntType,
01575        typename _CharT, typename _Traits>
01576     std::basic_istream<_CharT, _Traits>&
01577     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01578            binomial_distribution<_IntType>& __x)
01579     {
01580       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01581       typedef typename __istream_type::ios_base    __ios_base;
01582 
01583       const typename __ios_base::fmtflags __flags = __is.flags();
01584       __is.flags(__ios_base::dec | __ios_base::skipws);
01585 
01586       _IntType __t;
01587       double __p;
01588       __is >> __t >> __p >> __x._M_nd;
01589       __x.param(typename binomial_distribution<_IntType>::
01590         param_type(__t, __p));
01591 
01592       __is.flags(__flags);
01593       return __is;
01594     }
01595 
01596 
01597   template<typename _RealType, typename _CharT, typename _Traits>
01598     std::basic_ostream<_CharT, _Traits>&
01599     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01600            const exponential_distribution<_RealType>& __x)
01601     {
01602       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01603       typedef typename __ostream_type::ios_base    __ios_base;
01604 
01605       const typename __ios_base::fmtflags __flags = __os.flags();
01606       const _CharT __fill = __os.fill();
01607       const std::streamsize __precision = __os.precision();
01608       __os.flags(__ios_base::scientific | __ios_base::left);
01609       __os.fill(__os.widen(' '));
01610       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01611 
01612       __os << __x.lambda();
01613 
01614       __os.flags(__flags);
01615       __os.fill(__fill);
01616       __os.precision(__precision);
01617       return __os;
01618     }
01619 
01620   template<typename _RealType, typename _CharT, typename _Traits>
01621     std::basic_istream<_CharT, _Traits>&
01622     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01623            exponential_distribution<_RealType>& __x)
01624     {
01625       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01626       typedef typename __istream_type::ios_base    __ios_base;
01627 
01628       const typename __ios_base::fmtflags __flags = __is.flags();
01629       __is.flags(__ios_base::dec | __ios_base::skipws);
01630 
01631       _RealType __lambda;
01632       __is >> __lambda;
01633       __x.param(typename exponential_distribution<_RealType>::
01634         param_type(__lambda));
01635 
01636       __is.flags(__flags);
01637       return __is;
01638     }
01639 
01640 
01641   /**
01642    * Polar method due to Marsaglia.
01643    *
01644    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01645    * New York, 1986, Ch. V, Sect. 4.4.
01646    */
01647   template<typename _RealType>
01648     template<typename _UniformRandomNumberGenerator>
01649       typename normal_distribution<_RealType>::result_type
01650       normal_distribution<_RealType>::
01651       operator()(_UniformRandomNumberGenerator& __urng,
01652          const param_type& __param)
01653       {
01654     result_type __ret;
01655     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01656       __aurng(__urng);
01657 
01658     if (_M_saved_available)
01659       {
01660         _M_saved_available = false;
01661         __ret = _M_saved;
01662       }
01663     else
01664       {
01665         result_type __x, __y, __r2;
01666         do
01667           {
01668         __x = result_type(2.0) * __aurng() - 1.0;
01669         __y = result_type(2.0) * __aurng() - 1.0;
01670         __r2 = __x * __x + __y * __y;
01671           }
01672         while (__r2 > 1.0 || __r2 == 0.0);
01673 
01674         const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
01675         _M_saved = __x * __mult;
01676         _M_saved_available = true;
01677         __ret = __y * __mult;
01678       }
01679 
01680     __ret = __ret * __param.stddev() + __param.mean();
01681     return __ret;
01682       }
01683 
01684   template<typename _RealType>
01685     bool
01686     operator==(const std::normal_distribution<_RealType>& __d1,
01687            const std::normal_distribution<_RealType>& __d2)
01688     {
01689       if (__d1._M_param == __d2._M_param
01690       && __d1._M_saved_available == __d2._M_saved_available)
01691     {
01692       if (__d1._M_saved_available
01693           && __d1._M_saved == __d2._M_saved)
01694         return true;
01695       else if(!__d1._M_saved_available)
01696         return true;
01697       else
01698         return false;
01699     }
01700       else
01701     return false;
01702     }
01703 
01704   template<typename _RealType, typename _CharT, typename _Traits>
01705     std::basic_ostream<_CharT, _Traits>&
01706     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01707            const normal_distribution<_RealType>& __x)
01708     {
01709       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01710       typedef typename __ostream_type::ios_base    __ios_base;
01711 
01712       const typename __ios_base::fmtflags __flags = __os.flags();
01713       const _CharT __fill = __os.fill();
01714       const std::streamsize __precision = __os.precision();
01715       const _CharT __space = __os.widen(' ');
01716       __os.flags(__ios_base::scientific | __ios_base::left);
01717       __os.fill(__space);
01718       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01719 
01720       __os << __x.mean() << __space << __x.stddev()
01721        << __space << __x._M_saved_available;
01722       if (__x._M_saved_available)
01723     __os << __space << __x._M_saved;
01724 
01725       __os.flags(__flags);
01726       __os.fill(__fill);
01727       __os.precision(__precision);
01728       return __os;
01729     }
01730 
01731   template<typename _RealType, typename _CharT, typename _Traits>
01732     std::basic_istream<_CharT, _Traits>&
01733     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01734            normal_distribution<_RealType>& __x)
01735     {
01736       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01737       typedef typename __istream_type::ios_base    __ios_base;
01738 
01739       const typename __ios_base::fmtflags __flags = __is.flags();
01740       __is.flags(__ios_base::dec | __ios_base::skipws);
01741 
01742       double __mean, __stddev;
01743       __is >> __mean >> __stddev
01744        >> __x._M_saved_available;
01745       if (__x._M_saved_available)
01746     __is >> __x._M_saved;
01747       __x.param(typename normal_distribution<_RealType>::
01748         param_type(__mean, __stddev));
01749 
01750       __is.flags(__flags);
01751       return __is;
01752     }
01753 
01754 
01755   template<typename _RealType, typename _CharT, typename _Traits>
01756     std::basic_ostream<_CharT, _Traits>&
01757     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01758            const lognormal_distribution<_RealType>& __x)
01759     {
01760       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01761       typedef typename __ostream_type::ios_base    __ios_base;
01762 
01763       const typename __ios_base::fmtflags __flags = __os.flags();
01764       const _CharT __fill = __os.fill();
01765       const std::streamsize __precision = __os.precision();
01766       const _CharT __space = __os.widen(' ');
01767       __os.flags(__ios_base::scientific | __ios_base::left);
01768       __os.fill(__space);
01769       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01770 
01771       __os << __x.m() << __space << __x.s()
01772        << __space << __x._M_nd;
01773 
01774       __os.flags(__flags);
01775       __os.fill(__fill);
01776       __os.precision(__precision);
01777       return __os;
01778     }
01779 
01780   template<typename _RealType, typename _CharT, typename _Traits>
01781     std::basic_istream<_CharT, _Traits>&
01782     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01783            lognormal_distribution<_RealType>& __x)
01784     {
01785       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01786       typedef typename __istream_type::ios_base    __ios_base;
01787 
01788       const typename __ios_base::fmtflags __flags = __is.flags();
01789       __is.flags(__ios_base::dec | __ios_base::skipws);
01790 
01791       _RealType __m, __s;
01792       __is >> __m >> __s >> __x._M_nd;
01793       __x.param(typename lognormal_distribution<_RealType>::
01794         param_type(__m, __s));
01795 
01796       __is.flags(__flags);
01797       return __is;
01798     }
01799 
01800 
01801   template<typename _RealType, typename _CharT, typename _Traits>
01802     std::basic_ostream<_CharT, _Traits>&
01803     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01804            const chi_squared_distribution<_RealType>& __x)
01805     {
01806       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01807       typedef typename __ostream_type::ios_base    __ios_base;
01808 
01809       const typename __ios_base::fmtflags __flags = __os.flags();
01810       const _CharT __fill = __os.fill();
01811       const std::streamsize __precision = __os.precision();
01812       const _CharT __space = __os.widen(' ');
01813       __os.flags(__ios_base::scientific | __ios_base::left);
01814       __os.fill(__space);
01815       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01816 
01817       __os << __x.n() << __space << __x._M_gd;
01818 
01819       __os.flags(__flags);
01820       __os.fill(__fill);
01821       __os.precision(__precision);
01822       return __os;
01823     }
01824 
01825   template<typename _RealType, typename _CharT, typename _Traits>
01826     std::basic_istream<_CharT, _Traits>&
01827     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01828            chi_squared_distribution<_RealType>& __x)
01829     {
01830       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01831       typedef typename __istream_type::ios_base    __ios_base;
01832 
01833       const typename __ios_base::fmtflags __flags = __is.flags();
01834       __is.flags(__ios_base::dec | __ios_base::skipws);
01835 
01836       _RealType __n;
01837       __is >> __n >> __x._M_gd;
01838       __x.param(typename chi_squared_distribution<_RealType>::
01839         param_type(__n));
01840 
01841       __is.flags(__flags);
01842       return __is;
01843     }
01844 
01845 
01846   template<typename _RealType>
01847     template<typename _UniformRandomNumberGenerator>
01848       typename cauchy_distribution<_RealType>::result_type
01849       cauchy_distribution<_RealType>::
01850       operator()(_UniformRandomNumberGenerator& __urng,
01851          const param_type& __p)
01852       {
01853     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01854       __aurng(__urng);
01855     _RealType __u;
01856     do
01857       __u = __aurng();
01858     while (__u == 0.5);
01859 
01860     const _RealType __pi = 3.1415926535897932384626433832795029L;
01861     return __p.a() + __p.b() * std::tan(__pi * __u);
01862       }
01863 
01864   template<typename _RealType, typename _CharT, typename _Traits>
01865     std::basic_ostream<_CharT, _Traits>&
01866     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01867            const cauchy_distribution<_RealType>& __x)
01868     {
01869       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01870       typedef typename __ostream_type::ios_base    __ios_base;
01871 
01872       const typename __ios_base::fmtflags __flags = __os.flags();
01873       const _CharT __fill = __os.fill();
01874       const std::streamsize __precision = __os.precision();
01875       const _CharT __space = __os.widen(' ');
01876       __os.flags(__ios_base::scientific | __ios_base::left);
01877       __os.fill(__space);
01878       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01879 
01880       __os << __x.a() << __space << __x.b();
01881 
01882       __os.flags(__flags);
01883       __os.fill(__fill);
01884       __os.precision(__precision);
01885       return __os;
01886     }
01887 
01888   template<typename _RealType, typename _CharT, typename _Traits>
01889     std::basic_istream<_CharT, _Traits>&
01890     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01891            cauchy_distribution<_RealType>& __x)
01892     {
01893       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01894       typedef typename __istream_type::ios_base    __ios_base;
01895 
01896       const typename __ios_base::fmtflags __flags = __is.flags();
01897       __is.flags(__ios_base::dec | __ios_base::skipws);
01898 
01899       _RealType __a, __b;
01900       __is >> __a >> __b;
01901       __x.param(typename cauchy_distribution<_RealType>::
01902         param_type(__a, __b));
01903 
01904       __is.flags(__flags);
01905       return __is;
01906     }
01907 
01908 
01909   template<typename _RealType, typename _CharT, typename _Traits>
01910     std::basic_ostream<_CharT, _Traits>&
01911     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01912            const fisher_f_distribution<_RealType>& __x)
01913     {
01914       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01915       typedef typename __ostream_type::ios_base    __ios_base;
01916 
01917       const typename __ios_base::fmtflags __flags = __os.flags();
01918       const _CharT __fill = __os.fill();
01919       const std::streamsize __precision = __os.precision();
01920       const _CharT __space = __os.widen(' ');
01921       __os.flags(__ios_base::scientific | __ios_base::left);
01922       __os.fill(__space);
01923       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01924 
01925       __os << __x.m() << __space << __x.n()
01926        << __space << __x._M_gd_x << __space << __x._M_gd_y;
01927 
01928       __os.flags(__flags);
01929       __os.fill(__fill);
01930       __os.precision(__precision);
01931       return __os;
01932     }
01933 
01934   template<typename _RealType, typename _CharT, typename _Traits>
01935     std::basic_istream<_CharT, _Traits>&
01936     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01937            fisher_f_distribution<_RealType>& __x)
01938     {
01939       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01940       typedef typename __istream_type::ios_base    __ios_base;
01941 
01942       const typename __ios_base::fmtflags __flags = __is.flags();
01943       __is.flags(__ios_base::dec | __ios_base::skipws);
01944 
01945       _RealType __m, __n;
01946       __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
01947       __x.param(typename fisher_f_distribution<_RealType>::
01948         param_type(__m, __n));
01949 
01950       __is.flags(__flags);
01951       return __is;
01952     }
01953 
01954 
01955   template<typename _RealType, typename _CharT, typename _Traits>
01956     std::basic_ostream<_CharT, _Traits>&
01957     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01958            const student_t_distribution<_RealType>& __x)
01959     {
01960       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01961       typedef typename __ostream_type::ios_base    __ios_base;
01962 
01963       const typename __ios_base::fmtflags __flags = __os.flags();
01964       const _CharT __fill = __os.fill();
01965       const std::streamsize __precision = __os.precision();
01966       const _CharT __space = __os.widen(' ');
01967       __os.flags(__ios_base::scientific | __ios_base::left);
01968       __os.fill(__space);
01969       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01970 
01971       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
01972 
01973       __os.flags(__flags);
01974       __os.fill(__fill);
01975       __os.precision(__precision);
01976       return __os;
01977     }
01978 
01979   template<typename _RealType, typename _CharT, typename _Traits>
01980     std::basic_istream<_CharT, _Traits>&
01981     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01982            student_t_distribution<_RealType>& __x)
01983     {
01984       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01985       typedef typename __istream_type::ios_base    __ios_base;
01986 
01987       const typename __ios_base::fmtflags __flags = __is.flags();
01988       __is.flags(__ios_base::dec | __ios_base::skipws);
01989 
01990       _RealType __n;
01991       __is >> __n >> __x._M_nd >> __x._M_gd;
01992       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
01993 
01994       __is.flags(__flags);
01995       return __is;
01996     }
01997 
01998 
01999   template<typename _RealType>
02000     void
02001     gamma_distribution<_RealType>::param_type::
02002     _M_initialize()
02003     {
02004       _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
02005 
02006       const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
02007       _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
02008     }
02009 
02010   /**
02011    * Marsaglia, G. and Tsang, W. W.
02012    * "A Simple Method for Generating Gamma Variables"
02013    * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
02014    */
02015   template<typename _RealType>
02016     template<typename _UniformRandomNumberGenerator>
02017       typename gamma_distribution<_RealType>::result_type
02018       gamma_distribution<_RealType>::
02019       operator()(_UniformRandomNumberGenerator& __urng,
02020          const param_type& __param)
02021       {
02022     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02023       __aurng(__urng);
02024 
02025     result_type __u, __v, __n;
02026     const result_type __a1 = (__param._M_malpha
02027                   - _RealType(1.0) / _RealType(3.0));
02028 
02029     do
02030       {
02031         do
02032           {
02033         __n = _M_nd(__urng);
02034         __v = result_type(1.0) + __param._M_a2 * __n; 
02035           }
02036         while (__v <= 0.0);
02037 
02038         __v = __v * __v * __v;
02039         __u = __aurng();
02040       }
02041     while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
02042            && (std::log(__u) > (0.5 * __n * __n + __a1
02043                     * (1.0 - __v + std::log(__v)))));
02044 
02045     if (__param.alpha() == __param._M_malpha)
02046       return __a1 * __v * __param.beta();
02047     else
02048       {
02049         do
02050           __u = __aurng();
02051         while (__u == 0.0);
02052         
02053         return (std::pow(__u, result_type(1.0) / __param.alpha())
02054             * __a1 * __v * __param.beta());
02055       }
02056       }
02057 
02058   template<typename _RealType, typename _CharT, typename _Traits>
02059     std::basic_ostream<_CharT, _Traits>&
02060     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02061            const gamma_distribution<_RealType>& __x)
02062     {
02063       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02064       typedef typename __ostream_type::ios_base    __ios_base;
02065 
02066       const typename __ios_base::fmtflags __flags = __os.flags();
02067       const _CharT __fill = __os.fill();
02068       const std::streamsize __precision = __os.precision();
02069       const _CharT __space = __os.widen(' ');
02070       __os.flags(__ios_base::scientific | __ios_base::left);
02071       __os.fill(__space);
02072       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02073 
02074       __os << __x.alpha() << __space << __x.beta()
02075        << __space << __x._M_nd;
02076 
02077       __os.flags(__flags);
02078       __os.fill(__fill);
02079       __os.precision(__precision);
02080       return __os;
02081     }
02082 
02083   template<typename _RealType, typename _CharT, typename _Traits>
02084     std::basic_istream<_CharT, _Traits>&
02085     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02086            gamma_distribution<_RealType>& __x)
02087     {
02088       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02089       typedef typename __istream_type::ios_base    __ios_base;
02090 
02091       const typename __ios_base::fmtflags __flags = __is.flags();
02092       __is.flags(__ios_base::dec | __ios_base::skipws);
02093 
02094       _RealType __alpha_val, __beta_val;
02095       __is >> __alpha_val >> __beta_val >> __x._M_nd;
02096       __x.param(typename gamma_distribution<_RealType>::
02097         param_type(__alpha_val, __beta_val));
02098 
02099       __is.flags(__flags);
02100       return __is;
02101     }
02102 
02103 
02104   template<typename _RealType>
02105     template<typename _UniformRandomNumberGenerator>
02106       typename weibull_distribution<_RealType>::result_type
02107       weibull_distribution<_RealType>::
02108       operator()(_UniformRandomNumberGenerator& __urng,
02109          const param_type& __p)
02110       {
02111     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02112       __aurng(__urng);
02113     return __p.b() * std::pow(-std::log(__aurng()),
02114                   result_type(1) / __p.a());
02115       }
02116 
02117   template<typename _RealType, typename _CharT, typename _Traits>
02118     std::basic_ostream<_CharT, _Traits>&
02119     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02120            const weibull_distribution<_RealType>& __x)
02121     {
02122       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02123       typedef typename __ostream_type::ios_base    __ios_base;
02124 
02125       const typename __ios_base::fmtflags __flags = __os.flags();
02126       const _CharT __fill = __os.fill();
02127       const std::streamsize __precision = __os.precision();
02128       const _CharT __space = __os.widen(' ');
02129       __os.flags(__ios_base::scientific | __ios_base::left);
02130       __os.fill(__space);
02131       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02132 
02133       __os << __x.a() << __space << __x.b();
02134 
02135       __os.flags(__flags);
02136       __os.fill(__fill);
02137       __os.precision(__precision);
02138       return __os;
02139     }
02140 
02141   template<typename _RealType, typename _CharT, typename _Traits>
02142     std::basic_istream<_CharT, _Traits>&
02143     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02144            weibull_distribution<_RealType>& __x)
02145     {
02146       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02147       typedef typename __istream_type::ios_base    __ios_base;
02148 
02149       const typename __ios_base::fmtflags __flags = __is.flags();
02150       __is.flags(__ios_base::dec | __ios_base::skipws);
02151 
02152       _RealType __a, __b;
02153       __is >> __a >> __b;
02154       __x.param(typename weibull_distribution<_RealType>::
02155         param_type(__a, __b));
02156 
02157       __is.flags(__flags);
02158       return __is;
02159     }
02160 
02161 
02162   template<typename _RealType>
02163     template<typename _UniformRandomNumberGenerator>
02164       typename extreme_value_distribution<_RealType>::result_type
02165       extreme_value_distribution<_RealType>::
02166       operator()(_UniformRandomNumberGenerator& __urng,
02167          const param_type& __p)
02168       {
02169     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02170       __aurng(__urng);
02171     return __p.a() - __p.b() * std::log(-std::log(__aurng()));
02172       }
02173 
02174   template<typename _RealType, typename _CharT, typename _Traits>
02175     std::basic_ostream<_CharT, _Traits>&
02176     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02177            const extreme_value_distribution<_RealType>& __x)
02178     {
02179       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02180       typedef typename __ostream_type::ios_base    __ios_base;
02181 
02182       const typename __ios_base::fmtflags __flags = __os.flags();
02183       const _CharT __fill = __os.fill();
02184       const std::streamsize __precision = __os.precision();
02185       const _CharT __space = __os.widen(' ');
02186       __os.flags(__ios_base::scientific | __ios_base::left);
02187       __os.fill(__space);
02188       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02189 
02190       __os << __x.a() << __space << __x.b();
02191 
02192       __os.flags(__flags);
02193       __os.fill(__fill);
02194       __os.precision(__precision);
02195       return __os;
02196     }
02197 
02198   template<typename _RealType, typename _CharT, typename _Traits>
02199     std::basic_istream<_CharT, _Traits>&
02200     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02201            extreme_value_distribution<_RealType>& __x)
02202     {
02203       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02204       typedef typename __istream_type::ios_base    __ios_base;
02205 
02206       const typename __ios_base::fmtflags __flags = __is.flags();
02207       __is.flags(__ios_base::dec | __ios_base::skipws);
02208 
02209       _RealType __a, __b;
02210       __is >> __a >> __b;
02211       __x.param(typename extreme_value_distribution<_RealType>::
02212         param_type(__a, __b));
02213 
02214       __is.flags(__flags);
02215       return __is;
02216     }
02217 
02218 
02219   template<typename _IntType>
02220     void
02221     discrete_distribution<_IntType>::param_type::
02222     _M_initialize()
02223     {
02224       if (_M_prob.size() < 2)
02225     {
02226       _M_prob.clear();
02227       return;
02228     }
02229 
02230       const double __sum = std::accumulate(_M_prob.begin(),
02231                        _M_prob.end(), 0.0);
02232       // Now normalize the probabilites.
02233       __detail::__transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
02234               std::bind2nd(std::divides<double>(), __sum));
02235       // Accumulate partial sums.
02236       _M_cp.reserve(_M_prob.size());
02237       std::partial_sum(_M_prob.begin(), _M_prob.end(),
02238                std::back_inserter(_M_cp));
02239       // Make sure the last cumulative probability is one.
02240       _M_cp[_M_cp.size() - 1] = 1.0;
02241     }
02242 
02243   template<typename _IntType>
02244     template<typename _Func>
02245       discrete_distribution<_IntType>::param_type::
02246       param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
02247       : _M_prob(), _M_cp()
02248       {
02249     const size_t __n = __nw == 0 ? 1 : __nw;
02250     const double __delta = (__xmax - __xmin) / __n;
02251 
02252     _M_prob.reserve(__n);
02253     for (size_t __k = 0; __k < __nw; ++__k)
02254       _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
02255 
02256     _M_initialize();
02257       }
02258 
02259   template<typename _IntType>
02260     template<typename _UniformRandomNumberGenerator>
02261       typename discrete_distribution<_IntType>::result_type
02262       discrete_distribution<_IntType>::
02263       operator()(_UniformRandomNumberGenerator& __urng,
02264          const param_type& __param)
02265       {
02266     if (__param._M_cp.empty())
02267       return result_type(0);
02268 
02269     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02270       __aurng(__urng);
02271 
02272     const double __p = __aurng();
02273     auto __pos = std::lower_bound(__param._M_cp.begin(),
02274                       __param._M_cp.end(), __p);
02275 
02276     return __pos - __param._M_cp.begin();
02277       }
02278 
02279   template<typename _IntType, typename _CharT, typename _Traits>
02280     std::basic_ostream<_CharT, _Traits>&
02281     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02282            const discrete_distribution<_IntType>& __x)
02283     {
02284       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02285       typedef typename __ostream_type::ios_base    __ios_base;
02286 
02287       const typename __ios_base::fmtflags __flags = __os.flags();
02288       const _CharT __fill = __os.fill();
02289       const std::streamsize __precision = __os.precision();
02290       const _CharT __space = __os.widen(' ');
02291       __os.flags(__ios_base::scientific | __ios_base::left);
02292       __os.fill(__space);
02293       __os.precision(std::numeric_limits<double>::max_digits10);
02294 
02295       std::vector<double> __prob = __x.probabilities();
02296       __os << __prob.size();
02297       for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
02298     __os << __space << *__dit;
02299 
02300       __os.flags(__flags);
02301       __os.fill(__fill);
02302       __os.precision(__precision);
02303       return __os;
02304     }
02305 
02306   template<typename _IntType, typename _CharT, typename _Traits>
02307     std::basic_istream<_CharT, _Traits>&
02308     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02309            discrete_distribution<_IntType>& __x)
02310     {
02311       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02312       typedef typename __istream_type::ios_base    __ios_base;
02313 
02314       const typename __ios_base::fmtflags __flags = __is.flags();
02315       __is.flags(__ios_base::dec | __ios_base::skipws);
02316 
02317       size_t __n;
02318       __is >> __n;
02319 
02320       std::vector<double> __prob_vec;
02321       __prob_vec.reserve(__n);
02322       for (; __n != 0; --__n)
02323     {
02324       double __prob;
02325       __is >> __prob;
02326       __prob_vec.push_back(__prob);
02327     }
02328 
02329       __x.param(typename discrete_distribution<_IntType>::
02330         param_type(__prob_vec.begin(), __prob_vec.end()));
02331 
02332       __is.flags(__flags);
02333       return __is;
02334     }
02335 
02336 
02337   template<typename _RealType>
02338     void
02339     piecewise_constant_distribution<_RealType>::param_type::
02340     _M_initialize()
02341     {
02342       if (_M_int.size() < 2
02343       || (_M_int.size() == 2
02344           && _M_int[0] == _RealType(0)
02345           && _M_int[1] == _RealType(1)))
02346     {
02347       _M_int.clear();
02348       _M_den.clear();
02349       return;
02350     }
02351 
02352       const double __sum = std::accumulate(_M_den.begin(),
02353                        _M_den.end(), 0.0);
02354 
02355       __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
02356                 std::bind2nd(std::divides<double>(), __sum));
02357 
02358       _M_cp.reserve(_M_den.size());
02359       std::partial_sum(_M_den.begin(), _M_den.end(),
02360                std::back_inserter(_M_cp));
02361 
02362       // Make sure the last cumulative probability is one.
02363       _M_cp[_M_cp.size() - 1] = 1.0;
02364 
02365       for (size_t __k = 0; __k < _M_den.size(); ++__k)
02366     _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
02367     }
02368 
02369   template<typename _RealType>
02370     template<typename _InputIteratorB, typename _InputIteratorW>
02371       piecewise_constant_distribution<_RealType>::param_type::
02372       param_type(_InputIteratorB __bbegin,
02373          _InputIteratorB __bend,
02374          _InputIteratorW __wbegin)
02375       : _M_int(), _M_den(), _M_cp()
02376       {
02377     if (__bbegin != __bend)
02378       {
02379         for (;;)
02380           {
02381         _M_int.push_back(*__bbegin);
02382         ++__bbegin;
02383         if (__bbegin == __bend)
02384           break;
02385 
02386         _M_den.push_back(*__wbegin);
02387         ++__wbegin;
02388           }
02389       }
02390 
02391     _M_initialize();
02392       }
02393 
02394   template<typename _RealType>
02395     template<typename _Func>
02396       piecewise_constant_distribution<_RealType>::param_type::
02397       param_type(initializer_list<_RealType> __bl, _Func __fw)
02398       : _M_int(), _M_den(), _M_cp()
02399       {
02400     _M_int.reserve(__bl.size());
02401     for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
02402       _M_int.push_back(*__biter);
02403 
02404     _M_den.reserve(_M_int.size() - 1);
02405     for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
02406       _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
02407 
02408     _M_initialize();
02409       }
02410 
02411   template<typename _RealType>
02412     template<typename _Func>
02413       piecewise_constant_distribution<_RealType>::param_type::
02414       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
02415       : _M_int(), _M_den(), _M_cp()
02416       {
02417     const size_t __n = __nw == 0 ? 1 : __nw;
02418     const _RealType __delta = (__xmax - __xmin) / __n;
02419 
02420     _M_int.reserve(__n + 1);
02421     for (size_t __k = 0; __k <= __nw; ++__k)
02422       _M_int.push_back(__xmin + __k * __delta);
02423 
02424     _M_den.reserve(__n);
02425     for (size_t __k = 0; __k < __nw; ++__k)
02426       _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
02427 
02428     _M_initialize();
02429       }
02430 
02431   template<typename _RealType>
02432     template<typename _UniformRandomNumberGenerator>
02433       typename piecewise_constant_distribution<_RealType>::result_type
02434       piecewise_constant_distribution<_RealType>::
02435       operator()(_UniformRandomNumberGenerator& __urng,
02436          const param_type& __param)
02437       {
02438     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02439       __aurng(__urng);
02440 
02441     const double __p = __aurng();
02442     if (__param._M_cp.empty())
02443       return __p;
02444 
02445     auto __pos = std::lower_bound(__param._M_cp.begin(),
02446                       __param._M_cp.end(), __p);
02447     const size_t __i = __pos - __param._M_cp.begin();
02448 
02449     const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
02450 
02451     return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
02452       }
02453 
02454   template<typename _RealType, typename _CharT, typename _Traits>
02455     std::basic_ostream<_CharT, _Traits>&
02456     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02457            const piecewise_constant_distribution<_RealType>& __x)
02458     {
02459       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02460       typedef typename __ostream_type::ios_base    __ios_base;
02461 
02462       const typename __ios_base::fmtflags __flags = __os.flags();
02463       const _CharT __fill = __os.fill();
02464       const std::streamsize __precision = __os.precision();
02465       const _CharT __space = __os.widen(' ');
02466       __os.flags(__ios_base::scientific | __ios_base::left);
02467       __os.fill(__space);
02468       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02469 
02470       std::vector<_RealType> __int = __x.intervals();
02471       __os << __int.size() - 1;
02472 
02473       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
02474     __os << __space << *__xit;
02475 
02476       std::vector<double> __den = __x.densities();
02477       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
02478     __os << __space << *__dit;
02479 
02480       __os.flags(__flags);
02481       __os.fill(__fill);
02482       __os.precision(__precision);
02483       return __os;
02484     }
02485 
02486   template<typename _RealType, typename _CharT, typename _Traits>
02487     std::basic_istream<_CharT, _Traits>&
02488     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02489            piecewise_constant_distribution<_RealType>& __x)
02490     {
02491       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02492       typedef typename __istream_type::ios_base    __ios_base;
02493 
02494       const typename __ios_base::fmtflags __flags = __is.flags();
02495       __is.flags(__ios_base::dec | __ios_base::skipws);
02496 
02497       size_t __n;
02498       __is >> __n;
02499 
02500       std::vector<_RealType> __int_vec;
02501       __int_vec.reserve(__n + 1);
02502       for (size_t __i = 0; __i <= __n; ++__i)
02503     {
02504       _RealType __int;
02505       __is >> __int;
02506       __int_vec.push_back(__int);
02507     }
02508 
02509       std::vector<double> __den_vec;
02510       __den_vec.reserve(__n);
02511       for (size_t __i = 0; __i < __n; ++__i)
02512     {
02513       double __den;
02514       __is >> __den;
02515       __den_vec.push_back(__den);
02516     }
02517 
02518       __x.param(typename piecewise_constant_distribution<_RealType>::
02519       param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
02520 
02521       __is.flags(__flags);
02522       return __is;
02523     }
02524 
02525 
02526   template<typename _RealType>
02527     void
02528     piecewise_linear_distribution<_RealType>::param_type::
02529     _M_initialize()
02530     {
02531       if (_M_int.size() < 2
02532       || (_M_int.size() == 2
02533           && _M_int[0] == _RealType(0)
02534           && _M_int[1] == _RealType(1)
02535           && _M_den[0] == _M_den[1]))
02536     {
02537       _M_int.clear();
02538       _M_den.clear();
02539       return;
02540     }
02541 
02542       double __sum = 0.0;
02543       _M_cp.reserve(_M_int.size() - 1);
02544       _M_m.reserve(_M_int.size() - 1);
02545       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
02546     {
02547       const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
02548       __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
02549       _M_cp.push_back(__sum);
02550       _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
02551     }
02552 
02553       //  Now normalize the densities...
02554       __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
02555               std::bind2nd(std::divides<double>(), __sum));
02556       //  ... and partial sums... 
02557       __detail::__transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
02558                 std::bind2nd(std::divides<double>(), __sum));
02559       //  ... and slopes.
02560       __detail::__transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
02561                 std::bind2nd(std::divides<double>(), __sum));
02562       //  Make sure the last cumulative probablility is one.
02563       _M_cp[_M_cp.size() - 1] = 1.0;
02564      }
02565 
02566   template<typename _RealType>
02567     template<typename _InputIteratorB, typename _InputIteratorW>
02568       piecewise_linear_distribution<_RealType>::param_type::
02569       param_type(_InputIteratorB __bbegin,
02570          _InputIteratorB __bend,
02571          _InputIteratorW __wbegin)
02572       : _M_int(), _M_den(), _M_cp(), _M_m()
02573       {
02574     for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
02575       {
02576         _M_int.push_back(*__bbegin);
02577         _M_den.push_back(*__wbegin);
02578       }
02579 
02580     _M_initialize();
02581       }
02582 
02583   template<typename _RealType>
02584     template<typename _Func>
02585       piecewise_linear_distribution<_RealType>::param_type::
02586       param_type(initializer_list<_RealType> __bl, _Func __fw)
02587       : _M_int(), _M_den(), _M_cp(), _M_m()
02588       {
02589     _M_int.reserve(__bl.size());
02590     _M_den.reserve(__bl.size());
02591     for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
02592       {
02593         _M_int.push_back(*__biter);
02594         _M_den.push_back(__fw(*__biter));
02595       }
02596 
02597     _M_initialize();
02598       }
02599 
02600   template<typename _RealType>
02601     template<typename _Func>
02602       piecewise_linear_distribution<_RealType>::param_type::
02603       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
02604       : _M_int(), _M_den(), _M_cp(), _M_m()
02605       {
02606     const size_t __n = __nw == 0 ? 1 : __nw;
02607     const _RealType __delta = (__xmax - __xmin) / __n;
02608 
02609     _M_int.reserve(__n + 1);
02610     _M_den.reserve(__n + 1);
02611     for (size_t __k = 0; __k <= __nw; ++__k)
02612       {
02613         _M_int.push_back(__xmin + __k * __delta);
02614         _M_den.push_back(__fw(_M_int[__k] + __delta));
02615       }
02616 
02617     _M_initialize();
02618       }
02619 
02620   template<typename _RealType>
02621     template<typename _UniformRandomNumberGenerator>
02622       typename piecewise_linear_distribution<_RealType>::result_type
02623       piecewise_linear_distribution<_RealType>::
02624       operator()(_UniformRandomNumberGenerator& __urng,
02625          const param_type& __param)
02626       {
02627     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02628       __aurng(__urng);
02629 
02630     const double __p = __aurng();
02631     if (__param._M_cp.empty())
02632       return __p;
02633 
02634     auto __pos = std::lower_bound(__param._M_cp.begin(),
02635                       __param._M_cp.end(), __p);
02636     const size_t __i = __pos - __param._M_cp.begin();
02637 
02638     const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
02639 
02640     const double __a = 0.5 * __param._M_m[__i];
02641     const double __b = __param._M_den[__i];
02642     const double __cm = __p - __pref;
02643 
02644     _RealType __x = __param._M_int[__i];
02645     if (__a == 0)
02646       __x += __cm / __b;
02647     else
02648       {
02649         const double __d = __b * __b + 4.0 * __a * __cm;
02650         __x += 0.5 * (std::sqrt(__d) - __b) / __a;
02651           }
02652 
02653         return __x;
02654       }
02655 
02656   template<typename _RealType, typename _CharT, typename _Traits>
02657     std::basic_ostream<_CharT, _Traits>&
02658     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02659            const piecewise_linear_distribution<_RealType>& __x)
02660     {
02661       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02662       typedef typename __ostream_type::ios_base    __ios_base;
02663 
02664       const typename __ios_base::fmtflags __flags = __os.flags();
02665       const _CharT __fill = __os.fill();
02666       const std::streamsize __precision = __os.precision();
02667       const _CharT __space = __os.widen(' ');
02668       __os.flags(__ios_base::scientific | __ios_base::left);
02669       __os.fill(__space);
02670       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02671 
02672       std::vector<_RealType> __int = __x.intervals();
02673       __os << __int.size() - 1;
02674 
02675       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
02676     __os << __space << *__xit;
02677 
02678       std::vector<double> __den = __x.densities();
02679       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
02680     __os << __space << *__dit;
02681 
02682       __os.flags(__flags);
02683       __os.fill(__fill);
02684       __os.precision(__precision);
02685       return __os;
02686     }
02687 
02688   template<typename _RealType, typename _CharT, typename _Traits>
02689     std::basic_istream<_CharT, _Traits>&
02690     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02691            piecewise_linear_distribution<_RealType>& __x)
02692     {
02693       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02694       typedef typename __istream_type::ios_base    __ios_base;
02695 
02696       const typename __ios_base::fmtflags __flags = __is.flags();
02697       __is.flags(__ios_base::dec | __ios_base::skipws);
02698 
02699       size_t __n;
02700       __is >> __n;
02701 
02702       std::vector<_RealType> __int_vec;
02703       __int_vec.reserve(__n + 1);
02704       for (size_t __i = 0; __i <= __n; ++__i)
02705     {
02706       _RealType __int;
02707       __is >> __int;
02708       __int_vec.push_back(__int);
02709     }
02710 
02711       std::vector<double> __den_vec;
02712       __den_vec.reserve(__n + 1);
02713       for (size_t __i = 0; __i <= __n; ++__i)
02714     {
02715       double __den;
02716       __is >> __den;
02717       __den_vec.push_back(__den);
02718     }
02719 
02720       __x.param(typename piecewise_linear_distribution<_RealType>::
02721       param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
02722 
02723       __is.flags(__flags);
02724       return __is;
02725     }
02726 
02727 
02728   template<typename _IntType>
02729     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
02730     {
02731       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
02732     _M_v.push_back(__detail::__mod<result_type,
02733                __detail::_Shift<result_type, 32>::__value>(*__iter));
02734     }
02735 
02736   template<typename _InputIterator>
02737     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
02738     {
02739       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
02740     _M_v.push_back(__detail::__mod<result_type,
02741                __detail::_Shift<result_type, 32>::__value>(*__iter));
02742     }
02743 
02744   template<typename _RandomAccessIterator>
02745     void
02746     seed_seq::generate(_RandomAccessIterator __begin,
02747                _RandomAccessIterator __end)
02748     {
02749       typedef typename iterator_traits<_RandomAccessIterator>::value_type
02750         _Type;
02751 
02752       if (__begin == __end)
02753     return;
02754 
02755       std::fill(__begin, __end, _Type(0x8b8b8b8bu));
02756 
02757       const size_t __n = __end - __begin;
02758       const size_t __s = _M_v.size();
02759       const size_t __t = (__n >= 623) ? 11
02760                : (__n >=  68) ? 7
02761                : (__n >=  39) ? 5
02762                : (__n >=   7) ? 3
02763                : (__n - 1) / 2;
02764       const size_t __p = (__n - __t) / 2;
02765       const size_t __q = __p + __t;
02766       const size_t __m = std::max(__s + 1, __n);
02767 
02768       for (size_t __k = 0; __k < __m; ++__k)
02769     {
02770       _Type __arg = (__begin[__k % __n]
02771              ^ __begin[(__k + __p) % __n]
02772              ^ __begin[(__k - 1) % __n]);
02773       _Type __r1 = __arg ^ (__arg >> 27);
02774       __r1 = __detail::__mod<_Type,
02775             __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
02776       _Type __r2 = __r1;
02777       if (__k == 0)
02778         __r2 += __s;
02779       else if (__k <= __s)
02780         __r2 += __k % __n + _M_v[__k - 1];
02781       else
02782         __r2 += __k % __n;
02783       __r2 = __detail::__mod<_Type,
02784                __detail::_Shift<_Type, 32>::__value>(__r2);
02785       __begin[(__k + __p) % __n] += __r1;
02786       __begin[(__k + __q) % __n] += __r2;
02787       __begin[__k % __n] = __r2;
02788     }
02789 
02790       for (size_t __k = __m; __k < __m + __n; ++__k)
02791     {
02792       _Type __arg = (__begin[__k % __n]
02793              + __begin[(__k + __p) % __n]
02794              + __begin[(__k - 1) % __n]);
02795       _Type __r3 = __arg ^ (__arg >> 27);
02796       __r3 = __detail::__mod<_Type,
02797            __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
02798       _Type __r4 = __r3 - __k % __n;
02799       __r4 = __detail::__mod<_Type,
02800                __detail::_Shift<_Type, 32>::__value>(__r4);
02801       __begin[(__k + __p) % __n] ^= __r3;
02802       __begin[(__k + __q) % __n] ^= __r4;
02803       __begin[__k % __n] = __r4;
02804     }
02805     }
02806 
02807   template<typename _RealType, size_t __bits,
02808        typename _UniformRandomNumberGenerator>
02809     _RealType
02810     generate_canonical(_UniformRandomNumberGenerator& __urng)
02811     {
02812       const size_t __b
02813     = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
02814                    __bits);
02815       const long double __r = static_cast<long double>(__urng.max())
02816                 - static_cast<long double>(__urng.min()) + 1.0L;
02817       const size_t __log2r = std::log(__r) / std::log(2.0L);
02818       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
02819       _RealType __sum = _RealType(0);
02820       _RealType __tmp = _RealType(1);
02821       for (; __k != 0; --__k)
02822     {
02823       __sum += _RealType(__urng() - __urng.min()) * __tmp;
02824       __tmp *= __r;
02825     }
02826       return __sum / __tmp;
02827     }
02828 
02829 _GLIBCXX_END_NAMESPACE_VERSION
02830 } // namespace
02831 
02832 #endif