Home | History | Annotate | Download | only in ext
      1 // Random number extensions -*- C++ -*-
      2 
      3 // Copyright (C) 2012-2013 Free Software Foundation, Inc.
      4 //
      5 // This file is part of the GNU ISO C++ Library.  This library is free
      6 // software; you can redistribute it and/or modify it under the
      7 // terms of the GNU General Public License as published by the
      8 // Free Software Foundation; either version 3, or (at your option)
      9 // any later version.
     10 
     11 // This library is distributed in the hope that it will be useful,
     12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
     13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     14 // GNU General Public License for more details.
     15 
     16 // Under Section 7 of GPL version 3, you are granted additional
     17 // permissions described in the GCC Runtime Library Exception, version
     18 // 3.1, as published by the Free Software Foundation.
     19 
     20 // You should have received a copy of the GNU General Public License and
     21 // a copy of the GCC Runtime Library Exception along with this program;
     22 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
     23 // <http://www.gnu.org/licenses/>.
     24 
     25 /** @file ext/random
     26  *  This file is a GNU extension to the Standard C++ Library.
     27  */
     28 
     29 #ifndef _EXT_RANDOM
     30 #define _EXT_RANDOM 1
     31 
     32 #pragma GCC system_header
     33 
     34 #if __cplusplus < 201103L
     35 # include <bits/c++0x_warning.h>
     36 #else
     37 
     38 #include <random>
     39 #include <array>
     40 #include <ext/cmath>
     41 #ifdef __SSE2__
     42 # include <x86intrin.h>
     43 #endif
     44 
     45 #ifdef _GLIBCXX_USE_C99_STDINT_TR1
     46 
     47 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
     48 {
     49 _GLIBCXX_BEGIN_NAMESPACE_VERSION
     50 
     51 #if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
     52 
     53   /* Mersenne twister implementation optimized for vector operations.
     54    *
     55    * Reference: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/
     56    */
     57   template<typename _UIntType, size_t __m,
     58 	   size_t __pos1, size_t __sl1, size_t __sl2,
     59 	   size_t __sr1, size_t __sr2,
     60 	   uint32_t __msk1, uint32_t __msk2,
     61 	   uint32_t __msk3, uint32_t __msk4,
     62 	   uint32_t __parity1, uint32_t __parity2,
     63 	   uint32_t __parity3, uint32_t __parity4>
     64     class simd_fast_mersenne_twister_engine
     65     {
     66       static_assert(std::is_unsigned<_UIntType>::value, "template argument "
     67 		    "substituting _UIntType not an unsigned integral type");
     68       static_assert(__sr1 < 32, "first right shift too large");
     69       static_assert(__sr2 < 16, "second right shift too large");
     70       static_assert(__sl1 < 32, "first left shift too large");
     71       static_assert(__sl2 < 16, "second left shift too large");
     72 
     73     public:
     74       typedef _UIntType result_type;
     75 
     76     private:
     77       static constexpr size_t m_w = sizeof(result_type) * 8;
     78       static constexpr size_t _M_nstate = __m / 128 + 1;
     79       static constexpr size_t _M_nstate32 = _M_nstate * 4;
     80 
     81       static_assert(std::is_unsigned<_UIntType>::value, "template argument "
     82 		    "substituting _UIntType not an unsigned integral type");
     83       static_assert(__pos1 < _M_nstate, "POS1 not smaller than state size");
     84       static_assert(16 % sizeof(_UIntType) == 0,
     85 		    "UIntType size must divide 16");
     86 
     87     public:
     88       static constexpr size_t state_size = _M_nstate * (16
     89 							/ sizeof(result_type));
     90       static constexpr result_type default_seed = 5489u;
     91 
     92       // constructors and member function
     93       explicit
     94       simd_fast_mersenne_twister_engine(result_type __sd = default_seed)
     95       { seed(__sd); }
     96 
     97       template<typename _Sseq, typename = typename
     98 	std::enable_if<!std::is_same<_Sseq,
     99 				     simd_fast_mersenne_twister_engine>::value>
    100 	       ::type>
    101 	explicit
    102 	simd_fast_mersenne_twister_engine(_Sseq& __q)
    103 	{ seed(__q); }
    104 
    105       void
    106       seed(result_type __sd = default_seed);
    107 
    108       template<typename _Sseq>
    109 	typename std::enable_if<std::is_class<_Sseq>::value>::type
    110 	seed(_Sseq& __q);
    111 
    112       static constexpr result_type
    113       min()
    114       { return 0; };
    115 
    116       static constexpr result_type
    117       max()
    118       { return std::numeric_limits<result_type>::max(); }
    119 
    120       void
    121       discard(unsigned long long __z);
    122 
    123       result_type
    124       operator()()
    125       {
    126 	if (__builtin_expect(_M_pos >= state_size, 0))
    127 	  _M_gen_rand();
    128 
    129 	return _M_stateT[_M_pos++];
    130       }
    131 
    132       template<typename _UIntType_2, size_t __m_2,
    133 	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
    134 	       size_t __sr1_2, size_t __sr2_2,
    135 	       uint32_t __msk1_2, uint32_t __msk2_2,
    136 	       uint32_t __msk3_2, uint32_t __msk4_2,
    137 	       uint32_t __parity1_2, uint32_t __parity2_2,
    138 	       uint32_t __parity3_2, uint32_t __parity4_2>
    139 	friend bool
    140 	operator==(const simd_fast_mersenne_twister_engine<_UIntType_2,
    141 		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
    142 		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
    143 		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __lhs,
    144 		   const simd_fast_mersenne_twister_engine<_UIntType_2,
    145 		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
    146 		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
    147 		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __rhs);
    148 
    149       template<typename _UIntType_2, size_t __m_2,
    150 	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
    151 	       size_t __sr1_2, size_t __sr2_2,
    152 	       uint32_t __msk1_2, uint32_t __msk2_2,
    153 	       uint32_t __msk3_2, uint32_t __msk4_2,
    154 	       uint32_t __parity1_2, uint32_t __parity2_2,
    155 	       uint32_t __parity3_2, uint32_t __parity4_2,
    156 	       typename _CharT, typename _Traits>
    157 	friend std::basic_ostream<_CharT, _Traits>&
    158 	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    159 		   const __gnu_cxx::simd_fast_mersenne_twister_engine
    160 		   <_UIntType_2,
    161 		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
    162 		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
    163 		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
    164 
    165       template<typename _UIntType_2, size_t __m_2,
    166 	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
    167 	       size_t __sr1_2, size_t __sr2_2,
    168 	       uint32_t __msk1_2, uint32_t __msk2_2,
    169 	       uint32_t __msk3_2, uint32_t __msk4_2,
    170 	       uint32_t __parity1_2, uint32_t __parity2_2,
    171 	       uint32_t __parity3_2, uint32_t __parity4_2,
    172 	       typename _CharT, typename _Traits>
    173 	friend std::basic_istream<_CharT, _Traits>&
    174 	operator>>(std::basic_istream<_CharT, _Traits>& __is,
    175 		   __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType_2,
    176 		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
    177 		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
    178 		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
    179 
    180     private:
    181       union
    182       {
    183 #ifdef __SSE2__
    184 	__m128i _M_state[_M_nstate];
    185 #endif
    186 	uint32_t _M_state32[_M_nstate32];
    187 	result_type _M_stateT[state_size];
    188       } __attribute__ ((__aligned__ (16)));
    189       size_t _M_pos;
    190 
    191       void _M_gen_rand(void);
    192       void _M_period_certification();
    193   };
    194 
    195 
    196   template<typename _UIntType, size_t __m,
    197 	   size_t __pos1, size_t __sl1, size_t __sl2,
    198 	   size_t __sr1, size_t __sr2,
    199 	   uint32_t __msk1, uint32_t __msk2,
    200 	   uint32_t __msk3, uint32_t __msk4,
    201 	   uint32_t __parity1, uint32_t __parity2,
    202 	   uint32_t __parity3, uint32_t __parity4>
    203     inline bool
    204     operator!=(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
    205 	       __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
    206 	       __msk4, __parity1, __parity2, __parity3, __parity4>& __lhs,
    207 	       const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
    208 	       __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
    209 	       __msk4, __parity1, __parity2, __parity3, __parity4>& __rhs)
    210     { return !(__lhs == __rhs); }
    211 
    212 
    213   /* Definitions for the SIMD-oriented Fast Mersenne Twister as defined
    214    * in the C implementation by Daito and Matsumoto, as both a 32-bit
    215    * and 64-bit version.
    216    */
    217   typedef simd_fast_mersenne_twister_engine<uint32_t, 607, 2,
    218 					    15, 3, 13, 3,
    219 					    0xfdff37ffU, 0xef7f3f7dU,
    220 					    0xff777b7dU, 0x7ff7fb2fU,
    221 					    0x00000001U, 0x00000000U,
    222 					    0x00000000U, 0x5986f054U>
    223     sfmt607;
    224 
    225   typedef simd_fast_mersenne_twister_engine<uint64_t, 607, 2,
    226 					    15, 3, 13, 3,
    227 					    0xfdff37ffU, 0xef7f3f7dU,
    228 					    0xff777b7dU, 0x7ff7fb2fU,
    229 					    0x00000001U, 0x00000000U,
    230 					    0x00000000U, 0x5986f054U>
    231     sfmt607_64;
    232 
    233 
    234   typedef simd_fast_mersenne_twister_engine<uint32_t, 1279, 7,
    235 					    14, 3, 5, 1,
    236 					    0xf7fefffdU, 0x7fefcfffU,
    237 					    0xaff3ef3fU, 0xb5ffff7fU,
    238 					    0x00000001U, 0x00000000U,
    239 					    0x00000000U, 0x20000000U>
    240     sfmt1279;
    241 
    242   typedef simd_fast_mersenne_twister_engine<uint64_t, 1279, 7,
    243 					    14, 3, 5, 1,
    244 					    0xf7fefffdU, 0x7fefcfffU,
    245 					    0xaff3ef3fU, 0xb5ffff7fU,
    246 					    0x00000001U, 0x00000000U,
    247 					    0x00000000U, 0x20000000U>
    248     sfmt1279_64;
    249 
    250 
    251   typedef simd_fast_mersenne_twister_engine<uint32_t, 2281, 12,
    252 					    19, 1, 5, 1,
    253 					    0xbff7ffbfU, 0xfdfffffeU,
    254 					    0xf7ffef7fU, 0xf2f7cbbfU,
    255 					    0x00000001U, 0x00000000U,
    256 					    0x00000000U, 0x41dfa600U>
    257     sfmt2281;
    258 
    259   typedef simd_fast_mersenne_twister_engine<uint64_t, 2281, 12,
    260 					    19, 1, 5, 1,
    261 					    0xbff7ffbfU, 0xfdfffffeU,
    262 					    0xf7ffef7fU, 0xf2f7cbbfU,
    263 					    0x00000001U, 0x00000000U,
    264 					    0x00000000U, 0x41dfa600U>
    265     sfmt2281_64;
    266 
    267 
    268   typedef simd_fast_mersenne_twister_engine<uint32_t, 4253, 17,
    269 					    20, 1, 7, 1,
    270 					    0x9f7bffffU, 0x9fffff5fU,
    271 					    0x3efffffbU, 0xfffff7bbU,
    272 					    0xa8000001U, 0xaf5390a3U,
    273 					    0xb740b3f8U, 0x6c11486dU>
    274     sfmt4253;
    275 
    276   typedef simd_fast_mersenne_twister_engine<uint64_t, 4253, 17,
    277 					    20, 1, 7, 1,
    278 					    0x9f7bffffU, 0x9fffff5fU,
    279 					    0x3efffffbU, 0xfffff7bbU,
    280 					    0xa8000001U, 0xaf5390a3U,
    281 					    0xb740b3f8U, 0x6c11486dU>
    282     sfmt4253_64;
    283 
    284 
    285   typedef simd_fast_mersenne_twister_engine<uint32_t, 11213, 68,
    286 					    14, 3, 7, 3,
    287 					    0xeffff7fbU, 0xffffffefU,
    288 					    0xdfdfbfffU, 0x7fffdbfdU,
    289 					    0x00000001U, 0x00000000U,
    290 					    0xe8148000U, 0xd0c7afa3U>
    291     sfmt11213;
    292 
    293   typedef simd_fast_mersenne_twister_engine<uint64_t, 11213, 68,
    294 					    14, 3, 7, 3,
    295 					    0xeffff7fbU, 0xffffffefU,
    296 					    0xdfdfbfffU, 0x7fffdbfdU,
    297 					    0x00000001U, 0x00000000U,
    298 					    0xe8148000U, 0xd0c7afa3U>
    299     sfmt11213_64;
    300 
    301 
    302   typedef simd_fast_mersenne_twister_engine<uint32_t, 19937, 122,
    303 					    18, 1, 11, 1,
    304 					    0xdfffffefU, 0xddfecb7fU,
    305 					    0xbffaffffU, 0xbffffff6U,
    306 					    0x00000001U, 0x00000000U,
    307 					    0x00000000U, 0x13c9e684U>
    308     sfmt19937;
    309 
    310   typedef simd_fast_mersenne_twister_engine<uint64_t, 19937, 122,
    311 					    18, 1, 11, 1,
    312 					    0xdfffffefU, 0xddfecb7fU,
    313 					    0xbffaffffU, 0xbffffff6U,
    314 					    0x00000001U, 0x00000000U,
    315 					    0x00000000U, 0x13c9e684U>
    316     sfmt19937_64;
    317 
    318 
    319   typedef simd_fast_mersenne_twister_engine<uint32_t, 44497, 330,
    320 					    5, 3, 9, 3,
    321 					    0xeffffffbU, 0xdfbebfffU,
    322 					    0xbfbf7befU, 0x9ffd7bffU,
    323 					    0x00000001U, 0x00000000U,
    324 					    0xa3ac4000U, 0xecc1327aU>
    325     sfmt44497;
    326 
    327   typedef simd_fast_mersenne_twister_engine<uint64_t, 44497, 330,
    328 					    5, 3, 9, 3,
    329 					    0xeffffffbU, 0xdfbebfffU,
    330 					    0xbfbf7befU, 0x9ffd7bffU,
    331 					    0x00000001U, 0x00000000U,
    332 					    0xa3ac4000U, 0xecc1327aU>
    333     sfmt44497_64;
    334 
    335 
    336   typedef simd_fast_mersenne_twister_engine<uint32_t, 86243, 366,
    337 					    6, 7, 19, 1,
    338 					    0xfdbffbffU, 0xbff7ff3fU,
    339 					    0xfd77efffU, 0xbf9ff3ffU,
    340 					    0x00000001U, 0x00000000U,
    341 					    0x00000000U, 0xe9528d85U>
    342     sfmt86243;
    343 
    344   typedef simd_fast_mersenne_twister_engine<uint64_t, 86243, 366,
    345 					    6, 7, 19, 1,
    346 					    0xfdbffbffU, 0xbff7ff3fU,
    347 					    0xfd77efffU, 0xbf9ff3ffU,
    348 					    0x00000001U, 0x00000000U,
    349 					    0x00000000U, 0xe9528d85U>
    350     sfmt86243_64;
    351 
    352 
    353   typedef simd_fast_mersenne_twister_engine<uint32_t, 132049, 110,
    354 					    19, 1, 21, 1,
    355 					    0xffffbb5fU, 0xfb6ebf95U,
    356 					    0xfffefffaU, 0xcff77fffU,
    357 					    0x00000001U, 0x00000000U,
    358 					    0xcb520000U, 0xc7e91c7dU>
    359     sfmt132049;
    360 
    361   typedef simd_fast_mersenne_twister_engine<uint64_t, 132049, 110,
    362 					    19, 1, 21, 1,
    363 					    0xffffbb5fU, 0xfb6ebf95U,
    364 					    0xfffefffaU, 0xcff77fffU,
    365 					    0x00000001U, 0x00000000U,
    366 					    0xcb520000U, 0xc7e91c7dU>
    367     sfmt132049_64;
    368 
    369 
    370   typedef simd_fast_mersenne_twister_engine<uint32_t, 216091, 627,
    371 					    11, 3, 10, 1,
    372 					    0xbff7bff7U, 0xbfffffffU,
    373 					    0xbffffa7fU, 0xffddfbfbU,
    374 					    0xf8000001U, 0x89e80709U,
    375 					    0x3bd2b64bU, 0x0c64b1e4U>
    376     sfmt216091;
    377 
    378   typedef simd_fast_mersenne_twister_engine<uint64_t, 216091, 627,
    379 					    11, 3, 10, 1,
    380 					    0xbff7bff7U, 0xbfffffffU,
    381 					    0xbffffa7fU, 0xffddfbfbU,
    382 					    0xf8000001U, 0x89e80709U,
    383 					    0x3bd2b64bU, 0x0c64b1e4U>
    384     sfmt216091_64;
    385 
    386 #endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
    387 
    388   /**
    389    * @brief A beta continuous distribution for random numbers.
    390    *
    391    * The formula for the beta probability density function is:
    392    * @f[
    393    *     p(x|\alpha,\beta) = \frac{1}{B(\alpha,\beta)}
    394    *                         x^{\alpha - 1} (1 - x)^{\beta - 1}
    395    * @f]
    396    */
    397   template<typename _RealType = double>
    398     class beta_distribution
    399     {
    400       static_assert(std::is_floating_point<_RealType>::value,
    401 		    "template argument not a floating point type");
    402 
    403     public:
    404       /** The type of the range of the distribution. */
    405       typedef _RealType result_type;
    406       /** Parameter type. */
    407       struct param_type
    408       {
    409 	typedef beta_distribution<_RealType> distribution_type;
    410 	friend class beta_distribution<_RealType>;
    411 
    412 	explicit
    413 	param_type(_RealType __alpha_val = _RealType(1),
    414 		   _RealType __beta_val = _RealType(1))
    415 	: _M_alpha(__alpha_val), _M_beta(__beta_val)
    416 	{
    417 	  _GLIBCXX_DEBUG_ASSERT(_M_alpha > _RealType(0));
    418 	  _GLIBCXX_DEBUG_ASSERT(_M_beta > _RealType(0));
    419 	}
    420 
    421 	_RealType
    422 	alpha() const
    423 	{ return _M_alpha; }
    424 
    425 	_RealType
    426 	beta() const
    427 	{ return _M_beta; }
    428 
    429 	friend bool
    430 	operator==(const param_type& __p1, const param_type& __p2)
    431 	{ return (__p1._M_alpha == __p2._M_alpha
    432 		  && __p1._M_beta == __p2._M_beta); }
    433 
    434       private:
    435 	void
    436 	_M_initialize();
    437 
    438 	_RealType _M_alpha;
    439 	_RealType _M_beta;
    440       };
    441 
    442     public:
    443       /**
    444        * @brief Constructs a beta distribution with parameters
    445        * @f$\alpha@f$ and @f$\beta@f$.
    446        */
    447       explicit
    448       beta_distribution(_RealType __alpha_val = _RealType(1),
    449 			_RealType __beta_val = _RealType(1))
    450       : _M_param(__alpha_val, __beta_val)
    451       { }
    452 
    453       explicit
    454       beta_distribution(const param_type& __p)
    455       : _M_param(__p)
    456       { }
    457 
    458       /**
    459        * @brief Resets the distribution state.
    460        */
    461       void
    462       reset()
    463       { }
    464 
    465       /**
    466        * @brief Returns the @f$\alpha@f$ of the distribution.
    467        */
    468       _RealType
    469       alpha() const
    470       { return _M_param.alpha(); }
    471 
    472       /**
    473        * @brief Returns the @f$\beta@f$ of the distribution.
    474        */
    475       _RealType
    476       beta() const
    477       { return _M_param.beta(); }
    478 
    479       /**
    480        * @brief Returns the parameter set of the distribution.
    481        */
    482       param_type
    483       param() const
    484       { return _M_param; }
    485 
    486       /**
    487        * @brief Sets the parameter set of the distribution.
    488        * @param __param The new parameter set of the distribution.
    489        */
    490       void
    491       param(const param_type& __param)
    492       { _M_param = __param; }
    493 
    494       /**
    495        * @brief Returns the greatest lower bound value of the distribution.
    496        */
    497       result_type
    498       min() const
    499       { return result_type(0); }
    500 
    501       /**
    502        * @brief Returns the least upper bound value of the distribution.
    503        */
    504       result_type
    505       max() const
    506       { return result_type(1); }
    507 
    508       /**
    509        * @brief Generating functions.
    510        */
    511       template<typename _UniformRandomNumberGenerator>
    512 	result_type
    513 	operator()(_UniformRandomNumberGenerator& __urng)
    514 	{ return this->operator()(__urng, _M_param); }
    515 
    516       template<typename _UniformRandomNumberGenerator>
    517 	result_type
    518 	operator()(_UniformRandomNumberGenerator& __urng,
    519 		   const param_type& __p);
    520 
    521       template<typename _ForwardIterator,
    522 	       typename _UniformRandomNumberGenerator>
    523 	void
    524 	__generate(_ForwardIterator __f, _ForwardIterator __t,
    525 		   _UniformRandomNumberGenerator& __urng)
    526 	{ this->__generate(__f, __t, __urng, _M_param); }
    527 
    528       template<typename _ForwardIterator,
    529 	       typename _UniformRandomNumberGenerator>
    530 	void
    531 	__generate(_ForwardIterator __f, _ForwardIterator __t,
    532 		   _UniformRandomNumberGenerator& __urng,
    533 		   const param_type& __p)
    534 	{ this->__generate_impl(__f, __t, __urng, __p); }
    535 
    536       template<typename _UniformRandomNumberGenerator>
    537 	void
    538 	__generate(result_type* __f, result_type* __t,
    539 		   _UniformRandomNumberGenerator& __urng,
    540 		   const param_type& __p)
    541 	{ this->__generate_impl(__f, __t, __urng, __p); }
    542 
    543       /**
    544        * @brief Return true if two beta distributions have the same
    545        *        parameters and the sequences that would be generated
    546        *        are equal.
    547        */
    548       friend bool
    549       operator==(const beta_distribution& __d1,
    550 		 const beta_distribution& __d2)
    551       { return __d1._M_param == __d2._M_param; }
    552 
    553       /**
    554        * @brief Inserts a %beta_distribution random number distribution
    555        * @p __x into the output stream @p __os.
    556        *
    557        * @param __os An output stream.
    558        * @param __x  A %beta_distribution random number distribution.
    559        *
    560        * @returns The output stream with the state of @p __x inserted or in
    561        * an error state.
    562        */
    563       template<typename _RealType1, typename _CharT, typename _Traits>
    564 	friend std::basic_ostream<_CharT, _Traits>&
    565 	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    566 		   const __gnu_cxx::beta_distribution<_RealType1>& __x);
    567 
    568       /**
    569        * @brief Extracts a %beta_distribution random number distribution
    570        * @p __x from the input stream @p __is.
    571        *
    572        * @param __is An input stream.
    573        * @param __x  A %beta_distribution random number generator engine.
    574        *
    575        * @returns The input stream with @p __x extracted or in an error state.
    576        */
    577       template<typename _RealType1, typename _CharT, typename _Traits>
    578 	friend std::basic_istream<_CharT, _Traits>&
    579 	operator>>(std::basic_istream<_CharT, _Traits>& __is,
    580 		   __gnu_cxx::beta_distribution<_RealType1>& __x);
    581 
    582     private:
    583       template<typename _ForwardIterator,
    584 	       typename _UniformRandomNumberGenerator>
    585 	void
    586 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
    587 			_UniformRandomNumberGenerator& __urng,
    588 			const param_type& __p);
    589 
    590       param_type _M_param;
    591     };
    592 
    593   /**
    594    * @brief Return true if two beta distributions are different.
    595    */
    596   template<typename _RealType>
    597     inline bool
    598     operator!=(const __gnu_cxx::beta_distribution<_RealType>& __d1,
    599 	       const __gnu_cxx::beta_distribution<_RealType>& __d2)
    600    { return !(__d1 == __d2); }
    601 
    602 
    603   /**
    604    * @brief A multi-variate normal continuous distribution for random numbers.
    605    *
    606    * The formula for the normal probability density function is
    607    * @f[
    608    *     p(\overrightarrow{x}|\overrightarrow{\mu },\Sigma) =
    609    *       \frac{1}{\sqrt{(2\pi )^k\det(\Sigma))}}
    610    *       e^{-\frac{1}{2}(\overrightarrow{x}-\overrightarrow{\mu})^\text{T}
    611    *          \Sigma ^{-1}(\overrightarrow{x}-\overrightarrow{\mu})}
    612    * @f]
    613    *
    614    * where @f$\overrightarrow{x}@f$ and @f$\overrightarrow{\mu}@f$ are
    615    * vectors of dimension @f$k@f$ and @f$\Sigma@f$ is the covariance
    616    * matrix (which must be positive-definite).
    617    */
    618   template<std::size_t _Dimen, typename _RealType = double>
    619     class normal_mv_distribution
    620     {
    621       static_assert(std::is_floating_point<_RealType>::value,
    622 		    "template argument not a floating point type");
    623       static_assert(_Dimen != 0, "dimension is zero");
    624 
    625     public:
    626       /** The type of the range of the distribution. */
    627       typedef std::array<_RealType, _Dimen> result_type;
    628       /** Parameter type. */
    629       class param_type
    630       {
    631 	static constexpr size_t _M_t_size = _Dimen * (_Dimen + 1) / 2;
    632 
    633       public:
    634 	typedef normal_mv_distribution<_Dimen, _RealType> distribution_type;
    635 	friend class normal_mv_distribution<_Dimen, _RealType>;
    636 
    637 	param_type()
    638 	{
    639 	  std::fill(_M_mean.begin(), _M_mean.end(), _RealType(0));
    640 	  auto __it = _M_t.begin();
    641 	  for (size_t __i = 0; __i < _Dimen; ++__i)
    642 	    {
    643 	      std::fill_n(__it, __i, _RealType(0));
    644 	      __it += __i;
    645 	      *__it++ = _RealType(1);
    646 	    }
    647 	}
    648 
    649 	template<typename _ForwardIterator1, typename _ForwardIterator2>
    650 	  param_type(_ForwardIterator1 __meanbegin,
    651 		     _ForwardIterator1 __meanend,
    652 		     _ForwardIterator2 __varcovbegin,
    653 		     _ForwardIterator2 __varcovend)
    654 	{
    655 	  __glibcxx_function_requires(_ForwardIteratorConcept<
    656 				      _ForwardIterator1>)
    657 	  __glibcxx_function_requires(_ForwardIteratorConcept<
    658 				      _ForwardIterator2>)
    659 	  _GLIBCXX_DEBUG_ASSERT(std::distance(__meanbegin, __meanend)
    660 				<= _Dimen);
    661 	  const auto __dist = std::distance(__varcovbegin, __varcovend);
    662 	  _GLIBCXX_DEBUG_ASSERT(__dist == _Dimen * _Dimen
    663 				|| __dist == _Dimen * (_Dimen + 1) / 2
    664 				|| __dist == _Dimen);
    665 
    666 	  if (__dist == _Dimen * _Dimen)
    667 	    _M_init_full(__meanbegin, __meanend, __varcovbegin, __varcovend);
    668 	  else if (__dist == _Dimen * (_Dimen + 1) / 2)
    669 	    _M_init_lower(__meanbegin, __meanend, __varcovbegin, __varcovend);
    670 	  else
    671 	    _M_init_diagonal(__meanbegin, __meanend,
    672 			     __varcovbegin, __varcovend);
    673 	}
    674 
    675 	param_type(std::initializer_list<_RealType> __mean,
    676 		   std::initializer_list<_RealType> __varcov)
    677 	{
    678 	  _GLIBCXX_DEBUG_ASSERT(__mean.size() <= _Dimen);
    679 	  _GLIBCXX_DEBUG_ASSERT(__varcov.size() == _Dimen * _Dimen
    680 				|| __varcov.size() == _Dimen * (_Dimen + 1) / 2
    681 				|| __varcov.size() == _Dimen);
    682 
    683 	  if (__varcov.size() == _Dimen * _Dimen)
    684 	    _M_init_full(__mean.begin(), __mean.end(),
    685 			 __varcov.begin(), __varcov.end());
    686 	  else if (__varcov.size() == _Dimen * (_Dimen + 1) / 2)
    687 	    _M_init_lower(__mean.begin(), __mean.end(),
    688 			  __varcov.begin(), __varcov.end());
    689 	  else
    690 	    _M_init_diagonal(__mean.begin(), __mean.end(),
    691 			     __varcov.begin(), __varcov.end());
    692 	}
    693 
    694 	std::array<_RealType, _Dimen>
    695 	mean() const
    696 	{ return _M_mean; }
    697 
    698 	std::array<_RealType, _M_t_size>
    699 	varcov() const
    700 	{ return _M_t; }
    701 
    702 	friend bool
    703 	operator==(const param_type& __p1, const param_type& __p2)
    704 	{ return __p1._M_mean == __p2._M_mean && __p1._M_t == __p2._M_t; }
    705 
    706       private:
    707 	template <typename _InputIterator1, typename _InputIterator2>
    708 	  void _M_init_full(_InputIterator1 __meanbegin,
    709 			    _InputIterator1 __meanend,
    710 			    _InputIterator2 __varcovbegin,
    711 			    _InputIterator2 __varcovend);
    712 	template <typename _InputIterator1, typename _InputIterator2>
    713 	  void _M_init_lower(_InputIterator1 __meanbegin,
    714 			     _InputIterator1 __meanend,
    715 			     _InputIterator2 __varcovbegin,
    716 			     _InputIterator2 __varcovend);
    717 	template <typename _InputIterator1, typename _InputIterator2>
    718 	  void _M_init_diagonal(_InputIterator1 __meanbegin,
    719 				_InputIterator1 __meanend,
    720 				_InputIterator2 __varbegin,
    721 				_InputIterator2 __varend);
    722 
    723 	std::array<_RealType, _Dimen> _M_mean;
    724 	std::array<_RealType, _M_t_size> _M_t;
    725       };
    726 
    727     public:
    728       normal_mv_distribution()
    729       : _M_param(), _M_nd()
    730       { }
    731 
    732       template<typename _ForwardIterator1, typename _ForwardIterator2>
    733 	normal_mv_distribution(_ForwardIterator1 __meanbegin,
    734 			       _ForwardIterator1 __meanend,
    735 			       _ForwardIterator2 __varcovbegin,
    736 			       _ForwardIterator2 __varcovend)
    737 	: _M_param(__meanbegin, __meanend, __varcovbegin, __varcovend),
    738 	  _M_nd()
    739 	{ }
    740 
    741       normal_mv_distribution(std::initializer_list<_RealType> __mean,
    742 			     std::initializer_list<_RealType> __varcov)
    743       : _M_param(__mean, __varcov), _M_nd()
    744       { }
    745 
    746       explicit
    747       normal_mv_distribution(const param_type& __p)
    748       : _M_param(__p), _M_nd()
    749       { }
    750 
    751       /**
    752        * @brief Resets the distribution state.
    753        */
    754       void
    755       reset()
    756       { _M_nd.reset(); }
    757 
    758       /**
    759        * @brief Returns the mean of the distribution.
    760        */
    761       result_type
    762       mean() const
    763       { return _M_param.mean(); }
    764 
    765       /**
    766        * @brief Returns the compact form of the variance/covariance
    767        * matrix of the distribution.
    768        */
    769       std::array<_RealType, _Dimen * (_Dimen + 1) / 2>
    770       varcov() const
    771       { return _M_param.varcov(); }
    772 
    773       /**
    774        * @brief Returns the parameter set of the distribution.
    775        */
    776       param_type
    777       param() const
    778       { return _M_param; }
    779 
    780       /**
    781        * @brief Sets the parameter set of the distribution.
    782        * @param __param The new parameter set of the distribution.
    783        */
    784       void
    785       param(const param_type& __param)
    786       { _M_param = __param; }
    787 
    788       /**
    789        * @brief Returns the greatest lower bound value of the distribution.
    790        */
    791       result_type
    792       min() const
    793       { result_type __res;
    794 	__res.fill(std::numeric_limits<_RealType>::lowest());
    795 	return __res; }
    796 
    797       /**
    798        * @brief Returns the least upper bound value of the distribution.
    799        */
    800       result_type
    801       max() const
    802       { result_type __res;
    803 	__res.fill(std::numeric_limits<_RealType>::max());
    804 	return __res; }
    805 
    806       /**
    807        * @brief Generating functions.
    808        */
    809       template<typename _UniformRandomNumberGenerator>
    810 	result_type
    811 	operator()(_UniformRandomNumberGenerator& __urng)
    812 	{ return this->operator()(__urng, _M_param); }
    813 
    814       template<typename _UniformRandomNumberGenerator>
    815 	result_type
    816 	operator()(_UniformRandomNumberGenerator& __urng,
    817 		   const param_type& __p);
    818 
    819       template<typename _ForwardIterator,
    820 	       typename _UniformRandomNumberGenerator>
    821 	void
    822 	__generate(_ForwardIterator __f, _ForwardIterator __t,
    823 		   _UniformRandomNumberGenerator& __urng)
    824 	{ return this->__generate_impl(__f, __t, __urng, _M_param); }
    825 
    826       template<typename _ForwardIterator,
    827 	       typename _UniformRandomNumberGenerator>
    828 	void
    829 	__generate(_ForwardIterator __f, _ForwardIterator __t,
    830 		   _UniformRandomNumberGenerator& __urng,
    831 		   const param_type& __p)
    832 	{ return this->__generate_impl(__f, __t, __urng, __p); }
    833 
    834       /**
    835        * @brief Return true if two multi-variant normal distributions have
    836        *        the same parameters and the sequences that would
    837        *        be generated are equal.
    838        */
    839       template<size_t _Dimen1, typename _RealType1>
    840 	friend bool
    841 	operator==(const
    842 		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
    843 		   __d1,
    844 		   const
    845 		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
    846 		   __d2);
    847 
    848       /**
    849        * @brief Inserts a %normal_mv_distribution random number distribution
    850        * @p __x into the output stream @p __os.
    851        *
    852        * @param __os An output stream.
    853        * @param __x  A %normal_mv_distribution random number distribution.
    854        *
    855        * @returns The output stream with the state of @p __x inserted or in
    856        * an error state.
    857        */
    858       template<size_t _Dimen1, typename _RealType1,
    859 	       typename _CharT, typename _Traits>
    860 	friend std::basic_ostream<_CharT, _Traits>&
    861 	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
    862 		   const
    863 		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
    864 		   __x);
    865 
    866       /**
    867        * @brief Extracts a %normal_mv_distribution random number distribution
    868        * @p __x from the input stream @p __is.
    869        *
    870        * @param __is An input stream.
    871        * @param __x  A %normal_mv_distribution random number generator engine.
    872        *
    873        * @returns The input stream with @p __x extracted or in an error
    874        *          state.
    875        */
    876       template<size_t _Dimen1, typename _RealType1,
    877 	       typename _CharT, typename _Traits>
    878 	friend std::basic_istream<_CharT, _Traits>&
    879 	operator>>(std::basic_istream<_CharT, _Traits>& __is,
    880 		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
    881 		   __x);
    882 
    883     private:
    884       template<typename _ForwardIterator,
    885 	       typename _UniformRandomNumberGenerator>
    886 	void
    887 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
    888 			_UniformRandomNumberGenerator& __urng,
    889 			const param_type& __p);
    890 
    891       param_type _M_param;
    892       std::normal_distribution<_RealType> _M_nd;
    893   };
    894 
    895   /**
    896    * @brief Return true if two multi-variate normal distributions are
    897    * different.
    898    */
    899   template<size_t _Dimen, typename _RealType>
    900     inline bool
    901     operator!=(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
    902 	       __d1,
    903 	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
    904 	       __d2)
    905     { return !(__d1 == __d2); }
    906 
    907 
    908   /**
    909    * @brief A Rice continuous distribution for random numbers.
    910    *
    911    * The formula for the Rice probability density function is
    912    * @f[
    913    *     p(x|\nu,\sigma) = \frac{x}{\sigma^2}
    914    *                       \exp\left(-\frac{x^2+\nu^2}{2\sigma^2}\right)
    915    *                       I_0\left(\frac{x \nu}{\sigma^2}\right)
    916    * @f]
    917    * where @f$I_0(z)@f$ is the modified Bessel function of the first kind
    918    * of order 0 and @f$\nu >= 0@f$ and @f$\sigma > 0@f$.
    919    *
    920    * <table border=1 cellpadding=10 cellspacing=0>
    921    * <caption align=top>Distribution Statistics</caption>
    922    * <tr><td>Mean</td><td>@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
    923    * <tr><td>Variance</td><td>@f$2\sigma^2 + \nu^2
    924    *                   + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
    925    * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
    926    * </table>
    927    * where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2.
    928    */
    929   template<typename _RealType = double>
    930     class
    931     rice_distribution
    932     {
    933       static_assert(std::is_floating_point<_RealType>::value,
    934 		    "template argument not a floating point type");
    935     public:
    936       /** The type of the range of the distribution. */
    937       typedef _RealType result_type;
    938       /** Parameter type. */
    939       struct param_type
    940       {
    941 	typedef rice_distribution<result_type> distribution_type;
    942 
    943 	param_type(result_type __nu_val = result_type(0),
    944 		   result_type __sigma_val = result_type(1))
    945 	: _M_nu(__nu_val), _M_sigma(__sigma_val)
    946 	{
    947 	  _GLIBCXX_DEBUG_ASSERT(_M_nu >= result_type(0));
    948 	  _GLIBCXX_DEBUG_ASSERT(_M_sigma > result_type(0));
    949 	}
    950 
    951 	result_type
    952 	nu() const
    953 	{ return _M_nu; }
    954 
    955 	result_type
    956 	sigma() const
    957 	{ return _M_sigma; }
    958 
    959 	friend bool
    960 	operator==(const param_type& __p1, const param_type& __p2)
    961 	{ return __p1._M_nu == __p2._M_nu
    962 	      && __p1._M_sigma == __p2._M_sigma; }
    963 
    964       private:
    965 	void _M_initialize();
    966 
    967 	result_type _M_nu;
    968 	result_type _M_sigma;
    969       };
    970 
    971       /**
    972        * @brief Constructors.
    973        */
    974       explicit
    975       rice_distribution(result_type __nu_val = result_type(0),
    976 			result_type __sigma_val = result_type(1))
    977       : _M_param(__nu_val, __sigma_val),
    978 	_M_ndx(__nu_val, __sigma_val),
    979 	_M_ndy(result_type(0), __sigma_val)
    980       { }
    981 
    982       explicit
    983       rice_distribution(const param_type& __p)
    984       : _M_param(__p),
    985 	_M_ndx(__p.nu(), __p.sigma()),
    986 	_M_ndy(result_type(0), __p.sigma())
    987       { }
    988 
    989       /**
    990        * @brief Resets the distribution state.
    991        */
    992       void
    993       reset()
    994       {
    995 	_M_ndx.reset();
    996 	_M_ndy.reset();
    997       }
    998 
    999       /**
   1000        * @brief Return the parameters of the distribution.
   1001        */
   1002       result_type
   1003       nu() const
   1004       { return _M_param.nu(); }
   1005 
   1006       result_type
   1007       sigma() const
   1008       { return _M_param.sigma(); }
   1009 
   1010       /**
   1011        * @brief Returns the parameter set of the distribution.
   1012        */
   1013       param_type
   1014       param() const
   1015       { return _M_param; }
   1016 
   1017       /**
   1018        * @brief Sets the parameter set of the distribution.
   1019        * @param __param The new parameter set of the distribution.
   1020        */
   1021       void
   1022       param(const param_type& __param)
   1023       { _M_param = __param; }
   1024 
   1025       /**
   1026        * @brief Returns the greatest lower bound value of the distribution.
   1027        */
   1028       result_type
   1029       min() const
   1030       { return result_type(0); }
   1031 
   1032       /**
   1033        * @brief Returns the least upper bound value of the distribution.
   1034        */
   1035       result_type
   1036       max() const
   1037       { return std::numeric_limits<result_type>::max(); }
   1038 
   1039       /**
   1040        * @brief Generating functions.
   1041        */
   1042       template<typename _UniformRandomNumberGenerator>
   1043 	result_type
   1044 	operator()(_UniformRandomNumberGenerator& __urng)
   1045 	{
   1046 	  result_type __x = this->_M_ndx(__urng);
   1047 	  result_type __y = this->_M_ndy(__urng);
   1048 #if _GLIBCXX_USE_C99_MATH_TR1
   1049 	  return std::hypot(__x, __y);
   1050 #else
   1051 	  return std::sqrt(__x * __x + __y * __y);
   1052 #endif
   1053 	}
   1054 
   1055       template<typename _UniformRandomNumberGenerator>
   1056 	result_type
   1057 	operator()(_UniformRandomNumberGenerator& __urng,
   1058 		   const param_type& __p)
   1059 	{
   1060 	  typename std::normal_distribution<result_type>::param_type
   1061 	    __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
   1062 	  result_type __x = this->_M_ndx(__px, __urng);
   1063 	  result_type __y = this->_M_ndy(__py, __urng);
   1064 #if _GLIBCXX_USE_C99_MATH_TR1
   1065 	  return std::hypot(__x, __y);
   1066 #else
   1067 	  return std::sqrt(__x * __x + __y * __y);
   1068 #endif
   1069 	}
   1070 
   1071       template<typename _ForwardIterator,
   1072 	       typename _UniformRandomNumberGenerator>
   1073 	void
   1074 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1075 		   _UniformRandomNumberGenerator& __urng)
   1076 	{ this->__generate(__f, __t, __urng, _M_param); }
   1077 
   1078       template<typename _ForwardIterator,
   1079 	       typename _UniformRandomNumberGenerator>
   1080 	void
   1081 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1082 		   _UniformRandomNumberGenerator& __urng,
   1083 		   const param_type& __p)
   1084 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1085 
   1086       template<typename _UniformRandomNumberGenerator>
   1087 	void
   1088 	__generate(result_type* __f, result_type* __t,
   1089 		   _UniformRandomNumberGenerator& __urng,
   1090 		   const param_type& __p)
   1091 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1092 
   1093       /**
   1094        * @brief Return true if two Rice distributions have
   1095        *        the same parameters and the sequences that would
   1096        *        be generated are equal.
   1097        */
   1098       friend bool
   1099       operator==(const rice_distribution& __d1,
   1100 		 const rice_distribution& __d2)
   1101       { return (__d1._M_param == __d2._M_param
   1102 		&& __d1._M_ndx == __d2._M_ndx
   1103 		&& __d1._M_ndy == __d2._M_ndy); }
   1104 
   1105       /**
   1106        * @brief Inserts a %rice_distribution random number distribution
   1107        * @p __x into the output stream @p __os.
   1108        *
   1109        * @param __os An output stream.
   1110        * @param __x  A %rice_distribution random number distribution.
   1111        *
   1112        * @returns The output stream with the state of @p __x inserted or in
   1113        * an error state.
   1114        */
   1115       template<typename _RealType1, typename _CharT, typename _Traits>
   1116 	friend std::basic_ostream<_CharT, _Traits>&
   1117 	operator<<(std::basic_ostream<_CharT, _Traits>&,
   1118 		   const rice_distribution<_RealType1>&);
   1119 
   1120       /**
   1121        * @brief Extracts a %rice_distribution random number distribution
   1122        * @p __x from the input stream @p __is.
   1123        *
   1124        * @param __is An input stream.
   1125        * @param __x A %rice_distribution random number
   1126        *            generator engine.
   1127        *
   1128        * @returns The input stream with @p __x extracted or in an error state.
   1129        */
   1130       template<typename _RealType1, typename _CharT, typename _Traits>
   1131 	friend std::basic_istream<_CharT, _Traits>&
   1132 	operator>>(std::basic_istream<_CharT, _Traits>&,
   1133 		   rice_distribution<_RealType1>&);
   1134 
   1135     private:
   1136       template<typename _ForwardIterator,
   1137 	       typename _UniformRandomNumberGenerator>
   1138 	void
   1139 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   1140 			_UniformRandomNumberGenerator& __urng,
   1141 			const param_type& __p);
   1142 
   1143       param_type _M_param;
   1144 
   1145       std::normal_distribution<result_type> _M_ndx;
   1146       std::normal_distribution<result_type> _M_ndy;
   1147     };
   1148 
   1149   /**
   1150    * @brief Return true if two Rice distributions are not equal.
   1151    */
   1152   template<typename _RealType1>
   1153     inline bool
   1154     operator!=(const rice_distribution<_RealType1>& __d1,
   1155 	       const rice_distribution<_RealType1>& __d2)
   1156     { return !(__d1 == __d2); }
   1157 
   1158 
   1159   /**
   1160    * @brief A Nakagami continuous distribution for random numbers.
   1161    *
   1162    * The formula for the Nakagami probability density function is
   1163    * @f[
   1164    *     p(x|\mu,\omega) = \frac{2\mu^\mu}{\Gamma(\mu)\omega^\mu}
   1165    *                       x^{2\mu-1}e^{-\mu x / \omega}
   1166    * @f]
   1167    * where @f$\Gamma(z)@f$ is the gamma function and @f$\mu >= 0.5@f$
   1168    * and @f$\omega > 0@f$.
   1169    */
   1170   template<typename _RealType = double>
   1171     class
   1172     nakagami_distribution
   1173     {
   1174       static_assert(std::is_floating_point<_RealType>::value,
   1175 		    "template argument not a floating point type");
   1176 
   1177     public:
   1178       /** The type of the range of the distribution. */
   1179       typedef _RealType result_type;
   1180       /** Parameter type. */
   1181       struct param_type
   1182       {
   1183 	typedef nakagami_distribution<result_type> distribution_type;
   1184 
   1185 	param_type(result_type __mu_val = result_type(1),
   1186 		   result_type __omega_val = result_type(1))
   1187 	: _M_mu(__mu_val), _M_omega(__omega_val)
   1188 	{
   1189 	  _GLIBCXX_DEBUG_ASSERT(_M_mu >= result_type(0.5L));
   1190 	  _GLIBCXX_DEBUG_ASSERT(_M_omega > result_type(0));
   1191 	}
   1192 
   1193 	result_type
   1194 	mu() const
   1195 	{ return _M_mu; }
   1196 
   1197 	result_type
   1198 	omega() const
   1199 	{ return _M_omega; }
   1200 
   1201 	friend bool
   1202 	operator==(const param_type& __p1, const param_type& __p2)
   1203 	{ return __p1._M_mu == __p2._M_mu
   1204 	      && __p1._M_omega == __p2._M_omega; }
   1205 
   1206       private:
   1207 	void _M_initialize();
   1208 
   1209 	result_type _M_mu;
   1210 	result_type _M_omega;
   1211       };
   1212 
   1213       /**
   1214        * @brief Constructors.
   1215        */
   1216       explicit
   1217       nakagami_distribution(result_type __mu_val = result_type(1),
   1218 			    result_type __omega_val = result_type(1))
   1219       : _M_param(__mu_val, __omega_val),
   1220 	_M_gd(__mu_val, __omega_val / __mu_val)
   1221       { }
   1222 
   1223       explicit
   1224       nakagami_distribution(const param_type& __p)
   1225       : _M_param(__p),
   1226 	_M_gd(__p.mu(), __p.omega() / __p.mu())
   1227       { }
   1228 
   1229       /**
   1230        * @brief Resets the distribution state.
   1231        */
   1232       void
   1233       reset()
   1234       { _M_gd.reset(); }
   1235 
   1236       /**
   1237        * @brief Return the parameters of the distribution.
   1238        */
   1239       result_type
   1240       mu() const
   1241       { return _M_param.mu(); }
   1242 
   1243       result_type
   1244       omega() const
   1245       { return _M_param.omega(); }
   1246 
   1247       /**
   1248        * @brief Returns the parameter set of the distribution.
   1249        */
   1250       param_type
   1251       param() const
   1252       { return _M_param; }
   1253 
   1254       /**
   1255        * @brief Sets the parameter set of the distribution.
   1256        * @param __param The new parameter set of the distribution.
   1257        */
   1258       void
   1259       param(const param_type& __param)
   1260       { _M_param = __param; }
   1261 
   1262       /**
   1263        * @brief Returns the greatest lower bound value of the distribution.
   1264        */
   1265       result_type
   1266       min() const
   1267       { return result_type(0); }
   1268 
   1269       /**
   1270        * @brief Returns the least upper bound value of the distribution.
   1271        */
   1272       result_type
   1273       max() const
   1274       { return std::numeric_limits<result_type>::max(); }
   1275 
   1276       /**
   1277        * @brief Generating functions.
   1278        */
   1279       template<typename _UniformRandomNumberGenerator>
   1280 	result_type
   1281 	operator()(_UniformRandomNumberGenerator& __urng)
   1282 	{ return std::sqrt(this->_M_gd(__urng)); }
   1283 
   1284       template<typename _UniformRandomNumberGenerator>
   1285 	result_type
   1286 	operator()(_UniformRandomNumberGenerator& __urng,
   1287 		   const param_type& __p)
   1288 	{
   1289 	  typename std::gamma_distribution<result_type>::param_type
   1290 	    __pg(__p.mu(), __p.omega() / __p.mu());
   1291 	  return std::sqrt(this->_M_gd(__pg, __urng));
   1292 	}
   1293 
   1294       template<typename _ForwardIterator,
   1295 	       typename _UniformRandomNumberGenerator>
   1296 	void
   1297 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1298 		   _UniformRandomNumberGenerator& __urng)
   1299 	{ this->__generate(__f, __t, __urng, _M_param); }
   1300 
   1301       template<typename _ForwardIterator,
   1302 	       typename _UniformRandomNumberGenerator>
   1303 	void
   1304 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1305 		   _UniformRandomNumberGenerator& __urng,
   1306 		   const param_type& __p)
   1307 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1308 
   1309       template<typename _UniformRandomNumberGenerator>
   1310 	void
   1311 	__generate(result_type* __f, result_type* __t,
   1312 		   _UniformRandomNumberGenerator& __urng,
   1313 		   const param_type& __p)
   1314 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1315 
   1316       /**
   1317        * @brief Return true if two Nakagami distributions have
   1318        *        the same parameters and the sequences that would
   1319        *        be generated are equal.
   1320        */
   1321       friend bool
   1322       operator==(const nakagami_distribution& __d1,
   1323 		 const nakagami_distribution& __d2)
   1324       { return (__d1._M_param == __d2._M_param
   1325 		&& __d1._M_gd == __d2._M_gd); }
   1326 
   1327       /**
   1328        * @brief Inserts a %nakagami_distribution random number distribution
   1329        * @p __x into the output stream @p __os.
   1330        *
   1331        * @param __os An output stream.
   1332        * @param __x  A %nakagami_distribution random number distribution.
   1333        *
   1334        * @returns The output stream with the state of @p __x inserted or in
   1335        * an error state.
   1336        */
   1337       template<typename _RealType1, typename _CharT, typename _Traits>
   1338 	friend std::basic_ostream<_CharT, _Traits>&
   1339 	operator<<(std::basic_ostream<_CharT, _Traits>&,
   1340 		   const nakagami_distribution<_RealType1>&);
   1341 
   1342       /**
   1343        * @brief Extracts a %nakagami_distribution random number distribution
   1344        * @p __x from the input stream @p __is.
   1345        *
   1346        * @param __is An input stream.
   1347        * @param __x A %nakagami_distribution random number
   1348        *            generator engine.
   1349        *
   1350        * @returns The input stream with @p __x extracted or in an error state.
   1351        */
   1352       template<typename _RealType1, typename _CharT, typename _Traits>
   1353 	friend std::basic_istream<_CharT, _Traits>&
   1354 	operator>>(std::basic_istream<_CharT, _Traits>&,
   1355 		   nakagami_distribution<_RealType1>&);
   1356 
   1357     private:
   1358       template<typename _ForwardIterator,
   1359 	       typename _UniformRandomNumberGenerator>
   1360 	void
   1361 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   1362 			_UniformRandomNumberGenerator& __urng,
   1363 			const param_type& __p);
   1364 
   1365       param_type _M_param;
   1366 
   1367       std::gamma_distribution<result_type> _M_gd;
   1368     };
   1369 
   1370   /**
   1371    * @brief Return true if two Nakagami distributions are not equal.
   1372    */
   1373   template<typename _RealType>
   1374     inline bool
   1375     operator!=(const nakagami_distribution<_RealType>& __d1,
   1376 	       const nakagami_distribution<_RealType>& __d2)
   1377     { return !(__d1 == __d2); }
   1378 
   1379 
   1380   /**
   1381    * @brief A Pareto continuous distribution for random numbers.
   1382    *
   1383    * The formula for the Pareto cumulative probability function is
   1384    * @f[
   1385    *     P(x|\alpha,\mu) = 1 - \left(\frac{\mu}{x}\right)^\alpha
   1386    * @f]
   1387    * The formula for the Pareto probability density function is
   1388    * @f[
   1389    *     p(x|\alpha,\mu) = \frac{\alpha + 1}{\mu}
   1390    *                       \left(\frac{\mu}{x}\right)^{\alpha + 1}
   1391    * @f]
   1392    * where @f$x >= \mu@f$ and @f$\mu > 0@f$, @f$\alpha > 0@f$.
   1393    *
   1394    * <table border=1 cellpadding=10 cellspacing=0>
   1395    * <caption align=top>Distribution Statistics</caption>
   1396    * <tr><td>Mean</td><td>@f$\alpha \mu / (\alpha - 1)@f$
   1397    *              for @f$\alpha > 1@f$</td></tr>
   1398    * <tr><td>Variance</td><td>@f$\alpha \mu^2 / [(\alpha - 1)^2(\alpha - 2)]@f$
   1399    *              for @f$\alpha > 2@f$</td></tr>
   1400    * <tr><td>Range</td><td>@f$[\mu, \infty)@f$</td></tr>
   1401    * </table>
   1402    */
   1403   template<typename _RealType = double>
   1404     class
   1405     pareto_distribution
   1406     {
   1407       static_assert(std::is_floating_point<_RealType>::value,
   1408 		    "template argument not a floating point type");
   1409 
   1410     public:
   1411       /** The type of the range of the distribution. */
   1412       typedef _RealType result_type;
   1413       /** Parameter type. */
   1414       struct param_type
   1415       {
   1416 	typedef pareto_distribution<result_type> distribution_type;
   1417 
   1418 	param_type(result_type __alpha_val = result_type(1),
   1419 		   result_type __mu_val = result_type(1))
   1420 	: _M_alpha(__alpha_val), _M_mu(__mu_val)
   1421 	{
   1422 	  _GLIBCXX_DEBUG_ASSERT(_M_alpha > result_type(0));
   1423 	  _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0));
   1424 	}
   1425 
   1426 	result_type
   1427 	alpha() const
   1428 	{ return _M_alpha; }
   1429 
   1430 	result_type
   1431 	mu() const
   1432 	{ return _M_mu; }
   1433 
   1434 	friend bool
   1435 	operator==(const param_type& __p1, const param_type& __p2)
   1436 	{ return __p1._M_alpha == __p2._M_alpha && __p1._M_mu == __p2._M_mu; }
   1437 
   1438       private:
   1439 	void _M_initialize();
   1440 
   1441 	result_type _M_alpha;
   1442 	result_type _M_mu;
   1443       };
   1444 
   1445       /**
   1446        * @brief Constructors.
   1447        */
   1448       explicit
   1449       pareto_distribution(result_type __alpha_val = result_type(1),
   1450 			  result_type __mu_val = result_type(1))
   1451       : _M_param(__alpha_val, __mu_val),
   1452 	_M_ud()
   1453       { }
   1454 
   1455       explicit
   1456       pareto_distribution(const param_type& __p)
   1457       : _M_param(__p),
   1458 	_M_ud()
   1459       { }
   1460 
   1461       /**
   1462        * @brief Resets the distribution state.
   1463        */
   1464       void
   1465       reset()
   1466       {
   1467 	_M_ud.reset();
   1468       }
   1469 
   1470       /**
   1471        * @brief Return the parameters of the distribution.
   1472        */
   1473       result_type
   1474       alpha() const
   1475       { return _M_param.alpha(); }
   1476 
   1477       result_type
   1478       mu() const
   1479       { return _M_param.mu(); }
   1480 
   1481       /**
   1482        * @brief Returns the parameter set of the distribution.
   1483        */
   1484       param_type
   1485       param() const
   1486       { return _M_param; }
   1487 
   1488       /**
   1489        * @brief Sets the parameter set of the distribution.
   1490        * @param __param The new parameter set of the distribution.
   1491        */
   1492       void
   1493       param(const param_type& __param)
   1494       { _M_param = __param; }
   1495 
   1496       /**
   1497        * @brief Returns the greatest lower bound value of the distribution.
   1498        */
   1499       result_type
   1500       min() const
   1501       { return this->mu(); }
   1502 
   1503       /**
   1504        * @brief Returns the least upper bound value of the distribution.
   1505        */
   1506       result_type
   1507       max() const
   1508       { return std::numeric_limits<result_type>::max(); }
   1509 
   1510       /**
   1511        * @brief Generating functions.
   1512        */
   1513       template<typename _UniformRandomNumberGenerator>
   1514 	result_type
   1515 	operator()(_UniformRandomNumberGenerator& __urng)
   1516 	{
   1517 	  return this->mu() * std::pow(this->_M_ud(__urng),
   1518 				       -result_type(1) / this->alpha());
   1519 	}
   1520 
   1521       template<typename _UniformRandomNumberGenerator>
   1522 	result_type
   1523 	operator()(_UniformRandomNumberGenerator& __urng,
   1524 		   const param_type& __p)
   1525 	{
   1526 	  return __p.mu() * std::pow(this->_M_ud(__urng),
   1527 					   -result_type(1) / __p.alpha());
   1528 	}
   1529 
   1530       template<typename _ForwardIterator,
   1531 	       typename _UniformRandomNumberGenerator>
   1532 	void
   1533 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1534 		   _UniformRandomNumberGenerator& __urng)
   1535 	{ this->__generate(__f, __t, __urng, _M_param); }
   1536 
   1537       template<typename _ForwardIterator,
   1538 	       typename _UniformRandomNumberGenerator>
   1539 	void
   1540 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1541 		   _UniformRandomNumberGenerator& __urng,
   1542 		   const param_type& __p)
   1543 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1544 
   1545       template<typename _UniformRandomNumberGenerator>
   1546 	void
   1547 	__generate(result_type* __f, result_type* __t,
   1548 		   _UniformRandomNumberGenerator& __urng,
   1549 		   const param_type& __p)
   1550 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1551 
   1552       /**
   1553        * @brief Return true if two Pareto distributions have
   1554        *        the same parameters and the sequences that would
   1555        *        be generated are equal.
   1556        */
   1557       friend bool
   1558       operator==(const pareto_distribution& __d1,
   1559 		 const pareto_distribution& __d2)
   1560       { return (__d1._M_param == __d2._M_param
   1561 		&& __d1._M_ud == __d2._M_ud); }
   1562 
   1563       /**
   1564        * @brief Inserts a %pareto_distribution random number distribution
   1565        * @p __x into the output stream @p __os.
   1566        *
   1567        * @param __os An output stream.
   1568        * @param __x  A %pareto_distribution random number distribution.
   1569        *
   1570        * @returns The output stream with the state of @p __x inserted or in
   1571        * an error state.
   1572        */
   1573       template<typename _RealType1, typename _CharT, typename _Traits>
   1574 	friend std::basic_ostream<_CharT, _Traits>&
   1575 	operator<<(std::basic_ostream<_CharT, _Traits>&,
   1576 		   const pareto_distribution<_RealType1>&);
   1577 
   1578       /**
   1579        * @brief Extracts a %pareto_distribution random number distribution
   1580        * @p __x from the input stream @p __is.
   1581        *
   1582        * @param __is An input stream.
   1583        * @param __x A %pareto_distribution random number
   1584        *            generator engine.
   1585        *
   1586        * @returns The input stream with @p __x extracted or in an error state.
   1587        */
   1588       template<typename _RealType1, typename _CharT, typename _Traits>
   1589 	friend std::basic_istream<_CharT, _Traits>&
   1590 	operator>>(std::basic_istream<_CharT, _Traits>&,
   1591 		   pareto_distribution<_RealType1>&);
   1592 
   1593     private:
   1594       template<typename _ForwardIterator,
   1595 	       typename _UniformRandomNumberGenerator>
   1596 	void
   1597 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   1598 			_UniformRandomNumberGenerator& __urng,
   1599 			const param_type& __p);
   1600 
   1601       param_type _M_param;
   1602 
   1603       std::uniform_real_distribution<result_type> _M_ud;
   1604     };
   1605 
   1606   /**
   1607    * @brief Return true if two Pareto distributions are not equal.
   1608    */
   1609   template<typename _RealType>
   1610     inline bool
   1611     operator!=(const pareto_distribution<_RealType>& __d1,
   1612 	       const pareto_distribution<_RealType>& __d2)
   1613     { return !(__d1 == __d2); }
   1614 
   1615 
   1616   /**
   1617    * @brief A K continuous distribution for random numbers.
   1618    *
   1619    * The formula for the K probability density function is
   1620    * @f[
   1621    *     p(x|\lambda, \mu, \nu) = \frac{2}{x}
   1622    *             \left(\frac{\lambda\nu x}{\mu}\right)^{\frac{\lambda + \nu}{2}}
   1623    *             \frac{1}{\Gamma(\lambda)\Gamma(\nu)}
   1624    *             K_{\nu - \lambda}\left(2\sqrt{\frac{\lambda\nu x}{\mu}}\right)
   1625    * @f]
   1626    * where @f$I_0(z)@f$ is the modified Bessel function of the second kind
   1627    * of order @f$\nu - \lambda@f$ and @f$\lambda > 0@f$, @f$\mu > 0@f$
   1628    * and @f$\nu > 0@f$.
   1629    *
   1630    * <table border=1 cellpadding=10 cellspacing=0>
   1631    * <caption align=top>Distribution Statistics</caption>
   1632    * <tr><td>Mean</td><td>@f$\mu@f$</td></tr>
   1633    * <tr><td>Variance</td><td>@f$\mu^2\frac{\lambda + \nu + 1}{\lambda\nu}@f$</td></tr>
   1634    * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
   1635    * </table>
   1636    */
   1637   template<typename _RealType = double>
   1638     class
   1639     k_distribution
   1640     {
   1641       static_assert(std::is_floating_point<_RealType>::value,
   1642 		    "template argument not a floating point type");
   1643 
   1644     public:
   1645       /** The type of the range of the distribution. */
   1646       typedef _RealType result_type;
   1647       /** Parameter type. */
   1648       struct param_type
   1649       {
   1650 	typedef k_distribution<result_type> distribution_type;
   1651 
   1652 	param_type(result_type __lambda_val = result_type(1),
   1653 		   result_type __mu_val = result_type(1),
   1654 		   result_type __nu_val = result_type(1))
   1655 	: _M_lambda(__lambda_val), _M_mu(__mu_val), _M_nu(__nu_val)
   1656 	{
   1657 	  _GLIBCXX_DEBUG_ASSERT(_M_lambda > result_type(0));
   1658 	  _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0));
   1659 	  _GLIBCXX_DEBUG_ASSERT(_M_nu > result_type(0));
   1660 	}
   1661 
   1662 	result_type
   1663 	lambda() const
   1664 	{ return _M_lambda; }
   1665 
   1666 	result_type
   1667 	mu() const
   1668 	{ return _M_mu; }
   1669 
   1670 	result_type
   1671 	nu() const
   1672 	{ return _M_nu; }
   1673 
   1674 	friend bool
   1675 	operator==(const param_type& __p1, const param_type& __p2)
   1676 	{ return __p1._M_lambda == __p2._M_lambda
   1677 	      && __p1._M_mu == __p2._M_mu
   1678 	      && __p1._M_nu == __p2._M_nu; }
   1679 
   1680       private:
   1681 	void _M_initialize();
   1682 
   1683 	result_type _M_lambda;
   1684 	result_type _M_mu;
   1685 	result_type _M_nu;
   1686       };
   1687 
   1688       /**
   1689        * @brief Constructors.
   1690        */
   1691       explicit
   1692       k_distribution(result_type __lambda_val = result_type(1),
   1693 		     result_type __mu_val = result_type(1),
   1694 		     result_type __nu_val = result_type(1))
   1695       : _M_param(__lambda_val, __mu_val, __nu_val),
   1696 	_M_gd1(__lambda_val, result_type(1) / __lambda_val),
   1697 	_M_gd2(__nu_val, __mu_val / __nu_val)
   1698       { }
   1699 
   1700       explicit
   1701       k_distribution(const param_type& __p)
   1702       : _M_param(__p),
   1703 	_M_gd1(__p.lambda(), result_type(1) / __p.lambda()),
   1704 	_M_gd2(__p.nu(), __p.mu() / __p.nu())
   1705       { }
   1706 
   1707       /**
   1708        * @brief Resets the distribution state.
   1709        */
   1710       void
   1711       reset()
   1712       {
   1713 	_M_gd1.reset();
   1714 	_M_gd2.reset();
   1715       }
   1716 
   1717       /**
   1718        * @brief Return the parameters of the distribution.
   1719        */
   1720       result_type
   1721       lambda() const
   1722       { return _M_param.lambda(); }
   1723 
   1724       result_type
   1725       mu() const
   1726       { return _M_param.mu(); }
   1727 
   1728       result_type
   1729       nu() const
   1730       { return _M_param.nu(); }
   1731 
   1732       /**
   1733        * @brief Returns the parameter set of the distribution.
   1734        */
   1735       param_type
   1736       param() const
   1737       { return _M_param; }
   1738 
   1739       /**
   1740        * @brief Sets the parameter set of the distribution.
   1741        * @param __param The new parameter set of the distribution.
   1742        */
   1743       void
   1744       param(const param_type& __param)
   1745       { _M_param = __param; }
   1746 
   1747       /**
   1748        * @brief Returns the greatest lower bound value of the distribution.
   1749        */
   1750       result_type
   1751       min() const
   1752       { return result_type(0); }
   1753 
   1754       /**
   1755        * @brief Returns the least upper bound value of the distribution.
   1756        */
   1757       result_type
   1758       max() const
   1759       { return std::numeric_limits<result_type>::max(); }
   1760 
   1761       /**
   1762        * @brief Generating functions.
   1763        */
   1764       template<typename _UniformRandomNumberGenerator>
   1765 	result_type
   1766 	operator()(_UniformRandomNumberGenerator&);
   1767 
   1768       template<typename _UniformRandomNumberGenerator>
   1769 	result_type
   1770 	operator()(_UniformRandomNumberGenerator&, const param_type&);
   1771 
   1772       template<typename _ForwardIterator,
   1773 	       typename _UniformRandomNumberGenerator>
   1774 	void
   1775 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1776 		   _UniformRandomNumberGenerator& __urng)
   1777 	{ this->__generate(__f, __t, __urng, _M_param); }
   1778 
   1779       template<typename _ForwardIterator,
   1780 	       typename _UniformRandomNumberGenerator>
   1781 	void
   1782 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   1783 		   _UniformRandomNumberGenerator& __urng,
   1784 		   const param_type& __p)
   1785 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1786 
   1787       template<typename _UniformRandomNumberGenerator>
   1788 	void
   1789 	__generate(result_type* __f, result_type* __t,
   1790 		   _UniformRandomNumberGenerator& __urng,
   1791 		   const param_type& __p)
   1792 	{ this->__generate_impl(__f, __t, __urng, __p); }
   1793 
   1794       /**
   1795        * @brief Return true if two K distributions have
   1796        *        the same parameters and the sequences that would
   1797        *        be generated are equal.
   1798        */
   1799       friend bool
   1800       operator==(const k_distribution& __d1,
   1801 		 const k_distribution& __d2)
   1802       { return (__d1._M_param == __d2._M_param
   1803 		&& __d1._M_gd1 == __d2._M_gd1
   1804 		&& __d1._M_gd2 == __d2._M_gd2); }
   1805 
   1806       /**
   1807        * @brief Inserts a %k_distribution random number distribution
   1808        * @p __x into the output stream @p __os.
   1809        *
   1810        * @param __os An output stream.
   1811        * @param __x  A %k_distribution random number distribution.
   1812        *
   1813        * @returns The output stream with the state of @p __x inserted or in
   1814        * an error state.
   1815        */
   1816       template<typename _RealType1, typename _CharT, typename _Traits>
   1817 	friend std::basic_ostream<_CharT, _Traits>&
   1818 	operator<<(std::basic_ostream<_CharT, _Traits>&,
   1819 		   const k_distribution<_RealType1>&);
   1820 
   1821       /**
   1822        * @brief Extracts a %k_distribution random number distribution
   1823        * @p __x from the input stream @p __is.
   1824        *
   1825        * @param __is An input stream.
   1826        * @param __x A %k_distribution random number
   1827        *            generator engine.
   1828        *
   1829        * @returns The input stream with @p __x extracted or in an error state.
   1830        */
   1831       template<typename _RealType1, typename _CharT, typename _Traits>
   1832 	friend std::basic_istream<_CharT, _Traits>&
   1833 	operator>>(std::basic_istream<_CharT, _Traits>&,
   1834 		   k_distribution<_RealType1>&);
   1835 
   1836     private:
   1837       template<typename _ForwardIterator,
   1838 	       typename _UniformRandomNumberGenerator>
   1839 	void
   1840 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   1841 			_UniformRandomNumberGenerator& __urng,
   1842 			const param_type& __p);
   1843 
   1844       param_type _M_param;
   1845 
   1846       std::gamma_distribution<result_type> _M_gd1;
   1847       std::gamma_distribution<result_type> _M_gd2;
   1848     };
   1849 
   1850   /**
   1851    * @brief Return true if two K distributions are not equal.
   1852    */
   1853   template<typename _RealType>
   1854     inline bool
   1855     operator!=(const k_distribution<_RealType>& __d1,
   1856 	       const k_distribution<_RealType>& __d2)
   1857     { return !(__d1 == __d2); }
   1858 
   1859 
   1860   /**
   1861    * @brief An arcsine continuous distribution for random numbers.
   1862    *
   1863    * The formula for the arcsine probability density function is
   1864    * @f[
   1865    *     p(x|a,b) = \frac{1}{\pi \sqrt{(x - a)(b - x)}}
   1866    * @f]
   1867    * where @f$x >= a@f$ and @f$x <= b@f$.
   1868    *
   1869    * <table border=1 cellpadding=10 cellspacing=0>
   1870    * <caption align=top>Distribution Statistics</caption>
   1871    * <tr><td>Mean</td><td>@f$ (a + b) / 2 @f$</td></tr>
   1872    * <tr><td>Variance</td><td>@f$ (b - a)^2 / 8 @f$</td></tr>
   1873    * <tr><td>Range</td><td>@f$[a, b]@f$</td></tr>
   1874    * </table>
   1875    */
   1876   template<typename _RealType = double>
   1877     class
   1878     arcsine_distribution
   1879     {
   1880       static_assert(std::is_floating_point<_RealType>::value,
   1881 		    "template argument not a floating point type");
   1882 
   1883     public:
   1884       /** The type of the range of the distribution. */
   1885       typedef _RealType result_type;
   1886       /** Parameter type. */
   1887       struct param_type
   1888       {
   1889 	typedef arcsine_distribution<result_type> distribution_type;
   1890 
   1891 	param_type(result_type __a = result_type(0),
   1892 		   result_type __b = result_type(1))
   1893 	: _M_a(__a), _M_b(__b)
   1894 	{
   1895 	  _GLIBCXX_DEBUG_ASSERT(_M_a <= _M_b);
   1896 	}
   1897 
   1898 	result_type
   1899 	a() const
   1900 	{ return _M_a; }
   1901 
   1902 	result_type
   1903 	b() const
   1904 	{ return _M_b; }
   1905 
   1906 	friend bool
   1907 	operator==(const param_type& __p1, const param_type& __p2)
   1908 	{ return __p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b; }
   1909 
   1910       private:
   1911 	void _M_initialize();
   1912 
   1913 	result_type _M_a;
   1914 	result_type _M_b;
   1915       };
   1916 
   1917       /**
   1918        * @brief Constructors.
   1919        */
   1920       explicit
   1921       arcsine_distribution(result_type __a = result_type(0),
   1922 			   result_type __b = result_type(1))
   1923       : _M_param(__a, __b),
   1924 	_M_ud(-1.5707963267948966192313216916397514L,
   1925 	      +1.5707963267948966192313216916397514L)
   1926       { }
   1927 
   1928       explicit
   1929       arcsine_distribution(const param_type& __p)
   1930       : _M_param(__p),
   1931 	_M_ud(-1.5707963267948966192313216916397514L,
   1932 	      +1.5707963267948966192313216916397514L)
   1933       { }
   1934 
   1935       /**
   1936        * @brief Resets the distribution state.
   1937        */
   1938       void
   1939       reset()
   1940       { _M_ud.reset(); }
   1941 
   1942       /**
   1943        * @brief Return the parameters of the distribution.
   1944        */
   1945       result_type
   1946       a() const
   1947       { return _M_param.a(); }
   1948 
   1949       result_type
   1950       b() const
   1951       { return _M_param.b(); }
   1952 
   1953       /**
   1954        * @brief Returns the parameter set of the distribution.
   1955        */
   1956       param_type
   1957       param() const
   1958       { return _M_param; }
   1959 
   1960       /**
   1961        * @brief Sets the parameter set of the distribution.
   1962        * @param __param The new parameter set of the distribution.
   1963        */
   1964       void
   1965       param(const param_type& __param)
   1966       { _M_param = __param; }
   1967 
   1968       /**
   1969        * @brief Returns the greatest lower bound value of the distribution.
   1970        */
   1971       result_type
   1972       min() const
   1973       { return this->a(); }
   1974 
   1975       /**
   1976        * @brief Returns the least upper bound value of the distribution.
   1977        */
   1978       result_type
   1979       max() const
   1980       { return this->b(); }
   1981 
   1982       /**
   1983        * @brief Generating functions.
   1984        */
   1985       template<typename _UniformRandomNumberGenerator>
   1986 	result_type
   1987 	operator()(_UniformRandomNumberGenerator& __urng)
   1988 	{
   1989 	  result_type __x = std::sin(this->_M_ud(__urng));
   1990 	  return (__x * (this->b() - this->a())
   1991 		  + this->a() + this->b()) / result_type(2);
   1992 	}
   1993 
   1994       template<typename _UniformRandomNumberGenerator>
   1995 	result_type
   1996 	operator()(_UniformRandomNumberGenerator& __urng,
   1997 		   const param_type& __p)
   1998 	{
   1999 	  result_type __x = std::sin(this->_M_ud(__urng));
   2000 	  return (__x * (__p.b() - __p.a())
   2001 		  + __p.a() + __p.b()) / result_type(2);
   2002 	}
   2003 
   2004       template<typename _ForwardIterator,
   2005 	       typename _UniformRandomNumberGenerator>
   2006 	void
   2007 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2008 		   _UniformRandomNumberGenerator& __urng)
   2009 	{ this->__generate(__f, __t, __urng, _M_param); }
   2010 
   2011       template<typename _ForwardIterator,
   2012 	       typename _UniformRandomNumberGenerator>
   2013 	void
   2014 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2015 		   _UniformRandomNumberGenerator& __urng,
   2016 		   const param_type& __p)
   2017 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2018 
   2019       template<typename _UniformRandomNumberGenerator>
   2020 	void
   2021 	__generate(result_type* __f, result_type* __t,
   2022 		   _UniformRandomNumberGenerator& __urng,
   2023 		   const param_type& __p)
   2024 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2025 
   2026       /**
   2027        * @brief Return true if two arcsine distributions have
   2028        *        the same parameters and the sequences that would
   2029        *        be generated are equal.
   2030        */
   2031       friend bool
   2032       operator==(const arcsine_distribution& __d1,
   2033 		 const arcsine_distribution& __d2)
   2034       { return (__d1._M_param == __d2._M_param
   2035 		&& __d1._M_ud == __d2._M_ud); }
   2036 
   2037       /**
   2038        * @brief Inserts a %arcsine_distribution random number distribution
   2039        * @p __x into the output stream @p __os.
   2040        *
   2041        * @param __os An output stream.
   2042        * @param __x  A %arcsine_distribution random number distribution.
   2043        *
   2044        * @returns The output stream with the state of @p __x inserted or in
   2045        * an error state.
   2046        */
   2047       template<typename _RealType1, typename _CharT, typename _Traits>
   2048 	friend std::basic_ostream<_CharT, _Traits>&
   2049 	operator<<(std::basic_ostream<_CharT, _Traits>&,
   2050 		   const arcsine_distribution<_RealType1>&);
   2051 
   2052       /**
   2053        * @brief Extracts a %arcsine_distribution random number distribution
   2054        * @p __x from the input stream @p __is.
   2055        *
   2056        * @param __is An input stream.
   2057        * @param __x A %arcsine_distribution random number
   2058        *            generator engine.
   2059        *
   2060        * @returns The input stream with @p __x extracted or in an error state.
   2061        */
   2062       template<typename _RealType1, typename _CharT, typename _Traits>
   2063 	friend std::basic_istream<_CharT, _Traits>&
   2064 	operator>>(std::basic_istream<_CharT, _Traits>&,
   2065 		   arcsine_distribution<_RealType1>&);
   2066 
   2067     private:
   2068       template<typename _ForwardIterator,
   2069 	       typename _UniformRandomNumberGenerator>
   2070 	void
   2071 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   2072 			_UniformRandomNumberGenerator& __urng,
   2073 			const param_type& __p);
   2074 
   2075       param_type _M_param;
   2076 
   2077       std::uniform_real_distribution<result_type> _M_ud;
   2078     };
   2079 
   2080   /**
   2081    * @brief Return true if two arcsine distributions are not equal.
   2082    */
   2083   template<typename _RealType>
   2084     inline bool
   2085     operator!=(const arcsine_distribution<_RealType>& __d1,
   2086 	       const arcsine_distribution<_RealType>& __d2)
   2087     { return !(__d1 == __d2); }
   2088 
   2089 
   2090   /**
   2091    * @brief A Hoyt continuous distribution for random numbers.
   2092    *
   2093    * The formula for the Hoyt probability density function is
   2094    * @f[
   2095    *     p(x|q,\omega) = \frac{(1 + q^2)x}{q\omega}
   2096    *                     \exp\left(-\frac{(1 + q^2)^2 x^2}{4 q^2 \omega}\right)
   2097    *                       I_0\left(\frac{(1 - q^4) x^2}{4 q^2 \omega}\right)
   2098    * @f]
   2099    * where @f$I_0(z)@f$ is the modified Bessel function of the first kind
   2100    * of order 0 and @f$0 < q < 1@f$.
   2101    *
   2102    * <table border=1 cellpadding=10 cellspacing=0>
   2103    * <caption align=top>Distribution Statistics</caption>
   2104    * <tr><td>Mean</td><td>@f$ \sqrt{\frac{2}{\pi}} \sqrt{\frac{\omega}{1 + q^2}}
   2105    *                       E(1 - q^2) @f$</td></tr>
   2106    * <tr><td>Variance</td><td>@f$ \omega \left(1 - \frac{2E^2(1 - q^2)}
   2107    *                                      {\pi (1 + q^2)}\right) @f$</td></tr>
   2108    * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
   2109    * </table>
   2110    * where @f$E(x)@f$ is the elliptic function of the second kind.
   2111    */
   2112   template<typename _RealType = double>
   2113     class
   2114     hoyt_distribution
   2115     {
   2116       static_assert(std::is_floating_point<_RealType>::value,
   2117 		    "template argument not a floating point type");
   2118 
   2119     public:
   2120       /** The type of the range of the distribution. */
   2121       typedef _RealType result_type;
   2122       /** Parameter type. */
   2123       struct param_type
   2124       {
   2125 	typedef hoyt_distribution<result_type> distribution_type;
   2126 
   2127 	param_type(result_type __q = result_type(0.5L),
   2128 		   result_type __omega = result_type(1))
   2129 	: _M_q(__q), _M_omega(__omega)
   2130 	{
   2131 	  _GLIBCXX_DEBUG_ASSERT(_M_q > result_type(0));
   2132 	  _GLIBCXX_DEBUG_ASSERT(_M_q < result_type(1));
   2133 	}
   2134 
   2135 	result_type
   2136 	q() const
   2137 	{ return _M_q; }
   2138 
   2139 	result_type
   2140 	omega() const
   2141 	{ return _M_omega; }
   2142 
   2143 	friend bool
   2144 	operator==(const param_type& __p1, const param_type& __p2)
   2145 	{ return __p1._M_q == __p2._M_q
   2146 	      && __p1._M_omega == __p2._M_omega; }
   2147 
   2148       private:
   2149 	void _M_initialize();
   2150 
   2151 	result_type _M_q;
   2152 	result_type _M_omega;
   2153       };
   2154 
   2155       /**
   2156        * @brief Constructors.
   2157        */
   2158       explicit
   2159       hoyt_distribution(result_type __q = result_type(0.5L),
   2160 			result_type __omega = result_type(1))
   2161       : _M_param(__q, __omega),
   2162 	_M_ad(result_type(0.5L) * (result_type(1) + __q * __q),
   2163 	      result_type(0.5L) * (result_type(1) + __q * __q)
   2164 				/ (__q * __q)),
   2165 	_M_ed(result_type(1))
   2166       { }
   2167 
   2168       explicit
   2169       hoyt_distribution(const param_type& __p)
   2170       : _M_param(__p),
   2171 	_M_ad(result_type(0.5L) * (result_type(1) + __p.q() * __p.q()),
   2172 	      result_type(0.5L) * (result_type(1) + __p.q() * __p.q())
   2173 				/ (__p.q() * __p.q())),
   2174 	_M_ed(result_type(1))
   2175       { }
   2176 
   2177       /**
   2178        * @brief Resets the distribution state.
   2179        */
   2180       void
   2181       reset()
   2182       {
   2183 	_M_ad.reset();
   2184 	_M_ed.reset();
   2185       }
   2186 
   2187       /**
   2188        * @brief Return the parameters of the distribution.
   2189        */
   2190       result_type
   2191       q() const
   2192       { return _M_param.q(); }
   2193 
   2194       result_type
   2195       omega() const
   2196       { return _M_param.omega(); }
   2197 
   2198       /**
   2199        * @brief Returns the parameter set of the distribution.
   2200        */
   2201       param_type
   2202       param() const
   2203       { return _M_param; }
   2204 
   2205       /**
   2206        * @brief Sets the parameter set of the distribution.
   2207        * @param __param The new parameter set of the distribution.
   2208        */
   2209       void
   2210       param(const param_type& __param)
   2211       { _M_param = __param; }
   2212 
   2213       /**
   2214        * @brief Returns the greatest lower bound value of the distribution.
   2215        */
   2216       result_type
   2217       min() const
   2218       { return result_type(0); }
   2219 
   2220       /**
   2221        * @brief Returns the least upper bound value of the distribution.
   2222        */
   2223       result_type
   2224       max() const
   2225       { return std::numeric_limits<result_type>::max(); }
   2226 
   2227       /**
   2228        * @brief Generating functions.
   2229        */
   2230       template<typename _UniformRandomNumberGenerator>
   2231 	result_type
   2232 	operator()(_UniformRandomNumberGenerator& __urng);
   2233 
   2234       template<typename _UniformRandomNumberGenerator>
   2235 	result_type
   2236 	operator()(_UniformRandomNumberGenerator& __urng,
   2237 		   const param_type& __p);
   2238 
   2239       template<typename _ForwardIterator,
   2240 	       typename _UniformRandomNumberGenerator>
   2241 	void
   2242 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2243 		   _UniformRandomNumberGenerator& __urng)
   2244 	{ this->__generate(__f, __t, __urng, _M_param); }
   2245 
   2246       template<typename _ForwardIterator,
   2247 	       typename _UniformRandomNumberGenerator>
   2248 	void
   2249 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2250 		   _UniformRandomNumberGenerator& __urng,
   2251 		   const param_type& __p)
   2252 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2253 
   2254       template<typename _UniformRandomNumberGenerator>
   2255 	void
   2256 	__generate(result_type* __f, result_type* __t,
   2257 		   _UniformRandomNumberGenerator& __urng,
   2258 		   const param_type& __p)
   2259 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2260 
   2261       /**
   2262        * @brief Return true if two Hoyt distributions have
   2263        *        the same parameters and the sequences that would
   2264        *        be generated are equal.
   2265        */
   2266       friend bool
   2267       operator==(const hoyt_distribution& __d1,
   2268 		 const hoyt_distribution& __d2)
   2269       { return (__d1._M_param == __d2._M_param
   2270 		&& __d1._M_ad == __d2._M_ad
   2271 		&& __d1._M_ed == __d2._M_ed); }
   2272 
   2273       /**
   2274        * @brief Inserts a %hoyt_distribution random number distribution
   2275        * @p __x into the output stream @p __os.
   2276        *
   2277        * @param __os An output stream.
   2278        * @param __x  A %hoyt_distribution random number distribution.
   2279        *
   2280        * @returns The output stream with the state of @p __x inserted or in
   2281        * an error state.
   2282        */
   2283       template<typename _RealType1, typename _CharT, typename _Traits>
   2284 	friend std::basic_ostream<_CharT, _Traits>&
   2285 	operator<<(std::basic_ostream<_CharT, _Traits>&,
   2286 		   const hoyt_distribution<_RealType1>&);
   2287 
   2288       /**
   2289        * @brief Extracts a %hoyt_distribution random number distribution
   2290        * @p __x from the input stream @p __is.
   2291        *
   2292        * @param __is An input stream.
   2293        * @param __x A %hoyt_distribution random number
   2294        *            generator engine.
   2295        *
   2296        * @returns The input stream with @p __x extracted or in an error state.
   2297        */
   2298       template<typename _RealType1, typename _CharT, typename _Traits>
   2299 	friend std::basic_istream<_CharT, _Traits>&
   2300 	operator>>(std::basic_istream<_CharT, _Traits>&,
   2301 		   hoyt_distribution<_RealType1>&);
   2302 
   2303     private:
   2304       template<typename _ForwardIterator,
   2305 	       typename _UniformRandomNumberGenerator>
   2306 	void
   2307 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   2308 			_UniformRandomNumberGenerator& __urng,
   2309 			const param_type& __p);
   2310 
   2311       param_type _M_param;
   2312 
   2313       __gnu_cxx::arcsine_distribution<result_type> _M_ad;
   2314       std::exponential_distribution<result_type> _M_ed;
   2315     };
   2316 
   2317   /**
   2318    * @brief Return true if two Hoyt distributions are not equal.
   2319    */
   2320   template<typename _RealType>
   2321     inline bool
   2322     operator!=(const hoyt_distribution<_RealType>& __d1,
   2323 	       const hoyt_distribution<_RealType>& __d2)
   2324     { return !(__d1 == __d2); }
   2325 
   2326 
   2327   /**
   2328    * @brief A triangular distribution for random numbers.
   2329    *
   2330    * The formula for the triangular probability density function is
   2331    * @f[
   2332    *                  / 0                          for x < a
   2333    *     p(x|a,b,c) = | \frac{2(x-a)}{(c-a)(b-a)}  for a <= x <= b
   2334    *                  | \frac{2(c-x)}{(c-a)(c-b)}  for b < x <= c
   2335    *                  \ 0                          for c < x
   2336    * @f]
   2337    *
   2338    * <table border=1 cellpadding=10 cellspacing=0>
   2339    * <caption align=top>Distribution Statistics</caption>
   2340    * <tr><td>Mean</td><td>@f$ \frac{a+b+c}{2} @f$</td></tr>
   2341    * <tr><td>Variance</td><td>@f$ \frac{a^2+b^2+c^2-ab-ac-bc}
   2342    *                                   {18}@f$</td></tr>
   2343    * <tr><td>Range</td><td>@f$[a, c]@f$</td></tr>
   2344    * </table>
   2345    */
   2346   template<typename _RealType = double>
   2347     class triangular_distribution
   2348     {
   2349       static_assert(std::is_floating_point<_RealType>::value,
   2350 		    "template argument not a floating point type");
   2351 
   2352     public:
   2353       /** The type of the range of the distribution. */
   2354       typedef _RealType result_type;
   2355       /** Parameter type. */
   2356       struct param_type
   2357       {
   2358 	friend class triangular_distribution<_RealType>;
   2359 
   2360 	explicit
   2361 	param_type(_RealType __a = _RealType(0),
   2362 		   _RealType __b = _RealType(0.5),
   2363 		   _RealType __c = _RealType(1))
   2364 	: _M_a(__a), _M_b(__b), _M_c(__c)
   2365 	{
   2366 	  _GLIBCXX_DEBUG_ASSERT(_M_a <= _M_b);
   2367 	  _GLIBCXX_DEBUG_ASSERT(_M_b <= _M_c);
   2368 	  _GLIBCXX_DEBUG_ASSERT(_M_a < _M_c);
   2369 
   2370 	  _M_r_ab = (_M_b - _M_a) / (_M_c - _M_a);
   2371 	  _M_f_ab_ac = (_M_b - _M_a) * (_M_c - _M_a);
   2372 	  _M_f_bc_ac = (_M_c - _M_b) * (_M_c - _M_a);
   2373 	}
   2374 
   2375 	_RealType
   2376 	a() const
   2377 	{ return _M_a; }
   2378 
   2379 	_RealType
   2380 	b() const
   2381 	{ return _M_b; }
   2382 
   2383 	_RealType
   2384 	c() const
   2385 	{ return _M_c; }
   2386 
   2387 	friend bool
   2388 	operator==(const param_type& __p1, const param_type& __p2)
   2389 	{ return (__p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b
   2390 		  && __p1._M_c == __p2._M_c); }
   2391 
   2392       private:
   2393 
   2394 	_RealType _M_a;
   2395 	_RealType _M_b;
   2396 	_RealType _M_c;
   2397 	_RealType _M_r_ab;
   2398 	_RealType _M_f_ab_ac;
   2399 	_RealType _M_f_bc_ac;
   2400       };
   2401 
   2402       /**
   2403        * @brief Constructs a triangle distribution with parameters
   2404        * @f$ a @f$, @f$ b @f$ and @f$ c @f$.
   2405        */
   2406       explicit
   2407       triangular_distribution(result_type __a = result_type(0),
   2408 			      result_type __b = result_type(0.5),
   2409 			      result_type __c = result_type(1))
   2410       : _M_param(__a, __b, __c)
   2411       { }
   2412 
   2413       explicit
   2414       triangular_distribution(const param_type& __p)
   2415       : _M_param(__p)
   2416       { }
   2417 
   2418       /**
   2419        * @brief Resets the distribution state.
   2420        */
   2421       void
   2422       reset()
   2423       { }
   2424 
   2425       /**
   2426        * @brief Returns the @f$ a @f$ of the distribution.
   2427        */
   2428       result_type
   2429       a() const
   2430       { return _M_param.a(); }
   2431 
   2432       /**
   2433        * @brief Returns the @f$ b @f$ of the distribution.
   2434        */
   2435       result_type
   2436       b() const
   2437       { return _M_param.b(); }
   2438 
   2439       /**
   2440        * @brief Returns the @f$ c @f$ of the distribution.
   2441        */
   2442       result_type
   2443       c() const
   2444       { return _M_param.c(); }
   2445 
   2446       /**
   2447        * @brief Returns the parameter set of the distribution.
   2448        */
   2449       param_type
   2450       param() const
   2451       { return _M_param; }
   2452 
   2453       /**
   2454        * @brief Sets the parameter set of the distribution.
   2455        * @param __param The new parameter set of the distribution.
   2456        */
   2457       void
   2458       param(const param_type& __param)
   2459       { _M_param = __param; }
   2460 
   2461       /**
   2462        * @brief Returns the greatest lower bound value of the distribution.
   2463        */
   2464       result_type
   2465       min() const
   2466       { return _M_param._M_a; }
   2467 
   2468       /**
   2469        * @brief Returns the least upper bound value of the distribution.
   2470        */
   2471       result_type
   2472       max() const
   2473       { return _M_param._M_c; }
   2474 
   2475       /**
   2476        * @brief Generating functions.
   2477        */
   2478       template<typename _UniformRandomNumberGenerator>
   2479 	result_type
   2480 	operator()(_UniformRandomNumberGenerator& __urng)
   2481 	{ return this->operator()(__urng, _M_param); }
   2482 
   2483       template<typename _UniformRandomNumberGenerator>
   2484 	result_type
   2485 	operator()(_UniformRandomNumberGenerator& __urng,
   2486 		   const param_type& __p)
   2487 	{
   2488 	  std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
   2489 	    __aurng(__urng);
   2490 	  result_type __rnd = __aurng();
   2491 	  if (__rnd <= __p._M_r_ab)
   2492 	    return __p.a() + std::sqrt(__rnd * __p._M_f_ab_ac);
   2493 	  else
   2494 	    return __p.c() - std::sqrt((result_type(1) - __rnd)
   2495 				       * __p._M_f_bc_ac);
   2496 	}
   2497 
   2498       template<typename _ForwardIterator,
   2499 	       typename _UniformRandomNumberGenerator>
   2500 	void
   2501 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2502 		   _UniformRandomNumberGenerator& __urng)
   2503 	{ this->__generate(__f, __t, __urng, _M_param); }
   2504 
   2505       template<typename _ForwardIterator,
   2506 	       typename _UniformRandomNumberGenerator>
   2507 	void
   2508 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2509 		   _UniformRandomNumberGenerator& __urng,
   2510 		   const param_type& __p)
   2511 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2512 
   2513       template<typename _UniformRandomNumberGenerator>
   2514 	void
   2515 	__generate(result_type* __f, result_type* __t,
   2516 		   _UniformRandomNumberGenerator& __urng,
   2517 		   const param_type& __p)
   2518 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2519 
   2520       /**
   2521        * @brief Return true if two triangle distributions have the same
   2522        *        parameters and the sequences that would be generated
   2523        *        are equal.
   2524        */
   2525       friend bool
   2526       operator==(const triangular_distribution& __d1,
   2527 		 const triangular_distribution& __d2)
   2528       { return __d1._M_param == __d2._M_param; }
   2529 
   2530       /**
   2531        * @brief Inserts a %triangular_distribution random number distribution
   2532        * @p __x into the output stream @p __os.
   2533        *
   2534        * @param __os An output stream.
   2535        * @param __x  A %triangular_distribution random number distribution.
   2536        *
   2537        * @returns The output stream with the state of @p __x inserted or in
   2538        * an error state.
   2539        */
   2540       template<typename _RealType1, typename _CharT, typename _Traits>
   2541 	friend std::basic_ostream<_CharT, _Traits>&
   2542 	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   2543 		   const __gnu_cxx::triangular_distribution<_RealType1>& __x);
   2544 
   2545       /**
   2546        * @brief Extracts a %triangular_distribution random number distribution
   2547        * @p __x from the input stream @p __is.
   2548        *
   2549        * @param __is An input stream.
   2550        * @param __x  A %triangular_distribution random number generator engine.
   2551        *
   2552        * @returns The input stream with @p __x extracted or in an error state.
   2553        */
   2554       template<typename _RealType1, typename _CharT, typename _Traits>
   2555 	friend std::basic_istream<_CharT, _Traits>&
   2556 	operator>>(std::basic_istream<_CharT, _Traits>& __is,
   2557 		   __gnu_cxx::triangular_distribution<_RealType1>& __x);
   2558 
   2559     private:
   2560       template<typename _ForwardIterator,
   2561 	       typename _UniformRandomNumberGenerator>
   2562 	void
   2563 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   2564 			_UniformRandomNumberGenerator& __urng,
   2565 			const param_type& __p);
   2566 
   2567       param_type _M_param;
   2568     };
   2569 
   2570   /**
   2571    * @brief Return true if two triangle distributions are different.
   2572    */
   2573   template<typename _RealType>
   2574     inline bool
   2575     operator!=(const __gnu_cxx::triangular_distribution<_RealType>& __d1,
   2576 	       const __gnu_cxx::triangular_distribution<_RealType>& __d2)
   2577    { return !(__d1 == __d2); }
   2578 
   2579 
   2580   /**
   2581    * @brief A von Mises distribution for random numbers.
   2582    *
   2583    * The formula for the von Mises probability density function is
   2584    * @f[
   2585    *     p(x|\mu,\kappa) = \frac{e^{\kappa \cos(x-\mu)}}
   2586    *                            {2\pi I_0(\kappa)}
   2587    * @f]
   2588    *
   2589    * The generating functions use the method according to:
   2590    *
   2591    * D. J. Best and N. I. Fisher, 1979. "Efficient Simulation of the
   2592    * von Mises Distribution", Journal of the Royal Statistical Society.
   2593    * Series C (Applied Statistics), Vol. 28, No. 2, pp. 152-157.
   2594    *
   2595    * <table border=1 cellpadding=10 cellspacing=0>
   2596    * <caption align=top>Distribution Statistics</caption>
   2597    * <tr><td>Mean</td><td>@f$ \mu @f$</td></tr>
   2598    * <tr><td>Variance</td><td>@f$ 1-I_1(\kappa)/I_0(\kappa) @f$</td></tr>
   2599    * <tr><td>Range</td><td>@f$[-\pi, \pi]@f$</td></tr>
   2600    * </table>
   2601    */
   2602   template<typename _RealType = double>
   2603     class von_mises_distribution
   2604     {
   2605       static_assert(std::is_floating_point<_RealType>::value,
   2606 		    "template argument not a floating point type");
   2607 
   2608     public:
   2609       /** The type of the range of the distribution. */
   2610       typedef _RealType result_type;
   2611       /** Parameter type. */
   2612       struct param_type
   2613       {
   2614 	friend class von_mises_distribution<_RealType>;
   2615 
   2616 	explicit
   2617 	param_type(_RealType __mu = _RealType(0),
   2618 		   _RealType __kappa = _RealType(1))
   2619 	: _M_mu(__mu), _M_kappa(__kappa)
   2620 	{
   2621 	  const _RealType __pi = __gnu_cxx::__math_constants<_RealType>::__pi;
   2622 	  _GLIBCXX_DEBUG_ASSERT(_M_mu >= -__pi && _M_mu <= __pi);
   2623 	  _GLIBCXX_DEBUG_ASSERT(_M_kappa >= _RealType(0));
   2624 
   2625 	  auto __tau = std::sqrt(_RealType(4) * _M_kappa * _M_kappa
   2626 				 + _RealType(1)) + _RealType(1);
   2627 	  auto __rho = ((__tau - std::sqrt(_RealType(2) * __tau))
   2628 			/ (_RealType(2) * _M_kappa));
   2629 	  _M_r = (_RealType(1) + __rho * __rho) / (_RealType(2) * __rho);
   2630 	}
   2631 
   2632 	_RealType
   2633 	mu() const
   2634 	{ return _M_mu; }
   2635 
   2636 	_RealType
   2637 	kappa() const
   2638 	{ return _M_kappa; }
   2639 
   2640 	friend bool
   2641 	operator==(const param_type& __p1, const param_type& __p2)
   2642 	{ return (__p1._M_mu == __p2._M_mu
   2643 		  && __p1._M_kappa == __p2._M_kappa); }
   2644 
   2645       private:
   2646 	_RealType _M_mu;
   2647 	_RealType _M_kappa;
   2648 	_RealType _M_r;
   2649       };
   2650 
   2651       /**
   2652        * @brief Constructs a von Mises distribution with parameters
   2653        * @f$\mu@f$ and @f$\kappa@f$.
   2654        */
   2655       explicit
   2656       von_mises_distribution(result_type __mu = result_type(0),
   2657 			     result_type __kappa = result_type(1))
   2658 	: _M_param(__mu, __kappa)
   2659       { }
   2660 
   2661       explicit
   2662       von_mises_distribution(const param_type& __p)
   2663       : _M_param(__p)
   2664       { }
   2665 
   2666       /**
   2667        * @brief Resets the distribution state.
   2668        */
   2669       void
   2670       reset()
   2671       { }
   2672 
   2673       /**
   2674        * @brief Returns the @f$ \mu @f$ of the distribution.
   2675        */
   2676       result_type
   2677       mu() const
   2678       { return _M_param.mu(); }
   2679 
   2680       /**
   2681        * @brief Returns the @f$ \kappa @f$ of the distribution.
   2682        */
   2683       result_type
   2684       kappa() const
   2685       { return _M_param.kappa(); }
   2686 
   2687       /**
   2688        * @brief Returns the parameter set of the distribution.
   2689        */
   2690       param_type
   2691       param() const
   2692       { return _M_param; }
   2693 
   2694       /**
   2695        * @brief Sets the parameter set of the distribution.
   2696        * @param __param The new parameter set of the distribution.
   2697        */
   2698       void
   2699       param(const param_type& __param)
   2700       { _M_param = __param; }
   2701 
   2702       /**
   2703        * @brief Returns the greatest lower bound value of the distribution.
   2704        */
   2705       result_type
   2706       min() const
   2707       {
   2708 	return -__gnu_cxx::__math_constants<result_type>::__pi;
   2709       }
   2710 
   2711       /**
   2712        * @brief Returns the least upper bound value of the distribution.
   2713        */
   2714       result_type
   2715       max() const
   2716       {
   2717 	return __gnu_cxx::__math_constants<result_type>::__pi;
   2718       }
   2719 
   2720       /**
   2721        * @brief Generating functions.
   2722        */
   2723       template<typename _UniformRandomNumberGenerator>
   2724 	result_type
   2725 	operator()(_UniformRandomNumberGenerator& __urng)
   2726 	{ return this->operator()(__urng, _M_param); }
   2727 
   2728       template<typename _UniformRandomNumberGenerator>
   2729 	result_type
   2730 	operator()(_UniformRandomNumberGenerator& __urng,
   2731 		   const param_type& __p)
   2732 	{
   2733 	  const result_type __pi
   2734 	    = __gnu_cxx::__math_constants<result_type>::__pi;
   2735 	  std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
   2736 	    __aurng(__urng);
   2737 
   2738 	  result_type __f;
   2739 	  while (1)
   2740 	    {
   2741 	      result_type __rnd = std::cos(__pi * __aurng());
   2742 	      __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd);
   2743 	      result_type __c = __p._M_kappa * (__p._M_r - __f);
   2744 
   2745 	      result_type __rnd2 = __aurng();
   2746 	      if (__c * (result_type(2) - __c) > __rnd2)
   2747 		break;
   2748 	      if (std::log(__c / __rnd2) >= __c - result_type(1))
   2749 		break;
   2750 	    }
   2751 
   2752 	  result_type __res = std::acos(__f);
   2753 #if _GLIBCXX_USE_C99_MATH_TR1
   2754 	  __res = std::copysign(__res, __aurng() - result_type(0.5));
   2755 #else
   2756 	  if (__aurng() < result_type(0.5))
   2757 	    __res = -__res;
   2758 #endif
   2759 	  __res += __p._M_mu;
   2760 	  if (__res > __pi)
   2761 	    __res -= result_type(2) * __pi;
   2762 	  else if (__res < -__pi)
   2763 	    __res += result_type(2) * __pi;
   2764 	  return __res;
   2765 	}
   2766 
   2767       template<typename _ForwardIterator,
   2768 	       typename _UniformRandomNumberGenerator>
   2769 	void
   2770 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2771 		   _UniformRandomNumberGenerator& __urng)
   2772 	{ this->__generate(__f, __t, __urng, _M_param); }
   2773 
   2774       template<typename _ForwardIterator,
   2775 	       typename _UniformRandomNumberGenerator>
   2776 	void
   2777 	__generate(_ForwardIterator __f, _ForwardIterator __t,
   2778 		   _UniformRandomNumberGenerator& __urng,
   2779 		   const param_type& __p)
   2780 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2781 
   2782       template<typename _UniformRandomNumberGenerator>
   2783 	void
   2784 	__generate(result_type* __f, result_type* __t,
   2785 		   _UniformRandomNumberGenerator& __urng,
   2786 		   const param_type& __p)
   2787 	{ this->__generate_impl(__f, __t, __urng, __p); }
   2788 
   2789       /**
   2790        * @brief Return true if two von Mises distributions have the same
   2791        *        parameters and the sequences that would be generated
   2792        *        are equal.
   2793        */
   2794       friend bool
   2795       operator==(const von_mises_distribution& __d1,
   2796 		 const von_mises_distribution& __d2)
   2797       { return __d1._M_param == __d2._M_param; }
   2798 
   2799       /**
   2800        * @brief Inserts a %von_mises_distribution random number distribution
   2801        * @p __x into the output stream @p __os.
   2802        *
   2803        * @param __os An output stream.
   2804        * @param __x  A %von_mises_distribution random number distribution.
   2805        *
   2806        * @returns The output stream with the state of @p __x inserted or in
   2807        * an error state.
   2808        */
   2809       template<typename _RealType1, typename _CharT, typename _Traits>
   2810 	friend std::basic_ostream<_CharT, _Traits>&
   2811 	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
   2812 		   const __gnu_cxx::von_mises_distribution<_RealType1>& __x);
   2813 
   2814       /**
   2815        * @brief Extracts a %von_mises_distribution random number distribution
   2816        * @p __x from the input stream @p __is.
   2817        *
   2818        * @param __is An input stream.
   2819        * @param __x  A %von_mises_distribution random number generator engine.
   2820        *
   2821        * @returns The input stream with @p __x extracted or in an error state.
   2822        */
   2823       template<typename _RealType1, typename _CharT, typename _Traits>
   2824 	friend std::basic_istream<_CharT, _Traits>&
   2825 	operator>>(std::basic_istream<_CharT, _Traits>& __is,
   2826 		   __gnu_cxx::von_mises_distribution<_RealType1>& __x);
   2827 
   2828     private:
   2829       template<typename _ForwardIterator,
   2830 	       typename _UniformRandomNumberGenerator>
   2831 	void
   2832 	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
   2833 			_UniformRandomNumberGenerator& __urng,
   2834 			const param_type& __p);
   2835 
   2836       param_type _M_param;
   2837     };
   2838 
   2839   /**
   2840    * @brief Return true if two von Mises distributions are different.
   2841    */
   2842   template<typename _RealType>
   2843     inline bool
   2844     operator!=(const __gnu_cxx::von_mises_distribution<_RealType>& __d1,
   2845 	       const __gnu_cxx::von_mises_distribution<_RealType>& __d2)
   2846    { return !(__d1 == __d2); }
   2847 
   2848 _GLIBCXX_END_NAMESPACE_VERSION
   2849 } // namespace __gnu_cxx
   2850 
   2851 #include "ext/opt_random.h"
   2852 #include "random.tcc"
   2853 
   2854 #endif // _GLIBCXX_USE_C99_STDINT_TR1
   2855 
   2856 #endif // C++11
   2857 
   2858 #endif // _EXT_RANDOM
   2859