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

   1 // Random number extensions -*- C++ -*-
   2
   3 // Copyright (C) 2012-2014 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       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 #endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
 406
 407   /**
 408    * Iteration method due to M.D. J<o:>hnk.
 409    *
 410    * M.D. J<o:>hnk, Erzeugung von betaverteilten und gammaverteilten
 411    * Zufallszahlen, Metrika, Volume 8, 1964
 412    */
 413   template<typename _RealType>
 414     template<typename _UniformRandomNumberGenerator>
 415       typename beta_distribution<_RealType>::result_type
 416       beta_distribution<_RealType>::
 417       operator()(_UniformRandomNumberGenerator& __urng,
 418                  const param_type& __param)
 419       {
 420         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
 421           __aurng(__urng);
 422
 423         result_type __x, __y;
 424         do
 425           {
 426             __x = std::exp(std::log(__aurng()) / __param.alpha());
 427             __y = std::exp(std::log(__aurng()) / __param.beta());
 428           }
 429         while (__x + __y > result_type(1));
 430
 431         return __x / (__x + __y);
 432       }
 433
 434   template<typename _RealType>
 435     template<typename _OutputIterator,
 436              typename _UniformRandomNumberGenerator>
 437       void
 438       beta_distribution<_RealType>::
 439       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 440                       _UniformRandomNumberGenerator& __urng,
 441                       const param_type& __param)
 442       {
 443         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 444
 445         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
 446           __aurng(__urng);
 447
 448         while (__f != __t)
 449           {
 450             result_type __x, __y;
 451             do
 452               {
 453                 __x = std::exp(std::log(__aurng()) / __param.alpha());
 454                 __y = std::exp(std::log(__aurng()) / __param.beta());
 455               }
 456             while (__x + __y > result_type(1));
 457
 458             *__f++ = __x / (__x + __y);
 459           }
 460       }
 461
 462   template<typename _RealType, typename _CharT, typename _Traits>
 463     std::basic_ostream<_CharT, _Traits>&
 464     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 465                const __gnu_cxx::beta_distribution<_RealType>& __x)
 466     {
 467       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 468       typedef typename __ostream_type::ios_base    __ios_base;
 469
 470       const typename __ios_base::fmtflags __flags = __os.flags();
 471       const _CharT __fill = __os.fill();
 472       const std::streamsize __precision = __os.precision();
 473       const _CharT __space = __os.widen(' ');
 474       __os.flags(__ios_base::scientific | __ios_base::left);
 475       __os.fill(__space);
 476       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 477
 478       __os << __x.alpha() << __space << __x.beta();
 479
 480       __os.flags(__flags);
 481       __os.fill(__fill);
 482       __os.precision(__precision);
 483       return __os;
 484     }
 485
 486   template<typename _RealType, typename _CharT, typename _Traits>
 487     std::basic_istream<_CharT, _Traits>&
 488     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 489                __gnu_cxx::beta_distribution<_RealType>& __x)
 490     {
 491       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 492       typedef typename __istream_type::ios_base    __ios_base;
 493
 494       const typename __ios_base::fmtflags __flags = __is.flags();
 495       __is.flags(__ios_base::dec | __ios_base::skipws);
 496
 497       _RealType __alpha_val, __beta_val;
 498       __is >> __alpha_val >> __beta_val;
 499       __x.param(typename __gnu_cxx::beta_distribution<_RealType>::
 500                 param_type(__alpha_val, __beta_val));
 501
 502       __is.flags(__flags);
 503       return __is;
 504     }
 505
 506
 507   template<std::size_t _Dimen, typename _RealType>
 508     template<typename _InputIterator1, typename _InputIterator2>
 509       void
 510       normal_mv_distribution<_Dimen, _RealType>::param_type::
 511       _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
 512                    _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
 513       {
 514         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
 515         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
 516         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
 517                   _M_mean.end(), _RealType(0));
 518
 519         // Perform the Cholesky decomposition
 520         auto __w = _M_t.begin();
 521         for (size_t __j = 0; __j < _Dimen; ++__j)
 522           {
 523             _RealType __sum = _RealType(0);
 524
 525             auto __slitbegin = __w;
 526             auto __cit = _M_t.begin();
 527             for (size_t __i = 0; __i < __j; ++__i)
 528               {
 529                 auto __slit = __slitbegin;
 530                 _RealType __s = *__varcovbegin++;
 531                 for (size_t __k = 0; __k < __i; ++__k)
 532                   __s -= *__slit++ * *__cit++;
 533
 534                 *__w++ = __s /= *__cit++;
 535                 __sum += __s * __s;
 536               }
 537
 538             __sum = *__varcovbegin - __sum;
 539             if (__builtin_expect(__sum <= _RealType(0), 0))
 540               std::__throw_runtime_error(__N("normal_mv_distribution::"
 541                                              "param_type::_M_init_full"));
 542             *__w++ = std::sqrt(__sum);
 543
 544             std::advance(__varcovbegin, _Dimen - __j);
 545           }
 546       }
 547
 548   template<std::size_t _Dimen, typename _RealType>
 549     template<typename _InputIterator1, typename _InputIterator2>
 550       void
 551       normal_mv_distribution<_Dimen, _RealType>::param_type::
 552       _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
 553                     _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
 554       {
 555         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
 556         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
 557         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
 558                   _M_mean.end(), _RealType(0));
 559
 560         // Perform the Cholesky decomposition
 561         auto __w = _M_t.begin();
 562         for (size_t __j = 0; __j < _Dimen; ++__j)
 563           {
 564             _RealType __sum = _RealType(0);
 565
 566             auto __slitbegin = __w;
 567             auto __cit = _M_t.begin();
 568             for (size_t __i = 0; __i < __j; ++__i)
 569               {
 570                 auto __slit = __slitbegin;
 571                 _RealType __s = *__varcovbegin++;
 572                 for (size_t __k = 0; __k < __i; ++__k)
 573                   __s -= *__slit++ * *__cit++;
 574
 575                 *__w++ = __s /= *__cit++;
 576                 __sum += __s * __s;
 577               }
 578
 579             __sum = *__varcovbegin++ - __sum;
 580             if (__builtin_expect(__sum <= _RealType(0), 0))
 581               std::__throw_runtime_error(__N("normal_mv_distribution::"
 582                                              "param_type::_M_init_full"));
 583             *__w++ = std::sqrt(__sum);
 584           }
 585       }
 586
 587   template<std::size_t _Dimen, typename _RealType>
 588     template<typename _InputIterator1, typename _InputIterator2>
 589       void
 590       normal_mv_distribution<_Dimen, _RealType>::param_type::
 591       _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
 592                        _InputIterator2 __varbegin, _InputIterator2 __varend)
 593       {
 594         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
 595         __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
 596         std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
 597                   _M_mean.end(), _RealType(0));
 598
 599         auto __w = _M_t.begin();
 600         size_t __step = 0;
 601         while (__varbegin != __varend)
 602           {
 603             std::fill_n(__w, __step, _RealType(0));
 604             __w += __step++;
 605             if (__builtin_expect(*__varbegin < _RealType(0), 0))
 606               std::__throw_runtime_error(__N("normal_mv_distribution::"
 607                                              "param_type::_M_init_diagonal"));
 608             *__w++ = std::sqrt(*__varbegin++);
 609           }
 610       }
 611
 612   template<std::size_t _Dimen, typename _RealType>
 613     template<typename _UniformRandomNumberGenerator>
 614       typename normal_mv_distribution<_Dimen, _RealType>::result_type
 615       normal_mv_distribution<_Dimen, _RealType>::
 616       operator()(_UniformRandomNumberGenerator& __urng,
 617                  const param_type& __param)
 618       {
 619         result_type __ret;
 620
 621         _M_nd.__generate(__ret.begin(), __ret.end(), __urng);
 622
 623         auto __t_it = __param._M_t.crbegin();
 624         for (size_t __i = _Dimen; __i > 0; --__i)
 625           {
 626             _RealType __sum = _RealType(0);
 627             for (size_t __j = __i; __j > 0; --__j)
 628               __sum += __ret[__j - 1] * *__t_it++;
 629             __ret[__i - 1] = __sum;
 630           }
 631
 632         return __ret;
 633       }
 634
 635   template<std::size_t _Dimen, typename _RealType>
 636     template<typename _ForwardIterator, typename _UniformRandomNumberGenerator>
 637       void
 638       normal_mv_distribution<_Dimen, _RealType>::
 639       __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
 640                       _UniformRandomNumberGenerator& __urng,
 641                       const param_type& __param)
 642       {
 643         __glibcxx_function_requires(_Mutable_ForwardIteratorConcept<
 644                                     _ForwardIterator>)
 645         while (__f != __t)
 646           *__f++ = this->operator()(__urng, __param);
 647       }
 648
 649   template<size_t _Dimen, typename _RealType>
 650     bool
 651     operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
 652                __d1,
 653                const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
 654                __d2)
 655     {
 656       return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd;
 657     }
 658
 659   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
 660     std::basic_ostream<_CharT, _Traits>&
 661     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 662                const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
 663     {
 664       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 665       typedef typename __ostream_type::ios_base    __ios_base;
 666
 667       const typename __ios_base::fmtflags __flags = __os.flags();
 668       const _CharT __fill = __os.fill();
 669       const std::streamsize __precision = __os.precision();
 670       const _CharT __space = __os.widen(' ');
 671       __os.flags(__ios_base::scientific | __ios_base::left);
 672       __os.fill(__space);
 673       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 674
 675       auto __mean = __x._M_param.mean();
 676       for (auto __it : __mean)
 677         __os << __it << __space;
 678       auto __t = __x._M_param.varcov();
 679       for (auto __it : __t)
 680         __os << __it << __space;
 681
 682       __os << __x._M_nd;
 683
 684       __os.flags(__flags);
 685       __os.fill(__fill);
 686       __os.precision(__precision);
 687       return __os;
 688     }
 689
 690   template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
 691     std::basic_istream<_CharT, _Traits>&
 692     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 693                __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
 694     {
 695       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 696       typedef typename __istream_type::ios_base    __ios_base;
 697
 698       const typename __ios_base::fmtflags __flags = __is.flags();
 699       __is.flags(__ios_base::dec | __ios_base::skipws);
 700
 701       std::array<_RealType, _Dimen> __mean;
 702       for (auto& __it : __mean)
 703         __is >> __it;
 704       std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov;
 705       for (auto& __it : __varcov)
 706         __is >> __it;
 707
 708       __is >> __x._M_nd;
 709
 710       __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
 711                 param_type(__mean.begin(), __mean.end(),
 712                            __varcov.begin(), __varcov.end()));
 713
 714       __is.flags(__flags);
 715       return __is;
 716     }
 717
 718
 719   template<typename _RealType>
 720     template<typename _OutputIterator,
 721              typename _UniformRandomNumberGenerator>
 722       void
 723       rice_distribution<_RealType>::
 724       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 725                       _UniformRandomNumberGenerator& __urng,
 726                       const param_type& __p)
 727       {
 728         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 729
 730         while (__f != __t)
 731           {
 732             typename std::normal_distribution<result_type>::param_type
 733               __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
 734             result_type __x = this->_M_ndx(__px, __urng);
 735             result_type __y = this->_M_ndy(__py, __urng);
 736 #if _GLIBCXX_USE_C99_MATH_TR1
 737             *__f++ = std::hypot(__x, __y);
 738 #else
 739             *__f++ = std::sqrt(__x * __x + __y * __y);
 740 #endif
 741           }
 742       }
 743
 744   template<typename _RealType, typename _CharT, typename _Traits>
 745     std::basic_ostream<_CharT, _Traits>&
 746     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 747                const rice_distribution<_RealType>& __x)
 748     {
 749       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 750       typedef typename __ostream_type::ios_base    __ios_base;
 751
 752       const typename __ios_base::fmtflags __flags = __os.flags();
 753       const _CharT __fill = __os.fill();
 754       const std::streamsize __precision = __os.precision();
 755       const _CharT __space = __os.widen(' ');
 756       __os.flags(__ios_base::scientific | __ios_base::left);
 757       __os.fill(__space);
 758       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 759
 760       __os << __x.nu() << __space << __x.sigma();
 761       __os << __space << __x._M_ndx;
 762       __os << __space << __x._M_ndy;
 763
 764       __os.flags(__flags);
 765       __os.fill(__fill);
 766       __os.precision(__precision);
 767       return __os;
 768     }
 769
 770   template<typename _RealType, typename _CharT, typename _Traits>
 771     std::basic_istream<_CharT, _Traits>&
 772     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 773                rice_distribution<_RealType>& __x)
 774     {
 775       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 776       typedef typename __istream_type::ios_base    __ios_base;
 777
 778       const typename __ios_base::fmtflags __flags = __is.flags();
 779       __is.flags(__ios_base::dec | __ios_base::skipws);
 780
 781       _RealType __nu_val, __sigma_val;
 782       __is >> __nu_val >> __sigma_val;
 783       __is >> __x._M_ndx;
 784       __is >> __x._M_ndy;
 785       __x.param(typename rice_distribution<_RealType>::
 786                 param_type(__nu_val, __sigma_val));
 787
 788       __is.flags(__flags);
 789       return __is;
 790     }
 791
 792
 793   template<typename _RealType>
 794     template<typename _OutputIterator,
 795              typename _UniformRandomNumberGenerator>
 796       void
 797       nakagami_distribution<_RealType>::
 798       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 799                       _UniformRandomNumberGenerator& __urng,
 800                       const param_type& __p)
 801       {
 802         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 803
 804         typename std::gamma_distribution<result_type>::param_type
 805           __pg(__p.mu(), __p.omega() / __p.mu());
 806         while (__f != __t)
 807           *__f++ = std::sqrt(this->_M_gd(__pg, __urng));
 808       }
 809
 810   template<typename _RealType, typename _CharT, typename _Traits>
 811     std::basic_ostream<_CharT, _Traits>&
 812     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 813                const nakagami_distribution<_RealType>& __x)
 814     {
 815       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 816       typedef typename __ostream_type::ios_base    __ios_base;
 817
 818       const typename __ios_base::fmtflags __flags = __os.flags();
 819       const _CharT __fill = __os.fill();
 820       const std::streamsize __precision = __os.precision();
 821       const _CharT __space = __os.widen(' ');
 822       __os.flags(__ios_base::scientific | __ios_base::left);
 823       __os.fill(__space);
 824       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 825
 826       __os << __x.mu() << __space << __x.omega();
 827       __os << __space << __x._M_gd;
 828
 829       __os.flags(__flags);
 830       __os.fill(__fill);
 831       __os.precision(__precision);
 832       return __os;
 833     }
 834
 835   template<typename _RealType, typename _CharT, typename _Traits>
 836     std::basic_istream<_CharT, _Traits>&
 837     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 838                nakagami_distribution<_RealType>& __x)
 839     {
 840       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 841       typedef typename __istream_type::ios_base    __ios_base;
 842
 843       const typename __ios_base::fmtflags __flags = __is.flags();
 844       __is.flags(__ios_base::dec | __ios_base::skipws);
 845
 846       _RealType __mu_val, __omega_val;
 847       __is >> __mu_val >> __omega_val;
 848       __is >> __x._M_gd;
 849       __x.param(typename nakagami_distribution<_RealType>::
 850                 param_type(__mu_val, __omega_val));
 851
 852       __is.flags(__flags);
 853       return __is;
 854     }
 855
 856
 857   template<typename _RealType>
 858     template<typename _OutputIterator,
 859              typename _UniformRandomNumberGenerator>
 860       void
 861       pareto_distribution<_RealType>::
 862       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 863                       _UniformRandomNumberGenerator& __urng,
 864                       const param_type& __p)
 865       {
 866         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 867
 868         result_type __mu_val = __p.mu();
 869         result_type __malphinv = -result_type(1) / __p.alpha();
 870         while (__f != __t)
 871           *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv);
 872       }
 873
 874   template<typename _RealType, typename _CharT, typename _Traits>
 875     std::basic_ostream<_CharT, _Traits>&
 876     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 877                const pareto_distribution<_RealType>& __x)
 878     {
 879       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 880       typedef typename __ostream_type::ios_base    __ios_base;
 881
 882       const typename __ios_base::fmtflags __flags = __os.flags();
 883       const _CharT __fill = __os.fill();
 884       const std::streamsize __precision = __os.precision();
 885       const _CharT __space = __os.widen(' ');
 886       __os.flags(__ios_base::scientific | __ios_base::left);
 887       __os.fill(__space);
 888       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 889
 890       __os << __x.alpha() << __space << __x.mu();
 891       __os << __space << __x._M_ud;
 892
 893       __os.flags(__flags);
 894       __os.fill(__fill);
 895       __os.precision(__precision);
 896       return __os;
 897     }
 898
 899   template<typename _RealType, typename _CharT, typename _Traits>
 900     std::basic_istream<_CharT, _Traits>&
 901     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 902                pareto_distribution<_RealType>& __x)
 903     {
 904       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
 905       typedef typename __istream_type::ios_base    __ios_base;
 906
 907       const typename __ios_base::fmtflags __flags = __is.flags();
 908       __is.flags(__ios_base::dec | __ios_base::skipws);
 909
 910       _RealType __alpha_val, __mu_val;
 911       __is >> __alpha_val >> __mu_val;
 912       __is >> __x._M_ud;
 913       __x.param(typename pareto_distribution<_RealType>::
 914                 param_type(__alpha_val, __mu_val));
 915
 916       __is.flags(__flags);
 917       return __is;
 918     }
 919
 920
 921   template<typename _RealType>
 922     template<typename _UniformRandomNumberGenerator>
 923       typename k_distribution<_RealType>::result_type
 924       k_distribution<_RealType>::
 925       operator()(_UniformRandomNumberGenerator& __urng)
 926       {
 927         result_type __x = this->_M_gd1(__urng);
 928         result_type __y = this->_M_gd2(__urng);
 929         return std::sqrt(__x * __y);
 930       }
 931
 932   template<typename _RealType>
 933     template<typename _UniformRandomNumberGenerator>
 934       typename k_distribution<_RealType>::result_type
 935       k_distribution<_RealType>::
 936       operator()(_UniformRandomNumberGenerator& __urng,
 937                  const param_type& __p)
 938       {
 939         typename std::gamma_distribution<result_type>::param_type
 940           __p1(__p.lambda(), result_type(1) / __p.lambda()),
 941           __p2(__p.nu(), __p.mu() / __p.nu());
 942         result_type __x = this->_M_gd1(__p1, __urng);
 943         result_type __y = this->_M_gd2(__p2, __urng);
 944         return std::sqrt(__x * __y);
 945       }
 946
 947   template<typename _RealType>
 948     template<typename _OutputIterator,
 949              typename _UniformRandomNumberGenerator>
 950       void
 951       k_distribution<_RealType>::
 952       __generate_impl(_OutputIterator __f, _OutputIterator __t,
 953                       _UniformRandomNumberGenerator& __urng,
 954                       const param_type& __p)
 955       {
 956         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
 957
 958         typename std::gamma_distribution<result_type>::param_type
 959           __p1(__p.lambda(), result_type(1) / __p.lambda()),
 960           __p2(__p.nu(), __p.mu() / __p.nu());
 961         while (__f != __t)
 962           {
 963             result_type __x = this->_M_gd1(__p1, __urng);
 964             result_type __y = this->_M_gd2(__p2, __urng);
 965             *__f++ = std::sqrt(__x * __y);
 966           }
 967       }
 968
 969   template<typename _RealType, typename _CharT, typename _Traits>
 970     std::basic_ostream<_CharT, _Traits>&
 971     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
 972                const k_distribution<_RealType>& __x)
 973     {
 974       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
 975       typedef typename __ostream_type::ios_base    __ios_base;
 976
 977       const typename __ios_base::fmtflags __flags = __os.flags();
 978       const _CharT __fill = __os.fill();
 979       const std::streamsize __precision = __os.precision();
 980       const _CharT __space = __os.widen(' ');
 981       __os.flags(__ios_base::scientific | __ios_base::left);
 982       __os.fill(__space);
 983       __os.precision(std::numeric_limits<_RealType>::max_digits10);
 984
 985       __os << __x.lambda() << __space << __x.mu() << __space << __x.nu();
 986       __os << __space << __x._M_gd1;
 987       __os << __space << __x._M_gd2;
 988
 989       __os.flags(__flags);
 990       __os.fill(__fill);
 991       __os.precision(__precision);
 992       return __os;
 993     }
 994
 995   template<typename _RealType, typename _CharT, typename _Traits>
 996     std::basic_istream<_CharT, _Traits>&
 997     operator>>(std::basic_istream<_CharT, _Traits>& __is,
 998                k_distribution<_RealType>& __x)
 999     {
1000       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1001       typedef typename __istream_type::ios_base    __ios_base;
1002
1003       const typename __ios_base::fmtflags __flags = __is.flags();
1004       __is.flags(__ios_base::dec | __ios_base::skipws);
1005
1006       _RealType __lambda_val, __mu_val, __nu_val;
1007       __is >> __lambda_val >> __mu_val >> __nu_val;
1008       __is >> __x._M_gd1;
1009       __is >> __x._M_gd2;
1010       __x.param(typename k_distribution<_RealType>::
1011                 param_type(__lambda_val, __mu_val, __nu_val));
1012
1013       __is.flags(__flags);
1014       return __is;
1015     }
1016
1017
1018   template<typename _RealType>
1019     template<typename _OutputIterator,
1020              typename _UniformRandomNumberGenerator>
1021       void
1022       arcsine_distribution<_RealType>::
1023       __generate_impl(_OutputIterator __f, _OutputIterator __t,
1024                       _UniformRandomNumberGenerator& __urng,
1025                       const param_type& __p)
1026       {
1027         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
1028
1029         result_type __dif = __p.b() - __p.a();
1030         result_type __sum = __p.a() + __p.b();
1031         while (__f != __t)
1032           {
1033             result_type __x = std::sin(this->_M_ud(__urng));
1034             *__f++ = (__x * __dif + __sum) / result_type(2);
1035           }
1036       }
1037
1038   template<typename _RealType, typename _CharT, typename _Traits>
1039     std::basic_ostream<_CharT, _Traits>&
1040     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1041                const arcsine_distribution<_RealType>& __x)
1042     {
1043       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1044       typedef typename __ostream_type::ios_base    __ios_base;
1045
1046       const typename __ios_base::fmtflags __flags = __os.flags();
1047       const _CharT __fill = __os.fill();
1048       const std::streamsize __precision = __os.precision();
1049       const _CharT __space = __os.widen(' ');
1050       __os.flags(__ios_base::scientific | __ios_base::left);
1051       __os.fill(__space);
1052       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1053
1054       __os << __x.a() << __space << __x.b();
1055       __os << __space << __x._M_ud;
1056
1057       __os.flags(__flags);
1058       __os.fill(__fill);
1059       __os.precision(__precision);
1060       return __os;
1061     }
1062
1063   template<typename _RealType, typename _CharT, typename _Traits>
1064     std::basic_istream<_CharT, _Traits>&
1065     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1066                arcsine_distribution<_RealType>& __x)
1067     {
1068       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1069       typedef typename __istream_type::ios_base    __ios_base;
1070
1071       const typename __ios_base::fmtflags __flags = __is.flags();
1072       __is.flags(__ios_base::dec | __ios_base::skipws);
1073
1074       _RealType __a, __b;
1075       __is >> __a >> __b;
1076       __is >> __x._M_ud;
1077       __x.param(typename arcsine_distribution<_RealType>::
1078                 param_type(__a, __b));
1079
1080       __is.flags(__flags);
1081       return __is;
1082     }
1083
1084
1085   template<typename _RealType>
1086     template<typename _UniformRandomNumberGenerator>
1087       typename hoyt_distribution<_RealType>::result_type
1088       hoyt_distribution<_RealType>::
1089       operator()(_UniformRandomNumberGenerator& __urng)
1090       {
1091         result_type __x = this->_M_ad(__urng);
1092         result_type __y = this->_M_ed(__urng);
1093         return (result_type(2) * this->q()
1094                   / (result_type(1) + this->q() * this->q()))
1095                * std::sqrt(this->omega() * __x * __y);
1096       }
1097
1098   template<typename _RealType>
1099     template<typename _UniformRandomNumberGenerator>
1100       typename hoyt_distribution<_RealType>::result_type
1101       hoyt_distribution<_RealType>::
1102       operator()(_UniformRandomNumberGenerator& __urng,
1103                  const param_type& __p)
1104       {
1105         result_type __q2 = __p.q() * __p.q();
1106         result_type __num = result_type(0.5L) * (result_type(1) + __q2);
1107         typename __gnu_cxx::arcsine_distribution<result_type>::param_type
1108           __pa(__num, __num / __q2);
1109         result_type __x = this->_M_ad(__pa, __urng);
1110         result_type __y = this->_M_ed(__urng);
1111         return (result_type(2) * __p.q() / (result_type(1) + __q2))
1112                * std::sqrt(__p.omega() * __x * __y);
1113       }
1114
1115   template<typename _RealType>
1116     template<typename _OutputIterator,
1117              typename _UniformRandomNumberGenerator>
1118       void
1119       hoyt_distribution<_RealType>::
1120       __generate_impl(_OutputIterator __f, _OutputIterator __t,
1121                       _UniformRandomNumberGenerator& __urng,
1122                       const param_type& __p)
1123       {
1124         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
1125
1126         result_type __2q = result_type(2) * __p.q();
1127         result_type __q2 = __p.q() * __p.q();
1128         result_type __q2p1 = result_type(1) + __q2;
1129         result_type __num = result_type(0.5L) * __q2p1;
1130         result_type __omega = __p.omega();
1131         typename __gnu_cxx::arcsine_distribution<result_type>::param_type
1132           __pa(__num, __num / __q2);
1133         while (__f != __t)
1134           {
1135             result_type __x = this->_M_ad(__pa, __urng);
1136             result_type __y = this->_M_ed(__urng);
1137             *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y);
1138           }
1139       }
1140
1141   template<typename _RealType, typename _CharT, typename _Traits>
1142     std::basic_ostream<_CharT, _Traits>&
1143     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1144                const hoyt_distribution<_RealType>& __x)
1145     {
1146       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1147       typedef typename __ostream_type::ios_base    __ios_base;
1148
1149       const typename __ios_base::fmtflags __flags = __os.flags();
1150       const _CharT __fill = __os.fill();
1151       const std::streamsize __precision = __os.precision();
1152       const _CharT __space = __os.widen(' ');
1153       __os.flags(__ios_base::scientific | __ios_base::left);
1154       __os.fill(__space);
1155       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1156
1157       __os << __x.q() << __space << __x.omega();
1158       __os << __space << __x._M_ad;
1159       __os << __space << __x._M_ed;
1160
1161       __os.flags(__flags);
1162       __os.fill(__fill);
1163       __os.precision(__precision);
1164       return __os;
1165     }
1166
1167   template<typename _RealType, typename _CharT, typename _Traits>
1168     std::basic_istream<_CharT, _Traits>&
1169     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1170                hoyt_distribution<_RealType>& __x)
1171     {
1172       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1173       typedef typename __istream_type::ios_base    __ios_base;
1174
1175       const typename __ios_base::fmtflags __flags = __is.flags();
1176       __is.flags(__ios_base::dec | __ios_base::skipws);
1177
1178       _RealType __q, __omega;
1179       __is >> __q >> __omega;
1180       __is >> __x._M_ad;
1181       __is >> __x._M_ed;
1182       __x.param(typename hoyt_distribution<_RealType>::
1183                 param_type(__q, __omega));
1184
1185       __is.flags(__flags);
1186       return __is;
1187     }
1188
1189
1190   template<typename _RealType>
1191     template<typename _OutputIterator,
1192              typename _UniformRandomNumberGenerator>
1193       void
1194       triangular_distribution<_RealType>::
1195       __generate_impl(_OutputIterator __f, _OutputIterator __t,
1196                       _UniformRandomNumberGenerator& __urng,
1197                       const param_type& __param)
1198       {
1199         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
1200
1201         while (__f != __t)
1202           *__f++ = this->operator()(__urng, __param);
1203       }
1204
1205   template<typename _RealType, typename _CharT, typename _Traits>
1206     std::basic_ostream<_CharT, _Traits>&
1207     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1208                const __gnu_cxx::triangular_distribution<_RealType>& __x)
1209     {
1210       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1211       typedef typename __ostream_type::ios_base    __ios_base;
1212
1213       const typename __ios_base::fmtflags __flags = __os.flags();
1214       const _CharT __fill = __os.fill();
1215       const std::streamsize __precision = __os.precision();
1216       const _CharT __space = __os.widen(' ');
1217       __os.flags(__ios_base::scientific | __ios_base::left);
1218       __os.fill(__space);
1219       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1220
1221       __os << __x.a() << __space << __x.b() << __space << __x.c();
1222
1223       __os.flags(__flags);
1224       __os.fill(__fill);
1225       __os.precision(__precision);
1226       return __os;
1227     }
1228
1229   template<typename _RealType, typename _CharT, typename _Traits>
1230     std::basic_istream<_CharT, _Traits>&
1231     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1232                __gnu_cxx::triangular_distribution<_RealType>& __x)
1233     {
1234       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1235       typedef typename __istream_type::ios_base    __ios_base;
1236
1237       const typename __ios_base::fmtflags __flags = __is.flags();
1238       __is.flags(__ios_base::dec | __ios_base::skipws);
1239
1240       _RealType __a, __b, __c;
1241       __is >> __a >> __b >> __c;
1242       __x.param(typename __gnu_cxx::triangular_distribution<_RealType>::
1243                 param_type(__a, __b, __c));
1244
1245       __is.flags(__flags);
1246       return __is;
1247     }
1248
1249
1250   template<typename _RealType>
1251     template<typename _UniformRandomNumberGenerator>
1252       typename von_mises_distribution<_RealType>::result_type
1253       von_mises_distribution<_RealType>::
1254       operator()(_UniformRandomNumberGenerator& __urng,
1255                  const param_type& __p)
1256       {
1257         const result_type __pi
1258           = __gnu_cxx::__math_constants<result_type>::__pi;
1259         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1260           __aurng(__urng);
1261
1262         result_type __f;
1263         while (1)
1264           {
1265             result_type __rnd = std::cos(__pi * __aurng());
1266             __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd);
1267             result_type __c = __p._M_kappa * (__p._M_r - __f);
1268
1269             result_type __rnd2 = __aurng();
1270             if (__c * (result_type(2) - __c) > __rnd2)
1271               break;
1272             if (std::log(__c / __rnd2) >= __c - result_type(1))
1273               break;
1274           }
1275
1276         result_type __res = std::acos(__f);
1277 #if _GLIBCXX_USE_C99_MATH_TR1
1278         __res = std::copysign(__res, __aurng() - result_type(0.5));
1279 #else
1280         if (__aurng() < result_type(0.5))
1281           __res = -__res;
1282 #endif
1283         __res += __p._M_mu;
1284         if (__res > __pi)
1285           __res -= result_type(2) * __pi;
1286         else if (__res < -__pi)
1287           __res += result_type(2) * __pi;
1288         return __res;
1289       }
1290
1291   template<typename _RealType>
1292     template<typename _OutputIterator,
1293              typename _UniformRandomNumberGenerator>
1294       void
1295       von_mises_distribution<_RealType>::
1296       __generate_impl(_OutputIterator __f, _OutputIterator __t,
1297                       _UniformRandomNumberGenerator& __urng,
1298                       const param_type& __param)
1299       {
1300         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
1301
1302         while (__f != __t)
1303           *__f++ = this->operator()(__urng, __param);
1304       }
1305
1306   template<typename _RealType, typename _CharT, typename _Traits>
1307     std::basic_ostream<_CharT, _Traits>&
1308     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1309                const __gnu_cxx::von_mises_distribution<_RealType>& __x)
1310     {
1311       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1312       typedef typename __ostream_type::ios_base    __ios_base;
1313
1314       const typename __ios_base::fmtflags __flags = __os.flags();
1315       const _CharT __fill = __os.fill();
1316       const std::streamsize __precision = __os.precision();
1317       const _CharT __space = __os.widen(' ');
1318       __os.flags(__ios_base::scientific | __ios_base::left);
1319       __os.fill(__space);
1320       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1321
1322       __os << __x.mu() << __space << __x.kappa();
1323
1324       __os.flags(__flags);
1325       __os.fill(__fill);
1326       __os.precision(__precision);
1327       return __os;
1328     }
1329
1330   template<typename _RealType, typename _CharT, typename _Traits>
1331     std::basic_istream<_CharT, _Traits>&
1332     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1333                __gnu_cxx::von_mises_distribution<_RealType>& __x)
1334     {
1335       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1336       typedef typename __istream_type::ios_base    __ios_base;
1337
1338       const typename __ios_base::fmtflags __flags = __is.flags();
1339       __is.flags(__ios_base::dec | __ios_base::skipws);
1340
1341       _RealType __mu, __kappa;
1342       __is >> __mu >> __kappa;
1343       __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>::
1344                 param_type(__mu, __kappa));
1345
1346       __is.flags(__flags);
1347       return __is;
1348     }
1349
1350
1351   template<typename _UIntType>
1352     template<typename _UniformRandomNumberGenerator>
1353       typename hypergeometric_distribution<_UIntType>::result_type
1354       hypergeometric_distribution<_UIntType>::
1355       operator()(_UniformRandomNumberGenerator& __urng,
1356                  const param_type& __param)
1357       {
1358         std::__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1359           __aurng(__urng);
1360
1361         result_type __a = __param.successful_size();
1362         result_type __b = __param.total_size();
1363         result_type __k = 0;
1364
1365         if (__param.total_draws() < __param.total_size() / 2)
1366           {
1367             for (result_type __i = 0; __i < __param.total_draws(); ++__i)
1368               {
1369                 if (__b * __aurng() < __a)
1370                   {
1371                     ++__k;
1372                     if (__k == __param.successful_size())
1373                       return __k;
1374                    --__a;
1375                   }
1376                 --__b;
1377               }
1378             return __k;
1379           }
1380         else
1381           {
1382             for (result_type __i = 0; __i < __param.unsuccessful_size(); ++__i)
1383               {
1384                 if (__b * __aurng() < __a)
1385                   {
1386                     ++__k;
1387                     if (__k == __param.successful_size())
1388                       return __param.successful_size() - __k;
1389                     --__a;
1390                   }
1391                 --__b;
1392               }
1393             return __param.successful_size() - __k;
1394           }
1395       }
1396
1397   template<typename _UIntType>
1398     template<typename _OutputIterator,
1399              typename _UniformRandomNumberGenerator>
1400       void
1401       hypergeometric_distribution<_UIntType>::
1402       __generate_impl(_OutputIterator __f, _OutputIterator __t,
1403                       _UniformRandomNumberGenerator& __urng,
1404                       const param_type& __param)
1405       {
1406         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
1407
1408         while (__f != __t)
1409           *__f++ = this->operator()(__urng);
1410       }
1411
1412   template<typename _UIntType, typename _CharT, typename _Traits>
1413     std::basic_ostream<_CharT, _Traits>&
1414     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1415                const __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
1416     {
1417       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1418       typedef typename __ostream_type::ios_base    __ios_base;
1419
1420       const typename __ios_base::fmtflags __flags = __os.flags();
1421       const _CharT __fill = __os.fill();
1422       const std::streamsize __precision = __os.precision();
1423       const _CharT __space = __os.widen(' ');
1424       __os.flags(__ios_base::scientific | __ios_base::left);
1425       __os.fill(__space);
1426       __os.precision(std::numeric_limits<_UIntType>::max_digits10);
1427
1428       __os << __x.total_size() << __space << __x.successful_size() << __space
1429            << __x.total_draws();
1430
1431       __os.flags(__flags);
1432       __os.fill(__fill);
1433       __os.precision(__precision);
1434       return __os;
1435     }
1436
1437   template<typename _UIntType, typename _CharT, typename _Traits>
1438     std::basic_istream<_CharT, _Traits>&
1439     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1440                __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
1441     {
1442       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1443       typedef typename __istream_type::ios_base    __ios_base;
1444
1445       const typename __ios_base::fmtflags __flags = __is.flags();
1446       __is.flags(__ios_base::dec | __ios_base::skipws);
1447
1448       _UIntType __total_size, __successful_size, __total_draws;
1449       __is >> __total_size >> __successful_size >> __total_draws;
1450       __x.param(typename __gnu_cxx::hypergeometric_distribution<_UIntType>::
1451                 param_type(__total_size, __successful_size, __total_draws));
1452
1453       __is.flags(__flags);
1454       return __is;
1455     }
1456
1457
1458   template<typename _RealType>
1459     template<typename _UniformRandomNumberGenerator>
1460       typename logistic_distribution<_RealType>::result_type
1461       logistic_distribution<_RealType>::
1462       operator()(_UniformRandomNumberGenerator& __urng,
1463                  const param_type& __p)
1464       {
1465         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1466           __aurng(__urng);
1467
1468         result_type __arg = result_type(1);
1469         while (__arg == result_type(1) || __arg == result_type(0))
1470           __arg = __aurng();
1471         return __p.a()
1472              + __p.b() * std::log(__arg / (result_type(1) - __arg));
1473       }
1474
1475   template<typename _RealType>
1476     template<typename _OutputIterator,
1477              typename _UniformRandomNumberGenerator>
1478       void
1479       logistic_distribution<_RealType>::
1480       __generate_impl(_OutputIterator __f, _OutputIterator __t,
1481                       _UniformRandomNumberGenerator& __urng,
1482                       const param_type& __p)
1483       {
1484         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
1485         std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1486           __aurng(__urng);
1487
1488         while (__f != __t)
1489           {
1490             result_type __arg = result_type(1);
1491             while (__arg == result_type(1) || __arg == result_type(0))
1492               __arg = __aurng();
1493             *__f++ = __p.a()
1494                    + __p.b() * std::log(__arg / (result_type(1) - __arg));
1495           }
1496       }
1497
1498   template<typename _RealType, typename _CharT, typename _Traits>
1499     std::basic_ostream<_CharT, _Traits>&
1500     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1501                const logistic_distribution<_RealType>& __x)
1502     {
1503       typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
1504       typedef typename __ostream_type::ios_base    __ios_base;
1505
1506       const typename __ios_base::fmtflags __flags = __os.flags();
1507       const _CharT __fill = __os.fill();
1508       const std::streamsize __precision = __os.precision();
1509       const _CharT __space = __os.widen(' ');
1510       __os.flags(__ios_base::scientific | __ios_base::left);
1511       __os.fill(__space);
1512       __os.precision(std::numeric_limits<_RealType>::max_digits10);
1513
1514       __os << __x.a() << __space << __x.b();
1515
1516       __os.flags(__flags);
1517       __os.fill(__fill);
1518       __os.precision(__precision);
1519       return __os;
1520     }
1521
1522   template<typename _RealType, typename _CharT, typename _Traits>
1523     std::basic_istream<_CharT, _Traits>&
1524     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1525                logistic_distribution<_RealType>& __x)
1526     {
1527       typedef std::basic_istream<_CharT, _Traits>  __istream_type;
1528       typedef typename __istream_type::ios_base    __ios_base;
1529
1530       const typename __ios_base::fmtflags __flags = __is.flags();
1531       __is.flags(__ios_base::dec | __ios_base::skipws);
1532
1533       _RealType __a, __b;
1534       __is >> __a >> __b;
1535       __x.param(typename logistic_distribution<_RealType>::
1536                 param_type(__a, __b));
1537
1538       __is.flags(__flags);
1539       return __is;
1540     }
1541
1542
1543   template<std::size_t _Dimen, typename _RealType>
1544     template<typename _UniformRandomNumberGenerator>
1545       typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type
1546       uniform_on_sphere_distribution<_Dimen, _RealType>::
1547       operator()(_UniformRandomNumberGenerator& __urng,
1548                  const param_type& __p)
1549       {
1550         result_type __ret;
1551         _RealType __sum = _RealType(0);
1552
1553         std::generate(__ret.begin(), __ret.end(),
1554                       [&__urng, &__sum, this](){ _RealType __t = _M_nd(__urng);
1555                                                  __sum += __t * __t;
1556                                                  return __t; });
1557         auto __norm = std::sqrt(__sum);
1558         std::transform(__ret.begin(), __ret.end(), __ret.begin(),
1559                        [__norm](_RealType __val){ return __val / __norm; });
1560
1561         return __ret;
1562       }
1563
1564   template<std::size_t _Dimen, typename _RealType>
1565     template<typename _OutputIterator,
1566              typename _UniformRandomNumberGenerator>
1567       void
1568       uniform_on_sphere_distribution<_Dimen, _RealType>::
1569       __generate_impl(_OutputIterator __f, _OutputIterator __t,
1570                       _UniformRandomNumberGenerator& __urng,
1571                       const param_type& __param)
1572       {
1573         __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator>)
1574
1575         while (__f != __t)
1576           *__f++ = this->operator()(__urng, __param);
1577       }
1578
1579   template<std::size_t _Dimen, typename _RealType, typename _CharT,
1580            typename _Traits>
1581     std::basic_ostream<_CharT, _Traits>&
1582     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1583                const __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
1584                                                                _RealType>& __x)
1585     {
1586       return __os << __x._M_nd;
1587     }
1588
1589   template<std::size_t _Dimen, typename _RealType, typename _CharT,
1590            typename _Traits>
1591     std::basic_istream<_CharT, _Traits>&
1592     operator>>(std::basic_istream<_CharT, _Traits>& __is,
1593                __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
1594                                                          _RealType>& __x)
1595     {
1596       return __is >> __x._M_nd;
1597     }
1598
1599 _GLIBCXX_END_NAMESPACE_VERSION
1600 } // namespace
1601
1602
1603 #endif // _EXT_RANDOM_TCC