OpenE2K/gcc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
19ece7ec02
gcc/libstdc++-v3/include/ext/random
Edward Smith-Rowland 19ece7ec02 Implement the Nakagami statistical distribution as an extension. 
From-SVN: r191434
2012-09-18 16:14:09 +00:00

1370 lines
41 KiB
C++
Raw Blame History

// Random number extensions -*- C++ -*-
// Copyright (C) 2012 Free Software Foundation, Inc.
//
// This file is part of the GNU ISO C++ Library. This library is free
// software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the
// Free Software Foundation; either version 3, or (at your option)
// any later version.
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// Under Section 7 of GPL version 3, you are granted additional
// permissions described in the GCC Runtime Library Exception, version
// 3.1, as published by the Free Software Foundation.
// You should have received a copy of the GNU General Public License and
// a copy of the GCC Runtime Library Exception along with this program;
// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
// <http://www.gnu.org/licenses/>.
/** @file ext/random
* This file is a GNU extension to the Standard C++ Library.
*/
#ifndef _EXT_RANDOM
#define _EXT_RANDOM 1
#pragma GCC system_header
#include <random>
#include <array>
#ifdef __SSE2__
# include <x86intrin.h>
#endif
namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
{
_GLIBCXX_BEGIN_NAMESPACE_VERSION
/* Mersenne twister implementation optimized for vector operations.
*
* Reference: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/
*/
template<typename _UIntType, size_t __m,
size_t __pos1, size_t __sl1, size_t __sl2,
size_t __sr1, size_t __sr2,
uint32_t __msk1, uint32_t __msk2,
uint32_t __msk3, uint32_t __msk4,
uint32_t __parity1, uint32_t __parity2,
uint32_t __parity3, uint32_t __parity4>
class simd_fast_mersenne_twister_engine
{
static_assert(std::is_unsigned<_UIntType>::value, "template argument "
"substituting _UIntType not an unsigned integral type");
static_assert(__sr1 < 32, "first right shift too large");
static_assert(__sr2 < 16, "second right shift too large");
static_assert(__sl1 < 32, "first left shift too large");
static_assert(__sl2 < 16, "second left shift too large");
public:
typedef _UIntType result_type;
private:
static constexpr size_t m_w = sizeof(result_type) * 8;
static constexpr size_t _M_nstate = __m / 128 + 1;
static constexpr size_t _M_nstate32 = _M_nstate * 4;
static_assert(std::is_unsigned<_UIntType>::value, "template argument "
"substituting _UIntType not an unsigned integral type");
static_assert(__pos1 < _M_nstate, "POS1 not smaller than state size");
static_assert(16 % sizeof(_UIntType) == 0,
"UIntType size must divide 16");
public:
static constexpr size_t state_size = _M_nstate * (16
/ sizeof(result_type));
static constexpr result_type default_seed = 5489u;
// constructors and member function
explicit
simd_fast_mersenne_twister_engine(result_type __sd = default_seed)
{ seed(__sd); }
template<typename _Sseq, typename = typename
std::enable_if<!std::is_same<_Sseq,
simd_fast_mersenne_twister_engine>::value>
::type>
explicit
simd_fast_mersenne_twister_engine(_Sseq& __q)
{ seed(__q); }
void
seed(result_type __sd = default_seed);
template<typename _Sseq>
typename std::enable_if<std::is_class<_Sseq>::value>::type
seed(_Sseq& __q);
static constexpr result_type
min()
{ return 0; };
static constexpr result_type
max()
{ return std::numeric_limits<result_type>::max(); }
void
discard(unsigned long long __z);
result_type
operator()()
{
if (__builtin_expect(_M_pos >= state_size, 0))
_M_gen_rand();
return _M_stateT[_M_pos++];
}
template<typename _UIntType_2, size_t __m_2,
size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
size_t __sr1_2, size_t __sr2_2,
uint32_t __msk1_2, uint32_t __msk2_2,
uint32_t __msk3_2, uint32_t __msk4_2,
uint32_t __parity1_2, uint32_t __parity2_2,
uint32_t __parity3_2, uint32_t __parity4_2>
friend bool
operator==(const simd_fast_mersenne_twister_engine<_UIntType_2,
__m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
__msk1_2, __msk2_2, __msk3_2, __msk4_2,
__parity1_2, __parity2_2, __parity3_2, __parity4_2>& __lhs,
const simd_fast_mersenne_twister_engine<_UIntType_2,
__m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
__msk1_2, __msk2_2, __msk3_2, __msk4_2,
__parity1_2, __parity2_2, __parity3_2, __parity4_2>& __rhs);
template<typename _UIntType_2, size_t __m_2,
size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
size_t __sr1_2, size_t __sr2_2,
uint32_t __msk1_2, uint32_t __msk2_2,
uint32_t __msk3_2, uint32_t __msk4_2,
uint32_t __parity1_2, uint32_t __parity2_2,
uint32_t __parity3_2, uint32_t __parity4_2,
typename _CharT, typename _Traits>
friend std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
const __gnu_cxx::simd_fast_mersenne_twister_engine
<_UIntType_2,
__m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
__msk1_2, __msk2_2, __msk3_2, __msk4_2,
__parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
template<typename _UIntType_2, size_t __m_2,
size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
size_t __sr1_2, size_t __sr2_2,
uint32_t __msk1_2, uint32_t __msk2_2,
uint32_t __msk3_2, uint32_t __msk4_2,
uint32_t __parity1_2, uint32_t __parity2_2,
uint32_t __parity3_2, uint32_t __parity4_2,
typename _CharT, typename _Traits>
friend std::basic_istream<_CharT, _Traits>&
operator>>(std::basic_istream<_CharT, _Traits>& __is,
__gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType_2,
__m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
__msk1_2, __msk2_2, __msk3_2, __msk4_2,
__parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
private:
union
{
#ifdef __SSE2__
__m128i _M_state[_M_nstate];
#endif
uint32_t _M_state32[_M_nstate32];
result_type _M_stateT[state_size];
} __attribute__ ((__aligned__ (16)));
size_t _M_pos;
void _M_gen_rand(void);
void _M_period_certification();
};
template<typename _UIntType, size_t __m,
size_t __pos1, size_t __sl1, size_t __sl2,
size_t __sr1, size_t __sr2,
uint32_t __msk1, uint32_t __msk2,
uint32_t __msk3, uint32_t __msk4,
uint32_t __parity1, uint32_t __parity2,
uint32_t __parity3, uint32_t __parity4>
inline bool
operator!=(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
__m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
__msk4, __parity1, __parity2, __parity3, __parity4>& __lhs,
const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
__m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
__msk4, __parity1, __parity2, __parity3, __parity4>& __rhs)
{ return !(__lhs == __rhs); }
/* Definitions for the SIMD-oriented Fast Mersenne Twister as defined
* in the C implementation by Daito and Matsumoto, as both a 32-bit
* and 64-bit version.
*/
typedef simd_fast_mersenne_twister_engine<uint32_t, 607, 2,
15, 3, 13, 3,
0xfdff37ffU, 0xef7f3f7dU,
0xff777b7dU, 0x7ff7fb2fU,
0x00000001U, 0x00000000U,
0x00000000U, 0x5986f054U>
sfmt607;
typedef simd_fast_mersenne_twister_engine<uint64_t, 607, 2,
15, 3, 13, 3,
0xfdff37ffU, 0xef7f3f7dU,
0xff777b7dU, 0x7ff7fb2fU,
0x00000001U, 0x00000000U,
0x00000000U, 0x5986f054U>
sfmt607_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 1279, 7,
14, 3, 5, 1,
0xf7fefffdU, 0x7fefcfffU,
0xaff3ef3fU, 0xb5ffff7fU,
0x00000001U, 0x00000000U,
0x00000000U, 0x20000000U>
sfmt1279;
typedef simd_fast_mersenne_twister_engine<uint64_t, 1279, 7,
14, 3, 5, 1,
0xf7fefffdU, 0x7fefcfffU,
0xaff3ef3fU, 0xb5ffff7fU,
0x00000001U, 0x00000000U,
0x00000000U, 0x20000000U>
sfmt1279_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 2281, 12,
19, 1, 5, 1,
0xbff7ffbfU, 0xfdfffffeU,
0xf7ffef7fU, 0xf2f7cbbfU,
0x00000001U, 0x00000000U,
0x00000000U, 0x41dfa600U>
sfmt2281;
typedef simd_fast_mersenne_twister_engine<uint64_t, 2281, 12,
19, 1, 5, 1,
0xbff7ffbfU, 0xfdfffffeU,
0xf7ffef7fU, 0xf2f7cbbfU,
0x00000001U, 0x00000000U,
0x00000000U, 0x41dfa600U>
sfmt2281_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 4253, 17,
20, 1, 7, 1,
0x9f7bffffU, 0x9fffff5fU,
0x3efffffbU, 0xfffff7bbU,
0xa8000001U, 0xaf5390a3U,
0xb740b3f8U, 0x6c11486dU>
sfmt4253;
typedef simd_fast_mersenne_twister_engine<uint64_t, 4253, 17,
20, 1, 7, 1,
0x9f7bffffU, 0x9fffff5fU,
0x3efffffbU, 0xfffff7bbU,
0xa8000001U, 0xaf5390a3U,
0xb740b3f8U, 0x6c11486dU>
sfmt4253_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 11213, 68,
14, 3, 7, 3,
0xeffff7fbU, 0xffffffefU,
0xdfdfbfffU, 0x7fffdbfdU,
0x00000001U, 0x00000000U,
0xe8148000U, 0xd0c7afa3U>
sfmt11213;
typedef simd_fast_mersenne_twister_engine<uint64_t, 11213, 68,
14, 3, 7, 3,
0xeffff7fbU, 0xffffffefU,
0xdfdfbfffU, 0x7fffdbfdU,
0x00000001U, 0x00000000U,
0xe8148000U, 0xd0c7afa3U>
sfmt11213_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 19937, 122,
18, 1, 11, 1,
0xdfffffefU, 0xddfecb7fU,
0xbffaffffU, 0xbffffff6U,
0x00000001U, 0x00000000U,
0x00000000U, 0x13c9e684U>
sfmt19937;
typedef simd_fast_mersenne_twister_engine<uint64_t, 19937, 122,
18, 1, 11, 1,
0xdfffffefU, 0xddfecb7fU,
0xbffaffffU, 0xbffffff6U,
0x00000001U, 0x00000000U,
0x00000000U, 0x13c9e684U>
sfmt19937_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 44497, 330,
5, 3, 9, 3,
0xeffffffbU, 0xdfbebfffU,
0xbfbf7befU, 0x9ffd7bffU,
0x00000001U, 0x00000000U,
0xa3ac4000U, 0xecc1327aU>
sfmt44497;
typedef simd_fast_mersenne_twister_engine<uint64_t, 44497, 330,
5, 3, 9, 3,
0xeffffffbU, 0xdfbebfffU,
0xbfbf7befU, 0x9ffd7bffU,
0x00000001U, 0x00000000U,
0xa3ac4000U, 0xecc1327aU>
sfmt44497_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 86243, 366,
6, 7, 19, 1,
0xfdbffbffU, 0xbff7ff3fU,
0xfd77efffU, 0xbf9ff3ffU,
0x00000001U, 0x00000000U,
0x00000000U, 0xe9528d85U>
sfmt86243;
typedef simd_fast_mersenne_twister_engine<uint64_t, 86243, 366,
6, 7, 19, 1,
0xfdbffbffU, 0xbff7ff3fU,
0xfd77efffU, 0xbf9ff3ffU,
0x00000001U, 0x00000000U,
0x00000000U, 0xe9528d85U>
sfmt86243_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 132049, 110,
19, 1, 21, 1,
0xffffbb5fU, 0xfb6ebf95U,
0xfffefffaU, 0xcff77fffU,
0x00000001U, 0x00000000U,
0xcb520000U, 0xc7e91c7dU>
sfmt132049;
typedef simd_fast_mersenne_twister_engine<uint64_t, 132049, 110,
19, 1, 21, 1,
0xffffbb5fU, 0xfb6ebf95U,
0xfffefffaU, 0xcff77fffU,
0x00000001U, 0x00000000U,
0xcb520000U, 0xc7e91c7dU>
sfmt132049_64;
typedef simd_fast_mersenne_twister_engine<uint32_t, 216091, 627,
11, 3, 10, 1,
0xbff7bff7U, 0xbfffffffU,
0xbffffa7fU, 0xffddfbfbU,
0xf8000001U, 0x89e80709U,
0x3bd2b64bU, 0x0c64b1e4U>
sfmt216091;
typedef simd_fast_mersenne_twister_engine<uint64_t, 216091, 627,
11, 3, 10, 1,
0xbff7bff7U, 0xbfffffffU,
0xbffffa7fU, 0xffddfbfbU,
0xf8000001U, 0x89e80709U,
0x3bd2b64bU, 0x0c64b1e4U>
sfmt216091_64;
/**
* @brief A beta continuous distribution for random numbers.
*
* The formula for the beta probability density function is:
* @f[
* p(x|\alpha,\beta) = \frac{1}{B(\alpha,\beta)}
* x^{\alpha - 1} (1 - x)^{\beta - 1}
* @f]
*/
template<typename _RealType = double>
class beta_distribution
{
static_assert(std::is_floating_point<_RealType>::value,
"template argument not a floating point type");
public:
/** The type of the range of the distribution. */
typedef _RealType result_type;
/** Parameter type. */
struct param_type
{
typedef beta_distribution<_RealType> distribution_type;
friend class beta_distribution<_RealType>;
explicit
param_type(_RealType __alpha_val = _RealType(1),
_RealType __beta_val = _RealType(1))
: _M_alpha(__alpha_val), _M_beta(__beta_val)
{
_GLIBCXX_DEBUG_ASSERT(_M_alpha > _RealType(0));
_GLIBCXX_DEBUG_ASSERT(_M_beta > _RealType(0));
}
_RealType
alpha() const
{ return _M_alpha; }
_RealType
beta() const
{ return _M_beta; }
friend bool
operator==(const param_type& __p1, const param_type& __p2)
{ return (__p1._M_alpha == __p2._M_alpha
&& __p1._M_beta == __p2._M_beta); }
private:
void
_M_initialize();
_RealType _M_alpha;
_RealType _M_beta;
};
public:
/**
* @brief Constructs a beta distribution with parameters
* @f$\alpha@f$ and @f$\beta@f$.
*/
explicit
beta_distribution(_RealType __alpha_val = _RealType(1),
_RealType __beta_val = _RealType(1))
: _M_param(__alpha_val, __beta_val)
{ }
explicit
beta_distribution(const param_type& __p)
: _M_param(__p)
{ }
/**
* @brief Resets the distribution state.
*/
void
reset()
{ }
/**
* @brief Returns the @f$\alpha@f$ of the distribution.
*/
_RealType
alpha() const
{ return _M_param.alpha(); }
/**
* @brief Returns the @f$\beta@f$ of the distribution.
*/
_RealType
beta() const
{ return _M_param.beta(); }
/**
* @brief Returns the parameter set of the distribution.
*/
param_type
param() const
{ return _M_param; }
/**
* @brief Sets the parameter set of the distribution.
* @param __param The new parameter set of the distribution.
*/
void
param(const param_type& __param)
{ _M_param = __param; }
/**
* @brief Returns the greatest lower bound value of the distribution.
*/
result_type
min() const
{ return result_type(0); }
/**
* @brief Returns the least upper bound value of the distribution.
*/
result_type
max() const
{ return result_type(1); }
/**
* @brief Generating functions.
*/
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{ return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p);
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng)
{ this->__generate(__f, __t, __urng, this->param()); }
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{ this->__generate_impl(__f, __t, __urng, __p); }
template<typename _UniformRandomNumberGenerator>
void
__generate(result_type* __f, result_type* __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{ this->__generate_impl(__f, __t, __urng, __p); }
/**
* @brief Inserts a %beta_distribution random number distribution
* @p __x into the output stream @p __os.
*
* @param __os An output stream.
* @param __x A %beta_distribution random number distribution.
*
* @returns The output stream with the state of @p __x inserted or in
* an error state.
*/
template<typename _RealType1, typename _CharT, typename _Traits>
friend std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
const __gnu_cxx::beta_distribution<_RealType1>& __x);
/**
* @brief Extracts a %beta_distribution random number distribution
* @p __x from the input stream @p __is.
*
* @param __is An input stream.
* @param __x A %beta_distribution random number generator engine.
*
* @returns The input stream with @p __x extracted or in an error state.
*/
template<typename _RealType1, typename _CharT, typename _Traits>
friend std::basic_istream<_CharT, _Traits>&
operator>>(std::basic_istream<_CharT, _Traits>& __is,
__gnu_cxx::beta_distribution<_RealType1>& __x);
private:
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p);
param_type _M_param;
};
/**
* @brief Return true if two beta distributions have the same
* parameters and the sequences that would be generated
* are equal.
*/
template<typename _RealType>
inline bool
operator==(const __gnu_cxx::beta_distribution<_RealType>& __d1,
const __gnu_cxx::beta_distribution<_RealType>& __d2)
{ return __d1.param() == __d2.param(); }
/**
* @brief Return true if two beta distributions are different.
*/
template<typename _RealType>
inline bool
operator!=(const __gnu_cxx::beta_distribution<_RealType>& __d1,
const __gnu_cxx::beta_distribution<_RealType>& __d2)
{ return !(__d1 == __d2); }
/**
* @brief A multi-variate normal continuous distribution for random numbers.
*
* The formula for the normal probability density function is
* @f[
* p(\overrightarrow{x}|\overrightarrow{\mu },\Sigma) =
* \frac{1}{\sqrt{(2\pi )^k\det(\Sigma))}}
* e^{-\frac{1}{2}(\overrightarrow{x}-\overrightarrow{\mu})^\text{T}
* \Sigma ^{-1}(\overrightarrow{x}-\overrightarrow{\mu})}
* @f]
*
* where @f$\overrightarrow{x}@f$ and @f$\overrightarrow{\mu}@f$ are
* vectors of dimension @f$k@f$ and @f$\Sigma@f$ is the covariance
* matrix (which must be positive-definite).
*/
template<std::size_t _Dimen, typename _RealType = double>
class normal_mv_distribution
{
static_assert(std::is_floating_point<_RealType>::value,
"template argument not a floating point type");
static_assert(_Dimen != 0, "dimension is zero");
public:
/** The type of the range of the distribution. */
typedef std::array<_RealType, _Dimen> result_type;
/** Parameter type. */
class param_type
{
static constexpr size_t _M_t_size = _Dimen * (_Dimen + 1) / 2;
public:
typedef normal_mv_distribution<_Dimen, _RealType> distribution_type;
friend class normal_mv_distribution<_Dimen, _RealType>;
param_type()
{
std::fill(_M_mean.begin(), _M_mean.end(), _RealType(0));
auto __it = _M_t.begin();
for (size_t __i = 0; __i < _Dimen; ++__i)
{
std::fill_n(__it, __i, _RealType(0));
__it += __i;
*__it++ = _RealType(1);
}
}
template<typename _ForwardIterator1, typename _ForwardIterator2>
param_type(_ForwardIterator1 __meanbegin,
_ForwardIterator1 __meanend,
_ForwardIterator2 __varcovbegin,
_ForwardIterator2 __varcovend)
{
__glibcxx_function_requires(_ForwardIteratorConcept<
_ForwardIterator1>)
__glibcxx_function_requires(_ForwardIteratorConcept<
_ForwardIterator2>)
_GLIBCXX_DEBUG_ASSERT(std::distance(__meanbegin, __meanend)
<= _Dimen);
const auto __dist = std::distance(__varcovbegin, __varcovend);
_GLIBCXX_DEBUG_ASSERT(__dist == _Dimen * _Dimen
|| __dist == _Dimen * (_Dimen + 1) / 2
|| __dist == _Dimen);
if (__dist == _Dimen * _Dimen)
_M_init_full(__meanbegin, __meanend, __varcovbegin, __varcovend);
else if (__dist == _Dimen * (_Dimen + 1) / 2)
_M_init_lower(__meanbegin, __meanend, __varcovbegin, __varcovend);
else
_M_init_diagonal(__meanbegin, __meanend,
__varcovbegin, __varcovend);
}
param_type(std::initializer_list<_RealType> __mean,
std::initializer_list<_RealType> __varcov)
{
_GLIBCXX_DEBUG_ASSERT(__mean.size() <= _Dimen);
_GLIBCXX_DEBUG_ASSERT(__varcov.size() == _Dimen * _Dimen
|| __varcov.size() == _Dimen * (_Dimen + 1) / 2
|| __varcov.size() == _Dimen);
if (__varcov.size() == _Dimen * _Dimen)
_M_init_full(__mean.begin(), __mean.end(),
__varcov.begin(), __varcov.end());
else if (__varcov.size() == _Dimen * (_Dimen + 1) / 2)
_M_init_lower(__mean.begin(), __mean.end(),
__varcov.begin(), __varcov.end());
else
_M_init_diagonal(__mean.begin(), __mean.end(),
__varcov.begin(), __varcov.end());
}
std::array<_RealType, _Dimen>
mean() const
{ return _M_mean; }
std::array<_RealType, _M_t_size>
varcov() const
{ return _M_t; }
friend bool
operator==(const param_type& __p1, const param_type& __p2)
{ return __p1._M_mean == __p2._M_mean && __p1._M_t == __p2._M_t; }
private:
template <typename _InputIterator1, typename _InputIterator2>
void _M_init_full(_InputIterator1 __meanbegin,
_InputIterator1 __meanend,
_InputIterator2 __varcovbegin,
_InputIterator2 __varcovend);
template <typename _InputIterator1, typename _InputIterator2>
void _M_init_lower(_InputIterator1 __meanbegin,
_InputIterator1 __meanend,
_InputIterator2 __varcovbegin,
_InputIterator2 __varcovend);
template <typename _InputIterator1, typename _InputIterator2>
void _M_init_diagonal(_InputIterator1 __meanbegin,
_InputIterator1 __meanend,
_InputIterator2 __varbegin,
_InputIterator2 __varend);
std::array<_RealType, _Dimen> _M_mean;
std::array<_RealType, _M_t_size> _M_t;
};
public:
normal_mv_distribution()
: _M_param(), _M_nd()
{ }
template<typename _ForwardIterator1, typename _ForwardIterator2>
normal_mv_distribution(_ForwardIterator1 __meanbegin,
_ForwardIterator1 __meanend,
_ForwardIterator2 __varcovbegin,
_ForwardIterator2 __varcovend)
: _M_param(__meanbegin, __meanend, __varcovbegin, __varcovend),
_M_nd()
{ }
normal_mv_distribution(std::initializer_list<_RealType> __mean,
std::initializer_list<_RealType> __varcov)
: _M_param(__mean, __varcov), _M_nd()
{ }
explicit
normal_mv_distribution(const param_type& __p)
: _M_param(__p), _M_nd()
{ }
/**
* @brief Resets the distribution state.
*/
void
reset()
{ _M_nd.reset(); }
/**
* @brief Returns the mean of the distribution.
*/
result_type
mean() const
{ return _M_param.mean(); }
/**
* @brief Returns the compact form of the variance/covariance
* matrix of the distribution.
*/
std::array<_RealType, _Dimen * (_Dimen + 1) / 2>
varcov() const
{ return _M_param.varcov(); }
/**
* @brief Returns the parameter set of the distribution.
*/
param_type
param() const
{ return _M_param; }
/**
* @brief Sets the parameter set of the distribution.
* @param __param The new parameter set of the distribution.
*/
void
param(const param_type& __param)
{ _M_param = __param; }
/**
* @brief Returns the greatest lower bound value of the distribution.
*/
result_type
min() const
{ result_type __res;
__res.fill(std::numeric_limits<_RealType>::min());
return __res; }
/**
* @brief Returns the least upper bound value of the distribution.
*/
result_type
max() const
{ result_type __res;
__res.fill(std::numeric_limits<_RealType>::max());
return __res; }
/**
* @brief Generating functions.
*/
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{ return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p);
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng)
{ return this->__generate_impl(__f, __t, __urng, this->param()); }
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{ return this->__generate_impl(__f, __t, __urng, __p); }
/**
* @brief Return true if two multi-variant normal distributions have
* the same parameters and the sequences that would
* be generated are equal.
*/
template<size_t _Dimen1, typename _RealType1>
friend bool
operator==(const
__gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
__d1,
const
__gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
__d2);
/**
* @brief Inserts a %normal_mv_distribution random number distribution
* @p __x into the output stream @p __os.
*
* @param __os An output stream.
* @param __x A %normal_mv_distribution random number distribution.
*
* @returns The output stream with the state of @p __x inserted or in
* an error state.
*/
template<size_t _Dimen1, typename _RealType1,
typename _CharT, typename _Traits>
friend std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
const
__gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
__x);
/**
* @brief Extracts a %normal_mv_distribution random number distribution
* @p __x from the input stream @p __is.
*
* @param __is An input stream.
* @param __x A %normal_mv_distribution random number generator engine.
*
* @returns The input stream with @p __x extracted or in an error
* state.
*/
template<size_t _Dimen1, typename _RealType1,
typename _CharT, typename _Traits>
friend std::basic_istream<_CharT, _Traits>&
operator>>(std::basic_istream<_CharT, _Traits>& __is,
__gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
__x);
private:
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p);
param_type _M_param;
std::normal_distribution<_RealType> _M_nd;
};
/**
* @brief Return true if two multi-variate normal distributions are
* different.
*/
template<size_t _Dimen, typename _RealType>
inline bool
operator!=(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
__d1,
const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
__d2)
{ return !(__d1 == __d2); }
/**
* @brief A Rice continuous distribution for random numbers.
*
* The formula for the Rice probability density function is
* @f[
* p(x|\nu,\sigma) = \frac{x}{\sigma^2}
* \exp\left(-\frac{x^2+\nu^2}{2\sigma^2}\right)
* I_0\left(\frac{x \nu}{\sigma^2}\right)
* @f]
* where @f$I_0(z)@f$ is the modified Bessel function of the first kind
* of order 0 and @f$\nu >= 0@f$ and @f$\sigma > 0@f$.
*
* <table border=1 cellpadding=10 cellspacing=0>
* <caption align=top>Distribution Statistics</caption>
* <tr><td>Mean</td><td>@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
* <tr><td>Variance</td><td>@f$2\sigma^2 + \nu^2
* + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
* <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
* </table>
* where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2.
*/
template<typename _RealType = double>
class
rice_distribution
{
static_assert(std::is_floating_point<_RealType>::value,
"template argument not a floating point type");
public:
/** The type of the range of the distribution. */
typedef _RealType result_type;
/** Parameter type. */
struct param_type
{
typedef rice_distribution<result_type> distribution_type;
param_type(result_type __nu = result_type(0),
result_type __sigma = result_type(1))
: _M_nu(__nu), _M_sigma(__sigma)
{
_GLIBCXX_DEBUG_ASSERT(_M_nu >= result_type(0));
_GLIBCXX_DEBUG_ASSERT(_M_sigma > result_type(0));
}
result_type
nu() const
{ return _M_nu; }
result_type
sigma() const
{ return _M_sigma; }
friend bool
operator==(const param_type& __p1, const param_type& __p2)
{ return __p1._M_nu == __p2._M_nu
&& __p1._M_sigma == __p2._M_sigma; }
private:
void _M_initialize();
result_type _M_nu;
result_type _M_sigma;
};
/**
* @brief Constructors.
*/
explicit
rice_distribution(result_type __nu = result_type(0),
result_type __sigma = result_type(1))
: _M_param(__nu, __sigma),
_M_ndx(__nu, __sigma),
_M_ndy(result_type(0), __sigma)
{ }
explicit
rice_distribution(const param_type& __p)
: _M_param(__p),
_M_ndx(__p.nu(), __p.sigma()),
_M_ndy(result_type(0), __p.sigma())
{ }
/**
* @brief Resets the distribution state.
*/
void
reset()
{
_M_ndx.reset();
_M_ndy.reset();
}
/**
* @brief Return the parameters of the distribution.
*/
result_type
nu() const
{ return _M_param.nu(); }
result_type
sigma() const
{ return _M_param.sigma(); }
/**
* @brief Returns the parameter set of the distribution.
*/
param_type
param() const
{ return _M_param; }
/**
* @brief Sets the parameter set of the distribution.
* @param __param The new parameter set of the distribution.
*/
void
param(const param_type& __param)
{ _M_param = __param; }
/**
* @brief Returns the greatest lower bound value of the distribution.
*/
result_type
min() const
{ return result_type(0); }
/**
* @brief Returns the least upper bound value of the distribution.
*/
result_type
max() const
{ return std::numeric_limits<result_type>::max(); }
/**
* @brief Generating functions.
*/
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{
result_type __x = this->_M_ndx(__urng);
result_type __y = this->_M_ndy(__urng);
return std::hypot(__x, __y);
}
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
typename std::normal_distribution<result_type>::param_type
__px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
result_type __x = this->_M_ndx(__px, __urng);
result_type __y = this->_M_ndy(__py, __urng);
return std::hypot(__x, __y);
}
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng)
{ this->__generate(__f, __t, __urng, this->param()); }
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{ this->__generate_impl(__f, __t, __urng, __p); }
template<typename _UniformRandomNumberGenerator>
void
__generate(result_type* __f, result_type* __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{ this->__generate_impl(__f, __t, __urng, __p); }
/**
* @brief Return true if two Rice distributions have
* the same parameters and the sequences that would
* be generated are equal.
*/
friend bool
operator==(const rice_distribution& __d1,
const rice_distribution& __d2)
{ return (__d1.param() == __d2.param()
&& __d1._M_ndx == __d2._M_ndx
&& __d1._M_ndy == __d2._M_ndy); }
/**
* @brief Inserts a %rice_distribution random number distribution
* @p __x into the output stream @p __os.
*
* @param __os An output stream.
* @param __x A %rice_distribution random number distribution.
*
* @returns The output stream with the state of @p __x inserted or in
* an error state.
*/
template<typename _RealType1, typename _CharT, typename _Traits>
friend std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>&,
const rice_distribution<_RealType1>&);
/**
* @brief Extracts a %rice_distribution random number distribution
* @p __x from the input stream @p __is.
*
* @param __is An input stream.
* @param __x A %rice_distribution random number
* generator engine.
*
* @returns The input stream with @p __x extracted or in an error state.
*/
template<typename _RealType1, typename _CharT, typename _Traits>
friend std::basic_istream<_CharT, _Traits>&
operator>>(std::basic_istream<_CharT, _Traits>&,
rice_distribution<_RealType1>&);
private:
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p);
param_type _M_param;
std::normal_distribution<result_type> _M_ndx;
std::normal_distribution<result_type> _M_ndy;
};
/**
* @brief Return true if two Rice distributions are not equal.
*/
template<typename _RealType1>
inline bool
operator!=(const rice_distribution<_RealType1>& __d1,
const rice_distribution<_RealType1>& __d2)
{ return !(__d1 == __d2); }
/**
* @brief A Nakagami continuous distribution for random numbers.
*
* The formula for the Nakagami probability density function is
* @f[
* p(x|\mu,\omega) = \frac{2\mu^\mu}{\Gamma(\mu)\omega^\mu}
* x^{2\mu-1}e^{-\mu x / \omega}
* @f]
* where @f$\Gamma(z)@f$ is the gamma function and @f$\mu >= 0.5@f$
* and @f$\omega > 0@f$.
*/
template<typename _RealType = double>
class
nakagami_distribution
{
static_assert(std::is_floating_point<_RealType>::value,
"template argument not a floating point type");
public:
/** The type of the range of the distribution. */
typedef _RealType result_type;
/** Parameter type. */
struct param_type
{
typedef nakagami_distribution<result_type> distribution_type;
param_type(result_type __mu = result_type(1),
result_type __omega = result_type(1))
: _M_mu(__mu), _M_omega(__omega)
{
_GLIBCXX_DEBUG_ASSERT(_M_mu >= result_type(0.5L));
_GLIBCXX_DEBUG_ASSERT(_M_omega > result_type(0));
}
result_type
mu() const
{ return _M_mu; }
result_type
omega() const
{ return _M_omega; }
friend bool
operator==(const param_type& __p1, const param_type& __p2)
{ return __p1._M_mu == __p2._M_mu
&& __p1._M_omega == __p2._M_omega; }
private:
void _M_initialize();
result_type _M_mu;
result_type _M_omega;
};
/**
* @brief Constructors.
*/
explicit
nakagami_distribution(result_type __mu = result_type(1),
result_type __omega = result_type(1))
: _M_param(__mu, __omega),
_M_gd(__mu, __omega / __mu)
{ }
explicit
nakagami_distribution(const param_type& __p)
: _M_param(__p),
_M_gd(__p.mu(), __p.omega() / __p.mu())
{ }
/**
* @brief Resets the distribution state.
*/
void
reset()
{ _M_gd.reset(); }
/**
* @brief Return the parameters of the distribution.
*/
result_type
mu() const
{ return _M_param.mu(); }
result_type
omega() const
{ return _M_param.omega(); }
/**
* @brief Returns the parameter set of the distribution.
*/
param_type
param() const
{ return _M_param; }
/**
* @brief Sets the parameter set of the distribution.
* @param __param The new parameter set of the distribution.
*/
void
param(const param_type& __param)
{ _M_param = __param; }
/**
* @brief Returns the greatest lower bound value of the distribution.
*/
result_type
min() const
{ return result_type(0); }
/**
* @brief Returns the least upper bound value of the distribution.
*/
result_type
max() const
{ return std::numeric_limits<result_type>::max(); }
/**
* @brief Generating functions.
*/
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{ return std::sqrt(this->_M_gd(__urng)); }
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
typename std::gamma_distribution<result_type>::param_type
__pg(__p.mu(), __p.omega() / __p.mu());
return std::sqrt(this->_M_gd(__pg, __urng));
}
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng)
{ this->__generate(__f, __t, __urng, this->param()); }
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{ this->__generate_impl(__f, __t, __urng, __p); }
template<typename _UniformRandomNumberGenerator>
void
__generate(result_type* __f, result_type* __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{ this->__generate_impl(__f, __t, __urng, __p); }
/**
* @brief Return true if two Nakagami distributions have
* the same parameters and the sequences that would
* be generated are equal.
*/
friend bool
operator==(const nakagami_distribution& __d1,
const nakagami_distribution& __d2)
{ return (__d1.param() == __d2.param()
&& __d1._M_gd == __d2._M_gd); }
/**
* @brief Inserts a %nakagami_distribution random number distribution
* @p __x into the output stream @p __os.
*
* @param __os An output stream.
* @param __x A %nakagami_distribution random number distribution.
*
* @returns The output stream with the state of @p __x inserted or in
* an error state.
*/
template<typename _RealType1, typename _CharT, typename _Traits>
friend std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>&,
const nakagami_distribution<_RealType1>&);
/**
* @brief Extracts a %nakagami_distribution random number distribution
* @p __x from the input stream @p __is.
*
* @param __is An input stream.
* @param __x A %nakagami_distribution random number
* generator engine.
*
* @returns The input stream with @p __x extracted or in an error state.
*/
template<typename _RealType1, typename _CharT, typename _Traits>
friend std::basic_istream<_CharT, _Traits>&
operator>>(std::basic_istream<_CharT, _Traits>&,
nakagami_distribution<_RealType1>&);
private:
template<typename _ForwardIterator,
typename _UniformRandomNumberGenerator>
void
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
_UniformRandomNumberGenerator& __urng,
const param_type& __p);
param_type _M_param;
std::gamma_distribution<result_type> _M_gd;
};
/**
* @brief Return true if two Nakagami distributions are not equal.
*/
template<typename _RealType>
inline bool
operator!=(const nakagami_distribution<_RealType>& __d1,
const nakagami_distribution<_RealType>& __d2)
{ return !(__d1 == __d2); }
_GLIBCXX_END_NAMESPACE_VERSION
} // namespace __gnu_cxx
#include "opt_random.h"
#include "random.tcc"
#endif /* _EXT_RANDOM */
Reference in New Issue View Git Blame Copy Permalink