Skip to content

Commit

Permalink
fix hypot with multiple arguments
Browse files Browse the repository at this point in the history
  • Loading branch information
Jorge Fernandez-de-Cossio-Diaz committed Jan 12, 2020
1 parent aa35ee2 commit 7699f7d
Show file tree
Hide file tree
Showing 3 changed files with 109 additions and 24 deletions.
57 changes: 40 additions & 17 deletions base/math.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -610,37 +610,38 @@ julia> hypot(3, 4im)
5.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))
function hypot(x::T, y::T) where T<:AbstractFloat
hypot(x::Number) = abs(float(x))
hypot(x::Number, xs::Number...) = hypot(promote(x, xs...)...)
function hypot(x::T, y::T) where T <: Number
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(x) || isinf(y)
return T(Inf)
if isinf(ax) || isinf(ay)
return oftype(ax, Inf) * oneunit(x)
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
if ay <= ax*sqrt(eps(typeof(ax))/2) #Note: This also gets ay == 0
return ax * oneunit(x)
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 @@ -652,15 +653,16 @@ function hypot(x::T, y::T) where T<:AbstractFloat
else
if h <= 2*ay
delta = h-ay
h -= muladd(delta,delta-2*(ax-ay),ax*(2*delta - ax))/(2*h)
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)
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(x)
end


"""
hypot(x...)
Expand All @@ -675,7 +677,28 @@ julia> hypot(3, 4im, 12.0)
13.0
```
"""
hypot(x::Number...) = sqrt(sum(abs2(y) for y in x))
function hypot(x1::T, x2::T, x3::T, xs::T...) where {T<:Number}
v = (float(x1), float(x2), float(x3), float.(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
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)
Expand Down
67 changes: 61 additions & 6 deletions test/math.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,6 @@ 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
Expand Down Expand Up @@ -1032,8 +1031,64 @@ 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
end
9 changes: 8 additions & 1 deletion test/testhelpers/Furlongs.jl
View file
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 @@ -46,7 +53,7 @@ for f in (:real,:imag,:complex,:+,:-)
end

import Base: +, -, ==, !=, <, <=, isless, isequal, *, /, //, div, rem, mod, ^, hypot
for op in (:+, :-, :hypot)
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

0 comments on commit 7699f7d

Please sign in to comment.