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