Skip to content

Commit

Permalink
simplify hypot
Browse files Browse the repository at this point in the history
  • Loading branch information
cossio committed Jan 14, 2020
1 parent d3addb5 commit b8a74b1
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 106 deletions.
114 changes: 9 additions & 105 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]+\\].*)*"
Expand All @@ -608,107 +602,22 @@ 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
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)

Expand Down Expand Up @@ -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)))

Expand Down
2 changes: 1 addition & 1 deletion test/testhelpers/Furlongs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit b8a74b1

Please sign in to comment.