Fix rational power function[WIP] (#286)

* Fix rational power function

* Modify rational power function

* Removing Big interval

* Add rounding

* Add rational power function

* improve rational power function

* Using ccall function

* wrapping mpfr function

* add tests

* add version

* Add project.toml file

* Remove older julia versions

* correcting in julia versions

Author: yashrajgupta

.travis.yml

@@ -10,8 +10,6 @@ os:
   - osx
 
 julia:
-  - 0.7
-  - 1.0
   - 1.1
   - nightly
 

Project.toml

@@ -0,0 +1,23 @@
+name = "IntervalArithmetic"
+uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+
+[compat]
+CRlibm = "≥ 0.7.0"
+StaticArrays = "≥ 0.8.0"
+FastRounding = "≥ 0.1.2"
+SetRounding = "≥ 0.2.0"
+RecipesBase = "≥ 0.5.0"
+julia = "≥ 1.1.0"
+
+[deps]
+CRlibm = "96374032-68de-5a5b-8d9e-752f78720389"
+FastRounding = "fa42c844-2597-5d31-933b-ebd51ab2693f"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SetRounding = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]

REQUIRE

@@ -1,4 +1,4 @@
-julia 0.7
+julia 1.1
 CRlibm 0.7
 StaticArrays 0.8
 FastRounding 0.1.2

appveyor.yml

@@ -1,7 +1,6 @@
 environment:
   matrix:
-  - julia_version: 0.7
-  - julia_version: 1
+  - julia_version: 1.1
   - julia_version: nightly
 
 platform:

src/IntervalArithmetic.jl

@@ -43,6 +43,9 @@ import Base:  # for IntervalBox
     getindex, setindex,
     iterate, eltype
 
+import Base.MPFR: MPFRRoundingMode
+import Base.MPFR: MPFRRoundUp, MPFRRoundDown, MPFRRoundNearest, MPFRRoundToZero, MPFRRoundFromZero
+
 import .Broadcast: broadcasted
 
 export

src/intervals/functions.jl

@@ -6,7 +6,7 @@
 # Use the BigFloat version from MPFR instead, which is correctly-rounded:
 
 # Write explicitly like this to avoid ambiguity warnings:
-for T in (:Integer, :Rational, :Float64, :BigFloat, :Interval)
+for T in (:Integer, :Float64, :BigFloat, :Interval)
     @eval ^(a::Interval{Float64}, x::$T) = atomic(Interval{Float64}, big53(a)^x)
 end
 
@@ -142,25 +142,61 @@ function ^(a::Interval{Rational{T}}, x::AbstractFloat) where T<:Integer
 end
 
 # Rational power
-function ^(a::Interval{BigFloat}, r::Rational{S}) where S<:Integer
-    T = BigFloat
+function ^(a::Interval{T}, x::Rational) where T
     domain = Interval{T}(0, Inf)
 
+    p = x.num
+    q = x.den
+
+    isempty(a) && return emptyinterval(a)
+    x == 0 && return one(a)
+
     if a == zero(a)
-        a = a ∩ domain
-        r > zero(r) && return zero(a)
+        x > zero(x) && return zero(a)
         return emptyinterval(a)
     end
 
-    isinteger(r) && return atomic(Interval{T}, a^round(S,r))
-    r == one(S)//2 && return sqrt(a)
+    x == (1//2) && return sqrt(a)
+
+    if x >= 0
+        if a.lo ≥ 0
+            isinteger(x) && return a ^ Int64(x)
+            a = @biginterval(a)
+            ui = convert(Culong, q)
+            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)
+            b = convert(Interval{T}, b)
+            return b^p
+        end
 
-    a = a ∩ domain
-    (isempty(r) || isempty(a)) && return emptyinterval(a)
+        if a.lo < 0 && a.hi ≥ 0
+            isinteger(x) && return a ^ Int64(x)
+            a = a ∩ Interval{T}(0, Inf)
+            a = @biginterval(a)
+            ui = convert(Culong, q)
+            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)
+            b = convert(Interval{T}, b)
+            return b^p
+        end
+
+        if a.hi < 0
+            isinteger(x) && return a ^ Int64(x)
+            return emptyinterval(a)
+        end
 
-    y = atomic(Interval{BigFloat}, r)
+    end
+
+    if x < 0
+        return inv(a^(-x))
+    end
 
-    a^y
 end
 
 # Interval power of an interval:

test/interval_tests/numeric.jl

@@ -114,9 +114,15 @@ end
     setprecision(Interval, 256)
     x = @biginterval(27)
     y = x^(1//3)
-    @test (0 < diam(y) < 1e-76)
+    @test diam(y) == 0
     y = x^(1/3)
-    @test (0 < diam(y) < 1e-76)
+    @test (0 <= diam(y) < 1e-76)
+    x = @biginterval(9.595703125)
+    y = x^(1//3)
+    @test diam(y) == 0
+    x = @biginterval(0.1)
+    y = x^(1//3)
+    @test (0 <= diam(y) < 1e-76)
 
 end
 

test/interval_tests/power.jl

@@ -0,0 +1,60 @@
+#Test library imports
+using Test
+
+#Arithmetic library imports
+using IntervalArithmetic
+
+#Preamble
+setprecision(53)
+setprecision(Interval, Float64)
+setrounding(Interval, :tight)
+# Set full format, and show decorations
+@format full
+@testset "rational_power_test" begin
+        @test ^(∅, 1//3) == ∅
+        @test ^(1 .. 8, 1//3) == Interval(1, 2)
+        @test ^(2 .. 8, 1//3) ⊇  Interval(2^(1//3), 2)
+        @test ^(1 .. 9, 1//3) ⊇  Interval(1, 9^(1//3))
+        @test ^(2 .. 9, 1//3) ⊇  Interval(2^(1//3), 9^(1//3))
+        @test ^(-1 .. 8, 1//3) == Interval(0, 2)
+        @test ^(-2 .. 8, 1//3) ⊇  Interval(0, 2)
+        @test ^(-1 .. 9, 1//3) ⊇  Interval(0, 9^(1//3))
+        @test ^(-2 .. 9, 1//3) ⊇  Interval(0, 9^(1//3))
+        @test ^(1 .. 8, -1//3) == Interval(0.5, 1)
+        @test ^(2 .. 8, -1//3) ⊇  Interval(0.5, 2^(-1//3))
+        @test ^(1 .. 9, -1//3) ⊇  Interval(9^(-1//3), 1)
+        @test ^(2 .. 9, -1//3) ⊇  Interval(9^(-1//3), 2^(-1//3))
+        @test ^(-1 .. 8, -1//3) == Interval(0.5, Inf)
+        @test ^(-2 .. 8, -1//3) ⊇  Interval(0.5, Inf)
+        @test ^(-1 .. 9, -1//3) ⊇  Interval(9^(-1//3), Inf)
+        @test ^(-2 .. 9, -1//3) ⊇  Interval(9^(-1//3), Inf)
+        @test ^(-2 .. 4 , 1//2) == Interval(0, 2)
+        @test ^(-2 .. 8 , 1//3) == Interval(0, 2)
+        @test ^(-8 .. -2 , 1//3) == ∅
+        @test ^(-8 .. -2 , 1//2) == ∅
+        @test ^(-8 .. -2 , -1//3) == ∅
+        @test ^(-8 .. -2 , -1//2) == ∅
+        @test ^(∅, 2//3) == ∅
+        @test ^(1 .. 8, 2//3) == Interval(1, 4)
+        @test ^(2 .. 8, 2//3) ⊇  Interval(2^(2//3), 4)
+        @test ^(1 .. 9, 2//3) ⊇  Interval(1, 9^(2//3))
+        @test ^(2 .. 9, 2//3) ⊇  Interval(2^(2//3), 9^(2//3))
+        @test ^(-1 .. 8, 2//3) == Interval(0, 4)
+        @test ^(-2 .. 8, 2//3) ⊇  Interval(0, 4)
+        @test ^(-1 .. 9, 2//3) ⊇  Interval(0, 9^(2//3))
+        @test ^(-2 .. 9, 2//3) ⊇  Interval(0, 9^(2//3))
+        @test ^(1 .. 8, -2//3) == Interval(0.25, 1)
+        @test ^(2 .. 8, -2//3) ⊇  Interval(0.25, 2^(-2//3))
+        @test ^(1 .. 9, -2//3) ⊇  Interval(9^(-2//3), 1)
+        @test ^(2 .. 9, -2//3) ⊇  Interval(9^(-2//3), 2^(-2//3))
+        @test ^(-1 .. 8, -2//3) == Interval(0.25, Inf)
+        @test ^(-2 .. 8, -2//3) ⊇  Interval(0.25, Inf)
+        @test ^(-1 .. 9, -2//3) ⊇  Interval(9^(-2//3), Inf)
+        @test ^(-2 .. 9, -2//3) ⊇  Interval(9^(-2//3), Inf)
+        @test ^(-2 .. 4 , 3//2) == Interval(0, 8)
+        @test ^(-2 .. 8 , 2//3) == Interval(0, 4)
+        @test ^(-8 .. -2 , 2//3) == ∅
+        @test ^(-8 .. -2 , 3//2) == ∅
+        @test ^(-8 .. -2 , -2//3) == ∅
+        @test ^(-8 .. -2 , -3//2) == ∅
+end