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