improving special functions in GSoC 2017

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

improving special functions in GSoC 2017

Colin Macdonald-2
I added a rough blurb for a GSoC 2017 project:

http://wiki.octave.org/Summer_of_Code_Project_Ideas#Numerical

Motivation: I've seen various poor behaviour from our special functions
over the last couple months.  I think someone mentioned improving all
that would be a good GSoC project.  I agree.

Please link bugs there if you know of others.

Please also add your name as potential mentors: I'm not an expert on
special functions, but I'm willing to help out.

Colin

Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

Marco Caliari-4
> I added a rough blurb for a GSoC 2017 project:
>
> http://wiki.octave.org/Summer_of_Code_Project_Ideas#Numerical
>
> Motivation: I've seen various poor behaviour from our special functions over
> the last couple months.  I think someone mentioned improving all that would
> be a good GSoC project.  I agree.
>
> Please link bugs there if you know of others.
>
> Please also add your name as potential mentors: I'm not an expert on special
> functions, but I'm willing to help out.

Dear Colin,

I'm not an expert too, but please add my name as potential mentor. I
completely agree with you that "Special functions are expected by users to
just work".

Best regards,

Marco

Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

Colin Macdonald-2
On 25/06/16 02:54, Marco Caliari wrote:
> I'm not an expert too, but please add my name as potential mentor. I
> completely agree with you that "Special functions are expected by users
> to just work".

Great, added.  Thank you (or whoever) added other more bug links.

Colin

Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

Colin Macdonald-2
On 25/06/16 02:54, Marco Caliari wrote:
> "Special functions are expected by users to just work".

So I just found that our "besselj" gives NaNs for values like "1e10".

 >> besselj(1, 1e9)
ans = -5.21042264155388e-06
 >> besselj(1, 1e10)
ans = NaN + NaNi


https://savannah.gnu.org/bugs/index.php?48316


That's not some wild exotic function: its a common Bessel function!


Why are we not just calling some free-licensed library that nails these
things?  For example, SciPy does:

 >>> scipy.special.jn(1,1e9)
-5.2104226415538769e-06
 >>> scipy.special.jn(1,1e10)
-7.676506113846571e-06


Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

jbect
Le 26/06/2016 à 08:08, Colin Macdonald a écrit :

> On 25/06/16 02:54, Marco Caliari wrote:
>> "Special functions are expected by users to just work".
>
> So I just found that our "besselj" gives NaNs for values like "1e10".
>
> >> besselj(1, 1e9)
> ans = -5.21042264155388e-06
> >> besselj(1, 1e10)
> ans = NaN + NaNi
>
>
> https://savannah.gnu.org/bugs/index.php?48316
>
>
> That's not some wild exotic function: its a common Bessel function!
>
>
> Why are we not just calling some free-licensed library that nails
> these things?  For example, SciPy does:
>
> >>> scipy.special.jn(1,1e9)
> -5.2104226415538769e-06
> >>> scipy.special.jn(1,1e10)
> -7.676506113846571e-06
>
>

Another possibility is Boost :
http://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/special.html.


Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

Colin Macdonald-2
On 25/06/16 23:22, Julien Bect wrote:
> Another possibility is Boost :
> http://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/special.html.

And another possibility is math.h: "man jn" :)

Colin

Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

John W. Eaton
Administrator
In reply to this post by Colin Macdonald-2
On 06/26/2016 02:08 AM, Colin Macdonald wrote:

> On 25/06/16 02:54, Marco Caliari wrote:
>> "Special functions are expected by users to just work".
>
> So I just found that our "besselj" gives NaNs for values like "1e10".
>
>  >> besselj(1, 1e9)
> ans = -5.21042264155388e-06
>  >> besselj(1, 1e10)
> ans = NaN + NaNi
>
>
> https://savannah.gnu.org/bugs/index.php?48316
>
>
> That's not some wild exotic function: its a common Bessel function!
>
>
> Why are we not just calling some free-licensed library that nails these
> things?

Amos was the best thing I knew of at the time.  If there is something
better now, or a standard library solution, then perhaps we should be
using it.  Propose a patch and a test suite.

jwe



Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

Robert T. Short
On 06/25/2016 11:52 PM, John W. Eaton wrote:

> On 06/26/2016 02:08 AM, Colin Macdonald wrote:
>> On 25/06/16 02:54, Marco Caliari wrote:
>>> "Special functions are expected by users to just work".
>>
>> So I just found that our "besselj" gives NaNs for values like "1e10".
>>
>>  >> besselj(1, 1e9)
>> ans = -5.21042264155388e-06
>>  >> besselj(1, 1e10)
>> ans = NaN + NaNi
>>
>>
>> https://savannah.gnu.org/bugs/index.php?48316
>>
>>
>> That's not some wild exotic function: its a common Bessel function!
>>
>>
>> Why are we not just calling some free-licensed library that nails these
>> things?
>
> Amos was the best thing I knew of at the time.  If there is something
> better now, or a standard library solution, then perhaps we should be
> using it.  Propose a patch and a test suite.
>
> jwe
>
>
>
>
>
There was some discussion of this a few years ago (something like 2011
or 2012).  At the time Amos was still the best option that we could
find.  I don't remember all the details, but other libraries (including
boost) didn't handle things like complex arguments or other some other
crucial element.  Hopefully things have changed since then, but that was
the conclusion at that time.  I was not able to find that thread so my
memory is not refreshed.  I put in a bunch of tests at the time we had
the discussion.  It seems, though, that I didn't do any for complex
arguments.  I have no clue how I missed that.

I played around with python (scipy.special) a bit and didn't see any
obvious issues so it might be worth looking at the library they use.

Bob


Reply | Threaded
Open this post in threaded view
|

Re: improving special functions in GSoC 2017

Marco Caliari-4
On Sun, 26 Jun 2016, Robert T. Short wrote:

> On 06/25/2016 11:52 PM, John W. Eaton wrote:
>> On 06/26/2016 02:08 AM, Colin Macdonald wrote:
>>> On 25/06/16 02:54, Marco Caliari wrote:
>>>> "Special functions are expected by users to just work".
>>>
>>> So I just found that our "besselj" gives NaNs for values like "1e10".
>>>
>>>  >> besselj(1, 1e9)
>>> ans = -5.21042264155388e-06
>>>  >> besselj(1, 1e10)
>>> ans = NaN + NaNi
>>>
>>>
>>> https://savannah.gnu.org/bugs/index.php?48316
>>>
>>>
>>> That's not some wild exotic function: its a common Bessel function!
>>>
>>>
>>> Why are we not just calling some free-licensed library that nails these
>>> things?
>>
>> Amos was the best thing I knew of at the time.  If there is something
>> better now, or a standard library solution, then perhaps we should be using
>> it.  Propose a patch and a test suite.
>>
>> jwe
>>
>>
>>
>>
>>
> There was some discussion of this a few years ago (something like 2011 or
> 2012).  At the time Amos was still the best option that we could find.  I
> don't remember all the details, but other libraries (including boost) didn't
> handle things like complex arguments or other some other crucial element.

Yes, for instance boost, GSL, and Cephes require both the arguments of
gammainc to be nonnegative. Matlab allows x negative, and in this case the
result is complex.

> Hopefully things have changed since then, but that was the conclusion at that
> time.  I was not able to find that thread so my memory is not refreshed.  I
> put in a bunch of tests at the time we had the discussion.  It seems, though,
> that I didn't do any for complex arguments.  I have no clue how I missed
> that.

Should we start with *the* list of special functions with the domain of
application?

specfun1: A1 x B1 x ... x Z1 -> complex_field

where the domain is taken from the equivalent Matlab function.

Marco

Reply | Threaded
Open this post in threaded view
|

Special Functions [was: improving special functions in GSoC 2017]

Robert T. Short
On 06/28/2016 02:29 AM, Marco Caliari wrote:

> On Sun, 26 Jun 2016, Robert T. Short wrote:
>
>> On 06/25/2016 11:52 PM, John W. Eaton wrote:
>>> On 06/26/2016 02:08 AM, Colin Macdonald wrote:
>>>> On 25/06/16 02:54, Marco Caliari wrote:
>>>>> "Special functions are expected by users to just work".
>>>>
>>>> So I just found that our "besselj" gives NaNs for values like "1e10".
>>>>
>>>>  >> besselj(1, 1e9)
>>>> ans = -5.21042264155388e-06
>>>>  >> besselj(1, 1e10)
>>>> ans = NaN + NaNi
>>>>
>>>>
>>>> https://savannah.gnu.org/bugs/index.php?48316
>>>>
>>>>
>>>> That's not some wild exotic function: its a common Bessel function!
>>>>
>>>>
>>>> Why are we not just calling some free-licensed library that nails
>>>> these
>>>> things?
>>>
>>> Amos was the best thing I knew of at the time.  If there is
>>> something better now, or a standard library solution, then perhaps
>>> we should be using it.  Propose a patch and a test suite.
>>>
>>> jwe
>>>
>>>
>>>
>>>
>>>
>> There was some discussion of this a few years ago (something like
>> 2011 or 2012).  At the time Amos was still the best option that we
>> could find.  I don't remember all the details, but other libraries
>> (including boost) didn't handle things like complex arguments or
>> other some other crucial element.
>
> Yes, for instance boost, GSL, and Cephes require both the arguments of
> gammainc to be nonnegative. Matlab allows x negative, and in this case
> the result is complex.
>
>> Hopefully things have changed since then, but that was the conclusion
>> at that time.  I was not able to find that thread so my memory is not
>> refreshed.  I put in a bunch of tests at the time we had the
>> discussion.  It seems, though, that I didn't do any for complex
>> arguments.  I have no clue how I missed that.
>
> Should we start with *the* list of special functions with the domain
> of application?
>
> specfun1: A1 x B1 x ... x Z1 -> complex_field
>
> where the domain is taken from the equivalent Matlab function.
>
> Marco
>
>
Marco, jwe

Attached is a list of MATLAB special functions.  I started with

http://www.mathworks.com/help/matlab/special-functions-1.html

I think it is correct, but I could easily have some mistakes.

Bob


ListOfSpecialFunctions.pdf (38K) Download Attachment