Fixes in power with rationals and nthroot

JuliaIntervals · May 28, 2022 · 260b4eb · 260b4eb
1 parent b19f757
commit 260b4eb
Showing 1 changed file with 37 additions and 56 deletions.
diff --git a/src/intervals/arithmetic/power.jl b/src/intervals/arithmetic/power.jl
@@ -130,6 +130,8 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer}
 
     isempty(a) && return emptyinterval(a)
     iszero(x) && return one(a)
+    # x < 0 && return inv(a^(-x))
+    x < 0 && return F( inv(bigequiv(a^(-x)) ) )
 
     if isthinzero(a)
         x > zero(x) && return zero(a)
@@ -140,44 +142,26 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer}
 
     x == (1//2) && return sqrt(a)
 
-    alo = inf(a)
-    ahi = sup(a)
-
-    if x >= 0
-        if alo ≥ 0
-            abig = @biginterval(a)
-            ui = convert(Culong, q)
-            low = BigFloat()
-            high = BigFloat()
-            ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , inf(abig) , ui, MPFRRoundDown)
-            ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , sup(abig) , ui, MPFRRoundUp)
-            b = interval(low, high)
-            b = convert(F, b)
-            return b^p
-        end
-
-        if alo < 0 && ahi ≥ 0
-            a = a ∩ F(0, Inf)
-            abig = @biginterval(a)
-            ui = convert(Culong, q)
-            low = BigFloat()
-            high = BigFloat()
-            ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , inf(abig) , ui, MPFRRoundDown)
-            ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , sup(abig) , ui, MPFRRoundUp)
-            b = interval(low, high)
-            b = convert(F, b)
-            return b^p
-        end
-
-        if ahi < 0
-            return emptyinterval(a)
-        end
+    alo, ahi = bounds(a)
 
+    if ahi < 0
+        return emptyinterval(a)
     end
 
-    if x < 0
-        return inv(a^(-x))
+    if alo < 0 && ahi ≥ 0
+        a = a ∩ F(0, Inf)
     end
+
+    abig = bigequiv(a)
+    abiglo, abighi = bounds(abig)
+    ui = convert(Culong, q)
+    low = BigFloat()
+    high = BigFloat()
+    ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , abiglo , ui, MPFRRoundDown)
+    ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , abighi , ui, MPFRRoundUp)
+    b = interval(low, high) ^ p
+    return F(b)
+
 end
 
 # Interval power of an interval:
@@ -297,35 +281,32 @@ end
 Compute the real n-th root of Interval.
 """
 function nthroot(a::Interval{BigFloat}, n::Integer)
+    isempty(a) && return a
     n == 1 && return a
     n == 2 && return sqrt(a)
-    n < 0 && isthinzero(a) && return emptyinterval(a)
-    isempty(a) && return a
-    if n > 0
-        alo = inf(a)
-        ahi = sup(a)
-        ahi < 0 && iseven(n) && return emptyinterval(BigFloat)
-        if alo < 0 && ahi >= 0 && iseven(n)
-            a = a ∩ Interval{BigFloat}(0, Inf)
-        end
-        ui = convert(Culong, n)
-        low = BigFloat()
-        high = BigFloat()
-        ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , a.lo , ui, MPFRRoundDown)
-        ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , a.hi , ui, MPFRRoundUp)
-        b = interval(low , high)
-        return b
-    elseif n < 0
-        return inv(nthroot(a, -n))
-    elseif n == 0
-        return emptyinterval(a)
+    n == 0 && return emptyinterval(a)
+    # n < 0 && isthinzero(a) && return emptyinterval(a)
+    n < 0 && return inv(nthroot(a, -n))
+
+    alo, ahi = bounds(a)
+    ahi < 0 && iseven(n) && return emptyinterval(BigFloat)
+    if alo < 0 && ahi >= 0 && iseven(n)
+        a = a ∩ Interval{BigFloat}(0, Inf)
+        alo, ahi = bounds(a)
     end
+    ui = convert(Culong, n)
+    low = BigFloat()
+    high = BigFloat()
+    ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , alo , ui, MPFRRoundDown)
+    ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , ahi , ui, MPFRRoundUp)
+    return interval(low , high)
 end
 
-function nthroot(a::Interval{T}, n::Integer) where T
+function nthroot(a::F, n::Integer) where {T, F<:Interval{T}}
     n == 1 && return a
     n == 2 && return sqrt(a)
+    n < 0 && return inv(nthroot(a, -n))
     abig = @biginterval(a)
     b = nthroot(abig, n)
-    return convert(Interval{T}, b)
+    return F(b)
 end