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

Make sure exp(-Inf) is zero for vectorized expressions. This fixes #2385.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2021-12-08 17:57:23 +00:00
parent 39a6aff16c
commit f04fd8b168
2 changed files with 10 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
View File

@ -440,6 +440,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
EIGEN_UNUSED
Packet pexp_float(const Packet _x)
{
const Packet cst_zero = pset1<Packet>(0.0f);
const Packet cst_1 = pset1<Packet>(1.0f);
const Packet cst_half = pset1<Packet>(0.5f);
const Packet cst_exp_hi = pset1<Packet>( 88.723f);
@ -454,7 +455,8 @@ Packet pexp_float(const Packet _x)
const Packet cst_cephes_exp_p5 = pset1<Packet>(5.0000001201E-1f);
// Clamp x.
Packet x = pmax(pmin(_x, cst_exp_hi), cst_exp_lo);
Packet zero_mask = pcmp_lt(_x, cst_exp_lo);
Packet x = pmin(_x, cst_exp_hi);
// Express exp(x) as exp(m*ln(2) + r), start by extracting
// m = floor(x/ln(2) + 0.5).
@ -483,7 +485,7 @@ Packet pexp_float(const Packet _x)
// Return 2^m * exp(r).
// TODO: replace pldexp with faster implementation since y in [-1, 1).
return pmax(pldexp(y,m), _x);
return pselect(zero_mask, cst_zero, pmax(pldexp(y,m), _x));
}
template <typename Packet>
@ -492,7 +494,7 @@ EIGEN_UNUSED
Packet pexp_double(const Packet _x)
{
Packet x = _x;
const Packet cst_zero = pset1<Packet>(0.0f);
const Packet cst_1 = pset1<Packet>(1.0);
const Packet cst_2 = pset1<Packet>(2.0);
const Packet cst_half = pset1<Packet>(0.5);
@ -514,7 +516,8 @@ Packet pexp_double(const Packet _x)
Packet tmp, fx;
// clamp x
x = pmax(pmin(x, cst_exp_hi), cst_exp_lo);
Packet zero_mask = pcmp_lt(_x, cst_exp_lo);
x = pmin(x, cst_exp_hi);
// Express exp(x) as exp(g + n*log(2)).
fx = pmadd(cst_cephes_LOG2EF, x, cst_half);
@ -552,7 +555,7 @@ Packet pexp_double(const Packet _x)
// Construct the result 2^n * exp(g) = e * x. The max is used to catch
// non-finite values in the input.
// TODO: replace pldexp with faster implementation since x in [-1, 1).
return pmax(pldexp(x,fx), _x);
return pselect(zero_mask, cst_zero, pmax(pldexp(x,fx), _x));
}
// The following code is inspired by the following stack-overflow answer:

3
test/packetmath.cpp
View File

@ -48,7 +48,7 @@ inline bool REF_MUL(const bool& a, const bool& b) {
template <typename T>
inline T REF_FREXP(const T& x, T& exp) {
int iexp;
int iexp = 0;
EIGEN_USING_STD(frexp)
const T out = static_cast<T>(frexp(x, &iexp));
exp = static_cast<T>(iexp);
@ -713,6 +713,7 @@ void packetmath_real() {
for (int i = 0; i < size; ++i) {
data1[i] = Scalar(internal::random<double>(-87, 88));
data2[i] = Scalar(internal::random<double>(-87, 88));
data1[0] = -NumTraits<Scalar>::infinity();
}
CHECK_CWISE1_IF(PacketTraits::HasExp, std::exp, internal::pexp);