Internal change
PiperOrigin-RevId: 538300588
Change-Id: I9da178ce4621e072d87b4d1219aa97d4c7d86ce0
diff --git a/Eigen/src/Core/MathFunctionsImpl.h b/Eigen/src/Core/MathFunctionsImpl.h
index e5ae03d..86c06b2 100644
--- a/Eigen/src/Core/MathFunctionsImpl.h
+++ b/Eigen/src/Core/MathFunctionsImpl.h
@@ -80,26 +80,26 @@
using Scalar = typename unpacket_traits<Packet>::type;
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Packet
run(const Packet& a, const Packet& approx_rsqrt) {
+ constexpr Scalar kMinusHalf = Scalar(-1)/Scalar(2);
+ const Packet cst_minus_half = pset1<Packet>(kMinusHalf);
const Packet cst_minus_one = pset1<Packet>(Scalar(-1));
- const Packet cst_minus_half = pset1<Packet>(Scalar(-1)/Scalar(2));
-
- // Refine the approximation using one Newton-Raphson step:
- // The approximation is expressed this way to avoid over/under-flows.
- // x' = x - (x/2) * ( (a*x)*x - 1)
- Packet x = approx_rsqrt;
+ Packet inv_sqrt = approx_rsqrt;
for (int step = 0; step < Steps; ++step) {
- Packet minushalfx = pmul(cst_minus_half, x);
- Packet ax = pmul(a, x);
- Packet ax2m1 = pmadd(ax, x, cst_minus_one);
- x = pmadd(ax2m1, minushalfx, x);
+ // Refine the approximation using one Newton-Raphson step:
+ // h_n = (x * inv_sqrt) * inv_sqrt - 1 (so that h_n is nearly 0).
+ // inv_sqrt = inv_sqrt - 0.5 * inv_sqrt * h_n
+ Packet r2 = pmul(a, inv_sqrt);
+ Packet half_r = pmul(inv_sqrt, cst_minus_half);
+ Packet h_n = pmadd(r2, inv_sqrt, cst_minus_one);
+ inv_sqrt = pmadd(half_r, h_n, inv_sqrt);
}
// If x is NaN, then either:
// 1) the input is NaN
// 2) zero and infinity were multiplied
// In either of these cases, return approx_rsqrt
- return pselect(pisnan(x), approx_rsqrt, x);
+ return pselect(pisnan(inv_sqrt), approx_rsqrt, inv_sqrt);
}
};
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
index 13fafe3..1bfee7d 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
@@ -1970,258 +1970,204 @@
};
namespace unary_pow {
-template <typename ScalarExponent, bool IsIntegerAtCompileTime = NumTraits<ScalarExponent>::IsInteger>
-struct is_odd {
- static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ScalarExponent run(const ScalarExponent& x) {
- ScalarExponent xdiv2 = x / ScalarExponent(2);
- ScalarExponent floorxdiv2 = numext::floor(xdiv2);
- return xdiv2 != floorxdiv2;
+
+template <typename ScalarExponent, bool IsInteger = NumTraits<ScalarExponent>::IsInteger>
+struct exponent_helper {
+ using safe_abs_type = ScalarExponent;
+ static constexpr ScalarExponent one_half = ScalarExponent(0.5);
+ // these routines assume that exp is an integer stored as a floating point type
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ScalarExponent safe_abs(const ScalarExponent& exp) {
+ return numext::abs(exp);
+ }
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_odd(const ScalarExponent& exp) {
+ eigen_assert(((numext::isfinite)(exp) && exp == numext::floor(exp)) && "exp must be an integer");
+ ScalarExponent exp_div_2 = exp * one_half;
+ ScalarExponent floor_exp_div_2 = numext::floor(exp_div_2);
+ return exp_div_2 != floor_exp_div_2;
+ }
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ScalarExponent floor_div_two(const ScalarExponent& exp) {
+ ScalarExponent exp_div_2 = exp * one_half;
+ return numext::floor(exp_div_2);
}
};
+
template <typename ScalarExponent>
-struct is_odd<ScalarExponent, true> {
- static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ScalarExponent run(const ScalarExponent& x) {
- return x % ScalarExponent(2);
+struct exponent_helper<ScalarExponent, true> {
+ // if `exp` is a signed integer type, cast it to its unsigned counterpart to safely store its absolute value
+ // consider the (rare) case where `exp` is an int32_t: abs(-2147483648) != 2147483648
+ using safe_abs_type = typename numext::get_integer_by_size<sizeof(ScalarExponent)>::unsigned_type;
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE safe_abs_type safe_abs(const ScalarExponent& exp) {
+ ScalarExponent mask = exp ^ numext::abs(exp);
+ safe_abs_type result = static_cast<safe_abs_type>(exp);
+ return result ^ mask;
+ }
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_odd(const safe_abs_type& exp) {
+ return exp % safe_abs_type(2) != safe_abs_type(0);
+ }
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE safe_abs_type floor_div_two(const safe_abs_type& exp) {
+ return exp >> safe_abs_type(1);
}
};
template <typename Packet, typename ScalarExponent,
- bool BaseIsIntegerType = NumTraits<typename unpacket_traits<Packet>::type>::IsInteger>
-struct do_div {
+ bool ReciprocateIfExponentIsNegative =
+ !NumTraits<typename unpacket_traits<Packet>::type>::IsInteger && NumTraits<ScalarExponent>::IsSigned>
+struct reciprocate {
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent& exponent) {
- typedef typename unpacket_traits<Packet>::type Scalar;
+ using Scalar = typename unpacket_traits<Packet>::type;
const Packet cst_pos_one = pset1<Packet>(Scalar(1));
return exponent < 0 ? pdiv(cst_pos_one, x) : x;
}
};
template <typename Packet, typename ScalarExponent>
-struct do_div<Packet, ScalarExponent, true> {
+struct reciprocate<Packet, ScalarExponent, false> {
// pdiv not defined, nor necessary for integer base types
- static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent& exponent) {
- EIGEN_UNUSED_VARIABLE(exponent);
- return x;
- }
+ // if the exponent is unsigned, then the exponent cannot be negative
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent&) { return x; }
};
template <typename Packet, typename ScalarExponent>
-static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet int_pow(const Packet& x, const ScalarExponent& exponent) {
- typedef typename unpacket_traits<Packet>::type Scalar;
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet int_pow(const Packet& x, const ScalarExponent& exponent) {
+ using Scalar = typename unpacket_traits<Packet>::type;
+ using ExponentHelper = exponent_helper<ScalarExponent>;
+ using AbsExponentType = typename ExponentHelper::safe_abs_type;
const Packet cst_pos_one = pset1<Packet>(Scalar(1));
- if (exponent == 0) return cst_pos_one;
- Packet result = x;
+ if (exponent == ScalarExponent(0)) return cst_pos_one;
+
+ Packet result = reciprocate<Packet, ScalarExponent>::run(x, exponent);
Packet y = cst_pos_one;
- ScalarExponent m = numext::abs(exponent);
+ AbsExponentType m = ExponentHelper::safe_abs(exponent);
+
while (m > 1) {
- bool odd = is_odd<ScalarExponent>::run(m);
+ bool odd = ExponentHelper::is_odd(m);
if (odd) y = pmul(y, result);
result = pmul(result, result);
- m = numext::floor(m / ScalarExponent(2));
+ m = ExponentHelper::floor_div_two(m);
}
- result = pmul(y, result);
- result = do_div<Packet, ScalarExponent>::run(result, exponent);
- return result;
+
+ return pmul(y, result);
}
template <typename Packet>
-static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet gen_pow(const Packet& x,
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet gen_pow(const Packet& x,
const typename unpacket_traits<Packet>::type& exponent) {
const Packet exponent_packet = pset1<Packet>(exponent);
return generic_pow_impl(x, exponent_packet);
}
template <typename Packet, typename ScalarExponent>
-static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet handle_nonint_int_errors(const Packet& x, const Packet& powx,
- const ScalarExponent& exponent) {
- typedef typename unpacket_traits<Packet>::type Scalar;
-
- // non-integer base, integer exponent case
-
- const bool exponent_is_odd = is_odd<ScalarExponent>::run(exponent);
- const bool exponent_is_neg = exponent < 0;
-
- const Packet exp_is_odd = exponent_is_odd ? ptrue(x) : pzero(x);
- const Packet exp_is_neg = exponent_is_neg ? ptrue(x) : pzero(x);
-
- const Scalar pos_zero = Scalar(0);
- const Scalar neg_zero = -Scalar(0);
- const Scalar pos_one = Scalar(1);
- const Scalar pos_inf = NumTraits<Scalar>::infinity();
- const Scalar neg_inf = -NumTraits<Scalar>::infinity();
-
- const Packet cst_pos_zero = pset1<Packet>(pos_zero);
- const Packet cst_neg_zero = pset1<Packet>(neg_zero);
- const Packet cst_pos_one = pset1<Packet>(pos_one);
- const Packet cst_pos_inf = pset1<Packet>(pos_inf);
- const Packet cst_neg_inf = pset1<Packet>(neg_inf);
-
- const Packet abs_x = pabs(x);
- const Packet abs_x_is_zero = pcmp_eq(abs_x, cst_pos_zero);
- const Packet abs_x_is_one = pcmp_eq(abs_x, cst_pos_one);
- const Packet abs_x_is_inf = pcmp_eq(abs_x, cst_pos_inf);
-
- const Packet x_has_signbit = psignbit(x);
- const Packet x_is_neg = pandnot(x_has_signbit, abs_x_is_zero);
- const Packet x_is_neg_zero = pand(x_has_signbit, abs_x_is_zero);
-
- if (exponent == 0) {
- return cst_pos_one;
- }
-
- Packet pow_is_pos_inf = pand(pandnot(abs_x_is_zero, x_is_neg_zero), pand(exp_is_odd, exp_is_neg));
- pow_is_pos_inf = por(pow_is_pos_inf, pand(abs_x_is_zero, pandnot(exp_is_neg, exp_is_odd)));
- pow_is_pos_inf = por(pow_is_pos_inf, pand(pand(abs_x_is_inf, x_is_neg), pandnot(pnot(exp_is_neg), exp_is_odd)));
- pow_is_pos_inf = por(pow_is_pos_inf, pandnot(pandnot(abs_x_is_inf, x_is_neg), exp_is_neg));
-
- Packet pow_is_neg_inf = pand(x_is_neg_zero, pand(exp_is_neg, exp_is_odd));
- pow_is_neg_inf = por(pow_is_neg_inf, pand(pand(abs_x_is_inf, x_is_neg), pandnot(exp_is_odd, exp_is_neg)));
-
- Packet pow_is_pos_zero = pandnot(abs_x_is_zero, exp_is_neg);
- pow_is_pos_zero = por(pow_is_pos_zero, pand(pand(abs_x_is_inf, x_is_neg), pandnot(exp_is_neg, exp_is_odd)));
- pow_is_pos_zero = por(pow_is_pos_zero, pand(pandnot(abs_x_is_inf, x_is_neg), exp_is_neg));
-
- Packet pow_is_neg_zero = pand(x_is_neg_zero, pandnot(exp_is_odd, exp_is_neg));
- pow_is_neg_zero = por(pow_is_neg_zero, pand(pand(abs_x_is_inf, x_is_neg), pand(exp_is_odd, exp_is_neg)));
-
- Packet result = pselect(pow_is_neg_inf, cst_neg_inf, powx);
- result = pselect(pow_is_neg_zero, cst_neg_zero, result);
- result = pselect(pow_is_pos_zero, cst_pos_zero, result);
- result = pselect(pow_is_pos_inf, cst_pos_inf, result);
- result = pselect(pandnot(abs_x_is_one, x_is_neg), cst_pos_one, result);
- return result;
-}
-
-template <typename Packet, typename ScalarExponent>
-static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet handle_nonint_nonint_errors(const Packet& x, const Packet& powx,
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet handle_nonint_nonint_errors(const Packet& x, const Packet& powx,
const ScalarExponent& exponent) {
- typedef typename unpacket_traits<Packet>::type Scalar;
+ using Scalar = typename unpacket_traits<Packet>::type;
// non-integer base and exponent case
- const bool exponent_is_fin = (numext::isfinite)(exponent);
- const bool exponent_is_nan = (numext::isnan)(exponent);
- const bool exponent_is_neg = exponent < 0;
- const bool exponent_is_inf = !exponent_is_fin && !exponent_is_nan;
-
- const Packet exp_is_neg = exponent_is_neg ? ptrue(x) : pzero(x);
- const Packet exp_is_inf = exponent_is_inf ? ptrue(x) : pzero(x);
-
const Scalar pos_zero = Scalar(0);
+ const Scalar all_ones = ptrue<Scalar>(Scalar());
const Scalar pos_one = Scalar(1);
const Scalar pos_inf = NumTraits<Scalar>::infinity();
- const Scalar nan = NumTraits<Scalar>::quiet_NaN();
- const Packet cst_pos_zero = pset1<Packet>(pos_zero);
+ const Packet cst_pos_zero = pzero(x);
const Packet cst_pos_one = pset1<Packet>(pos_one);
const Packet cst_pos_inf = pset1<Packet>(pos_inf);
- const Packet cst_nan = pset1<Packet>(nan);
+
+ const bool exponent_is_not_fin = !(numext::isfinite)(exponent);
+ const bool exponent_is_neg = exponent < ScalarExponent(0);
+ const bool exponent_is_pos = exponent > ScalarExponent(0);
+
+ const Packet exp_is_not_fin = pset1<Packet>(exponent_is_not_fin ? all_ones : pos_zero);
+ const Packet exp_is_neg = pset1<Packet>(exponent_is_neg ? all_ones : pos_zero);
+ const Packet exp_is_pos = pset1<Packet>(exponent_is_pos ? all_ones : pos_zero);
+ const Packet exp_is_inf = pand(exp_is_not_fin, por(exp_is_neg, exp_is_pos));
+ const Packet exp_is_nan = pandnot(exp_is_not_fin, por(exp_is_neg, exp_is_pos));
+
+ const Packet x_is_le_zero = pcmp_le(x, cst_pos_zero);
+ const Packet x_is_ge_zero = pcmp_le(cst_pos_zero, x);
+ const Packet x_is_zero = pand(x_is_le_zero, x_is_ge_zero);
const Packet abs_x = pabs(x);
- const Packet abs_x_is_zero = pcmp_eq(abs_x, cst_pos_zero);
- const Packet abs_x_is_lt_one = pcmp_lt(abs_x, cst_pos_one);
- const Packet abs_x_is_gt_one = pcmp_lt(cst_pos_one, abs_x);
- const Packet abs_x_is_one = pcmp_eq(abs_x, cst_pos_one);
+ const Packet abs_x_is_le_one = pcmp_le(abs_x, cst_pos_one);
+ const Packet abs_x_is_ge_one = pcmp_le(cst_pos_one, abs_x);
const Packet abs_x_is_inf = pcmp_eq(abs_x, cst_pos_inf);
+ const Packet abs_x_is_one = pand(abs_x_is_le_one, abs_x_is_ge_one);
- const Packet x_has_signbit = psignbit(x);
- const Packet x_is_neg = pandnot(x_has_signbit, abs_x_is_zero);
+ Packet pow_is_inf_if_exp_is_neg = por(x_is_zero, pand(abs_x_is_le_one, exp_is_inf));
+ Packet pow_is_inf_if_exp_is_pos = por(abs_x_is_inf, pand(abs_x_is_ge_one, exp_is_inf));
+ Packet pow_is_one = pand(abs_x_is_one, por(exp_is_inf, x_is_ge_zero));
- if (exponent_is_nan) {
- return pselect(pandnot(abs_x_is_one, x_is_neg), cst_pos_one, cst_nan);
- }
-
- Packet pow_is_pos_zero = pandnot(abs_x_is_zero, exp_is_neg);
- pow_is_pos_zero = por(pow_is_pos_zero, pand(abs_x_is_gt_one, pand(exp_is_inf, exp_is_neg)));
- pow_is_pos_zero = por(pow_is_pos_zero, pand(abs_x_is_lt_one, pandnot(exp_is_inf, exp_is_neg)));
- pow_is_pos_zero = por(pow_is_pos_zero, pand(abs_x_is_inf, exp_is_neg));
-
- const Packet pow_is_pos_one = pand(abs_x_is_one, exp_is_inf);
-
- Packet pow_is_pos_inf = pand(abs_x_is_zero, exp_is_neg);
- pow_is_pos_inf = por(pow_is_pos_inf, pand(abs_x_is_lt_one, pand(exp_is_inf, exp_is_neg)));
- pow_is_pos_inf = por(pow_is_pos_inf, pand(abs_x_is_gt_one, pandnot(exp_is_inf, exp_is_neg)));
- pow_is_pos_inf = por(pow_is_pos_inf, pandnot(abs_x_is_inf, exp_is_neg));
-
- const Packet pow_is_nan = pandnot(pandnot(x_is_neg, abs_x_is_inf), exp_is_inf);
-
- Packet result = pselect(pow_is_pos_inf, cst_pos_inf, powx);
- result = pselect(pow_is_pos_one, cst_pos_one, result);
- result = pselect(pow_is_pos_zero, cst_pos_zero, result);
- result = pselect(pow_is_nan, cst_nan, result);
- result = pselect(pandnot(abs_x_is_one, x_is_neg), cst_pos_one, result);
- return result;
-}
-
-template <typename Packet, typename ScalarExponent, std::enable_if_t<NumTraits<typename unpacket_traits<Packet>::type>::IsSigned, bool> = true>
-static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet handle_int_int(const Packet& x, const ScalarExponent& exponent) {
- typedef typename unpacket_traits<Packet>::type Scalar;
-
- // signed integer base, integer exponent case
-
- // This routine handles negative and very large positive exponents
- // Signed integer overflow and divide by zero is undefined behavior
- // Unsigned integers do not overflow
-
- const bool exponent_is_odd = unary_pow::is_odd<ScalarExponent>::run(exponent);
-
- const Scalar zero = Scalar(0);
- const Scalar pos_one = Scalar(1);
-
- const Packet cst_zero = pset1<Packet>(zero);
- const Packet cst_pos_one = pset1<Packet>(pos_one);
-
- const Packet abs_x = pabs(x);
-
- const Packet pow_is_zero = exponent < 0 ? pcmp_lt(cst_pos_one, abs_x) : pzero(x);
- const Packet pow_is_one = pcmp_eq(cst_pos_one, abs_x);
- const Packet pow_is_neg = exponent_is_odd ? pcmp_lt(x, cst_zero) : pzero(x);
-
- Packet result = pselect(pow_is_zero, cst_zero, x);
- result = pselect(pandnot(pow_is_one, pow_is_neg), cst_pos_one, result);
- result = pselect(pand(pow_is_one, pow_is_neg), pnegate(cst_pos_one), result);
- return result;
-}
-
-template <typename Packet, typename ScalarExponent, std::enable_if_t<!NumTraits<typename unpacket_traits<Packet>::type>::IsSigned, bool> = true>
-static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet handle_int_int(const Packet& x, const ScalarExponent& exponent) {
- typedef typename unpacket_traits<Packet>::type Scalar;
-
- // unsigned integer base, integer exponent case
-
- // This routine handles negative and very large positive exponents
- // Signed integer overflow and divide by zero is undefined behavior
- // Unsigned integers do not overflow
-
- const Scalar zero = Scalar(0);
- const Scalar pos_one = Scalar(1);
-
- const Packet cst_zero = pset1<Packet>(zero);
- const Packet cst_pos_one = pset1<Packet>(pos_one);
-
- const Packet pow_is_zero = exponent < 0 ? pcmp_lt(cst_pos_one, x) : pzero(x);
- const Packet pow_is_one = pcmp_eq(cst_pos_one, x);
-
- Packet result = pselect(pow_is_zero, cst_zero, x);
+ Packet result = powx;
+ result = por(x_is_le_zero, result);
+ result = pselect(pow_is_inf_if_exp_is_neg, pand(cst_pos_inf, exp_is_neg), result);
+ result = pselect(pow_is_inf_if_exp_is_pos, pand(cst_pos_inf, exp_is_pos), result);
+ result = por(exp_is_nan, result);
result = pselect(pow_is_one, cst_pos_one, result);
return result;
}
+template <typename Packet, typename ScalarExponent,
+ std::enable_if_t<NumTraits<typename unpacket_traits<Packet>::type>::IsSigned, bool> = true>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet handle_negative_exponent(const Packet& x, const ScalarExponent& exponent) {
+ using Scalar = typename unpacket_traits<Packet>::type;
+
+ // singed integer base, signed integer exponent case
+
+ // This routine handles negative exponents.
+ // The return value is either 0, 1, or -1.
+
+ const Scalar pos_zero = Scalar(0);
+ const Scalar all_ones = ptrue<Scalar>(Scalar());
+ const Scalar pos_one = Scalar(1);
+
+ const Packet cst_pos_one = pset1<Packet>(pos_one);
+
+ const bool exponent_is_odd = exponent % ScalarExponent(2) != ScalarExponent(0);
+
+ const Packet exp_is_odd = pset1<Packet>(exponent_is_odd ? all_ones : pos_zero);
+
+ const Packet abs_x = pabs(x);
+ const Packet abs_x_is_one = pcmp_eq(abs_x, cst_pos_one);
+
+ Packet result = pselect(exp_is_odd, x, abs_x);
+ result = pand(abs_x_is_one, result);
+ return result;
+}
+
+template <typename Packet, typename ScalarExponent,
+ std::enable_if_t<!NumTraits<typename unpacket_traits<Packet>::type>::IsSigned, bool> = true>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet handle_negative_exponent(const Packet& x, const ScalarExponent&) {
+ using Scalar = typename unpacket_traits<Packet>::type;
+
+ // unsigned integer base, signed integer exponent case
+
+ // This routine handles negative exponents.
+ // The return value is either 0 or 1
+
+ const Scalar pos_one = Scalar(1);
+
+ const Packet cst_pos_one = pset1<Packet>(pos_one);
+
+ const Packet x_is_one = pcmp_eq(x, cst_pos_one);
+
+ return pand(x_is_one, x);
+}
+
+
} // end namespace unary_pow
template <typename Packet, typename ScalarExponent,
bool BaseIsIntegerType = NumTraits<typename unpacket_traits<Packet>::type>::IsInteger,
- bool ExponentIsIntegerType = NumTraits<ScalarExponent>::IsInteger>
+ bool ExponentIsIntegerType = NumTraits<ScalarExponent>::IsInteger,
+ bool ExponentIsSigned = NumTraits<ScalarExponent>::IsSigned>
struct unary_pow_impl;
-template <typename Packet, typename ScalarExponent>
-struct unary_pow_impl<Packet, ScalarExponent, false, false> {
+template <typename Packet, typename ScalarExponent, bool ExponentIsSigned>
+struct unary_pow_impl<Packet, ScalarExponent, false, false, ExponentIsSigned> {
typedef typename unpacket_traits<Packet>::type Scalar;
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent& exponent) {
const bool exponent_is_integer = (numext::isfinite)(exponent) && numext::round(exponent) == exponent;
if (exponent_is_integer) {
- Packet result = unary_pow::int_pow(x, exponent);
- result = unary_pow::handle_nonint_int_errors(x, result, exponent);
- return result;
+ return unary_pow::int_pow(x, exponent);
} else {
Packet result = unary_pow::gen_pow(x, exponent);
result = unary_pow::handle_nonint_nonint_errors(x, result, exponent);
@@ -2230,27 +2176,32 @@
}
};
-template <typename Packet, typename ScalarExponent>
-struct unary_pow_impl<Packet, ScalarExponent, false, true> {
+template <typename Packet, typename ScalarExponent, bool ExponentIsSigned>
+struct unary_pow_impl<Packet, ScalarExponent, false, true, ExponentIsSigned> {
typedef typename unpacket_traits<Packet>::type Scalar;
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent& exponent) {
- Packet result = unary_pow::int_pow(x, exponent);
- result = unary_pow::handle_nonint_int_errors(x, result, exponent);
- return result;
+ return unary_pow::int_pow(x, exponent);
}
};
template <typename Packet, typename ScalarExponent>
-struct unary_pow_impl<Packet, ScalarExponent, true, true> {
- typedef typename unpacket_traits<Packet>::type Scalar;
- static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent& exponent) {
- if (exponent < 0 || exponent > NumTraits<Scalar>::digits()) {
- return unary_pow::handle_int_int(x, exponent);
- }
- else {
- return unary_pow::int_pow(x, exponent);
- }
+struct unary_pow_impl<Packet, ScalarExponent, true, true, true> {
+ typedef typename unpacket_traits<Packet>::type Scalar;
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent& exponent) {
+ if (exponent < ScalarExponent(0)) {
+ return unary_pow::handle_negative_exponent(x, exponent);
+ } else {
+ return unary_pow::int_pow(x, exponent);
}
+ }
+};
+
+template <typename Packet, typename ScalarExponent>
+struct unary_pow_impl<Packet, ScalarExponent, true, true, false> {
+ typedef typename unpacket_traits<Packet>::type Scalar;
+ static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x, const ScalarExponent& exponent) {
+ return unary_pow::int_pow(x, exponent);
+ }
};
} // end namespace internal
diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp
index 7601125..989359f 100644
--- a/test/array_cwise.cpp
+++ b/test/array_cwise.cpp
@@ -191,53 +191,97 @@
*/
}
+template <typename Base, typename Exponent, bool ExpIsInteger = NumTraits<Exponent>::IsInteger>
+struct ref_pow {
+ static Base run(Base base, Exponent exponent) {
+ EIGEN_USING_STD(pow);
+ return static_cast<Base>(pow(base, static_cast<Base>(exponent)));
+ }
+};
-template <typename Scalar>
-void pow_scalar_exponent_test() {
- using Int_t = typename internal::make_integer<Scalar>::type;
- const Scalar tol = test_precision<Scalar>();
+template <typename Base, typename Exponent>
+struct ref_pow<Base, Exponent, true> {
+ static Base run(Base base, Exponent exponent) {
+ EIGEN_USING_STD(pow);
+ return static_cast<Base>(pow(base, exponent));
+ }
+};
- std::vector<Scalar> abs_vals = special_values<Scalar>();
- const Index num_vals = (Index)abs_vals.size();
- Map<Array<Scalar, Dynamic, 1>> bases(abs_vals.data(), num_vals);
+template <typename Exponent, bool ExpIsInteger = NumTraits<Exponent>::IsInteger>
+struct pow_helper {
+ static bool is_integer_impl(const Exponent& exp) { return (numext::isfinite)(exp) && exp == numext::floor(exp); }
+ static bool is_odd_impl(const Exponent& exp) {
+ Exponent exp_div_2 = exp / Exponent(2);
+ Exponent floor_exp_div_2 = numext::floor(exp_div_2);
+ return exp_div_2 != floor_exp_div_2;
+ }
+};
+template <typename Exponent>
+struct pow_helper<Exponent, true> {
+ static bool is_integer_impl(const Exponent&) { return true; }
+ static bool is_odd_impl(const Exponent& exp) { return exp % 2 != 0; }
+};
+template <typename Exponent>
+bool is_integer(const Exponent& exp) {
+ return pow_helper<Exponent>::is_integer_impl(exp);
+}
+template <typename Exponent>
+bool is_odd(const Exponent& exp) {
+ return pow_helper<Exponent>::is_odd_impl(exp);
+}
+template <typename Base, typename Exponent>
+void float_pow_test_impl() {
+ const Base tol = test_precision<Base>();
+ std::vector<Base> abs_base_vals = special_values<Base>();
+ std::vector<Exponent> abs_exponent_vals = special_values<Exponent>();
+ for (int i = 0; i < 100; i++) {
+ abs_base_vals.push_back(internal::random<Base>(Base(0), Base(10)));
+ abs_exponent_vals.push_back(internal::random<Exponent>(Exponent(0), Exponent(10)));
+ }
+ const Index num_repeats = internal::packet_traits<Base>::size + 1;
+ ArrayX<Base> bases(num_repeats), eigenPow(num_repeats);
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<Scalar>(NumTraits<Int_t>::highest()));
- if (exponent_is_integer) {
- Int_t exponent_as_int = static_cast<Int_t>(exponent);
- Array<Scalar, Dynamic, 1> eigenPow = bases.pow(exponent_as_int);
- for (Index j = 0; j < num_vals; j++) {
- Scalar e = static_cast<Scalar>(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));
- if ((a == a) && (e == e)) success &= (bool)numext::signbit(e) == (bool)numext::signbit(a);
- all_pass &= success;
- if (!success) {
- std::cout << "pow(" << bases(j) << "," << exponent << ") = " << a << " != " << e << std::endl;
- }
- }
- } else {
- // test floating point exponent code path
- Array<Scalar, Dynamic, 1> eigenPow = bases.pow(exponent);
- for (Index j = 0; j < num_vals; j++) {
- Scalar e = static_cast<Scalar>(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));
- if ((a == a) && (e == e)) success &= (bool)numext::signbit(e) == (bool)numext::signbit(a);
- all_pass &= success;
- if (!success) {
- std::cout << "pow(" << bases(j) << "," << exponent << ") = " << a << " != " << e << std::endl;
+ for (Base abs_base : abs_base_vals)
+ for (Base base : {negative_or_zero(abs_base), abs_base}) {
+ bases.setConstant(base);
+ for (Exponent abs_exponent : abs_exponent_vals) {
+ for (Exponent exponent : {negative_or_zero(abs_exponent), abs_exponent}) {
+ eigenPow = bases.pow(exponent);
+ for (Index j = 0; j < num_repeats; j++) {
+ Base e = ref_pow<Base, Exponent>::run(bases(j), exponent);
+ if (is_integer(exponent)) {
+ // std::pow may return an incorrect result for a very large integral exponent
+ // if base is negative and the exponent is odd, then the result must be negative
+ // if std::pow returns otherwise, flip the sign
+ bool exp_is_odd = is_odd(exponent);
+ bool base_is_neg = !(numext::isnan)(base) && (bool)numext::signbit(base);
+ bool result_is_neg = exp_is_odd && base_is_neg;
+ bool ref_is_neg = !(numext::isnan)(e) && (bool)numext::signbit(e);
+ bool flip_sign = result_is_neg != ref_is_neg;
+ if (flip_sign) e = -e;
+ }
+
+ Base a = eigenPow(j);
+ #ifdef EIGEN_COMP_MSVC
+ // Work around MSVC return value on underflow.
+ // if std::pow returns 0 and Eigen returns a denormalized value, then skip the test
+ int fpclass = std::fpclassify(a);
+ if (e == Base(0) && fpclass == FP_SUBNORMAL) continue;
+ #endif
+
+ bool both_nan = (numext::isnan)(a) && (numext::isnan)(e);
+ bool exact_or_approx = (a == e) || internal::isApprox(a, e, tol);
+ bool same_sign = (bool)numext::signbit(e) == (bool)numext::signbit(a);
+ bool success = both_nan || (exact_or_approx && same_sign);
+ all_pass &= success;
+ if (!success) {
+ std::cout << "pow(" << bases(j) << "," << exponent << ") = " << a << " != " << e << std::endl;
+ }
}
}
}
}
- }
VERIFY(all_pass);
}
@@ -259,24 +303,9 @@
}
}
-template <typename Base, typename Exponent, bool ExpIsInteger = NumTraits<Exponent>::IsInteger>
-struct ref_pow {
- static Base run(Base base, Exponent exponent) {
- EIGEN_USING_STD(pow);
- return static_cast<Base>(pow(base, static_cast<Base>(exponent)));
- }
-};
-
-template <typename Base, typename Exponent>
-struct ref_pow<Base, Exponent, true> {
- static Base run(Base base, Exponent exponent) {
- EIGEN_USING_STD(pow);
- return static_cast<Base>(pow(base, exponent));
- }
-};
-
template <typename Base, typename Exponent>
void test_exponent(Exponent exponent) {
+ EIGEN_STATIC_ASSERT(NumTraits<Base>::IsInteger,THIS TEST IS ONLY INTENDED FOR BASE INTEGER TYPES)
const Base max_abs_bases = static_cast<Base>(10000);
// avoid integer overflow in Base type
Base threshold = calc_overflow_threshold<Base, Exponent>(numext::abs(exponent));
@@ -300,10 +329,10 @@
for (Base a : y) {
Base e = ref_pow<Base, Exponent>::run(base, exponent);
bool pass = (a == e);
- if (!NumTraits<Base>::IsInteger) {
- pass = pass || (((numext::isfinite)(e) && internal::isApprox(a, e)) ||
- ((numext::isnan)(a) && (numext::isnan)(e)));
- }
+ //if (!NumTraits<Base>::IsInteger) {
+ // pass = pass || (((numext::isfinite)(e) && internal::isApprox(a, e)) ||
+ // ((numext::isnan)(a) && (numext::isnan)(e)));
+ //}
all_pass &= pass;
if (!pass) {
std::cout << "pow(" << base << "," << exponent << ") = " << a << " != " << e << std::endl;
@@ -314,7 +343,7 @@
}
template <typename Base, typename Exponent>
-void unary_pow_test() {
+void int_pow_test_impl() {
Exponent max_exponent = static_cast<Exponent>(NumTraits<Base>::digits());
Exponent min_exponent = negative_or_zero(max_exponent);
@@ -323,21 +352,26 @@
}
}
+void float_pow_test() {
+ float_pow_test_impl<float, float>();
+ float_pow_test_impl<double, double>();
+}
+
void mixed_pow_test() {
// The following cases will test promoting a smaller exponent type
// to a wider base type.
- unary_pow_test<double, int>();
- unary_pow_test<double, float>();
- unary_pow_test<float, half>();
- unary_pow_test<double, half>();
- unary_pow_test<float, bfloat16>();
- unary_pow_test<double, bfloat16>();
+ float_pow_test_impl<double, int>();
+ float_pow_test_impl<double, float>();
+ float_pow_test_impl<float, half>();
+ float_pow_test_impl<double, half>();
+ float_pow_test_impl<float, bfloat16>();
+ float_pow_test_impl<double, bfloat16>();
// Although in the following cases the exponent cannot be represented exactly
// in the base type, we do not perform a conversion, but implement
// the operation using repeated squaring.
- unary_pow_test<float, int>();
- unary_pow_test<double, long long>();
+ float_pow_test_impl<float, int>();
+ float_pow_test_impl<double, long long>();
// The following cases will test promoting a wider exponent type
// to a narrower base type. This should compile but would generate a
@@ -346,20 +380,20 @@
}
void int_pow_test() {
- unary_pow_test<int, int>();
- unary_pow_test<unsigned int, unsigned int>();
- unary_pow_test<long long, long long>();
- unary_pow_test<unsigned long long, unsigned long long>();
+ int_pow_test_impl<int, int>();
+ int_pow_test_impl<unsigned int, unsigned int>();
+ int_pow_test_impl<long long, long long>();
+ int_pow_test_impl<unsigned long long, unsigned long long>();
// Although in the following cases the exponent cannot be represented exactly
// in the base type, we do not perform a conversion, but implement the
// operation using repeated squaring.
- unary_pow_test<long long, int>();
- unary_pow_test<int, unsigned int>();
- unary_pow_test<unsigned int, int>();
- unary_pow_test<long long, unsigned long long>();
- unary_pow_test<unsigned long long, long long>();
- unary_pow_test<long long, int>();
+ int_pow_test_impl<long long, int>();
+ int_pow_test_impl<int, unsigned int>();
+ int_pow_test_impl<unsigned int, int>();
+ int_pow_test_impl<long long, unsigned long long>();
+ int_pow_test_impl<unsigned long long, long long>();
+ int_pow_test_impl<long long, int>();
}
namespace Eigen {
@@ -849,7 +883,6 @@
// Test pow and atan2 on special IEEE values.
unary_ops_test<Scalar>();
binary_ops_test<Scalar>();
- pow_scalar_exponent_test<Scalar>();
VERIFY_IS_APPROX(log10(m3), log(m3)/numext::log(Scalar(10)));
VERIFY_IS_APPROX(log2(m3), log(m3)/numext::log(Scalar(2)));
@@ -1223,6 +1256,7 @@
}
for(int i = 0; i < g_repeat; i++) {
+ CALL_SUBTEST_5( float_pow_test() );
CALL_SUBTEST_6( int_pow_test() );
CALL_SUBTEST_7( mixed_pow_test() );
CALL_SUBTEST_8( signbit_tests() );
