libstdc++
random.tcc
Go to the documentation of this file.
00001 // random number generation (out of line) -*- C++ -*-
00002 
00003 // Copyright (C) 2009-2019 Free Software Foundation, Inc.
00004 //
00005 // This file is part of the GNU ISO C++ Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the
00007 // terms of the GNU General Public License as published by the
00008 // Free Software Foundation; either version 3, or (at your option)
00009 // any later version.
00010 
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 
00016 // Under Section 7 of GPL version 3, you are granted additional
00017 // permissions described in the GCC Runtime Library Exception, version
00018 // 3.1, as published by the Free Software Foundation.
00019 
00020 // You should have received a copy of the GNU General Public License and
00021 // a copy of the GCC Runtime Library Exception along with this program;
00022 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
00023 // <http://www.gnu.org/licenses/>.
00024 
00025 /** @file bits/random.tcc
00026  *  This is an internal header file, included by other library headers.
00027  *  Do not attempt to use it directly. @headername{random}
00028  */
00029 
00030 #ifndef _RANDOM_TCC
00031 #define _RANDOM_TCC 1
00032 
00033 #include <numeric> // std::accumulate and std::partial_sum
00034 
00035 namespace std _GLIBCXX_VISIBILITY(default)
00036 {
00037 _GLIBCXX_BEGIN_NAMESPACE_VERSION
00038 
00039   /*
00040    * (Further) implementation-space details.
00041    */
00042   namespace __detail
00043   {
00044     // General case for x = (ax + c) mod m -- use Schrage's algorithm
00045     // to avoid integer overflow.
00046     //
00047     // Preconditions:  a > 0, m > 0.
00048     //
00049     // Note: only works correctly for __m % __a < __m / __a.
00050     template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
00051       _Tp
00052       _Mod<_Tp, __m, __a, __c, false, true>::
00053       __calc(_Tp __x)
00054       {
00055         if (__a == 1)
00056           __x %= __m;
00057         else
00058           {
00059             static const _Tp __q = __m / __a;
00060             static const _Tp __r = __m % __a;
00061 
00062             _Tp __t1 = __a * (__x % __q);
00063             _Tp __t2 = __r * (__x / __q);
00064             if (__t1 >= __t2)
00065               __x = __t1 - __t2;
00066             else
00067               __x = __m - __t2 + __t1;
00068           }
00069 
00070         if (__c != 0)
00071           {
00072             const _Tp __d = __m - __x;
00073             if (__d > __c)
00074               __x += __c;
00075             else
00076               __x = __c - __d;
00077           }
00078         return __x;
00079       }
00080 
00081     template<typename _InputIterator, typename _OutputIterator,
00082              typename _Tp>
00083       _OutputIterator
00084       __normalize(_InputIterator __first, _InputIterator __last,
00085                   _OutputIterator __result, const _Tp& __factor)
00086       {
00087         for (; __first != __last; ++__first, ++__result)
00088           *__result = *__first / __factor;
00089         return __result;
00090       }
00091 
00092   } // namespace __detail
00093 
00094   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00095     constexpr _UIntType
00096     linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
00097 
00098   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00099     constexpr _UIntType
00100     linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
00101 
00102   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00103     constexpr _UIntType
00104     linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
00105 
00106   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00107     constexpr _UIntType
00108     linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
00109 
00110   /**
00111    * Seeds the LCR with integral value @p __s, adjusted so that the
00112    * ring identity is never a member of the convergence set.
00113    */
00114   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00115     void
00116     linear_congruential_engine<_UIntType, __a, __c, __m>::
00117     seed(result_type __s)
00118     {
00119       if ((__detail::__mod<_UIntType, __m>(__c) == 0)
00120           && (__detail::__mod<_UIntType, __m>(__s) == 0))
00121         _M_x = 1;
00122       else
00123         _M_x = __detail::__mod<_UIntType, __m>(__s);
00124     }
00125 
00126   /**
00127    * Seeds the LCR engine with a value generated by @p __q.
00128    */
00129   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
00130     template<typename _Sseq>
00131       auto
00132       linear_congruential_engine<_UIntType, __a, __c, __m>::
00133       seed(_Sseq& __q)
00134       -> _If_seed_seq<_Sseq>
00135       {
00136         const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
00137                                         : std::__lg(__m);
00138         const _UIntType __k = (__k0 + 31) / 32;
00139         uint_least32_t __arr[__k + 3];
00140         __q.generate(__arr + 0, __arr + __k + 3);
00141         _UIntType __factor = 1u;
00142         _UIntType __sum = 0u;
00143         for (size_t __j = 0; __j < __k; ++__j)
00144           {
00145             __sum += __arr[__j + 3] * __factor;
00146             __factor *= __detail::_Shift<_UIntType, 32>::__value;
00147           }
00148         seed(__sum);
00149       }
00150 
00151   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
00152            typename _CharT, typename _Traits>
00153     std::basic_ostream<_CharT, _Traits>&
00154     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00155                const linear_congruential_engine<_UIntType,
00156                                                 __a, __c, __m>& __lcr)
00157     {
00158       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00159       typedef typename __ostream_type::ios_base    __ios_base;
00160 
00161       const typename __ios_base::fmtflags __flags = __os.flags();
00162       const _CharT __fill = __os.fill();
00163       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00164       __os.fill(__os.widen(' '));
00165 
00166       __os << __lcr._M_x;
00167 
00168       __os.flags(__flags);
00169       __os.fill(__fill);
00170       return __os;
00171     }
00172 
00173   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
00174            typename _CharT, typename _Traits>
00175     std::basic_istream<_CharT, _Traits>&
00176     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00177                linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
00178     {
00179       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00180       typedef typename __istream_type::ios_base    __ios_base;
00181 
00182       const typename __ios_base::fmtflags __flags = __is.flags();
00183       __is.flags(__ios_base::dec);
00184 
00185       __is >> __lcr._M_x;
00186 
00187       __is.flags(__flags);
00188       return __is;
00189     }
00190 
00191 
00192   template<typename _UIntType,
00193            size_t __w, size_t __n, size_t __m, size_t __r,
00194            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00195            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00196            _UIntType __f>
00197     constexpr size_t
00198     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00199                             __s, __b, __t, __c, __l, __f>::word_size;
00200 
00201   template<typename _UIntType,
00202            size_t __w, size_t __n, size_t __m, size_t __r,
00203            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00204            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00205            _UIntType __f>
00206     constexpr size_t
00207     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00208                             __s, __b, __t, __c, __l, __f>::state_size;
00209 
00210   template<typename _UIntType,
00211            size_t __w, size_t __n, size_t __m, size_t __r,
00212            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00213            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00214            _UIntType __f>
00215     constexpr size_t
00216     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00217                             __s, __b, __t, __c, __l, __f>::shift_size;
00218 
00219   template<typename _UIntType,
00220            size_t __w, size_t __n, size_t __m, size_t __r,
00221            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00222            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00223            _UIntType __f>
00224     constexpr size_t
00225     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00226                             __s, __b, __t, __c, __l, __f>::mask_bits;
00227 
00228   template<typename _UIntType,
00229            size_t __w, size_t __n, size_t __m, size_t __r,
00230            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00231            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00232            _UIntType __f>
00233     constexpr _UIntType
00234     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00235                             __s, __b, __t, __c, __l, __f>::xor_mask;
00236 
00237   template<typename _UIntType,
00238            size_t __w, size_t __n, size_t __m, size_t __r,
00239            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00240            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00241            _UIntType __f>
00242     constexpr size_t
00243     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00244                             __s, __b, __t, __c, __l, __f>::tempering_u;
00245    
00246   template<typename _UIntType,
00247            size_t __w, size_t __n, size_t __m, size_t __r,
00248            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00249            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00250            _UIntType __f>
00251     constexpr _UIntType
00252     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00253                             __s, __b, __t, __c, __l, __f>::tempering_d;
00254 
00255   template<typename _UIntType,
00256            size_t __w, size_t __n, size_t __m, size_t __r,
00257            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00258            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00259            _UIntType __f>
00260     constexpr size_t
00261     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00262                             __s, __b, __t, __c, __l, __f>::tempering_s;
00263 
00264   template<typename _UIntType,
00265            size_t __w, size_t __n, size_t __m, size_t __r,
00266            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00267            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00268            _UIntType __f>
00269     constexpr _UIntType
00270     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00271                             __s, __b, __t, __c, __l, __f>::tempering_b;
00272 
00273   template<typename _UIntType,
00274            size_t __w, size_t __n, size_t __m, size_t __r,
00275            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00276            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00277            _UIntType __f>
00278     constexpr size_t
00279     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00280                             __s, __b, __t, __c, __l, __f>::tempering_t;
00281 
00282   template<typename _UIntType,
00283            size_t __w, size_t __n, size_t __m, size_t __r,
00284            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00285            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00286            _UIntType __f>
00287     constexpr _UIntType
00288     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00289                             __s, __b, __t, __c, __l, __f>::tempering_c;
00290 
00291   template<typename _UIntType,
00292            size_t __w, size_t __n, size_t __m, size_t __r,
00293            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00294            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00295            _UIntType __f>
00296     constexpr size_t
00297     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00298                             __s, __b, __t, __c, __l, __f>::tempering_l;
00299 
00300   template<typename _UIntType,
00301            size_t __w, size_t __n, size_t __m, size_t __r,
00302            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00303            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00304            _UIntType __f>
00305     constexpr _UIntType
00306     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00307                             __s, __b, __t, __c, __l, __f>::
00308                                               initialization_multiplier;
00309 
00310   template<typename _UIntType,
00311            size_t __w, size_t __n, size_t __m, size_t __r,
00312            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00313            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00314            _UIntType __f>
00315     constexpr _UIntType
00316     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00317                             __s, __b, __t, __c, __l, __f>::default_seed;
00318 
00319   template<typename _UIntType,
00320            size_t __w, size_t __n, size_t __m, size_t __r,
00321            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00322            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00323            _UIntType __f>
00324     void
00325     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00326                             __s, __b, __t, __c, __l, __f>::
00327     seed(result_type __sd)
00328     {
00329       _M_x[0] = __detail::__mod<_UIntType,
00330         __detail::_Shift<_UIntType, __w>::__value>(__sd);
00331 
00332       for (size_t __i = 1; __i < state_size; ++__i)
00333         {
00334           _UIntType __x = _M_x[__i - 1];
00335           __x ^= __x >> (__w - 2);
00336           __x *= __f;
00337           __x += __detail::__mod<_UIntType, __n>(__i);
00338           _M_x[__i] = __detail::__mod<_UIntType,
00339             __detail::_Shift<_UIntType, __w>::__value>(__x);
00340         }
00341       _M_p = state_size;
00342     }
00343 
00344   template<typename _UIntType,
00345            size_t __w, size_t __n, size_t __m, size_t __r,
00346            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00347            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00348            _UIntType __f>
00349     template<typename _Sseq>
00350       auto
00351       mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00352                               __s, __b, __t, __c, __l, __f>::
00353       seed(_Sseq& __q)
00354       -> _If_seed_seq<_Sseq>
00355       {
00356         const _UIntType __upper_mask = (~_UIntType()) << __r;
00357         const size_t __k = (__w + 31) / 32;
00358         uint_least32_t __arr[__n * __k];
00359         __q.generate(__arr + 0, __arr + __n * __k);
00360 
00361         bool __zero = true;
00362         for (size_t __i = 0; __i < state_size; ++__i)
00363           {
00364             _UIntType __factor = 1u;
00365             _UIntType __sum = 0u;
00366             for (size_t __j = 0; __j < __k; ++__j)
00367               {
00368                 __sum += __arr[__k * __i + __j] * __factor;
00369                 __factor *= __detail::_Shift<_UIntType, 32>::__value;
00370               }
00371             _M_x[__i] = __detail::__mod<_UIntType,
00372               __detail::_Shift<_UIntType, __w>::__value>(__sum);
00373 
00374             if (__zero)
00375               {
00376                 if (__i == 0)
00377                   {
00378                     if ((_M_x[0] & __upper_mask) != 0u)
00379                       __zero = false;
00380                   }
00381                 else if (_M_x[__i] != 0u)
00382                   __zero = false;
00383               }
00384           }
00385         if (__zero)
00386           _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
00387         _M_p = state_size;
00388       }
00389 
00390   template<typename _UIntType, size_t __w,
00391            size_t __n, size_t __m, size_t __r,
00392            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00393            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00394            _UIntType __f>
00395     void
00396     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00397                             __s, __b, __t, __c, __l, __f>::
00398     _M_gen_rand(void)
00399     {
00400       const _UIntType __upper_mask = (~_UIntType()) << __r;
00401       const _UIntType __lower_mask = ~__upper_mask;
00402 
00403       for (size_t __k = 0; __k < (__n - __m); ++__k)
00404         {
00405           _UIntType __y = ((_M_x[__k] & __upper_mask)
00406                            | (_M_x[__k + 1] & __lower_mask));
00407           _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
00408                        ^ ((__y & 0x01) ? __a : 0));
00409         }
00410 
00411       for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
00412         {
00413           _UIntType __y = ((_M_x[__k] & __upper_mask)
00414                            | (_M_x[__k + 1] & __lower_mask));
00415           _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
00416                        ^ ((__y & 0x01) ? __a : 0));
00417         }
00418 
00419       _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
00420                        | (_M_x[0] & __lower_mask));
00421       _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
00422                        ^ ((__y & 0x01) ? __a : 0));
00423       _M_p = 0;
00424     }
00425 
00426   template<typename _UIntType, size_t __w,
00427            size_t __n, size_t __m, size_t __r,
00428            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00429            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00430            _UIntType __f>
00431     void
00432     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00433                             __s, __b, __t, __c, __l, __f>::
00434     discard(unsigned long long __z)
00435     {
00436       while (__z > state_size - _M_p)
00437         {
00438           __z -= state_size - _M_p;
00439           _M_gen_rand();
00440         }
00441       _M_p += __z;
00442     }
00443 
00444   template<typename _UIntType, size_t __w,
00445            size_t __n, size_t __m, size_t __r,
00446            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00447            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00448            _UIntType __f>
00449     typename
00450     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00451                             __s, __b, __t, __c, __l, __f>::result_type
00452     mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
00453                             __s, __b, __t, __c, __l, __f>::
00454     operator()()
00455     {
00456       // Reload the vector - cost is O(n) amortized over n calls.
00457       if (_M_p >= state_size)
00458         _M_gen_rand();
00459 
00460       // Calculate o(x(i)).
00461       result_type __z = _M_x[_M_p++];
00462       __z ^= (__z >> __u) & __d;
00463       __z ^= (__z << __s) & __b;
00464       __z ^= (__z << __t) & __c;
00465       __z ^= (__z >> __l);
00466 
00467       return __z;
00468     }
00469 
00470   template<typename _UIntType, size_t __w,
00471            size_t __n, size_t __m, size_t __r,
00472            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00473            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00474            _UIntType __f, typename _CharT, typename _Traits>
00475     std::basic_ostream<_CharT, _Traits>&
00476     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00477                const mersenne_twister_engine<_UIntType, __w, __n, __m,
00478                __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
00479     {
00480       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00481       typedef typename __ostream_type::ios_base    __ios_base;
00482 
00483       const typename __ios_base::fmtflags __flags = __os.flags();
00484       const _CharT __fill = __os.fill();
00485       const _CharT __space = __os.widen(' ');
00486       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00487       __os.fill(__space);
00488 
00489       for (size_t __i = 0; __i < __n; ++__i)
00490         __os << __x._M_x[__i] << __space;
00491       __os << __x._M_p;
00492 
00493       __os.flags(__flags);
00494       __os.fill(__fill);
00495       return __os;
00496     }
00497 
00498   template<typename _UIntType, size_t __w,
00499            size_t __n, size_t __m, size_t __r,
00500            _UIntType __a, size_t __u, _UIntType __d, size_t __s,
00501            _UIntType __b, size_t __t, _UIntType __c, size_t __l,
00502            _UIntType __f, typename _CharT, typename _Traits>
00503     std::basic_istream<_CharT, _Traits>&
00504     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00505                mersenne_twister_engine<_UIntType, __w, __n, __m,
00506                __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
00507     {
00508       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00509       typedef typename __istream_type::ios_base    __ios_base;
00510 
00511       const typename __ios_base::fmtflags __flags = __is.flags();
00512       __is.flags(__ios_base::dec | __ios_base::skipws);
00513 
00514       for (size_t __i = 0; __i < __n; ++__i)
00515         __is >> __x._M_x[__i];
00516       __is >> __x._M_p;
00517 
00518       __is.flags(__flags);
00519       return __is;
00520     }
00521 
00522 
00523   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00524     constexpr size_t
00525     subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
00526 
00527   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00528     constexpr size_t
00529     subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
00530 
00531   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00532     constexpr size_t
00533     subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
00534 
00535   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00536     constexpr _UIntType
00537     subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
00538 
00539   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00540     void
00541     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00542     seed(result_type __value)
00543     {
00544       std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
00545         __lcg(__value == 0u ? default_seed : __value);
00546 
00547       const size_t __n = (__w + 31) / 32;
00548 
00549       for (size_t __i = 0; __i < long_lag; ++__i)
00550         {
00551           _UIntType __sum = 0u;
00552           _UIntType __factor = 1u;
00553           for (size_t __j = 0; __j < __n; ++__j)
00554             {
00555               __sum += __detail::__mod<uint_least32_t,
00556                        __detail::_Shift<uint_least32_t, 32>::__value>
00557                          (__lcg()) * __factor;
00558               __factor *= __detail::_Shift<_UIntType, 32>::__value;
00559             }
00560           _M_x[__i] = __detail::__mod<_UIntType,
00561             __detail::_Shift<_UIntType, __w>::__value>(__sum);
00562         }
00563       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
00564       _M_p = 0;
00565     }
00566 
00567   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
00568     template<typename _Sseq>
00569       auto
00570       subtract_with_carry_engine<_UIntType, __w, __s, __r>::
00571       seed(_Sseq& __q)
00572       -> _If_seed_seq<_Sseq>
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<double>(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             // 1 / 78
01308             const double __178 = 0.0128205128205128205128205128205128L;
01309             // e^(1 / 78)
01310             const double __e178 = 1.0129030479320018583185514777512983L;
01311             const double __c5 = __c4 + __e178;
01312             const double __c = __param._M_cb + __c5;
01313             const double __2cx = 2 * (2 * __m + __param._M_d);
01314 
01315             bool __reject = true;
01316             do
01317               {
01318                 const double __u = __c * __aurng();
01319                 const double __e = -std::log(1.0 - __aurng());
01320 
01321                 double __w = 0.0;
01322 
01323                 if (__u <= __c1)
01324                   {
01325                     const double __n = _M_nd(__urng);
01326                     const double __y = -std::abs(__n) * __param._M_sm - 1;
01327                     __x = std::floor(__y);
01328                     __w = -__n * __n / 2;
01329                     if (__x < -__m)
01330                       continue;
01331                   }
01332                 else if (__u <= __c2)
01333                   {
01334                     const double __n = _M_nd(__urng);
01335                     const double __y = 1 + std::abs(__n) * __param._M_scx;
01336                     __x = std::ceil(__y);
01337                     __w = __y * (2 - __y) * __param._M_1cx;
01338                     if (__x > __param._M_d)
01339                       continue;
01340                   }
01341                 else if (__u <= __c3)
01342                   // NB: This case not in the book, nor in the Errata,
01343                   // but should be ok...
01344                   __x = -1;
01345                 else if (__u <= __c4)
01346                   __x = 0;
01347                 else if (__u <= __c5)
01348                   {
01349                     __x = 1;
01350                     // Only in the Errata, see libstdc++/83237.
01351                     __w = __178;
01352                   }
01353                 else
01354                   {
01355                     const double __v = -std::log(1.0 - __aurng());
01356                     const double __y = __param._M_d
01357                                      + __v * __2cx / __param._M_d;
01358                     __x = std::ceil(__y);
01359                     __w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
01360                   }
01361 
01362                 __reject = (__w - __e - __x * __param._M_lm_thr
01363                             > __param._M_lfm - std::lgamma(__x + __m + 1));
01364 
01365                 __reject |= __x + __m >= __thr;
01366 
01367               } while (__reject);
01368 
01369             return result_type(__x + __m + __naf);
01370           }
01371         else
01372 #endif
01373           {
01374             _IntType     __x = 0;
01375             double __prod = 1.0;
01376 
01377             do
01378               {
01379                 __prod *= __aurng();
01380                 __x += 1;
01381               }
01382             while (__prod > __param._M_lm_thr);
01383 
01384             return __x - 1;
01385           }
01386       }
01387 
01388   template<typename _IntType>
01389     template<typename _ForwardIterator,
01390              typename _UniformRandomNumberGenerator>
01391       void
01392       poisson_distribution<_IntType>::
01393       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01394                       _UniformRandomNumberGenerator& __urng,
01395                       const param_type& __param)
01396       {
01397         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01398         // We could duplicate everything from operator()...
01399         while (__f != __t)
01400           *__f++ = this->operator()(__urng, __param);
01401       }
01402 
01403   template<typename _IntType,
01404            typename _CharT, typename _Traits>
01405     std::basic_ostream<_CharT, _Traits>&
01406     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01407                const poisson_distribution<_IntType>& __x)
01408     {
01409       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01410       typedef typename __ostream_type::ios_base    __ios_base;
01411 
01412       const typename __ios_base::fmtflags __flags = __os.flags();
01413       const _CharT __fill = __os.fill();
01414       const std::streamsize __precision = __os.precision();
01415       const _CharT __space = __os.widen(' ');
01416       __os.flags(__ios_base::scientific | __ios_base::left);
01417       __os.fill(__space);
01418       __os.precision(std::numeric_limits<double>::max_digits10);
01419 
01420       __os << __x.mean() << __space << __x._M_nd;
01421 
01422       __os.flags(__flags);
01423       __os.fill(__fill);
01424       __os.precision(__precision);
01425       return __os;
01426     }
01427 
01428   template<typename _IntType,
01429            typename _CharT, typename _Traits>
01430     std::basic_istream<_CharT, _Traits>&
01431     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01432                poisson_distribution<_IntType>& __x)
01433     {
01434       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01435       typedef typename __istream_type::ios_base    __ios_base;
01436 
01437       const typename __ios_base::fmtflags __flags = __is.flags();
01438       __is.flags(__ios_base::skipws);
01439 
01440       double __mean;
01441       __is >> __mean >> __x._M_nd;
01442       __x.param(typename poisson_distribution<_IntType>::param_type(__mean));
01443 
01444       __is.flags(__flags);
01445       return __is;
01446     }
01447 
01448 
01449   template<typename _IntType>
01450     void
01451     binomial_distribution<_IntType>::param_type::
01452     _M_initialize()
01453     {
01454       const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
01455 
01456       _M_easy = true;
01457 
01458 #if _GLIBCXX_USE_C99_MATH_TR1
01459       if (_M_t * __p12 >= 8)
01460         {
01461           _M_easy = false;
01462           const double __np = std::floor(_M_t * __p12);
01463           const double __pa = __np / _M_t;
01464           const double __1p = 1 - __pa;
01465 
01466           const double __pi_4 = 0.7853981633974483096156608458198757L;
01467           const double __d1x =
01468             std::sqrt(__np * __1p * std::log(32 * __np
01469                                              / (81 * __pi_4 * __1p)));
01470           _M_d1 = std::round(std::max<double>(1.0, __d1x));
01471           const double __d2x =
01472             std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
01473                                              / (__pi_4 * __pa)));
01474           _M_d2 = std::round(std::max<double>(1.0, __d2x));
01475 
01476           // sqrt(pi / 2)
01477           const double __spi_2 = 1.2533141373155002512078826424055226L;
01478           _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
01479           _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
01480           _M_c = 2 * _M_d1 / __np;
01481           _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
01482           const double __a12 = _M_a1 + _M_s2 * __spi_2;
01483           const double __s1s = _M_s1 * _M_s1;
01484           _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
01485                              * 2 * __s1s / _M_d1
01486                              * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
01487           const double __s2s = _M_s2 * _M_s2;
01488           _M_s = (_M_a123 + 2 * __s2s / _M_d2
01489                   * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
01490           _M_lf = (std::lgamma(__np + 1)
01491                    + std::lgamma(_M_t - __np + 1));
01492           _M_lp1p = std::log(__pa / __1p);
01493 
01494           _M_q = -std::log(1 - (__p12 - __pa) / __1p);
01495         }
01496       else
01497 #endif
01498         _M_q = -std::log(1 - __p12);
01499     }
01500 
01501   template<typename _IntType>
01502     template<typename _UniformRandomNumberGenerator>
01503       typename binomial_distribution<_IntType>::result_type
01504       binomial_distribution<_IntType>::
01505       _M_waiting(_UniformRandomNumberGenerator& __urng,
01506                  _IntType __t, double __q)
01507       {
01508         _IntType __x = 0;
01509         double __sum = 0.0;
01510         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01511           __aurng(__urng);
01512 
01513         do
01514           {
01515             if (__t == __x)
01516               return __x;
01517             const double __e = -std::log(1.0 - __aurng());
01518             __sum += __e / (__t - __x);
01519             __x += 1;
01520           }
01521         while (__sum <= __q);
01522 
01523         return __x - 1;
01524       }
01525 
01526   /**
01527    * A rejection algorithm when t * p >= 8 and a simple waiting time
01528    * method - the second in the referenced book - otherwise.
01529    * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
01530    * is defined.
01531    *
01532    * Reference:
01533    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01534    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
01535    */
01536   template<typename _IntType>
01537     template<typename _UniformRandomNumberGenerator>
01538       typename binomial_distribution<_IntType>::result_type
01539       binomial_distribution<_IntType>::
01540       operator()(_UniformRandomNumberGenerator& __urng,
01541                  const param_type& __param)
01542       {
01543         result_type __ret;
01544         const _IntType __t = __param.t();
01545         const double __p = __param.p();
01546         const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
01547         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
01548           __aurng(__urng);
01549 
01550 #if _GLIBCXX_USE_C99_MATH_TR1
01551         if (!__param._M_easy)
01552           {
01553             double __x;
01554 
01555             // See comments above...
01556             const double __naf =
01557               (1 - std::numeric_limits<double>::epsilon()) / 2;
01558             const double __thr =
01559               std::numeric_limits<_IntType>::max() + __naf;
01560 
01561             const double __np = std::floor(__t * __p12);
01562 
01563             // sqrt(pi / 2)
01564             const double __spi_2 = 1.2533141373155002512078826424055226L;
01565             const double __a1 = __param._M_a1;
01566             const double __a12 = __a1 + __param._M_s2 * __spi_2;
01567             const double __a123 = __param._M_a123;
01568             const double __s1s = __param._M_s1 * __param._M_s1;
01569             const double __s2s = __param._M_s2 * __param._M_s2;
01570 
01571             bool __reject;
01572             do
01573               {
01574                 const double __u = __param._M_s * __aurng();
01575 
01576                 double __v;
01577 
01578                 if (__u <= __a1)
01579                   {
01580                     const double __n = _M_nd(__urng);
01581                     const double __y = __param._M_s1 * std::abs(__n);
01582                     __reject = __y >= __param._M_d1;
01583                     if (!__reject)
01584                       {
01585                         const double __e = -std::log(1.0 - __aurng());
01586                         __x = std::floor(__y);
01587                         __v = -__e - __n * __n / 2 + __param._M_c;
01588                       }
01589                   }
01590                 else if (__u <= __a12)
01591                   {
01592                     const double __n = _M_nd(__urng);
01593                     const double __y = __param._M_s2 * std::abs(__n);
01594                     __reject = __y >= __param._M_d2;
01595                     if (!__reject)
01596                       {
01597                         const double __e = -std::log(1.0 - __aurng());
01598                         __x = std::floor(-__y);
01599                         __v = -__e - __n * __n / 2;
01600                       }
01601                   }
01602                 else if (__u <= __a123)
01603                   {
01604                     const double __e1 = -std::log(1.0 - __aurng());
01605                     const double __e2 = -std::log(1.0 - __aurng());
01606 
01607                     const double __y = __param._M_d1
01608                                      + 2 * __s1s * __e1 / __param._M_d1;
01609                     __x = std::floor(__y);
01610                     __v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
01611                                                     -__y / (2 * __s1s)));
01612                     __reject = false;
01613                   }
01614                 else
01615                   {
01616                     const double __e1 = -std::log(1.0 - __aurng());
01617                     const double __e2 = -std::log(1.0 - __aurng());
01618 
01619                     const double __y = __param._M_d2
01620                                      + 2 * __s2s * __e1 / __param._M_d2;
01621                     __x = std::floor(-__y);
01622                     __v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
01623                     __reject = false;
01624                   }
01625 
01626                 __reject = __reject || __x < -__np || __x > __t - __np;
01627                 if (!__reject)
01628                   {
01629                     const double __lfx =
01630                       std::lgamma(__np + __x + 1)
01631                       + std::lgamma(__t - (__np + __x) + 1);
01632                     __reject = __v > __param._M_lf - __lfx
01633                              + __x * __param._M_lp1p;
01634                   }
01635 
01636                 __reject |= __x + __np >= __thr;
01637               }
01638             while (__reject);
01639 
01640             __x += __np + __naf;
01641 
01642             const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
01643                                             __param._M_q);
01644             __ret = _IntType(__x) + __z;
01645           }
01646         else
01647 #endif
01648           __ret = _M_waiting(__urng, __t, __param._M_q);
01649 
01650         if (__p12 != __p)
01651           __ret = __t - __ret;
01652         return __ret;
01653       }
01654 
01655   template<typename _IntType>
01656     template<typename _ForwardIterator,
01657              typename _UniformRandomNumberGenerator>
01658       void
01659       binomial_distribution<_IntType>::
01660       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01661                       _UniformRandomNumberGenerator& __urng,
01662                       const param_type& __param)
01663       {
01664         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01665         // We could duplicate everything from operator()...
01666         while (__f != __t)
01667           *__f++ = this->operator()(__urng, __param);
01668       }
01669 
01670   template<typename _IntType,
01671            typename _CharT, typename _Traits>
01672     std::basic_ostream<_CharT, _Traits>&
01673     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01674                const binomial_distribution<_IntType>& __x)
01675     {
01676       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01677       typedef typename __ostream_type::ios_base    __ios_base;
01678 
01679       const typename __ios_base::fmtflags __flags = __os.flags();
01680       const _CharT __fill = __os.fill();
01681       const std::streamsize __precision = __os.precision();
01682       const _CharT __space = __os.widen(' ');
01683       __os.flags(__ios_base::scientific | __ios_base::left);
01684       __os.fill(__space);
01685       __os.precision(std::numeric_limits<double>::max_digits10);
01686 
01687       __os << __x.t() << __space << __x.p()
01688            << __space << __x._M_nd;
01689 
01690       __os.flags(__flags);
01691       __os.fill(__fill);
01692       __os.precision(__precision);
01693       return __os;
01694     }
01695 
01696   template<typename _IntType,
01697            typename _CharT, typename _Traits>
01698     std::basic_istream<_CharT, _Traits>&
01699     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01700                binomial_distribution<_IntType>& __x)
01701     {
01702       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01703       typedef typename __istream_type::ios_base    __ios_base;
01704 
01705       const typename __ios_base::fmtflags __flags = __is.flags();
01706       __is.flags(__ios_base::dec | __ios_base::skipws);
01707 
01708       _IntType __t;
01709       double __p;
01710       __is >> __t >> __p >> __x._M_nd;
01711       __x.param(typename binomial_distribution<_IntType>::
01712                 param_type(__t, __p));
01713 
01714       __is.flags(__flags);
01715       return __is;
01716     }
01717 
01718 
01719   template<typename _RealType>
01720     template<typename _ForwardIterator,
01721              typename _UniformRandomNumberGenerator>
01722       void
01723       std::exponential_distribution<_RealType>::
01724       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01725                       _UniformRandomNumberGenerator& __urng,
01726                       const param_type& __p)
01727       {
01728         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01729         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01730           __aurng(__urng);
01731         while (__f != __t)
01732           *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
01733       }
01734 
01735   template<typename _RealType, typename _CharT, typename _Traits>
01736     std::basic_ostream<_CharT, _Traits>&
01737     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01738                const exponential_distribution<_RealType>& __x)
01739     {
01740       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01741       typedef typename __ostream_type::ios_base    __ios_base;
01742 
01743       const typename __ios_base::fmtflags __flags = __os.flags();
01744       const _CharT __fill = __os.fill();
01745       const std::streamsize __precision = __os.precision();
01746       __os.flags(__ios_base::scientific | __ios_base::left);
01747       __os.fill(__os.widen(' '));
01748       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01749 
01750       __os << __x.lambda();
01751 
01752       __os.flags(__flags);
01753       __os.fill(__fill);
01754       __os.precision(__precision);
01755       return __os;
01756     }
01757 
01758   template<typename _RealType, typename _CharT, typename _Traits>
01759     std::basic_istream<_CharT, _Traits>&
01760     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01761                exponential_distribution<_RealType>& __x)
01762     {
01763       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01764       typedef typename __istream_type::ios_base    __ios_base;
01765 
01766       const typename __ios_base::fmtflags __flags = __is.flags();
01767       __is.flags(__ios_base::dec | __ios_base::skipws);
01768 
01769       _RealType __lambda;
01770       __is >> __lambda;
01771       __x.param(typename exponential_distribution<_RealType>::
01772                 param_type(__lambda));
01773 
01774       __is.flags(__flags);
01775       return __is;
01776     }
01777 
01778 
01779   /**
01780    * Polar method due to Marsaglia.
01781    *
01782    * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
01783    * New York, 1986, Ch. V, Sect. 4.4.
01784    */
01785   template<typename _RealType>
01786     template<typename _UniformRandomNumberGenerator>
01787       typename normal_distribution<_RealType>::result_type
01788       normal_distribution<_RealType>::
01789       operator()(_UniformRandomNumberGenerator& __urng,
01790                  const param_type& __param)
01791       {
01792         result_type __ret;
01793         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01794           __aurng(__urng);
01795 
01796         if (_M_saved_available)
01797           {
01798             _M_saved_available = false;
01799             __ret = _M_saved;
01800           }
01801         else
01802           {
01803             result_type __x, __y, __r2;
01804             do
01805               {
01806                 __x = result_type(2.0) * __aurng() - 1.0;
01807                 __y = result_type(2.0) * __aurng() - 1.0;
01808                 __r2 = __x * __x + __y * __y;
01809               }
01810             while (__r2 > 1.0 || __r2 == 0.0);
01811 
01812             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
01813             _M_saved = __x * __mult;
01814             _M_saved_available = true;
01815             __ret = __y * __mult;
01816           }
01817 
01818         __ret = __ret * __param.stddev() + __param.mean();
01819         return __ret;
01820       }
01821 
01822   template<typename _RealType>
01823     template<typename _ForwardIterator,
01824              typename _UniformRandomNumberGenerator>
01825       void
01826       normal_distribution<_RealType>::
01827       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01828                       _UniformRandomNumberGenerator& __urng,
01829                       const param_type& __param)
01830       {
01831         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01832 
01833         if (__f == __t)
01834           return;
01835 
01836         if (_M_saved_available)
01837           {
01838             _M_saved_available = false;
01839             *__f++ = _M_saved * __param.stddev() + __param.mean();
01840 
01841             if (__f == __t)
01842               return;
01843           }
01844 
01845         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01846           __aurng(__urng);
01847 
01848         while (__f + 1 < __t)
01849           {
01850             result_type __x, __y, __r2;
01851             do
01852               {
01853                 __x = result_type(2.0) * __aurng() - 1.0;
01854                 __y = result_type(2.0) * __aurng() - 1.0;
01855                 __r2 = __x * __x + __y * __y;
01856               }
01857             while (__r2 > 1.0 || __r2 == 0.0);
01858 
01859             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
01860             *__f++ = __y * __mult * __param.stddev() + __param.mean();
01861             *__f++ = __x * __mult * __param.stddev() + __param.mean();
01862           }
01863 
01864         if (__f != __t)
01865           {
01866             result_type __x, __y, __r2;
01867             do
01868               {
01869                 __x = result_type(2.0) * __aurng() - 1.0;
01870                 __y = result_type(2.0) * __aurng() - 1.0;
01871                 __r2 = __x * __x + __y * __y;
01872               }
01873             while (__r2 > 1.0 || __r2 == 0.0);
01874 
01875             const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
01876             _M_saved = __x * __mult;
01877             _M_saved_available = true;
01878             *__f = __y * __mult * __param.stddev() + __param.mean();
01879           }
01880       }
01881 
01882   template<typename _RealType>
01883     bool
01884     operator==(const std::normal_distribution<_RealType>& __d1,
01885                const std::normal_distribution<_RealType>& __d2)
01886     {
01887       if (__d1._M_param == __d2._M_param
01888           && __d1._M_saved_available == __d2._M_saved_available)
01889         {
01890           if (__d1._M_saved_available
01891               && __d1._M_saved == __d2._M_saved)
01892             return true;
01893           else if(!__d1._M_saved_available)
01894             return true;
01895           else
01896             return false;
01897         }
01898       else
01899         return false;
01900     }
01901 
01902   template<typename _RealType, typename _CharT, typename _Traits>
01903     std::basic_ostream<_CharT, _Traits>&
01904     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01905                const normal_distribution<_RealType>& __x)
01906     {
01907       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01908       typedef typename __ostream_type::ios_base    __ios_base;
01909 
01910       const typename __ios_base::fmtflags __flags = __os.flags();
01911       const _CharT __fill = __os.fill();
01912       const std::streamsize __precision = __os.precision();
01913       const _CharT __space = __os.widen(' ');
01914       __os.flags(__ios_base::scientific | __ios_base::left);
01915       __os.fill(__space);
01916       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01917 
01918       __os << __x.mean() << __space << __x.stddev()
01919            << __space << __x._M_saved_available;
01920       if (__x._M_saved_available)
01921         __os << __space << __x._M_saved;
01922 
01923       __os.flags(__flags);
01924       __os.fill(__fill);
01925       __os.precision(__precision);
01926       return __os;
01927     }
01928 
01929   template<typename _RealType, typename _CharT, typename _Traits>
01930     std::basic_istream<_CharT, _Traits>&
01931     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01932                normal_distribution<_RealType>& __x)
01933     {
01934       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01935       typedef typename __istream_type::ios_base    __ios_base;
01936 
01937       const typename __ios_base::fmtflags __flags = __is.flags();
01938       __is.flags(__ios_base::dec | __ios_base::skipws);
01939 
01940       double __mean, __stddev;
01941       __is >> __mean >> __stddev
01942            >> __x._M_saved_available;
01943       if (__x._M_saved_available)
01944         __is >> __x._M_saved;
01945       __x.param(typename normal_distribution<_RealType>::
01946                 param_type(__mean, __stddev));
01947 
01948       __is.flags(__flags);
01949       return __is;
01950     }
01951 
01952 
01953   template<typename _RealType>
01954     template<typename _ForwardIterator,
01955              typename _UniformRandomNumberGenerator>
01956       void
01957       lognormal_distribution<_RealType>::
01958       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
01959                       _UniformRandomNumberGenerator& __urng,
01960                       const param_type& __p)
01961       {
01962         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
01963           while (__f != __t)
01964             *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
01965       }
01966 
01967   template<typename _RealType, typename _CharT, typename _Traits>
01968     std::basic_ostream<_CharT, _Traits>&
01969     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01970                const lognormal_distribution<_RealType>& __x)
01971     {
01972       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01973       typedef typename __ostream_type::ios_base    __ios_base;
01974 
01975       const typename __ios_base::fmtflags __flags = __os.flags();
01976       const _CharT __fill = __os.fill();
01977       const std::streamsize __precision = __os.precision();
01978       const _CharT __space = __os.widen(' ');
01979       __os.flags(__ios_base::scientific | __ios_base::left);
01980       __os.fill(__space);
01981       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01982 
01983       __os << __x.m() << __space << __x.s()
01984            << __space << __x._M_nd;
01985 
01986       __os.flags(__flags);
01987       __os.fill(__fill);
01988       __os.precision(__precision);
01989       return __os;
01990     }
01991 
01992   template<typename _RealType, typename _CharT, typename _Traits>
01993     std::basic_istream<_CharT, _Traits>&
01994     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01995                lognormal_distribution<_RealType>& __x)
01996     {
01997       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01998       typedef typename __istream_type::ios_base    __ios_base;
01999 
02000       const typename __ios_base::fmtflags __flags = __is.flags();
02001       __is.flags(__ios_base::dec | __ios_base::skipws);
02002 
02003       _RealType __m, __s;
02004       __is >> __m >> __s >> __x._M_nd;
02005       __x.param(typename lognormal_distribution<_RealType>::
02006                 param_type(__m, __s));
02007 
02008       __is.flags(__flags);
02009       return __is;
02010     }
02011 
02012   template<typename _RealType>
02013     template<typename _ForwardIterator,
02014              typename _UniformRandomNumberGenerator>
02015       void
02016       std::chi_squared_distribution<_RealType>::
02017       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02018                       _UniformRandomNumberGenerator& __urng)
02019       {
02020         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02021         while (__f != __t)
02022           *__f++ = 2 * _M_gd(__urng);
02023       }
02024 
02025   template<typename _RealType>
02026     template<typename _ForwardIterator,
02027              typename _UniformRandomNumberGenerator>
02028       void
02029       std::chi_squared_distribution<_RealType>::
02030       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02031                       _UniformRandomNumberGenerator& __urng,
02032                       const typename
02033                       std::gamma_distribution<result_type>::param_type& __p)
02034       {
02035         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02036         while (__f != __t)
02037           *__f++ = 2 * _M_gd(__urng, __p);
02038       }
02039 
02040   template<typename _RealType, typename _CharT, typename _Traits>
02041     std::basic_ostream<_CharT, _Traits>&
02042     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02043                const chi_squared_distribution<_RealType>& __x)
02044     {
02045       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02046       typedef typename __ostream_type::ios_base    __ios_base;
02047 
02048       const typename __ios_base::fmtflags __flags = __os.flags();
02049       const _CharT __fill = __os.fill();
02050       const std::streamsize __precision = __os.precision();
02051       const _CharT __space = __os.widen(' ');
02052       __os.flags(__ios_base::scientific | __ios_base::left);
02053       __os.fill(__space);
02054       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02055 
02056       __os << __x.n() << __space << __x._M_gd;
02057 
02058       __os.flags(__flags);
02059       __os.fill(__fill);
02060       __os.precision(__precision);
02061       return __os;
02062     }
02063 
02064   template<typename _RealType, typename _CharT, typename _Traits>
02065     std::basic_istream<_CharT, _Traits>&
02066     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02067                chi_squared_distribution<_RealType>& __x)
02068     {
02069       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02070       typedef typename __istream_type::ios_base    __ios_base;
02071 
02072       const typename __ios_base::fmtflags __flags = __is.flags();
02073       __is.flags(__ios_base::dec | __ios_base::skipws);
02074 
02075       _RealType __n;
02076       __is >> __n >> __x._M_gd;
02077       __x.param(typename chi_squared_distribution<_RealType>::
02078                 param_type(__n));
02079 
02080       __is.flags(__flags);
02081       return __is;
02082     }
02083 
02084 
02085   template<typename _RealType>
02086     template<typename _UniformRandomNumberGenerator>
02087       typename cauchy_distribution<_RealType>::result_type
02088       cauchy_distribution<_RealType>::
02089       operator()(_UniformRandomNumberGenerator& __urng,
02090                  const param_type& __p)
02091       {
02092         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02093           __aurng(__urng);
02094         _RealType __u;
02095         do
02096           __u = __aurng();
02097         while (__u == 0.5);
02098 
02099         const _RealType __pi = 3.1415926535897932384626433832795029L;
02100         return __p.a() + __p.b() * std::tan(__pi * __u);
02101       }
02102 
02103   template<typename _RealType>
02104     template<typename _ForwardIterator,
02105              typename _UniformRandomNumberGenerator>
02106       void
02107       cauchy_distribution<_RealType>::
02108       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02109                       _UniformRandomNumberGenerator& __urng,
02110                       const param_type& __p)
02111       {
02112         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02113         const _RealType __pi = 3.1415926535897932384626433832795029L;
02114         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02115           __aurng(__urng);
02116         while (__f != __t)
02117           {
02118             _RealType __u;
02119             do
02120               __u = __aurng();
02121             while (__u == 0.5);
02122 
02123             *__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
02124           }
02125       }
02126 
02127   template<typename _RealType, typename _CharT, typename _Traits>
02128     std::basic_ostream<_CharT, _Traits>&
02129     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02130                const cauchy_distribution<_RealType>& __x)
02131     {
02132       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02133       typedef typename __ostream_type::ios_base    __ios_base;
02134 
02135       const typename __ios_base::fmtflags __flags = __os.flags();
02136       const _CharT __fill = __os.fill();
02137       const std::streamsize __precision = __os.precision();
02138       const _CharT __space = __os.widen(' ');
02139       __os.flags(__ios_base::scientific | __ios_base::left);
02140       __os.fill(__space);
02141       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02142 
02143       __os << __x.a() << __space << __x.b();
02144 
02145       __os.flags(__flags);
02146       __os.fill(__fill);
02147       __os.precision(__precision);
02148       return __os;
02149     }
02150 
02151   template<typename _RealType, typename _CharT, typename _Traits>
02152     std::basic_istream<_CharT, _Traits>&
02153     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02154                cauchy_distribution<_RealType>& __x)
02155     {
02156       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02157       typedef typename __istream_type::ios_base    __ios_base;
02158 
02159       const typename __ios_base::fmtflags __flags = __is.flags();
02160       __is.flags(__ios_base::dec | __ios_base::skipws);
02161 
02162       _RealType __a, __b;
02163       __is >> __a >> __b;
02164       __x.param(typename cauchy_distribution<_RealType>::
02165                 param_type(__a, __b));
02166 
02167       __is.flags(__flags);
02168       return __is;
02169     }
02170 
02171 
02172   template<typename _RealType>
02173     template<typename _ForwardIterator,
02174              typename _UniformRandomNumberGenerator>
02175       void
02176       std::fisher_f_distribution<_RealType>::
02177       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02178                       _UniformRandomNumberGenerator& __urng)
02179       {
02180         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02181         while (__f != __t)
02182           *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
02183       }
02184 
02185   template<typename _RealType>
02186     template<typename _ForwardIterator,
02187              typename _UniformRandomNumberGenerator>
02188       void
02189       std::fisher_f_distribution<_RealType>::
02190       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02191                       _UniformRandomNumberGenerator& __urng,
02192                       const param_type& __p)
02193       {
02194         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02195         typedef typename std::gamma_distribution<result_type>::param_type
02196           param_type;
02197         param_type __p1(__p.m() / 2);
02198         param_type __p2(__p.n() / 2);
02199         while (__f != __t)
02200           *__f++ = ((_M_gd_x(__urng, __p1) * n())
02201                     / (_M_gd_y(__urng, __p2) * m()));
02202       }
02203 
02204   template<typename _RealType, typename _CharT, typename _Traits>
02205     std::basic_ostream<_CharT, _Traits>&
02206     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02207                const fisher_f_distribution<_RealType>& __x)
02208     {
02209       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02210       typedef typename __ostream_type::ios_base    __ios_base;
02211 
02212       const typename __ios_base::fmtflags __flags = __os.flags();
02213       const _CharT __fill = __os.fill();
02214       const std::streamsize __precision = __os.precision();
02215       const _CharT __space = __os.widen(' ');
02216       __os.flags(__ios_base::scientific | __ios_base::left);
02217       __os.fill(__space);
02218       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02219 
02220       __os << __x.m() << __space << __x.n()
02221            << __space << __x._M_gd_x << __space << __x._M_gd_y;
02222 
02223       __os.flags(__flags);
02224       __os.fill(__fill);
02225       __os.precision(__precision);
02226       return __os;
02227     }
02228 
02229   template<typename _RealType, typename _CharT, typename _Traits>
02230     std::basic_istream<_CharT, _Traits>&
02231     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02232                fisher_f_distribution<_RealType>& __x)
02233     {
02234       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02235       typedef typename __istream_type::ios_base    __ios_base;
02236 
02237       const typename __ios_base::fmtflags __flags = __is.flags();
02238       __is.flags(__ios_base::dec | __ios_base::skipws);
02239 
02240       _RealType __m, __n;
02241       __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y;
02242       __x.param(typename fisher_f_distribution<_RealType>::
02243                 param_type(__m, __n));
02244 
02245       __is.flags(__flags);
02246       return __is;
02247     }
02248 
02249 
02250   template<typename _RealType>
02251     template<typename _ForwardIterator,
02252              typename _UniformRandomNumberGenerator>
02253       void
02254       std::student_t_distribution<_RealType>::
02255       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02256                       _UniformRandomNumberGenerator& __urng)
02257       {
02258         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02259         while (__f != __t)
02260           *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
02261       }
02262 
02263   template<typename _RealType>
02264     template<typename _ForwardIterator,
02265              typename _UniformRandomNumberGenerator>
02266       void
02267       std::student_t_distribution<_RealType>::
02268       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02269                       _UniformRandomNumberGenerator& __urng,
02270                       const param_type& __p)
02271       {
02272         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02273         typename std::gamma_distribution<result_type>::param_type
02274           __p2(__p.n() / 2, 2);
02275         while (__f != __t)
02276           *__f++ =  _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
02277       }
02278 
02279   template<typename _RealType, typename _CharT, typename _Traits>
02280     std::basic_ostream<_CharT, _Traits>&
02281     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02282                const student_t_distribution<_RealType>& __x)
02283     {
02284       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02285       typedef typename __ostream_type::ios_base    __ios_base;
02286 
02287       const typename __ios_base::fmtflags __flags = __os.flags();
02288       const _CharT __fill = __os.fill();
02289       const std::streamsize __precision = __os.precision();
02290       const _CharT __space = __os.widen(' ');
02291       __os.flags(__ios_base::scientific | __ios_base::left);
02292       __os.fill(__space);
02293       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02294 
02295       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
02296 
02297       __os.flags(__flags);
02298       __os.fill(__fill);
02299       __os.precision(__precision);
02300       return __os;
02301     }
02302 
02303   template<typename _RealType, typename _CharT, typename _Traits>
02304     std::basic_istream<_CharT, _Traits>&
02305     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02306                student_t_distribution<_RealType>& __x)
02307     {
02308       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02309       typedef typename __istream_type::ios_base    __ios_base;
02310 
02311       const typename __ios_base::fmtflags __flags = __is.flags();
02312       __is.flags(__ios_base::dec | __ios_base::skipws);
02313 
02314       _RealType __n;
02315       __is >> __n >> __x._M_nd >> __x._M_gd;
02316       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
02317 
02318       __is.flags(__flags);
02319       return __is;
02320     }
02321 
02322 
02323   template<typename _RealType>
02324     void
02325     gamma_distribution<_RealType>::param_type::
02326     _M_initialize()
02327     {
02328       _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
02329 
02330       const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
02331       _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
02332     }
02333 
02334   /**
02335    * Marsaglia, G. and Tsang, W. W.
02336    * "A Simple Method for Generating Gamma Variables"
02337    * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
02338    */
02339   template<typename _RealType>
02340     template<typename _UniformRandomNumberGenerator>
02341       typename gamma_distribution<_RealType>::result_type
02342       gamma_distribution<_RealType>::
02343       operator()(_UniformRandomNumberGenerator& __urng,
02344                  const param_type& __param)
02345       {
02346         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02347           __aurng(__urng);
02348 
02349         result_type __u, __v, __n;
02350         const result_type __a1 = (__param._M_malpha
02351                                   - _RealType(1.0) / _RealType(3.0));
02352 
02353         do
02354           {
02355             do
02356               {
02357                 __n = _M_nd(__urng);
02358                 __v = result_type(1.0) + __param._M_a2 * __n; 
02359               }
02360             while (__v <= 0.0);
02361 
02362             __v = __v * __v * __v;
02363             __u = __aurng();
02364           }
02365         while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
02366                && (std::log(__u) > (0.5 * __n * __n + __a1
02367                                     * (1.0 - __v + std::log(__v)))));
02368 
02369         if (__param.alpha() == __param._M_malpha)
02370           return __a1 * __v * __param.beta();
02371         else
02372           {
02373             do
02374               __u = __aurng();
02375             while (__u == 0.0);
02376             
02377             return (std::pow(__u, result_type(1.0) / __param.alpha())
02378                     * __a1 * __v * __param.beta());
02379           }
02380       }
02381 
02382   template<typename _RealType>
02383     template<typename _ForwardIterator,
02384              typename _UniformRandomNumberGenerator>
02385       void
02386       gamma_distribution<_RealType>::
02387       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02388                       _UniformRandomNumberGenerator& __urng,
02389                       const param_type& __param)
02390       {
02391         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02392         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02393           __aurng(__urng);
02394 
02395         result_type __u, __v, __n;
02396         const result_type __a1 = (__param._M_malpha
02397                                   - _RealType(1.0) / _RealType(3.0));
02398 
02399         if (__param.alpha() == __param._M_malpha)
02400           while (__f != __t)
02401             {
02402               do
02403                 {
02404                   do
02405                     {
02406                       __n = _M_nd(__urng);
02407                       __v = result_type(1.0) + __param._M_a2 * __n;
02408                     }
02409                   while (__v <= 0.0);
02410 
02411                   __v = __v * __v * __v;
02412                   __u = __aurng();
02413                 }
02414               while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
02415                      && (std::log(__u) > (0.5 * __n * __n + __a1
02416                                           * (1.0 - __v + std::log(__v)))));
02417 
02418               *__f++ = __a1 * __v * __param.beta();
02419             }
02420         else
02421           while (__f != __t)
02422             {
02423               do
02424                 {
02425                   do
02426                     {
02427                       __n = _M_nd(__urng);
02428                       __v = result_type(1.0) + __param._M_a2 * __n;
02429                     }
02430                   while (__v <= 0.0);
02431 
02432                   __v = __v * __v * __v;
02433                   __u = __aurng();
02434                 }
02435               while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
02436                      && (std::log(__u) > (0.5 * __n * __n + __a1
02437                                           * (1.0 - __v + std::log(__v)))));
02438 
02439               do
02440                 __u = __aurng();
02441               while (__u == 0.0);
02442 
02443               *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
02444                         * __a1 * __v * __param.beta());
02445             }
02446       }
02447 
02448   template<typename _RealType, typename _CharT, typename _Traits>
02449     std::basic_ostream<_CharT, _Traits>&
02450     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02451                const gamma_distribution<_RealType>& __x)
02452     {
02453       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02454       typedef typename __ostream_type::ios_base    __ios_base;
02455 
02456       const typename __ios_base::fmtflags __flags = __os.flags();
02457       const _CharT __fill = __os.fill();
02458       const std::streamsize __precision = __os.precision();
02459       const _CharT __space = __os.widen(' ');
02460       __os.flags(__ios_base::scientific | __ios_base::left);
02461       __os.fill(__space);
02462       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02463 
02464       __os << __x.alpha() << __space << __x.beta()
02465            << __space << __x._M_nd;
02466 
02467       __os.flags(__flags);
02468       __os.fill(__fill);
02469       __os.precision(__precision);
02470       return __os;
02471     }
02472 
02473   template<typename _RealType, typename _CharT, typename _Traits>
02474     std::basic_istream<_CharT, _Traits>&
02475     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02476                gamma_distribution<_RealType>& __x)
02477     {
02478       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02479       typedef typename __istream_type::ios_base    __ios_base;
02480 
02481       const typename __ios_base::fmtflags __flags = __is.flags();
02482       __is.flags(__ios_base::dec | __ios_base::skipws);
02483 
02484       _RealType __alpha_val, __beta_val;
02485       __is >> __alpha_val >> __beta_val >> __x._M_nd;
02486       __x.param(typename gamma_distribution<_RealType>::
02487                 param_type(__alpha_val, __beta_val));
02488 
02489       __is.flags(__flags);
02490       return __is;
02491     }
02492 
02493 
02494   template<typename _RealType>
02495     template<typename _UniformRandomNumberGenerator>
02496       typename weibull_distribution<_RealType>::result_type
02497       weibull_distribution<_RealType>::
02498       operator()(_UniformRandomNumberGenerator& __urng,
02499                  const param_type& __p)
02500       {
02501         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02502           __aurng(__urng);
02503         return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
02504                                   result_type(1) / __p.a());
02505       }
02506 
02507   template<typename _RealType>
02508     template<typename _ForwardIterator,
02509              typename _UniformRandomNumberGenerator>
02510       void
02511       weibull_distribution<_RealType>::
02512       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02513                       _UniformRandomNumberGenerator& __urng,
02514                       const param_type& __p)
02515       {
02516         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02517         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02518           __aurng(__urng);
02519         auto __inv_a = result_type(1) / __p.a();
02520 
02521         while (__f != __t)
02522           *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
02523                                       __inv_a);
02524       }
02525 
02526   template<typename _RealType, typename _CharT, typename _Traits>
02527     std::basic_ostream<_CharT, _Traits>&
02528     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02529                const weibull_distribution<_RealType>& __x)
02530     {
02531       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02532       typedef typename __ostream_type::ios_base    __ios_base;
02533 
02534       const typename __ios_base::fmtflags __flags = __os.flags();
02535       const _CharT __fill = __os.fill();
02536       const std::streamsize __precision = __os.precision();
02537       const _CharT __space = __os.widen(' ');
02538       __os.flags(__ios_base::scientific | __ios_base::left);
02539       __os.fill(__space);
02540       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02541 
02542       __os << __x.a() << __space << __x.b();
02543 
02544       __os.flags(__flags);
02545       __os.fill(__fill);
02546       __os.precision(__precision);
02547       return __os;
02548     }
02549 
02550   template<typename _RealType, typename _CharT, typename _Traits>
02551     std::basic_istream<_CharT, _Traits>&
02552     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02553                weibull_distribution<_RealType>& __x)
02554     {
02555       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02556       typedef typename __istream_type::ios_base    __ios_base;
02557 
02558       const typename __ios_base::fmtflags __flags = __is.flags();
02559       __is.flags(__ios_base::dec | __ios_base::skipws);
02560 
02561       _RealType __a, __b;
02562       __is >> __a >> __b;
02563       __x.param(typename weibull_distribution<_RealType>::
02564                 param_type(__a, __b));
02565 
02566       __is.flags(__flags);
02567       return __is;
02568     }
02569 
02570 
02571   template<typename _RealType>
02572     template<typename _UniformRandomNumberGenerator>
02573       typename extreme_value_distribution<_RealType>::result_type
02574       extreme_value_distribution<_RealType>::
02575       operator()(_UniformRandomNumberGenerator& __urng,
02576                  const param_type& __p)
02577       {
02578         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02579           __aurng(__urng);
02580         return __p.a() - __p.b() * std::log(-std::log(result_type(1)
02581                                                       - __aurng()));
02582       }
02583 
02584   template<typename _RealType>
02585     template<typename _ForwardIterator,
02586              typename _UniformRandomNumberGenerator>
02587       void
02588       extreme_value_distribution<_RealType>::
02589       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02590                       _UniformRandomNumberGenerator& __urng,
02591                       const param_type& __p)
02592       {
02593         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02594         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
02595           __aurng(__urng);
02596 
02597         while (__f != __t)
02598           *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
02599                                                           - __aurng()));
02600       }
02601 
02602   template<typename _RealType, typename _CharT, typename _Traits>
02603     std::basic_ostream<_CharT, _Traits>&
02604     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02605                const extreme_value_distribution<_RealType>& __x)
02606     {
02607       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02608       typedef typename __ostream_type::ios_base    __ios_base;
02609 
02610       const typename __ios_base::fmtflags __flags = __os.flags();
02611       const _CharT __fill = __os.fill();
02612       const std::streamsize __precision = __os.precision();
02613       const _CharT __space = __os.widen(' ');
02614       __os.flags(__ios_base::scientific | __ios_base::left);
02615       __os.fill(__space);
02616       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02617 
02618       __os << __x.a() << __space << __x.b();
02619 
02620       __os.flags(__flags);
02621       __os.fill(__fill);
02622       __os.precision(__precision);
02623       return __os;
02624     }
02625 
02626   template<typename _RealType, typename _CharT, typename _Traits>
02627     std::basic_istream<_CharT, _Traits>&
02628     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02629                extreme_value_distribution<_RealType>& __x)
02630     {
02631       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02632       typedef typename __istream_type::ios_base    __ios_base;
02633 
02634       const typename __ios_base::fmtflags __flags = __is.flags();
02635       __is.flags(__ios_base::dec | __ios_base::skipws);
02636 
02637       _RealType __a, __b;
02638       __is >> __a >> __b;
02639       __x.param(typename extreme_value_distribution<_RealType>::
02640                 param_type(__a, __b));
02641 
02642       __is.flags(__flags);
02643       return __is;
02644     }
02645 
02646 
02647   template<typename _IntType>
02648     void
02649     discrete_distribution<_IntType>::param_type::
02650     _M_initialize()
02651     {
02652       if (_M_prob.size() < 2)
02653         {
02654           _M_prob.clear();
02655           return;
02656         }
02657 
02658       const double __sum = std::accumulate(_M_prob.begin(),
02659                                            _M_prob.end(), 0.0);
02660       // Now normalize the probabilites.
02661       __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
02662                             __sum);
02663       // Accumulate partial sums.
02664       _M_cp.reserve(_M_prob.size());
02665       std::partial_sum(_M_prob.begin(), _M_prob.end(),
02666                        std::back_inserter(_M_cp));
02667       // Make sure the last cumulative probability is one.
02668       _M_cp[_M_cp.size() - 1] = 1.0;
02669     }
02670 
02671   template<typename _IntType>
02672     template<typename _Func>
02673       discrete_distribution<_IntType>::param_type::
02674       param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
02675       : _M_prob(), _M_cp()
02676       {
02677         const size_t __n = __nw == 0 ? 1 : __nw;
02678         const double __delta = (__xmax - __xmin) / __n;
02679 
02680         _M_prob.reserve(__n);
02681         for (size_t __k = 0; __k < __nw; ++__k)
02682           _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
02683 
02684         _M_initialize();
02685       }
02686 
02687   template<typename _IntType>
02688     template<typename _UniformRandomNumberGenerator>
02689       typename discrete_distribution<_IntType>::result_type
02690       discrete_distribution<_IntType>::
02691       operator()(_UniformRandomNumberGenerator& __urng,
02692                  const param_type& __param)
02693       {
02694         if (__param._M_cp.empty())
02695           return result_type(0);
02696 
02697         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02698           __aurng(__urng);
02699 
02700         const double __p = __aurng();
02701         auto __pos = std::lower_bound(__param._M_cp.begin(),
02702                                       __param._M_cp.end(), __p);
02703 
02704         return __pos - __param._M_cp.begin();
02705       }
02706 
02707   template<typename _IntType>
02708     template<typename _ForwardIterator,
02709              typename _UniformRandomNumberGenerator>
02710       void
02711       discrete_distribution<_IntType>::
02712       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02713                       _UniformRandomNumberGenerator& __urng,
02714                       const param_type& __param)
02715       {
02716         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02717 
02718         if (__param._M_cp.empty())
02719           {
02720             while (__f != __t)
02721               *__f++ = result_type(0);
02722             return;
02723           }
02724 
02725         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02726           __aurng(__urng);
02727 
02728         while (__f != __t)
02729           {
02730             const double __p = __aurng();
02731             auto __pos = std::lower_bound(__param._M_cp.begin(),
02732                                           __param._M_cp.end(), __p);
02733 
02734             *__f++ = __pos - __param._M_cp.begin();
02735           }
02736       }
02737 
02738   template<typename _IntType, typename _CharT, typename _Traits>
02739     std::basic_ostream<_CharT, _Traits>&
02740     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02741                const discrete_distribution<_IntType>& __x)
02742     {
02743       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02744       typedef typename __ostream_type::ios_base    __ios_base;
02745 
02746       const typename __ios_base::fmtflags __flags = __os.flags();
02747       const _CharT __fill = __os.fill();
02748       const std::streamsize __precision = __os.precision();
02749       const _CharT __space = __os.widen(' ');
02750       __os.flags(__ios_base::scientific | __ios_base::left);
02751       __os.fill(__space);
02752       __os.precision(std::numeric_limits<double>::max_digits10);
02753 
02754       std::vector<double> __prob = __x.probabilities();
02755       __os << __prob.size();
02756       for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
02757         __os << __space << *__dit;
02758 
02759       __os.flags(__flags);
02760       __os.fill(__fill);
02761       __os.precision(__precision);
02762       return __os;
02763     }
02764 
02765   template<typename _IntType, typename _CharT, typename _Traits>
02766     std::basic_istream<_CharT, _Traits>&
02767     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02768                discrete_distribution<_IntType>& __x)
02769     {
02770       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02771       typedef typename __istream_type::ios_base    __ios_base;
02772 
02773       const typename __ios_base::fmtflags __flags = __is.flags();
02774       __is.flags(__ios_base::dec | __ios_base::skipws);
02775 
02776       size_t __n;
02777       __is >> __n;
02778 
02779       std::vector<double> __prob_vec;
02780       __prob_vec.reserve(__n);
02781       for (; __n != 0; --__n)
02782         {
02783           double __prob;
02784           __is >> __prob;
02785           __prob_vec.push_back(__prob);
02786         }
02787 
02788       __x.param(typename discrete_distribution<_IntType>::
02789                 param_type(__prob_vec.begin(), __prob_vec.end()));
02790 
02791       __is.flags(__flags);
02792       return __is;
02793     }
02794 
02795 
02796   template<typename _RealType>
02797     void
02798     piecewise_constant_distribution<_RealType>::param_type::
02799     _M_initialize()
02800     {
02801       if (_M_int.size() < 2
02802           || (_M_int.size() == 2
02803               && _M_int[0] == _RealType(0)
02804               && _M_int[1] == _RealType(1)))
02805         {
02806           _M_int.clear();
02807           _M_den.clear();
02808           return;
02809         }
02810 
02811       const double __sum = std::accumulate(_M_den.begin(),
02812                                            _M_den.end(), 0.0);
02813 
02814       __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
02815                             __sum);
02816 
02817       _M_cp.reserve(_M_den.size());
02818       std::partial_sum(_M_den.begin(), _M_den.end(),
02819                        std::back_inserter(_M_cp));
02820 
02821       // Make sure the last cumulative probability is one.
02822       _M_cp[_M_cp.size() - 1] = 1.0;
02823 
02824       for (size_t __k = 0; __k < _M_den.size(); ++__k)
02825         _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
02826     }
02827 
02828   template<typename _RealType>
02829     template<typename _InputIteratorB, typename _InputIteratorW>
02830       piecewise_constant_distribution<_RealType>::param_type::
02831       param_type(_InputIteratorB __bbegin,
02832                  _InputIteratorB __bend,
02833                  _InputIteratorW __wbegin)
02834       : _M_int(), _M_den(), _M_cp()
02835       {
02836         if (__bbegin != __bend)
02837           {
02838             for (;;)
02839               {
02840                 _M_int.push_back(*__bbegin);
02841                 ++__bbegin;
02842                 if (__bbegin == __bend)
02843                   break;
02844 
02845                 _M_den.push_back(*__wbegin);
02846                 ++__wbegin;
02847               }
02848           }
02849 
02850         _M_initialize();
02851       }
02852 
02853   template<typename _RealType>
02854     template<typename _Func>
02855       piecewise_constant_distribution<_RealType>::param_type::
02856       param_type(initializer_list<_RealType> __bl, _Func __fw)
02857       : _M_int(), _M_den(), _M_cp()
02858       {
02859         _M_int.reserve(__bl.size());
02860         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
02861           _M_int.push_back(*__biter);
02862 
02863         _M_den.reserve(_M_int.size() - 1);
02864         for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
02865           _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
02866 
02867         _M_initialize();
02868       }
02869 
02870   template<typename _RealType>
02871     template<typename _Func>
02872       piecewise_constant_distribution<_RealType>::param_type::
02873       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
02874       : _M_int(), _M_den(), _M_cp()
02875       {
02876         const size_t __n = __nw == 0 ? 1 : __nw;
02877         const _RealType __delta = (__xmax - __xmin) / __n;
02878 
02879         _M_int.reserve(__n + 1);
02880         for (size_t __k = 0; __k <= __nw; ++__k)
02881           _M_int.push_back(__xmin + __k * __delta);
02882 
02883         _M_den.reserve(__n);
02884         for (size_t __k = 0; __k < __nw; ++__k)
02885           _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
02886 
02887         _M_initialize();
02888       }
02889 
02890   template<typename _RealType>
02891     template<typename _UniformRandomNumberGenerator>
02892       typename piecewise_constant_distribution<_RealType>::result_type
02893       piecewise_constant_distribution<_RealType>::
02894       operator()(_UniformRandomNumberGenerator& __urng,
02895                  const param_type& __param)
02896       {
02897         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02898           __aurng(__urng);
02899 
02900         const double __p = __aurng();
02901         if (__param._M_cp.empty())
02902           return __p;
02903 
02904         auto __pos = std::lower_bound(__param._M_cp.begin(),
02905                                       __param._M_cp.end(), __p);
02906         const size_t __i = __pos - __param._M_cp.begin();
02907 
02908         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
02909 
02910         return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
02911       }
02912 
02913   template<typename _RealType>
02914     template<typename _ForwardIterator,
02915              typename _UniformRandomNumberGenerator>
02916       void
02917       piecewise_constant_distribution<_RealType>::
02918       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
02919                       _UniformRandomNumberGenerator& __urng,
02920                       const param_type& __param)
02921       {
02922         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
02923         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
02924           __aurng(__urng);
02925 
02926         if (__param._M_cp.empty())
02927           {
02928             while (__f != __t)
02929               *__f++ = __aurng();
02930             return;
02931           }
02932 
02933         while (__f != __t)
02934           {
02935             const double __p = __aurng();
02936 
02937             auto __pos = std::lower_bound(__param._M_cp.begin(),
02938                                           __param._M_cp.end(), __p);
02939             const size_t __i = __pos - __param._M_cp.begin();
02940 
02941             const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
02942 
02943             *__f++ = (__param._M_int[__i]
02944                       + (__p - __pref) / __param._M_den[__i]);
02945           }
02946       }
02947 
02948   template<typename _RealType, typename _CharT, typename _Traits>
02949     std::basic_ostream<_CharT, _Traits>&
02950     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
02951                const piecewise_constant_distribution<_RealType>& __x)
02952     {
02953       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
02954       typedef typename __ostream_type::ios_base    __ios_base;
02955 
02956       const typename __ios_base::fmtflags __flags = __os.flags();
02957       const _CharT __fill = __os.fill();
02958       const std::streamsize __precision = __os.precision();
02959       const _CharT __space = __os.widen(' ');
02960       __os.flags(__ios_base::scientific | __ios_base::left);
02961       __os.fill(__space);
02962       __os.precision(std::numeric_limits<_RealType>::max_digits10);
02963 
02964       std::vector<_RealType> __int = __x.intervals();
02965       __os << __int.size() - 1;
02966 
02967       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
02968         __os << __space << *__xit;
02969 
02970       std::vector<double> __den = __x.densities();
02971       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
02972         __os << __space << *__dit;
02973 
02974       __os.flags(__flags);
02975       __os.fill(__fill);
02976       __os.precision(__precision);
02977       return __os;
02978     }
02979 
02980   template<typename _RealType, typename _CharT, typename _Traits>
02981     std::basic_istream<_CharT, _Traits>&
02982     operator>>(std::basic_istream<_CharT, _Traits>& __is,
02983                piecewise_constant_distribution<_RealType>& __x)
02984     {
02985       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
02986       typedef typename __istream_type::ios_base    __ios_base;
02987 
02988       const typename __ios_base::fmtflags __flags = __is.flags();
02989       __is.flags(__ios_base::dec | __ios_base::skipws);
02990 
02991       size_t __n;
02992       __is >> __n;
02993 
02994       std::vector<_RealType> __int_vec;
02995       __int_vec.reserve(__n + 1);
02996       for (size_t __i = 0; __i <= __n; ++__i)
02997         {
02998           _RealType __int;
02999           __is >> __int;
03000           __int_vec.push_back(__int);
03001         }
03002 
03003       std::vector<double> __den_vec;
03004       __den_vec.reserve(__n);
03005       for (size_t __i = 0; __i < __n; ++__i)
03006         {
03007           double __den;
03008           __is >> __den;
03009           __den_vec.push_back(__den);
03010         }
03011 
03012       __x.param(typename piecewise_constant_distribution<_RealType>::
03013           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
03014 
03015       __is.flags(__flags);
03016       return __is;
03017     }
03018 
03019 
03020   template<typename _RealType>
03021     void
03022     piecewise_linear_distribution<_RealType>::param_type::
03023     _M_initialize()
03024     {
03025       if (_M_int.size() < 2
03026           || (_M_int.size() == 2
03027               && _M_int[0] == _RealType(0)
03028               && _M_int[1] == _RealType(1)
03029               && _M_den[0] == _M_den[1]))
03030         {
03031           _M_int.clear();
03032           _M_den.clear();
03033           return;
03034         }
03035 
03036       double __sum = 0.0;
03037       _M_cp.reserve(_M_int.size() - 1);
03038       _M_m.reserve(_M_int.size() - 1);
03039       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
03040         {
03041           const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
03042           __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
03043           _M_cp.push_back(__sum);
03044           _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
03045         }
03046 
03047       //  Now normalize the densities...
03048       __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
03049                             __sum);
03050       //  ... and partial sums... 
03051       __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
03052       //  ... and slopes.
03053       __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
03054 
03055       //  Make sure the last cumulative probablility is one.
03056       _M_cp[_M_cp.size() - 1] = 1.0;
03057      }
03058 
03059   template<typename _RealType>
03060     template<typename _InputIteratorB, typename _InputIteratorW>
03061       piecewise_linear_distribution<_RealType>::param_type::
03062       param_type(_InputIteratorB __bbegin,
03063                  _InputIteratorB __bend,
03064                  _InputIteratorW __wbegin)
03065       : _M_int(), _M_den(), _M_cp(), _M_m()
03066       {
03067         for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
03068           {
03069             _M_int.push_back(*__bbegin);
03070             _M_den.push_back(*__wbegin);
03071           }
03072 
03073         _M_initialize();
03074       }
03075 
03076   template<typename _RealType>
03077     template<typename _Func>
03078       piecewise_linear_distribution<_RealType>::param_type::
03079       param_type(initializer_list<_RealType> __bl, _Func __fw)
03080       : _M_int(), _M_den(), _M_cp(), _M_m()
03081       {
03082         _M_int.reserve(__bl.size());
03083         _M_den.reserve(__bl.size());
03084         for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
03085           {
03086             _M_int.push_back(*__biter);
03087             _M_den.push_back(__fw(*__biter));
03088           }
03089 
03090         _M_initialize();
03091       }
03092 
03093   template<typename _RealType>
03094     template<typename _Func>
03095       piecewise_linear_distribution<_RealType>::param_type::
03096       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
03097       : _M_int(), _M_den(), _M_cp(), _M_m()
03098       {
03099         const size_t __n = __nw == 0 ? 1 : __nw;
03100         const _RealType __delta = (__xmax - __xmin) / __n;
03101 
03102         _M_int.reserve(__n + 1);
03103         _M_den.reserve(__n + 1);
03104         for (size_t __k = 0; __k <= __nw; ++__k)
03105           {
03106             _M_int.push_back(__xmin + __k * __delta);
03107             _M_den.push_back(__fw(_M_int[__k] + __delta));
03108           }
03109 
03110         _M_initialize();
03111       }
03112 
03113   template<typename _RealType>
03114     template<typename _UniformRandomNumberGenerator>
03115       typename piecewise_linear_distribution<_RealType>::result_type
03116       piecewise_linear_distribution<_RealType>::
03117       operator()(_UniformRandomNumberGenerator& __urng,
03118                  const param_type& __param)
03119       {
03120         __detail::_Adaptor<_UniformRandomNumberGenerator, double>
03121           __aurng(__urng);
03122 
03123         const double __p = __aurng();
03124         if (__param._M_cp.empty())
03125           return __p;
03126 
03127         auto __pos = std::lower_bound(__param._M_cp.begin(),
03128                                       __param._M_cp.end(), __p);
03129         const size_t __i = __pos - __param._M_cp.begin();
03130 
03131         const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
03132 
03133         const double __a = 0.5 * __param._M_m[__i];
03134         const double __b = __param._M_den[__i];
03135         const double __cm = __p - __pref;
03136 
03137         _RealType __x = __param._M_int[__i];
03138         if (__a == 0)
03139           __x += __cm / __b;
03140         else
03141           {
03142             const double __d = __b * __b + 4.0 * __a * __cm;
03143             __x += 0.5 * (std::sqrt(__d) - __b) / __a;
03144           }
03145 
03146         return __x;
03147       }
03148 
03149   template<typename _RealType>
03150     template<typename _ForwardIterator,
03151              typename _UniformRandomNumberGenerator>
03152       void
03153       piecewise_linear_distribution<_RealType>::
03154       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
03155                       _UniformRandomNumberGenerator& __urng,
03156                       const param_type& __param)
03157       {
03158         __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
03159         // We could duplicate everything from operator()...
03160         while (__f != __t)
03161           *__f++ = this->operator()(__urng, __param);
03162       }
03163 
03164   template<typename _RealType, typename _CharT, typename _Traits>
03165     std::basic_ostream<_CharT, _Traits>&
03166     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
03167                const piecewise_linear_distribution<_RealType>& __x)
03168     {
03169       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
03170       typedef typename __ostream_type::ios_base    __ios_base;
03171 
03172       const typename __ios_base::fmtflags __flags = __os.flags();
03173       const _CharT __fill = __os.fill();
03174       const std::streamsize __precision = __os.precision();
03175       const _CharT __space = __os.widen(' ');
03176       __os.flags(__ios_base::scientific | __ios_base::left);
03177       __os.fill(__space);
03178       __os.precision(std::numeric_limits<_RealType>::max_digits10);
03179 
03180       std::vector<_RealType> __int = __x.intervals();
03181       __os << __int.size() - 1;
03182 
03183       for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
03184         __os << __space << *__xit;
03185 
03186       std::vector<double> __den = __x.densities();
03187       for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
03188         __os << __space << *__dit;
03189 
03190       __os.flags(__flags);
03191       __os.fill(__fill);
03192       __os.precision(__precision);
03193       return __os;
03194     }
03195 
03196   template<typename _RealType, typename _CharT, typename _Traits>
03197     std::basic_istream<_CharT, _Traits>&
03198     operator>>(std::basic_istream<_CharT, _Traits>& __is,
03199                piecewise_linear_distribution<_RealType>& __x)
03200     {
03201       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
03202       typedef typename __istream_type::ios_base    __ios_base;
03203 
03204       const typename __ios_base::fmtflags __flags = __is.flags();
03205       __is.flags(__ios_base::dec | __ios_base::skipws);
03206 
03207       size_t __n;
03208       __is >> __n;
03209 
03210       std::vector<_RealType> __int_vec;
03211       __int_vec.reserve(__n + 1);
03212       for (size_t __i = 0; __i <= __n; ++__i)
03213         {
03214           _RealType __int;
03215           __is >> __int;
03216           __int_vec.push_back(__int);
03217         }
03218 
03219       std::vector<double> __den_vec;
03220       __den_vec.reserve(__n + 1);
03221       for (size_t __i = 0; __i <= __n; ++__i)
03222         {
03223           double __den;
03224           __is >> __den;
03225           __den_vec.push_back(__den);
03226         }
03227 
03228       __x.param(typename piecewise_linear_distribution<_RealType>::
03229           param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin()));
03230 
03231       __is.flags(__flags);
03232       return __is;
03233     }
03234 
03235 
03236   template<typename _IntType>
03237     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
03238     {
03239       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
03240         _M_v.push_back(__detail::__mod<result_type,
03241                        __detail::_Shift<result_type, 32>::__value>(*__iter));
03242     }
03243 
03244   template<typename _InputIterator>
03245     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
03246     {
03247       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
03248         _M_v.push_back(__detail::__mod<result_type,
03249                        __detail::_Shift<result_type, 32>::__value>(*__iter));
03250     }
03251 
03252   template<typename _RandomAccessIterator>
03253     void
03254     seed_seq::generate(_RandomAccessIterator __begin,
03255                        _RandomAccessIterator __end)
03256     {
03257       typedef typename iterator_traits<_RandomAccessIterator>::value_type
03258         _Type;
03259 
03260       if (__begin == __end)
03261         return;
03262 
03263       std::fill(__begin, __end, _Type(0x8b8b8b8bu));
03264 
03265       const size_t __n = __end - __begin;
03266       const size_t __s = _M_v.size();
03267       const size_t __t = (__n >= 623) ? 11
03268                        : (__n >=  68) ? 7
03269                        : (__n >=  39) ? 5
03270                        : (__n >=   7) ? 3
03271                        : (__n - 1) / 2;
03272       const size_t __p = (__n - __t) / 2;
03273       const size_t __q = __p + __t;
03274       const size_t __m = std::max(size_t(__s + 1), __n);
03275 
03276       for (size_t __k = 0; __k < __m; ++__k)
03277         {
03278           _Type __arg = (__begin[__k % __n]
03279                          ^ __begin[(__k + __p) % __n]
03280                          ^ __begin[(__k - 1) % __n]);
03281           _Type __r1 = __arg ^ (__arg >> 27);
03282           __r1 = __detail::__mod<_Type,
03283                     __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
03284           _Type __r2 = __r1;
03285           if (__k == 0)
03286             __r2 += __s;
03287           else if (__k <= __s)
03288             __r2 += __k % __n + _M_v[__k - 1];
03289           else
03290             __r2 += __k % __n;
03291           __r2 = __detail::__mod<_Type,
03292                    __detail::_Shift<_Type, 32>::__value>(__r2);
03293           __begin[(__k + __p) % __n] += __r1;
03294           __begin[(__k + __q) % __n] += __r2;
03295           __begin[__k % __n] = __r2;
03296         }
03297 
03298       for (size_t __k = __m; __k < __m + __n; ++__k)
03299         {
03300           _Type __arg = (__begin[__k % __n]
03301                          + __begin[(__k + __p) % __n]
03302                          + __begin[(__k - 1) % __n]);
03303           _Type __r3 = __arg ^ (__arg >> 27);
03304           __r3 = __detail::__mod<_Type,
03305                    __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
03306           _Type __r4 = __r3 - __k % __n;
03307           __r4 = __detail::__mod<_Type,
03308                    __detail::_Shift<_Type, 32>::__value>(__r4);
03309           __begin[(__k + __p) % __n] ^= __r3;
03310           __begin[(__k + __q) % __n] ^= __r4;
03311           __begin[__k % __n] = __r4;
03312         }
03313     }
03314 
03315   template<typename _RealType, size_t __bits,
03316            typename _UniformRandomNumberGenerator>
03317     _RealType
03318     generate_canonical(_UniformRandomNumberGenerator& __urng)
03319     {
03320       static_assert(std::is_floating_point<_RealType>::value,
03321                     "template argument must be a floating point type");
03322 
03323       const size_t __b
03324         = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
03325                    __bits);
03326       const long double __r = static_cast<long double>(__urng.max())
03327                             - static_cast<long double>(__urng.min()) + 1.0L;
03328       const size_t __log2r = std::log(__r) / std::log(2.0L);
03329       const size_t __m = std::max<size_t>(1UL,
03330                                           (__b + __log2r - 1UL) / __log2r);
03331       _RealType __ret;
03332       _RealType __sum = _RealType(0);
03333       _RealType __tmp = _RealType(1);
03334       for (size_t __k = __m; __k != 0; --__k)
03335         {
03336           __sum += _RealType(__urng() - __urng.min()) * __tmp;
03337           __tmp *= __r;
03338         }
03339       __ret = __sum / __tmp;
03340       if (__builtin_expect(__ret >= _RealType(1), 0))
03341         {
03342 #if _GLIBCXX_USE_C99_MATH_TR1
03343           __ret = std::nextafter(_RealType(1), _RealType(0));
03344 #else
03345           __ret = _RealType(1)
03346             - std::numeric_limits<_RealType>::epsilon() / _RealType(2);
03347 #endif
03348         }
03349       return __ret;
03350     }
03351 
03352 _GLIBCXX_END_NAMESPACE_VERSION
03353 } // namespace
03354 
03355 #endif