From 7699f7d5adefb30e7c60f9f1605d44f0bc8aa8b3 Mon Sep 17 00:00:00 2001 From: Jorge Fernandez-de-Cossio-Diaz Date: Fri, 7 Dec 2018 12:50:27 -0500 Subject: [PATCH] fix hypot with multiple arguments --- base/math.jl | 57 +++++++++++++++++++++--------- test/math.jl | 67 ++++++++++++++++++++++++++++++++---- test/testhelpers/Furlongs.jl | 9 ++++- 3 files changed, 109 insertions(+), 24 deletions(-) diff --git a/base/math.jl b/base/math.jl index 20cb374c5c33bd..15e3811662dfff 100644 --- a/base/math.jl +++ b/base/math.jl @@ -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. @@ -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...) @@ -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) diff --git a/test/math.jl b/test/math.jl index 4515f8ca7b8625..3214e42cdadf35 100644 --- a/test/math.jl +++ b/test/math.jl @@ -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 @@ -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 diff --git a/test/testhelpers/Furlongs.jl b/test/testhelpers/Furlongs.jl index 3c2f7a035e91cf..b10036f556b98e 100644 --- a/test/testhelpers/Furlongs.jl +++ b/test/testhelpers/Furlongs.jl @@ -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 @@ -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)