From f167c50d945a21762e0897bbeb14e1e4d4261299 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Mon, 9 Apr 2018 14:24:27 -0500 Subject: [PATCH 1/4] Speed up trig by a few times --- src/intervals/trigonometric.jl | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/src/intervals/trigonometric.jl b/src/intervals/trigonometric.jl index c8cd07012..c9adcd327 100644 --- a/src/intervals/trigonometric.jl +++ b/src/intervals/trigonometric.jl @@ -1,9 +1,20 @@ # 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}) where {T} = multiply_by_positive_constant(0.5, pi_interval(T)) +half_pi(::Type{T}) where {T} = 0.5 * pi_interval(T) +half_pi(x::T) where {T<:AbstractFloat} = half_pi(T) + +two_pi(::Type{T}) where {T} = multiply_by_positive_constant(2.0, 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 +30,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 From 69d783ee432c313058cd53ade0269899c8ca22c3 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Tue, 10 Apr 2018 11:54:52 -0500 Subject: [PATCH 2/4] Fix sin(x::BigFloat) --- src/intervals/trigonometric.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/intervals/trigonometric.jl b/src/intervals/trigonometric.jl index c9adcd327..390429833 100644 --- a/src/intervals/trigonometric.jl +++ b/src/intervals/trigonometric.jl @@ -10,7 +10,7 @@ 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}) where {T} = multiply_by_positive_constant(0.5, pi_interval(T)) +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) @@ -50,6 +50,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)) From 564d06030a2bc254133bc9839a34ab00ea817877 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Tue, 10 Apr 2018 11:57:00 -0500 Subject: [PATCH 3/4] Fix sin for big interval --- src/intervals/trigonometric.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/intervals/trigonometric.jl b/src/intervals/trigonometric.jl index 390429833..59a2d7dce 100644 --- a/src/intervals/trigonometric.jl +++ b/src/intervals/trigonometric.jl @@ -14,7 +14,8 @@ half_pi(::Type{Float64}) = multiply_by_positive_constant(0.5, pi_interval(Float6 half_pi(::Type{T}) where {T} = 0.5 * pi_interval(T) half_pi(x::T) where {T<:AbstractFloat} = half_pi(T) -two_pi(::Type{T}) where {T} = multiply_by_positive_constant(2.0, pi_interval(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) ) From 776e41996b479622df5465bcdc35324bd71779e3 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Wed, 11 Apr 2018 16:33:16 -0500 Subject: [PATCH 4/4] Revert tan to use previous find_quadrants --- src/intervals/trigonometric.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/intervals/trigonometric.jl b/src/intervals/trigonometric.jl index 59a2d7dce..88cd9d05c 100644 --- a/src/intervals/trigonometric.jl +++ b/src/intervals/trigonometric.jl @@ -131,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)