path: root/libstdc++-v3/include/ext/random

// 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>
#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++];
      }

#ifdef __SSE2__
      friend bool
      operator==(const simd_fast_mersenne_twister_engine& __lhs,
		 const simd_fast_mersenne_twister_engine& __rhs)
      { __m128i __res = _mm_cmpeq_epi8(__lhs._M_state[0], __rhs._M_state[0]);
	for (size_t __i = 1; __i < __lhs._M_nstate; ++__i)
	  __res = _mm_and_si128(__res, _mm_cmpeq_epi8(__lhs._M_state[__i],
						      __rhs._M_state[__i]));
	return (_mm_movemask_epi8(__res) == 0xffff
		&& __lhs._M_pos == __rhs._M_pos); }
#else
      friend bool
      operator==(const simd_fast_mersenne_twister_engine& __lhs,
		 const simd_fast_mersenne_twister_engine& __rhs)
      { return (std::equal(__lhs._M_stateT, __lhs._M_stateT + state_size,
			   __rhs._M_stateT)
		&& __lhs._M_pos == __rhs._M_pos); }
#endif

      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 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.param() == __d2.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<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 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); }



_GLIBCXX_END_NAMESPACE_VERSION
} // namespace std

#include "random.tcc"

#endif /* _EXT_RANDOM */