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 <[email protected]> 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
