Fast and Accurate Bessel Function Computation

John Harrison Intel Corporation, JF1-13 2111 NE 25th Avenue Hillsboro OR 97124, USA

                                 Email: johnh@ichips.intel.com

Abstract

The Bessel functions are considered relatively difficult to compute. Although they have a simple power series expansion that is everywhere convergent, they exhibit approximately periodic behavior which makes the direct use of the power series impractically slow and numerically unstable. We describe an alternative method based on systematic expansion around the zeros, refining existing techniques based on Hankel expansions, which mostly avoids the use of multiprecision arithmetic while yielding accurate results.

1. Introduction and overview

Bessel functions are certain canonical solutions to the differential equations

x^{2} \frac{d^{2}y}{dx^{2}} + x\frac{dy}{dx} + (x^{2} − n^{2})y = 0

We will consider only the case where n is an integer. The canonical solutions considered are the Bessel functions of the first kind, J_{n}(x), nonsingular at x = 0, and those of the second kind, Y_{n}(x), which are singular there. In each case, the integer n is referred to as the order of the Bessel function. Figure 1 shows a plot of J_{0}(x) and J_{1}(x) near the origin, while Figure 2 is a similar plot for Y_{0}(x) and Y_{1}(x).

The Bessel functions J_{n}(x) have power series that are convergent everywhere, with better convergence than the familiar series for the exponential or trigonometric functions:

J_{n}(x) = ∑_{m=0}^{∞} \frac{(−1)^{m}(x/2)^{n+2m}}{m!(n + m)!}

However, the direct use of the power series would require too many terms for large x, and even for moderate x is likely to be quite numerically unstable close to the zeros. The difficulty is that the final value J_{n}(x) can be small for large x even when the intermediate terms of the power series are large. The trigonometric functions like sin(x) also have this property, but they are still quite easy because they are exactly periodic. The Bessel functions are not quite periodic, though they do start to look more and more like scaled trigonometric functions for large x, roughly speaking:^{1}

J_{n}(x) ≈ \sqrt{\frac{2}{πx}} cos(x − [n/2 + 1/4]π)

Y_{n}(x) ≈ \sqrt{\frac{2}{πx}} sin(x − [n/2 + 1/4]π)

For extensive detail on the theory of the Bessel functions, as well as a little history and explanation of how they arise in physical applications, the reader is referred to Watson’s monograph [9].

The not-quite periodicity has led to some pessimism about the prospects of computing the Bessel functions with the same kind of relative accuracy guarantees as for most elementary transcendental functions. For example Hart et al.

[3] say:

However, because of the large number of zeros of these functions, it is impractical to construct minimum relative error subroutines, and the relative error is likely to be unbounded in the neighborhood of the zeros.

However, it is important to remember that this was written at a time when giving good relative accuracy guarantees even on the basic trigonometric functions would have been considered impractical too. We shall see that, at least for specific order and specific floating-point precisions, producing results with good relative accuracy is not so very difficult. In this paper we focus exclusively on the functions J_{0}, J_{1}, Y_{0} and Y_{1} in double-precision (binary64) floating-point arithmetic. The results should generalize straightforwardly to other specific Bessel functions of integer order and other floating-point formats, but the various parameters, polynomial degrees, domain bounds, worst-case results etc. would need to be computed afresh for each such instance.

Our general approach

The main theme in our proposed approach is to expand each function about its zeros. More precisely, we want to

1. These are not, properly speaking, asymptotic results ∼, since the zeros of the two functions do not coincide exactly. The exact relationship will be made precise below.

-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0  5  10  15  20 J_{0}(x) J_{1}(x)

Figure 1. Bessel function of the first kind, J_{0} and J_{1}

-3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 0  5  10  15  20 Y_{0}(x) Y_{1}(x)

Figure 2. Bessel function of the second kind, Y_{0} and Y_{1}

formulate the algorithms to move the inevitable cancellation forward in the computation to a point before there are rounding errors to be magnified. For example, if the input argument x is close to a zero z, we want to, in effect, compute x − z accurately at once and use that value in subsequent stages.

If we can successfully realize this idea, we can expect accurate results even close to zeros without performing all the intermediate computations in very high precision. Still, in order to assess the accuracy required we need a rigorous examination of the zeros to see how close they may be to double-precision floating-point numbers, and this will also form part of the present document. We will consider in turn:

• Evaluation near 0 for the singular Y_{n}

• Evaluation for ‘small’ arguments, roughly |x| < 45, but away from the singularities of the Y_{n} at zero

• Evaluation for ‘large’ arguments, roughly |x| ≥ 45.

2. Y_{n} near singularities

The Y_{n} have various singularities at x = 0, including one of the form log(x), so it doesn’t seem practical to use polynomial or rational approximations in this neighborhood. However, if we incorporate the logarithm and perhaps reciprocal powers, tractable approximations are available. For example we have (see [9] §3.51):

Y_{0}(x) = \frac{2}{π} ( \begin{matrix}[γ + log(x/2)]J_{0}(x)− \\ ∑_{m=1}^{∞} \frac{(−1)^{m}(x/2)^{2m}}{(m!)^{2}} [1 + \frac{1}{2} + ⋯ + \frac{1}{m}]) \end{matrix}

If we use general floating-point coefficients, we can write this as follows, using two even power series W_{0}(x) and Z_{0}(x):

Y_{0}(x) = W_{0}(x) log(x) − Z_{0}(x)

A similar but slightly more complicated expansion also works for Y_{1}(x), in which case the corresponding power series W_{1}(x) and Z_{1}(x) are odd:

Y_{1}(x) = W_{1}(x) log(x) − Z_{1}(x) − \frac{2}{πx}

For the moderate ranges required, minimax economizations of the W_{i}(x) and Z_{i}(x) only require about degree 16 with only alternate terms present for double-precision results.

3. Small arguments

For fairly small arguments, say |x| < 45, we can take the idea of expanding around zeros at face value. We don’t need to do any kind of computation of zeros at runtime, but can just pre-tabulate the zeros z_{1}, z_{2}, … z_{N}, since there aren’t so many in range. Then at runtime, given an input x we start by reducing the input argument x to r = x − z_{k}, where z_{k} is approximately the closest zero to x, and then use an expansion in terms of r. Since even for small arguments the zeros are spaced fairly regularly at intervals of about π — see Table 1 — we can quite efficiently find an appropriate k by a simple division (or more precisely, reciprocal multiplication) and scaling to an integer. If we store z_{k} in two pieces and subtract the high part first, we get an accurate reduced argument

r = (x − z_{k}^{hi}) − z_{k}^{lo}

We use double-extended precision for the intermediate result, which ensures that the first subtraction is exact. (Since our estimate of the closest zero at the low end can be a little wayward, it is not immediately clear that we could rely on Sterbenz’s lemma here if we used double-precision throughout.)

In fact, we tabulate both zeros and extrema among the z_{k}. This roughly halves the size of each interval of approximation to about −π/4 ≤ r ≤ π/4 and so permits us to use polynomials of lower degree. It also removes potential concerns over monotonicity near the extrema where the successive intervals of approximation fit together.

We then use a separate polynomial p_{k}(r) for each k, with each one precomputed. About most z_{k}, we can achieve accuracy adequate for double-precision results using polynomials p_{k}(r) of degree 14. The exceptions are the small zeros of the Y_{n} where the nearness of the singularities causes problems. Even here, provided in the case of the smallest zero we cut the range off well away from the origin, the polynomials required have tractable degree. Table 2 summarizes the degree of polynomials required to achieve relative errors suitable for the various floating-point precisions for the four functions J_{0}, J_{1}, Y_{0} and Y_{1}. We cut off the range for the first zero of Y_{0} at r = −0.1.

4. Asymptotic expansions

If we turn to larger arguments with about |x| ≥ 45, the approach traditionally advocated is based on Hankel’s asymptotic expansions. The Bessel functions can be expressed as

J_{n}(x) = \sqrt{\frac{2}{πx}}( \begin{matrix}cos(x − [n/2 + 1/4]π) · P_{n}(x)− \\ sin(x − [n/2 + 1/4]π) · Q_{n}(x)) \end{matrix}

and

Y_{n}(x) = \sqrt{\frac{2}{πx}}( \begin{matrix}sin(x − [n/2 + 1/4]π) · P_{n}(x)+ \\ cos(x − [n/2 + 1/4]π) · Q_{n}(x)) \end{matrix}

where the auxiliary functions P_{n}(x) and Q_{n}(x) may, for example, be expressed as integrals — see [9] §7.2. For reasonably large values of x these auxiliary functions can be well approximated by asymptotic expansions:

P_{n}(x) ∼ ∑_{m=0}^{∞} \frac{(−1)^{m}(n, 2m)}{(2x)^{2m}}

Q_{n}(x) ∼ ∑_{m=0}^{∞} \frac{(−1)^{m}(n, 2m + 1)}{(2x)^{2m+1}}

where the notation (n, m) denotes:

(n, m) = \frac{(4n^{2} − 1^{2})(4n^{2} − 3^{2}) ⋯ (4n^{2} − [2m − 1]^{2})}{2^{2m}m!}

For example, we have

J_{0}(x) = \sqrt{\frac{2}{πx}} [cos(x − π/4)P_{0}(x) − sin(x − π/4)Q_{0}(x)]

where

P_{0}(x) ∼ 1−\frac{9}{128x^{2}} +\frac{3675}{32768x^{4}} −\frac{2401245}{4194304x^{6}} +\frac{13043905875}{2147483648x^{8}} −⋯

and

Q_{0}(x) ∼ −\frac{1}{8x} + \frac{75}{1024x^{3}} − \frac{59535}{262144x^{5}} + \frac{57972915}{33554432x^{7}} − ⋯

\begin{matrix}J_{0} & J_{1} & Y_{0} & Y_{1} \\ 2.404825 & 0.000000 & 0.893576 & 2.197141 \\ 5.520078 & 3.831705 & 3.957678 & 5.429681 \\ 8.653727 & 7.015586 & 7.086051 & 8.596005 \\ 11.791534 & 10.173468 & 10.222345 & 11.749154 \\ 14.930917 & 13.323691 & 13.361097 & 14.897442 \\ 18.071063 & 16.470630 & 16.500922 & 18.043402 \\ 21.211636 & 19.615858 & 19.641309 & 21.188068 \\ 24.352471 & 22.760084 & 22.782028 & 24.331942 \\ 27.493479 & 25.903672 & 25.922957 & 27.475294 \\ 30.634606 & 29.046828 & 29.064030 & 30.618286 \\ 33.775820 & 32.189679 & 32.205204 & 33.761017 \\ 36.917098 & 35.332307 & 35.346452 & 36.903555 \\ 40.058425 & 38.474766 & 38.487756 & 40.045944 \\ 43.199791 & 41.617094 & 41.629104 & 43.188218 \\ 46.341188 & 44.759318 & 44.770486 & 46.330399 \\ 49.482609 & 47.901460 & 47.911896 & 49.472505 \\ 52.624051 & 51.043535 & 51.053328 & 52.614550 \\ 55.765510 & 54.185553 & 54.194779 & 55.756544 \\ 58.906983 & 57.327525 & 57.336245 & 58.898496 \\ 62.048469 & 60.469457 & 60.477725 & 62.040411 \\ 65.189964 & 63.611356 & 63.619215 & 65.182295 \\ 68.331469 & 66.753226 & 66.760716 & 68.324152 \\ 71.472981 & 69.895071 & 69.902224 & 71.465986 \\ 74.614500 & 73.036895 & 73.043740 & 74.607799 \end{matrix}

Table 1. The approximate values of first zeros of J_{0}, J_{1}, Y_{0} and Y_{1}

Precision Relative error Typical degree Worst-case degree
Single  2^{−31}  8-9 (all but 3)  12
Double  2^{−60}  14 (all but 6)  23
Extended  2^{−71}  16 (all but 7)  28

Table 2. Degrees of polynomials needed around small zeros and extrema

Although the series P_{n}(x) and Q_{n}(x) actually diverge for all x, they are valid asymptotic expansions, meaning that at whatever point m we truncate the series to give P_{n}^{m}(x), we can make the relative error as small as we wish by making x sufficiently large:

lim_{x→∞} \frac{|P_{n}^{m}(x) − P_{n}(x)|}{|P_{n}(x)|} = 1

In fact, it turns out — see [9] §7.32 — that the error in truncating the series for P_{n}(x) and Q_{n}(x) at a given term is no more than the value of the first term neglected and has the same sign, provided 2p ≥ n − 1/2 where the term in x^{−p} is the last one included

The Hankel formulas are much more appropriate for large x because they allow us to exploit periodicity directly. However they are still numerically unstable near the zeros, because the two terms cos(⋯)P_{n}(x) and sin(⋯)Q_{n}(x) may cancel substantially. This means that we need to compute the trigonometric functions in significantly more than working precision in order to obtain accurate results. In accordance with our general theme, we will try to modify the asymptotic expansions to move the cancellation forward.

Modified expansions

The main novelty in our approach to computing the Bessel functions for large arguments is to expand either J_{n}(x) or Y_{n}(x) in terms of a single trigonometric function with a further perturbation of its input argument as follows:

J_{n}(x) = \sqrt{\frac{2}{πx}} · β_{n}(x) cos(x − [n/2 + 1/4]π − α_{n}(x))

Y_{n}(x) = \sqrt{\frac{2}{πx}} · β_{n}(x) sin(x − [n/2 + 1/4]π − α_{n}(x))

where α_{n}(x) and β_{n}(x) are chosen appropriately. These equations serve to define α_{n}(x) and β_{n}(x) in terms of the Bessel functions by the following, where the ≡ sign indicates congruence modulo π, i.e. ignoring integer multiples of π. (In some of the plots we eliminate spurious π multiples by successively applying tan then atan to such expressions.)

α_{n}(x) ≡ x − [n/2 + 1/4]π − atan(Y_{n}(x)/J_{n}(x))
≡ x − [n/2 + 3/4]π + atan(J_{n}(x)/Y_{n}(x))

and

β_{n}(x) = \sqrt{\frac{πx}{2}(J_{n}(x)^{2} + Y_{n}(x)^{2})}

We will next obtain asymptotic expansions for the α_{n}(x) and β_{n}(x) starting from those for P_{n}(x) and Q_{n}(x). Our derivations are purely formal, but we believe it should be possible to prove rigorously that the error on truncating the expansions is bounded by the first term neglected and has the same sign — see for example the expansion of the closely related ‘modulus’ and ‘phase’ functions in §9.2.31 of [1]. Using the addition formula we obtain:

J_{n}(x) = \sqrt{\frac{2}{πx}}·β_{n}(x)[ \begin{matrix}cos(x − [n/2 + 1/4]π) cos(α_{n}(x))+ \\ sin(x − [n/2 + 1/4]π) sin(α_{n}(x))] \end{matrix}

Comparing coefficients, we see

β_{n}(x) cos(α_{n}(x)) = P_{n}(x)

and β_{n}(x) sin(α_{n}(x)) = −Q_{n}(x)

These imply that we should choose β_{n}(x) as follows:

β_{n}(x)^{2} = β_{n}(x)^{2}[cos(α_{n})^{2} + sin(α_{n})^{2}]
= (β_{n}(x) cos(α_{n}))^{2} + (β_{n}(x) sin(α_{n}))^{2}
= P_{n}(x)^{2} + Q_{n}(x)^{2}

If we assume β_{n}(x) is nonzero even at the zeros (and it will be since Q_{n}(x) ≈ 1 for the large x we are interested in), we can also immediately obtain tan(α_{n}(x)) = −Q_{n}(x)/P_{n}(x). We can carry through the appropriate computations on formal power series — see §15.52 of [9] for more rigorous proofs of some related results. We have:

tan(α_{0}(x)) = −Q_{0}(x)/P_{0}(x)
= \frac{1}{8x} − \frac{33}{512x^{3}} + \frac{3417}{16384x^{5}} − \frac{3427317}{2097152x^{7}} + ⋯

and composing this with the series for the arctangent function, we obtain

α_{0}(x) = \frac{1}{8x}−\frac{25}{384x^{3}} +\frac{1073}{5120x^{5}} −\frac{375733}{229376x^{7}} +\frac{55384775}{2359296x^{9}} −⋯

Similarly, using β_{0}(x) = \sqrt{P_{0}(x)^{2} + Q_{0}(x)^{2}} we get

β_{0}(x) = 1− \frac{1}{16x^{2}} + \frac{53}{512x^{4}} − \frac{4447}{8192x^{6}} + \frac{3066403}{524288x^{8}} −· · · · · ·

In exactly the same way we get:

α_{1}(x) = −\frac{3}{8x}+\frac{21}{128x^{3}} −\frac{1899}{5120x^{5}} +\frac{543483}{229376x^{7}} −\frac{8027901}{262144x^{9}} +⋯

and

β_{1}(x) = 1 + \frac{3}{16x^{2}} − \frac{99}{512x^{4}} + \frac{6597}{8192x^{6}} − \frac{4057965}{524288x^{8}} + ⋯

For moderate x, the asymptotic series quickly become good approximations, as shown in Figures 3 and 4, though to get the high accuracy we demand, somewhat larger x and many more terms of the series are necessary, explaining our chosen cutoff point of |x| ≥ 45. Nevertheless, we can and do utilize a minimax economization of the expansion, which considerably reduces the required degree.

Watson [9] presents (around pp. 213–4) a discussion of some formulas due to Stieltjes for estimating the remainders on truncating the original expansions, and presumably similar analysis could be applied to the α_{n}(x) and β_{n}(x).

The modified expansions immediately give rise to relatively simple computation patterns by truncating the series appropriately, e.g.

J_{0}(x) ≈ \sqrt{\frac{2}{πx}}(1− \frac{1}{16x^{2}} + \frac{53}{512x^{4}} ) cos(x− \frac{π}{4} − \frac{1}{8x} + \frac{25}{384x^{3}} )

Y_{0}(x) ≈ \sqrt{\frac{2}{πx}}(1− \frac{1}{16x^{2}} + \frac{53}{512x^{4}} ) sin(x− \frac{π}{4} − \frac{1}{8x} + \frac{25}{384x^{3}} )

J_{1}(x) ≈ \sqrt{\frac{2}{πx}}(1+\frac{3}{16x^{2}} −\frac{99}{512x^{4}} ) cos(x−\frac{3π}{4} +\frac{3}{8x}−\frac{21}{128x^{3}} )

Y_{1}(x) ≈ \sqrt{\frac{2}{πx}}(1+\frac{3}{16x^{2}} −\frac{99}{512x^{4}} ) sin(x−\frac{3π}{4} +\frac{3}{8x}−\frac{21}{128x^{3}} )

This seems a much more promising approach because numerical cancellation only needs to be dealt with in the simple algebraic expression x − α_{n}(x), integrated with standard trigonometric range reduction. The other components of the final answer, a trigonometric function, inverse square root and β_{n}(x) are all well-behaved. The overall computational burden is therefore not much worse than with the trigonometric functions.

We now need to proceed with a careful analysis to show the range limits of this technique and the points at which we can afford to truncate the series. Since the value of β_{n}(x) is approximately 1, the absolute error in it contributes directly to a relative error in the overall result. In the case of α_{n}(x) however, a given absolute error may translate into a much larger relative error if the result is close to zero. In order to place a bound on how much the error can blow up, we need to know how close the large zeros may come to floatingpoint numbers.

5. Worst-case zeros

The zeros of the Bessel functions are spaced about π from each other. Assuming that their low-order bits are randomly distributed, we would expect the closest one of the zeros would come to a precision-p floating-point number would be about 2^{−(p+log_{2} p)}, say 2^{−60} for double precision [6]. Since at least the work of Kahan (see the nearpi program) it has been known that this is indeed the case for pure multiples of π, and we might expect similar results here. While from a practical point of view it seems safe to make this kind of naive assumption, and perhaps leave a little margin of safety, a more satisfactory approach is to verify this fact rigorously. We will do so for double precision, our main focus.

0 0.05 0.1 0.15 0.2 0.25 0.3 1  1.5  2  2.5  3  3.5  4  4.5  5 α_{0}(x) = atan(tan(x − π/4 − atan(Y_{0}(x)/J_{0}(x)))) 3 terms of asymptotic expansion

Figure 3. Plot of α_{0}(x) for moderate x versus approximation

0.8 0.85 0.9 0.95 1 1.05 1.1 1  1.5  2  2.5  3  3.5  4  4.5  5 β_{0}(x) = \sqrt{\frac{πx}{2}(J_{0}(x)^{2} + Y_{0}(x)^{2})} 3 terms of asymptotic expansion

Figure 4. Plot of β_{0}(x) for moderate x versus approximation

Locating the zeros

We have already accurately computed the small zeros of the Bessel functions of interest, since these were needed to produce the various expansions used at the low end. (All these were computed by straightforward use of the infinite series, which can be computationally lengthy but is quite straightforward.) We can therefore simply exhaustively check the small zeros to see how close they come to doubleprecision numbers.

For larger values of x we use a technique originally due to Stokes — see §15.52 of [9]. We can start with our functions α_{n}(x) (which were considered by Stokes in this context, though not as part of a computational method) and perform some further power series manipulations to obtain expansions for the various roots. For example, since

J_{0}(x) = \sqrt{\frac{2}{πx}} · β_{0}(x) cos(x − π/4 − α_{0}(x))

we seek values for which the cosine term here is zero, i.e. where for some integer m we have

x − π/4 − α_{0}(x) = (m − 1/2)π

or if we write p = (m − 1/4)π:

x − α_{0}(x) = p

Using the asymptotic expansion we can write:

x − \frac{1}{8x} + \frac{25}{384x^{3}} − \frac{1073}{5120x^{5}} + \frac{375733}{229376x^{7}} − ⋯ = p

or

x (1 − \frac{1}{8x^{2}} + \frac{25}{384x^{4}} − \frac{1073}{5120x^{6}} + \frac{375733}{229376x^{8}} − ⋯) = p

Multiplying both sides by the inverse of the series 1 − \frac{1}{8x^{2}} + \frac{25}{384x^{4}} − ⋯ we get:

x = p (1 + \frac{1}{8x^{2}} − \frac{19}{384x^{4}} + \frac{2999}{15360x^{6}} − \frac{16352423}{10321920x^{8}} + ⋯)

and therefore:

\frac{1}{p} = \frac{1}{x} + \frac{1}{8x^{3}} − \frac{19}{384x^{5}} + \frac{2999}{15360x^{7}} − \frac{16352423}{10321920x^{9}} + ⋯

By reversion of the power series we obtain:

\frac{1}{x} = \frac{1}{p} − \frac{1}{8p^{3}} + \frac{37}{384p^{5}} − \frac{1373}{5120p^{7}} + \frac{19575433}{10321920p^{9}} − ⋯

so

\frac{p}{x} = 1 − \frac{1}{8p^{2}} + \frac{37}{384p^{4}} − \frac{1373}{5120p^{6}} + \frac{19575433}{10321920p^{8}} − ⋯

Inverting the series again we get

\frac{x}{p} = 1 + \frac{1}{8p^{2}} − \frac{31}{384p^{4}} + \frac{3779}{15360p^{6}} − \frac{6277237}{3440640p^{8}} + ⋯

and so we finally obtain x in terms of p = (m − 1/4)π:

x = p + \frac{1}{8p} − \frac{31}{384p^{3}} + \frac{3779}{15360p^{5}} − \frac{6277237}{3440640p^{7}} + ⋯

These power series computations are quite tedious by hand but can be handled automatically in some computer algebra systems. We used our own implementation of power series in the OCaml language via higher-order functions [5], and with this we can obtain any reasonable truncation of these series automatically in an acceptable time. Needless to say, our results are consistent with the hand computations in [9], though often we need more terms than are given there. Zeros of the other Bessel functions are obtained in the same way. For Y_{0}(x) the expansion is the same except that we use p = (m + 1/4)π, and we have the following expansion for the zeros of J_{1}(x) and Y_{1}(x):

x = p−\frac{3}{8p}+\frac{3}{128p^{3}} −\frac{1179}{5120p^{5}} +\frac{1951209}{1146880p^{7}} −\frac{223791831}{9175040p^{9}}+⋯

where p = (m + 1/4)π for J_{1}(x) and p = (m − 1/4)π for Y_{1}(x).

For moderate m, say up to a few thousand, million, or even billion, we can examine the zeros exhaustively using a suitable truncation of this formula to approximate them. And for m > 2^{70} or so, all terms but the first are sufficiently small in magnitude that we can reduce matters to the well-known linear case, which has already been analyzed for trigonometric range reduction. In the middle range, we have applied the Lefèvre-Muller technique [4] to the Stokes expansions. The idea of this technique is to cover the range with a large (but quite feasible) number of linear approximations that are accurate enough for us to be able to read off local results from. We then find ‘worst cases’ for the linear approximations in a straightforward manner by scaling the function so its gradient is approximated by a suitable rational approximation a/b, then solving congruences mod b.

Results

Table 3 shows all the double-precision floating-point numbers in the range 0 ≤ x ≤ 2^{90} that are within 2^{−55} of a zero of one of our four basic functions J_{0}(x), J_{1}(x), Y_{0}(x) and Y_{1}(x). As noted above, there are no surprises beyond 2^{90}, since the zeros are so close to (n − 1/4)π that we know roughly what to expect from previous studies on the trigonometric functions — see e.g. [8] and [6]. (Note that if we have results for double-precision numbers close to nπ, the results for multiples nπ/4 and a fortiori (n − 1/4)π can be at worst 1/4 as large.) Indeed, the larger zeros of J_{0} and Y_{1}, and those of J_{1} and Y_{0}, are already becoming very close near the top of our range. Still, for completeness, Table 4 lists the closest zeros of the functions to double-precision numbers over the whole range. The overall lesson is that the results are in line with naive statistical expectations.

6. A reference implementation

To realize the ideas set out above, we have programmed a reference implementation of the double-precision Bessel functions J_{0}, J_{1}, Y_{0} and Y_{1} designed for the Intel® Itanium® architecture. This architecture allows us to use double-extended precision for internal calculations, so most of the computation can be done ‘naively’ with good final error bounds. The sole exception is the accurate computation of α_{n}(x) and its integration with trigonometric range reduction. Here we make good use of the fma (fused multiplyadd) that this architecture offers.

The higher-order part of the α_{n}(x) series, consisting of terms that are still below about 2^{−60} in size, can also be calculated naively. But for the low-order part of the series we want to compute the summation in two double-extended pieces to ensure absolute accuracy of order 2^{−120}. One reasonable way of doing this is to make repeated use of

Exact value of double Approximate decimal value Distance from zero Function
2^{14} × 6617649673795284  1.08423572255 × 10^{20}  2^{−58.438199858}  Y_{1}
2^{37} × 6311013172270677  8.67379045745 × 10^{26}  2^{−58.4361612221}  Y_{1}
2^{37} × 6311013172270677  8.67379045745 × 10^{26}  2^{−58.4361612218}  J_{0}
2^{14} × 6617649673795284  1.08423572255 × 10^{20}  2^{−58.4356039989}  J_{0}
2^{−14} × 4643410941512688  283411312348.0  2^{−57.4750503229}  J_{0}
2^{23} × 7209129755475690  6.04745635398 × 10^{22}  2^{−57.2168697228}  J_{1}
2^{23} × 7209129755475690  6.04745635398 × 10^{22}  2^{−57.216867725}  Y_{0}
2^{28} × 8451279557623493  2.26862308183 × 10^{24}  2^{−57.1898493381}  Y_{1}
2^{28} × 8451279557623493  2.26862308183 × 10^{24}  2^{−57.1898492858}  J_{0}
2^{−26} × 8368094255856943  124694321.392  2^{−56.9469283257}  Y_{0}
2^{16} × 4963237255346463  3.25270716766 × 10^{20}  2^{−56.8515062657}  J_{1}
2^{16} × 4963237255346463  3.25270716766 × 10^{20}  2^{−56.8512179523}  Y_{0}
2^{32} × 8754199225116346  3.75989993745 × 10^{25}  2^{−56.8148254699}  Y_{1}
2^{32} × 8754199225116346  3.75989993745 × 10^{25}  2^{−56.8148254675}  J_{0}
2^{−18} × 7757980709970194  29594347801.1  2^{−56.6251590484}  Y_{0}
2^{−42} × 5944707359537560  1351.66996177  2^{−56.5865796184}  J_{1}
2^{−22} × 5798262669118148  1382413546.83  2^{−56.3577856958}  J_{0}
2^{−30} × 4670568619103095  4349805.99126  2^{−56.1575895598}  J_{1}
2^{16} × 8272062092244105  5.42117861277 × 10^{20}  2^{−56.1144022739}  Y_{1}
2^{16} × 8272062092244105  5.42117861277 × 10^{20}  2^{−56.1142984844}  J_{0}
2^{0} × 6027843377079719  6.02784337708 × 10^{15}  2^{−55.9891793419}  J_{0}
2^{−10} × 5256649930600386  5.13344719785 × 10^{12}  2^{−55.7177747539}  J_{0}
2^{−16} × 6535297120514194  99720720222.7  2^{−55.7166597393}  Y_{0}
2^{−35} × 6458928246558283  187979.55262  2^{−55.5787196547}  Y_{1}
2^{−26} × 5245948062070016  78170717.6875  2^{−55.509211239}  Y_{0}
2^{−20} × 7547179409128835  7197551163.8  2^{−55.4766161411}  Y_{1}
2^{−26} × 8441237061651159  125784234.131  2^{−55.4755750981}  Y_{0}
2^{−45} × 7046625970325583  200.277155793  2^{−55.4323204082}  J_{0}
2^{−40} × 8936924570334870  8128.08554686  2^{−55.3841525656}  J_{1}
2^{−29} × 8114890393276829  15115161.2276  2^{−55.1139949895}  J_{1}
2^{−53} × 8048625784723434  0.893576966279  2^{−55.0615557705}  Y_{0}

Table 3. Doubles in range [0, 2^{90}] within 2^{−55} of Bessel zero

a special 2-part Horner step, which takes an existing 2part value (h, l) and produces a new one (h′, l′) using 2part coefficients (c_{H}, c_{L}) and dividing by a 2-part variable x_{H} + x_{L}. (In our application we have x_{H} + x_{L} = x^{2} and the (c_{H}, c_{L}) are 2-part floating-point approximations to the α_{n} series coefficients.) That is, we get to good accuracy:

(h′ + l′) = (c_{H} + c_{L}) + \frac{1}{x_{H} + x_{L}} (h + l)

In a preamble, which only needs to be executed once for a given x_{H} +x_{L}, we compute an approximate inverse y and correction e:

y = 1/x_{H} e_{1} = 1 − x_{H} · y e_{2} = e_{1} − x_{L} · y e = y · e_{2}

Each Horner step then consists of the following computation. One can show that provided the terms decrease in size (as they do in our application), this maintains a relative error of about twice the working precision.

h′ = c_{H} + y · h t = c_{L} + e · h
r_{1} = h′ − c_{H}
r = y · h − r_{1}
s = t + r
l′ = s + y · l

The latency may seem like 5 fmas, but note that we get h′ from h in just one latency and only use l to obtain l′ on the last step. Thus we can in fact pipeline a number of successive applications separated only by one fma latency. If we have a 2-part result (r, c) from argument reduction and (h, l) is a 2-part α_{n}^{(1)} +α_{n}^{(3)}/x^{2} +α_{n}^{(5)}/x^{4} +⋯ then we can make the final combination using the same step with only a 1-part x needed (this is just our input number x, whereas earlier we were dividing by x^{2} at each stage):

(h′ + l′) = (r + c) + \frac{1}{x}(h + l)

and then perform a conventional floating-point addition h′ + l′ to obtain the final result r as a double-extended number. Since this is bounded by approximately |r| ≤ π/4 and the sine and cosine functions are well-conditioned in this area, we can just evaluate sin(r) or cos(r) straightforwardly.

Exact value of double Approximate decimal value Distance from zero Function
2^{796} × 6381956970095103 2.65968632416 × 10^{255}  2^{−61.8879179362}  Y_{0}
2^{796} × 6381956970095103 2.65968632416 × 10^{255}  2^{−61.8879179362}  J_{1}
2^{78} × 5916243447979695  1.78807486485 × 10^{39}  2^{−59.9300883695}  Y_{1}
2^{78} × 5916243447979695  1.78807486485 × 10^{39}  2^{−59.9300883695}  J_{0}
2^{938} × 8444920710073313 1.96214685729 × 10^{298}  2^{−59.7839697826}  Y_{0}
2^{938} × 8444920710073313 1.96214685729 × 10^{298}  2^{−59.7839697826}  J_{1}
2^{524} × 5850965514341686 3.21325554979 × 10^{173}  2^{−59.4812472194}  Y_{1}
2^{524} × 5850965514341686 3.21325554979 × 10^{173}  2^{−59.4812472194}  J_{0}
2^{807} × 8160885118204141 6.96536316842 × 10^{258}  2^{−59.1402789847}  Y_{1}
2^{807} × 8160885118204141 6.96536316842 × 10^{258}  2^{−59.1402789847}  J_{0}
2^{627} × 8360820580228475 4.65646330711 × 10^{204}  2^{−59.0913203893}  Y_{0}
2^{627} × 8360820580228475 4.65646330711 × 10^{204}  2^{−59.0913203893}  J_{1}
2^{144} × 8583082635084172  1.91409138863 × 10^{59}  2^{−59.0538090794}  Y_{1}
2^{144} × 8583082635084172  1.91409138863 × 10^{59}  2^{−59.0538090794}  J_{0}
2^{242} × 7958046405119485  5.6242603729 × 10^{88}  2^{−58.9935827902}  Y_{0}
2^{242} × 7958046405119485  5.6242603729 × 10^{88}  2^{−58.9935827902}  J_{1}
2^{561} × 6808218460873451 5.13879213026 × 10^{184}  2^{−58.9110422712}  Y_{0}
2^{561} × 6808218460873451 5.13879213026 × 10^{184}  2^{−58.9110422712}  J_{1}
2^{200} × 7636753411044619  1.22717895908 × 10^{76}  2^{−58.8981116632}  Y_{1}
2^{200} × 7636753411044619  1.22717895908 × 10^{76}  2^{−58.8981116632}  J_{0}
2^{193} × 6366906923947931  7.99314450027 × 10^{73}  2^{−58.5326981594}  Y_{1}
2^{193} × 6366906923947931  7.99314450027 × 10^{73}  2^{−58.5326981594}  J_{0}
2^{14} × 6617649673795284  1.08423572255 × 10^{20}  2^{−58.438199858}  Y_{1}
2^{37} × 6311013172270677  8.67379045745 × 10^{26}  2^{−58.4361612221}  Y_{1}
2^{37} × 6311013172270677  8.67379045745 × 10^{26}  2^{−58.4361612218}  J_{0}
2^{14} × 6617649673795284  1.08423572255 × 10^{20}  2^{−58.4356039989}  J_{0}
2^{886} × 5648695676206402 2.91423343144 × 10^{282}  2^{−58.152019003}  Y_{0}
2^{886} × 5648695676206402 2.91423343144 × 10^{282}  2^{−58.152019003}  J_{1}
2^{968} × 6221301883130153 1.55209063452 × 10^{307}  2^{−58.1018949515}  Y_{1}
2^{968} × 6221301883130153 1.55209063452 × 10^{307}  2^{−58.1018949515}  J_{0}
2^{502} × 6951690029616219 9.10223874078 × 10^{166}  2^{−58.0556510652}  Y_{0}
2^{502} × 6951690029616219 9.10223874078 × 10^{166}  2^{−58.0556510652}  J_{1}
2^{525} × 8776448271512529 9.63976664937 × 10^{173}  2^{−57.8962847187}  Y_{0}
2^{525} × 8776448271512529 9.63976664937 × 10^{173}  2^{−57.8962847187}  J_{1}
2^{90} × 6101578227064009  7.55338799011 × 10^{42}  2^{−57.8318447312}  Y_{0}
2^{90} × 6101578227064009  7.55338799011 × 10^{42}  2^{−57.8318447312}  J_{1}
2^{636} × 8438258240500718 2.40619075865 × 10^{207}  2^{−57.7878018343}  Y_{0}
2^{636} × 8438258240500718 2.40619075865 × 10^{207}  2^{−57.7878018343}  J_{1}
2^{129} × 4595526034082901  3.12755295225 × 10^{54}  2^{−57.7272223193}  Y_{1}
2^{129} × 4595526034082901  3.12755295225 × 10^{54}  2^{−57.7272223193}  J_{0}
2^{434} × 7750232865711478 3.43821372626 × 10^{146}  2^{−57.7099830586}  Y_{1}
2^{434} × 7750232865711478 3.43821372626 × 10^{146}  2^{−57.7099830586}  J_{0}

Table 4. Doubles closest to Bessel zero

But in the interest of speed we use our own custom implementation of these trigonometric functions that bypasses the usual range reduction step, because this range reduction has already been performed, as described next.

One may be able to use an off-the-shelf argument reduction routine by reducing modulo π/4 and then modifying the result, provided that enough information about the quotient as well as the remainder is provided. Note that one can always adapt an argument reduction for 2^{a}π for 2^{b}π by multiplying the input and dividing the output(s) by 2^{a−b}. Instead, we programmed our own variant of the standard scheme [7], which optimizes for speed and simplicity at the cost of much greater storage. For an input x = 2^{a}m with m ∈ ℤ we observe that

(x/π) mod 1 = (2^{a}m/π) mod 1 = (m · 2^{a}/π) mod 1

We may apply the modulus twice within the calculation without changing the result, since m is an integer:

(2^{a}m/π) mod 1 = (m · [(2^{a}/π) mod 1]) mod 1

In IEEE double precision there are only about 2048 possibilities for a, so we can simply precompute and tabulate the values:

P_{a} = ((2^{a}/π) mod 1)/2^{a}

so we just need to calculate at runtime (P_{a} · x) mod 1. Each P_{a} is stored in three parts using floating-point numbers, so P_{a} = P_{a}^{hi} + P_{a}^{med} + P_{a}^{lo}. The final computation uses some straightforward fma-based techniques to ensure accuracy in the computations, also adding and subtracting a ‘shifter’ S = 2^{62} + 2^{61} to fix the binary point and force integer truncation. The final result x_{2} is an accurate reduced argument (x/π) mod 1, though it may be a little greater than 1/2 in magnitude, which can then be postmultiplied by π.

N_{S} = S + x · P_{a}^{hi}
N = S − N_{S}
x_{0} = N + x · P_{a}^{hi}
x_{1} = x_{0} + x · P_{a}^{med}
x_{2} = x_{1} + x · P_{a}^{lo}

The complete code (written in C) runs in about 160 cycles for all arguments, and with more aggressive scheduling and optimization, we believe it could be brought below 100 cycles.

7. Further work

Our α_{n} series were designed to achieve better than 2^{−120} absolute accuracy over the whole range. However, this very high degree of accuracy is only required in the neighborhood of the zeros. It would be interesting to investigate computing a minimax economization using a more refined weight function to take this into account; it may be that a significant further reduction in the degree is possible. At least we could explore more sophisticated methods of arriving at minimax approximations, taking into account the fact that the coefficients are machine numbers [2].

In general the above techniques should generalize straightforwardly to J_{n} and Y_{n} for fixed, moderate n, yielding algorithms that are fast and accurate at the cost of relatively large tables. The problem of implementing generic functions when n is another parameter is potentially more difficult, and has not yet been considered. Another more ambitious goal still would be to strive for perfect rounding, but for this we would need more comprehensive ‘worst case’ data for the functions in general, not merely for their zeros.

Our current code is just a simple prototype, and makes use of the somewhat uncommon combination, found on the Intel® Itanium® architecture, of the double-extended floatingpoint type and the fused multiply-add. It should not be too difficult to produce a portable C implementation that would run on any reasonable platform supporting IEEE floatingpoint, though the performance figures would probably not be quite as good as those reported here. This could be of fairly wide interest, since we are not aware of any other double-precision Bessel function implementations offering the combination of accuracy and speed we attain here.

Acknowledgements

The author is grateful to Peter Tang, both for suggesting the Bessel functions as an interesting line of research and for his explanation of the established theory. Thanks also to the anonymous reviewers for ARITH, whose comments were of great help in improving the paper.

References

[1] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, volume 55 of Applied Mathematics Series. US National Bureau of Standards, 1964.

[2] N. Brisebarre, J.-M. Muller, and A. Tisserand. Computing machine-efficient polynomial approximations. ACM Transactions on Mathematical Software, 32:236–256, 2006.

[3] J. F. Hart, E. W. Cheney, C. L. Lawson, H. J. Maehly, C. K. Mesztenyi, J. R. Rice, H. G. Thatcher, and C. Witzgall. Computer Approximations. Robert E. Krieger, 1978.

[4] V. Lefèvre and J.-M. Muller. Worst cases for correct rounding of the elementary functions in double precision. Research Report 4044, INRIA, 2000.

[5] M. D. McIlroy. Squinting at power series. Software — Practice and Experience, 20:661–683, 1990.

[6] J.-M. Muller. Elementary functions: Algorithms and Implementation. Birkhäuser, 2nd edition, 2006.

[7] M. Payne and R. Hanek. Radian reduction for trigonometric functions. SIGNUM Newsletter, 18(1):19–24, 1983.

[8] R. A. Smith. A continued fraction analysis of trigonometric argument reduction. IEEE Transactions on Computers, 44:1348– 1351, 1995.

[9] G. N. Watson. A treatise on the theory of Bessel functions. Cambridge University Press, 1922.