diff --git a/docs/src/functions_list.md b/docs/src/functions_list.md index dc595591..4f6f1203 100644 --- a/docs/src/functions_list.md +++ b/docs/src/functions_list.md @@ -151,3 +151,9 @@ ellipe eta zeta ``` + +## Special Trigonometric Functions +```@docs +sincu +sinhcu +``` \ No newline at end of file diff --git a/docs/src/functions_overview.md b/docs/src/functions_overview.md index ca22f5ac..5ac3d5b4 100644 --- a/docs/src/functions_overview.md +++ b/docs/src/functions_overview.md @@ -118,3 +118,11 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi |:-------- |:----------- | | [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` | | [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` | + +## Special Trigonometric Functions + +| Function | Description | +|:-------- |:----------- | +| [`sincu(x)`](@ref SpecialFunctions.sincu) | [Unnormalized sine cardinal function](https://en.wikipedia.org/wiki/Sinc_function) at `x` | +| [`sinhcu(x)`](@ref SpecialFunctions.sinhcu) | [Unnormalized hyperbolic sinc function](https://mathworld.wolfram.com/SinhcFunction.html) at `x` | + diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 0d296cc0..492a410f 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -90,7 +90,11 @@ export expintx, sinint, cosint, - lbinomial + lbinomial, + + # Cardinal trigonometric + sincu, + sinhcu include("bessel.jl") include("erf.jl") @@ -103,6 +107,8 @@ include("gamma.jl") include("gamma_inc.jl") include("betanc.jl") include("beta_inc.jl") +include("trig.jl") + if !isdefined(Base, :get_extension) include("../ext/SpecialFunctionsChainRulesCoreExt.jl") end diff --git a/src/trig.jl b/src/trig.jl new file mode 100644 index 00000000..bdadfa07 --- /dev/null +++ b/src/trig.jl @@ -0,0 +1,71 @@ +# u corresponds to unnormalized + +# sinc/sincu is zero when the real part is Inf and imag is finite +isinf_real(x::Real) = isinf(x) +isinf_real(x::Number) = isinf(real(x)) && isfinite(imag(x)) + +# sinhc/sinhcu is zero when the imag part is Inf and real is finite +isinf_imag(x::Real) = false +isinf_imag(x::Number) = isfinite(real(x)) && isinf(imag(x)) + +# sincu copied from Boost library and correct limit behavior added +# https://www.boost.org/doc/libs/1_87_1/boost/math/special_functions/sinc.hpp +""" + sincu(x) + +Compute the unnormalized sinc function ``\\operatorname{sincu}(x) = \\sin(x) / (x)`` +with accuracy near the origin. +""" +sincu(x) = _sinc(float(x)) +function _sinc(x::Union{T,Complex{T}}) where {T} + if isinf_real(x) + return zero(x) + end + + nrm = fastabs(x) + if nrm >= 3.3*sqrt(sqrt(eps(T))) + return sin(x)/x + else + # |x| < (eps*120)^(1/4) + return 1 - x*x/6 + end +end + +# sinhcu copied from Boost library and correct limit behavior added +# https://www.boost.org/doc/libs/1_87_1/boost/math/special_functions/sinhc.hpp + +""" + sinhcu(x) + +Compute the unnormalized sinhc function ``\\operatorname{sinhcu}(x) = \\sinh(x) / (x)`` +with accuracy accuracy near the origin. +""" +sinhcu(x) = _sinhcu(float(x)) +function _sinhcu(x::Union{T,Complex{T}}) where {T} + taylor_0_bound = eps(T) + taylor_2_bound = sqrt(taylor_0_bound) + taylor_n_bound = sqrt(taylor_2_bound) + + if isinf_imag(x) + return zero(x) + end + + nrm = fastabs(x) + + if nrm >= taylor_n_bound || isnan(nrm) + return sinh(x)/x + else + # approximation by taylor series in x at 0 up to order 0 + res = one(x) + if nrm >= taylor_0_bound + x2 = x*x + # approximation by taylor series in x at 0 up to order 2 + res += x2/6 + if nrm >= taylor_2_bound + # approximation by taylor series in x at 0 up to order 4 + res += (x2*x2)/120 + end + end + return res + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 69dcafe7..6d92a2e1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,7 +33,8 @@ tests = [ "logabsgamma", "sincosint", "other_tests", - "chainrules" + "chainrules", + "trig" ] const testdir = dirname(@__FILE__) diff --git a/test/trig.jl b/test/trig.jl new file mode 100644 index 00000000..73690671 --- /dev/null +++ b/test/trig.jl @@ -0,0 +1,57 @@ +@testset "Special trigonometric functions" begin + @testset "sincu (unnormalized sinc)" begin + a = 1.0 + @test sincu(a) ≈ sin(a)/a + @test sincu(-a) ≈ sin(-a)/-a + + b = complex(2.0, 3.0) + @test sincu(b) ≈ sin(b)/b + @test sincu(-b) ≈ sin(-b)/-b + + # Right below threshold + c = sqrt(sqrt(eps(Float64))) + @test sincu(c) == 1 - c^2/6 + d = complex(0,c) + @test sincu(d) == 1 - d^2/6 + + # Limits/nans + @test sincu(0) == 1.0 + @test sincu(Inf) == 0.0 + @test sincu(-Inf) == 0.0 + @test isnan(sincu(NaN)) + end + + @testset "sinhcu (unnormalized sinhc)" begin + a = 1.0 + @test sinhcu(a) ≈ sinh(a)/a + @test sinhcu(-a) ≈ sinh(-a)/-a + + b = complex(2.0, 3.0) + @test sinhcu(b) ≈ sinh(b)/b + @test sinhcu(-b) ≈ sinh(-b)/-b + + # 0th order approximation + c = sqrt(eps(Float64))-eps(Float64) + @test sinhcu(c) == 1.0 + d = complex(0,c) + @test sinhcu(d) == 1.0 + + # 2nd order approximation + e = sqrt(eps(Float64)) + @test sinhcu(e) == 1 + e^2/6 + f = complex(0,e) + @test sinhcu(f) == 1 + f^2/6 + + # 4th order approximation + g = sqrt(sqrt(eps(Float64))) + @test sinhcu(g) == 1 + g^2/6 + g^4/120 + h = complex(0,g) + @test sinhcu(h) == 1 + h^2/6 + h^4/120 + + # Limits/nans + @test sinhcu(0.0) == 1.0 + @test sinhcu(Inf*im) == 0.0 + @test sinhcu(-Inf*im) == 0.0 + @test isnan(sinhcu(NaN)) + end +end