random.tcc

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