diff --git a/Eigen/src/Core/arch/AVX/Complex.h b/Eigen/src/Core/arch/AVX/Complex.h index a5e6499c5..6a8bee890 100644 --- a/Eigen/src/Core/arch/AVX/Complex.h +++ b/Eigen/src/Core/arch/AVX/Complex.h @@ -40,6 +40,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -238,6 +239,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -431,6 +433,16 @@ EIGEN_STRONG_INLINE Packet4cf psqrt(const Packet4cf& a) { return psqrt_complex(a); } +template <> +EIGEN_STRONG_INLINE Packet2cd plog(const Packet2cd& a) { + return plog_complex(a); +} + +template <> +EIGEN_STRONG_INLINE Packet4cf plog(const Packet4cf& a) { + return plog_complex(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 f2c8ce630..c14b4a0b6 100644 --- a/Eigen/src/Core/arch/AVX512/Complex.h +++ b/Eigen/src/Core/arch/AVX512/Complex.h @@ -39,6 +39,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -223,6 +224,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -448,6 +450,16 @@ EIGEN_STRONG_INLINE Packet8cf psqrt(const Packet8cf& a) { return psqrt_complex(a); } +template <> +EIGEN_STRONG_INLINE Packet4cd plog(const Packet4cd& a) { + return plog_complex(a); +} + +template <> +EIGEN_STRONG_INLINE Packet8cf plog(const Packet8cf& a) { + return plog_complex(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 7bfc61d21..e3c4436c5 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h @@ -98,6 +98,7 @@ struct packet_traits > : default_packet_traits { HasMin = 0, HasMax = 0, HasSqrt = 1, + HasLog = 1, #ifdef EIGEN_VECTORIZE_VSX HasBlend = 1, #endif @@ -369,6 +370,11 @@ EIGEN_STRONG_INLINE Packet2cf psqrt(const Packet2cf& a) { return psqrt_complex(a); } +template <> +EIGEN_STRONG_INLINE Packet2cf plog(const Packet2cf& a) { + return plog_complex(a); +} + //---------- double ---------- #ifdef EIGEN_VECTORIZE_VSX struct Packet1cd { @@ -433,6 +439,7 @@ struct packet_traits > : default_packet_traits { HasMin = 0, HasMax = 0, HasSqrt = 1, + HasLog = 1, HasSetLinear = 0 }; }; @@ -618,6 +625,11 @@ EIGEN_STRONG_INLINE Packet1cd psqrt(const Packet1cd& a) { return psqrt_complex(a); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog(const Packet1cd& a) { + return plog_complex(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 d84b1cc59..118426fbe 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h @@ -964,6 +964,34 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pdiv_complex(const Pa return Packet(pdiv(result_scaled.v, y_max)); } +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_complex(const Packet& x) { + typedef typename unpacket_traits::type Scalar; + typedef typename Scalar::value_type RealScalar; + typedef typename unpacket_traits::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(NumTraits::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 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psqrt_complex(const Packet& a) { typedef typename unpacket_traits::type Scalar; @@ -1076,6 +1104,41 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psqrt_complex(const P 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 +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet phypot_complex(const Packet& a) { + typedef typename unpacket_traits::type Scalar; + typedef typename Scalar::value_type RealScalar; + typedef typename unpacket_traits::as_real RealPacket; + + const RealPacket cst_zero_rp = pset1(static_cast(0.0)); + const RealPacket cst_minus_one_rp = pset1(static_cast(-1.0)); + const RealPacket cst_two_rp = pset1(static_cast(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 struct psign_impl::type>::IsComplex && !NumTraits::type>::IsInteger>> { diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h index ade9f3f41..960bb67fb 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h @@ -113,6 +113,10 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pdiv_complex(const Pa template struct ppolevl; +/** \internal \returns log(x) for complex types */ +template +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 82408474d..22c776574 100644 --- a/Eigen/src/Core/arch/NEON/Complex.h +++ b/Eigen/src/Core/arch/NEON/Complex.h @@ -62,6 +62,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -436,6 +437,16 @@ EIGEN_STRONG_INLINE Packet2cf psqrt(const Packet2cf& a) { return psqrt_complex(a); } +template <> +EIGEN_STRONG_INLINE Packet1cf plog(const Packet1cf& a) { + return plog_complex(a); +} + +template <> +EIGEN_STRONG_INLINE Packet2cf plog(const Packet2cf& a) { + return plog_complex(a); +} + //---------- double ---------- #if EIGEN_ARCH_ARM64 && !EIGEN_APPLE_DOUBLE_NEON_BUG @@ -468,6 +479,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -657,6 +669,11 @@ EIGEN_STRONG_INLINE Packet1cd psqrt(const Packet1cd& a) { return psqrt_complex(a); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog(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 4c5c499e1..76c3a0547 100644 --- a/Eigen/src/Core/arch/SSE/Complex.h +++ b/Eigen/src/Core/arch/SSE/Complex.h @@ -42,6 +42,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -232,6 +233,7 @@ struct packet_traits > : default_packet_traits { HasDiv = 1, HasNegate = 1, HasSqrt = 1, + HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -412,6 +414,16 @@ EIGEN_STRONG_INLINE Packet2cf psqrt(const Packet2cf& a) { return psqrt_complex(a); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog(const Packet1cd& a) { + return plog_complex(a); +} + +template <> +EIGEN_STRONG_INLINE Packet2cf plog(const Packet2cf& a) { + return plog_complex(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 4000e0567..e8bd17da1 100644 --- a/Eigen/src/Core/arch/ZVector/Complex.h +++ b/Eigen/src/Core/arch/ZVector/Complex.h @@ -60,6 +60,7 @@ struct packet_traits > : default_packet_traits { HasSub = 1, HasMul = 1, HasDiv = 1, + HasLog = 1, HasNegate = 1, HasAbs = 0, HasAbs2 = 0, @@ -83,6 +84,7 @@ struct packet_traits > : default_packet_traits { HasSub = 1, HasMul = 1, HasDiv = 1, + HasLog = 1, HasNegate = 1, HasAbs = 0, HasAbs2 = 0, @@ -248,6 +250,11 @@ EIGEN_STRONG_INLINE Packet1cd pdiv(const Packet1cd& a, const Packet1c return pdiv_complex(a, b); } +template <> +EIGEN_STRONG_INLINE Packet1cd plog(const Packet1cd& a, const Packet1cd& b) { + return plog_complex(a, b); +} + EIGEN_STRONG_INLINE Packet1cd pcplxflip /**/ (const Packet1cd& x) { return Packet1cd(preverse(Packet2d(x.v))); } @@ -424,6 +431,11 @@ EIGEN_STRONG_INLINE Packet2cf pdiv(const Packet2cf& a, const Packet2c return pdiv_complex(a, b); } +template <> +EIGEN_STRONG_INLINE Packet2cf plog(const Packet2cf& a, const Packet2cf& b) { + return plog_complex(a, b); +} + EIGEN_STRONG_INLINE Packet2cf pcplxflip /**/ (const Packet2cf& x) { Packet2cf res; res.cd[0] = pcplxflip(x.cd[0]); diff --git a/test/packetmath.cpp b/test/packetmath.cpp index a97cce77c..b2cca731e 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -1340,6 +1340,8 @@ void packetmath_complex() { 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(1e2); @@ -1401,6 +1403,47 @@ void packetmath_complex() { 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(), internal::random()); + } + 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::infinity(); + const RealScalar nan = std::numeric_limits::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(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 diff --git a/test/packetmath_test_shared.h b/test/packetmath_test_shared.h index e1acb684d..a4bb34766 100644 --- a/test/packetmath_test_shared.h +++ b/test/packetmath_test_shared.h @@ -124,6 +124,19 @@ bool areEqual(const Scalar* a, const Scalar* b, int size) { return true; } +template +bool areApprox(const Scalar* a, const Scalar* b, int size, const typename NumTraits::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 @@ bool areEqual(const Scalar* a, const Scalar* b, int size) { 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::epsilon(); \ + for (int j = 0; j < N; j += PacketSize) \ + internal::pstore(data2 + j, internal::plog(internal::pload(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 struct packet_helper { template