Module jy 

The oscillatory Bessel functions $J_0$, $J_1$, $Y_0$, $Y_1$.

What “accuracy” means for a function with zeros

This is the one decision to make before any tolerance is written, and getting it wrong produces a test that is either impossible to pass or meaningless.

$J_\nu$ and $Y_\nu$ oscillate through zero forever. At a zero the relative error of any implementation is unbounded (the true value is 0 and the computed one is not), so a relative-error contract is not merely hard to meet, it is not a statement about anything. What every implementation actually delivers, and what this kernel promises, is accuracy relative to the envelope:

\left|\,\hat{f}(x) - f(x)\,\right| \;\lesssim\; C\,\varepsilon\,\sqrt{\frac{2}{\pi x}}

since $\sqrt{2/\pi x}$ is the amplitude the oscillation rides on. Equivalently: absolute error scaled by $\sqrt{x}$ is bounded. Tests here compare on that basis.

Below x = 8 there is a stronger guarantee. It is why the fits are shaped the way they are. Each sub-8 region carries one zero of the function, factored out as $(x + x_k)\left((x - x_{k1}/256) - x_{k2}\right)$: $x_{k1}/256$ is a power-of-two-scaled integer and therefore exact, so the subtraction near the root loses nothing and full relative accuracy survives at the first two or three zeros. Nobody does this above 8 (Boost included) because the number of zeros to factor grows without bound.

Regions

$J$ splits at 4 and 8, $Y_0$ at 3, 5.5 and 8, $Y_1$ at 4 and 8. Above 8 all four share the Hankel form: one amplitude pair in $(8/x)^2$ against sin x and cos x.

That last point is the whole reason this file exists rather than a port of fdlibm, whose j0f splits the asymptotic envelope alone into four sub-intervals with a rational apiece. A branch picks one and skips three. A vector unit evaluates all four and discards three. Boost’s single Hankel region is higher degree and strictly cheaper here.

Y calls J

$Y_\nu$ is singular at the origin, and the singularity is carried by a $\frac{2}{\pi}\ln(x/x_k)\,J_\nu(x)$ term rather than by the rational, so these kernels call the $J$ kernels, exactly as $K$ calls $I$. The log is taken about the region’s own root, not as a bare $\ln x$, which is what stops that term from swamping the rational near the zero.

Functions

bessel_j0_impl

$J_0(x)$. Even in x.

bessel_j1_impl

$J_1(x)$. Odd in x.

bessel_jn_pair_impl

(J_{N-1}, J_N) for N >= 2.

bessel_jv_impl

J_n(x) with a per-lane order.

bessel_y_impl

$Y_0(x)$ and $Y_1(x)$, selected by ORDER_ONE.

bessel_yn_recur

(Y_{N-1}, Y_N) by upward recurrence from the two closed forms.

bessel_yv_impl

Y_n(x) with a per-lane order.