Release 1.0 (#728)

* Release 1.0

* Polish configuration system

* Polish display of unbounded rational intervals

* Fix test

* Reintroduce `@interval` with two arguments and last cleanup

* Fix docstrings

* Fix `@interval`

* Use `exact` instead of `ExactReal` as a constructor

* Run tests on release candidates

Author: OlivierHnt

.github/workflows/CI.yml

@@ -13,6 +13,7 @@ jobs:
         version:
           - '1.10'
           - '1'
+          - 'pre'
           - 'nightly'
         os:
           - ubuntu-latest
@@ -25,8 +26,8 @@ jobs:
           - os: macOS-latest
             arch: x86
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

Project.toml

@@ -1,7 +1,7 @@
 name = "IntervalArithmetic"
 uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 repo = "https://github.com/JuliaIntervals/IntervalArithmetic.jl.git"
-version = "0.23.0"
+version = "1.0.0"
 
 [deps]
 CRlibm = "96374032-68de-5a5b-8d9e-752f78720389"

ext/IntervalArithmeticLinearAlgebraExt.jl

@@ -13,10 +13,12 @@ end
 
 # contructor for `UniformScaling`
 
-IntervalArithmetic.interval(::Type{T}, J::LinearAlgebra.UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) where {T} =
-    LinearAlgebra.UniformScaling(interval(T, J.λ, d; format = format))
-IntervalArithmetic.interval(J::LinearAlgebra.UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) =
-    LinearAlgebra.UniformScaling(interval(J.λ, d; format = format))
+IntervalArithmetic._interval_infsup(::Type{T}, J::LinearAlgebra.UniformScaling, H::LinearAlgebra.UniformScaling, d::IntervalArithmetic.Decoration) where {T<:IntervalArithmetic.NumTypes} =
+    LinearAlgebra.UniformScaling(IntervalArithmetic._interval_infsup(T, J.λ, H.λ, d))
+
+IntervalArithmetic._infer_numtype(J::LinearAlgebra.UniformScaling) = numtype(eltype(J))
+
+IntervalArithmetic.exact(J::LinearAlgebra.UniformScaling) = LinearAlgebra.UniformScaling(exact(J.λ))
 
 # by-pass generic `opnorm` from LinearAlgebra to prevent NG flag
 
@@ -139,32 +141,22 @@ end
 
 # matrix multiplication
 
-function IntervalArithmetic.configure_matmul(matmul::Symbol)
-    matmul ∈ (:slow, :fast) || return throw(ArgumentError("only the matrix multiplication mode `:slow` and `:fast` are available"))
-
-    @eval begin
-        function LinearAlgebra.mul!(C::AbstractVector{<:RealOrComplexI}, A::AbstractVecOrMat, B::AbstractVector, α::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 IntervalArithmetic._mul!(IntervalArithmetic.MatMulMode{$(QuoteNode(matmul))}(), C, A, B, α, β)
-        end
-
-        function LinearAlgebra.mul!(C::AbstractMatrix{<:RealOrComplexI}, A::AbstractVecOrMat, B::AbstractVecOrMat, α::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 IntervalArithmetic._mul!(IntervalArithmetic.MatMulMode{$(QuoteNode(matmul))}(), C, A, B, α, β)
-        end
-    end
-
-    return matmul
+function LinearAlgebra.mul!(C::AbstractVector{<:RealOrComplexI}, A::AbstractVecOrMat, B::AbstractVector, α::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
 
-IntervalArithmetic.configure_matmul(:fast)
+function LinearAlgebra.mul!(C::AbstractMatrix{<:RealOrComplexI}, A::AbstractVecOrMat, B::AbstractVecOrMat, α::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
 
 #
 
 LinearAlgebra.mul!(C::AbstractVecOrMat{<:RealOrComplexI}, A::AbstractMatrix{<:RealOrComplexI}, B::AbstractVecOrMat{<:RealOrComplexI}) =
     LinearAlgebra.mul!(C, A, B, interval(true), interval(false))
 
-function IntervalArithmetic._mul!(::IntervalArithmetic.MatMulMode{:slow}, C, A::AbstractMatrix, B::AbstractVecOrMat, α, β)
+function _mul!(::IntervalArithmetic.MatMulMode{:slow}, C, A::AbstractMatrix, B::AbstractVecOrMat, α, β)
     if iszero(α)
         if iszero(β)
             C .= zero(eltype(C))
@@ -221,14 +213,14 @@ end
 # fast matrix multiplication
 # Note: Rump's algorithm
 
-function IntervalArithmetic._mul!(::IntervalArithmetic.MatMulMode{:fast}, C, A, B, α, β)
+function _mul!(::IntervalArithmetic.MatMulMode{:fast}, C, A, B, α, β)
     if Int == Int32
         @info "Fast multiplication is not supported on 32-bit systems, using the slow version"
-        return IntervalArithmetic._mul!(IntervalArithmetic.MatMulMode{:slow}(), C, A, B, α, β)
+        return _mul!(IntervalArithmetic.MatMulMode{:slow}(), C, A, B, α, β)
     else
         numtype(eltype(C)) <: Union{Float16,Float32,Float64} && return _fastmul!(C, A, B, α, β)
         @info "Fast multiplication is only supported for `Union{Float16,Float32,Float64}`, using the slow version"
-        return IntervalArithmetic._mul!(IntervalArithmetic.MatMulMode{:slow}(), C, A, B, α, β)
+        return _mul!(IntervalArithmetic.MatMulMode{:slow}(), C, A, B, α, β)
     end
 end
 
@@ -413,11 +405,11 @@ function ___mul(A::AbstractMatrix{T}, B::AbstractVecOrMat{Interval{T}}) where {T
 
     stride_A  = _to_stride_64(A)
     stride_mB = _to_stride_64(mB)
-    C₁ = IntervalArithmetic._sub_round.(
+    C₁ = IntervalArithmetic._fround.(-,
         T.(_call_gem_openblas!(cache_2, stride_A, stride_mB, RoundDown), RoundDown),
         T.(rC, RoundUp),
         RoundDown)
-    C₂ = cache_2; C₂ .= IntervalArithmetic._add_round.(
+    C₂ = cache_2; C₂ .= IntervalArithmetic._fround.(+,
         T.(_call_gem_openblas!(cache_2, stride_A, stride_mB, RoundUp), RoundUp),
         T.(rC, RoundUp),
         RoundUp)
@@ -435,11 +427,11 @@ function ___mul(A::AbstractMatrix{Interval{T}}, B::AbstractVecOrMat{T}) where {T
 
     stride_mA = _to_stride_64(mA)
     stride_B  = _to_stride_64(B)
-    C₁ = IntervalArithmetic._sub_round.(
+    C₁ = IntervalArithmetic._fround.(-,
         T.(_call_gem_openblas!(cache_2, stride_mA, stride_B, RoundDown), RoundDown),
         T.(rC, RoundUp),
         RoundDown)
-    C₂ = cache_2; C₂ .= IntervalArithmetic._add_round.(
+    C₂ = cache_2; C₂ .= IntervalArithmetic._fround.(+,
         T.(_call_gem_openblas!(cache_2, stride_mA, stride_B, RoundUp), RoundUp),
         T.(rC, RoundUp),
         RoundUp)
@@ -455,76 +447,33 @@ function ___mul(A::AbstractMatrix{Interval{T}}, B::AbstractVecOrMat{Interval{T}}
     cache_2 = zeros(Float64, size(A, 1), size(B, 2))
 
     absmA_rB = _call_gem_openblas!(cache_1, _to_stride_64(abs.(mA)), _to_stride_64(rB), RoundUp)
-    U = rB; U .= IntervalArithmetic._add_round.(abs.(mB), rB, RoundUp)
+    U = rB; U .= IntervalArithmetic._fround.(+, abs.(mB), rB, RoundUp)
     rA_U = _call_gem_openblas!(cache_2, _to_stride_64(rA), _to_stride_64(U), RoundUp)
 
     cache_3 = zeros(Float64, size(A, 1), size(B, 2))
 
     stride_mA = _to_stride_64(mA)
     stride_mB = _to_stride_64(mB)
-    C₁ = IntervalArithmetic._sub_round.(
+    C₁ = IntervalArithmetic._fround.(-,
         T.(_call_gem_openblas!(cache_3, stride_mA, stride_mB, RoundDown), RoundDown),
-        IntervalArithmetic._add_round.(T.(absmA_rB, RoundUp), T.(rA_U, RoundUp), RoundUp),
+        IntervalArithmetic._fround.(+, T.(absmA_rB, RoundUp), T.(rA_U, RoundUp), RoundUp),
         RoundDown)
-    C₂ = cache_3; C₂ .= IntervalArithmetic._add_round.(
+    C₂ = cache_3; C₂ .= IntervalArithmetic._fround.(+,
         T.(_call_gem_openblas!(cache_3, stride_mA, stride_mB, RoundUp), RoundUp),
-        IntervalArithmetic._add_round.(T.(absmA_rB, RoundUp), T.(rA_U, RoundUp), RoundUp),
+        IntervalArithmetic._fround.(+, T.(absmA_rB, RoundUp), T.(rA_U, RoundUp), RoundUp),
         RoundUp)
 
     return C₁, C₂
 end
 
-# function ___mul(A::AbstractMatrix{Interval{T}}, B::AbstractVecOrMat{Interval{T}}) where {T<:AbstractFloat}
-#     k = size(A, 2)
-#     u2 = eps(T) # twice the unit roundoff
-#     @assert (2k + 2) * u2 ≤ 1
-
-#     mA, rA = _vec_or_mat_midradius(A)
-#     mB, rB = _vec_or_mat_midradius(B)
-
-#     cache_1 = zeros(T, size(A, 1), size(B, 2))
-#     cache_2 = zeros(T, size(A, 1), size(B, 2))
-#     mC, μ = _fused_matmul!(cache_1, cache_2, mA, rA, mB, rB)
-
-#     γ = IntervalArithmetic._add_round.(IntervalArithmetic._mul_round.(convert(T, k + 1), eps.(μ), RoundUp), IntervalArithmetic._mul_round(IntervalArithmetic._inv_round(u2, RoundUp), floatmin(T), RoundUp), RoundUp)
-
-#     U = mA; U .= IntervalArithmetic._add_round.(abs.(mA), rA, RoundUp)
-#     V = mB; V .= IntervalArithmetic._add_round.(abs.(mB), rB, RoundUp)
-
-#     cache_3 = zeros(Float64, size(A, 1), size(B, 2))
-#     rC = IntervalArithmetic._add_round.(IntervalArithmetic._sub_round.(T.(_call_gem_openblas!(cache_3, _to_stride_64(U), _to_stride_64(V), RoundUp), RoundUp), μ, RoundUp), 2 .* γ, RoundUp)
-
-#     C₁ = μ; C₁ .= IntervalArithmetic._sub_round.(mC, rC, RoundDown)
-#     C₂ = γ; C₂ .= IntervalArithmetic._add_round.(mC, rC, RoundUp)
-
-#     return C₁, C₂
-# end
-
-# function _fused_matmul!(mC, μ, mA, rA, mB, rB)
-#     Threads.@threads for j ∈ axes(mB, 2)
-#         for l ∈ axes(mA, 2)
-#             @inbounds for i ∈ axes(mA, 1)
-#                 a, c = mA[i,l], rA[i,l]
-#                 b, d = mB[l,j], rB[l,j]
-#                 e = sign(a) * min(abs(a), c)
-#                 f = sign(b) * min(abs(b), d)
-#                 p = a*b + e*f
-#                 mC[i,j] += p
-#                 μ[i,j] += abs(p)
-#             end
-#         end
-#     end
-#     return mC, μ
-# end
-
 _to_stride_64(A::StridedArray{Float64}) = A
 _to_stride_64(A::StridedArray{<:AbstractFloat}) = Float64.(A, RoundUp)
 _to_stride_64(A::AbstractVector) = _to_stride_64(Vector(A))
 _to_stride_64(A::AbstractMatrix) = _to_stride_64(Matrix(A))
 
 function _vec_or_mat_midradius(A::AbstractVecOrMat{Interval{T}}) where {T<:AbstractFloat}
-    mA = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp)
-    rA = IntervalArithmetic._sub_round.(mA, inf.(A), RoundUp)
+    mA = IntervalArithmetic._fround.(/, IntervalArithmetic._fround.(+, inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp)
+    rA = IntervalArithmetic._fround.(-, mA, inf.(A), RoundUp)
     return mA, rA
 end