Implement tanpi (#48575)

See issue: #48226

base/exports.jl

@@ -352,6 +352,7 @@ export
     tan,
     tand,
     tanh,
+    tanpi,
     trailing_ones,
     trailing_zeros,
     trunc,

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

base/mpfr.jl

@@ -17,7 +17,7 @@ import
         cbrt, typemax, typemin, unsafe_trunc, floatmin, floatmax, rounding,
         setrounding, maxintfloat, widen, significand, frexp, tryparse, iszero,
         isone, big, _string_n, decompose, minmax,
-        sinpi, cospi, sincospi, sind, cosd, tand, asind, acosd, atand
+        sinpi, cospi, sincospi, tanpi, sind, cosd, tand, asind, acosd, atand
 
 import ..Rounding: rounding_raw, setrounding_raw
 
@@ -790,7 +790,7 @@ function sum(arr::AbstractArray{BigFloat})
 end
 
 # Functions for which NaN results are converted to DomainError, following Base
-for f in (:sin, :cos, :tan, :sec, :csc, :acos, :asin, :atan, :acosh, :asinh, :atanh, :sinpi, :cospi)
+for f in (:sin, :cos, :tan, :sec, :csc, :acos, :asin, :atan, :acosh, :asinh, :atanh, :sinpi, :cospi, :tanpi)
     @eval begin
         function ($f)(x::BigFloat)
             isnan(x) && return x

base/special/trig.jl

@@ -725,43 +725,61 @@ end
 
 # Uses minimax polynomial of sin(π * x) for π * x in [0, .25]
 @inline function sinpi_kernel(x::Float64)
+    sinpi_kernel_wide(x)
+end
+@inline function sinpi_kernel_wide(x::Float64)
     x² = x*x
     x⁴ = x²*x²
     r  = evalpoly(x², (2.5501640398773415, -0.5992645293202981, 0.08214588658006512,
-                      -7.370429884921779e-3, 4.662827319453555e-4, -2.1717412523382308e-5))
+                       -7.370429884921779e-3, 4.662827319453555e-4, -2.1717412523382308e-5))
     return muladd(3.141592653589793, x, x*muladd(-5.16771278004997,
                   x², muladd(x⁴, r,  1.2245907532225998e-16)))
 end
 @inline function sinpi_kernel(x::Float32)
+    Float32(sinpi_kernel_wide(x))
+end
+@inline function sinpi_kernel_wide(x::Float32)
     x = Float64(x)
-    return Float32(x*evalpoly(x*x, (3.1415926535762266, -5.167712769188119,
-                                    2.5501626483206374, -0.5992021090314925, 0.08100185277841528)))
+    return x*evalpoly(x*x, (3.1415926535762266, -5.167712769188119,
+                            2.5501626483206374, -0.5992021090314925, 0.08100185277841528))
 end
 
 @inline function sinpi_kernel(x::Float16)
+    Float16(sinpi_kernel_wide(x))
+end
+@inline function sinpi_kernel_wide(x::Float16)
     x = Float32(x)
-    return Float16(x*evalpoly(x*x, (3.1415927f0, -5.1677127f0, 2.5501626f0, -0.5992021f0, 0.081001855f0)))
+    return x*evalpoly(x*x, (3.1415927f0, -5.1677127f0, 2.5501626f0, -0.5992021f0, 0.081001855f0))
 end
 
 # Uses minimax polynomial of cos(π * x) for π * x in [0, .25]
 @inline function cospi_kernel(x::Float64)
+    cospi_kernel_wide(x)
+end
+@inline function cospi_kernel_wide(x::Float64)
     x² = x*x
     r = x²*evalpoly(x², (4.058712126416765, -1.3352627688537357, 0.23533063027900392,
-                        -0.025806887811869204, 1.9294917136379183e-3, -1.0368935675474665e-4))
+                         -0.025806887811869204, 1.9294917136379183e-3, -1.0368935675474665e-4))
     a_x² = 4.934802200544679 * x²
     a_x²lo = muladd(3.109686485461973e-16, x², muladd(4.934802200544679, x², -a_x²))
 
     w  = 1.0-a_x²
     return w + muladd(x², r, ((1.0-w)-a_x²) - a_x²lo)
 end
 @inline function cospi_kernel(x::Float32)
+    Float32(cospi_kernel_wide(x))
+end
+@inline function cospi_kernel_wide(x::Float32)
     x = Float64(x)
-    return Float32(evalpoly(x*x, (1.0, -4.934802200541122, 4.058712123568637,
-                                 -1.3352624040152927, 0.23531426791507182, -0.02550710082498761)))
+    return evalpoly(x*x, (1.0, -4.934802200541122, 4.058712123568637,
+                          -1.3352624040152927, 0.23531426791507182, -0.02550710082498761))
 end
 @inline function cospi_kernel(x::Float16)
+    Float16(cospi_kernel_wide(x))
+end
+@inline function cospi_kernel_wide(x::Float16)
     x = Float32(x)
-    return Float16(evalpoly(x*x, (1.0f0, -4.934802f0, 4.058712f0, -1.3352624f0, 0.23531426f0, -0.0255071f0)))
+    return evalpoly(x*x, (1.0f0, -4.934802f0, 4.058712f0, -1.3352624f0, 0.23531426f0, -0.0255071f0))
 end
 
 """
@@ -867,12 +885,56 @@ function sincospi(_x::T) where T<:Union{IEEEFloat, Rational}
     return si, co
 end
 
+"""
+    tanpi(x)
+
+Compute ``\\tan(\\pi x)`` more accurately than `tan(pi*x)`, especially for large `x`.
+
+See also [`tand`](@ref), [`sinpi`](@ref), [`cospi`](@ref), [`sincospi`](@ref).
+"""
+
+function tanpi(_x::T) where T<:Union{IEEEFloat, Rational}
+    # This is modified from sincospi.
+    # Would it be faster or more accurate to make a tanpi_kernel?
+    x = abs(_x)
+    if !isfinite(x)
+        isnan(x) && return x
+        throw(DomainError(x, "`x` cannot be infinite."))
+    end
+    # For large x, answers are all zero.
+    # All integer values for floats larger than maxintfloat are even.
+    if T <: AbstractFloat
+        x >= maxintfloat(T) && return copysign(zero(T), _x)
+    end
+
+    # reduce to interval [0, 0.5]
+    n = round(2*x)
+    rx = float(muladd(T(-.5), n, x))
+    n = Int64(n) & 3
+    si, co = sinpi_kernel_wide(rx), cospi_kernel_wide(rx)
+    if n==0
+        si, co = si, co
+    elseif n==1
+        si, co  = co, zero(T)-si
+    elseif n==2
+        si, co  = zero(T)-si, zero(T)-co
+    else
+        si, co  = zero(T)-co, si
+    end
+    si = ifelse(signbit(_x), -si, si)
+    return float(T)(si / co)
+end
+
 sinpi(x::Integer) = x >= 0 ? zero(float(x)) : -zero(float(x))
 cospi(x::Integer) = isodd(x) ? -one(float(x)) : one(float(x))
+tanpi(x::Integer) = x >= 0 ? (isodd(x) ? -zero(float(x)) : zero(float(x))) :
+                             (isodd(x) ? zero(float(x)) : -zero(float(x)))
 sincospi(x::Integer) = (sinpi(x), cospi(x))
 sinpi(x::Real) = sin(pi*x)
 cospi(x::Real) = cos(pi*x)
 sincospi(x::Real) = sincos(pi*x)
+tanpi(x::Real) = tan(pi*x)
+tanpi(x::Complex) = sinpi(x) / cospi(x) # Is there a better way to do this?
 
 function sinpi(z::Complex{T}) where T
     F = float(T)

test/math.jl

@@ -69,8 +69,9 @@ end
     @test repr(Any[pi ℯ; ℯ pi]) == "Any[π ℯ; ℯ π]"
     @test string(pi) == "π"
 
-    @test sin(π) === sinpi(1) == tan(π) == sinpi(1 // 1) == 0
-    @test cos(π) === cospi(1) == sec(π) == cospi(1 // 1) == -1
+    @test sin(π) == sind(180) === sinpi(1) === sinpi(1//1) == tan(π) == 0
+    @test tan(π) == tand(180) === tanpi(1) === tanpi(1//1) === -0.0
+    @test cos(π) == cosd(180) === cospi(1) === cospi(1//1) == sec(π) == -1
     @test csc(π) == 1/0 && cot(π) == -1/0
     @test sincos(π) === sincospi(1) == (0, -1)
 end
@@ -181,6 +182,7 @@ end
             @test cbrt(x) ≈ cbrt(big(x))
             @test cos(x) ≈ cos(big(x))
             @test cosh(x) ≈ cosh(big(x))
+	    @test cospi(x) ≈ cospi(big(x))
             @test exp(x) ≈ exp(big(x))
             @test exp10(x) ≈ exp10(big(x))
             @test exp2(x) ≈ exp2(big(x))
@@ -194,9 +196,11 @@ end
             @test log2(x) ≈ log2(big(x))
             @test sin(x) ≈ sin(big(x))
             @test sinh(x) ≈ sinh(big(x))
+            @test sinpi(x) ≈ sinpi(big(x))
             @test sqrt(x) ≈ sqrt(big(x))
             @test tan(x) ≈ tan(big(x))
             @test tanh(x) ≈ tanh(big(x))
+	    @test tanpi(x) ≈ tanpi(big(x))
             @test sec(x) ≈ sec(big(x))
             @test csc(x) ≈ csc(big(x))
             @test secd(x) ≈ secd(big(x))
@@ -499,6 +503,22 @@ end
             @test cospi(convert(T,-1.5))::fT ⩲ zero(fT)
             @test_throws DomainError cospi(convert(T,Inf))
         end
+	@testset "trig pi functions accuracy" for numerator in -20:1:20
+	    for func in (sinpi, cospi, tanpi,
+	                 x -> sincospi(x)[1],
+			 x -> sincospi(x)[2])
+                x = numerator // 20
+		# Check that rational function works
+		@test func(x) ≈ func(BigFloat(x))
+		# Use short value so that wider values will be exactly equal
+		shortx = Float16(x)
+		# Compare to BigFloat value
+		bigvalue = func(BigFloat(shortx))
+		for T in (Float16,Float32,Float64)
+		    @test func(T(shortx)) ≈ T(bigvalue)
+                end
+            end
+	end
         @testset begin
             # If the machine supports fma (fused multiply add), we require exact equality.
             # Otherwise, we only require approximate equality.
@@ -529,14 +549,18 @@ end
     @test ismissing(scdm[2])
 end
 
-@testset "Integer and Inf args for sinpi/cospi/sinc/cosc" begin
+@testset "Integer and Inf args for sinpi/cospi/tanpi/sinc/cosc" begin
     for (sinpi, cospi) in ((sinpi, cospi), (x->sincospi(x)[1], x->sincospi(x)[2]))
-        @test sinpi(1) == 0
-        @test sinpi(-1) == -0
+        @test sinpi(1) === 0.0
+        @test sinpi(-1) === -0.0
         @test cospi(1) == -1
         @test cospi(2) == 1
     end
 
+    @test tanpi(1) === -0.0
+    @test tanpi(-1) === 0.0
+    @test tanpi(2) === 0.0
+    @test tanpi(-2) === -0.0
     @test sinc(1) == 0
     @test sinc(complex(1,0)) == 0
     @test sinc(0) == 1
@@ -589,14 +613,15 @@ end
     end
 end
 
-@testset "Irrational args to sinpi/cospi/sinc/cosc" begin
+@testset "Irrational args to sinpi/cospi/tanpi/sinc/cosc" begin
     for x in (pi, ℯ, Base.MathConstants.golden)
         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)) ≈ ComplexF64(sinpi(complex(big(x), big(x))))
             @test cospi(complex(x, x)) ≈ ComplexF64(cospi(complex(big(x), big(x))))
         end
+	@test tanpi(x) ≈ Float64(tanpi(big(x)))
         @test sinc(x)  ≈ Float64(sinc(big(x)))
         @test cosc(x)  ≈ Float64(cosc(big(x)))
         @test sinc(complex(x, x))  ≈ ComplexF64(sinc(complex(big(x),  big(x))))
@@ -626,7 +651,7 @@ end
 end
 
 @testset "trig function type stability" begin
-    @testset "$T $f" for T = (Float32,Float64,BigFloat,Rational{Int16},Complex{Int32},ComplexF16), f = (sind,cosd,sinpi,cospi)
+    @testset "$T $f" for T = (Float32,Float64,BigFloat,Rational{Int16},Complex{Int32},ComplexF16), f = (sind,cosd,sinpi,cospi,tanpi)
         @test Base.return_types(f,Tuple{T}) == [float(T)]
     end
     @testset "$T sincospi" for T = (Float32,Float64,BigFloat,Rational{Int16},Complex{Int32},ComplexF16)