Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Use sincos from libm in cis
Browse files Browse the repository at this point in the history
yuyichao committed Apr 27, 2017

Verified

This commit was signed with the committer’s verified signature.
yuyichao Yichao Yu
1 parent 19e3f1b commit d5a3d02
Showing 2 changed files with 41 additions and 3 deletions.
10 changes: 9 additions & 1 deletion base/complex.jl
Original file line number Diff line number Diff line change
@@ -420,8 +420,16 @@ sqrt(z::Complex) = sqrt(float(z))
# return Complex(abs(iz)/r/2, copysign(r,iz))
# end

function cis(v::Union{Float64,Float32})
res = FastMath.cis_fast(v) # Note: not yet defined at this point
return Math.nan_dom_err(res, v)
end

# compute exp(im*theta)
cis(theta::Real) = Complex(cos(theta),sin(theta))
cis(theta::AbstractFloat) = Complex(cos(theta),sin(theta))

# Convert integer/irrational to float so that they can also use the `sincos` path
cis(theta::Real) = cis(float(theta))

"""
cis(z)
34 changes: 32 additions & 2 deletions base/fastmath.jl
Original file line number Diff line number Diff line change
@@ -273,6 +273,38 @@ atan2_fast(x::Float64, y::Float64) =

# explicit implementations

# FIXME: Change to `ccall((:sincos, libm))` when `Ref` calling convention can be
# stack allocated.
@inline function cis_fast(v::Float64)
s, c = Base.llvmcall("""
%f = bitcast i8 *%1 to void (double, double *, double *)*
%pres = alloca [2 x double]
%p1 = getelementptr inbounds [2 x double], [2 x double]* %pres, i64 0, i64 0
%p2 = getelementptr inbounds [2 x double], [2 x double]* %pres, i64 0, i64 1
call void %f(double %0, double *nocapture noalias %p1, double *nocapture noalias %p2)
%res = load [2 x double], [2 x double]* %pres
ret [2 x double] %res
""", Tuple{Float64,Float64}, Tuple{Float64,Ptr{Void}}, v, cglobal((:sincos, libm)))
return Complex(c, s)
end

@inline function cis_fast(v::Float32)
s, c = Base.llvmcall("""
%f = bitcast i8 *%1 to void (float, float *, float *)*
%pres = alloca [2 x float]
%p1 = getelementptr inbounds [2 x float], [2 x float]* %pres, i64 0, i64 0
%p2 = getelementptr inbounds [2 x float], [2 x float]* %pres, i64 0, i64 1
call void %f(float %0, float *nocapture noalias %p1, float *nocapture noalias %p2)
%res = load [2 x float], [2 x float]* %pres
ret [2 x float] %res
""", Tuple{Float32,Float32}, Tuple{Float32,Ptr{Void}}, v, cglobal((:sincosf, libm)))
return Complex(c, s)
end

# Fallback for `::AbstractFloat` that's more specific than `::Real`
cis_fast(x::AbstractFloat) = cis(x)
cis_fast(v::Real) = cis_fast(float(v)::AbstractFloat)

@fastmath begin
exp10_fast(x::T) where {T<:FloatTypes} = exp2(log2(T(10))*x)
exp10_fast(x::Integer) = exp10(float(x))
@@ -287,8 +319,6 @@ atan2_fast(x::Float64, y::Float64) =

# complex numbers

cis_fast(x::T) where {T<:FloatTypes} = Complex{T}(cos(x), sin(x))

# See <http://en.cppreference.com/w/cpp/numeric/complex>
pow_fast(x::T, y::T) where {T<:ComplexTypes} = exp(y*log(x))
pow_fast(x::T, y::Complex{T}) where {T<:FloatTypes} = exp(y*log(x))

0 comments on commit d5a3d02

Please sign in to comment.