libstdc++
random.tcc
Go to the documentation of this file.
00001 // Random number extensions -*- C++ -*-
00002 
00003 // Copyright (C) 2012-2015 Free Software Foundation, Inc.
00004 //
00005 // This file is part of the GNU ISO C++ Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the
00007 // terms of the GNU General Public License as published by the
00008 // Free Software Foundation; either version 3, or (at your option)
00009 // any later version.
00010 
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 
00016 // Under Section 7 of GPL version 3, you are granted additional
00017 // permissions described in the GCC Runtime Library Exception, version
00018 // 3.1, as published by the Free Software Foundation.
00019 
00020 // You should have received a copy of the GNU General Public License and
00021 // a copy of the GCC Runtime Library Exception along with this program;
00022 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
00023 // <http://www.gnu.org/licenses/>.
00024 
00025 /** @file ext/random.tcc
00026  *  This is an internal header file, included by other library headers.
00027  *  Do not attempt to use it directly. @headername{ext/random}
00028  */
00029 
00030 #ifndef _EXT_RANDOM_TCC
00031 #define _EXT_RANDOM_TCC 1
00032 
00033 #pragma GCC system_header
00034 
00035 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
00036 {
00037 _GLIBCXX_BEGIN_NAMESPACE_VERSION
00038 
00039 #if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
00040 
00041   template<typename _UIntType, size_t __m,
00042            size_t __pos1, size_t __sl1, size_t __sl2,
00043            size_t __sr1, size_t __sr2,
00044            uint32_t __msk1, uint32_t __msk2,
00045            uint32_t __msk3, uint32_t __msk4,
00046            uint32_t __parity1, uint32_t __parity2,
00047            uint32_t __parity3, uint32_t __parity4>
00048     void simd_fast_mersenne_twister_engine<_UIntType, __m,
00049                                            __pos1, __sl1, __sl2, __sr1, __sr2,
00050                                            __msk1, __msk2, __msk3, __msk4,
00051                                            __parity1, __parity2, __parity3,
00052                                            __parity4>::
00053     seed(_UIntType __seed)
00054     {
00055       _M_state32[0] = static_cast<uint32_t>(__seed);
00056       for (size_t __i = 1; __i < _M_nstate32; ++__i)
00057         _M_state32[__i] = (1812433253UL
00058                            * (_M_state32[__i - 1] ^ (_M_state32[__i - 1] >> 30))
00059                            + __i);
00060       _M_pos = state_size;
00061       _M_period_certification();
00062     }
00063 
00064 
00065   namespace {
00066 
00067     inline uint32_t _Func1(uint32_t __x)
00068     {
00069       return (__x ^ (__x >> 27)) * UINT32_C(1664525);
00070     }
00071 
00072     inline uint32_t _Func2(uint32_t __x)
00073     {
00074       return (__x ^ (__x >> 27)) * UINT32_C(1566083941);
00075     }
00076 
00077   }
00078 
00079 
00080   template<typename _UIntType, size_t __m,
00081            size_t __pos1, size_t __sl1, size_t __sl2,
00082            size_t __sr1, size_t __sr2,
00083            uint32_t __msk1, uint32_t __msk2,
00084            uint32_t __msk3, uint32_t __msk4,
00085            uint32_t __parity1, uint32_t __parity2,
00086            uint32_t __parity3, uint32_t __parity4>
00087     template<typename _Sseq>
00088       typename std::enable_if<std::is_class<_Sseq>::value>::type
00089       simd_fast_mersenne_twister_engine<_UIntType, __m,
00090                                         __pos1, __sl1, __sl2, __sr1, __sr2,
00091                                         __msk1, __msk2, __msk3, __msk4,
00092                                         __parity1, __parity2, __parity3,
00093                                         __parity4>::
00094       seed(_Sseq& __q)
00095       {
00096         size_t __lag;
00097 
00098         if (_M_nstate32 >= 623)
00099           __lag = 11;
00100         else if (_M_nstate32 >= 68)
00101           __lag = 7;
00102         else if (_M_nstate32 >= 39)
00103           __lag = 5;
00104         else
00105           __lag = 3;
00106         const size_t __mid = (_M_nstate32 - __lag) / 2;
00107 
00108         std::fill(_M_state32, _M_state32 + _M_nstate32, UINT32_C(0x8b8b8b8b));
00109         uint32_t __arr[_M_nstate32];
00110         __q.generate(__arr + 0, __arr + _M_nstate32);
00111 
00112         uint32_t __r = _Func1(_M_state32[0] ^ _M_state32[__mid]
00113                               ^ _M_state32[_M_nstate32  - 1]);
00114         _M_state32[__mid] += __r;
00115         __r += _M_nstate32;
00116         _M_state32[__mid + __lag] += __r;
00117         _M_state32[0] = __r;
00118 
00119         for (size_t __i = 1, __j = 0; __j < _M_nstate32; ++__j)
00120           {
00121             __r = _Func1(_M_state32[__i]
00122                          ^ _M_state32[(__i + __mid) % _M_nstate32]
00123                          ^ _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
00124             _M_state32[(__i + __mid) % _M_nstate32] += __r;
00125             __r += __arr[__j] + __i;
00126             _M_state32[(__i + __mid + __lag) % _M_nstate32] += __r;
00127             _M_state32[__i] = __r;
00128             __i = (__i + 1) % _M_nstate32;
00129           }
00130         for (size_t __j = 0; __j < _M_nstate32; ++__j)
00131           {
00132             const size_t __i = (__j + 1) % _M_nstate32;
00133             __r = _Func2(_M_state32[__i]
00134                          + _M_state32[(__i + __mid) % _M_nstate32]
00135                          + _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
00136             _M_state32[(__i + __mid) % _M_nstate32] ^= __r;
00137             __r -= __i;
00138             _M_state32[(__i + __mid + __lag) % _M_nstate32] ^= __r;
00139             _M_state32[__i] = __r;
00140           }
00141 
00142         _M_pos = state_size;
00143         _M_period_certification();
00144       }
00145 
00146 
00147   template<typename _UIntType, size_t __m,
00148            size_t __pos1, size_t __sl1, size_t __sl2,
00149            size_t __sr1, size_t __sr2,
00150            uint32_t __msk1, uint32_t __msk2,
00151            uint32_t __msk3, uint32_t __msk4,
00152            uint32_t __parity1, uint32_t __parity2,
00153            uint32_t __parity3, uint32_t __parity4>
00154     void simd_fast_mersenne_twister_engine<_UIntType, __m,
00155                                            __pos1, __sl1, __sl2, __sr1, __sr2,
00156                                            __msk1, __msk2, __msk3, __msk4,
00157                                            __parity1, __parity2, __parity3,
00158                                            __parity4>::
00159     _M_period_certification(void)
00160     {
00161       static const uint32_t __parity[4] = { __parity1, __parity2,
00162                                             __parity3, __parity4 };
00163       uint32_t __inner = 0;
00164       for (size_t __i = 0; __i < 4; ++__i)
00165         if (__parity[__i] != 0)
00166           __inner ^= _M_state32[__i] & __parity[__i];
00167 
00168       if (__builtin_parity(__inner) & 1)
00169         return;
00170       for (size_t __i = 0; __i < 4; ++__i)
00171         if (__parity[__i] != 0)
00172           {
00173             _M_state32[__i] ^= 1 << (__builtin_ffs(__parity[__i]) - 1);
00174             return;
00175           }
00176       __builtin_unreachable();
00177     }
00178 
00179 
00180   template<typename _UIntType, size_t __m,
00181            size_t __pos1, size_t __sl1, size_t __sl2,
00182            size_t __sr1, size_t __sr2,
00183            uint32_t __msk1, uint32_t __msk2,
00184            uint32_t __msk3, uint32_t __msk4,
00185            uint32_t __parity1, uint32_t __parity2,
00186            uint32_t __parity3, uint32_t __parity4>
00187     void simd_fast_mersenne_twister_engine<_UIntType, __m,
00188                                            __pos1, __sl1, __sl2, __sr1, __sr2,
00189                                            __msk1, __msk2, __msk3, __msk4,
00190                                            __parity1, __parity2, __parity3,
00191                                            __parity4>::
00192     discard(unsigned long long __z)
00193     {
00194       while (__z > state_size - _M_pos)
00195         {
00196           __z -= state_size - _M_pos;
00197 
00198           _M_gen_rand();
00199         }
00200 
00201       _M_pos += __z;
00202     }
00203 
00204 
00205 #ifndef  _GLIBCXX_OPT_HAVE_RANDOM_SFMT_GEN_READ
00206 
00207   namespace {
00208 
00209     template<size_t __shift>
00210       inline void __rshift(uint32_t *__out, const uint32_t *__in)
00211       {
00212         uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
00213                          | static_cast<uint64_t>(__in[2]));
00214         uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
00215                          | static_cast<uint64_t>(__in[0]));
00216 
00217         uint64_t __oh = __th >> (__shift * 8);
00218         uint64_t __ol = __tl >> (__shift * 8);
00219         __ol |= __th << (64 - __shift * 8);
00220         __out[1] = static_cast<uint32_t>(__ol >> 32);
00221         __out[0] = static_cast<uint32_t>(__ol);
00222         __out[3] = static_cast<uint32_t>(__oh >> 32);
00223         __out[2] = static_cast<uint32_t>(__oh);
00224       }
00225 
00226 
00227     template<size_t __shift>
00228       inline void __lshift(uint32_t *__out, const uint32_t *__in)
00229       {
00230         uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
00231                          | static_cast<uint64_t>(__in[2]));
00232         uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
00233                          | static_cast<uint64_t>(__in[0]));
00234 
00235         uint64_t __oh = __th << (__shift * 8);
00236         uint64_t __ol = __tl << (__shift * 8);
00237         __oh |= __tl >> (64 - __shift * 8);
00238         __out[1] = static_cast<uint32_t>(__ol >> 32);
00239         __out[0] = static_cast<uint32_t>(__ol);
00240         __out[3] = static_cast<uint32_t>(__oh >> 32);
00241         __out[2] = static_cast<uint32_t>(__oh);
00242       }
00243 
00244 
00245     template<size_t __sl1, size_t __sl2, size_t __sr1, size_t __sr2,
00246              uint32_t __msk1, uint32_t __msk2, uint32_t __msk3, uint32_t __msk4>
00247       inline void __recursion(uint32_t *__r,
00248                               const uint32_t *__a, const uint32_t *__b,
00249                               const uint32_t *__c, const uint32_t *__d)
00250       {
00251         uint32_t __x[4];
00252         uint32_t __y[4];
00253 
00254         __lshift<__sl2>(__x, __a);
00255         __rshift<__sr2>(__y, __c);
00256         __r[0] = (__a[0] ^ __x[0] ^ ((__b[0] >> __sr1) & __msk1)
00257                   ^ __y[0] ^ (__d[0] << __sl1));
00258         __r[1] = (__a[1] ^ __x[1] ^ ((__b[1] >> __sr1) & __msk2)
00259                   ^ __y[1] ^ (__d[1] << __sl1));
00260         __r[2] = (__a[2] ^ __x[2] ^ ((__b[2] >> __sr1) & __msk3)
00261                   ^ __y[2] ^ (__d[2] << __sl1));
00262         __r[3] = (__a[3] ^ __x[3] ^ ((__b[3] >> __sr1) & __msk4)
00263                   ^ __y[3] ^ (__d[3] << __sl1));
00264       }
00265 
00266   }
00267 
00268 
00269   template<typename _UIntType, size_t __m,
00270            size_t __pos1, size_t __sl1, size_t __sl2,
00271            size_t __sr1, size_t __sr2,
00272            uint32_t __msk1, uint32_t __msk2,
00273            uint32_t __msk3, uint32_t __msk4,
00274            uint32_t __parity1, uint32_t __parity2,
00275            uint32_t __parity3, uint32_t __parity4>
00276     void simd_fast_mersenne_twister_engine<_UIntType, __m,
00277                                            __pos1, __sl1, __sl2, __sr1, __sr2,
00278                                            __msk1, __msk2, __msk3, __msk4,
00279                                            __parity1, __parity2, __parity3,
00280                                            __parity4>::
00281     _M_gen_rand(void)
00282     {
00283       const uint32_t *__r1 = &_M_state32[_M_nstate32 - 8];
00284       const uint32_t *__r2 = &_M_state32[_M_nstate32 - 4];
00285       static constexpr size_t __pos1_32 = __pos1 * 4;
00286 
00287       size_t __i;
00288       for (__i = 0; __i < _M_nstate32 - __pos1_32; __i += 4)
00289         {
00290           __recursion<__sl1, __sl2, __sr1, __sr2,
00291                       __msk1, __msk2, __msk3, __msk4>
00292             (&_M_state32[__i], &_M_state32[__i],
00293              &_M_state32[__i + __pos1_32], __r1, __r2);
00294           __r1 = __r2;
00295           __r2 = &_M_state32[__i];
00296         }
00297 
00298       for (; __i < _M_nstate32; __i += 4)
00299         {
00300           __recursion<__sl1, __sl2, __sr1, __sr2,
00301                       __msk1, __msk2, __msk3, __msk4>
00302             (&_M_state32[__i], &_M_state32[__i],
00303              &_M_state32[__i + __pos1_32 - _M_nstate32], __r1, __r2);
00304           __r1 = __r2;
00305           __r2 = &_M_state32[__i];
00306         }
00307 
00308       _M_pos = 0;
00309     }
00310 
00311 #endif
00312 
00313 #ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_OPERATOREQUAL
00314   template<typename _UIntType, size_t __m,
00315            size_t __pos1, size_t __sl1, size_t __sl2,
00316            size_t __sr1, size_t __sr2,
00317            uint32_t __msk1, uint32_t __msk2,
00318            uint32_t __msk3, uint32_t __msk4,
00319            uint32_t __parity1, uint32_t __parity2,
00320            uint32_t __parity3, uint32_t __parity4>
00321     bool
00322     operator==(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
00323                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
00324                __msk1, __msk2, __msk3, __msk4,
00325                __parity1, __parity2, __parity3, __parity4>& __lhs,
00326                const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
00327                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
00328                __msk1, __msk2, __msk3, __msk4,
00329                __parity1, __parity2, __parity3, __parity4>& __rhs)
00330     {
00331       typedef __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
00332                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
00333                __msk1, __msk2, __msk3, __msk4,
00334                __parity1, __parity2, __parity3, __parity4> __engine;
00335       return (std::equal(__lhs._M_stateT,
00336                          __lhs._M_stateT + __engine::state_size,
00337                          __rhs._M_stateT)
00338               && __lhs._M_pos == __rhs._M_pos);
00339     }
00340 #endif
00341 
00342   template<typename _UIntType, size_t __m,
00343            size_t __pos1, size_t __sl1, size_t __sl2,
00344            size_t __sr1, size_t __sr2,
00345            uint32_t __msk1, uint32_t __msk2,
00346            uint32_t __msk3, uint32_t __msk4,
00347            uint32_t __parity1, uint32_t __parity2,
00348            uint32_t __parity3, uint32_t __parity4,
00349            typename _CharT, typename _Traits>
00350     std::basic_ostream<_CharT, _Traits>&
00351     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00352                const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
00353                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
00354                __msk1, __msk2, __msk3, __msk4,
00355                __parity1, __parity2, __parity3, __parity4>& __x)
00356     {
00357       typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
00358       typedef typename __ostream_type::ios_base __ios_base;
00359 
00360       const typename __ios_base::fmtflags __flags = __os.flags();
00361       const _CharT __fill = __os.fill();
00362       const _CharT __space = __os.widen(' ');
00363       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
00364       __os.fill(__space);
00365 
00366       for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
00367         __os << __x._M_state32[__i] << __space;
00368       __os << __x._M_pos;
00369 
00370       __os.flags(__flags);
00371       __os.fill(__fill);
00372       return __os;
00373     }
00374 
00375 
00376   template<typename _UIntType, size_t __m,
00377            size_t __pos1, size_t __sl1, size_t __sl2,
00378            size_t __sr1, size_t __sr2,
00379            uint32_t __msk1, uint32_t __msk2,
00380            uint32_t __msk3, uint32_t __msk4,
00381            uint32_t __parity1, uint32_t __parity2,
00382            uint32_t __parity3, uint32_t __parity4,
00383            typename _CharT, typename _Traits>
00384     std::basic_istream<_CharT, _Traits>&
00385     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00386                __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
00387                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
00388                __msk1, __msk2, __msk3, __msk4,
00389                __parity1, __parity2, __parity3, __parity4>& __x)
00390     {
00391       typedef std::basic_istream<_CharT, _Traits> __istream_type;
00392       typedef typename __istream_type::ios_base __ios_base;
00393 
00394       const typename __ios_base::fmtflags __flags = __is.flags();
00395       __is.flags(__ios_base::dec | __ios_base::skipws);
00396 
00397       for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
00398         __is >> __x._M_state32[__i];
00399       __is >> __x._M_pos;
00400 
00401       __is.flags(__flags);
00402       return __is;
00403     }
00404 
00405 #endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
00406 
00407   /**
00408    * Iteration method due to M.D. J<o:>hnk.
00409    *
00410    * M.D. J<o:>hnk, Erzeugung von betaverteilten und gammaverteilten
00411    * Zufallszahlen, Metrika, Volume 8, 1964
00412    */
00413   template<typename _RealType>
00414     template<typename _UniformRandomNumberGenerator>
00415       typename beta_distribution<_RealType>::result_type
00416       beta_distribution<_RealType>::
00417       operator()(_UniformRandomNumberGenerator& __urng,
00418                  const param_type& __param)
00419       {
00420         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
00421           __aurng(__urng);
00422 
00423         result_type __x, __y;
00424         do
00425           {
00426             __x = std::exp(std::log(__aurng()) / __param.alpha());
00427             __y = std::exp(std::log(__aurng()) / __param.beta());
00428           }
00429         while (__x + __y > result_type(1));
00430 
00431         return __x / (__x + __y);
00432       }
00433 
00434   template<typename _RealType>
00435     template<typename _OutputIterator,
00436              typename _UniformRandomNumberGenerator>
00437       void
00438       beta_distribution<_RealType>::
00439       __generate_impl(_OutputIterator __f, _OutputIterator __t,
00440                       _UniformRandomNumberGenerator& __urng,
00441                       const param_type& __param)
00442       {
00443         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
00444 
00445         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
00446           __aurng(__urng);
00447 
00448         while (__f != __t)
00449           {
00450             result_type __x, __y;
00451             do
00452               {
00453                 __x = std::exp(std::log(__aurng()) / __param.alpha());
00454                 __y = std::exp(std::log(__aurng()) / __param.beta());
00455               }
00456             while (__x + __y > result_type(1));
00457 
00458             *__f++ = __x / (__x + __y);
00459           }
00460       }
00461 
00462   template<typename _RealType, typename _CharT, typename _Traits>
00463     std::basic_ostream<_CharT, _Traits>&
00464     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00465                const __gnu_cxx::beta_distribution<_RealType>& __x)
00466     {
00467       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00468       typedef typename __ostream_type::ios_base    __ios_base;
00469 
00470       const typename __ios_base::fmtflags __flags = __os.flags();
00471       const _CharT __fill = __os.fill();
00472       const std::streamsize __precision = __os.precision();
00473       const _CharT __space = __os.widen(' ');
00474       __os.flags(__ios_base::scientific | __ios_base::left);
00475       __os.fill(__space);
00476       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00477 
00478       __os << __x.alpha() << __space << __x.beta();
00479 
00480       __os.flags(__flags);
00481       __os.fill(__fill);
00482       __os.precision(__precision);
00483       return __os;
00484     }
00485 
00486   template<typename _RealType, typename _CharT, typename _Traits>
00487     std::basic_istream<_CharT, _Traits>&
00488     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00489                __gnu_cxx::beta_distribution<_RealType>& __x)
00490     {
00491       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00492       typedef typename __istream_type::ios_base    __ios_base;
00493 
00494       const typename __ios_base::fmtflags __flags = __is.flags();
00495       __is.flags(__ios_base::dec | __ios_base::skipws);
00496 
00497       _RealType __alpha_val, __beta_val;
00498       __is >> __alpha_val >> __beta_val;
00499       __x.param(typename __gnu_cxx::beta_distribution<_RealType>::
00500                 param_type(__alpha_val, __beta_val));
00501 
00502       __is.flags(__flags);
00503       return __is;
00504     }
00505 
00506 
00507   template<std::size_t _Dimen, typename _RealType>
00508     template<typename _InputIterator1, typename _InputIterator2>
00509       void
00510       normal_mv_distribution<_Dimen, _RealType>::param_type::
00511       _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
00512                    _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
00513       {
00514         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
00515         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
00516         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
00517                   _M_mean.end(), _RealType(0));
00518 
00519         // Perform the Cholesky decomposition
00520         auto __w = _M_t.begin();
00521         for (size_t __j = 0; __j < _Dimen; ++__j)
00522           {
00523             _RealType __sum = _RealType(0);
00524 
00525             auto __slitbegin = __w;
00526             auto __cit = _M_t.begin();
00527             for (size_t __i = 0; __i < __j; ++__i)
00528               {
00529                 auto __slit = __slitbegin;
00530                 _RealType __s = *__varcovbegin++;
00531                 for (size_t __k = 0; __k < __i; ++__k)
00532                   __s -= *__slit++ * *__cit++;
00533 
00534                 *__w++ = __s /= *__cit++;
00535                 __sum += __s * __s;
00536               }
00537 
00538             __sum = *__varcovbegin - __sum;
00539             if (__builtin_expect(__sum <= _RealType(0), 0))
00540               std::__throw_runtime_error(__N("normal_mv_distribution::"
00541                                              "param_type::_M_init_full"));
00542             *__w++ = std::sqrt(__sum);
00543 
00544             std::advance(__varcovbegin, _Dimen - __j);
00545           }
00546       }
00547 
00548   template<std::size_t _Dimen, typename _RealType>
00549     template<typename _InputIterator1, typename _InputIterator2>
00550       void
00551       normal_mv_distribution<_Dimen, _RealType>::param_type::
00552       _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
00553                     _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
00554       {
00555         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
00556         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
00557         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
00558                   _M_mean.end(), _RealType(0));
00559 
00560         // Perform the Cholesky decomposition
00561         auto __w = _M_t.begin();
00562         for (size_t __j = 0; __j < _Dimen; ++__j)
00563           {
00564             _RealType __sum = _RealType(0);
00565 
00566             auto __slitbegin = __w;
00567             auto __cit = _M_t.begin();
00568             for (size_t __i = 0; __i < __j; ++__i)
00569               {
00570                 auto __slit = __slitbegin;
00571                 _RealType __s = *__varcovbegin++;
00572                 for (size_t __k = 0; __k < __i; ++__k)
00573                   __s -= *__slit++ * *__cit++;
00574 
00575                 *__w++ = __s /= *__cit++;
00576                 __sum += __s * __s;
00577               }
00578 
00579             __sum = *__varcovbegin++ - __sum;
00580             if (__builtin_expect(__sum <= _RealType(0), 0))
00581               std::__throw_runtime_error(__N("normal_mv_distribution::"
00582                                              "param_type::_M_init_full"));
00583             *__w++ = std::sqrt(__sum);
00584           }
00585       }
00586 
00587   template<std::size_t _Dimen, typename _RealType>
00588     template<typename _InputIterator1, typename _InputIterator2>
00589       void
00590       normal_mv_distribution<_Dimen, _RealType>::param_type::
00591       _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
00592                        _InputIterator2 __varbegin, _InputIterator2 __varend)
00593       {
00594         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
00595         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
00596         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
00597                   _M_mean.end(), _RealType(0));
00598 
00599         auto __w = _M_t.begin();
00600         size_t __step = 0;
00601         while (__varbegin != __varend)
00602           {
00603             std::fill_n(__w, __step, _RealType(0));
00604             __w += __step++;
00605             if (__builtin_expect(*__varbegin < _RealType(0), 0))
00606               std::__throw_runtime_error(__N("normal_mv_distribution::"
00607                                              "param_type::_M_init_diagonal"));
00608             *__w++ = std::sqrt(*__varbegin++);
00609           }
00610       }
00611 
00612   template<std::size_t _Dimen, typename _RealType>
00613     template<typename _UniformRandomNumberGenerator>
00614       typename normal_mv_distribution<_Dimen, _RealType>::result_type
00615       normal_mv_distribution<_Dimen, _RealType>::
00616       operator()(_UniformRandomNumberGenerator& __urng,
00617                  const param_type& __param)
00618       {
00619         result_type __ret;
00620 
00621         _M_nd.__generate(__ret.begin(), __ret.end(), __urng);
00622 
00623         auto __t_it = __param._M_t.crbegin();
00624         for (size_t __i = _Dimen; __i > 0; --__i)
00625           {
00626             _RealType __sum = _RealType(0);
00627             for (size_t __j = __i; __j > 0; --__j)
00628               __sum += __ret[__j - 1] * *__t_it++;
00629             __ret[__i - 1] = __sum;
00630           }
00631 
00632         return __ret;
00633       }
00634 
00635   template<std::size_t _Dimen, typename _RealType>
00636     template<typename _ForwardIterator, typename _UniformRandomNumberGenerator>
00637       void
00638       normal_mv_distribution<_Dimen, _RealType>::
00639       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
00640                       _UniformRandomNumberGenerator& __urng,
00641                       const param_type& __param)
00642       {
00643         __glibcxx_function_requires(_Mutable_ForwardIteratorConcept<
00644                                     _ForwardIterator>)
00645         while (__f != __t)
00646           *__f++ = this->operator()(__urng, __param);
00647       }
00648 
00649   template<size_t _Dimen, typename _RealType>
00650     bool
00651     operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
00652                __d1,
00653                const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
00654                __d2)
00655     {
00656       return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd;
00657     }
00658 
00659   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
00660     std::basic_ostream<_CharT, _Traits>&
00661     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00662                const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
00663     {
00664       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00665       typedef typename __ostream_type::ios_base    __ios_base;
00666 
00667       const typename __ios_base::fmtflags __flags = __os.flags();
00668       const _CharT __fill = __os.fill();
00669       const std::streamsize __precision = __os.precision();
00670       const _CharT __space = __os.widen(' ');
00671       __os.flags(__ios_base::scientific | __ios_base::left);
00672       __os.fill(__space);
00673       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00674 
00675       auto __mean = __x._M_param.mean();
00676       for (auto __it : __mean)
00677         __os << __it << __space;
00678       auto __t = __x._M_param.varcov();
00679       for (auto __it : __t)
00680         __os << __it << __space;
00681 
00682       __os << __x._M_nd;
00683 
00684       __os.flags(__flags);
00685       __os.fill(__fill);
00686       __os.precision(__precision);
00687       return __os;
00688     }
00689 
00690   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
00691     std::basic_istream<_CharT, _Traits>&
00692     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00693                __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
00694     {
00695       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00696       typedef typename __istream_type::ios_base    __ios_base;
00697 
00698       const typename __ios_base::fmtflags __flags = __is.flags();
00699       __is.flags(__ios_base::dec | __ios_base::skipws);
00700 
00701       std::array<_RealType, _Dimen> __mean;
00702       for (auto& __it : __mean)
00703         __is >> __it;
00704       std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov;
00705       for (auto& __it : __varcov)
00706         __is >> __it;
00707 
00708       __is >> __x._M_nd;
00709 
00710       __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
00711                 param_type(__mean.begin(), __mean.end(),
00712                            __varcov.begin(), __varcov.end()));
00713 
00714       __is.flags(__flags);
00715       return __is;
00716     }
00717 
00718 
00719   template<typename _RealType>
00720     template<typename _OutputIterator,
00721              typename _UniformRandomNumberGenerator>
00722       void
00723       rice_distribution<_RealType>::
00724       __generate_impl(_OutputIterator __f, _OutputIterator __t,
00725                       _UniformRandomNumberGenerator& __urng,
00726                       const param_type& __p)
00727       {
00728         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
00729 
00730         while (__f != __t)
00731           {
00732             typename std::normal_distribution<result_type>::param_type
00733               __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
00734             result_type __x = this->_M_ndx(__px, __urng);
00735             result_type __y = this->_M_ndy(__py, __urng);
00736 #if _GLIBCXX_USE_C99_MATH_TR1
00737             *__f++ = std::hypot(__x, __y);
00738 #else
00739             *__f++ = std::sqrt(__x * __x + __y * __y);
00740 #endif
00741           }
00742       }
00743 
00744   template<typename _RealType, typename _CharT, typename _Traits>
00745     std::basic_ostream<_CharT, _Traits>&
00746     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00747                const rice_distribution<_RealType>& __x)
00748     {
00749       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00750       typedef typename __ostream_type::ios_base    __ios_base;
00751 
00752       const typename __ios_base::fmtflags __flags = __os.flags();
00753       const _CharT __fill = __os.fill();
00754       const std::streamsize __precision = __os.precision();
00755       const _CharT __space = __os.widen(' ');
00756       __os.flags(__ios_base::scientific | __ios_base::left);
00757       __os.fill(__space);
00758       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00759 
00760       __os << __x.nu() << __space << __x.sigma();
00761       __os << __space << __x._M_ndx;
00762       __os << __space << __x._M_ndy;
00763 
00764       __os.flags(__flags);
00765       __os.fill(__fill);
00766       __os.precision(__precision);
00767       return __os;
00768     }
00769 
00770   template<typename _RealType, typename _CharT, typename _Traits>
00771     std::basic_istream<_CharT, _Traits>&
00772     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00773                rice_distribution<_RealType>& __x)
00774     {
00775       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00776       typedef typename __istream_type::ios_base    __ios_base;
00777 
00778       const typename __ios_base::fmtflags __flags = __is.flags();
00779       __is.flags(__ios_base::dec | __ios_base::skipws);
00780 
00781       _RealType __nu_val, __sigma_val;
00782       __is >> __nu_val >> __sigma_val;
00783       __is >> __x._M_ndx;
00784       __is >> __x._M_ndy;
00785       __x.param(typename rice_distribution<_RealType>::
00786                 param_type(__nu_val, __sigma_val));
00787 
00788       __is.flags(__flags);
00789       return __is;
00790     }
00791 
00792 
00793   template<typename _RealType>
00794     template<typename _OutputIterator,
00795              typename _UniformRandomNumberGenerator>
00796       void
00797       nakagami_distribution<_RealType>::
00798       __generate_impl(_OutputIterator __f, _OutputIterator __t,
00799                       _UniformRandomNumberGenerator& __urng,
00800                       const param_type& __p)
00801       {
00802         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
00803 
00804         typename std::gamma_distribution<result_type>::param_type
00805           __pg(__p.mu(), __p.omega() / __p.mu());
00806         while (__f != __t)
00807           *__f++ = std::sqrt(this->_M_gd(__pg, __urng));
00808       }
00809 
00810   template<typename _RealType, typename _CharT, typename _Traits>
00811     std::basic_ostream<_CharT, _Traits>&
00812     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00813                const nakagami_distribution<_RealType>& __x)
00814     {
00815       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00816       typedef typename __ostream_type::ios_base    __ios_base;
00817 
00818       const typename __ios_base::fmtflags __flags = __os.flags();
00819       const _CharT __fill = __os.fill();
00820       const std::streamsize __precision = __os.precision();
00821       const _CharT __space = __os.widen(' ');
00822       __os.flags(__ios_base::scientific | __ios_base::left);
00823       __os.fill(__space);
00824       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00825 
00826       __os << __x.mu() << __space << __x.omega();
00827       __os << __space << __x._M_gd;
00828 
00829       __os.flags(__flags);
00830       __os.fill(__fill);
00831       __os.precision(__precision);
00832       return __os;
00833     }
00834 
00835   template<typename _RealType, typename _CharT, typename _Traits>
00836     std::basic_istream<_CharT, _Traits>&
00837     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00838                nakagami_distribution<_RealType>& __x)
00839     {
00840       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00841       typedef typename __istream_type::ios_base    __ios_base;
00842 
00843       const typename __ios_base::fmtflags __flags = __is.flags();
00844       __is.flags(__ios_base::dec | __ios_base::skipws);
00845 
00846       _RealType __mu_val, __omega_val;
00847       __is >> __mu_val >> __omega_val;
00848       __is >> __x._M_gd;
00849       __x.param(typename nakagami_distribution<_RealType>::
00850                 param_type(__mu_val, __omega_val));
00851 
00852       __is.flags(__flags);
00853       return __is;
00854     }
00855 
00856 
00857   template<typename _RealType>
00858     template<typename _OutputIterator,
00859              typename _UniformRandomNumberGenerator>
00860       void
00861       pareto_distribution<_RealType>::
00862       __generate_impl(_OutputIterator __f, _OutputIterator __t,
00863                       _UniformRandomNumberGenerator& __urng,
00864                       const param_type& __p)
00865       {
00866         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
00867 
00868         result_type __mu_val = __p.mu();
00869         result_type __malphinv = -result_type(1) / __p.alpha();
00870         while (__f != __t)
00871           *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv);
00872       }
00873 
00874   template<typename _RealType, typename _CharT, typename _Traits>
00875     std::basic_ostream<_CharT, _Traits>&
00876     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00877                const pareto_distribution<_RealType>& __x)
00878     {
00879       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00880       typedef typename __ostream_type::ios_base    __ios_base;
00881 
00882       const typename __ios_base::fmtflags __flags = __os.flags();
00883       const _CharT __fill = __os.fill();
00884       const std::streamsize __precision = __os.precision();
00885       const _CharT __space = __os.widen(' ');
00886       __os.flags(__ios_base::scientific | __ios_base::left);
00887       __os.fill(__space);
00888       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00889 
00890       __os << __x.alpha() << __space << __x.mu();
00891       __os << __space << __x._M_ud;
00892 
00893       __os.flags(__flags);
00894       __os.fill(__fill);
00895       __os.precision(__precision);
00896       return __os;
00897     }
00898 
00899   template<typename _RealType, typename _CharT, typename _Traits>
00900     std::basic_istream<_CharT, _Traits>&
00901     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00902                pareto_distribution<_RealType>& __x)
00903     {
00904       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
00905       typedef typename __istream_type::ios_base    __ios_base;
00906 
00907       const typename __ios_base::fmtflags __flags = __is.flags();
00908       __is.flags(__ios_base::dec | __ios_base::skipws);
00909 
00910       _RealType __alpha_val, __mu_val;
00911       __is >> __alpha_val >> __mu_val;
00912       __is >> __x._M_ud;
00913       __x.param(typename pareto_distribution<_RealType>::
00914                 param_type(__alpha_val, __mu_val));
00915 
00916       __is.flags(__flags);
00917       return __is;
00918     }
00919 
00920 
00921   template<typename _RealType>
00922     template<typename _UniformRandomNumberGenerator>
00923       typename k_distribution<_RealType>::result_type
00924       k_distribution<_RealType>::
00925       operator()(_UniformRandomNumberGenerator& __urng)
00926       {
00927         result_type __x = this->_M_gd1(__urng);
00928         result_type __y = this->_M_gd2(__urng);
00929         return std::sqrt(__x * __y);
00930       }
00931 
00932   template<typename _RealType>
00933     template<typename _UniformRandomNumberGenerator>
00934       typename k_distribution<_RealType>::result_type
00935       k_distribution<_RealType>::
00936       operator()(_UniformRandomNumberGenerator& __urng,
00937                  const param_type& __p)
00938       {
00939         typename std::gamma_distribution<result_type>::param_type
00940           __p1(__p.lambda(), result_type(1) / __p.lambda()),
00941           __p2(__p.nu(), __p.mu() / __p.nu());
00942         result_type __x = this->_M_gd1(__p1, __urng);
00943         result_type __y = this->_M_gd2(__p2, __urng);
00944         return std::sqrt(__x * __y);
00945       }
00946 
00947   template<typename _RealType>
00948     template<typename _OutputIterator,
00949              typename _UniformRandomNumberGenerator>
00950       void
00951       k_distribution<_RealType>::
00952       __generate_impl(_OutputIterator __f, _OutputIterator __t,
00953                       _UniformRandomNumberGenerator& __urng,
00954                       const param_type& __p)
00955       {
00956         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
00957 
00958         typename std::gamma_distribution<result_type>::param_type
00959           __p1(__p.lambda(), result_type(1) / __p.lambda()),
00960           __p2(__p.nu(), __p.mu() / __p.nu());
00961         while (__f != __t)
00962           {
00963             result_type __x = this->_M_gd1(__p1, __urng);
00964             result_type __y = this->_M_gd2(__p2, __urng);
00965             *__f++ = std::sqrt(__x * __y);
00966           }
00967       }
00968 
00969   template<typename _RealType, typename _CharT, typename _Traits>
00970     std::basic_ostream<_CharT, _Traits>&
00971     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
00972                const k_distribution<_RealType>& __x)
00973     {
00974       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
00975       typedef typename __ostream_type::ios_base    __ios_base;
00976 
00977       const typename __ios_base::fmtflags __flags = __os.flags();
00978       const _CharT __fill = __os.fill();
00979       const std::streamsize __precision = __os.precision();
00980       const _CharT __space = __os.widen(' ');
00981       __os.flags(__ios_base::scientific | __ios_base::left);
00982       __os.fill(__space);
00983       __os.precision(std::numeric_limits<_RealType>::max_digits10);
00984 
00985       __os << __x.lambda() << __space << __x.mu() << __space << __x.nu();
00986       __os << __space << __x._M_gd1;
00987       __os << __space << __x._M_gd2;
00988 
00989       __os.flags(__flags);
00990       __os.fill(__fill);
00991       __os.precision(__precision);
00992       return __os;
00993     }
00994 
00995   template<typename _RealType, typename _CharT, typename _Traits>
00996     std::basic_istream<_CharT, _Traits>&
00997     operator>>(std::basic_istream<_CharT, _Traits>& __is,
00998                k_distribution<_RealType>& __x)
00999     {
01000       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01001       typedef typename __istream_type::ios_base    __ios_base;
01002 
01003       const typename __ios_base::fmtflags __flags = __is.flags();
01004       __is.flags(__ios_base::dec | __ios_base::skipws);
01005 
01006       _RealType __lambda_val, __mu_val, __nu_val;
01007       __is >> __lambda_val >> __mu_val >> __nu_val;
01008       __is >> __x._M_gd1;
01009       __is >> __x._M_gd2;
01010       __x.param(typename k_distribution<_RealType>::
01011                 param_type(__lambda_val, __mu_val, __nu_val));
01012 
01013       __is.flags(__flags);
01014       return __is;
01015     }
01016 
01017 
01018   template<typename _RealType>
01019     template<typename _OutputIterator,
01020              typename _UniformRandomNumberGenerator>
01021       void
01022       arcsine_distribution<_RealType>::
01023       __generate_impl(_OutputIterator __f, _OutputIterator __t,
01024                       _UniformRandomNumberGenerator& __urng,
01025                       const param_type& __p)
01026       {
01027         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
01028 
01029         result_type __dif = __p.b() - __p.a();
01030         result_type __sum = __p.a() + __p.b();
01031         while (__f != __t)
01032           {
01033             result_type __x = std::sin(this->_M_ud(__urng));
01034             *__f++ = (__x * __dif + __sum) / result_type(2);
01035           }
01036       }
01037 
01038   template<typename _RealType, typename _CharT, typename _Traits>
01039     std::basic_ostream<_CharT, _Traits>&
01040     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01041                const arcsine_distribution<_RealType>& __x)
01042     {
01043       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01044       typedef typename __ostream_type::ios_base    __ios_base;
01045 
01046       const typename __ios_base::fmtflags __flags = __os.flags();
01047       const _CharT __fill = __os.fill();
01048       const std::streamsize __precision = __os.precision();
01049       const _CharT __space = __os.widen(' ');
01050       __os.flags(__ios_base::scientific | __ios_base::left);
01051       __os.fill(__space);
01052       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01053 
01054       __os << __x.a() << __space << __x.b();
01055       __os << __space << __x._M_ud;
01056 
01057       __os.flags(__flags);
01058       __os.fill(__fill);
01059       __os.precision(__precision);
01060       return __os;
01061     }
01062 
01063   template<typename _RealType, typename _CharT, typename _Traits>
01064     std::basic_istream<_CharT, _Traits>&
01065     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01066                arcsine_distribution<_RealType>& __x)
01067     {
01068       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01069       typedef typename __istream_type::ios_base    __ios_base;
01070 
01071       const typename __ios_base::fmtflags __flags = __is.flags();
01072       __is.flags(__ios_base::dec | __ios_base::skipws);
01073 
01074       _RealType __a, __b;
01075       __is >> __a >> __b;
01076       __is >> __x._M_ud;
01077       __x.param(typename arcsine_distribution<_RealType>::
01078                 param_type(__a, __b));
01079 
01080       __is.flags(__flags);
01081       return __is;
01082     }
01083 
01084 
01085   template<typename _RealType>
01086     template<typename _UniformRandomNumberGenerator>
01087       typename hoyt_distribution<_RealType>::result_type
01088       hoyt_distribution<_RealType>::
01089       operator()(_UniformRandomNumberGenerator& __urng)
01090       {
01091         result_type __x = this->_M_ad(__urng);
01092         result_type __y = this->_M_ed(__urng);
01093         return (result_type(2) * this->q()
01094                   / (result_type(1) + this->q() * this->q()))
01095                * std::sqrt(this->omega() * __x * __y);
01096       }
01097 
01098   template<typename _RealType>
01099     template<typename _UniformRandomNumberGenerator>
01100       typename hoyt_distribution<_RealType>::result_type
01101       hoyt_distribution<_RealType>::
01102       operator()(_UniformRandomNumberGenerator& __urng,
01103                  const param_type& __p)
01104       {
01105         result_type __q2 = __p.q() * __p.q();
01106         result_type __num = result_type(0.5L) * (result_type(1) + __q2);
01107         typename __gnu_cxx::arcsine_distribution<result_type>::param_type
01108           __pa(__num, __num / __q2);
01109         result_type __x = this->_M_ad(__pa, __urng);
01110         result_type __y = this->_M_ed(__urng);
01111         return (result_type(2) * __p.q() / (result_type(1) + __q2))
01112                * std::sqrt(__p.omega() * __x * __y);
01113       }
01114 
01115   template<typename _RealType>
01116     template<typename _OutputIterator,
01117              typename _UniformRandomNumberGenerator>
01118       void
01119       hoyt_distribution<_RealType>::
01120       __generate_impl(_OutputIterator __f, _OutputIterator __t,
01121                       _UniformRandomNumberGenerator& __urng,
01122                       const param_type& __p)
01123       {
01124         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
01125 
01126         result_type __2q = result_type(2) * __p.q();
01127         result_type __q2 = __p.q() * __p.q();
01128         result_type __q2p1 = result_type(1) + __q2;
01129         result_type __num = result_type(0.5L) * __q2p1;
01130         result_type __omega = __p.omega();
01131         typename __gnu_cxx::arcsine_distribution<result_type>::param_type
01132           __pa(__num, __num / __q2);
01133         while (__f != __t)
01134           {
01135             result_type __x = this->_M_ad(__pa, __urng);
01136             result_type __y = this->_M_ed(__urng);
01137             *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y);
01138           }
01139       }
01140 
01141   template<typename _RealType, typename _CharT, typename _Traits>
01142     std::basic_ostream<_CharT, _Traits>&
01143     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01144                const hoyt_distribution<_RealType>& __x)
01145     {
01146       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01147       typedef typename __ostream_type::ios_base    __ios_base;
01148 
01149       const typename __ios_base::fmtflags __flags = __os.flags();
01150       const _CharT __fill = __os.fill();
01151       const std::streamsize __precision = __os.precision();
01152       const _CharT __space = __os.widen(' ');
01153       __os.flags(__ios_base::scientific | __ios_base::left);
01154       __os.fill(__space);
01155       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01156 
01157       __os << __x.q() << __space << __x.omega();
01158       __os << __space << __x._M_ad;
01159       __os << __space << __x._M_ed;
01160 
01161       __os.flags(__flags);
01162       __os.fill(__fill);
01163       __os.precision(__precision);
01164       return __os;
01165     }
01166 
01167   template<typename _RealType, typename _CharT, typename _Traits>
01168     std::basic_istream<_CharT, _Traits>&
01169     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01170                hoyt_distribution<_RealType>& __x)
01171     {
01172       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01173       typedef typename __istream_type::ios_base    __ios_base;
01174 
01175       const typename __ios_base::fmtflags __flags = __is.flags();
01176       __is.flags(__ios_base::dec | __ios_base::skipws);
01177 
01178       _RealType __q, __omega;
01179       __is >> __q >> __omega;
01180       __is >> __x._M_ad;
01181       __is >> __x._M_ed;
01182       __x.param(typename hoyt_distribution<_RealType>::
01183                 param_type(__q, __omega));
01184 
01185       __is.flags(__flags);
01186       return __is;
01187     }
01188 
01189 
01190   template<typename _RealType>
01191     template<typename _OutputIterator,
01192              typename _UniformRandomNumberGenerator>
01193       void
01194       triangular_distribution<_RealType>::
01195       __generate_impl(_OutputIterator __f, _OutputIterator __t,
01196                       _UniformRandomNumberGenerator& __urng,
01197                       const param_type& __param)
01198       {
01199         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
01200 
01201         while (__f != __t)
01202           *__f++ = this->operator()(__urng, __param);
01203       }
01204 
01205   template<typename _RealType, typename _CharT, typename _Traits>
01206     std::basic_ostream<_CharT, _Traits>&
01207     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01208                const __gnu_cxx::triangular_distribution<_RealType>& __x)
01209     {
01210       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01211       typedef typename __ostream_type::ios_base    __ios_base;
01212 
01213       const typename __ios_base::fmtflags __flags = __os.flags();
01214       const _CharT __fill = __os.fill();
01215       const std::streamsize __precision = __os.precision();
01216       const _CharT __space = __os.widen(' ');
01217       __os.flags(__ios_base::scientific | __ios_base::left);
01218       __os.fill(__space);
01219       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01220 
01221       __os << __x.a() << __space << __x.b() << __space << __x.c();
01222 
01223       __os.flags(__flags);
01224       __os.fill(__fill);
01225       __os.precision(__precision);
01226       return __os;
01227     }
01228 
01229   template<typename _RealType, typename _CharT, typename _Traits>
01230     std::basic_istream<_CharT, _Traits>&
01231     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01232                __gnu_cxx::triangular_distribution<_RealType>& __x)
01233     {
01234       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01235       typedef typename __istream_type::ios_base    __ios_base;
01236 
01237       const typename __ios_base::fmtflags __flags = __is.flags();
01238       __is.flags(__ios_base::dec | __ios_base::skipws);
01239 
01240       _RealType __a, __b, __c;
01241       __is >> __a >> __b >> __c;
01242       __x.param(typename __gnu_cxx::triangular_distribution<_RealType>::
01243                 param_type(__a, __b, __c));
01244 
01245       __is.flags(__flags);
01246       return __is;
01247     }
01248 
01249 
01250   template<typename _RealType>
01251     template<typename _UniformRandomNumberGenerator>
01252       typename von_mises_distribution<_RealType>::result_type
01253       von_mises_distribution<_RealType>::
01254       operator()(_UniformRandomNumberGenerator& __urng,
01255                  const param_type& __p)
01256       {
01257         const result_type __pi
01258           = __gnu_cxx::__math_constants<result_type>::__pi;
01259         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01260           __aurng(__urng);
01261 
01262         result_type __f;
01263         while (1)
01264           {
01265             result_type __rnd = std::cos(__pi * __aurng());
01266             __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd);
01267             result_type __c = __p._M_kappa * (__p._M_r - __f);
01268 
01269             result_type __rnd2 = __aurng();
01270             if (__c * (result_type(2) - __c) > __rnd2)
01271               break;
01272             if (std::log(__c / __rnd2) >= __c - result_type(1))
01273               break;
01274           }
01275 
01276         result_type __res = std::acos(__f);
01277 #if _GLIBCXX_USE_C99_MATH_TR1
01278         __res = std::copysign(__res, __aurng() - result_type(0.5));
01279 #else
01280         if (__aurng() < result_type(0.5))
01281           __res = -__res;
01282 #endif
01283         __res += __p._M_mu;
01284         if (__res > __pi)
01285           __res -= result_type(2) * __pi;
01286         else if (__res < -__pi)
01287           __res += result_type(2) * __pi;
01288         return __res;
01289       }
01290 
01291   template<typename _RealType>
01292     template<typename _OutputIterator,
01293              typename _UniformRandomNumberGenerator>
01294       void
01295       von_mises_distribution<_RealType>::
01296       __generate_impl(_OutputIterator __f, _OutputIterator __t,
01297                       _UniformRandomNumberGenerator& __urng,
01298                       const param_type& __param)
01299       {
01300         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
01301 
01302         while (__f != __t)
01303           *__f++ = this->operator()(__urng, __param);
01304       }
01305 
01306   template<typename _RealType, typename _CharT, typename _Traits>
01307     std::basic_ostream<_CharT, _Traits>&
01308     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01309                const __gnu_cxx::von_mises_distribution<_RealType>& __x)
01310     {
01311       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01312       typedef typename __ostream_type::ios_base    __ios_base;
01313 
01314       const typename __ios_base::fmtflags __flags = __os.flags();
01315       const _CharT __fill = __os.fill();
01316       const std::streamsize __precision = __os.precision();
01317       const _CharT __space = __os.widen(' ');
01318       __os.flags(__ios_base::scientific | __ios_base::left);
01319       __os.fill(__space);
01320       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01321 
01322       __os << __x.mu() << __space << __x.kappa();
01323 
01324       __os.flags(__flags);
01325       __os.fill(__fill);
01326       __os.precision(__precision);
01327       return __os;
01328     }
01329 
01330   template<typename _RealType, typename _CharT, typename _Traits>
01331     std::basic_istream<_CharT, _Traits>&
01332     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01333                __gnu_cxx::von_mises_distribution<_RealType>& __x)
01334     {
01335       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01336       typedef typename __istream_type::ios_base    __ios_base;
01337 
01338       const typename __ios_base::fmtflags __flags = __is.flags();
01339       __is.flags(__ios_base::dec | __ios_base::skipws);
01340 
01341       _RealType __mu, __kappa;
01342       __is >> __mu >> __kappa;
01343       __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>::
01344                 param_type(__mu, __kappa));
01345 
01346       __is.flags(__flags);
01347       return __is;
01348     }
01349 
01350 
01351   template<typename _UIntType>
01352     template<typename _UniformRandomNumberGenerator>
01353       typename hypergeometric_distribution<_UIntType>::result_type
01354       hypergeometric_distribution<_UIntType>::
01355       operator()(_UniformRandomNumberGenerator& __urng,
01356                  const param_type& __param)
01357       {
01358         std::__detail::_Adaptor<_UniformRandomNumberGenerator, double>
01359           __aurng(__urng);
01360 
01361         result_type __a = __param.successful_size();
01362         result_type __b = __param.total_size();
01363         result_type __k = 0;
01364 
01365         if (__param.total_draws() < __param.total_size() / 2)
01366           {
01367             for (result_type __i = 0; __i < __param.total_draws(); ++__i)
01368               {
01369                 if (__b * __aurng() < __a)
01370                   {
01371                     ++__k;
01372                     if (__k == __param.successful_size())
01373                       return __k;
01374                    --__a;
01375                   }
01376                 --__b;
01377               }
01378             return __k;
01379           }
01380         else
01381           {
01382             for (result_type __i = 0; __i < __param.unsuccessful_size(); ++__i)
01383               {
01384                 if (__b * __aurng() < __a)
01385                   {
01386                     ++__k;
01387                     if (__k == __param.successful_size())
01388                       return __param.successful_size() - __k;
01389                     --__a;
01390                   }
01391                 --__b;
01392               }
01393             return __param.successful_size() - __k;
01394           }
01395       }
01396 
01397   template<typename _UIntType>
01398     template<typename _OutputIterator,
01399              typename _UniformRandomNumberGenerator>
01400       void
01401       hypergeometric_distribution<_UIntType>::
01402       __generate_impl(_OutputIterator __f, _OutputIterator __t,
01403                       _UniformRandomNumberGenerator& __urng,
01404                       const param_type& __param)
01405       {
01406         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
01407 
01408         while (__f != __t)
01409           *__f++ = this->operator()(__urng);
01410       }
01411 
01412   template<typename _UIntType, typename _CharT, typename _Traits>
01413     std::basic_ostream<_CharT, _Traits>&
01414     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01415                const __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
01416     {
01417       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01418       typedef typename __ostream_type::ios_base    __ios_base;
01419 
01420       const typename __ios_base::fmtflags __flags = __os.flags();
01421       const _CharT __fill = __os.fill();
01422       const std::streamsize __precision = __os.precision();
01423       const _CharT __space = __os.widen(' ');
01424       __os.flags(__ios_base::scientific | __ios_base::left);
01425       __os.fill(__space);
01426       __os.precision(std::numeric_limits<_UIntType>::max_digits10);
01427 
01428       __os << __x.total_size() << __space << __x.successful_size() << __space
01429            << __x.total_draws();
01430 
01431       __os.flags(__flags);
01432       __os.fill(__fill);
01433       __os.precision(__precision);
01434       return __os;
01435     }
01436 
01437   template<typename _UIntType, typename _CharT, typename _Traits>
01438     std::basic_istream<_CharT, _Traits>&
01439     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01440                __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
01441     {
01442       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01443       typedef typename __istream_type::ios_base    __ios_base;
01444 
01445       const typename __ios_base::fmtflags __flags = __is.flags();
01446       __is.flags(__ios_base::dec | __ios_base::skipws);
01447 
01448       _UIntType __total_size, __successful_size, __total_draws;
01449       __is >> __total_size >> __successful_size >> __total_draws;
01450       __x.param(typename __gnu_cxx::hypergeometric_distribution<_UIntType>::
01451                 param_type(__total_size, __successful_size, __total_draws));
01452 
01453       __is.flags(__flags);
01454       return __is;
01455     }
01456 
01457 
01458   template<typename _RealType>
01459     template<typename _UniformRandomNumberGenerator>
01460       typename logistic_distribution<_RealType>::result_type
01461       logistic_distribution<_RealType>::
01462       operator()(_UniformRandomNumberGenerator& __urng,
01463                  const param_type& __p)
01464       {
01465         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01466           __aurng(__urng);
01467 
01468         result_type __arg = result_type(1);
01469         while (__arg == result_type(1) || __arg == result_type(0))
01470           __arg = __aurng();
01471         return __p.a()
01472              + __p.b() * std::log(__arg / (result_type(1) - __arg));
01473       }
01474 
01475   template<typename _RealType>
01476     template<typename _OutputIterator,
01477              typename _UniformRandomNumberGenerator>
01478       void
01479       logistic_distribution<_RealType>::
01480       __generate_impl(_OutputIterator __f, _OutputIterator __t,
01481                       _UniformRandomNumberGenerator& __urng,
01482                       const param_type& __p)
01483       {
01484         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
01485         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
01486           __aurng(__urng);
01487 
01488         while (__f != __t)
01489           {
01490             result_type __arg = result_type(1);
01491             while (__arg == result_type(1) || __arg == result_type(0))
01492               __arg = __aurng();
01493             *__f++ = __p.a()
01494                    + __p.b() * std::log(__arg / (result_type(1) - __arg));
01495           }
01496       }
01497 
01498   template<typename _RealType, typename _CharT, typename _Traits>
01499     std::basic_ostream<_CharT, _Traits>&
01500     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01501                const logistic_distribution<_RealType>& __x)
01502     {
01503       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
01504       typedef typename __ostream_type::ios_base    __ios_base;
01505 
01506       const typename __ios_base::fmtflags __flags = __os.flags();
01507       const _CharT __fill = __os.fill();
01508       const std::streamsize __precision = __os.precision();
01509       const _CharT __space = __os.widen(' ');
01510       __os.flags(__ios_base::scientific | __ios_base::left);
01511       __os.fill(__space);
01512       __os.precision(std::numeric_limits<_RealType>::max_digits10);
01513 
01514       __os << __x.a() << __space << __x.b();
01515 
01516       __os.flags(__flags);
01517       __os.fill(__fill);
01518       __os.precision(__precision);
01519       return __os;
01520     }
01521 
01522   template<typename _RealType, typename _CharT, typename _Traits>
01523     std::basic_istream<_CharT, _Traits>&
01524     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01525                logistic_distribution<_RealType>& __x)
01526     {
01527       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
01528       typedef typename __istream_type::ios_base    __ios_base;
01529 
01530       const typename __ios_base::fmtflags __flags = __is.flags();
01531       __is.flags(__ios_base::dec | __ios_base::skipws);
01532 
01533       _RealType __a, __b;
01534       __is >> __a >> __b;
01535       __x.param(typename logistic_distribution<_RealType>::
01536                 param_type(__a, __b));
01537 
01538       __is.flags(__flags);
01539       return __is;
01540     }
01541 
01542 
01543   namespace {
01544 
01545     // Helper class for the uniform_on_sphere_distribution generation
01546     // function.
01547     template<std::size_t _Dimen, typename _RealType>
01548       class uniform_on_sphere_helper
01549       {
01550         typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>::
01551           result_type result_type;
01552 
01553       public:
01554         template<typename _NormalDistribution,
01555                  typename _UniformRandomNumberGenerator>
01556         result_type operator()(_NormalDistribution& __nd,
01557                                _UniformRandomNumberGenerator& __urng)
01558         {
01559           result_type __ret;
01560           typename result_type::value_type __norm;
01561 
01562           do
01563             {
01564               auto __sum = _RealType(0);
01565 
01566               std::generate(__ret.begin(), __ret.end(),
01567                             [&__nd, &__urng, &__sum](){
01568                               _RealType __t = __nd(__urng);
01569                               __sum += __t * __t;
01570                               return __t; });
01571               __norm = std::sqrt(__sum);
01572             }
01573           while (__norm == _RealType(0) || ! std::isfinite(__norm));
01574 
01575           std::transform(__ret.begin(), __ret.end(), __ret.begin(),
01576                          [__norm](_RealType __val){ return __val / __norm; });
01577 
01578           return __ret;
01579         }
01580       };
01581 
01582 
01583     template<typename _RealType>
01584       class uniform_on_sphere_helper<2, _RealType>
01585       {
01586         typedef typename uniform_on_sphere_distribution<2, _RealType>::
01587           result_type result_type;
01588 
01589       public:
01590         template<typename _NormalDistribution,
01591                  typename _UniformRandomNumberGenerator>
01592         result_type operator()(_NormalDistribution&,
01593                                _UniformRandomNumberGenerator& __urng)
01594         {
01595           result_type __ret;
01596           _RealType __sq;
01597           std::__detail::_Adaptor<_UniformRandomNumberGenerator,
01598                                   _RealType> __aurng(__urng);
01599 
01600           do
01601             {
01602               __ret[0] = _RealType(2) * __aurng() - _RealType(1);
01603               __ret[1] = _RealType(2) * __aurng() - _RealType(1);
01604 
01605               __sq = __ret[0] * __ret[0] + __ret[1] * __ret[1];
01606             }
01607           while (__sq == _RealType(0) || __sq > _RealType(1));
01608 
01609 #if _GLIBCXX_USE_C99_MATH_TR1
01610           // Yes, we do not just use sqrt(__sq) because hypot() is more
01611           // accurate.
01612           auto __norm = std::hypot(__ret[0], __ret[1]);
01613 #else
01614           auto __norm = std::sqrt(__sq);
01615 #endif
01616           __ret[0] /= __norm;
01617           __ret[1] /= __norm;
01618 
01619           return __ret;
01620         }
01621       };
01622 
01623   }
01624 
01625 
01626   template<std::size_t _Dimen, typename _RealType>
01627     template<typename _UniformRandomNumberGenerator>
01628       typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type
01629       uniform_on_sphere_distribution<_Dimen, _RealType>::
01630       operator()(_UniformRandomNumberGenerator& __urng,
01631                  const param_type& __p)
01632       {
01633         uniform_on_sphere_helper<_Dimen, _RealType> __helper;
01634         return __helper(_M_nd, __urng);
01635       }
01636 
01637   template<std::size_t _Dimen, typename _RealType>
01638     template<typename _OutputIterator,
01639              typename _UniformRandomNumberGenerator>
01640       void
01641       uniform_on_sphere_distribution<_Dimen, _RealType>::
01642       __generate_impl(_OutputIterator __f, _OutputIterator __t,
01643                       _UniformRandomNumberGenerator& __urng,
01644                       const param_type& __param)
01645       {
01646         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
01647 
01648         while (__f != __t)
01649           *__f++ = this->operator()(__urng, __param);
01650       }
01651 
01652   template<std::size_t _Dimen, typename _RealType, typename _CharT,
01653            typename _Traits>
01654     std::basic_ostream<_CharT, _Traits>&
01655     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
01656                const __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
01657                                                                _RealType>& __x)
01658     {
01659       return __os << __x._M_nd;
01660     }
01661 
01662   template<std::size_t _Dimen, typename _RealType, typename _CharT,
01663            typename _Traits>
01664     std::basic_istream<_CharT, _Traits>&
01665     operator>>(std::basic_istream<_CharT, _Traits>& __is,
01666                __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
01667                                                          _RealType>& __x)
01668     {
01669       return __is >> __x._M_nd;
01670     }
01671 
01672 _GLIBCXX_END_NAMESPACE_VERSION
01673 } // namespace
01674 
01675 
01676 #endif // _EXT_RANDOM_TCC