diff --git a/docs/src/manual/configuration.md b/docs/src/manual/configuration.md index 08fd2a8d..781ecda0 100644 --- a/docs/src/manual/configuration.md +++ b/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(π) diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index 7978085b..f0ae3c1e 100644 --- a/src/IntervalArithmetic.jl +++ b/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} diff --git a/src/intervals/arithmetic/hyperbolic.jl b/src/intervals/arithmetic/hyperbolic.jl index 31c1796f..4b2ef59d 100644 --- a/src/intervals/arithmetic/hyperbolic.jl +++ b/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)) diff --git a/src/intervals/arithmetic/trigonometric.jl b/src/intervals/arithmetic/trigonometric.jl index dcdd9bf3..d598749d 100644 --- a/src/intervals/arithmetic/trigonometric.jl +++ b/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}) diff --git a/src/intervals/construction.jl b/src/intervals/construction.jl index 2ae2e874..3b051600 100644 --- a/src/intervals/construction.jl +++ b/src/intervals/construction.jl @@ -99,7 +99,7 @@ _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()`: @@ -107,10 +107,10 @@ end # 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 diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl index 9f8aae53..4adbe6ff 100644 --- a/src/intervals/rounding.jl +++ b/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) diff --git a/test/intervals/arithmetic/hyperbolic.jl b/test/intervals/arithmetic/hyperbolic.jl index e6efbdc8..63354ba5 100644 --- a/test/intervals/arithmetic/hyperbolic.jl +++ b/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 diff --git a/test/intervals/arithmetic/trigonometric.jl b/test/intervals/arithmetic/trigonometric.jl index b0f9617d..e108eb5d 100644 --- a/test/intervals/arithmetic/trigonometric.jl +++ b/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 diff --git a/test/intervals/construction.jl b/test/intervals/construction.jl index 9f89f653..ab61022f 100644 --- a/test/intervals/construction.jl +++ b/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}