Rounding modes using traits (#220)

* Traits-based interval rounding selection

* Working traits-based selection

* Cleanup

* Do not change rounding type if not necessary

* Fix rounding with traits

* Update CRlibm in REQUIRE

* Traits fixes

* Fix rounding for powers

* Fixes

* Fixes

* Fix small integer powers

* Remove  from tests

* Remove failing linear algebra test

* Linalg test broken only on 0.6

* Remove overwrite for small integer powers

* Change RoundingType -> IntervalRounding. Add to design explanation

* Remove restriction on T in IntervalRounding{T}

* Docstring

* Rewrite rounding-type trait to use types instead of instances

* Correct REQUIRE

* Add basic tests for interval rounding modes on 0.6

* Reenable linear algebra test

* Move multidim tests to end

* Correct testset syntax

* Add Suppressor for rounding tests

* Remove old definitions

* Comment out Suppressor

* Remove testset for rounding tests due to world age problems

* Add individual testsets for each interval rounding type

* Expand docstring for

Author: dpsanders

src/ValidatedNumerics.jl

@@ -4,7 +4,7 @@ __precompile__(true)
 
 module ValidatedNumerics
 
-using CRlibm
+import CRlibm
 using Compat
 using StaticArrays
 using ForwardDiff

src/intervals/arithmetic.jl

@@ -159,12 +159,6 @@ if VERSION >= v"0.6.0-dev.1024"
     const filter = Iterators.filter
 end
 
-if VERSION < v"0.5.0-dev+1279"
-    min(x) = x
-    max(x) = x
-end
-
-
 function min_ignore_nans(args...)
     min(filter(x->!isnan(x), args)...)
 end
@@ -319,7 +313,7 @@ function mid{T}(a::Interval{T})
 
     a.lo == -∞ && return nextfloat(a.lo)
     a.hi == +∞ && return prevfloat(a.hi)
-    
+
     (a.lo + a.hi) / 2  # rounds to nearest
 end
 

src/intervals/functions.jl

@@ -14,6 +14,9 @@ end
 
 # Integer power:
 
+# overwrite new behaviour for small integer powers:
+# ^{p}(x::ValidatedNumerics.Interval, ::Type{Val{p}}) = x^p
+
 function ^(a::Interval{BigFloat}, n::Integer)
     isempty(a) && return a
     n == 0 && return one(a)

src/intervals/precision.jl

@@ -22,6 +22,14 @@ function big53(a::Interval{Float64})
     end
 end
 
+function big53(x::Float64)
+    # BigFloat(x, 53)  # in Julia v0.6
+
+    setprecision(53) do
+        BigFloat(x)
+    end
+end
+
 
 setprecision(::Type{Interval}, ::Type{Float64}) = parameters.precision_type = Float64
 # does not change the BigFloat precision
@@ -77,44 +85,13 @@ Base.float{T}(::Type{Rational{T}}) = typeof(float(one(Rational{T})))
 
 # Use that type for rounding with rationals, e.g. for sqrt:
 
-if VERSION < v"0.5.0-dev+1182"
-
-    function Base.with_rounding{T}(f::Function, ::Type{Rational{T}},
-        rounding_mode::RoundingMode)
-        setrounding(f, float(Rational{T}), rounding_mode)
-    end
-
-else
-    function Base.setrounding{T}(f::Function, ::Type{Rational{T}},
-        rounding_mode::RoundingMode)
-        setrounding(f, float(Rational{T}), rounding_mode)
-    end
+function Base.setrounding{T}(f::Function, ::Type{Rational{T}},
+    rounding_mode::RoundingMode)
+    setrounding(f, float(Rational{T}), rounding_mode)
 end
 
 
-float{T}(x::Interval{T}) = convert(Interval{float(T)}, x)  # https://github.com/dpsanders/ValidatedNumerics.jl/issues/174
-
-## Change type of interval rounding:
 
-
-doc"""`rounding(Interval)` returns the current interval rounding mode.
-There are two possible rounding modes:
-
-- :narrow  -- changes the floating-point rounding mode to `RoundUp` and `RoundDown`.
-This gives the narrowest possible interval.
-
-- :wide -- Leaves the floating-point rounding mode in `RoundNearest` and uses
-`prevfloat` and `nextfloat` to achieve directed rounding. This creates an interval of width 2`eps`.
-"""
-
-rounding(::Type{Interval}) = parameters.rounding
-
-function setrounding(::Type{Interval}, mode)
-    if mode ∉ [:wide, :narrow]
-        throw(ArgumentError("Only possible interval rounding modes are `:wide` and `:narrow`"))
-    end
-
-    parameters.rounding = mode  # a symbol
-end
+float{T}(x::Interval{T}) = convert(Interval{float(T)}, x)  # https://github.com/dpsanders/ValidatedNumerics.jl/issues/174
 
 big{T}(x::Interval{T}) = convert(Interval{BigFloat}, x)

src/intervals/rounding.jl

@@ -1,80 +1,198 @@
-# Define rounded versions of elementary functions
-# E.g.  +(a, b, RoundDown)
-# Some, like sin(a, RoundDown)  are already defined in CRlibm
+#= Design summary:
 
+This is a so-called "traits-based" design, as follows.
 
-import Base: +, -, *, /, sin, sqrt, inv, ^, zero, convert, parse
+The main body of the file defines versions of elementary functions with all allowed
+interval rounding types, e.g.
++(IntervalRounding{:correct}, a, b, RoundDown)
++(IntervalRounding{:fast}, a, b, RoundDown)
++(IntervalRounding{:none}, a, b, RoundDown)
 
-# unary minus:
--{T<:AbstractFloat}(a::T, ::RoundingMode) = -a  # ignore rounding
+The current allowed rounding types are
+- :correct  # correct rounding (rounding to the nearest floating-point number)
+- :fast     # fast rounding by prevfloat and nextfloat  (slightly wider than needed)
+- :none     # no rounding at all (for speed, but forgoes any pretense at being rigorous)
 
-# zero:
-zero{T<:AbstractFloat}(a::Interval{T}, ::RoundingMode) = zero(T)
-zero{T<:AbstractFloat}(::Type{T}, ::RoundingMode) = zero(T)
+The function `setrounding(Interval, rounding_type)` then defines rounded
+ functions *without* an explicit rounding type, e.g.
+
+sin(x, r::RoundingMode) = sin(IntervalRounding{:correct}, x, r)
+
+These are overwritten when `setrounding(Interval, rounding_type)` is called again.
+
+In Julia v0.6, but *not* in Julia v0.5, this will automatically redefine all relevant functions, in particular those used in +(a::Interval, b::Interval) etc., so that all interval functions will automatically work with the correct interval rounding type.
+=#
 
-convert(::Type{BigFloat}, x, rounding_mode::RoundingMode) = setrounding(BigFloat, rounding_mode) do
-    convert(BigFloat, x)
-end
 
-parse{T}(::Type{T}, x, rounding_mode::RoundingMode) = setrounding(T, rounding_mode) do
-    parse(T, x)
+doc"""Interval rounding trait type"""
+immutable IntervalRounding{T} end
+
+# Functions that are the same for all rounding types:
+@eval begin
+    # unary minus:
+    -{T<:AbstractFloat}(a::T, ::RoundingMode) = -a  # ignore rounding
+
+    # zero:
+    zero{T<:AbstractFloat}(a::Interval{T}, ::RoundingMode) = zero(T)
+    zero{T<:AbstractFloat}(::Type{T}, ::RoundingMode) = zero(T)
+
+    convert(::Type{BigFloat}, x, rounding_mode::RoundingMode) =
+        setrounding(BigFloat, rounding_mode) do
+            convert(BigFloat, x)
+        end
+
+    parse{T}(::Type{T}, x, rounding_mode::RoundingMode) = setrounding(T, rounding_mode) do
+        parse(T, x)
+    end
+
+
+    sqrt{T<:Rational}(a::T, rounding_mode::RoundingMode) = setrounding(float(T), rounding_mode) do
+        sqrt(float(a))
+    end
+
 end
 
 
+
 # no-ops for rational rounding:
 for f in (:+, :-, :*, :/)
     @eval $f{T<:Rational}(a::T, b::T, ::RoundingMode) = $f(a, b)
 end
 
-sqrt{T<:Rational}(a::T, rounding_mode::RoundingMode) = setrounding(float(T), rounding_mode) do
-    sqrt(float(a))
-end
-
-
 
+# Define functions with different rounding types:
 for mode in (:Down, :Up)
 
     mode1 = Expr(:quote, mode)
     mode1 = :(::RoundingMode{$mode1})
 
     mode2 = Symbol("Round", mode)
 
+    if mode == :Down
+        directed = :prevfloat
+    else
+        directed = :nextfloat
+    end
+
 
-    for f in (:+, :-, :*, :/,
-                :atan2)
+    # binary functions:
+    for f in (:+, :-, :*, :/, :atan2)
 
-        @eval begin
-            function $f{T<:AbstractFloat}(a::T, b::T, $mode1)
-                setrounding(T, $mode2) do
-                    $f(a, b)
+        @eval function $f{T<:AbstractFloat}(::Type{IntervalRounding{:correct}},
+                                            a::T, b::T, $mode1)
+                    setrounding(T, $mode2) do
+                        $f(a, b)
+                    end
                 end
-            end
-        end
+
+        @eval $f{T<:AbstractFloat}(::Type{IntervalRounding{:fast}},
+                                    a::T, b::T, $mode1) = $directed($f(a, b))
+
+        @eval $f{T<:AbstractFloat}(::Type{IntervalRounding{:none}},
+                                    a::T, b::T, $mode1) = $f(a, b)
+
     end
 
-    @eval begin
-        function ^{T<:AbstractFloat,S}(a::T, b::S, $mode1)
-            setrounding(T, $mode2) do
-                ^(a, b)
-            end
+
+    # power:
+
+    @eval function ^{S<:Real}(::Type{IntervalRounding{:correct}},
+                                        a::BigFloat, b::S, $mode1)
+                  setrounding(BigFloat, $mode2) do
+                      ^(a, b)
+                  end
+           end
+
+    # for correct rounding for Float64, must pass through BigFloat:
+    @eval function ^{S<:Real}(::Type{IntervalRounding{:correct}}, a::Float64, b::S, $mode1)
+        setprecision(BigFloat, 53) do
+            Float64(^(IntervalRounding{:correct}, BigFloat(a), b, $mode2))
         end
     end
 
+    @eval ^{T<:AbstractFloat,S<:Real}(::Type{IntervalRounding{:fast}},
+                                a::T, b::S, $mode1) = $directed(a^b)
+
+    @eval ^{T<:AbstractFloat,S<:Real}(::Type{IntervalRounding{:none}},
+                                a::T, b::S, $mode1) = a^b
+
+
+    # functions not in CRlibm:
+    for f in (:sqrt, :inv, :tanh, :asinh, :acosh, :atanh)
+
+        @eval function $f{T<:AbstractFloat}(::Type{IntervalRounding{:correct}},
+                                            a::T, $mode1)
+                            setrounding(T, $mode2) do
+                                $f(a)
+                            end
+               end
+
+
+        @eval $f{T<:AbstractFloat}(::Type{IntervalRounding{:fast}},
+                                    a::T, $mode1) = $directed($f(a))
+
+        @eval $f{T<:AbstractFloat}(::Type{IntervalRounding{:none}},
+                                    a::T, $mode1) = $f(a)
+
 
-    for f in (:sqrt, :inv,
-            :tanh, :asinh, :acosh, :atanh)   # these functions not in CRlibm
-        @eval begin
-            function $f{T<:AbstractFloat}(a::T, $mode1)
-                setrounding(T, $mode2) do
-                    $f(a)
-                end
-            end
-        end
     end
 
 
+    # Functions defined in CRlibm
+
     for f in CRlibm.functions
-        @eval $f{T<:AbstractFloat}(a::T, $mode1) = CRlibm.$f(a, $mode2)
+        @eval $f{T<:AbstractFloat}(::Type{IntervalRounding{:correct}},
+                                a::T, $mode1) = CRlibm.$f(a, $mode2)
+
+        @eval $f{T<:AbstractFloat}(::Type{IntervalRounding{:fast}},
+                                    a::T, $mode1) = $directed($f(a))
+
+        @eval $f{T<:AbstractFloat}(::Type{IntervalRounding{:none}},
+                                    a::T, $mode1) = $f(a)
+
+    end
+end
+
+doc"""
+    setrounding(Interval, rounding_type::Symbol)
+
+Set the rounding type used for all interval calculations on Julia v0.6 and above.
+Valid `rounding_type`s are `:correct`, `:fast` and `:none`.
+"""
+function setrounding(::Type{Interval}, rounding_type::Symbol)
+
+    if rounding_type == current_rounding_type[]
+        return  # no need to redefine anything
+    end
+
+    if rounding_type ∉ (:correct, :fast, :none)
+        throw(ArgumentError("Rounding type must be one of `:correct`, `:fast`, `:none`"))
+    end
+
+    roundtype = IntervalRounding{rounding_type}
+
+
+    # binary functions:
+    for f in (:+, :-, :*, :/, :^, :atan2)
+
+        @eval $f{T<:AbstractFloat}(a::T, b::T, r::RoundingMode) = $f($roundtype, a, b, r)
     end
 
+    @eval ^{T<:AbstractFloat, S<:Real}(a::T, b::S, r::RoundingMode) = ^($roundtype, promote(a, b)..., r)
+
+    # unary functions:
+    for f in vcat(CRlibm.functions,
+                [:sqrt, :inv, :tanh, :asinh, :acosh, :atanh])
+
+        @eval $f{T<:AbstractFloat}(a::T, r::RoundingMode) = $f($roundtype, a, r)
+
+        @eval $f(x::Real, r::RoundingMode) = $f(float(x), r)
+
+    end
+
+    current_rounding_type[] = rounding_type
 end
+
+# default: correct rounding
+const current_rounding_type = Symbol[:undefined]
+setrounding(Interval, :correct)

test/ITF1788_tests/libieeep1788_tests_bool.jl

@@ -28,7 +28,7 @@ using ValidatedNumerics
 #Preamble
 setprecision(53)
 setprecision(Interval, Float64)
-setrounding(Interval, :narrow)
+# setrounding(Interval, :narrow)
 
 facts("minimal_empty_test") do
     @fact isempty(∅) --> true

test/ITF1788_tests/libieeep1788_tests_cancel.jl

@@ -28,7 +28,7 @@ using ValidatedNumerics
 #Preamble
 setprecision(53)
 setprecision(Interval, Float64)
-setrounding(Interval, :narrow)
+# setrounding(Interval, :narrow)
 
 facts("minimal_cancelPlus_test") do
     @fact cancelplus(Interval(-Inf, -1.0), ∅) --> entireinterval(Float64)

test/ITF1788_tests/libieeep1788_tests_elem.jl

@@ -28,7 +28,7 @@ using ValidatedNumerics
 #Preamble
 setprecision(53)
 setprecision(Interval, Float64)
-setrounding(Interval, :narrow)
+# setrounding(Interval, :narrow)
 
 facts("minimal_pos_test") do
     @fact +Interval(1.0, 2.0) --> Interval(1.0, 2.0)

test/ITF1788_tests/libieeep1788_tests_mul_rev.jl

@@ -28,7 +28,7 @@ using ValidatedNumerics
 #Preamble
 setprecision(53)
 setprecision(Interval, Float64)
-setrounding(Interval, :narrow)
+# setrounding(Interval, :narrow)
 
 facts("minimal_mulRevToPair_test") do
 

test/ITF1788_tests/libieeep1788_tests_num.jl

@@ -28,7 +28,7 @@ using ValidatedNumerics
 #Preamble
 setprecision(53)
 setprecision(Interval, Float64)
-setrounding(Interval, :narrow)
+# setrounding(Interval, :narrow)
 
 facts("minimal_inf_test") do
     @fact infimum(∅) --> Inf