Implement plog_complex
diff --git a/Eigen/src/Core/arch/AVX/Complex.h b/Eigen/src/Core/arch/AVX/Complex.h index a5e6499..6a8bee8 100644 --- a/Eigen/src/Core/arch/AVX/Complex.h +++ b/Eigen/src/Core/arch/AVX/Complex.h
@@ -40,6 +40,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -238,6 +239,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -431,6 +433,16 @@ return psqrt_complex<Packet4cf>(a); } +template <> +EIGEN_STRONG_INLINE Packet2cd plog<Packet2cd>(const Packet2cd& a) { + return plog_complex<Packet2cd>(a); +} + +template <> +EIGEN_STRONG_INLINE Packet4cf plog<Packet4cf>(const Packet4cf& a) { + return plog_complex<Packet4cf>(a); +} + } // end namespace internal } // end namespace Eigen
diff --git a/Eigen/src/Core/arch/AVX512/Complex.h b/Eigen/src/Core/arch/AVX512/Complex.h index f2c8ce6..c14b4a0 100644 --- a/Eigen/src/Core/arch/AVX512/Complex.h +++ b/Eigen/src/Core/arch/AVX512/Complex.h
@@ -39,6 +39,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -223,6 +224,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -448,6 +450,16 @@ return psqrt_complex<Packet8cf>(a); } +template <> +EIGEN_STRONG_INLINE Packet4cd plog<Packet4cd>(const Packet4cd& a) { + return plog_complex<Packet4cd>(a); +} + +template <> +EIGEN_STRONG_INLINE Packet8cf plog<Packet8cf>(const Packet8cf& a) { + return plog_complex<Packet8cf>(a); +} + } // end namespace internal } // end namespace Eigen
diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h index 7bfc61d..e3c4436 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h
@@ -98,6 +98,7 @@ HasMin = 0, HasMax = 0, HasSqrt = 1, + HasLog = 1, #ifdef EIGEN_VECTORIZE_VSX HasBlend = 1, #endif @@ -369,6 +370,11 @@ return psqrt_complex<Packet2cf>(a); } +template <> +EIGEN_STRONG_INLINE Packet2cf plog<Packet2cf>(const Packet2cf& a) { + return plog_complex<Packet2cf>(a); +} + //---------- double ---------- #ifdef EIGEN_VECTORIZE_VSX struct Packet1cd { @@ -433,6 +439,7 @@ HasMin = 0, HasMax = 0, HasSqrt = 1, + HasLog = 1, HasSetLinear = 0 }; }; @@ -618,6 +625,11 @@ return psqrt_complex<Packet1cd>(a); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog<Packet1cd>(const Packet1cd& a) { + return plog_complex<Packet1cd>(a); +} + #endif // __VSX__ } // end namespace internal
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index d84b1cc..118426f 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
@@ -965,6 +965,34 @@ } template <typename Packet> +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_complex(const Packet& x) { + typedef typename unpacket_traits<Packet>::type Scalar; + typedef typename Scalar::value_type RealScalar; + typedef typename unpacket_traits<Packet>::as_real RealPacket; + + RealPacket real_mask_rp = peven_mask(x.v); + Packet real_mask(real_mask_rp); + + // Real part + RealPacket x_flip = pcplxflip(x).v; // b, a + Packet x_norm = phypot_complex(x); // sqrt(a^2 + b^2), sqrt(a^2 + b^2) + RealPacket xlogr = plog(x_norm.v); // log(sqrt(a^2 + b^2)), log(sqrt(a^2 + b^2)) + + // Imag part + RealPacket ximg = patan2(x.v, x_flip); // atan2(a, b), atan2(b, a) + + const RealPacket cst_pos_inf = pset1<RealPacket>(NumTraits<RealScalar>::infinity()); + RealPacket x_abs = pabs(x.v); + RealPacket is_x_pos_inf = pcmp_eq(x_abs, cst_pos_inf); + RealPacket is_y_pos_inf = pcplxflip(Packet(is_x_pos_inf)).v; + RealPacket is_any_inf = por(is_x_pos_inf, is_y_pos_inf); + RealPacket xreal = pselect(is_any_inf, cst_pos_inf, xlogr); + + Packet xres = pselect(real_mask, Packet(xreal), Packet(ximg)); // log(sqrt(a^2 + b^2)), atan2(b, a) + return xres; +} + +template <typename Packet> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psqrt_complex(const Packet& a) { typedef typename unpacket_traits<Packet>::type Scalar; typedef typename Scalar::value_type RealScalar; @@ -1076,6 +1104,41 @@ return pselect(is_imag_inf, imag_inf_result, pselect(is_real_inf, real_inf_result, result)); } +// \internal \returns the norm of a complex number z = x + i*y, defined as sqrt(x^2 + y^2). +// Implemented using the hypot(a,b) algorithm from https://doi.org/10.48550/arXiv.1904.09481 +template <typename Packet> +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet phypot_complex(const Packet& a) { + typedef typename unpacket_traits<Packet>::type Scalar; + typedef typename Scalar::value_type RealScalar; + typedef typename unpacket_traits<Packet>::as_real RealPacket; + + const RealPacket cst_zero_rp = pset1<RealPacket>(static_cast<RealScalar>(0.0)); + const RealPacket cst_minus_one_rp = pset1<RealPacket>(static_cast<RealScalar>(-1.0)); + const RealPacket cst_two_rp = pset1<RealPacket>(static_cast<RealScalar>(2.0)); + const RealPacket evenmask = peven_mask(a.v); + + RealPacket a_abs = pabs(a.v); + RealPacket a_flip = pcplxflip(Packet(a_abs)).v; // |b|, |a| + RealPacket a_all = pselect(evenmask, a_abs, a_flip); // |a|, |a| + RealPacket b_all = pselect(evenmask, a_flip, a_abs); // |b|, |b| + + RealPacket a2 = pmul(a.v, a.v); // |a^2, b^2| + RealPacket a2_flip = pcplxflip(Packet(a2)).v; // |b^2, a^2| + RealPacket h = psqrt(padd(a2, a2_flip)); // |sqrt(a^2 + b^2), sqrt(a^2 + b^2)| + RealPacket h_sq = pmul(h, h); // |a^2 + b^2, a^2 + b^2| + RealPacket a_sq = pselect(evenmask, a2, a2_flip); // |a^2, a^2| + RealPacket m_h_sq = pmul(h_sq, cst_minus_one_rp); + RealPacket m_a_sq = pmul(a_sq, cst_minus_one_rp); + RealPacket x = psub(psub(pmadd(h, h, m_h_sq), pmadd(b_all, b_all, psub(a_sq, h_sq))), pmadd(a_all, a_all, m_a_sq)); + h = psub(h, pdiv(x, pmul(cst_two_rp, h))); // |h - x/(2*h), h - x/(2*h)| + + // handle zero-case + RealPacket iszero = pcmp_eq(por(a_abs, a_flip), cst_zero_rp); + + h = pandnot(h, iszero); // |sqrt(a^2+b^2), sqrt(a^2+b^2)| + return Packet(h); // |sqrt(a^2+b^2), sqrt(a^2+b^2)| +} + template <typename Packet> struct psign_impl<Packet, std::enable_if_t<!NumTraits<typename unpacket_traits<Packet>::type>::IsComplex && !NumTraits<typename unpacket_traits<Packet>::type>::IsInteger>> {
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h index ade9f3f..960bb67 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h
@@ -113,6 +113,10 @@ template <typename Packet, int N> struct ppolevl; +/** \internal \returns log(x) for complex types */ +template <typename Packet> +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_complex(const Packet& x); + // Macros for instantiating these generic functions for different backends. #define EIGEN_PACKET_FUNCTION(METHOD, SCALAR, PACKET) \ template <> \
diff --git a/Eigen/src/Core/arch/NEON/Complex.h b/Eigen/src/Core/arch/NEON/Complex.h index 8240847..22c7765 100644 --- a/Eigen/src/Core/arch/NEON/Complex.h +++ b/Eigen/src/Core/arch/NEON/Complex.h
@@ -62,6 +62,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -436,6 +437,16 @@ return psqrt_complex<Packet2cf>(a); } +template <> +EIGEN_STRONG_INLINE Packet1cf plog<Packet1cf>(const Packet1cf& a) { + return plog_complex(a); +} + +template <> +EIGEN_STRONG_INLINE Packet2cf plog<Packet2cf>(const Packet2cf& a) { + return plog_complex(a); +} + //---------- double ---------- #if EIGEN_ARCH_ARM64 && !EIGEN_APPLE_DOUBLE_NEON_BUG @@ -468,6 +479,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -657,6 +669,11 @@ return psqrt_complex<Packet1cd>(a); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog<Packet1cd>(const Packet1cd& a) { + return plog_complex(a); +} + #endif // EIGEN_ARCH_ARM64 } // end namespace internal
diff --git a/Eigen/src/Core/arch/SSE/Complex.h b/Eigen/src/Core/arch/SSE/Complex.h index 4c5c499..76c3a05 100644 --- a/Eigen/src/Core/arch/SSE/Complex.h +++ b/Eigen/src/Core/arch/SSE/Complex.h
@@ -42,6 +42,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -232,6 +233,7 @@ HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -412,6 +414,16 @@ return psqrt_complex<Packet2cf>(a); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog<Packet1cd>(const Packet1cd& a) { + return plog_complex<Packet1cd>(a); +} + +template <> +EIGEN_STRONG_INLINE Packet2cf plog<Packet2cf>(const Packet2cf& a) { + return plog_complex<Packet2cf>(a); +} + } // end namespace internal } // end namespace Eigen
diff --git a/Eigen/src/Core/arch/ZVector/Complex.h b/Eigen/src/Core/arch/ZVector/Complex.h index 4000e05..e8bd17d 100644 --- a/Eigen/src/Core/arch/ZVector/Complex.h +++ b/Eigen/src/Core/arch/ZVector/Complex.h
@@ -60,6 +60,7 @@ HasSub = 1, HasMul = 1, HasDiv = 1, + HasLog = 1, HasNegate = 1, HasAbs = 0, HasAbs2 = 0, @@ -83,6 +84,7 @@ HasSub = 1, HasMul = 1, HasDiv = 1, + HasLog = 1, HasNegate = 1, HasAbs = 0, HasAbs2 = 0, @@ -248,6 +250,11 @@ return pdiv_complex(a, b); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog<Packet1cd>(const Packet1cd& a, const Packet1cd& b) { + return plog_complex(a, b); +} + EIGEN_STRONG_INLINE Packet1cd pcplxflip /*<Packet1cd>*/ (const Packet1cd& x) { return Packet1cd(preverse(Packet2d(x.v))); } @@ -424,6 +431,11 @@ return pdiv_complex(a, b); } +template <> +EIGEN_STRONG_INLINE Packet2cf plog<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { + return plog_complex(a, b); +} + EIGEN_STRONG_INLINE Packet2cf pcplxflip /*<Packet2cf>*/ (const Packet2cf& x) { Packet2cf res; res.cd[0] = pcplxflip(x.cd[0]);
diff --git a/test/packetmath.cpp b/test/packetmath.cpp index a97cce7..b2cca73 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp
@@ -1340,6 +1340,8 @@ EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4]; EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4]; EIGEN_ALIGN_MAX Scalar pval[PacketSize * 4]; + EIGEN_ALIGN_MAX RealScalar realdata[PacketSize * 4]; + EIGEN_ALIGN_MAX RealScalar realref[PacketSize * 4]; for (int i = 0; i < size; ++i) { data1[i] = internal::random<Scalar>() * Scalar(1e2); @@ -1401,6 +1403,47 @@ data1[3] = Scalar(-inf, nan); CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4); } + if (PacketTraits::HasLog) { + for (int i = 0; i < size; ++i) { + data1[i] = Scalar(internal::random<RealScalar>(), internal::random<RealScalar>()); + } + CHECK_CWISE1_N(std::log, internal::plog, 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_IM1ULP_N(std::log, internal::plog, 4); + } + } + // Set reference results to nan. + // Some architectures don't handle IEEE edge cases correctly + ref[0] = Scalar(nan, nan); + ref[1] = Scalar(nan, nan); + ref[2] = Scalar(nan, nan); + ref[3] = Scalar(nan, nan); + for (RealScalar x : {zero, one}) { + data1[0] = Scalar(x, nan); + data1[1] = Scalar(-x, nan); + data1[2] = Scalar(nan, x); + data1[3] = Scalar(nan, -x); + for (int j = 0; j < size; j += PacketSize) + internal::pstore(data2 + j, internal::plog(internal::pload<Packet>(data1 + j))); + VERIFY(test::areApprox(ref, data2, 4)); + } + data1[0] = Scalar(inf, nan); + data1[1] = Scalar(-inf, nan); + data1[2] = Scalar(nan, inf); + data1[3] = Scalar(nan, -inf); + CHECK_CWISE1_IM1ULP_N(std::log, internal::plog, 4); + } } template <typename Scalar, typename Packet>
diff --git a/test/packetmath_test_shared.h b/test/packetmath_test_shared.h index e1acb68..a4bb347 100644 --- a/test/packetmath_test_shared.h +++ b/test/packetmath_test_shared.h
@@ -124,6 +124,19 @@ return true; } +template <typename Scalar> +bool areApprox(const Scalar* a, const Scalar* b, int size, const typename NumTraits<Scalar>::Real& precision) { + for (int i = 0; i < size; ++i) { + if (numext::not_equal_strict(a[i], b[i]) && !internal::isApprox(a[i], b[i], precision) && + !((numext::isnan)(a[i]) && (numext::isnan)(b[i]))) { + print_mismatch(a, b, size); + std::cout << "Values differ in position " << i << ": " << a[i] << " vs " << b[i] << std::endl; + return false; + } + } + return true; +} + #define CHECK_CWISE1(REFOP, POP) \ { \ for (int i = 0; i < PacketSize; ++i) ref[i] = REFOP(data1[i]); \ @@ -141,6 +154,29 @@ VERIFY(test::areApprox(ref, data2, N) && #POP); \ } +// Checks component-wise for input of complex type of size N. The real and +// the imaginary part are compared separately, with 1ULP relaxed condition +// for the imaginary part. All of data1 data2, ref, realdata1 and realref +// should have size at least ceil(N/PacketSize)*PacketSize to avoid +// memory access errors. +#define CHECK_CWISE1_IM1ULP_N(REFOP, POP, N) \ + { \ + RealScalar eps_1ulp = RealScalar(1e1) * std::numeric_limits<RealScalar>::epsilon(); \ + for (int j = 0; j < N; j += PacketSize) \ + internal::pstore(data2 + j, internal::plog(internal::pload<Packet>(data1 + j))); \ + for (int i = 0; i < N; ++i) { \ + ref[i] = REFOP(data1[i]); \ + realref[i] = ref[i].imag(); \ + realdata[i] = data2[i].imag(); \ + } \ + VERIFY(test::areApprox(realdata, realref, N, eps_1ulp)); \ + for (int i = 0; i < N; ++i) { \ + realdata[i] = data2[i].real(); \ + realref[i] = ref[i].real(); \ + } \ + VERIFY(test::areApprox(realdata, realref, N)); \ + } + template <bool Cond, typename Packet> struct packet_helper { template <typename T>