## Bessel functions

The topic of this month were Bessel functions. On the bug tracker it is reported only a problem regarding the "J" one (see bug 48316), but the same problem is present in every type of Bessel function (they return`NaN + NaNi`when the argument is too big).

### Amos, Cephes, C++ and GSL

Actually, Bessel functions in Octave are computed via the Amos library, written in Fortran. Studying the implementation I discovered that the reported bug follows from the fact that if the input is too large in module, the function`zbesj.f`set IERR to 4 (IERR is a variable which describe how the algorithm has terminate) and set the output to zero, then

`lo-specfunc.cc`return NaN when IERR=4. Obviusly, the same happen for other Bessel functions.

What I initially did was to "unlock" these .f files in such a way to give still IERR=4, but computing anyway the output, and modify

`lo-specfun.cc`in order to return the value even if IERR is 4. Then I tested the accuracy, together with other libraries.

On the bug report where suggested the Cephes library, so I tested also them, the C++ Special mathematical functions library, and the GSL (Gnu Scientific library). Unfortunately, both these alternatives work worse than Amos. I also tried to study and implement some asymptotic expansions by myself to use in the cases which give inaccurate results, unfortunately without success.

For completeness, in the following table there are some results of the tests (ERR in Cephes refers to

`cos total loss of precision error`):

1e09 | 1e10 | 1e11 | 1e12 | 1e13 | |
---|---|---|---|---|---|

Amos | 1.6257e-16 | 0.0000e+00 | 0.0000e+00 | 1.3379e-16 | 1.1905e-16 |

Cephes | 2.825060e-08 | ERR | ERR | ERR | ERR |

GSL | 4.8770e-16 | 2.2068e-16 | 4.2553e-16 | 1.3379e-16 | 1.1905e-16 |

C++ | 2.82506e-08 | 2.68591e-07 | 1.55655e-05 | 8.58396e-07 | 0.000389545 |

1e15 | 1e20 | 1e25 | 1e30 | |
---|---|---|---|---|

Amos | 1.3522e-16 | 1.6256e-16 | 15.22810 | 2.04092 |

GSL | 1.3522e-16 | 0 | 15.22810 | 2.04092 |

### The problem with the double precision

As I explained in my last post the problem seems to be that there are no efficient algorithms in double precision for arguments so large. In fact, the error done by Amos is quite small if compared with the value computed with SageMath (which I'm using as reference one), but only if we use the double precision also in Sage: using more digits, one can see that even the first digit change. Here the tests:`sage: bessel_J(1,10^40).n()`

7.95201138653707e-21

sage: bessel_J(1,10^40).n(digits=16)

-7.522035532561300e-21

sage: bessel_J(1,10^40).n(digits=20)

-5.7641641376196142462e-21

sage: bessel_J(1,10^40).n(digits=25)

1.620655257788838811378880e-21

sage: bessel_J(1,10^40).n(digits=30)

1.42325562284462027823158499095e-21

sage: bessel_J(1,10^40).n(digits=35)

1.4232556228446202782315849909515310e-21

sage: bessel_J(1,10^40).n(digits=40)

1.423255622844620278231584990951531026335e-21

The value stabilizes only when we use more than 30 digits (two times the digits used in double precision).

7.95201138653707e-21

sage: bessel_J(1,10^40).n(digits=16)

-7.522035532561300e-21

sage: bessel_J(1,10^40).n(digits=20)

-5.7641641376196142462e-21

sage: bessel_J(1,10^40).n(digits=25)

1.620655257788838811378880e-21

sage: bessel_J(1,10^40).n(digits=30)

1.42325562284462027823158499095e-21

sage: bessel_J(1,10^40).n(digits=35)

1.4232556228446202782315849909515310e-21

sage: bessel_J(1,10^40).n(digits=40)

1.423255622844620278231584990951531026335e-21

### The decision

Even if we are aware of the fact that the result is not always accurate, for MATLAB compatibility we decided to unlock the Amos, since they are still the most accurate and, even more important, the type of the inputs is the same as in MATLAB (while, for example, GSL doesn't accept complex $x$ value). Moreover, in Octave is possible to obtain in output the value of IERR, thing that is not possible in MATLAB.You can find the work on the bookmark "bessel" of my repository.