Core: Fix numext::sign for complex values at the ends of the range

libeigen/eigen!2989

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/src/Core/MathFunctions.h

diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 26fb2b7..e8499d1 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -838,8 +838,10 @@
     EIGEN_USING_STD(abs);
     real_type aa = abs(a);
     if (aa == real_type(0)) return Scalar(0);
-    aa = real_type(1) / aa;
-    return Scalar(numext::real(a) * aa, numext::imag(a) * aa);
+    // Divide rather than multiply by 1/aa: the reciprocal overflows for subnormal aa and is itself
+    // subnormal, hence inexact, for aa near the top of the range. psign_impl for complex packets
+    // divides for the same reason.
+    return Scalar(numext::real(a) / aa, numext::imag(a) / aa);
   }
 };
 

test/numext.cpp

diff --git a/test/numext.cpp b/test/numext.cpp
index d1578b3..94b4f03 100644
--- a/test/numext.cpp
+++ b/test/numext.cpp
@@ -234,6 +234,33 @@
   }
 }
 
+// numext::sign(z) = z / |z| for complex z. The interesting inputs are the ends of the range, where
+// forming 1/|z| first would overflow (subnormal |z|) or itself be subnormal, hence inexact (|z| near max).
+template <typename T>
+void check_complex_sign() {
+  typedef typename NumTraits<T>::Real Real;
+  const Real zero(0), one(1);
+
+  VERIFY_IS_EQUAL(numext::sign(T(zero, zero)), T(zero, zero));
+  VERIFY_IS_EQUAL(numext::sign(T(one, zero)), T(one, zero));
+  VERIFY_IS_EQUAL(numext::sign(T(zero, -one)), T(zero, -one));
+
+  for (Real r : {std::numeric_limits<Real>::denorm_min(), (std::numeric_limits<Real>::min)(),
+                 (std::numeric_limits<Real>::max)()}) {
+    VERIFY_IS_EQUAL(numext::sign(T(r, zero)), T(one, zero));
+    VERIFY_IS_EQUAL(numext::sign(T(-r, zero)), T(-one, zero));
+    VERIFY_IS_EQUAL(numext::sign(T(zero, r)), T(zero, one));
+  }
+
+  // Off the axes the magnitude is only as accurate as abs() itself, which cannot resolve |z| for a z
+  // whose components are at the bottom of the subnormal range; from the smallest normal upwards it can.
+  for (Real r : {(std::numeric_limits<Real>::min)(), one, (std::numeric_limits<Real>::max)() / Real(2)}) {
+    const T s = numext::sign(T(r, r));
+    VERIFY_IS_APPROX(numext::abs(s), one);
+    VERIFY_IS_EQUAL(numext::real(s), numext::imag(s));
+  }
+}
+
 template <typename T>
 void check_arg() {
   typedef typename NumTraits<T>::Real Real;
@@ -570,6 +597,9 @@
     CALL_SUBTEST(check_abs<std::complex<float>>());
     CALL_SUBTEST(check_abs<std::complex<double>>());
 
+    CALL_SUBTEST(check_complex_sign<std::complex<float>>());
+    CALL_SUBTEST(check_complex_sign<std::complex<double>>());
+
     CALL_SUBTEST(check_arg<std::complex<float>>());
     CALL_SUBTEST(check_arg<std::complex<double>>());