Extend configure for flavour, power mode, matmul mode and default numtype

Author: OlivierHnt

src/IntervalArithmetic.jl

@@ -1,7 +1,57 @@
 """
     IntervalArithmetic
 
-Library for validated numerics using interval arithmetic.
+Library for validated numerics using interval arithmetic, offering tools to
+rigorously bound errors in numerical computations by representing values as
+intervals and ensuring that computed results contain the true value. It provides
+accurate, and efficient methods, ideal for scientific computing,
+computer-assisted proofs, and any domain requiring certified numerical results.
+
+Key features:
+
+- **Bound Type**: The default numerical type used to represent the bounds of the
+  intervals. The default is `Float64`, but other subtypes of
+  [`IntervalArithmetic.NumTypes`](@ref) can be used to adjust precision.
+
+- **Flavor**: The interval representation that adhere to the IEEE Standard
+  1788-2015. By default, it uses the set-based flavor, which excludes infinity
+  to be part of an interval. Learn more: [`IntervalArithmetic.Flavor`](@ref).
+
+- **Interval Rounding**: Controls the rounding behavior for interval arithmetic
+  operations. By default, the library employs correct rounding to ensure that
+  bounds are computed as tightly as possible. Learn more:
+  [`IntervalArithmetic.IntervalRounding`](@ref).
+
+- **Power mode**: A performance setting for power operations. The default mode
+  uses an efficient algorithm prioritizing fast computation, but it can be
+  adjusted for more precise, slower calculations if needed. Learn more:
+  [`IntervalArithmetic.PowerMode`](@ref).
+
+- **Matrix Multiplication mode**: A performance setting for matrix
+  multiplication operations. The default mode uses an efficient algorithm
+  prioritizing fast computation, but it can be changed to use the standard
+  definition of matrix multiplication. Learn more:
+  [`IntervalArithmetic.MatMulMode`](@ref).
+
+The default behaviors described above can be configured via
+[`IntervalArithmetic.configure`](@ref).
+
+**Display Settings**: controls how intervals are displayed. By default,
+intervals are shown using the standard mathematical notation ``[a, b]``, along
+with decorations and up to 6 significant digits. Learn more:
+[`setdisplay`](@ref).
+
+# Usage
+
+```julia
+using IntervalArithmetic
+
+# Create an interval
+x = interval(1.0, 2.0)
+
+# Perform a rigorous computation
+x + exp(interval(π))
+```
 
 Learn more: https://github.com/JuliaIntervals/IntervalArithmetic.jl
 """
@@ -32,10 +82,6 @@ include("display.jl")
 
 #
 
-include("symbols.jl")
-
-#
-
 import LinearAlgebra
 import OpenBLASConsistentFPCSR_jll # 32-bit systems are not supported
 
@@ -50,6 +96,27 @@ include("matmul.jl")
 
 #
 
+function configure_numtype(numtype::Type{<:NumTypes})
+    @eval promote_numtype(::Type{T}, ::Type{S}) where {T,S} = promote_type($numtype, numtype(T), numtype(S))
+    @eval macro interval(expr)
+        return _wrap_interval($numtype, expr)
+    end
+    @eval _parse(str::AbstractString) = parse(Interval{$numtype}, str)
+    @eval emptyinterval() = emptyinterval(Interval{$numtype})
+    @eval entireinterval() = entireinterval(Interval{$numtype})
+    @eval nai() = nai(Interval{$numtype})
+    return numtype
+end
+
+function configure_flavor(flavor::Symbol)
+    flavor == :set_based || return throw(ArgumentError("only the interval flavor `:set_based` is supported and implemented"))
+    @eval zero_times_infinity(::Type{T}) where {T<:NumTypes} = zero_times_infinity(Flavor{$(QuoteNode(flavor))}(), T)
+    @eval div_by_thin_zero(x::BareInterval) = div_by_thin_zero(Flavor{$(QuoteNode(flavor))}(), x)
+    @eval contains_infinity(x::BareInterval) = contains_infinity(Flavor{$(QuoteNode(flavor))}(), x)
+    @eval is_valid_interval(a::Real, b::Real) = is_valid_interval(Flavor{$(QuoteNode(flavor))}(), a, b)
+    return flavor
+end
+
 function configure_rounding(rounding::Symbol)
     rounding ∈ (:correct, :none) || return throw(ArgumentError("only the rounding mode `:correct` and `:none` are available"))
 
@@ -86,15 +153,80 @@ function configure_rounding(rounding::Symbol)
     return rounding
 end
 
-function configure(; rounding::Symbol=:correct)
+function configure_power(power::Symbol)
+    power ∈ (:slow, :fast) || return throw(ArgumentError("only the power mode `:slow` and `:fast` are available"))
+
+    for f ∈ (:_select_pow, :_select_pown)
+        @eval $f(x, y) = $f(PowerMode{$(QuoteNode(power))}(), x, y)
+    end
+
+    return power
+end
+
+function 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 _mul!(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 _mul!(MatMulMode{$(QuoteNode(matmul))}(), C, A, B, α, β)
+        end
+    end
+
+    return matmul
+end
+
+"""
+    configure(; numtype=Float64, flavor=:set_based, rounding=:correct, power=:fast, matmul=:fast)
+
+Configure the default behavior for:
+
+- **Bound Type**: The default numerical type used to represent the bounds of the
+  intervals. The default is `Float64`, but other subtypes of
+  [`IntervalArithmetic.NumTypes`](@ref) can be used to adjust precision.
+
+- **Flavor**: The interval representation that adhere to the IEEE Standard
+  1788-2015. By default, it uses the set-based flavor, which excludes infinity
+  to be part of an interval. Learn more: [`IntervalArithmetic.Flavor`](@ref).
+
+- **Interval Rounding**: Controls the rounding behavior for interval arithmetic
+  operations. By default, the library employs correct rounding to ensure that
+  bounds are computed as tightly as possible. Learn more:
+  [`IntervalArithmetic.IntervalRounding`](@ref).
+
+- **Power mode**: A performance setting for power operations. The default mode
+  uses an efficient algorithm prioritizing fast computation, but it can be
+  adjusted for more precise, slower calculations if needed. Learn more:
+  [`IntervalArithmetic.PowerMode`](@ref).
+
+- **Matrix Multiplication mode**: A performance setting for matrix
+  multiplication operations. The default mode uses an efficient algorithm
+  prioritizing fast computation, but it can be changed to use the standard
+  definition of matrix multiplication. Learn more:
+  [`IntervalArithmetic.MatMulMode`](@ref).
+"""
+function configure(; numtype::Type{<:NumTypes}=Float64, flavor::Symbol=:set_based, rounding::Symbol=:correct, power::Symbol=:fast, matmul::Symbol=:fast)
+    configure_numtype(numtype)
+    configure_flavor(flavor)
     configure_rounding(rounding)
-    return rounding
+    configure_power(power)
+    configure_matmul(matmul)
+    return numtype, flavor, rounding, power, matmul
 end
 
 configure()
 
 #
 
+include("symbols.jl")
+
+#
+
 bareinterval(::Type{BigFloat}, a::AbstractIrrational) = _unsafe_bareinterval(BigFloat, a, a)
 
 # Note: generated functions must be defined after all the methods they use

src/intervals/arithmetic/absmax.jl

@@ -27,7 +27,7 @@ Implement the square absolute value.
     This function calls `^` internally, hence it depends on
     `IntervalArithmetic.power_mode()`.
 """
-Base.abs2(x::BareInterval) = _select_pown(power_mode(), x, 2) # not in the IEEE Standard 1788-2015
+Base.abs2(x::BareInterval) = _select_pown(x, 2) # not in the IEEE Standard 1788-2015
 
 function Base.abs2(x::Interval)
     r = abs2(bareinterval(x))

src/intervals/arithmetic/power.jl

@@ -20,7 +20,7 @@ Available mode types:
 """
 struct PowerMode{T} end
 
-power_mode() = PowerMode{:fast}()
+#
 
 """
     ^(x::BareInterval, y::BareInterval)
@@ -49,13 +49,13 @@ Interval{Float64}(-Inf, Inf, trv)
 ```
 """
 function Base.:^(x::BareInterval, y::BareInterval)
-    isthininteger(y) || return _select_pow(power_mode(), x, y)
-    return _select_pown(power_mode(), x, Integer(sup(y)))
+    isthininteger(y) || return _select_pow(x, y)
+    return _select_pown(x, Integer(sup(y)))
 end
 
 function Base.:^(x::Interval, y::Interval)
-    isthininteger(y) || return _select_pow(power_mode(), x, y)
-    r = _select_pown(power_mode(), x, Integer(sup(y)))
+    isthininteger(y) || return _select_pow(x, y)
+    r = _select_pown(x, Integer(sup(y)))
     d = min(decoration(r), decoration(y))
     t = isguaranteed(r) & isguaranteed(y)
     return _unsafe_interval(bareinterval(r), d, t)
@@ -65,7 +65,7 @@ Base.:^(x::Interval, n::Integer) = ^(x, n//one(n))
 Base.:^(x::Interval, y::Rational) = ^(x, convert(Interval{typeof(y)}, y))
 
 function Base.:^(x::Complex{<:Interval}, y::Complex{<:Interval})
-    if isreal(x) && isthininteger(y)
+    if isthinzero(imag(x)) && isthininteger(y)
         r = real(x) ^ real(y)
         d = min(decoration(x), decoration(y), decoration(r))
         t = isguaranteed(x) & isguaranteed(y)
@@ -96,7 +96,7 @@ Base.:^(x::Complex{<:Interval}, y::Interval) = ^(promote(x, y)...)
 Base.:^(x::Interval, y::Complex{<:Interval}) = ^(promote(x, y)...)
 
 # overwrite behaviour for small integer powers from https://github.com/JuliaLang/julia/pull/24240
-Base.literal_pow(::typeof(^), x::Interval, ::Val{n}) where {n} = _select_pown(power_mode(), x, n)
+Base.literal_pow(::typeof(^), x::Interval, ::Val{n}) where {n} = _select_pown(x, n)
 Base.literal_pow(::typeof(^), x::Complex{<:Interval}, ::Val{n}) where {n} = ^(x, interval(n))
 
 # helper functions for power
@@ -310,9 +310,9 @@ Compute the hypotenuse.
 
 Implement the `hypot` function of the IEEE Standard 1788-2015 (Table 10.5).
 """
-Base.hypot(x::BareInterval, y::BareInterval) = sqrt(_select_pown(power_mode(), x, 2) + _select_pown(power_mode(), y, 2))
+Base.hypot(x::BareInterval, y::BareInterval) = sqrt(_select_pown(x, 2) + _select_pown(y, 2))
 
-Base.hypot(x::Interval, y::Interval) = sqrt(_select_pown(power_mode(), x, 2) + _select_pown(power_mode(), y, 2))
+Base.hypot(x::Interval, y::Interval) = sqrt(_select_pown(x, 2) + _select_pown(y, 2))
 
 """
     fastpow(x, y)

src/intervals/construction.jl

@@ -8,35 +8,6 @@ Constant for the supported types of interval bounds. This is set to
 """
 const NumTypes = Union{Rational,AbstractFloat}
 
-"""
-    default_numtype()
-
-Return the default bound type used in [`promote_numtype`](@ref). By default,
-`default_numtype()` is set to `Float64`. It can be modified by redefining the
-function, however it should be set to a concrete subtype of [`NumTypes`](@ref).
-
-# Examples
-
-```jldoctest
-julia> IntervalArithmetic.default_numtype() = Float32
-
-julia> typeof(interval(1, 2))
-Interval{Float32}
-
-julia> typeof(interval(1, big(2)))
-Interval{BigFloat}
-
-julia> IntervalArithmetic.default_numtype() = Float64
-
-julia> typeof(interval(1, 2))
-Interval{Float64}
-
-julia> typeof(interval(1, big(2)))
-Interval{BigFloat}
-```
-"""
-default_numtype() = Float64
-
 """
     promote_numtype(T, S)
 
@@ -45,12 +16,11 @@ Return the bound type used to construct intervals. The bound type is given by
 when `T` is a `Rational{R}` and `S` is an `AbstractIrrational` (or vice-versa),
 in which case the bound type is given by `Rational{promote_type(R, Int64)}`. In
 all other cases, the bound type is given by
-`promote_type(default_numtype(), T, S)`.
+`promote_type([DEFAULT BOUND TYPE], T, S)`.
 """
 promote_numtype(::Type{T}, ::Type{S}) where {T<:NumTypes,S<:NumTypes} = promote_type(T, S)
 promote_numtype(::Type{T}, ::Type{S}) where {T<:NumTypes,S} = promote_type(numtype(T), numtype(S))
 promote_numtype(::Type{T}, ::Type{S}) where {T,S<:NumTypes} = promote_type(numtype(T), numtype(S))
-promote_numtype(::Type{T}, ::Type{S}) where {T,S} = promote_type(default_numtype(), numtype(T), numtype(S))
 
 promote_numtype(::Type{Rational{T}}, ::Type{<:AbstractIrrational}) where {T<:Integer} = Rational{promote_type(T, Int64)}
 promote_numtype(::Type{<:AbstractIrrational}, ::Type{Rational{T}}) where {T<:Integer} = Rational{promote_type(T, Int64)}
@@ -659,10 +629,6 @@ julia> @interval Float64 sin(1) exp(1)
 Interval{Float64}(0.8414709848078965, 2.7182818284590455, com)
 ```
 """
-macro interval(expr)
-    return _wrap_interval(default_numtype(), expr)
-end
-
 macro interval(T, expr)
     return _wrap_interval(T, expr)
 end

src/intervals/flavor.jl

@@ -28,9 +28,6 @@ Some flavors `F` include:
 # Examples
 
 ```jldoctest
-julia> IntervalArithmetic.default_flavor()
-IntervalArithmetic.Flavor{:set_based}()
-
 julia> IntervalArithmetic.is_valid_interval(Inf, Inf)
 false
 
@@ -49,50 +46,37 @@ true
 """
 struct Flavor{F} end
 
-"""
-    default_flavor()
-
-Return the default flavor used to handle edge cases.
-"""
-default_flavor() = Flavor{:set_based}()
+#
 
 """
-    zero_times_infinity([F::Flavor=default_flavor()], T<:NumTypes)
+    zero_times_infinity([F::Flavor,] T<:NumTypes)
 
 For the given flavor `F`, return ``0 \\times \\infty`` as an instance of type
 `T`.
 """
 zero_times_infinity(::Flavor{:set_based}, ::Type{T}) where {T<:NumTypes} = zero(T)
 
-zero_times_infinity(::Type{T}) where {T<:NumTypes} = zero_times_infinity(default_flavor(), T)
-
 """
-    div_by_thin_zero([F::Flavor=default_flavor()], x::BareInterval)
+    div_by_thin_zero([F::Flavor,] x::BareInterval)
 
 For the given flavor `F`, divide `x` by the interval containing only ``0``.
 """
 div_by_thin_zero(::Flavor{:set_based}, ::BareInterval{T}) where {T<:NumTypes} =
     emptyinterval(BareInterval{T})
 
-div_by_thin_zero(x::BareInterval) = div_by_thin_zero(default_flavor(), x)
-
 """
-    contains_infinity([F::Flavor=default_flavor()], x::BareInterval)
+    contains_infinity([F::Flavor,] x::BareInterval)
 
 For the given flavor `F`, test whether `x` contains infinity.
 """
 contains_infinity(::Flavor{:set_based}, ::BareInterval) = false
 
-contains_infinity(x::BareInterval) = contains_infinity(default_flavor(), x)
-
 """
-    is_valid_interval([F::Flavor=default_flavor()], a::Real, b::Real)
+    is_valid_interval([F::Flavor,] a::Real, b::Real)
 
 For the given flavor `F`, test whether ``[a, b]`` is a valid interval.
 """
 is_valid_interval(::Flavor{:set_based}, a::Real, b::Real) = b - a ≥ 0
 # to prevent issues with division by zero, e.g. `is_valid_interval(1//0, 1//0)`
 is_valid_interval(::Flavor{:set_based}, a::Rational, b::Rational) =
     !((a > b) | (a == typemax(typeof(a))) | (b == typemin(typeof(b))))
-
-is_valid_interval(a::Real, b::Real) = is_valid_interval(default_flavor(), a, b)