Re: [Bug-gsl] problem in cubic solver

Re: [Bug-gsl] problem in cubic solver

[Andrew W. Steiner]

On the topic of cubics, have you tried translating CERNLIB's rrteq3
into C for root? I have found it to give better results in general
than GSL's cubic routine, and thus it might be worth adding a
translation of the CERNLIB code to GSL.

My attempt at this translation is given at line 517 of
http://o2scl.svn.sourceforge.net/viewvc/o2scl/trunk/src/other/poly.cpp?revision=35&view=markup

Take care,
Andrew

On Tue, Apr 28, 2009 at 11:11 AM, Lorenzo Moneta <Lorenzo.Moneta@cern.ch> wrote:
>
> Hello,
>
>  I have found a problem with solver for cubic equation (both complex and real one)
>  In some particular case a NaN is returned, as shown in the attached text example.
> The problem is observed mainly on 64 bit architectures (for example Linux 64-bit  with gcc 4.3)  and not on 32 bit architectures.
>
> This is due to a problem in a sqrt. A patch is attached fixing the problem for the complex routine (gsl_poly_complex_solve_cubic ).
> A similar patch can be probably done also for the real routine
>
>  This patch fails the current test for polynomial in gsl, however in my opinion, this is acceptable, because the test condition is too strict. With a slight change of the test coefficient, you will have a similar failure also with the current version.
>
>  Best Regards
>
>  Lorenzo Moneta
>  ROOT project  (http://root.cern.ch )
>  CERN
>
>
>
>
>
>
> _______________________________________________
> Bug-gsl mailing list
> Bug-gsl@gnu.org
> http://lists.gnu.org/mailman/listinfo/bug-gsl
>

Re: [Bug-gsl] problem in cubic solver

[Lorenzo Moneta]

Hi Andrew,

  Thank you for the link, we have not in Root performed this  
translation. We have currently an implementation very similar to the  
GSL one
(but only finding real roots).
  Thank you for the link, I will try your code and check  as well if  
it gives better numerical stability. In that case it wold be good also  
to have it in GSL.

I remember also, you were working on importing the GSL quartic in GSL.  
Is this still planned ?
We are using presently  in ROOT this code,   it would be good if we  
could use instead GSL directly.

Best Regards

  Lorenzo




On Apr 28, 2009, at 6:31 PM, Andrew W. Steiner wrote:

> On the topic of cubics, have you tried translating CERNLIB's rrteq3
> into C for root? I have found it to give better results in general
> than GSL's cubic routine, and thus it might be worth adding a
> translation of the CERNLIB code to GSL.
>
> My attempt at this translation is given at line 517 of
> http://o2scl.svn.sourceforge.net/viewvc/o2scl/trunk/src/other/poly.cpp?revision=35&view=markup
>
> Take care,
> Andrew
>
> On Tue, Apr 28, 2009 at 11:11 AM, Lorenzo Moneta <Lorenzo.Moneta@cern.ch 
> > wrote:
>>
>> Hello,
>>
>>  I have found a problem with solver for cubic equation (both  
>> complex and real one)
>>  In some particular case a NaN is returned, as shown in the  
>> attached text example.
>> The problem is observed mainly on 64 bit architectures (for example  
>> Linux 64-bit  with gcc 4.3)  and not on 32 bit architectures.
>>
>> This is due to a problem in a sqrt. A patch is attached fixing the  
>> problem for the complex routine (gsl_poly_complex_solve_cubic ).
>> A similar patch can be probably done also for the real routine
>>
>>  This patch fails the current test for polynomial in gsl, however  
>> in my opinion, this is acceptable, because the test condition is  
>> too strict. With a slight change of the test coefficient, you will  
>> have a similar failure also with the current version.
>>
>>  Best Regards
>>
>>  Lorenzo Moneta
>>  ROOT project  (http://root.cern.ch )
>>  CERN
>>
>>
>>
>>
>>
>>
>> _______________________________________________
>> Bug-gsl mailing list
>> Bug-gsl@gnu.org
>> http://lists.gnu.org/mailman/listinfo/bug-gsl
>>

Re: [Bug-gsl] problem in cubic solver

[Andrew W. Steiner]

Ah yes. I seem to remember having some difficulty at the time because
I didn't understand cvs at all. Now that it's a couple years later I'm
older and wiser and should probably return to it this summer. I've
made many improvements since then in all of my polynomial routines,
and it shouldn't be too difficult to just take them from the code I've
already written for o2scl.

On Tue, Apr 28, 2009 at 4:05 PM, Lorenzo Moneta <Lorenzo.Moneta@cern.ch> wrote:
> Hi Andrew,
>
>  Thank you for the link, we have not in Root performed this translation. We
> have currently an implementation very similar to the GSL one
> (but only finding real roots).
>  Thank you for the link, I will try your code and check  as well if it gives
> better numerical stability. In that case it wold be good also to have it in
> GSL.
>
> I remember also, you were working on importing the GSL quartic in GSL. Is
> this still planned ?
> We are using presently  in ROOT this code,   it would be good if we could
> use instead GSL directly.
>
> Best Regards
>
>  Lorenzo
>
>
>
>
> On Apr 28, 2009, at 6:31 PM, Andrew W. Steiner wrote:
>
>> On the topic of cubics, have you tried translating CERNLIB's rrteq3
>> into C for root? I have found it to give better results in general
>> than GSL's cubic routine, and thus it might be worth adding a
>> translation of the CERNLIB code to GSL.
>>
>> My attempt at this translation is given at line 517 of
>>
>> http://o2scl.svn.sourceforge.net/viewvc/o2scl/trunk/src/other/poly.cpp?revision=35&view=markup
>>
>> Take care,
>> Andrew
>>
>> On Tue, Apr 28, 2009 at 11:11 AM, Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
>> wrote:
>>>
>>> Hello,
>>>
>>>  I have found a problem with solver for cubic equation (both complex and
>>> real one)
>>>  In some particular case a NaN is returned, as shown in the attached text
>>> example.
>>> The problem is observed mainly on 64 bit architectures (for example Linux
>>> 64-bit  with gcc 4.3)  and not on 32 bit architectures.
>>>
>>> This is due to a problem in a sqrt. A patch is attached fixing the
>>> problem for the complex routine (gsl_poly_complex_solve_cubic ).
>>> A similar patch can be probably done also for the real routine
>>>
>>>  This patch fails the current test for polynomial in gsl, however in my
>>> opinion, this is acceptable, because the test condition is too strict. With
>>> a slight change of the test coefficient, you will have a similar failure
>>> also with the current version.
>>>
>>>  Best Regards
>>>
>>>  Lorenzo Moneta
>>>  ROOT project  (http://root.cern.ch )
>>>  CERN
>>>
>>>
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> Bug-gsl mailing list
>>> Bug-gsl@gnu.org
>>> http://lists.gnu.org/mailman/listinfo/bug-gsl
>>>
>
>