lunedì 28 agosto 2017

Final Resume

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

sabato 19 agosto 2017

Integral functions

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.

sabato 12 agosto 2017

betaincinv

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

mercoledì 2 agosto 2017

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.

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.

martedì 11 luglio 2017

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.

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.

martedì 20 giugno 2017

Timetable: modification

Timetable: modification

Timetable: modification

According to my timetable (that you can find here), during this last week of June, I should've work on the input validation of betainc. Since a new bug related to this function has been found and, moreover, the actual implementation doesn't accept the "lower" or "upper" tail (as MATLAB do), me and my mentor decided to use this week to start studying how to rewrite betainc (main references will be [1] and [2]) and to use the last part fo the GSoC to actually implement it. In this way, my timetable remain almost identical (I will use July to work on Bessel functions) and I will be able to fix also this problem.

[1] Abramowitz, Stegun "Handbook of Mathematical Functions"
[2] Cuyt, Brevik Petersen, Vendonk, Waadeland "Handbook of Continued Fractions for Special Functions"

venerdì 16 giugno 2017

gammaincinv

The gammaincinv function

The inverse of the incomplete gamma function

Definition

Given two real values $y\in[0,1]$ and $a>0$, gammaincinv(y,a,tail) gives the value $x$ such that \[P(x,a) = y,\quad\text{ or }\quad Q(x,a) = y\] depending on the value of tail ("lower" or "upper").

Computation

The computation is carried out via standard Newton method, solving, for $x$ \[ y - P(x,a) = 0 \quad\text{ or }\quad y - Q(x,a) = 0. \] The Newton method uses as exit criterion a control on the residual and on the maximum number of iterations.
For numerical stability, the initial guess and which gammainc function is inverted depends on the values of $y$ and $a$. We use the following property: \[\text{gammaincinv}(y,a,\text{`upper'})=\text{gammaincinv}(1-y,a,\text{`lower'})\] which follows immediately by the property \[P(x,a) + Q(x,a) = 1.\] We start computing the trivial values, i.e. $a=1,y=0$ and $y=1$. Then we rename the variables as follows (in order to obtain a more classical notation): if tail = "lower", then $p=y,q=1-y$; while if tail = "upper", $q=y,p=1-y$. We will invert with $p$ as $y$ in the cases in which we invert the lower incomplete gamma function, and with $q$ as $y$ in the cases in which we invert the upper one.
Now, the initial guesses (and the relative version of the incomplete gamma function to invert) are the following:
  1. $p<\dfrac{(0.2(1+a))^a}{\Gamma(1+a)}$. We define \[r = (p\,\Gamma(1+a))^{\frac{1}{a}}\] and the parameters \begin{align*} c_2 &= \dfrac{1}{a+1}\\ c_3 &= \dfrac{3a+5}{2(a+1)^2(a+2)}\\ c_4 &= \dfrac{8a^2+33a+31}{3(a+1)^3(a+2)(a+3)}\\ c_5 &= \dfrac{125a^4+1179a^3+3971a^2+5661a+2888}{24(1+a)^4(a+2)^2(a+3)(a+4)}\\ \end{align*} Then the initial guess is \[x_0 = r + c_2\,r^2+c_3\,r^3+c_4\,r^4+c_5\,r^5\] and we invert the lower incomplete gamma function.
  2. $q<\dfrac{e^{-a/2}}{\Gamma(a+1)}, a\in(0,10)$. The initial guess is \[x_0 = - \log (q) - \log(\Gamma(a))\] and we invert the upper incomplete gamma function.
  3. $a\in(0,1)$. The initial guess is \[x_0 = (p\Gamma(a+1))^{\frac{1}{a}}\] and we invert the lower incomplete gamma function.
  4. Other cases. The initial guess is \[x_0 = a\,t^3\] where \[t = 1-d-\sqrt{d}\,\text{norminv}(q),\quad d=\dfrac{1}{9a}\] (where norminv is the inverse of the normal distribution function) and we invert the upper incomplete gamma function.
The first three initial guesses can be found in the article "Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function rations" by Gil, Segura ad Temme; while the last one can be found in the documentation of the function igami.c of the cephes library.
You can find my work on my repository on bitbucket here, bookmark "gammainc".

Comments

Working on this function I found some problems in gammainc. To solve them, I just changed the flag for the condition 5 and added separate functions to treat Inf cases.

martedì 6 giugno 2017

gammainc

The gammainc function

The incomplete gamma function

Definition

There are actually four types of incomplete gamma functions:
  • lower (which is the default one) $$P(a,x) = \frac{1}{\Gamma(a)}\int_{0}^{x}e^{-t}t^{a-1}$$
  • upper $$Q(a,x) = \frac{1}{\Gamma(a)}\int_{x}^{+\infty}t^{a-1}e^{-t}dt$$
  • scaled lower and scaled upper which return the corresponding incomplete gamma function scaled by $$\Gamma(a+1)\frac{e^x}{x^a}$$

Computation

To compute the incomplete gamma function, there is a total of seven strategies depending on the position of the arguments in the $x-a$ plane.
  1. $x=0,a=0$. In this case, the output is given as follows:
    • for "lower" and "scaledlower" is $1$
    • for "upper" and "scaledupper" is $0$
  2. $x=0,a\neq0$. In this case:
    • for "lower" is $0$
    • for "upper" and "scalelower" is $1$
    • for "scaledupper" is $+\infty$
  3. $x\neq0,a=0$. In this case:
    • for "lower" is $1$
    • for "scalelower" is $e^x$
    • for "upper" and "scaledupper" is $0$
  4. $x\neq0,a=1$. In this case:
    • for "lower" is $1-e^x$
    • for "upper" is $e^{-x}$
    • for "scaledlower" is $\frac{e^x-1}{x}$
    • for "scaledupper" is $\frac{1}{x}$
  5. $0< x\leq36,a\in\{2,3,4,\ldots,17,18\}$. In this case we used the following expansion (which holds true only for $n\in\mathbb{N}$): denote $\gamma(a,x):=\Gamma(a)P(a,x)$, then $$\gamma(n,x) = (n-1)!\left( 1-e^{-x}\sum_{k=0}^{n-1} \frac{x^k}{k!}\right).$$
  6. $x+0.25< a $ or $x<0$, or $x$ not real or $|x|<1$ or $a<5$. Here we used the expansion $$\gamma(a,x) = e^{-x} x^a \sum_{k=0}^{+\infty} \frac{\Gamma(a)}{\Gamma(a+1+k)}x^k.$$
  7. for all the remaining cases was used a continued fraction expansion. Call $\Gamma(a,x) = \Gamma(a) Q(a,x)$, the expansion is $$\Gamma(a,x) = e^{-x}x^a\left( \frac{1}{x+1-a-}\frac{1(1-a)}{x+3-a-}\frac{2(2-a)}{x+5-a-}\cdots \right).$$ This case is handled by an external .cc function.

My contribution

Actually I dind't participate from the start in the (re)making of this function: most of the work done on it is due to Nir and Marco (see the discussion on the bug tracker).
My contribution was, before the GSoC, to adapt the codes in such a way to make the type of the output (and the tolerances inside the algorithms) coherent on the type of the input (single, int or double). Then I corrected few small bugs and made a patch helped by Carne during the OctConf in Geneve.
During this first week of GSoC I fixed the input validation and added tests to this function, finding some problem in the implementation. Most of them depended on the fact that the continued fractions were used in cases in which they are not optimal. To solve these problems I changed the conditions to use the series expansion instead of the continued fractions.
Now gammainc works properly also with complex argument (for the $x$ variable), while Matlab doesn't accept non-real input.
You can find my work on my repository on bitbucket here, bookmark "gammainc". I will work on this same bookmark the next two weeks while I will work on gammaincinv.

martedì 30 maggio 2017

Timetable

Timetable

Timetable

During this first period I searched for the bugs related to special functions in the bug tracker. I found these four bugs that need to be fixed:
  1. gammainc: the upper version rounds down to zero if the output is small. For this bug there is a lot of work already done by Marco and Nir. I will just fix the last few things (like the input validation).
  2. gammaincinv: after finishing gammainc, I should start working on the inverse (which is missing in Octave).
  3. betainc: it is necessary to fix the input validation.
  4. besselj: it gives NaN if the input is too large. Actually, also other versions of the Bessel function have the same problem.
My idea of the timetable is the following
  • First part (May 30 - June 26: 4 weeks):
    • 1st week: finish to fix gammainc (input validation, documentation, add tests)
    • 2nd and 3rd weeks: write gammaincinv (via Newton's method)
    • 4th week: fix the input validation for betainc (and, if needed, add tests and fix the documentation)
  • Second part (July 1 - July 24: 3 weeks): fix the Bessel functions. Some ideas comprehend to test other libraries in addition to the amos (which is the library currently in use) and try to implement new algorithms. Add tests and revise the documentation.
  • Final part (July 29 - August 21: 3 weeks): add tests to the other special functions to make sure they work properly, trying to fix the problems that, eventually, will be found and revise the documentation if needed.

sabato 13 maggio 2017

expint

The expint function

The exponential integral

This is the first function I worked on in Octave: you can find the discussion here. Even if it will not be part of my GSoC project, I think it would be interesting to show how the function has been rewritten.

Definition and main properties

The canonical definition for the exponential integral is $$E_i (z) = -\int_{-z}^{+\infty} \frac{e^{-t}}{t}dt$$ but, for Matlab compatibility, expint compute the function $$E_1(z) =\int_{z}^{+\infty}\frac{e^{-t}}{t}dt.$$ The two definitions are related, for positive real values $z$, by $$E_1(-z) = -E_i(z) - \mathbb{i}\pi.$$ More in general, the function $E_n$ is defined as $$E_n(z) = \int_1^{+\infty}\frac{e^{-zt}}{t^n}$$ and the following recurrence relation holds true: $$E_{n+1}(z) = \frac{1}{n}\left(e^{-z}-zE_n(z)\right).$$ Moreover $$E_n(\bar{z}) = \overline{E_n(z)}.$$

Computation

To compute $E_1$, I implemented three different strategies:
  • Taylor series (Abramowitz, Stegun, "Handbook of Mathematical Functions", formula 5.1.11, p 229): $$E_1(z) = -\gamma -\log z -\sum_{n=1}^{\infty}\frac{(-1)^nz^n}{nn!}$$ where $\gamma\approx 0.57721$ is the Euler's constant.
  • Continued fractions (Abramowitz, Stegun, "Handbook of Mathematical Functions", formula 5.1.22, p 229): $$E_1(z) = e^{-z}\left(\frac{1}{z+}\frac{1}{1+}\frac{1}{z+}\frac{2}{1+}\frac{2}{z+}\cdots\right)$$ or, in a more explicit notation $$E_1(z)=e^{-z}\frac{1}{z+\frac{1}{1+\frac{1}{z+\frac{2}{1+\frac{2}{z+\ddots}}}}}$$ This formulation has been implemented using the modified Lentz algorithm ("Numerical recipes in Fortran 77" p.165).
  • Asymptotic series ("Selected Asymptotic Methods with Application to Magnetics and Antennas" formula A.10 p 161): $$E_1(z)\approx \frac{e^{-z}}{z}\sum_{n=0}^{N}\frac{(-1)^n n!}{z^n}.$$ A difficulty about this approximation is that the series is divergent, and that is the reason why the sum goes up to a certain $N$ "big" and not up to infinity.
Then I tested them in the square $[-100,100]\times[-100,100]$ of the complex plane, comparing the result with the ones given by the Octave symbolic package. With these informations, I identified three zones of the complex plane: in each zone one strategy is better than the others.
Then the final function divides the input into three parts, accordingly to the position of the same in the complex plane, and compute the function using the opportune strategy.

martedì 9 maggio 2017

Introduction

Introduction to the project

Introduction

About me

My name is Michele Ginesi, I obtained a Bachelor's degree in Applied Mathematics in Verona, Italy. Now I am getting a Master degree in Mathematics in the same university.
I was selected to partecipate to the Google Summer of Code under GNU Octave for the project Make specfuns special again.

About the project

Special functions are an interesting and important topic in Mathematics, and so it is fundamental to have a way to compute them in an accurate way.
Some examples of special functions (whith some important application of the same) are:
  • Gamma function $\Gamma$. This is one of the most important: it can be viewed as the extension of the factorial function ($\Gamma(n)=(n-1)!$, for $n\in\mathbb{N}$) and it is a component in various distribution functions in probability theory.
  • Beta function $B$. This was the first known Scattering amplitude in String theory.
  • Bessel functions. These are the canonical solutions $u(r)$ of the Bessel's differential equation $$r^2\frac{d^2u}{dr^2}+r\frac{du}{dr}+(r^2-\alpha^2)u=0 $$ for any complex number $\alpha$. At a first approach may seems that they are an end in themselves, but actually the Bessel equation describes the radial component of the two dimensional wave equation ($u_{tt}-c^2\Delta u = 0$) applied on a disc.
The most common strategies used to approximate special functions are Taylor series and continued fractions, sometimes asymptotic series. Some particular functions can be implemented via recurrence formula or other type relations with other functions (for example, $B(z,w) = \Gamma(z)\Gamma(w)/\Gamma(z+w)$).
This project will be divided into three main parts:
  1. Fix the already known bugs related to special functions (e.g. #48316, #47800, #48036).
  2. When the known bugs will be fixed, I will proceed to add new tests to make sure that all the functions are accurate.
  3. Fix the new problems/bugs that, eventually, will be found during the second phase.
The main reference for this work will be Handbook of Mathematical Functions by Irene Stegun and Milton Abramowitz which contains all the functions I will work on completed with (almost) all the expansions that are needed to implement them; and Handbook of Continued Fractions for Special Functions by Annie Cuyt, Vigdis Brevik Petersen, Brigitte Verdonk, Haakon Waadeland and William B. Jones.
To test the functions I will use, in addition to the Handbook of Mathematical Functions, Sage and the Octave symbolic package to get a reference value.
Here you can find my repository on which I will work.