fix: Correctness fixes in interval_operations

src/intervals/interval_operations/bisect.jl

@@ -58,6 +58,7 @@ _check_minceable(x) = isbounded(x) ||
 In-place version of [`mince`](@ref).
 """
 function mince!(v::AbstractVector{<:BareInterval}, x::BareInterval{T}, n::Integer) where {T<:NumTypes}
+    resize!(v, n)
     if isempty_interval(x) # an empty interval has no nodes to interpolate
         for i ∈ 1:n
             v[i] = emptyinterval(x)
@@ -66,13 +67,18 @@ function mince!(v::AbstractVector{<:BareInterval}, x::BareInterval{T}, n::Intege
     end
     _check_minceable(x)
     nodes = LinRange(inf(x), sup(x), n+1)
+    lo = inf(x)
     @inbounds for i ∈ 1:n
-        v[i] = _unsafe_bareinterval(T, nodes[i], nodes[i+1])
+        # `LinRange` may interpolate non-monotonically, clamp the nodes to keep every piece valid
+        hi = i == n ? sup(x) : min(max(nodes[i+1], lo), sup(x))
+        v[i] = _unsafe_bareinterval(T, lo, hi)
+        lo = hi
     end
     return v
 end
 
 function mince!(v::AbstractVector{<:Interval}, x::Interval{T}, n::Integer) where {T<:NumTypes}
+    resize!(v, n)
     if isnai(x) | isempty_interval(x) # neither has nodes to interpolate
         y = isnai(x) ? nai(x) : emptyinterval(x)
         for i ∈ 1:n
@@ -84,10 +90,14 @@ function mince!(v::AbstractVector{<:Interval}, x::Interval{T}, n::Integer) where
     nodes = LinRange(inf(x), sup(x), n+1)
     d = decoration(x)
     t = isguaranteed(x)
+    lo = inf(x)
     @inbounds for i ∈ 1:n
-        rᵢ = _unsafe_bareinterval(T, nodes[i], nodes[i+1])
+        # `LinRange` may interpolate non-monotonically, clamp the nodes to keep every piece valid
+        hi = i == n ? sup(x) : min(max(nodes[i+1], lo), sup(x))
+        rᵢ = _unsafe_bareinterval(T, lo, hi)
         dᵢ = min(d, decoration(rᵢ))
         v[i] = _unsafe_interval(rᵢ, dᵢ, t)
+        lo = hi
     end
     return v
 end

src/intervals/interval_operations/boolean.jl

@@ -88,10 +88,9 @@ See also: [`issubset_interval`](@ref) and [`isinterior`](@ref).
 isstrictsubset(x::BareInterval, y::BareInterval) = issubset_interval(x, y) & !isequal_interval(x, y)
 
 isstrictsubset(x::Interval, y::Interval) = issubset_interval(x, y) & !isequal_interval(x, y)
-isstrictsubset(x::Complex{<:Interval}, y::Complex{<:Interval}) =
-    (isstrictsubset(real(x), real(y)) & issubset_interval(imag(x), imag(y))) | (issubset_interval(real(x), real(y)) & isstrictsubset(imag(x), imag(y)))
-isstrictsubset(x::Complex{<:Interval}, y::Interval) = isstrictsubset(real(x), y) & isthinzero(imag(x))
-isstrictsubset(x::Interval, y::Complex{<:Interval}) = isstrictsubset(x, real(y)) & in_interval(0, imag(y))
+isstrictsubset(x::Complex{<:Interval}, y::Complex{<:Interval}) = issubset_interval(x, y) & !isequal_interval(x, y)
+isstrictsubset(x::Complex{<:Interval}, y::Interval) = issubset_interval(x, y) & !isequal_interval(x, y)
+isstrictsubset(x::Interval, y::Complex{<:Interval}) = issubset_interval(x, y) & !isequal_interval(x, y)
 
 isstrictsubset(x::AbstractVector, y::AbstractVector) = issubset_interval(x, y) & any(t -> isstrictsubset(t[1], t[2]), zip(x, y))
 
@@ -145,8 +144,14 @@ function isdisjoint_interval(x::Interval, y::Interval)
     return isdisjoint_interval(bareinterval(x), bareinterval(y))
 end
 isdisjoint_interval(x::Complex{<:Interval}, y::Complex{<:Interval}) = isdisjoint_interval(real(x), real(y)) | isdisjoint_interval(imag(x), imag(y))
-isdisjoint_interval(x::Complex{<:Interval}, y::Interval) = isdisjoint_interval(real(x), y) | !in_interval(0, imag(x))
-isdisjoint_interval(x::Interval, y::Complex{<:Interval}) = isdisjoint_interval(x, real(y)) | !in_interval(0, imag(y))
+function isdisjoint_interval(x::Complex{<:Interval}, y::Interval)
+    isnai(x) | isnai(y) && return false
+    return isdisjoint_interval(real(x), y) | !in_interval(0, imag(x))
+end
+function isdisjoint_interval(x::Interval, y::Complex{<:Interval})
+    isnai(x) | isnai(y) && return false
+    return isdisjoint_interval(x, real(y)) | !in_interval(0, imag(y))
+end
 
 function isdisjoint_interval(x::AbstractVector, y::AbstractVector)
     n = length(x)
@@ -164,8 +169,7 @@ _isdisjoint_interval(x, y, z, w...) = _isdisjoint_interval(x, y) && _isdisjoint_
 """
     isweakless(x, y)
 
-Test whether `inf(x) ≤ inf(y)` and `sup(x) ≤ sup(y)`, where `<` is replaced by
-`≤` for infinite values.
+Test whether `inf(x) ≤ inf(y)` and `sup(x) ≤ sup(y)`.
 
 Implement the `less` function of the IEEE Standard 1788-2015
 (Table 10.3, and Sections 10.5.10 and 12.12.9).
@@ -197,7 +201,7 @@ end
 """
     precedes(x, y)
 
-Test whether any element of `x` is lesser or equal to every elements of `y`.
+Test whether every element of `x` is lesser or equal to every element of `y`.
 
 Implement the `precedes` function of the IEEE Standard 1788-2015
 (Table 10.3, and Sections 10.5.10 and 12.12.9).
@@ -212,7 +216,7 @@ end
 """
     strictprecedes(x, y)
 
-Test whether any element of `x` is strictly lesser than every elements of `y`.
+Test whether every element of `x` is strictly lesser than every element of `y`.
 
 Implement the `strictPrecedes` function of the IEEE Standard 1788-2015
 (Table 10.3, and Sections 10.5.10 and 12.12.9).
@@ -296,7 +300,7 @@ Implement the `isNaI` function of the IEEE Standard 1788-2015 (Section 12.12.9).
 isnai(::BareInterval) = false
 
 isnai(x::Interval) = decoration(x) == ill
-isnai(x::Complex{<:Interval}) = isnai(real(x)) & isnai(imag(x))
+isnai(x::Complex{<:Interval}) = isnai(real(x)) | isnai(imag(x))
 
 """
     isbounded(x)

src/intervals/interval_operations/cancellative.jl

@@ -7,11 +7,11 @@
 
 Compute the unique interval `z` such that `y + z == x`.
 
-The keywork `dec` argument controls the decoration of the result. It can be
+The keyword `dec` argument controls the decoration of the result. It can be
 either a specific decoration, or one of two following options:
     - `:default`: if at least one of the input intervals is `ill`,
         then the result is `ill`, otherwise it is `trv` (Section 11.7.1).
-    - `:auto`: the ouptut has the minimal decoration of the inputs.
+    - `:auto`: the output has the minimal decoration of the inputs.
 
 Implement the `cancelMinus` function of the IEEE Standard 1788-2015 (Section 9.2).
 """
@@ -32,32 +32,36 @@ function cancelminus(x::BareInterval{T}, y::BareInterval{T}) where {T<:NumTypes}
     z = _unsafe_bareinterval(T, z_lo, z_hi)
     isunbounded(z) && return z
 
-    # corner case 1 (page 62) involving finite precision for diam(x) and diam(y)
-    w_lo, w_hi = bounds(@round(T, inf(y) + z_lo, sup(y) + z_hi))
-    (diam(x) == diam(y)) & ((prevfloat(inf(x)) > w_lo) | (nextfloat(sup(x)) < w_hi)) && return entireinterval(BareInterval{T})
+    if T <: AbstractFloat # rational arithmetic is exact, so corner case 1 cannot occur
+        # corner case 1 (page 62) involving finite precision for diam(x) and diam(y)
+        w_lo, w_hi = bounds(@round(T, inf(y) + z_lo, sup(y) + z_hi))
+        (diam(x) == diam(y)) & ((prevfloat(inf(x)) > w_lo) | (nextfloat(sup(x)) < w_hi)) && return entireinterval(BareInterval{T})
+    end
 
     return z
 end
 cancelminus(x::BareInterval, y::BareInterval) = cancelminus(promote(x, y)...)
 
 function cancelminus(x::Interval, y::Interval; dec = :default)
-    r = cancelminus(bareinterval(x), bareinterval(y))
-    d = min(decoration(x), decoration(y), decoration(r))
+    d = min(decoration(x), decoration(y))
+    d == ill && return nai(promote_type(numtype(x), numtype(y))) # one of the inputs is an NaI
+    r = cancelminus(_bareinterval(x), _bareinterval(y))
+    d = min(d, decoration(r))
     t = isguaranteed(x) & isguaranteed(y)
     return _set_decoration(_unsafe_interval(r, d, t), dec)
 end
 
 """
     cancelplus(x, y; dec = :default)
 
-Compute the unique interval `z` such that `y - z == x`; this is semantically
+Compute the unique interval `z` such that `z - y == x`; this is semantically
 equivalent to `cancelminus(x, -y)`.
 
-The keywork `dec` argument controls the decoration of the result. It can be
+The keyword `dec` argument controls the decoration of the result. It can be
 either a specific decoration, or one of two following options:
     - `:default`: if at least one of the input intervals is `ill`,
         then the result is `ill`, otherwise it is `trv` (Section 11.7.1).
-    - `:auto`: the ouptut has the minimal decoration of the inputs.
+    - `:auto`: the output has the minimal decoration of the inputs.
 
 Implement the `cancelPlus` function of the IEEE Standard 1788-2015 (Section 9.2).
 """

src/intervals/interval_operations/constants.jl

@@ -5,7 +5,7 @@
     emptyinterval(T=[default_numtype()])
 
 Create an empty interval. This interval is an exception to the fact that the
-lower bound is larger than the upper one.
+lower bound is smaller than or equal to the upper one.
 
 Implement the `empty` function of the IEEE Standard 1788-2015 (Section 10.5.2).
 """

src/intervals/interval_operations/extended_div.jl

@@ -7,6 +7,8 @@
 Two-output division.
 
 Implement the `mulRevToPair` function of the IEEE Standard 1788-2015 (Section 10.5.5).
+`extended_div(x, y)` corresponds to `mulRevToPair(y, x)` of the standard (the
+numerator comes first).
 """
 function extended_div(x::BareInterval{T}, y::BareInterval{T}) where {T<:NumTypes}
     ylo, yhi = bounds(y)
@@ -21,6 +23,8 @@ function extended_div(x::BareInterval{T}, y::BareInterval{T}) where {T<:NumTypes
     end
 end
 
+extended_div(x::BareInterval, y::BareInterval) = extended_div(promote(x, y)...)
+
 function extended_div(x::Interval, y::Interval)
     bx = bareinterval(x)
     by = bareinterval(y)

src/intervals/interval_operations/numeric.jl

@@ -78,12 +78,25 @@ end
 
 bounds(x::Real) = bounds(interval(x))
 
+# analogues of `nextfloat(typemin(T))` and `prevfloat(typemax(T))` for
+# `Rational` bounds (cf. Section 12.12.8); unbounded integer types such as
+# `BigInt` have no smallest or largest finite value
+_finitemin(::Type{Rational{T}}) where {T<:Integer} = convert(Rational{T}, typemin(T))
+_finitemax(::Type{Rational{T}}) where {T<:Integer} = convert(Rational{T}, typemax(T))
+_finitemin(::Type{Rational{BigInt}}) =
+    throw(ArgumentError("cannot compute the midpoint of unbounded intervals with `Rational{BigInt}` bounds; there is no smallest finite `Rational{BigInt}`"))
+_finitemax(::Type{Rational{BigInt}}) =
+    throw(ArgumentError("cannot compute the midpoint of unbounded intervals with `Rational{BigInt}` bounds; there is no largest finite `Rational{BigInt}`"))
+
 """
     mid(x, α = 0.5)
 
 Relative midpoint of `x`, for `α` between 0 and 1 such that `mid(x, 0)` is the
 lower bound of the interval, `mid(x, 1)` its upper bound, and `mid(x, 0.5)` its
-midpoint.
+midpoint. For an unbounded interval, the finite bound is returned whenever `α`
+selects its side, and the infinite bound is replaced by the largest finite
+value of the bound type otherwise (cf. Section 12.12.8 of the IEEE Standard
+1788-2015).
 
 Implement the `mid` function of the IEEE Standard 1788-2015 (Table 9.2).
 
@@ -99,8 +112,14 @@ function mid(x::BareInterval{T}, α = 0.5) where {T<:AbstractFloat}
         return nextfloat(typemin(T))
     else
         lo, hi = bounds(x)
-        lo == typemin(T) && return nextfloat(lo) # cf. Section 12.12.8
-        hi == typemax(T) && return prevfloat(hi) # cf. Section 12.12.8
+        if lo == typemin(T) # cf. Section 12.12.8
+            α > 0.5 && return _normalisezero(hi)
+            return nextfloat(lo)
+        end
+        if hi == typemax(T) # cf. Section 12.12.8
+            α < 0.5 && return _normalisezero(lo)
+            return prevfloat(hi)
+        end
         β = convert(T, α)
         midpoint = β * (hi + lo * (1/β - 1)) # exactly 0.5 * (hi + lo) for β = 0.5
         midpoint = ifelse(isfinite(midpoint), midpoint, (1 - β) * lo + β * hi)
@@ -114,12 +133,18 @@ function mid(x::BareInterval{Rational{T}}, α = 1//2) where {T<:Integer}
     isempty_interval(x) && return throw(ArgumentError("cannot compute the midpoint of empty intervals; cannot return a `Rational` NaN"))
     if isentire_interval(x)
         α == 0.5 && return zero(Rational{T})
-        α > 0.5 && return convert(Rational{T}, typemax(T))
-        return convert(Rational{T}, typemin(T))
+        α > 0.5 && return _finitemax(Rational{T})
+        return _finitemin(Rational{T})
     else
         lo, hi = bounds(x)
-        lo == typemin(Rational{T}) && return convert(Rational{T}, typemin(T)) # cf. Section 12.12.8
-        hi == typemax(Rational{T}) && return convert(Rational{T}, typemax(T)) # cf. Section 12.12.8
+        if lo == typemin(Rational{T}) # cf. Section 12.12.8
+            α > 0.5 && return _normalisezero(hi)
+            return _finitemin(Rational{T})
+        end
+        if hi == typemax(Rational{T}) # cf. Section 12.12.8
+            α < 0.5 && return _normalisezero(lo)
+            return _finitemax(Rational{T})
+        end
         β = convert(Rational{T}, α)
         return _normalisezero((1 - β) * lo + β * hi)
     end
@@ -215,7 +240,7 @@ See also: [`inf`](@ref), [`sup`](@ref), [`bounds`](@ref), [`mid`](@ref),
 """
 function midradius(x::BareInterval)
     m = mid(x)
-    return m, max(m - inf(x), sup(x) - m)
+    return m, max(_fround(-, m, inf(x), RoundUp), _fround(-, sup(x), m, RoundUp)) # cf. Section 12.12.8
 end
 function midradius(x::BareInterval{<:Rational}) # needed to avoid integer overflow error
     m = mid(x)
@@ -263,7 +288,11 @@ function mag(x::Interval{<:Rational})
 end
 
 mag(x::Real) = mag(interval(x))
-mag(x::Complex) = sup(abs(interval(x)))
+function mag(x::Complex)
+    z = interval(x)
+    isempty_interval(z) && return mag(emptyinterval(real(z)))
+    return sup(abs(z))
+end
 
 """
     mig(x)
@@ -295,20 +324,26 @@ function mig(x::Interval{<:Rational})
 end
 
 mig(x::Real) = mig(interval(x))
-mig(x::Complex) = inf(abs(interval(x)))
+function mig(x::Complex)
+    z = interval(x)
+    isempty_interval(z) && return mig(emptyinterval(real(z)))
+    return inf(abs(z))
+end
 
 """
     dist(x, y)
 
-Distance between `x` and `y`.
+Upper bound of the Hausdorff distance between `x` and `y`: the bound
+differences are rounded upward, and equal bounds, including infinite ones, are
+at distance zero.
 """
 function dist(x::BareInterval{T}, y::BareInterval{T}) where {T<:AbstractFloat}
     isempty_interval(x) | isempty_interval(y) && return convert(T, NaN)
-    return max(abs(inf(x) - inf(y)), abs(sup(x) - sup(y)))
+    return max(_dist(inf(x), inf(y)), _dist(sup(x), sup(y)))
 end
 function dist(x::BareInterval{T}, y::BareInterval{T}) where {T<:Rational}
     isempty_interval(x) | isempty_interval(y) && return throw(ArgumentError("cannot compute the distance of empty intervals; cannot return a `Rational` NaN"))
-    return max(abs(inf(x) - inf(y)), abs(sup(x) - sup(y)))
+    return max(_dist(inf(x), inf(y)), _dist(sup(x), sup(y)))
 end
 dist(x::BareInterval, y::BareInterval) = dist(promote(x, y)...)
 
@@ -321,3 +356,11 @@ function dist(x::Interval{T}, y::Interval{T}) where {T<:Rational}
     return dist(bareinterval(x), bareinterval(y))
 end
 dist(x::Interval, y::Interval) = dist(promote(x, y)...)
+
+# used internally, upper bound of `abs(a - b)` with equal, possibly infinite,
+# bounds at distance zero
+function _dist(a::T, b::T) where {T<:NumTypes}
+    a == b && return zero(T)
+    isinf(a) & isinf(b) && return typemax(T) # opposite infinities
+    return _fround(-, max(a, b), min(a, b), RoundUp)
+end