From 686bf94bac386df7045493c1ec94790827cf28e1 Mon Sep 17 00:00:00 2001 From: mattsignorelli Date: Sun, 6 Apr 2025 11:30:24 -0400 Subject: [PATCH 1/3] Add (hyperbolic) unnormalized sinus cardinal functions --- src/SpecialFunctions.jl | 8 +++++- src/trigcardinal.jl | 52 +++++++++++++++++++++++++++++++++++++ test/runtests.jl | 3 ++- test/trigcardinal.jl | 57 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 118 insertions(+), 2 deletions(-) create mode 100644 src/trigcardinal.jl create mode 100644 test/trigcardinal.jl diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 0d296cc0..4cb86d38 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("trigcardinal.jl") + if !isdefined(Base, :get_extension) include("../ext/SpecialFunctionsChainRulesCoreExt.jl") end diff --git a/src/trigcardinal.jl b/src/trigcardinal.jl new file mode 100644 index 00000000..17237278 --- /dev/null +++ b/src/trigcardinal.jl @@ -0,0 +1,52 @@ +# u corresponds to unnormalized + +# sincu copied exactly from Boost library +# https://www.boost.org/doc/libs/1_87_1/boost/math/special_functions/sinc.hpp +sincu(x) = _sinc(float(x)) +function _sinc(x::Union{T,Complex{T}}) where {T} + if Base.Math.isinf_real(x) + return zero(x) + else + nrm = Base.Math.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 +end + +# sinhcu copied exactly from Boost library +# https://www.boost.org/doc/libs/1_87_1/boost/math/special_functions/sinhc.hpp +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 Base.Math.isinf_real(x) + return typemax(x) + elseif isnan(x) + return NaN + end + + nrm = Base.Math.fastabs(x) + + if nrm >= taylor_n_bound + 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..53c4e694 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,7 +33,8 @@ tests = [ "logabsgamma", "sincosint", "other_tests", - "chainrules" + "chainrules", + "trigcardinal" ] const testdir = dirname(@__FILE__) diff --git a/test/trigcardinal.jl b/test/trigcardinal.jl new file mode 100644 index 00000000..d38df4a3 --- /dev/null +++ b/test/trigcardinal.jl @@ -0,0 +1,57 @@ +@testset "Cardinal trigonometric" 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) == Inf + @test sinhcu(-Inf) == Inf + @test isnan(sinhcu(NaN)) + end +end From fab1241a9e08f08d5006f2d46ae05622c58a30a7 Mon Sep 17 00:00:00 2001 From: mattsignorelli Date: Sun, 6 Apr 2025 12:11:42 -0400 Subject: [PATCH 2/3] fix sinhcu limit, work for Julia min --- src/trigcardinal.jl | 34 ++++++++++++++++++++-------------- test/trigcardinal.jl | 4 ++-- 2 files changed, 22 insertions(+), 16 deletions(-) diff --git a/src/trigcardinal.jl b/src/trigcardinal.jl index 17237278..0ee44d3e 100644 --- a/src/trigcardinal.jl +++ b/src/trigcardinal.jl @@ -1,19 +1,27 @@ # 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 exactly from Boost library # https://www.boost.org/doc/libs/1_87_1/boost/math/special_functions/sinc.hpp sincu(x) = _sinc(float(x)) function _sinc(x::Union{T,Complex{T}}) where {T} - if Base.Math.isinf_real(x) + if isinf_real(x) return zero(x) + end + + nrm = fastabs(x) + if nrm >= 3.3*sqrt(sqrt(eps(T))) + return sin(x)/x else - nrm = Base.Math.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 + # |x| < (eps*120)^(1/4) + return 1 - x*x/6 end end @@ -25,15 +33,13 @@ function _sinhcu(x::Union{T,Complex{T}}) where {T} taylor_2_bound = sqrt(taylor_0_bound) taylor_n_bound = sqrt(taylor_2_bound) - if Base.Math.isinf_real(x) - return typemax(x) - elseif isnan(x) - return NaN + if isinf_imag(x) + return zero(x) end - nrm = Base.Math.fastabs(x) + nrm = fastabs(x) - if nrm >= taylor_n_bound + if nrm >= taylor_n_bound || isnan(nrm) return sinh(x)/x else # approximation by taylor series in x at 0 up to order 0 diff --git a/test/trigcardinal.jl b/test/trigcardinal.jl index d38df4a3..23a4bff4 100644 --- a/test/trigcardinal.jl +++ b/test/trigcardinal.jl @@ -50,8 +50,8 @@ # Limits/nans @test sinhcu(0.0) == 1.0 - @test sinhcu(Inf) == Inf - @test sinhcu(-Inf) == Inf + @test sinhcu(Inf*im) == 0.0 + @test sinhcu(-Inf*im) == 0.0 @test isnan(sinhcu(NaN)) end end From 4377c7a332382aaae069f65be6eacdaf8eedb6c3 Mon Sep 17 00:00:00 2001 From: mattsignorelli Date: Sun, 6 Apr 2025 12:41:44 -0400 Subject: [PATCH 3/3] update docs + remove abuses of word cardinal --- docs/src/functions_list.md | 6 ++++++ docs/src/functions_overview.md | 8 ++++++++ src/SpecialFunctions.jl | 2 +- src/{trigcardinal.jl => trig.jl} | 19 ++++++++++++++++--- test/runtests.jl | 2 +- test/{trigcardinal.jl => trig.jl} | 2 +- 6 files changed, 33 insertions(+), 6 deletions(-) rename src/{trigcardinal.jl => trig.jl} (78%) rename test/{trigcardinal.jl => trig.jl} (96%) 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 4cb86d38..492a410f 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -107,7 +107,7 @@ include("gamma.jl") include("gamma_inc.jl") include("betanc.jl") include("beta_inc.jl") -include("trigcardinal.jl") +include("trig.jl") if !isdefined(Base, :get_extension) include("../ext/SpecialFunctionsChainRulesCoreExt.jl") diff --git a/src/trigcardinal.jl b/src/trig.jl similarity index 78% rename from src/trigcardinal.jl rename to src/trig.jl index 0ee44d3e..bdadfa07 100644 --- a/src/trigcardinal.jl +++ b/src/trig.jl @@ -8,14 +8,20 @@ isinf_real(x::Number) = isinf(real(x)) && isfinite(imag(x)) isinf_imag(x::Real) = false isinf_imag(x::Number) = isfinite(real(x)) && isinf(imag(x)) -# sincu copied exactly from Boost library +# 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 @@ -25,8 +31,15 @@ function _sinc(x::Union{T,Complex{T}}) where {T} end end -# sinhcu copied exactly from Boost library +# 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) diff --git a/test/runtests.jl b/test/runtests.jl index 53c4e694..6d92a2e1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,7 +34,7 @@ tests = [ "sincosint", "other_tests", "chainrules", - "trigcardinal" + "trig" ] const testdir = dirname(@__FILE__) diff --git a/test/trigcardinal.jl b/test/trig.jl similarity index 96% rename from test/trigcardinal.jl rename to test/trig.jl index 23a4bff4..73690671 100644 --- a/test/trigcardinal.jl +++ b/test/trig.jl @@ -1,4 +1,4 @@ -@testset "Cardinal trigonometric" begin +@testset "Special trigonometric functions" begin @testset "sincu (unnormalized sinc)" begin a = 1.0 @test sincu(a) ≈ sin(a)/a