Rational Approximations for the Modified Bessel Function of the First Kind – I1(x) for Computations with Double Precision

by on November 12, 2015

In connection with our previous study, in this post we analyse existing rational approximations for modified Bessel function of the first kind of order one – .

At first, we consider well known and widely used methods, then propose new approximations which deliver the lowest relative error and have favorable computation structure (pure polynomials).

We study properties of the following rational approximations: Blair & Edwards[1], Numerical Recipes[2], SPECFUN 2.5[3], and Boost 1.59.0[4]. Similar to , all of them have the structure and intervals, proposed by Blair & Edwards[1]:

(1)

where and are polynomials of corresponding degrees.

Actual accuracy of each approximation in terms of machine epsilon, eps = 2-52 ˜ 2.22e-16, is plotted below:
(Please see our report on I0(x) for technical details).

Blair & Edwards[1]

Numerical Recipes, 3rd edition, 2007[2]

SPECFUN 2.5[3]

SPECFUN approximation gives arbitrary large errors near 1/x = 0.0033705641 (the plot is trimmed). As it turned out, polynomial in numerator has a root in this point, which means that approximation cannot be used at all. DON’T use the SPECFUN for computing the !

Boost 1.59.0[4]

GNU GSL 2.0[5]
(GNU GSL library relies on Chebyshev series)

NAG Toolbox for MATLAB[7]
(Most likely NAG library uses Chebyshev expansion)

NAG routine has serious accuracy degradation for . Accuracy plot for the range:

MATLAB 2015a[6]
(MATLAB uses Chebyshev series?)

As we see, in all cases, relative errors go way beyond the limits we naturally expect from double precision (red lines show corridor where error ideally should be).

Peak relative error in terms of machine epsilon –

Method/Software
Minimax Rational Approximations:
Blair & Edwards 4.62 2.55
Numerical Recipes, 2007 5.28 4.34
SPECFUN 2.5 4.62
Boost 1.59.0 4.39 2.55
Refined – (3), (4) 2.06 0.92
Chebyshev Expansions:
GNU GSL 2.0 2.12 3.26
NAG Mark 24 5.32 257.01
MATLAB 2015a 3.41 2.82

computed over 50K uniformly distributed random samples in each interval

Refined Approximations

The main issue with approximations above is that they were built only by minimizing theoretical relative error – computed in extended precision. The rounding and cancellation effects (e.g. conditioning of approximations) with respect to double precision were not included in the search process.

We built software which rigorously iterates through all possible rational approximations (1), and selects optimal by two criteria: (a) minimum relative error (with respect to double precision) and (b) high approximation order. The latter was assigned with lower priority though. The split point is also subject to search.

Similar to I0(x) we use extended approximation form for small arguments:

(2)

The has beautiful power series at 0, where the first few coefficients do not introduce any additional rounding errors (division by 2 do not alter mantissa, only the exponent). They are perfect candidates to be included in approximation form, as this allows us to free several degrees of freedom which can be spent to attain higher accuracy.

Found approximations have the following properties:

Relative error is up to 2.06eps on and 0.92eps on . Although we didn’t manage to push error lower the machine epsilon completely, but still we reached the highest accuracy level among all others mentioned above.

Final expressions for optimal approximations:

• For small arguments, :

(3)

(use Horner scheme to compute the outer polynomial as well)

• For large arguments, :

(4)

Coefficients rounded to 25 decimal digits:

n for
08.3333333333333333311567967e-02
16.9444444444444450632369337e-03
23.4722222222221933634809047e-04
31.1574074074079326676719210e-05
42.7557319223490964726712181e-07
54.9209498642000488498034902e-09
66.8346524852208360284643288e-11
77.5940608652484019265663823e-13
86.9036483611746228414130722e-15
95.2305536429160706017626195e-17
103.3486060280590464196327437e-19
111.8645262719811753663834319e-21
127.9611250107842314599760659e-24
135.3251032089995165438568695e-26
n for
0 3.9894228040143270388374079e-01
1-1.4960335515072058522575487e-01
2-4.6751048269476797374239762e-02
3-4.0907267094886972971863462e-02
4-5.7501487840859800117669379e-02
5-1.1428156617865937773864845e-01
6 6.7988447242260666801129937e-02
7-2.2694203870019250176636896e+01
8 9.7548286270114208672947525e+02
9-2.9286459257939415083570152e+04
10 4.9934855620495985742805154e+05
11 5.7682364160056137069002930e+05
12-3.1576840778898356890175020e+08
13 1.0484906321376589515223174e+10
14-2.0918193917759394367113655e+11
15 2.9320804098307168426392082e+12
16-3.0147278411132255281401004e+13
17 2.2950466603697814797615042e+14
18-1.2816007548999035598180100e+15
19 5.1086996139908353110844064e+15
20-1.3774917783425787550429723e+16
21 2.2531580094188348024267027e+16
22-1.6895178303473738478791245e+16

text

Conclusion

Obtained approximations for deliver the lowest relative error compared to all commonly used approximations in mainstream software libraries: Blair & Edwards[1], Numerical Recipes[2], SPECFUN 2.5[3], Boost 1.59.0[4] and even GNU GSL[5].

3. Cody W.J., Algorithm 715: SPECFUN – A Portable FORTRAN Package of Special Function Routines and Test Drivers, ACM Transactions on Mathematical Software, Volume 19, Number 1, March 1993, pages 22-32. SPECFUN: Modified Bessel function of the first kind, I1(x).