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