From 0f6f75bd8a0445edc3361659e065f15a29e2743c Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 23 Dec 2018 17:26:21 +0100 Subject: [PATCH] Implement a faster fix for sin/cos of large entries that also correctly handle INF input. --- .../arch/Default/GenericPacketMathFunctions.h | 13 +++++++------ test/packetmath.cpp | 16 ++++++++++++++++ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index 80bcc077d..7ceaea894 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h @@ -289,15 +289,9 @@ Packet psincos_float(const Packet& _x) const Packet cst_coscof_p1 = pset1(-1.388731625493765E-003f); const Packet cst_coscof_p2 = pset1( 4.166664568298827E-002f); const Packet cst_cephes_FOPI = pset1( 1.27323954473516f); // 4 / M_PI - const Packet cst_sincos_max_arg = pset1( 13176795.0f); // Approx. (2**24) / (4/Pi). Packet x = pabs(_x); - // Prevent sin/cos from generating values larger than 1.0 in magnitude - // for very large arguments by setting x to 0.0. - Packet small_or_nan_mask = pcmp_lt_or_nan(x, cst_sincos_max_arg); - x = pand(x, small_or_nan_mask); - // Scale x by 4/Pi to find x's octant. Packet y = pmul(x, cst_cephes_FOPI); @@ -348,6 +342,13 @@ Packet psincos_float(const Packet& _x) y = ComputeSine ? pselect(poly_mask,y2,y1) : pselect(poly_mask,y1,y2); + // For very large arguments the the reduction to the [-Pi/4,+Pi/4] range + // does not work thus leading to sine/cosine out of the [-1:1] range. + // Since computing the sine/cosine for very large entry entries makes little + // sense in term of accuracy, we simply clamp to [-1,1]: + y = pmin(y,pset1( 1.f)); + y = pmax(y,pset1(-1.f)); + // Update the sign return pxor(y, sign_bit); } diff --git a/test/packetmath.cpp b/test/packetmath.cpp index bf0312a73..9f647530b 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -579,6 +579,22 @@ template void packetmath_real() VERIFY(data2[0]<=Scalar(1.) && data2[0]>=Scalar(-1.)); VERIFY(data2[1]<=Scalar(1.) && data2[1]>=Scalar(-1.)); } + + data1[0] = std::numeric_limits::infinity(); + data1[1] = -std::numeric_limits::infinity(); + h.store(data2, internal::psin(h.load(data1))); + VERIFY((numext::isnan)(data2[0])); + VERIFY((numext::isnan)(data2[1])); + + h.store(data2, internal::pcos(h.load(data1))); + VERIFY((numext::isnan)(data2[0])); + VERIFY((numext::isnan)(data2[1])); + + data1[0] = std::numeric_limits::quiet_NaN(); + h.store(data2, internal::psin(h.load(data1))); + VERIFY((numext::isnan)(data2[0])); + h.store(data2, internal::pcos(h.load(data1))); + VERIFY((numext::isnan)(data2[0])); } } }