From 35c009007e0be2723d14b4bc1e224718d4e589e0 Mon Sep 17 00:00:00 2001
From: Simon Byrne <simonbyrne@gmail.com>
Date: Thu, 5 Nov 2020 13:41:26 -0800
Subject: [PATCH] Fix overflow for hypot in Float16

Fixes #38311
---
 base/math.jl | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 981f152ed383d..17f7081cc10b4 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -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.