From bbcea6fba6d33dfe643936d444b31d6575ddd9c8 Mon Sep 17 00:00:00 2001 From: OlivierHnt <38465572+OlivierHnt@users.noreply.github.com> Date: Wed, 2 Oct 2024 12:37:21 -0400 Subject: [PATCH] Rework fix --- src/intervals/arithmetic/trigonometric.jl | 54 +++++++++++------------ 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/src/intervals/arithmetic/trigonometric.jl b/src/intervals/arithmetic/trigonometric.jl index 1c9ab3b1..1b89b5f9 100644 --- a/src/intervals/arithmetic/trigonometric.jl +++ b/src/intervals/arithmetic/trigonometric.jl @@ -9,17 +9,17 @@ function _quadrant(f, x::T) where {T<:AbstractFloat} PI_LO, PI_HI = bounds(PI) if abs(x) ≤ PI_LO # (-π, π) r2 = 2x # should be exact for floats - r2 ≤ -PI_HI && return 2, floor(x) # (-π, -π/2) - r2 < -PI_LO && return f(2, 3), floor(x) # (-π, -π/2) or [-π/2, 0) - r2 < 0 && return 3, floor(x) # [-π/2, 0) - r2 ≤ PI_LO && return 0, floor(x) # [0, π/2) - r2 < PI_HI && return f(0, 1), floor(x) # [0, π/2) or [π/2, π) - return 1, floor(x) # [π/2, π) + r2 ≤ -PI_HI && return 2 # (-π, -π/2) + r2 < -PI_LO && return f(2, 3) # (-π, -π/2) or [-π/2, 0) + r2 < 0 && return 3 # [-π/2, 0) + r2 ≤ PI_LO && return 0 # [0, π/2) + r2 < PI_HI && return f(0, 1) # [0, π/2) or [π/2, π) + return 1 # [π/2, π) else k = _unsafe_scale(bareinterval(x) / PI, convert(T, 2)) fk = floor(k) lfk, hfk = bounds(fk) - return f(Int(mod(lfk, 4)), Int(mod(hfk, 4))), f(lfk, hfk) + return f(Int(mod(lfk, 4)), Int(mod(hfk, 4))) end end @@ -69,20 +69,20 @@ function Base.sin(x::BareInterval{T}) where {T<:AbstractFloat} lo, hi = bounds(x) - lo_quadrant, lq = _quadrant(min, lo) - hi_quadrant, hq = _quadrant(max, hi) + lo_quadrant = _quadrant(min, lo) + hi_quadrant = _quadrant(max, hi) - if hq - lq > 3 - return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 2π - - elseif lo_quadrant == hi_quadrant + if lo_quadrant == hi_quadrant + d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 2π (lo_quadrant == 1) | (lo_quadrant == 2) && return @round(T, sin(hi), sin(lo)) # decreasing return @round(T, sin(lo), sin(hi)) elseif lo_quadrant == 3 && hi_quadrant == 0 + d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 3π/2 return @round(T, sin(lo), sin(hi)) # increasing elseif lo_quadrant == 1 && hi_quadrant == 2 + d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 3π/2 return @round(T, sin(hi), sin(lo)) # decreasing elseif (lo_quadrant == 0 || lo_quadrant == 3) && (hi_quadrant == 1 || hi_quadrant == 2) @@ -173,20 +173,20 @@ function Base.cos(x::BareInterval{T}) where {T<:AbstractFloat} lo, hi = bounds(x) - lo_quadrant, lq = _quadrant(min, lo) - hi_quadrant, hq = _quadrant(max, hi) + lo_quadrant = _quadrant(min, lo) + hi_quadrant = _quadrant(max, hi) - if hq - lq > 3 - return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 2π - - elseif lo_quadrant == hi_quadrant + if lo_quadrant == hi_quadrant + d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 2π (lo_quadrant == 2) | (lo_quadrant == 3) && return @round(T, cos(lo), cos(hi)) # increasing return @round(T, cos(hi), cos(lo)) elseif lo_quadrant == 2 && hi_quadrant == 3 + d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 3π/2 return @round(T, cos(lo), cos(hi)) elseif lo_quadrant == 0 && hi_quadrant == 1 + d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 3π/2 return @round(T, cos(hi), cos(lo)) elseif (lo_quadrant == 2 || lo_quadrant == 3) && (hi_quadrant == 0 || hi_quadrant == 1) @@ -277,8 +277,8 @@ function Base.tan(x::BareInterval{T}) where {T<:AbstractFloat} lo, hi = bounds(x) - lo_quadrant, _ = _quadrant(min, lo) - hi_quadrant, _ = _quadrant(max, hi) + lo_quadrant = _quadrant(min, lo) + hi_quadrant = _quadrant(max, hi) lo_quadrant_mod = mod(lo_quadrant, 2) hi_quadrant_mod = mod(hi_quadrant, 2) @@ -317,8 +317,8 @@ function Base.cot(x::BareInterval{T}) where {T<:AbstractFloat} lo, hi = bounds(x) - lo_quadrant, _ = _quadrant(min, lo) - hi_quadrant, _ = _quadrant(max, hi) + lo_quadrant = _quadrant(min, lo) + hi_quadrant = _quadrant(max, hi) if (lo_quadrant == 2 || lo_quadrant == 3) && hi == 0 return @round(T, typemin(T), cot(lo)) # singularity from the left @@ -349,8 +349,8 @@ function Base.sec(x::BareInterval{T}) where {T<:AbstractFloat} lo, hi = bounds(x) - lo_quadrant, _ = _quadrant(min, lo) - hi_quadrant, _ = _quadrant(max, hi) + lo_quadrant = _quadrant(min, lo) + hi_quadrant = _quadrant(max, hi) if lo_quadrant == hi_quadrant (lo_quadrant == 0) | (lo_quadrant == 1) && return @round(T, sec(lo), sec(hi)) # increasing @@ -387,8 +387,8 @@ function Base.csc(x::BareInterval{T}) where {T<:AbstractFloat} lo, hi = bounds(x) - lo_quadrant, _ = _quadrant(min, lo) - hi_quadrant, _ = _quadrant(max, hi) + lo_quadrant = _quadrant(min, lo) + hi_quadrant = _quadrant(max, hi) if (lo_quadrant == 2 || lo_quadrant == 3) && hi == 0 # singularity from the left