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 <address@hidden> 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 <address@hidden>
> 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