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.