fixed bugs towards ITF1788 compliance

.github/workflows/CI.yml

@@ -11,8 +11,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
-          - '1.6'
+          - '1'
           - 'nightly'
         os:
           - ubuntu-latest

.gitignore

@@ -1,3 +1,4 @@
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+Manifest.toml

Project.toml

@@ -1,12 +1,12 @@
 name = "IntervalContractors"
 uuid = "15111844-de3b-5229-b4ba-526f2f385dc9"
-version = "0.4.3"
+version = "0.4.6"
 
 [deps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 
 [compat]
-IntervalArithmetic = "0.16, 0.17, 0.18"
+IntervalArithmetic = "0.16, 0.17, 0.18, 0.19, 0.20"
 julia = "1.3, 1.4, 1.5"
 
 [extras]

src/IntervalContractors.jl

@@ -12,7 +12,8 @@ export plus_rev, minus_rev, inv_rev,
         asin_rev, acos_rev, atan_rev,
         sinh_rev, cosh_rev, tanh_rev,
         asinh_rev, acosh_rev, atanh_rev,
-        mul_rev_IEEE1788
+        mul_rev_IEEE1788, mul_rev_to_pair,
+        pow_rev1, pow_rev2
 
 using IntervalArithmetic
 
@@ -32,6 +33,7 @@ include("inverse_trig.jl")
 include("hyperbolic.jl")
 include("inverse_hyperbolic.jl")
 include("extrema.jl")
+include("decorated.jl")
 
 """
 Dictionary mapping functions to their reverse functions.

src/arithmetic.jl

@@ -1,6 +1,16 @@
 
 """
-Reverse plus
+    plus_rev(a::Interval, b::Interval[, c::Interval])
+
+Reverse addition. Calculates the preimage of `a = b + c` for `b` and `c`.
+
+### Output
+
+The triplet `(a, b_new, c_new)` where
+
+- `a` remains unchanged
+- `b_new` is the interval hull of the set ``{x ∈ b : ∃ y ∈ c, x + y ∈ a}``
+- `c_new` is the interval hull of the set ``{y ∈ c : ∃ x ∈ b, x + y ∈ a}``
 """
 function plus_rev(a::Interval, b::Interval, c::Interval)  # a = b + c
     # a = a ∩ (b + c)  # add this line for plus contractor (as opposed to reverse function)
@@ -13,7 +23,17 @@ end
 plus_rev(a,b,c) = plus_rev(promote(a,b,c)...)
 
 """
-Reverse minus
+    minus_rev(a::Interval, b::Interval[, c::Interval])
+
+Reverse subtraction. Calculates the preimage of `a = b - c` for `b` and `c`.
+
+### Output
+
+The triplet `(a, b_new, c_new)` where
+
+- `a` remains unchanged
+- `b_new` is the interval hull of the set ``{x ∈ b : ∃ y ∈ c, x - y ∈ a}``
+- `c_new` is the interval hull of the set ``{y ∈ c : ∃ x ∈ b, x - y ∈ a}``
 """
 function minus_rev(a::Interval, b::Interval, c::Interval)  # a = b - c
 
@@ -76,7 +96,16 @@ end
 div_rev(a,b,c) = div_rev(promote(a,b,c)...)
 
 """
-Reverse inverse
+    inv_rev(a::Interval, b::Interval)
+
+Reverse inverse. Calculates the interval hull of the preimage of a = b⁻¹
+
+### Output
+
+Pair `(a, b_new)` where
+
+- `a` is unchanged
+- `b_new` is the interval hull of the set ``{x ∈ b : x⁻¹ ∈ a}``
 """
 function inv_rev(a::Interval, b::Interval)  # a = inv(b)
 
@@ -88,9 +117,24 @@ end
 inv_rev(a,b) = inv_rev(promote(a,b)...)
 
 """
-Reverse power
+    power_rev(a::Interval, b::Interval, n::Integer)
+
+Reverse power. Calculates the preimage of `a = bⁿ`.  See section 10.5.4 of the
+IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+The triplet `(a, b_new, n)` where
+
+- `a` and `n` are unchanged
+- `b_new` is the interval hull of the set ``{x ∈ b : xⁿ ∈ a}``
 """
-function power_rev(a::Interval, b::Interval, n::Integer)  # a = b^n,  log(a) = n.log(b),  b = a^(1/n)
+function power_rev(a::Interval{T}, b::Interval{T}, n::Integer) where T  # a = b^n,  log(a) = n.log(b),  b = a^(1/n)
+
+    if iszero(n)
+        1 ∈ a && return (a, entireinterval(T) ∩ b, n)
+        return (a, emptyinterval(T), n)
+    end
 
     if n == 2  # a = b^2
         root = √a
@@ -117,6 +161,7 @@ function power_rev(a::Interval, b::Interval, n::Integer)  # a = b^n,  log(a) = n
     return (a, b, n)
 end
 
+power_rev(a::Interval{T}, n::Integer) where T = power_rev(a, entireinterval(T), n)
 
 function power_rev(a::Interval, b::Interval, c::Interval)  # a = b^c
 
@@ -135,7 +180,16 @@ power_rev(a, b, c) = power_rev(promote(a, b, c)...)
 
 
 """
-Reverse square root
+    sqrt_rev(a::Interval, b::Interval)
+
+Reverse square root. Calculates the preimage of `a = √b`.
+
+### Output
+
+The pair `(a, b_new)` where
+
+- `a` is unchanged
+- `b_new` is the interval hull of the set ``{x ∈ b : √x ∈ a}``
 """
 function sqrt_rev(a::Interval, b::Interval)  # a = sqrt(b)
 
@@ -150,7 +204,17 @@ sqrt_rev(a,b) = sqrt_rev(promote(a,b)...)
 # IEEE-1788 style
 
 """
-Reverse sqr
+    sqrt_rev(c::Interval[, x::Interval])
+
+Reverse square. Calculates the preimage of `a = x²`. If `x` is not provided, then
+byt default ``[-∞, ∞]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+The pair `(c, x_new)` where
+
+- `c` is unchanged
+- `x_new` is the interval hull of the set ``{x ∈ b : x² ∈ a}``
 """
 function sqr_rev(c, x)   # c = x^2;  refine x
 
@@ -162,10 +226,18 @@ function sqr_rev(c, x)   # c = x^2;  refine x
     return (c, hull(x1, x2))
 end
 
-sqr_rev(c) = sqr_rev(c, -∞..∞)
-
 """
-Reverse abs
+    abs_rev(c::Interval[, x::Interval])
+
+Reverse absolute value. Calculates the preimage of `a = |x|`. If `x` is not provided, then
+byt default ``[-∞, ∞]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+The pair `(c, x_new)` where
+
+- `c` is unchanged
+- `x_new` is the interval hull of the set ``{x ∈ b : |x| ∈ a}``
 """
 function abs_rev(y, x)   # y = abs(x); refine x
 
@@ -193,22 +265,83 @@ sign_rev(a,b) = sign_rev(promote(a,b)...)
 ## IEEE-1788 versions:
 
 """
-According to the IEEE-1788 standard:
+    mul_rev_IEEE1788(b::Interval, c::Interval[, x::Interval])
+
+Reverse multiplication. Computes the preimage of ``c=x * b`` with respect to `x`. If `x` is not provided,
+then byt default ``[-∞, ∞]`` is used.. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic.
 
-- `∘_rev1(b, c, x)` is the subset of `x` such that `x ∘ b` is defined and in `c`;
-- `∘_rev2(a, c, x)` is the subset of `x` such that `a ∘ x` is defined and in `c`
+### Output
 
-When `∘` is commutative, these agree and we write `∘_rev(b, c, x)`.
+- `x_new` the interval hull of the set ``{t ∈ x : ∃ y ∈ b, t*y ∈ c}
 """
+mul_rev_IEEE1788(b, c, x) = mul_rev(c, x, b)[2]
 
-function mul_rev_IEEE1788(b, c, x)   # c = b*x
-    return x ∩ (c / b)
-end
+"""
+    pow_rev1(b::Interval, c::Interval[, x::Interval])
 
+Reverse power 1. Computes the preimage of ``c=xᵇ`` with respect to `x`. If `x` is not provided,
+then byt default ``[-∞, ∞]`` is used.. See section 10.5.4 of the
+IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+- `x_new` the interval hull of the set ``{t ∈ x : ∃ y ∈ b, tʸ ∈ c}
+"""
 function pow_rev1(b, c, x)   # c = x^b
     return x ∩ c^(1/b)  # replace by 1//b
 end
 
+"""
+    pow_rev2(b::Interval, c::Interval[, x::Interval])
+
+Reverse power 2. Computes the preimage of ``c = aˣ`` with respect to `x`. If `x` is not provided, then
+byt default ``[-∞, ∞]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+- `x_new` the interval hull of the set ``{t ∈ x : ∃ y ∈ b, tʸ ∈ c}
+"""
 function pow_rev2(a, c, x)   # c = a^x
-    return x ∩ (log(c) / lob(a))
+    return x ∩ (log(c) / log(a))
+end
+
+"""
+    mul_rev_to_pair(b::Interval, c::Interval)
+
+Computes the division c/b, but returns a pair of intervals instead of a single interval.
+If the set corresponding to c/b is composed by two disjoint intervals, then it returns the
+two intervals. If c/b is a single or empty interval, then the second interval in the pair
+is set to empty. See section 10.5.5 of the IEEE 1788-2015 standard for interval arithmetic.
+
+### Example
+
+```jldoctest
+julia> mul_rev_to_pair(-1..1, 1..2)
+([-∞, -1], [1, ∞])
+
+julia> mul_rev_to_pair(1..2, 3..4)
+([1.5, 4], ∅)
+```
+
+"""
+mul_rev_to_pair(b::Interval, c::Interval) = extended_div(c, b)
+
+function mul_rev_to_pair(b::DecoratedInterval{T}, c::DecoratedInterval{T}) where T
+    (isnai(b) || isnai(c)) && return (nai(T), nai(T))
+
+    0 ∉ b && return (c/b, DecoratedInterval(emptyinterval(T), trv))
+
+    x1, x2 = extended_div(interval(c), interval(b))
+    return (DecoratedInterval(x1, trv), DecoratedInterval(x2, trv))
+end
+
+mul_rev_to_pair(b::Interval, c::Interval) = extended_div(c, b)
+
+function mul_rev_to_pair(b::DecoratedInterval{T}, c::DecoratedInterval{T}) where T
+    (isnai(b) || isnai(c)) && return (nai(T), nai(T))
+
+    0 ∉ b && return (c/b, DecoratedInterval(emptyinterval(T), trv))
+
+    x1, x2 = extended_div(interval(c), interval(b))
+    return (DecoratedInterval(x1, trv), DecoratedInterval(x2, trv))
 end

src/decorated.jl

@@ -0,0 +1,35 @@
+entiredecorated(T) = DecoratedInterval(entireinterval(T))
+
+for op in (:sqr_rev, :abs_rev, :sin_rev, :cos_rev, :tan_rev, :cosh_rev, :sinh_rev, :tanh_rev)
+    @eval begin 
+        function $op(a::DecoratedInterval{T}, x::DecoratedInterval{T}) where T 
+            ( isnai(a) || isnai(x) ) && return nai(T)
+            bare = $op(interval(a), interval(x))
+            return (DecoratedInterval(bare[1], trv), DecoratedInterval(bare[2], trv))
+        end
+    end 
+    @eval $op(a::Interval{T}) where T = $op(a, entireinterval(T))
+    @eval $op(a::DecoratedInterval{T}) where T = $op(a, entiredecorated(T)) 
+end
+
+function power_rev(a::DecoratedInterval{T}, x::DecoratedInterval{T}, n::Integer) where T
+    ( isnai(a) || isnai(x) ) && return nai(T)
+    bare = power_rev(interval(a), interval(x), n)
+    return (DecoratedInterval(bare[1], trv), DecoratedInterval(bare[2], trv), n)
+end
+
+power_rev(a::DecoratedInterval{T}, n::Integer) where T = power_rev(a, entiredecorated(T), n)
+
+for op in (:mul_rev_IEEE1788, :pow_rev1, :pow_rev2)
+    @eval begin
+        function $op(b::DecoratedInterval{T}, c::DecoratedInterval{T}, x::DecoratedInterval{T}) where T
+            (isnai(b) || isnai(c) || isnai(x) ) && return nai(T)
+            bare = $op(interval(b), interval(c), interval(x))
+            return DecoratedInterval(bare, trv)
+        end
+    end
+
+    @eval $op(a::Interval{T}, b::Interval{T}) where T = $op(a, b, entireinterval(T))
+    @eval $op(a::DecoratedInterval{T}, b::DecoratedInterval{T}) where T = $op(a, b, entiredecorated(T))
+
+end

src/hyperbolic.jl

@@ -1,5 +1,15 @@
 """
-Reverse function for `sinh`.
+    sinh_rev(c::Interval[, x::Interval])
+
+Reverse hyperbolic sine. Calculates the preimage of `a = sinh(x)`. If `x` is not provided, then
+byt default ``[-∞, ∞]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+The pair `(c, x_new)` where
+
+- `c` is unchanged
+- `x_new` is the interval hull of the set ``{x ∈ b : sinh(x) ∈ a}``
 """
 function sinh_rev(y::Interval, x::Interval)
     x = x ∩ asinh(y)
@@ -8,17 +18,37 @@ function sinh_rev(y::Interval, x::Interval)
 end
 
 """
-Reverse function for `cosh`.
+    cosh_rev(c::Interval[, x::Interval])
+
+Reverse square root. Calculates the preimage of `a = cosh(x)`. If `x` is not provided, then
+byt default ``[-∞, ∞]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+The pair `(c, x_new)` where
+
+- `c` is unchanged
+- `x_new` is the interval hull of the set ``{x ∈ b : cosh(x) ∈ a}``
 """
 function cosh_rev(y::Interval,x::Interval)
     y_new = y ∩ Interval(1.,∞)
-    x = x ∩ acosh(y)
+    x = (x ∩ acosh(y)) ∪ (x ∩ -acosh(y))
 
     return y_new, x
 end
 
 """
-Reverse function for `tanh`.
+    tanh_rev(c::Interval[, x::Interval])
+
+Reverse square root. Calculates the preimage of `a = tanh(x)`. If `x` is not provided, then
+byt default ``[-∞, ∞]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic.
+
+### Output
+
+The pair `(c, x_new)` where
+
+- `c` is unchanged
+- `x_new` is the interval hull of the set ``{x ∈ b : tanh(x) ∈ a}``
 """
 function tanh_rev(y::Interval,x::Interval)
     y_new = y ∩ Interval(-1.,1.)