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