diff options
author | Eugene Brevdo <ebrevdo@gmail.com> | 2016-03-03 19:39:41 -0800 |
---|---|---|
committer | Eugene Brevdo <ebrevdo@gmail.com> | 2016-03-03 19:39:41 -0800 |
commit | 7ea35bfa1c0b4950feae65d49c0e6f2cbf3691d9 (patch) | |
tree | 9e5471a21d70b96b8084b9716f1abcf8614521dd /Eigen/src | |
parent | ab3dc0b0fe64c34fab110f15914b0b9fcc0329da (diff) |
Initial implementation of igamma and igammac.
Diffstat (limited to 'Eigen/src')
-rw-r--r-- | Eigen/src/Core/GenericPacketMath.h | 10 | ||||
-rw-r--r-- | Eigen/src/Core/GlobalFunctions.h | 30 | ||||
-rw-r--r-- | Eigen/src/Core/NumTraits.h | 5 | ||||
-rw-r--r-- | Eigen/src/Core/SpecialFunctions.h | 291 | ||||
-rw-r--r-- | Eigen/src/Core/arch/CUDA/MathFunctions.h | 18 | ||||
-rw-r--r-- | Eigen/src/Core/arch/CUDA/PacketMath.h | 4 | ||||
-rw-r--r-- | Eigen/src/Core/functors/BinaryFunctors.h | 49 | ||||
-rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/util/Meta.h | 5 |
9 files changed, 413 insertions, 1 deletions
diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h index 02882bdea..ead0253df 100644 --- a/Eigen/src/Core/GenericPacketMath.h +++ b/Eigen/src/Core/GenericPacketMath.h @@ -78,6 +78,8 @@ struct default_packet_traits HasDiGamma = 0, HasErf = 0, HasErfc = 0, + HasIGamma = 0, + HasIGammac = 0, HasRound = 0, HasFloor = 0, @@ -457,6 +459,14 @@ Packet perf(const Packet& a) { using numext::erf; return erf(a); } 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_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +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_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS +Packet pigammac(const Packet& a, const Packet& x) { using numext::igammac; return igammac(a, 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 396da8e71..7df0fdda9 100644 --- a/Eigen/src/Core/GlobalFunctions.h +++ b/Eigen/src/Core/GlobalFunctions.h @@ -129,6 +129,36 @@ namespace Eigen ); } + /** \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. + * + */ + 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() + ); + } + + /** \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. + * + */ + 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() + ); + } + namespace internal { EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(real,scalar_real_op) diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h index 6a596bb7d..7ddb4a867 100644 --- a/Eigen/src/Core/NumTraits.h +++ b/Eigen/src/Core/NumTraits.h @@ -95,6 +95,11 @@ template<typename T> struct GenericNumTraits static inline T infinity() { return numext::numeric_limits<T>::infinity(); } + + EIGEN_DEVICE_FUNC + static inline T quiet_NaN() { + return numext::numeric_limits<T>::quiet_NaN(); + } }; template<typename T> struct NumTraits : GenericNumTraits<T> diff --git a/Eigen/src/Core/SpecialFunctions.h b/Eigen/src/Core/SpecialFunctions.h index b02ad9a1f..ff2146afc 100644 --- a/Eigen/src/Core/SpecialFunctions.h +++ b/Eigen/src/Core/SpecialFunctions.h @@ -283,7 +283,7 @@ struct digamma_impl { Scalar p, q, nz, s, w, y; bool negative; - const Scalar maxnum = numext::numeric_limits<Scalar>::infinity(); + const Scalar maxnum = NumTraits<Scalar>::infinity(); const Scalar m_pi = 3.14159265358979323846; negative = 0; @@ -401,6 +401,282 @@ struct erfc_impl<double> { }; #endif // EIGEN_HAS_C99_MATH +/**************************************************************************** + * Implementation of igammac (complemented incomplete gamma integral) * + ****************************************************************************/ + +template <typename Scalar> +struct igammac_retval { + typedef Scalar type; +}; + +#ifndef EIGEN_HAS_C99_MATH + +template <typename Scalar> +struct igammac_impl { + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +#else + +template <typename Scalar> struct igamma_impl; // predeclare igamma_impl + +template <typename Scalar> +struct igammac_impl { + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + /* igamc() + * + * Incomplete gamma integral (modified for Eigen) + * + * + * + * SYNOPSIS: + * + * double a, x, y, igamc(); + * + * y = igamc( a, x ); + * + * DESCRIPTION: + * + * The function is defined by + * + * + * igamc(a,x) = 1 - igam(a,x) + * + * inf. + * - + * 1 | | -t a-1 + * = ----- | e t dt. + * - | | + * | (a) - + * x + * + * + * In this implementation both arguments must be positive. + * The integral is evaluated by either a power series or + * continued fraction expansion, depending on the relative + * values of a and x. + * + * ACCURACY (float): + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 30000 7.8e-6 5.9e-7 + * + * + * ACCURACY (double): + * + * Tested at random a, x. + * a x Relative error: + * arithmetic domain domain # trials peak rms + * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15 + * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15 + * + */ + /* + Cephes Math Library Release 2.2: June, 1992 + Copyright 1985, 1987, 1992 by Stephen L. Moshier + Direct inquiries to 30 Frost Street, Cambridge, MA 02140 + */ + const Scalar zero = 0; + const Scalar one = 1; + const Scalar two = 2; + const Scalar machep = NumTraits<Scalar>::epsilon(); + const Scalar maxlog = ::log(NumTraits<Scalar>::highest()); + const Scalar big = one / machep; + + Scalar ans, ax, c, yc, r, t, y, z; + Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; + + if ((x <= zero) || ( a <= zero)) { + return one; + } + + if ((x < one) || (x < a)) { + return (one - igamma_impl<Scalar>::run(a, x)); + } + + ax = a * ::log(x) - x - lgamma_impl<Scalar>::run(a); + if( ax < -maxlog ) { // underflow + return zero; + } + ax = ::exp(ax); + + // continued fraction + y = one - a; + z = x + y + one; + c = zero; + pkm2 = one; + qkm2 = x; + pkm1 = x + one; + qkm1 = z * x; + ans = pkm1/qkm1; + + do { + c += one; + y += one; + z += two; + yc = y * c; + pk = pkm1 * z - pkm2 * yc; + qk = qkm1 * z - qkm2 * yc; + if( qk != zero ) { + r = pk/qk; + t = ::abs( (ans - r)/r ); + ans = r; + } else { + t = one; + } + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + if (::abs(pk) > big) { + pkm2 *= machep; + pkm1 *= machep; + qkm2 *= machep; + qkm1 *= machep; + } + } while( t > machep ); + + return ( ans * ax ); + } +}; + +#endif // EIGEN_HAS_C99_MATH + +/**************************************************************************** + * Implementation of igamma (incomplete gamma integral) * + ****************************************************************************/ + +template <typename Scalar> +struct igamma_retval { + typedef Scalar type; +}; + +#ifndef EIGEN_HAS_C99_MATH + +template <typename Scalar> +struct igamma_impl { + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +#else + +template <typename Scalar> +struct igamma_impl { + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + /* igam() + * Incomplete gamma integral + * + * + * + * SYNOPSIS: + * + * double a, x, y, igam(); + * + * y = igam( a, x ); + * + * DESCRIPTION: + * + * The function is defined by + * + * x + * - + * 1 | | -t a-1 + * igam(a,x) = ----- | e t dt. + * - | | + * | (a) - + * 0 + * + * + * In this implementation both arguments must be positive. + * The integral is evaluated by either a power series or + * continued fraction expansion, depending on the relative + * values of a and x. + * + * ACCURACY (double): + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 200000 3.6e-14 2.9e-15 + * IEEE 0,100 300000 9.9e-14 1.5e-14 + * + * + * ACCURACY (float): + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 20000 7.8e-6 5.9e-7 + * + */ + /* + Cephes Math Library Release 2.2: June, 1992 + Copyright 1985, 1987, 1992 by Stephen L. Moshier + Direct inquiries to 30 Frost Street, Cambridge, MA 02140 + */ + + + /* left tail of incomplete gamma function: + * + * inf. k + * a -x - x + * x e > ---------- + * - - + * k=0 | (a+k+1) + * + */ + const Scalar zero = 0; + const Scalar one = 1; + const Scalar machep = NumTraits<Scalar>::epsilon(); + const Scalar maxlog = ::log(NumTraits<Scalar>::highest()); + + double ans, ax, c, r; + + if( (x <= zero) || ( a <= zero) ) { + return zero; + } + + if( (x > one) && (x > a ) ) { + return (one - igammac_impl<Scalar>::run(a,x)); + } + + /* Compute x**a * exp(-x) / gamma(a) */ + ax = a * ::log(x) - x - lgamma_impl<Scalar>::run(a); + if( ax < -maxlog ) { + // underflow + return zero; + } + ax = ::exp(ax); + + /* power series */ + r = a; + c = one; + ans = one; + + do { + r += one; + c *= x/r; + ans += c; + } while( c/ans > machep ); + + return( ans * ax/a ); + } +}; + +#endif // EIGEN_HAS_C99_MATH + } // end namespace internal namespace numext { @@ -429,8 +705,21 @@ EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar) return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x); } +template <typename Scalar> +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma, Scalar) + igamma(const Scalar& a, const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(igamma, Scalar)::run(a, x); +} + +template <typename Scalar> +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igammac, Scalar) + igammac(const Scalar& a, const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(igammac, Scalar)::run(a, x); +} + } // end namespace numext + } // end namespace Eigen #endif // EIGEN_SPECIAL_FUNCTIONS_H diff --git a/Eigen/src/Core/arch/CUDA/MathFunctions.h b/Eigen/src/Core/arch/CUDA/MathFunctions.h index a2c06a817..6e84d3af8 100644 --- a/Eigen/src/Core/arch/CUDA/MathFunctions.h +++ b/Eigen/src/Core/arch/CUDA/MathFunctions.h @@ -116,6 +116,24 @@ double2 perfc<double2>(const double2& a) return make_double2(erfc(a.x), erfc(a.y)); } +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 pigamma<float4>(const float4& a, const float4& x) +{ + using numext::pigamma; + return make_float4( + pigamma(a.x, x.x), + pigamma(a.y, x.y), + pigamma(a.z, x.z), + pigamma(a.w, x.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 pigammac<double2>(const double2& a, const double& x) +{ + using numext::pigammac; + return make_double2(pigammac(a.x, x.x), pigammac(a.y, x.y)); +} + #endif diff --git a/Eigen/src/Core/arch/CUDA/PacketMath.h b/Eigen/src/Core/arch/CUDA/PacketMath.h index d3d9f910e..d2563030b 100644 --- a/Eigen/src/Core/arch/CUDA/PacketMath.h +++ b/Eigen/src/Core/arch/CUDA/PacketMath.h @@ -43,6 +43,8 @@ template<> struct packet_traits<float> : default_packet_traits HasDiGamma = 1, HasErf = 1, HasErfc = 1, + HasIgamma = 1, + HasIGammac = 1, HasBlend = 0, }; @@ -67,6 +69,8 @@ template<> struct packet_traits<double> : default_packet_traits HasDiGamma = 1, HasErf = 1, HasErfc = 1, + HasIGamma = 1, + HasIGammac = 1, HasBlend = 0, }; diff --git a/Eigen/src/Core/functors/BinaryFunctors.h b/Eigen/src/Core/functors/BinaryFunctors.h index 4962d625c..5cdfff845 100644 --- a/Eigen/src/Core/functors/BinaryFunctors.h +++ b/Eigen/src/Core/functors/BinaryFunctors.h @@ -337,6 +337,55 @@ template<> struct functor_traits<scalar_boolean_or_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 { + 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::pigammac(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 { + 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/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index f09632375..a102e5457 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -206,6 +206,8 @@ template<typename Scalar> struct scalar_add_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 LhsScalar,typename RhsScalar=LhsScalar> struct scalar_product_op; template<typename LhsScalar,typename RhsScalar> struct scalar_multiple2_op; diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 6b35179f2..24e8a6d8a 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -148,6 +148,7 @@ template<typename T> struct numeric_limits static T (max)() { assert(false && "Highest not supported for this type"); } static T (min)() { assert(false && "Lowest not supported for this type"); } static T infinity() { assert(false && "Infinity not supported for this type"); } + static T quiet_NaN() { assert(false && "quiet_NaN not supported for this type"); } }; template<> struct numeric_limits<float> { @@ -159,6 +160,8 @@ template<> struct numeric_limits<float> static float (min)() { return FLT_MIN; } EIGEN_DEVICE_FUNC static float infinity() { return CUDART_INF_F; } + EIGEN_DEVICE_FUNC + static float quiet_NaN() { return CUDART_NAN_F; } }; template<> struct numeric_limits<double> { @@ -170,6 +173,8 @@ template<> struct numeric_limits<double> static double (min)() { return DBL_MIN; } EIGEN_DEVICE_FUNC static double infinity() { return CUDART_INF; } + EIGEN_DEVICE_FUNC + static double quiet_NaN() { return CUDART_NAN; } }; template<> struct numeric_limits<int> { |