Add unnormalized sinc/sinhc functions

docs/src/functions_list.md

@@ -151,3 +151,9 @@ ellipe
 eta
 zeta
 ```
+
+## Special Trigonometric Functions
+```@docs
+sincu
+sinhcu
+```

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` |
+

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

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

test/runtests.jl

@@ -33,7 +33,8 @@ tests = [
     "logabsgamma",
     "sincosint",
     "other_tests",
-    "chainrules"
+    "chainrules",
+    "trig"
 ]
 
 const testdir = dirname(@__FILE__)

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