Optimization: mpfr_cmp2 now accepts any combination of real arguments

and computes the sign of |b| - |c| (in addition to the number of
cancelled bits); mpfr_add, mpfr_sub, mpfr_sub1, mpfr_agm and some
tests updated to take this change into account.


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

8 files changed, 284 insertions, 238 deletions

diff --git a/add.c b/add.c
index f35f72b18..aa91104e0 100644
--- a/add.c
+++ b/add.c
@@ -87,33 +87,7 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
 
   if (MPFR_SIGN(b) != MPFR_SIGN(c))
     { /* signs differ, it's a subtraction */
-      if (MPFR_EXP(b) < MPFR_EXP(c))
-        {
-          return mpfr_sub1(a, c, b, rnd_mode,
-                           (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b));
-        }
-      else if (MPFR_EXP(b) > MPFR_EXP(c))
-        {
-          return mpfr_sub1(a, b, c, rnd_mode,
-                           (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c));
-        }
-      else
-        { /* MPFR_EXP(b) == MPFR_EXP(c) */
-          int d = mpfr_cmp_abs(b,c);
-          if (d == 0)
-            {
-              if (rnd_mode == GMP_RNDD)
-                MPFR_SET_NEG(a);
-              else
-                MPFR_SET_POS(a);
-              MPFR_SET_ZERO(a);
-              MPFR_RET(0);
-            }
-          else if (d > 0)
-            return mpfr_sub1(a, b, c, rnd_mode, 0);
-          else
-            return mpfr_sub1(a, c, b, rnd_mode, 0);
-        }
+      return mpfr_sub1(a, b, c, rnd_mode, 0);
     }
   else
     { /* signs are equal, it's an addition */
diff --git a/agm.c b/agm.c
index 5ce3c7fb8..7f8a35213 100644
--- a/agm.c
+++ b/agm.c
@@ -162,9 +162,8 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
   /* Main loop */
 
   while (go_on) {
-    int err,  eq, can_round;
-    
-    eq=0;
+    int err, can_round;
+    mp_prec_t eq;
     
     err=1 + (int) ((3.0/2.0*(double)_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
 	     +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
@@ -176,18 +175,16 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
     }
 
     /* Calculus of un and vn */
-    while (eq<=p-2) {
-      mpfr_mul(tmp,u,v,GMP_RNDN);
-      mpfr_sqrt (tmpu, tmp, GMP_RNDN); 
-      mpfr_add(tmp,u,v,GMP_RNDN);
-      mpfr_div_2exp(tmpv,tmp,1,GMP_RNDN);
-      mpfr_set(u,tmpu,GMP_RNDN);
-      mpfr_set(v,tmpv,GMP_RNDN);
-      if (mpfr_cmp(v,u) >= 0)
-	eq = mpfr_cmp2(v,u);
-      else
-	eq = mpfr_cmp2(u,v);
-    }
+    do
+      {
+        mpfr_mul(tmp, u, v, GMP_RNDN);
+        mpfr_sqrt (tmpu, tmp, GMP_RNDN); 
+        mpfr_add(tmp, u, v, GMP_RNDN);
+        mpfr_div_2exp(tmpv, tmp, 1, GMP_RNDN);
+        mpfr_set(u, tmpu, GMP_RNDN);
+        mpfr_set(v, tmpv, GMP_RNDN);
+      }
+    while (mpfr_cmp2(u, v, &eq) != 0 && eq <= p - 2);
 
     /* Roundability of the result */
       can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q);
diff --git a/cmp2.c b/cmp2.c
index 725cf5867..45bbb8a39 100644
--- a/cmp2.c
+++ b/cmp2.c
@@ -25,72 +25,124 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-/* returns the number of cancelled bits when one subtracts abs(c) from abs(b). 
-   Assumes |b| >= |c|, which implies MPFR_EXP(b)>=MPFR_EXP(c).
-   if |b| = |c|, returns prec(b).
+/* Returns the number of cancelled bits when one subtracts |c| from |b|
+   if |b| != |c|, and the sign.
 
    Assumes neither of b or c is NaN or +/- infinity.
 
-   In other terms mpfr_cmp2 (b, c) returns EXP(b) - EXP(b-c).
+   In other terms, if |b| != |c|, mpfr_cmp2 (b, c) returns
+   EXP(max(|b|,|c|)) - EXP(|b| - |c|).
 */
-mp_prec_t
-mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c)
+
+int
+mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mp_prec_t *cancel)
 {
   mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0;
   mp_size_t bn, cn;
   mp_exp_unsigned_t diff_exp;
   mp_prec_t res = 0;
+  int sign;
 
   MPFR_ASSERTN(MPFR_IS_FP(b));
   MPFR_ASSERTN(MPFR_IS_FP(c));
 
-  if (MPFR_IS_ZERO(c))
+  /* Optimized case x - x */
+  if (b == c)
     return 0;
 
-  MPFR_ASSERTN(MPFR_NOTZERO(b));
+  if (MPFR_IS_ZERO(b))
+    {
+      if (MPFR_IS_ZERO(c))
+        return 0;
 
-  bp = MPFR_MANT(b);
-  cp = MPFR_MANT(c);
+      *cancel = 0;
+      return -1;
+    }
 
-  bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
-  cn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
+  if (MPFR_IS_ZERO(c))
+    {
+      *cancel = 0;
+      return 1;
+    }
 
-  MPFR_ASSERTN(MPFR_EXP(b) >= MPFR_EXP(c));
-  diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
-  
-  if (diff_exp == 0) /* otherwise the shifted most significant limb of c
-			cannot match bp[bn] */
+  if (MPFR_EXP(b) >= MPFR_EXP(c))
     {
-      while (bn>=0 && cn>=0 && bp[bn] == cp[cn])
-	{
-	  bn--;
-	  cn--;
-	  res += BITS_PER_MP_LIMB;
-	}
+      sign = 1;
+      diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
 
-      if (bn < 0) /* b = c */
-	return MPFR_PREC(b);
+      bp = MPFR_MANT(b);
+      cp = MPFR_MANT(c);
 
-      if (cn < 0) /* c discards exactly the upper part of b */
-	{
-          unsigned int z;
+      bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
+      cn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
 
-	  while (bn>=0 && bp[bn]==0)
-	    {
-	      bn--;
-	      res += BITS_PER_MP_LIMB;
-	    }
+      if (diff_exp == 0)
+        {
+          while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn])
+            {
+              bn--;
+              cn--;
+              res += BITS_PER_MP_LIMB;
+            }
+
+          if (bn < 0)
+            {
+              if (cn < 0) /* b = c */
+                return 0;
 
-	  if (bn < 0) /* b = c */
-	    return MPFR_PREC(b);
+              bp = cp;
+              bn = cn;
+              cn = -1;
+              sign = -1;
+            }
 
-	  count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
-	  return res + z;
-	}
+          if (cn < 0) /* c discards exactly the upper part of b */
+            {
+              unsigned int z;
+
+              MPFR_ASSERTN(bn >= 0);
+
+              while (bp[bn] == 0)
+                {
+                  if (--bn < 0) /* b = c */
+                    return 0;
+                  res += BITS_PER_MP_LIMB;
+                }
+
+              count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
+              *cancel = res + z;
+              return sign;
+            }
+
+          MPFR_ASSERTN(bn >= 0);
+          MPFR_ASSERTN(cn >= 0);
+          MPFR_ASSERTN(bp[bn] != cp[cn]);
+          if (bp[bn] < cp[cn])
+            {
+              mp_limb_t *tp;
+              mp_size_t tn;
+
+              tp = bp; bp = cp; cp = tp;
+              tn = bn; bn = cn; cn = tn;
+              sign = -1;
+            }
+        }
+    } /* MPFR_EXP(b) >= MPFR_EXP(c) */
+  else /* MPFR_EXP(b) < MPFR_EXP(c) */
+    {
+      sign = -1;
+      diff_exp = (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b);
+
+      bp = MPFR_MANT(c);
+      cp = MPFR_MANT(b);
+
+      bn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
+      cn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
     }
 
   /* now we have removed the identical upper limbs of b and c
-     (can happen only when diff_exp = 0): bp[bn] > cc, bn>=0, cn>=0 */
+     (can happen only when diff_exp = 0), and after the possible
+     swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0 */
 
   if (diff_exp < BITS_PER_MP_LIMB)
     {
@@ -135,7 +187,10 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c)
     {
       res--;
       if (dif != 0)
-	return res;
+        {
+          *cancel = res;
+          return sign;
+        }
     }
   else /* high_dif = 0 */
     {
@@ -144,7 +199,10 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c)
       count_leading_zeros(z, dif); /* dif > 1 here */
       res += z;
       if (dif != (MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - z - 1)))
-	return res; /* dif is not a power of two */
+        { /* dif is not a power of two */
+          *cancel = res;
+          return sign;
+        }
     }
 
   /* now result is res + (low(b) < low(c)) */
@@ -166,18 +224,26 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c)
 	  cn--;
 	}
       if (bp[bn] != cc)
-	return res + (bp[bn] < cc);
+        {
+          *cancel = res + (bp[bn] < cc);
+          return sign;
+        }
       bn--;
     }
 
   if (bn < 0)
     {
       if (lastc != 0)
-	return res + 1;
-      while (cn >= 0 && cp[cn] == 0)
-        cn--;
-      return res + (cn >= 0);
+        res++;
+      else
+        {
+          while (cn >= 0 && cp[cn] == 0)
+            cn--;
+          if (cn >= 0)
+            res++;
+        }
     }
 
-  return res; /* remainder from c is 0 */
+  *cancel = res;
+  return sign;
 }
diff --git a/mpfr-impl.h b/mpfr-impl.h
index b9b96cdb1..fe824c55b 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -146,7 +146,7 @@ void mpfr_restore_emin_emax _PROTO ((void));
 int mpfr_add1 _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
                        mp_rnd_t, mp_exp_unsigned_t));
 int mpfr_sub1 _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
-                       mp_rnd_t, mp_exp_unsigned_t));
+                       mp_rnd_t, int));
 int mpfr_round_raw_generic _PROTO ((mp_limb_t *, mp_limb_t *, mp_prec_t, int,
 				    mp_prec_t, mp_rnd_t, int *, int));
 int mpfr_can_round_raw _PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t,
@@ -154,7 +154,7 @@ int mpfr_can_round_raw _PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t,
 double mpfr_get_d2 _PROTO ((mpfr_srcptr, long)); 
 mp_size_t mpn_sqrtrem_new _PROTO ((mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t));
 int mpfr_cmp_abs _PROTO ((mpfr_srcptr, mpfr_srcptr));
-mp_prec_t mpfr_cmp2 _PROTO ((mpfr_srcptr, mpfr_srcptr));
+int mpfr_cmp2 _PROTO ((mpfr_srcptr, mpfr_srcptr, mp_prec_t *));
 long _mpfr_ceil_log2 _PROTO ((double));
 long _mpfr_floor_log2 _PROTO ((double));
 double _mpfr_ceil_exp2 _PROTO ((double));
diff --git a/sub.c b/sub.c
index e830bf692..95f05c9b6 100644
--- a/sub.c
+++ b/sub.c
@@ -29,124 +29,83 @@ int
 mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
 {
   if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
-  {
-    MPFR_SET_NAN(a);
-    MPFR_RET_NAN;
-  }
+    {
+      MPFR_SET_NAN(a);
+      MPFR_RET_NAN;
+    }
 
   MPFR_CLEAR_NAN(a);
 
   if (MPFR_IS_INF(b))
-  {
-    if (!MPFR_IS_INF(c) || MPFR_SIGN(b) != MPFR_SIGN(c))
-    {
-      MPFR_SET_INF(a);
-      MPFR_SET_SAME_SIGN(a, b);
-      MPFR_RET(0); /* exact */
-    }
-    else
     {
-      MPFR_SET_NAN(a);
-      MPFR_RET_NAN;
+      if (!MPFR_IS_INF(c) || MPFR_SIGN(b) != MPFR_SIGN(c))
+        {
+          MPFR_SET_INF(a);
+          MPFR_SET_SAME_SIGN(a, b);
+          MPFR_RET(0); /* exact */
+        }
+      else
+        {
+          MPFR_SET_NAN(a);
+          MPFR_RET_NAN;
+        }
     }
-  }
   else
     if (MPFR_IS_INF(c))
-    {
-      MPFR_SET_INF(a);
-      if (MPFR_SIGN(c) == MPFR_SIGN(a))
-        MPFR_CHANGE_SIGN(a);
-      MPFR_RET(0); /* exact */
-    }
+      {
+        MPFR_SET_INF(a);
+        if (MPFR_SIGN(c) == MPFR_SIGN(a))
+          MPFR_CHANGE_SIGN(a);
+        MPFR_RET(0); /* exact */
+      }
 
   MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_IS_FP(c));
 
   if (MPFR_IS_ZERO(b))
-  {
-    if (MPFR_IS_ZERO(c))
     {
-      if (MPFR_SIGN(a) !=
-          (rnd_mode != GMP_RNDD ?
-           ((MPFR_SIGN(b) < 0 && MPFR_SIGN(c) > 0) ? -1 : 1) :
-           ((MPFR_SIGN(b) > 0 && MPFR_SIGN(c) < 0) ? 1 : -1)))
-        MPFR_CHANGE_SIGN(a);
-      MPFR_CLEAR_INF(a);
-      MPFR_SET_ZERO(a);
-      MPFR_RET(0); /* 0 - 0 is exact */
+      if (MPFR_IS_ZERO(c))
+        {
+          if (MPFR_SIGN(a) !=
+              (rnd_mode != GMP_RNDD ?
+               ((MPFR_SIGN(b) < 0 && MPFR_SIGN(c) > 0) ? -1 : 1) :
+               ((MPFR_SIGN(b) > 0 && MPFR_SIGN(c) < 0) ? 1 : -1)))
+            MPFR_CHANGE_SIGN(a);
+          MPFR_CLEAR_INF(a);
+          MPFR_SET_ZERO(a);
+          MPFR_RET(0); /* 0 - 0 is exact */
+        }
+      return mpfr_neg (a, c, rnd_mode);
     }
-    return mpfr_neg (a, c, rnd_mode);
-  }
 
   if (MPFR_IS_ZERO(c))
-  {
-    return mpfr_set (a, b, rnd_mode);
-  }
+    {
+      return mpfr_set (a, b, rnd_mode);
+    }
 
   MPFR_CLEAR_INF(a);
 
   if (MPFR_SIGN(b) == MPFR_SIGN(c))
-  { /* signs are equal, it's a real subtraction */
-    if (MPFR_EXP(b) < MPFR_EXP(c))
-    { /* exchange rounding modes towards +/- infinity */
-      int inexact;
-      if (rnd_mode == GMP_RNDU)
-        rnd_mode = GMP_RNDD;
-      else if (rnd_mode == GMP_RNDD)
-        rnd_mode = GMP_RNDU;
-      inexact = - mpfr_sub1(a, c, b, rnd_mode,
-                            (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b));
-      MPFR_CHANGE_SIGN(a);
-      return inexact;
-    }
-    else if (MPFR_EXP(b) > MPFR_EXP(c))
-      return mpfr_sub1(a, b, c, rnd_mode,
-                       (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c));
-    else
-    { /* MPFR_EXP(b) == MPFR_EXP(c) */
-      int d = mpfr_cmp_abs (b, c);
-
-      if (d == 0)
-      {
-        if (rnd_mode == GMP_RNDD)
-          MPFR_SET_NEG(a);
-        else
-          MPFR_SET_POS(a);
-        MPFR_SET_ZERO(a);
-        MPFR_RET(0);
-      }
-      else if (d > 0)
-        return mpfr_sub1 (a, b, c, rnd_mode, 0);
-      else
-      { /* exchange rounding modes towards +/- infinity */
-        int inexact;
-        if (rnd_mode == GMP_RNDU)
-          rnd_mode = GMP_RNDD;
-        else if (rnd_mode == GMP_RNDD)
-          rnd_mode = GMP_RNDU;
-	inexact = - mpfr_sub1 (a, c, b, rnd_mode, 0);
-	MPFR_CHANGE_SIGN(a);
-        return inexact;
-      }
+    { /* signs are equal, it's a real subtraction */
+      return mpfr_sub1(a, b, c, rnd_mode, 1);
     }
-  }
   else
-  { /* signs differ, it's an addition */
-    if (MPFR_EXP(b) < MPFR_EXP(c))
-    { /* exchange rounding modes towards +/- infinity */
-      int inexact;
-      if (rnd_mode == GMP_RNDU)
-        rnd_mode = GMP_RNDD;
-      else if (rnd_mode == GMP_RNDD)
-        rnd_mode = GMP_RNDU;
-      inexact = mpfr_add1(a, c, b, rnd_mode,
-                          (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b));
-      MPFR_CHANGE_SIGN(a);
-      return inexact;
-    }
-    else
-    {
-      return mpfr_add1(a, b, c, rnd_mode,
-                       (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c));
+    { /* signs differ, it's an addition */
+      if (MPFR_EXP(b) < MPFR_EXP(c))
+        { /* exchange rounding modes towards +/- infinity */
+          int inexact;
+          if (rnd_mode == GMP_RNDU)
+            rnd_mode = GMP_RNDD;
+          else if (rnd_mode == GMP_RNDD)
+            rnd_mode = GMP_RNDU;
+          inexact = mpfr_add1(a, c, b, rnd_mode,
+                              (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b));
+          MPFR_CHANGE_SIGN(a);
+          return inexact;
+        }
+      else
+        {
+          return mpfr_add1(a, b, c, rnd_mode,
+                           (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c));
+        }
     }
-  }
 }
diff --git a/sub1.c b/sub1.c
index d381ee52d..99e119d8e 100644
--- a/sub1.c
+++ b/sub1.c
@@ -32,10 +32,12 @@ MA 02111-1307, USA. */
 */
 
 int
-mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
-	   mp_rnd_t rnd_mode, mp_exp_unsigned_t diff_exp)
+mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode,
+           int sub)
 {
-  unsigned long cancel, cancel1;
+  int sign;
+  mp_exp_unsigned_t diff_exp;
+  mp_prec_t cancel, cancel1;
   mp_size_t cancel2, an, bn, cn, cn0;
   mp_limb_t *ap, *bp, *cp;
   mp_limb_t carry, bb, cc, borrow = 0;
@@ -43,18 +45,36 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
   int sh, k;
   TMP_DECL(marker);
 
-#ifdef DEBUG
-  printf("\nenter mpfr_sub, rnd_mode=%s:\n", mpfr_print_rnd_mode(rnd_mode));
-  printf("b="); if (MPFR_SIGN(b)>0) putchar(' '); mpfr_print_binary(b); putchar('\n');
-  printf("c="); if (MPFR_SIGN(c)>0) putchar(' '); for (k=0; k<diff_exp; k++) putchar(' '); mpfr_print_binary(c); putchar('\n');
-  printf("PREC(a)=%u PREC(b)=%u PREC(c)=%u\n", MPFR_PREC(a), MPFR_PREC(b),
-	 MPFR_PREC(c));
-#endif
   TMP_MARK(marker);
   ap = MPFR_MANT(a);
   an = 1 + (MPFR_PREC(a) - 1) / BITS_PER_MP_LIMB;
 
-  cancel = mpfr_cmp2 (b, c);
+  sign = mpfr_cmp2 (b, c, &cancel);
+  if (sign == 0)
+    {
+      if (rnd_mode == GMP_RNDD)
+        MPFR_SET_NEG(a);
+      else
+        MPFR_SET_POS(a);
+      MPFR_SET_ZERO(a);
+      MPFR_RET(0);
+    }
+
+  /* If subtraction: sign(a) = sign * sign(b) */
+  if (sub && MPFR_SIGN(a) != sign * MPFR_SIGN(b))
+    MPFR_CHANGE_SIGN(a);
+
+  if (sign < 0) /* swap b and c so that |b| > |c| */
+    {
+      mpfr_srcptr t;
+      t = b; b = c; c = t;
+    }
+
+  /* If addition: sign(a) = sign of the larger argument in absolute value */
+  if (!sub)
+    MPFR_SET_SAME_SIGN(a, b);
+
+  diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
 
   /* reserve a space to store b aligned with the result, i.e. shifted by
      (-cancel) % BITS_PER_MP_LIMB to the right */
@@ -112,10 +132,6 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
   printf("cancel=%u cancel1=%u cancel2=%d\n", cancel, cancel1, cancel2);
 #endif
 
-  /* adjust sign of result */
-  if (MPFR_SIGN(a)*MPFR_SIGN(b) < 0)
-    MPFR_CHANGE_SIGN(a);
-
   /*               ap[an-1]        ap[0]
              <----------------+-----------|---->
              <----------PREC(a)----------><-sh->
@@ -214,8 +230,8 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
     }
   else /* directed rounding: set rnd_mode to RNDZ iff towards zero */
     {
-      if (((rnd_mode == GMP_RNDD) && (MPFR_SIGN(b) > 0)) ||
-	  ((rnd_mode == GMP_RNDU) && (MPFR_SIGN(b) < 0)))
+      if (((rnd_mode == GMP_RNDD) && (MPFR_SIGN(a) > 0)) ||
+	  ((rnd_mode == GMP_RNDU) && (MPFR_SIGN(a) < 0)))
 	rnd_mode = GMP_RNDZ;
 
       if (carry)
@@ -403,5 +419,5 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
 #endif
   /* check that result is msb-normalized */
   MPFR_ASSERTN(ap[an-1] > ~ap[an-1]);
-  return inexact * MPFR_SIGN(b);
+  return inexact * MPFR_SIGN(a);
 }
diff --git a/tests/tcmp.c b/tests/tcmp.c
index 98b7b4412..df6f23bdd 100644
--- a/tests/tcmp.c
+++ b/tests/tcmp.c
@@ -36,6 +36,7 @@ main (void)
   double x, y;
   mpfr_t xx, yy;
   int i, c;
+  mp_prec_t p;
 
   mpfr_init (xx);
   mpfr_init (yy);
@@ -54,18 +55,30 @@ main (void)
   mpfr_set_prec (yy, 65);
   mpfr_set_str_raw(xx, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
   mpfr_set_str_raw(yy, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
-  if (mpfr_cmp2(xx,yy)!=64) { printf("Error (1) in mpfr_cmp\n"); exit(1); }
+  p = 0;
+  if (mpfr_cmp2(xx, yy, &p) <= 0 || p != 64)
+    {
+      printf("Error (1) in mpfr_cmp2\n");
+      exit(1);
+    }
   mpfr_set_str_raw(xx, "0.10100010001110110111000010001000010011111101000100011101000011100");
   mpfr_set_str_raw(yy, "0.10100010001110110111000010001000010011111101000100011101000011011");
-  if (mpfr_cmp2(xx,yy)!=64) { printf("Error (2) in mpfr_cmp\n"); exit(1); }
+  p = 0;
+  if (mpfr_cmp2(xx, yy, &p) <= 0 || p != 64)
+    {
+      printf("Error (2) in mpfr_cmp2\n");
+      exit(1);
+    }
 
   mpfr_set_prec (xx, 160); mpfr_set_prec (yy, 160);
   mpfr_set_str_raw (xx, "0.1E1");
   mpfr_set_str_raw (yy, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100");
-  if (mpfr_cmp2 (xx, yy) != 144) {
-    printf("Error (3) in mpfr_cmp\n");
-    exit(1);
-  }
+  p = 0;
+  if (mpfr_cmp2 (xx, yy, &p) <= 0 || p != 144)
+    {
+      printf("Error (3) in mpfr_cmp2\n");
+      exit(1);
+    }
 
   mpfr_set_prec(xx, 53); mpfr_set_prec(yy, 200);
   mpfr_set_d(xx, 1.0, 0);
diff --git a/tests/tcmp2.c b/tests/tcmp2.c
index 61bfaaa1c..70349fd40 100644
--- a/tests/tcmp2.c
+++ b/tests/tcmp2.c
@@ -57,7 +57,8 @@ void
 worst_cases ()
 {
   mpfr_t x, y;
-  unsigned int i, j, k, l, b, expected;
+  unsigned int i, j, k, b, expected;
+  mp_prec_t l;
 
   mpfr_init2 (x, 200);
   mpfr_init2 (y, 200);
@@ -68,24 +69,26 @@ worst_cases ()
       mpfr_set_ui (x, 1, GMP_RNDN);
       mpfr_div_2exp (y, y, 1, GMP_RNDN); /* y = 1/2^i */
 
-      if ((l = mpfr_cmp2 (x, y)) != 1)
+      l = 0;
+      if (mpfr_cmp2 (x, y, &l) <= 0 || l != 1)
 	{
 	  fprintf (stderr, "Error in mpfr_cmp2:\nx=");
 	  mpfr_out_str (stderr, 2, 0, x, GMP_RNDN);
 	  fprintf (stderr, "\ny=");
 	  mpfr_out_str (stderr, 2, 0, y, GMP_RNDN);
-	  fprintf (stderr, "\ngot %u instead of %u\n", l, 1);
+	  fprintf (stderr, "\ngot %lu instead of %u\n", l, 1);
 	  exit(1);
 	}
 
       mpfr_add (x, x, y, GMP_RNDN); /* x = 1 + 1/2^i */
-      if ((l = mpfr_cmp2 (x, y)) != 0)
+      l = 0;
+      if (mpfr_cmp2 (x, y, &l) <= 0 || l != 0)
 	{
 	  fprintf (stderr, "Error in mpfr_cmp2:\nx=");
 	  mpfr_out_str (stderr, 2, 0, x, GMP_RNDN);
 	  fprintf (stderr, "\ny=");
 	  mpfr_out_str (stderr, 2, 0, y, GMP_RNDN);
-	  fprintf (stderr, "\ngot %u instead of %u\n", l, 0);
+	  fprintf (stderr, "\ngot %lu instead of %u\n", l, 0);
 	  exit(1);
 	}
     }
@@ -110,7 +113,7 @@ worst_cases ()
 
 	      set_bit (x, i + j + k + 2, 1);
 	      set_bit (y, i + j + k + 2, 0);
-	      l = mpfr_cmp2 (x, y);
+	      l = 0; mpfr_cmp2 (x, y, &l);
 	      expected = i + j + k + 1;
 	      if (l != expected)
 		{
@@ -118,14 +121,14 @@ worst_cases ()
 		  mpfr_out_str (stderr, 2, 0, x, GMP_RNDN);
 		  fprintf (stderr, "\ny=");
 		  mpfr_out_str (stderr, 2, 0, y, GMP_RNDN);
-		  fprintf (stderr, "\ngot %u instead of %u\n", l, expected);
+		  fprintf (stderr, "\ngot %lu instead of %u\n", l, expected);
 		  exit(1);
 		}
 
 	      set_bit (x, i + j + k + 2, 0);
 	      set_bit (x, i + j + k + 3, 0);
 	      set_bit (y, i + j + k + 3, 1);
-	      l = mpfr_cmp2 (x, y);
+	      l = 0; mpfr_cmp2 (x, y, &l);
 	      expected = i + j + k + 2;
 	      if (l != expected)
 		{
@@ -133,7 +136,7 @@ worst_cases ()
 		  mpfr_out_str (stderr, 2, 0, x, GMP_RNDN);
 		  fprintf (stderr, "\ny=");
 		  mpfr_out_str (stderr, 2, 0, y, GMP_RNDN);
-		  fprintf (stderr, "\ngot %u instead of %u\n", l, expected);
+		  fprintf (stderr, "\ngot %lu instead of %u\n", l, expected);
 		  exit(1);
 		}
 	    }
@@ -148,7 +151,7 @@ void
 tcmp2 (double x, double y, int i)
 {
   mpfr_t xx, yy;
-  int j;
+  mp_prec_t j;
 
   if (i==-1) {
     if (x==y) i=53;
@@ -157,13 +160,25 @@ tcmp2 (double x, double y, int i)
   mpfr_init2(xx, 53); mpfr_init2(yy, 53);
   mpfr_set_d (xx, x, GMP_RNDN);
   mpfr_set_d (yy, y, GMP_RNDN);
-  j = mpfr_cmp2 (xx, yy);
-  if (j != i) {
+  j = 0;
+  if (mpfr_cmp2 (xx, yy, &j) == 0)
+    {
+      if (x != y)
+        {
+          fprintf (stderr, "Error in mpfr_cmp2 for\nx=");
+          mpfr_out_str (stderr, 2, 0, xx, GMP_RNDN);
+          fprintf (stderr, "\ny=");
+          mpfr_out_str (stderr, 2, 0, yy, GMP_RNDN);
+          fprintf (stderr, "\ngot sign 0 for x != y\n");
+          exit(1);
+        }
+    }
+  else if (j != i) {
     fprintf (stderr, "Error in mpfr_cmp2 for\nx=");
     mpfr_out_str (stderr, 2, 0, xx, GMP_RNDN);
     fprintf (stderr, "\ny=");
     mpfr_out_str (stderr, 2, 0, yy, GMP_RNDN);
-    fprintf (stderr, "\ngot %u instead of %u\n", j, i);
+    fprintf (stderr, "\ngot %lu instead of %u\n", j, i);
     exit(1);
   }
   mpfr_clear(xx); mpfr_clear(yy);
@@ -172,7 +187,7 @@ tcmp2 (double x, double y, int i)
 void special ()
 {
   mpfr_t x, y;
-  int j;
+  mp_prec_t j;
 
   mpfr_init (x); mpfr_init (y);
 
@@ -181,7 +196,8 @@ void special ()
   mpfr_set_prec (y, 65);
   mpfr_set_str_raw (x, "0.10000000000000000000000000000000000001110010010110100110011110000E1");
   mpfr_set_str_raw (y, "0.11100100101101001100111011111111110001101001000011101001001010010E-35");
-  if ((j = mpfr_cmp2 (x, y)) != 1) {
+  j = 0;
+  if (mpfr_cmp2 (x, y, &j) <= 0 || j != 1) {
     printf ("Error in mpfr_cmp2:\n");
     printf ("x=");
     mpfr_print_binary (x);
@@ -189,29 +205,31 @@ void special ()
     printf ("y=");
     mpfr_print_binary (y);
     putchar ('\n');
-    printf ("got %d, expected 1\n", j);
+    printf ("got %lu, expected 1\n", j);
     exit (1);
   }
 
   mpfr_set_prec(x, 127); mpfr_set_prec(y, 127);
   mpfr_set_str_raw(x, "0.1011010000110111111000000101011110110001000101101011011110010010011110010000101101000010011001100110010000000010110000101000101E6");
   mpfr_set_str_raw(y, "0.1011010000110111111000000101011011111100011101000011001111000010100010100110110100110010011001100110010000110010010110000010110E6");
-  if ((j=mpfr_cmp2(x, y)) != 32) {
+  j = 0;
+  if (mpfr_cmp2(x, y, &j) <= 0 || j != 32) {
     printf("Error in mpfr_cmp2:\n");
     printf("x="); mpfr_print_binary(x); putchar('\n');
     printf("y="); mpfr_print_binary(y); putchar('\n');
-    printf("got %d, expected 32\n", j);
+    printf("got %lu, expected 32\n", j);
     exit(1);
   }
 
   mpfr_set_prec (x, 128); mpfr_set_prec (y, 239);
   mpfr_set_str_raw (x, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100100000000000E167");
   mpfr_set_str_raw (y, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100011111111111111111111111111111111111111111111111011111100101011100011001101000100111000010000000000101100110000111111000101E167");
-  if ((j=mpfr_cmp2(x, y)) != 164) {
+  j = 0;
+  if (mpfr_cmp2(x, y, &j) <= 0 || j != 164) {
     printf("Error in mpfr_cmp2:\n");
     printf("x="); mpfr_print_binary(x); putchar('\n');
     printf("y="); mpfr_print_binary(y); putchar('\n');
-    printf("got %d, expected 164\n", j);
+    printf("got %lu, expected 164\n", j);
     exit(1);
   }
 
@@ -219,11 +237,12 @@ void special ()
   mpfr_set_prec (x, 130); mpfr_set_prec (y, 130);
   mpfr_set_str_raw (x, "0.1100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2");
   mpfr_set_str_raw (y, "0.1011111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100E2");
-  if ((j=mpfr_cmp2(x, y)) != 127) {
+  j = 0;
+  if (mpfr_cmp2(x, y, &j) <= 0 || j != 127) {
     printf("Error in mpfr_cmp2:\n");
     printf("x="); mpfr_print_binary(x); putchar('\n');
     printf("y="); mpfr_print_binary(y); putchar('\n');
-    printf("got %d, expected 127\n", j);
+    printf("got %lu, expected 127\n", j);
     exit(1);
   }
 
@@ -231,11 +250,12 @@ void special ()
   mpfr_set_prec (x, 65); mpfr_set_prec (y, 65);
   mpfr_set_ui (x, 5, GMP_RNDN);
   mpfr_set_str_raw (y, "0.10011111111111111111111111111111111111111111111111111111111111101E3");
-  if ((j=mpfr_cmp2(x, y)) != 63) {
+  j = 0;
+  if (mpfr_cmp2(x, y, &j) <= 0 || j != 63) {
     printf("Error in mpfr_cmp2:\n");
     printf("x="); mpfr_print_binary(x); putchar('\n');
     printf("y="); mpfr_print_binary(y); putchar('\n');
-    printf("got %d, expected 63\n", j);
+    printf("got %lu, expected 63\n", j);
     exit(1);
   }
 
@@ -243,11 +263,12 @@ void special ()
   mpfr_set_prec (x, 65); mpfr_set_prec (y, 65);
   mpfr_set_str_raw (x, "0.10011011111000101001110000000000000000000000000000000000000000000E-69");
   mpfr_set_str_raw (y, "0.10011011111000101001101111111111111111111111111111111111111111101E-69");
-  if ((j=mpfr_cmp2(x, y)) != 63) {
+  j = 0;
+  if (mpfr_cmp2(x, y, &j) <= 0 || j != 63) {
     printf("Error in mpfr_cmp2:\n");
     printf("x="); mpfr_print_binary(x); putchar('\n');
     printf("y="); mpfr_print_binary(y); putchar('\n');
-    printf("got %d, expected 63\n", j);
+    printf("got %lu, expected 63\n", j);
     exit(1);
   }