diff --git a/src/intervals/arithmetic/basic.jl b/src/intervals/arithmetic/basic.jl index 7f4887541..5b056b4c4 100644 --- a/src/intervals/arithmetic/basic.jl +++ b/src/intervals/arithmetic/basic.jl @@ -8,17 +8,9 @@ +(a::Interval) = a # Not in the IEEE standard -""" - -(a::Interval) - -Implement the `neg` function of the IEEE Std 1788-2015 (Table 9.1). -""" --(a::F) where {F<:Interval} = F(-sup(a), -inf(a)) - - """ +(a::Interval, b::Real) - +(a::Real, a::Interval) + +(a::Real, b::Interval) +(a::Interval, b::Interval) Implement the `add` function of the IEEE Std 1788-2015 (Table 9.1). @@ -34,10 +26,18 @@ function +(a::F, b::F) where {F<:Interval} (isempty(a) || isempty(b)) && return emptyinterval(F) return @round(F, inf(a) + inf(b), sup(a) + sup(b)) end ++(a::Interval, b::Interval) = +(promote(a, b)...) + +""" + -(a::Interval) + +Implement the `neg` function of the IEEE Std 1788-2015 (Table 9.1). +""" +-(a::F) where {F<:Interval} = F(-sup(a), -inf(a)) """ -(a::Interval, b::Real) - -(a::Real, a::Interval) + -(a::Real, b::Interval) -(a::Interval, b::Interval) Implement the `sub` function of the IEEE Std 1788-2015 (Table 9.1). @@ -46,19 +46,18 @@ function -(a::F, b::T) where {T<:Real, F<:Interval{T}} isempty(a) && return emptyinterval(F) return @round(F, inf(a) - b, sup(a) - b) end - function -(b::T, a::F) where {T, F<:Interval{T}} isempty(a) && return emptyinterval(F) return @round(F, b - sup(a), b - inf(a)) end +-(a::F, b::Real) where {F<:Interval} = a - F(b) +-(a::Real, b::F) where {F<:Interval} = F(a) - b function -(a::F, b::F) where {F<:Interval} (isempty(a) || isempty(b)) && return emptyinterval(F) return @round(F, inf(a) - sup(b), sup(a) - inf(b)) end - --(a::F, b::Real) where {F<:Interval} = a - F(b) --(a::Real, b::F) where {F<:Interval} = F(a) - b +-(a::Interval, b::Interval) = -(promote(a, b)...) """ scale(α, a::Interval) @@ -71,7 +70,7 @@ For efficiency, does not check that the constant is positive. """ *(a::Interval, b::Real) - *(a::Real, a::Interval) + *(a::Real, b::Interval) *(a::Interval, b::Interval) Implement the `mul` function of the IEEE Std 1788-2015 (Table 9.1). @@ -88,25 +87,21 @@ function *(x::T, a::F) where {T<:Real, F<:Interval{T}} return @round(F, sup(a)*x, inf(a)*x) end end - *(x::T, a::F) where {T<:Real, S, F<:Interval{S}} = Interval{S}(x) * a *(a::F, x::T) where {T<:Real, S, F<:Interval{S}} = x*a function *(a::F, b::F) where {F<:Interval} (isempty(a) || isempty(b)) && return emptyinterval(F) - (isthinzero(a) || isthinzero(b)) && return zero(F) - (isbounded(a) && isbounded(b)) && return mult(*, a, b) - return mult((x, y, r) -> unbounded_mult(F, x, y, r), a, b) end - +*(a::Interval, b::Interval) = *(promote(a, b)...) # Helper functions for multiplication function unbounded_mult(::Type{F}, x::T, y::T, r::RoundingMode) where {T, F<:Interval{T}} - iszero(x) && return sign(y)*zero_times_infinity(T) - iszero(y) && return sign(x)*zero_times_infinity(T) + iszero(x) && return sign(y) * zero_times_infinity(T) + iszero(y) && return sign(x) * zero_times_infinity(T) return *(x, y, r) end @@ -114,11 +109,11 @@ function mult(op, a::F, b::F) where {T, F<:Interval{T}} if inf(b) >= zero(T) inf(a) >= zero(T) && return @round(F, op(inf(a), inf(b)), op(sup(a), sup(b))) sup(a) <= zero(T) && return @round(F, op(inf(a), sup(b)), op(sup(a), inf(b))) - return @round(F, inf(a)*sup(b), sup(a)*sup(b)) # when zero(T) ∈ a + return @round(F, inf(a)*sup(b), sup(a)*sup(b)) # zero(T) ∈ a elseif sup(b) <= zero(T) inf(a) >= zero(T) && return @round(F, op(sup(a), inf(b)), op(inf(a), sup(b))) sup(a) <= zero(T) && return @round(F, op(sup(a), sup(b)), op(inf(a), inf(b))) - return @round(F, sup(a)*inf(b), inf(a)*inf(b)) # when zero(T) ∈ a + return @round(F, sup(a)*inf(b), inf(a)*inf(b)) # zero(T) ∈ a else inf(a) > zero(T) && return @round(F, op(sup(a), inf(b)), op(sup(a), sup(b))) sup(a) < zero(T) && return @round(F, op(inf(a), sup(b)), op(inf(a), inf(b))) @@ -129,7 +124,7 @@ end """ /(a::Interval, b::Real) - /(a::Real, a::Interval) + /(a::Real, b::Interval) /(a::Interval, b::Interval) Implement the `div` function of the IEEE Std 1788-2015 (Table 9.1). @@ -182,6 +177,7 @@ function /(a::F, b::F) where {T, F<:Interval{T}} end end end +/(a::Interval, b::Interval) = /(promote(a, b)...) """ inv(a::Interval) diff --git a/src/intervals/construction.jl b/src/intervals/construction.jl index 51fdd2739..bf5f36404 100644 --- a/src/intervals/construction.jl +++ b/src/intervals/construction.jl @@ -115,6 +115,11 @@ Interval(x::Irrational) = Interval{default_bound()}(x) return :(return $res) # Set body of the function to return the precomputed result end +# promotion + +Base.promote_rule(::Type{Interval{T}}, ::Type{Interval{S}}) where {T,S} = + Interval{promote_type(T, S)} + """ interval(a, b) @@ -123,15 +128,13 @@ If so, then an `Interval(a, b)` object is returned; if not, a warning is printed and the empty interval is returned. """ function interval(a::T, b::S) where {T<:Real, S<:Real} - if !is_valid_interval(a, b) - @warn "Invalid input, empty interval is returned" - return emptyinterval(promote_type(T, S)) - end - - return Interval(a, b) + is_valid_interval(a, b) && return Interval(a, b) + @warn "Invalid input, empty interval is returned" + return emptyinterval(promote_type(T, S)) end interval(a::Real) = interval(a, a) +interval(a::Interval) = interval(inf(a), sup(a)) # Check the validity of the interval const checked_interval = interval @@ -240,4 +243,4 @@ function bigequiv(x::AbstractFloat) end float(x::Interval{T}) where T = atomic(Interval{float(T)}, x) -big(x::Interval) = atomic(Interval{BigFloat}, x) \ No newline at end of file +big(x::Interval) = atomic(Interval{BigFloat}, x) diff --git a/test/interval_tests/construction.jl b/test/interval_tests/construction.jl index d8c9a6e91..32d46979c 100644 --- a/test/interval_tests/construction.jl +++ b/test/interval_tests/construction.jl @@ -39,6 +39,9 @@ using Test @test big(ℯ) in Interval{Float32}(0, ℯ) @test big(π) in Interval{Float32}(π, 4) + @test interval(Interval(pi)) ≛ Interval(pi) + @test interval(Interval(NaN, -Inf)) ≛ emptyinterval() + # a < Inf and b > -Inf @test @interval("1e300") ≛ Interval(9.999999999999999e299, 1.0e300) @test @interval("-1e307") ≛ Interval(-1.0000000000000001e307, -1.0e307) @@ -164,6 +167,16 @@ end @test convert(Interval{BigFloat}, x) === x end +@testset "Promotion between intervals" begin + x = Interval{Float64}(π) + y = Interval{BigFloat}(π) + x_, y_ = promote(x, y) + + @test promote_type(typeof(x), typeof(y)) == Interval{BigFloat} + @test bounds(x_) == (BigFloat(inf(x), RoundDown), BigFloat(sup(x), RoundUp)) + @test y_ ≛ y +end + @testset "Typed intervals" begin @test typeof(@interval Float64 1 2) == Interval{Float64} @test typeof(@interval 1 2) == Interval{Float64} @@ -185,6 +198,11 @@ end a = convert(Interval{Float64}, @biginterval(3, 4)) @test typeof(a) == Interval{Float64} + + pi64, pi32 = Interval{Float64}(pi), Interval{Float32}(pi) + x, y = promote(pi64, pi32) + @test x ≛ pi64 + @test y ≛ Interval{Float64}(pi32) end @testset "Interval{T} constructor" begin diff --git a/test/interval_tests/numeric.jl b/test/interval_tests/numeric.jl index 22c5a0d4c..2d43835d7 100644 --- a/test/interval_tests/numeric.jl +++ b/test/interval_tests/numeric.jl @@ -40,6 +40,12 @@ end @test a + b ≛ Interval(+(a.lo, b.lo, RoundDown), +(a.hi, b.hi, RoundUp)) @test -a ≛ Interval(-a.hi, -a.lo) @test a - b ≛ Interval(-(a.lo, b.hi, RoundDown), -(a.hi, b.lo, RoundUp)) + for f in (:+, :-, :*, :/) + @eval begin + @test $f(Interval{Float64}(pi), Interval{Float32}(pi)) ≛ + $f(Interval{Float64}(pi), Interval{Float64}(Interval{Float32}(pi))) + end + end @test Interval(1//4,1//2) + Interval(2//3) ≛ Interval(11//12, 7//6) @test_broken Interval(1//4,1//2) - Interval(2//3) ≛ Interval(-5//12, -1//6)