Speed up sin and cos by a factor 6.5 (#117)

* Speed up trig by a few times

* Fix sin(x::BigFloat)

* Fix sin for big interval

* Revert tan to use previous find_quadrants

Author: dpsanders

src/intervals/trigonometric.jl

@@ -1,9 +1,21 @@
 # This file is part of the IntervalArithmetic.jl package; MIT licensed
 
-half_pi(::Type{T}) where T = pi_interval(T) / 2
-half_pi(x::T) where T<:AbstractFloat = half_pi(T)
 
-two_pi(::Type{T})  where T = pi_interval(T) * 2
+# hard code this for efficiency:
+const one_over_half_pi_interval = inv(0.5 * pi_interval(Float64))
+
+"""
+Multiply an interval by a positive constant.
+For efficiency, does not check that the constant is positive.
+"""
+multiply_by_positive_constant(α, x::Interval) = @round(α*x.lo, α*x.hi)
+
+half_pi(::Type{Float64}) = multiply_by_positive_constant(0.5, pi_interval(Float64))
+half_pi(::Type{T}) where {T} = 0.5 * pi_interval(T)
+half_pi(x::T) where {T<:AbstractFloat} = half_pi(T)
+
+two_pi(::Type{Float64}) = multiply_by_positive_constant(2.0, pi_interval(Float64))
+two_pi(::Type{T}) where {T} = 2 * pi_interval(T)
 
 range_atan2(::Type{T}) where {T<:Real} = Interval(-(pi_interval(T).hi), pi_interval(T).hi)
 half_range_atan2(::Type{T}) where {T} = (temp = half_pi(T); Interval(-(temp.hi), temp.hi) )
@@ -19,9 +31,17 @@ This is a rather indirect way to determine if π/2 and 3π/2 are contained
 in the interval; cf. the formula for sine of an interval in
 Tucker, *Validated Numerics*."""
 
-function find_quadrants(x::T) where {T<:AbstractFloat}
+function find_quadrants(x::T) where {T}
     temp = atomic(Interval{T}, x) / half_pi(x)
-    (floor(temp.lo), floor(temp.hi))
+
+    return SVector(floor(temp.lo), floor(temp.hi))
+end
+
+function find_quadrants(x::Float64)
+    temp = multiply_by_positive_constant(x, one_over_half_pi_interval)
+    # x / half_pi(Float64)
+
+    return SVector(floor(temp.lo), floor(temp.hi))
 end
 
 function sin(a::Interval{T}) where T
@@ -31,6 +51,8 @@ function sin(a::Interval{T}) where T
 
     diam(a) > two_pi(T).lo && return whole_range
 
+    # The following is equiavlent to doing temp = a / half_pi  and
+    # taking floor(a.lo), floor(a.hi)
     lo_quadrant = minimum(find_quadrants(a.lo))
     hi_quadrant = maximum(find_quadrants(a.hi))
 
@@ -109,14 +131,19 @@ function cos(a::Interval{T}) where T
     end
 end
 
+function find_quadrants_tan(x::T) where {T}
+    temp = atomic(Interval{T}, x) / half_pi(x)
+
+    return SVector(floor(temp.lo), floor(temp.hi))
+end
 
 function tan(a::Interval{T}) where T
     isempty(a) && return a
 
     diam(a) > pi_interval(T).lo && return entireinterval(a)
 
-    lo_quadrant = minimum(find_quadrants(a.lo))
-    hi_quadrant = maximum(find_quadrants(a.hi))
+    lo_quadrant = minimum(find_quadrants_tan(a.lo))
+    hi_quadrant = maximum(find_quadrants_tan(a.hi))
 
     lo_quadrant_mod = mod(lo_quadrant, 2)
     hi_quadrant_mod = mod(hi_quadrant, 2)