## lunedì 24 luglio 2017

### Second period resume

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 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).

### 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.

## Bessel functions

During the second part of GSoC I have to work on Bessel functions. The bug is related in particular on the Bessel function of first kind and regard the fact that the function is not computed if the $x$ argument is too big.
During this week I studied the actual implementation i.e. the Amos library, written in Fortran. The problem is related to the fact that zbesj.f simply refuse to compute the result when the input argument is too large (the same problem happens for all bessel functions, both for real and complex arguments). What I did was to "unlock" the amos in such a way to compute the value in any case (which seems to be the strategy used by MATLAB). Then I compared the relative errors with the relative errors done by the Gnu Scientific Library (GSL).
At first, I would precise some limitation present in GSL:
• The parameter $x$ must be real, while it can be complex in Amos and MATLAB. The same holds for the parameter $\alpha$.
• The class of the output does not adapt to the class of the input, returning always a double.
Doing some tests, it seems that the amos work better in terms of accuracy (even if we are talking about errors which do not differ in the order). I concentrated on values in which Amos usually refuse to compute, since in every other zone of the $x-\alpha$ plane it is known the error is in the order of 1e-15. For $\alpha\in\{-1,0,1\}$ they return the same result, while for other values of $\alpha$, Amos are in general more accurate.
Anyway, I would remark the fact that there are not, as far as I know, accurate algorithms in double precision for very large magnitude. In fact, for such arguments, both Amos and GSL make a relative error of the order of the unity. This problem is evident when using SageMath to compute accurate values, e.g.
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 values "stabilize" only when the number of digits is bigger than 30, far away from double precision.

## Incomplete gamma function

I also gave a look on how to use GSL to eventually improve the incomplete gamma function. It is not possible, however, to use only GSL functions gsl_sf_gamma_inc_P and gsl_sf_gamma_inc_Q due to their limitations:
• There is no the "scaled" option.
• The value of $x$ must be real. This may be not a problem since even MATLAB does not accept complex value (while the version I worked on does).
• The parameter $x$ must be positive. This is actually a problem of MATLAB compatibility.
• The class of the output does not adapt to the class of the input, returning always a double.
I tested the accuracy of the GSL and they work better when $a<1$ and $x\ll1$ so I think I will fix gammainc.m using the algorithm of gsl_sf_gamma_inc_P and gsl_sf_gamma_inc_Q for that values of the input arguments.

## Incomplete Beta function

The actual incomplete beta function needs to be replaced. I've already studied the continued fraction exapnsion (which seems to be the best one to use). In GSL is implemented in a good way but still present two limitations:
• Does not adapt to the input class (it always work in double precision).
• There is no "upper" version. This is missing also in the current betainc function, but it is present in MATLAB. Unfortunately, it is not sufficient to compute it as $1 - I_x(a,b)$, it is necessary to find an accurate way to compute directly it.
So I will take inspiration to GSL version of the function but I think I will write beatinc as a single .m file.

## sabato 1 luglio 2017

### First period resume

First period resume

# First period resume

Here I present a brief resume of the work done in this first month.

## 1st week: Gammainc

The incomplete gamma function presented a series of inaccurate results in the previous implementation, so it has been decided to rewrite it from scratch. A big part of the work was done by Marco and Nir (here the discussion). I gave my contribution by fixing some problems in the implementation and by making a functioning commit (thanks to the suggestions given by Carne during the OctConf in Geneve).

## 2nd and 3rd week: Gammaincinv

The inverse of the incomplete gamma function was completely missing in Octave (and this was a problem of compatibility with Matlab). Now the function is present, written as a single .m file. The implementation consist in a simple, but efficient, Newton's method.

## Last week: Betainc

When I first submitted my timeline, there was only a bug related to betainc, involving the input validation. During the first period of GSoC emerged a new bug of "incorrect result" type. From this, me and my mentor decided that it is necessary to rewrite the function, so I used this last week to study efficient algorithms to evaluate it. In particular, seems that the most efficient way is to use continude fractions, maybe after a rescale of the function. I'm already working on it, but I will complete it during the first week of august, after finishing the work on Bessel functions.