Reduce the cost of hull and intersect_interval

src/intervals/construction.jl

@@ -44,6 +44,13 @@ Fields:
 - `lo :: T`
 - `hi :: T`
 
+The stored bounds satisfy two representation invariants: neither is `NaN`, the
+empty interval being stored as `(typemax(T), typemin(T))`; and a zero bound is
+stored as `+0`, with [`inf`](@ref) restoring the `-0` required by the standard.
+The inner constructor normalizes zero bounds, while keeping `NaN` out is left to
+the callers. [`sup`](@ref) and [`bounds`](@ref) rely on both invariants to read
+the fields directly.
+
 Constructor compliant with the IEEE Standard 1788-2015: [`bareinterval`](@ref).
 
 See also: [`Interval`](@ref).
@@ -295,6 +302,8 @@ _unsafe_interval
 # used only to construct intervals
 _inf(x::Interval) = x.bareinterval.lo
 _sup(x::Interval) = x.bareinterval.hi
+# Access the bare interval after ruling out NaI, without repeating the warning check.
+_bareinterval(x::Interval) = x.bareinterval
 #
 
 # avoid inlining the expanded code from @warn

src/intervals/interval_operations/constants.jl

@@ -9,11 +9,7 @@ lower bound is larger than the upper one.
 
 Implement the `empty` function of the IEEE Standard 1788-2015 (Section 10.5.2).
 """
-emptyinterval(::Type{BareInterval{T}}) where {T<:AbstractFloat} = _unsafe_bareinterval(T, convert(T, NaN), convert(T, NaN))
-# note: `using Base.unsafe_rational(Int, 0, 0)` as an equivalent to `NaN` for `Rational`
-# does not work well since most codes for `Rational` assume that the denominator cannot be zero
-# e.g. `iszero(Base.unsafe_rational(Int, 0, 0)) == true`
-emptyinterval(::Type{BareInterval{T}}) where {T<:Rational} = _unsafe_bareinterval(T, typemax(T), typemin(T))
+emptyinterval(::Type{BareInterval{T}}) where {T<:NumTypes} = _unsafe_bareinterval(T, typemax(T), typemin(T))
 emptyinterval(::BareInterval{T}) where {T<:NumTypes} = emptyinterval(BareInterval{T})
 
 emptyinterval(::Type{Interval{T}}) where {T<:NumTypes} = _unsafe_interval(emptyinterval(BareInterval{T}), trv, true)

src/intervals/interval_operations/numeric.jl

@@ -14,9 +14,8 @@ Implement the `inf` function of the IEEE Standard 1788-2015 (Table 9.2).
 See also: [`sup`](@ref), [`bounds`](@ref), [`mid`](@ref), [`diam`](@ref),
 [`radius`](@ref) and [`midradius`](@ref).
 """
-inf(x::BareInterval{T}) where {T<:AbstractFloat} = ifelse(isnan(x.lo), typemax(T), ifelse(iszero(x.lo), copysign(x.lo, -1), x.lo))
-function inf(x::BareInterval{T}) where {T<:BigFloat}
-    isnan(x.lo) && return typemax(T) # typemax(x.lo)
+inf(x::BareInterval{<:AbstractFloat}) = ifelse(iszero(x.lo), copysign(x.lo, -1), x.lo)
+function inf(x::BareInterval{BigFloat})
     iszero(x.lo) && return copysign(x.lo, -1)
     return x.lo
 end
@@ -44,13 +43,7 @@ Implement the `sup` function of the IEEE Standard 1788-2015 (Table 9.2).
 See also: [`inf`](@ref), [`bounds`](@ref), [`mid`](@ref), [`diam`](@ref),
 [`radius`](@ref) and [`midradius`](@ref).
 """
-sup(x::BareInterval{T}) where {T<:AbstractFloat} = ifelse(isnan(x.hi), typemin(T), x.hi)
-function sup(x::BareInterval{T}) where {T<:BigFloat}
-    isnan(x.hi) && return typemin(T) # typemin(x.hi)
-    return x.hi
-end
-
-sup(x::BareInterval{<:Rational}) = x.hi
+sup(x::BareInterval) = x.hi
 
 function sup(x::Interval{T}) where {T<:AbstractFloat}
     isnai(x) && return convert(T, NaN)
@@ -72,8 +65,7 @@ not normalize the infimum of the interval.
 See also: [`inf`](@ref), [`sup`](@ref), [`mid`](@ref), [`diam`](@ref),
 [`radius`](@ref) and [`midradius`](@ref).
 """
-bounds(x::BareInterval{T}) where {T<:AbstractFloat} = (ifelse(isnan(x.lo), typemax(T), x.lo), sup(x))
-bounds(x::BareInterval{<:Rational}) = (inf(x), sup(x))
+bounds(x::BareInterval) = (x.lo, x.hi)
 
 function bounds(x::Interval{T}) where {T<:AbstractFloat}
     isnai(x) && return (convert(T, NaN), convert(T, NaN))

src/intervals/interval_operations/set_operations.jl

@@ -13,39 +13,57 @@ end
 """
     intersect_interval(x, y; dec = :default)
 
-Returns the intersection of the intervals `x` and `y`, considered as (extended)
-sets of real numbers. That is, the set that contains the points common in `x`
-and `y`.
+Return the intersection of `x` and `y`, considered as extended sets of real
+numbers.
 
-The keywork `dec` argument controls the decoration of the result. It can be
-either a specific decoration, or one of two following options:
+The keyword argument `dec` controls the decoration of the result. It can be a
+specific decoration or one of the following two 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 `intersection` function of the IEEE Standard 1788-2015 (Section 9.3).
 """
 function intersect_interval(x::BareInterval{T}, y::BareInterval{T}) where {T<:NumTypes}
-    lo = max(inf(x), inf(y))
-    hi = min(sup(x), sup(y))
-    if lo > hi
-        return emptyinterval(BareInterval{T})
-    else
-        return _unsafe_bareinterval(T, lo, hi)
-    end
+    # An empty interval has bounds `(typemax(T), typemin(T))`, so an empty
+    # operand makes `lo > hi` without requiring a separate check.
+    lo = max(x.lo, y.lo)
+    hi = min(x.hi, y.hi)
+    lo > hi && return emptyinterval(BareInterval{T})
+    return _unsafe_bareinterval(T, lo, hi)
 end
 intersect_interval(x::BareInterval, y::BareInterval) = intersect_interval(promote(x, y)...)
 
 function intersect_interval(x::Interval{T}, y::Interval{S}; dec = :default) where {T<:NumTypes,S<:NumTypes}
-    isnai(x) | isnai(y) && return nai(promote_type(T, S))
-    r = intersect_interval(bareinterval(x), bareinterval(y))
     d = min(decoration(x), decoration(y))
+    d == ill && return nai(promote_type(T, S)) # one of the inputs is an NaI
+    r = intersect_interval(_bareinterval(x), _bareinterval(y))
     t = isguaranteed(x) & isguaranteed(y)
+    dec === :default && return _unsafe_interval(r, trv, t)
     return _set_decoration(_unsafe_interval(r, d, t), dec)
 end
 
 intersect_interval(x, y, z, w...; dec = :default) =
     reduce((a, b) -> intersect_interval(a, b; dec = dec), (x, y, z, w...))
+
+# Bare intervals carry no decoration, so they take no `dec` keyword.
+intersect_interval(x::BareInterval, y::BareInterval, z::BareInterval, w::Vararg{BareInterval,N}) where {N} =
+    reduce(intersect_interval, (x, y, z, w...))
+
+# Share the decoration checks across all arguments. `Vararg{Interval,N}` is
+# required over `Interval...`: without the `N` the method is not specialized on
+# the number of arguments past six, and the reductions below then allocate and
+# dispatch dynamically.
+function intersect_interval(x::Interval, y::Interval, z::Interval, w::Vararg{Interval,N}; dec = :default) where {N}
+    xs = (x, y, z, w...)
+    r = mapreduce(_bareinterval, intersect_interval, xs)
+    d = mapreduce(decoration, min, xs)
+    d == ill && return nai(numtype(r)) # one of the inputs is an NaI
+    t = mapreduce(isguaranteed, &, xs)
+    dec === :default && return _unsafe_interval(r, trv, t)
+    return _set_decoration(_unsafe_interval(r, d, t), dec)
+end
+
 intersect_interval(x::Complex, y::Complex; dec = :default) =
     complex(intersect_interval(real(x), real(y); dec = dec), intersect_interval(imag(x), imag(y); dec = dec))
 intersect_interval(x::Real, y::Complex; dec = :default) =
@@ -60,30 +78,51 @@ Return the interval hull of the intervals `x` and `y`, considered as (extended)
 sets of real numbers, i.e. the smallest interval that contains all of `x` and
 `y`.
 
-The keywork `dec` argument controls the decoration of the result. It can be
-either a specific decoration, or one of two following options:
+The keyword argument `dec` controls the decoration of the result. It can be a
+specific decoration or one of the following two 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 `convexHull` function of the IEEE Standard 1788-2015 (Section 9.3).
 """
 function hull(x::BareInterval{T}, y::BareInterval{T}) where {T<:NumTypes}
-    isempty_interval(x) & isempty_interval(y) && return x
-    return _unsafe_bareinterval(T, min(inf(x), inf(y)), max(sup(x), sup(y)))
+    # The bounds of an empty interval are neutral under `min` and `max`, so no
+    # separate emptiness check is required.
+    return _unsafe_bareinterval(T, min(x.lo, y.lo), max(x.hi, y.hi))
 end
 hull(x::BareInterval, y::BareInterval) = hull(promote(x, y)...)
 
 function hull(x::Interval{T}, y::Interval{S}; dec = :default) where {T<:NumTypes,S<:NumTypes}
-    isnai(x) | isnai(y) && return nai(promote_type(T, S))
-    r = hull(bareinterval(x), bareinterval(y))
     d = min(decoration(x), decoration(y))
+    d == ill && return nai(promote_type(T, S)) # one of the inputs is an NaI
+    r = hull(_bareinterval(x), _bareinterval(y))
     t = isguaranteed(x) & isguaranteed(y)
+    dec === :default && return _unsafe_interval(r, trv, t)
     return _set_decoration(_unsafe_interval(r, d, t), dec)
 end
 
 hull(x, y, z, w...; dec = :default) =
     reduce((a, b) -> hull(a, b; dec = dec), (x, y, z, w...))
+
+# Bare intervals carry no decoration, so they take no `dec` keyword.
+hull(x::BareInterval, y::BareInterval, z::BareInterval, w::Vararg{BareInterval,N}) where {N} =
+    reduce(hull, (x, y, z, w...))
+
+# Share the decoration checks across all arguments. `Vararg{Interval,N}` is
+# required over `Interval...`: without the `N` the method is not specialized on
+# the number of arguments past six, and the reductions below then allocate and
+# dispatch dynamically.
+function hull(x::Interval, y::Interval, z::Interval, w::Vararg{Interval,N}; dec = :default) where {N}
+    xs = (x, y, z, w...)
+    r = mapreduce(_bareinterval, hull, xs)
+    d = mapreduce(decoration, min, xs)
+    d == ill && return nai(numtype(r)) # one of the inputs is an NaI
+    t = mapreduce(isguaranteed, &, xs)
+    dec === :default && return _unsafe_interval(r, trv, t)
+    return _set_decoration(_unsafe_interval(r, d, t), dec)
+end
+
 hull(x::Complex, y::Complex; dec = :default) =
     complex(hull(real(x), real(y); dec = dec), hull(imag(x), imag(y); dec = dec))
 hull(x::Real, y::Complex; dec = :default) =

test/interval_tests/construction.jl

@@ -1,3 +1,30 @@
+@testset "Representation invariants of the bounds" begin
+    # `inf` and `sup` read the fields directly, so no operation may store a
+    # `NaN` bound; the empty interval is `(typemax(T), typemin(T))` instead
+    for T ∈ (Float16, Float32, Float64, BigFloat, Rational{Int}, Rational{BigInt})
+        e = emptyinterval(BareInterval{T})
+        @test (e.lo == typemax(T)) & (e.hi == typemin(T))
+        @test isempty_interval(e)
+        @test (inf(e) == typemax(T)) & (sup(e) == typemin(T))
+
+        f = bareinterval(T, 2, 1) # ill-formed, hence empty
+        @test (f.lo == typemax(T)) & (f.hi == typemin(T))
+
+        g = nai(Interval{T})
+        @test (g.bareinterval.lo == typemax(T)) & (g.bareinterval.hi == typemin(T))
+    end
+
+    # a zero bound is stored as `+0`, `inf` restoring the `-0` of the standard
+    x = bareinterval(0.0, 1.0)
+    @test (x.lo === 0.0) & (inf(x) === -0.0) & (sup(x) === 1.0)
+    y = bareinterval(-1.0, -0.0)
+    @test (y.hi === 0.0) & (sup(y) === 0.0)
+
+    # `bounds` reports what is stored, without normalizing the infimum
+    @test bounds(x) === (0.0, 1.0)
+    @test bounds(emptyinterval(BareInterval{Float64})) === (Inf, -Inf)
+end
+
 @testset "Difference between checked and unchecked bare intervals" begin
     @test IntervalArithmetic._unsafe_bareinterval(Float64, 1, 2) === bareinterval(1, 2)
 

test/interval_tests/set_operations.jl

@@ -1,5 +1,122 @@
 using IntervalArithmetic: interval_diff
 
+# the arguments must be locals of a function for `@allocated` to be meaningful,
+# and the call must be compiled before it is measured
+function alloc_hull8(x)
+    hull(x, x, x, x, x, x, x, x)
+    return @allocated hull(x, x, x, x, x, x, x, x)
+end
+
+function alloc_intersect8(x)
+    intersect_interval(x, x, x, x, x, x, x, x)
+    return @allocated intersect_interval(x, x, x, x, x, x, x, x)
+end
+
+@testset "hull and intersect_interval" begin
+    # the empty interval is neutral for `hull` and absorbing for
+    # `intersect_interval`, in either position and for both interval types
+    for T ∈ (Float64, Float32, BigFloat, Rational{Int})
+        x = bareinterval(T, 1, 2)
+        e = emptyinterval(BareInterval{T})
+
+        @test isequal_interval(hull(e, e), e)
+        @test isequal_interval(hull(x, e), x)
+        @test isequal_interval(hull(e, x), x)
+        @test isequal_interval(hull(x, bareinterval(T, 5, 6)), bareinterval(T, 1, 6))
+
+        @test isempty_interval(intersect_interval(e, e))
+        @test isempty_interval(intersect_interval(x, e))
+        @test isempty_interval(intersect_interval(e, x))
+        @test isempty_interval(intersect_interval(x, bareinterval(T, 5, 6)))
+        @test isequal_interval(intersect_interval(x, bareinterval(T, 0, 3)), x)
+
+        y = interval(T, 1, 2)
+        f = emptyinterval(Interval{T})
+
+        @test isempty_interval(hull(f, f))
+        @test isequal_interval(hull(y, f), y)
+        @test isequal_interval(hull(f, y), y)
+        @test isempty_interval(intersect_interval(y, f))
+        @test isempty_interval(intersect_interval(f, y))
+    end
+
+    # an NaI operand poisons the result, in either position and at any arity
+    n = nai(Float64)
+    x = interval(1, 2)
+
+    @test isnai(hull(n, x)) & isnai(hull(x, n)) & isnai(hull(n, n))
+    @test isnai(intersect_interval(n, x)) & isnai(intersect_interval(x, n))
+    @test isnai(hull(x, x, n)) & isnai(hull(x, n, x)) & isnai(hull(n, x, x))
+    @test isnai(hull(x, x, x, n))
+    @test isnai(intersect_interval(x, x, n)) & isnai(intersect_interval(n, x, x))
+
+    # `dec = :default` yields `trv`, `dec = :auto` the minimal input decoration
+    y = interval(3, 4)
+    z = interval(5, Inf)
+
+    @test decoration(hull(x, y)) == trv
+    @test decoration(hull(x, y; dec = :auto)) == com
+    @test decoration(hull(x, z; dec = :auto)) == dac
+    @test decoration(hull(x, y; dec = def)) == def
+    @test decoration(intersect_interval(x, y)) == trv
+    @test decoration(intersect_interval(x, interval(2, 3); dec = :auto)) == com
+    @test_throws ArgumentError hull(x, y; dec = :nonsense)
+    @test_throws ArgumentError intersect_interval(x, y; dec = :nonsense)
+
+    # `dec = :auto` cannot promise `com` for an unbounded hull
+    @test decoration(hull(interval(-Inf, 0), interval(1, 2); dec = :auto)) == dac
+
+    # bare intervals accept the variadic forms, without a `dec` keyword
+    bs = (bareinterval(1, 2), bareinterval(-3, 0), bareinterval(5, 6))
+
+    @test isequal_interval(hull(bs...), bareinterval(-3, 6))
+    @test isequal_interval(hull(bs..., bareinterval(7, 9)), bareinterval(-3, 9))
+    @test isempty_interval(intersect_interval(bs...))
+    @test isequal_interval(
+        intersect_interval(bareinterval(0, 4), bareinterval(1, 5), bareinterval(2, 6)),
+        bareinterval(2, 4))
+    @test_throws MethodError hull(bs...; dec = :auto)
+
+    # the variadic forms agree with the pairwise reduction they replace
+    args = (interval(1, 2), interval(-3, 0), interval(5, 6), interval(-1, 8),
+        interval(0, 1), interval(-2, 2), interval(4, 7), interval(-5, 5))
+    for dec ∈ (:default, :auto, trv, def, com)
+        for k ∈ 3:8
+            xs = args[1:k]
+            @test isequal_interval(hull(xs...; dec = dec), reduce((a, b) -> hull(a, b; dec = dec), xs))
+            @test decoration(hull(xs...; dec = dec)) == decoration(reduce((a, b) -> hull(a, b; dec = dec), xs))
+            @test isequal_interval(intersect_interval(xs...; dec = dec), reduce((a, b) -> intersect_interval(a, b; dec = dec), xs))
+            @test decoration(intersect_interval(xs...; dec = dec)) == decoration(reduce((a, b) -> intersect_interval(a, b; dec = dec), xs))
+        end
+    end
+
+    # The variadic methods must specialize on the number of arguments: without
+    # that, the reductions over the arguments allocate and dispatch dynamically
+    # past six arguments, which costs two orders of magnitude.
+    @test alloc_hull8(interval(1, 2)) == 0
+    @test alloc_intersect8(interval(1, 2)) == 0
+    @test alloc_hull8(bareinterval(1, 2)) == 0
+    @test alloc_intersect8(bareinterval(1, 2)) == 0
+
+    # `isguaranteed` is the conjunction over all the arguments
+    ng = convert(Interval{Float64}, 1)
+
+    @test !isguaranteed(hull(x, ng))
+    @test !isguaranteed(hull(x, y, ng))
+    @test !isguaranteed(hull(x, y, x, ng))
+    @test isguaranteed(hull(x, y, x))
+    @test !isguaranteed(intersect_interval(x, y, ng))
+
+    # mixed bound types promote, at any arity
+    @test numtype(hull(interval(Float32, 1, 2), interval(Float64, 3, 4))) === Float64
+    @test numtype(hull(interval(Float32, 1, 2), interval(Float32, 0, 1), interval(Float64, 3, 4))) === Float64
+    @test numtype(intersect_interval(interval(Float32, 1, 2), interval(Float32, 0, 3), interval(Float64, 1, 4))) === Float64
+
+    # a zero lower bound of the result still reports as `-0.0`
+    @test inf(hull(interval(0, 1), interval(2, 3))) === -0.0
+    @test inf(intersect_interval(interval(0, 1), interval(-1, 3))) === -0.0
+end
+
 @testset "removed interval" begin
     @test_throws ArgumentError intersect(interval(1))
     @test_throws ArgumentError intersect(interval(1), 2, [1], 4., 5)