Provide relative tolerances for besselj1(x) and besselj0(x) near zeros

[heltonmc]

besselj0(x) and besselj1(x) are less accurate than in some cases compared to SpecialFunctions. It is necessary to have as high as accuracy as possible here to use in recurrence.

A particularly egregious example is...

julia> x = 36.917175397041404
36.917175397041404
julia> b = Bessels.besselj0
besselj0 (generic function with 3 methods)
julia> b = Bessels.besselj0(x)
1.011768220172076e-5
julia> s = SpecialFunctions.besselj0(BigFloat(x))
1.011768220171316454254168807990064731641395604156791812684039845214461189126123e-05
julia> (1 - s/b)/eps(Float64)
3380.747393649221098313930426393567621377521727019041833762438852447240038176093

#SpecialFunctions
julia> sf = SpecialFunctions.besselj0(x)
1.0117682201713196e-5
julia> (1 - s/sf)/eps(Float64)
13.84518197098556328359063545573609602974659455973659014095473715553658825833736

besselj0 error

besselj1 error

[heltonmc]

It looks like this issue is limited to the asymptotic expansions for large arguments.
https://github.com/heltonmc/Bessels.jl/blob/ded7c2d1a401377df6752e5b4e2e0a2958f3974c/src/besselj.jl#L107-L129

[oscardssmith]

yeah. Accuracy near zeros of the function is somewhat fundamentally hard. I think the best answer might be to detect whether besselj(x) is near 0, and if it is, use an expansion around the zero.

[heltonmc]

Hmmm well we can approximate their zeros fairly quickly....

x - pi/4 - evalpoly(inv(x)^2, (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296, 24713030909/46137344, -7780757249041/436207616)) = (m-1/2)*pi

Where m is some integer... will need some care though...

[heltonmc]

Looking at this a bit more. It is probably best to absorb the pi/4 within the evaluation.

This problem can then be reduced to simply computing
(cos(x) * (cos(xn) + sin(xn)) + sin(x)*(cos(xn) - sin(xn)))
as accurately as possible.

An example close to a root is

(x, xn) = (46.341188371661815, -0.002696731212363751)

So then

julia> (cos(x) * (cos(xn) + sin(xn)) + sin(x)*(cos(xn) - sin(xn)))
-1.2212453270876722e-15

[heltonmc]

Or more simply we must calculate

cos(x + xn) + sin(x + xn)

more accurately. Near roots cos(x + xn) ≈ -sin(x + xn) is true and we get some cancellation. We might want to check this condition then do the cos(x + xn) + sin(x + xn) with higher precision or some expansion? Or have a better routine in general for that calculation.

[oscardssmith]

Isn't this just sin(x + pi/4 + xn)?

[heltonmc]

I've been trying to reconcile why the SpecialFunctions.jl implementation is better near the zeros. They suggest using this version. https://github.com/JuliaMath/openlibm/blob/0edf8d6929efc09c33e1e3342483a33dcd7517d1/src/e_j0f.c#L116-L124

But yes essentially we want to compute cos(x - pi/4 + xn) which can be expressed as inv(sqrt(2)*(cos(x+xn) + sin(x + xn) and then was trying to employ the suggestion in the the link for adding sin and cos. I haven't found anything to work much better but using inv(sqrt(2)*(cos(x+xn) + sin(x + xn) does improve accuracy slightly compared to cos(x - pi/4 + xn) but still struggles with a lot of cancellation

[heltonmc]

https://github.com/heltonmc/Bessels.jl/blob/ded7c2d1a401377df6752e5b4e2e0a2958f3974c/src/besselj.jl#L58-L63

So for our code replacing these lines with

a = SQ2OPI(T) * sqrt(xinv) * p
xn = xinv * q
b = 0.7071067811865475 * (Float64(cos(BigFloat(x) + BigFloat(xn)) + sin(BigFloat(x) + BigFloat(xn))))

fixes this issue but obviously that's cheating 😄

[heltonmc]

This is improved by #88. Which I think closes this issue.