From 20fce70e5a27c5924035ec592bad6e063b8aef17 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Antonio=20S=C3=A1nchez?= <cantonios@google.com>
Date: Fri, 6 Mar 2026 21:37:26 +0000
Subject: [PATCH] Fix another complex div edge case.

libeigen/eigen!2257
---
 .../Core/arch/Default/GenericPacketMathComplex.h | 11 ++++++-----
 test/packetmath.cpp                              | 16 ++++++++++++++++
 2 files changed, 22 insertions(+), 5 deletions(-)

diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h b/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h
index 6363d9ec1..8ceaf967b 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h
@@ -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
   // numbers x = a + i*b, y = c + i*d.
   const RealPacket one = pset1<RealPacket>(RealScalar(1));
-  const RealPacket y_flip = pcplxflip(y).v;
-  // We need to avoid dividing by  Inf/Inf, so use a mask to carefully
-  // apply the scale.
-  const RealPacket mask = pcmp_lt(pabs(y.v), pabs(y_flip));  // |c| < |d|
-  const RealPacket y_scaled = pselect(mask, pdiv(y.v, y_flip), one);
+  const RealPacket abs_y = pabs(y.v);
+  const RealPacket abs_y_flip = pcplxflip(Packet(abs_y)).v;
+
+  const RealPacket mask = pcmp_lt(abs_y, abs_y_flip);  // |c| < |d|
+  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);
   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
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index 37437aee0..33fb852a2 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -1624,6 +1624,22 @@ void packetmath_complex() {
   const RealScalar inf = std::numeric_limits<RealScalar>::infinity();
   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.
   {
     std::array<RealScalar, 8> special_values = {zero, one, inf, nan, -zero, -one, -inf, -nan};