simplify hypot

JuliaLang · Jan 14, 2020 · 768f33c · 768f33c
1 parent d3addb5
commit 768f33c
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 106 deletions.
diff --git a/base/math.jl b/base/math.jl
@@ -583,15 +583,9 @@ julia> sqrt(big(complex(-81)))
 sqrt(x::Real) = sqrt(float(x))
 
 """
-    hypot(x, y)
-
-Compute the hypotenuse ``\\sqrt{|x|^2+|y|^2}`` avoiding overflow and underflow.
+    hypot(x...)
 
-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
+Compute the hypotenuse ``\\sqrt{\\sum |x_i|^2}`` avoiding overflow and underflow.
 
 # Examples
 ```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
@@ -608,107 +602,22 @@ Stacktrace:
 
 julia> hypot(3, 4im)
 5.0
-```
-"""
-hypot(x::Number) = abs(float(x))
-hypot(x::Number, xs::Number...) = hypot(promote(x, xs...)...)
-function hypot(x::T, y::T) where T<:Number
-	# preserves unit
-	axu = abs(float(x))
-	ayu = abs(float(y))
-
-	# unitless
-	ax = axu / oneunit(axu)
-	ay = ayu / oneunit(ayu)
-
-
-    #ax = abs(float(x / oneunit(x)))
-    #ay = abs(float(y / oneunit(y)))
-
-    # Return Inf if either or both imputs is Inf (Compliance with IEEE754)
-    if isinf(ax) || isinf(ay)
-    	return oftype(axu, Inf)
-    end
-
-    # Order the operands
-    if ay > ax
-    	axu, ayu = ayu, axu
-        ax, ay = ay, ax
-    end
-
-    # Widely varying operands
-    if ay ≤ ax * sqrt(eps(typeof(ax)) / 2)  # Note: This also gets ay == 0
-        return axu
-    end
 
-    # Operands do not vary widely
-    scale = eps(sqrt(floatmin(ax)))  # rescaling constant
-    if ax > sqrt(floatmax(ax) / 2)
-        ax = ax * scale
-        ay = ay * scale
-        scale = inv(scale)
-    elseif ay < sqrt(floatmin(ax))
-        ax = ax / scale
-        ay = ay / scale
-    else
-        scale = oneunit(scale)
-    end
-
-    h = sqrt(muladd(ax, ax, ay * ay))
-
-    if Base.Math.FMA_NATIVE # This branch is correctly rounded but requires a native hardware fma.
-        hsquared = h * h
-        axsquared = ax * ax
-        h -= (fma(-ay, ay, hsquared - axsquared) + fma(h, h, -hsquared) - fma(ax, ax, -axsquared))/(2*h)
-    else # This branch is within one ulp of correctly rounded.
-        if h ≤ 2ay
-            delta = h - ay
-            h -= muladd(delta, delta - 2 * (ax - ay), ax * (2delta - ax)) / (2h)
-        else
-            delta = h-ax
-            h -= muladd(delta, delta, muladd(ay, (4delta - ay), 2delta * (ax - 2ay))) / (2h)
-        end
-    end
-    return h * scale * oneunit(axu)
-end
-
-"""
-    hypot(x...)
-
-Compute the hypotenuse ``\\sqrt{\\sum |x_i|^2}`` avoiding overflow and underflow.
-
-# Examples
-```jldoctest
 julia> hypot(-5.7)
 5.7
 
 julia> hypot(3, 4im, 12.0)
 13.0
 ```
 """
-function hypot(x1::T, x2::T, x3::T, xs::T...) where {T<:Number}
-    v = float.((x1, x2, x3, xs...))
-    # speculatively try fast naive approach
-    s = sum(abs2, v)
-    if isnan(s)
-        return any(isinf, v) ? oftype(sqrt(s), Inf) : sqrt(s) # IEEE 754
-    elseif isinf(s) || s ≤ oneunit(s) * _floatmin(one(s)) * (length(v) - 1)
-        # if overflow/underflow, try normalization
-        ma = maximum(abs, v)
-	    if iszero(ma) || isinf(ma)
-	        return ma
-	    else
-	        return ma * sqrt(sum(y -> abs2(y / ma), v))
-	    end
-    else
-        return sqrt(s)
-    end
+hypot(x::Number, xs::Number...) = _hypot(float.(promote(x, xs...))...)
+function _hypot(x::T...) where {T<:Number}
+    maxabs = maximum(abs, x)
+    isnan(maxabs) && any(isinf, x) && return T(Inf)
+    (iszero(maxabs) || isinf(maxabs)) && return maxabs
+    return maxabs * sqrt(sum(y -> abs2(y / maxabs), x))
 end
 
-_floatmin(x::AbstractFloat) = _floatmin(typeof(x))
-_floatmin(::Type{T}) where {T<:AbstractFloat} = nextfloat(zero(T)) / eps(T)
-_floatmin(::Type{T}) where {T<:IEEEFloat} = floatmin(T)
-
 atan(y::Real, x::Real) = atan(promote(float(y),float(x))...)
 atan(y::T, x::T) where {T<:AbstractFloat} = Base.no_op_err("atan", T)
 
@@ -1137,12 +1046,7 @@ for func in (:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh,
     end
 end
 
-for func in (:atan,:hypot)
-    @eval begin
-        $func(a::Float16,b::Float16) = Float16($func(Float32(a),Float32(b)))
-    end
-end
-
+atan(a::Float16,b::Float16) = Float16(atan(Float32(a),Float32(b)))
 cbrt(a::Float16) = Float16(cbrt(Float32(a)))
 sincos(a::Float16) = Float16.(sincos(Float32(a)))
 

diff --git a/test/math.jl b/test/math.jl
@@ -276,6 +276,7 @@ end
             @test hypot(T(0), T(0)) === T(0)
             @test hypot(T(Inf), T(Inf)) === T(Inf)
             @test hypot(T(Inf), T(x)) === T(Inf)
+            @test hypot(T(Inf), T(NaN)) === T(Inf)
             @test isnan_type(T, hypot(T(x), T(NaN)))
         end
     end

diff --git a/test/testhelpers/Furlongs.jl b/test/testhelpers/Furlongs.jl
@@ -52,7 +52,7 @@ for f in (:real,:imag,:complex,:+,:-)
     @eval Base.$f(x::Furlong{p}) where {p} = Furlong{p}($f(x.val))
 end
 
-import Base: +, -, ==, !=, <, <=, isless, isequal, *, /, //, div, rem, mod, ^, hypot
+import Base: +, -, ==, !=, <, <=, isless, isequal, *, /, //, div, rem, mod, ^
 for op in (:+, :-)
     @eval function $op(x::Furlong{p}, y::Furlong{p}) where {p}
         v = $op(x.val, y.val)