From 6cca8a66a3daedb551f4f80170966d74b6143ba6 Mon Sep 17 00:00:00 2001 From: Ramakrishnan Muthukrishnan 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