flintlib / arb

Arb has been merged into FLINT -- use https://github.com/flintlib/flint/ instead
http://arblib.org/
GNU Lesser General Public License v2.1
456 stars 137 forks source link

Sudden Loss/Reset of Precision in Bessel Functions #360

Open fingolfin opened 3 years ago

fingolfin commented 3 years ago

On https://github.com/Nemocas/Nemo.jl/issues/954 there is a report for some loss/reset of precision in the besselk, bessely and besselj (while besseli seems to be unaffected).

In short, besselk(2.0, 57.9668555791947111700112064539202300037282162) is evaluated for various precisions; notice the jump from precision 110 to 120:

(p, besselk) = (100, 1.13687286812830688844382559e-26)
(p, besselk) = (110, 1.13687286812830688844382558995e-26)
(p, besselk) = (120, [+/- 1.19e-9])
(p, besselk) = (130, [+/- 1.17e-12])
(p, besselk) = (140, [+/- 1.18e-15])
(p, besselk) = (150, [+/- 1.15e-18])
(p, besselk) = (160, [+/- 1.14e-21])
(p, besselk) = (170, [+/- 1.13e-24])
(p, besselk) = (180, 1e-26)
(p, besselk) = (190, 1.137e-26)
(p, besselk) = (200, 1.136873e-26)
(p, besselk) = (210, 1.136872868e-26)
(p, besselk) = (220, 1.136872868128e-26)
(p, besselk) = (230, 1.136872868128307e-26)
(p, besselk) = (240, 1.136872868128306888e-26)
(p, besselk) = (250, 1.136872868128306888444e-26)

Further examples for bessely and besselj can be found at the Nemo issue I linked above.

fredrik-johansson commented 3 years ago

Same problem as #228 (and others): there is a suboptimal cutoff between the expansion at zero and the asymptotic expansion at infinity.

tthsqe12 commented 2 years ago

Honestly, I think you will end up implementing the Laplace integral for U(z) when z is too small for the reciprocal transform, but not so small that the asymptotic series can be usefully truncated.

fredrik-johansson commented 2 years ago

This example now works for Bessel K thanks to falling back to numerical integration:

10   1e-26
20   1.137e-26
30   1.136873e-26
40   1.136872868e-26
50   1.136872868128e-26
60   1.136872868128307e-26
70   1.136872868128306888e-26
80   1.136872868128306888444e-26
90   1.136872868128306888443826e-26
100   1.136872868128306888443825590e-26
110   1.136872868128306888443825589950e-26
120   1.136872868128306888443825589950094e-26
130   1.136872868128306888443825589950094325e-26
140   1.136872868128306888443825589950094324651e-26
150   1.13687286812830688844382558995009432465094e-26
160   1.136872868128306888443825589950094324650935557e-26
170   1.136872868128306888443825589950094324650935557413e-26
180   1.136872868128306888443825589950094324650935557413225e-26
190   1.136872868128306888443825589950094324650935557413224860e-26
200   1.136872868128306888443825589950094324650935557413224859798e-26
210   1.136872868128306888443825589950094324650935557413224859797745e-26
220   1.136872868128306888443825589950094324650935557413224859797744589e-26
230   1.136872868128306888443825589950094324650935557413224859797744589545e-26
240   1.136872868128306888443825589950094324650935557413224859797744589544662e-26
250   1.136872868128306888443825589950094324650935557413224859797744589544661843e-26

The general problem is far from solved though (other Bessel functions; nonreal arguments).