Home | History | Annotate | Line # | Download | only in ext
      1 // Random number extensions -*- C++ -*-
      2 
      3 // Copyright (C) 2012-2022 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 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
     36 {
     37 _GLIBCXX_BEGIN_NAMESPACE_VERSION
     38 
     39 #if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
     40 
     41   template<typename _UIntType, size_t __m,
     42 	   size_t __pos1, size_t __sl1, size_t __sl2,
     43 	   size_t __sr1, size_t __sr2,
     44 	   uint32_t __msk1, uint32_t __msk2,
     45 	   uint32_t __msk3, uint32_t __msk4,
     46 	   uint32_t __parity1, uint32_t __parity2,
     47 	   uint32_t __parity3, uint32_t __parity4>
     48     void simd_fast_mersenne_twister_engine<_UIntType, __m,
     49 					   __pos1, __sl1, __sl2, __sr1, __sr2,
     50 					   __msk1, __msk2, __msk3, __msk4,
     51 					   __parity1, __parity2, __parity3,
     52 					   __parity4>::
     53     seed(_UIntType __seed)
     54     {
     55       _M_state32[0] = static_cast<uint32_t>(__seed);
     56       for (size_t __i = 1; __i < _M_nstate32; ++__i)
     57 	_M_state32[__i] = (1812433253UL
     58 			   * (_M_state32[__i - 1] ^ (_M_state32[__i - 1] >> 30))
     59 			   + __i);
     60       _M_pos = state_size;
     61       _M_period_certification();
     62     }
     63 
     64 
     65   namespace {
     66 
     67     inline uint32_t _Func1(uint32_t __x)
     68     {
     69       return (__x ^ (__x >> 27)) * UINT32_C(1664525);
     70     }
     71 
     72     inline uint32_t _Func2(uint32_t __x)
     73     {
     74       return (__x ^ (__x >> 27)) * UINT32_C(1566083941);
     75     }
     76 
     77   }
     78 
     79 
     80   template<typename _UIntType, size_t __m,
     81 	   size_t __pos1, size_t __sl1, size_t __sl2,
     82 	   size_t __sr1, size_t __sr2,
     83 	   uint32_t __msk1, uint32_t __msk2,
     84 	   uint32_t __msk3, uint32_t __msk4,
     85 	   uint32_t __parity1, uint32_t __parity2,
     86 	   uint32_t __parity3, uint32_t __parity4>
     87     template<typename _Sseq>
     88       auto
     89       simd_fast_mersenne_twister_engine<_UIntType, __m,
     90 					__pos1, __sl1, __sl2, __sr1, __sr2,
     91 					__msk1, __msk2, __msk3, __msk4,
     92 					__parity1, __parity2, __parity3,
     93 					__parity4>::
     94       seed(_Sseq& __q)
     95       -> _If_seed_seq<_Sseq>
     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 	    result_type>)
    446 
    447 	std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
    448 	  __aurng(__urng);
    449 
    450 	while (__f != __t)
    451 	  {
    452 	    result_type __x, __y;
    453 	    do
    454 	      {
    455 		__x = std::exp(std::log(__aurng()) / __param.alpha());
    456 		__y = std::exp(std::log(__aurng()) / __param.beta());
    457 	      }
    458 	    while (__x + __y > result_type(1));
    459 
    460 	    *__f++ = __x / (__x + __y);
    461 	  }
    462       }
    463 
    464   template<typename _RealType, typename _CharT, typename _Traits>
    465     std::basic_ostream<_CharT, _Traits>&
    466     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    467 	       const __gnu_cxx::beta_distribution<_RealType>& __x)
    468     {
    469       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    470       typedef typename __ostream_type::ios_base    __ios_base;
    471 
    472       const typename __ios_base::fmtflags __flags = __os.flags();
    473       const _CharT __fill = __os.fill();
    474       const std::streamsize __precision = __os.precision();
    475       const _CharT __space = __os.widen(' ');
    476       __os.flags(__ios_base::scientific | __ios_base::left);
    477       __os.fill(__space);
    478       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    479 
    480       __os << __x.alpha() << __space << __x.beta();
    481 
    482       __os.flags(__flags);
    483       __os.fill(__fill);
    484       __os.precision(__precision);
    485       return __os;
    486     }
    487 
    488   template<typename _RealType, typename _CharT, typename _Traits>
    489     std::basic_istream<_CharT, _Traits>&
    490     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    491 	       __gnu_cxx::beta_distribution<_RealType>& __x)
    492     {
    493       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    494       typedef typename __istream_type::ios_base    __ios_base;
    495 
    496       const typename __ios_base::fmtflags __flags = __is.flags();
    497       __is.flags(__ios_base::dec | __ios_base::skipws);
    498 
    499       _RealType __alpha_val, __beta_val;
    500       __is >> __alpha_val >> __beta_val;
    501       __x.param(typename __gnu_cxx::beta_distribution<_RealType>::
    502 		param_type(__alpha_val, __beta_val));
    503 
    504       __is.flags(__flags);
    505       return __is;
    506     }
    507 
    508 
    509   template<std::size_t _Dimen, typename _RealType>
    510     template<typename _InputIterator1, typename _InputIterator2>
    511       void
    512       normal_mv_distribution<_Dimen, _RealType>::param_type::
    513       _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
    514 		   _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
    515       {
    516 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
    517 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
    518 	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
    519 		  _M_mean.end(), _RealType(0));
    520 
    521 	// Perform the Cholesky decomposition
    522 	auto __w = _M_t.begin();
    523 	for (size_t __j = 0; __j < _Dimen; ++__j)
    524 	  {
    525 	    _RealType __sum = _RealType(0);
    526 
    527 	    auto __slitbegin = __w;
    528 	    auto __cit = _M_t.begin();
    529 	    for (size_t __i = 0; __i < __j; ++__i)
    530 	      {
    531 		auto __slit = __slitbegin;
    532 		_RealType __s = *__varcovbegin++;
    533 		for (size_t __k = 0; __k < __i; ++__k)
    534 		  __s -= *__slit++ * *__cit++;
    535 
    536 		*__w++ = __s /= *__cit++;
    537 		__sum += __s * __s;
    538 	      }
    539 
    540 	    __sum = *__varcovbegin - __sum;
    541 	    if (__builtin_expect(__sum <= _RealType(0), 0))
    542 	      std::__throw_runtime_error(__N("normal_mv_distribution::"
    543 					     "param_type::_M_init_full"));
    544 	    *__w++ = std::sqrt(__sum);
    545 
    546 	    std::advance(__varcovbegin, _Dimen - __j);
    547 	  }
    548       }
    549 
    550   template<std::size_t _Dimen, typename _RealType>
    551     template<typename _InputIterator1, typename _InputIterator2>
    552       void
    553       normal_mv_distribution<_Dimen, _RealType>::param_type::
    554       _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
    555 		    _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
    556       {
    557 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
    558 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
    559 	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
    560 		  _M_mean.end(), _RealType(0));
    561 
    562 	// Perform the Cholesky decomposition
    563 	auto __w = _M_t.begin();
    564 	for (size_t __j = 0; __j < _Dimen; ++__j)
    565 	  {
    566 	    _RealType __sum = _RealType(0);
    567 
    568 	    auto __slitbegin = __w;
    569 	    auto __cit = _M_t.begin();
    570 	    for (size_t __i = 0; __i < __j; ++__i)
    571 	      {
    572 		auto __slit = __slitbegin;
    573 		_RealType __s = *__varcovbegin++;
    574 		for (size_t __k = 0; __k < __i; ++__k)
    575 		  __s -= *__slit++ * *__cit++;
    576 
    577 		*__w++ = __s /= *__cit++;
    578 		__sum += __s * __s;
    579 	      }
    580 
    581 	    __sum = *__varcovbegin++ - __sum;
    582 	    if (__builtin_expect(__sum <= _RealType(0), 0))
    583 	      std::__throw_runtime_error(__N("normal_mv_distribution::"
    584 					     "param_type::_M_init_lower"));
    585 	    *__w++ = std::sqrt(__sum);
    586 	  }
    587       }
    588 
    589   template<std::size_t _Dimen, typename _RealType>
    590     template<typename _InputIterator1, typename _InputIterator2>
    591       void
    592       normal_mv_distribution<_Dimen, _RealType>::param_type::
    593       _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
    594 		       _InputIterator2 __varbegin, _InputIterator2 __varend)
    595       {
    596 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
    597 	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
    598 	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
    599 		  _M_mean.end(), _RealType(0));
    600 
    601 	auto __w = _M_t.begin();
    602 	size_t __step = 0;
    603 	while (__varbegin != __varend)
    604 	  {
    605 	    std::fill_n(__w, __step, _RealType(0));
    606 	    __w += __step++;
    607 	    if (__builtin_expect(*__varbegin < _RealType(0), 0))
    608 	      std::__throw_runtime_error(__N("normal_mv_distribution::"
    609 					     "param_type::_M_init_diagonal"));
    610 	    *__w++ = std::sqrt(*__varbegin++);
    611 	  }
    612       }
    613 
    614   template<std::size_t _Dimen, typename _RealType>
    615     template<typename _UniformRandomNumberGenerator>
    616       typename normal_mv_distribution<_Dimen, _RealType>::result_type
    617       normal_mv_distribution<_Dimen, _RealType>::
    618       operator()(_UniformRandomNumberGenerator& __urng,
    619 		 const param_type& __param)
    620       {
    621 	result_type __ret;
    622 
    623 	_M_nd.__generate(__ret.begin(), __ret.end(), __urng);
    624 
    625 	auto __t_it = __param._M_t.crbegin();
    626 	for (size_t __i = _Dimen; __i > 0; --__i)
    627 	  {
    628 	    _RealType __sum = _RealType(0);
    629 	    for (size_t __j = __i; __j > 0; --__j)
    630 	      __sum += __ret[__j - 1] * *__t_it++;
    631 	    __ret[__i - 1] = __sum;
    632 	  }
    633 
    634 	return __ret;
    635       }
    636 
    637   template<std::size_t _Dimen, typename _RealType>
    638     template<typename _ForwardIterator, typename _UniformRandomNumberGenerator>
    639       void
    640       normal_mv_distribution<_Dimen, _RealType>::
    641       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
    642 		      _UniformRandomNumberGenerator& __urng,
    643 		      const param_type& __param)
    644       {
    645 	__glibcxx_function_requires(_Mutable_ForwardIteratorConcept<
    646 				    _ForwardIterator>)
    647 	while (__f != __t)
    648 	  *__f++ = this->operator()(__urng, __param);
    649       }
    650 
    651   template<size_t _Dimen, typename _RealType>
    652     bool
    653     operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
    654 	       __d1,
    655 	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
    656 	       __d2)
    657     {
    658       return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd;
    659     }
    660 
    661   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
    662     std::basic_ostream<_CharT, _Traits>&
    663     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    664 	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
    665     {
    666       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    667       typedef typename __ostream_type::ios_base    __ios_base;
    668 
    669       const typename __ios_base::fmtflags __flags = __os.flags();
    670       const _CharT __fill = __os.fill();
    671       const std::streamsize __precision = __os.precision();
    672       const _CharT __space = __os.widen(' ');
    673       __os.flags(__ios_base::scientific | __ios_base::left);
    674       __os.fill(__space);
    675       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    676 
    677       auto __mean = __x._M_param.mean();
    678       for (auto __it : __mean)
    679 	__os << __it << __space;
    680       auto __t = __x._M_param.varcov();
    681       for (auto __it : __t)
    682 	__os << __it << __space;
    683 
    684       __os << __x._M_nd;
    685 
    686       __os.flags(__flags);
    687       __os.fill(__fill);
    688       __os.precision(__precision);
    689       return __os;
    690     }
    691 
    692   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
    693     std::basic_istream<_CharT, _Traits>&
    694     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    695 	       __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
    696     {
    697       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    698       typedef typename __istream_type::ios_base    __ios_base;
    699 
    700       const typename __ios_base::fmtflags __flags = __is.flags();
    701       __is.flags(__ios_base::dec | __ios_base::skipws);
    702 
    703       std::array<_RealType, _Dimen> __mean;
    704       for (auto& __it : __mean)
    705 	__is >> __it;
    706       std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov;
    707       for (auto& __it : __varcov)
    708 	__is >> __it;
    709 
    710       __is >> __x._M_nd;
    711 
    712       // The param_type temporary is built with a private constructor,
    713       // to skip the Cholesky decomposition that would be performed
    714       // otherwise.
    715       __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
    716 		param_type(__mean, __varcov));
    717 
    718       __is.flags(__flags);
    719       return __is;
    720     }
    721 
    722 
    723   template<typename _RealType>
    724     template<typename _OutputIterator,
    725 	     typename _UniformRandomNumberGenerator>
    726       void
    727       rice_distribution<_RealType>::
    728       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    729 		      _UniformRandomNumberGenerator& __urng,
    730 		      const param_type& __p)
    731       {
    732 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
    733 	    result_type>)
    734 
    735 	while (__f != __t)
    736 	  {
    737 	    typename std::normal_distribution<result_type>::param_type
    738 	      __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
    739 	    result_type __x = this->_M_ndx(__px, __urng);
    740 	    result_type __y = this->_M_ndy(__py, __urng);
    741 #if _GLIBCXX_USE_C99_MATH_TR1
    742 	    *__f++ = std::hypot(__x, __y);
    743 #else
    744 	    *__f++ = std::sqrt(__x * __x + __y * __y);
    745 #endif
    746 	  }
    747       }
    748 
    749   template<typename _RealType, typename _CharT, typename _Traits>
    750     std::basic_ostream<_CharT, _Traits>&
    751     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    752 	       const rice_distribution<_RealType>& __x)
    753     {
    754       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    755       typedef typename __ostream_type::ios_base    __ios_base;
    756 
    757       const typename __ios_base::fmtflags __flags = __os.flags();
    758       const _CharT __fill = __os.fill();
    759       const std::streamsize __precision = __os.precision();
    760       const _CharT __space = __os.widen(' ');
    761       __os.flags(__ios_base::scientific | __ios_base::left);
    762       __os.fill(__space);
    763       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    764 
    765       __os << __x.nu() << __space << __x.sigma();
    766       __os << __space << __x._M_ndx;
    767       __os << __space << __x._M_ndy;
    768 
    769       __os.flags(__flags);
    770       __os.fill(__fill);
    771       __os.precision(__precision);
    772       return __os;
    773     }
    774 
    775   template<typename _RealType, typename _CharT, typename _Traits>
    776     std::basic_istream<_CharT, _Traits>&
    777     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    778 	       rice_distribution<_RealType>& __x)
    779     {
    780       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    781       typedef typename __istream_type::ios_base    __ios_base;
    782 
    783       const typename __ios_base::fmtflags __flags = __is.flags();
    784       __is.flags(__ios_base::dec | __ios_base::skipws);
    785 
    786       _RealType __nu_val, __sigma_val;
    787       __is >> __nu_val >> __sigma_val;
    788       __is >> __x._M_ndx;
    789       __is >> __x._M_ndy;
    790       __x.param(typename rice_distribution<_RealType>::
    791 		param_type(__nu_val, __sigma_val));
    792 
    793       __is.flags(__flags);
    794       return __is;
    795     }
    796 
    797 
    798   template<typename _RealType>
    799     template<typename _OutputIterator,
    800 	     typename _UniformRandomNumberGenerator>
    801       void
    802       nakagami_distribution<_RealType>::
    803       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    804 		      _UniformRandomNumberGenerator& __urng,
    805 		      const param_type& __p)
    806       {
    807 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
    808 	    result_type>)
    809 
    810 	typename std::gamma_distribution<result_type>::param_type
    811 	  __pg(__p.mu(), __p.omega() / __p.mu());
    812 	while (__f != __t)
    813 	  *__f++ = std::sqrt(this->_M_gd(__pg, __urng));
    814       }
    815 
    816   template<typename _RealType, typename _CharT, typename _Traits>
    817     std::basic_ostream<_CharT, _Traits>&
    818     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    819 	       const nakagami_distribution<_RealType>& __x)
    820     {
    821       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    822       typedef typename __ostream_type::ios_base    __ios_base;
    823 
    824       const typename __ios_base::fmtflags __flags = __os.flags();
    825       const _CharT __fill = __os.fill();
    826       const std::streamsize __precision = __os.precision();
    827       const _CharT __space = __os.widen(' ');
    828       __os.flags(__ios_base::scientific | __ios_base::left);
    829       __os.fill(__space);
    830       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    831 
    832       __os << __x.mu() << __space << __x.omega();
    833       __os << __space << __x._M_gd;
    834 
    835       __os.flags(__flags);
    836       __os.fill(__fill);
    837       __os.precision(__precision);
    838       return __os;
    839     }
    840 
    841   template<typename _RealType, typename _CharT, typename _Traits>
    842     std::basic_istream<_CharT, _Traits>&
    843     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    844 	       nakagami_distribution<_RealType>& __x)
    845     {
    846       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    847       typedef typename __istream_type::ios_base    __ios_base;
    848 
    849       const typename __ios_base::fmtflags __flags = __is.flags();
    850       __is.flags(__ios_base::dec | __ios_base::skipws);
    851 
    852       _RealType __mu_val, __omega_val;
    853       __is >> __mu_val >> __omega_val;
    854       __is >> __x._M_gd;
    855       __x.param(typename nakagami_distribution<_RealType>::
    856 		param_type(__mu_val, __omega_val));
    857 
    858       __is.flags(__flags);
    859       return __is;
    860     }
    861 
    862 
    863   template<typename _RealType>
    864     template<typename _OutputIterator,
    865 	     typename _UniformRandomNumberGenerator>
    866       void
    867       pareto_distribution<_RealType>::
    868       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    869 		      _UniformRandomNumberGenerator& __urng,
    870 		      const param_type& __p)
    871       {
    872 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
    873 	    result_type>)
    874 
    875 	result_type __mu_val = __p.mu();
    876 	result_type __malphinv = -result_type(1) / __p.alpha();
    877 	while (__f != __t)
    878 	  *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv);
    879       }
    880 
    881   template<typename _RealType, typename _CharT, typename _Traits>
    882     std::basic_ostream<_CharT, _Traits>&
    883     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    884 	       const pareto_distribution<_RealType>& __x)
    885     {
    886       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    887       typedef typename __ostream_type::ios_base    __ios_base;
    888 
    889       const typename __ios_base::fmtflags __flags = __os.flags();
    890       const _CharT __fill = __os.fill();
    891       const std::streamsize __precision = __os.precision();
    892       const _CharT __space = __os.widen(' ');
    893       __os.flags(__ios_base::scientific | __ios_base::left);
    894       __os.fill(__space);
    895       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    896 
    897       __os << __x.alpha() << __space << __x.mu();
    898       __os << __space << __x._M_ud;
    899 
    900       __os.flags(__flags);
    901       __os.fill(__fill);
    902       __os.precision(__precision);
    903       return __os;
    904     }
    905 
    906   template<typename _RealType, typename _CharT, typename _Traits>
    907     std::basic_istream<_CharT, _Traits>&
    908     operator>>(std::basic_istream<_CharT, _Traits>& __is,
    909 	       pareto_distribution<_RealType>& __x)
    910     {
    911       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
    912       typedef typename __istream_type::ios_base    __ios_base;
    913 
    914       const typename __ios_base::fmtflags __flags = __is.flags();
    915       __is.flags(__ios_base::dec | __ios_base::skipws);
    916 
    917       _RealType __alpha_val, __mu_val;
    918       __is >> __alpha_val >> __mu_val;
    919       __is >> __x._M_ud;
    920       __x.param(typename pareto_distribution<_RealType>::
    921 		param_type(__alpha_val, __mu_val));
    922 
    923       __is.flags(__flags);
    924       return __is;
    925     }
    926 
    927 
    928   template<typename _RealType>
    929     template<typename _UniformRandomNumberGenerator>
    930       typename k_distribution<_RealType>::result_type
    931       k_distribution<_RealType>::
    932       operator()(_UniformRandomNumberGenerator& __urng)
    933       {
    934 	result_type __x = this->_M_gd1(__urng);
    935 	result_type __y = this->_M_gd2(__urng);
    936 	return std::sqrt(__x * __y);
    937       }
    938 
    939   template<typename _RealType>
    940     template<typename _UniformRandomNumberGenerator>
    941       typename k_distribution<_RealType>::result_type
    942       k_distribution<_RealType>::
    943       operator()(_UniformRandomNumberGenerator& __urng,
    944 		 const param_type& __p)
    945       {
    946 	typename std::gamma_distribution<result_type>::param_type
    947 	  __p1(__p.lambda(), result_type(1) / __p.lambda()),
    948 	  __p2(__p.nu(), __p.mu() / __p.nu());
    949 	result_type __x = this->_M_gd1(__p1, __urng);
    950 	result_type __y = this->_M_gd2(__p2, __urng);
    951 	return std::sqrt(__x * __y);
    952       }
    953 
    954   template<typename _RealType>
    955     template<typename _OutputIterator,
    956 	     typename _UniformRandomNumberGenerator>
    957       void
    958       k_distribution<_RealType>::
    959       __generate_impl(_OutputIterator __f, _OutputIterator __t,
    960 		      _UniformRandomNumberGenerator& __urng,
    961 		      const param_type& __p)
    962       {
    963 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
    964 	    result_type>)
    965 
    966 	typename std::gamma_distribution<result_type>::param_type
    967 	  __p1(__p.lambda(), result_type(1) / __p.lambda()),
    968 	  __p2(__p.nu(), __p.mu() / __p.nu());
    969 	while (__f != __t)
    970 	  {
    971 	    result_type __x = this->_M_gd1(__p1, __urng);
    972 	    result_type __y = this->_M_gd2(__p2, __urng);
    973 	    *__f++ = std::sqrt(__x * __y);
    974 	  }
    975       }
    976 
    977   template<typename _RealType, typename _CharT, typename _Traits>
    978     std::basic_ostream<_CharT, _Traits>&
    979     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    980 	       const k_distribution<_RealType>& __x)
    981     {
    982       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
    983       typedef typename __ostream_type::ios_base    __ios_base;
    984 
    985       const typename __ios_base::fmtflags __flags = __os.flags();
    986       const _CharT __fill = __os.fill();
    987       const std::streamsize __precision = __os.precision();
    988       const _CharT __space = __os.widen(' ');
    989       __os.flags(__ios_base::scientific | __ios_base::left);
    990       __os.fill(__space);
    991       __os.precision(std::numeric_limits<_RealType>::max_digits10);
    992 
    993       __os << __x.lambda() << __space << __x.mu() << __space << __x.nu();
    994       __os << __space << __x._M_gd1;
    995       __os << __space << __x._M_gd2;
    996 
    997       __os.flags(__flags);
    998       __os.fill(__fill);
    999       __os.precision(__precision);
   1000       return __os;
   1001     }
   1002 
   1003   template<typename _RealType, typename _CharT, typename _Traits>
   1004     std::basic_istream<_CharT, _Traits>&
   1005     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1006 	       k_distribution<_RealType>& __x)
   1007     {
   1008       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1009       typedef typename __istream_type::ios_base    __ios_base;
   1010 
   1011       const typename __ios_base::fmtflags __flags = __is.flags();
   1012       __is.flags(__ios_base::dec | __ios_base::skipws);
   1013 
   1014       _RealType __lambda_val, __mu_val, __nu_val;
   1015       __is >> __lambda_val >> __mu_val >> __nu_val;
   1016       __is >> __x._M_gd1;
   1017       __is >> __x._M_gd2;
   1018       __x.param(typename k_distribution<_RealType>::
   1019 		param_type(__lambda_val, __mu_val, __nu_val));
   1020 
   1021       __is.flags(__flags);
   1022       return __is;
   1023     }
   1024 
   1025 
   1026   template<typename _RealType>
   1027     template<typename _OutputIterator,
   1028 	     typename _UniformRandomNumberGenerator>
   1029       void
   1030       arcsine_distribution<_RealType>::
   1031       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1032 		      _UniformRandomNumberGenerator& __urng,
   1033 		      const param_type& __p)
   1034       {
   1035 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1036 	    result_type>)
   1037 
   1038 	result_type __dif = __p.b() - __p.a();
   1039 	result_type __sum = __p.a() + __p.b();
   1040 	while (__f != __t)
   1041 	  {
   1042 	    result_type __x = std::sin(this->_M_ud(__urng));
   1043 	    *__f++ = (__x * __dif + __sum) / result_type(2);
   1044 	  }
   1045       }
   1046 
   1047   template<typename _RealType, typename _CharT, typename _Traits>
   1048     std::basic_ostream<_CharT, _Traits>&
   1049     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1050 	       const arcsine_distribution<_RealType>& __x)
   1051     {
   1052       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1053       typedef typename __ostream_type::ios_base    __ios_base;
   1054 
   1055       const typename __ios_base::fmtflags __flags = __os.flags();
   1056       const _CharT __fill = __os.fill();
   1057       const std::streamsize __precision = __os.precision();
   1058       const _CharT __space = __os.widen(' ');
   1059       __os.flags(__ios_base::scientific | __ios_base::left);
   1060       __os.fill(__space);
   1061       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1062 
   1063       __os << __x.a() << __space << __x.b();
   1064       __os << __space << __x._M_ud;
   1065 
   1066       __os.flags(__flags);
   1067       __os.fill(__fill);
   1068       __os.precision(__precision);
   1069       return __os;
   1070     }
   1071 
   1072   template<typename _RealType, typename _CharT, typename _Traits>
   1073     std::basic_istream<_CharT, _Traits>&
   1074     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1075 	       arcsine_distribution<_RealType>& __x)
   1076     {
   1077       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1078       typedef typename __istream_type::ios_base    __ios_base;
   1079 
   1080       const typename __ios_base::fmtflags __flags = __is.flags();
   1081       __is.flags(__ios_base::dec | __ios_base::skipws);
   1082 
   1083       _RealType __a, __b;
   1084       __is >> __a >> __b;
   1085       __is >> __x._M_ud;
   1086       __x.param(typename arcsine_distribution<_RealType>::
   1087 		param_type(__a, __b));
   1088 
   1089       __is.flags(__flags);
   1090       return __is;
   1091     }
   1092 
   1093 
   1094   template<typename _RealType>
   1095     template<typename _UniformRandomNumberGenerator>
   1096       typename hoyt_distribution<_RealType>::result_type
   1097       hoyt_distribution<_RealType>::
   1098       operator()(_UniformRandomNumberGenerator& __urng)
   1099       {
   1100 	result_type __x = this->_M_ad(__urng);
   1101 	result_type __y = this->_M_ed(__urng);
   1102 	return (result_type(2) * this->q()
   1103 		  / (result_type(1) + this->q() * this->q()))
   1104 	       * std::sqrt(this->omega() * __x * __y);
   1105       }
   1106 
   1107   template<typename _RealType>
   1108     template<typename _UniformRandomNumberGenerator>
   1109       typename hoyt_distribution<_RealType>::result_type
   1110       hoyt_distribution<_RealType>::
   1111       operator()(_UniformRandomNumberGenerator& __urng,
   1112 		 const param_type& __p)
   1113       {
   1114 	result_type __q2 = __p.q() * __p.q();
   1115 	result_type __num = result_type(0.5L) * (result_type(1) + __q2);
   1116 	typename __gnu_cxx::arcsine_distribution<result_type>::param_type
   1117 	  __pa(__num, __num / __q2);
   1118 	result_type __x = this->_M_ad(__pa, __urng);
   1119 	result_type __y = this->_M_ed(__urng);
   1120 	return (result_type(2) * __p.q() / (result_type(1) + __q2))
   1121 	       * std::sqrt(__p.omega() * __x * __y);
   1122       }
   1123 
   1124   template<typename _RealType>
   1125     template<typename _OutputIterator,
   1126 	     typename _UniformRandomNumberGenerator>
   1127       void
   1128       hoyt_distribution<_RealType>::
   1129       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1130 		      _UniformRandomNumberGenerator& __urng,
   1131 		      const param_type& __p)
   1132       {
   1133 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1134 	    result_type>)
   1135 
   1136 	result_type __2q = result_type(2) * __p.q();
   1137 	result_type __q2 = __p.q() * __p.q();
   1138 	result_type __q2p1 = result_type(1) + __q2;
   1139 	result_type __num = result_type(0.5L) * __q2p1;
   1140 	result_type __omega = __p.omega();
   1141 	typename __gnu_cxx::arcsine_distribution<result_type>::param_type
   1142 	  __pa(__num, __num / __q2);
   1143 	while (__f != __t)
   1144 	  {
   1145 	    result_type __x = this->_M_ad(__pa, __urng);
   1146 	    result_type __y = this->_M_ed(__urng);
   1147 	    *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y);
   1148 	  }
   1149       }
   1150 
   1151   template<typename _RealType, typename _CharT, typename _Traits>
   1152     std::basic_ostream<_CharT, _Traits>&
   1153     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1154 	       const hoyt_distribution<_RealType>& __x)
   1155     {
   1156       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1157       typedef typename __ostream_type::ios_base    __ios_base;
   1158 
   1159       const typename __ios_base::fmtflags __flags = __os.flags();
   1160       const _CharT __fill = __os.fill();
   1161       const std::streamsize __precision = __os.precision();
   1162       const _CharT __space = __os.widen(' ');
   1163       __os.flags(__ios_base::scientific | __ios_base::left);
   1164       __os.fill(__space);
   1165       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1166 
   1167       __os << __x.q() << __space << __x.omega();
   1168       __os << __space << __x._M_ad;
   1169       __os << __space << __x._M_ed;
   1170 
   1171       __os.flags(__flags);
   1172       __os.fill(__fill);
   1173       __os.precision(__precision);
   1174       return __os;
   1175     }
   1176 
   1177   template<typename _RealType, typename _CharT, typename _Traits>
   1178     std::basic_istream<_CharT, _Traits>&
   1179     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1180 	       hoyt_distribution<_RealType>& __x)
   1181     {
   1182       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1183       typedef typename __istream_type::ios_base    __ios_base;
   1184 
   1185       const typename __ios_base::fmtflags __flags = __is.flags();
   1186       __is.flags(__ios_base::dec | __ios_base::skipws);
   1187 
   1188       _RealType __q, __omega;
   1189       __is >> __q >> __omega;
   1190       __is >> __x._M_ad;
   1191       __is >> __x._M_ed;
   1192       __x.param(typename hoyt_distribution<_RealType>::
   1193 		param_type(__q, __omega));
   1194 
   1195       __is.flags(__flags);
   1196       return __is;
   1197     }
   1198 
   1199 
   1200   template<typename _RealType>
   1201     template<typename _OutputIterator,
   1202 	     typename _UniformRandomNumberGenerator>
   1203       void
   1204       triangular_distribution<_RealType>::
   1205       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1206 		      _UniformRandomNumberGenerator& __urng,
   1207 		      const param_type& __param)
   1208       {
   1209 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1210 	    result_type>)
   1211 
   1212 	while (__f != __t)
   1213 	  *__f++ = this->operator()(__urng, __param);
   1214       }
   1215 
   1216   template<typename _RealType, typename _CharT, typename _Traits>
   1217     std::basic_ostream<_CharT, _Traits>&
   1218     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1219 	       const __gnu_cxx::triangular_distribution<_RealType>& __x)
   1220     {
   1221       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1222       typedef typename __ostream_type::ios_base    __ios_base;
   1223 
   1224       const typename __ios_base::fmtflags __flags = __os.flags();
   1225       const _CharT __fill = __os.fill();
   1226       const std::streamsize __precision = __os.precision();
   1227       const _CharT __space = __os.widen(' ');
   1228       __os.flags(__ios_base::scientific | __ios_base::left);
   1229       __os.fill(__space);
   1230       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1231 
   1232       __os << __x.a() << __space << __x.b() << __space << __x.c();
   1233 
   1234       __os.flags(__flags);
   1235       __os.fill(__fill);
   1236       __os.precision(__precision);
   1237       return __os;
   1238     }
   1239 
   1240   template<typename _RealType, typename _CharT, typename _Traits>
   1241     std::basic_istream<_CharT, _Traits>&
   1242     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1243 	       __gnu_cxx::triangular_distribution<_RealType>& __x)
   1244     {
   1245       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1246       typedef typename __istream_type::ios_base    __ios_base;
   1247 
   1248       const typename __ios_base::fmtflags __flags = __is.flags();
   1249       __is.flags(__ios_base::dec | __ios_base::skipws);
   1250 
   1251       _RealType __a, __b, __c;
   1252       __is >> __a >> __b >> __c;
   1253       __x.param(typename __gnu_cxx::triangular_distribution<_RealType>::
   1254 		param_type(__a, __b, __c));
   1255 
   1256       __is.flags(__flags);
   1257       return __is;
   1258     }
   1259 
   1260 
   1261   template<typename _RealType>
   1262     template<typename _UniformRandomNumberGenerator>
   1263       typename von_mises_distribution<_RealType>::result_type
   1264       von_mises_distribution<_RealType>::
   1265       operator()(_UniformRandomNumberGenerator& __urng,
   1266 		 const param_type& __p)
   1267       {
   1268 	const result_type __pi
   1269 	  = __gnu_cxx::__math_constants<result_type>::__pi;
   1270 	std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
   1271 	  __aurng(__urng);
   1272 
   1273 	result_type __f;
   1274 	while (1)
   1275 	  {
   1276 	    result_type __rnd = std::cos(__pi * __aurng());
   1277 	    __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd);
   1278 	    result_type __c = __p._M_kappa * (__p._M_r - __f);
   1279 
   1280 	    result_type __rnd2 = __aurng();
   1281 	    if (__c * (result_type(2) - __c) > __rnd2)
   1282 	      break;
   1283 	    if (std::log(__c / __rnd2) >= __c - result_type(1))
   1284 	      break;
   1285 	  }
   1286 
   1287 	result_type __res = std::acos(__f);
   1288 #if _GLIBCXX_USE_C99_MATH_TR1
   1289 	__res = std::copysign(__res, __aurng() - result_type(0.5));
   1290 #else
   1291 	if (__aurng() < result_type(0.5))
   1292 	  __res = -__res;
   1293 #endif
   1294 	__res += __p._M_mu;
   1295 	if (__res > __pi)
   1296 	  __res -= result_type(2) * __pi;
   1297 	else if (__res < -__pi)
   1298 	  __res += result_type(2) * __pi;
   1299 	return __res;
   1300       }
   1301 
   1302   template<typename _RealType>
   1303     template<typename _OutputIterator,
   1304 	     typename _UniformRandomNumberGenerator>
   1305       void
   1306       von_mises_distribution<_RealType>::
   1307       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1308 		      _UniformRandomNumberGenerator& __urng,
   1309 		      const param_type& __param)
   1310       {
   1311 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1312 	    result_type>)
   1313 
   1314 	while (__f != __t)
   1315 	  *__f++ = this->operator()(__urng, __param);
   1316       }
   1317 
   1318   template<typename _RealType, typename _CharT, typename _Traits>
   1319     std::basic_ostream<_CharT, _Traits>&
   1320     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1321 	       const __gnu_cxx::von_mises_distribution<_RealType>& __x)
   1322     {
   1323       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1324       typedef typename __ostream_type::ios_base    __ios_base;
   1325 
   1326       const typename __ios_base::fmtflags __flags = __os.flags();
   1327       const _CharT __fill = __os.fill();
   1328       const std::streamsize __precision = __os.precision();
   1329       const _CharT __space = __os.widen(' ');
   1330       __os.flags(__ios_base::scientific | __ios_base::left);
   1331       __os.fill(__space);
   1332       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1333 
   1334       __os << __x.mu() << __space << __x.kappa();
   1335 
   1336       __os.flags(__flags);
   1337       __os.fill(__fill);
   1338       __os.precision(__precision);
   1339       return __os;
   1340     }
   1341 
   1342   template<typename _RealType, typename _CharT, typename _Traits>
   1343     std::basic_istream<_CharT, _Traits>&
   1344     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1345 	       __gnu_cxx::von_mises_distribution<_RealType>& __x)
   1346     {
   1347       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1348       typedef typename __istream_type::ios_base    __ios_base;
   1349 
   1350       const typename __ios_base::fmtflags __flags = __is.flags();
   1351       __is.flags(__ios_base::dec | __ios_base::skipws);
   1352 
   1353       _RealType __mu, __kappa;
   1354       __is >> __mu >> __kappa;
   1355       __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>::
   1356 		param_type(__mu, __kappa));
   1357 
   1358       __is.flags(__flags);
   1359       return __is;
   1360     }
   1361 
   1362 
   1363   template<typename _UIntType>
   1364     template<typename _UniformRandomNumberGenerator>
   1365       typename hypergeometric_distribution<_UIntType>::result_type
   1366       hypergeometric_distribution<_UIntType>::
   1367       operator()(_UniformRandomNumberGenerator& __urng,
   1368 		 const param_type& __param)
   1369       {
   1370 	std::__detail::_Adaptor<_UniformRandomNumberGenerator, double>
   1371 	  __aurng(__urng);
   1372 
   1373 	result_type __a = __param.successful_size();
   1374 	result_type __b = __param.total_size();
   1375 	result_type __k = 0;
   1376 
   1377 	if (__param.total_draws() < __param.total_size() / 2)
   1378 	  {
   1379 	    for (result_type __i = 0; __i < __param.total_draws(); ++__i)
   1380 	      {
   1381 		if (__b * __aurng() < __a)
   1382 		  {
   1383 		    ++__k;
   1384 		    if (__k == __param.successful_size())
   1385 		      return __k;
   1386 		   --__a;
   1387 		  }
   1388 		--__b;
   1389 	      }
   1390 	    return __k;
   1391 	  }
   1392 	else
   1393 	  {
   1394 	    for (result_type __i = 0; __i < __param.unsuccessful_size(); ++__i)
   1395 	      {
   1396 		if (__b * __aurng() < __a)
   1397 		  {
   1398 		    ++__k;
   1399 		    if (__k == __param.successful_size())
   1400 		      return __param.successful_size() - __k;
   1401 		    --__a;
   1402 		  }
   1403 		--__b;
   1404 	      }
   1405 	    return __param.successful_size() - __k;
   1406 	  }
   1407       }
   1408 
   1409   template<typename _UIntType>
   1410     template<typename _OutputIterator,
   1411 	     typename _UniformRandomNumberGenerator>
   1412       void
   1413       hypergeometric_distribution<_UIntType>::
   1414       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1415 		      _UniformRandomNumberGenerator& __urng,
   1416 		      const param_type& __param)
   1417       {
   1418 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1419 	    result_type>)
   1420 
   1421 	while (__f != __t)
   1422 	  *__f++ = this->operator()(__urng);
   1423       }
   1424 
   1425   template<typename _UIntType, typename _CharT, typename _Traits>
   1426     std::basic_ostream<_CharT, _Traits>&
   1427     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1428 	       const __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
   1429     {
   1430       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1431       typedef typename __ostream_type::ios_base    __ios_base;
   1432 
   1433       const typename __ios_base::fmtflags __flags = __os.flags();
   1434       const _CharT __fill = __os.fill();
   1435       const std::streamsize __precision = __os.precision();
   1436       const _CharT __space = __os.widen(' ');
   1437       __os.flags(__ios_base::scientific | __ios_base::left);
   1438       __os.fill(__space);
   1439       __os.precision(std::numeric_limits<_UIntType>::max_digits10);
   1440 
   1441       __os << __x.total_size() << __space << __x.successful_size() << __space
   1442 	   << __x.total_draws();
   1443 
   1444       __os.flags(__flags);
   1445       __os.fill(__fill);
   1446       __os.precision(__precision);
   1447       return __os;
   1448     }
   1449 
   1450   template<typename _UIntType, typename _CharT, typename _Traits>
   1451     std::basic_istream<_CharT, _Traits>&
   1452     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1453 	       __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
   1454     {
   1455       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1456       typedef typename __istream_type::ios_base    __ios_base;
   1457 
   1458       const typename __ios_base::fmtflags __flags = __is.flags();
   1459       __is.flags(__ios_base::dec | __ios_base::skipws);
   1460 
   1461       _UIntType __total_size, __successful_size, __total_draws;
   1462       __is >> __total_size >> __successful_size >> __total_draws;
   1463       __x.param(typename __gnu_cxx::hypergeometric_distribution<_UIntType>::
   1464 		param_type(__total_size, __successful_size, __total_draws));
   1465 
   1466       __is.flags(__flags);
   1467       return __is;
   1468     }
   1469 
   1470 
   1471   template<typename _RealType>
   1472     template<typename _UniformRandomNumberGenerator>
   1473       typename logistic_distribution<_RealType>::result_type
   1474       logistic_distribution<_RealType>::
   1475       operator()(_UniformRandomNumberGenerator& __urng,
   1476 		 const param_type& __p)
   1477       {
   1478 	std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
   1479 	  __aurng(__urng);
   1480 
   1481 	result_type __arg = result_type(1);
   1482 	while (__arg == result_type(1) || __arg == result_type(0))
   1483 	  __arg = __aurng();
   1484 	return __p.a()
   1485 	     + __p.b() * std::log(__arg / (result_type(1) - __arg));
   1486       }
   1487 
   1488   template<typename _RealType>
   1489     template<typename _OutputIterator,
   1490 	     typename _UniformRandomNumberGenerator>
   1491       void
   1492       logistic_distribution<_RealType>::
   1493       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1494 		      _UniformRandomNumberGenerator& __urng,
   1495 		      const param_type& __p)
   1496       {
   1497 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1498 	    result_type>)
   1499 
   1500 	std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
   1501 	  __aurng(__urng);
   1502 
   1503 	while (__f != __t)
   1504 	  {
   1505 	    result_type __arg = result_type(1);
   1506 	    while (__arg == result_type(1) || __arg == result_type(0))
   1507 	      __arg = __aurng();
   1508 	    *__f++ = __p.a()
   1509 		   + __p.b() * std::log(__arg / (result_type(1) - __arg));
   1510 	  }
   1511       }
   1512 
   1513   template<typename _RealType, typename _CharT, typename _Traits>
   1514     std::basic_ostream<_CharT, _Traits>&
   1515     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1516 	       const logistic_distribution<_RealType>& __x)
   1517     {
   1518       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1519       typedef typename __ostream_type::ios_base    __ios_base;
   1520 
   1521       const typename __ios_base::fmtflags __flags = __os.flags();
   1522       const _CharT __fill = __os.fill();
   1523       const std::streamsize __precision = __os.precision();
   1524       const _CharT __space = __os.widen(' ');
   1525       __os.flags(__ios_base::scientific | __ios_base::left);
   1526       __os.fill(__space);
   1527       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1528 
   1529       __os << __x.a() << __space << __x.b();
   1530 
   1531       __os.flags(__flags);
   1532       __os.fill(__fill);
   1533       __os.precision(__precision);
   1534       return __os;
   1535     }
   1536 
   1537   template<typename _RealType, typename _CharT, typename _Traits>
   1538     std::basic_istream<_CharT, _Traits>&
   1539     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1540 	       logistic_distribution<_RealType>& __x)
   1541     {
   1542       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1543       typedef typename __istream_type::ios_base    __ios_base;
   1544 
   1545       const typename __ios_base::fmtflags __flags = __is.flags();
   1546       __is.flags(__ios_base::dec | __ios_base::skipws);
   1547 
   1548       _RealType __a, __b;
   1549       __is >> __a >> __b;
   1550       __x.param(typename logistic_distribution<_RealType>::
   1551 		param_type(__a, __b));
   1552 
   1553       __is.flags(__flags);
   1554       return __is;
   1555     }
   1556 
   1557 
   1558   namespace {
   1559 
   1560     // Helper class for the uniform_on_sphere_distribution generation
   1561     // function.
   1562     template<std::size_t _Dimen, typename _RealType>
   1563       class uniform_on_sphere_helper
   1564       {
   1565 	typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>::
   1566 	  result_type result_type;
   1567 
   1568       public:
   1569 	template<typename _NormalDistribution,
   1570 		 typename _UniformRandomNumberGenerator>
   1571 	result_type operator()(_NormalDistribution& __nd,
   1572 			       _UniformRandomNumberGenerator& __urng)
   1573         {
   1574 	  result_type __ret;
   1575 	  typename result_type::value_type __norm;
   1576 
   1577 	  do
   1578 	    {
   1579 	      auto __sum = _RealType(0);
   1580 
   1581 	      std::generate(__ret.begin(), __ret.end(),
   1582 			    [&__nd, &__urng, &__sum](){
   1583 			      _RealType __t = __nd(__urng);
   1584 			      __sum += __t * __t;
   1585 			      return __t; });
   1586 	      __norm = std::sqrt(__sum);
   1587 	    }
   1588 	  while (__norm == _RealType(0) || ! __builtin_isfinite(__norm));
   1589 
   1590 	  std::transform(__ret.begin(), __ret.end(), __ret.begin(),
   1591 			 [__norm](_RealType __val){ return __val / __norm; });
   1592 
   1593 	  return __ret;
   1594         }
   1595       };
   1596 
   1597 
   1598     template<typename _RealType>
   1599       class uniform_on_sphere_helper<2, _RealType>
   1600       {
   1601 	typedef typename uniform_on_sphere_distribution<2, _RealType>::
   1602 	  result_type result_type;
   1603 
   1604       public:
   1605 	template<typename _NormalDistribution,
   1606 		 typename _UniformRandomNumberGenerator>
   1607 	result_type operator()(_NormalDistribution&,
   1608 			       _UniformRandomNumberGenerator& __urng)
   1609         {
   1610 	  result_type __ret;
   1611 	  _RealType __sq;
   1612 	  std::__detail::_Adaptor<_UniformRandomNumberGenerator,
   1613 				  _RealType> __aurng(__urng);
   1614 
   1615 	  do
   1616 	    {
   1617 	      __ret[0] = _RealType(2) * __aurng() - _RealType(1);
   1618 	      __ret[1] = _RealType(2) * __aurng() - _RealType(1);
   1619 
   1620 	      __sq = __ret[0] * __ret[0] + __ret[1] * __ret[1];
   1621 	    }
   1622 	  while (__sq == _RealType(0) || __sq > _RealType(1));
   1623 
   1624 #if _GLIBCXX_USE_C99_MATH_TR1
   1625 	  // Yes, we do not just use sqrt(__sq) because hypot() is more
   1626 	  // accurate.
   1627 	  auto __norm = std::hypot(__ret[0], __ret[1]);
   1628 #else
   1629 	  auto __norm = std::sqrt(__sq);
   1630 #endif
   1631 	  __ret[0] /= __norm;
   1632 	  __ret[1] /= __norm;
   1633 
   1634 	  return __ret;
   1635         }
   1636       };
   1637 
   1638   }
   1639 
   1640 
   1641   template<std::size_t _Dimen, typename _RealType>
   1642     template<typename _UniformRandomNumberGenerator>
   1643       typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type
   1644       uniform_on_sphere_distribution<_Dimen, _RealType>::
   1645       operator()(_UniformRandomNumberGenerator& __urng,
   1646 		 const param_type& __p)
   1647       {
   1648         uniform_on_sphere_helper<_Dimen, _RealType> __helper;
   1649         return __helper(_M_nd, __urng);
   1650       }
   1651 
   1652   template<std::size_t _Dimen, typename _RealType>
   1653     template<typename _OutputIterator,
   1654 	     typename _UniformRandomNumberGenerator>
   1655       void
   1656       uniform_on_sphere_distribution<_Dimen, _RealType>::
   1657       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1658 		      _UniformRandomNumberGenerator& __urng,
   1659 		      const param_type& __param)
   1660       {
   1661 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1662 	    result_type>)
   1663 
   1664 	while (__f != __t)
   1665 	  *__f++ = this->operator()(__urng, __param);
   1666       }
   1667 
   1668   template<std::size_t _Dimen, typename _RealType, typename _CharT,
   1669 	   typename _Traits>
   1670     std::basic_ostream<_CharT, _Traits>&
   1671     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1672 	       const __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
   1673 							       _RealType>& __x)
   1674     {
   1675       return __os << __x._M_nd;
   1676     }
   1677 
   1678   template<std::size_t _Dimen, typename _RealType, typename _CharT,
   1679 	   typename _Traits>
   1680     std::basic_istream<_CharT, _Traits>&
   1681     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1682 	       __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
   1683 							 _RealType>& __x)
   1684     {
   1685       return __is >> __x._M_nd;
   1686     }
   1687 
   1688 
   1689   namespace {
   1690 
   1691     // Helper class for the uniform_inside_sphere_distribution generation
   1692     // function.
   1693     template<std::size_t _Dimen, bool _SmallDimen, typename _RealType>
   1694       class uniform_inside_sphere_helper;
   1695 
   1696     template<std::size_t _Dimen, typename _RealType>
   1697       class uniform_inside_sphere_helper<_Dimen, false, _RealType>
   1698       {
   1699 	using result_type
   1700 	  = typename uniform_inside_sphere_distribution<_Dimen, _RealType>::
   1701 	    result_type;
   1702 
   1703       public:
   1704 	template<typename _UniformOnSphereDistribution,
   1705 		 typename _UniformRandomNumberGenerator>
   1706 	result_type
   1707 	operator()(_UniformOnSphereDistribution& __uosd,
   1708 		   _UniformRandomNumberGenerator& __urng,
   1709 		   _RealType __radius)
   1710         {
   1711 	  std::__detail::_Adaptor<_UniformRandomNumberGenerator,
   1712 				  _RealType> __aurng(__urng);
   1713 
   1714 	  _RealType __pow = 1 / _RealType(_Dimen);
   1715 	  _RealType __urt = __radius * std::pow(__aurng(), __pow);
   1716 	  result_type __ret = __uosd(__aurng);
   1717 
   1718 	  std::transform(__ret.begin(), __ret.end(), __ret.begin(),
   1719 			 [__urt](_RealType __val)
   1720 			 { return __val * __urt; });
   1721 
   1722 	  return __ret;
   1723         }
   1724       };
   1725 
   1726     // Helper class for the uniform_inside_sphere_distribution generation
   1727     // function specialized for small dimensions.
   1728     template<std::size_t _Dimen, typename _RealType>
   1729       class uniform_inside_sphere_helper<_Dimen, true, _RealType>
   1730       {
   1731 	using result_type
   1732 	  = typename uniform_inside_sphere_distribution<_Dimen, _RealType>::
   1733 	    result_type;
   1734 
   1735       public:
   1736 	template<typename _UniformOnSphereDistribution,
   1737 		 typename _UniformRandomNumberGenerator>
   1738 	result_type
   1739 	operator()(_UniformOnSphereDistribution&,
   1740 		   _UniformRandomNumberGenerator& __urng,
   1741 		   _RealType __radius)
   1742         {
   1743 	  result_type __ret;
   1744 	  _RealType __sq;
   1745 	  _RealType __radsq = __radius * __radius;
   1746 	  std::__detail::_Adaptor<_UniformRandomNumberGenerator,
   1747 				  _RealType> __aurng(__urng);
   1748 
   1749 	  do
   1750 	    {
   1751 	      __sq = _RealType(0);
   1752 	      for (int i = 0; i < _Dimen; ++i)
   1753 		{
   1754 		  __ret[i] = _RealType(2) * __aurng() - _RealType(1);
   1755 		  __sq += __ret[i] * __ret[i];
   1756 		}
   1757 	    }
   1758 	  while (__sq > _RealType(1));
   1759 
   1760 	  for (int i = 0; i < _Dimen; ++i)
   1761             __ret[i] *= __radius;
   1762 
   1763 	  return __ret;
   1764         }
   1765       };
   1766   } // namespace
   1767 
   1768   //
   1769   //  Experiments have shown that rejection is more efficient than transform
   1770   //  for dimensions less than 8.
   1771   //
   1772   template<std::size_t _Dimen, typename _RealType>
   1773     template<typename _UniformRandomNumberGenerator>
   1774       typename uniform_inside_sphere_distribution<_Dimen, _RealType>::result_type
   1775       uniform_inside_sphere_distribution<_Dimen, _RealType>::
   1776       operator()(_UniformRandomNumberGenerator& __urng,
   1777 		 const param_type& __p)
   1778       {
   1779         uniform_inside_sphere_helper<_Dimen, _Dimen < 8, _RealType> __helper;
   1780         return __helper(_M_uosd, __urng, __p.radius());
   1781       }
   1782 
   1783   template<std::size_t _Dimen, typename _RealType>
   1784     template<typename _OutputIterator,
   1785 	     typename _UniformRandomNumberGenerator>
   1786       void
   1787       uniform_inside_sphere_distribution<_Dimen, _RealType>::
   1788       __generate_impl(_OutputIterator __f, _OutputIterator __t,
   1789 		      _UniformRandomNumberGenerator& __urng,
   1790 		      const param_type& __param)
   1791       {
   1792 	__glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
   1793 	    result_type>)
   1794 
   1795 	while (__f != __t)
   1796 	  *__f++ = this->operator()(__urng, __param);
   1797       }
   1798 
   1799   template<std::size_t _Dimen, typename _RealType, typename _CharT,
   1800 	   typename _Traits>
   1801     std::basic_ostream<_CharT, _Traits>&
   1802     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   1803 	       const __gnu_cxx::uniform_inside_sphere_distribution<_Dimen,
   1804 								_RealType>& __x)
   1805     {
   1806       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
   1807       typedef typename __ostream_type::ios_base    __ios_base;
   1808 
   1809       const typename __ios_base::fmtflags __flags = __os.flags();
   1810       const _CharT __fill = __os.fill();
   1811       const std::streamsize __precision = __os.precision();
   1812       const _CharT __space = __os.widen(' ');
   1813       __os.flags(__ios_base::scientific | __ios_base::left);
   1814       __os.fill(__space);
   1815       __os.precision(std::numeric_limits<_RealType>::max_digits10);
   1816 
   1817       __os << __x.radius() << __space << __x._M_uosd;
   1818 
   1819       __os.flags(__flags);
   1820       __os.fill(__fill);
   1821       __os.precision(__precision);
   1822 
   1823       return __os;
   1824     }
   1825 
   1826   template<std::size_t _Dimen, typename _RealType, typename _CharT,
   1827 	   typename _Traits>
   1828     std::basic_istream<_CharT, _Traits>&
   1829     operator>>(std::basic_istream<_CharT, _Traits>& __is,
   1830 	       __gnu_cxx::uniform_inside_sphere_distribution<_Dimen,
   1831 							     _RealType>& __x)
   1832     {
   1833       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
   1834       typedef typename __istream_type::ios_base    __ios_base;
   1835 
   1836       const typename __ios_base::fmtflags __flags = __is.flags();
   1837       __is.flags(__ios_base::dec | __ios_base::skipws);
   1838 
   1839       _RealType __radius_val;
   1840       __is >> __radius_val >> __x._M_uosd;
   1841       __x.param(typename uniform_inside_sphere_distribution<_Dimen, _RealType>::
   1842 		param_type(__radius_val));
   1843 
   1844       __is.flags(__flags);
   1845 
   1846       return __is;
   1847     }
   1848 
   1849 _GLIBCXX_END_NAMESPACE_VERSION
   1850 } // namespace __gnu_cxx
   1851 
   1852 
   1853 #endif // _EXT_RANDOM_TCC
   1854