mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-27 07:29:52 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add generic PacketMath implementation of the Error Function (erf).

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2019-09-19 12:48:30 -07:00
parent e02d429637
commit 6de5ed08d8
20 changed files with 265 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/Core/GenericPacketMath.h
View File

@ -531,6 +531,10 @@ Packet pcosh(const Packet& a) { EIGEN_USING_STD_MATH(cosh); return cosh(a); }
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet ptanh(const Packet& a) { EIGEN_USING_STD_MATH(tanh); return tanh(a); }
/** \internal \returns the error function of \a a (coeff-wise). */
template <typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet perf(const Packet& a) { return numext::erf(a); }
/** \internal \returns the exp of \a a (coeff-wise) */
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet pexp(const Packet& a) { EIGEN_USING_STD_MATH(exp); return exp(a); }

26
Eigen/src/Core/MathFunctions.h
View File

@ -891,7 +891,7 @@ template<typename T> EIGEN_DEVICE_FUNC bool isnan_impl(const std::complex<T>& x)
template<typename T> EIGEN_DEVICE_FUNC bool isinf_impl(const std::complex<T>& x);
template<typename T> T generic_fast_tanh_float(const T& a_x);
template<typename T> T generic_fast_erf_float(const T& a_x);
} // end namespace internal
/****************************************************************************
@ -1579,6 +1579,30 @@ template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double tanh(const double &x) { return ::tanh(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T erf(const T &x) {
EIGEN_USING_STD_MATH(tanh);
return erf(x);
}
#if (!defined(EIGEN_GPUCC)) && EIGEN_FAST_MATH && !defined(SYCL_DEVICE_ONLY)
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float erf(float x) { return internal::generic_fast_erf_float(x); }
#endif
#if defined(SYCL_DEVICE_ONLY)
SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(erf, erf)
#endif
#if !EIGEN_HAS_CXX11 || defined(EIGEN_GPUCC)
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float erf(const float &x) { return ::erff(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double erf(const double &x) { return ::erf(x); }
#endif
template <typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T fmod(const T& a, const T& b) {

52
Eigen/src/Core/MathFunctionsImpl.h
View File

@ -66,6 +66,58 @@ T generic_fast_tanh_float(const T& a_x)
return pdiv(p, q);
}
/** \internal \returns the error function of \a a (coeff-wise)
Doesn't do anything fancy, just a 13/8-degree rational interpolant which
is accurate up to a couple of ulp in the range [-4, 4], outside of which
fl(erf(x)) = +/-1.
This implementation works on both scalars and Ts.
*/
template <typename T>
T generic_fast_erf_float(const T& a_x) {
// Clamp the inputs to the range [-4, 4] since anything outside
// this range is +/-1.0f in single-precision.
const T plus_4 = pset1<T>(4.f);
const T minus_4 = pset1<T>(-4.f);
const T x = pmax(pmin(a_x, plus_4), minus_4);
// The monomial coefficients of the numerator polynomial (odd).
const T alpha_1 = pset1<T>(-1.60960333262415e-02f);
const T alpha_3 = pset1<T>(-2.95459980854025e-03f);
const T alpha_5 = pset1<T>(-7.34990630326855e-04f);
const T alpha_7 = pset1<T>(-5.69250639462346e-05f);
const T alpha_9 = pset1<T>(-2.10102402082508e-06f);
const T alpha_11 = pset1<T>(2.77068142495902e-08f);
const T alpha_13 = pset1<T>(-2.72614225801306e-10f);
// The monomial coefficients of the denominator polynomial (even).
const T beta_0 = pset1<T>(-1.42647390514189e-02f);
const T beta_2 = pset1<T>(-7.37332916720468e-03f);
const T beta_4 = pset1<T>(-1.68282697438203e-03f);
const T beta_6 = pset1<T>(-2.13374055278905e-04f);
const T beta_8 = pset1<T>(-1.45660718464996e-05f);
// Since the polynomials are odd/even, we need x^2.
const T x2 = pmul(x, x);
// Evaluate the numerator polynomial p.
T p = pmadd(x2, alpha_13, alpha_11);
p = pmadd(x2, p, alpha_9);
p = pmadd(x2, p, alpha_7);
p = pmadd(x2, p, alpha_5);
p = pmadd(x2, p, alpha_3);
p = pmadd(x2, p, alpha_1);
p = pmul(x, p);
// Evaluate the denominator polynomial p.
T q = pmadd(x2, beta_8, beta_6);
q = pmadd(x2, q, beta_4);
q = pmadd(x2, q, beta_2);
q = pmadd(x2, q, beta_0);
// Divide the numerator by the denominator.
return pdiv(p, q);
}
template<typename RealScalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
RealScalar positive_real_hypot(const RealScalar& x, const RealScalar& y)

8
Eigen/src/Core/arch/AVX/MathFunctions.h
View File

@ -62,6 +62,14 @@ ptanh<Packet8f>(const Packet8f& x) {
return internal::generic_fast_tanh_float(x);
}
// Error function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
perf<Packet8f>(const Packet8f& x) {
return internal::generic_fast_erf_float(x);
}
// Exponential function for doubles.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4d
pexp<Packet4d>(const Packet4d& x) {

21
Eigen/src/Core/arch/AVX/PacketMath.h
View File

@ -62,21 +62,22 @@ template<> struct packet_traits<float> : default_packet_traits
enum {
Vectorizable = 1,
AlignedOnScalar = 1,
size=8,
size = 8,
HasHalfPacket = 1,
HasDiv = 1,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasLog1p = 1,
HasExpm1 = 1,
HasExp = 1,
HasDiv = 1,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasLog1p = 1,
HasExpm1 = 1,
HasExp = 1,
HasNdtri = 1,
HasBessel = 1,
HasBessel = 1,
HasSqrt = 1,
HasRsqrt = 1,
HasTanh = EIGEN_FAST_MATH,
HasTanh = EIGEN_FAST_MATH,
HasErf = EIGEN_FAST_MATH,
HasBlend = 1,
HasRound = 1,
HasFloor = 1,

12
Eigen/src/Core/arch/AVX512/MathFunctions.h
View File

@ -405,6 +405,18 @@ Packet16f pexpm1<Packet16f>(const Packet16f& _x) {
}
#endif
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
ptanh<Packet16f>(const Packet16f& _x) {
return internal::generic_fast_tanh_float(_x);
}
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet16f
perf<Packet16f>(const Packet16f& _x) {
return internal::generic_fast_erf_float(_x);
}
} // end namespace internal
} // end namespace Eigen

4
Eigen/src/Core/arch/AVX512/PacketMath.h
View File

@ -98,11 +98,13 @@ template<> struct packet_traits<float> : default_packet_traits
HasLog1p = 1,
HasExpm1 = 1,
HasNdtri = 1,
#endif
HasBessel = 1,
#endif
HasExp = 1,
HasSqrt = EIGEN_FAST_MATH,
HasRsqrt = EIGEN_FAST_MATH,
HasTanh = EIGEN_FAST_MATH,
HasErf = EIGEN_FAST_MATH,
#endif
HasDiv = 1
};

14
Eigen/src/Core/arch/AltiVec/MathFunctions.h
View File

@ -76,6 +76,20 @@ Packet2d pexp<Packet2d>(const Packet2d& _x)
}
#endif
// Hyperbolic Tangent function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
ptanh<Packet4f>(const Packet4f& x) {
return internal::generic_fast_tanh_float(x);
}
// Error function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
perf<Packet4f>(const Packet4f& x) {
return internal::generic_fast_erf_float(x);
}
} // end namespace internal
} // end namespace Eigen

35
Eigen/src/Core/arch/AltiVec/PacketMath.h
View File

@ -120,27 +120,27 @@ static Packet16uc p16uc_COMPLEX32_REV2 = vec_sld(p16uc_PSET64_HI, p16uc_PSET64_L
#define EIGEN_PPC_PREFETCH(ADDR) asm( " dcbt [%[addr]]\n" :: [addr] "r" (ADDR) : "cc" );
#endif
template<> struct packet_traits<float> : default_packet_traits
{
template <>
struct packet_traits<float> : default_packet_traits {
typedef Packet4f type;
typedef Packet4f half;
enum {
Vectorizable = 1,
AlignedOnScalar = 1,
size=4,
size = 4,
HasHalfPacket = 1,
HasAdd = 1,
HasSub = 1,
HasMul = 1,
HasDiv = 1,
HasMin = 1,
HasMax = 1,
HasAbs = 1,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasExp = 1,
HasAdd = 1,
HasSub = 1,
HasMul = 1,
HasDiv = 1,
HasMin = 1,
HasMax = 1,
HasAbs = 1,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasExp = 1,
#ifdef __VSX__
HasSqrt = 1,
#if !EIGEN_COMP_CLANG
@ -151,6 +151,8 @@ template<> struct packet_traits<float> : default_packet_traits
#else
HasSqrt = 0,
HasRsqrt = 0,
HasTanh = EIGEN_FAST_MATH,
HasErf = EIGEN_FAST_MATH,
#endif
HasRound = 1,
HasFloor = 1,
@ -159,8 +161,8 @@ template<> struct packet_traits<float> : default_packet_traits
HasBlend = 1
};
};
template<> struct packet_traits<int> : default_packet_traits
{
template <>
struct packet_traits<int> : default_packet_traits {
typedef Packet4i type;
typedef Packet4i half;
enum {
@ -177,7 +179,6 @@ template<> struct packet_traits<int> : default_packet_traits
};
};
template<> struct unpacket_traits<Packet4f>
{
typedef float type;

7
Eigen/src/Core/arch/MSA/MathFunctions.h
View File

@ -321,6 +321,13 @@ pcos<Packet4f>(const Packet4f& x) {
return psincos_inner_msa_float</* sine */ false>(x);
}
// Error function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
perf<Packet4f>(const Packet4f& x) {
return internal::generic_fast_erf_float(x);
}
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet2d
pexp<Packet2d>(const Packet2d& _x) {

1
Eigen/src/Core/arch/MSA/PacketMath.h
View File

@ -88,6 +88,7 @@ struct packet_traits<float> : default_packet_traits {
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasTanh = EIGEN_FAST_MATH,
HasErf = EIGEN_FAST_MATH,
HasLog = 1,
HasExp = 1,
HasSqrt = 1,

14
Eigen/src/Core/arch/NEON/MathFunctions.h
View File

@ -36,6 +36,20 @@ Packet4f pcos<Packet4f>(const Packet4f& x)
return pcos_float(x);
}
// Hyperbolic Tangent function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
ptanh<Packet4f>(const Packet4f& x) {
return internal::generic_fast_tanh_float(x);
}
// Error function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
perf<Packet4f>(const Packet4f& x) {
return internal::generic_fast_erf_float(x);
}
} // end namespace internal
} // end namespace Eigen

28
Eigen/src/Core/arch/NEON/PacketMath.h
View File

@ -98,28 +98,30 @@ typedef uint32x4_t Packet4ui;
#define EIGEN_ARM_PREFETCH(ADDR)
#endif
template<> struct packet_traits<float> : default_packet_traits
{
template <>
struct packet_traits<float> : default_packet_traits {
typedef Packet4f type;
typedef Packet4f half; // Packet2f intrinsics not implemented yet
typedef Packet4f half; // Packet2f intrinsics not implemented yet
enum {
Vectorizable = 1,
AlignedOnScalar = 1,
size = 4,
HasHalfPacket=0, // Packet2f intrinsics not implemented yet
HasDiv = 1,
HasHalfPacket = 0, // Packet2f intrinsics not implemented yet
HasDiv = 1,
HasFloor = 1,
// FIXME check the Has*
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasExp = 1,
HasSqrt = 0
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasExp = 1,
HasSqrt = 0,
HasTanh = EIGEN_FAST_MATH,
HasErf = EIGEN_FAST_MATH
};
};
template<> struct packet_traits<int32_t> : default_packet_traits
{
template <>
struct packet_traits<int32_t> : default_packet_traits {
typedef Packet4i type;
typedef Packet4i half; // Packet2i intrinsics not implemented yet
enum {

7
Eigen/src/Core/arch/SSE/MathFunctions.h
View File

@ -147,6 +147,13 @@ ptanh<Packet4f>(const Packet4f& x) {
return internal::generic_fast_tanh_float(x);
}
// Error function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
perf<Packet4f>(const Packet4f& a) {
return internal::generic_fast_erf_float(a);
}
} // end namespace internal
namespace numext {

29
Eigen/src/Core/arch/SSE/PacketMath.h
View File

@ -96,28 +96,29 @@ template<> struct is_arithmetic<__m128d> { enum { value = true }; };
// Use the packet_traits defined in AVX/PacketMath.h instead if we're going
// to leverage AVX instructions.
#ifndef EIGEN_VECTORIZE_AVX
template<> struct packet_traits<float> : default_packet_traits
{
template <>
struct packet_traits<float> : default_packet_traits {
typedef Packet4f type;
typedef Packet4f half;
enum {
Vectorizable = 1,
AlignedOnScalar = 1,
size=4,
size = 4,
HasHalfPacket = 0,
HasDiv = 1,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasLog1p = 1,
HasExpm1 = 1,
HasDiv = 1,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
HasLog1p = 1,
HasExpm1 = 1,
HasNdtri = 1,
HasExp = 1,
HasBessel = 1,
HasExp = 1,
HasBessel = 1,
HasSqrt = 1,
HasRsqrt = 1,
HasTanh = EIGEN_FAST_MATH,
HasTanh = EIGEN_FAST_MATH,
HasErf = EIGEN_FAST_MATH,
HasBlend = 1,
HasFloor = 1
@ -128,8 +129,8 @@ template<> struct packet_traits<float> : default_packet_traits
#endif
};
};
template<> struct packet_traits<double> : default_packet_traits
{
template <>
struct packet_traits<double> : default_packet_traits {
typedef Packet2d type;
typedef Packet2d half;
enum {

14
Eigen/src/Core/arch/ZVector/MathFunctions.h
View File

@ -225,6 +225,20 @@ Packet4f prsqrt<Packet4f>(const Packet4f& x) {
return res;
}
// Hyperbolic Tangent function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
ptanh<Packet4f>(const Packet4f& x) {
return internal::generic_fast_tanh_float(x);
}
// Error function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
perf<Packet4f>(const Packet4f& x) {
return internal::generic_fast_erf_float(x);
}
} // end namespace internal
} // end namespace Eigen

32
Eigen/src/Core/arch/ZVector/PacketMath.h
View File

@ -173,33 +173,35 @@ template<> struct packet_traits<int> : default_packet_traits
};
};
template<> struct packet_traits<float> : default_packet_traits
{
template <>
struct packet_traits<float> : default_packet_traits {
typedef Packet4f type;
typedef Packet4f half;
enum {
Vectorizable = 1,
AlignedOnScalar = 1,
size=4,
size = 4,
HasHalfPacket = 0,
HasAdd = 1,
HasSub = 1,
HasMul = 1,
HasDiv = 1,
HasMin = 1,
HasMax = 1,
HasAbs = 1,
HasSin = 0,
HasCos = 0,
HasLog = 0,
HasAdd = 1,
HasSub = 1,
HasMul = 1,
HasDiv = 1,
HasMin = 1,
HasMax = 1,
HasAbs = 1,
HasSin = 0,
HasCos = 0,
HasLog = 0,
#if !defined(__ARCH__) || (defined(__ARCH__) && __ARCH__ >= 12)
HasExp = 0,
HasExp = 0,
#else
HasExp = 1,
HasExp = 1,
#endif
HasSqrt = 1,
HasRsqrt = 1,
HasTanh = 1,
HasErf = 1,
HasRound = 1,
HasFloor = 1,
HasCeil = 1,

2
test/packetmath.cpp
View File

@ -578,7 +578,7 @@ template<typename Scalar,typename Packet> void packetmath_real()
h.store(data2, internal::plgamma(h.load(data1)));
VERIFY((numext::isnan)(data2[0]));
}
{
if (internal::packet_traits<Scalar>::HasErf) {
data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
packet_helper<internal::packet_traits<Scalar>::HasErf,Packet> h;
h.store(data2, internal::perf(h.load(data1)));

40
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h
View File

@ -238,25 +238,39 @@ struct functor_traits<scalar_polygamma_op<Scalar> >
};
/** \internal
* \brief Template functor to compute the Gauss error function of a
* scalar
* \sa class CwiseUnaryOp, Cwise::erf()
* \brief Template functor to compute the error function of a scalar
* \sa class CwiseUnaryOp, ArrayBase::erf()
*/
template<typename Scalar> struct scalar_erf_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_erf_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const {
using numext::erf; return erf(a);
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar
operator()(const Scalar& a) const {
return numext::erf(a);
}
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
return perf(x);
}
typedef typename packet_traits<Scalar>::type Packet;
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& a) const { return internal::perf(a); }
};
template<typename Scalar>
struct functor_traits<scalar_erf_op<Scalar> >
{
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
PacketAccess = packet_traits<Scalar>::HasErf,
Cost =
(PacketAccess
#ifdef EIGEN_VECTORIZE_FMA
// Haswell can issue 2 add/mul/madd per cycle.
// 10 pmadd, 2 pmul, 1 div, 2 other
? (2 * NumTraits<Scalar>::AddCost +
7 * NumTraits<Scalar>::MulCost +
scalar_div_cost<Scalar, packet_traits<Scalar>::HasDiv>::value)
#else
? (12 * NumTraits<Scalar>::AddCost +
12 * NumTraits<Scalar>::MulCost +
scalar_div_cost<Scalar, packet_traits<Scalar>::HasDiv>::value)
#endif
// Assume for simplicity that this is as expensive as an exp().
: (functor_traits<scalar_exp_op<Scalar> >::Cost))
};
};

5
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h
View File

@ -30,10 +30,6 @@ Packet pzeta(const Packet& x, const Packet& q) { using numext::zeta; return zeta
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); }
@ -77,4 +73,3 @@ Packet pbetainc(const Packet& a, const Packet& b,const Packet& x) { using numext
} // end namespace Eigen
#endif // EIGEN_SPECIALFUNCTIONS_PACKETMATH_H