Subversion Repositories altor32

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

// .

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

#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>::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

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

   *

   * 

      
      

         
920
         

         

         
   * 
Distribution Statistics
      
      

         
921
         

         

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

         
922
         

         

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

         
923
         

         

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

         
924
         

         

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

         
925
         

         

         
   * 

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

   */

  template