Adding jinc function (#266)

* Add jinc (also called sombrero, besinc) function

* Add tests for jinc function

* Add more jinc tests to verify that polynomial works as well

* Add redefinition of fastabs for jinc and clean jinc body

* Some code comments regarded cleaning of jinc function

* Add fastabs to SpecialFunctions

* Add relative error comparison to tests

* Add proper relative error testing

* Adapt jinc function to get type stability

* Add test cases for jinc function

* Use fastabs in expint and gamma for complex number checking

* Adapt \cong for complex numbers

* Add tests for fastabs

* Replace relerrc with complex relerr in tests for gamma and expint

* Fix test cases for jinc with correct big expression

* Revert \cong threshold back to 1e-6 for Float32

* Relax accuracy threshold for Float32 slightly

* Use \isapprox for all complex number tests of jinc

* Remove complex argument handling of \cong in runtest

* Remove explicit zero handling in jinc because evalpoly already does

* Add complex handling for relerr in runtests

* Add jinc function to documentation

docs/src/functions_list.md

@@ -50,6 +50,7 @@ SpecialFunctions.besseli
 SpecialFunctions.besselix
 SpecialFunctions.besselk
 SpecialFunctions.besselkx
+SpecialFunctions.jinc
 SpecialFunctions.ellipk
 SpecialFunctions.ellipe
 SpecialFunctions.eta

docs/src/functions_overview.md

@@ -78,6 +78,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 | [`besselix(nu,z)`](@ref SpecialFunctions.besselix)            | scaled modified Bessel function of the first kind of order `nu` at `z`                                                                                          |
 | [`besselk(nu,z)`](@ref SpecialFunctions.besselk)              | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z`                                               |
 | [`besselkx(nu,z)`](@ref SpecialFunctions.besselkx)            | scaled modified Bessel function of the second kind of order `nu` at `z`                                                                                         |
+| [`jinc(x)`](@ref SpecialFunctions.jinc)                       | scaled [Bessel function of the first kind divided by `x`](https://en.wikipedia.org/wiki/Sombrero_function). A.k.a. sombrero or besinc                           |
 
 
 ## [Elliptic Integrals](https://dlmf.nist.gov/19)

src/SpecialFunctions.jl

@@ -27,6 +27,7 @@ export
     bessely0,
     bessely1,
     besselyx,
+    jinc,
     dawson,
     ellipk,
     ellipe,
@@ -83,4 +84,13 @@ for f in (:beta, :lbeta)
 end
 polygamma(m::Integer, x::Missing) = missing
 
+ # In future just use `fastabs` from Base.Math
+ # https://github.com/JuliaLang/julia/blob/93fb785831dcfcc442f82fab8746f0244c5274ae/base/special/trig.jl#L1057
+if isdefined(Base.Math, :fastabs)
+    import Base.Math: fastabs
+else
+    fastabs(x::Number) = abs(x)
+    fastabs(x::Complex) = abs(real(x)) + abs(imag(x))
+end
+
 end # module

src/bessel.jl

@@ -758,3 +758,33 @@ External links: [DLMF](https://dlmf.nist.gov/10.2.6), [Wikipedia](https://en.wik
 See also: [`hankelh2(nu,x)`](@ref SpecialFunctions.hankelh2)
 """
 hankelh2x(nu, z) = besselhx(nu, 2, z)
+
+"""
+    jinc(x)
+
+Bessel function of the first kind divided by `x`.
+Following convention: ``\\operatorname{jinc}{x} = \\frac{2 \\cdot J_1{\\pi x}}{\\pi x}``.
+Sometimes known as sombrero or besinc function.
+
+External links: [Wikipedia](https://en.wikipedia.org/wiki/Sombrero_function)
+"""
+jinc(x::Number) = _jinc(float(x))
+_jinc(x::Number) = iszero(x) ? one(x) : _jinc_core(x)
+_jinc(x::Float16) = Float16(_jinc(Float32(x)))
+_jinc(x::ComplexF16) = ComplexF16(_jinc(ComplexF32(x)))
+
+# below these thresholds we evaluate a fourth order Taylor polynomial
+_jinc_threshold(::Type{Float64}) = 0.002
+_jinc_threshold(::Type{Float32}) = 0.05f0
+
+
+ # for small arguments, a Taylor series is faster
+@inline function _jinc(x::Union{T,Complex{T}}) where {T<:Union{Float32,Float64}}
+    if fastabs(x) < _jinc_threshold(T)
+        return @evalpoly(x^2, T(1), -T(π)^2/8, T(π)^4/192)
+    else
+        return _jinc_core(x)
+    end
+end
+ # the actual definition of jinc
+_jinc_core(x::Number) =  2 * besselj1(π*x) / (π*x)

src/expint.jl

@@ -178,7 +178,7 @@ function En_cf_nogamma(ν::Number, z::Number, n::Int=1000)
         conv && i > 4 && break
 
         # rescale
-        if max(abs(real(A)), abs(imag(A))) > scale
+        if fastabs(A) > scale
             A     /= scale
             Aprev /= scale
             B     /= scale

src/gamma.jl

@@ -474,7 +474,7 @@ Dirichlet eta function ``\\eta(s) = \\sum^\\infty_{n=1}(-1)^{n-1}/n^{s}``.
 """
 function eta(z::ComplexOrReal{Float64})
     δz = 1 - z
-    if abs(real(δz)) + abs(imag(δz)) < 7e-3 # Taylor expand around z==1
+    if fastabs(δz) < 7e-3 # Taylor expand around z==1
         return 0.6931471805599453094172321214581765 *
                @evalpoly(δz,
                          1.0,

test/bessel.jl

@@ -271,4 +271,99 @@ end
             @test f(0,1) ≈ f(0,Complex{Float16}(1))
         end
     end
+
+    @testset "jinc function" begin
+
+        x1 = big"1.8"
+        res1 = big"-0.116401402395234261196060221416100814277048989838340270527643894539781608507904800871249915340068849141209639297793338757544621191431028495364540423087719279432811811"
+        x2 = big"0.00001"
+        res2 = big"0.9999999998766299449914564074224855053704639002961349883799411022244488435839201448375732357069707105107902868318079100001712780899305533775708875568993566765287247356"
+        x3 = big"0.02"
+        res3 = big"0.9995066009475120760121935544369072906921634576001177613216005027724298234311657036029970452969291362640681823594159438475599412434802474217121139032145806264666152465837468052082455102325311416327533"
+
+        root1 = big"1.21966989126650445492653884746525517787935933077511212945638126557694328028076014425087191879391333148570410142424433731699"
+        for T in [Float16, Float32, Float64]
+            @test jinc(T(x1)) ≅ T(res1)
+            @test jinc(T(0)) ≅ T(1)
+            x = root1 
+            # check the first root for negative and positive
+            @test jinc(T(x)) + T(1) ≅ T(1)
+            @test jinc(T(-x)) + T(1) ≅ T(1)
+
+            # check for symmetry
+            for x = -5.0:1.0:3.0
+                @test jinc(T(10^x)) ≅ jinc(T(- 10^x))
+
+            end
+
+            x = T(x2)
+            res = T(res2)
+            @test jinc(x) ≅ res
+
+            x = T(x3)
+            res = T(res3)
+            @test jinc(x) ≅ res
+
+            # type stability for 0
+            @test typeof(jinc(T(0))) === T
+        end
+
+        # test big float
+        begin
+            T = BigFloat
+            @test T(jinc(T(x1))) ≈  T(res1)
+            #= @test jinc(T(0)) ≈ T(1) =#
+            x = root1
+            # check the first root for negative and positive
+            @test jinc(T(x)) + T(1) ≈ T(1)
+            @test jinc(T(-x)) + T(1) ≈ T(1)
+
+            # check for symmetry
+            for x = -5.0:1.0:3.0
+                @test jinc(T(10^x)) ≈ jinc(T(- 10^x))
+
+            end
+
+            x = T(x2)
+            res = T(res2)
+            @test T(jinc(x)) ≈ res
+
+            x = T(x3)
+            res = T(res3)
+            @test T(jinc(x)) ≈ res
+
+            # type stability for 0
+            @test typeof(jinc(T(0))) === T
+        end
+
+
+        x1 = big"1.0" + big"0.0001"*1im
+        res1 = big"0.1811917493658519438199345187518953631147425189362440669941643986059307908067224704314102961998438798566229144099942611910610729246007355473812330821120666591291397178705295331963573655500130149361351" - big"0.0000970867870757251899047609572999879575662488432474075980361828257989869563037501527523384396635449088435640500216646530541742421182815911923875286404692082399771721882711636718746219759416544606077"*1im
+
+        x2 = big"0.00001" + big"10"*1im
+        res2 = big"1.971105539196620210542904558564976891448851717386016064466205835417372529649732623778146798009632930245313240967710133352055225553347549825667896762467087506209685330050397030027682762563950964962993e11" - big"5.89917589196246805438832769814196828580119495400655210987269285701113079909408122244782544425726443894463237062342568390611353700647229983489787030738912900107228049806443988000825445873073676204e6"*1im
+
+        x3 = big"0.00002"*1im
+        res3 = big"1.0000000004934802201356421734767362282021102443681623763651414603818274528461323973841772525693725599437763705213335500965076723534424132751557185487068200219257886657409530057472519875335446749379510"
+
+        # check also some complex numbers
+        # besselj1 is not type stable and returns Complex{Float64}
+        # the stricter ≅ doesn't hold here
+        # besselj1 seems to produce inaccurate results for some Float32 cases
+        for T in [Complex{Float16}, Complex{Float32}, Complex{Float64}]
+            x = T(x1)
+            res = T(res1)
+            @test T(jinc(x)) ≈ res
+
+            x = T(x2)
+            res = T(res2)
+            @test T(jinc(x)) ≈ res
+
+            x = T(x3)
+            res = T(res3)
+            @test T(jinc(x)) + T(1im) ≅ res + T(1im) 
+
+        end
+
+    end
 end

test/expint.jl

@@ -32,7 +32,7 @@ using Base.MathConstants
 , 18.922784338495287 - 0.029203673105875738im, 18.422784340698875 - 0.029203673041504553im, 17.92278434433198 - 0.029203672935374417im, 17.422784350321955 - 0.029203672760395406im, 16.922784360197753 - 0.029203672471903787im, 16.42278437648019 - 0.029203671996261525im, 15.922784403325396 - 0.02920367121206002im, 15.422784447585654 - 0.029203669919130378im, 14.922784520558482 - 0.029203667787449886im, 14.42278464087033 - 0.029203664272903267im, 13.922784839231014 - 0.029203658478396427im, 13.422785166272453 - 0.02920364892487225im, 12.922785705472512 - 0.029203633173780545im, 12.422786594462801 - 0.029203607204639166im, 11.922788060159139 - 0.029203564388813802im, 11.422790476681525 - 0.029203493797588825im, 10.922794460847014 - 0.029203377412707156im, 10.422801029608085 - 0.029203185527489635im, 9.92281185961712 - 0.029202869165002185im, 9.422829715155434 - 0.029202347578921022im, 8.922859153613475 - 0.02920148764919001im, 8.422907688480091 - 0.02920006992002703im, 7.922987706377455 - 0.02919773263004735im, 7.4231196266023955 - 0.029193879498721455im, 6.923337107299312 - 0.02918752786907706im, 6.423695620754588 - 0.029177058818573274im, 5.924286569298801 - 0.029159806468544838im, 5.425260497751445 - 0.02913138441743693im, 4.926865199419432 - 0.029084584836757257im, 4.429508095337361 - 0.029007589432432816im, 3.9338578706686858 - 0.02888109015248326im, 3.4410087804823353 - 0.028673731153572076im, 2.9527427760716067 - 0.028335098763862032im, 2.4719383098140226 - 0.02778549583478644im, 2.0031828640491396 - 0.026902526581439178im, 1.5536245933117467 - 0.025507584670758247im, 1.1339838503079331 - 0.023363827088483065im, 0.7593072088540906 - 0.020215469144043718im, 0.4483165720604514 - 0.015923527983522815im, 0.21922704942196994 - 0.010743431043978405im, 0.08026140877749337 - 0.0056148853587022im, 0.01865598084904698 - 0.0019258134104280294im, 0.0020975029211462598 - 0.00032968586251067337im, 7.260678405534787e-5 - 1.7926254192086143e-5im, 3.6394029588254425e-7 - 1.4661826414313507e-7im, 7.41111524749911e-11 - 5.228655363875185e-11im, 6.721230570877014e-17 - 1.0361363222368277e-16im, -1.864169823362996e-27 - 3.5720968579640937e-26im, -8.139775526992696e-42 - 4.2816798952068195e-42im, -8.5705278203875e-68 + 2.3494757832291492e-67im, 1.5313896916407397e-109 - 1.893701952982838e-109im, 1.3253076410048056e-178 + 1.2528740964916723e-178im, 2.2308653435430693e-292 - 1.5308786810615063e-292im, 0.0 - 0.0im, -0.0 - 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, -0.0 + 0.0im, 0.0 - 0.0im, -0.0 - 0.0im, 19.42278433697267 - 0.42920367234736084im, 18.922784338188503 - 0.42920367179092506im, 18.422784340193076 - 0.42920367087351746im, 17.922784343498055 - 0.42920366936096793im, 17.422784348947047 - 0.42920366686719563im, 16.92278435793092 - 0.42920366275565996im, 16.422784372742814 - 0.42920365597688387im, 15.922784397163502 - 0.4292036448005718im, 15.42278443742641 - 0.4292036263739489im, 14.922784503808721 - 0.4292035959935854im, 14.422784613254645 - 0.42920354590483817im, 13.922784793700458 - 0.4292034633224673im, 13.422785091205279 - 0.42920332716718834im, 12.92278558170773 - 0.4292031026851722im, 12.422786390409346 - 0.42920273257713776im, 11.922787723732338 - 0.4292021223728026im, 11.42278992200878 - 0.4292011163177121im, 10.922793546349741 - 0.4291994576181139im, 10.422799521866436 - 0.4291967228979313im, 9.922809373797097 - 0.4291922141422782im, 9.422825616800965 - 0.4291847805579199im, 8.922852396759225 - 0.4291725249129162im, 8.422896548827053 - 0.42915231948695526im, 7.92296934159787 - 0.4291190083193409im, 7.4230893520139265 - 0.4290640927834337im, 6.923287203307911 - 0.4289735667620085im, 6.423613371154225 - 0.4288243536951455im, 5.924151039178316 - 0.42857844921551824im, 5.425037254279619 - 0.42817331000498887im, 4.926497697671408 - 0.4275061325334096im, 4.428903718874863 - 0.42640827127518255im, 3.93286557186195 - 0.4246039756763536im, 3.4393839757047893 - 0.42164482116812363im, 2.95009415537207 - 0.41680819592206714im, 2.467652457856141 - 0.4089474035324931im, 1.9963316101738395 - 0.39629004411950314im, 1.5428903310988162 - 0.37622082443650384im, 1.117715990264714 - 0.345200709719881im, 0.7359758618950316 - 0.29924118518393067im, 0.41778887902644557 - 0.23578449257048073im, 0.1849728717554206 - 0.15794164682101264im, 0.050624891660021015 - 0.07977458336698368im, 0.002271877217438789 - 0.02412005211237192im, -0.001975490551566147 - 0.0025306436416951427im, -0.0001392658852773198 + 4.621052039983524e-5im, 8.145281857187584e-7 + 8.59924915803196e-7im, -4.4236745140561634e-10 - 3.40054561329938e-10im, -1.5072272094878889e-16 - 2.45592289751446e-15im, -2.0153450436042867e-24 + 4.5231338334005084e-24im, 3.0603115497994336e-38 - 5.757093007467674e-39im, 1.3186180582171105e-61 + 9.894484316990687e-62im, -1.4612609528592482e-100 - 9.442607791153564e-100im, 2.828756024728357e-163 + 1.1630396096583835e-162im, 2.335195331306485e-266 - 2.241880282600362e-266im, -0.0 + 0.0im, 0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 19.4227843364907 - 0.8292036716852221im, 18.92278433739387 - 0.8292036706992425im, 18.42278433888295 - 0.829203669073637im, 17.922784341338023 - 0.8292036663934665im, 17.422784345385757 - 0.8292036619746125im, 16.92278435205934 - 0.829203654689154im, 16.42278436306222 - 0.8292036426774636im, 15.922784381202906 - 0.8292036228735347im, 15.422784411111834 - 0.8292035902223763im, 14.922784460423323 - 0.8292035363897192im, 14.422784541724223 - 0.8292034476346783im, 13.92278467576675 - 0.8292033013023703im, 13.422784896765515 - 0.8292030600412242im, 12.922785261130894 - 0.8292026622689573im, 12.422785861867869 - 0.8292020064536764im, 11.922786852315776 - 0.8292009251979334im, 11.42278848528851 - 0.8291991425109292im, 10.922791177605953 - 0.8291962033633026im, 10.422795616488505 - 0.8291913575453713im, 9.922802934972712 - 0.8291833681892412im, 9.422815001124505 - 0.8291701960955198im, 8.922834894876095 - 0.8291484793314605im, 8.422867694210316 - 0.8291126753839623im, 7.92292177139537 - 0.8290536472183477im, 7.42301093021208 - 0.8289563331961032im, 6.9231579299103165 - 0.8287959084429407im, 6.4234002959472605 - 0.8285314642336306im, 5.923799902211866 - 0.8280956093926821im, 5.424458774659457 - 0.8273773865745692im, 4.925545161071208 - 0.8261942707329614im, 4.427336549853862 - 0.8242464495407521im, 3.9302906924978545 - 0.8210426567295104im, 3.4351629431935944 - 0.8157811723325826im, 2.9432002972868325 - 0.8071624764467379im, 2.4564625100037736 - 0.7931041315254844im, 1.9783531955862408 - 0.7703339671393372im, 1.5144930910715089 - 0.7338853118696992im, 1.0741247206259925 - 0.6766887196739264im, 0.6722202469513102 - 0.58992874694513im, 0.3319947929095889 - 0.4658866726525998im, 0.08535470889369866 - 0.30643466528170005im, -0.03733900572449072 - 0.1386367619198214im, -0.042067480553635336 - 0.02120741184151642im, -0.005991077985344485 + 0.007220357363965503im, 0.0008386668749499684 + 7.043588351797712e-5im, -1.9573199354809665e-5 + 6.832060917729574e-6im, -6.140424743950199e-8 - 6.232131309598241e-9im, 5.6850215965543594e-12 - 5.394718311174526e-13im, -1.6938945500157091e-18 + 4.0701113104347157e-19im, -1.462310657221111e-29 + 4.0692910148679145e-29im, -1.856502255835392e-46 + 5.895551561127694e-47im, 3.920757198528301e-75 + 5.504396946923525e-75im, -1.1054632599207572e-121 - 1.5359750079774523e-122im, 4.0192755581761105e-199 - 1.0681629965033397e-198im, 0.0 + 0.0im, -0.0 - 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 0.0 - 0.0im, 19.42278433578893 - 1.2292036712630388im, 18.922784336236848 - 1.2292036700031799im, 18.422784336975337 - 1.2292036679260239im, 17.922784338192905 - 1.2292036645013722im, 17.42278434020033 - 1.2292036588550765im, 16.922784343510017 - 1.2292036495459089im, 16.42278434896677 - 1.229203634197686im, 15.92278435796343 - 1.2292036088927445im, 15.42278437279642 - 1.2292035671719501im, 14.922784397251887 - 1.2292034983859899im, 14.422784437572137 - 1.229203384977118im, 13.922784504049007 - 1.2292031979975087im, 13.422784613650869 - 1.2292028897202767im, 12.92278479435388 - 1.229202381457121im, 12.422785092283025 - 1.2292015434730454im, 11.922785583485812 - 1.2292001618714208im, 11.422786393344117 - 1.2291978839969164im, 10.92278772857968 - 1.229194128420797im, 10.422789930024404 - 1.2291879365335145im, 9.922793559629715 - 1.2291777278670097im, 9.422799543936591 - 1.2291608967023038im, 8.922809410660811 - 1.2291331470229703im, 8.422825678873323 - 1.2290873961342346im, 7.922852502617603 - 1.229011967295957im, 7.422896732921453 - 1.2288876105836488im, 6.9229696711123205 - 1.2286825930366383im, 6.4230899659415925 - 1.2283446088960583im, 5.923288407509832 - 1.2277874560252655im, 5.423615878292022 - 1.2268691076253266im, 4.924156590178831 - 1.2253556636275011im, 4.4250502556307625 - 1.2228622026816482im, 3.926529580285248 - 1.2187560440033898im, 3.428984608482909 - 1.211999386761129im, 2.9330755977603307 - 1.2008957200141752im, 2.4399371502299614 - 1.182687871199021im, 1.9515619183211097 - 1.1529402079332955im, 1.4715529768315847 - 1.1046463591926041im, 1.006651810570212 - 1.0271180359240253im, 0.5698682817468509 - 0.9051891431813432im, 0.18654339164699657 - 0.7208595063511073im, -0.09640051243934544 - 0.46385927728183574im, -0.2087824175708163 - 0.16483426811518412im, -0.11405381653019356 + 0.05081491251114011im, 0.024538092462979437 + 0.037639224571643835im, -0.0024341445399757335 - 0.010461008300708324im, 0.0013375970782092976 - 9.523009242235569e-5im, 1.7854684906309634e-5 - 5.564494645760924e-5im, 2.501216238190226e-7 - 3.793280250670107e-7im, -1.5160273071553551e-10 - 1.415734833400537e-10im, -3.115103669354081e-16 + 8.306464570043674e-16im, -1.6389552556716997e-24 - 5.260446033283102e-25im, 7.876560123250489e-39 + 6.611239210483262e-39im, -1.7581312513562253e-62 + 4.713557266756052e-62im, 2.386089044922975e-100 + 9.514717594582031e-101im, -1.5923279422546296e-163 + 2.07656250727342e-163im, -2.3832441726055056e-267 - 4.44002680709751e-267im, 0.0 - 0.0im, 0.0 - 0.0im, -0.0 - 0.0im, -0.0 + 0.0im, -0.0 - 0.0im]
         for (x, y) in zip(inputs, sympy_outputs)
             if y == 5.7333170913660523e-42 + 4.019319913618439e-39im
-                @test relerrc(expint(x), y) ≤ 5e-12
+                @test relerr(expint(x), y) ≤ 5e-12
             else
                 @test expint(x) ≅ y
             end
@@ -46,7 +46,7 @@ using Base.MathConstants
 
     @testset "En" begin
         # ≅ but worse
-        ⩭(x, y) = relerrc(x, y) ≤ 5e-12
+        ⩭(x, y) = relerr(x, y) ≤ 5e-12
 
         @test expint(0, 1) ⩭ 1/e
         @test expint(0, 2) ⩭ 1/(2e^2)