libstdc++-v3/include/ext/random.tcc

   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 _GLIBCXX_END_NAMESPACE_VERSION
1191 } // namespace
1192
1193
1194 #endif // _EXT_RANDOM_TCC