added a test for exact powers for mpfr_root

and fixed mpfr_root for negative x (and odd k)


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@9970 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 23 insertions, 16 deletions

diff --git a/src/root.c b/src/root.c
index 7946c0dce..3d78fe758 100644
--- a/src/root.c
+++ b/src/root.c
@@ -186,22 +186,27 @@ mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode)
 
 /* Compute y <- x^(1/k) using exp(log(x)/k).
    Assume all special cases have been eliminated before.
-   In the extended exponent range, overflows/underflows are not possible. */
+   In the extended exponent range, overflows/underflows are not possible.
+   Assume x > 0, or x < 0 and k odd.
+*/
 static int
 mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode)
 {
-  int inexact, exact_root = 0;
+  int neg, inexact, exact_root = 0;
   mpfr_prec_t w; /* working precision */
-  mpfr_t t;
+  mpfr_t absx, t;
   MPFR_GROUP_DECL(group);
   MPFR_TMP_DECL(marker);
   MPFR_ZIV_DECL(loop);
   MPFR_SAVE_EXPO_DECL (expo);
+  
+  neg = MPFR_IS_NEG (x);
+  MPFR_TMP_INIT_ABS (absx, x);
 
   MPFR_TMP_MARK(marker);
   w = MPFR_PREC(y) + 10;
   /* Take some guard bits to prepare for the 'expt' lost bits below.
-     If x < 2^k, then log(x) < k, thus taking log2(k) bits should be fine. */
+     If |x| < 2^k, then log|x| < k, thus taking log2(k) bits should be fine. */
   if (MPFR_GET_EXP(x) > 0)
     w += MPFR_INT_CEIL_LOG2 (MPFR_GET_EXP(x));
   MPFR_GROUP_INIT_1(group, w, t);
@@ -212,30 +217,30 @@ mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode)
       mpfr_exp_t expt;
       unsigned int err;
 
-      mpfr_log (t, x, MPFR_RNDN);
-      /* t = log(x) * (1 + theta) with |theta| <= 2^(-w) */
+      mpfr_log (t, absx, MPFR_RNDN);
+      /* t = log|x| * (1 + theta) with |theta| <= 2^(-w) */
       mpfr_div_ui (t, t, k, MPFR_RNDN);
       expt = MPFR_GET_EXP (t);
-      /* t = log(x)/k * (1 + theta) + eps with |theta| <= 2^(-w)
+      /* t = log|x|/k * (1 + theta) + eps with |theta| <= 2^(-w)
          and |eps| <= 1/2 ulp(t), thus the total error is bounded
          by 1.5 * 2^(expt - w) */
       mpfr_exp (t, t, MPFR_RNDN);
-      /* t = x^(1/k) * exp(tau) * (1 + theta1) with
+      /* t = |x|^(1/k) * exp(tau) * (1 + theta1) with
          |tau| <= 1.5 * 2^(expt - w) and |theta1| <= 2^(-w).
          For |tau| <= 0.5 we have |exp(tau)-1| < 4/3*tau, thus
          for w >= expt + 2 we have:
-         t = x^(1/k) * (1 + 2^(expt+2)*theta2) * (1 + theta1) with
+         t = |x|^(1/k) * (1 + 2^(expt+2)*theta2) * (1 + theta1) with
          |theta1|, |theta2| <= 2^(-w).
          If expt+2 > 0, as long as w >= 1, we have:
-         t = x^(1/k) * (1 + 2^(expt+3)*theta3) with |theta3| < 2^(-w).
+         t = |x|^(1/k) * (1 + 2^(expt+3)*theta3) with |theta3| < 2^(-w).
          For expt+2 = 0, we have:
-         t = x^(1/k) * (1 + 2^2*theta3) with |theta3| < 2^(-w).
+         t = |x|^(1/k) * (1 + 2^2*theta3) with |theta3| < 2^(-w).
          Finally for expt+2 < 0 we have:
-         t = x^(1/k) * (1 + 2*theta3) with |theta3| < 2^(-w).
+         t = |x|^(1/k) * (1 + 2*theta3) with |theta3| < 2^(-w).
       */
       err = (expt + 2 > 0) ? expt + 3
         : (expt + 2 == 0) ? 2 : 1;
-      /* now t = x^(1/k) * (1 + 2^(err-w)) thus the error is at most
+      /* now t = |x|^(1/k) * (1 + 2^(err-w)) thus the error is at most
          2^(EXP(t) - w + err) */
       if (MPFR_LIKELY (MPFR_CAN_ROUND(t, w - err, MPFR_PREC(y), rnd_mode)))
         break;
@@ -250,9 +255,10 @@ mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode)
         mpfr_init2 (zk, MPFR_PREC(x));
         mpfr_set (z, t, MPFR_RNDN);
         inexact = mpfr_pow_ui (zk, z, k, MPFR_RNDN);
-        exact_root = !inexact && mpfr_equal_p (zk, x);
+        exact_root = !inexact && mpfr_equal_p (zk, absx);
         if (exact_root) /* z is the exact root, thus round z directly */
-          inexact = mpfr_set (y, z, rnd_mode);
+          inexact = (neg == 0) ? mpfr_set (y, z, rnd_mode)
+            : mpfr_neg (y, z, rnd_mode);
         mpfr_clear (zk);
         mpfr_clear (z);
         if (exact_root)
@@ -265,7 +271,8 @@ mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode)
   MPFR_ZIV_FREE (loop);
 
   if (!exact_root)
-    inexact = mpfr_set (y, t, rnd_mode);
+    inexact = (neg == 0) ? mpfr_set (y, t, rnd_mode)
+      : mpfr_neg (y, t, rnd_mode);
 
   MPFR_GROUP_CLEAR(group);
   MPFR_TMP_FREE(marker);