Skip to content

Commit

Permalink
Use sincos from libm in cis
Browse files Browse the repository at this point in the history
  • Loading branch information
yuyichao committed Apr 27, 2017
1 parent b838f2e commit bd84708
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 1 deletion.
39 changes: 38 additions & 1 deletion base/complex.jl
Original file line number Diff line number Diff line change
Expand Up @@ -420,8 +420,45 @@ sqrt(z::Complex) = sqrt(float(z))
# return Complex(abs(iz)/r/2, copysign(r,iz))
# end

# FIXME: Change to `ccall((:sincos, libm_name))` when `Ref` calling convention can be
# stack allocated. Also in `fastmath.jl`
function cis(v::Float64)
s, c = 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_name)))
if !isnan(v) && (isnan(s) || isnan(c))
throw(DomainError())
end
return Complex(c, s)
end

function cis(v::Float32)
s, c = 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_name)))
if !isnan(v) && (isnan(s) || isnan(c))
throw(DomainError())
end
return Complex(c, s)
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)
Expand Down
30 changes: 30 additions & 0 deletions base/fastmath.jl
Original file line number Diff line number Diff line change
Expand Up @@ -273,6 +273,36 @@ atan2_fast(x::Float64, y::Float64) =

# explicit implementations

# FIXME: Change to `ccall((:sincos, libm_name))` when `Ref` calling convention can be
# stack allocated. Also in `complex.jl`
@inline function cis_fast(v::Float64)
s, c = 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_name)))
return Complex(c, s)
end

@inline function cis_fast(v::Float32)
s, c = 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_name)))
return Complex(c, s)
end

cis_fast(v::Real) = cis_fast(float(v))

@fastmath begin
exp10_fast(x::T) where {T<:FloatTypes} = exp2(log2(T(10))*x)
exp10_fast(x::Integer) = exp10(float(x))
Expand Down

0 comments on commit bd84708

Please sign in to comment.