Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement tanpi #48575

Merged
merged 1 commit into from
Feb 7, 2023
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Implement tanpi
See issue: #48226
bobcassels committed Feb 7, 2023
commit 7ab07f04b8b8c878543abd18cc45f8ac8f03c54a
1 change: 1 addition & 0 deletions base/exports.jl
Original file line number Diff line number Diff line change
@@ -352,6 +352,7 @@ export
tan,
tand,
tanh,
tanpi,
trailing_ones,
trailing_zeros,
trunc,
2 changes: 1 addition & 1 deletion base/math.jl
Original file line number Diff line number Diff line change
@@ -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,
4 changes: 2 additions & 2 deletions base/mpfr.jl
Original file line number Diff line number Diff line change
@@ -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
78 changes: 70 additions & 8 deletions base/special/trig.jl
Original file line number Diff line number Diff line change
@@ -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⁴ =*
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
r =*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 *
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)
39 changes: 32 additions & 7 deletions test/math.jl
Original file line number Diff line number Diff line change
@@ -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)