Second period resume
Here I present a brief resume of the work done in this second month.
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
lospecfunc.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
lospecfun.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.6257e16 
0.0000e+00 
0.0000e+00 
1.3379e16 
1.1905e16 
Cephes 
2.825060e08 
ERR 
ERR 
ERR 
ERR 
GSL 
4.8770e16 
2.2068e16 
4.2553e16 
1.3379e16 
1.1905e16 
C++ 
2.82506e08 
2.68591e07 
1.55655e05 
8.58396e07 
0.000389545 

1e15 
1e20 
1e25 
1e30 
Amos 
1.3522e16 
1.6256e16 
15.22810 
2.04092 
GSL 
1.3522e16 
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.95201138653707e21
sage: bessel_J(1,10^40).n(digits=16)
7.522035532561300e21
sage: bessel_J(1,10^40).n(digits=20)
5.7641641376196142462e21
sage: bessel_J(1,10^40).n(digits=25)
1.620655257788838811378880e21
sage: bessel_J(1,10^40).n(digits=30)
1.42325562284462027823158499095e21
sage: bessel_J(1,10^40).n(digits=35)
1.4232556228446202782315849909515310e21
sage: bessel_J(1,10^40).n(digits=40)
1.423255622844620278231584990951531026335e21
The value stabilizes only when we use more than 30 digits (two times the digits used in double precision).
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.
Betainc
During these last days I also started to implement betainc from scratch, I think it will be ready for the first days of August. Then, it will be necessary to rewrite also betaincinv, since the actual version doesn't have the "upper" version. This should not be too difficult. I think we can use a simple Newton method (as for gammaincinv), the only problem will be, as for gammaincinv, to find good initial guesses.