fix for expt bug

fix for expt bug

[Ramakrishnan Muthukrishnan]

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

Hi,

As reported here:

<https://savannah.gnu.org/bugs/index.php?31464>

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
-- 
  Ramakrishnan

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

From c23939a1c2b7bfe1b3cf20abb6a7b431699281a1 Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan <vu3rdd@gmail.com>
Date: Sun, 31 Oct 2010 23:22:52 +0530
Subject: [PATCH] Fix for bug #31464. expt needs to treat negative bases specially.
 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 = 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 "(= 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)"
+           (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))))
 
 ;;;
 ;;; asinh
-- 
1.7.2.3

Re: fix for expt bug

[Mark H Weaver]

Ramakrishnan,

Your fix is incorrect.  You have assumed that (-a)^b = -(a^b),
but this is true only if b is an odd integer.  A correct relation
is (-a)^b = (-1)^b * a^b.

As reported by Ludovic:

scheme@(guile-user)> (expt -2742638075.5 2)
$8 = 7.52206361318235e18-1842.31337891184i
scheme@(guile-user)> (* -2742638075.5 -2742638075.5)
$9 = 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 = -(a^b)   (where a is positive real and b is an odd integer)
  (-a)^b =   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 <vu3rdd@gmail.com> wrote:
> <https://savannah.gnu.org/bugs/index.php?31464>
>
> 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
> -- 
>   Ramakrishnan
>
> From c23939a1c2b7bfe1b3cf20abb6a7b431699281a1 Mon Sep 17 00:00:00 2001
> From: Ramakrishnan Muthukrishnan <vu3rdd@gmail.com>
> Date: Sun, 31 Oct 2010 23:22:52 +0530
> Subject: [PATCH] Fix for bug #31464. expt needs to treat negative bases specially.
>  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 = 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 "(= 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)"
> +           (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))))
>  
>  ;;;
>  ;;; asinh

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

On Mon, Nov 1, 2010 at 5:33 AM, Mark H Weaver <mhw@netris.org> wrote:
> Ramakrishnan,
>
> Your fix is incorrect.  You have assumed that (-a)^b = -(a^b),

Mark,

Yes, you are right. Not sure what I was thinking when I made the patch. :-(

> 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 = -(a^b)   (where a is positive real and b is an odd integer)
>  (-a)^b =   a^b    (where a is positive real and b is an even integer)

Thanks. I will take a dig at doing this solution.

-- 
  Ramakrishnan

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

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

I reworked the patch a bit as per Mark's comments. It takes care of
the sign part and handles integer power values properly. It also fixes
another bug with the current expt implementation where base is an
integer and the power value is not (eg: (expt 2 2.0) gives an error in
the current git master. The patch fixes this too.)

Also adds some test cases. Please review.

Thanks
Ramakrishnan

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

From 025bde78d4c199dee1d2857e913d69ce4d7c2e59 Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan <vu3rdd@gmail.com>
Date: Sun, 31 Oct 2010 23:22:52 +0530
Subject: [PATCH] Fix for bug #31464. expt needs to treat negative bases specially.

* libguile/numbers.c: 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.

* 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            |   20 +++++++++++++++++++-
 test-suite/tests/numbers.test |   15 ++++++++++++++-
 2 files changed, 33 insertions(+), 2 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index fbc6cc8..5bbf4b0 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5445,12 +5445,30 @@ 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 (x)) &&
+      scm_is_true (scm_exact_p (y)) &&
+      !SCM_FRACTIONP (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_real (y) &&
+           (scm_to_double (x) < 0) &&
+           !SCM_FRACTIONP (y))
+    {
+      SCM x_abs, result;
+
+      x_abs = scm_abs (x);
+      result = scm_from_double (pow (scm_to_double (x_abs), scm_to_double (y)));
+
+      if (scm_is_true (scm_equal_p (y, scm_floor (y))) &&
+          scm_is_true (scm_odd_p (y)))
+        return scm_product (scm_from_double (-1.0), result);
+      else
+        return 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..637449f 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -2892,7 +2892,20 @@
   (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))"
+            (eqv-loosely? -2742638075.5 (expt -2742638075.5 1)))
+  (pass-if "(eqv-loosely? 4.0 (expt -2.0 2.0))"
+           (eqv-loosely? 4.0 (expt -2.0 2.0)))
+  (pass-if "(eqv-loosely? -0.125 (expt -2.0 -3.0))"
+           (eqv-loosely? -0.125 (expt -2.0 -3.0)))
+  (pass-if "(eqv-loosely? -0.125 (expt -2 -3.0))"
+           (eqv-loosely? -0.125 (expt -2 -3.0)))
+  (pass-if "(eqv-loosely? 0.25 (expt 2.0 -2.0))"
+           (eqv-loosely? 0.25 (expt 2.0 -2.0)))
+  (pass-if "(eqv-loosely? 1.0000000000000002+1.7320508075688772i (expt -8 1/3))"
+            (eqv-loosely? 1.0000000000000002+1.7320508075688772i (expt -8 1/3))))
 
 ;;;
 ;;; asinh
-- 
1.7.2.3

Re: fix for expt bug

[Mark H Weaver]

Hi Ramakrishnan,

Thanks for the revised patch.  There are still some problems:

> diff --git a/libguile/numbers.c b/libguile/numbers.c
> index fbc6cc8..5bbf4b0 100644
> --- a/libguile/numbers.c
> +++ b/libguile/numbers.c
> @@ -5445,12 +5445,30 @@ 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 (x)) &&
> +      scm_is_true (scm_exact_p (y)) &&
> +      !SCM_FRACTIONP (y))
>      return scm_integer_expt (x, y);

This is an improvement, but isn't quite right.  scm_integer_expt
requires that y be an exact integer, but x is allowed to be any scheme
number whatsoever.  So the "scm_exact_p (x)" doesn't belong.  It should
simply be changed to "scm_exact_p (y)" instead.

The other problem is that !SCM_FRACTIONP is not the right test.
Although it is currently true that the only exact numbers in guile are
integers and rationals, there's no guarantee that other exact numbers
won't be added in the future.

The test above should be this:

    if (scm_is_true (scm_exact_p (y)) && scm_is_integer (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_real (y) &&
> +           (scm_to_double (x) < 0) &&
> +           !SCM_FRACTIONP (y))

The !SCM_FRACTIONP (y) does not belong here.  Instead, you should test
scm_is_integer (y), which will also make the scm_is_real (y) test
redundant.  As is, your code will produce the wrong answer if y is an
inexact floating-point number such as 0.5.  SCM_FRACTIONP returns true
only for exact rational numbers.

So the above test should be:

  else if (scm_is_real (x) && scm_is_integer (y) && scm_to_double (x) < 0)

> +    {
> +      SCM x_abs, result;
> +
> +      x_abs = scm_abs (x);

You've already established that x is a negative real, so it is wasteful
to do another test in scm_abs.  Just negate x.

Also, it's inefficient to create so many intermediate floating-point SCM
values, since it involves allocating heap cells.  You don't need x_abs,
and "result" can be a C double.  See my suggested code below.

> +      result = scm_from_double (pow (scm_to_double (x_abs), scm_to_double (y)));
> +
> +      if (scm_is_true (scm_equal_p (y, scm_floor (y))) &&

The test above is better replaced by scm_is_integer (y), but it's needed
in the outer conditional anyway, so you can drop it here.

> +          scm_is_true (scm_odd_p (y)))
> +        return scm_product (scm_from_double (-1.0), result);
> +      else
> +        return result;
> +    }

In summary, I think the new clause should be something like: (UNTESTED!)

  else if (scm_is_real (x) && scm_is_integer (y) && scm_to_double (x) < 0)
  {
    double result;

    result = pow (-scm_to_double (x), scm_to_double (y));
    if (scm_is_true (scm_odd_p (y)))
      result = -result;
    return scm_from_double (result);
  }

     Best,
      Mark

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

On Tue, Nov 2, 2010 at 10:25 AM, Mark H Weaver <mhw@netris.org> wrote:
> Hi Ramakrishnan,
>
> Thanks for the revised patch.  There are still some problems:
> [...]
>
> This is an improvement, but isn't quite right.  scm_integer_expt
> requires that y be an exact integer, but x is allowed to be any scheme
> number whatsoever.  So the "scm_exact_p (x)" doesn't belong.  It should
> simply be changed to "scm_exact_p (y)" instead.
>
> The other problem is that !SCM_FRACTIONP is not the right test.
> Although it is currently true that the only exact numbers in guile are
> integers and rationals, there's no guarantee that other exact numbers
> won't be added in the future.
>
> The test above should be this:
>
>    if (scm_is_true (scm_exact_p (y)) && scm_is_integer (y))

Thanks Mark, for the comments. I am learning a lot. I will make the
changes and test them.

-- 
  Ramakrishnan

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

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

Here is the updated patch.

thanks
Ramakrishnan

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

From e320d79c8f3cd8b7ddcb9c2d13356e34a3346cfe Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan <vu3rdd@gmail.com>
Date: Sun, 31 Oct 2010 23:22:52 +0530
Subject: [PATCH] Fix for bug #31464. expt needs to treat negative bases specially.

* libguile/numbers.c: 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.

* 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            |   13 ++++++++++++-
 test-suite/tests/numbers.test |   15 ++++++++++++++-
 2 files changed, 26 insertions(+), 2 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index fbc6cc8..4b81d98 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5445,12 +5445,23 @@ 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;
+      
+      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..637449f 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -2892,7 +2892,20 @@
   (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))"
+            (eqv-loosely? -2742638075.5 (expt -2742638075.5 1)))
+  (pass-if "(eqv-loosely? 4.0 (expt -2.0 2.0))"
+           (eqv-loosely? 4.0 (expt -2.0 2.0)))
+  (pass-if "(eqv-loosely? -0.125 (expt -2.0 -3.0))"
+           (eqv-loosely? -0.125 (expt -2.0 -3.0)))
+  (pass-if "(eqv-loosely? -0.125 (expt -2 -3.0))"
+           (eqv-loosely? -0.125 (expt -2 -3.0)))
+  (pass-if "(eqv-loosely? 0.25 (expt 2.0 -2.0))"
+           (eqv-loosely? 0.25 (expt 2.0 -2.0)))
+  (pass-if "(eqv-loosely? 1.0000000000000002+1.7320508075688772i (expt -8 1/3))"
+            (eqv-loosely? 1.0000000000000002+1.7320508075688772i (expt -8 1/3))))
 
 ;;;
 ;;; asinh
-- 
1.7.2.3

Re: fix for expt bug

[Mark H Weaver]

Hi Ramakrishnan,

The code in your latest patch looks good to me, though the commit
message has some problems, and I'd add more test cases:

> * libguile/numbers.c: 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.

In the equation above, -x^n parses as -(x^n), but it should be (-x)^n.
Same problem with (-1^n).  Also, the absolute value should be removed.
It is superfluous for the case you are handling (x>0), and for other
cases it is erroneous.

The equation above should be:  (-x)^n = (-1)^n * x^n

Also, I would add a couple more test cases:

* Test the case from Ludovic's original bug report:
  (expt -2742638075.5 2) should loosely equal
  7.52206361318235e18

* Test that (expt -1 0.5) is loosely equal to +i.
  (I believe this would have failed with your second patch)

    Thanks,
      Mark

Re: fix for expt bug

[Ludovic Courtès]

Hello!

Thank you Ramakrishnan for the patch, and thanks Mark for the detailed
review!  It looks like we’re getting close.  :-)

In addition to Mark’s latest comments, just a few cosmetic notes:

Ramakrishnan Muthukrishnan <vu3rdd@gmail.com> writes:

> * libguile/numbers.c: 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.

Please follow GNU change log conventions, i.e., describe the changes,
not the rationale (info "(standards) Change Logs").  Provide appropriate
explanations in ‘scm_expt’, though.

> * 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.

Ditto: list the new test case names.

> +  (pass-if "(eqv-loosely? -2742638075.5 (expt -2742638075.5 1))"
> +            (eqv-loosely? -2742638075.5 (expt -2742638075.5 1)))

The first argument to ‘pass-if’ can be omitted altogether.

Thanks,
Ludo’.

Re: fix for expt bug

[Mark H Weaver]

Hi Ramakrishnan,

I noticed one more problem with your commit message.  It should mention
the other included bug fix, for when the exponent is an inexact integer,
e.g. (expt 2 2.0).

   Thanks,
     Mark

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

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

Thanks again, Mark and Ludovic.

Attached is an updated patch.

thanks
-- 
  Ramakrishnan

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

From a1dd2da8562ddeb2052f2994ad0302bcc8d5d1a2 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..e4b03aa 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

Re: fix for expt bug

[Mark H Weaver]

Hi Ramakrishnan,

We're almost there, but you neglected one of the comments I made about
your previous patch.

> +      /* 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.
> +      */

I repeat: In the equation above, -x^n parses as -(x^n), but it should be
(-x)^n.  Same problem with (-1^n).  Exponentiation has higher precedence
than negation.

The equation above should be:  (-x)^n = (-1)^n * x^n

Feel free to add more parentheses if you like, but the ones I have
included are essential and must not be removed.

Everything else looks good to me.

    Thanks,
      Mark

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

[-- 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

Re: fix for expt bug

[Mark H Weaver]

Ramakrishnan and others,

I just realized that there is a better way to fix these bugs.  We don't
need a new top-level case in expt after all.  Instead, we generalize the
scm_integer_expt case to support inexact integer exponents.

Within that case, if the exponent is an inexact integer, then we make it
exact and make the base inexact, and then call scm_integer_expt.

I think this strategy is better for these reasons:

* Fewer top-level cases to check, thus making expt faster.
* More precise answers in the case of inexact integer exponents.
* More cases handled without spurious imaginary parts turning up, for
  example (expt +234234234.5i 4.0)
* Simpler code

Something like the code below.  (WARNING: UNTESTED!)

What do you all think?  Can you see any problems with this strategy?

    Best,
     Mark


index fbc6cc8..f38772e 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5445,8 +5445,14 @@ 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
+	return 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)));

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

On Wed, Nov 3, 2010 at 9:52 PM, Mark H Weaver <mhw@netris.org> wrote:
> Ramakrishnan and others,
>
> I just realized that there is a better way to fix these bugs.  We don't
> need a new top-level case in expt after all.  Instead, we generalize the
> scm_integer_expt case to support inexact integer exponents.
>
> Within that case, if the exponent is an inexact integer, then we make it
> exact and make the base inexact, and then call scm_integer_expt.

Mark,

Why do we need to convert the base to inexact? is there any problem if
they are just as it is and we convert only the exponent to exact when
they are exact?

-- 
  Ramakrishnan

Re: fix for expt bug

[Ramakrishnan Muthukrishnan]

On Wed, Nov 3, 2010 at 11:23 PM, Ramakrishnan Muthukrishnan
<vu3rdd@gmail.com> wrote:
> On Wed, Nov 3, 2010 at 9:52 PM, Mark H Weaver <mhw@netris.org> wrote:
>> Ramakrishnan and others,
>>
>> I just realized that there is a better way to fix these bugs.  We don't
>> need a new top-level case in expt after all.  Instead, we generalize the
>> scm_integer_expt case to support inexact integer exponents.
>>
>> Within that case, if the exponent is an inexact integer, then we make it
>> exact and make the base inexact, and then call scm_integer_expt.
>
> Mark,
>
> Why do we need to convert the base to inexact? is there any problem if
> they are just as it is and we convert only the exponent to exact when
> they are exact?

scheme@(guile-user)> (integer-expt (exact->inexact 3/2) (inexact->exact 4.0))
$17 = 5.0625
scheme@(guile-user)> (integer-expt 3/2 (inexact->exact 4.0))
$18 = 81/16
scheme@(guile-user)> (integer-expt (exact->inexact 3/2) (inexact->exact 4.0))
$19 = 5.0625

We want an output representation in an inexact form. For that reason,
we would want to do exact to inexact conversion of base. Is that
correct?

Mark, I guess your patch solves the problem in a much more efficient way.

thanks
-- 
  Ramakrishnan

Re: fix for expt bug

[Ludovic Courtès]

Hi Mark,

Mark H Weaver <mhw@netris.org> writes:

> I just realized that there is a better way to fix these bugs.  We don't
> need a new top-level case in expt after all.  Instead, we generalize the
> scm_integer_expt case to support inexact integer exponents.

You mean “inexact number”, right?

The doc for ‘integer-expt’ makes the same mistake:

  "Return @var{n} raised to the power @var{k}.  @var{k} must be an\n"
  "exact integer, @var{n} can be any number.\n"

However, an integer is necessarily an exact number, whereas the converse
isn’t true:

--8<---------------cut here---------------start------------->8---
scheme@(guile-user) [25]> (exact? 2/3)
$4 = #t
scheme@(guile-user) [25]> (integer? 2/3)
$5 = #f
--8<---------------cut here---------------end--------------->8---

In actuality, ‘integer-expt’ expects K to be an integer (fixnum or
bignum).

So I suspect that your patch doesn’t work because ‘inexact->exact’
returns something that’s obviously not necessarily an integer:

--8<---------------cut here---------------start------------->8---
scheme@(guile-user) [25]> (inexact->exact 2.3)
$6 = 2589569785738035/1125899906842624
--8<---------------cut here---------------end--------------->8---

Does that make sense?

Thanks,
Ludo’.

Re: fix for expt bug

[Mark H Weaver]

ludo@gnu.org (Ludovic Courtès) writes:
>> I just realized that there is a better way to fix these bugs.  We don't
>> need a new top-level case in expt after all.  Instead, we generalize the
>> scm_integer_expt case to support inexact integer exponents.
>
> You mean “inexact number”, right?

No, I meant "inexact integer", for example 3.0.

> However, an integer is necessarily an exact number,

No, (integer? 3.0) returns #t, as it should, according to R5RS.
R5RS's description of "integer?" gives this precise example, and
guile's implementation agrees.

> So I suspect that your patch doesn’t work because ‘inexact->exact’
> returns something that’s obviously not necessarily an integer:

I'm only talking about inexact integers such as 3.0, so that

  (expt 2 3.0) ==> (integer-expt 2.0 3) ==> 8.0

However, I have since realized that it is not enough to convert the base
to inexact; I must also convert integer-expt's result to inexact,
because there is one case where making the base inexact is not enough:

  (expt 3 0.0) ==> (integer-expt 3.0 0) ==> 1

Though (expt 3 0.0) should reduce to 1.0.  So the code needs to coerce
the result of integer-expt to inexact.

I am working on a patch to fix these and some other problems.

    Mark

Re: fix for expt bug

[Ludovic Courtès]

Hi Mark!

Mark H Weaver <mhw@netris.org> writes:

> No, (integer? 3.0) returns #t, as it should, according to R5RS.
> R5RS's description of "integer?" gives this precise example, and
> guile's implementation agrees.

Damn, I had never realized that, shame on me.

Ludo’.

Re: fix for expt bug

[Andy Wingo]

Hi,

On Mon 08 Nov 2010 22:08, ludo@gnu.org (Ludovic Courtès) writes:

> Mark H Weaver <mhw@netris.org> writes:
>
>> No, (integer? 3.0) returns #t, as it should, according to R5RS.
>> R5RS's description of "integer?" gives this precise example, and
>> guile's implementation agrees.
>
> Damn, I had never realized that, shame on me.

Bill Schottstaedt has a nice rant on
https://ccrma.stanford.edu/software/snd/snd/s7.html. I think all of his
examples are taken from Guile...

    I can't find the right tone for this section; this is the 400-th
    revision; I wish I were a better writer! I think the exact/inexact
    distinction in Scheme is confused and useless, and leads to verbose
    and buggy code. In some Schemes, "rational" means "could possibly be
    expressed equally well as a ratio (floats are approximations)". In
    s7 it's: "is actually expressed as a ratio (or an integer of
    course)"; otherwise "rational?" is the same as "real?":

    (not-s7-scheme)> (rational? (sqrt 2))
    #t

    As I understand it, "inexact" originally meant "floating point", and
    "exact" meant integer or ratio of integers. But words have a life of
    their own. 0.0 somehow became an "inexact" integer (although it can
    be represented exactly in floating point). +inf.0 must be an integer
    — its fractional part is explicitly zero! But +nan.0... And then
    there's:

    (not-s7-scheme)> (integer? 9007199254740993.1)
    #t

    When does this matter? I often need to index into a vector, but the
    index is a float (a "real" in Scheme-speak: its fractional part can
    be non-zero). In one scheme:

    (not-s7-scheme)> (vector-ref #(0) (floor 0.1))
    ERROR: Wrong type (expecting exact integer): 0.0   ; [why?  "it's probably a programmer mistake"!]

    Not to worry, I'll use inexact->exact:

    (not-s7-scheme)> (inexact->exact 0.1)              ; [why? "floats are ratios"!]
    3602879701896397/36028797018963968

    So I end up using the verbose (floor (inexact->exact ...))
    everywhere, and even then I have no guarantee that I'll get a legal
    vector index. When I started work on s7, I thought perhaps "exact"
    could mean "is represented exactly in the computer". We'd have
    integers and ratios exact; reals and complex exact if they are
    exactly represented in the current floating point
    implementation. 0.0 and 0.5 might be exact if the printout isn't
    misleading, and 0.1 is inexact. "integer?" and friends would refer
    instead to the programmer's point of view. That is, if the
    programmer uses 1 or if the thing prints as 1, it is the integer 1,
    whereas 1.0 means floating point (not integer!). And to keep
    exactness in view, we'd have to monitor which operations introduce
    inexactness — a kind of interval arithmetic. But then what would
    inexact->exact do? If we discard the exact/inexact distinction, we
    can maintain backwards compatibility via:

        (define exact? rational?)
        (define (inexact? x) (not (rational? x)))
        (define inexact->exact rationalize) ; or floor
        (define (exact->inexact x) (* x 1.0))

There is also Mike Sperber's paper a few years ago about Scheme's
numeric tower being borked.

Anyway, just to say that you're in good company :)

Andy
-- 
http://wingolog.org/

Re: fix for expt bug

[Ludovic Courtès]

Hi,

Andy Wingo <wingo@pobox.com> writes:

> There is also Mike Sperber's paper a few years ago about Scheme's
> numeric tower being borked.
>
> Anyway, just to say that you're in good company :)

Heh, good to know, and interesting link!

Ludo’.

[PATCH] Fix bugs in expt and integer-expt

[Mark H Weaver]

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

Here's my patch to fix some bugs in expt and integer-expt.
Please review.

     Mark


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: Patch to fix bugs in expt and integer-expt --]
[-- Type: text/x-diff, Size: 6321 bytes --]

From 19a24107af814e34fe63d02478bfb9441354625e Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
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

Re: [PATCH] Fix bugs in expt and integer-expt

[Ludovic Courtès]

Hello!

Looks good to me, now that I know about inexact integers.  ;-)

Just a couple of details:

Mark H Weaver <mhw@netris.org> writes:

> From 19a24107af814e34fe63d02478bfb9441354625e Mon Sep 17 00:00:00 2001
> From: Mark H Weaver <mhw@netris.org>
> 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.

Please use the C function names for C code, in the log.

Also, can you add a reference to bug #31464 in the log and next to the
corresponding test cases?

Finally, can you document the change for 0^x with x != 0 in NEWS and
explain the departure from R5RS in doc/ref/api-data.texi?

Thanks!

Ludo’.

Re: [PATCH] Fix bugs in expt and integer-expt

[Andy Wingo]

On Fri 05 Nov 2010 03:18, Mark H Weaver <mhw@netris.org> writes:

> Here's my patch to fix some bugs in expt and integer-expt.
> Please review.

Applied with Ludovic's suggestions, thanks to you and Ramakrishnan!

Andy
-- 
http://wingolog.org/