Merge pull request #414 from mfherbst/stackoverflowfixes

Extend BigFloat fallback to other floating-point types

src/IntervalArithmetic.jl

@@ -52,7 +52,7 @@ import .Broadcast: broadcasted
 export
     AbstractInterval, Interval,
     interval,
-    @interval, @biginterval, @floatinterval, @make_interval,
+    @interval, @biginterval, @floatinterval,
     diam, radius, mid, mag, mig, hull,
     emptyinterval, ∅, ∞, isempty, isinterior, ⪽, nthroot,
     precedes, strictprecedes, ≼, ≺, ⊂, ⊃, ⊇, contains_zero,

src/intervals/functions.jl

@@ -7,7 +7,7 @@
 
 # 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 ^(a::Interval{Float64}, x::$T) = atomic(Interval{Float64}, bigequiv(a)^x)
 end
 
 
@@ -294,7 +294,7 @@ for f in (:exp2, :exp10, :cbrt)
             end
         end
 
-    @eval ($f)(a::Interval{Float64}) = atomic(Interval{Float64}, $f(big53(a)))  # no CRlibm version
+    @eval ($f)(a::Interval{T}) where T = atomic(Interval{T}, $f(bigequiv(a)))  # no CRlibm version
 
     @eval function ($f)(a::Interval{BigFloat})
             isempty(a) && return a
@@ -327,7 +327,7 @@ end
 
 hypot(a::Interval{BigFloat}, b::Interval{BigFloat}) = sqrt(a^2 + b^2)
 
-hypot(a::Interval{T}, b::Interval{T}) where T= atomic(Interval{T}, hypot(big53(a), big53(b)))
+hypot(a::Interval{T}, b::Interval{T}) where T= atomic(Interval{T}, hypot(bigequiv(a), bigequiv(b)))
 
 """
 nthroot(a::Interval{BigFloat}, n::Integer)
@@ -357,6 +357,6 @@ function nthroot(a::Interval{BigFloat}, n::Integer)
 end
 
 function nthroot(a::Interval{T}, n::Integer) where T
-    b = nthroot(big53(a), n)
+    b = nthroot(bigequiv(a), n)
     return convert(Interval{T}, b)
-end 
+end 

src/intervals/hyperbolic.jl

@@ -16,7 +16,7 @@ end
 # function tanh(a::Interval{Float64})
 #     isempty(a) && return a
 #
-#     float(tanh(big53(a)))
+#     float(tanh(bigequiv(a)))
 # end
 
 function tanh(a::Interval{BigFloat})
@@ -66,9 +66,9 @@ sech(a::Interval{BigFloat}) = 1/cosh(a)
 # Float64 versions of functions missing from CRlibm:
 for f in (:tanh, :asinh, :acosh, :atanh, :coth, :csch, :sech)
 
-    @eval function ($f)(a::Interval{Float64})
+    @eval function ($f)(a::Interval{T}) where T
         isempty(a) && return a
 
-        atomic(Interval{Float64}, ($f)(big53(a)) )
+        atomic(Interval{T}, ($f)(bigequiv(a)) )
     end
 end

src/intervals/macros.jl

@@ -16,8 +16,14 @@ Examples:
     @interval(1/3^2)
 ```
 """
-macro interval(expr1, expr2...)
-    make_interval(:(parameters.precision_type), expr1, expr2)
+macro interval(expr1)
+    make_interval(:(parameters.precision_type), expr1, ())
+end
+macro interval(expr1, expr2)
+    make_interval(:(parameters.precision_type), expr1, (expr2, ))
+end
+macro interval(T, expr1, expr2)
+    make_interval(T, expr1, (expr2, ))
 end
 
 "The `@floatinterval` macro constructs an interval with `Float64` entries."
@@ -78,7 +84,7 @@ end
 and making each literal (0.1, 1, etc.) into a corresponding interval construction,
 by calling `transform`."""
 function make_interval(T, expr1, expr2)
-    # @show expr1, expr2
+    # @show T, expr1, expr2
 
     expr1 = transform(expr1, :atomic, :(Interval{$T}))
 

src/intervals/precision.jl

@@ -15,17 +15,15 @@ const parameters = IntervalParameters()
 
 ## Precision:
 
-"`big53` creates an equivalent `BigFloat` interval to a given `Float64` interval."
-function big53(a::Interval{Float64})
-    setprecision(Interval, 53) do  # precision of Float64
+"`bigequiv` creates an equivalent `BigFloat` interval to a given `AbstractFloat` interval."
+function bigequiv(a::Interval{T}) where {T <: AbstractFloat}
+    setprecision(Interval, precision(T)) do  # precision of T
         atomic(Interval{BigFloat}, a)
     end
 end
 
-function big53(x::Float64)
-    # BigFloat(x, 53)  # in Julia v0.6
-
-    setprecision(53) do
+function bigequiv(x::AbstractFloat)
+    setprecision(precision(x)) do
         BigFloat(x)
     end
 end

src/intervals/rounding.jl

@@ -53,15 +53,21 @@ parse(::Type{T}, x::AbstractString, rounding_mode::RoundingMode) where {T} = set
     parse(T, x)
 end
 
-# use BigFloat parser to get round issues on Windows:
-function parse(::Type{Float64}, s::AbstractString, r::RoundingMode)
-    a = setprecision(BigFloat, 53) do
-            setrounding(BigFloat, r) do
-                parse(BigFloat, s)   # correctly takes account of rounding mode
+# use higher precision float parser to get round issues on Windows
+@static if Sys.iswindows()
+    function parse(::Type{Float64}, s::AbstractString, r::RoundingMode)
+        a = setprecision(BigFloat, 53) do
+                setrounding(BigFloat, r) do
+                    parse(BigFloat, s)   # correctly takes account of rounding mode
+                end
             end
-        end
 
-    return Float64(a, r)
+        return Float64(a, r)
+    end
+
+    function parse(::Type{T}, s::AbstractString, r::RoundingMode) where {T <: Union{Float16, Float32}}
+        return T(parse(Float64, s, r), r)
+    end
 end
 
 

src/intervals/trigonometric.jl

@@ -323,10 +323,10 @@ end
 
 #atan{T<:Real, S<:Real}(y::Interval{T}, x::Interval{S}) = atan(promote(y, x)...)
 
-function atan(y::Interval{Float64}, x::Interval{Float64})
-    (isempty(y) || isempty(x)) && return emptyinterval(Float64)
+function atan(y::Interval{T}, x::Interval{T}) where T
+    (isempty(y) || isempty(x)) && return emptyinterval(T)
 
-    atomic(Interval{Float64}, atan(big53(y), big53(x)))
+    atomic(Interval{T}, atan(bigequiv(y), bigequiv(x)))
 end
 
 
@@ -430,6 +430,6 @@ for f in (:cot, :csc, :sec)
     @eval function ($f)(a::Interval{T}) where T
         isempty(a) && return a
 
-        atomic(Interval{T}, ($f)(big53(a)) )
+        atomic(Interval{T}, ($f)(bigequiv(a)) )
     end
-end
+end

test/interval_tests/construction.jl

@@ -216,6 +216,15 @@ end
     # @test imag(b) == Interval(-15.200784463067956, -15.20078446306795)
 end
 
+@testset "Typed intervals" begin
+    setprecision(Interval, Float64)
+
+    @test typeof(@interval Float64 1 2) == Interval{Float64}
+    @test typeof(@interval         1 2) == Interval{Float64}
+    @test typeof(@interval Float32 1 2) == Interval{Float32}
+    @test typeof(@interval Float16 1 2) == Interval{Float16}
+end
+
 @testset ".. tests" begin
 
     a = 0.1..0.3
@@ -337,7 +346,7 @@ end
     @test isequal(hash(x), hash(y))
 
     x = @interval(0.1)
-    y = IntervalArithmetic.big53(x)
+    y = IntervalArithmetic.bigequiv(x)
     @test isequal(x, y)
     @test isequal(hash(x), hash(y))
 

test/interval_tests/hyperbolic.jl

@@ -46,6 +46,14 @@ setprecision(Interval, Float64)
     @test tanh(@biginterval(-4.5, 0.1)) ⊆ tanh(@interval(-4.5, 0.1))
     @test tanh(@biginterval(1.3, 6.3)) ⊆ tanh(@interval(1.3, 6.3))
 
+    @test tanh(@interval(0.5))       ⊆ tanh(@interval(Float32, 0.5, 0.5))
+    @test tanh(@interval(0.5, 1.67)) ⊆ tanh(@interval(Float32, 0.5, 1.67))
+    @test tanh(@interval(1.67, 3.2)) ⊆ tanh(@interval(Float32, 1.67, 3.2))
+    @test tanh(@interval(2.1, 5.6))  ⊆ tanh(@interval(Float32, 2.1, 5.6))
+    @test tanh(@interval(0.5, 8.5))  ⊆ tanh(@interval(Float32, 0.5, 8.5))
+    @test tanh(@interval(-4.5, 0.1)) ⊆ tanh(@interval(Float32, -4.5, 0.1))
+    @test tanh(@interval(1.3, 6.3))  ⊆ tanh(@interval(Float32, 1.3, 6.3))
+
     @test asinh(@biginterval(1)) ⊆ asinh(@interval(1))
     @test asinh(@biginterval(0.9, 2)) ⊆ asinh(@interval(0.9, 2))
     @test asinh(@biginterval(3, 4)) ⊆ asinh(@interval(3, 4))

test/interval_tests/numeric.jl

@@ -361,6 +361,9 @@ end
         @test cbrt(big(2..3)) ⊆ cbrt(2..3)
 
         @test cbrt(big(3..4)) != cbrt(3..4)
+
+        @test cbrt(2f0..3f0) == Interval(1.259921f0, 1.4422497f0)
+        @test cbrt(2..3) ⊆ cbrt(2f0..3f0)
     end
 
     @testset "inv" begin

test/interval_tests/trig.jl

@@ -188,6 +188,21 @@ end
         Interval(-1.5707963267948968, 1.5707963267948968)
     @test atan(@biginterval(-0.1, 0.1), @biginterval(-0.0, 0.1)) ⊆
         atan(@interval(-0.1, 0.1), @interval(0.0, 0.1))
+
+    @test atan(@interval(Float32, -0.1, 0.1), @interval(Float32, 0.1, 0.1)) ==
+        Interval(-0.78539824f0, 0.78539824f0)
+    @test atan(@interval(-0.1, 0.1), @interval(0.1, 0.1)) ⊆
+        atan(@interval(Float32, -0.1, 0.1), @interval(Float32, 0.1, 0.1))
+    @test atan(@interval(Float32, 0.0, 0.0), @interval(Float32, -0.0, 0.1)) ==
+        @interval(Float32, 0.0, 0.0)
+    @test atan(@interval(Float32, 0.0, 0.1), @interval(Float32, -0.0, 0.1)) ==
+        Interval(0.0, 1.5707964f0)
+    @test atan(@interval(Float32, -0.1, 0.0), @interval(Float32, 0.0, 0.1)) ==
+        Interval(-1.5707964f0, 0.0)
+    @test atan(@interval(Float32, -0.1, 0.1), @interval(Float32, -0.0, 0.1)) ==
+        Interval(-1.5707964f0, 1.5707964f0)
+    @test atan(@interval(-0.1, 0.1), @interval(-0.0, 0.1)) ⊆
+        atan(@interval(Float32, -0.1, 0.1), @interval(Float32, 0.0, 0.1))
 end
 
 @testset "Trig" begin