Fix cancelminus (#509)

* Fix cancelminus

* Bump patch version

* Add few tests related with one corner case mentioned in the standard

Project.toml

@@ -1,7 +1,7 @@
 name = "IntervalArithmetic"
 uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 repo = "https://github.com/JuliaIntervals/IntervalArithmetic.jl.git"
-version = "0.20.3"
+version = "0.20.4"
 
 [deps]
 CRlibm = "96374032-68de-5a5b-8d9e-752f78720389"

src/intervals/arithmetic.jl

@@ -531,21 +531,21 @@ function cancelminus(a::Interval{T}, b::Interval{T}) where T<:Real
     c_lo = @round_down(a.lo - b.lo)
     c_hi = @round_up(a.hi - b.hi)
 
+    # Interval(c_lo, c_hi) is not a proper interval
     c_lo > c_hi && return entireinterval(T)
+    c = Interval(c_lo, c_hi)
 
-    c_lo == Inf && return Interval(prevfloat(c_lo), c_hi)
-    c_hi == -Inf && return Interval(c_lo, nextfloat(c_hi))
+    # Corner case 2 (page 62), involving unbounded c
+    (c_lo == Inf || c_hi == -Inf) && return widen(c)
+    isunbounded(c) && return c
 
+    # Corner case 1 (page 62) involving finite precision for diam(a) and diam(b)
     a_lo = @round_down(b.lo + c_lo)
     a_hi = @round_up(b.hi + c_hi)
+    (diam(a) == diam(b)) && (nextfloat(a.hi) < a_hi || prevfloat(a.lo) > a_lo) &&
+        return entireinterval(T)
 
-    if a_lo ≤ a.lo ≤ a.hi ≤ a_hi
-        (nextfloat(a.hi) < a_hi || prevfloat(a.lo) > a_hi) &&
-            return entireinterval(T)
-        return Interval(c_lo, c_hi)
-     end
-
-    return entireinterval(T)
+    return c
 end
 cancelminus(a::Interval, b::Interval) = cancelminus(promote(a, b)...)
 

test/interval_tests/consistency.jl

@@ -262,6 +262,16 @@ setprecision(Interval, Float64)
         @test cancelplus(Interval(0.0), Interval(1.0)) == Interval(1.0)
         @test cancelminus(Interval(-5.0, 0.0), Interval(0.0, 5.0)) == Interval(-5.0)
         @test cancelplus(Interval(-5.0, 0.0), Interval(0.0, 5.0)) == Interval(0.0)
+        @test cancelminus(Interval(1e308), -Interval(1e308)) == widen(Interval(Inf))
+        @test cancelplus(Interval(1e308), Interval(1e308)) == widen(Interval(Inf))
+        @test cancelminus(Interval(nextfloat(1e308)), -Interval(nextfloat(1e308))) ==
+                widen(Interval(Inf))
+        @test cancelplus(Interval(nextfloat(1e308)), Interval(nextfloat(1e308))) ==
+                widen(Interval(Inf))
+        @test cancelminus(Interval(prevfloat(big(Inf))), -Interval(prevfloat(big(Inf)))) ==
+                widen(Interval(big(Inf)))
+        @test cancelplus(Interval(prevfloat(big(Inf))), Interval(prevfloat(big(Inf)))) ==
+                widen(Interval(big(Inf)))
     end
 
     @testset "mid and radius" begin