GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 03:39:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add more tests for corner cases of log1p and expm1. Add handling of infinite arguments to log1p such that log1p(inf) = inf.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2019-08-28 12:20:21 -07:00
parent 6e77f9bef3
commit 1187bb65ad
3 changed files with 17 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -551,7 +551,8 @@ namespace std_fallback {
Scalar x1p = RealScalar(1) + x; Scalar x1p = RealScalar(1) + x;
Scalar log_1p = log(x1p); Scalar log_1p = log(x1p);
const bool is_small = numext::equal_strict(x1p, Scalar(1)); const bool is_small = numext::equal_strict(x1p, Scalar(1));
return is_small ? x : x * (log_1p / (x1p - RealScalar(1))); const bool is_inf = numext::equal_strict(x1p, log_1p);
return (is_small || is_inf) ? x : x * (log_1p / (x1p - RealScalar(1)));
} }
} }

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

@ -137,8 +137,9 @@ Packet generic_plog1p(const Packet& x)
Packet xp1 = padd(x, one); Packet xp1 = padd(x, one);
Packet small_mask = pcmp_eq(xp1, one); Packet small_mask = pcmp_eq(xp1, one);
Packet log1 = plog(xp1); Packet log1 = plog(xp1);
Packet inf_mask = pcmp_eq(xp1, log1);
Packet log_large = pmul(x, pdiv(log1, psub(xp1, one))); Packet log_large = pmul(x, pdiv(log1, psub(xp1, one)));
return pselect(small_mask, x, log_large); return pselect(por(small_mask, inf_mask), x, log_large);
} }
/** \internal \returns exp(x)-1 computed using W. Kahan's formula. /** \internal \returns exp(x)-1 computed using W. Kahan's formula.

14
test/packetmath.cpp
View File

@ -607,8 +607,12 @@ template<typename Scalar,typename Packet> void packetmath_real()
CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma); CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma);
CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf); CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf);
CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc); CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc);
CHECK_CWISE1_IF(PacketTraits::HasExpm1, std::expm1, internal::pexpm1); data1[0] = std::numeric_limits<Scalar>::infinity();
data1[1] = Scalar(-1);
CHECK_CWISE1_IF(PacketTraits::HasLog1p, std::log1p, internal::plog1p); CHECK_CWISE1_IF(PacketTraits::HasLog1p, std::log1p, internal::plog1p);
data1[0] = std::numeric_limits<Scalar>::infinity();
data1[1] = -std::numeric_limits<Scalar>::infinity();
CHECK_CWISE1_IF(PacketTraits::HasExpm1, std::expm1, internal::pexpm1);
#endif #endif
if(PacketSize>=2) if(PacketSize>=2)
@ -648,6 +652,14 @@ template<typename Scalar,typename Packet> void packetmath_real()
h.store(data2, internal::plog(h.load(data1))); h.store(data2, internal::plog(h.load(data1)));
VERIFY((numext::isinf)(data2[0])); VERIFY((numext::isinf)(data2[0]));
} }
if(PacketTraits::HasLog1p) {
packet_helper<PacketTraits::HasLog1p,Packet> h;
data1[0] = Scalar(-2);
data1[1] = -std::numeric_limits<Scalar>::infinity();
h.store(data2, internal::plog1p(h.load(data1)));
VERIFY((numext::isnan)(data2[0]));
VERIFY((numext::isnan)(data2[1]));
}
if(PacketTraits::HasSqrt) if(PacketTraits::HasSqrt)
{ {
packet_helper<PacketTraits::HasSqrt,Packet> h; packet_helper<PacketTraits::HasSqrt,Packet> h;