Fix #790

docs/src/manual/configuration.md

@@ -22,6 +22,9 @@ Each keyword argument sets a specific configuration option:
 - `matmul`: control the matrix multiplication algorithm.
 - `nthreads`: control the number of threads used by the `:fast` matrix multiplication algorithm.
 
+!!! warning
+    Interval results are only correct when subnormal numbers are not flushed to zero. Do not call `set_zero_subnormals(true)`.
+
 ```@repl
 using IntervalArithmetic
 x = interval(π)

src/IntervalArithmetic.jl

@@ -285,7 +285,7 @@ _value_max(::Type{T}) where {T<:AbstractFloat} = floatmax(T)
 
 #
 
-bareinterval(::Type{BigFloat}, a::AbstractIrrational) = _unsafe_bareinterval(BigFloat, a, a)
+bareinterval(::Type{T}, a::AbstractIrrational) where {T<:Union{BigFloat,Rational{BigInt}}} = _unsafe_bareinterval(T, a, a)
 
 # Note: generated functions must be defined after all the methods they use
 @generated function bareinterval(::Type{T}, a::AbstractIrrational) where {T<:NumTypes}

src/intervals/arithmetic/hyperbolic.jl

@@ -199,13 +199,11 @@ function Base.acoth(x::BareInterval{T}) where {T<:AbstractFloat}
     singular_domain = _unsafe_bareinterval(T, -one(T), one(T))
     issubset_interval(x, singular_domain) && return emptyinterval(BareInterval{T})
     lo, hi = bounds(x)
-    if lo ≤ -1 && 1 ≤ hi
+    if lo < -1 && 1 < hi
         return entireinterval(BareInterval{T})
-    elseif lo < -1 < hi || lo < 1 < hi
-        return emptyinterval(BareInterval{T})
-    elseif lo == 1
+    elseif lo ≤ 1 < hi
         return @round(T, acoth(hi), typemax(T))
-    elseif hi == -1
+    elseif lo < -1 ≤ hi
         return @round(T, typemin(T), acoth(lo))
     else
         return @round(T, acoth(hi), acoth(lo))

src/intervals/arithmetic/trigonometric.jl

@@ -4,8 +4,8 @@
 
 # helper functions
 
-_big_pi(::T) where {T<:AbstractFloat} =
-    _unsafe_bareinterval(BigFloat, BigFloat(π, RoundDown; precision = 256), BigFloat(π, RoundUp; precision = 256))
+_big_pi(x::T) where {T<:AbstractFloat} =
+    _unsafe_bareinterval(BigFloat, BigFloat(π, RoundDown; precision = exponent(x) + precision(x) + 64), BigFloat(π, RoundUp; precision = exponent(x) + precision(x) + 64))
 
 _big_pi(x::BigFloat) =
     _unsafe_bareinterval(BigFloat, BigFloat(π, RoundDown; precision = precision(x)+32), BigFloat(π, RoundUp; precision = precision(x)+32))
@@ -22,7 +22,8 @@ function _quadrant(f, x::T) where {T<:AbstractFloat}
         x2 <  PI_HI && return f(0, 1) # [0, π/2) or [π/2, π)
         return 1 # [π/2, π)
     else
-        k = bareinterval(BigFloat, x2) / _big_pi(x)
+        X = BigFloat(x; precision = precision(x))
+        k = _unsafe_scale(_unsafe_bareinterval(BigFloat, X, X), BigFloat(2)) / _big_pi(x)
         fk = floor(k)
         lfk, hfk = bounds(fk)
         return Int(mod(f(lfk, hfk), 4))
@@ -279,7 +280,7 @@ Implement the `tan` function of the IEEE Standard 1788-2015 (Table 9.1).
 function Base.tan(x::BareInterval{T}) where {T<:AbstractFloat}
     isempty_interval(x) && return x
 
-    diam(x) > sup(bareinterval(T, π)) && return entireinterval(BareInterval{T})
+    diam(x) ≥ sup(bareinterval(T, π)) && return entireinterval(BareInterval{T})
 
     lo, hi = bounds(x)
 
@@ -317,7 +318,7 @@ This function is not part of the IEEE Standard 1788-2015.
 function Base.cot(x::BareInterval{T}) where {T<:AbstractFloat}
     isempty_interval(x) && return x
 
-    diam(x) > sup(bareinterval(T, π)) && return entireinterval(BareInterval{T})
+    diam(x) ≥ sup(bareinterval(T, π)) && return entireinterval(BareInterval{T})
 
     isthinzero(x) && return emptyinterval(BareInterval{T})
 

src/intervals/construction.jl

@@ -99,18 +99,18 @@ _round(::Type{T}, a, r::RoundingMode) where {T<:NumTypes} = __round(T, a, r)
 _round(::Type{<:NumTypes}, ::AbstractIrrational, ::RoundingMode) = throw(ArgumentError("only irrationals from MathConstants or IrrationalConstants.jl are supported"))
 for irr ∈ (:(:π), :(:γ), :(:catalan)) # irrationals supported by MPFR
     @eval _round(::Type{T}, a::Irrational{$irr}, r::RoundingMode) where {T<:NumTypes} =
-        __round(T, BigFloat(a, r), r)
+        __round(T, BigFloat(a, r; precision = precision(float(T))), r)
 end
 # irrationals not supported by MPFR, use their exact formula ℯ = exp(1), φ = (1+sqrt(5))/2.
 # The rounding mode is fixed to `:correct` rather than taken from `default_rounding()`:
 # the enclosure of a constant must be rigorous whatever the configuration, and the
 # generated `bareinterval(::Type{T}, ::AbstractIrrational)` bakes the result in at its
 # first generation, which method invalidation does not track.
 _round(::Type{T}, ::Irrational{:ℯ}, r::RoundingMode) where {T<:NumTypes} =
-    __round(T, _fround(exp, IntervalRounding{:correct}(), one(BigFloat), r), r)
+    __round(T, _fround(exp, IntervalRounding{:correct}(), BigFloat(1; precision = precision(float(T))), r), r)
 function _round(::Type{T}, ::Irrational{:φ}, r::RoundingMode) where {T<:NumTypes}
     ir = IntervalRounding{:correct}()
-    sqrt5 = _fround(sqrt, ir, BigFloat(5), r)
+    sqrt5 = _fround(sqrt, ir, BigFloat(5; precision = precision(float(T))), r)
     return __round(T, _fround(/, ir, _fround(+, ir, one(BigFloat), sqrt5, r), BigFloat(2), r), r)
 end
 # floats

src/intervals/rounding.jl

@@ -106,7 +106,7 @@ for f ∈ (:cbrt, :exp, :exp2, :exp10, :expm1, # exponential
                 bigz::Ref{BigFloat},
                 one(bigx)::Ref{BigFloat},
                 bigx::Ref{BigFloat},
-                RoundUp::Base.MPFR.MPFRRoundingMode
+                RoundDown::Base.MPFR.MPFRRoundingMode
             )::Int32
             bigw = _fround($g, ir, bigz, r)
             return BigFloat(bigw, r; precision = prec)
@@ -119,7 +119,7 @@ for f ∈ (:cbrt, :exp, :exp2, :exp10, :expm1, # exponential
                 bigz::Ref{BigFloat},
                 one(bigx)::Ref{BigFloat},
                 bigx::Ref{BigFloat},
-                RoundDown::Base.MPFR.MPFRRoundingMode
+                RoundUp::Base.MPFR.MPFRRoundingMode
             )::Int32
             bigw = _fround($g, ir, bigz, r)
             return BigFloat(bigw, r; precision = prec)

test/intervals/arithmetic/hyperbolic.jl

@@ -178,11 +178,16 @@ end
 @testset "acoth" begin
     @test isempty_interval(acoth(bareinterval(-1.0, 1.0)))
     @test isentire_interval(acoth(bareinterval(-2.0, 2.0)))
-    @test isempty_interval(acoth(bareinterval(0.5, 2.0)))
-    @test isempty_interval(acoth(bareinterval(-2.0, 0.5)))
+    @test isequal_interval(acoth(bareinterval(0.5, 2.0)), bareinterval(0.5493061443340548, Inf))
+    @test isequal_interval(acoth(bareinterval(-2.0, 0.5)), bareinterval(-Inf, -0.5493061443340548))
+    @test isequal_interval(acoth(bareinterval(-1.0, 2.0)), bareinterval(0.5493061443340548, Inf))
+    @test isequal_interval(acoth(bareinterval(-2.0, 1.0)), bareinterval(-Inf, -0.5493061443340548))
     @test isequal_interval(acoth(bareinterval(1.0, 2.0)), bareinterval(0.5493061443340548, Inf))
     @test isequal_interval(acoth(bareinterval(-2.0, -1.0)), bareinterval(-Inf, -0.5493061443340548))
     @test isequal_interval(acoth(bareinterval(2.0, 3.0)), bareinterval(0.34657359027997264, 0.5493061443340549))
+    x = 1.0122020710883537
+    r = acoth(bareinterval(x))
+    @test inf(r) ≤ acoth(big(x)) ≤ sup(r)
     @test isequal_interval(acoth(bareinterval(-3.0, -2.0)), bareinterval(-0.5493061443340549, -0.34657359027997264))
     @test isempty_interval(acoth(emptyinterval(BareInterval{Float64})))
     @test isentire_interval(acoth(entireinterval(BareInterval{Float64})))
@@ -193,7 +198,7 @@ end
     @test isempty_interval(r)
     @test decoration(r) == trv
     r = acoth(interval(0.5, 2.0))
-    @test isempty_interval(r)
+    @test isequal_interval(r, interval(0.5493061443340548, Inf))
     @test decoration(r) == trv
 end
 

test/intervals/arithmetic/trigonometric.jl

@@ -11,7 +11,7 @@ using Logging
     @test IntervalArithmetic._quadrant(min, 100.0) == 3
     @test IntervalArithmetic._quadrant(max, 100.0) == 3
     Logging.with_logger(Logging.NullLogger()) do
-        @test_throws InexactError IntervalArithmetic._quadrant(min, Inf)
+        @test_throws DomainError IntervalArithmetic._quadrant(min, Inf)
     end
 
     for x ∈ (1.5707963267948966, 1.5707963267948968, 3.141592653589793, 3.1415926535897936, -1.5707963267948966, -3.141592653589793)
@@ -24,7 +24,7 @@ using Logging
 
     bp = IntervalArithmetic._big_pi(1.0)
     @test bp isa BareInterval{BigFloat}
-    @test precision(inf(bp)) == 256
+    @test precision(inf(bp)) == precision(1.0) + 64
     @test inf(bp) < π < sup(bp)
     bp = IntervalArithmetic._big_pi(BigFloat(1))
     @test precision(inf(bp)) == precision(BigFloat(1)) + 32
@@ -263,6 +263,7 @@ end
     @test isequal_interval(tan(interval(0.5, 1.67)), entireinterval())
     @test isequal_interval(tan(interval(1.67, 3.2)), interval(-10.047182299210307, 0.05847385445957865))
     @test isequal_interval(tan(interval(6.638314112824137, 8.38263151220128)), entireinterval())
+    @test isentire_interval(tan(bareinterval(-3.1415926535897936, 0.0)))
 
     @test issubset_interval(tan(interval(BigFloat, 0.5, 0.5)), tan(interval(0.5)))
     @test isequal_interval(tan(interval(BigFloat, 0.5, 1.67)), entireinterval(BigFloat))
@@ -281,6 +282,7 @@ end
     @test isequal_interval(cot(bareinterval(-1.0, 0.0)), bareinterval(-Inf, -0.6420926159343306))
     @test isequal_interval(cot(bareinterval(2.0, 3.0)), bareinterval(-7.015252551434534, -0.4576575543602857))
     @test isentire_interval(cot(bareinterval(3.0, 3.5)))
+    @test isentire_interval(cot(bareinterval(-3.1415926535897936, 0.0)))
     @test isentire_interval(cot(entireinterval(BareInterval{Float64})))
     @test isempty_interval(cot(emptyinterval(BareInterval{Float64})))
     @test numtype(cot(bareinterval(1//2))) === Float64
@@ -373,6 +375,9 @@ end
     r = acot(bareinterval(1.0, 2.0))
     @test bounds(r) == (0.4636476090008061, 0.7853981633974484)
     @test inf(r) ≤ acot(2.0) && acot(1.0) ≤ sup(r)
+    x = 0x1.9c14c6477256cp-1
+    r = acot(bareinterval(x))
+    @test inf(r) ≤ acot(big(x)) ≤ sup(r)
     @test isempty_interval(acot(emptyinterval(BareInterval{Float64})))
 
     # acot(0) = π/2; the jump at 0 is decorated like the branch cut of two-argument atan
@@ -493,6 +498,15 @@ end
     @test isequal_interval(sin(x), interval(-1, 1))
     @test isequal_interval(cos(x), interval(-1, 1))
     @test isequal_interval(tan(x), interval(-Inf, Inf))
+
+    for f ∈ (sin, cos, tan, cot, sec, csc)
+        @test in_interval(f(big(1e308)), f(interval(1e308)))
+    end
+
+    @test isbounded(tan(interval(0x1.3c0ca428c59fbp+255)))
+    setprecision(BigFloat, 1000) do
+        @test isbounded(tan(interval(0x1.117203831ce2fp+626)))
+    end
 end
 
 @testset "inverse roots of unity" begin

test/intervals/construction.jl

@@ -147,6 +147,13 @@ end
     setprecision(BigFloat, 53) do
         @test bounds(bareinterval(BigFloat, 1, π)) == (1.0, 3.1415926535897936)
     end
+    setprecision(BigFloat, 20) do
+        for c ∈ (π, ℯ, MathConstants.φ)
+            @test nextfloat(IntervalArithmetic._round(Float64, c, RoundDown)) == IntervalArithmetic._round(Float64, c, RoundUp)
+        end
+    end
+    p20 = setprecision(() -> bareinterval(Rational{BigInt}, π), BigFloat, 20)
+    @test diam(bareinterval(Rational{BigInt}, π)) < diam(p20)
 
     @test bareinterval(1) === bareinterval(1.0, 1.0)
     @test typeof(bareinterval(1)) == BareInterval{Float64}