Improve support for complex intervals (#636)

Author: OlivierHnt

src/display.jl

@@ -369,7 +369,6 @@ function _round_string_down(s::String)
                     idx == 1 ? "." : "",
                     '9'^(len - (idx + 1)))
             return string(new_mantissa, 'e', exponent)
-
         end
     end
 end

src/intervals/arithmetic/hyperbolic.jl

@@ -25,6 +25,10 @@ for f ∈ (:sinh, :tanh, :asinh)
     end
 end
 
+Base.sinh(x::Complex{Interval{T}}) where {T<:NumTypes} = (exp(x) - exp(-x)) / interval(T, 2)
+
+Base.tanh(x::Complex{<:Interval}) = sinh(x) / cosh(x)
+
 """
     cosh(::BareInterval)
     cosh(::Interval)
@@ -44,6 +48,8 @@ function Base.cosh(x::Interval)
     return _unsafe_interval(r, d, isguaranteed(x))
 end
 
+Base.cosh(x::Complex{Interval{T}}) where {T<:NumTypes} = (exp(x) + exp(-x)) / interval(T, 2)
+
 """
     coth(::BareInterval)
     coth(::Interval)

src/intervals/arithmetic/power.jl

@@ -48,8 +48,8 @@ Interval{Float64}(-Inf, Inf, trv)
 ```
 """
 function Base.:^(x::BareInterval, y::BareInterval)
-    isthininteger(y) && return _select_pown(power_mode(), x, Integer(sup(y)))
-    return _select_pow(power_mode(), x, y)
+    isthininteger(y) || return _select_pow(power_mode(), x, y)
+    return _select_pown(power_mode(), x, Integer(sup(y)))
 end
 
 function Base.:^(x::Interval, y::Interval)
@@ -65,9 +65,17 @@ Base.:^(x::Interval, n::Integer) = ^(x, n//one(n))
 Base.:^(x::Rational, y::Interval) = ^(convert(Interval{typeof(x)}, x), y)
 Base.:^(x::Interval, y::Rational) = ^(x, convert(Interval{typeof(y)}, y))
 
+Base.:^(x::Complex{Interval{T}}, y::Complex{Interval{T}}) where {T<:NumTypes} = exp(y * log(x))
+Base.:^(x::Complex{<:Interval}, y::Complex{<:Interval}) = ^(promote(x, y)...)
+Base.:^(x::Complex{<:Interval}, y::Real) = ^(promote(x, y)...)
+Base.:^(x::Real, y::Complex{<:Interval}) = ^(promote(x, y)...)
+# needed to avoid method ambiguities
+Base.:^(x::Complex{<:Interval}, n::Integer) = ^(promote(x, n)...)
+
 # overwrite behaviour for small integer powers from https://github.com/JuliaLang/julia/pull/24240
 # Base.literal_pow(::typeof(^), x::Interval, ::Val{n}) where {n} = x^n
 Base.literal_pow(::typeof(^), x::Interval, ::Val{n}) where {n} = _select_pown(power_mode(), x, n)
+Base.literal_pow(::typeof(^), x::Complex{Interval{T}}, ::Val{n}) where {T<:NumTypes,n} = exp(interval(T, n) * log(x))
 
 # helper functions for power
 
@@ -413,6 +421,14 @@ for f ∈ (:cbrt, :exp, :exp2, :exp10, :expm1)
     end
 end
 
+Base.exp(x::Complex{<:Interval}) = exp(real(x)) * cis(imag(x))
+
+Base.exp2(x::Complex{<:Interval}) = exp2(real(x)) * cis(imag(x) * log(interval(numtype(x), 2)))
+
+Base.exp10(x::Complex{<:Interval}) = exp10(real(x)) * cis(imag(x) * log(interval(numtype(x), 10)))
+
+Base.expm1(x::Complex{<:Interval}) = exp(x) - interval(numtype(x), 1)
+
 #
 
 for f ∈ (:log, :log2, :log10)
@@ -437,6 +453,12 @@ for f ∈ (:log, :log2, :log10)
     end
 end
 
+Base.log(x::Complex{<:Interval}) = complex(log(abs(x)), angle(x))
+
+Base.log2(x::Complex{<:Interval}) = complex(log2(abs(x)), angle(x)/log(interval(numtype(x), 2)))
+
+Base.log10(x::Complex{<:Interval}) = complex(log10(abs(x)), angle(x)/log(interval(numtype(x), 10)))
+
 function Base.log1p(x::BareInterval{T}) where {T<:AbstractFloat}
     domain = _unsafe_bareinterval(T, -one(T), typemax(T))
     x = intersect_interval(x, domain)
@@ -454,3 +476,5 @@ function Base.log1p(x::Interval{T}) where {T<:NumTypes}
     d = min(d, ifelse(isinterior(bx, domain), d, trv))
     return _unsafe_interval(r, d, isguaranteed(x))
 end
+
+Base.log1p(x::Complex{<:Interval}) = log(interval(numtype(x), 1) + x)

src/intervals/arithmetic/trigonometric.jl

@@ -99,6 +99,8 @@ function Base.sin(x::Interval)
     return _unsafe_interval(r, d, isguaranteed(x))
 end
 
+Base.sin(x::Complex{Interval{T}}) where {T<:NumTypes} = complex(sin(real(x)) * cosh(imag(x)), cos(real(x)) * sinh(imag(x)))
+
 # not in the IEEE Standard 1788-2015
 
 if Int == Int32 && VERSION < v"1.10"
@@ -147,6 +149,9 @@ function Base.sinpi(x::Interval)
     return _unsafe_interval(r, d, isguaranteed(x))
 end
 
+Base.sinpi(x::Complex{Interval{T}}) where {T<:NumTypes} =
+    complex(sinpi(real(x)) * cosh(imag(x) * interval(T, π)), cospi(real(x)) * sinh(imag(x) * interval(T, π)))
+
 """
     cos(::BareInterval)
     cos(::Interval)
@@ -196,6 +201,8 @@ function Base.cos(x::Interval)
     return _unsafe_interval(r, d, isguaranteed(x))
 end
 
+Base.cos(x::Complex{Interval{T}}) where {T<:NumTypes} = complex(cos(real(x)) * cosh(imag(x)), -sin(real(x)) * sinh(imag(x)))
+
 # not in the IEEE Standard 1788-2015
 
 if Int == Int32 && VERSION < v"1.10"
@@ -246,6 +253,9 @@ function Base.cospi(x::Interval)
     return _unsafe_interval(r, d, isguaranteed(x))
 end
 
+Base.cospi(x::Complex{Interval{T}}) where {T<:NumTypes} =
+    complex(cospi(real(x)) * cosh(imag(x) * interval(T, π)), -sinpi(real(x)) * sinh(imag(x) * interval(T, π)))
+
 """
     tan(::BareInterval)
     tan(::Interval)
@@ -282,6 +292,8 @@ function Base.tan(x::Interval)
     return _unsafe_interval(r, d, isguaranteed(x))
 end
 
+Base.tan(x::Complex{<:Interval}) = sin(x) / cos(x)
+
 """
     cot(::BareInterval)
     cot(::Interval)

src/intervals/interval_operations/numeric.jl

@@ -88,7 +88,7 @@ See also: [`inf`](@ref), [`sup`](@ref), [`bounds`](@ref), [`diam`](@ref),
 [`radius`](@ref) and [`midradius`](@ref).
 """
 function mid(x::BareInterval{T}, α = 0.5) where {T<:AbstractFloat}
-    0 ≤ α ≤ 1 || throw(DomainError(α, "α must be between 0 and 1"))
+    0 ≤ α ≤ 1 || return throw(DomainError(α, "α must be between 0 and 1"))
     isempty_interval(x) && return convert(T, NaN)
     if isentire_interval(x)
         α == 0.5 && return zero(T)
@@ -105,7 +105,7 @@ function mid(x::BareInterval{T}, α = 0.5) where {T<:AbstractFloat}
     end
 end
 function mid(x::BareInterval{T}, α = 1//2) where {T<:Rational}
-    0 ≤ α ≤ 1 || throw(DomainError(α, "α must be between 0 and 1"))
+    0 ≤ α ≤ 1 || return throw(DomainError(α, "α must be between 0 and 1"))
     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(T)
@@ -123,12 +123,12 @@ function mid(x::BareInterval{T}, α = 1//2) where {T<:Rational}
 end
 
 function mid(x::Interval{T}, α = 0.5) where {T<:AbstractFloat}
-    0 ≤ α ≤ 1 || throw(DomainError(α, "α must be between 0 and 1"))
+    0 ≤ α ≤ 1 || return throw(DomainError(α, "α must be between 0 and 1"))
     isnai(x) && return convert(T, NaN)
     return mid(bareinterval(x), α)
 end
 function mid(x::Interval{<:Rational}, α = 1//2)
-    0 ≤ α ≤ 1 || throw(DomainError(α, "α must be between 0 and 1"))
+    0 ≤ α ≤ 1 || return throw(DomainError(α, "α must be between 0 and 1"))
     isnai(x) && return throw(ArgumentError("cannot compute the midpoint of an NaI; cannot return a `Rational` NaN"))
     return mid(bareinterval(x), α)
 end

test/interval_tests/consistency.jl

@@ -290,6 +290,20 @@
 
         # Test putting functions in interval:
         @test issubset_interval(log(interval(-2, 5)), interval(-Inf, log(interval(5))))
+
+        #
+
+        x = interval(0, 3) + interval(0, 4)*interval(im)
+        @test isequal_interval(abs2(x), interval(0, 25))
+        @test isequal_interval(abs(x), interval(0, 5))
+
+        y = interval(-1, 1) + interval(-2, 2)*interval(im)
+        @test inf(abs(y)) == 0
+        @test inf(abs2(y)) == 0
+
+        @test mag(x) == 4
+        @test mig(x) == 0
+        @test mid(x) == 1.5 + 2im
     end
 
     @testset "Interval power of an interval" begin

test/interval_tests/power.jl

@@ -84,12 +84,16 @@ end
     @test isequal_interval(pow(interval(-3, 2), interval(-1, 1)), interval(0, Inf, IntervalArithmetic.trv))
 end
 
-# @testset "Complex{<:Interval}" begin
-#     a = interval(3 + 4im)
-#     b = exp(a)
-#     @test real(b) == interval(-13.12878308146216, -13.128783081462153)
-#     @test imag(b) == interval(-15.200784463067956, -15.20078446306795)
-# end
+@testset "Complex{<:Interval}" begin
+    a = interval(3 + 4im)
+    b = exp(a)
+    @test isequal_interval(real(b), interval(-13.12878308146216, -13.128783081462153))
+    @test isequal_interval(imag(b), interval(-15.200784463067956, -15.20078446306795))
+
+    z = exp(-im * interval(π))
+    @test in_interval(-1, real(z))
+    @test in_interval(0, imag(z))
+end
 
 @testset "Literal powers" begin
     x = interval(-1,1)
@@ -100,4 +104,4 @@ end
     @test_broken isguaranteed(x ^ 2305843009213693952)
     @test isequal_interval(x^2, interval(0,1))
     @test isequal_interval(x^3, x)
-end
+end

test/interval_tests/trigonometric.jl

@@ -19,6 +19,11 @@ end
     @test issubset_interval(sin(interval(BigFloat, 0.5, 8.5)), sin(interval(0.5, 8.5)))
     @test issubset_interval(sin(interval(BigFloat, -4.5, 0.1)), sin(interval(-4.5, 0.1)))
     @test issubset_interval(sin(interval(BigFloat, 1.3, 6.3)), sin(interval(1.3, 6.3)))
+
+    #
+
+    z = interval(3, 1e-7; format = :midpoint) + interval(4, 1e-7; format = :midpoint) * interval(im)
+    @test issubset_interval(sin(z), complex(sin(real(z)) * cosh(imag(z)), sinh(imag(z)) * cos(real(z))))
 end
 
 @testset "cos" begin