* include/ext/random: Add __gnu_cxx:normal_mv_distribution<> class.

	* include/ext/random.tccAdd out-of-line functions for
	__gnu_cxx::normal_mv_distribution<>.
	* testsuite/26_numerics/random/normal_mv_distribution/
	operators/equal.cc: New file.
	* testsuite/26_numerics/random/normal_mv_distribution/
	operators/serialize.cc: New file.
	* testsuite/26_numerics/random/normal_mv_distribution/
	operators/inequal.cc: New file.
	* testsuite/26_numerics/random/normal_mv_distribution/
	cons/default.cc: New file.
	* testsuite/26_numerics/random/normal_mv_distribution/
	cons/parms.cc: New file.
	* testsuite/26_numerics/random/normal_mv_distribution/
	requirements/explicit_instantiation/1.cc: New file.
	* testsuite/26_numerics/random/normal_mv_distribution/
	requirements/typedefs.cc: New file.



git-svn-id: https://gcc.gnu.org/svn/gcc/trunk@190960 138bc75d-0d04-0410-961f-82ee72b054a4

2 files changed, 519 insertions, 0 deletions

diff --git a/libstdc++-v3/include/ext/random b/libstdc++-v3/include/ext/random
index 9563e6a0500..6bb438a8558 100644
--- a/libstdc++-v3/include/ext/random
+++ b/libstdc++-v3/include/ext/random
@@ -32,6 +32,7 @@
 #pragma GCC system_header
 
 #include <random>
+#include <array>
 #ifdef __SSE2__
 # include <x86intrin.h>
 #endif
@@ -590,6 +591,310 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
     { 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); }
+
 
 _GLIBCXX_END_NAMESPACE_VERSION
 } // namespace std
diff --git a/libstdc++-v3/include/ext/random.tcc b/libstdc++-v3/include/ext/random.tcc
index 1776b0df163..0fa006af0bd 100644
--- a/libstdc++-v3/include/ext/random.tcc
+++ b/libstdc++-v3/include/ext/random.tcc
@@ -538,6 +538,220 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
       return __is;
     }
 
+
+  template<std::size_t _Dimen, typename _RealType>
+    template<typename _InputIterator1, typename _InputIterator2>
+      void
+      normal_mv_distribution<_Dimen, _RealType>::param_type::
+      _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
+		   _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
+      {
+	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
+	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
+	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
+		  _M_mean.end(), _RealType(0));
+
+	// Perform the Cholesky decomposition
+	auto __w = _M_t.begin();
+	for (size_t __j = 0; __j < _Dimen; ++__j)
+	  {
+	    _RealType __sum = _RealType(0);
+
+	    auto __slitbegin = __w;
+	    auto __cit = _M_t.begin();
+	    for (size_t __i = 0; __i < __j; ++__i)
+	      {
+		auto __slit = __slitbegin;
+		_RealType __s = *__varcovbegin++;
+		for (size_t __k = 0; __k < __i; ++__k)
+		  __s -= *__slit++ * *__cit++;
+
+		*__w++ = __s /= *__cit++;
+		__sum += __s * __s;
+	      }
+
+	    __sum = *__varcovbegin - __sum;
+	    if (__builtin_expect(__sum <= _RealType(0), 0))
+	      std::__throw_runtime_error(__N("normal_mv_distribution::"
+					     "param_type::_M_init_full"));
+	    *__w++ = std::sqrt(__sum);
+
+	    std::advance(__varcovbegin, _Dimen - __j);
+	  }
+      }
+
+  template<std::size_t _Dimen, typename _RealType>
+    template<typename _InputIterator1, typename _InputIterator2>
+      void
+      normal_mv_distribution<_Dimen, _RealType>::param_type::
+      _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
+		    _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
+      {
+	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
+	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
+	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
+		  _M_mean.end(), _RealType(0));
+
+	// Perform the Cholesky decomposition
+	auto __w = _M_t.begin();
+	for (size_t __j = 0; __j < _Dimen; ++__j)
+	  {
+	    _RealType __sum = _RealType(0);
+
+	    auto __slitbegin = __w;
+	    auto __cit = _M_t.begin();
+	    for (size_t __i = 0; __i < __j; ++__i)
+	      {
+		auto __slit = __slitbegin;
+		_RealType __s = *__varcovbegin++;
+		for (size_t __k = 0; __k < __i; ++__k)
+		  __s -= *__slit++ * *__cit++;
+
+		*__w++ = __s /= *__cit++;
+		__sum += __s * __s;
+	      }
+
+	    __sum = *__varcovbegin++ - __sum;
+	    if (__builtin_expect(__sum <= _RealType(0), 0))
+	      std::__throw_runtime_error(__N("normal_mv_distribution::"
+					     "param_type::_M_init_full"));
+	    *__w++ = std::sqrt(__sum);
+	  }
+      }
+
+  template<std::size_t _Dimen, typename _RealType>
+    template<typename _InputIterator1, typename _InputIterator2>
+      void
+      normal_mv_distribution<_Dimen, _RealType>::param_type::
+      _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
+		       _InputIterator2 __varbegin, _InputIterator2 __varend)
+      {
+	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
+	__glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
+	std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
+		  _M_mean.end(), _RealType(0));
+
+	auto __w = _M_t.begin();
+	size_t __step = 0;
+	while (__varbegin != __varend)
+	  {
+	    std::fill_n(__w, __step, _RealType(0));
+	    __w += __step++;
+	    if (__builtin_expect(*__varbegin < _RealType(0), 0))
+	      std::__throw_runtime_error(__N("normal_mv_distribution::"
+					     "param_type::_M_init_diagonal"));
+	    *__w++ = std::sqrt(*__varbegin++);
+	  }
+      }
+
+  template<std::size_t _Dimen, typename _RealType>
+    template<typename _UniformRandomNumberGenerator>
+      typename normal_mv_distribution<_Dimen, _RealType>::result_type
+      normal_mv_distribution<_Dimen, _RealType>::
+      operator()(_UniformRandomNumberGenerator& __urng,
+		 const param_type& __param)
+      {
+	result_type __ret;
+
+	for (size_t __i = 0; __i < _Dimen; ++__i)
+	  __ret[__i] = _M_nd(__urng);
+
+	auto __t_it = __param._M_t.crbegin();
+	for (size_t __i = _Dimen; __i > 0; --__i)
+	  {
+	    _RealType __sum = _RealType(0);
+	    for (size_t __j = __i; __j > 0; --__j)
+	      __sum += __ret[__j - 1] * *__t_it++;
+	    __ret[__i - 1] = __sum;
+	  }
+
+	return __ret;
+      }
+
+  template<std::size_t _Dimen, typename _RealType>
+    template<typename _ForwardIterator, typename _UniformRandomNumberGenerator>
+      void
+      normal_mv_distribution<_Dimen, _RealType>::
+      __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
+		      _UniformRandomNumberGenerator& __urng,
+		      const param_type& __param)
+      {
+	__glibcxx_function_requires(_Mutable_ForwardIteratorConcept<
+				    _ForwardIterator>)
+	while (__f != __t)
+	  *__f++ = this->operator()(__urng, __param);
+      }
+
+  template<size_t _Dimen, typename _RealType>
+    bool
+    operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
+	       __d1,
+	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
+	       __d2)
+    {
+      return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd;
+    }
+
+  template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
+    std::basic_ostream<_CharT, _Traits>&
+    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
+	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
+    {
+      typedef std::basic_ostream<_CharT, _Traits>  __ostream_type;
+      typedef typename __ostream_type::ios_base    __ios_base;
+
+      const typename __ios_base::fmtflags __flags = __os.flags();
+      const _CharT __fill = __os.fill();
+      const std::streamsize __precision = __os.precision();
+      const _CharT __space = __os.widen(' ');
+      __os.flags(__ios_base::scientific | __ios_base::left);
+      __os.fill(__space);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
+
+      auto __mean = __x._M_param.mean();
+      for (auto __it : __mean)
+	__os << __it << __space;
+      auto __t = __x._M_param.varcov();
+      for (auto __it : __t)
+	__os << __it << __space;
+
+      __os << __x._M_nd;
+
+      __os.flags(__flags);
+      __os.fill(__fill);
+      __os.precision(__precision);
+      return __os;
+    }
+
+  template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
+    std::basic_istream<_CharT, _Traits>&
+    operator>>(std::basic_istream<_CharT, _Traits>& __is,
+	       __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
+    {
+      typedef std::basic_istream<_CharT, _Traits>  __istream_type;
+      typedef typename __istream_type::ios_base    __ios_base;
+
+      const typename __ios_base::fmtflags __flags = __is.flags();
+      __is.flags(__ios_base::dec | __ios_base::skipws);
+
+      std::array<_RealType, _Dimen> __mean;
+      for (auto& __it : __mean)
+	__is >> __it;
+      std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov;
+      for (auto& __it : __varcov)
+	__is >> __it;
+
+      __is >> __x._M_nd;
+
+      __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
+		param_type(__mean.begin(), __mean.end(),
+			   __varcov.begin(), __varcov.end()));
+
+      __is.flags(__flags);
+      return __is;
+    }
+
+
 _GLIBCXX_END_NAMESPACE_VERSION
 } // namespace