diff --git a/base/exports.jl b/base/exports.jl index 316025db9ce6c..3b49e3b96e622 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -329,6 +329,7 @@ export sinc, sincos, sincosd, + sincospi, sind, sinh, sinpi, diff --git a/base/math.jl b/base/math.jl index ce9780269ba1c..515f634bee383 100644 --- a/base/math.jl +++ b/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, diff --git a/base/special/trig.jl b/base/special/trig.jl index c848ebfe492d1..0ebf2c7f6c358 100644 --- a/base/special/trig.jl +++ b/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::Integer) = x >= 0 ? zero(float(x)) : -zero(float(x)) cospi(x::Integer) = isodd(x) ? -one(float(x)) : one(float(x)) +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 diff --git a/test/math.jl b/test/math.jl index 79deaac247a7e..d36883c344680 100644 --- a/test/math.jl +++ b/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