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.

expint

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.

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.