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,
01652          _IntType __t, double __q)
01653       {
01654     _IntType __x = 0;
01655     double __sum = 0.0;
01656     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01657       __aurng(__urng);
01658 
01659     do
01660       {
01661         if (__t == __x)
01662           return __x;
01663         const double __e = -std::log(1.0 - __aurng());
01664         __sum += __e / (__t - __x);
01665         __x += 1;
01666       }
01667     while (__sum <= __q);
01668 
01669     return __x - 1;
01670       }
01671 
01672   /**
01673    * A rejection algorithm when t * p >= 8 and a simple waiting time
01674    * method - the second in the referenced book - otherwise.
01675    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
01676    * is defined.
01677    *
01678    * Reference:
01679    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01680    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
01681    */
01682   template<typename _IntType>
01683     template<typename _UniformRandomNumberGenerator>
01684       typename binomial_distribution<_IntType>::result_type
01685       binomial_distribution<_IntType>::
01686       operator()(_UniformRandomNumberGenerator& __urng,
01687          const param_type& __param)
01688       {
01689     result_type __ret;
01690     const _IntType __t = __param.t();
01691     const double __p = __param.p();
01692     const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
01693     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01694       __aurng(__urng);
01695 
01696 #if _GLIBCXX_USE_C99_MATH_TR1
01697     if (!__param._M_easy)
01698       {
01699         double __x;
01700 
01701         // See comments above...
01702         const double __naf =
01703           (1 - std::numeric_limits<double>::epsilon()) / 2;
01704         const double __thr =
01705           std::numeric_limits<_IntType>::max() + __naf;
01706 
01707         const double __np = std::floor(__t * __p12);
01708 
01709         // sqrt(pi / 2)
01710         const double __spi_2 = 1.2533141373155002512078826424055226L;
01711         const double __a1 = __param._M_a1;
01712         const double __a12 = __a1 + __param._M_s2 * __spi_2;
01713         const double __a123 = __param._M_a123;
01714         const double __s1s = __param._M_s1 * __param._M_s1;
01715         const double __s2s = __param._M_s2 * __param._M_s2;
01716 
01717         bool __reject;
01718         do
01719           {
01720         const double __u = __param._M_s * __aurng();
01721 
01722         double __v;
01723 
01724         if (__u <= __a1)
01725           {
01726             const double __n = _M_nd(__urng);
01727             const double __y = __param._M_s1 * std::abs(__n);
01728             __reject = __y >= __param._M_d1;
01729             if (!__reject)
01730               {
01731             const double __e = -std::log(1.0 - __aurng());
01732             __x = std::floor(__y);
01733             __v = -__e - __n * __n / 2 + __param._M_c;
01734               }
01735           }
01736         else if (__u <= __a12)
01737           {
01738             const double __n = _M_nd(__urng);
01739             const double __y = __param._M_s2 * std::abs(__n);
01740             __reject = __y >= __param._M_d2;
01741             if (!__reject)
01742               {
01743             const double __e = -std::log(1.0 - __aurng());
01744             __x = std::floor(-__y);
01745             __v = -__e - __n * __n / 2;
01746               }
01747           }
01748         else if (__u <= __a123)
01749           {
01750             const double __e1 = -std::log(1.0 - __aurng());
01751             const double __e2 = -std::log(1.0 - __aurng());
01752 
01753             const double __y = __param._M_d1
01754                      + 2 * __s1s * __e1 / __param._M_d1;
01755             __x = std::floor(__y);
01756             __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
01757                             -__y / (2 * __s1s)));
01758             __reject = false;
01759           }
01760         else
01761           {
01762             const double __e1 = -std::log(1.0 - __aurng());
01763             const double __e2 = -std::log(1.0 - __aurng());
01764 
01765             const double __y = __param._M_d2
01766                      + 2 * __s2s * __e1 / __param._M_d2;
01767             __x = std::floor(-__y);
01768             __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
01769             __reject = false;
01770           }
01771 
01772         __reject = __reject || __x < -__np || __x > __t - __np;
01773         if (!__reject)
01774           {
01775             const double __lfx =
01776               std::lgamma(__np + __x + 1)
01777               + std::lgamma(__t - (__np + __x) + 1);
01778             __reject = __v > __param._M_lf - __lfx
01779                  + __x * __param._M_lp1p;
01780           }
01781 
01782         __reject |= __x + __np >= __thr;
01783           }
01784         while (__reject);
01785 
01786         __x += __np + __naf;
01787 
01788         const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
01789                         __param._M_q);
01790         __ret = _IntType(__x) + __z;
01791       }
01792     else
01793 #endif
01794       __ret = _M_waiting(__urng, __t, __param._M_q);
01795 
01796     if (__p12 != __p)
01797       __ret = __t - __ret;
01798     return __ret;
01799       }
01800 
01801   template<typename _IntType>
01802     template<typename _ForwardIterator,
01803          typename _UniformRandomNumberGenerator>
01804       void
01805       binomial_distribution<_IntType>::
01806       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01807               _UniformRandomNumberGenerator& __urng,
01808               const param_type& __param)
01809       {
01810     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01811     // We could duplicate everything from operator()...
01812     while (__f != __t)
01813       *__f++ = this->operator()(__urng, __param);
01814       }
01815 
01816   template<typename _IntType,
01817        typename _CharT, typename _Traits>
01818     std::basic_ostream<_CharT, _Traits>&
01819     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01820            const binomial_distribution<_IntType>& __x)
01821     {
01822       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01823       typedef typename __ostream_type::ios_base    __ios_base;
01824 
01825       const typename __ios_base::fmtflags __flags = __os.flags();
01826       const _CharT __fill = __os.fill();
01827       const std::streamsize __precision = __os.precision();
01828       const _CharT __space = __os.widen(' ');
01829       __os.flags(__ios_base::scientific | __ios_base::left);
01830       __os.fill(__space);
01831       __os.precision(std::numeric_limits<double>::max_digits10);
01832 
01833       __os << __x.t() << __space << __x.p()
01834        << __space << __x._M_nd;
01835 
01836       __os.flags(__flags);
01837       __os.fill(__fill);
01838       __os.precision(__precision);
01839       return __os;
01840     }
01841 
01842   template<typename _IntType,
01843        typename _CharT, typename _Traits>
01844     std::basic_istream<_CharT, _Traits>&
01845     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01846            binomial_distribution<_IntType>& __x)
01847     {
01848       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01849       typedef typename __istream_type::ios_base    __ios_base;
01850 
01851       const typename __ios_base::fmtflags __flags = __is.flags();
01852       __is.flags(__ios_base::dec | __ios_base::skipws);
01853 
01854       _IntType __t;
01855       double __p;
01856       __is >> __t >> __p >> __x._M_nd;
01857       __x.param(typename binomial_distribution<_IntType>::
01858         param_type(__t, __p));
01859 
01860       __is.flags(__flags);
01861       return __is;
01862     }
01863 
01864 
01865   template<typename _RealType>
01866     template<typename _ForwardIterator,
01867          typename _UniformRandomNumberGenerator>
01868       void
01869       std::exponential_distribution<_RealType>::
01870       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01871               _UniformRandomNumberGenerator& __urng,
01872               const param_type& __p)
01873       {
01874     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01875     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01876       __aurng(__urng);
01877     while (__f != __t)
01878       *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
01879       }
01880 
01881   template<typename _RealType, typename _CharT, typename _Traits>
01882     std::basic_ostream<_CharT, _Traits>&
01883     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01884            const exponential_distribution<_RealType>& __x)
01885     {
01886       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01887       typedef typename __ostream_type::ios_base    __ios_base;
01888 
01889       const typename __ios_base::fmtflags __flags = __os.flags();
01890       const _CharT __fill = __os.fill();
01891       const std::streamsize __precision = __os.precision();
01892       __os.flags(__ios_base::scientific | __ios_base::left);
01893       __os.fill(__os.widen(' '));
01894       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01895 
01896       __os << __x.lambda();
01897 
01898       __os.flags(__flags);
01899       __os.fill(__fill);
01900       __os.precision(__precision);
01901       return __os;
01902     }
01903 
01904   template<typename _RealType, typename _CharT, typename _Traits>
01905     std::basic_istream<_CharT, _Traits>&
01906     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01907            exponential_distribution<_RealType>& __x)
01908     {
01909       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01910       typedef typename __istream_type::ios_base    __ios_base;
01911 
01912       const typename __ios_base::fmtflags __flags = __is.flags();
01913       __is.flags(__ios_base::dec | __ios_base::skipws);
01914 
01915       _RealType __lambda;
01916       __is >> __lambda;
01917       __x.param(typename exponential_distribution<_RealType>::
01918         param_type(__lambda));
01919 
01920       __is.flags(__flags);
01921       return __is;
01922     }
01923 
01924 
01925   /**
01926    * Polar method due to Marsaglia.
01927    *
01928    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01929    * New York, 1986, Ch. V, Sect. 4.4.
01930    */
01931   template<typename _RealType>
01932     template<typename _UniformRandomNumberGenerator>
01933       typename normal_distribution<_RealType>::result_type
01934       normal_distribution<_RealType>::
01935       operator()(_UniformRandomNumberGenerator& __urng,
01936          const param_type& __param)
01937       {
01938     result_type __ret;
01939     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01940       __aurng(__urng);
01941 
01942     if (_M_saved_available)
01943       {
01944         _M_saved_available = false;
01945         __ret = _M_saved;
01946       }
01947     else
01948       {
01949         result_type __x, __y, __r2;
01950         do
01951           {
01952         __x = result_type(2.0) * __aurng() - 1.0;
01953         __y = result_type(2.0) * __aurng() - 1.0;
01954         __r2 = __x * __x + __y * __y;
01955           }
01956         while (__r2 > 1.0 || __r2 == 0.0);
01957 
01958         const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
01959         _M_saved = __x * __mult;
01960         _M_saved_available = true;
01961         __ret = __y * __mult;
01962       }
01963 
01964     __ret = __ret * __param.stddev() + __param.mean();
01965     return __ret;
01966       }
01967 
01968   template<typename _RealType>
01969     template<typename _ForwardIterator,
01970          typename _UniformRandomNumberGenerator>
01971       void
01972       normal_distribution<_RealType>::
01973       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01974               _UniformRandomNumberGenerator& __urng,
01975               const param_type& __param)
01976       {
01977     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01978 
01979     if (__f == __t)
01980       return;
01981 
01982     if (_M_saved_available)
01983       {
01984         _M_saved_available = false;
01985         *__f++ = _M_saved * __param.stddev() + __param.mean();
01986 
01987         if (__f == __t)
01988           return;
01989       }
01990 
01991     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01992       __aurng(__urng);
01993 
01994     while (__f + 1 < __t)
01995       {
01996         result_type __x, __y, __r2;
01997         do
01998           {
01999         __x = result_type(2.0) * __aurng() - 1.0;
02000         __y = result_type(2.0) * __aurng() - 1.0;
02001         __r2 = __x * __x + __y * __y;
02002           }
02003         while (__r2 > 1.0 || __r2 == 0.0);
02004 
02005         const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
02006         *__f++ = __y * __mult * __param.stddev() + __param.mean();
02007         *__f++ = __x * __mult * __param.stddev() + __param.mean();
02008       }
02009 
02010     if (__f != __t)
02011       {
02012         result_type __x, __y, __r2;
02013         do
02014           {
02015         __x = result_type(2.0) * __aurng() - 1.0;
02016         __y = result_type(2.0) * __aurng() - 1.0;
02017         __r2 = __x * __x + __y * __y;
02018           }
02019         while (__r2 > 1.0 || __r2 == 0.0);
02020 
02021         const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
02022         _M_saved = __x * __mult;
02023         _M_saved_available = true;
02024         *__f = __y * __mult * __param.stddev() + __param.mean();
02025       }
02026       }
02027 
02028   template<typename _RealType>
02029     bool
02030     operator==(const std::normal_distribution<_RealType>& __d1,
02031            const std::normal_distribution<_RealType>& __d2)
02032     {
02033       if (__d1._M_param == __d2._M_param
02034       && __d1._M_saved_available == __d2._M_saved_available)
02035     {
02036       if (__d1._M_saved_available
02037           && __d1._M_saved == __d2._M_saved)
02038         return true;
02039       else if(!__d1._M_saved_available)
02040         return true;
02041       else
02042         return false;
02043     }
02044       else
02045     return false;
02046     }
02047 
02048   template<typename _RealType, typename _CharT, typename _Traits>
02049     std::basic_ostream<_CharT, _Traits>&
02050     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02051            const normal_distribution<_RealType>& __x)
02052     {
02053       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02054       typedef typename __ostream_type::ios_base    __ios_base;
02055 
02056       const typename __ios_base::fmtflags __flags = __os.flags();
02057       const _CharT __fill = __os.fill();
02058       const std::streamsize __precision = __os.precision();
02059       const _CharT __space = __os.widen(' ');
02060       __os.flags(__ios_base::scientific | __ios_base::left);
02061       __os.fill(__space);
02062       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02063 
02064       __os << __x.mean() << __space << __x.stddev()
02065        << __space << __x._M_saved_available;
02066       if (__x._M_saved_available)
02067     __os << __space << __x._M_saved;
02068 
02069       __os.flags(__flags);
02070       __os.fill(__fill);
02071       __os.precision(__precision);
02072       return __os;
02073     }
02074 
02075   template<typename _RealType, typename _CharT, typename _Traits>
02076     std::basic_istream<_CharT, _Traits>&
02077     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02078            normal_distribution<_RealType>& __x)
02079     {
02080       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02081       typedef typename __istream_type::ios_base    __ios_base;
02082 
02083       const typename __ios_base::fmtflags __flags = __is.flags();
02084       __is.flags(__ios_base::dec | __ios_base::skipws);
02085 
02086       double __mean, __stddev;
02087       __is >> __mean >> __stddev
02088        >> __x._M_saved_available;
02089       if (__x._M_saved_available)
02090     __is >> __x._M_saved;
02091       __x.param(typename normal_distribution<_RealType>::
02092         param_type(__mean, __stddev));
02093 
02094       __is.flags(__flags);
02095       return __is;
02096     }
02097 
02098 
02099   template<typename _RealType>
02100     template<typename _ForwardIterator,
02101          typename _UniformRandomNumberGenerator>
02102       void
02103       lognormal_distribution<_RealType>::
02104       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02105               _UniformRandomNumberGenerator& __urng,
02106               const param_type& __p)
02107       {
02108     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02109       while (__f != __t)
02110         *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
02111       }
02112 
02113   template<typename _RealType, typename _CharT, typename _Traits>
02114     std::basic_ostream<_CharT, _Traits>&
02115     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02116            const lognormal_distribution<_RealType>& __x)
02117     {
02118       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02119       typedef typename __ostream_type::ios_base    __ios_base;
02120 
02121       const typename __ios_base::fmtflags __flags = __os.flags();
02122       const _CharT __fill = __os.fill();
02123       const std::streamsize __precision = __os.precision();
02124       const _CharT __space = __os.widen(' ');
02125       __os.flags(__ios_base::scientific | __ios_base::left);
02126       __os.fill(__space);
02127       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02128 
02129       __os << __x.m() << __space << __x.s()
02130        << __space << __x._M_nd;
02131 
02132       __os.flags(__flags);
02133       __os.fill(__fill);
02134       __os.precision(__precision);
02135       return __os;
02136     }
02137 
02138   template<typename _RealType, typename _CharT, typename _Traits>
02139     std::basic_istream<_CharT, _Traits>&
02140     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02141            lognormal_distribution<_RealType>& __x)
02142     {
02143       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02144       typedef typename __istream_type::ios_base    __ios_base;
02145 
02146       const typename __ios_base::fmtflags __flags = __is.flags();
02147       __is.flags(__ios_base::dec | __ios_base::skipws);
02148 
02149       _RealType __m, __s;
02150       __is >> __m >> __s >> __x._M_nd;
02151       __x.param(typename lognormal_distribution<_RealType>::
02152         param_type(__m, __s));
02153 
02154       __is.flags(__flags);
02155       return __is;
02156     }
02157 
02158   template<typename _RealType>
02159     template<typename _ForwardIterator,
02160          typename _UniformRandomNumberGenerator>
02161       void
02162       std::chi_squared_distribution<_RealType>::
02163       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02164               _UniformRandomNumberGenerator& __urng)
02165       {
02166     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02167     while (__f != __t)
02168       *__f++ = 2 * _M_gd(__urng);
02169       }
02170 
02171   template<typename _RealType>
02172     template<typename _ForwardIterator,
02173          typename _UniformRandomNumberGenerator>
02174       void
02175       std::chi_squared_distribution<_RealType>::
02176       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02177               _UniformRandomNumberGenerator& __urng,
02178               const typename
02179               std::gamma_distribution<result_type>::param_type& __p)
02180       {
02181     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02182     while (__f != __t)
02183       *__f++ = 2 * _M_gd(__urng, __p);
02184       }
02185 
02186   template<typename _RealType, typename _CharT, typename _Traits>
02187     std::basic_ostream<_CharT, _Traits>&
02188     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02189            const chi_squared_distribution<_RealType>& __x)
02190     {
02191       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02192       typedef typename __ostream_type::ios_base    __ios_base;
02193 
02194       const typename __ios_base::fmtflags __flags = __os.flags();
02195       const _CharT __fill = __os.fill();
02196       const std::streamsize __precision = __os.precision();
02197       const _CharT __space = __os.widen(' ');
02198       __os.flags(__ios_base::scientific | __ios_base::left);
02199       __os.fill(__space);
02200       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02201 
02202       __os << __x.n() << __space << __x._M_gd;
02203 
02204       __os.flags(__flags);
02205       __os.fill(__fill);
02206       __os.precision(__precision);
02207       return __os;
02208     }
02209 
02210   template<typename _RealType, typename _CharT, typename _Traits>
02211     std::basic_istream<_CharT, _Traits>&
02212     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02213            chi_squared_distribution<_RealType>& __x)
02214     {
02215       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02216       typedef typename __istream_type::ios_base    __ios_base;
02217 
02218       const typename __ios_base::fmtflags __flags = __is.flags();
02219       __is.flags(__ios_base::dec | __ios_base::skipws);
02220 
02221       _RealType __n;
02222       __is >> __n >> __x._M_gd;
02223       __x.param(typename chi_squared_distribution<_RealType>::
02224         param_type(__n));
02225 
02226       __is.flags(__flags);
02227       return __is;
02228     }
02229 
02230 
02231   template<typename _RealType>
02232     template<typename _UniformRandomNumberGenerator>
02233       typename cauchy_distribution<_RealType>::result_type
02234       cauchy_distribution<_RealType>::
02235       operator()(_UniformRandomNumberGenerator& __urng,
02236          const param_type& __p)
02237       {
02238     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02239       __aurng(__urng);
02240     _RealType __u;
02241     do
02242       __u = __aurng();
02243     while (__u == 0.5);
02244 
02245     const _RealType __pi = 3.1415926535897932384626433832795029L;
02246     return __p.a() + __p.b() * std::tan(__pi * __u);
02247       }
02248 
02249   template<typename _RealType>
02250     template<typename _ForwardIterator,
02251          typename _UniformRandomNumberGenerator>
02252       void
02253       cauchy_distribution<_RealType>::
02254       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02255               _UniformRandomNumberGenerator& __urng,
02256               const param_type& __p)
02257       {
02258     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02259     const _RealType __pi = 3.1415926535897932384626433832795029L;
02260     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02261       __aurng(__urng);
02262     while (__f != __t)
02263       {
02264         _RealType __u;
02265         do
02266           __u = __aurng();
02267         while (__u == 0.5);
02268 
02269         *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
02270       }
02271       }
02272 
02273   template<typename _RealType, typename _CharT, typename _Traits>
02274     std::basic_ostream<_CharT, _Traits>&
02275     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02276            const cauchy_distribution<_RealType>& __x)
02277     {
02278       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02279       typedef typename __ostream_type::ios_base    __ios_base;
02280 
02281       const typename __ios_base::fmtflags __flags = __os.flags();
02282       const _CharT __fill = __os.fill();
02283       const std::streamsize __precision = __os.precision();
02284       const _CharT __space = __os.widen(' ');
02285       __os.flags(__ios_base::scientific | __ios_base::left);
02286       __os.fill(__space);
02287       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02288 
02289       __os << __x.a() << __space << __x.b();
02290 
02291       __os.flags(__flags);
02292       __os.fill(__fill);
02293       __os.precision(__precision);
02294       return __os;
02295     }
02296 
02297   template<typename _RealType, typename _CharT, typename _Traits>
02298     std::basic_istream<_CharT, _Traits>&
02299     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02300            cauchy_distribution<_RealType>& __x)
02301     {
02302       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02303       typedef typename __istream_type::ios_base    __ios_base;
02304 
02305       const typename __ios_base::fmtflags __flags = __is.flags();
02306       __is.flags(__ios_base::dec | __ios_base::skipws);
02307 
02308       _RealType __a, __b;
02309       __is >> __a >> __b;
02310       __x.param(typename cauchy_distribution<_RealType>::
02311         param_type(__a, __b));
02312 
02313       __is.flags(__flags);
02314       return __is;
02315     }
02316 
02317 
02318   template<typename _RealType>
02319     template<typename _ForwardIterator,
02320          typename _UniformRandomNumberGenerator>
02321       void
02322       std::fisher_f_distribution<_RealType>::
02323       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02324               _UniformRandomNumberGenerator& __urng)
02325       {
02326     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02327     while (__f != __t)
02328       *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
02329       }
02330 
02331   template<typename _RealType>
02332     template<typename _ForwardIterator,
02333          typename _UniformRandomNumberGenerator>
02334       void
02335       std::fisher_f_distribution<_RealType>::
02336       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02337               _UniformRandomNumberGenerator& __urng,
02338               const param_type& __p)
02339       {
02340     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02341     typedef typename std::gamma_distribution<result_type>::param_type
02342       param_type;
02343     param_type __p1(__p.m() / 2);
02344     param_type __p2(__p.n() / 2);
02345     while (__f != __t)
02346       *__f++ = ((_M_gd_x(__urng, __p1) * n())
02347             / (_M_gd_y(__urng, __p2) * m()));
02348       }
02349 
02350   template<typename _RealType, typename _CharT, typename _Traits>
02351     std::basic_ostream<_CharT, _Traits>&
02352     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02353            const fisher_f_distribution<_RealType>& __x)
02354     {
02355       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02356       typedef typename __ostream_type::ios_base    __ios_base;
02357 
02358       const typename __ios_base::fmtflags __flags = __os.flags();
02359       const _CharT __fill = __os.fill();
02360       const std::streamsize __precision = __os.precision();
02361       const _CharT __space = __os.widen(' ');
02362       __os.flags(__ios_base::scientific | __ios_base::left);
02363       __os.fill(__space);
02364       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02365 
02366       __os << __x.m() << __space << __x.n()
02367        << __space << __x._M_gd_x << __space << __x._M_gd_y;
02368 
02369       __os.flags(__flags);
02370       __os.fill(__fill);
02371       __os.precision(__precision);
02372       return __os;
02373     }
02374 
02375   template<typename _RealType, typename _CharT, typename _Traits>
02376     std::basic_istream<_CharT, _Traits>&
02377     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02378            fisher_f_distribution<_RealType>& __x)
02379     {
02380       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02381       typedef typename __istream_type::ios_base    __ios_base;
02382 
02383       const typename __ios_base::fmtflags __flags = __is.flags();
02384       __is.flags(__ios_base::dec | __ios_base::skipws);
02385 
02386       _RealType __m, __n;
02387       __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
02388       __x.param(typename fisher_f_distribution<_RealType>::
02389         param_type(__m, __n));
02390 
02391       __is.flags(__flags);
02392       return __is;
02393     }
02394 
02395 
02396   template<typename _RealType>
02397     template<typename _ForwardIterator,
02398          typename _UniformRandomNumberGenerator>
02399       void
02400       std::student_t_distribution<_RealType>::
02401       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02402               _UniformRandomNumberGenerator& __urng)
02403       {
02404     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02405     while (__f != __t)
02406       *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
02407       }
02408 
02409   template<typename _RealType>
02410     template<typename _ForwardIterator,
02411          typename _UniformRandomNumberGenerator>
02412       void
02413       std::student_t_distribution<_RealType>::
02414       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02415               _UniformRandomNumberGenerator& __urng,
02416               const param_type& __p)
02417       {
02418     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02419     typename std::gamma_distribution<result_type>::param_type
02420       __p2(__p.n() / 2, 2);
02421     while (__f != __t)
02422       *__f++ =  _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
02423       }
02424 
02425   template<typename _RealType, typename _CharT, typename _Traits>
02426     std::basic_ostream<_CharT, _Traits>&
02427     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02428            const student_t_distribution<_RealType>& __x)
02429     {
02430       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02431       typedef typename __ostream_type::ios_base    __ios_base;
02432 
02433       const typename __ios_base::fmtflags __flags = __os.flags();
02434       const _CharT __fill = __os.fill();
02435       const std::streamsize __precision = __os.precision();
02436       const _CharT __space = __os.widen(' ');
02437       __os.flags(__ios_base::scientific | __ios_base::left);
02438       __os.fill(__space);
02439       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02440 
02441       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
02442 
02443       __os.flags(__flags);
02444       __os.fill(__fill);
02445       __os.precision(__precision);
02446       return __os;
02447     }
02448 
02449   template<typename _RealType, typename _CharT, typename _Traits>
02450     std::basic_istream<_CharT, _Traits>&
02451     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02452            student_t_distribution<_RealType>& __x)
02453     {
02454       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02455       typedef typename __istream_type::ios_base    __ios_base;
02456 
02457       const typename __ios_base::fmtflags __flags = __is.flags();
02458       __is.flags(__ios_base::dec | __ios_base::skipws);
02459 
02460       _RealType __n;
02461       __is >> __n >> __x._M_nd >> __x._M_gd;
02462       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
02463 
02464       __is.flags(__flags);
02465       return __is;
02466     }
02467 
02468 
02469   template<typename _RealType>
02470     void
02471     gamma_distribution<_RealType>::param_type::
02472     _M_initialize()
02473     {
02474       _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
02475 
02476       const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
02477       _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
02478     }
02479 
02480   /**
02481    * Marsaglia, G. and Tsang, W. W.
02482    * "A Simple Method for Generating Gamma Variables"
02483    * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
02484    */
02485   template<typename _RealType>
02486     template<typename _UniformRandomNumberGenerator>
02487       typename gamma_distribution<_RealType>::result_type
02488       gamma_distribution<_RealType>::
02489       operator()(_UniformRandomNumberGenerator& __urng,
02490          const param_type& __param)
02491       {
02492     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02493       __aurng(__urng);
02494 
02495     result_type __u, __v, __n;
02496     const result_type __a1 = (__param._M_malpha
02497                   - _RealType(1.0) / _RealType(3.0));
02498 
02499     do
02500       {
02501         do
02502           {
02503         __n = _M_nd(__urng);
02504         __v = result_type(1.0) + __param._M_a2 * __n; 
02505           }
02506         while (__v <= 0.0);
02507 
02508         __v = __v * __v * __v;
02509         __u = __aurng();
02510       }
02511     while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
02512            && (std::log(__u) > (0.5 * __n * __n + __a1
02513                     * (1.0 - __v + std::log(__v)))));
02514 
02515     if (__param.alpha() == __param._M_malpha)
02516       return __a1 * __v * __param.beta();
02517     else
02518       {
02519         do
02520           __u = __aurng();
02521         while (__u == 0.0);
02522         
02523         return (std::pow(__u, result_type(1.0) / __param.alpha())
02524             * __a1 * __v * __param.beta());
02525       }
02526       }
02527 
02528   template<typename _RealType>
02529     template<typename _ForwardIterator,
02530          typename _UniformRandomNumberGenerator>
02531       void
02532       gamma_distribution<_RealType>::
02533       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02534               _UniformRandomNumberGenerator& __urng,
02535               const param_type& __param)
02536       {
02537     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02538     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02539       __aurng(__urng);
02540 
02541     result_type __u, __v, __n;
02542     const result_type __a1 = (__param._M_malpha
02543                   - _RealType(1.0) / _RealType(3.0));
02544 
02545     if (__param.alpha() == __param._M_malpha)
02546       while (__f != __t)
02547         {
02548           do
02549         {
02550           do
02551             {
02552               __n = _M_nd(__urng);
02553               __v = result_type(1.0) + __param._M_a2 * __n;
02554             }
02555           while (__v <= 0.0);
02556 
02557           __v = __v * __v * __v;
02558           __u = __aurng();
02559         }
02560           while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
02561              && (std::log(__u) > (0.5 * __n * __n + __a1
02562                       * (1.0 - __v + std::log(__v)))));
02563 
02564           *__f++ = __a1 * __v * __param.beta();
02565         }
02566     else
02567       while (__f != __t)
02568         {
02569           do
02570         {
02571           do
02572             {
02573               __n = _M_nd(__urng);
02574               __v = result_type(1.0) + __param._M_a2 * __n;
02575             }
02576           while (__v <= 0.0);
02577 
02578           __v = __v * __v * __v;
02579           __u = __aurng();
02580         }
02581           while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
02582              && (std::log(__u) > (0.5 * __n * __n + __a1
02583                       * (1.0 - __v + std::log(__v)))));
02584 
02585           do
02586         __u = __aurng();
02587           while (__u == 0.0);
02588 
02589           *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
02590             * __a1 * __v * __param.beta());
02591         }
02592       }
02593 
02594   template<typename _RealType, typename _CharT, typename _Traits>
02595     std::basic_ostream<_CharT, _Traits>&
02596     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02597            const gamma_distribution<_RealType>& __x)
02598     {
02599       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02600       typedef typename __ostream_type::ios_base    __ios_base;
02601 
02602       const typename __ios_base::fmtflags __flags = __os.flags();
02603       const _CharT __fill = __os.fill();
02604       const std::streamsize __precision = __os.precision();
02605       const _CharT __space = __os.widen(' ');
02606       __os.flags(__ios_base::scientific | __ios_base::left);
02607       __os.fill(__space);
02608       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02609 
02610       __os << __x.alpha() << __space << __x.beta()
02611        << __space << __x._M_nd;
02612 
02613       __os.flags(__flags);
02614       __os.fill(__fill);
02615       __os.precision(__precision);
02616       return __os;
02617     }
02618 
02619   template<typename _RealType, typename _CharT, typename _Traits>
02620     std::basic_istream<_CharT, _Traits>&
02621     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02622            gamma_distribution<_RealType>& __x)
02623     {
02624       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02625       typedef typename __istream_type::ios_base    __ios_base;
02626 
02627       const typename __ios_base::fmtflags __flags = __is.flags();
02628       __is.flags(__ios_base::dec | __ios_base::skipws);
02629 
02630       _RealType __alpha_val, __beta_val;
02631       __is >> __alpha_val >> __beta_val >> __x._M_nd;
02632       __x.param(typename gamma_distribution<_RealType>::
02633         param_type(__alpha_val, __beta_val));
02634 
02635       __is.flags(__flags);
02636       return __is;
02637     }
02638 
02639 
02640   template<typename _RealType>
02641     template<typename _UniformRandomNumberGenerator>
02642       typename weibull_distribution<_RealType>::result_type
02643       weibull_distribution<_RealType>::
02644       operator()(_UniformRandomNumberGenerator& __urng,
02645          const param_type& __p)
02646       {
02647     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02648       __aurng(__urng);
02649     return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
02650                   result_type(1) / __p.a());
02651       }
02652 
02653   template<typename _RealType>
02654     template<typename _ForwardIterator,
02655          typename _UniformRandomNumberGenerator>
02656       void
02657       weibull_distribution<_RealType>::
02658       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02659               _UniformRandomNumberGenerator& __urng,
02660               const param_type& __p)
02661       {
02662     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02663     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02664       __aurng(__urng);
02665     auto __inv_a = result_type(1) / __p.a();
02666 
02667     while (__f != __t)
02668       *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
02669                       __inv_a);
02670       }
02671 
02672   template<typename _RealType, typename _CharT, typename _Traits>
02673     std::basic_ostream<_CharT, _Traits>&
02674     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02675            const weibull_distribution<_RealType>& __x)
02676     {
02677       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02678       typedef typename __ostream_type::ios_base    __ios_base;
02679 
02680       const typename __ios_base::fmtflags __flags = __os.flags();
02681       const _CharT __fill = __os.fill();
02682       const std::streamsize __precision = __os.precision();
02683       const _CharT __space = __os.widen(' ');
02684       __os.flags(__ios_base::scientific | __ios_base::left);
02685       __os.fill(__space);
02686       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02687 
02688       __os << __x.a() << __space << __x.b();
02689 
02690       __os.flags(__flags);
02691       __os.fill(__fill);
02692       __os.precision(__precision);
02693       return __os;
02694     }
02695 
02696   template<typename _RealType, typename _CharT, typename _Traits>
02697     std::basic_istream<_CharT, _Traits>&
02698     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02699            weibull_distribution<_RealType>& __x)
02700     {
02701       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02702       typedef typename __istream_type::ios_base    __ios_base;
02703 
02704       const typename __ios_base::fmtflags __flags = __is.flags();
02705       __is.flags(__ios_base::dec | __ios_base::skipws);
02706 
02707       _RealType __a, __b;
02708       __is >> __a >> __b;
02709       __x.param(typename weibull_distribution<_RealType>::
02710         param_type(__a, __b));
02711 
02712       __is.flags(__flags);
02713       return __is;
02714     }
02715 
02716 
02717   template<typename _RealType>
02718     template<typename _UniformRandomNumberGenerator>
02719       typename extreme_value_distribution<_RealType>::result_type
02720       extreme_value_distribution<_RealType>::
02721       operator()(_UniformRandomNumberGenerator& __urng,
02722          const param_type& __p)
02723       {
02724     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02725       __aurng(__urng);
02726     return __p.a() - __p.b() * std::log(-std::log(result_type(1)
02727                               - __aurng()));
02728       }
02729 
02730   template<typename _RealType>
02731     template<typename _ForwardIterator,
02732          typename _UniformRandomNumberGenerator>
02733       void
02734       extreme_value_distribution<_RealType>::
02735       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02736               _UniformRandomNumberGenerator& __urng,
02737               const param_type& __p)
02738       {
02739     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02740     __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02741       __aurng(__urng);
02742 
02743     while (__f != __t)
02744       *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
02745                               - __aurng()));
02746       }
02747 
02748   template<typename _RealType, typename _CharT, typename _Traits>
02749     std::basic_ostream<_CharT, _Traits>&
02750     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02751            const extreme_value_distribution<_RealType>& __x)
02752     {
02753       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02754       typedef typename __ostream_type::ios_base    __ios_base;
02755 
02756       const typename __ios_base::fmtflags __flags = __os.flags();
02757       const _CharT __fill = __os.fill();
02758       const std::streamsize __precision = __os.precision();
02759       const _CharT __space = __os.widen(' ');
02760       __os.flags(__ios_base::scientific | __ios_base::left);
02761       __os.fill(__space);
02762       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02763 
02764       __os << __x.a() << __space << __x.b();
02765 
02766       __os.flags(__flags);
02767       __os.fill(__fill);
02768       __os.precision(__precision);
02769       return __os;
02770     }
02771 
02772   template<typename _RealType, typename _CharT, typename _Traits>
02773     std::basic_istream<_CharT, _Traits>&
02774     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02775            extreme_value_distribution<_RealType>& __x)
02776     {
02777       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02778       typedef typename __istream_type::ios_base    __ios_base;
02779 
02780       const typename __ios_base::fmtflags __flags = __is.flags();
02781       __is.flags(__ios_base::dec | __ios_base::skipws);
02782 
02783       _RealType __a, __b;
02784       __is >> __a >> __b;
02785       __x.param(typename extreme_value_distribution<_RealType>::
02786         param_type(__a, __b));
02787 
02788       __is.flags(__flags);
02789       return __is;
02790     }
02791 
02792 
02793   template<typename _IntType>
02794     void
02795     discrete_distribution<_IntType>::param_type::
02796     _M_initialize()
02797     {
02798       if (_M_prob.size() < 2)
02799     {
02800       _M_prob.clear();
02801       return;
02802     }
02803 
02804       const double __sum = std::accumulate(_M_prob.begin(),
02805                        _M_prob.end(), 0.0);
02806       // Now normalize the probabilites.
02807       __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
02808                 __sum);
02809       // Accumulate partial sums.
02810       _M_cp.reserve(_M_prob.size());
02811       std::partial_sum(_M_prob.begin(), _M_prob.end(),
02812                std::back_inserter(_M_cp));
02813       // Make sure the last cumulative probability is one.
02814       _M_cp[_M_cp.size() - 1] = 1.0;
02815     }
02816 
02817   template<typename _IntType>
02818     template<typename _Func>
02819       discrete_distribution<_IntType>::param_type::
02820       param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
02821       : _M_prob(), _M_cp()
02822       {
02823     const size_t __n = __nw == 0 ? 1 : __nw;
02824     const double __delta = (__xmax - __xmin) / __n;
02825 
02826     _M_prob.reserve(__n);
02827     for (size_t __k = 0; __k < __nw; ++__k)
02828       _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
02829 
02830     _M_initialize();
02831       }
02832 
02833   template<typename _IntType>
02834     template<typename _UniformRandomNumberGenerator>
02835       typename discrete_distribution<_IntType>::result_type
02836       discrete_distribution<_IntType>::
02837       operator()(_UniformRandomNumberGenerator& __urng,
02838          const param_type& __param)
02839       {
02840     if (__param._M_cp.empty())
02841       return result_type(0);
02842 
02843     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02844       __aurng(__urng);
02845 
02846     const double __p = __aurng();
02847     auto __pos = std::lower_bound(__param._M_cp.begin(),
02848                       __param._M_cp.end(), __p);
02849 
02850     return __pos - __param._M_cp.begin();
02851       }
02852 
02853   template<typename _IntType>
02854     template<typename _ForwardIterator,
02855          typename _UniformRandomNumberGenerator>
02856       void
02857       discrete_distribution<_IntType>::
02858       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02859               _UniformRandomNumberGenerator& __urng,
02860               const param_type& __param)
02861       {
02862     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02863 
02864     if (__param._M_cp.empty())
02865       {
02866         while (__f != __t)
02867           *__f++ = result_type(0);
02868         return;
02869       }
02870 
02871     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02872       __aurng(__urng);
02873 
02874     while (__f != __t)
02875       {
02876         const double __p = __aurng();
02877         auto __pos = std::lower_bound(__param._M_cp.begin(),
02878                       __param._M_cp.end(), __p);
02879 
02880         *__f++ = __pos - __param._M_cp.begin();
02881       }
02882       }
02883 
02884   template<typename _IntType, typename _CharT, typename _Traits>
02885     std::basic_ostream<_CharT, _Traits>&
02886     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02887            const discrete_distribution<_IntType>& __x)
02888     {
02889       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02890       typedef typename __ostream_type::ios_base    __ios_base;
02891 
02892       const typename __ios_base::fmtflags __flags = __os.flags();
02893       const _CharT __fill = __os.fill();
02894       const std::streamsize __precision = __os.precision();
02895       const _CharT __space = __os.widen(' ');
02896       __os.flags(__ios_base::scientific | __ios_base::left);
02897       __os.fill(__space);
02898       __os.precision(std::numeric_limits<double>::max_digits10);
02899 
02900       std::vector<double> __prob = __x.probabilities();
02901       __os << __prob.size();
02902       for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
02903     __os << __space << *__dit;
02904 
02905       __os.flags(__flags);
02906       __os.fill(__fill);
02907       __os.precision(__precision);
02908       return __os;
02909     }
02910 
02911   template<typename _IntType, typename _CharT, typename _Traits>
02912     std::basic_istream<_CharT, _Traits>&
02913     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02914            discrete_distribution<_IntType>& __x)
02915     {
02916       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02917       typedef typename __istream_type::ios_base    __ios_base;
02918 
02919       const typename __ios_base::fmtflags __flags = __is.flags();
02920       __is.flags(__ios_base::dec | __ios_base::skipws);
02921 
02922       size_t __n;
02923       __is >> __n;
02924 
02925       std::vector<double> __prob_vec;
02926       __prob_vec.reserve(__n);
02927       for (; __n != 0; --__n)
02928     {
02929       double __prob;
02930       __is >> __prob;
02931       __prob_vec.push_back(__prob);
02932     }
02933 
02934       __x.param(typename discrete_distribution<_IntType>::
02935         param_type(__prob_vec.begin(), __prob_vec.end()));
02936 
02937       __is.flags(__flags);
02938       return __is;
02939     }
02940 
02941 
02942   template<typename _RealType>
02943     void
02944     piecewise_constant_distribution<_RealType>::param_type::
02945     _M_initialize()
02946     {
02947       if (_M_int.size() < 2
02948       || (_M_int.size() == 2
02949           && _M_int[0] == _RealType(0)
02950           && _M_int[1] == _RealType(1)))
02951     {
02952       _M_int.clear();
02953       _M_den.clear();
02954       return;
02955     }
02956 
02957       const double __sum = std::accumulate(_M_den.begin(),
02958                        _M_den.end(), 0.0);
02959 
02960       __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
02961                 __sum);
02962 
02963       _M_cp.reserve(_M_den.size());
02964       std::partial_sum(_M_den.begin(), _M_den.end(),
02965                std::back_inserter(_M_cp));
02966 
02967       // Make sure the last cumulative probability is one.
02968       _M_cp[_M_cp.size() - 1] = 1.0;
02969 
02970       for (size_t __k = 0; __k < _M_den.size(); ++__k)
02971     _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
02972     }
02973 
02974   template<typename _RealType>
02975     template<typename _InputIteratorB, typename _InputIteratorW>
02976       piecewise_constant_distribution<_RealType>::param_type::
02977       param_type(_InputIteratorB __bbegin,
02978          _InputIteratorB __bend,
02979          _InputIteratorW __wbegin)
02980       : _M_int(), _M_den(), _M_cp()
02981       {
02982     if (__bbegin != __bend)
02983       {
02984         for (;;)
02985           {
02986         _M_int.push_back(*__bbegin);
02987         ++__bbegin;
02988         if (__bbegin == __bend)
02989           break;
02990 
02991         _M_den.push_back(*__wbegin);
02992         ++__wbegin;
02993           }
02994       }
02995 
02996     _M_initialize();
02997       }
02998 
02999   template<typename _RealType>
03000     template<typename _Func>
03001       piecewise_constant_distribution<_RealType>::param_type::
03002       param_type(initializer_list<_RealType> __bl, _Func __fw)
03003       : _M_int(), _M_den(), _M_cp()
03004       {
03005     _M_int.reserve(__bl.size());
03006     for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
03007       _M_int.push_back(*__biter);
03008 
03009     _M_den.reserve(_M_int.size() - 1);
03010     for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
03011       _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
03012 
03013     _M_initialize();
03014       }
03015 
03016   template<typename _RealType>
03017     template<typename _Func>
03018       piecewise_constant_distribution<_RealType>::param_type::
03019       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
03020       : _M_int(), _M_den(), _M_cp()
03021       {
03022     const size_t __n = __nw == 0 ? 1 : __nw;
03023     const _RealType __delta = (__xmax - __xmin) / __n;
03024 
03025     _M_int.reserve(__n + 1);
03026     for (size_t __k = 0; __k <= __nw; ++__k)
03027       _M_int.push_back(__xmin + __k * __delta);
03028 
03029     _M_den.reserve(__n);
03030     for (size_t __k = 0; __k < __nw; ++__k)
03031       _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
03032 
03033     _M_initialize();
03034       }
03035 
03036   template<typename _RealType>
03037     template<typename _UniformRandomNumberGenerator>
03038       typename piecewise_constant_distribution<_RealType>::result_type
03039       piecewise_constant_distribution<_RealType>::
03040       operator()(_UniformRandomNumberGenerator& __urng,
03041          const param_type& __param)
03042       {
03043     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
03044       __aurng(__urng);
03045 
03046     const double __p = __aurng();
03047     if (__param._M_cp.empty())
03048       return __p;
03049 
03050     auto __pos = std::lower_bound(__param._M_cp.begin(),
03051                       __param._M_cp.end(), __p);
03052     const size_t __i = __pos - __param._M_cp.begin();
03053 
03054     const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
03055 
03056     return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
03057       }
03058 
03059   template<typename _RealType>
03060     template<typename _ForwardIterator,
03061          typename _UniformRandomNumberGenerator>
03062       void
03063       piecewise_constant_distribution<_RealType>::
03064       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
03065               _UniformRandomNumberGenerator& __urng,
03066               const param_type& __param)
03067       {
03068     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
03069     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
03070       __aurng(__urng);
03071 
03072     if (__param._M_cp.empty())
03073       {
03074         while (__f != __t)
03075           *__f++ = __aurng();
03076         return;
03077       }
03078 
03079     while (__f != __t)
03080       {
03081         const double __p = __aurng();
03082 
03083         auto __pos = std::lower_bound(__param._M_cp.begin(),
03084                       __param._M_cp.end(), __p);
03085         const size_t __i = __pos - __param._M_cp.begin();
03086 
03087         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
03088 
03089         *__f++ = (__param._M_int[__i]
03090               + (__p - __pref) / __param._M_den[__i]);
03091       }
03092       }
03093 
03094   template<typename _RealType, typename _CharT, typename _Traits>
03095     std::basic_ostream<_CharT, _Traits>&
03096     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
03097            const piecewise_constant_distribution<_RealType>& __x)
03098     {
03099       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
03100       typedef typename __ostream_type::ios_base    __ios_base;
03101 
03102       const typename __ios_base::fmtflags __flags = __os.flags();
03103       const _CharT __fill = __os.fill();
03104       const std::streamsize __precision = __os.precision();
03105       const _CharT __space = __os.widen(' ');
03106       __os.flags(__ios_base::scientific | __ios_base::left);
03107       __os.fill(__space);
03108       __os.precision(std::numeric_limits<_RealType>::max_digits10);
03109 
03110       std::vector<_RealType> __int = __x.intervals();
03111       __os << __int.size() - 1;
03112 
03113       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
03114     __os << __space << *__xit;
03115 
03116       std::vector<double> __den = __x.densities();
03117       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
03118     __os << __space << *__dit;
03119 
03120       __os.flags(__flags);
03121       __os.fill(__fill);
03122       __os.precision(__precision);
03123       return __os;
03124     }
03125 
03126   template<typename _RealType, typename _CharT, typename _Traits>
03127     std::basic_istream<_CharT, _Traits>&
03128     operator>>(std::basic_istream<_CharT, _Traits>& __is,
03129            piecewise_constant_distribution<_RealType>& __x)
03130     {
03131       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
03132       typedef typename __istream_type::ios_base    __ios_base;
03133 
03134       const typename __ios_base::fmtflags __flags = __is.flags();
03135       __is.flags(__ios_base::dec | __ios_base::skipws);
03136 
03137       size_t __n;
03138       __is >> __n;
03139 
03140       std::vector<_RealType> __int_vec;
03141       __int_vec.reserve(__n + 1);
03142       for (size_t __i = 0; __i <= __n; ++__i)
03143     {
03144       _RealType __int;
03145       __is >> __int;
03146       __int_vec.push_back(__int);
03147     }
03148 
03149       std::vector<double> __den_vec;
03150       __den_vec.reserve(__n);
03151       for (size_t __i = 0; __i < __n; ++__i)
03152     {
03153       double __den;
03154       __is >> __den;
03155       __den_vec.push_back(__den);
03156     }
03157 
03158       __x.param(typename piecewise_constant_distribution<_RealType>::
03159       param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
03160 
03161       __is.flags(__flags);
03162       return __is;
03163     }
03164 
03165 
03166   template<typename _RealType>
03167     void
03168     piecewise_linear_distribution<_RealType>::param_type::
03169     _M_initialize()
03170     {
03171       if (_M_int.size() < 2
03172       || (_M_int.size() == 2
03173           && _M_int[0] == _RealType(0)
03174           && _M_int[1] == _RealType(1)
03175           && _M_den[0] == _M_den[1]))
03176     {
03177       _M_int.clear();
03178       _M_den.clear();
03179       return;
03180     }
03181 
03182       double __sum = 0.0;
03183       _M_cp.reserve(_M_int.size() - 1);
03184       _M_m.reserve(_M_int.size() - 1);
03185       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
03186     {
03187       const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
03188       __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
03189       _M_cp.push_back(__sum);
03190       _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
03191     }
03192 
03193       //  Now normalize the densities...
03194       __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
03195                 __sum);
03196       //  ... and partial sums... 
03197       __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
03198       //  ... and slopes.
03199       __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
03200 
03201       //  Make sure the last cumulative probablility is one.
03202       _M_cp[_M_cp.size() - 1] = 1.0;
03203      }
03204 
03205   template<typename _RealType>
03206     template<typename _InputIteratorB, typename _InputIteratorW>
03207       piecewise_linear_distribution<_RealType>::param_type::
03208       param_type(_InputIteratorB __bbegin,
03209          _InputIteratorB __bend,
03210          _InputIteratorW __wbegin)
03211       : _M_int(), _M_den(), _M_cp(), _M_m()
03212       {
03213     for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
03214       {
03215         _M_int.push_back(*__bbegin);
03216         _M_den.push_back(*__wbegin);
03217       }
03218 
03219     _M_initialize();
03220       }
03221 
03222   template<typename _RealType>
03223     template<typename _Func>
03224       piecewise_linear_distribution<_RealType>::param_type::
03225       param_type(initializer_list<_RealType> __bl, _Func __fw)
03226       : _M_int(), _M_den(), _M_cp(), _M_m()
03227       {
03228     _M_int.reserve(__bl.size());
03229     _M_den.reserve(__bl.size());
03230     for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
03231       {
03232         _M_int.push_back(*__biter);
03233         _M_den.push_back(__fw(*__biter));
03234       }
03235 
03236     _M_initialize();
03237       }
03238 
03239   template<typename _RealType>
03240     template<typename _Func>
03241       piecewise_linear_distribution<_RealType>::param_type::
03242       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
03243       : _M_int(), _M_den(), _M_cp(), _M_m()
03244       {
03245     const size_t __n = __nw == 0 ? 1 : __nw;
03246     const _RealType __delta = (__xmax - __xmin) / __n;
03247 
03248     _M_int.reserve(__n + 1);
03249     _M_den.reserve(__n + 1);
03250     for (size_t __k = 0; __k <= __nw; ++__k)
03251       {
03252         _M_int.push_back(__xmin + __k * __delta);
03253         _M_den.push_back(__fw(_M_int[__k] + __delta));
03254       }
03255 
03256     _M_initialize();
03257       }
03258 
03259   template<typename _RealType>
03260     template<typename _UniformRandomNumberGenerator>
03261       typename piecewise_linear_distribution<_RealType>::result_type
03262       piecewise_linear_distribution<_RealType>::
03263       operator()(_UniformRandomNumberGenerator& __urng,
03264          const param_type& __param)
03265       {
03266     __detail::_Adaptor<_UniformRandomNumberGenerator, double>
03267       __aurng(__urng);
03268 
03269     const double __p = __aurng();
03270     if (__param._M_cp.empty())
03271       return __p;
03272 
03273     auto __pos = std::lower_bound(__param._M_cp.begin(),
03274                       __param._M_cp.end(), __p);
03275     const size_t __i = __pos - __param._M_cp.begin();
03276 
03277     const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
03278 
03279     const double __a = 0.5 * __param._M_m[__i];
03280     const double __b = __param._M_den[__i];
03281     const double __cm = __p - __pref;
03282 
03283     _RealType __x = __param._M_int[__i];
03284     if (__a == 0)
03285       __x += __cm / __b;
03286     else
03287       {
03288         const double __d = __b * __b + 4.0 * __a * __cm;
03289         __x += 0.5 * (std::sqrt(__d) - __b) / __a;
03290           }
03291 
03292         return __x;
03293       }
03294 
03295   template<typename _RealType>
03296     template<typename _ForwardIterator,
03297          typename _UniformRandomNumberGenerator>
03298       void
03299       piecewise_linear_distribution<_RealType>::
03300       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
03301               _UniformRandomNumberGenerator& __urng,
03302               const param_type& __param)
03303       {
03304     __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
03305     // We could duplicate everything from operator()...
03306     while (__f != __t)
03307       *__f++ = this->operator()(__urng, __param);
03308       }
03309 
03310   template<typename _RealType, typename _CharT, typename _Traits>
03311     std::basic_ostream<_CharT, _Traits>&
03312     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
03313            const piecewise_linear_distribution<_RealType>& __x)
03314     {
03315       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
03316       typedef typename __ostream_type::ios_base    __ios_base;
03317 
03318       const typename __ios_base::fmtflags __flags = __os.flags();
03319       const _CharT __fill = __os.fill();
03320       const std::streamsize __precision = __os.precision();
03321       const _CharT __space = __os.widen(' ');
03322       __os.flags(__ios_base::scientific | __ios_base::left);
03323       __os.fill(__space);
03324       __os.precision(std::numeric_limits<_RealType>::max_digits10);
03325 
03326       std::vector<_RealType> __int = __x.intervals();
03327       __os << __int.size() - 1;
03328 
03329       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
03330     __os << __space << *__xit;
03331 
03332       std::vector<double> __den = __x.densities();
03333       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
03334     __os << __space << *__dit;
03335 
03336       __os.flags(__flags);
03337       __os.fill(__fill);
03338       __os.precision(__precision);
03339       return __os;
03340     }
03341 
03342   template<typename _RealType, typename _CharT, typename _Traits>
03343     std::basic_istream<_CharT, _Traits>&
03344     operator>>(std::basic_istream<_CharT, _Traits>& __is,
03345            piecewise_linear_distribution<_RealType>& __x)
03346     {
03347       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
03348       typedef typename __istream_type::ios_base    __ios_base;
03349 
03350       const typename __ios_base::fmtflags __flags = __is.flags();
03351       __is.flags(__ios_base::dec | __ios_base::skipws);
03352 
03353       size_t __n;
03354       __is >> __n;
03355 
03356       std::vector<_RealType> __int_vec;
03357       __int_vec.reserve(__n + 1);
03358       for (size_t __i = 0; __i <= __n; ++__i)
03359     {
03360       _RealType __int;
03361       __is >> __int;
03362       __int_vec.push_back(__int);
03363     }
03364 
03365       std::vector<double> __den_vec;
03366       __den_vec.reserve(__n + 1);
03367       for (size_t __i = 0; __i <= __n; ++__i)
03368     {
03369       double __den;
03370       __is >> __den;
03371       __den_vec.push_back(__den);
03372     }
03373 
03374       __x.param(typename piecewise_linear_distribution<_RealType>::
03375       param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
03376 
03377       __is.flags(__flags);
03378       return __is;
03379     }
03380 
03381 
03382   template<typename _IntType>
03383     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
03384     {
03385       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
03386     _M_v.push_back(__detail::__mod<result_type,
03387                __detail::_Shift<result_type, 32>::__value>(*__iter));
03388     }
03389 
03390   template<typename _InputIterator>
03391     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
03392     {
03393       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
03394     _M_v.push_back(__detail::__mod<result_type,
03395                __detail::_Shift<result_type, 32>::__value>(*__iter));
03396     }
03397 
03398   template<typename _RandomAccessIterator>
03399     void
03400     seed_seq::generate(_RandomAccessIterator __begin,
03401                _RandomAccessIterator __end)
03402     {
03403       typedef typename iterator_traits<_RandomAccessIterator>::value_type
03404         _Type;
03405 
03406       if (__begin == __end)
03407     return;
03408 
03409       std::fill(__begin, __end, _Type(0x8b8b8b8bu));
03410 
03411       const size_t __n = __end - __begin;
03412       const size_t __s = _M_v.size();
03413       const size_t __t = (__n >= 623) ? 11
03414                : (__n >=  68) ? 7
03415                : (__n >=  39) ? 5
03416                : (__n >=   7) ? 3
03417                : (__n - 1) / 2;
03418       const size_t __p = (__n - __t) / 2;
03419       const size_t __q = __p + __t;
03420       const size_t __m = std::max(size_t(__s + 1), __n);
03421 
03422       for (size_t __k = 0; __k < __m; ++__k)
03423     {
03424       _Type __arg = (__begin[__k % __n]
03425              ^ __begin[(__k + __p) % __n]
03426              ^ __begin[(__k - 1) % __n]);
03427       _Type __r1 = __arg ^ (__arg >> 27);
03428       __r1 = __detail::__mod<_Type,
03429             __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
03430       _Type __r2 = __r1;
03431       if (__k == 0)
03432         __r2 += __s;
03433       else if (__k <= __s)
03434         __r2 += __k % __n + _M_v[__k - 1];
03435       else
03436         __r2 += __k % __n;
03437       __r2 = __detail::__mod<_Type,
03438                __detail::_Shift<_Type, 32>::__value>(__r2);
03439       __begin[(__k + __p) % __n] += __r1;
03440       __begin[(__k + __q) % __n] += __r2;
03441       __begin[__k % __n] = __r2;
03442     }
03443 
03444       for (size_t __k = __m; __k < __m + __n; ++__k)
03445     {
03446       _Type __arg = (__begin[__k % __n]
03447              + __begin[(__k + __p) % __n]
03448              + __begin[(__k - 1) % __n]);
03449       _Type __r3 = __arg ^ (__arg >> 27);
03450       __r3 = __detail::__mod<_Type,
03451            __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
03452       _Type __r4 = __r3 - __k % __n;
03453       __r4 = __detail::__mod<_Type,
03454                __detail::_Shift<_Type, 32>::__value>(__r4);
03455       __begin[(__k + __p) % __n] ^= __r3;
03456       __begin[(__k + __q) % __n] ^= __r4;
03457       __begin[__k % __n] = __r4;
03458     }
03459     }
03460 
03461   template<typename _RealType, size_t __bits,
03462        typename _UniformRandomNumberGenerator>
03463     _RealType
03464     generate_canonical(_UniformRandomNumberGenerator& __urng)
03465     {
03466       const size_t __b
03467     = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
03468                    __bits);
03469       const long double __r = static_cast<long double>(__urng.max())
03470                 - static_cast<long double>(__urng.min()) + 1.0L;
03471       const size_t __log2r = std::log(__r) / std::log(2.0L);
03472       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
03473       _RealType __sum = _RealType(0);
03474       _RealType __tmp = _RealType(1);
03475       for (; __k != 0; --__k)
03476     {
03477       __sum += _RealType(__urng() - __urng.min()) * __tmp;
03478       __tmp *= __r;
03479     }
03480       return __sum / __tmp;
03481     }
03482 
03483 _GLIBCXX_END_NAMESPACE_VERSION
03484 } // namespace
03485 
03486 #endif