diff --git a/base/math.jl b/base/math.jl index 0a9fac443c602..064c3e6058f37 100644 --- 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; @@ -538,6 +544,47 @@ Stacktrace: ``` """ 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 index 9887f7c7a1289..cc5ed20412984 100644 --- a/test/math.jl +++ b/test/math.jl @@ -191,6 +191,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)