From 710e43de653ccc3d0937134b35f4fdf97417c343 Mon Sep 17 00:00:00 2001
From: cfborges <50250495+cfborges@users.noreply.github.com>
Date: Mon, 1 Jul 2019 15:09:06 -0700
Subject: [PATCH] Correctly rounded variant of the hypot code (#32345)

* Implements the correctly rounded variant of the hypot code only in the case where there is a native FMA

* Adds comments explaining the functioning of the two branches in the computation.
---
 base/math.jl | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index d1a11fbdbad63..4aa2d6acbca39 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -579,12 +579,20 @@ function hypot(x::T,y::T) where T<:AbstractFloat
         scale = one(scale)
     end
     h = sqrt(muladd(ax,ax,ay*ay))
-    if h <= 2*ay
-        delta = h-ay
-        h -= muladd(delta,delta-2*(ax-ay),ax*(2*delta - ax))/(2*h)
+    # This branch is correctly rounded but requires a native hardware fma.
+    if Base.Math.FMA_NATIVE
+        hsquared = h*h
+        axsquared = ax*ax
+        h -= (fma(-ay,ay,hsquared-axsquared) + fma(h,h,-hsquared) - fma(ax,ax,-axsquared))/(2*h)
+    # This branch is within one ulp of correctly rounded.
     else
-        delta = h-ax
-        h -= muladd(delta,delta,muladd(ay,(4*delta-ay),2*delta*(ax-2*ay)))/(2*h)
+        if h <= 2*ay
+            delta = h-ay
+            h -= muladd(delta,delta-2*(ax-ay),ax*(2*delta - ax))/(2*h)
+        else
+            delta = h-ax
+            h -= muladd(delta,delta,muladd(ay,(4*delta-ay),2*delta*(ax-2*ay)))/(2*h)
+        end
     end
     h*scale
 end