Fix quadrant

src/intervals/arithmetic/trigonometric.jl

@@ -9,16 +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 # (-π, -π/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, π)
+        r2 ≤ -PI_HI && return 2, 2 # (-π, -π/2)
+        r2 < -PI_LO && return f(2, 3), f(2, 3) # (-π, -π/2) or [-π/2, 0)
+        r2 <  0     && return 3, 3 # [-π/2, 0)
+        r2 ≤  PI_LO && return 0, 0 # [0, π/2)
+        r2 <  PI_HI && return f(0, 1), f(0, 1) # [0, π/2) or [π/2, π)
+        return 1, 1 # [π/2, π)
     else
         k = _unsafe_scale(bareinterval(x) / PI, convert(T, 2))
         fk = floor(k)
-        return f(mod(inf(fk), 4), mod(sup(fk), 4))
+        lfk, hfk = Int(inf(fk)), Int(sup(fk))
+        return f(mod(lfk, 4), mod(hfk, 4)), f(lfk, hfk)
     end
 end
 
@@ -68,11 +69,13 @@ function Base.sin(x::BareInterval{T}) where {T<:AbstractFloat}
 
     lo, hi = bounds(x)
 
-    lo_quadrant = _quadrant(min, lo)
-    hi_quadrant = _quadrant(max, hi)
+    lo_quadrant, lq = _quadrant(min, lo)
+    hi_quadrant, hq = _quadrant(max, hi)
 
-    if lo_quadrant == hi_quadrant
-        d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T))
+    if hq - lq > 3
+        return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 2π
+
+    elseif lo_quadrant == hi_quadrant
         (lo_quadrant == 1) | (lo_quadrant == 2) && return @round(T, sin(hi), sin(lo)) # decreasing
         return @round(T, sin(lo), sin(hi))
 
@@ -170,11 +173,13 @@ function Base.cos(x::BareInterval{T}) where {T<:AbstractFloat}
 
     lo, hi = bounds(x)
 
-    lo_quadrant = _quadrant(min, lo)
-    hi_quadrant = _quadrant(max, hi)
+    lo_quadrant, lq = _quadrant(min, lo)
+    hi_quadrant, hq = _quadrant(max, hi)
 
-    if lo_quadrant == hi_quadrant
-        d ≥ PI_HI && return _unsafe_bareinterval(T, -one(T), one(T))
+    if hq - lq > 3
+        return _unsafe_bareinterval(T, -one(T), one(T)) # diameter ≥ 2π
+
+    elseif lo_quadrant == hi_quadrant
         (lo_quadrant == 2) | (lo_quadrant == 3) && return @round(T, cos(lo), cos(hi)) # increasing
         return @round(T, cos(hi), cos(lo))
 
@@ -272,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)
 
@@ -312,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
@@ -344,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
@@ -382,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