From mboxrd@z Thu Jan 1 00:00:00 1970 Path: news.gmane.org!not-for-mail From: Mark H Weaver Newsgroups: gmane.lisp.guile.devel Subject: Re: fix for expt bug Date: Sun, 31 Oct 2010 20:03:46 -0400 Message-ID: <8762whnc4d.fsf@yeeloong.netris.org> References: NNTP-Posting-Host: lo.gmane.org Mime-Version: 1.0 Content-Type: text/plain; charset=utf-8 Content-Transfer-Encoding: quoted-printable X-Trace: dough.gmane.org 1288569865 11964 80.91.229.12 (1 Nov 2010 00:04:25 GMT) X-Complaints-To: usenet@dough.gmane.org NNTP-Posting-Date: Mon, 1 Nov 2010 00:04:25 +0000 (UTC) Cc: guile-devel To: Ramakrishnan Muthukrishnan Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Mon Nov 01 01:04:21 2010 Return-path: Envelope-to: guile-devel@m.gmane.org Original-Received: from lists.gnu.org ([199.232.76.165]) by lo.gmane.org with esmtp (Exim 4.69) (envelope-from ) id 1PChsp-0000QT-KB for guile-devel@m.gmane.org; Mon, 01 Nov 2010 01:04:19 +0100 Original-Received: from localhost ([127.0.0.1]:57939 helo=lists.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PChso-00010T-On for guile-devel@m.gmane.org; Sun, 31 Oct 2010 20:04:18 -0400 Original-Received: from [140.186.70.92] (port=42187 helo=eggs.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PChsf-000109-2A for guile-devel@gnu.org; Sun, 31 Oct 2010 20:04:10 -0400 Original-Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1PChsd-0002Ld-Ck for guile-devel@gnu.org; Sun, 31 Oct 2010 20:04:08 -0400 Original-Received: from world.peace.net ([216.204.32.208]:37511) by eggs.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1PChsd-0002G0-4W for guile-devel@gnu.org; Sun, 31 Oct 2010 20:04:07 -0400 Original-Received: from c-98-217-70-13.hsd1.ma.comcast.net ([98.217.70.13] helo=freedomincluded) by world.peace.net with esmtpa (Exim 4.69) (envelope-from ) id 1PChsL-0003j2-6J; Sun, 31 Oct 2010 20:03:49 -0400 Original-Received: from mhw by freedomincluded with local (Exim 4.69) (envelope-from ) id 1PChsJ-0000Y7-1W; Sun, 31 Oct 2010 20:03:47 -0400 In-Reply-To: (Ramakrishnan Muthukrishnan's message of "Sun, 31 Oct 2010 23:42:49 +0530") User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/23.1 (gnu/linux) X-detected-operating-system: by eggs.gnu.org: GNU/Linux 2.6 (newer, 3) X-BeenThere: guile-devel@gnu.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Developers list for Guile, the GNU extensibility library" List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Original-Sender: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Errors-To: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Xref: news.gmane.org gmane.lisp.guile.devel:11093 Archived-At: Ramakrishnan, Your fix is incorrect. You have assumed that (-a)^b =3D -(a^b), but this is true only if b is an odd integer. A correct relation is (-a)^b =3D (-1)^b * a^b. As reported by Ludovic: scheme@(guile-user)> (expt -2742638075.5 2) $8 =3D 7.52206361318235e18-1842.31337891184i scheme@(guile-user)> (* -2742638075.5 -2742638075.5) $9 =3D 7.52206361318235e18 The reported bug is simply due to roundoff error. You must know something about complex numbers and specifically complex exponentials to understand what's happening here. The spurious imaginary part is actually quite miniscule compared with the real part. If one looks at these results in polar form, the complex argument (angle) should ideally be an exact zero, but is instead an inexact number very close to 0 (about 2.45e-16 radians in this case). The most straightforward solution to this bug is to support a few special cases, such as: (-a)^b =3D -(a^b) (where a is positive real and b is an odd integer) (-a)^b =3D a^b (where a is positive real and b is an even integer) A more elaborate solution would involve supporting a polar representation of complex numbers where the angle can be an exact rational times pi, but this is almost certainly not worth it. Best, Mark Ramakrishnan Muthukrishnan wrote: > > > expt currently wrongly prints the results for any negative base > numbers. The attached patch (also attached with the bug tracker) has a > fix and also adds a test case. > > This is my first patch to guile, so please correct me if I have done > something wrong. I will be happy to rework it. I am yet to recieve the > assignment papers from the FSF (but they have been despatched already > by the FSF), I will mail them back as soon as I get it. > > thanks > --=20 > =C2=A0 Ramakrishnan > > From c23939a1c2b7bfe1b3cf20abb6a7b431699281a1 Mon Sep 17 00:00:00 2001 > From: Ramakrishnan Muthukrishnan > Date: Sun, 31 Oct 2010 23:22:52 +0530 > Subject: [PATCH] Fix for bug #31464. expt needs to treat negative bases s= pecially. > Also adding test-suite cases for expt. > > * libguile/numbers.c: if base is negative, find the absolute value > of the base, find log and then later reapply the negative sign. > The bug was caused by the fact that log of a negative number is > undefined. > > * test-suite/tests/numbers.test: Two new test cases for expt. For > cases where the base is negative and the power to be raised is > not an integer, the result should be a complex number. > --- > libguile/numbers.c | 11 +++++++++++ > test-suite/tests/numbers.test | 7 ++++++- > 2 files changed, 17 insertions(+), 1 deletions(-) > > diff --git a/libguile/numbers.c b/libguile/numbers.c > index fbc6cc8..07ae95d 100644 > --- a/libguile/numbers.c > +++ b/libguile/numbers.c > @@ -5451,6 +5451,17 @@ SCM_DEFINE (scm_expt, "expt", 2, 0, 0, > { > return scm_from_double (pow (scm_to_double (x), scm_to_double (y))= ); > } > + else if (scm_is_real (x) && > + scm_is_real (y) && > + (scm_to_double (x) < 0) && > + !SCM_FRACTIONP (y)) > + { > + SCM x_abs; > + > + x_abs =3D scm_abs (x); > + return scm_product(scm_from_double (-1.0), > + scm_exp (scm_product (scm_log (x_abs), y))); > + } > else > return scm_exp (scm_product (scm_log (x), y)); > } > diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test > index 3c3e14f..3d53bbe 100644 > --- a/test-suite/tests/numbers.test > +++ b/test-suite/tests/numbers.test > @@ -2892,7 +2892,12 @@ > (pass-if "(=3D 1 (expt 0 0))" (=3D 1 (expt 0 0))) > (pass-if "(=3D 1 (expt 0 0.0))" (=3D 1 (expt 0 0.0))) > (pass-if "(=3D 1 (expt 0.0 0))" (=3D 1 (expt 0.0 0))) > - (pass-if "(=3D 1 (expt 0.0 0.0))" (=3D 1 (expt 0.0 0.0)))) > + (pass-if "(=3D 1 (expt 0.0 0.0))" (=3D 1 (expt 0.0 0.0))) > + ;; tests for non zero values of base and exponent. > + (pass-if "(eqv-loosely? -2742638075.5 (expt -2742638075.5 1)" > + (eqv-loosely? -2742638075.5 (expt -2742638075.5 1))) > + (pass-if "(eqv-loosely? 1.0000000000000002+1.7320508075688772i (expt -= 8 1/3)" > + (eqv-loosely? 1.0000000000000002+1.7320508075688772i (expt -8= 1/3)))) >=20=20 > ;;; > ;;; asinh