diff --git a/base/math.jl b/base/math.jl index a08ba8df5b1c3d..eb90defa22d8b9 100644 --- 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,77 +602,7 @@ 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 @@ -686,29 +610,14 @@ 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 index 02437d5a62628a..8bb6faab205d08 100644 --- 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 index b10036f556b98e..d1c180144a57b9 100644 --- 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)