Fix pexp complex test edge-cases.
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index bb553e7..3f28068 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -164,7 +164,7 @@
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC static inline RealScalar& run(Scalar& x) { return reinterpret_cast<RealScalar*>(&x)[1]; }
EIGEN_DEVICE_FUNC static inline const RealScalar& run(const Scalar& x) {
- return reinterpret_cast<RealScalar*>(&x)[1];
+ return reinterpret_cast<const RealScalar*>(&x)[1];
}
};
@@ -1541,6 +1541,25 @@
return exp(x);
}
+// MSVC screws up some edge-cases for std::exp(complex).
+#ifdef EIGEN_COMP_MSVC
+template <typename RealScalar>
+EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex<RealScalar> exp(const std::complex<RealScalar>& x) {
+ EIGEN_USING_STD(exp);
+ // If z is (x,±∞) (for any finite x), the result is (NaN,NaN) and FE_INVALID is raised.
+ // If z is (x,NaN) (for any finite x), the result is (NaN,NaN) and FE_INVALID may be raised.
+ if ((isfinite)(real_ref(x)) && !(isfinite)(imag_ref(x))) {
+ return std::complex<RealScalar>(NumTraits<RealScalar>::quiet_NaN(), NumTraits<RealScalar>::quiet_NaN());
+ }
+ // If z is (+∞,±∞), the result is (±∞,NaN) and FE_INVALID is raised (the sign of the real part is unspecified)
+ // If z is (+∞,NaN), the result is (±∞,NaN) (the sign of the real part is unspecified)
+ if ((real_ref(x) == NumTraits<RealScalar>::infinity() && !(isfinite)(imag_ref(x)))) {
+ return std::complex<RealScalar>(NumTraits<RealScalar>::infinity(), NumTraits<RealScalar>::quiet_NaN());
+ }
+ return exp(x);
+}
+#endif
+
#if defined(SYCL_DEVICE_ONLY)
SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(exp, exp)
#endif
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
index b44cb78..78dbf20 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
@@ -1068,40 +1068,39 @@
typedef typename unpacket_traits<Packet>::type Scalar;
typedef typename Scalar::value_type RealScalar;
const RealPacket even_mask = peven_mask(a.v);
- const Packet even_maskp = Packet(even_mask);
const RealPacket odd_mask = pcplxflip(Packet(even_mask)).v;
- Packet p0y = Packet(pand(odd_mask, a.v));
- Packet py0 = pcplxflip(p0y);
- Packet pyy = padd(p0y, py0);
+ // Let a = x + iy.
+ // exp(a) = exp(x) * cis(y), plus some special edge-case handling.
- RealPacket sincos = psincos_float<false, RealPacket, true>(pyy.v);
- RealPacket cossin = pcplxflip(Packet(sincos)).v;
+ // exp(x):
+ RealPacket x = pand(a.v, even_mask);
+ x = por(x, pcplxflip(Packet(x)).v);
+ RealPacket expx = pexp(x); // exp(x);
+
+ // cis(y):
+ RealPacket y = pand(odd_mask, a.v);
+ y = por(y, pcplxflip(Packet(y)).v);
+ RealPacket cisy = psincos_float<false, RealPacket, true>(y);
+ cisy = pcplxflip(Packet(cisy)).v; // cos(y) + i * sin(y)
const RealPacket cst_pos_inf = pset1<RealPacket>(NumTraits<RealScalar>::infinity());
const RealPacket cst_neg_inf = pset1<RealPacket>(-NumTraits<RealScalar>::infinity());
- Packet x_is_inf = Packet(pcmp_eq(a.v, cst_pos_inf));
- Packet x_is_minf = Packet(pcmp_eq(a.v, cst_neg_inf));
- Packet x_is_zero = Packet(pcmp_eq(pzero(a).v, a.v));
- Packet x_real_is_inf = pand(even_maskp, x_is_inf);
- Packet x_real_is_minf = pand(even_maskp, x_is_minf);
- Packet inf0 = pset1<Packet>(Scalar(NumTraits<RealScalar>::infinity(), RealScalar(0)));
- Packet x_is_inf0 = pand(x_real_is_inf, pcplxflip(x_is_zero));
- x_is_inf0 = por(x_is_inf0, pcplxflip(x_is_inf0));
- Packet x_imag_goes_zero = pand(por(x_is_minf, x_is_inf), pcplxflip(x_real_is_minf));
- Packet x_is_nan = Packet(pisnan(a.v));
- Packet x_real_goes_zero = pand(x_is_nan, pcplxflip(x_real_is_minf));
- RealPacket pexp_real = pexp(a.v);
- Packet pexp_half = Packet(pand(even_mask, pexp_real));
- RealPacket xexp_flip_rp = pcplxflip(pexp_half).v;
- RealPacket xexp = padd(pexp_half.v, xexp_flip_rp);
- Packet result(pmul(cossin, xexp));
+ // If x is -inf, we know that cossin(y) is bounded,
+ // so the result is (0, +/-0), where the sign of the imaginary part comes
+ // from the sign of cossin(y).
+ RealPacket cisy_sign = por(pandnot(cisy, pabs(cisy)), pset1<RealPacket>(RealScalar(1)));
+ cisy = pselect(pcmp_eq(x, cst_neg_inf), cisy_sign, cisy);
- result = pselect(x_is_inf0, inf0, result);
- result = pselect(x_real_is_minf, pzero(a), result);
- result = pselect(x_imag_goes_zero, pzero(a), result);
- result = pselect(x_real_goes_zero, pzero(a), result);
+ // If x is inf, and cos(y) has unknown sign (y is inf or NaN), the result
+ // is (+/-inf, NaN), where the signs are undetermined (take the sign of y).
+ RealPacket y_sign = por(pandnot(y, pabs(y)), pset1<RealPacket>(RealScalar(1)));
+ cisy = pselect(pand(pcmp_eq(x, cst_pos_inf), pisnan(cisy)), pand(y_sign, even_mask), cisy);
+ Packet result = Packet(pmul(expx, cisy));
+
+ // If y is +/- 0, the input is real, so take the real result for consistency.
+ result = pselect(Packet(pcmp_eq(y, pzero(y))), Packet(por(pand(expx, even_mask), pand(y, odd_mask))), result);
return result;
}
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index c05f873..92e72cb 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -277,6 +277,7 @@
template <typename Scalar, typename Packet>
void packetmath_boolean_mask_ops() {
+ using RealScalar = typename NumTraits<Scalar>::Real;
const int PacketSize = internal::unpacket_traits<Packet>::size;
const int size = 2 * PacketSize;
EIGEN_ALIGN_MAX Scalar data1[size];
@@ -289,7 +290,7 @@
CHECK_CWISE1(internal::ptrue, internal::ptrue);
CHECK_CWISE2_IF(true, internal::pandnot, internal::pandnot);
for (int i = 0; i < PacketSize; ++i) {
- data1[i] = Scalar(i);
+ data1[i] = Scalar(RealScalar(i));
data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
}
@@ -1335,6 +1336,62 @@
template <typename Scalar, typename Packet, bool HasExp = internal::packet_traits<Scalar>::HasExp>
struct exp_complex_test_impl {
typedef typename Scalar::value_type RealScalar;
+
+ static Scalar pexp1(const Scalar& x) {
+ Packet px = internal::pset1<Packet>(x);
+ Packet py = internal::pexp(px);
+ return internal::pfirst(py);
+ }
+
+ static Scalar cis(const RealScalar& x) { return Scalar(numext::cos(x), numext::sin(x)); }
+
+ // Verify equality with signed zero.
+ static bool is_exactly_equal(const RealScalar& a, const RealScalar& b) {
+ // NaNs are always unsigned, and always compare not equal directly.
+ if ((numext::isnan)(a)) {
+ return (numext::isnan)(b);
+ }
+ // Signed zero.
+ RealScalar zero(0);
+ if (a == zero) {
+ // Signs are either 0 or NaN, so verify that their comparisons to zero are equal.
+ return (a == b) && ((numext::signbit(a) == zero) == (numext::signbit(b) == zero));
+ }
+ // Allow _some_ tolerance.
+ return verifyIsApprox(a, b);
+ }
+
+ // Verify equality with signed zero.
+ static bool is_exactly_equal(const Scalar& a, const Scalar& b) {
+ bool result = is_exactly_equal(numext::real_ref(a), numext::real_ref(b)) &&
+ is_exactly_equal(numext::imag_ref(a), numext::imag_ref(b));
+ if (!result) {
+ std::cout << a << " != " << b << std::endl;
+ }
+ return result;
+ }
+
+ static bool is_sign_exp_unspecified(const Scalar& z) {
+ const RealScalar inf = std::numeric_limits<RealScalar>::infinity();
+ // If z is (-∞,±∞), the result is (±0,±0) (signs are unspecified)
+ if (numext::real_ref(z) == -inf && (numext::isinf)(numext::imag_ref(z))) {
+ return true;
+ }
+ // If z is (+∞,±∞), the result is (±∞,NaN) and FE_INVALID is raised (the sign of the real part is unspecified)
+ if (numext::real_ref(z) == +inf && (numext::isinf)(numext::imag_ref(z))) {
+ return true;
+ }
+ // If z is (-∞,NaN), the result is (±0,±0) (signs are unspecified)
+ if (numext::real_ref(z) == -inf && (numext::isnan)(numext::imag_ref(z))) {
+ return true;
+ }
+ // If z is (+∞,NaN), the result is (±∞,NaN) (the sign of the real part is unspecified)
+ if (numext::real_ref(z) == +inf && (numext::isnan)(numext::imag_ref(z))) {
+ return true;
+ }
+ return false;
+ }
+
static void run(Scalar* data1, Scalar* data2, Scalar* ref, int size) {
const int PacketSize = internal::unpacket_traits<Packet>::size;
@@ -1343,27 +1400,45 @@
}
CHECK_CWISE1_N(std::exp, internal::pexp, size);
- // Test misc. corner cases.
- const RealScalar zero = RealScalar(0);
- const RealScalar one = RealScalar(1);
- const RealScalar inf = std::numeric_limits<RealScalar>::infinity();
- const RealScalar nan = std::numeric_limits<RealScalar>::quiet_NaN();
- for (RealScalar x : {zero, one, inf}) {
- for (RealScalar y : {zero, one, inf}) {
- data1[0] = Scalar(x, y);
- data1[1] = Scalar(-x, y);
- data1[2] = Scalar(x, -y);
- data1[3] = Scalar(-x, -y);
- CHECK_CWISE1_N(std::exp, internal::pexp, 4);
+ // Test all corner cases (and more).
+ const RealScalar edges[] = {RealScalar(0),
+ RealScalar(1),
+ RealScalar(2),
+ RealScalar(EIGEN_PI / 2),
+ RealScalar(EIGEN_PI),
+ RealScalar(3 * EIGEN_PI / 2),
+ RealScalar(2 * EIGEN_PI),
+ numext::log(NumTraits<RealScalar>::highest()) - 1,
+ NumTraits<RealScalar>::highest(),
+ std::numeric_limits<RealScalar>::infinity(),
+ std::numeric_limits<RealScalar>::quiet_NaN(),
+ -RealScalar(0),
+ -RealScalar(1),
+ -RealScalar(2),
+ -RealScalar(EIGEN_PI / 2),
+ -RealScalar(EIGEN_PI),
+ -RealScalar(3 * EIGEN_PI / 2),
+ -RealScalar(2 * EIGEN_PI),
+ -numext::log(NumTraits<RealScalar>::highest()) + 1,
+ -NumTraits<RealScalar>::highest(),
+ -std::numeric_limits<RealScalar>::infinity(),
+ -std::numeric_limits<RealScalar>::quiet_NaN()};
+
+ for (RealScalar x : edges) {
+ for (RealScalar y : edges) {
+ Scalar z = Scalar(x, y);
+ Scalar w = pexp1(z);
+ if (is_sign_exp_unspecified(z)) {
+ Scalar abs_w = Scalar(numext::abs(numext::real_ref(w)), numext::abs(numext::imag_ref(w)));
+ Scalar expected = numext::exp(z);
+ Scalar abs_expected =
+ Scalar(numext::abs(numext::real_ref(expected)), numext::abs(numext::imag_ref(expected)));
+ VERIFY(is_exactly_equal(abs_w, abs_expected));
+ } else {
+ VERIFY(is_exactly_equal(w, numext::exp(z)));
+ }
}
}
- for (RealScalar x : {zero, one, inf}) {
- data1[0] = Scalar(x, nan);
- data1[1] = Scalar(-x, nan);
- data1[2] = Scalar(nan, x);
- data1[3] = Scalar(nan, -x);
- CHECK_CWISE1_N(std::exp, internal::pexp, 4);
- }
}
};
