Fix overflow when rounding (JuliaLang#42033)

LilithHafner · Mar 8, 2022 · 0854304 · 0854304
1 parent 33553f2
commit 0854304
Show file tree

Hide file tree

Showing 2 changed files with 34 additions and 2 deletions.
diff --git a/base/floatfuncs.jl b/base/floatfuncs.jl
@@ -166,6 +166,16 @@ function _round_invstep(x, invstep, r::RoundingMode)
     return y
 end
 
+# round x to multiples of 1/(invstepsqrt^2)
+# Using square root of step prevents overflowing
+function _round_invstepsqrt(x, invstepsqrt, r::RoundingMode)
+    y = round((x * invstepsqrt) * invstepsqrt, r) / invstepsqrt / invstepsqrt
+    if !isfinite(y)
+        return x
+    end
+    return y
+end
+
 # round x to multiples of step
 function _round_step(x, step, r::RoundingMode)
     # TODO: use div with rounding mode
@@ -186,10 +196,15 @@ function _round_digits(x, r::RoundingMode, digits::Integer, base)
     fx = float(x)
     if digits >= 0
         invstep = oftype(fx, base)^digits
-        _round_invstep(fx, invstep, r)
+        if isfinite(invstep)
+            return _round_invstep(fx, invstep, r)
+        else
+            invstepsqrt = oftype(fx, base)^oftype(fx, digits/2)
+            return _round_invstepsqrt(fx, invstepsqrt, r)
+        end
     else
         step = oftype(fx, base)^-digits
-        _round_step(fx, step, r)
+        return _round_step(fx, step, r)
     end
 end
 

diff --git a/test/floatfuncs.jl b/test/floatfuncs.jl
@@ -119,6 +119,23 @@ end
     @test round(Float32(1.2), sigdigits=5) === Float32(1.2)
     @test round(Float16(0.6), sigdigits=2) === Float16(0.6)
     @test round(Float16(1.1), sigdigits=70) === Float16(1.1)
+
+    # issue 37171
+    @test round(9.87654321e-308, sigdigits = 1) ≈ 1.0e-307
+    @test round(9.87654321e-308, sigdigits = 2) ≈ 9.9e-308
+    @test round(9.87654321e-308, sigdigits = 3) ≈ 9.88e-308
+    @test round(9.87654321e-308, sigdigits = 4) ≈ 9.877e-308
+    @test round(9.87654321e-308, sigdigits = 5) ≈ 9.8765e-308
+    @test round(9.87654321e-308, sigdigits = 6) ≈ 9.87654e-308
+    @test round(9.87654321e-308, sigdigits = 7) ≈ 9.876543e-308
+    @test round(9.87654321e-308, sigdigits = 8) ≈ 9.8765432e-308
+    @test round(9.87654321e-308, sigdigits = 9) ≈ 9.87654321e-308
+    @test round(9.87654321e-308, sigdigits = 10) ≈ 9.87654321e-308
+    @test round(9.87654321e-308, sigdigits = 11) ≈ 9.87654321e-308
+
+    @inferred round(Float16(1.), sigdigits=2)
+    @inferred round(Float32(1.), sigdigits=2)
+    @inferred round(Float64(1.), sigdigits=2)
 end
 
 @testset "literal pow matches runtime pow matches optimized pow" begin