GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-10 18:59:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implement plog_complex

Browse Source
This commit is contained in:
Damiano Franzò 2023-12-07 13:41:09 +01:00 committed by Rasmus Munk Larsen
parent 043442e21b
commit 7fd7a3f946
10 changed files with 223 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Eigen/src/Core/arch/AVX/Complex.h
View File

@ -40,6 +40,7 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
HasDiv = 1,
HasNegate = 1,
HasSqrt = 1,
HasLog = 1,
HasAbs = 0,
HasAbs2 = 0,
HasMin = 0,
@ -238,6 +239,7 @@ struct packet_traits<std::complex<double> > : 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<Packet4cf>(const Packet4cf& a) {
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

12
Eigen/src/Core/arch/AVX512/Complex.h
View File

@ -39,6 +39,7 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
HasDiv = 1,
HasNegate = 1,
HasSqrt = 1,
HasLog = 1,
HasAbs = 0,
HasAbs2 = 0,
HasMin = 0,
@ -223,6 +224,7 @@ struct packet_traits<std::complex<double> > : 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<Packet8cf>(const Packet8cf& a) {
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

12
Eigen/src/Core/arch/AltiVec/Complex.h
View File

@ -98,6 +98,7 @@ struct packet_traits<std::complex<float> > : 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<Packet2cf>(const Packet2cf& a) {
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 @@ struct packet_traits<std::complex<double> > : default_packet_traits {
HasMin = 0,
HasMax = 0,
HasSqrt = 1,
HasLog = 1,
HasSetLinear = 0
};
};
@ -618,6 +625,11 @@ EIGEN_STRONG_INLINE Packet1cd psqrt<Packet1cd>(const Packet1cd& a) {
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

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

@ -964,6 +964,34 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pdiv_complex(const Pa
return Packet(pdiv(result_scaled.v, y_max));
}
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;
@ -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 <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>> {

4
Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h
View File

@ -113,6 +113,10 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pdiv_complex(const Pa
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 <> \

17
Eigen/src/Core/arch/NEON/Complex.h
View File

@ -62,6 +62,7 @@ struct packet_traits<std::complex<float> > : 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<Packet2cf>(const Packet2cf& a) {
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 @@ struct packet_traits<std::complex<double> > : 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<Packet1cd>(const Packet1cd& a) {
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

12
Eigen/src/Core/arch/SSE/Complex.h
View File

@ -42,6 +42,7 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
HasDiv = 1,
HasNegate = 1,
HasSqrt = 1,
HasLog = 1,
HasAbs = 0,
HasAbs2 = 0,
HasMin = 0,
@ -232,6 +233,7 @@ struct packet_traits<std::complex<double> > : 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<Packet2cf>(const Packet2cf& a) {
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

12
Eigen/src/Core/arch/ZVector/Complex.h
View File

@ -60,6 +60,7 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
HasSub = 1,
HasMul = 1,
HasDiv = 1,
HasLog = 1,
HasNegate = 1,
HasAbs = 0,
HasAbs2 = 0,
@ -83,6 +84,7 @@ struct packet_traits<std::complex<double> > : 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<Packet1cd>(const Packet1cd& a, const Packet1c
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 @@ EIGEN_STRONG_INLINE Packet2cf pdiv<Packet2cf>(const Packet2cf& a, const Packet2c
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]);

43
test/packetmath.cpp
View File

@ -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>() * 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<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>

36
test/packetmath_test_shared.h
View File

@ -124,6 +124,19 @@ bool areEqual(const Scalar* a, const Scalar* b, int size) {
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 @@ 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<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>