add sincospi (JuliaLang#35816)

Author: simeonschaub

base/exports.jl

@@ -329,6 +329,7 @@ export
     sinc,
     sincos,
     sincosd,
+    sincospi,
     sind,
     sinh,
     sinpi,

base/math.jl

@@ -5,7 +5,7 @@ module Math
 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,
+       sinpi, cospi, sincospi, sinc, cosc,
        cosd, cotd, cscd, secd, sind, tand, sincosd,
        acosd, acotd, acscd, asecd, asind, atand,
        rad2deg, deg2rad,

base/special/trig.jl

@@ -777,8 +777,8 @@ function sinpi(x::T) where T<:AbstractFloat
     end
 end
 
-# Integers and Rationals
-function sinpi(x::T) where T<:Union{Integer,Rational}
+# Rationals
+function sinpi(x::T) where T<:Rational
     Tf = float(T)
     if !isfinite(x)
         throw(DomainError(x, "`x` must be finite."))
@@ -836,8 +836,8 @@ function cospi(x::T) where T<:AbstractFloat
     end
 end
 
-# Integers and Rationals
-function cospi(x::T) where T<:Union{Integer,Rational}
+# Rationals
+function cospi(x::T) where T<:Rational
     if !isfinite(x)
         throw(DomainError(x, "`x` must be finite."))
     end
@@ -860,10 +860,79 @@ function cospi(x::T) where T<:Union{Integer,Rational}
     end
 end
 
+"""
+    sincospi(x)
+
+Simultaneously compute `sinpi(x)` and `cospi(x)`, where the `x` is in radians.
+"""
+function sincospi(x::T) where T<:AbstractFloat
+    if !isfinite(x)
+        isnan(x) && return x, x
+        throw(DomainError(x, "`x` cannot be infinite."))
+    end
+
+    ax = abs(x)
+    s = maxintfloat(T)
+    ax >= s && return (copysign(zero(T), x), one(T)) # even integer-valued
+
+    # reduce to interval [-1,1]
+    # assumes RoundNearest rounding mode
+    t = 3*(s/2)
+    rx = x-((x+t)-t) # zeros may be incorrectly signed
+    arx = abs(rx)
+
+    # same selection scheme as sinpi and cospi
+    if (arx == 0) | (arx == 1)
+        return copysign(zero(T), x), ifelse(ax % 2 == 0, one(T), -one(T))
+    elseif arx < 0.25
+        return sincos_kernel(mulpi_ext(rx))
+    elseif arx < 0.75
+        y = mulpi_ext(T(0.5) - arx)
+        return copysign(cos_kernel(y), rx), sin_kernel(y)
+    else
+        y_si = mulpi_ext(copysign(one(T), rx) - rx)
+        y_co = mulpi_ext(one(T) - arx)
+        return sin_kernel(y_si), -cos_kernel(y_co)
+    end
+end
+
+# Rationals
+function sincospi(x::T) where T<:Rational
+    Tf = float(T)
+    if !isfinite(x)
+        throw(DomainError(x, "`x` must be finite."))
+    end
+
+    # until we get an IEEE remainder function (#9283)
+    rx = rem(x,2)
+    if rx > 1
+        rx -= 2
+    elseif rx < -1
+        rx += 2
+    end
+    arx = abs(rx)
+
+    # same selection scheme as sinpi and cospi
+    if (arx == 0) | (arx == 1)
+        return copysign(zero(Tf),x), ifelse(iseven(numerator(x)), one(Tf), -one(Tf))
+    elseif arx < 0.25
+        return sincos_kernel(mulpi_ext(rx))
+    elseif arx < 0.75
+        y = mulpi_ext(T(0.5) - arx)
+        return copysign(cos_kernel(y), rx), sin_kernel(y)
+    else
+        y_si = mulpi_ext(copysign(one(T), rx) - rx)
+        y_co = mulpi_ext(one(T) - arx)
+        return sin_kernel(y_si), -cos_kernel(y_co)
+    end
+end
+
 sinpi(x::T) where {T<:Integer} = zero(signed(T))
 cospi(x::T) where {T<:Integer} = isodd(x) ? -one(signed(T)) : one(signed(T))
+sincospi(x::Integer) = (sinpi(x), cospi(x))
 sinpi(x::Real) = sinpi(float(x))
 cospi(x::Real) = cospi(float(x))
+sincospi(x::Real) = sincospi(float(x))
 
 function sinpi(z::Complex{T}) where T
     F = float(T)
@@ -893,7 +962,8 @@ function sinpi(z::Complex{T}) where T
         end
     else
         pizi = pi*zi
-        Complex(sinpi(zr)*cosh(pizi), cospi(zr)*sinh(pizi))
+        sipi, copi = sincospi(zr)
+        Complex(sipi*cosh(pizi), copi*sinh(pizi))
     end
 end
 
@@ -928,7 +998,55 @@ function cospi(z::Complex{T}) where T
         end
     else
         pizi = pi*zi
-        Complex(cospi(zr)*cosh(pizi), -sinpi(zr)*sinh(pizi))
+        sipi, copi = sincospi(zr)
+        Complex(copi*cosh(pizi), -sipi*sinh(pizi))
+    end
+end
+
+function sincospi(z::Complex{T}) where T
+    F = float(T)
+    zr, zi = reim(z)
+    if isinteger(zr)
+        # zr = ...,-2,-1,0,1,2,...
+        # sin(pi*zr) == ±0
+        # cos(pi*zr) == ±1
+        # cosh(pi*zi) > 0
+        s = copysign(zero(F),zr)
+        c_pos = isa(zr,Integer) ? iseven(zr) : isinteger(zr/2)
+        pizi = pi*zi
+        sh, ch = sinh(pizi), cosh(pizi)
+        (
+            Complex(s, c_pos ? sh : -sh),
+            Complex(c_pos ? ch : -ch, isnan(zi) ? s : -flipsign(s,zi)),
+        )
+    elseif isinteger(2*zr)
+        # zr = ...,-1.5,-0.5,0.5,1.5,2.5,...
+        # sin(pi*zr) == ±1
+        # cos(pi*zr) == +0
+        # sign(sinh(pi*zi)) == sign(zi)
+        s_pos = isinteger((2*zr-1)/4)
+        pizi = pi*zi
+        sh, ch = sinh(pizi), cosh(pizi)
+        (
+            Complex(s_pos ? ch : -ch, isnan(zi) ? zero(F) : copysign(zero(F),zi)),
+            Complex(zero(F), s_pos ? -sh : sh),
+        )
+    elseif !isfinite(zr)
+        if zi == 0
+            Complex(F(NaN), F(zi)), Complex(F(NaN), isnan(zr) ? zero(F) : -flipsign(F(zi),zr))
+        elseif isinf(zi)
+            Complex(F(NaN), F(zi)), Complex(F(Inf), F(NaN))
+        else
+            Complex(F(NaN), F(NaN)), Complex(F(NaN), F(NaN))
+        end
+    else
+        pizi = pi*zi
+        sipi, copi = sincospi(zr)
+        sihpi, cohpi = sinh(pizi), cosh(pizi)
+        (
+            Complex(sipi*cohpi, copi*sihpi),
+            Complex(copi*cohpi, -sipi*sihpi),
+        )
     end
 end
 

test/math.jl

@@ -408,8 +408,11 @@ end
             @test sincosd(convert(T, 270))::fTsc === (  -one(fT), zero(fT) )
         end
 
-        @testset "sinpi and cospi" begin
-            for x = -3:0.3:3
+        @testset "$name" for (name, (sinpi, cospi)) in (
+            "sinpi and cospi" => (sinpi, cospi),
+            "sincospi" => (x->sincospi(x)[1], x->sincospi(x)[2])
+        )
+            @testset "pi * $x" for x = -3:0.3:3
                 @test sinpi(convert(T,x))::fT ≈ convert(fT,sin(pi*x)) atol=eps(pi*convert(fT,x))
                 @test cospi(convert(T,x))::fT ≈ convert(fT,cos(pi*x)) atol=eps(pi*convert(fT,x))
             end
@@ -433,10 +436,13 @@ end
             @test cosd(convert(T,60)) == 0.5
             @test sind(convert(T,150)) == 0.5
             @test sinpi(one(T)/convert(T,6)) == 0.5
+            @test sincospi(one(T)/convert(T,6))[1] == 0.5
             @test_throws DomainError sind(convert(T,Inf))
             @test_throws DomainError cosd(convert(T,Inf))
             T != Float32 && @test cospi(one(T)/convert(T,3)) == 0.5
+            T != Float32 && @test sincospi(one(T)/convert(T,3))[2] == 0.5
             T == Rational{Int} && @test sinpi(5//6) == 0.5
+            T == Rational{Int} && @test sincospi(5//6)[1] == 0.5
         end
     end
     scdm = sincosd(missing)
@@ -445,10 +451,12 @@ end
 end
 
 @testset "Integer args to sinpi/cospi/sinc/cosc" begin
-    @test sinpi(1) == 0
-    @test sinpi(-1) == -0
-    @test cospi(1) == -1
-    @test cospi(2) == 1
+    for (sinpi, cospi) in ((sinpi, cospi), (x->sincospi(x)[1], x->sincospi(x)[2]))
+        @test sinpi(1) == 0
+        @test sinpi(-1) == -0
+        @test cospi(1) == -1
+        @test cospi(2) == 1
+    end
 
     @test sinc(1) == 0
     @test sinc(complex(1,0)) == 0
@@ -462,20 +470,25 @@ end
 
 @testset "Irrational args to sinpi/cospi/sinc/cosc" begin
     for x in (pi, ℯ, Base.MathConstants.golden)
-        @test sinpi(x) ≈ Float64(sinpi(big(x)))
-        @test cospi(x) ≈ Float64(cospi(big(x)))
+        for (sinpi, cospi) in ((sinpi, cospi), (x->sincospi(x)[1], x->sincospi(x)[2]))
+            @test sinpi(x) ≈ Float64(sinpi(big(x)))
+            @test cospi(x) ≈ Float64(cospi(big(x)))
+            @test sinpi(complex(x, x)) ≈ Complex{Float64}(sinpi(complex(big(x), big(x))))
+            @test cospi(complex(x, x)) ≈ Complex{Float64}(cospi(complex(big(x), big(x))))
+        end
         @test sinc(x)  ≈ Float64(sinc(big(x)))
         @test cosc(x)  ≈ Float64(cosc(big(x)))
-        @test sinpi(complex(x, x)) ≈ Complex{Float64}(sinpi(complex(big(x), big(x))))
-        @test cospi(complex(x, x)) ≈ Complex{Float64}(cospi(complex(big(x), big(x))))
         @test sinc(complex(x, x))  ≈ Complex{Float64}(sinc(complex(big(x),  big(x))))
         @test cosc(complex(x, x))  ≈ Complex{Float64}(cosc(complex(big(x),  big(x))))
     end
 end
 
 @testset "trig function type stability" begin
-    @testset "$T $f" for T = (Float32,Float64,BigFloat), f = (sind,cosd,sinpi,cospi)
-        @test Base.return_types(f,Tuple{T}) == [T]
+    @testset "$T $f" for T = (Float32,Float64,BigFloat,Rational{Int16},Complex{Int32},ComplexF16), f = (sind,cosd,sinpi,cospi)
+        @test Base.return_types(f,Tuple{T}) == [float(T)]
+    end
+    @testset "$T sincospi" for T = (Float32,Float64,BigFloat,Rational{Int16},Complex{Int32},ComplexF16)
+        @test Base.return_types(sincospi,Tuple{T}) == [Tuple{float(T),float(T)}]
     end
 end