diff --git a/src/intervals/interval_operations/bisect.jl b/src/intervals/interval_operations/bisect.jl
index 7b73e734..dd80d2a4 100644
--- a/src/intervals/interval_operations/bisect.jl
+++ b/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
diff --git a/src/intervals/interval_operations/boolean.jl b/src/intervals/interval_operations/boolean.jl
index 16dd3b8b..52585b39 100644
--- a/src/intervals/interval_operations/boolean.jl
+++ b/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)
diff --git a/src/intervals/interval_operations/cancellative.jl b/src/intervals/interval_operations/cancellative.jl
index 1181b31f..ea1590c1 100644
--- a/src/intervals/interval_operations/cancellative.jl
+++ b/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,17 +32,21 @@ 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
@@ -50,14 +54,14 @@ 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).
 """
diff --git a/src/intervals/interval_operations/constants.jl b/src/intervals/interval_operations/constants.jl
index 51fa1f39..88429219 100644
--- a/src/intervals/interval_operations/constants.jl
+++ b/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).
 """
diff --git a/src/intervals/interval_operations/extended_div.jl b/src/intervals/interval_operations/extended_div.jl
index ae3928d3..fdb6df5b 100644
--- a/src/intervals/interval_operations/extended_div.jl
+++ b/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)
diff --git a/src/intervals/interval_operations/numeric.jl b/src/intervals/interval_operations/numeric.jl
index dd6439fd..7dfe16f3 100644
--- a/src/intervals/interval_operations/numeric.jl
+++ b/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
diff --git a/src/intervals/interval_operations/set_operations.jl b/src/intervals/interval_operations/set_operations.jl
index 89dd5790..ec3503dd 100644
--- a/src/intervals/interval_operations/set_operations.jl
+++ b/src/intervals/interval_operations/set_operations.jl
@@ -5,7 +5,7 @@
 _set_decoration(x::Interval, d::Decoration) = setdecoration(x, d)
 
 function _set_decoration(x::Interval, d::Symbol)
-    d === :auto    && return x
+    d === :auto    && return setdecoration(x, decoration(x))
     d === :default && return setdecoration(x, min(trv, decoration(x)))
     return throw(ArgumentError("decoration option must be `:default`, `:auto` or a specific decoration"))
 end
@@ -39,7 +39,6 @@ function intersect_interval(x::Interval{T}, y::Interval{S}; dec = :default) wher
     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
 
@@ -60,7 +59,6 @@ function intersect_interval(x::Interval, y::Interval, z::Interval, w::Vararg{Int
     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
 
@@ -98,7 +96,6 @@ function hull(x::Interval{T}, y::Interval{S}; dec = :default) where {T<:NumTypes
     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
 
@@ -119,7 +116,6 @@ function hull(x::Interval, y::Interval, z::Interval, w::Vararg{Interval,N}; dec
     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
 
@@ -143,11 +139,11 @@ const union_interval = hull
 Remove the interior of `y` from `x`. If `x` and `y` are vectors, then they are
 treated as multi-dimensional intervals.
 
-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.
 """
 interiordiff(x::BareInterval, y::BareInterval) =
     interiordiff!(Vector{promote_type(typeof(x), typeof(y))}(undef, 0), x, y)
@@ -162,30 +158,15 @@ interiordiff(x::AbstractVector{<:Interval}, y::AbstractVector{<:Interval}; dec =
     interiordiff!(Vector{promote_type(typeof(x), typeof(y))}(undef, 0), x, y; dec = dec)
 
 """
-    interiordiff!(x, y; dec = :default)
+    interiordiff!(v, x, y; dec = :default)
 
 In-place version of [`interiordiff`](@ref).
 """
 function interiordiff!(v::AbstractVector, x::BareInterval{T}, y::BareInterval{T}) where {T<:NumTypes}
-    if isinterior(x, y)
-        empty!(v)
-    elseif isdisjoint_interval(x, y)
-        resize!(v, 1)
-        @inbounds v[begin] = x
-    else
-        inter = intersect_interval(x, y)
-        if issubset_interval(y, x)
-            resize!(v, 2)
-            @inbounds v[begin] = _unsafe_bareinterval(T, inf(x), inf(inter))
-            @inbounds v[end] = _unsafe_bareinterval(T, sup(inter), sup(x))
-        elseif isweakless(x, inter)
-            resize!(v, 1)
-            @inbounds v[begin] = _unsafe_bareinterval(T, inf(x), inf(inter))
-        else
-            resize!(v, 1)
-            @inbounds v[begin] = _unsafe_bareinterval(T, sup(inter), sup(x))
-        end
-    end
+    h₁, h₂, _ = _interiordiff(x, y, nothing)
+    empty!(v)
+    isempty_interval(h₁) || push!(v, h₁)
+    isempty_interval(h₂) || push!(v, h₂)
     return v
 end
 interiordiff!(v::AbstractVector, x::BareInterval, y::BareInterval) = interiordiff!(v, promote(x, y)...)
@@ -193,31 +174,17 @@ interiordiff!(v::AbstractVector{<:AbstractVector}, x::AbstractVector{<:BareInter
     _interiordiff!(v, x, y, nothing)
 
 function interiordiff!(v::AbstractVector, x::Interval{T}, y::Interval{T}; dec = :default) where {T<:NumTypes}
-    if isinterior(x, y)
-        empty!(v)
-    elseif isdisjoint_interval(x, y)
+    d = min(decoration(x), decoration(y))
+    if d == ill # one of the inputs is an NaI
         resize!(v, 1)
-        @inbounds v[begin] = x
-    else
-        d = min(decoration(x), decoration(y))
-        t = isguaranteed(x) & isguaranteed(y)
-        inter = intersect_interval(x, y)
-        if issubset_interval(y, x)
-            resize!(v, 2)
-            r1 = _unsafe_bareinterval(T, inf(x), inf(inter))
-            @inbounds v[begin] = _set_decoration(_unsafe_interval(r1, d, t), dec)
-            r2 = _unsafe_bareinterval(T, sup(inter), sup(x))
-            @inbounds v[end] = _set_decoration(_unsafe_interval(r2, d, t), dec)
-        elseif isweakless(x, inter)
-            resize!(v, 1)
-            r1 = _unsafe_bareinterval(T, inf(x), inf(inter))
-            @inbounds v[begin] = _set_decoration(_unsafe_interval(r1, d, t), dec)
-        else
-            resize!(v, 1)
-            r2 = _unsafe_bareinterval(T, sup(inter), sup(x))
-            @inbounds v[begin] = _set_decoration(_unsafe_interval(r2, d, t), dec)
-        end
+        @inbounds v[begin] = nai(T)
+        return v
     end
+    t = isguaranteed(x) & isguaranteed(y)
+    h₁, h₂, _ = _interiordiff(_bareinterval(x), _bareinterval(y), nothing)
+    empty!(v)
+    isempty_interval(h₁) || push!(v, _set_decoration(_unsafe_interval(h₁, d, t), dec))
+    isempty_interval(h₂) || push!(v, _set_decoration(_unsafe_interval(h₂, d, t), dec))
     return v
 end
 interiordiff!(v::AbstractVector, x::Interval, y::Interval; dec = :default) = interiordiff!(v, promote(x, y)...; dec = dec)
@@ -232,8 +199,14 @@ function _interiordiff!(v::AbstractVector{<:AbstractVector}, x::AbstractVector,
     N == len || return throw(DimensionMismatch("x has length $N, y has length $len"))
 
     if any(t -> isdisjoint_interval(t[1], t[2]), zip(x, y))
-        resize!(v, 1)
-        @inbounds v[begin] = copy(x)
+        # mirror the filter of the general case: an empty or NaI component
+        # means that `x` does not describe a box
+        if any(t -> isempty_interval(t) | isnai(t), x)
+            empty!(v)
+        else
+            resize!(v, 1)
+            @inbounds v[begin] = copy(x)
+        end
     else
         resize!(v, 2N)
 
@@ -261,16 +234,16 @@ function _interiordiff!(v::AbstractVector{<:AbstractVector}, x::AbstractVector,
 end
 
 function _interiordiff(x::BareInterval{T}, y::BareInterval{T}, ::Nothing) where {T<:NumTypes}
-    isdisjoint_interval(x, y) && return (x, emptyinterval(BareInterval{T}), emptyinterval(BareInterval{T}))
-
     inter = intersect_interval(x, y)
+    isdisjoint_interval(x, y) && return (x, emptyinterval(BareInterval{T}), inter)
     isinterior(x, y) && return (emptyinterval(BareInterval{T}), emptyinterval(BareInterval{T}), inter)
-    isequal_interval(x, y) && return (_unsafe_bareinterval(T, inf(x), inf(x)), _unsafe_bareinterval(T, sup(x), sup(x)), inter)
-
-    inf(x) == inf(inter) && return (_unsafe_bareinterval(T, sup(inter), sup(x)), emptyinterval(BareInterval{T}), inter)
-    sup(x) == sup(inter) && return (_unsafe_bareinterval(T, inf(x), inf(inter)), emptyinterval(BareInterval{T}), inter)
 
-    return (_unsafe_bareinterval(T, inf(x), inf(y)), _unsafe_bareinterval(T, sup(y), sup(x)), inter)
+    # a piece of `x` survives on a given side only if the corresponding bound
+    # of `y` is finite, since an infinite bound means that the interior of `y`
+    # is unbounded on that side
+    h₁ = (inf(x) ≤ inf(y)) & (inf(y) != typemin(T)) ? _unsafe_bareinterval(T, inf(x), inf(y)) : emptyinterval(BareInterval{T})
+    h₂ = (sup(y) ≤ sup(x)) & (sup(y) != typemax(T)) ? _unsafe_bareinterval(T, sup(y), sup(x)) : emptyinterval(BareInterval{T})
+    return (h₁, h₂, inter)
 end
 _interiordiff(x::BareInterval, y::BareInterval, dec::Nothing) = _interiordiff(promote(x, y)..., dec)
 
@@ -284,15 +257,20 @@ end
 #
 
 function interval_diff(x::Interval{T}, y::Interval) where {T<:NumTypes}
-    isdisjoint_interval(x, y) && return [x]
+    d = min(decoration(x), decoration(y))
+    d == ill && return [nai(T)] # one of the inputs is an NaI
+    t = isguaranteed(x) & isguaranteed(y)
+    # the pieces are decorated `trv` (Section 11.7.1), including in the
+    # disjoint branch where `x` is returned with its decoration downgraded
+    isdisjoint_interval(x, y) && return [_unsafe_interval(_bareinterval(x), trv, t)]
     issubset_interval(x, y) && return Interval{T}[]
 
     intersection = intersect_interval(x, y)
-    inf(x) == inf(intersection) && return [interval(sup(intersection), sup(x))]
-    sup(x) == sup(intersection) && return [interval(inf(x), inf(intersection))]
+    inf(x) == inf(intersection) && return [_unsafe_interval(_unsafe_bareinterval(T, sup(intersection), sup(x)), trv, t)]
+    sup(x) == sup(intersection) && return [_unsafe_interval(_unsafe_bareinterval(T, inf(x), inf(intersection)), trv, t)]
 
     return [
-        interval(inf(x), inf(intersection)),
-        interval(sup(intersection), sup(x))
+        _unsafe_interval(_unsafe_bareinterval(T, inf(x), inf(intersection)), trv, t),
+        _unsafe_interval(_unsafe_bareinterval(T, sup(intersection), sup(x)), trv, t)
     ]
 end
diff --git a/test/interval_tests/display.jl b/test/interval_tests/display.jl
index bb035c0f..d9e7c095 100644
--- a/test/interval_tests/display.jl
+++ b/test/interval_tests/display.jl
@@ -50,7 +50,7 @@ setprecision(BigFloat, 256) do
 
             @test sprint(show, MIME("text/plain"), emptyinterval(BareInterval{Float64})) == "∅"
 
-            @test sprint(show, MIME("text/plain"), a) == "0.65 ± 0.65"
+            @test sprint(show, MIME("text/plain"), a) == "0.65 ± 0.650001"
             @test sprint(show, MIME("text/plain"), large_expo) ==
                 "(5.0e+123456788 ± 5.0e+123456788)₂₅₆"
         end
@@ -156,9 +156,9 @@ setprecision(BigFloat, 256) do
 
                 @test sprint(show, MIME("text/plain"), a)    == "(1.5 ± 0.5)_com"
                 @test sprint(show, MIME("text/plain"), a_NG) == "(1.5 ± 0.5)_com_NG"
-                @test sprint(show, MIME("text/plain"), b)    == "(0.65 ± 0.65)_com"
-                @test sprint(show, MIME("text/plain"), b32)  == "(0.65f0 ± 0.65f0)_com"
-                @test sprint(show, MIME("text/plain"), b16)  == "(Float16(0.65) ± Float16(0.65))_com"
+                @test sprint(show, MIME("text/plain"), b)    == "(0.65 ± 0.650001)_com"
+                @test sprint(show, MIME("text/plain"), b32)  == "(0.65f0 ± 0.650001f0)_com"
+                @test sprint(show, MIME("text/plain"), b16)  == "(Float16(0.65) ± Float16(0.6504))_com"
                 @test sprint(show, MIME("text/plain"), br)   == "(1//10 ± 6//5)_com"
                 @test sprint(show, MIME("text/plain"), c)    == "(1.79769e+308 ± ∞)_dac"
                 @test sprint(show, MIME("text/plain"), cr)   == "(9223372036854775807//1 ± ∞)_dac"
@@ -174,9 +174,9 @@ setprecision(BigFloat, 256) do
 
                 @test sprint(show, MIME("text/plain"), a)    == "1.5 ± 0.5"
                 @test sprint(show, MIME("text/plain"), a_NG) == "(1.5 ± 0.5)_NG"
-                @test sprint(show, MIME("text/plain"), b)    == "0.65 ± 0.65"
-                @test sprint(show, MIME("text/plain"), b32)  == "0.65f0 ± 0.65f0"
-                @test sprint(show, MIME("text/plain"), b16)  == "Float16(0.65) ± Float16(0.65)"
+                @test sprint(show, MIME("text/plain"), b)    == "0.65 ± 0.650001"
+                @test sprint(show, MIME("text/plain"), b32)  == "0.65f0 ± 0.650001f0"
+                @test sprint(show, MIME("text/plain"), b16)  == "Float16(0.65) ± Float16(0.6504)"
                 @test sprint(show, MIME("text/plain"), br)   == "1//10 ± 6//5"
                 @test sprint(show, MIME("text/plain"), c)    == "1.79769e+308 ± ∞"
                 @test sprint(show, MIME("text/plain"), cr)   == "9223372036854775807//1 ± ∞"