Final Resume

During the GSoC I worked on different special functions that needed to be improved or implemented from scratch. Discussing with my mentors and the community, we decided that my work should be pushed on a copy of the scource code of Octave on my repository [1] and then I should have work with different bookmarks for each function I had to work on. When different functions happened to be related (e.g. gammainc and gammaincinv), I worked on these on the same bookmark. I present now a summary and the bookmarks related to the functions.

Incomplete gamma function

bookmark: gammainc
first commit: d1e03faf080b
last commit: 107dc1d24c1b
added files: /libinterp/corefcn/__gammainc_lentz__.cc, /scripts/specfun/gammainc.m, /scripts/specfun/gammaincinv.m
removed files:/libinterp/corefcn/gammainc.cc, /liboctave/external/slatec-fn/dgami.f, /liboctave/external/slatec-fn/dgamit.f, /liboctave/external/slatec-fn/gami.f, /liboctave/external/slatec-fn/gamit.f, /liboctave/external/slatec-fn/xdgami.f, /liboctave/external/slatec-fn/xdgamit.f, /liboctave/external/slatec-fn/xgmainc.f, /liboctave/external/slatec-fn/xsgmainc.f
modified files: NEWS, /doc/interpreter/arith.txi, /libinterp/corefcn/module.mk, /liboctave/external/slatec-fn/module.mk, /liboctave/numeric/lo-specfun.cc, /scripts/specfun/module.mk

Summary of the work

On this bookmark I worked on the incomplete gamma function and its inverse.
The incomplete gamma function gammainc had both missing features (it were missed the "scaled" options) and some problem of inaccurate result type (see bug # 47800). Part of the work was already been done by Marco and Nir, I had to finish it. We decided to implement it as a single .m file (gammainc.m) which call (for some inputs) a subfunction written in C++ (__gammainc_lentz__.cc).
The inverse of the incomplete gamma function was missing in Octave (see bug # 48036). I implemented it as a single .m file (gammaincinv.m) which uses a Newton method.

Bessel functions

bookmark: bessel
first commit: aef0656026cc
last commit: e9468092daf9
modified files: /liboctave/external/amos/README, /liboctave/external/amos/cbesh.f, /liboctave/external/amos/cbesi.f, /liboctave/external/amos/cbesj.f, /liboctave/external/amos/cbesk.f, /liboctave/external/amos/zbesh.f, /liboctave/external/amos/zbesi.f, /liboctave/external/amos/zbesj.f, /liboctave/external/amos/zbesk.f, /liboctave/numeric/lo-specfun.cc, /scripts/specfun/bessel.m

Summary of the work

On this bookmark I worked on Bessel functions.
There was a bug reporting NaN as output when the argument $x$ was too large in magnitude (see bug # 48316). The problem was given by Amos library, which refuses to compute the output in such cases. I started "unlocking" this library, in such a way to compute the output even when the argument was over the limit setted by the library. Then I compared the results with other libraries (e.g. Cephes [2], Gnu Scientific library [3] and C++ special function library [4]) and some implementations I made. In the end, I discovered that the "unlocked" Amos were the best one to use, so we decided to maintain them (in the "unlocked" form), modifying the error variable to explain the loss of accuracy.

Incomplete beta function

bookmark: betainc
first commit: 712a069d2860
last commit: e0c0dd40f096
added files: /libinterp/corefcn/__betainc_lentz__.cc, /scripts/specfun/betainc.m, /scripts/specfun/betaincinv.m
removed files: /libinterp/corefcn/betainc.cc, /liboctave/external/slatec-fn/betai.f, /liboctave/external/slatec-fn/dbetai.f, /liboctave/external/slatec-fn/xbetai.f, /liboctave/external/slatec-fn/xdbetai.f
modified files: /libinterp/corefcn/module.mk, /liboctave/external/slatec-fn/module.mk, /liboctave/numeric/lo-specfun.cc, /liboctave/numeric/lo-specfun.h, /scripts/specfun/module.mk, /scripts/statistics/distributions/betainv.m, /scripts/statistics/distributions/binocdf.m

Summary of the work

On this bookmark I worked on the incomplete beta function and its inverse.
The incomplete beta function missed the "upper" version and had reported bugs on input validation (see bug # 34405) and inaccurate result (see bug # 51157). We decided to rewrite it from scratch. It is now implemented ad a single .m file (betainc.m) which make the input validation part, then the output is computed using a continued fraction evaluation, done by a C++ function (__betainc_lentz__.cc).
The inverse was present in Octave but missed the "upper" version (since it was missing also in betainc itself). The function is now written as a single .m file (betaincinv.m) which implement a Newton method where the initial guess is computed by few steps of bisection method.

Integral functions

bookmark: expint
first commit: 61d533c7d2d8
last commit: d5222cffb1a5
added files:/libinterp/corefcn/__expint_lentz__.cc, /scripts/specfun/cosint.m, /scripts/specfun/sinint.m
modified files: /doc/interpreter/arith.txi, /libinterp/corefcn/module.mk, /scripts/specfun/expint.m, /scripts/specfun/module.mk

Summary of the work

On this bookmark I worked on exponential integral, sine integral and cosine integral. I already rewrote the exponential integral before the GSoC. Here I just moved the Lentz algorithm to an external C++ function (__expint_lentz__.cc), accordingly to gammainc and betainc. I've also modified the exit criterion for the asymptotic expansion using [5] (pages 1 -- 4) as reference.
The functions sinint and cosint were present only in the symbolic package of Octave but was missing a numerical implementation in the core. I wrote them as .m files (sinint.m and cosint.m). Both codes use the series expansion near the origin and relations with expint for the other values.

To do

There is still room for improvement for some of the functions I wrote. In particular, gammainc can be improved in accuracy for certain couple of values, and I would like to make a template version for the various Lentz algorithms in C++ so to avoid code duplication in the functions.
In October I will start a PhD in Computer Science, still here in Verona. This will permit me to remain in contact with my mentor Marco Caliari, so that we will work on these aspects.

[1] https://bitbucket.org/M_Ginesi/octave
[2] http://www.netlib.org/cephes/
[3] https://www.gnu.org/software/gsl/
[4] http://en.cppreference.com/w/cpp/numeric/special_math
[5] N. Bleistein and R.A. Handelsman, "Asymptotic Expansions of Integrals", Dover Publications, 1986.

Integral functions

During the last week I made few modifications to expint.m and wrote sinint.m and cosint.m from scratch. All the work done can be found on the bookmark expint of my repository.

expint

As I mentioned here I rewrote expint.m from scratch before the GSoC. During the last week I moved the Lentz algorithm to a .cc function (in order to remain coherent with the implementations of gammainc and betainc) and added few tests.

sinint

The sinint function is present in the symbolic package, but is not present a numerical implementation in the core.
The sine integral is defined as $$ \text{Si} (z) = \int_0^z \frac{\sin(t)}{t}\,dt. $$ To compute it we use the series expansion $$ \text{Si}(z) = \sum_{n=0}^\infty \frac{(-1)^n z^{2n+1}}{(2n+1)(2n+1)!} $$ when the module of the argument is smaller than 2. For bigger values we use the following relation with the exponential integral $$ \text{Si} = \frac{1}{2i} (E_1(iz)-E_1(-iz)) + \frac{\pi}{2},\quad |\text{arg}(z)| < \frac{\pi}{2}$$ and the following simmetry relations $$ \text{Si}(-z) = -\text{Si}(z), $$ $$ \text{Si}(\bar{z}) = \overline {\text{Si}(z)}. $$ The function is write as a single .m file.

cosint

As the sinint function, also cosint is present in the symbolic package, but there is not a numerical implementation in the core.
The cosine integral is defined as $$ \text{Ci} (z) = -\int_z^\infty \frac{\cos(t)}{t}\,dt. $$ An equivalent definition is $$ \text{Ci} (z) = \gamma + \log z + \int_0^z \frac{\cos t - 1}{t}\,dt. $$ To compute it we use the series expansion $$ \text{Ci}(z) = \gamma + \log z + \sum_{n=1}^\infty \frac{(-1)^n z^{2n}}{(2n)(2n)!} $$ when the module of the argument is smaller than 2. For bigger values we use the following relation with the exponential integral $$ \text{Ci} = -\frac{1}{2} (E_1(iz)+E_1(-iz)),\quad |\text{arg}(z)| < \frac{\pi}{2}$$ and the following simmetry relations $$ \text{Ci}(-z) = \text{Ci}(z) -i\pi,\quad 0<\text{arg}(z)<\pi, $$ $$ \text{Ci}(\bar{z}) = \overline{\text{Ci}(z)} .$$ As for sinint, also cosint is written as a single .m file.

betaincinv

The inverse of the incomplete beta function was present in Octave, but without the "upper" option (since it was missing in betainc itself). We decided to rewrite it from scratch using Newton method, as for gammaincinv (see my post on it if you are interested).
To make the code numerically more accurate, we decide which version ("lower" or "upper") invert depending on the inputs.
At first we compute the trivial values (0 and 1). Then the remaining terms are divided in two sets: those that will be inverted with the "lower" version, and those that will be inverted with the "upper" one. For both cases, we perform 10 iterations of bisection method and then we perform a Newton method.
The implementation (together with the new implementation of betainc) can be found on my repository, bookmark "betainc".

betainc

The betainc function has two bugs reported: #34405 on the input validation and #51157 on inaccurate result. Moreover, it is missing the "upper" version, which is present in MATLAB.

The function

The incomplete beta function ratio is defined as $$I_x(a,b) = \dfrac{B_x(a,b)}{B(a,b)},\quad 0\le x \le 1,\,a>0,\,b>0,$$ where $B(a,b)$ is the classical beta function and $$B_x(a,b)=\int_0^x t^{a-1}(1-t)^{b-1}\,dt.$$ In the "upper" version the integral goes from $x$ to $1$. To compute this we will use the fact that $$\begin{array}{rcl} I_x(a,b) + I_x^U(a,b) &=& \dfrac{1}{B(a,b)}\left( \int_0^x t^{a-1}(1-t)^{b-1}\,dt + \int_x^1 t^{a-1}(1-t)^{b-1}\,dt\right)\\ &=&\dfrac{1}{B(a,b)}\int_0^1 t^{a-1}(1-t)^{b-1}\,dt\\ &=&\dfrac{B(a,b)}{B(a,b)}\\ &=&1 \end{array}$$ and the relation $$I_x(a,b) + I_{1-x}(b,a) = 1$$ so that $$I_x^U(a,b) = I_{1-x}(b,a).$$

The implementation

Even if it is possible to obtain a Taylor series representation of the incomplete beta function, it seems to not be used. Indeed the MATLAB help cite only the continuous fraction representation present in "Handbook of Mathematical Functions" by Abramowitz and Stegun: $$I_x(a,b) = \dfrac{x^a(1-x)^b}{aB(a,b)}\left(\dfrac{1}{1+} \dfrac{d_1}{1+} \dfrac{d_2}{1+}\ldots\right)$$ with $$d_{2m+1} = -\dfrac{(a+m)(a+b+m)}{(a+2m)(a+2m+1)}x$$ and $$d_{2m} = \dfrac{m(b-m)}{(a+2m-1)(a+2m)}x$$ which seems to be the same strategy used by GSL. To be more precise, this continued fraction is computed directly when $$x<\dfrac{a-1}{a+b-2}$$ otherwise, the computed fraction is used to compute $I_{1-x}(b,a)$ and then it is used the fact that $$I_x(a,b) = 1-I_{1-x}(b,a).$$ In my implementation I use a continued fraction present in "Handboob of Continued Fractions for Special Functions" by Cuyt, Petersen, Verdonk, Waadeland and Jones, which is more complicated but converges in fewer steps: $$\dfrac{B(a,b)I_x(a,b)}{x^a(1-x)^b} = \mathop{\huge{\text{K}}}_{m=1}^\infty \left(\dfrac{\alpha_m(x)}{\beta_m(x)}\right),$$ where $$\begin{array}{rcl} \alpha_1(x) &=&1,\\ \alpha_{m+1}(x) &=&\dfrac{(a+m-1)(a+b+m-1)(b-m)m}{(a+2m-1)^2}x^2,\quad m\geq 1,\\ \beta_{m+1}(x) &=&a + 2m + \left( \dfrac{m(b-m)}{a+2m-1} - \dfrac{(a+m)(a+b+m)}{a+2m+1} \right)x,\quad m\geq 0. \end{array}$$ This is most useful when $$x\leq\dfrac{a}{a+b},$$ thus, the continued fraction is computed directly when this condition is satisfied, while it is used to evaluate $I_{1-x}(b,a)$ otherwise.
The function is now written as a .m file, which check the validity of the inputs and divide the same in the values which need to be rescaled and in those wo doesn't need it. Then the continued fraction is computed by an external .c function. Finally, the .m file explicit $I_x(a,b)$.

betaincinv

Next step will be to write the inverse. It was already present in Octave, but is missing the upper version, so it has to be rewritten.

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.

Gnu Scientific Library

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.

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.