desired precision for Jinc approximation

[NicolasRobidoux]

Anthony:

I'm going to break the approximation of Jinc in two pieces: |x|<=r3 (minimax approx. custom built by Chantal Racette and I), and |x|>r3, for which I'll use the standard (Hart's) method. (Here, r3 is the true zero of BesselJ1 before pi scaling, which is just a bit larger than 10.)

Do you really care about getting close to full double precision for Jinc evaluation?

For example, my preliminary results give max absolute error on [-r3,r3] less than 8e-09 (when the rational function is evaluated in double precision; in single precision, it's less than 9e-7 <- not of interest to IM anyway, but interesting to me because of SSE/Altivec/MMX/... and GPUs) at a cost of 23 flops (including one division, but excluding the squaring back of the square rooted distance, which adds one flop to the total).

I would guess that this is good enough.

Agreed?

(I can, of course, produce more precise approximations, but I think it's a waste of flops for IM.)

[NicolasRobidoux]

A nice reference about floating point approximation of the Bessel functions http://www.systomath.com/include/Boost- ... essel.html.

See, also, http://www.netlib.org/specfun/j1y1.

(I'm not using this for the [-r3,r3] approximation: I'm using minimax software.)

[NicolasRobidoux]

Another question for Anthony:

Are you OK with me getting rid of BesselOrderOne and approximating Jinc directly. For example, right now, BesselOrderOne is computed from J1 (Jinc itself) for small x by multiplying J1 by x, and then Jinc is computed by redividing by x.

Is BesselOrderOne used for anything else? (We could recover it from Jinc anyway

[anthony]

I don't think ultra high precision is needed.