From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 8797 invoked by alias); 28 Apr 2009 20:36:29 -0000 Received: (qmail 8785 invoked by uid 22791); 28 Apr 2009 20:36:28 -0000 X-SWARE-Spam-Status: No, hits=-2.0 required=5.0 tests=BAYES_00,SARE_MSGID_LONG40,SPF_PASS X-Spam-Check-By: sourceware.org Received: from ey-out-1920.google.com (HELO ey-out-1920.google.com) (74.125.78.148) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 28 Apr 2009 20:36:20 +0000 Received: by ey-out-1920.google.com with SMTP id 13so178165eye.24 for ; Tue, 28 Apr 2009 13:36:17 -0700 (PDT) MIME-Version: 1.0 Received: by 10.210.43.11 with SMTP id q11mr558755ebq.59.1240950977325; Tue, 28 Apr 2009 13:36:17 -0700 (PDT) In-Reply-To: <99A23E5B-6A5F-476F-B98E-516EA7DB80C4@cern.ch> References: <254EE974-5E07-4EE0-87F2-B2A4FAF11472@cern.ch> <63c059b10904280931o6600dacg80f040ff32a1ac71@mail.gmail.com> <99A23E5B-6A5F-476F-B98E-516EA7DB80C4@cern.ch> Date: Tue, 28 Apr 2009 20:36:00 -0000 Message-ID: <63c059b10904281336m24e2bdc0y5f3d0a94c8cd1d2d@mail.gmail.com> Subject: Re: [Bug-gsl] problem in cubic solver From: "Andrew W. Steiner" To: Lorenzo Moneta Cc: GSL Discuss Mailing List Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2009-q2/txt/msg00003.txt.bz2 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 wr= ote: > Hi Andrew, > > =A0Thank 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). > =A0Thank you for the link, I will try your code and check =A0as well if i= t 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 =A0in ROOT this code, =A0 it would be good if we c= ould > use instead GSL directly. > > Best Regards > > =A0Lorenzo > > > > > 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?r= evision=3D35&view=3Dmarkup >> >> Take care, >> Andrew >> >> On Tue, Apr 28, 2009 at 11:11 AM, Lorenzo Moneta >> wrote: >>> >>> Hello, >>> >>> =A0I have found a problem with solver for cubic equation (both complex = and >>> real one) >>> =A0In 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 Lin= ux >>> 64-bit =A0with gcc 4.3) =A0and 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 >>> >>> =A0This 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. >>> >>> =A0Best Regards >>> >>> =A0Lorenzo Moneta >>> =A0ROOT project =A0(http://root.cern.ch ) >>> =A0CERN >>> >>> >>> >>> >>> >>> >>> _______________________________________________ >>> Bug-gsl mailing list >>> Bug-gsl@gnu.org >>> http://lists.gnu.org/mailman/listinfo/bug-gsl >>> > >