devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix another complex div edge case.

Browse Source 
libeigen/eigen!2257
This commit is contained in:
Antonio Sánchez
2026-03-06 21:37:26 +00:00
committed by Rasmus Munk Larsen
parent 5bacb5be9a
commit 20fce70e5a
2 changed files with 22 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
Eigen/src/Core/arch/Default/GenericPacketMathComplex.h
View File

@@ -29,11 +29,12 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pdiv_complex(const Pa
// are a pair length-2 SIMD vectors representing a single pair of complex // are a pair length-2 SIMD vectors representing a single pair of complex
// numbers x = a + i*b, y = c + i*d. // numbers x = a + i*b, y = c + i*d.
const RealPacket one = pset1<RealPacket>(RealScalar(1)); const RealPacket one = pset1<RealPacket>(RealScalar(1));
const RealPacket y_flip = pcplxflip(y).v; const RealPacket abs_y = pabs(y.v);
// We need to avoid dividing by Inf/Inf, so use a mask to carefully const RealPacket abs_y_flip = pcplxflip(Packet(abs_y)).v;
// apply the scale.
const RealPacket mask = pcmp_lt(pabs(y.v), pabs(y_flip)); // |c| < |d| const RealPacket mask = pcmp_lt(abs_y, abs_y_flip); // |c| < |d|
const RealPacket y_scaled = pselect(mask, pdiv(y.v, y_flip), one); RealPacket y_scaled = pselect(mask, pdiv(abs_y, abs_y_flip), one);
y_scaled = por(y_scaled, pandnot(y.v, abs_y)); // copy signs in case |c| == |d|
RealPacket denom = pmul(y.v, y_scaled); RealPacket denom = pmul(y.v, y_scaled);
denom = padd(denom, pcplxflip(Packet(denom)).v); // c * c' + d * d' denom = padd(denom, pcplxflip(Packet(denom)).v); // c * c' + d * d'
Packet num = pmul(x, pconj(Packet(y_scaled))); // a * c' + b * d', -a * d + b * c Packet num = pmul(x, pconj(Packet(y_scaled))); // a * c' + b * d', -a * d + b * c

16
test/packetmath.cpp
View File

@@ -1624,6 +1624,22 @@ void packetmath_complex() {
const RealScalar inf = std::numeric_limits<RealScalar>::infinity(); const RealScalar inf = std::numeric_limits<RealScalar>::infinity();
const RealScalar nan = std::numeric_limits<RealScalar>::quiet_NaN(); const RealScalar nan = std::numeric_limits<RealScalar>::quiet_NaN();
// Test division by a denominator with equal real and imaginary magnitudes
// to ensure pdiv scaling avoids division by zero (e.g. 1.0 - 1.0i).
if (PacketTraits::HasDiv) {
for (int i = 0; i < PacketSize; ++i) {
data1[i] = Scalar(one, zero);
RealScalar sign_re = (i & 1) ? -one : one;
RealScalar sign_im = (i & 2) ? -one : one;
data2[i] = Scalar(sign_re, sign_im);
}
internal::pstore(pval, internal::pdiv(internal::pload<Packet>(data1), internal::pload<Packet>(data2)));
for (int i = 0; i < PacketSize; ++i) {
Scalar expected = data1[i] / data2[i];
VERIFY_IS_APPROX(pval[i], expected);
}
}
// Multiplication and Division. // Multiplication and Division.
{ {
std::array<RealScalar, 8> special_values = {zero, one, inf, nan, -zero, -one, -inf, -nan}; std::array<RealScalar, 8> special_values = {zero, one, inf, nan, -zero, -one, -inf, -nan};