From fa7aab4105b4e097fc982f9fa6ae8e40f1077dda Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 31 May 2020 12:39:04 -0400
Subject: [PATCH 01/29] Add configuration.jl

---
 src/configuration.jl | 66 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 66 insertions(+)
 create mode 100644 src/configuration.jl

diff --git a/src/configuration.jl b/src/configuration.jl
new file mode 100644
index 000000000..f84d13a7a
--- /dev/null
+++ b/src/configuration.jl
@@ -0,0 +1,66 @@
+"""
+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
+"""
+mutable struct Configuration
+    directed_rounding::Symbol
+    powers::Symbol
+end
+
+function Base.show(io::IO, config::Configuration)
+    println(io, "IntervalArithmetic.Configuration:")
+    println(io, "- directed_rounding: $(config.directed_rounding)")
+    println(io, "- powers: $(config.powers)")
+end
+
+
+const configuration = Configuration(:tight, :fast)
+
+const directed_rounding_values = (:tight, :fast, :none)
+const power_values = (:tight, :fast)
+
+function configure!(; directed_rounding=nothing, powers=nothing)
+    if directed_rounding != nothing
+        if directed_rounding ∉ directed_rounding_values
+            throw(ArgumentError("directed_rounding must be one of $directed_rounding_values"))
+        end
+
+        configuration.directed_rounding = directed_rounding
+
+        # TODO: Implement!
+    end
+
+    if powers != nothing
+        if powers ∉ power_values
+            throw(ArgumentError("powers must be one of $directed_rounding_values"))
+        end
+
+        configuration.powers = powers
+
+        # TODO: Implement!
+    end
+
+    return configuration
+
+end
+
+configuration
+
+configure!(powers=:tight)
+
+
+ulp(f::F, x::T) where {F, T<:AbstractFloat} = 1

From ce61769fdf2e925beebfdcdf3c2849cd3904cc25 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sat, 6 Jun 2020 00:41:36 -0400
Subject: [PATCH 02/29] IntervalRounding -> DirectedRounding

---
 src/intervals/rounding.jl | 74 +++++++++++++++++++--------------------
 1 file changed, 37 insertions(+), 37 deletions(-)

diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index 557f98f3e..a4a6fdd41 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -4,11 +4,11 @@ This is a so-called "traits-based" design, as follows.
 
 The main body of the file defines versions of elementary functions with all allowed
 interval rounding types, e.g.
-+(IntervalRounding{:fast}, a, b, RoundDown)
-+(IntervalRounding{:tight}, a, b, RoundDown)
-+(IntervalRounding{:accurate}, a, b, RoundDown)
-+(IntervalRounding{:slow}, a, b, RoundDown)
-+(IntervalRounding{:none}, a, b, RoundDown)
++(DirectedRounding{:fast}, a, b, RoundDown)
++(DirectedRounding{:tight}, a, b, RoundDown)
++(DirectedRounding{:accurate}, a, b, RoundDown)
++(DirectedRounding{:slow}, a, b, RoundDown)
++(DirectedRounding{:none}, a, b, RoundDown)
 
 The current allowed rounding types are
 - :fast     # fast, tight (correct) rounding with errorfree arithmetic via FastRounding.jl
@@ -20,7 +20,7 @@ The current allowed rounding types are
 The function `setrounding(Interval, rounding_type)` then defines rounded
  functions *without* an explicit rounding type, e.g.
 
-sin(x, r::RoundingMode) = sin(IntervalRounding{:slow}, x, r)
+sin(x, r::RoundingMode) = sin(DirectedRounding{:slow}, x, r)
 
 These are overwritten when `setrounding(Interval, rounding_type)` is called again.
 
@@ -29,7 +29,7 @@ In Julia v0.6 and later (but *not* in Julia v0.5), this automatically redefines
 
 
 """Interval rounding trait type"""
-struct IntervalRounding{T} end
+struct DirectedRounding{T} end
 
 
 
@@ -89,7 +89,7 @@ for (op, f) in ( (:+, :add), (:-, :sub), (:*, :mul), (:/, :div) )
 
             mode2 = Symbol("Round", mode)
 
-            @eval $op(::IntervalRounding{:fast},
+            @eval $op(::DirectedRounding{:fast},
                                 a::$T, b::$T, $mode1) = $ff(a, b, $mode2)
         end
     end
@@ -97,9 +97,9 @@ end
 
 # inv and sqrt:
 
-# inv(::IntervalRounding{:fast}, a::T, r::RoundingMode) where T<:Union{Float32, Float64} = inv_round(a, r)
+# inv(::DirectedRounding{:fast}, a::T, r::RoundingMode) where T<:Union{Float32, Float64} = inv_round(a, r)
 #
-# sqrt(::IntervalRounding{:fast}, a::T, r::RoundingMode) where T<:Union{Float32, Float64} = sqrt_round(a, r)
+# sqrt(::DirectedRounding{:fast}, a::T, r::RoundingMode) where T<:Union{Float32, Float64} = sqrt_round(a, r)
 
 
 for T in (Float32, Float64)
@@ -110,10 +110,10 @@ for T in (Float32, Float64)
 
         mode2 = Symbol("Round", mode)
 
-        @eval inv(::IntervalRounding{:fast},
+        @eval inv(::DirectedRounding{:fast},
                             a::$T, $mode1) = inv_round(a, $mode2)
 
-        @eval sqrt(::IntervalRounding{:fast},
+        @eval sqrt(::DirectedRounding{:fast},
                             a::$T, $mode1) = sqrt_round(a, $mode2)
         end
 end
@@ -125,13 +125,13 @@ for T in (Float32, Float64)
             mode1 = Expr(:quote, mode)
             mode1 = :(::RoundingMode{$mode1})
             ff = Symbol(f, suffix)
-            @eval $op(::IntervalRounding{:tight},
+            @eval $op(::DirectedRounding{:tight},
                                 a::$T, b::$T, $mode1) = $ff(a, b)
         end
     end
 
-    @eval inv(::IntervalRounding{:tight}, a::$T, ::RoundingMode{:Down}) = div_down(one($T), a)
-    @eval inv(::IntervalRounding{:tight}, a::$T, ::RoundingMode{:Up}) = div_up(one($T), a)
+    @eval inv(::DirectedRounding{:tight}, a::$T, ::RoundingMode{:Down}) = div_down(one($T), a)
+    @eval inv(::DirectedRounding{:tight}, a::$T, ::RoundingMode{:Up}) = div_up(one($T), a)
 end
 
 
@@ -152,31 +152,31 @@ for mode in (:Down, :Up)
     # binary functions:
     for f in (:+, :-, :*, :/, :atan)
 
-        @eval function $f(::IntervalRounding{:slow},
+        @eval function $f(::DirectedRounding{:slow},
                           a::T, b::T, $mode1) where T<:AbstractFloat
                     setrounding(T, $mode2) do
                         $f(a, b)
                     end
                 end
 
-        @eval function $f(::IntervalRounding{:fast},
+        @eval function $f(::DirectedRounding{:fast},
                                   a::T, b::T, $mode1) where T<:AbstractFloat
                             setrounding(T, $mode2) do
                                 $f(a, b)
                             end
                         end
         
-        @eval function $f(::IntervalRounding{:tight},
+        @eval function $f(::DirectedRounding{:tight},
                             a::T, b::T, $mode1) where T<:AbstractFloat
                       setrounding(T, $mode2) do
                           $f(a, b)
                       end
                   end
 
-        @eval $f(::IntervalRounding{:accurate},
+        @eval $f(::DirectedRounding{:accurate},
                   a::T, b::T, $mode1) where {T<:AbstractFloat} = $directed($f(a, b))
 
-        @eval $f(::IntervalRounding{:none},
+        @eval $f(::DirectedRounding{:none},
                   a::T, b::T, $mode1) where {T<:AbstractFloat} = $f(a, b)
 
     end
@@ -184,7 +184,7 @@ for mode in (:Down, :Up)
 
     # power:
 
-    @eval function ^(::IntervalRounding{:slow},
+    @eval function ^(::DirectedRounding{:slow},
                                a::BigFloat, b::S, $mode1) where S<:Real
                   setrounding(BigFloat, $mode2) do
                       ^(a, b)
@@ -192,16 +192,16 @@ for mode in (:Down, :Up)
            end
 
     # for correct rounding for Float64, must pass through BigFloat:
-    @eval function ^(::IntervalRounding{:slow}, a::Float64, b::S, $mode1) where S<:Real
+    @eval function ^(::DirectedRounding{:slow}, a::Float64, b::S, $mode1) where S<:Real
         setprecision(BigFloat, 53) do
-            Float64(^(IntervalRounding{:slow}, BigFloat(a), b, $mode2))
+            Float64(^(DirectedRounding{:slow}, BigFloat(a), b, $mode2))
         end
     end
 
-    @eval ^(::IntervalRounding{:accurate},
+    @eval ^(::DirectedRounding{:accurate},
       a::T, b::S, $mode1) where {T<:AbstractFloat,S<:Real} = $directed(a^b)
 
-    @eval ^(::IntervalRounding{:none},
+    @eval ^(::DirectedRounding{:none},
       a::T, b::S, $mode1) where {T<:AbstractFloat,S<:Real} = a^b
 
 
@@ -209,27 +209,27 @@ for mode in (:Down, :Up)
     for f in (:sqrt, :inv, :tanh, :asinh, :acosh, :atanh, :cbrt)
 
 
-        @eval function $f(::IntervalRounding{:slow},
+        @eval function $f(::DirectedRounding{:slow},
                           a::T, $mode1) where T<:AbstractFloat
                             setrounding(T, $mode2) do
                                 $f(a)
                             end
                end
 
-        @eval function $f(::IntervalRounding{:fast},
+        @eval function $f(::DirectedRounding{:fast},
                                  a::T, $mode1) where T<:AbstractFloat
-                                   $f(IntervalRounding{:slow}(), a, $mode2)
+                                   $f(DirectedRounding{:slow}(), a, $mode2)
                       end
 
-        @eval function $f(::IntervalRounding{:tight},
+        @eval function $f(::DirectedRounding{:tight},
                         a::T, $mode1) where T<:AbstractFloat
-                          $f(IntervalRounding{:slow}(), a, $mode2)
+                          $f(DirectedRounding{:slow}(), a, $mode2)
                 end
 
-        @eval $f(::IntervalRounding{:accurate},
+        @eval $f(::DirectedRounding{:accurate},
                   a::T, $mode1) where {T<:AbstractFloat} = $directed($f(a))
 
-        @eval $f(::IntervalRounding{:none},
+        @eval $f(::DirectedRounding{:none},
                   a::T, $mode1) where {T<:AbstractFloat} = $f(a)
 
 
@@ -239,13 +239,13 @@ for mode in (:Down, :Up)
     # Functions defined in CRlibm
 
     for f in CRlibm.functions
-        @eval $f(::IntervalRounding{:slow},
+        @eval $f(::DirectedRounding{:slow},
               a::T, $mode1) where {T<:AbstractFloat} = CRlibm.$f(a, $mode2)
 
-        @eval $f(::IntervalRounding{:accurate},
+        @eval $f(::DirectedRounding{:accurate},
                   a::T, $mode1) where {T<:AbstractFloat} = $directed($f(a))
 
-        @eval $f(::IntervalRounding{:none},
+        @eval $f(::DirectedRounding{:none},
                   a::T, $mode1) where {T<:AbstractFloat} = $f(a)
 
     end
@@ -264,7 +264,7 @@ function _setrounding(::Type{Interval}, rounding_type::Symbol)
         throw(ArgumentError("Rounding type must be one of $rounding_types"))
     end
 
-    roundtype = IntervalRounding{rounding_type}()
+    roundtype = DirectedRounding{rounding_type}()
 
 
     # binary functions:
@@ -280,7 +280,7 @@ function _setrounding(::Type{Interval}, rounding_type::Symbol)
 
 
     if rounding_type in (:fast, :tight)   # for remaining functions, use CRlibm
-        roundtype = IntervalRounding{:slow}()
+        roundtype = DirectedRounding{:slow}()
     end
 
 

From 4f1a4b06cf30701019cdd58dc8a11936593a8b1b Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sat, 6 Jun 2020 00:48:06 -0400
Subject: [PATCH 03/29] Move power functions to powers.jl

---
 src/intervals/functions.jl | 281 -------------------------------------
 src/intervals/intervals.jl |   1 +
 2 files changed, 1 insertion(+), 281 deletions(-)

diff --git a/src/intervals/functions.jl b/src/intervals/functions.jl
index 3f3aadd8a..74483c2a4 100644
--- a/src/intervals/functions.jl
+++ b/src/intervals/functions.jl
@@ -1,286 +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}, big53(a)^x)
-end
-
-
-# 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} = 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)
-    a = a ∩ domain
-
-    isempty(a) && return a
-
-    @round(sqrt(a.lo), sqrt(a.hi))  # `sqrt` is correctly-rounded
-end
-
-
-function cbrt(a::Interval{T}) where T
-    isempty(a) && return a
-
-    @round(cbrt(a.lo), cbrt(a.hi))  
-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
         function ($f)(a::Interval{T}) where T
diff --git a/src/intervals/intervals.jl b/src/intervals/intervals.jl
index 78149832a..e89939e5b 100644
--- a/src/intervals/intervals.jl
+++ b/src/intervals/intervals.jl
@@ -133,6 +133,7 @@ include("conversion.jl")
 include("precision.jl")
 include("set_operations.jl")
 include("arithmetic.jl")
+include("powers.jl")
 include("functions.jl")
 include("trigonometric.jl")
 include("hyperbolic.jl")

From 0a5e2ac8acc52bb9955deea7512bb166f41e6e2f Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sat, 6 Jun 2020 00:48:31 -0400
Subject: [PATCH 04/29] Add powers.jl

---
 src/intervals/powers.jl | 278 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 278 insertions(+)
 create mode 100644 src/intervals/powers.jl

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
new file mode 100644
index 000000000..7a52a5064
--- /dev/null
+++ b/src/intervals/powers.jl
@@ -0,0 +1,278 @@
+
+## 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}, big53(a)^x)
+end
+
+
+# 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} = 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)
+    a = a ∩ domain
+
+    isempty(a) && return a
+
+    @round(sqrt(a.lo), sqrt(a.hi))  # `sqrt` is correctly-rounded
+end
+
+
+function cbrt(a::Interval{T}) where T
+    isempty(a) && return a
+
+    @round(cbrt(a.lo), cbrt(a.hi))
+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

From 089a7e6d94785e05b01bd0563e46dbb2b35d51ef Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sat, 6 Jun 2020 01:05:23 -0400
Subject: [PATCH 05/29] De-export setrounding and rename to
 set_directed_rounding

---
 src/IntervalArithmetic.jl |  2 +-
 src/intervals/rounding.jl | 17 +++++++++--------
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 8eb5657f1..6eab5efa8 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -80,7 +80,7 @@ export
 
 export showfull
 
-import Base: rounding, setrounding, setprecision
+import Base: setprecision
 
 
 
diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index a4a6fdd41..ee03b1853 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -165,7 +165,7 @@ for mode in (:Down, :Up)
                                 $f(a, b)
                             end
                         end
-        
+
         @eval function $f(::DirectedRounding{:tight},
                             a::T, b::T, $mode1) where T<:AbstractFloat
                       setrounding(T, $mode2) do
@@ -254,7 +254,8 @@ end
 
 const rounding_types = (:fast, :tight, :accurate, :slow, :none)
 
-function _setrounding(::Type{Interval}, rounding_type::Symbol)
+"Redefines all relevant functions to use the specified directed rounding mode"
+function _set_directed_rounding(rounding_type::Symbol)
 
     if rounding_type == current_rounding_type[]
         return  # no need to redefine anything
@@ -305,12 +306,12 @@ function _setrounding(::Type{Interval}, rounding_type::Symbol)
 end
 
 """
-    setrounding(Interval, rounding_type::Symbol)
+    set_directed_rounding(rounding_type::Symbol)
 
-Set the rounding type used for all interval calculations on Julia v0.6 and above.
+Set the directed rounding mode used for all interval calculations and above.
 Valid `rounding_type`s are $rounding_types.
 """
-function setrounding(::Type{Interval}, rounding_type::Symbol)
+function set_directed_rounding(rounding_type::Symbol)
 
     # suppress redefinition warnings:
     # modified from OhMyREPL.jl
@@ -323,7 +324,7 @@ function setrounding(::Type{Interval}, rounding_type::Symbol)
     old_stderr = stderr
     redirect_stderr(io)
 
-    _setrounding(Interval, rounding_type)
+    _set_directed_rounding(rounding_type)
 
     redirect_stderr(old_stderr)
 
@@ -356,10 +357,10 @@ function setrounding(::Type{Interval}, rounding_type::Symbol)
 
 end
 
-rounding(Interval) = current_rounding_type[]
+directed_rounding(Interval) = current_rounding_type[]
 
 
 
 # default: correct rounding
 const current_rounding_type = Symbol[:undefined]
-setrounding(Interval, :tight)
+set_directed_rounding(Interval, :tight)

From c6feb3f06d5a25bad9c9a32fcf46c05e374def35 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sat, 6 Jun 2020 01:17:34 -0400
Subject: [PATCH 06/29] Move sqrt and cbrt back to functions.jl. Add PowerType
 dispatch for powers

---
 src/intervals/functions.jl | 17 ++++++++++++++
 src/intervals/powers.jl    | 46 +++++++++++++++-----------------------
 2 files changed, 35 insertions(+), 28 deletions(-)

diff --git a/src/intervals/functions.jl b/src/intervals/functions.jl
index 74483c2a4..0a1865fbf 100644
--- a/src/intervals/functions.jl
+++ b/src/intervals/functions.jl
@@ -1,5 +1,22 @@
 # This file is part of the IntervalArithmetic.jl package; MIT licensed
 
+function sqrt(a::Interval{T}) where T
+    domain = Interval{T}(0, Inf)
+    a = a ∩ domain
+
+    isempty(a) && return a
+
+    @round(sqrt(a.lo), sqrt(a.hi))  # `sqrt` is correctly-rounded
+end
+
+
+function cbrt(a::Interval{T}) where T
+    isempty(a) && return a
+
+    @round(cbrt(a.lo), cbrt(a.hi))
+end
+
+
 for f in (:exp, :expm1)
     @eval begin
         function ($f)(a::Interval{T}) where T
diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 7a52a5064..3c5f03cc1 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -1,12 +1,18 @@
 
 ## Powers
 
+struct PowerType{T} end
+
+## Design:
+# PowerType{:tight} gives tight powers of intervals using BigFloat
+# PowerType{:fast} gives slightly less-tight but faster powers, using repeated multiplication (power_by_squaring), using the current directed rounding method
+
 # 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}, big53(a)^x)
+    @eval ^(::PowerType{:tight}, a::Interval{Float64}, x::$T) = atomic(Interval{Float64}, big53(a)^x)
 end
 
 
@@ -22,7 +28,7 @@ Base.eltype(x::Interval{T}) where {T<:Real} = T
 
 
 
-function ^(a::Interval{BigFloat}, n::Integer)
+function ^(::PowerType{:tight}, a::Interval{BigFloat}, n::Integer)
     isempty(a) && return a
     n == 0 && return one(a)
     n == 1 && return a
@@ -88,7 +94,7 @@ end
 ^(a::Interval{BigFloat}, x::AbstractFloat) = a^big(x)
 
 # Floating-point power of a BigFloat interval:
-function ^(a::Interval{BigFloat}, x::BigFloat)
+function ^(::PowerType{:tight}, a::Interval{BigFloat}, x::BigFloat)
 
     domain = Interval{BigFloat}(0, Inf)
 
@@ -145,7 +151,7 @@ function ^(a::Interval{Rational{T}}, x::AbstractFloat) where T<:Integer
 end
 
 # Rational power
-function ^(a::Interval{T}, x::Rational) where T
+function ^(::PowerType{:tight}, a::Interval{T}, x::Rational) where T
     domain = Interval{T}(0, Inf)
 
     p = x.num
@@ -203,7 +209,7 @@ function ^(a::Interval{T}, x::Rational) where T
 end
 
 # Interval power of an interval:
-function ^(a::Interval{BigFloat}, x::Interval)
+function ^(::PowerType{:tight}, a::Interval{BigFloat}, x::Interval)
     T = BigFloat
     domain = Interval{T}(0, Inf)
 
@@ -223,32 +229,16 @@ function ^(a::Interval{BigFloat}, x::Interval)
 end
 
 
-function sqrt(a::Interval{T}) where T
-    domain = Interval{T}(0, Inf)
-    a = a ∩ domain
-
-    isempty(a) && return a
-
-    @round(sqrt(a.lo), sqrt(a.hi))  # `sqrt` is correctly-rounded
-end
-
-
-function cbrt(a::Interval{T}) where T
-    isempty(a) && return a
-
-    @round(cbrt(a.lo), cbrt(a.hi))
-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
+A fast (?) implementation of `x^n` using `power_by_squaring`.
+Usually return an interval that is slightly larger than that calculated
+using PowerType{:tight}, but is guaranteed to be a correct
 enclosure when using multiplication with correct rounding.
 """
-function pow(x::Interval, n::Integer)  # fast integer power
+function ^(::PowerType{:fast}, x::Interval, n::Integer)  # fast integer power
     if n < 0
-        return 1/pow(x, -n)
+        return 1 / pow(x, -n)
     end
 
     isempty(x) && return x
@@ -269,7 +259,7 @@ function pow(x::Interval, n::Integer)  # fast integer power
 
 end
 
-function pow(x::Interval, y::Real)  # fast real power, including for y an Interval
+function pow(::PowerType{:fast}, x::Interval, y::Real)  # fast real power, including for y an Interval
 
     isempty(x) && return x
     isinteger(y) && return pow(x, Int(y.lo))

From 81134b960d8e87bef27441deabd67133bd222da1 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sat, 6 Jun 2020 01:43:28 -0400
Subject: [PATCH 07/29] Configuration seems to be working

---
 src/IntervalArithmetic.jl |  2 ++
 src/configuration.jl      | 40 +++++++++++++--------------------------
 src/intervals/powers.jl   |  8 ++++++++
 src/intervals/rounding.jl | 14 +++++++-------
 4 files changed, 30 insertions(+), 34 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 6eab5efa8..30945e6a3 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -116,6 +116,8 @@ include("multidim/multidim.jl")
 include("bisect.jl")
 include("decorations/decorations.jl")
 
+include("configuration.jl")
+
 include("rand.jl")
 include("parsing.jl")
 include("display.jl")
diff --git a/src/configuration.jl b/src/configuration.jl
index f84d13a7a..66a6af5aa 100644
--- a/src/configuration.jl
+++ b/src/configuration.jl
@@ -16,51 +16,37 @@ Holds configuration information for the package:
     - `: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
 """
-mutable struct Configuration
-    directed_rounding::Symbol
-    powers::Symbol
-end
-
-function Base.show(io::IO, config::Configuration)
-    println(io, "IntervalArithmetic.Configuration:")
-    println(io, "- directed_rounding: $(config.directed_rounding)")
-    println(io, "- powers: $(config.powers)")
-end
 
+const configuration = Dict(:directed_rounding => :tight,
+                           :powers => :fast)
 
-const configuration = Configuration(:tight, :fast)
-
-const directed_rounding_values = (:tight, :fast, :none)
-const power_values = (:tight, :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 ∉ directed_rounding_values
-            throw(ArgumentError("directed_rounding must be one of $directed_rounding_values"))
+        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
+        configuration[:directed_rounding] = directed_rounding
 
-        # TODO: Implement!
+        set_directed_rounding(directed_rounding)
     end
 
     if powers != nothing
-        if powers ∉ power_values
-            throw(ArgumentError("powers must be one of $directed_rounding_values"))
+        if powers ∉ allowed_values[:powers]
+            throw(ArgumentError("powers must be one of $(allowed_values[:powers])"))
         end
 
-        configuration.powers = powers
+        configuration[:powers] = powers
 
-        # TODO: Implement!
+        set_power_type(powers)
     end
 
     return configuration
 
 end
 
-configuration
-
-configure!(powers=:tight)
-
 
-ulp(f::F, x::T) where {F, T<:AbstractFloat} = 1
+# configure!(directed_rounding=:tight, powers=:tight)
diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 3c5f03cc1..1432d83dc 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -266,3 +266,11 @@ function pow(::PowerType{:fast}, x::Interval, y::Real)  # fast real power, inclu
     return exp(y * log(x))
 
 end
+
+
+function set_power_type(power_type)
+
+    type = PowerType{power_type}()
+
+    @eval ^(x::Interval{T}, n::S) where {T, S<:Integer} = ^($type, x::Interval{T}, n)
+end
diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index ee03b1853..1bcbe2528 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -257,9 +257,9 @@ const rounding_types = (:fast, :tight, :accurate, :slow, :none)
 "Redefines all relevant functions to use the specified directed rounding mode"
 function _set_directed_rounding(rounding_type::Symbol)
 
-    if rounding_type == current_rounding_type[]
-        return  # no need to redefine anything
-    end
+    # if rounding_type == current_rounding_type[]
+    #     return  # no need to redefine anything
+    # end
 
     if rounding_type ∉ rounding_types
         throw(ArgumentError("Rounding type must be one of $rounding_types"))
@@ -302,7 +302,7 @@ function _set_directed_rounding(rounding_type::Symbol)
 
     end
 
-    current_rounding_type[] = rounding_type
+    # current_rounding_type[] = rounding_type
 end
 
 """
@@ -357,10 +357,10 @@ function set_directed_rounding(rounding_type::Symbol)
 
 end
 
-directed_rounding(Interval) = current_rounding_type[]
+# directed_rounding(Interval) = current_rounding_type[]
 
 
 
 # default: correct rounding
-const current_rounding_type = Symbol[:undefined]
-set_directed_rounding(Interval, :tight)
+#const current_rounding_type = Symbol[:undefined]
+set_directed_rounding(:tight)

From cfddaef63555d4f9f301cb9aecde34d97084e6ff Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 15:50:43 -0400
Subject: [PATCH 08/29] Make powers non-recursive (incorporates
 nonrecursive_powers branch)

---
 src/intervals/powers.jl | 93 ++++++++++++++++++++++++++++++++++-------
 1 file changed, 79 insertions(+), 14 deletions(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 1432d83dc..27f061360 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -236,33 +236,61 @@ Usually return an interval that is slightly larger than that calculated
 using PowerType{:tight}, but is guaranteed to be a correct
 enclosure when using multiplication with correct rounding.
 """
-function ^(::PowerType{:fast}, x::Interval, n::Integer)  # fast integer power
+function ^(::PowerType{:fast}, x::Interval{T}, n::Integer) where {T}  # fast integer power
+
+    n == 0 && return one(x)
+
+    isempty(x) && return x
+
     if n < 0
-        return 1 / pow(x, -n)
+        negative_power = true
+        n = -n
+    else
+        negative_power = false
     end
 
-    isempty(x) && return x
+    if iseven(n)
+        if 0 ∈ x
 
-    if iseven(n) && 0 ∈ x
+            result =  Interval(zero(T),
+                        power_by_squaring(mag(x), n, RoundUp))
 
-        return hull(zero(x),
-                    hull(Base.power_by_squaring(Interval(mig(x)), n),
-                        Base.power_by_squaring(Interval(mag(x)), n))
-            )
+        elseif x.lo > 0
+            result = Interval(power_by_squaring(x.lo, n, RoundDown),
+                            power_by_squaring(x.hi, n, RoundUp))
 
-    else
+        else  # x.lo < x.hi < 0
+            result = Interval(power_by_squaring(-x.hi, n, RoundDown),
+                            power_by_squaring(-x.lo, n, RoundUp))
+        end
+
+    else  # odd n
 
-      return hull( Base.power_by_squaring(Interval(x.lo), n),
-                    Base.power_by_squaring(Interval(x.hi), n) )
+         a = power_by_squaring(x.lo, n, RoundDown)
+         b = power_by_squaring(x.hi, n, RoundUp)
+
+        result = Interval(a, b)
 
     end
+    #
+    # else  # completely negative interval
+    #     a = power_by_squaring(x.lo, n, RoundDown)
+    #     b = power_by_squaring(x.hi, n, RoundUp)
+    #
+    #     return Interval(a, b)
+    #end
 
-end
+    if negative_power
+        return inv(result)
+    else
+        return result
+    end
 
-function pow(::PowerType{:fast}, x::Interval, y::Real)  # fast real power, including for y an Interval
+end
+function ^(::PowerType{:fast}, x::Interval{T}, y::Real) where {T}  # fast real power, including for y an Interval
 
     isempty(x) && return x
-    isinteger(y) && return pow(x, Int(y.lo))
+    isinteger(y) && return x^(Int(y.lo))
     return exp(y * log(x))
 
 end
@@ -274,3 +302,40 @@ function set_power_type(power_type)
 
     @eval ^(x::Interval{T}, n::S) where {T, S<:Integer} = ^($type, x::Interval{T}, n)
 end
+
+
+
+# power_by_squaring adapted from Base Julia to add rounding mode:
+function power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
+
+    # assumes p is positive
+
+    if p == 1
+        return x
+    elseif p == 0
+        return one(x)
+    elseif p == 2
+        return *(x, x, r)  # multiplication with directed rounding
+    end
+    # elseif p < 0
+    #     isone(x) && return copy(x)
+    #     isone(-x) && return iseven(p) ? one(x) : copy(x)
+    #     Base.throw_domerr_powbysq(x, p)
+    # end
+
+    t = trailing_zeros(p) + 1
+    p >>= t
+    while (t -= 1) > 0
+        x = *(x, x, r)
+    end
+    y = x
+    while p > 0
+        t = trailing_zeros(p) + 1
+        p >>= t
+        while (t -= 1) >= 0
+            x = *(x, x, r)
+        end
+        y = *(y, x, r)
+    end
+    return y
+end

From 044c947cf7997a0d994668edfa12479f53a07ffb Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 15:51:48 -0400
Subject: [PATCH 09/29] Bump version of FastRounding in Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 948ab1a7b..58c461512 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,7 +16,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 CRlibm = "0.7, 0.8"
-FastRounding = "0.2"
+FastRounding = "0.2, 0.3"
 Polynomials = "0.7"
 RecipesBase = "1.0"
 RoundingEmulator = "0.2"

From 524cd4c82bdc87602aa04b79e3a18c6e976acfc1 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 16:29:15 -0400
Subject: [PATCH 10/29] Add Base annotation

---
 src/intervals/powers.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 27f061360..886e14e57 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -306,7 +306,7 @@ end
 
 
 # power_by_squaring adapted from Base Julia to add rounding mode:
-function power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
+function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
 
     # assumes p is positive
 

From de103ee12e0f77ef4d1d0bcfb5aea98d4f61fbe3 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 16:50:24 -0400
Subject: [PATCH 11/29] Import power_by_squaring

---
 src/IntervalArithmetic.jl | 3 ++-
 src/configuration.jl      | 2 +-
 src/intervals/powers.jl   | 2 +-
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 30945e6a3..8bf759bf3 100644
--- a/src/IntervalArithmetic.jl
+++ b/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,
diff --git a/src/configuration.jl b/src/configuration.jl
index 66a6af5aa..02810daf5 100644
--- a/src/configuration.jl
+++ b/src/configuration.jl
@@ -49,4 +49,4 @@ function configure!(; directed_rounding=nothing, powers=nothing)
 end
 
 
-# configure!(directed_rounding=:tight, powers=:tight)
+configure!(directed_rounding=:tight, powers=:fast)
diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 886e14e57..1fee29a23 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -300,7 +300,7 @@ function set_power_type(power_type)
 
     type = PowerType{power_type}()
 
-    @eval ^(x::Interval{T}, n::S) where {T, S<:Integer} = ^($type, x::Interval{T}, n)
+    @eval ^(x::Interval{T}, n::S) where {T, S<:Real} = ^($type, x::Interval{T}, n)
 end
 
 

From b2fb8ecc32a8a02de457a21acaca6cbf042ab6ad Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 16:52:20 -0400
Subject: [PATCH 12/29] Move configuration call to init

---
 src/IntervalArithmetic.jl | 2 ++
 src/configuration.jl      | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 8bf759bf3..e9bd7c7f6 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -107,6 +107,8 @@ function __init__()
 
     setprecision(Interval, 256)  # set up pi
     setprecision(Interval, Float64)
+
+    configure!(directed_rounding=:tight, powers=:fast)
 end
 
 
diff --git a/src/configuration.jl b/src/configuration.jl
index 02810daf5..f42e5c5b2 100644
--- a/src/configuration.jl
+++ b/src/configuration.jl
@@ -49,4 +49,4 @@ function configure!(; directed_rounding=nothing, powers=nothing)
 end
 
 
-configure!(directed_rounding=:tight, powers=:fast)
+# configure!(directed_rounding=:tight, powers=:fast)

From 7a8838c0e105c716ff2275d26dcd48ed2462288b Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 17:00:18 -0400
Subject: [PATCH 13/29] Fix power ambiguity

---
 src/intervals/powers.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 1fee29a23..93e8d0543 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -287,7 +287,7 @@ function ^(::PowerType{:fast}, x::Interval{T}, n::Integer) where {T}  # fast int
     end
 
 end
-function ^(::PowerType{:fast}, x::Interval{T}, y::Real) where {T}  # fast real power, including for y an Interval
+function ^(::PowerType{:fast}, x::Interval{T}, y::S) where {T, S<:Real}  # fast real power, including for y an Interval
 
     isempty(x) && return x
     isinteger(y) && return x^(Int(y.lo))
@@ -300,7 +300,9 @@ function set_power_type(power_type)
 
     type = PowerType{power_type}()
 
-    @eval ^(x::Interval{T}, n::S) where {T, S<:Real} = ^($type, x::Interval{T}, n)
+    for S in (Integer, Real)
+        @eval ^(x::Interval{T}, n::$S) where {T} = ^($type, x::Interval{T}, n)
+    end
 end
 
 

From c3a07c6d9b7f7f917c056cce855155b29884ae4c Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 17:07:26 -0400
Subject: [PATCH 14/29] Move negative power check

---
 src/intervals/powers.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 93e8d0543..abae0ed15 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -242,13 +242,13 @@ function ^(::PowerType{:fast}, x::Interval{T}, n::Integer) where {T}  # fast int
 
     isempty(x) && return x
 
+    negative_power = false
+
     if n < 0
         negative_power = true
         n = -n
-    else
-        negative_power = false
     end
-
+        
     if iseven(n)
         if 0 ∈ x
 

From 6e9ed32a5569acf9baf4deb9e90b39d269f99d50 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 17:12:17 -0400
Subject: [PATCH 15/29] Try some inlining

---
 src/intervals/rounding.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index 1bcbe2528..5de91c099 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -130,8 +130,8 @@ for T in (Float32, Float64)
         end
     end
 
-    @eval inv(::DirectedRounding{:tight}, a::$T, ::RoundingMode{:Down}) = div_down(one($T), a)
-    @eval inv(::DirectedRounding{:tight}, a::$T, ::RoundingMode{:Up}) = div_up(one($T), a)
+    @eval @inline inv(::DirectedRounding{:tight}, a::$T, ::RoundingMode{:Down}) = div_down(one($T), a)
+    @eval @inline inv(::DirectedRounding{:tight}, a::$T, ::RoundingMode{:Up}) = div_up(one($T), a)
 end
 
 

From 1e714b5d3218540c7b30505cfcafe19ddf493c02 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 17:58:21 -0400
Subject: [PATCH 16/29] Fix name mismatch

---
 src/configuration.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/configuration.jl b/src/configuration.jl
index f42e5c5b2..e8a4343b6 100644
--- a/src/configuration.jl
+++ b/src/configuration.jl
@@ -31,6 +31,11 @@ function configure!(; directed_rounding=nothing, powers=nothing)
 
         configuration[:directed_rounding] = directed_rounding
 
+        # name mismatch
+        if directed_rounding == :fast
+            directed_rounding = :accurate
+        end
+
         set_directed_rounding(directed_rounding)
     end
 

From 42f6c74d34d5235bc64d7b1e9e05bbd7ac70c671 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 19:08:48 -0400
Subject: [PATCH 17/29] More inlining

---
 src/multidim/intervalbox.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/multidim/intervalbox.jl b/src/multidim/intervalbox.jl
index 9c1d26b0e..97b5c59da 100644
--- a/src/multidim/intervalbox.jl
+++ b/src/multidim/intervalbox.jl
@@ -10,7 +10,11 @@ end
 # IntervalBox(x::Interval) = IntervalBox( SVector(x) )  # single interval treated as tuple with one element
 
 IntervalBox(x::Interval...) = IntervalBox(SVector(x))
-IntervalBox(x::SVector) = IntervalBox(interval.(x))
+
+IntervalBox(x::SVector{N,Interval{T}}) where {N,T}) = IntervalBox{N,T}(x)
+
+IntervalBox(x::SVector{N,T}) where {N,T}) = IntervalBox(interval.(x))
+
 IntervalBox(x::Tuple) = IntervalBox(SVector(x))
 IntervalBox(x::Real) = IntervalBox(interval.(x))
 IntervalBox(x...) = IntervalBox(x)

From 8989bb39c0f53b267f8b4f295b2f151f24714549 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 7 Jun 2020 19:13:36 -0400
Subject: [PATCH 18/29] Typos

---
 src/multidim/intervalbox.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/multidim/intervalbox.jl b/src/multidim/intervalbox.jl
index 97b5c59da..1bdfe149a 100644
--- a/src/multidim/intervalbox.jl
+++ b/src/multidim/intervalbox.jl
@@ -11,9 +11,7 @@ end
 
 IntervalBox(x::Interval...) = IntervalBox(SVector(x))
 
-IntervalBox(x::SVector{N,Interval{T}}) where {N,T}) = IntervalBox{N,T}(x)
-
-IntervalBox(x::SVector{N,T}) where {N,T}) = IntervalBox(interval.(x))
+IntervalBox(x::SVector{N,T}) where {N,T} = IntervalBox(interval.(x))
 
 IntervalBox(x::Tuple) = IntervalBox(SVector(x))
 IntervalBox(x::Real) = IntervalBox(interval.(x))

From d3aa7c0e179b092f6d62e2c648f6ef4bef7e87b8 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Tue, 9 Jun 2020 16:55:40 -0400
Subject: [PATCH 19/29] Fix warn -> @warn

---
 src/intervals/rounding.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index 5de91c099..239e134a7 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -344,11 +344,11 @@ function set_directed_rounding(rounding_type::Symbol)
         if !startswith(line, "WARNING: Method definition")
 
             if startswith(line, "WARNING") || startswith(line, "Use")
-                warn(line)
+                @warn(line)
 
             else
                 println("Error on line: ", line)
-                error(line)
+                @error(line)
             end
         end
     end

From 73fc00173f57b21dd20331ca56c402981d6c4cbd Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Tue, 9 Jun 2020 17:52:30 -0400
Subject: [PATCH 20/29] Fix (apparently) issues with precompilation

---
 src/IntervalArithmetic.jl |  2 +-
 src/intervals/rounding.jl | 10 ++++++----
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index e9bd7c7f6..e3af4be31 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -108,7 +108,7 @@ function __init__()
     setprecision(Interval, 256)  # set up pi
     setprecision(Interval, Float64)
 
-    configure!(directed_rounding=:tight, powers=:fast)
+    # configure!(directed_rounding=:tight, powers=:fast)
 end
 
 
diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index 239e134a7..2fbcaf5d3 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -319,17 +319,16 @@ function set_directed_rounding(rounding_type::Symbol)
     # dump the warnings to a file, and check the file to make
     # sure they are only redefinition warnings
 
-    path, io = mktemp()
+    path, io = mktemp(cleanup=false)
 
     old_stderr = stderr
-    redirect_stderr(io)
+    # redirect_stderr(io)
 
     _set_directed_rounding(rounding_type)
 
+    @show path
     redirect_stderr(old_stderr)
 
-    close(io)
-
     # check
     lines = readlines(path)
 
@@ -353,6 +352,9 @@ function set_directed_rounding(rounding_type::Symbol)
         end
     end
 
+    close(io)
+
+
     return rounding_type
 
 end

From d585a4de1b1e8254186542c9a4080aec04774966 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Tue, 9 Jun 2020 17:54:20 -0400
Subject: [PATCH 21/29] Comment out unnecessary lines

---
 src/intervals/rounding.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index 2fbcaf5d3..f6fa798db 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -321,13 +321,13 @@ function set_directed_rounding(rounding_type::Symbol)
 
     path, io = mktemp(cleanup=false)
 
-    old_stderr = stderr
+    # old_stderr = stderr
     # redirect_stderr(io)
 
     _set_directed_rounding(rounding_type)
 
-    @show path
-    redirect_stderr(old_stderr)
+    # @show path
+    # redirect_stderr(old_stderr)
 
     # check
     lines = readlines(path)

From 997e32cdde65d2d66e95abfb198e88fc8b2c549d Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Wed, 10 Jun 2020 23:59:06 -0400
Subject: [PATCH 22/29] Fix rational and real powers

---
 src/intervals/powers.jl | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index abae0ed15..622993a11 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -91,7 +91,7 @@ function sqr(a::Interval{T}) where T<:Real
     # return @round(mig(a)^2, mag(a)^2)
 end
 
-^(a::Interval{BigFloat}, x::AbstractFloat) = a^big(x)
+^(::PowerType{:tight}, a::Interval{BigFloat}, x::AbstractFloat) = a^big(x)
 
 # Floating-point power of a BigFloat interval:
 function ^(::PowerType{:tight}, a::Interval{BigFloat}, x::BigFloat)
@@ -144,7 +144,7 @@ function ^(::PowerType{:tight}, a::Interval{BigFloat}, x::BigFloat)
     return hull(lo, hi)
 end
 
-function ^(a::Interval{Rational{T}}, x::AbstractFloat) where T<:Integer
+function ^(::PowerType{:tight}, 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)
@@ -248,7 +248,7 @@ function ^(::PowerType{:fast}, x::Interval{T}, n::Integer) where {T}  # fast int
         negative_power = true
         n = -n
     end
-        
+
     if iseven(n)
         if 0 ∈ x
 
@@ -296,16 +296,27 @@ function ^(::PowerType{:fast}, x::Interval{T}, y::S) where {T, S<:Real}  # fast
 end
 
 
+# ^(::PowerType{:fast}, x::Interval{T}, n::Integer) where {T} = ^($type, x::Interval{T}, n)
+function ^(::PowerType{:fast}, x::Interval{T}, y::Rational) where {T}
+    ^(PowerType{:fast}(), x, Interval(y.num) / y.den)
+end
+
+
 function set_power_type(power_type)
 
     type = PowerType{power_type}()
 
-    for S in (Integer, Real)
-        @eval ^(x::Interval{T}, n::$S) where {T} = ^($type, x::Interval{T}, n)
+    for S in (Integer, Real, Rational)
+        @eval ^(x::Interval{T}, y::$S) where {T} = ^($type, x::Interval{T}, y)
     end
 end
 
 
+set_power_type(:tight)
+
+
+pow(x, y) = ^(PowerType{:fast}(), x, y)
+
 
 # power_by_squaring adapted from Base Julia to add rounding mode:
 function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)

From 3c00c9965ad9c1b0344f76355ccc6a18669801b5 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sat, 13 Jun 2020 12:53:00 -0400
Subject: [PATCH 23/29] Simplify real power

---
 src/intervals/powers.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 622993a11..9bb5ae2a0 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -290,8 +290,10 @@ end
 function ^(::PowerType{:fast}, x::Interval{T}, y::S) where {T, S<:Real}  # fast real power, including for y an Interval
 
     isempty(x) && return x
-    isinteger(y) && return x^(Int(y.lo))
-    return exp(y * log(x))
+
+    y2 = Interval(y)
+
+    return exp(y2 * log(x))
 
 end
 

From d35b4156d0eba4f881400aecb3b3b73a94b118bd Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 21 Jun 2020 21:17:50 -0400
Subject: [PATCH 24/29] Reinstate integer power

---
 src/intervals/powers.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 9bb5ae2a0..6cefed9df 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -291,6 +291,8 @@ function ^(::PowerType{:fast}, x::Interval{T}, y::S) where {T, S<:Real}  # fast
 
     isempty(x) && return x
 
+    isinteger(y) && return ^(PowerType{:fast}(), x, convert(Int, y.lo))
+
     y2 = Interval(y)
 
     return exp(y2 * log(x))

From 5751634550cbb46746a91dc7583e393d0c0f80a2 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Sun, 21 Jun 2020 22:27:22 -0400
Subject: [PATCH 25/29] Move configuration.jl into intervals; add call to
 configure

---
 src/IntervalArithmetic.jl            | 2 --
 src/{ => intervals}/configuration.jl | 2 +-
 src/intervals/intervals.jl           | 3 ++-
 src/intervals/powers.jl              | 2 +-
 src/intervals/rounding.jl            | 2 +-
 5 files changed, 5 insertions(+), 6 deletions(-)
 rename src/{ => intervals}/configuration.jl (97%)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index e3af4be31..f24852648 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -119,8 +119,6 @@ include("multidim/multidim.jl")
 include("bisect.jl")
 include("decorations/decorations.jl")
 
-include("configuration.jl")
-
 include("rand.jl")
 include("parsing.jl")
 include("display.jl")
diff --git a/src/configuration.jl b/src/intervals/configuration.jl
similarity index 97%
rename from src/configuration.jl
rename to src/intervals/configuration.jl
index e8a4343b6..ab21d62be 100644
--- a/src/configuration.jl
+++ b/src/intervals/configuration.jl
@@ -54,4 +54,4 @@ function configure!(; directed_rounding=nothing, powers=nothing)
 end
 
 
-# configure!(directed_rounding=:tight, powers=:fast)
+configure!(directed_rounding=:tight, powers=:fast)
diff --git a/src/intervals/intervals.jl b/src/intervals/intervals.jl
index e89939e5b..57f2ffbfc 100644
--- a/src/intervals/intervals.jl
+++ b/src/intervals/intervals.jl
@@ -128,12 +128,13 @@ 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")
 include("arithmetic.jl")
-include("powers.jl")
 include("functions.jl")
 include("trigonometric.jl")
 include("hyperbolic.jl")
diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 6cefed9df..d6bf28114 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -316,7 +316,7 @@ function set_power_type(power_type)
 end
 
 
-set_power_type(:tight)
+# set_power_type(:tight)
 
 
 pow(x, y) = ^(PowerType{:fast}(), x, y)
diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index f6fa798db..113626f00 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -365,4 +365,4 @@ end
 
 # default: correct rounding
 #const current_rounding_type = Symbol[:undefined]
-set_directed_rounding(:tight)
+# set_directed_rounding(:tight)

From b669b090a0168f0df381b57c38cac779d47cb233 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Mon, 21 Jun 2021 10:37:52 -0400
Subject: [PATCH 26/29] Update deps

---
 Project.toml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 58c461512..1093be033 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,14 +15,14 @@ SetRounding = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-CRlibm = "0.7, 0.8"
+CRlibm = "0.7, 0.8, 1"
 FastRounding = "0.2, 0.3"
-Polynomials = "0.7"
+Polynomials = "0.7, 1, 2"
 RecipesBase = "1.0"
 RoundingEmulator = "0.2"
 SetRounding = "0.2"
-StaticArrays = "0.8, 0.9, 0.10, 0.11, 0.12"
-julia = "1.3, 1.4, 1.5"
+StaticArrays = "0.8, 0.9, 0.10, 0.11, 0.12, 1"
+julia = "1.3, 1.4, 1.5, 1.6"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 57f401f4a9789154dcc44c17d752a5179b4de032 Mon Sep 17 00:00:00 2001
From: David Sanders <dsanders@Davids-MacBook-Pro.local>
Date: Fri, 17 Sep 2021 11:45:17 -0400
Subject: [PATCH 27/29] Fix power_by_squaring so that (-Inf..(-1e300))^3 is
 correct

---
 src/intervals/intervals.jl | 10 ++++++++++
 src/intervals/powers.jl    | 25 +++++++++++++++++++++++--
 2 files changed, 33 insertions(+), 2 deletions(-)

diff --git a/src/intervals/intervals.jl b/src/intervals/intervals.jl
index 57f2ffbfc..382d115cd 100644
--- a/src/intervals/intervals.jl
+++ b/src/intervals/intervals.jl
@@ -19,6 +19,10 @@ struct Interval{T<:Real} <: AbstractInterval{T}
 
     function Interval{T}(a::Real, b::Real) where T<:Real
 
+        if isinf(a) && a == b
+            return emptyinterval(T)
+        end
+    
         if validity_check
 
             if is_valid_interval(a, b)
@@ -86,6 +90,10 @@ function is_valid_interval(a::Real, b::Real)
         end
     end
 
+    if isinf(a) && a == b 
+        return true 
+    end
+
     if a > b
         if isinf(a) && isinf(b)
             return true  # empty interval = [∞,-∞]
@@ -107,6 +115,8 @@ end
 `interval(a, b)` checks whether [a, b] is a valid `Interval`, which is the case if `-∞ <= a <= b <= ∞`, using the (non-exported) `is_valid_interval` function. If so, then an `Interval(a, b)` object is returned; if not, then an error is thrown.
 """
 function interval(a::Real, b::Real)
+
+    
     if !is_valid_interval(a, b)
         throw(ArgumentError("`[$a, $b]` is not a valid interval. Need `a ≤ b` to construct `interval(a, b)`."))
     end
diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index d6bf28114..979aef695 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -327,8 +327,24 @@ function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
 
     # assumes p is positive
 
+    p_orig = p
+
+    # handle sign explicitly by calculating with positive
+    s = sign(x)
+    x = abs(x)
+
+    if s < 0 && isodd(p)
+        # we need to reverse the rounding mode:
+        if r == RoundDown 
+            r = RoundUp 
+        
+        elseif r == RoundUp
+            r = RoundDown 
+        end
+    end
+
     if p == 1
-        return x
+        return copysign(x, s)
     elseif p == 0
         return one(x)
     elseif p == 2
@@ -345,7 +361,9 @@ function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
     while (t -= 1) > 0
         x = *(x, x, r)
     end
+
     y = x
+    
     while p > 0
         t = trailing_zeros(p) + 1
         p >>= t
@@ -354,5 +372,8 @@ function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
         end
         y = *(y, x, r)
     end
-    return y
+
+    if isodd(p_orig)
+        return copysign(y, s)   
+    end
 end

From d40067e6e029b45d28fc5a4ba5dbfbf681258489 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Thu, 10 Feb 2022 22:36:30 -0500
Subject: [PATCH 28/29] Fix power_by_squaring

---
 src/intervals/powers.jl | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 979aef695..438f0a1b1 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -238,6 +238,8 @@ enclosure when using multiplication with correct rounding.
 """
 function ^(::PowerType{:fast}, x::Interval{T}, n::Integer) where {T}  # fast integer power
 
+    @show x, n
+
     n == 0 && return one(x)
 
     isempty(x) && return x
@@ -335,11 +337,11 @@ function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
 
     if s < 0 && isodd(p)
         # we need to reverse the rounding mode:
-        if r == RoundDown 
-            r = RoundUp 
-        
+        if r == RoundDown
+            r = RoundUp
+
         elseif r == RoundUp
-            r = RoundDown 
+            r = RoundDown
         end
     end
 
@@ -363,7 +365,7 @@ function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
     end
 
     y = x
-    
+
     while p > 0
         t = trailing_zeros(p) + 1
         p >>= t
@@ -374,6 +376,9 @@ function Base.power_by_squaring(x::AbstractFloat, p::Integer, r::RoundingMode)
     end
 
     if isodd(p_orig)
-        return copysign(y, s)   
+        return copysign(y, s)
     end
+
+    return y
+
 end

From 9e36c1409f5d5586ccc3358396d7f88ca53cce32 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Thu, 10 Feb 2022 22:36:59 -0500
Subject: [PATCH 29/29] Fix power_by_squaring

---
 src/intervals/powers.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/intervals/powers.jl b/src/intervals/powers.jl
index 438f0a1b1..45d58ae01 100644
--- a/src/intervals/powers.jl
+++ b/src/intervals/powers.jl
@@ -238,8 +238,6 @@ enclosure when using multiplication with correct rounding.
 """
 function ^(::PowerType{:fast}, x::Interval{T}, n::Integer) where {T}  # fast integer power
 
-    @show x, n
-
     n == 0 && return one(x)
 
     isempty(x) && return x