From 14c847dc0eda98c4ad848766225222a26676e60e Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Wed, 12 Oct 2022 20:12:08 +0000 Subject: [PATCH] Refactor special values test for pow, and add a similar test for atan2 --- Eigen/src/Core/functors/BinaryFunctors.h | 2 +- test/array_cwise.cpp | 134 +++++++++++++---------- 2 files changed, 76 insertions(+), 60 deletions(-) diff --git a/Eigen/src/Core/functors/BinaryFunctors.h b/Eigen/src/Core/functors/BinaryFunctors.h index 5bb4e0f9e..38f123cff 100644 --- a/Eigen/src/Core/functors/BinaryFunctors.h +++ b/Eigen/src/Core/functors/BinaryFunctors.h @@ -523,7 +523,7 @@ struct scalar_atan2_op { // See https://en.cppreference.com/w/cpp/numeric/math/atan2 // for how corner cases are supposed to be handled according to the // IEEE floating-point standard (IEC 60559). - const Packet kSignMask = pnegate(pzero(x)); + const Packet kSignMask = pset1(-Scalar(0.0)); const Packet kPi = pset1(Scalar(EIGEN_PI)); const Packet kPiO2 = pset1(Scalar(EIGEN_PI / 2)); const Packet kPiO4 = pset1(Scalar(EIGEN_PI / 4)); diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp index 73393ff76..a7e0ff48c 100644 --- a/test/array_cwise.cpp +++ b/test/array_cwise.cpp @@ -7,12 +7,11 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include #include "main.h" - -// Test the corner cases of pow(x, y) for real types. -template -void pow_test() { +template +std::vector special_values() { const Scalar zero = Scalar(0); const Scalar eps = Eigen::NumTraits::epsilon(); const Scalar one = Scalar(1); @@ -27,25 +26,19 @@ void pow_test() { const Scalar max = (std::numeric_limits::max)(); const Scalar max_exp = (static_cast(int(Eigen::NumTraits::max_exponent())) * Scalar(EIGEN_LN2)) / eps; - const static Scalar abs_vals[] = {zero, - denorm_min, - min, - eps, - sqrt_half, - one, - sqrt2, - two, - three, - max_exp, - max, - inf, - nan}; - const int abs_cases = 13; + return {zero, denorm_min, min, eps, sqrt_half, one, sqrt2, two, three, max_exp, max, inf, nan}; +} + +template +void special_value_pairs(Array& x, + Array& y) { + std::vector abs_vals = special_values(); + const int abs_cases = abs_vals.size(); const int num_cases = 2*abs_cases * 2*abs_cases; - // Repeat the same value to make sure we hit the vectorized path. - const int num_repeats = 32; - Array x(num_repeats, num_cases); - Array y(num_repeats, num_cases); + // ensure both vectorized and non-vectorized paths taken + const int num_repeats = 2 * internal::packet_traits::size + 1; + x.resize(num_repeats, num_cases); + y.resize(num_repeats, num_cases); int count = 0; for (int i = 0; i < abs_cases; ++i) { const Scalar abs_x = abs_vals[i]; @@ -64,65 +57,85 @@ void pow_test() { } } } +} - Array actual = x.pow(y); +template +void binary_op_test(std::string name, Fn fun, RefFn ref) { const Scalar tol = test_precision(); + Array x; + Array y; + special_value_pairs(x, y); + + Array actual = fun(x, y); bool all_pass = true; - for (int i = 0; i < 1; ++i) { - for (int j = 0; j < num_cases; ++j) { - Scalar e = static_cast(std::pow(x(i,j), y(i,j))); + for (int i = 0; i < x.rows(); ++i) { + for (int j = 0; j < x.cols(); ++j) { + Scalar e = static_cast(ref(x(i,j), y(i,j))); Scalar a = actual(i, j); bool success = (a==e) || ((numext::isfinite)(e) && internal::isApprox(a, e, tol)) || ((numext::isnan)(a) && (numext::isnan)(e)); all_pass &= success; if (!success) { - std::cout << "pow(" << x(i,j) << "," << y(i,j) << ") = " << a << " != " << e << std::endl; + std::cout << name << "(" << x(i,j) << "," << y(i,j) << ") = " << a << " != " << e << std::endl; } } } + VERIFY(all_pass); +} - typedef typename internal::make_integer::type Int_t; +template +void binary_ops_test() { + binary_op_test("pow", + [](auto x, auto y) { return Eigen::pow(x, y); }, + [](auto x, auto y) { return std::pow(x, y); }); + binary_op_test("atan2", + [](auto x, auto y) { return Eigen::atan2(x, y); }, + [](auto x, auto y) { return std::atan2(x, y); }); +} - // ensure both vectorized and non-vectorized paths taken - Index test_size = 2 * internal::packet_traits::size + 1; - - Array eigenPow(test_size); - for (int i = 0; i < num_cases; ++i) { - Array bases = x.col(i); - for (Scalar abs_exponent : abs_vals){ - for (Scalar exponent : {-abs_exponent, abs_exponent}){ - // test floating point exponent code path - eigenPow.setZero(); - eigenPow = bases.pow(exponent); - for (int j = 0; j < num_repeats; j++){ +template +void pow_scalar_exponent_test() { + using Int_t = typename internal::make_integer::type; + const Scalar tol = test_precision(); + + std::vector abs_vals = special_values(); + const int num_vals = abs_vals.size(); + Map> bases(abs_vals.data(), num_vals); + + bool all_pass = true; + for (Scalar abs_exponent : abs_vals) { + for (Scalar exponent : {-abs_exponent, abs_exponent}) { + // test integer exponent code path + bool exponent_is_integer = (numext::isfinite)(exponent) && (numext::round(exponent) == exponent) && + (numext::abs(exponent) < static_cast(NumTraits::highest())); + if (exponent_is_integer) { + Int_t exponent_as_int = static_cast(exponent); + Array eigenPow = bases.pow(exponent_as_int); + for (int j = 0; j < num_vals; j++) { Scalar e = static_cast(std::pow(bases(j), exponent)); Scalar a = eigenPow(j); - bool success = (a == e) || ((numext::isfinite)(e) && internal::isApprox(a, e, tol)) || ((numext::isnan)(a) && (numext::isnan)(e)); + bool success = (a == e) || ((numext::isfinite)(e) && internal::isApprox(a, e, tol)) || + ((numext::isnan)(a) && (numext::isnan)(e)); all_pass &= success; if (!success) { - std::cout << "pow(" << x(i, j) << "," << y(i, j) << ") = " << a << " != " << e << std::endl; + std::cout << "pow(" << bases(j) << "," << exponent << ") = " << a << " != " << e << std::endl; } } - // test integer exponent code path - bool exponent_is_integer = (numext::isfinite)(exponent) && (numext::round(exponent) == exponent) && (numext::abs(exponent) < static_cast(NumTraits::highest())); - if (exponent_is_integer) - { - Int_t exponent_as_int = static_cast(exponent); - eigenPow.setZero(); - eigenPow = bases.pow(exponent_as_int); - for (int j = 0; j < num_repeats; j++){ - Scalar e = static_cast(std::pow(bases(j), exponent)); - Scalar a = eigenPow(j); - bool success = (a == e) || ((numext::isfinite)(e) && internal::isApprox(a, e, tol)) || ((numext::isnan)(a) && (numext::isnan)(e)); - all_pass &= success; - if (!success) { - std::cout << "pow(" << x(i, j) << "," << y(i, j) << ") = " << a << " != " << e << std::endl; - } + } else { + // test floating point exponent code path + Array eigenPow = bases.pow(exponent); + for (int j = 0; j < num_vals; j++) { + Scalar e = static_cast(std::pow(bases(j), exponent)); + Scalar a = eigenPow(j); + bool success = (a == e) || ((numext::isfinite)(e) && internal::isApprox(a, e, tol)) || + ((numext::isnan)(a) && (numext::isnan)(e)); + all_pass &= success; + if (!success) { + std::cout << "pow(" << bases(j) << "," << exponent << ") = " << a << " != " << e << std::endl; } } } } } - VERIFY(all_pass); } @@ -626,7 +639,10 @@ template void array_real(const ArrayType& m) // Avoid inf and NaN. m3 = (m1.square()::epsilon()).select(Scalar(1),m3); VERIFY_IS_APPROX(m3.pow(RealScalar(-2)), m3.square().inverse()); - pow_test(); + + // Test pow and atan2 on special IEEE values. + binary_ops_test(); + pow_scalar_exponent_test(); VERIFY_IS_APPROX(log10(m3), log(m3)/numext::log(Scalar(10))); VERIFY_IS_APPROX(log2(m3), log(m3)/numext::log(Scalar(2)));