-
Notifications
You must be signed in to change notification settings - Fork 72
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Speed up sin and cos by a factor 6.5 #117
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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)) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. What about returning a tuple here, instead of a There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It used to be a tuple; changing it to an In any case all of this code can be simplified a lot when #119 is possible, I believe. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Fair enough. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. OK, I found an alternative solution: return a tuple and do
Too late... |
||
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) | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Perhaps instead of defining
multiply_by_positive_constant
we should rather define the methods*(a::T, x::Interval)
and*(x::Interval, a::T)
(inarithmetic.jl
), whereT<:Real
(or maybe<:Number
), or even unspecified, as it is here. Currently we use promotion to an interval and rounding, which may slow down things.Question: Below you use this function for to define
half_pi(::Type{Float64})
, but not for other types; can you explain why?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I use it for
two_pi
as well.Actually the
half_pi
function is no longer used, I believe.We may well want to define explicit methods for
*(constant, Interval)
. I've opened an issue: #120There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That was actually the suggestion.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Right, but I'd prefer to leave it for a separate PR.
Also, if you actually know that the constant is positive then this specialised version is even slightly quicker.