Rework test suite

ext/IntervalArithmeticLinearAlgebraExt.jl

@@ -235,6 +235,20 @@ function LinearAlgebra.mul!(C::AbstractMatrix{<:RealOrComplexI}, A::AbstractVecO
     return _mul!(IntervalArithmetic.default_matmul(), C, A, B, α, β)
 end
 
+# disambiguate from the LinearAlgebra methods for structured matrices, needed for Julia v1.10
+const _StructuredMatrix = Union{LinearAlgebra.Bidiagonal,LinearAlgebra.Diagonal,LinearAlgebra.SymTridiagonal,LinearAlgebra.Tridiagonal}
+
+for (CT, AT, BT) ∈ (
+        (:(AbstractMatrix{<:RealOrComplexI}), :AbstractMatrix,    :_StructuredMatrix),
+        (:(AbstractMatrix{<:RealOrComplexI}), :_StructuredMatrix, :AbstractMatrix),
+        (:(AbstractMatrix{<:RealOrComplexI}), :_StructuredMatrix, :_StructuredMatrix),
+        (:(AbstractVector{<:RealOrComplexI}), :_StructuredMatrix, :AbstractVector))
+    @eval function LinearAlgebra.mul!(C::$CT, A::$AT, B::$BT, α::Number, β::Number)
+        size(A, 2) == size(B, 1) || return throw(DimensionMismatch("The number of columns of A must match the number of rows of B."))
+        return _mul!(IntervalArithmetic.default_matmul(), C, A, B, α, β)
+    end
+end
+
 #
 
 LinearAlgebra.mul!(C::AbstractVecOrMat{<:RealOrComplexI}, A::AbstractMatrix{<:RealOrComplexI}, B::AbstractVecOrMat{<:RealOrComplexI}) =

src/IntervalArithmetic.jl

@@ -272,7 +272,7 @@ function sample(rng::Random.AbstractRNG, x::Interval{T}) where {T<:NumTypes}
     β = rand(rng, float(T))
     lo = ifelse(lo == typemin(T), _value_min(T), lo)
     hi = ifelse(hi == typemax(T), _value_max(T), hi)
-    val = convert(T, (1 - β) * lo + β * hi)
+    val = _value_round(T, (1 - β) * lo + β * hi)
     val = ifelse(val < lo, lo, val)
     val = ifelse(val > hi, hi, val)
     return val
@@ -284,6 +284,10 @@ _value_max(::Type{T}) where {T<:AbstractFloat} = floatmax(T)
 _value_min(::Type{Rational{T}}) where {T<:Integer} = convert(Rational{T}, typemin(T))
 _value_max(::Type{Rational{T}}) where {T<:Integer} = convert(Rational{T}, typemax(T))
 
+_value_round(::Type{T}, val) where {T<:AbstractFloat} = convert(T, val)
+# the exact fraction of a float may overflow the integer type
+_value_round(::Type{Rational{T}}, val) where {T<:Integer} = rationalize(T, val)
+
     export sample
 
 #

src/display.jl

@@ -328,8 +328,11 @@ function _round_string(x::AbstractFloat, sigdigits::Int)
     return str
 end
 
-_round_string(x::AbstractFloat, sigdigits::Int, mode::RoundingMode) =
-    _round_string(round(x, mode; sigdigits = sigdigits), sigdigits)
+# print as Julia Base whenever the shortest representation fits within `sigdigits`
+function _round_string(x::AbstractFloat, sigdigits::Int, mode::RoundingMode)
+    _count_sigdigits(string(x)) ≤ sigdigits && return _round_string(x, sigdigits)
+    return _round_string(round(x, mode; sigdigits = sigdigits), sigdigits)
+end
 
 _count_sigdigits(s::AbstractString) = length(replace(split(s, r"[eE]")[1], '-' => "", '.' => "", r"^0+" => ""))
 

src/intervals/arithmetic/trigonometric.jl

@@ -501,26 +501,34 @@ Base.atan(x::Complex{Interval{T}}) where {T<:NumTypes} =
     acot(::BareInterval)
     acot(::Interval)
 
-This function is not part of the IEEE Standard 1788-2015.
+This function is not part of the IEEE Standard 1788-2015. It uses the convention
+`acot(0) = π/2` and the jump at 0 is decorated like the branch cut of the
+two-argument `atan` described as `atan2` in the IEEE Standard 1788-2015
+(Table 9.1).
 """
 function Base.acot(x::BareInterval{T}) where {T<:AbstractFloat}
     isempty_interval(x) && return x
     lo, hi = bounds(x)
-    if lo < 0 < hi || lo == hi == 0
+    if lo < 0 && hi ≥ 0
         HALF_PI_HI = sup(_half_pi(T))
         return _unsafe_bareinterval(T, -HALF_PI_HI, HALF_PI_HI)
     elseif lo == 0
+        hi == 0 && return _half_pi(T)
         return @round(T, acot(hi), +sup(_half_pi(T)))
-    elseif hi == 0
-        return @round(T, -sup(_half_pi(T)), acot(lo))
     else
         return @round(T, acot(hi), acot(lo))
     end
 end
 
 Base.acot(x::BareInterval{<:Rational}) = acot(float(x))
 
-# automatically defined for `Interval` since it is a subtype of `Real`
+function Base.acot(x::Interval)
+    bx = bareinterval(x)
+    r = acot(bx)
+    d = min(decoration(x), decoration(r))
+    d = min(d, ifelse(in_interval(0, bx), ifelse(inf(bx) < 0, def, dac), d))
+    return _unsafe_interval(r, d, isguaranteed(x))
+end
 
 """
     atan(::BareInterval, ::BareInterval)
@@ -588,7 +596,7 @@ function Base.atan(y::Interval, x::Interval)
             ifelse(in_interval(0, by),
                    ifelse(in_interval(0, bx),
                           trv,
-                          ifelse(sup(bx) < 0, ifelse(inf(by) < 0, def, d), d)),
+                          ifelse(sup(bx) < 0, ifelse(inf(by) < 0, def, dac), d)),
                    d))
     t = isguaranteed(y) & isguaranteed(x)
     return _unsafe_interval(r, d, t)

src/intervals/rounding.jl

@@ -138,15 +138,19 @@ for f ∈ (:cbrt, :exp, :exp2, :exp10, :expm1, # exponential
         end
     end
 
+    # if CRlibm has no routine, dispatch falls back to the shared `:correct` methods above
     if f ∈ CRlibm.functions
         @eval _fround(::typeof($f), ::IntervalRounding{:correct}, x::Float16, r::RoundingMode) = Float16(_fround($f, Float64(x), r), r)
         @eval _fround(::typeof($f), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, r::RoundingMode) = CRlibm.$f(x, r)
     end
 
-    @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Down}) = prevfloat(CoreMath.$coremath_f(x))
-    @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Up}) = nextfloat(CoreMath.$coremath_f(x))
+    # if CoreMath has no routine, dispatch falls back to the shared `:correct` methods above
+    if f ∈ CoreMath.univariate_functions
+        @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Down}) = prevfloat(CoreMath.$coremath_f(x))
+        @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Up}) = nextfloat(CoreMath.$coremath_f(x))
+    end
 
-    @eval _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, r::RoundingMode) = $f(x)
+    @eval _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x)
 end
 
 # 2-argument functions