Skip to content

Commit

Permalink
Fix overflow for hypot in Float16
Browse files Browse the repository at this point in the history
Fixes #38311
  • Loading branch information
simonbyrne committed Dec 1, 2020
1 parent 59d676c commit 35c0090
Showing 1 changed file with 10 additions and 7 deletions.
17 changes: 10 additions & 7 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -673,16 +673,19 @@ function _hypot(x, y)
end

# Operands do not vary widely
scale = eps(typeof(ax))*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))
L = sqrt(floatmin(ax))/2 # a power of 2, so inv(L) is exact
if ax > inv(L)
ax = ax*L
ay = ay*L
scale = inv(L)
elseif ay < L
scale = L*eps(typeof(ax))
# inv(scale) would underflow in Float16
# is a power of 2, so should be optimized to a multplication by the compiler for other types
ax = ax/scale
ay = ay/scale
else
scale = oneunit(scale)
scale = oneunit(L)
end
h = sqrt(muladd(ax, ax, ay*ay))
# This branch is correctly rounded but requires a native hardware fma.
Expand Down

0 comments on commit 35c0090

Please sign in to comment.