inexact flag was wrong in some cases

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

1 files changed, 54 insertions, 37 deletions

diff --git a/sqrt.c b/sqrt.c
index b92779f67..99d8b11e7 100644
--- a/sqrt.c
+++ b/sqrt.c
@@ -1,6 +1,7 @@
 /* mpfr_sqrt -- square root of a floating-point number
 
-Copyright (C) 1999 Free Software Foundation.
+Copyright (C) 1999-2001 Free Software Foundation.
+Contributed by the Spaces project.
 
 This file is part of the MPFR Library.
 
@@ -44,7 +45,7 @@ mpfr_sqrt (r, u, rnd_mode)
   mp_size_t prec, err;
   mp_limb_t q_limb;
   long rw, nw, k; 
-  int exact = 0, t;
+  int inexact = 0, t;
   unsigned long cc = 0; 
   char can_round = 0; 
   TMP_DECL(marker0);
@@ -53,18 +54,18 @@ mpfr_sqrt (r, u, rnd_mode)
 
   if (MPFR_IS_NAN(u)) {
     MPFR_SET_NAN(r);
-    return 1;
+    return 1; /* NaN is always inexact */
   }
 
   if (MPFR_SIGN(u) < 0) {
     if (MPFR_IS_INF(u) || MPFR_NOTZERO(u)) {
       MPFR_SET_NAN(r);
-      return 1;
+      return 1; /* NaN is always inexact */
     }
     else { /* sqrt(-0) = -0 */
       MPFR_SET_ZERO(r);
       if (MPFR_SIGN(r) > 0) MPFR_CHANGE_SIGN(r);
-      return 0;
+      return 0; /* zero is exact */
     }
   }
 
@@ -74,7 +75,7 @@ mpfr_sqrt (r, u, rnd_mode)
   if (MPFR_IS_INF(u)) 
     { 
       MPFR_SET_INF(r);
-      return 1;
+      return 0; /* infinity is exact */
     }
 
   MPFR_CLEAR_INF(r);
@@ -86,11 +87,10 @@ mpfr_sqrt (r, u, rnd_mode)
       MPFR_EXP(r) = 0;
       rsize = (prec-1)/BITS_PER_MP_LIMB + 1;
       MPN_ZERO(MPFR_MANT(r), rsize);
-      return 0;
+      return 0; /* zero is exact */
     }
 
   up = MPFR_MANT(u);
-
   usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; 
 
 #ifdef DEBUG
@@ -164,14 +164,15 @@ mpfr_sqrt (r, u, rnd_mode)
       q_limb = mpn_sqrtrem_new (rp, remp, tmp, rsize);
 
 #ifdef DEBUG
-      printf("The result is : \n"); 
-      printf("sqrt : "); 
-      for(k = rrsize - 1; k >= 0; k--) { printf("%lu ", rp[k]); }
-      printf("(q_limb = %lu)\n", q_limb); 
+      printf ("The result is : \n"); 
+      printf ("sqrt : "); 
+      for (k = rrsize - 1; k >= 0; k--)
+	printf ("%lu ", rp[k]);
+      printf ("(inexact = %lu)\n", q_limb); 
 #endif
       
-      can_round = (mpfr_can_round_raw(rp, rrsize, 1, err, 
-				      GMP_RNDZ, rnd_mode, MPFR_PREC(r))); 
+      can_round = mpfr_can_round_raw(rp, rrsize, 1, err, 
+				     GMP_RNDZ, rnd_mode, MPFR_PREC(r));
 
       /* If we used all the limbs of both the dividend and the divisor, 
 	 then we have the correct RNDZ rounding */
@@ -186,60 +187,76 @@ mpfr_sqrt (r, u, rnd_mode)
     }
   while (!can_round && (rsize < 2*usize) 
 	 && (rsize += 2) && (rrsize ++)); 
-
+#ifdef DEBUG
+  printf ("can_round = %d\n", can_round);
+#endif
 
   /* This part may be deplaced upper to avoid a few mpfr_can_round_raw */
   /* when the square root is exact. It is however very unprobable that */
   /* it would improve the behaviour of the present code on average.    */
 
-  if (!q_limb) /* possibly exact */
+  if (!q_limb) /* the sqrtrem call was exact, possible exact square root */
     {
       /* if we have taken into account the whole of up */
       for (k = usize - rsize - 1; k >= 0; k ++)
 	if (up[k]) break; 
       
-      if (k < 0) { exact = 1; goto fin; }
+      if (k < 0)
+	goto fin; /* exact square root ==> inexact = 0 */
     }
 
   if (can_round) 
     {
-      cc = mpfr_round_raw(rp, rp, err, 0, MPFR_PREC(r), rnd_mode, NULL);
+      cc = mpfr_round_raw (rp, rp, err, 0, MPFR_PREC(r), rnd_mode, &inexact);
+      if (!inexact) /* exact high part: inexact flag depends from remainder */
+	inexact = -q_limb;
       rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1; 
     }
   else
     /* Use the return value of sqrtrem to decide of the rounding         */
     /* Note that at this point the sqrt has been computed                */
-    /* EXACTLY. If rounding = GMP_RNDZ, do nothing [comes from           */
-    /* the fact that the exact square root can end with a bunch of ones, */
-    /* and in that case we indeed cannot round if we do not know that    */
-    /* the computation was exact.                                        */
+    /* EXACTLY.                                                          */
     switch (rnd_mode)
       {
       case GMP_RNDZ : 
-      case GMP_RNDD : break; 
+      case GMP_RNDD :
+	inexact = -1; /* result is truncated */
+	break; 
 
       case GMP_RNDN : 
 	/* Not in the situation ...0 111111 */
 	rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1);
-	if (rw) { rw = BITS_PER_MP_LIMB - rw; nw = 0; } else nw = 1; 
+	if (rw)
+	  {
+	    rw = BITS_PER_MP_LIMB - rw;
+	    nw = 0;
+	  }
+	else
+	  nw = 1; 
 	if ((rp[nw] >> rw) & 1 &&                     /* Not 0111111111 */
 	    (q_limb ||                                /* Nonzero remainder */
 	    (rw ? (rp[nw] >> (rw + 1)) & 1 : 
 	     (rp[nw] >> (BITS_PER_MP_LIMB - 1)) & 1))) /* or even rounding */ 
-	  cc = mpn_add_1(rp + nw, rp + nw, rrsize, ((mp_limb_t)1) << rw); 
+	  {
+	    cc = mpn_add_1 (rp + nw, rp + nw, rrsize, ((mp_limb_t)1) << rw);
+	    inexact = 1;
+	  }
+	else
+	  inexact = -1;
 	break;
  
-      case GMP_RNDU : 
-	if (q_limb)
-	  {
-	    t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1); 
-	    if (t) 
-	      {
-		cc = mpn_add_1(rp, rp, rrsize, 1 << (BITS_PER_MP_LIMB - t)); 
-	      }
-	    else
-	      cc = mpn_add_1 (rp, rp, rrsize, 1); 
-	  }	      
+      case GMP_RNDU:
+	/* we should arrive here only when the result is inexact,
+	   i.e. either q_limb > 0 (the remainder from mpn_sqrtrem is non-zero)
+	   or up[0..usize-rsize-1] is non zero, thus we have to add one
+	   ulp, and inexact = 1 */
+	inexact = 1;
+	t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1); 
+	rsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
+	if (t) 
+	    cc = mpn_add_1 (rp + rrsize - rsize, rp + rrsize - rsize, rsize, 1 << (BITS_PER_MP_LIMB - t));
+	else
+	    cc = mpn_add_1 (rp + rrsize - rsize, rp + rrsize - rsize, rsize, 1);
       }
 
   if (cc) {
@@ -260,5 +277,5 @@ mpfr_sqrt (r, u, rnd_mode)
   TMP_FREE (marker);
   }
   TMP_FREE(marker0);
-  return exact;
+  return inexact;
 }