Skip to content
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.

Sign up for GitHub

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

Add Julia port of atan2 from Openlibm and add some tests. #23468

Merged
merged 3 commits into from
Oct 3, 2017
Merged

Add Julia port of atan2 from Openlibm and add some tests. #23468

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
da0f952
Add Julia port of atan2 from Openlibm and add some tests.
pkofod Aug 22, 2017
07cc6fc
Merge branch 'master' into atan2
pkofod Sep 28, 2017
4198f1a
Whitespace fix.
pkofod Sep 29, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 0 additions & 3 deletions base/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -536,9 +536,6 @@ quadrant of the return value.
atan2(y::Real, x::Real) = atan2(promote(float(y),float(x))...)
atan2(y::T, x::T) where {T<:AbstractFloat} = Base.no_op_err("atan2", T)

atan2(y::Float64, x::Float64) = ccall((:atan2,libm), Float64, (Float64, Float64,), y, x)
atan2(y::Float32, x::Float32) = ccall((:atan2f,libm), Float32, (Float32, Float32), y, x)

max(x::T, y::T) where {T<:AbstractFloat} = ifelse((y > x) | (signbit(y) < signbit(x)),
ifelse(isnan(x), x, y), ifelse(isnan(y), y, x))

Expand Down
1 change: 1 addition & 0 deletions base/special/rem_pio2.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ Return the high word of `x` as a `UInt32`.
Return positive part of the high word of `x` as a `UInt32`.
"""
@inline poshighword(x::UInt64) = unsafe_trunc(UInt32,x >> 32)&0x7fffffff
@inline poshighword(x::Float32) = reinterpret(UInt32, x)&0x7fffffff
@inline poshighword(x::Float64) = poshighword(reinterpret(UInt64, x))

@inline function cody_waite_2c_pio2(x::Float64, fn, n)
Expand Down
101 changes: 101 additions & 0 deletions base/special/trig.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,6 +175,107 @@ function asin(x::T) where T<:Union{Float32, Float64}
return asin_kernel(t, x)
end

# atan2 methods
ATAN2_PI_LO(::Type{Float32}) = -8.7422776573f-08
ATAN2_RATIO_BIT_SHIFT(::Type{Float32}) = 23
ATAN2_RATIO_THRESHOLD(::Type{Float32}) = 26

ATAN2_PI_LO(::Type{Float64}) = 1.2246467991473531772E-16
ATAN2_RATIO_BIT_SHIFT(::Type{Float64}) = 20
ATAN2_RATIO_THRESHOLD(::Type{Float64}) = 60

function atan2(y::T, x::T) where T<:Union{Float32, Float64}
# Method :
# M1) Reduce y to positive by atan2(y,x)=-atan2(-y,x).
# M2) Reduce x to positive by (if x and y are unexceptional):
# ARG (x+iy) = arctan(y/x) ... if x > 0,
# ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
#
# Special cases:
#
# S1) ATAN2((anything), NaN ) is NaN;
# S2) ATAN2(NAN , (anything) ) is NaN;
# S3) ATAN2(+-0, +(anything but NaN)) is +-0 ;
# S4) ATAN2(+-0, -(anything but NaN)) is +-pi ;
# S5) ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
# S6) ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
# S7) ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
# S8) ATAN2(+-INF,+INF ) is +-pi/4 ;
# S9) ATAN2(+-INF,-INF ) is +-3pi/4;
# S10) ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
if isnan(x) || isnan(y) # S1 or S2
return T(NaN)
end

if x == T(1.0) # then y/x = y and x > 0, see M2
return atan(y)
end
# generate an m ∈ {0, 1, 2, 3} to branch off of
m = 2*signbit(x) + 1*signbit(y)

if iszero(y)
if m == 0 || m == 1
return y # atan(+-0, +anything) = +-0
elseif m == 2
return T(pi) # atan(+0, -anything) = pi
elseif m == 3
return -T(pi) # atan(-0, -anything) =-pi
end
elseif iszero(x)
return flipsign(T(pi)/2, y)
end

if isinf(x)
if isinf(y)
if m == 0
return T(pi)/4 # atan(+Inf), +Inf))
elseif m == 1
return -T(pi)/4 # atan(-Inf), +Inf))
elseif m == 2
return 3*T(pi)/4 # atan(+Inf, -Inf)
elseif m == 3
return -3*T(pi)/4 # atan(-Inf,-Inf)
end
else
if m == 0
return zero(T) # atan(+...,+Inf) */
elseif m == 1
return -zero(T) # atan(-...,+Inf) */
elseif m == 2
return T(pi) # atan(+...,-Inf) */
elseif m == 3
return -T(pi) # atan(-...,-Inf) */
end
end
end

# x wasn't Inf, but y is
isinf(y) && return copysign(T(pi)/2, y)

ypw = poshighword(y)
xpw = poshighword(x)
# compute y/x for Float32
k = reinterpret(Int32, ypw-xpw)>>ATAN2_RATIO_BIT_SHIFT(T)

if k > ATAN2_RATIO_THRESHOLD(T) # |y/x| > threshold
z=T(pi)/2+T(0.5)*ATAN2_PI_LO(T)
m&=1;
elseif x<0 && k < -ATAN2_RATIO_THRESHOLD(T) # 0 > |y|/x > threshold
z = zero(T)
else #safe to do y/x
z = atan(abs(y/x))
end

if m == 0
return z # atan(+,+)
elseif m == 1
return -z # atan(-,+)
elseif m == 2
return T(pi)-(z-ATAN2_PI_LO(T)) # atan(+,-)
else # default case m == 3
return (z-ATAN2_PI_LO(T))-T(pi) # atan(-,-)
end
end
# acos methods
ACOS_X_MIN_THRESHOLD(::Type{Float32}) = 2.0f0^-26
ACOS_X_MIN_THRESHOLD(::Type{Float64}) = 2.0^-57
Expand Down
54 changes: 54 additions & 0 deletions test/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -704,6 +704,60 @@ end
end
end

@testset "atan2" begin
for T in (Float32, Float64)
@test atan2(T(NaN), T(NaN)) === T(NaN)
@test atan2(T(NaN), T(0.1)) === T(NaN)
@test atan2(T(0.1), T(NaN)) === T(NaN)
r = T(randn())
absr = abs(r)
# y zero
@test atan2(T(r), one(T)) === atan(T(r))
@test atan2(zero(T), absr) === zero(T)
@test atan2(-zero(T), absr) === -zero(T)
@test atan2(zero(T), -absr) === T(pi)
@test atan2(-zero(T), -absr) === -T(pi)
# x zero and y not zero
@test atan2(one(T), zero(T)) === T(pi)/2
@test atan2(-one(T), zero(T)) === -T(pi)/2
# isinf(x) == true && isinf(y) == true
@test atan2(T(Inf), T(Inf)) === T(pi)/4 # m == 0 (see atan2 code)
@test atan2(-T(Inf), T(Inf)) === -T(pi)/4 # m == 1
@test atan2(T(Inf), -T(Inf)) === 3*T(pi)/4 # m == 2
@test atan2(-T(Inf), -T(Inf)) === -3*T(pi)/4 # m == 3
# isinf(x) == true && isinf(y) == false
@test atan2(absr, T(Inf)) === zero(T) # m == 0
@test atan2(-absr, T(Inf)) === -zero(T) # m == 1
@test atan2(absr, -T(Inf)) === T(pi) # m == 2
@test atan2(-absr, -T(Inf)) === -T(pi) # m == 3
# isinf(y) == true && isinf(x) == false
@test atan2(T(Inf), absr) === T(pi)/2
@test atan2(-T(Inf), absr) === -T(pi)/2
@test atan2(T(Inf), -absr) === T(pi)/2
@test atan2(-T(Inf), -absr) === -T(pi)/2
# |y/x| above high threshold
atanpi = T(1.5707963267948966)
@test atan2(T(2.0^61), T(1.0)) === atanpi # m==0
@test atan2(-T(2.0^61), T(1.0)) === -atanpi # m==1
@test atan2(T(2.0^61), -T(1.0)) === atanpi # m==2
@test atan2(-T(2.0^61), -T(1.0)) === -atanpi # m==3
@test atan2(-T(Inf), -absr) === -T(pi)/2
# |y|/x between 0 and low threshold
@test atan2(T(2.0^-61), -T(1.0)) === T(pi) # m==2
@test atan2(-T(2.0^-61), -T(1.0)) === -T(pi) # m==3
# y/x is "safe" ("arbitrary values", just need to hit the branch)
_ATAN2_PI_LO(::Type{Float32}) = -8.7422776573f-08
_ATAN2_PI_LO(::Type{Float64}) = 1.2246467991473531772E-16
@test atan2(T(5.0), T(2.5)) === atan(abs(T(5.0)/T(2.5)))
@test atan2(-T(5.0), T(2.5)) === -atan(abs(-T(5.0)/T(2.5)))
@test atan2(T(5.0), -T(2.5)) === T(pi)-(atan(abs(T(5.0)/-T(2.5)))-_ATAN2_PI_LO(T))
@test atan2(-T(5.0), -T(2.5)) === -(T(pi)-atan(abs(-T(5.0)/-T(2.5)))-_ATAN2_PI_LO(T))
@test atan2(T(1235.2341234), T(2.5)) === atan(abs(T(1235.2341234)/T(2.5)))
@test atan2(-T(1235.2341234), T(2.5)) === -atan(abs(-T(1235.2341234)/T(2.5)))
@test atan2(T(1235.2341234), -T(2.5)) === T(pi)-(atan(abs(T(1235.2341234)/-T(2.5)))-_ATAN2_PI_LO(T))
@test atan2(-T(1235.2341234), -T(2.5)) === -(T(pi)-(atan(abs(-T(1235.2341234)/T(2.5)))-_ATAN2_PI_LO(T)))
end
end
@testset "acos #23283" begin
for T in (Float32, Float64)
@test acos(zero(T)) === T(pi)/2
Expand Down