Re: About special functions reogranization

[Rik]

On 07/19/2017 01:20 AM, Michele wrote:

On 07/14/2017 09:04 PM, Rik wrote:

On 07/13/2017 09:00 AM, address@hidden wrote:

Subject:

About special functions reogranization

From:

Michele <address@hidden>

Date:

07/13/2017 01:49 AM

To:

address@hidden

List-Post:

<mailto:address@hidden>

Content-Transfer-Encoding:

7bit

Precedence:

list

MIME-Version:

1.0

Message-ID:

<address@hidden>

Content-Type:

text/plain; charset=utf-8; format=flowed

Message:

3

Dear all,
I'm Michele Ginesi, the student working on special functions for GSoC2017.
During the last period I saw on this mailing list some discussions about the reorganization of special functions in Octave. I saw there is the idea to use the standard math library of C++ or GSL library.
Unfortunately, as I explained in the last post on my blog [1], I think it is not the case to fully substitute the old (but still the more accurate) Fortran code present in the Amos library.
For example, for what concern Bessel functions, I would  suggest to "unlock" the Amos, since the bug related to them [2] is caused by the fact that if the magnitude of the input is too large, the output is simply not computed (you can find more details on [1], as well as other thoughts relative to incomplete gamma and beta functions).
I hope you can give me some opinions and/or suggestions.
Thank you for your kind attention.

[1] http://gsocspecfun.blogspot.it/2017/07/gnu-scientific-library.html
[2] https://savannah.gnu.org/bugs/?func=detailitem&item_id=48316

Michele,

I read your summary of the state of the bessel, gammainc, and betainc functions with interest.

For the Bessel functions, you state that Matlab allows a complex alpha parameter, but this is not the case (http://www.mathworks.com/help/matlab/ref/besselj.html).  This doesn't detract much from your argument.  The real issue which prevents using the GSL library is that the input itself must be real, while AMOS and Matlab both support a complex input Z.

For gammainc, you assert that there is no "scaled" option.  But, there is an unnormalized version of the function gsl_sf_gamma_inc.  It wouldn't be too difficult to transform the unnormalized output to the scaled output defined by Matlab.  It does involve multiplication by a scaling factor, but underflow potential, which is what the scaled option in Matlab is seeking to prevent, looks much more unlikely in this case.  However, you're still probably right to create your own implementation given that the GSL library requires X to be positive for all versions of the function.  The current implementation of gammainc is in Fortran.  Will you be writing in Fortran, C++, or Octave?  I doubt you can get the performance out of an Octave m-file that you can out of the other two compiled languages, unless there are very few iterative loops in the algorithm.

I agree with your conclusion that betainc is insufficient.  You mention writing an m-file to replace what is now Fortran code.  Again, my concern here is that the performance will not be high enough unless the algorithm does not use iterative loops.

Happy coding,
Rik

Dear Rik,

thanks for your comments. I already worked on gammainc, you can see the details in [1]. The important thing is that Marco, Nir and I implemented it as an .m file, with a subroutine in C++ (you can find this implementation in the bookmark "gammainc" of my repository [2]). Every special function has more or less the same structure: the domain is split into subdomains and different subfunctions are called. Each subfunction implements an algorithm, the most used ones being series expansions and continued fractions. I think that I could continue my work (e.g. betainc) in this way, finding algorithms that are accurate in their domains. Then, with maybe the help of the community, start to translate the subfunctions from .m files to C++ (I'm not autonomous with C++). In this way, at the end of GSoC we will have accurate algorithms, that in the next time can become also fast.

This is a good strategy.  Concentrate on correctness first, and performance second.  It is usually not too difficult to translate an algorithm, once it has been written, from one programming language to another.  And the C++ used by Octave has a lot of features of the m-file language itself which will make the transition even easier.

--Rik