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: [PATCH] Fix bugs in expt and integer-expt Date: Thu, 04 Nov 2010 22:18:34 -0400 Message-ID: <87k4ksv7gl.fsf_-_@yeeloong.netris.org> References: <8762whnc4d.fsf@yeeloong.netris.org> <87r5f4l3yt.fsf@yeeloong.netris.org> <87aalrxim7.fsf@yeeloong.netris.org> <87wrouwf6b.fsf@yeeloong.netris.org> <87bp66rns2.fsf@gnu.org> NNTP-Posting-Host: lo.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="=-=-=" X-Trace: dough.gmane.org 1288923551 17436 80.91.229.12 (5 Nov 2010 02:19:11 GMT) X-Complaints-To: usenet@dough.gmane.org NNTP-Posting-Date: Fri, 5 Nov 2010 02:19:11 +0000 (UTC) To: guile-devel@gnu.org Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Fri Nov 05 03:19:07 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 1PEBtT-0003tq-7X for guile-devel@m.gmane.org; Fri, 05 Nov 2010 03:19:07 +0100 Original-Received: from localhost ([127.0.0.1]:58700 helo=lists.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PEBtS-0007Ej-7H for guile-devel@m.gmane.org; Thu, 04 Nov 2010 22:19:06 -0400 Original-Received: from [140.186.70.92] (port=53469 helo=eggs.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PEBtM-0007DU-Ct for guile-devel@gnu.org; Thu, 04 Nov 2010 22:19:01 -0400 Original-Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1PEBtK-0006mC-V3 for guile-devel@gnu.org; Thu, 04 Nov 2010 22:19:00 -0400 Original-Received: from world.peace.net ([216.204.32.208]:33508) by eggs.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1PEBtK-0006m7-PJ for guile-devel@gnu.org; Thu, 04 Nov 2010 22:18:58 -0400 Original-Received: from turntable.mit.edu ([18.160.0.29] helo=freedomincluded) by world.peace.net with esmtpa (Exim 4.69) (envelope-from ) id 1PEBtD-0002aL-M9; Thu, 04 Nov 2010 22:18:51 -0400 Original-Received: from mhw by freedomincluded with local (Exim 4.69) (envelope-from ) id 1PEBsx-000757-B2; Thu, 04 Nov 2010 22:18:35 -0400 In-Reply-To: <87bp66rns2.fsf@gnu.org> ("Ludovic =?utf-8?Q?Court=C3=A8s=22'?= =?utf-8?Q?s?= message of "Thu, 04 Nov 2010 00:27:25 +0100") 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:11119 Archived-At: --=-=-= Here's my patch to fix some bugs in expt and integer-expt. Please review. Mark --=-=-= Content-Type: text/x-diff Content-Disposition: inline; filename=0001-Fix-bugs-in-expt-and-integer-expt.patch Content-Description: Patch to fix bugs in expt and integer-expt >From 19a24107af814e34fe63d02478bfb9441354625e Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Thu, 4 Nov 2010 22:10:02 -0400 Subject: [PATCH] Fix bugs in expt and integer-expt * libguile/numbers.c (expt): Fix bug that caused expt to throw an exception whenever the base was exact and the exponent was an inexact integer, e.g. (expt 5 6.0). (expt): Fix bug that caused expt to introduce spurious imaginary parts in the result when the base was an inexact negative real and the exponent was an integer, e.g. (expt -1.0 2) (integer-expt, expt): Change behavior of (integer-expt 0 -1), and therefore also (expt 0 -1), to return NaN, per R6RS (actually, R6RS says we should throw an exception or return an "unspecified number object", but for now we use NaN). Formerly we returned 0, per R5RS. R5RS claims that 0^x=0 for all non-zero x, but that's mathematically incorrect, and probably an oversight. (integer-expt): Consistently throw a wrong-argument-type exception when the exponent is inexact. Formerly, it didn't always check this if the base was 0, 1, or -1. * test-suite/tests/numbers.test ("integer-expt", "expt"): Add tests. --- libguile/numbers.c | 41 ++++++++++++++++++++++++++++++---- test-suite/tests/numbers.test | 48 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 79 insertions(+), 10 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index fbc6cc8..b699843 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -1780,10 +1780,20 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0, SCM acc = SCM_I_MAKINUM (1L); SCM_VALIDATE_NUMBER (SCM_ARG1, n); + if (!SCM_I_INUMP (k) && !SCM_BIGP (k)) + SCM_WRONG_TYPE_ARG (2, k); - /* 0^0 == 1 according to R5RS */ - if (scm_is_eq (n, SCM_INUM0) || scm_is_eq (n, acc)) - return scm_is_false (scm_zero_p(k)) ? n : acc; + if (scm_is_true (scm_zero_p (n))) + { + if (scm_is_true (scm_zero_p (k))) /* 0^0 == 1 per R5RS */ + return acc; /* return exact 1, regardless of n */ + else if (scm_is_true (scm_positive_p (k))) + return n; + else /* return NaN for (0 ^ k) for negative k per R6RS */ + return scm_nan (); + } + else if (scm_is_eq (n, acc)) + return acc; else if (scm_is_eq (n, SCM_I_MAKINUM (-1L))) return scm_is_false (scm_even_p (k)) ? n : acc; @@ -5445,8 +5455,29 @@ SCM_DEFINE (scm_expt, "expt", 2, 0, 0, "Return @var{x} raised to the power of @var{y}.") #define FUNC_NAME s_scm_expt { - if (scm_is_true (scm_exact_p (x)) && scm_is_integer (y)) - return scm_integer_expt (x, y); + if (scm_is_integer (y)) + { + if (scm_is_true (scm_exact_p (y))) + return scm_integer_expt (x, y); + else + { + /* Here we handle the case where the exponent is an inexact + integer. We make the exponent exact in order to use + scm_integer_expt, and thus avoid the spurious imaginary + parts that may result from round-off errors in the general + e^(y log x) method below (for example when squaring a large + negative number). In this case, we must return an inexact + result for correctness. We also make the base inexact so + that scm_integer_expt will use fast inexact arithmetic + internally. Note that making the base inexact is not + sufficient to guarantee an inexact result, because + scm_integer_expt will return an exact 1 when the exponent + is 0, even if the base is inexact. */ + return scm_exact_to_inexact + (scm_integer_expt (scm_exact_to_inexact (x), + scm_inexact_to_exact (y))); + } + } else if (scm_is_real (x) && scm_is_real (y) && scm_to_double (x) >= 0.0) { return scm_from_double (pow (scm_to_double (x), scm_to_double (y))); diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index 3c3e14f..6ee4bbc 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -2889,10 +2889,32 @@ (with-test-prefix "expt" (pass-if-exception "non-numeric base" exception:wrong-type-arg (expt #t 0)) - (pass-if "(= 1 (expt 0 0))" (= 1 (expt 0 0))) - (pass-if "(= 1 (expt 0 0.0))" (= 1 (expt 0 0.0))) - (pass-if "(= 1 (expt 0.0 0))" (= 1 (expt 0.0 0))) - (pass-if "(= 1 (expt 0.0 0.0))" (= 1 (expt 0.0 0.0)))) + (pass-if (eqv? 1 (expt 0 0))) + (pass-if (eqv? 1 (expt 0.0 0))) + (pass-if (eqv? 1.0 (expt 0 0.0))) + (pass-if (eqv? 1.0 (expt 0.0 0.0))) + (pass-if (nan? (expt 0 -1))) + (pass-if (nan? (expt 0 -1.0))) + (pass-if (nan? (expt 0.0 -1))) + (pass-if (nan? (expt 0.0 -1.0))) + (pass-if (eqv? 0 (expt 0 3))) + (pass-if (= 0 (expt 0 4.0))) + (pass-if (eqv? 0.0 (expt 0.0 5))) + (pass-if (eqv? 0.0 (expt 0.0 6.0))) + (pass-if (eqv? -2742638075.5 (expt -2742638075.5 1))) + (pass-if (eqv? (* -2742638075.5 -2742638075.5) + (expt -2742638075.5 2))) + (pass-if (eqv? 4.0 (expt -2.0 2.0))) + (pass-if (eqv? -1/8 (expt -2 -3))) + (pass-if (eqv? -0.125 (expt -2.0 -3))) + (pass-if (eqv? -0.125 (expt -2 -3.0))) + (pass-if (eqv? -0.125 (expt -2.0 -3.0))) + (pass-if (eqv? 0.25 (expt 2.0 -2.0))) + (pass-if (eqv? (* -1.0 12398 12398) (expt +12398i 2.0))) + (pass-if (eqv-loosely? +i (expt -1 0.5))) + (pass-if (eqv-loosely? +i (expt -1 1/2))) + (pass-if (eqv-loosely? 1.0+1.7320508075688i (expt -8 1/3)))) + ;;; ;;; asinh @@ -3014,7 +3036,23 @@ (pass-if-exception "2^-inf" exception:wrong-type-arg (integer-expt 2 -inf.0)) (pass-if-exception "2^nan" exception:wrong-type-arg - (integer-expt 2 +nan.0))) + (integer-expt 2 +nan.0)) + + (pass-if (eqv? 1 (integer-expt 0 0))) + (pass-if (eqv? 1 (integer-expt 0.0 0))) + (pass-if (nan? (integer-expt 0 -1))) + (pass-if (nan? (integer-expt 0.0 -1))) + (pass-if (eqv? 0 (integer-expt 0 3))) + (pass-if (eqv? 0.0 (integer-expt 0.0 5))) + (pass-if (eqv? -2742638075.5 (integer-expt -2742638075.5 1))) + (pass-if (eqv? (* -2742638075.5 -2742638075.5) + (integer-expt -2742638075.5 2))) + (pass-if (eqv? 4.0 (integer-expt -2.0 2))) + (pass-if (eqv? -1/8 (integer-expt -2 -3))) + (pass-if (eqv? -0.125 (integer-expt -2.0 -3))) + (pass-if (eqv? 0.25 (integer-expt 2.0 -2))) + (pass-if (eqv? (* -1.0 12398 12398) (integer-expt +12398.0i 2)))) + ;;; ;;; integer-length -- 1.5.6.5 --=-=-=--