diff options
24 files changed, 982 insertions, 787 deletions
diff --git a/Eigen/Core b/Eigen/Core index 1f3a6504e..946ed0677 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -328,7 +328,6 @@ using std::ptrdiff_t; #include "src/Core/NumTraits.h" #include "src/Core/MathFunctions.h" -#include "src/Core/SpecialFunctions.h" #include "src/Core/GenericPacketMath.h" #if defined EIGEN_VECTORIZE_AVX diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h index 76a75dee1..16a122370 100644 --- a/Eigen/src/Core/GenericPacketMath.h +++ b/Eigen/src/Core/GenericPacketMath.h @@ -435,42 +435,6 @@ Packet pfloor(const Packet& a) { using numext::floor; return floor(a); } template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pceil(const Packet& a) { using numext::ceil; return ceil(a); } -/** \internal \returns the ln(|gamma(\a a)|) (coeff-wise) */ -template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS -Packet plgamma(const Packet& a) { using numext::lgamma; return lgamma(a); } - -/** \internal \returns the derivative of lgamma, psi(\a a) (coeff-wise) */ -template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS -Packet pdigamma(const Packet& a) { using numext::digamma; return digamma(a); } - -/** \internal \returns the zeta function of two arguments (coeff-wise) */ -template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS -Packet pzeta(const Packet& x, const Packet& q) { using numext::zeta; return zeta(x, q); } - -/** \internal \returns the polygamma function (coeff-wise) */ -template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS -Packet ppolygamma(const Packet& n, const Packet& x) { using numext::polygamma; return polygamma(n, x); } - -/** \internal \returns the erf(\a a) (coeff-wise) */ -template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS -Packet perf(const Packet& a) { using numext::erf; return erf(a); } - -/** \internal \returns the erfc(\a a) (coeff-wise) */ -template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS -Packet perfc(const Packet& a) { using numext::erfc; return erfc(a); } - -/** \internal \returns the incomplete gamma function igamma(\a a, \a x) */ -template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -Packet pigamma(const Packet& a, const Packet& x) { using numext::igamma; return igamma(a, x); } - -/** \internal \returns the complementary incomplete gamma function igammac(\a a, \a x) */ -template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -Packet pigammac(const Packet& a, const Packet& x) { using numext::igammac; return igammac(a, x); } - -/** \internal \returns the complementary incomplete gamma function betainc(\a a, \a b, \a x) */ -template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -Packet pbetainc(const Packet& a, const Packet& b,const Packet& x) { using numext::betainc; return betainc(a, b, x); } - /*************************************************************************** * The following functions might not have to be overwritten for vectorized types ***************************************************************************/ diff --git a/Eigen/src/Core/GlobalFunctions.h b/Eigen/src/Core/GlobalFunctions.h index b9c3ec25b..879a93e6b 100644 --- a/Eigen/src/Core/GlobalFunctions.h +++ b/Eigen/src/Core/GlobalFunctions.h @@ -174,111 +174,6 @@ namespace Eigen } #endif - /** \cpp11 \returns an expression of the coefficient-wise igamma(\a a, \a x) to the given arrays. - * - * This function computes the coefficient-wise incomplete gamma function. - * - * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, - * or float/double in non c++11 mode, the user has to provide implementations of igammac(T,T) for any scalar - * type T to be supported. - * - * \sa Eigen::igammac(), Eigen::lgamma() - */ - template<typename Derived,typename ExponentDerived> - inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_op<typename Derived::Scalar>, const Derived, const ExponentDerived> - igamma(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) - { - return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_op<typename Derived::Scalar>, const Derived, const ExponentDerived>( - a.derived(), - x.derived() - ); - } - - /** \cpp11 \returns an expression of the coefficient-wise igammac(\a a, \a x) to the given arrays. - * - * This function computes the coefficient-wise complementary incomplete gamma function. - * - * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, - * or float/double in non c++11 mode, the user has to provide implementations of igammac(T,T) for any scalar - * type T to be supported. - * - * \sa Eigen::igamma(), Eigen::lgamma() - */ - template<typename Derived,typename ExponentDerived> - inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igammac_op<typename Derived::Scalar>, const Derived, const ExponentDerived> - igammac(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) - { - return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igammac_op<typename Derived::Scalar>, const Derived, const ExponentDerived>( - a.derived(), - x.derived() - ); - } - - /** \cpp11 \returns an expression of the coefficient-wise polygamma(\a n, \a x) to the given arrays. - * - * It returns the \a n -th derivative of the digamma(psi) evaluated at \c x. - * - * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, - * or float/double in non c++11 mode, the user has to provide implementations of polygamma(T,T) for any scalar - * type T to be supported. - * - * \sa Eigen::digamma() - */ - // * \warning Be careful with the order of the parameters: x.polygamma(n) is equivalent to polygamma(n,x) - // * \sa ArrayBase::polygamma() - template<typename DerivedN,typename DerivedX> - inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_polygamma_op<typename DerivedX::Scalar>, const DerivedN, const DerivedX> - polygamma(const Eigen::ArrayBase<DerivedN>& n, const Eigen::ArrayBase<DerivedX>& x) - { - return Eigen::CwiseBinaryOp<Eigen::internal::scalar_polygamma_op<typename DerivedX::Scalar>, const DerivedN, const DerivedX>( - n.derived(), - x.derived() - ); - } - - /** \cpp11 \returns an expression of the coefficient-wise betainc(\a x, \a a, \a b) to the given arrays. - * - * This function computes the regularized incomplete beta function (integral). - * - * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, - * or float/double in non c++11 mode, the user has to provide implementations of betainc(T,T,T) for any scalar - * type T to be supported. - * - * \sa Eigen::betainc(), Eigen::lgamma() - */ - template<typename ArgADerived, typename ArgBDerived, typename ArgXDerived> - inline const Eigen::CwiseTernaryOp<Eigen::internal::scalar_betainc_op<typename ArgXDerived::Scalar>, const ArgADerived, const ArgBDerived, const ArgXDerived> - betainc(const Eigen::ArrayBase<ArgADerived>& a, const Eigen::ArrayBase<ArgBDerived>& b, const Eigen::ArrayBase<ArgXDerived>& x) - { - return Eigen::CwiseTernaryOp<Eigen::internal::scalar_betainc_op<typename ArgXDerived::Scalar>, const ArgADerived, const ArgBDerived, const ArgXDerived>( - a.derived(), - b.derived(), - x.derived() - ); - } - - - /** \returns an expression of the coefficient-wise zeta(\a x, \a q) to the given arrays. - * - * It returns the Riemann zeta function of two arguments \a x and \a q: - * - * \param x is the exposent, it must be > 1 - * \param q is the shift, it must be > 0 - * - * \note This function supports only float and double scalar types. To support other scalar types, the user has - * to provide implementations of zeta(T,T) for any scalar type T to be supported. - * - * \sa ArrayBase::zeta() - */ - template<typename DerivedX,typename DerivedQ> - inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_zeta_op<typename DerivedX::Scalar>, const DerivedX, const DerivedQ> - zeta(const Eigen::ArrayBase<DerivedX>& x, const Eigen::ArrayBase<DerivedQ>& q) - { - return Eigen::CwiseBinaryOp<Eigen::internal::scalar_zeta_op<typename DerivedX::Scalar>, const DerivedX, const DerivedQ>( - x.derived(), - q.derived() - ); - } namespace internal { diff --git a/Eigen/src/Core/arch/CUDA/Half.h b/Eigen/src/Core/arch/CUDA/Half.h index 87bdbfd1e..61692b83d 100644 --- a/Eigen/src/Core/arch/CUDA/Half.h +++ b/Eigen/src/Core/arch/CUDA/Half.h @@ -454,36 +454,6 @@ template <> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half maxi(const Eigen:: return f1 < f2 ? b : a; #endif } - -#if EIGEN_HAS_C99_MATH -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half lgamma(const Eigen::half& a) { - return Eigen::half(Eigen::numext::lgamma(static_cast<float>(a))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half digamma(const Eigen::half& a) { - return Eigen::half(Eigen::numext::digamma(static_cast<float>(a))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half zeta(const Eigen::half& x, const Eigen::half& q) { - return Eigen::half(Eigen::numext::zeta(static_cast<float>(x), static_cast<float>(q))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half polygamma(const Eigen::half& n, const Eigen::half& x) { - return Eigen::half(Eigen::numext::polygamma(static_cast<float>(n), static_cast<float>(x))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erf(const Eigen::half& a) { - return Eigen::half(Eigen::numext::erf(static_cast<float>(a))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erfc(const Eigen::half& a) { - return Eigen::half(Eigen::numext::erfc(static_cast<float>(a))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igamma(const Eigen::half& a, const Eigen::half& x) { - return Eigen::half(Eigen::numext::igamma(static_cast<float>(a), static_cast<float>(x))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igammac(const Eigen::half& a, const Eigen::half& x) { - return Eigen::half(Eigen::numext::igammac(static_cast<float>(a), static_cast<float>(x))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half betainc(const Eigen::half& a, const Eigen::half& b, const Eigen::half& x) { - return Eigen::half(Eigen::numext::betainc(static_cast<float>(a), static_cast<float>(b), static_cast<float>(x))); -} -#endif } // end namespace numext } // end namespace Eigen diff --git a/Eigen/src/Core/functors/BinaryFunctors.h b/Eigen/src/Core/functors/BinaryFunctors.h index 2c1331208..dc3690444 100644 --- a/Eigen/src/Core/functors/BinaryFunctors.h +++ b/Eigen/src/Core/functors/BinaryFunctors.h @@ -429,57 +429,6 @@ template<> struct functor_traits<scalar_boolean_xor_op> { }; }; -/** \internal - * \brief Template functor to compute the incomplete gamma function igamma(a, x) - * - * \sa class CwiseBinaryOp, Cwise::igamma - */ -template<typename Scalar> struct scalar_igamma_op : binary_op_base<Scalar,Scalar> -{ - EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const { - using numext::igamma; return igamma(a, x); - } - template<typename Packet> - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const { - return internal::pigamma(a, x); - } -}; -template<typename Scalar> -struct functor_traits<scalar_igamma_op<Scalar> > { - enum { - // Guesstimate - Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasIGamma - }; -}; - - -/** \internal - * \brief Template functor to compute the complementary incomplete gamma function igammac(a, x) - * - * \sa class CwiseBinaryOp, Cwise::igammac - */ -template<typename Scalar> struct scalar_igammac_op : binary_op_base<Scalar,Scalar> -{ - EIGEN_EMPTY_STRUCT_CTOR(scalar_igammac_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const { - using numext::igammac; return igammac(a, x); - } - template<typename Packet> - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const - { - return internal::pigammac(a, x); - } -}; -template<typename Scalar> -struct functor_traits<scalar_igammac_op<Scalar> > { - enum { - // Guesstimate - Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasIGammac - }; -}; //---------- binary functors bound to a constant, thus appearing as a unary functor ---------- diff --git a/Eigen/src/Core/functors/TernaryFunctors.h b/Eigen/src/Core/functors/TernaryFunctors.h index 8b9e53062..b254e96c6 100644 --- a/Eigen/src/Core/functors/TernaryFunctors.h +++ b/Eigen/src/Core/functors/TernaryFunctors.h @@ -16,29 +16,7 @@ namespace internal { //---------- associative ternary functors ---------- -/** \internal - * \brief Template functor to compute the incomplete beta integral betainc(a, b, x) - * - */ -template<typename Scalar> struct scalar_betainc_op { - EIGEN_EMPTY_STRUCT_CTOR(scalar_betainc_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& x, const Scalar& a, const Scalar& b) const { - using numext::betainc; return betainc(x, a, b); - } - template<typename Packet> - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& x, const Packet& a, const Packet& b) const - { - return internal::pbetainc(x, a, b); - } -}; -template<typename Scalar> -struct functor_traits<scalar_betainc_op<Scalar> > { - enum { - // Guesstimate - Cost = 400 * NumTraits<Scalar>::MulCost + 400 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasBetaInc - }; -}; + } // end namespace internal diff --git a/Eigen/src/Core/functors/UnaryFunctors.h b/Eigen/src/Core/functors/UnaryFunctors.h index a7d8c3b52..04208c9fe 100644 --- a/Eigen/src/Core/functors/UnaryFunctors.h +++ b/Eigen/src/Core/functors/UnaryFunctors.h @@ -473,142 +473,6 @@ struct functor_traits<scalar_asin_op<Scalar> > /** \internal - * \brief Template functor to compute the natural log of the absolute - * value of Gamma of a scalar - * \sa class CwiseUnaryOp, Cwise::lgamma() - */ -template<typename Scalar> struct scalar_lgamma_op { - EIGEN_EMPTY_STRUCT_CTOR(scalar_lgamma_op) - EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { - using numext::lgamma; return lgamma(a); - } - typedef typename packet_traits<Scalar>::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::plgamma(a); } -}; -template<typename Scalar> -struct functor_traits<scalar_lgamma_op<Scalar> > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasLGamma - }; -}; - -/** \internal - * \brief Template functor to compute psi, the derivative of lgamma of a scalar. - * \sa class CwiseUnaryOp, Cwise::digamma() - */ -template<typename Scalar> struct scalar_digamma_op { - EIGEN_EMPTY_STRUCT_CTOR(scalar_digamma_op) - EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { - using numext::digamma; return digamma(a); - } - typedef typename packet_traits<Scalar>::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::pdigamma(a); } -}; -template<typename Scalar> -struct functor_traits<scalar_digamma_op<Scalar> > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasDiGamma - }; -}; - -/** \internal - * \brief Template functor to compute the Riemann Zeta function of two arguments. - * \sa class CwiseUnaryOp, Cwise::zeta() - */ -template<typename Scalar> struct scalar_zeta_op { - EIGEN_EMPTY_STRUCT_CTOR(scalar_zeta_op) - EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& x, const Scalar& q) const { - using numext::zeta; return zeta(x, q); - } - typedef typename packet_traits<Scalar>::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& x, const Packet& q) const { return internal::pzeta(x, q); } -}; -template<typename Scalar> -struct functor_traits<scalar_zeta_op<Scalar> > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasZeta - }; -}; - -/** \internal - * \brief Template functor to compute the polygamma function. - * \sa class CwiseUnaryOp, Cwise::polygamma() - */ -template<typename Scalar> struct scalar_polygamma_op { - EIGEN_EMPTY_STRUCT_CTOR(scalar_polygamma_op) - EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& n, const Scalar& x) const { - using numext::polygamma; return polygamma(n, x); - } - typedef typename packet_traits<Scalar>::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& n, const Packet& x) const { return internal::ppolygamma(n, x); } -}; -template<typename Scalar> -struct functor_traits<scalar_polygamma_op<Scalar> > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasPolygamma - }; -}; - -/** \internal - * \brief Template functor to compute the Gauss error function of a - * scalar - * \sa class CwiseUnaryOp, Cwise::erf() - */ -template<typename Scalar> struct scalar_erf_op { - EIGEN_EMPTY_STRUCT_CTOR(scalar_erf_op) - EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { - using numext::erf; return erf(a); - } - typedef typename packet_traits<Scalar>::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perf(a); } -}; -template<typename Scalar> -struct functor_traits<scalar_erf_op<Scalar> > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasErf - }; -}; - -/** \internal - * \brief Template functor to compute the Complementary Error Function - * of a scalar - * \sa class CwiseUnaryOp, Cwise::erfc() - */ -template<typename Scalar> struct scalar_erfc_op { - EIGEN_EMPTY_STRUCT_CTOR(scalar_erfc_op) - EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { - using numext::erfc; return erfc(a); - } - typedef typename packet_traits<Scalar>::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perfc(a); } -}; -template<typename Scalar> -struct functor_traits<scalar_erfc_op<Scalar> > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, - PacketAccess = packet_traits<Scalar>::HasErfc - }; -}; - - -/** \internal * \brief Template functor to compute the atan of a scalar * \sa class CwiseUnaryOp, ArrayBase::atan() */ diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 1c90c0e2b..ea107393a 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -203,15 +203,21 @@ template<typename Scalar> struct scalar_random_op; template<typename Scalar> struct scalar_constant_op; template<typename Scalar> struct scalar_identity_op; template<typename Scalar,bool iscpx> struct scalar_sign_op; -template<typename Scalar> struct scalar_igamma_op; -template<typename Scalar> struct scalar_igammac_op; -template<typename Scalar> struct scalar_betainc_op; - template<typename Scalar,typename ScalarExponent> struct scalar_pow_op; template<typename LhsScalar,typename RhsScalar=LhsScalar> struct scalar_hypot_op; template<typename LhsScalar,typename RhsScalar=LhsScalar> struct scalar_product_op; template<typename LhsScalar,typename RhsScalar=LhsScalar> struct scalar_quotient_op; +// SpecialFunctions module +template<typename Scalar> struct scalar_lgamma_op; +template<typename Scalar> struct scalar_digamma_op; +template<typename Scalar> struct scalar_erf_op; +template<typename Scalar> struct scalar_erfc_op; +template<typename Scalar> struct scalar_igamma_op; +template<typename Scalar> struct scalar_igammac_op; +template<typename Scalar> struct scalar_zeta_op; +template<typename Scalar> struct scalar_betainc_op; + } // end namespace internal struct IOFormat; diff --git a/Eigen/src/plugins/ArrayCwiseBinaryOps.h b/Eigen/src/plugins/ArrayCwiseBinaryOps.h index 19e25ab62..62fb303d9 100644 --- a/Eigen/src/plugins/ArrayCwiseBinaryOps.h +++ b/Eigen/src/plugins/ArrayCwiseBinaryOps.h @@ -330,6 +330,8 @@ operator^(const EIGEN_CURRENT_STORAGE_BASE_CLASS<OtherDerived> &other) const #if 0 /** \cpp11 \returns an expression of the coefficient-wise polygamma function. * + * \specialfunctions_module + * * It returns the \a n -th derivative of the digamma(psi) evaluated at \c *this. * * \warning Be careful with the order of the parameters: x.polygamma(n) is equivalent to polygamma(n,x) @@ -346,6 +348,8 @@ polygamma(const EIGEN_CURRENT_STORAGE_BASE_CLASS<DerivedN> &n) const /** \returns an expression of the coefficient-wise zeta function. * + * \specialfunctions_module + * * It returns the Riemann zeta function of two arguments \c *this and \a q: * * \param *this is the exposent, it must be > 1 diff --git a/Eigen/src/plugins/ArrayCwiseUnaryOps.h b/Eigen/src/plugins/ArrayCwiseUnaryOps.h index 9e42bb540..db02e299c 100644 --- a/Eigen/src/plugins/ArrayCwiseUnaryOps.h +++ b/Eigen/src/plugins/ArrayCwiseUnaryOps.h @@ -22,10 +22,6 @@ typedef CwiseUnaryOp<internal::scalar_atan_op<Scalar>, const Derived> AtanReturn typedef CwiseUnaryOp<internal::scalar_tanh_op<Scalar>, const Derived> TanhReturnType; typedef CwiseUnaryOp<internal::scalar_sinh_op<Scalar>, const Derived> SinhReturnType; typedef CwiseUnaryOp<internal::scalar_cosh_op<Scalar>, const Derived> CoshReturnType; -typedef CwiseUnaryOp<internal::scalar_lgamma_op<Scalar>, const Derived> LgammaReturnType; -typedef CwiseUnaryOp<internal::scalar_digamma_op<Scalar>, const Derived> DigammaReturnType; -typedef CwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived> ErfReturnType; -typedef CwiseUnaryOp<internal::scalar_erfc_op<Scalar>, const Derived> ErfcReturnType; typedef CwiseUnaryOp<internal::scalar_square_op<Scalar>, const Derived> SquareReturnType; typedef CwiseUnaryOp<internal::scalar_cube_op<Scalar>, const Derived> CubeReturnType; typedef CwiseUnaryOp<internal::scalar_round_op<Scalar>, const Derived> RoundReturnType; @@ -324,77 +320,6 @@ cosh() const return CoshReturnType(derived()); } -/** \cpp11 \returns an expression of the coefficient-wise ln(|gamma(*this)|). - * - * Example: \include Cwise_lgamma.cpp - * Output: \verbinclude Cwise_lgamma.out - * - * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, - * or float/double in non c++11 mode, the user has to provide implementations of lgamma(T) for any scalar - * type T to be supported. - * - * \sa digamma() - */ -EIGEN_DEVICE_FUNC -inline const LgammaReturnType -lgamma() const -{ - return LgammaReturnType(derived()); -} - -/** \returns an expression of the coefficient-wise digamma (psi, derivative of lgamma). - * - * \note This function supports only float and double scalar types. To support other scalar types, - * the user has to provide implementations of digamma(T) for any scalar - * type T to be supported. - * - * \sa Eigen::digamma(), Eigen::polygamma(), lgamma() - */ -EIGEN_DEVICE_FUNC -inline const DigammaReturnType -digamma() const -{ - return DigammaReturnType(derived()); -} - -/** \cpp11 \returns an expression of the coefficient-wise Gauss error - * function of *this. - * - * Example: \include Cwise_erf.cpp - * Output: \verbinclude Cwise_erf.out - * - * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, - * or float/double in non c++11 mode, the user has to provide implementations of erf(T) for any scalar - * type T to be supported. - * - * \sa erfc() - */ -EIGEN_DEVICE_FUNC -inline const ErfReturnType -erf() const -{ - return ErfReturnType(derived()); -} - -/** \cpp11 \returns an expression of the coefficient-wise Complementary error - * function of *this. - * - * Example: \include Cwise_erfc.cpp - * Output: \verbinclude Cwise_erfc.out - * - * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, - * or float/double in non c++11 mode, the user has to provide implementations of erfc(T) for any scalar - * type T to be supported. - * - * \sa erf() - */ -EIGEN_DEVICE_FUNC -inline const ErfcReturnType -erfc() const -{ - return ErfcReturnType(derived()); -} - /** \returns an expression of the coefficient-wise inverse of *this. * * Example: \include Cwise_inverse.cpp @@ -538,3 +463,90 @@ operator!() const THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_OF_BOOL); return BooleanNotReturnType(derived()); } + + +// --- SpecialFunctions module --- + +typedef CwiseUnaryOp<internal::scalar_lgamma_op<Scalar>, const Derived> LgammaReturnType; +typedef CwiseUnaryOp<internal::scalar_digamma_op<Scalar>, const Derived> DigammaReturnType; +typedef CwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived> ErfReturnType; +typedef CwiseUnaryOp<internal::scalar_erfc_op<Scalar>, const Derived> ErfcReturnType; + +/** \cpp11 \returns an expression of the coefficient-wise ln(|gamma(*this)|). + * + * \specialfunctions_module + * + * Example: \include Cwise_lgamma.cpp + * Output: \verbinclude Cwise_lgamma.out + * + * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations of lgamma(T) for any scalar + * type T to be supported. + * + * \sa digamma() + */ +EIGEN_DEVICE_FUNC +inline const LgammaReturnType +lgamma() const +{ + return LgammaReturnType(derived()); +} + +/** \returns an expression of the coefficient-wise digamma (psi, derivative of lgamma). + * + * \specialfunctions_module + * + * \note This function supports only float and double scalar types. To support other scalar types, + * the user has to provide implementations of digamma(T) for any scalar + * type T to be supported. + * + * \sa Eigen::digamma(), Eigen::polygamma(), lgamma() + */ +EIGEN_DEVICE_FUNC +inline const DigammaReturnType +digamma() const +{ + return DigammaReturnType(derived()); +} + +/** \cpp11 \returns an expression of the coefficient-wise Gauss error + * function of *this. + * + * \specialfunctions_module + * + * Example: \include Cwise_erf.cpp + * Output: \verbinclude Cwise_erf.out + * + * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations of erf(T) for any scalar + * type T to be supported. + * + * \sa erfc() + */ +EIGEN_DEVICE_FUNC +inline const ErfReturnType +erf() const +{ + return ErfReturnType(derived()); +} + +/** \cpp11 \returns an expression of the coefficient-wise Complementary error + * function of *this. + * + * \specialfunctions_module + * + * Example: \include Cwise_erfc.cpp + * Output: \verbinclude Cwise_erfc.out + * + * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations of erfc(T) for any scalar + * type T to be supported. + * + * \sa erf() + */ +EIGEN_DEVICE_FUNC +inline const ErfcReturnType +erfc() const +{ + return ErfcReturnType(derived()); +} diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index 0c3673f89..cfc15a0f9 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -216,6 +216,7 @@ ALIASES = "only_for_vectors=This is only for vectors (either row- "lu_module=This is defined in the %LU module. \code #include <Eigen/LU> \endcode" \ "qr_module=This is defined in the %QR module. \code #include <Eigen/QR> \endcode" \ "svd_module=This is defined in the %SVD module. \code #include <Eigen/SVD> \endcode" \ + "specialfunctions_module=This is defined in the SpecialFunctions module. \code #include <Eigen/SpecialFunctions> \endcode" \ "label=\bug" \ "matrixworld=<a href='#matrixonly' style='color:green;text-decoration: none;'>*</a>" \ "arrayworld=<a href='#arrayonly' style='color:blue;text-decoration: none;'>*</a>" \ diff --git a/test/array.cpp b/test/array.cpp index 6347c8067..d734e604a 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -234,12 +234,7 @@ template<typename ArrayType> void array_real(const ArrayType& m) VERIFY_IS_APPROX(m1.sinh(), sinh(m1)); VERIFY_IS_APPROX(m1.cosh(), cosh(m1)); VERIFY_IS_APPROX(m1.tanh(), tanh(m1)); -#if EIGEN_HAS_C99_MATH - VERIFY_IS_APPROX(m1.lgamma(), lgamma(m1)); - VERIFY_IS_APPROX(m1.digamma(), digamma(m1)); - VERIFY_IS_APPROX(m1.erf(), erf(m1)); - VERIFY_IS_APPROX(m1.erfc(), erfc(m1)); -#endif // EIGEN_HAS_C99_MATH + VERIFY_IS_APPROX(m1.arg(), arg(m1)); VERIFY_IS_APPROX(m1.round(), round(m1)); VERIFY_IS_APPROX(m1.floor(), floor(m1)); @@ -313,88 +308,6 @@ template<typename ArrayType> void array_real(const ArrayType& m) m1 += ArrayType::Constant(rows,cols,Scalar(tiny)); VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse()); - - -#if EIGEN_HAS_C99_MATH - // check special functions (comparing against numpy implementation) - if (!NumTraits<Scalar>::IsComplex) - { - - { - // Test various propreties of igamma & igammac. These are normalized - // gamma integrals where - // igammac(a, x) = Gamma(a, x) / Gamma(a) - // igamma(a, x) = gamma(a, x) / Gamma(a) - // where Gamma and gamma are considered the standard unnormalized - // upper and lower incomplete gamma functions, respectively. - ArrayType a = m1.abs() + 2; - ArrayType x = m2.abs() + 2; - ArrayType zero = ArrayType::Zero(rows, cols); - ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0)); - ArrayType a_m1 = a - one; - ArrayType Gamma_a_x = Eigen::igammac(a, x) * a.lgamma().exp(); - ArrayType Gamma_a_m1_x = Eigen::igammac(a_m1, x) * a_m1.lgamma().exp(); - ArrayType gamma_a_x = Eigen::igamma(a, x) * a.lgamma().exp(); - ArrayType gamma_a_m1_x = Eigen::igamma(a_m1, x) * a_m1.lgamma().exp(); - - // Gamma(a, 0) == Gamma(a) - VERIFY_IS_APPROX(Eigen::igammac(a, zero), one); - - // Gamma(a, x) + gamma(a, x) == Gamma(a) - VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp()); - - // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x) - VERIFY_IS_APPROX(Gamma_a_x, (a - 1) * Gamma_a_m1_x + x.pow(a-1) * (-x).exp()); - - // gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x) - VERIFY_IS_APPROX(gamma_a_x, (a - 1) * gamma_a_m1_x - x.pow(a-1) * (-x).exp()); - } - - // Check exact values of igamma and igammac against a third party calculation. - Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; - Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; - - // location i*6+j corresponds to a_s[i], x_s[j]. - Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); - Scalar igamma_s[][6] = {{0.0, nan, nan, nan, nan, nan}, - {0.0, 0.6321205588285578, 0.7768698398515702, - 0.9816843611112658, 9.999500016666262e-05, 1.0}, - {0.0, 0.4275932955291202, 0.608374823728911, - 0.9539882943107686, 7.522076445089201e-07, 1.0}, - {0.0, 0.01898815687615381, 0.06564245437845008, - 0.5665298796332909, 4.166333347221828e-18, 1.0}, - {0.0, 0.9999780593618628, 0.9999899967080838, - 0.9999996219837988, 0.9991370418689945, 1.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.5042041932513908}}; - Scalar igammac_s[][6] = {{nan, nan, nan, nan, nan, nan}, - {1.0, 0.36787944117144233, 0.22313016014842982, - 0.018315638888734182, 0.9999000049998333, 0.0}, - {1.0, 0.5724067044708798, 0.3916251762710878, - 0.04601170568923136, 0.9999992477923555, 0.0}, - {1.0, 0.9810118431238462, 0.9343575456215499, - 0.4334701203667089, 1.0, 0.0}, - {1.0, 2.1940638138146658e-05, 1.0003291916285e-05, - 3.7801620118431334e-07, 0.0008629581310054535, - 0.0}, - {1.0, 1.0, 1.0, 1.0, 1.0, 0.49579580674813944}}; - for (int i = 0; i < 6; ++i) { - for (int j = 0; j < 6; ++j) { - if ((std::isnan)(igamma_s[i][j])) { - VERIFY((std::isnan)(numext::igamma(a_s[i], x_s[j]))); - } else { - VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]); - } - - if ((std::isnan)(igammac_s[i][j])) { - VERIFY((std::isnan)(numext::igammac(a_s[i], x_s[j]))); - } else { - VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][j]); - } - } - } - } -#endif // EIGEN_HAS_C99_MATH - // check inplace transpose m3 = m1; m3.transposeInPlace(); @@ -537,242 +450,8 @@ template<typename ArrayType> void min_max(const ArrayType& m) } -template<typename X, typename Y> -void verify_component_wise(const X& x, const Y& y) -{ - for(Index i=0; i<x.size(); ++i) - { - if((numext::isfinite)(y(i))) - VERIFY_IS_APPROX( x(i), y(i) ); - else if((numext::isnan)(y(i))) - VERIFY((numext::isnan)(x(i))); - else - VERIFY_IS_EQUAL( x(i), y(i) ); - } -} - -// check special functions (comparing against numpy implementation) -template<typename ArrayType> void array_special_functions() -{ - using std::abs; - using std::sqrt; - typedef typename ArrayType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; - - Scalar plusinf = std::numeric_limits<Scalar>::infinity(); - Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); - - // Check the zeta function against scipy.special.zeta - { - ArrayType x(7), q(7), res(7), ref(7); - x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9; - q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345; - ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan; - CALL_SUBTEST( verify_component_wise(ref, ref); ); - CALL_SUBTEST( res = x.zeta(q); verify_component_wise(res, ref); ); - CALL_SUBTEST( res = zeta(x,q); verify_component_wise(res, ref); ); - } - - // digamma - { - ArrayType x(7), res(7), ref(7); - x << 1, 1.5, 4, -10.5, 10000.5, 0, -1; - ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, plusinf, plusinf; - CALL_SUBTEST( verify_component_wise(ref, ref); ); - - CALL_SUBTEST( res = x.digamma(); verify_component_wise(res, ref); ); - CALL_SUBTEST( res = digamma(x); verify_component_wise(res, ref); ); - } - - -#if EIGEN_HAS_C99_MATH - { - ArrayType n(11), x(11), res(11), ref(11); - n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170; - x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64; - ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927; - CALL_SUBTEST( verify_component_wise(ref, ref); ); - - if(sizeof(RealScalar)>=8) { // double - // Reason for commented line: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 - // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res, ref); ); - CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res, ref); ); - } - else { - // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res.head(8), ref.head(8)); ); - CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res.head(8), ref.head(8)); ); - } - } -#endif - -#if EIGEN_HAS_C99_MATH - { - // Inputs and ground truth generated with scipy via: - // a = np.logspace(-3, 3, 5) - 1e-3 - // b = np.logspace(-3, 3, 5) - 1e-3 - // x = np.linspace(-0.1, 1.1, 5) - // (full_a, full_b, full_x) = np.vectorize(lambda a, b, x: (a, b, x))(*np.ix_(a, b, x)) - // full_a = full_a.flatten().tolist() # same for full_b, full_x - // v = scipy.special.betainc(full_a, full_b, full_x).flatten().tolist() - // - // Note in Eigen, we call betainc with arguments in the order (x, a, b). - ArrayType a(125); - ArrayType b(125); - ArrayType x(125); - ArrayType v(125); - ArrayType res(125); - - a << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, - 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, - 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, - 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, - 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, - 999.999, 999.999, 999.999; - - b << 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999, - 0.999, 0.999, 0.999, 0.999, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999, - 999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, - 0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999; - - x << -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, - 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, - 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, - 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, - -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, - 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, - 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, - 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, - 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, - 0.8, 1.1; - - v << nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, - nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, - nan, nan, nan, 0.47972119876364683, 0.5, 0.5202788012363533, nan, nan, - 0.9518683957740043, 0.9789663010413743, 0.9931729188073435, nan, nan, - 0.999995949033062, 0.9999999999993698, 0.9999999999999999, nan, nan, - 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, - nan, nan, nan, nan, nan, 0.006827081192655869, 0.0210336989586256, - 0.04813160422599567, nan, nan, 0.20014344256217678, 0.5000000000000001, - 0.7998565574378232, nan, nan, 0.9991401428435834, 0.999999999698403, - 0.9999999999999999, nan, nan, 0.9999999999999999, 0.9999999999999999, - 0.9999999999999999, nan, nan, nan, nan, nan, nan, nan, - 1.0646600232370887e-25, 6.301722877826246e-13, 4.050966937974938e-06, - nan, nan, 7.864342668429763e-23, 3.015969667594166e-10, - 0.0008598571564165444, nan, nan, 6.031987710123844e-08, - 0.5000000000000007, 0.9999999396801229, nan, nan, 0.9999999999999999, - 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan, nan, nan, - nan, 0.0, 7.029920380986636e-306, 2.2450728208591345e-101, nan, nan, - 0.0, 9.275871147869727e-302, 1.2232913026152827e-97, nan, nan, 0.0, - 3.0891393081932924e-252, 2.9303043666183996e-60, nan, nan, - 2.248913486879199e-196, 0.5000000000004947, 0.9999999999999999, nan; - - CALL_SUBTEST(res = betainc(a, b, x); - verify_component_wise(res, v);); - } - - // Test various properties of betainc - { - ArrayType m1 = ArrayType::Random(32); - ArrayType m2 = ArrayType::Random(32); - ArrayType m3 = ArrayType::Random(32); - ArrayType one = ArrayType::Constant(32, Scalar(1.0)); - const Scalar eps = std::numeric_limits<Scalar>::epsilon(); - ArrayType a = (m1 * 4.0).exp(); - ArrayType b = (m2 * 4.0).exp(); - ArrayType x = m3.abs(); - - // betainc(a, 1, x) == x**a - CALL_SUBTEST( - ArrayType test = betainc(a, one, x); - ArrayType expected = x.pow(a); - verify_component_wise(test, expected);); - - // betainc(1, b, x) == 1 - (1 - x)**b - CALL_SUBTEST( - ArrayType test = betainc(one, b, x); - ArrayType expected = one - (one - x).pow(b); - verify_component_wise(test, expected);); - - // betainc(a, b, x) == 1 - betainc(b, a, 1-x) - CALL_SUBTEST( - ArrayType test = betainc(a, b, x) + betainc(b, a, one - x); - ArrayType expected = one; - verify_component_wise(test, expected);); - - // betainc(a+1, b, x) = betainc(a, b, x) - x**a * (1 - x)**b / (a * beta(a, b)) - CALL_SUBTEST( - ArrayType num = x.pow(a) * (one - x).pow(b); - ArrayType denom = a * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); - // Add eps to rhs and lhs so that component-wise test doesn't result in - // nans when both outputs are zeros. - ArrayType expected = betainc(a, b, x) - num / denom + eps; - ArrayType test = betainc(a + one, b, x) + eps; - if (sizeof(Scalar) >= 8) { // double - verify_component_wise(test, expected); - } else { - // Reason for limited test: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 - verify_component_wise(test.head(8), expected.head(8)); - }); - - // betainc(a, b+1, x) = betainc(a, b, x) + x**a * (1 - x)**b / (b * beta(a, b)) - CALL_SUBTEST( - // Add eps to rhs and lhs so that component-wise test doesn't result in - // nans when both outputs are zeros. - ArrayType num = x.pow(a) * (one - x).pow(b); - ArrayType denom = b * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); - ArrayType expected = betainc(a, b, x) + num / denom + eps; - ArrayType test = betainc(a, b + one, x) + eps; - verify_component_wise(test, expected);); - } -#endif -} - void test_array() { -#ifndef EIGEN_HAS_C99_MATH - std::cerr << "WARNING: testing of special math functions disabled" << std::endl; -#endif - for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( array(Array<float, 1, 1>()) ); CALL_SUBTEST_2( array(Array22f()) ); @@ -812,7 +491,4 @@ void test_array() VERIFY((internal::is_same< internal::global_math_functions_filtering_base<Xpr>::type, ArrayBase<Xpr> >::value)); - - CALL_SUBTEST_7(array_special_functions<ArrayXf>()); - CALL_SUBTEST_7(array_special_functions<ArrayXd>()); } diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt index 6d0cf4f9d..7478b6b0d 100644 --- a/unsupported/Eigen/CMakeLists.txt +++ b/unsupported/Eigen/CMakeLists.txt @@ -17,6 +17,7 @@ set(Eigen_HEADERS Polynomials Skyline SparseExtra + SpecialFunctions Splines ) diff --git a/unsupported/Eigen/CXX11/Tensor b/unsupported/Eigen/CXX11/Tensor index 79bac2f67..f7b94cee1 100644 --- a/unsupported/Eigen/CXX11/Tensor +++ b/unsupported/Eigen/CXX11/Tensor @@ -15,6 +15,7 @@ #include <Eigen/src/Core/util/DisableStupidWarnings.h> +#include "../SpecialFunctions" #include "src/util/CXX11Meta.h" #include "src/util/MaxSizeVector.h" diff --git a/unsupported/Eigen/SpecialFunctions b/unsupported/Eigen/SpecialFunctions new file mode 100644 index 000000000..b4ed83b30 --- /dev/null +++ b/unsupported/Eigen/SpecialFunctions @@ -0,0 +1,57 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SPECIALFUNCTIONS_MODULE +#define EIGEN_SPECIALFUNCTIONS_MODULE + +#include "../../Eigen/Core" + +#include "../../Eigen/src/Core/util/DisableStupidWarnings.h" + +namespace Eigen { + +/** + * \defgroup SpecialFunctions_Module Special math functions module + * + * This module features additional coefficient-wise math functions available + * within the numext:: namespace for the scalar version, and as method and/or free + * functions of Array. Those include: + * + * - erf + * - erfc + * - lgamma + * - igamma + * - igammac + * - digamma + * - polygamma + * - zeta + * - betainc + * + * \code + * #include <unsupported/Eigen/SpecialFunctions> + * \endcode + */ +//@{ + +} + +#include "src/SpecialFunctions/SpecialFunctionsImpl.h" +#include "src/SpecialFunctions/SpecialFunctionsPacketMath.h" +#include "src/SpecialFunctions/SpecialFunctionsHalf.h" +#include "src/SpecialFunctions/SpecialFunctionsFunctors.h" +#include "src/SpecialFunctions/SpecialFunctionsArrayAPI.h" + +namespace Eigen { +//@} +} + + +#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h" + +#endif // EIGEN_SPECIALFUNCTIONS_MODULE diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index a7e8c7553..f42946793 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -11,5 +11,6 @@ ADD_SUBDIRECTORY(NumericalDiff) ADD_SUBDIRECTORY(Polynomials) ADD_SUBDIRECTORY(Skyline) ADD_SUBDIRECTORY(SparseExtra) +ADD_SUBDIRECTORY(SpecialFunctions) ADD_SUBDIRECTORY(KroneckerProduct) ADD_SUBDIRECTORY(Splines) diff --git a/unsupported/Eigen/src/SpecialFunctions/CMakeLists.txt b/unsupported/Eigen/src/SpecialFunctions/CMakeLists.txt new file mode 100644 index 000000000..00c7afd7b --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_SpecialFunctions_SRCS "*.h") + +INSTALL(FILES + ${Eigen_SpecialFunctions_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/SpecialFunctions COMPONENT Devel + ) diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h new file mode 100644 index 000000000..ed415db99 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h @@ -0,0 +1,124 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + + +#ifndef EIGEN_SPECIALFUNCTIONS_ARRAYAPI_H +#define EIGEN_SPECIALFUNCTIONS_ARRAYAPI_H + +namespace Eigen { + +/** \cpp11 \returns an expression of the coefficient-wise igamma(\a a, \a x) to the given arrays. + * + * This function computes the coefficient-wise incomplete gamma function. + * + * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations of igammac(T,T) for any scalar + * type T to be supported. + * + * \sa Eigen::igammac(), Eigen::lgamma() + */ +template<typename Derived,typename ExponentDerived> +inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_op<typename Derived::Scalar>, const Derived, const ExponentDerived> +igamma(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) +{ + return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_op<typename Derived::Scalar>, const Derived, const ExponentDerived>( + a.derived(), + x.derived() + ); +} + +/** \cpp11 \returns an expression of the coefficient-wise igammac(\a a, \a x) to the given arrays. + * + * This function computes the coefficient-wise complementary incomplete gamma function. + * + * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations of igammac(T,T) for any scalar + * type T to be supported. + * + * \sa Eigen::igamma(), Eigen::lgamma() + */ +template<typename Derived,typename ExponentDerived> +inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igammac_op<typename Derived::Scalar>, const Derived, const ExponentDerived> +igammac(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) +{ + return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igammac_op<typename Derived::Scalar>, const Derived, const ExponentDerived>( + a.derived(), + x.derived() + ); +} + +/** \cpp11 \returns an expression of the coefficient-wise polygamma(\a n, \a x) to the given arrays. + * + * It returns the \a n -th derivative of the digamma(psi) evaluated at \c x. + * + * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations of polygamma(T,T) for any scalar + * type T to be supported. + * + * \sa Eigen::digamma() + */ +// * \warning Be careful with the order of the parameters: x.polygamma(n) is equivalent to polygamma(n,x) +// * \sa ArrayBase::polygamma() +template<typename DerivedN,typename DerivedX> +inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_polygamma_op<typename DerivedX::Scalar>, const DerivedN, const DerivedX> +polygamma(const Eigen::ArrayBase<DerivedN>& n, const Eigen::ArrayBase<DerivedX>& x) +{ + return Eigen::CwiseBinaryOp<Eigen::internal::scalar_polygamma_op<typename DerivedX::Scalar>, const DerivedN, const DerivedX>( + n.derived(), + x.derived() + ); +} + +/** \cpp11 \returns an expression of the coefficient-wise betainc(\a x, \a a, \a b) to the given arrays. + * + * This function computes the regularized incomplete beta function (integral). + * + * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations of betainc(T,T,T) for any scalar + * type T to be supported. + * + * \sa Eigen::betainc(), Eigen::lgamma() + */ +template<typename ArgADerived, typename ArgBDerived, typename ArgXDerived> +inline const Eigen::CwiseTernaryOp<Eigen::internal::scalar_betainc_op<typename ArgXDerived::Scalar>, const ArgADerived, const ArgBDerived, const ArgXDerived> +betainc(const Eigen::ArrayBase<ArgADerived>& a, const Eigen::ArrayBase<ArgBDerived>& b, const Eigen::ArrayBase<ArgXDerived>& x) +{ + return Eigen::CwiseTernaryOp<Eigen::internal::scalar_betainc_op<typename ArgXDerived::Scalar>, const ArgADerived, const ArgBDerived, const ArgXDerived>( + a.derived(), + b.derived(), + x.derived() + ); +} + + +/** \returns an expression of the coefficient-wise zeta(\a x, \a q) to the given arrays. + * + * It returns the Riemann zeta function of two arguments \a x and \a q: + * + * \param x is the exposent, it must be > 1 + * \param q is the shift, it must be > 0 + * + * \note This function supports only float and double scalar types. To support other scalar types, the user has + * to provide implementations of zeta(T,T) for any scalar type T to be supported. + * + * \sa ArrayBase::zeta() + */ +template<typename DerivedX,typename DerivedQ> +inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_zeta_op<typename DerivedX::Scalar>, const DerivedX, const DerivedQ> +zeta(const Eigen::ArrayBase<DerivedX>& x, const Eigen::ArrayBase<DerivedQ>& q) +{ + return Eigen::CwiseBinaryOp<Eigen::internal::scalar_zeta_op<typename DerivedX::Scalar>, const DerivedX, const DerivedQ>( + x.derived(), + q.derived() + ); +} + +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_ARRAYAPI_H diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h new file mode 100644 index 000000000..d8f2363be --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h @@ -0,0 +1,236 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Eugene Brevdo <ebrevdo@gmail.com> +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SPECIALFUNCTIONS_FUNCTORS_H +#define EIGEN_SPECIALFUNCTIONS_FUNCTORS_H + +namespace Eigen { + +namespace internal { + + +/** \internal + * \brief Template functor to compute the incomplete gamma function igamma(a, x) + * + * \sa class CwiseBinaryOp, Cwise::igamma + */ +template<typename Scalar> struct scalar_igamma_op : binary_op_base<Scalar,Scalar> +{ + EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_op) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const { + using numext::igamma; return igamma(a, x); + } + template<typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const { + return internal::pigamma(a, x); + } +}; +template<typename Scalar> +struct functor_traits<scalar_igamma_op<Scalar> > { + enum { + // Guesstimate + Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasIGamma + }; +}; + + +/** \internal + * \brief Template functor to compute the complementary incomplete gamma function igammac(a, x) + * + * \sa class CwiseBinaryOp, Cwise::igammac + */ +template<typename Scalar> struct scalar_igammac_op : binary_op_base<Scalar,Scalar> +{ + EIGEN_EMPTY_STRUCT_CTOR(scalar_igammac_op) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const { + using numext::igammac; return igammac(a, x); + } + template<typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const + { + return internal::pigammac(a, x); + } +}; +template<typename Scalar> +struct functor_traits<scalar_igammac_op<Scalar> > { + enum { + // Guesstimate + Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasIGammac + }; +}; + + +/** \internal + * \brief Template functor to compute the incomplete beta integral betainc(a, b, x) + * + */ +template<typename Scalar> struct scalar_betainc_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_betainc_op) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& x, const Scalar& a, const Scalar& b) const { + using numext::betainc; return betainc(x, a, b); + } + template<typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& x, const Packet& a, const Packet& b) const + { + return internal::pbetainc(x, a, b); + } +}; +template<typename Scalar> +struct functor_traits<scalar_betainc_op<Scalar> > { + enum { + // Guesstimate + Cost = 400 * NumTraits<Scalar>::MulCost + 400 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasBetaInc + }; +}; + + +/** \internal + * \brief Template functor to compute the natural log of the absolute + * value of Gamma of a scalar + * \sa class CwiseUnaryOp, Cwise::lgamma() + */ +template<typename Scalar> struct scalar_lgamma_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_lgamma_op) + EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { + using numext::lgamma; return lgamma(a); + } + typedef typename packet_traits<Scalar>::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::plgamma(a); } +}; +template<typename Scalar> +struct functor_traits<scalar_lgamma_op<Scalar> > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasLGamma + }; +}; + +/** \internal + * \brief Template functor to compute psi, the derivative of lgamma of a scalar. + * \sa class CwiseUnaryOp, Cwise::digamma() + */ +template<typename Scalar> struct scalar_digamma_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_digamma_op) + EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { + using numext::digamma; return digamma(a); + } + typedef typename packet_traits<Scalar>::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::pdigamma(a); } +}; +template<typename Scalar> +struct functor_traits<scalar_digamma_op<Scalar> > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasDiGamma + }; +}; + +/** \internal + * \brief Template functor to compute the Riemann Zeta function of two arguments. + * \sa class CwiseUnaryOp, Cwise::zeta() + */ +template<typename Scalar> struct scalar_zeta_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_zeta_op) + EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& x, const Scalar& q) const { + using numext::zeta; return zeta(x, q); + } + typedef typename packet_traits<Scalar>::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& x, const Packet& q) const { return internal::pzeta(x, q); } +}; +template<typename Scalar> +struct functor_traits<scalar_zeta_op<Scalar> > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasZeta + }; +}; + +/** \internal + * \brief Template functor to compute the polygamma function. + * \sa class CwiseUnaryOp, Cwise::polygamma() + */ +template<typename Scalar> struct scalar_polygamma_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_polygamma_op) + EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& n, const Scalar& x) const { + using numext::polygamma; return polygamma(n, x); + } + typedef typename packet_traits<Scalar>::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& n, const Packet& x) const { return internal::ppolygamma(n, x); } +}; +template<typename Scalar> +struct functor_traits<scalar_polygamma_op<Scalar> > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasPolygamma + }; +}; + +/** \internal + * \brief Template functor to compute the Gauss error function of a + * scalar + * \sa class CwiseUnaryOp, Cwise::erf() + */ +template<typename Scalar> struct scalar_erf_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_erf_op) + EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { + using numext::erf; return erf(a); + } + typedef typename packet_traits<Scalar>::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perf(a); } +}; +template<typename Scalar> +struct functor_traits<scalar_erf_op<Scalar> > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasErf + }; +}; + +/** \internal + * \brief Template functor to compute the Complementary Error Function + * of a scalar + * \sa class CwiseUnaryOp, Cwise::erfc() + */ +template<typename Scalar> struct scalar_erfc_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_erfc_op) + EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const { + using numext::erfc; return erfc(a); + } + typedef typename packet_traits<Scalar>::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perfc(a); } +}; +template<typename Scalar> +struct functor_traits<scalar_erfc_op<Scalar> > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasErfc + }; +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_FUNCTORS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h new file mode 100644 index 000000000..553bcda6a --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h @@ -0,0 +1,47 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SPECIALFUNCTIONS_HALF_H +#define EIGEN_SPECIALFUNCTIONS_HALF_H + +namespace Eigen { +namespace numext { + +#if EIGEN_HAS_C99_MATH +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half lgamma(const Eigen::half& a) { + return Eigen::half(Eigen::numext::lgamma(static_cast<float>(a))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half digamma(const Eigen::half& a) { + return Eigen::half(Eigen::numext::digamma(static_cast<float>(a))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half zeta(const Eigen::half& x, const Eigen::half& q) { + return Eigen::half(Eigen::numext::zeta(static_cast<float>(x), static_cast<float>(q))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half polygamma(const Eigen::half& n, const Eigen::half& x) { + return Eigen::half(Eigen::numext::polygamma(static_cast<float>(n), static_cast<float>(x))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erf(const Eigen::half& a) { + return Eigen::half(Eigen::numext::erf(static_cast<float>(a))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erfc(const Eigen::half& a) { + return Eigen::half(Eigen::numext::erfc(static_cast<float>(a))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igamma(const Eigen::half& a, const Eigen::half& x) { + return Eigen::half(Eigen::numext::igamma(static_cast<float>(a), static_cast<float>(x))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igammac(const Eigen::half& a, const Eigen::half& x) { + return Eigen::half(Eigen::numext::igammac(static_cast<float>(a), static_cast<float>(x))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half betainc(const Eigen::half& a, const Eigen::half& b, const Eigen::half& x) { + return Eigen::half(Eigen::numext::betainc(static_cast<float>(a), static_cast<float>(b), static_cast<float>(x))); +} +#endif + +} // end namespace numext +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_HALF_H diff --git a/Eigen/src/Core/SpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index a657cb854..52619fc0c 100644 --- a/Eigen/src/Core/SpecialFunctions.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -1382,7 +1382,7 @@ struct betainc_helper<double> { */ t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) - lgamma_impl<double>::run(b) + u + numext::log(s); - return s = exp(t); + return s = numext::exp(t); } }; diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h new file mode 100644 index 000000000..46d60d323 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h @@ -0,0 +1,58 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SPECIALFUNCTIONS_PACKETMATH_H +#define EIGEN_SPECIALFUNCTIONS_PACKETMATH_H + +namespace Eigen { + +namespace internal { + +/** \internal \returns the ln(|gamma(\a a)|) (coeff-wise) */ +template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +Packet plgamma(const Packet& a) { using numext::lgamma; return lgamma(a); } + +/** \internal \returns the derivative of lgamma, psi(\a a) (coeff-wise) */ +template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +Packet pdigamma(const Packet& a) { using numext::digamma; return digamma(a); } + +/** \internal \returns the zeta function of two arguments (coeff-wise) */ +template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +Packet pzeta(const Packet& x, const Packet& q) { using numext::zeta; return zeta(x, q); } + +/** \internal \returns the polygamma function (coeff-wise) */ +template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +Packet ppolygamma(const Packet& n, const Packet& x) { using numext::polygamma; return polygamma(n, x); } + +/** \internal \returns the erf(\a a) (coeff-wise) */ +template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +Packet perf(const Packet& a) { using numext::erf; return erf(a); } + +/** \internal \returns the erfc(\a a) (coeff-wise) */ +template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +Packet perfc(const Packet& a) { using numext::erfc; return erfc(a); } + +/** \internal \returns the incomplete gamma function igamma(\a a, \a x) */ +template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +Packet pigamma(const Packet& a, const Packet& x) { using numext::igamma; return igamma(a, x); } + +/** \internal \returns the complementary incomplete gamma function igammac(\a a, \a x) */ +template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +Packet pigammac(const Packet& a, const Packet& x) { using numext::igammac; return igammac(a, x); } + +/** \internal \returns the complementary incomplete gamma function betainc(\a a, \a b, \a x) */ +template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +Packet pbetainc(const Packet& a, const Packet& b,const Packet& x) { using numext::betainc; return betainc(a, b, x); } + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_PACKETMATH_H + diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index c9a70d7a7..5137b51cf 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -109,6 +109,7 @@ ei_add_test(gmres) ei_add_test(minres) ei_add_test(levenberg_marquardt) ei_add_test(kronecker_product) +ei_add_test(special_functions) # TODO: The following test names are prefixed with the cxx11 string, since historically # the tests depended on c++11. This isn't the case anymore so we ought to rename them. diff --git a/unsupported/test/special_functions.cpp b/unsupported/test/special_functions.cpp new file mode 100644 index 000000000..057fb3e92 --- /dev/null +++ b/unsupported/test/special_functions.cpp @@ -0,0 +1,345 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include "../Eigen/SpecialFunctions" + +template<typename X, typename Y> +void verify_component_wise(const X& x, const Y& y) +{ + for(Index i=0; i<x.size(); ++i) + { + if((numext::isfinite)(y(i))) + VERIFY_IS_APPROX( x(i), y(i) ); + else if((numext::isnan)(y(i))) + VERIFY((numext::isnan)(x(i))); + else + VERIFY_IS_EQUAL( x(i), y(i) ); + } +} + +template<typename ArrayType> void array_special_functions() +{ + using std::abs; + using std::sqrt; + typedef typename ArrayType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + Scalar plusinf = std::numeric_limits<Scalar>::infinity(); + Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); + + Index rows = internal::random<Index>(1,30); + Index cols = 1; + + // API + { + ArrayType m1 = ArrayType::Random(rows,cols); +#if EIGEN_HAS_C99_MATH + VERIFY_IS_APPROX(m1.lgamma(), lgamma(m1)); + VERIFY_IS_APPROX(m1.digamma(), digamma(m1)); + VERIFY_IS_APPROX(m1.erf(), erf(m1)); + VERIFY_IS_APPROX(m1.erfc(), erfc(m1)); +#endif // EIGEN_HAS_C99_MATH + } + + +#if EIGEN_HAS_C99_MATH + // check special functions (comparing against numpy implementation) + if (!NumTraits<Scalar>::IsComplex) + { + + { + ArrayType m1 = ArrayType::Random(rows,cols); + ArrayType m2 = ArrayType::Random(rows,cols); + + // Test various propreties of igamma & igammac. These are normalized + // gamma integrals where + // igammac(a, x) = Gamma(a, x) / Gamma(a) + // igamma(a, x) = gamma(a, x) / Gamma(a) + // where Gamma and gamma are considered the standard unnormalized + // upper and lower incomplete gamma functions, respectively. + ArrayType a = m1.abs() + 2; + ArrayType x = m2.abs() + 2; + ArrayType zero = ArrayType::Zero(rows, cols); + ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0)); + ArrayType a_m1 = a - one; + ArrayType Gamma_a_x = Eigen::igammac(a, x) * a.lgamma().exp(); + ArrayType Gamma_a_m1_x = Eigen::igammac(a_m1, x) * a_m1.lgamma().exp(); + ArrayType gamma_a_x = Eigen::igamma(a, x) * a.lgamma().exp(); + ArrayType gamma_a_m1_x = Eigen::igamma(a_m1, x) * a_m1.lgamma().exp(); + + // Gamma(a, 0) == Gamma(a) + VERIFY_IS_APPROX(Eigen::igammac(a, zero), one); + + // Gamma(a, x) + gamma(a, x) == Gamma(a) + VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp()); + + // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x) + VERIFY_IS_APPROX(Gamma_a_x, (a - 1) * Gamma_a_m1_x + x.pow(a-1) * (-x).exp()); + + // gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x) + VERIFY_IS_APPROX(gamma_a_x, (a - 1) * gamma_a_m1_x - x.pow(a-1) * (-x).exp()); + } + + { + // Check exact values of igamma and igammac against a third party calculation. + Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; + Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; + + // location i*6+j corresponds to a_s[i], x_s[j]. + Scalar igamma_s[][6] = {{0.0, nan, nan, nan, nan, nan}, + {0.0, 0.6321205588285578, 0.7768698398515702, + 0.9816843611112658, 9.999500016666262e-05, 1.0}, + {0.0, 0.4275932955291202, 0.608374823728911, + 0.9539882943107686, 7.522076445089201e-07, 1.0}, + {0.0, 0.01898815687615381, 0.06564245437845008, + 0.5665298796332909, 4.166333347221828e-18, 1.0}, + {0.0, 0.9999780593618628, 0.9999899967080838, + 0.9999996219837988, 0.9991370418689945, 1.0}, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.5042041932513908}}; + Scalar igammac_s[][6] = {{nan, nan, nan, nan, nan, nan}, + {1.0, 0.36787944117144233, 0.22313016014842982, + 0.018315638888734182, 0.9999000049998333, 0.0}, + {1.0, 0.5724067044708798, 0.3916251762710878, + 0.04601170568923136, 0.9999992477923555, 0.0}, + {1.0, 0.9810118431238462, 0.9343575456215499, + 0.4334701203667089, 1.0, 0.0}, + {1.0, 2.1940638138146658e-05, 1.0003291916285e-05, + 3.7801620118431334e-07, 0.0008629581310054535, + 0.0}, + {1.0, 1.0, 1.0, 1.0, 1.0, 0.49579580674813944}}; + for (int i = 0; i < 6; ++i) { + for (int j = 0; j < 6; ++j) { + if ((std::isnan)(igamma_s[i][j])) { + VERIFY((std::isnan)(numext::igamma(a_s[i], x_s[j]))); + } else { + VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]); + } + + if ((std::isnan)(igammac_s[i][j])) { + VERIFY((std::isnan)(numext::igammac(a_s[i], x_s[j]))); + } else { + VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][j]); + } + } + } + } + } +#endif // EIGEN_HAS_C99_MATH + + // Check the zeta function against scipy.special.zeta + { + ArrayType x(7), q(7), res(7), ref(7); + x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9; + q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345; + ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan; + CALL_SUBTEST( verify_component_wise(ref, ref); ); + CALL_SUBTEST( res = x.zeta(q); verify_component_wise(res, ref); ); + CALL_SUBTEST( res = zeta(x,q); verify_component_wise(res, ref); ); + } + + // digamma + { + ArrayType x(7), res(7), ref(7); + x << 1, 1.5, 4, -10.5, 10000.5, 0, -1; + ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, plusinf, plusinf; + CALL_SUBTEST( verify_component_wise(ref, ref); ); + + CALL_SUBTEST( res = x.digamma(); verify_component_wise(res, ref); ); + CALL_SUBTEST( res = digamma(x); verify_component_wise(res, ref); ); + } + + +#if EIGEN_HAS_C99_MATH + { + ArrayType n(11), x(11), res(11), ref(11); + n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170; + x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64; + ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927; + CALL_SUBTEST( verify_component_wise(ref, ref); ); + + if(sizeof(RealScalar)>=8) { // double + // Reason for commented line: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 + // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res, ref); ); + CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res, ref); ); + } + else { + // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res.head(8), ref.head(8)); ); + CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res.head(8), ref.head(8)); ); + } + } +#endif + +#if EIGEN_HAS_C99_MATH + { + // Inputs and ground truth generated with scipy via: + // a = np.logspace(-3, 3, 5) - 1e-3 + // b = np.logspace(-3, 3, 5) - 1e-3 + // x = np.linspace(-0.1, 1.1, 5) + // (full_a, full_b, full_x) = np.vectorize(lambda a, b, x: (a, b, x))(*np.ix_(a, b, x)) + // full_a = full_a.flatten().tolist() # same for full_b, full_x + // v = scipy.special.betainc(full_a, full_b, full_x).flatten().tolist() + // + // Note in Eigen, we call betainc with arguments in the order (x, a, b). + ArrayType a(125); + ArrayType b(125); + ArrayType x(125); + ArrayType v(125); + ArrayType res(125); + + a << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, + 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, + 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, + 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, + 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, + 999.999, 999.999, 999.999; + + b << 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999, + 0.999, 0.999, 0.999, 0.999, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999, + 999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, + 0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999; + + x << -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, + 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, + 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, + 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, + -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, + 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, + 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, + 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, + 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, + 0.8, 1.1; + + v << nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, + nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, + nan, nan, nan, 0.47972119876364683, 0.5, 0.5202788012363533, nan, nan, + 0.9518683957740043, 0.9789663010413743, 0.9931729188073435, nan, nan, + 0.999995949033062, 0.9999999999993698, 0.9999999999999999, nan, nan, + 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, + nan, nan, nan, nan, nan, 0.006827081192655869, 0.0210336989586256, + 0.04813160422599567, nan, nan, 0.20014344256217678, 0.5000000000000001, + 0.7998565574378232, nan, nan, 0.9991401428435834, 0.999999999698403, + 0.9999999999999999, nan, nan, 0.9999999999999999, 0.9999999999999999, + 0.9999999999999999, nan, nan, nan, nan, nan, nan, nan, + 1.0646600232370887e-25, 6.301722877826246e-13, 4.050966937974938e-06, + nan, nan, 7.864342668429763e-23, 3.015969667594166e-10, + 0.0008598571564165444, nan, nan, 6.031987710123844e-08, + 0.5000000000000007, 0.9999999396801229, nan, nan, 0.9999999999999999, + 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan, nan, nan, + nan, 0.0, 7.029920380986636e-306, 2.2450728208591345e-101, nan, nan, + 0.0, 9.275871147869727e-302, 1.2232913026152827e-97, nan, nan, 0.0, + 3.0891393081932924e-252, 2.9303043666183996e-60, nan, nan, + 2.248913486879199e-196, 0.5000000000004947, 0.9999999999999999, nan; + + CALL_SUBTEST(res = betainc(a, b, x); + verify_component_wise(res, v);); + } + + // Test various properties of betainc + { + ArrayType m1 = ArrayType::Random(32); + ArrayType m2 = ArrayType::Random(32); + ArrayType m3 = ArrayType::Random(32); + ArrayType one = ArrayType::Constant(32, Scalar(1.0)); + const Scalar eps = std::numeric_limits<Scalar>::epsilon(); + ArrayType a = (m1 * 4.0).exp(); + ArrayType b = (m2 * 4.0).exp(); + ArrayType x = m3.abs(); + + // betainc(a, 1, x) == x**a + CALL_SUBTEST( + ArrayType test = betainc(a, one, x); + ArrayType expected = x.pow(a); + verify_component_wise(test, expected);); + + // betainc(1, b, x) == 1 - (1 - x)**b + CALL_SUBTEST( + ArrayType test = betainc(one, b, x); + ArrayType expected = one - (one - x).pow(b); + verify_component_wise(test, expected);); + + // betainc(a, b, x) == 1 - betainc(b, a, 1-x) + CALL_SUBTEST( + ArrayType test = betainc(a, b, x) + betainc(b, a, one - x); + ArrayType expected = one; + verify_component_wise(test, expected);); + + // betainc(a+1, b, x) = betainc(a, b, x) - x**a * (1 - x)**b / (a * beta(a, b)) + CALL_SUBTEST( + ArrayType num = x.pow(a) * (one - x).pow(b); + ArrayType denom = a * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); + // Add eps to rhs and lhs so that component-wise test doesn't result in + // nans when both outputs are zeros. + ArrayType expected = betainc(a, b, x) - num / denom + eps; + ArrayType test = betainc(a + one, b, x) + eps; + if (sizeof(Scalar) >= 8) { // double + verify_component_wise(test, expected); + } else { + // Reason for limited test: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 + verify_component_wise(test.head(8), expected.head(8)); + }); + + // betainc(a, b+1, x) = betainc(a, b, x) + x**a * (1 - x)**b / (b * beta(a, b)) + CALL_SUBTEST( + // Add eps to rhs and lhs so that component-wise test doesn't result in + // nans when both outputs are zeros. + ArrayType num = x.pow(a) * (one - x).pow(b); + ArrayType denom = b * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); + ArrayType expected = betainc(a, b, x) + num / denom + eps; + ArrayType test = betainc(a, b + one, x) + eps; + verify_component_wise(test, expected);); + } +#endif +} + +void test_special_functions() +{ + CALL_SUBTEST_1(array_special_functions<ArrayXf>()); + CALL_SUBTEST_2(array_special_functions<ArrayXd>()); +} |