From 91aa5feeadea5b0ce4ba4b07863ace36a36a42c7 Mon Sep 17 00:00:00 2001
From: Inkrementator <70717315+Inkrementator@users.noreply.github.com>
Date: Mon, 17 Mar 2025 21:19:44 +0100
Subject: [PATCH] Unify `hypot` and abs of complex values

This ports the fixes from abs(T)(ComplexT) regarding NaN and skipping
potentially inaccuracy inducing mathematical no-ops when one arg is 0 to
hypot. Both functions do the same thing and should be deduplicated.
hypot also has some logic regarding under and overflows, and while I
don't fully understand it, it should probably not be removed for complex
numbers.
---
 std/complex.d        | 20 ++------------------
 std/math/algebraic.d | 30 +++++++++++++++++++++++-------
 2 files changed, 25 insertions(+), 25 deletions(-)

diff --git a/std/complex.d b/std/complex.d
index f372b399e..256d0dc2b 100644
--- a/std/complex.d
+++ b/std/complex.d
@@ -767,24 +767,8 @@ if (is(T R == Complex!R))
  */
 T abs(T)(Complex!T z) @safe pure nothrow @nogc
 {
-    import std.math.algebraic : sqrt;
-    import std.math.traits : isInfinity, isNaN;
-    import std.math : fabs;
-    import std.algorithm.comparison : max;
-
-    if (isNaN(z.re) || isNaN(z.im))
-        return T.nan;
-    if (isInfinity(z.re) || isInfinity(z.im))
-        return T.infinity;
-
-    const T maxabs = max(fabs(z.re), fabs(z.im));
-    if (maxabs == 0.0)
-        return maxabs;
-
-    const T qx = z.re / maxabs;
-    const T qy = z.im / maxabs;
-
-    return maxabs * sqrt(qx * qx + qy * qy);
+    import std.math.algebraic : hypot;
+    return hypot(z.re, z.im);
 }
 
 ///
diff --git a/std/math/algebraic.d b/std/math/algebraic.d
index 5870ecdc8..929e90ea2 100644
--- a/std/math/algebraic.d
+++ b/std/math/algebraic.d
@@ -294,10 +294,11 @@ real cbrt(real x) @trusted pure nothrow @nogc
  * hypot(x, -y) are equivalent.
  *
  *  $(TABLE_SV
- *  $(TR $(TH x)            $(TH y)            $(TH hypot(x, y)) $(TH invalid?))
- *  $(TR $(TD x)            $(TD $(PLUSMN)0.0) $(TD |x|)         $(TD no))
- *  $(TR $(TD $(PLUSMNINF)) $(TD y)            $(TD +$(INFIN))   $(TD no))
- *  $(TR $(TD $(PLUSMNINF)) $(TD $(NAN))       $(TD +$(INFIN))   $(TD no))
+ *  $(TR $(TH x)            $(TH y)                 $(TH hypot(x, y)) $(TH invalid?))
+ *  $(TR $(TD x)            $(TD $(PLUSMN)0.0)      $(TD |x|)         $(TD no))
+ *  $(TR $(TD $(PLUSMNINF)) $(TD y)                 $(TD +$(INFIN))   $(TD no))
+ *  $(TR $(TD $(PLUSMNINF)) $(TD $(NAN))            $(TD +$(INFIN))   $(TD no))
+ *  $(TR $(TD $(NAN))       $(TD y != $(PLUSMNINF)) $(TD $(NAN))      $(TD no))
  *  )
  */
 T hypot(T)(const T x, const T y) @safe pure nothrow @nogc
@@ -308,13 +309,14 @@ if (isFloatingPoint!T)
     // If both are huge, avoid overflow by scaling by 2^^-N.
     // If both are tiny, avoid underflow by scaling by 2^^N.
     import core.math : fabs, sqrt;
-    import std.math.traits : floatTraits, RealFormat;
+    import std.math.traits : floatTraits, RealFormat, isNaN;
+    import std.algorithm.comparison : max;
 
     alias F = floatTraits!T;
 
     T u = fabs(x);
     T v = fabs(y);
-    if (!(u >= v))  // check for NaN as well.
+    if (!(u >= v))
     {
         v = u;
         u = fabs(y);
@@ -322,6 +324,15 @@ if (isFloatingPoint!T)
         if (v == T.infinity) return v; // hypot(nan, inf) == inf
     }
 
+    if (u.isNaN || v.isNaN)
+        return T.nan;
+
+    const maxabs = max(u,v);
+    if (v == 0.0)
+    {
+        return u;
+    }
+
     static if (F.realFormat == RealFormat.ieeeSingle)
     {
         enum SQRTMIN = 0x1p-60f;
@@ -375,18 +386,23 @@ if (isFloatingPoint!T)
     }
 
     // both are in the normal range
-    return ratio * sqrt(u*u + v*v);
+    return ratio * sqrt(u^^2 + v^^2);
 }
 
 ///
 @safe unittest
 {
     import std.math.operations : feqrel;
+    import std.math.traits : isNaN;
 
     assert(hypot(1.0, 1.0).feqrel(1.4142) > 16);
     assert(hypot(3.0, 4.0).feqrel(5.0) > 16);
     assert(hypot(real.infinity, 1.0L) == real.infinity);
+    assert(hypot(1.0L, real.infinity) == real.infinity);
     assert(hypot(real.infinity, real.nan) == real.infinity);
+    assert(hypot(real.nan, real.infinity) == real.infinity);
+    assert(hypot(real.nan, 1.0L).isNaN);
+    assert(hypot(1.0L, real.nan).isNaN);
 }
 
 @safe unittest