Re: fix for expt bug
Hi, Andy Wingo wi...@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’.
Re: fix for expt bug
Hi, On Mon 08 Nov 2010 22:08, l...@gnu.org (Ludovic Courtès) writes: Mark H Weaver m...@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
Hi Mark! Mark H Weaver m...@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
l...@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
Thanks again, Mark and Ludovic.
Attached is an updated patch.
thanks
--
Ramakrishnan
From a1dd2da8562ddeb2052f2994ad0302bcc8d5d1a2 Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan vu3...@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.0002+1.7320508075688772i (expt -8 1/3
;;;
;;; asinh
--
1.7.2.3
Re: fix for expt bug
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 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
On Wed, Nov 3, 2010 at 9:02 PM, Mark H Weaver m...@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
From 6cca8a66a3daedb551f4f80170966d74b6143ba6 Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan vu3...@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.0002+1.7320508075688772i (expt -8 1/3
;;;
;;; asinh
--
1.7.2.3
Re: fix for expt bug
On Wed, Nov 3, 2010 at 9:52 PM, Mark H Weaver m...@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
On Wed, Nov 3, 2010 at 11:23 PM, Ramakrishnan Muthukrishnan vu3...@gmail.com wrote: On Wed, Nov 3, 2010 at 9:52 PM, Mark H Weaver m...@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
Hi Mark,
Mark H Weaver m...@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
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 (x0), 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
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
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
From 025bde78d4c199dee1d2857e913d69ce4d7c2e59 Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan vu3...@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.0002+1.7320508075688772i (expt -8 1/3))
+(eqv-loosely? 1.0002+1.7320508075688772i (expt -8 1/3
;;;
;;; asinh
--
1.7.2.3
Re: fix for expt bug
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
Here is the updated patch.
thanks
Ramakrishnan
From e320d79c8f3cd8b7ddcb9c2d13356e34a3346cfe Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan vu3...@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.0002+1.7320508075688772i (expt -8 1/3))
+(eqv-loosely? 1.0002+1.7320508075688772i (expt -8 1/3
;;;
;;; asinh
--
1.7.2.3
fix for expt bug
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
From c23939a1c2b7bfe1b3cf20abb6a7b431699281a1 Mon Sep 17 00:00:00 2001
From: Ramakrishnan Muthukrishnan vu3...@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.0002+1.7320508075688772i (expt -8 1/3)
+ (eqv-loosely? 1.0002+1.7320508075688772i (expt -8 1/3
;;;
;;; asinh
--
1.7.2.3
Re: fix for expt bug
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 vu3...@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 vu3...@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.0002+1.7320508075688772i (expt -8
1/3)
+ (eqv-loosely? 1.0002+1.7320508075688772i (expt -8
1/3
;;;
;;; asinh
Re: fix for expt bug
On Mon, Nov 1, 2010 at 5:33 AM, Mark H Weaver m...@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
