improved implementation of hypot(a,b) (#31922)

* Replacement for hypot(a,b) Provides a fast and accurate implementation of hypot() that leverages the fused multiply add where available. The approach is explained and tested in detail in the paper: An Improved Algorithm for hypot(a,b) by Carlos F. Borges The article is available online at ArXiv at the link https://arxiv.org/abs/1904.09481 This version gives correctly rounded results more than 99% of the time. It is 10 times better than the C math library hypot function in this respect.
JuliaLang · Jun 13, 2019 · 4a04600 · 4a04600
1 parent b9ea972
commit 4a04600
Show file tree

Hide file tree

Showing 2 changed files with 50 additions and 0 deletions.
diff --git a/base/math.jl b/base/math.jl
@@ -523,6 +523,12 @@ sqrt(x::Real) = sqrt(float(x))
 
 Compute the hypotenuse ``\\sqrt{|x|^2+|y|^2}`` avoiding overflow and underflow.
 
+This code is an implementation of the algorithm described in:
+An Improved Algorithm for `hypot(a,b)`
+by Carlos F. Borges
+The article is available online at ArXiv at the link
+  https://arxiv.org/abs/1904.09481
+
 # Examples
 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
 julia> a = 10^10;
@@ -541,6 +547,47 @@ julia> hypot(3, 4im)
 ```
 """
 hypot(x::Number, y::Number) = hypot(promote(x, y)...)
+hypot(x::Complex, y::Complex) = hypot(promote(abs(x),abs(y))...)
+hypot(x::Integer, y::Integer) = hypot(promote(float(x), float(y))...)
+function hypot(x::T,y::T) where T<:AbstractFloat
+    #Return Inf if either or both imputs is Inf (Compliance with IEEE754)
+    if isinf(x) || isinf(y)
+        return convert(T,Inf)
+    end
+
+    # Order the operands
+    ax,ay = abs(x), abs(y)
+    if ay > ax
+        ax,ay = ay,ax
+    end
+
+    # Widely varying operands
+    if ay <= ax*sqrt(eps(T)/2)  #Note: This also gets ay == 0
+        return ax
+    end
+
+    # Operands do not vary widely
+    scale = eps(sqrt(floatmin(T)))  #Rescaling constant
+    if ax > sqrt(floatmax(T)/2)
+        ax = ax*scale
+        ay = ay*scale
+        scale = inv(scale)
+    elseif ay < sqrt(floatmin(T))
+        ax = ax/scale
+        ay = ay/scale
+    else
+        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)
+    else
+        delta = h-ax
+        h -= muladd(delta,delta,muladd(ay,(4*delta-ay),2*delta*(ax-2*ay)))/(2*h)
+    end
+    h*scale
+end
 function hypot(x::T, y::T) where T<:Number
     ax = abs(x)
     ay = abs(y)

diff --git a/test/math.jl b/test/math.jl
@@ -197,6 +197,9 @@ end
             @test isequal(expm1(T(0)), T(0))
             @test expm1(T(1)) ≈ T(ℯ)-1 atol=10*eps(T)
             @test isequal(hypot(T(3),T(4)), T(5))
+            @test isequal(hypot(floatmax(T),T(1)),floatmax(T))
+            @test isequal(hypot(floatmin(T)*sqrt(eps(T)),T(0)),floatmin(T)*sqrt(eps(T)))
+            @test isequal(floatmin(T)*hypot(1.368423059742933,1.3510496552495361),hypot(floatmin(T)*1.368423059742933,floatmin(T)*1.3510496552495361))
             @test isequal(log(T(1)), T(0))
             @test isequal(log(ℯ,T(1)), T(0))
             @test log(T(ℯ)) ≈ T(1) atol=eps(T)