Home | History | Annotate | Download | only in ext
      1 // Random number extensions -*- C++ -*-
      2 
      3 // Copyright (C) 2012-2013 Free Software Foundation, Inc.
      4 //
      5 // This file is part of the GNU ISO C++ Library.  This library is free
      6 // software; you can redistribute it and/or modify it under the
      7 // terms of the GNU General Public License as published by the
      8 // Free Software Foundation; either version 3, or (at your option)
      9 // any later version.
     10 
     11 // This library is distributed in the hope that it will be useful,
     12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
     13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     14 // GNU General Public License for more details.
     15 
     16 // Under Section 7 of GPL version 3, you are granted additional
     17 // permissions described in the GCC Runtime Library Exception, version
     18 // 3.1, as published by the Free Software Foundation.
     19 
     20 // You should have received a copy of the GNU General Public License and
     21 // a copy of the GCC Runtime Library Exception along with this program;
     22 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
     23 // <http://www.gnu.org/licenses/>.
     24 
     25 /** @file ext/random.tcc
     26  *  This is an internal header file, included by other library headers.
     27  *  Do not attempt to use it directly. @headername{ext/random}
     28  */
     29 
     30 #ifndef _EXT_RANDOM_TCC
     31 #define _EXT_RANDOM_TCC 1
     32 
     33 #pragma GCC system_header
     34 
     35 
     36 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
     37 {
     38 _GLIBCXX_BEGIN_NAMESPACE_VERSION
     39 
     40 #if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
     41 
     42   template<typename _UIntType, size_t __m,
     43 	   size_t __pos1, size_t __sl1, size_t __sl2,
     44 	   size_t __sr1, size_t __sr2,
     45 	   uint32_t __msk1, uint32_t __msk2,
     46 	   uint32_t __msk3, uint32_t __msk4,
     47 	   uint32_t __parity1, uint32_t __parity2,
     48 	   uint32_t __parity3, uint32_t __parity4>
     49     void simd_fast_mersenne_twister_engine<_UIntType, __m,
     50 					   __pos1, __sl1, __sl2, __sr1, __sr2,
     51 					   __msk1, __msk2, __msk3, __msk4,
     52 					   __parity1, __parity2, __parity3,
     53 					   __parity4>::
     54     seed(_UIntType __seed)
     55     {
     56       _M_state32[0] = static_cast<uint32_t>(__seed);
     57       for (size_t __i = 1; __i < _M_nstate32; ++__i)
     58 	_M_state32[__i] = (1812433253UL
     59 			   * (_M_state32[__i - 1] ^ (_M_state32[__i - 1] >> 30))
     60 			   + __i);
     61       _M_pos = state_size;
     62       _M_period_certification();
     63     }
     64 
     65 
     66   namespace {
     67 
     68     inline uint32_t _Func1(uint32_t __x)
     69     {
     70       return (__x ^ (__x >> 27)) * UINT32_C(1664525);
     71     }
     72 
     73     inline uint32_t _Func2(uint32_t __x)
     74     {
     75       return (__x ^ (__x >> 27)) * UINT32_C(1566083941);
     76     }
     77 
     78   }
     79 
     80 
     81   template<typename _UIntType, size_t __m,
     82 	   size_t __pos1, size_t __sl1, size_t __sl2,
     83 	   size_t __sr1, size_t __sr2,
     84 	   uint32_t __msk1, uint32_t __msk2,
     85 	   uint32_t __msk3, uint32_t __msk4,
     86 	   uint32_t __parity1, uint32_t __parity2,
     87 	   uint32_t __parity3, uint32_t __parity4>
     88     template<typename _Sseq>
     89       typename std::enable_if<std::is_class<_Sseq>::value>::type
     90       simd_fast_mersenne_twister_engine<_UIntType, __m,
     91 					__pos1, __sl1, __sl2, __sr1, __sr2,
     92 					__msk1, __msk2, __msk3, __msk4,
     93 					__parity1, __parity2, __parity3,
     94 					__parity4>::
     95       seed(_Sseq& __q)
     96       {
     97 	size_t __lag;
     98 
     99 	if (_M_nstate32 >= 623)
    100 	  __lag = 11;
    101 	else if (_M_nstate32 >= 68)
    102 	  __lag = 7;
    103 	else if (_M_nstate32 >= 39)
    104 	  __lag = 5;
    105 	else
    106 	  __lag = 3;
    107 	const size_t __mid = (_M_nstate32 - __lag) / 2;
    108 
    109 	std::fill(_M_state32, _M_state32 + _M_nstate32, UINT32_C(0x8b8b8b8b));
    110 	uint32_t __arr[_M_nstate32];
    111 	__q.generate(__arr + 0, __arr + _M_nstate32);
    112 
    113 	uint32_t __r = _Func1(_M_state32[0] ^ _M_state32[__mid]
    114 			      ^ _M_state32[_M_nstate32  - 1]);
    115 	_M_state32[__mid] += __r;
    116 	__r += _M_nstate32;
    117 	_M_state32[__mid + __lag] += __r;
    118 	_M_state32[0] = __r;
    119 
    120 	for (size_t __i = 1, __j = 0; __j < _M_nstate32; ++__j)
    121 	  {
    122 	    __r = _Func1(_M_state32[__i]
    123 			 ^ _M_state32[(__i + __mid) % _M_nstate32]
    124 			 ^ _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
    125 	    _M_state32[(__i + __mid) % _M_nstate32] += __r;
    126 	    __r += __arr[__j] + __i;
    127 	    _M_state32[(__i + __mid + __lag) % _M_nstate32] += __r;
    128 	    _M_state32[__i] = __r;
    129 	    __i = (__i + 1) % _M_nstate32;
    130 	  }
    131 	for (size_t __j = 0; __j < _M_nstate32; ++__j)
    132 	  {
    133 	    const size_t __i = (__j + 1) % _M_nstate32;
    134 	    __r = _Func2(_M_state32[__i]
    135 			 + _M_state32[(__i + __mid) % _M_nstate32]
    136 			 + _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
    137 	    _M_state32[(__i + __mid) % _M_nstate32] ^= __r;
    138 	    __r -= __i;
    139 	    _M_state32[(__i + __mid + __lag) % _M_nstate32] ^= __r;
    140 	    _M_state32[__i] = __r;
    141 	  }
    142 
    143 	_M_pos = state_size;
    144 	_M_period_certification();
    145       }
    146 
    147 
    148   template<typename _UIntType, size_t __m,
    149 	   size_t __pos1, size_t __sl1, size_t __sl2,
    150 	   size_t __sr1, size_t __sr2,
    151 	   uint32_t __msk1, uint32_t __msk2,
    152 	   uint32_t __msk3, uint32_t __msk4,
    153 	   uint32_t __parity1, uint32_t __parity2,
    154 	   uint32_t __parity3, uint32_t __parity4>
    155     void simd_fast_mersenne_twister_engine<_UIntType, __m,
    156 					   __pos1, __sl1, __sl2, __sr1, __sr2,
    157 					   __msk1, __msk2, __msk3, __msk4,
    158 					   __parity1, __parity2, __parity3,
    159 					   __parity4>::
    160     _M_period_certification(void)
    161     {
    162       static const uint32_t __parity[4] = { __parity1, __parity2,
    163 					    __parity3, __parity4 };
    164       uint32_t __inner = 0;
    165       for (size_t __i = 0; __i < 4; ++__i)
    166 	if (__parity[__i] != 0)
    167 	  __inner ^= _M_state32[__i] & __parity[__i];
    168 
    169       if (__builtin_parity(__inner) & 1)
    170 	return;
    171       for (size_t __i = 0; __i < 4; ++__i)
    172 	if (__parity[__i] != 0)
    173 	  {
    174 	    _M_state32[__i] ^= 1 << (__builtin_ffs(__parity[__i]) - 1);
    175 	    return;
    176 	  }
    177       __builtin_unreachable();
    178     }
    179 
    180 
    181   template<typename _UIntType, size_t __m,
    182 	   size_t __pos1, size_t __sl1, size_t __sl2,
    183 	   size_t __sr1, size_t __sr2,
    184 	   uint32_t __msk1, uint32_t __msk2,
    185 	   uint32_t __msk3, uint32_t __msk4,
    186 	   uint32_t __parity1, uint32_t __parity2,
    187 	   uint32_t __parity3, uint32_t __parity4>
    188     void simd_fast_mersenne_twister_engine<_UIntType, __m,
    189 					   __pos1, __sl1, __sl2, __sr1, __sr2,
    190 					   __msk1, __msk2, __msk3, __msk4,
    191 					   __parity1, __parity2, __parity3,
    192 					   __parity4>::
    193     discard(unsigned long long __z)
    194     {
    195       while (__z > state_size - _M_pos)
    196 	{
    197 	  __z -= state_size - _M_pos;
    198 
    199 	  _M_gen_rand();
    200 	}
    201 
    202       _M_pos += __z;
    203     }
    204 
    205 
    206 #ifndef  _GLIBCXX_OPT_HAVE_RANDOM_SFMT_GEN_READ
    207 
    208   namespace {
    209 
    210     template<size_t __shift>
    211       inline void __rshift(uint32_t *__out, const uint32_t *__in)
    212       {
    213 	uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
    214 			 | static_cast<uint64_t>(__in[2]));
    215 	uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
    216 			 | static_cast<uint64_t>(__in[0]));
    217 
    218 	uint64_t __oh = __th >> (__shift * 8);
    219 	uint64_t __ol = __tl >> (__shift * 8);
    220 	__ol |= __th << (64 - __shift * 8);
    221 	__out[1] = static_cast<uint32_t>(__ol >> 32);
    222 	__out[0] = static_cast<uint32_t>(__ol);
    223 	__out[3] = static_cast<uint32_t>(__oh >> 32);
    224 	__out[2] = static_cast<uint32_t>(__oh);
    225       }
    226 
    227 
    228     template<size_t __shift>
    229       inline void __lshift(uint32_t *__out, const uint32_t *__in)
    230       {
    231 	uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
    232 			 | static_cast<uint64_t>(__in[2]));
    233 	uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
    234 			 | static_cast<uint64_t>(__in[0]));
    235 
    236 	uint64_t __oh = __th << (__shift * 8);
    237 	uint64_t __ol = __tl << (__shift * 8);
    238 	__oh |= __tl >> (64 - __shift * 8);
    239 	__out[1] = static_cast<uint32_t>(__ol >> 32);
    240 	__out[0] = static_cast<uint32_t>(__ol);
    241 	__out[3] = static_cast<uint32_t>(__oh >> 32);
    242 	__out[2] = static_cast<uint32_t>(__oh);
    243       }
    244 
    245 
    246     template<size_t __sl1, size_t __sl2, size_t __sr1, size_t __sr2,
    247 	     uint32_t __msk1, uint32_t __msk2, uint32_t __msk3, uint32_t __msk4>
    248       inline void __recursion(uint32_t *__r,
    249 			      const uint32_t *__a, const uint32_t *__b,
    250 			      const uint32_t *__c, const uint32_t *__d)
    251       {
    252 	uint32_t __x[4];
    253 	uint32_t __y[4];
    254 
    255 	__lshift<__sl2>(__x, __a);
    256 	__rshift<__sr2>(__y, __c);
    257 	__r[0] = (__a[0] ^ __x[0] ^ ((__b[0] >> __sr1) & __msk1)
    258 		  ^ __y[0] ^ (__d[0] << __sl1));
    259 	__r[1] = (__a[1] ^ __x[1] ^ ((__b[1] >> __sr1) & __msk2)
    260 		  ^ __y[1] ^ (__d[1] << __sl1));
    261 	__r[2] = (__a[2] ^ __x[2] ^ ((__b[2] >> __sr1) & __msk3)
    262 		  ^ __y[2] ^ (__d[2] << __sl1));
    263 	__r[3] = (__a[3] ^ __x[3] ^ ((__b[3] >> __sr1) & __msk4)
    264 		  ^ __y[3] ^ (__d[3] << __sl1));
    265       }
    266 
    267   }
    268 
    269 
    270   template<typename _UIntType, size_t __m,
    271 	   size_t __pos1, size_t __sl1, size_t __sl2,
    272 	   size_t __sr1, size_t __sr2,
    273 	   uint32_t __msk1, uint32_t __msk2,
    274 	   uint32_t __msk3, uint32_t __msk4,
    275 	   uint32_t __parity1, uint32_t __parity2,
    276 	   uint32_t __parity3, uint32_t __parity4>
    277     void simd_fast_mersenne_twister_engine<_UIntType, __m,
    278 					   __pos1, __sl1, __sl2, __sr1, __sr2,
    279 					   __msk1, __msk2, __msk3, __msk4,
    280 					   __parity1, __parity2, __parity3,
    281 					   __parity4>::
    282     _M_gen_rand(void)
    283     {
    284       const uint32_t *__r1 = &_M_state32[_M_nstate32 - 8];
    285       const uint32_t *__r2 = &_M_state32[_M_nstate32 - 4];
    286       static constexpr size_t __pos1_32 = __pos1 * 4;
    287 
    288       size_t __i;
    289       for (__i = 0; __i < _M_nstate32 - __pos1_32; __i += 4)
    290 	{
    291 	  __recursion<__sl1, __sl2, __sr1, __sr2,
    292 		      __msk1, __msk2, __msk3, __msk4>
    293 	    (&_M_state32[__i], &_M_state32[__i],
    294 	     &_M_state32[__i + __pos1_32], __r1, __r2);
    295 	  __r1 = __r2;
    296 	  __r2 = &_M_state32[__i];
    297 	}
    298 
    299       for (; __i < _M_nstate32; __i += 4)
    300 	{
    301 	  __recursion<__sl1, __sl2, __sr1, __sr2,
    302 		      __msk1, __msk2, __msk3, __msk4>
    303 	    (&_M_state32[__i], &_M_state32[__i],
    304 	     &_M_state32[__i + __pos1_32 - _M_nstate32], __r1, __r2);
    305 	  __r1 = __r2;
    306 	  __r2 = &_M_state32[__i];
    307 	}
    308 
    309       _M_pos = 0;
    310     }
    311 
    312 #endif
    313 
    314 #ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_OPERATOREQUAL
    315   template<typename _UIntType, size_t __m,
    316 	   size_t __pos1, size_t __sl1, size_t __sl2,
    317 	   size_t __sr1, size_t __sr2,
    318 	   uint32_t __msk1, uint32_t __msk2,
    319 	   uint32_t __msk3, uint32_t __msk4,
    320 	   uint32_t __parity1, uint32_t __parity2,
    321 	   uint32_t __parity3, uint32_t __parity4>
    322     bool
    323     operator==(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
    324 	       __m, __pos1, __sl1, __sl2, __sr1, __sr2,
    325 	       __msk1, __msk2, __msk3, __msk4,
    326 	       __parity1, __parity2, __parity3, __parity4>& __lhs,
    327 	       const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
    328 	       __m, __pos1, __sl1, __sl2, __sr1, __sr2,
    329 	       __msk1, __msk2, __msk3, __msk4,
    330 	       __parity1, __parity2, __parity3, __parity4>& __rhs)
    331     {
    332       typedef __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
    333 	       __m, __pos1, __sl1, __sl2, __sr1, __sr2,
    334 	       __msk1, __msk2, __msk3, __msk4,
    335 	       __parity1, __parity2, __parity3, __parity4> __engine;
    336       return (std::equal(__lhs._M_stateT,
    337 			 __lhs._M_stateT + __engine::state_size,
    338 			 __rhs._M_stateT)
    339 	      && __lhs._M_pos == __rhs._M_pos);
    340     }
    341 #endif
    342 
    343   template<typename _UIntType, size_t __m,
    344 	   size_t __pos1, size_t __sl1, size_t __sl2,
    345 	   size_t __sr1, size_t __sr2,
    346 	   uint32_t __msk1, uint32_t __msk2,
    347 	   uint32_t __msk3, uint32_t __msk4,
    348 	   uint32_t __parity1, uint32_t __parity2,
    349 	   uint32_t __parity3, uint32_t __parity4,
    350 	   typename _CharT, typename _Traits>
    351     std::basic_ostream<_CharT, _Traits>&
    352     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    353 	       const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
    354 	       __m, __pos1, __sl1, __sl2, __sr1, __sr2,
    355 	       __msk1, __msk2, __msk3, __msk4,
    356 	       __parity1, __parity2, __parity3, __parity4>& __x)
    357     {
    358       typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
    359       typedef typename __ostream_type::ios_base __ios_base;
    360 
    361       const typename __ios_base::fmtflags __flags = __os.flags();
    362       const _CharT __fill = __os.fill();
    363       const _CharT __space = __os.widen(' ');
    364       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
    365       __os.fill(__space);
    366 
    367       for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
    368 	__os << __x._M_state32[__i] << __space;
    369       __os << __x._M_pos;
    370 
    371       __os.flags(__flags);
    372       __os.fill(__fill);
    373       return __os;
    374     }
    375 
    376 
    377   template<typename _UIntType, size_t __m,
    378 	   size_t __pos1, size_t __sl1, size_t __sl2,
    379 	   size_t __sr1, size_t __sr2,
    380 	   uint32_t __msk1, uint32_t __msk2,
    381 	   uint32_t __msk3, uint32_t __msk4,
    382 	   uint32_t __parity1, uint32_t __parity2,
    383 	   uint32_t __parity3, uint32_t __parity4,
    384 	   typename _CharT, typename _Traits>
    385     std::basic_istream<_CharT, _Traits>&
    386     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    387 	       __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
    388 	       __m, __pos1, __sl1, __sl2, __sr1, __sr2,
    389 	       __msk1, __msk2, __msk3, __msk4,
    390 	       __parity1, __parity2, __parity3, __parity4>& __x)
    391     {
    392       typedef std::basic_istream<_CharT, _Traits> __istream_type;
    393       typedef typename __istream_type::ios_base __ios_base;
    394 
    395       const typename __ios_base::fmtflags __flags = __is.flags();
    396       __is.flags(__ios_base::dec | __ios_base::skipws);
    397 
    398       for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
    399 	__is >> __x._M_state32[__i];
    400       __is >> __x._M_pos;
    401 
    402       __is.flags(__flags);
    403       return __is;
    404     }
    405 
    406 #endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
    407 
    408   /**
    409    * Iteration method due to M.D. J<o:>hnk.
    410    *
    411    * M.D. J<o:>hnk, Erzeugung von betaverteilten und gammaverteilten
    412    * Zufallszahlen, Metrika, Volume 8, 1964
    413    */
    414   template<typename _RealType>
    415     template<typename _UniformRandomNumberGenerator>
    416       typename beta_distribution<_RealType>::result_type
    417       beta_distribution<_RealType>::
    418       operator()(_UniformRandomNumberGenerator& __urng,
    419 		 const param_type& __param)
    420       {
    421 	std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
    422 	  __aurng(__urng);
    423 
    424 	result_type __x, __y;
    425 	do
    426 	  {
    427 	    __x = std::exp(std::log(__aurng()) / __param.alpha());
    428 	    __y = std::exp(std::log(__aurng()) / __param.beta());
    429 	  }
    430 	while (__x + __y > result_type(1));
    431 
    432 	return __x / (__x + __y);
    433       }
    434 
    435   template<typename _RealType>
    436     template<typename _OutputIterator,
    437 	     typename _UniformRandomNumberGenerator>
    438       void
    439       beta_distribution<_RealType>::
    440       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    441 		      _UniformRandomNumberGenerator& __urng,
    442 		      const param_type& __param)
    443       {
    444 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
    445 
    446 	std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
    447 	  __aurng(__urng);
    448 
    449 	while (__f != __t)
    450 	  {
    451 	    result_type __x, __y;
    452 	    do
    453 	      {
    454 		__x = std::exp(std::log(__aurng()) / __param.alpha());
    455 		__y = std::exp(std::log(__aurng()) / __param.beta());
    456 	      }
    457 	    while (__x + __y > result_type(1));
    458 
    459 	    *__f++ = __x / (__x + __y);
    460 	  }
    461       }
    462 
    463   template<typename _RealType, typename _CharT, typename _Traits>
    464     std::basic_ostream<_CharT, _Traits>&
    465     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    466 	       const __gnu_cxx::beta_distribution<_RealType>& __x)
    467     {
    468       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    469       typedef typename __ostream_type::ios_base    __ios_base;
    470 
    471       const typename __ios_base::fmtflags __flags = __os.flags();
    472       const _CharT __fill = __os.fill();
    473       const std::streamsize __precision = __os.precision();
    474       const _CharT __space = __os.widen(' ');
    475       __os.flags(__ios_base::scientific | __ios_base::left);
    476       __os.fill(__space);
    477       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    478 
    479       __os << __x.alpha() << __space << __x.beta();
    480 
    481       __os.flags(__flags);
    482       __os.fill(__fill);
    483       __os.precision(__precision);
    484       return __os;
    485     }
    486 
    487   template<typename _RealType, typename _CharT, typename _Traits>
    488     std::basic_istream<_CharT, _Traits>&
    489     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    490 	       __gnu_cxx::beta_distribution<_RealType>& __x)
    491     {
    492       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    493       typedef typename __istream_type::ios_base    __ios_base;
    494 
    495       const typename __ios_base::fmtflags __flags = __is.flags();
    496       __is.flags(__ios_base::dec | __ios_base::skipws);
    497 
    498       _RealType __alpha_val, __beta_val;
    499       __is >> __alpha_val >> __beta_val;
    500       __x.param(typename __gnu_cxx::beta_distribution<_RealType>::
    501 		param_type(__alpha_val, __beta_val));
    502 
    503       __is.flags(__flags);
    504       return __is;
    505     }
    506 
    507 
    508   template<std::size_t _Dimen, typename _RealType>
    509     template<typename _InputIterator1, typename _InputIterator2>
    510       void
    511       normal_mv_distribution<_Dimen, _RealType>::param_type::
    512       _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
    513 		   _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
    514       {
    515 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
    516 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
    517 	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
    518 		  _M_mean.end(), _RealType(0));
    519 
    520 	// Perform the Cholesky decomposition
    521 	auto __w = _M_t.begin();
    522 	for (size_t __j = 0; __j < _Dimen; ++__j)
    523 	  {
    524 	    _RealType __sum = _RealType(0);
    525 
    526 	    auto __slitbegin = __w;
    527 	    auto __cit = _M_t.begin();
    528 	    for (size_t __i = 0; __i < __j; ++__i)
    529 	      {
    530 		auto __slit = __slitbegin;
    531 		_RealType __s = *__varcovbegin++;
    532 		for (size_t __k = 0; __k < __i; ++__k)
    533 		  __s -= *__slit++ * *__cit++;
    534 
    535 		*__w++ = __s /= *__cit++;
    536 		__sum += __s * __s;
    537 	      }
    538 
    539 	    __sum = *__varcovbegin - __sum;
    540 	    if (__builtin_expect(__sum <= _RealType(0), 0))
    541 	      std::__throw_runtime_error(__N("normal_mv_distribution::"
    542 					     "param_type::_M_init_full"));
    543 	    *__w++ = std::sqrt(__sum);
    544 
    545 	    std::advance(__varcovbegin, _Dimen - __j);
    546 	  }
    547       }
    548 
    549   template<std::size_t _Dimen, typename _RealType>
    550     template<typename _InputIterator1, typename _InputIterator2>
    551       void
    552       normal_mv_distribution<_Dimen, _RealType>::param_type::
    553       _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
    554 		    _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
    555       {
    556 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
    557 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
    558 	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
    559 		  _M_mean.end(), _RealType(0));
    560 
    561 	// Perform the Cholesky decomposition
    562 	auto __w = _M_t.begin();
    563 	for (size_t __j = 0; __j < _Dimen; ++__j)
    564 	  {
    565 	    _RealType __sum = _RealType(0);
    566 
    567 	    auto __slitbegin = __w;
    568 	    auto __cit = _M_t.begin();
    569 	    for (size_t __i = 0; __i < __j; ++__i)
    570 	      {
    571 		auto __slit = __slitbegin;
    572 		_RealType __s = *__varcovbegin++;
    573 		for (size_t __k = 0; __k < __i; ++__k)
    574 		  __s -= *__slit++ * *__cit++;
    575 
    576 		*__w++ = __s /= *__cit++;
    577 		__sum += __s * __s;
    578 	      }
    579 
    580 	    __sum = *__varcovbegin++ - __sum;
    581 	    if (__builtin_expect(__sum <= _RealType(0), 0))
    582 	      std::__throw_runtime_error(__N("normal_mv_distribution::"
    583 					     "param_type::_M_init_full"));
    584 	    *__w++ = std::sqrt(__sum);
    585 	  }
    586       }
    587 
    588   template<std::size_t _Dimen, typename _RealType>
    589     template<typename _InputIterator1, typename _InputIterator2>
    590       void
    591       normal_mv_distribution<_Dimen, _RealType>::param_type::
    592       _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
    593 		       _InputIterator2 __varbegin, _InputIterator2 __varend)
    594       {
    595 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
    596 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
    597 	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
    598 		  _M_mean.end(), _RealType(0));
    599 
    600 	auto __w = _M_t.begin();
    601 	size_t __step = 0;
    602 	while (__varbegin != __varend)
    603 	  {
    604 	    std::fill_n(__w, __step, _RealType(0));
    605 	    __w += __step++;
    606 	    if (__builtin_expect(*__varbegin < _RealType(0), 0))
    607 	      std::__throw_runtime_error(__N("normal_mv_distribution::"
    608 					     "param_type::_M_init_diagonal"));
    609 	    *__w++ = std::sqrt(*__varbegin++);
    610 	  }
    611       }
    612 
    613   template<std::size_t _Dimen, typename _RealType>
    614     template<typename _UniformRandomNumberGenerator>
    615       typename normal_mv_distribution<_Dimen, _RealType>::result_type
    616       normal_mv_distribution<_Dimen, _RealType>::
    617       operator()(_UniformRandomNumberGenerator& __urng,
    618 		 const param_type& __param)
    619       {
    620 	result_type __ret;
    621 
    622 	_M_nd.__generate(__ret.begin(), __ret.end(), __urng);
    623 
    624 	auto __t_it = __param._M_t.crbegin();
    625 	for (size_t __i = _Dimen; __i > 0; --__i)
    626 	  {
    627 	    _RealType __sum = _RealType(0);
    628 	    for (size_t __j = __i; __j > 0; --__j)
    629 	      __sum += __ret[__j - 1] * *__t_it++;
    630 	    __ret[__i - 1] = __sum;
    631 	  }
    632 
    633 	return __ret;
    634       }
    635 
    636   template<std::size_t _Dimen, typename _RealType>
    637     template<typename _ForwardIterator, typename _UniformRandomNumberGenerator>
    638       void
    639       normal_mv_distribution<_Dimen, _RealType>::
    640       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
    641 		      _UniformRandomNumberGenerator& __urng,
    642 		      const param_type& __param)
    643       {
    644 	__glibcxx_function_requires(_Mutable_ForwardIteratorConcept<
    645 				    _ForwardIterator>)
    646 	while (__f != __t)
    647 	  *__f++ = this->operator()(__urng, __param);
    648       }
    649 
    650   template<size_t _Dimen, typename _RealType>
    651     bool
    652     operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
    653 	       __d1,
    654 	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
    655 	       __d2)
    656     {
    657       return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd;
    658     }
    659 
    660   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
    661     std::basic_ostream<_CharT, _Traits>&
    662     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    663 	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
    664     {
    665       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    666       typedef typename __ostream_type::ios_base    __ios_base;
    667 
    668       const typename __ios_base::fmtflags __flags = __os.flags();
    669       const _CharT __fill = __os.fill();
    670       const std::streamsize __precision = __os.precision();
    671       const _CharT __space = __os.widen(' ');
    672       __os.flags(__ios_base::scientific | __ios_base::left);
    673       __os.fill(__space);
    674       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    675 
    676       auto __mean = __x._M_param.mean();
    677       for (auto __it : __mean)
    678 	__os << __it << __space;
    679       auto __t = __x._M_param.varcov();
    680       for (auto __it : __t)
    681 	__os << __it << __space;
    682 
    683       __os << __x._M_nd;
    684 
    685       __os.flags(__flags);
    686       __os.fill(__fill);
    687       __os.precision(__precision);
    688       return __os;
    689     }
    690 
    691   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
    692     std::basic_istream<_CharT, _Traits>&
    693     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    694 	       __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
    695     {
    696       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    697       typedef typename __istream_type::ios_base    __ios_base;
    698 
    699       const typename __ios_base::fmtflags __flags = __is.flags();
    700       __is.flags(__ios_base::dec | __ios_base::skipws);
    701 
    702       std::array<_RealType, _Dimen> __mean;
    703       for (auto& __it : __mean)
    704 	__is >> __it;
    705       std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov;
    706       for (auto& __it : __varcov)
    707 	__is >> __it;
    708 
    709       __is >> __x._M_nd;
    710 
    711       __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
    712 		param_type(__mean.begin(), __mean.end(),
    713 			   __varcov.begin(), __varcov.end()));
    714 
    715       __is.flags(__flags);
    716       return __is;
    717     }
    718 
    719 
    720   template<typename _RealType>
    721     template<typename _OutputIterator,
    722 	     typename _UniformRandomNumberGenerator>
    723       void
    724       rice_distribution<_RealType>::
    725       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    726 		      _UniformRandomNumberGenerator& __urng,
    727 		      const param_type& __p)
    728       {
    729 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
    730 
    731 	while (__f != __t)
    732 	  {
    733 	    typename std::normal_distribution<result_type>::param_type
    734 	      __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
    735 	    result_type __x = this->_M_ndx(__px, __urng);
    736 	    result_type __y = this->_M_ndy(__py, __urng);
    737 #if _GLIBCXX_USE_C99_MATH_TR1
    738 	    *__f++ = std::hypot(__x, __y);
    739 #else
    740 	    *__f++ = std::sqrt(__x * __x + __y * __y);
    741 #endif
    742 	  }
    743       }
    744 
    745   template<typename _RealType, typename _CharT, typename _Traits>
    746     std::basic_ostream<_CharT, _Traits>&
    747     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    748 	       const rice_distribution<_RealType>& __x)
    749     {
    750       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    751       typedef typename __ostream_type::ios_base    __ios_base;
    752 
    753       const typename __ios_base::fmtflags __flags = __os.flags();
    754       const _CharT __fill = __os.fill();
    755       const std::streamsize __precision = __os.precision();
    756       const _CharT __space = __os.widen(' ');
    757       __os.flags(__ios_base::scientific | __ios_base::left);
    758       __os.fill(__space);
    759       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    760 
    761       __os << __x.nu() << __space << __x.sigma();
    762       __os << __space << __x._M_ndx;
    763       __os << __space << __x._M_ndy;
    764 
    765       __os.flags(__flags);
    766       __os.fill(__fill);
    767       __os.precision(__precision);
    768       return __os;
    769     }
    770 
    771   template<typename _RealType, typename _CharT, typename _Traits>
    772     std::basic_istream<_CharT, _Traits>&
    773     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    774 	       rice_distribution<_RealType>& __x)
    775     {
    776       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    777       typedef typename __istream_type::ios_base    __ios_base;
    778 
    779       const typename __ios_base::fmtflags __flags = __is.flags();
    780       __is.flags(__ios_base::dec | __ios_base::skipws);
    781 
    782       _RealType __nu_val, __sigma_val;
    783       __is >> __nu_val >> __sigma_val;
    784       __is >> __x._M_ndx;
    785       __is >> __x._M_ndy;
    786       __x.param(typename rice_distribution<_RealType>::
    787 		param_type(__nu_val, __sigma_val));
    788 
    789       __is.flags(__flags);
    790       return __is;
    791     }
    792 
    793 
    794   template<typename _RealType>
    795     template<typename _OutputIterator,
    796 	     typename _UniformRandomNumberGenerator>
    797       void
    798       nakagami_distribution<_RealType>::
    799       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    800 		      _UniformRandomNumberGenerator& __urng,
    801 		      const param_type& __p)
    802       {
    803 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
    804 
    805 	typename std::gamma_distribution<result_type>::param_type
    806 	  __pg(__p.mu(), __p.omega() / __p.mu());
    807 	while (__f != __t)
    808 	  *__f++ = std::sqrt(this->_M_gd(__pg, __urng));
    809       }
    810 
    811   template<typename _RealType, typename _CharT, typename _Traits>
    812     std::basic_ostream<_CharT, _Traits>&
    813     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    814 	       const nakagami_distribution<_RealType>& __x)
    815     {
    816       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    817       typedef typename __ostream_type::ios_base    __ios_base;
    818 
    819       const typename __ios_base::fmtflags __flags = __os.flags();
    820       const _CharT __fill = __os.fill();
    821       const std::streamsize __precision = __os.precision();
    822       const _CharT __space = __os.widen(' ');
    823       __os.flags(__ios_base::scientific | __ios_base::left);
    824       __os.fill(__space);
    825       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    826 
    827       __os << __x.mu() << __space << __x.omega();
    828       __os << __space << __x._M_gd;
    829 
    830       __os.flags(__flags);
    831       __os.fill(__fill);
    832       __os.precision(__precision);
    833       return __os;
    834     }
    835 
    836   template<typename _RealType, typename _CharT, typename _Traits>
    837     std::basic_istream<_CharT, _Traits>&
    838     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    839 	       nakagami_distribution<_RealType>& __x)
    840     {
    841       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    842       typedef typename __istream_type::ios_base    __ios_base;
    843 
    844       const typename __ios_base::fmtflags __flags = __is.flags();
    845       __is.flags(__ios_base::dec | __ios_base::skipws);
    846 
    847       _RealType __mu_val, __omega_val;
    848       __is >> __mu_val >> __omega_val;
    849       __is >> __x._M_gd;
    850       __x.param(typename nakagami_distribution<_RealType>::
    851 		param_type(__mu_val, __omega_val));
    852 
    853       __is.flags(__flags);
    854       return __is;
    855     }
    856 
    857 
    858   template<typename _RealType>
    859     template<typename _OutputIterator,
    860 	     typename _UniformRandomNumberGenerator>
    861       void
    862       pareto_distribution<_RealType>::
    863       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    864 		      _UniformRandomNumberGenerator& __urng,
    865 		      const param_type& __p)
    866       {
    867 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
    868 
    869 	result_type __mu_val = __p.mu();
    870 	result_type __malphinv = -result_type(1) / __p.alpha();
    871 	while (__f != __t)
    872 	  *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv);
    873       }
    874 
    875   template<typename _RealType, typename _CharT, typename _Traits>
    876     std::basic_ostream<_CharT, _Traits>&
    877     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    878 	       const pareto_distribution<_RealType>& __x)
    879     {
    880       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    881       typedef typename __ostream_type::ios_base    __ios_base;
    882 
    883       const typename __ios_base::fmtflags __flags = __os.flags();
    884       const _CharT __fill = __os.fill();
    885       const std::streamsize __precision = __os.precision();
    886       const _CharT __space = __os.widen(' ');
    887       __os.flags(__ios_base::scientific | __ios_base::left);
    888       __os.fill(__space);
    889       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    890 
    891       __os << __x.alpha() << __space << __x.mu();
    892       __os << __space << __x._M_ud;
    893 
    894       __os.flags(__flags);
    895       __os.fill(__fill);
    896       __os.precision(__precision);
    897       return __os;
    898     }
    899 
    900   template<typename _RealType, typename _CharT, typename _Traits>
    901     std::basic_istream<_CharT, _Traits>&
    902     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    903 	       pareto_distribution<_RealType>& __x)
    904     {
    905       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    906       typedef typename __istream_type::ios_base    __ios_base;
    907 
    908       const typename __ios_base::fmtflags __flags = __is.flags();
    909       __is.flags(__ios_base::dec | __ios_base::skipws);
    910 
    911       _RealType __alpha_val, __mu_val;
    912       __is >> __alpha_val >> __mu_val;
    913       __is >> __x._M_ud;
    914       __x.param(typename pareto_distribution<_RealType>::
    915 		param_type(__alpha_val, __mu_val));
    916 
    917       __is.flags(__flags);
    918       return __is;
    919     }
    920 
    921 
    922   template<typename _RealType>
    923     template<typename _UniformRandomNumberGenerator>
    924       typename k_distribution<_RealType>::result_type
    925       k_distribution<_RealType>::
    926       operator()(_UniformRandomNumberGenerator& __urng)
    927       {
    928 	result_type __x = this->_M_gd1(__urng);
    929 	result_type __y = this->_M_gd2(__urng);
    930 	return std::sqrt(__x * __y);
    931       }
    932 
    933   template<typename _RealType>
    934     template<typename _UniformRandomNumberGenerator>
    935       typename k_distribution<_RealType>::result_type
    936       k_distribution<_RealType>::
    937       operator()(_UniformRandomNumberGenerator& __urng,
    938 		 const param_type& __p)
    939       {
    940 	typename std::gamma_distribution<result_type>::param_type
    941 	  __p1(__p.lambda(), result_type(1) / __p.lambda()),
    942 	  __p2(__p.nu(), __p.mu() / __p.nu());
    943 	result_type __x = this->_M_gd1(__p1, __urng);
    944 	result_type __y = this->_M_gd2(__p2, __urng);
    945 	return std::sqrt(__x * __y);
    946       }
    947 
    948   template<typename _RealType>
    949     template<typename _OutputIterator,
    950 	     typename _UniformRandomNumberGenerator>
    951       void
    952       k_distribution<_RealType>::
    953       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    954 		      _UniformRandomNumberGenerator& __urng,
    955 		      const param_type& __p)
    956       {
    957 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
    958 
    959 	typename std::gamma_distribution<result_type>::param_type
    960 	  __p1(__p.lambda(), result_type(1) / __p.lambda()),
    961 	  __p2(__p.nu(), __p.mu() / __p.nu());
    962 	while (__f != __t)
    963 	  {
    964 	    result_type __x = this->_M_gd1(__p1, __urng);
    965 	    result_type __y = this->_M_gd2(__p2, __urng);
    966 	    *__f++ = std::sqrt(__x * __y);
    967 	  }
    968       }
    969 
    970   template<typename _RealType, typename _CharT, typename _Traits>
    971     std::basic_ostream<_CharT, _Traits>&
    972     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    973 	       const k_distribution<_RealType>& __x)
    974     {
    975       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    976       typedef typename __ostream_type::ios_base    __ios_base;
    977 
    978       const typename __ios_base::fmtflags __flags = __os.flags();
    979       const _CharT __fill = __os.fill();
    980       const std::streamsize __precision = __os.precision();
    981       const _CharT __space = __os.widen(' ');
    982       __os.flags(__ios_base::scientific | __ios_base::left);
    983       __os.fill(__space);
    984       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    985 
    986       __os << __x.lambda() << __space << __x.mu() << __space << __x.nu();
    987       __os << __space << __x._M_gd1;
    988       __os << __space << __x._M_gd2;
    989 
    990       __os.flags(__flags);
    991       __os.fill(__fill);
    992       __os.precision(__precision);
    993       return __os;
    994     }
    995 
    996   template<typename _RealType, typename _CharT, typename _Traits>
    997     std::basic_istream<_CharT, _Traits>&
    998     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    999 	       k_distribution<_RealType>& __x)
   1000     {
   1001       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1002       typedef typename __istream_type::ios_base    __ios_base;
   1003 
   1004       const typename __ios_base::fmtflags __flags = __is.flags();
   1005       __is.flags(__ios_base::dec | __ios_base::skipws);
   1006 
   1007       _RealType __lambda_val, __mu_val, __nu_val;
   1008       __is >> __lambda_val >> __mu_val >> __nu_val;
   1009       __is >> __x._M_gd1;
   1010       __is >> __x._M_gd2;
   1011       __x.param(typename k_distribution<_RealType>::
   1012 		param_type(__lambda_val, __mu_val, __nu_val));
   1013 
   1014       __is.flags(__flags);
   1015       return __is;
   1016     }
   1017 
   1018 
   1019   template<typename _RealType>
   1020     template<typename _OutputIterator,
   1021 	     typename _UniformRandomNumberGenerator>
   1022       void
   1023       arcsine_distribution<_RealType>::
   1024       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1025 		      _UniformRandomNumberGenerator& __urng,
   1026 		      const param_type& __p)
   1027       {
   1028 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
   1029 
   1030 	result_type __dif = __p.b() - __p.a();
   1031 	result_type __sum = __p.a() + __p.b();
   1032 	while (__f != __t)
   1033 	  {
   1034 	    result_type __x = std::sin(this->_M_ud(__urng));
   1035 	    *__f++ = (__x * __dif + __sum) / result_type(2);
   1036 	  }
   1037       }
   1038 
   1039   template<typename _RealType, typename _CharT, typename _Traits>
   1040     std::basic_ostream<_CharT, _Traits>&
   1041     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1042 	       const arcsine_distribution<_RealType>& __x)
   1043     {
   1044       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1045       typedef typename __ostream_type::ios_base    __ios_base;
   1046 
   1047       const typename __ios_base::fmtflags __flags = __os.flags();
   1048       const _CharT __fill = __os.fill();
   1049       const std::streamsize __precision = __os.precision();
   1050       const _CharT __space = __os.widen(' ');
   1051       __os.flags(__ios_base::scientific | __ios_base::left);
   1052       __os.fill(__space);
   1053       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1054 
   1055       __os << __x.a() << __space << __x.b();
   1056       __os << __space << __x._M_ud;
   1057 
   1058       __os.flags(__flags);
   1059       __os.fill(__fill);
   1060       __os.precision(__precision);
   1061       return __os;
   1062     }
   1063 
   1064   template<typename _RealType, typename _CharT, typename _Traits>
   1065     std::basic_istream<_CharT, _Traits>&
   1066     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1067 	       arcsine_distribution<_RealType>& __x)
   1068     {
   1069       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1070       typedef typename __istream_type::ios_base    __ios_base;
   1071 
   1072       const typename __ios_base::fmtflags __flags = __is.flags();
   1073       __is.flags(__ios_base::dec | __ios_base::skipws);
   1074 
   1075       _RealType __a, __b;
   1076       __is >> __a >> __b;
   1077       __is >> __x._M_ud;
   1078       __x.param(typename arcsine_distribution<_RealType>::
   1079 		param_type(__a, __b));
   1080 
   1081       __is.flags(__flags);
   1082       return __is;
   1083     }
   1084 
   1085 
   1086   template<typename _RealType>
   1087     template<typename _UniformRandomNumberGenerator>
   1088       typename hoyt_distribution<_RealType>::result_type
   1089       hoyt_distribution<_RealType>::
   1090       operator()(_UniformRandomNumberGenerator& __urng)
   1091       {
   1092 	result_type __x = this->_M_ad(__urng);
   1093 	result_type __y = this->_M_ed(__urng);
   1094 	return (result_type(2) * this->q()
   1095 		  / (result_type(1) + this->q() * this->q()))
   1096 	       * std::sqrt(this->omega() * __x * __y);
   1097       }
   1098 
   1099   template<typename _RealType>
   1100     template<typename _UniformRandomNumberGenerator>
   1101       typename hoyt_distribution<_RealType>::result_type
   1102       hoyt_distribution<_RealType>::
   1103       operator()(_UniformRandomNumberGenerator& __urng,
   1104 		 const param_type& __p)
   1105       {
   1106 	result_type __q2 = __p.q() * __p.q();
   1107 	result_type __num = result_type(0.5L) * (result_type(1) + __q2);
   1108 	typename __gnu_cxx::arcsine_distribution<result_type>::param_type
   1109 	  __pa(__num, __num / __q2);
   1110 	result_type __x = this->_M_ad(__pa, __urng);
   1111 	result_type __y = this->_M_ed(__urng);
   1112 	return (result_type(2) * __p.q() / (result_type(1) + __q2))
   1113 	       * std::sqrt(__p.omega() * __x * __y);
   1114       }
   1115 
   1116   template<typename _RealType>
   1117     template<typename _OutputIterator,
   1118 	     typename _UniformRandomNumberGenerator>
   1119       void
   1120       hoyt_distribution<_RealType>::
   1121       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1122 		      _UniformRandomNumberGenerator& __urng,
   1123 		      const param_type& __p)
   1124       {
   1125 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
   1126 
   1127 	result_type __2q = result_type(2) * __p.q();
   1128 	result_type __q2 = __p.q() * __p.q();
   1129 	result_type __q2p1 = result_type(1) + __q2;
   1130 	result_type __num = result_type(0.5L) * __q2p1;
   1131 	result_type __omega = __p.omega();
   1132 	typename __gnu_cxx::arcsine_distribution<result_type>::param_type
   1133 	  __pa(__num, __num / __q2);
   1134 	while (__f != __t)
   1135 	  {
   1136 	    result_type __x = this->_M_ad(__pa, __urng);
   1137 	    result_type __y = this->_M_ed(__urng);
   1138 	    *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y);
   1139 	  }
   1140       }
   1141 
   1142   template<typename _RealType, typename _CharT, typename _Traits>
   1143     std::basic_ostream<_CharT, _Traits>&
   1144     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1145 	       const hoyt_distribution<_RealType>& __x)
   1146     {
   1147       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1148       typedef typename __ostream_type::ios_base    __ios_base;
   1149 
   1150       const typename __ios_base::fmtflags __flags = __os.flags();
   1151       const _CharT __fill = __os.fill();
   1152       const std::streamsize __precision = __os.precision();
   1153       const _CharT __space = __os.widen(' ');
   1154       __os.flags(__ios_base::scientific | __ios_base::left);
   1155       __os.fill(__space);
   1156       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1157 
   1158       __os << __x.q() << __space << __x.omega();
   1159       __os << __space << __x._M_ad;
   1160       __os << __space << __x._M_ed;
   1161 
   1162       __os.flags(__flags);
   1163       __os.fill(__fill);
   1164       __os.precision(__precision);
   1165       return __os;
   1166     }
   1167 
   1168   template<typename _RealType, typename _CharT, typename _Traits>
   1169     std::basic_istream<_CharT, _Traits>&
   1170     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1171 	       hoyt_distribution<_RealType>& __x)
   1172     {
   1173       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1174       typedef typename __istream_type::ios_base    __ios_base;
   1175 
   1176       const typename __ios_base::fmtflags __flags = __is.flags();
   1177       __is.flags(__ios_base::dec | __ios_base::skipws);
   1178 
   1179       _RealType __q, __omega;
   1180       __is >> __q >> __omega;
   1181       __is >> __x._M_ad;
   1182       __is >> __x._M_ed;
   1183       __x.param(typename hoyt_distribution<_RealType>::
   1184 		param_type(__q, __omega));
   1185 
   1186       __is.flags(__flags);
   1187       return __is;
   1188     }
   1189 
   1190 
   1191   template<typename _RealType>
   1192     template<typename _OutputIterator,
   1193 	     typename _UniformRandomNumberGenerator>
   1194       void
   1195       triangular_distribution<_RealType>::
   1196       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1197 		      _UniformRandomNumberGenerator& __urng,
   1198 		      const param_type& __param)
   1199       {
   1200 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
   1201 
   1202 	while (__f != __t)
   1203 	  *__f++ = this->operator()(__urng, __param);
   1204       }
   1205 
   1206   template<typename _RealType, typename _CharT, typename _Traits>
   1207     std::basic_ostream<_CharT, _Traits>&
   1208     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1209 	       const __gnu_cxx::triangular_distribution<_RealType>& __x)
   1210     {
   1211       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1212       typedef typename __ostream_type::ios_base    __ios_base;
   1213 
   1214       const typename __ios_base::fmtflags __flags = __os.flags();
   1215       const _CharT __fill = __os.fill();
   1216       const std::streamsize __precision = __os.precision();
   1217       const _CharT __space = __os.widen(' ');
   1218       __os.flags(__ios_base::scientific | __ios_base::left);
   1219       __os.fill(__space);
   1220       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1221 
   1222       __os << __x.a() << __space << __x.b() << __space << __x.c();
   1223 
   1224       __os.flags(__flags);
   1225       __os.fill(__fill);
   1226       __os.precision(__precision);
   1227       return __os;
   1228     }
   1229 
   1230   template<typename _RealType, typename _CharT, typename _Traits>
   1231     std::basic_istream<_CharT, _Traits>&
   1232     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1233 	       __gnu_cxx::triangular_distribution<_RealType>& __x)
   1234     {
   1235       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1236       typedef typename __istream_type::ios_base    __ios_base;
   1237 
   1238       const typename __ios_base::fmtflags __flags = __is.flags();
   1239       __is.flags(__ios_base::dec | __ios_base::skipws);
   1240 
   1241       _RealType __a, __b, __c;
   1242       __is >> __a >> __b >> __c;
   1243       __x.param(typename __gnu_cxx::triangular_distribution<_RealType>::
   1244 		param_type(__a, __b, __c));
   1245 
   1246       __is.flags(__flags);
   1247       return __is;
   1248     }
   1249 
   1250 
   1251   template<typename _RealType>
   1252     template<typename _OutputIterator,
   1253 	     typename _UniformRandomNumberGenerator>
   1254       void
   1255       von_mises_distribution<_RealType>::
   1256       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1257 		      _UniformRandomNumberGenerator& __urng,
   1258 		      const param_type& __param)
   1259       {
   1260 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
   1261 
   1262 	while (__f != __t)
   1263 	  *__f++ = this->operator()(__urng, __param);
   1264       }
   1265 
   1266   template<typename _RealType, typename _CharT, typename _Traits>
   1267     std::basic_ostream<_CharT, _Traits>&
   1268     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1269 	       const __gnu_cxx::von_mises_distribution<_RealType>& __x)
   1270     {
   1271       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1272       typedef typename __ostream_type::ios_base    __ios_base;
   1273 
   1274       const typename __ios_base::fmtflags __flags = __os.flags();
   1275       const _CharT __fill = __os.fill();
   1276       const std::streamsize __precision = __os.precision();
   1277       const _CharT __space = __os.widen(' ');
   1278       __os.flags(__ios_base::scientific | __ios_base::left);
   1279       __os.fill(__space);
   1280       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1281 
   1282       __os << __x.mu() << __space << __x.kappa();
   1283 
   1284       __os.flags(__flags);
   1285       __os.fill(__fill);
   1286       __os.precision(__precision);
   1287       return __os;
   1288     }
   1289 
   1290   template<typename _RealType, typename _CharT, typename _Traits>
   1291     std::basic_istream<_CharT, _Traits>&
   1292     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1293 	       __gnu_cxx::von_mises_distribution<_RealType>& __x)
   1294     {
   1295       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1296       typedef typename __istream_type::ios_base    __ios_base;
   1297 
   1298       const typename __ios_base::fmtflags __flags = __is.flags();
   1299       __is.flags(__ios_base::dec | __ios_base::skipws);
   1300 
   1301       _RealType __mu, __kappa;
   1302       __is >> __mu >> __kappa;
   1303       __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>::
   1304 		param_type(__mu, __kappa));
   1305 
   1306       __is.flags(__flags);
   1307       return __is;
   1308     }
   1309 
   1310 _GLIBCXX_END_NAMESPACE_VERSION
   1311 } // namespace
   1312 
   1313 
   1314 #endif // _EXT_RANDOM_TCC
   1315