diff --git a/src/intervals/trigonometric.jl b/src/intervals/trigonometric.jl index c8cd07012..88cd9d05c 100644 --- a/src/intervals/trigonometric.jl +++ b/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)