mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add more tests for pow and fix a corner case for huge exponent where the result is always zero or infinite unless x is one.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2021-02-05 16:58:49 -08:00
parent abcde69a79
commit 6e3b795f81
2 changed files with 32 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -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);

25
test/array_cwise.cpp
View File

@ -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) {