`cancelminus` specification

[Chris00]

The following code

a = @interval(1e308, 1e308)
b = @interval(1e-308, 1e-308)
@show cancelminus(a, b)

prints

cancelminus(a, b) = [9.99999e+307, 1.00001e+308]

Is it really what the standard mandates? Shouldn't it be Interval(prevfloat(Inf), Inf)?

[lucaferranti]

according to the standard (section 10.5.6) c = cancelminus(a, b) is the tightest interval c so that a ⊆ b + c. While checking for tightness is always a beast of itself, that property is satisfied with the current behavior

julia> a = 1e308..1e308
[9.99999e+307, 1.00001e+308]

julia> b = 1e-308..1e-308
[9.99999e-309, 1.00001e-308]

julia> c = cancelminus(a, b)
[9.99999e+307, 1.00001e+308]

julia> a ⊆ b + c
true

While wouldn't with the output you suggested

 julia> cc = Interval(prevfloat(Inf), Inf)
[1.79769e+308, ∞]

julia> a ⊆ b + cc
false

The standard also gives a numerical receipt on how to compute it

This specification leads to the following Level 1 algorithm. If x = ∅ and y is bounded, then z = ∅. If x != ∅
and y = ∅, then z has no value. If x = [xlo, xhi] and y = [ylo, yhi] are both nonempty and bounded, define zlo = xlo−ylo
and zhi = xhi − yhi. Then z is defined to be [zlo, zhi] if zlo ≤ zhi (equivalently if width(x) ≥ width(y)), and has no value
otherwise. If either x or y is unbounded, z has no value

edit: lower and upper bars were lost when copy-pasting from the standard, clarified notation

[Chris00]

Sorry, I copied the wrong example. Here is it

a = @interval(1e308, 1e308)
b = @interval(-1e-308, -1e-308)
@show cancelminus(a, b)

which yields

cancelminus(a, b) = [9.99999e+307, 1.00001e+308]

[lucaferranti]

sorry, but I still fail to see your point

julia> a = @interval(1e308, 1e308)
[9.99999e+307, 1.00001e+308]

julia> b = @interval(-1e-308, -1e-308)
[-1.00001e-308, -9.99999e-309]

julia> c = cancelminus(a, b)
[9.99999e+307, 1.00001e+308]

julia> cc = Interval(prevfloat(Inf), Inf)
[1.79769e+308, ∞]

julia> a ⊆ b + c
true

julia> a ⊆ b + cc
false

The standard also gives a numerical receipt on how to compute the cancel minus, I just noticed that when I copy-pasted the description in my previous comment I lost the upper and lower bars on the characters, making it harder to read, I've now edited to clarify.

[lucaferranti]

A-haa I think you meant the following example

julia> a = @interval(1e308, 1e308)
[9.99999e+307, 1.00001e+308]

julia> b = @interval(-1e308, -1e308)
[-1.00001e+308, -9.99999e+307]

julia> c = cancelminus(a, b)
[-∞, ∞]

julia> cc = Interval(prevfloat(Inf), Inf)
[1.79769e+308, ∞]

julia> a ⊆ b + c
true

julia> a ⊆ b + cc
true

Good point, it might be a bug in cancel minus, because the behavior in IntervalArithmetic.jl differs e.g. from the octave package

 X>> pkg load interval
>> a = infsup(1e308, 1e308)
a = [1e+308, 1.0001e+308]
>> b = infsup(-1e308, -1e308)
b = [-1.0001e+308, -1e+308]
>> c = cancelminus(a, b)
c = [1.7976e+308, inf]

[lucaferranti]

A quick solution that seems to fix this without breaking other tests is changing this line to

isfinite(a_hi) && nextfloat(a.hi) < a_hi

just removing the above mentioned line would lead to this ITF test to fail, @dpsanders do you have comments?

edit, for better comfort, the following is the critical test in ITF1788, displayed is the current behavior in IA.jl (which matches what is expected by ITF1788)

julia> a = interval(-1.0, 2.2204460492503126e-16)
[-1, 2.22045e-16]

julia> b = interval(-1.0, 2.2204460492503128e-16)
[-1, 2.22045e-16]

julia> cancelplus(a, b)
[-∞, ∞]

[Chris00]

A-haa I think you meant the following example

Yes, thanks. (Doing too many things at the same time... 😀)

[lucaferranti]

@lbenet do you have comments on this?

[lbenet]

Sorry for the long silence... I'll take a look on this this afternoon and over the weekend.

[lbenet]

@lucaferranti The problem is indeed related to this line. The subtle cases are described in the standard (sect 12.12.5, page 62), and that line was supposed to solve one of them.

I just pushed #509 which should solve the problem spotted by @Chris00, as well as the case identified by you. They correspond to two corner cases described in the standard.

Using #509, I obtain the following:

# Chris broken case
julia> a = @interval(1e308, 1e308)
[9.99999e+307, 1.00001e+308]

julia> b = -a
[-1.00001e+308, -9.99999e+307]

julia> cancelminus(a,b)
[1.79769e+308, ∞]

# Luca's example, slightly changed
julia> a = interval(-1.0, 2.2204460492503126e-16)
[-1, 2.22045e-16]

julia> b = interval(-2.2204460492503128e-16, 1.0)
[-2.22045e-16, 1]

julia> cancelminus(a,b)
[-∞, ∞]

[lucaferranti]

a-haa thank you for spotting that!