Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix hypot with more than two arguments #30301

Closed
wants to merge 3 commits into from
Closed

fix hypot with more than two arguments #30301

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d3addb5
fix hypot with multiple arguments
Dec 7, 2018
e9f7b06
simplify hypot
Jan 14, 2020
7d05851
combine with correctly rounded 2 arg version
Jan 15, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 43 additions & 38 deletions base/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -593,6 +593,10 @@ by Carlos F. Borges
The article is available online at ArXiv at the link
https://arxiv.org/abs/1904.09481

hypot(x...)

Compute the hypotenuse ``\\sqrt{\\sum |x_i|^2}`` avoiding overflow and underflow.

# Examples
```jldoctest; filter = r"Stacktrace:(\\n \\[[0-9]+\\].*)*"
julia> a = Int64(10)^10;
Expand All @@ -608,40 +612,52 @@ Stacktrace:

julia> hypot(3, 4im)
5.0

julia> hypot(-5.7)
5.7

julia> hypot(3, 4im, 12.0)
13.0
```
"""
hypot(x::Number, y::Number) = hypot(promote(x, y)...)
hypot(x::Complex, y::Complex) = hypot(abs(x), abs(y))
hypot(x::T, y::T) where {T<:Real} = hypot(float(x), float(y))
hypot(x::T, y::T) where {T<:Number} = (z = y/x; abs(x) * sqrt(one(z) + z*z))
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 T(Inf)
hypot(x::Number) = abs(float(x))
hypot(x::Number,y::Number,xs::Number...) = _hypot(float.(promote(x,y,xs...))...)
function _hypot(x::T, y::T) where {T<:Number}
# preserves unit
axu = abs(x)
ayu = abs(y)

# unitless
ax = axu / oneunit(axu)
ay = ayu / oneunit(ayu)

# 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
ax,ay = abs(x), abs(y)
if ay > ax
axu,ayu = ayu,axu
ax,ay = ay,ax
end

# Widely varying operands
if ay <= ax*sqrt(eps(T)/2) #Note: This also gets ay == 0
return ax
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(T))) #Rescaling constant
if ax > sqrt(floatmax(T)/2)
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(T))
elseif ay < sqrt(floatmin(ax))
ax = ax/scale
ay = ay/scale
else
scale = one(scale)
scale = oneunit(scale)
end
h = sqrt(muladd(ax,ax,ay*ay))
# This branch is correctly rounded but requires a native hardware fma.
Expand All @@ -659,24 +675,18 @@ function hypot(x::T, y::T) where T<:AbstractFloat
h -= muladd(delta,delta,muladd(ay,(4*delta-ay),2*delta*(ax-2*ay)))/(2*h)
end
end
return h*scale
return h*scale*oneunit(axu)
end
function _hypot(x::T...) where {T<:Number}
maxabs = maximum(abs, x)
if isnan(maxabs) && any(isinf, x)
return oftype(maxabs, Inf)
elseif (iszero(maxabs) || isinf(maxabs))
return maxabs
else
return maxabs * sqrt(sum(y -> abs2(y / maxabs), x))
end
end
cossio marked this conversation as resolved.
Show resolved Hide resolved

"""
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
```
"""
hypot(x::Number...) = sqrt(sum(abs2(y) for y in x))

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)
Expand Down Expand Up @@ -1106,12 +1116,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)))

Expand Down
69 changes: 64 additions & 5 deletions test/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1032,8 +1032,67 @@ end

isdefined(Main, :Furlongs) || @eval Main include("testhelpers/Furlongs.jl")
using .Main.Furlongs
@test hypot(Furlong(0), Furlong(0)) == Furlong(0.0)
@test hypot(Furlong(3), Furlong(4)) == Furlong(5.0)
@test hypot(Furlong(NaN), Furlong(Inf)) == Furlong(Inf)
@test hypot(Furlong(Inf), Furlong(NaN)) == Furlong(Inf)
@test hypot(Furlong(Inf), Furlong(Inf)) == Furlong(Inf)
@test (@inferred hypot(Furlong(0), Furlong(0))) == Furlong(0.0)
@test (@inferred hypot(Furlong(3), Furlong(4))) == Furlong(5.0)
@test (@inferred hypot(Furlong(NaN), Furlong(Inf))) == Furlong(Inf)
@test (@inferred hypot(Furlong(Inf), Furlong(NaN))) == Furlong(Inf)
@test (@inferred hypot(Furlong(0), Furlong(0), Furlong(0))) == Furlong(0.0)
@test (@inferred hypot(Furlong(Inf), Furlong(Inf))) == Furlong(Inf)
@test (@inferred hypot(Furlong(1), Furlong(1), Furlong(1))) == Furlong(sqrt(3))
@test (@inferred hypot(Furlong(Inf), Furlong(NaN), Furlong(0))) == Furlong(Inf)
@test (@inferred hypot(Furlong(Inf), Furlong(Inf), Furlong(Inf))) == Furlong(Inf)
@test isnan(hypot(Furlong(NaN), Furlong(0), Furlong(1)))

@testset "hypot" begin
@test_throws MethodError hypot()
@test (@inferred hypot(floatmax())) == floatmax()
@test (@inferred hypot(floatmax(), floatmax())) == Inf
@test (@inferred hypot(floatmin(), floatmin())) == √2floatmin()
@test (@inferred hypot(floatmin(), floatmin(), floatmin())) == √3floatmin()
@test (@inferred hypot(1e-162)) ≈ 1e-162
@test (@inferred hypot(2e-162, 1e-162, 1e-162)) ≈ hypot(2,1,1)*1e-162
@test (@inferred hypot(1e162)) ≈ 1e162
@test hypot(-2) === 2.0
@test hypot(-2, 0) === 2.0
let i = typemax(Int)
@test (@inferred hypot(i, i)) ≈ i * √2
@test (@inferred hypot(i, i, i)) ≈ i * √3
@test (@inferred hypot(i, i, i, i)) ≈ 2.0i
@test (@inferred hypot(i//1, 1//i, 1//i)) ≈ i
end
let i = typemin(Int)
@test (@inferred hypot(i, i)) ≈ -√2i
@test (@inferred hypot(i, i, i)) ≈ -√3i
@test (@inferred hypot(i, i, i, i)) ≈ -2.0i
end
@testset "$T" for T in (Float32, Float64)
@test (@inferred hypot(T(Inf), T(NaN))) == T(Inf) # IEEE754 says so
@test (@inferred hypot(T(Inf), T(3//2), T(NaN))) == T(Inf)
@test (@inferred hypot(T(1e10), T(1e10), T(1e10), T(1e10))) ≈ 2e10
@test isnan_type(T, hypot(T(3), T(3//4), T(NaN)))
@test hypot(T(1), T(0)) === T(1)
@test hypot(T(1), T(0), T(0)) === T(1)
@test (@inferred hypot(T(Inf), T(Inf), T(Inf))) == T(Inf)
for s in (zero(T), floatmin(T)*1e3, floatmax(T)*1e-3, T(Inf))
@test hypot(1s, 2s) ≈ s * hypot(1,2) rtol=8eps(T)
@test hypot(1s, 2s, 3s) ≈ s * hypot(1,2,3) rtol=8eps(T)
end
end
@testset "$T" for T in (Float16, Float32, Float64, BigFloat)
let x = 1.1sqrt(floatmin(T))
@test (@inferred hypot(x,x/4)) ≈ x * sqrt(17/BigFloat(16))
@test (@inferred hypot(x,x/4,x/4)) ≈ x * sqrt(9/BigFloat(8))
end
let x = 2sqrt(nextfloat(zero(T)))
@test (@inferred hypot(x,x/4)) ≈ x * sqrt(17/BigFloat(16))
@test (@inferred hypot(x,x/4,x/4)) ≈ x * sqrt(9/BigFloat(8))
end
let x = sqrt(nextfloat(zero(T))/eps(T))/8, f = sqrt(4eps(T))
@test hypot(x,x*f) ≈ x * hypot(one(f),f) rtol=eps(T)
@test hypot(x,x*f,x*f) ≈ x * hypot(one(f),f,f) rtol=eps(T)
end
end
cossio marked this conversation as resolved.
Show resolved Hide resolved
# hypot on Complex returns Real
@test (@inferred hypot(3, 4im)) === 5.0
@test (@inferred hypot(3, 4im, 12)) == 13.0
end
11 changes: 9 additions & 2 deletions test/testhelpers/Furlongs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,13 @@ Base.oneunit(x::Type{Furlong{p,T}}) where {p,T} = Furlong{p,T}(one(T))
Base.zero(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(zero(T))
Base.zero(::Type{Furlong{p,T}}) where {p,T} = Furlong{p,T}(zero(T))
Base.iszero(x::Furlong) = iszero(x.val)
Base.float(x::Furlong{p}) where {p} = Furlong{p}(float(x.val))
Base.eps(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(eps(T))
Base.eps(::Furlong{p,T}) where {p,T<:AbstractFloat} = eps(Furlong{p,T})
Base.floatmin(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floatmin(T))
Base.floatmin(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmin(Furlong{p,T})
Base.floatmax(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floatmax(T))
Base.floatmax(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmax(Furlong{p,T})

# convert Furlong exponent p to a canonical form. This
# is not type stable, but it doesn't matter since it is used
Expand All @@ -45,8 +52,8 @@ 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
for op in (:+, :-, :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)
Furlong{p}(v)
Expand Down