From 24086de7a606441af1a4db00f3b0524e9c1a434b Mon Sep 17 00:00:00 2001
From: Yichao Yu <yyc1992@gmail.com>
Date: Sat, 6 May 2017 12:23:54 -0400
Subject: [PATCH] Define `sincos`

Fix #10442
---
 base/exports.jl        |  1 +
 base/fastmath.jl       | 40 ++++++++++++++++++++++++++++++++++++++++
 base/math.jl           | 15 ++++++++++++++-
 base/mpfr.jl           | 13 ++++++++++++-
 base/sysimg.jl         |  8 ++++----
 doc/src/stdlib/math.md |  1 +
 test/fastmath.jl       |  1 +
 test/math.jl           |  9 +++++++++
 8 files changed, 82 insertions(+), 6 deletions(-)

diff --git a/base/exports.jl b/base/exports.jl
index 039c5328f9bb1..ce570f899b3de 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -405,6 +405,7 @@ export
     significand,
     sin,
     sinc,
+    sincos,
     sind,
     sinh,
     sinpi,
diff --git a/base/fastmath.jl b/base/fastmath.jl
index fbef3beb0f836..d05a80131d971 100644
--- a/base/fastmath.jl
+++ b/base/fastmath.jl
@@ -76,6 +76,7 @@ const fast_op =
          :min => :min_fast,
          :minmax => :minmax_fast,
          :sin => :sin_fast,
+         :sincos => :sincos_fast,
          :sinh => :sinh_fast,
          :sqrt => :sqrt_fast,
          :tan => :tan_fast,
@@ -273,6 +274,45 @@ atan2_fast(x::Float64, y::Float64) =
 
 # explicit implementations
 
+# FIXME: Change to `ccall((:sincos, libm))` when `Ref` calling convention can be
+#        stack allocated.
+@inline function sincos_fast(v::Float64)
+    return Base.llvmcall("""
+    %f = bitcast i8 *%1 to void (double, double *, double *)*
+    %ps = alloca double
+    %pc = alloca double
+    call void %f(double %0, double *%ps, double *%pc)
+    %s = load double, double* %ps
+    %c = load double, double* %pc
+    %res0 = insertvalue [2 x double] undef, double %s, 0
+    %res = insertvalue [2 x double] %res0, double %c, 1
+    ret [2 x double] %res
+    """, Tuple{Float64,Float64}, Tuple{Float64,Ptr{Void}}, v, cglobal((:sincos, libm)))
+end
+
+@inline function sincos_fast(v::Float32)
+    return Base.llvmcall("""
+    %f = bitcast i8 *%1 to void (float, float *, float *)*
+    %ps = alloca float
+    %pc = alloca float
+    call void %f(float %0, float *%ps, float *%pc)
+    %s = load float, float* %ps
+    %c = load float, float* %pc
+    %res0 = insertvalue [2 x float] undef, float %s, 0
+    %res = insertvalue [2 x float] %res0, float %c, 1
+    ret [2 x float] %res
+    """, Tuple{Float32,Float32}, Tuple{Float32,Ptr{Void}}, v, cglobal((:sincosf, libm)))
+end
+
+@inline function sincos_fast(v::Float16)
+    s, c = sincos_fast(Float32(v))
+    return Float16(s), Float16(c)
+end
+
+sincos_fast(v::AbstractFloat) = (sin_fast(v), cos_fast(v))
+sincos_fast(v::Real) = sincos_fast(float(v)::AbstractFloat)
+sincos_fast(v) = (sin_fast(v), cos_fast(v))
+
 @fastmath begin
     exp10_fast(x::T) where {T<:FloatTypes} = exp2(log2(T(10))*x)
     exp10_fast(x::Integer) = exp10(float(x))
diff --git a/base/math.jl b/base/math.jl
index 4d8eb64fc8d32..b5f691990574b 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -2,7 +2,7 @@
 
 module Math
 
-export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
+export sin, cos, sincos, tan, sinh, cosh, tanh, asin, acos, atan,
        asinh, acosh, atanh, sec, csc, cot, asec, acsc, acot,
        sech, csch, coth, asech, acsch, acoth,
        sinpi, cospi, sinc, cosc,
@@ -419,6 +419,19 @@ for f in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10,
     end
 end
 
+"""
+    sincos(x)
+
+Compute sine and cosine of `x`, where `x` is in radians.
+"""
+@inline function sincos(x)
+    res = Base.FastMath.sincos_fast(x)
+    if (isnan(res[1]) | isnan(res[2])) & !isnan(x)
+        throw(DomainError())
+    end
+    return res
+end
+
 sqrt(x::Float64) = sqrt_llvm(x)
 sqrt(x::Float32) = sqrt_llvm(x)
 
diff --git a/base/mpfr.jl b/base/mpfr.jl
index 5ff7007370d10..4855ca4909678 100644
--- a/base/mpfr.jl
+++ b/base/mpfr.jl
@@ -13,7 +13,7 @@ import
         nextfloat, prevfloat, promote_rule, rem, rem2pi, round, show, float,
         sum, sqrt, string, print, trunc, precision, exp10, expm1,
         gamma, lgamma, log1p,
-        eps, signbit, sin, cos, tan, sec, csc, cot, acos, asin, atan,
+        eps, signbit, sin, cos, sincos, tan, sec, csc, cot, acos, asin, atan,
         cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, atan2,
         cbrt, typemax, typemin, unsafe_trunc, realmin, realmax, rounding,
         setrounding, maxintfloat, widen, significand, frexp, tryparse, iszero, big
@@ -24,6 +24,8 @@ import Base.GMP: ClongMax, CulongMax, CdoubleMax, Limb
 
 import Base.Math.lgamma_r
 
+import Base.FastMath.sincos_fast
+
 function __init__()
     try
         # set exponent to full range by default
@@ -515,6 +517,15 @@ for f in (:exp, :exp2, :exp10, :expm1, :cosh, :sinh, :tanh, :sech, :csch, :coth,
     end
 end
 
+function sincos_fast(v::BigFloat)
+    s = BigFloat()
+    c = BigFloat()
+    ccall((:mpfr_sin_cos, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32),
+          &s, &c, &v, ROUNDING_MODE[])
+    return (s, c)
+end
+sincos(v::BigFloat) = sincos_fast(v)
+
 # return log(2)
 function big_ln2()
     c = BigFloat()
diff --git a/base/sysimg.jl b/base/sysimg.jl
index acfb4c7b68ee4..03b0e3e3b457e 100644
--- a/base/sysimg.jl
+++ b/base/sysimg.jl
@@ -258,6 +258,10 @@ importall .Order
 include("sort.jl")
 importall .Sort
 
+# Fast math
+include("fastmath.jl")
+importall .FastMath
+
 function deepcopy_internal end
 
 # BigInts and BigFloats
@@ -343,10 +347,6 @@ importall .DFT
 include("dsp.jl")
 importall .DSP
 
-# Fast math
-include("fastmath.jl")
-importall .FastMath
-
 # libgit2 support
 include("libgit2/libgit2.jl")
 
diff --git a/doc/src/stdlib/math.md b/doc/src/stdlib/math.md
index a195a354e9a0f..2e106fe4952ae 100644
--- a/doc/src/stdlib/math.md
+++ b/doc/src/stdlib/math.md
@@ -58,6 +58,7 @@ Base.:(!)
 Base.isapprox
 Base.sin
 Base.cos
+Base.sincos
 Base.tan
 Base.Math.sind
 Base.Math.cosd
diff --git a/test/fastmath.jl b/test/fastmath.jl
index 1a89a83503eda..80120e7c7b2fd 100644
--- a/test/fastmath.jl
+++ b/test/fastmath.jl
@@ -9,6 +9,7 @@
 @test macroexpand(:(@fastmath min(1))) == :(Base.FastMath.min_fast(1))
 @test macroexpand(:(@fastmath min)) == :(Base.FastMath.min_fast)
 @test macroexpand(:(@fastmath x.min)) == :(x.min)
+@test macroexpand(:(@fastmath sincos(x))) == :(Base.FastMath.sincos_fast(x))
 
 # basic arithmetic
 
diff --git a/test/math.jl b/test/math.jl
index 136c50f98463c..59cb22f12446d 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -627,3 +627,12 @@ end
 @testset "promote Float16 irrational #15359" begin
     @test typeof(Float16(.5) * pi) == Float16
 end
+
+@testset "sincos" begin
+    @test sincos(1.0) === (sin(1.0), cos(1.0))
+    @test sincos(1f0) === (sin(1f0), cos(1f0))
+    @test sincos(Float16(1)) === (sin(Float16(1)), cos(Float16(1)))
+    @test sincos(1) === (sin(1), cos(1))
+    @test sincos(big(1)) == (sin(big(1)), cos(big(1)))
+    @test sincos(big(1.0)) == (sin(big(1.0)), cos(big(1.0)))
+end