From 7ab07f04b8b8c878543abd18cc45f8ac8f03c54a Mon Sep 17 00:00:00 2001 From: Bob Cassels Date: Tue, 7 Feb 2023 09:22:53 -0500 Subject: [PATCH] Implement tanpi See issue: https://github.com/JuliaLang/julia/issues/48226 --- base/exports.jl | 1 + base/math.jl | 2 +- base/mpfr.jl | 4 +-- base/special/trig.jl | 78 +++++++++++++++++++++++++++++++++++++++----- test/math.jl | 39 ++++++++++++++++++---- 5 files changed, 106 insertions(+), 18 deletions(-) diff --git a/base/exports.jl b/base/exports.jl index 600b36b6c37c6..26b24b85651a2 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -352,6 +352,7 @@ export tan, tand, tanh, + tanpi, trailing_ones, trailing_zeros, trunc, diff --git a/base/math.jl b/base/math.jl index f41057c76cfc2..0faa6a00c3f53 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, sincospi, sinc, cosc, + sinpi, cospi, sincospi, tanpi, sinc, cosc, cosd, cotd, cscd, secd, sind, tand, sincosd, acosd, acotd, acscd, asecd, asind, atand, rad2deg, deg2rad, diff --git a/base/mpfr.jl b/base/mpfr.jl index d42beb0c59190..4485265b580de 100644 --- a/base/mpfr.jl +++ b/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 diff --git a/base/special/trig.jl b/base/special/trig.jl index 929e259913104..ed92f83bb52e2 100644 --- a/base/special/trig.jl +++ b/base/special/trig.jl @@ -725,29 +725,41 @@ 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²)) @@ -755,13 +767,19 @@ end 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) diff --git a/test/math.jl b/test/math.jl index f9af521de61ca..f08cc6c32ce15 100644 --- a/test/math.jl +++ b/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,7 +613,7 @@ 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))) @@ -597,6 +621,7 @@ end @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)