diff options
| -rw-r--r-- | Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h | 17 | ||||
| -rw-r--r-- | test/array_cwise.cpp | 25 |
2 files changed, 32 insertions, 10 deletions
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index e3e91f4ab..b4fa0489b 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h @@ -929,6 +929,11 @@ Packet generic_pow(const Packet& x, const Packet& y) { Packet y_is_pos = pandnot(ptrue(y), por(y_is_zero, y_is_neg)); Packet y_is_nan = pandnot(ptrue(y), pcmp_eq(y, y)); Packet abs_y_is_inf = pcmp_eq(pabs(y), cst_pos_inf); + // |y| is so large that (1+eps)^y over- or underflows. + EIGEN_CONSTEXPR Scalar huge_exponent = + (std::numeric_limits<Scalar>::digits * Scalar(EIGEN_LOG2E)) / + std::numeric_limits<Scalar>::epsilon(); + Packet abs_y_is_huge = pcmp_lt(pset1<Packet>(huge_exponent), pabs(y)); // Predicates for whether y is integer and/or even. Packet y_is_int = pcmp_eq(pfloor(y), y); @@ -937,14 +942,16 @@ Packet generic_pow(const Packet& x, const Packet& y) { // Predicates encoding special cases for the value of pow(x,y) Packet invalid_negative_x = pandnot(pandnot(pandnot(x_is_neg, abs_x_is_inf), y_is_int), abs_y_is_inf); + Packet pow_is_one = por(por(x_is_one, y_is_zero), + pand(x_is_neg_one, + por(abs_y_is_inf, pandnot(y_is_even, invalid_negative_x)))); Packet pow_is_nan = por(invalid_negative_x, por(x_is_nan, y_is_nan)); - Packet pow_is_one = por(por(y_is_zero, x_is_one), pand(x_is_neg_one, abs_y_is_inf)); Packet pow_is_zero = por(por(por(pand(x_is_zero, y_is_pos), pand(abs_x_is_inf, y_is_neg)), - pand(pand(abs_x_is_lt_one, abs_y_is_inf), y_is_pos)), - pand(pand(abs_x_is_gt_one, abs_y_is_inf), y_is_neg)); + pand(pand(abs_x_is_lt_one, abs_y_is_huge), y_is_pos)), + pand(pand(abs_x_is_gt_one, abs_y_is_huge), y_is_neg)); Packet pow_is_inf = por(por(por(pand(x_is_zero, y_is_neg), pand(abs_x_is_inf, y_is_pos)), - pand(pand(abs_x_is_lt_one, abs_y_is_inf), y_is_neg)), - pand(pand(abs_x_is_gt_one, abs_y_is_inf), y_is_pos)); + pand(pand(abs_x_is_lt_one, abs_y_is_huge), y_is_neg)), + pand(pand(abs_x_is_gt_one, abs_y_is_huge), y_is_pos)); // General computation of pow(x,y) for positive x or negative x and integer y. Packet negate_pow_abs = pandnot(x_is_neg, y_is_even); diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp index fa087a060..a1529bc96 100644 --- a/test/array_cwise.cpp +++ b/test/array_cwise.cpp @@ -15,18 +15,31 @@ template<typename Scalar> void pow_test() { const Scalar zero = Scalar(0); const Scalar one = Scalar(1); + const Scalar two = Scalar(2); + const Scalar three = Scalar(3); const Scalar sqrt_half = Scalar(std::sqrt(0.5)); const Scalar sqrt2 = Scalar(std::sqrt(2)); const Scalar inf = std::numeric_limits<Scalar>::infinity(); const Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); - const static Scalar abs_vals[] = {zero, sqrt_half, one, sqrt2, inf, nan}; - const int abs_cases = 6; + const Scalar min = (std::numeric_limits<Scalar>::min)(); + const Scalar max = (std::numeric_limits<Scalar>::max)(); + const static Scalar abs_vals[] = {zero, + sqrt_half, + one, + sqrt2, + two, + three, + min, + max, + inf, + nan}; + + const int abs_cases = 10; 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<Scalar, Dynamic, Dynamic> x(num_repeats, num_cases); Array<Scalar, Dynamic, Dynamic> y(num_repeats, num_cases); - Array<Scalar, Dynamic, Dynamic> expected(num_repeats, num_cases); int count = 0; for (int i = 0; i < abs_cases; ++i) { const Scalar abs_x = abs_vals[i]; @@ -39,7 +52,6 @@ void pow_test() { for (int repeat = 0; repeat < num_repeats; ++repeat) { x(repeat, count) = x_case; y(repeat, count) = y_case; - expected(repeat, count) = numext::pow(x_case, y_case); } ++count; } @@ -52,8 +64,11 @@ void pow_test() { bool all_pass = true; for (int i = 0; i < 1; ++i) { for (int j = 0; j < num_cases; ++j) { + // TODO(rmlarsen): Skip tests that trigger a known bug in pldexp for now. + if (std::abs(x(i,j)) == max || std::abs(x(i,j)) == min) continue; + + Scalar e = numext::pow(x(i,j), y(i,j)); Scalar a = actual(i, j); - Scalar e = expected(i, j); bool fail = !(a==e) && !internal::isApprox(a, e, tol) && !((numext::isnan)(a) && (numext::isnan)(e)); all_pass &= !fail; if (fail) { |