libstdc++
random.tcc
Go to the documentation of this file.
00001 // random number generation (out of line) -*- C++ -*-
00002 
00003 // Copyright (C) 2009-2015 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<double>::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       static_assert(std::is_floating_point<_RealType>::value,
03467                     "template argument not a floating point type");
03468 
03469       const size_t __b
03470         = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
03471                    __bits);
03472       const long double __r = static_cast<long double>(__urng.max())
03473                             - static_cast<long double>(__urng.min()) + 1.0L;
03474       const size_t __log2r = std::log(__r) / std::log(2.0L);
03475       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
03476       _RealType __sum = _RealType(0);
03477       _RealType __tmp = _RealType(1);
03478       for (; __k != 0; --__k)
03479         {
03480           __sum += _RealType(__urng() - __urng.min()) * __tmp;
03481           __tmp *= __r;
03482         }
03483       return __sum / __tmp;
03484     }
03485 
03486 _GLIBCXX_END_NAMESPACE_VERSION
03487 } // namespace
03488 
03489 #endif