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