[WIP / RFC] Configuration mechanism for directed rounding and fast powers

src/IntervalArithmetic.jl

@@ -37,7 +37,8 @@ import Base:
     abs, abs2,
     show,
     isinteger, setdiff,
-    parse, hash
+    parse, hash,
+    power_by_squaring
 
 import Base:  # for IntervalBox
     broadcast, length,
@@ -80,7 +81,7 @@ export
 
 export showfull
 
-import Base: rounding, setrounding, setprecision
+import Base: setprecision
 
 
 
@@ -106,6 +107,8 @@ function __init__()
 
     setprecision(Interval, 256)  # set up pi
     setprecision(Interval, Float64)
+
+    # configure!(directed_rounding=:tight, powers=:fast)
 end
 
 

src/intervals/configuration.jl

@@ -0,0 +1,57 @@
+"""
+Holds configuration information for the package:
+
+- `directed_rounding`:
+
+    Which method to use to implement **directed rounding**.
+    Possible values are:
+
+    - `:tight`: correct rounding using error-free arithmetic and CRlibm. Gives
+        the tightest resulting intervals, with a width of 1 ulp for each basic operation
+    - `:fast`: uses `prevfloat` and `nextfloat` to give a faster result, with width 2 ulps
+    - `:none`: no rounding; does *not* guarantee that the results are correct.
+
+- `powers`:
+    Method to implement powers. Possible values are
+    - `:tight`: gives tightest possible result by using `BigFloat` internally, and hence is slow
+    - `:fast`: uses repeated multiplication based on `Base.power_by_squaring`, and hence is faster but gives slightly looser results
+"""
+
+const configuration = Dict(:directed_rounding => :tight,
+                           :powers => :fast)
+
+const allowed_values = Dict(:directed_rounding => (:tight, :fast, :none),
+                            :powers => (:tight, :fast))
+
+function configure!(; directed_rounding=nothing, powers=nothing)
+    if directed_rounding != nothing
+        if directed_rounding ∉ allowed_values[:directed_rounding]
+            throw(ArgumentError("directed_rounding must be one of $(allowed_values[:directed_rounding])"))
+        end
+
+        configuration[:directed_rounding] = directed_rounding
+
+        # name mismatch
+        if directed_rounding == :fast
+            directed_rounding = :accurate
+        end
+
+        set_directed_rounding(directed_rounding)
+    end
+
+    if powers != nothing
+        if powers ∉ allowed_values[:powers]
+            throw(ArgumentError("powers must be one of $(allowed_values[:powers])"))
+        end
+
+        configuration[:powers] = powers
+
+        set_power_type(powers)
+    end
+
+    return configuration
+
+end
+
+
+configure!(directed_rounding=:tight, powers=:fast)

src/intervals/functions.jl

@@ -1,23 +1,5 @@
 # This file is part of the IntervalArithmetic.jl package; MIT licensed
 
-## Powers
-
-# CRlibm does not contain a correctly-rounded ^ function for Float64
-# Use the BigFloat version from MPFR instead, which is correctly-rounded:
-
-# Write explicitly like this to avoid ambiguity warnings:
-for T in (:Integer, :Float64, :BigFloat, :Interval)
-    @eval ^(a::Interval{Float64}, x::$T) = atomic(Interval{Float64}, bigequiv(a)^x)
-end
-^(a::Interval{Float32}, x::Interval) = atomic(Interval{Float32}, bigequiv(a)^x)
-
-# Integer power:
-
-# overwrite new behaviour for small integer powers from
-# https://github.com/JuliaLang/julia/pull/24240:
-
-Base.literal_pow(::typeof(^), x::Interval{T}, ::Val{p}) where {T,p} = x^p
-
 
 Base.eltype(x::Interval{T}) where {T<:Real} = Interval{T}
 
@@ -36,206 +18,6 @@ Float64
 numtype(::Interval{T}) where {T<:Real} = T
 
 
-function ^(a::Interval{BigFloat}, n::Integer)
-    isempty(a) && return a
-    n == 0 && return one(a)
-    n == 1 && return a
-    # n == 2 && return sqr(a)
-    n < 0 && a == zero(a) && return emptyinterval(a)
-
-    T = BigFloat
-
-    if isodd(n) # odd power
-        isentire(a) && return a
-        if n > 0
-            a.lo == 0 && return @round(0, a.hi^n)
-            a.hi == 0 && return @round(a.lo^n, 0)
-            return @round(a.lo^n, a.hi^n)
-        else
-            if a.lo ≥ 0
-                a.lo == 0 && return @round(a.hi^n, Inf)
-                return @round(a.hi^n, a.lo^n)
-
-            elseif a.hi ≤ 0
-                a.hi == 0 && return @round(-Inf, a.lo^n)
-                return @round(a.hi^n, a.lo^n)
-            else
-                return entireinterval(a)
-            end
-        end
-
-    else # even power
-        if n > 0
-            if a.lo ≥ 0
-                return @round(a.lo^n, a.hi^n)
-            elseif a.hi ≤ 0
-                return @round(a.hi^n, a.lo^n)
-            else
-                return @round(mig(a)^n, mag(a)^n)
-            end
-
-        else
-            if a.lo ≥ 0
-                return @round(a.hi^n, a.lo^n)
-            elseif a.hi ≤ 0
-                return @round(a.lo^n, a.hi^n)
-            else
-                return @round(mag(a)^n, mig(a)^n)
-            end
-        end
-    end
-end
-
-function sqr(a::Interval{T}) where T<:Real
-    return a^2
-    # isempty(a) && return a
-    # if a.lo ≥ zero(T)
-    #     return @round(a.lo^2, a.hi^2)
-    #
-    # elseif a.hi ≤ zero(T)
-    #     return @round(a.hi^2, a.lo^2)
-    # end
-    #
-    # return @round(mig(a)^2, mag(a)^2)
-end
-
-^(a::Interval{BigFloat}, x::AbstractFloat) = a^big(x)
-
-# Floating-point power of a BigFloat interval:
-function ^(a::Interval{BigFloat}, x::BigFloat)
-
-    domain = Interval{BigFloat}(0, Inf)
-
-    if a == zero(a)
-        a = a ∩ domain
-        x > zero(x) && return zero(a)
-        return emptyinterval(a)
-    end
-
-    isinteger(x) && return a^(round(Int, x))
-    x == 0.5 && return sqrt(a)
-
-    a = a ∩ domain
-    (isempty(x) || isempty(a)) && return emptyinterval(a)
-
-    xx = atomic(Interval{BigFloat}, x)
-
-    # @round() can't be used directly, because both arguments may
-    # Inf or -Inf, which throws an error
-    # lo = @round(a.lo^xx.lo, a.lo^xx.lo)
-    lolod = @round_down(a.lo^xx.lo)
-    lolou = @round_up(a.lo^xx.lo)
-    lo = (lolod == Inf || lolou == -Inf) ?
-        wideinterval(lolod) : Interval(lolod, lolou)
-
-    # lo1 = @round(a.lo^xx.hi, a.lo^xx.hi)
-    lohid = @round_down(a.lo^xx.hi)
-    lohiu = @round_up(a.lo^xx.hi)
-    lo1 = (lohid == Inf || lohiu == -Inf) ?
-        wideinterval(lohid) : Interval(lohid, lohiu)
-
-    # hi = @round(a.hi^xx.lo, a.hi^xx.lo)
-    hilod = @round_down(a.hi^xx.lo)
-    hilou = @round_up(a.hi^xx.lo)
-    hi = (hilod == Inf || hilou == -Inf) ?
-        wideinterval(hilod) : Interval(hilod, hilou)
-
-    # hi1 = @round(a.hi^xx.hi, a.hi^xx.hi)
-    hihid = @round_down(a.hi^xx.hi)
-    hihiu = @round_up(a.hi^xx.hi)
-    hi1 = (hihid == Inf || hihiu == -Inf) ?
-        wideinterval(hihid) : Interval(hihid, hihiu)
-
-    lo = hull(lo, lo1)
-    hi = hull(hi, hi1)
-
-    return hull(lo, hi)
-end
-
-function ^(a::Interval{Rational{T}}, x::AbstractFloat) where T<:Integer
-    a = Interval(a.lo.num/a.lo.den, a.hi.num/a.hi.den)
-    a = a^x
-    atomic(Interval{Rational{T}}, a)
-end
-
-# Rational power
-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)
-        x > zero(x) && return zero(a)
-        return emptyinterval(a)
-    end
-
-    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
-
-        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
-
-    end
-
-    if x < 0
-        return inv(a^(-x))
-    end
-
-end
-
-# Interval power of an interval:
-function ^(a::Interval{BigFloat}, x::Interval)
-    T = BigFloat
-    domain = Interval{T}(0, Inf)
-
-    a = a ∩ domain
-
-    (isempty(x) || isempty(a)) && return emptyinterval(a)
-
-    lo1 = a^x.lo
-    lo2 = a^x.hi
-    lo1 = hull(lo1, lo2)
-
-    hi1 = a^x.lo
-    hi2 = a^x.hi
-    hi1 = hull(hi1, hi2)
-
-    hull(lo1, hi1)
-end
-
 
 function sqrt(a::Interval{T}) where T
     domain = Interval{T}(0, Inf)
@@ -248,46 +30,6 @@ end
 
 
 
-"""
-    pow(x::Interval, n::Integer)
-
-A faster implementation of `x^n`, currently using `power_by_squaring`.
-`pow(x, n)` will usually return an interval that is slightly larger than that calculated by `x^n`, but is guaranteed to be a correct
-enclosure when using multiplication with correct rounding.
-"""
-function pow(x::Interval, n::Integer)  # fast integer power
-    if n < 0
-        return 1/pow(x, -n)
-    end
-
-    isempty(x) && return x
-
-    if iseven(n) && 0 ∈ x
-
-        return hull(zero(x),
-                    hull(Base.power_by_squaring(Interval(mig(x)), n),
-                        Base.power_by_squaring(Interval(mag(x)), n))
-            )
-
-    else
-
-      return hull( Base.power_by_squaring(Interval(x.lo), n),
-                    Base.power_by_squaring(Interval(x.hi), n) )
-
-    end
-
-end
-
-function pow(x::Interval, y::Real)  # fast real power, including for y an Interval
-
-    isempty(x) && return x
-    isinteger(y) && return pow(x, Int(y.lo))
-    return exp(y * log(x))
-
-end
-
-
-
 
 for f in (:exp, :expm1)
     @eval begin

src/intervals/intervals.jl

@@ -124,7 +124,9 @@ end
 include("special.jl")
 include("macros.jl")
 include("rounding_macros.jl")
+include("powers.jl")
 include("rounding.jl")
+include("configuration.jl")
 include("conversion.jl")
 include("precision.jl")
 include("set_operations.jl")