Re: fix for expt bug

[Ramakrishnan Muthukrishnan <vu3rdd@gmail.com>]

[-- Attachment #1: Type: text/plain, Size: 290 bytes --]

On Wed, Nov 3, 2010 at 9:02 PM, Mark H Weaver <mhw@netris.org> wrote:
> Hi Ramakrishnan,
>
> We're almost there, but you neglected one of the comments I made about
> your previous patch.

Sorry, I should pay more attention. :-(

Attaching the modified patch.
-- 
  Ramakrishnan

[-- Attachment #2: expt-fix.patch --]
[-- Type: text/x-patch, Size: 3102 bytes --]

From 6cca8a66a3daedb551f4f80170966d74b6143ba6 Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan <vu3rdd@gmail.com>
Date: Sun, 31 Oct 2010 23:22:52 +0530
Subject: [PATCH] Adding a case for `expt' when base is negative.

* libguile/numbers.c (scm_expt): Add a case when base is a
  negative number. Also fix the bug for the case when exponent
  is an inexact number.

* test-suite/tests/numbers.test ("expt"): Add tests for the cases
  when base is an integer and exponent is an inexact number, base
  is an integer and exponent is a rational number, base and
  exponent are both inexact and so on.
---
 libguile/numbers.c            |   18 +++++++++++++++++-
 test-suite/tests/numbers.test |   13 ++++++++++++-
 2 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index fbc6cc8..f07f53d 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5445,12 +5445,28 @@ 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))
+  if (scm_is_true (scm_exact_p (y)) && scm_is_integer (y))
     return scm_integer_expt (x, 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)));
     }
+  else if (scm_is_real (x) && scm_is_integer (y) && (scm_to_double (x) < 0))
+    {
+      double result;
+
+      /* If base is negative, expt needs to find (-x)^n = (-1^n) * (x^n).
+         We find x^n and then if n is odd, we also multiply the result
+         with -1. These changes apply only for cases where n is an
+         integer.
+      */
+      result = pow (-scm_to_double (x), scm_to_double (y));
+
+      if (scm_is_true (scm_odd_p (y)))
+        return scm_from_double (-result);
+      else
+        return scm_from_double (result);
+    }
   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..97c58c9 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -2892,7 +2892,18 @@
   (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 "(= 1 (expt 0.0 0.0))" (= 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)))
+  (pass-if (eqv-loosely? (* -2742638075.5 -2742638075.5)
+                         (expt -2742638075.5 2)))
+  (pass-if (eqv-loosely? 4.0 (expt -2.0 2.0)))
+  (pass-if (eqv-loosely? -0.125 (expt -2.0 -3.0)))
+  (pass-if (eqv-loosely? -0.125 (expt -2 -3.0)))
+  (pass-if (eqv-loosely? 0.25 (expt 2.0 -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.0000000000000002+1.7320508075688772i (expt -8 1/3))))
 
 ;;;
 ;;; asinh
-- 
1.7.2.3