Adding jinc function

[roflmaostc]

Hey,

we already discussed a bit in #264 about adding a jinc function.
We decided on the following convention, which is similar to sinc in terms of the peak height and argument scaling:

I wanted to use fastabs which is only available in the current master of julia. This is why I copied the definition to the if statement.

I would be happy to receive some feedback!

Thanks,

Felix

(Closes #264.)

[roflmaostc]

Error for 1.3: UndefVarError: evalpoly not defined

evalpoly needs julia >= 1.4. What should we do about that?

[codecov]

Codecov Report

Merging #266 (eb4473f) into master (9116108) will increase coverage by 0.08%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #266      +/-   ##
==========================================
+ Coverage   87.59%   87.68%   +0.08%     
==========================================
  Files          11       11              
  Lines        2508     2526      +18     
==========================================
+ Hits         2197     2215      +18     
  Misses        311      311              

Flag	Coverage Δ
unittests	87.68% <100.00%> (+0.08%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
src/SpecialFunctions.jl	100.00% <100.00%> (ø)
src/bessel.jl	93.48% <100.00%> (+0.25%)	⬆️
src/expint.jl	95.25% <100.00%> (+0.09%)	⬆️
src/gamma.jl	93.83% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 9116108...eb4473f. Read the comment docs.

[roflmaostc]

Just found out, that www.wolframalpha seems to struggle provide high precision results...
However, Mathematica produces similar results as Julia :D
Seems to be Julia :/

[roflmaostc]

I could track down the problem

Mathematica:

BesselJ[1, 1.8`200]

Julia:

besselj1(big(1.8))

Outputs:

0.58151695173116518347008479715619083566653627328448241718624528992169
0.581516951731165184221547593970460741118046297653496168575818368562848497914668

What should we do? besselj1 seems to be inaccurate 😿

[roflmaostc]

Maybe not good style, but the test cases I pushed fail. I'm not sure what's going on but the problem is probably somewhere in besselj1.

Does someone know why Mathematica (and some webpages) show different results and is it expected?

[stevengj]

You want

julia> besselj1(big"1.8")
0.5815169517311651834700847971561908356665362732844824171862452899216979041922363

which matches your Mathematica result.

The problem is that big(1.8) is not the number you are thinking of:

julia> big(1.8)
1.8000000000000000444089209850062616169452667236328125

That is because 1.8 is a Float64 literal, but 1.8 is not exactly representable as a binary floating-point number so it is rounded to something that differs from 1.8 in the 17th decimal place. In other words, big(1.8) means round 1.8 to the nearest Float64 and then convert to BigFloat.

[roflmaostc]

Sorry to bother. Your hint fixed many cases, but not all.
I tried to make sure that it's parsed correctly

julia> res2 = big"1.971105539196620210542904558564976891448851717386016064466205835417372529649732623778146798009632930245313240967710133352055225553347549825667896762467087506209685330050397030027682762563950964962993e11" - big"5.89917589196246805438832769814196828580119495400655210987269285701113079909408122244782544425726443894463237062342568390611353700647229983489787030738912900107228049806443988000825445873073676204e6"*1im
1.971105539196620210542904558564976891448851717386016064466205835417372529649727e+11 - 5.899175891962468054388327698141968285801194954006552109872692857011130799094111e+06im

julia> Complex{Float64}(res2)
1.9711055391966202e11 - 5.899175891962468e6im

julia> x2 = big"0.00001" + big"10"*1im
9.99999999999999999999999999999999999999999999999999999999999999999999999999994e-06 + 10.0im

julia> jinc(Complex{Float64}(x2))
1.9711055391966183e11 - 5.899175891974533e6im

However, the result still mismatches. Any further idea?

Again, Mathematica vs julia for besselj1:

BesselJ[1, 0.00001`200 + I*10.0`200]

0.02548617798056962930423261011063522051286601522259212166836850560273\
5224000250169633431820066537728081112783926189116335069422934386705411\
0177843126480606543290913529010319815371330493594661136321 + 
 2670.9883035791128339955440167296286794317890303094196461232014839508\
3287658223018716591008155410056811500340138735313841171781910267075143\
88244430315094098180814058360987122761051848662031490035143561 I

julia> besselj1(Complex{Float64}(x2))
0.025486177980733188 + 2670.988303579112im

[stevengj]

However, the result still mismatches. Any further idea?

z = besselj1(Complex{Float64}(x2)) matches the Mathematica result z₀ with relative error 1.81e-16, i.e. nearly to machine precision, as measured in the complex plane as abs(z-z₀)/abs(z₀). The relative error in the real part by itself, however is 6.41e-12 while the relative error in the (much larger) imaginary part by itself is 1.70e-16.

Because the error in the real part by itself is > 1e-13, it will fail the ≅ equality test. But that's not a big cause for concern — ≅ equality is sometimes too stringent. It's quite common for complex-valued special functions to be accurate only in the sense of abs(z-z₀)/abs(z₀) being small, not for the real and imaginary parts individually, and ≅ equality sometimes fails when the one of the two parts is much smaller than the other (i.e. close to the real or imaginary axes). Apparently the Bessel functions that we are using fall into this category.

[roflmaostc]

Ok, got that.

I relaxed \cong for complex number to check for abs accuracy if real and imag accuracy fails.

[stevengj]

I relaxed \cong for complex number to check for abs accuracy if real and imag accuracy fails.

No, please don't do that! Please revert this change.

In cases where the stronger accuracy condition holds, we want to ensure that it continues to hold!

Only change the tests for the cases that don't currently support this level of accuracy. In that case, you don't need a new operator at all — you can just use ≈ with an rtol argument, e.g.

@test z ≈ z₀ rtol=1e-13

[stevengj]

Looks good. Now you just need to add it to the documentation (at the end of the list of Bessel functions) here and here.