Subversion Repositories Kolibri OS

// Random number extensions -*- C++ -*-

// Copyright (C) 2012-2015 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

// .

/** @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

#if __cplusplus < 201103L

# include 

#else

#include 

#include 

#include 

#include 

#ifdef __SSE2__

# include 

#endif

#ifdef _GLIBCXX_USE_C99_STDINT_TR1

namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)

{

_GLIBCXX_BEGIN_NAMESPACE_VERSION

#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__

  /* Mersenne twister implementation optimized for vector operations.

   *

   * Reference: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/

   */

  template

	   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

	std::enable_if

				     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 std::enable_if::value>::type

	seed(_Sseq& __q);

      static constexpr result_type

      min()

      { return 0; };

      static constexpr result_type

      max()

      { return std::numeric_limits::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

	       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

	       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

	       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

	   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

					    15, 3, 13, 3,

					    0xfdff37ffU, 0xef7f3f7dU,

					    0xff777b7dU, 0x7ff7fb2fU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x5986f054U>

    sfmt607;

  typedef simd_fast_mersenne_twister_engine

					    15, 3, 13, 3,

					    0xfdff37ffU, 0xef7f3f7dU,

					    0xff777b7dU, 0x7ff7fb2fU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x5986f054U>

    sfmt607_64;

  typedef simd_fast_mersenne_twister_engine

					    14, 3, 5, 1,

					    0xf7fefffdU, 0x7fefcfffU,

					    0xaff3ef3fU, 0xb5ffff7fU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x20000000U>

    sfmt1279;

  typedef simd_fast_mersenne_twister_engine

					    14, 3, 5, 1,

					    0xf7fefffdU, 0x7fefcfffU,

					    0xaff3ef3fU, 0xb5ffff7fU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x20000000U>

    sfmt1279_64;

  typedef simd_fast_mersenne_twister_engine

					    19, 1, 5, 1,

					    0xbff7ffbfU, 0xfdfffffeU,

					    0xf7ffef7fU, 0xf2f7cbbfU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x41dfa600U>

    sfmt2281;

  typedef simd_fast_mersenne_twister_engine

					    19, 1, 5, 1,

					    0xbff7ffbfU, 0xfdfffffeU,

					    0xf7ffef7fU, 0xf2f7cbbfU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x41dfa600U>

    sfmt2281_64;

  typedef simd_fast_mersenne_twister_engine

					    20, 1, 7, 1,

					    0x9f7bffffU, 0x9fffff5fU,

					    0x3efffffbU, 0xfffff7bbU,

					    0xa8000001U, 0xaf5390a3U,

					    0xb740b3f8U, 0x6c11486dU>

    sfmt4253;

  typedef simd_fast_mersenne_twister_engine

					    20, 1, 7, 1,

					    0x9f7bffffU, 0x9fffff5fU,

					    0x3efffffbU, 0xfffff7bbU,

					    0xa8000001U, 0xaf5390a3U,

					    0xb740b3f8U, 0x6c11486dU>

    sfmt4253_64;

  typedef simd_fast_mersenne_twister_engine

					    14, 3, 7, 3,

					    0xeffff7fbU, 0xffffffefU,

					    0xdfdfbfffU, 0x7fffdbfdU,

					    0x00000001U, 0x00000000U,

					    0xe8148000U, 0xd0c7afa3U>

    sfmt11213;

  typedef simd_fast_mersenne_twister_engine

					    14, 3, 7, 3,

					    0xeffff7fbU, 0xffffffefU,

					    0xdfdfbfffU, 0x7fffdbfdU,

					    0x00000001U, 0x00000000U,

					    0xe8148000U, 0xd0c7afa3U>

    sfmt11213_64;

  typedef simd_fast_mersenne_twister_engine

					    18, 1, 11, 1,

					    0xdfffffefU, 0xddfecb7fU,

					    0xbffaffffU, 0xbffffff6U,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x13c9e684U>

    sfmt19937;

  typedef simd_fast_mersenne_twister_engine

					    18, 1, 11, 1,

					    0xdfffffefU, 0xddfecb7fU,

					    0xbffaffffU, 0xbffffff6U,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0x13c9e684U>

    sfmt19937_64;

  typedef simd_fast_mersenne_twister_engine

					    5, 3, 9, 3,

					    0xeffffffbU, 0xdfbebfffU,

					    0xbfbf7befU, 0x9ffd7bffU,

					    0x00000001U, 0x00000000U,

					    0xa3ac4000U, 0xecc1327aU>

    sfmt44497;

  typedef simd_fast_mersenne_twister_engine

					    5, 3, 9, 3,

					    0xeffffffbU, 0xdfbebfffU,

					    0xbfbf7befU, 0x9ffd7bffU,

					    0x00000001U, 0x00000000U,

					    0xa3ac4000U, 0xecc1327aU>

    sfmt44497_64;

  typedef simd_fast_mersenne_twister_engine

					    6, 7, 19, 1,

					    0xfdbffbffU, 0xbff7ff3fU,

					    0xfd77efffU, 0xbf9ff3ffU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0xe9528d85U>

    sfmt86243;

  typedef simd_fast_mersenne_twister_engine

					    6, 7, 19, 1,

					    0xfdbffbffU, 0xbff7ff3fU,

					    0xfd77efffU, 0xbf9ff3ffU,

					    0x00000001U, 0x00000000U,

					    0x00000000U, 0xe9528d85U>

    sfmt86243_64;

  typedef simd_fast_mersenne_twister_engine

					    19, 1, 21, 1,

					    0xffffbb5fU, 0xfb6ebf95U,

					    0xfffefffaU, 0xcff77fffU,

					    0x00000001U, 0x00000000U,

					    0xcb520000U, 0xc7e91c7dU>

    sfmt132049;

  typedef simd_fast_mersenne_twister_engine

					    19, 1, 21, 1,

					    0xffffbb5fU, 0xfb6ebf95U,

					    0xfffefffaU, 0xcff77fffU,

					    0x00000001U, 0x00000000U,

					    0xcb520000U, 0xc7e91c7dU>

    sfmt132049_64;

  typedef simd_fast_mersenne_twister_engine

					    11, 3, 10, 1,

					    0xbff7bff7U, 0xbfffffffU,

					    0xbffffa7fU, 0xffddfbfbU,

					    0xf8000001U, 0x89e80709U,

					    0x3bd2b64bU, 0x0c64b1e4U>

    sfmt216091;

  typedef simd_fast_mersenne_twister_engine

					    11, 3, 10, 1,

					    0xbff7bff7U, 0xbfffffffU,

					    0xbffffa7fU, 0xffddfbfbU,

					    0xf8000001U, 0x89e80709U,

					    0x3bd2b64bU, 0x0c64b1e4U>

    sfmt216091_64;

#endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__

  /**

   * @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

    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

	result_type

	operator()(_UniformRandomNumberGenerator& __urng)

	{ return this->operator()(__urng, _M_param); }

      template

	result_type

	operator()(_UniformRandomNumberGenerator& __urng,

		   const param_type& __p);

      template

	       typename _UniformRandomNumberGenerator>

	void

	__generate(_ForwardIterator __f, _ForwardIterator __t,

		   _UniformRandomNumberGenerator& __urng)

	{ this->__generate(__f, __t, __urng, _M_param); }

      template

	       typename _UniformRandomNumberGenerator>

	void

	__generate(_ForwardIterator __f, _ForwardIterator __t,

		   _UniformRandomNumberGenerator& __urng,

		   const param_type& __p)

	{ this->__generate_impl(__f, __t, __urng, __p); }

      template

	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 beta distributions have the same

       *        parameters and the sequences that would be generated

       *        are equal.

       */

      friend bool

      operator==(const beta_distribution& __d1,

		 const beta_distribution& __d2)

      { return __d1._M_param == __d2._M_param; }

      /**

       * @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

	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

	friend std::basic_istream<_CharT, _Traits>&

	operator>>(std::basic_istream<_CharT, _Traits>& __is,

		   __gnu_cxx::beta_distribution<_RealType1>& __x);

    private:

      template

	       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 are different.

   */

  template

    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

    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

	  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 

	  void _M_init_full(_InputIterator1 __meanbegin,

			    _InputIterator1 __meanend,

			    _InputIterator2 __varcovbegin,

			    _InputIterator2 __varcovend);

	template 

	  void _M_init_lower(_InputIterator1 __meanbegin,

			     _InputIterator1 __meanend,

			     _InputIterator2 __varcovbegin,

			     _InputIterator2 __varcovend);

	template 

	  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

	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>::lowest());

	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

	result_type

	operator()(_UniformRandomNumberGenerator& __urng)

	{ return this->operator()(__urng, _M_param); }

      template

	result_type

	operator()(_UniformRandomNumberGenerator& __urng,

		   const param_type& __p);

      template

	       typename _UniformRandomNumberGenerator>

	void

	__generate(_ForwardIterator __f, _ForwardIterator __t,

		   _UniformRandomNumberGenerator& __urng)

	{ return this->__generate_impl(__f, __t, __urng, _M_param); }

      template

	       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

	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

	       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

	       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 _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

    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$.

   *

   * 

        
        

           

           

           

           
   * 
Distribution Statistics
        
        

           

           

           

           
   * 
Mean
@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$
        
        

           

           

           

           
   * 
Variance
@f$2\sigma^2 + \nu^2
        
        

           

           

           

           
   *                   + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$
        
        

           

           

           

           
   * 
Range
@f$[0, \infty)@f$
        
        

           

           

           

           
   * 

   * where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2.

   */

  template

    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 distribution_type;

	param_type(result_type __nu_val = result_type(0),

		   result_type __sigma_val = result_type(1))

	: _M_nu(__nu_val), _M_sigma(__sigma_val)

	{

	  _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_val = result_type(0),

			result_type __sigma_val = result_type(1))

      : _M_param(__nu_val, __sigma_val),

	_M_ndx(__nu_val, __sigma_val),

	_M_ndy(result_type(0), __sigma_val)

      { }

      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::max(); }

      /**

       * @brief Generating functions.

       */

      template

	result_type

	operator()(_UniformRandomNumberGenerator& __urng)

	{

	  result_type __x = this->_M_ndx(__urng);

	  result_type __y = this->_M_ndy(__urng);

#if _GLIBCXX_USE_C99_MATH_TR1

	  return std::hypot(__x, __y);

#else

	  return std::sqrt(__x * __x + __y * __y);

#endif

	}

      template

	result_type

	operator()(_UniformRandomNumberGenerator& __urng,

		   const param_type& __p)

	{

	  typename std::normal_distribution::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);

#if _GLIBCXX_USE_C99_MATH_TR1

	  return std::hypot(__x, __y);

#else

	  return std::sqrt(__x * __x + __y * __y);

#endif

	}

      template

	       typename _UniformRandomNumberGenerator>

	void

	__generate(_ForwardIterator __f, _ForwardIterator __t,

		   _UniformRandomNumberGenerator& __urng)

	{ this->__generate(__f, __t, __urng, _M_param); }

      template

	       typename _UniformRandomNumberGenerator>

	void

	__generate(_ForwardIterator __f, _ForwardIterator __t,

		   _UniformRandomNumberGenerator& __urng,

		   const param_type& __p)

	{ this->__generate_impl(__f, __t, __urng, __p); }

      template

	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._M_param == __d2._M_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

	friend std::basic_ostream<_CharT, _Traits>&