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

   1 // Random number extensions -*- C++ -*-
   2
   3 // Copyright (C) 2012 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
  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       typename std::enable_if<std::is_class<_Sseq>::value>::type
  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       {
  96         size_t __lag;
  97
  98         if (_M_nstate32 >= 623)
  99           __lag = 11;
 100         else if (_M_nstate32 >= 68)
 101           __lag = 7;
 102         else if (_M_nstate32 >= 39)
 103           __lag = 5;
 104         else
 105           __lag = 3;
 106         const size_t __mid = (_M_nstate32 - __lag) / 2;
 107
 108         std::fill(_M_state32, _M_state32 + _M_nstate32, UINT32_C(0x8b8b8b8b));
 109         uint32_t __arr[_M_nstate32];
 110         __q.generate(__arr + 0, __arr + _M_nstate32);
 111
 112         uint32_t __r = _Func1(_M_state32[0] ^ _M_state32[__mid]
 113                               ^ _M_state32[_M_nstate32  - 1]);
 114         _M_state32[__mid] += __r;
 115         __r += _M_nstate32;
 116         _M_state32[__mid + __lag] += __r;
 117         _M_state32[0] = __r;
 118
 119         for (size_t __i = 1, __j = 0; __j < _M_nstate32; ++__j)
 120           {
 121             __r = _Func1(_M_state32[__i]
 122                          ^ _M_state32[(__i + __mid) % _M_nstate32]
 123                          ^ _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
 124             _M_state32[(__i + __mid) % _M_nstate32] += __r;
 125             __r += __arr[__j] + __i;
 126             _M_state32[(__i + __mid + __lag) % _M_nstate32] += __r;
 127             _M_state32[__i] = __r;
 128             __i = (__i + 1) % _M_nstate32;
 129           }
 130         for (size_t __j = 0; __j < _M_nstate32; ++__j)
 131           {
 132             const size_t __i = (__j + 1) % _M_nstate32;
 133             __r = _Func2(_M_state32[__i]
 134                          + _M_state32[(__i + __mid) % _M_nstate32]
 135                          + _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
 136             _M_state32[(__i + __mid) % _M_nstate32] ^= __r;
 137             __r -= __i;
 138             _M_state32[(__i + __mid + __lag) % _M_nstate32] ^= __r;
 139             _M_state32[__i] = __r;
 140           }
 141
 142         _M_pos = state_size;
 143         _M_period_certification();
 144       }
 145
 146
 147   template<typename _UIntType, size_t __m,
 148            size_t __pos1, size_t __sl1, size_t __sl2,
 149            size_t __sr1, size_t __sr2,
 150            uint32_t __msk1, uint32_t __msk2,
 151            uint32_t __msk3, uint32_t __msk4,
 152            uint32_t __parity1, uint32_t __parity2,
 153            uint32_t __parity3, uint32_t __parity4>
 154     void simd_fast_mersenne_twister_engine<_UIntType, __m,
 155                                            __pos1, __sl1, __sl2, __sr1, __sr2,
 156                                            __msk1, __msk2, __msk3, __msk4,
 157                                            __parity1, __parity2, __parity3,
 158                                            __parity4>::
 159     _M_period_certification(void)
 160     {
 161       static const uint32_t __parity[4] = { __parity1, __parity2,
 162                                             __parity3, __parity4 };
 163       uint32_t __inner = 0;
 164       for (size_t __i = 0; __i < 4; ++__i)
 165         if (__parity[__i] != 0)
 166           __inner ^= _M_state32[__i] & __parity[__i];
 167
 168       if (__builtin_parity(__inner) & 1)
 169         return;
 170       for (size_t __i = 0; __i < 4; ++__i)
 171         if (__parity[__i] != 0)
 172           {
 173             _M_state32[__i] ^= 1 << (__builtin_ffs(__parity[__i]) - 1);
 174             return;
 175           }
 176       __builtin_unreachable();
 177     }
 178
 179
 180   template<typename _UIntType, size_t __m,
 181            size_t __pos1, size_t __sl1, size_t __sl2,
 182            size_t __sr1, size_t __sr2,
 183            uint32_t __msk1, uint32_t __msk2,
 184            uint32_t __msk3, uint32_t __msk4,
 185            uint32_t __parity1, uint32_t __parity2,
 186            uint32_t __parity3, uint32_t __parity4>
 187     void simd_fast_mersenne_twister_engine<_UIntType, __m,
 188                                            __pos1, __sl1, __sl2, __sr1, __sr2,
 189                                            __msk1, __msk2, __msk3, __msk4,
 190                                            __parity1, __parity2, __parity3,
 191                                            __parity4>::
 192     discard(unsigned long long __z)
 193     {
 194       while (__z > state_size - _M_pos)
 195         {
 196           __z -= state_size - _M_pos;
 197
 198           _M_gen_rand();
 199         }
 200
 201       _M_pos += __z;
 202     }
 203
 204
 205 #ifndef  _GLIBCXX_OPT_HAVE_RANDOM_SFMT_GEN_READ
 206
 207   namespace {
 208
 209     template<size_t __shift>
 210       inline void __rshift(uint32_t *__out, const uint32_t *__in)
 211       {
 212         uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
 213                          | static_cast<uint64_t>(__in[2]));
 214         uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
 215                          | static_cast<uint64_t>(__in[0]));
 216
 217         uint64_t __oh = __th >> (__shift * 8);
 218         uint64_t __ol = __tl >> (__shift * 8);
 219         __ol |= __th << (64 - __shift * 8);
 220         __out[1] = static_cast<uint32_t>(__ol >> 32);
 221         __out[0] = static_cast<uint32_t>(__ol);
 222         __out[3] = static_cast<uint32_t>(__oh >> 32);
 223         __out[2] = static_cast<uint32_t>(__oh);
 224       }
 225
 226
 227     template<size_t __shift>
 228       inline void __lshift(uint32_t *__out, const uint32_t *__in)
 229       {
 230         uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
 231                          | static_cast<uint64_t>(__in[2]));
 232         uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
 233                          | static_cast<uint64_t>(__in[0]));
 234
 235         uint64_t __oh = __th << (__shift * 8);
 236         uint64_t __ol = __tl << (__shift * 8);
 237         __oh |= __tl >> (64 - __shift * 8);
 238         __out[1] = static_cast<uint32_t>(__ol >> 32);
 239         __out[0] = static_cast<uint32_t>(__ol);
 240         __out[3] = static_cast<uint32_t>(__oh >> 32);
 241         __out[2] = static_cast<uint32_t>(__oh);
 242       }
 243
 244
 245     template<size_t __sl1, size_t __sl2, size_t __sr1, size_t __sr2,
 246              uint32_t __msk1, uint32_t __msk2, uint32_t __msk3, uint32_t __msk4>
 247       inline void __recursion(uint32_t *__r,
 248                               const uint32_t *__a, const uint32_t *__b,
 249                               const uint32_t *__c, const uint32_t *__d)
 250       {
 251         uint32_t __x[4];
 252         uint32_t __y[4];
 253
 254         __lshift<__sl2>(__x, __a);
 255         __rshift<__sr2>(__y, __c);
 256         __r[0] = (__a[0] ^ __x[0] ^ ((__b[0] >> __sr1) & __msk1)
 257                   ^ __y[0] ^ (__d[0] << __sl1));
 258         __r[1] = (__a[1] ^ __x[1] ^ ((__b[1] >> __sr1) & __msk2)
 259                   ^ __y[1] ^ (__d[1] << __sl1));
 260         __r[2] = (__a[2] ^ __x[2] ^ ((__b[2] >> __sr1) & __msk3)
 261                   ^ __y[2] ^ (__d[2] << __sl1));
 262         __r[3] = (__a[3] ^ __x[3] ^ ((__b[3] >> __sr1) & __msk4)
 263                   ^ __y[3] ^ (__d[3] << __sl1));
 264       }
 265
 266   }
 267
 268
 269   template<typename _UIntType, size_t __m,
 270            size_t __pos1, size_t __sl1, size_t __sl2,
 271            size_t __sr1, size_t __sr2,
 272            uint32_t __msk1, uint32_t __msk2,
 273            uint32_t __msk3, uint32_t __msk4,
 274            uint32_t __parity1, uint32_t __parity2,
 275            uint32_t __parity3, uint32_t __parity4>
 276     void simd_fast_mersenne_twister_engine<_UIntType, __m,
 277                                            __pos1, __sl1, __sl2, __sr1, __sr2,
 278                                            __msk1, __msk2, __msk3, __msk4,
 279                                            __parity1, __parity2, __parity3,
 280                                            __parity4>::
 281     _M_gen_rand(void)
 282     {
 283       const uint32_t *__r1 = &_M_state32[_M_nstate32 - 8];
 284       const uint32_t *__r2 = &_M_state32[_M_nstate32 - 4];
 285       static constexpr size_t __pos1_32 = __pos1 * 4;
 286
 287       size_t __i;
 288       for (__i = 0; __i < _M_nstate32 - __pos1_32; __i += 4)
 289         {
 290           __recursion<__sl1, __sl2, __sr1, __sr2,
 291                       __msk1, __msk2, __msk3, __msk4>
 292             (&_M_state32[__i], &_M_state32[__i],
 293              &_M_state32[__i + __pos1_32], __r1, __r2);
 294           __r1 = __r2;
 295           __r2 = &_M_state32[__i];
 296         }
 297
 298       for (; __i < _M_nstate32; __i += 4)
 299         {
 300           __recursion<__sl1, __sl2, __sr1, __sr2,
 301                       __msk1, __msk2, __msk3, __msk4>
 302             (&_M_state32[__i], &_M_state32[__i],
 303              &_M_state32[__i + __pos1_32 - _M_nstate32], __r1, __r2);
 304           __r1 = __r2;
 305           __r2 = &_M_state32[__i];
 306         }
 307
 308       _M_pos = 0;
 309     }
 310
 311 #endif
 312
 313 #ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_OPERATOREQUAL
 314   template<typename _UIntType, size_t __m,
 315            size_t __pos1, size_t __sl1, size_t __sl2,
 316            size_t __sr1, size_t __sr2,
 317            uint32_t __msk1, uint32_t __msk2,
 318            uint32_t __msk3, uint32_t __msk4,
 319            uint32_t __parity1, uint32_t __parity2,
 320            uint32_t __parity3, uint32_t __parity4>
 321     bool
 322     operator==(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
 323                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
 324                __msk1, __msk2, __msk3, __msk4,
 325                __parity1, __parity2, __parity3, __parity4>& __lhs,
 326                const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
 327                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
 328                __msk1, __msk2, __msk3, __msk4,
 329                __parity1, __parity2, __parity3, __parity4>& __rhs)
 330     {
 331       typedef __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
 332                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
 333                __msk1, __msk2, __msk3, __msk4,
 334                __parity1, __parity2, __parity3, __parity4> __engine;
 335       return (std::equal(__lhs._M_stateT,
 336                          __lhs._M_stateT + __engine::state_size,
 337                          __rhs._M_stateT)
 338               && __lhs._M_pos == __rhs._M_pos);
 339     }
 340 #endif
 341
 342   template<typename _UIntType, size_t __m,
 343            size_t __pos1, size_t __sl1, size_t __sl2,
 344            size_t __sr1, size_t __sr2,
 345            uint32_t __msk1, uint32_t __msk2,
 346            uint32_t __msk3, uint32_t __msk4,
 347            uint32_t __parity1, uint32_t __parity2,
 348            uint32_t __parity3, uint32_t __parity4,
 349            typename _CharT, typename _Traits>
 350     std::basic_ostream<_CharT, _Traits>&
 351     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 352                const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
 353                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
 354                __msk1, __msk2, __msk3, __msk4,
 355                __parity1, __parity2, __parity3, __parity4>& __x)
 356     {
 357       typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
 358       typedef typename __ostream_type::ios_base __ios_base;
 359
 360       const typename __ios_base::fmtflags __flags = __os.flags();
 361       const _CharT __fill = __os.fill();
 362       const _CharT __space = __os.widen(' ');
 363       __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
 364       __os.fill(__space);
 365
 366       for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
 367         __os << __x._M_state32[__i] << __space;
 368       __os << __x._M_pos;
 369
 370       __os.flags(__flags);
 371       __os.fill(__fill);
 372       return __os;
 373     }
 374
 375
 376   template<typename _UIntType, size_t __m,
 377            size_t __pos1, size_t __sl1, size_t __sl2,
 378            size_t __sr1, size_t __sr2,
 379            uint32_t __msk1, uint32_t __msk2,
 380            uint32_t __msk3, uint32_t __msk4,
 381            uint32_t __parity1, uint32_t __parity2,
 382            uint32_t __parity3, uint32_t __parity4,
 383            typename _CharT, typename _Traits>
 384     std::basic_istream<_CharT, _Traits>&
 385     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 386                __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
 387                __m, __pos1, __sl1, __sl2, __sr1, __sr2,
 388                __msk1, __msk2, __msk3, __msk4,
 389                __parity1, __parity2, __parity3, __parity4>& __x)
 390     {
 391       typedef std::basic_istream<_CharT, _Traits> __istream_type;
 392       typedef typename __istream_type::ios_base __ios_base;
 393
 394       const typename __ios_base::fmtflags __flags = __is.flags();
 395       __is.flags(__ios_base::dec | __ios_base::skipws);
 396
 397       for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
 398         __is >> __x._M_state32[__i];
 399       __is >> __x._M_pos;
 400
 401       __is.flags(__flags);
 402       return __is;
 403     }
 404
 405
 406   /**
 407    * Iteration method due to M.D. J<o:>hnk.
 408    *
 409    * M.D. J<o:>hnk, Erzeugung von betaverteilten und gammaverteilten
 410    * Zufallszahlen, Metrika, Volume 8, 1964
 411    */
 412   template<typename _RealType>
 413     template<typename _UniformRandomNumberGenerator>
 414       typename beta_distribution<_RealType>::result_type
 415       beta_distribution<_RealType>::
 416       operator()(_UniformRandomNumberGenerator& __urng,
 417                  const param_type& __param)
 418       {
 419         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
 420           __aurng(__urng);
 421
 422         result_type __x, __y;
 423         do
 424           {
 425             __x = std::exp(std::log(__aurng()) / __param.alpha());
 426             __y = std::exp(std::log(__aurng()) / __param.beta());
 427           }
 428         while (__x + __y > result_type(1));
 429
 430         return __x / (__x + __y);
 431       }
 432
 433   template<typename _RealType>
 434     template<typename _OutputIterator,
 435              typename _UniformRandomNumberGenerator>
 436       void
 437       beta_distribution<_RealType>::
 438       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 439                       _UniformRandomNumberGenerator& __urng,
 440                       const param_type& __param)
 441       {
 442         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 443
 444         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
 445           __aurng(__urng);
 446
 447         while (__f != __t)
 448           {
 449             result_type __x, __y;
 450             do
 451               {
 452                 __x = std::exp(std::log(__aurng()) / __param.alpha());
 453                 __y = std::exp(std::log(__aurng()) / __param.beta());
 454               }
 455             while (__x + __y > result_type(1));
 456
 457             *__f++ = __x / (__x + __y);
 458           }
 459       }
 460
 461   template<typename _RealType, typename _CharT, typename _Traits>
 462     std::basic_ostream<_CharT, _Traits>&
 463     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 464                const __gnu_cxx::beta_distribution<_RealType>& __x)
 465     {
 466       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 467       typedef typename __ostream_type::ios_base    __ios_base;
 468
 469       const typename __ios_base::fmtflags __flags = __os.flags();
 470       const _CharT __fill = __os.fill();
 471       const std::streamsize __precision = __os.precision();
 472       const _CharT __space = __os.widen(' ');
 473       __os.flags(__ios_base::scientific | __ios_base::left);
 474       __os.fill(__space);
 475       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 476
 477       __os << __x.alpha() << __space << __x.beta();
 478
 479       __os.flags(__flags);
 480       __os.fill(__fill);
 481       __os.precision(__precision);
 482       return __os;
 483     }
 484
 485   template<typename _RealType, typename _CharT, typename _Traits>
 486     std::basic_istream<_CharT, _Traits>&
 487     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 488                __gnu_cxx::beta_distribution<_RealType>& __x)
 489     {
 490       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 491       typedef typename __istream_type::ios_base    __ios_base;
 492
 493       const typename __ios_base::fmtflags __flags = __is.flags();
 494       __is.flags(__ios_base::dec | __ios_base::skipws);
 495
 496       _RealType __alpha_val, __beta_val;
 497       __is >> __alpha_val >> __beta_val;
 498       __x.param(typename __gnu_cxx::beta_distribution<_RealType>::
 499                 param_type(__alpha_val, __beta_val));
 500
 501       __is.flags(__flags);
 502       return __is;
 503     }
 504
 505
 506   template<std::size_t _Dimen, typename _RealType>
 507     template<typename _InputIterator1, typename _InputIterator2>
 508       void
 509       normal_mv_distribution<_Dimen, _RealType>::param_type::
 510       _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
 511                    _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
 512       {
 513         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
 514         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
 515         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
 516                   _M_mean.end(), _RealType(0));
 517
 518         // Perform the Cholesky decomposition
 519         auto __w = _M_t.begin();
 520         for (size_t __j = 0; __j < _Dimen; ++__j)
 521           {
 522             _RealType __sum = _RealType(0);
 523
 524             auto __slitbegin = __w;
 525             auto __cit = _M_t.begin();
 526             for (size_t __i = 0; __i < __j; ++__i)
 527               {
 528                 auto __slit = __slitbegin;
 529                 _RealType __s = *__varcovbegin++;
 530                 for (size_t __k = 0; __k < __i; ++__k)
 531                   __s -= *__slit++ * *__cit++;
 532
 533                 *__w++ = __s /= *__cit++;
 534                 __sum += __s * __s;
 535               }
 536
 537             __sum = *__varcovbegin - __sum;
 538             if (__builtin_expect(__sum <= _RealType(0), 0))
 539               std::__throw_runtime_error(__N("normal_mv_distribution::"
 540                                              "param_type::_M_init_full"));
 541             *__w++ = std::sqrt(__sum);
 542
 543             std::advance(__varcovbegin, _Dimen - __j);
 544           }
 545       }
 546
 547   template<std::size_t _Dimen, typename _RealType>
 548     template<typename _InputIterator1, typename _InputIterator2>
 549       void
 550       normal_mv_distribution<_Dimen, _RealType>::param_type::
 551       _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
 552                     _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
 553       {
 554         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
 555         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
 556         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
 557                   _M_mean.end(), _RealType(0));
 558
 559         // Perform the Cholesky decomposition
 560         auto __w = _M_t.begin();
 561         for (size_t __j = 0; __j < _Dimen; ++__j)
 562           {
 563             _RealType __sum = _RealType(0);
 564
 565             auto __slitbegin = __w;
 566             auto __cit = _M_t.begin();
 567             for (size_t __i = 0; __i < __j; ++__i)
 568               {
 569                 auto __slit = __slitbegin;
 570                 _RealType __s = *__varcovbegin++;
 571                 for (size_t __k = 0; __k < __i; ++__k)
 572                   __s -= *__slit++ * *__cit++;
 573
 574                 *__w++ = __s /= *__cit++;
 575                 __sum += __s * __s;
 576               }
 577
 578             __sum = *__varcovbegin++ - __sum;
 579             if (__builtin_expect(__sum <= _RealType(0), 0))
 580               std::__throw_runtime_error(__N("normal_mv_distribution::"
 581                                              "param_type::_M_init_full"));
 582             *__w++ = std::sqrt(__sum);
 583           }
 584       }
 585
 586   template<std::size_t _Dimen, typename _RealType>
 587     template<typename _InputIterator1, typename _InputIterator2>
 588       void
 589       normal_mv_distribution<_Dimen, _RealType>::param_type::
 590       _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
 591                        _InputIterator2 __varbegin, _InputIterator2 __varend)
 592       {
 593         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
 594         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
 595         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
 596                   _M_mean.end(), _RealType(0));
 597
 598         auto __w = _M_t.begin();
 599         size_t __step = 0;
 600         while (__varbegin != __varend)
 601           {
 602             std::fill_n(__w, __step, _RealType(0));
 603             __w += __step++;
 604             if (__builtin_expect(*__varbegin < _RealType(0), 0))
 605               std::__throw_runtime_error(__N("normal_mv_distribution::"
 606                                              "param_type::_M_init_diagonal"));
 607             *__w++ = std::sqrt(*__varbegin++);
 608           }
 609       }
 610
 611   template<std::size_t _Dimen, typename _RealType>
 612     template<typename _UniformRandomNumberGenerator>
 613       typename normal_mv_distribution<_Dimen, _RealType>::result_type
 614       normal_mv_distribution<_Dimen, _RealType>::
 615       operator()(_UniformRandomNumberGenerator& __urng,
 616                  const param_type& __param)
 617       {
 618         result_type __ret;
 619
 620         _M_nd.__generate(__ret.begin(), __ret.end(), __urng);
 621
 622         auto __t_it = __param._M_t.crbegin();
 623         for (size_t __i = _Dimen; __i > 0; --__i)
 624           {
 625             _RealType __sum = _RealType(0);
 626             for (size_t __j = __i; __j > 0; --__j)
 627               __sum += __ret[__j - 1] * *__t_it++;
 628             __ret[__i - 1] = __sum;
 629           }
 630
 631         return __ret;
 632       }
 633
 634   template<std::size_t _Dimen, typename _RealType>
 635     template<typename _ForwardIterator, typename _UniformRandomNumberGenerator>
 636       void
 637       normal_mv_distribution<_Dimen, _RealType>::
 638       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
 639                       _UniformRandomNumberGenerator& __urng,
 640                       const param_type& __param)
 641       {
 642         __glibcxx_function_requires(_Mutable_ForwardIteratorConcept<
 643                                     _ForwardIterator>)
 644         while (__f != __t)
 645           *__f++ = this->operator()(__urng, __param);
 646       }
 647
 648   template<size_t _Dimen, typename _RealType>
 649     bool
 650     operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
 651                __d1,
 652                const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
 653                __d2)
 654     {
 655       return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd;
 656     }
 657
 658   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
 659     std::basic_ostream<_CharT, _Traits>&
 660     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 661                const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
 662     {
 663       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 664       typedef typename __ostream_type::ios_base    __ios_base;
 665
 666       const typename __ios_base::fmtflags __flags = __os.flags();
 667       const _CharT __fill = __os.fill();
 668       const std::streamsize __precision = __os.precision();
 669       const _CharT __space = __os.widen(' ');
 670       __os.flags(__ios_base::scientific | __ios_base::left);
 671       __os.fill(__space);
 672       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 673
 674       auto __mean = __x._M_param.mean();
 675       for (auto __it : __mean)
 676         __os << __it << __space;
 677       auto __t = __x._M_param.varcov();
 678       for (auto __it : __t)
 679         __os << __it << __space;
 680
 681       __os << __x._M_nd;
 682
 683       __os.flags(__flags);
 684       __os.fill(__fill);
 685       __os.precision(__precision);
 686       return __os;
 687     }
 688
 689   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
 690     std::basic_istream<_CharT, _Traits>&
 691     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 692                __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
 693     {
 694       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 695       typedef typename __istream_type::ios_base    __ios_base;
 696
 697       const typename __ios_base::fmtflags __flags = __is.flags();
 698       __is.flags(__ios_base::dec | __ios_base::skipws);
 699
 700       std::array<_RealType, _Dimen> __mean;
 701       for (auto& __it : __mean)
 702         __is >> __it;
 703       std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov;
 704       for (auto& __it : __varcov)
 705         __is >> __it;
 706
 707       __is >> __x._M_nd;
 708
 709       __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
 710                 param_type(__mean.begin(), __mean.end(),
 711                            __varcov.begin(), __varcov.end()));
 712
 713       __is.flags(__flags);
 714       return __is;
 715     }
 716
 717
 718   template<typename _RealType>
 719     template<typename _OutputIterator,
 720              typename _UniformRandomNumberGenerator>
 721       void
 722       rice_distribution<_RealType>::
 723       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 724                       _UniformRandomNumberGenerator& __urng,
 725                       const param_type& __p)
 726       {
 727         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 728
 729         while (__f != __t)
 730           {
 731             typename std::normal_distribution<result_type>::param_type
 732               __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
 733             result_type __x = this->_M_ndx(__px, __urng);
 734             result_type __y = this->_M_ndy(__py, __urng);
 735             *__f++ = std::hypot(__x, __y);
 736           }
 737       }
 738
 739   template<typename _RealType, typename _CharT, typename _Traits>
 740     std::basic_ostream<_CharT, _Traits>&
 741     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 742                const rice_distribution<_RealType>& __x)
 743     {
 744       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 745       typedef typename __ostream_type::ios_base    __ios_base;
 746
 747       const typename __ios_base::fmtflags __flags = __os.flags();
 748       const _CharT __fill = __os.fill();
 749       const std::streamsize __precision = __os.precision();
 750       const _CharT __space = __os.widen(' ');
 751       __os.flags(__ios_base::scientific | __ios_base::left);
 752       __os.fill(__space);
 753       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 754
 755       __os << __x.nu() << __space << __x.sigma();
 756       __os << __space << __x._M_ndx;
 757       __os << __space << __x._M_ndy;
 758
 759       __os.flags(__flags);
 760       __os.fill(__fill);
 761       __os.precision(__precision);
 762       return __os;
 763     }
 764
 765   template<typename _RealType, typename _CharT, typename _Traits>
 766     std::basic_istream<_CharT, _Traits>&
 767     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 768                rice_distribution<_RealType>& __x)
 769     {
 770       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 771       typedef typename __istream_type::ios_base    __ios_base;
 772
 773       const typename __ios_base::fmtflags __flags = __is.flags();
 774       __is.flags(__ios_base::dec | __ios_base::skipws);
 775
 776       _RealType __nu, __sigma;
 777       __is >> __nu >> __sigma;
 778       __is >> __x._M_ndx;
 779       __is >> __x._M_ndy;
 780       __x.param(typename rice_distribution<_RealType>::
 781                 param_type(__nu, __sigma));
 782
 783       __is.flags(__flags);
 784       return __is;
 785     }
 786
 787
 788   template<typename _RealType>
 789     template<typename _OutputIterator,
 790              typename _UniformRandomNumberGenerator>
 791       void
 792       nakagami_distribution<_RealType>::
 793       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 794                       _UniformRandomNumberGenerator& __urng,
 795                       const param_type& __p)
 796       {
 797         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 798
 799         typename std::gamma_distribution<result_type>::param_type
 800           __pg(__p.mu(), __p.omega() / __p.mu());
 801         while (__f != __t)
 802           *__f++ = std::sqrt(this->_M_gd(__pg, __urng));
 803       }
 804
 805   template<typename _RealType, typename _CharT, typename _Traits>
 806     std::basic_ostream<_CharT, _Traits>&
 807     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 808                const nakagami_distribution<_RealType>& __x)
 809     {
 810       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 811       typedef typename __ostream_type::ios_base    __ios_base;
 812
 813       const typename __ios_base::fmtflags __flags = __os.flags();
 814       const _CharT __fill = __os.fill();
 815       const std::streamsize __precision = __os.precision();
 816       const _CharT __space = __os.widen(' ');
 817       __os.flags(__ios_base::scientific | __ios_base::left);
 818       __os.fill(__space);
 819       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 820
 821       __os << __x.mu() << __space << __x.omega();
 822       __os << __space << __x._M_gd;
 823
 824       __os.flags(__flags);
 825       __os.fill(__fill);
 826       __os.precision(__precision);
 827       return __os;
 828     }
 829
 830   template<typename _RealType, typename _CharT, typename _Traits>
 831     std::basic_istream<_CharT, _Traits>&
 832     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 833                nakagami_distribution<_RealType>& __x)
 834     {
 835       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 836       typedef typename __istream_type::ios_base    __ios_base;
 837
 838       const typename __ios_base::fmtflags __flags = __is.flags();
 839       __is.flags(__ios_base::dec | __ios_base::skipws);
 840
 841       _RealType __mu, __omega;
 842       __is >> __mu >> __omega;
 843       __is >> __x._M_gd;
 844       __x.param(typename nakagami_distribution<_RealType>::
 845                 param_type(__mu, __omega));
 846
 847       __is.flags(__flags);
 848       return __is;
 849     }
 850
 851 _GLIBCXX_END_NAMESPACE_VERSION
 852 } // namespace
 853
 854
 855 #endif // _EXT_RANDOM_TCC