Bug fixed and update to check the exponents.

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

1 files changed, 289 insertions, 181 deletions

diff --git a/add.c b/add.c
index 082e8e1fd..37b3eac85 100644
--- a/add.c
+++ b/add.c
@@ -25,184 +25,263 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-/* #define DEBUG */
-
-extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, 
-			      mp_rnd_t, int));
-void mpfr_add1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t, int));
-
 #define ONE ((mp_limb_t) 1)
 
+extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
+                              mp_rnd_t, mp_exp_unsigned_t));
+void mpfr_add1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
+                       mp_rnd_t, mp_exp_unsigned_t));
+
 /* signs of b and c are supposed equal,
    diff_exp is the difference between the exponents of b and c,
    which is supposed >= 0 */
 
-void 
+void
 #if __STDC__
-mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, 
-	   mp_rnd_t rnd_mode, int diff_exp) 
+mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
+	   mp_rnd_t rnd_mode, mp_exp_unsigned_t diff_exp)
 #else
-mpfr_add1 (a, b, c, rnd_mode, diff_exp) 
+mpfr_add1 (a, b, c, rnd_mode, diff_exp)
      mpfr_ptr a;
      mpfr_srcptr b;
      mpfr_srcptr c;
      mp_rnd_t rnd_mode;
-     int diff_exp;
+     mp_exp_unsigned_t diff_exp;
 #endif
 {
-  mp_limb_t *ap, *bp, *cp, cc, c2=0, c3=0;
-  unsigned int an, bn, cn;
-  int sh, dif, k;
-  TMP_DECL(marker); 
-
-  TMP_MARK(marker); 
+  mp_limb_t *ap, *bp, *cp;
+  mp_limb_t cc, c2 = 0, c3 = 0;  /* should be removed later */
+  mp_exp_t dif;  /* should be removed later */
+  mp_prec_t aq, bq, cq;
+  mp_size_t an, bn, cn, k;
+  int sh;
+  int nulp = 1;
+  TMP_DECL(marker);
+
+  MPFR_ASSERTN(MPFR_NOTZERO(c));
+
+  TMP_MARK(marker);
   ap = MPFR_MANT(a);
   bp = MPFR_MANT(b);
   cp = MPFR_MANT(c);
-  if (ap == bp) {
-    bp = (mp_ptr) TMP_ALLOC(MPFR_ABSSIZE(b) * BYTES_PER_MP_LIMB); 
+
+  if (ap == bp)
+  {
+    bp = (mp_ptr) TMP_ALLOC(MPFR_ABSSIZE(b) * BYTES_PER_MP_LIMB);
     MPN_COPY (bp, ap, MPFR_ABSSIZE(b));
     if (ap == cp) { cp = bp; }
   }
   else if (ap == cp)
-    {
-      cp = (mp_ptr) TMP_ALLOC (MPFR_ABSSIZE(c) * BYTES_PER_MP_LIMB);
-      MPN_COPY(cp, ap, MPFR_ABSSIZE(c)); 
-    }
+  {
+    cp = (mp_ptr) TMP_ALLOC (MPFR_ABSSIZE(c) * BYTES_PER_MP_LIMB);
+    MPN_COPY(cp, ap, MPFR_ABSSIZE(c));
+  }
 
-  an = (MPFR_PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */
+  aq = MPFR_PREC(a);
+  bq = MPFR_PREC(b);
+  cq = MPFR_PREC(c);
 
-  sh = an*BITS_PER_MP_LIMB-MPFR_PREC(a); /* non-significant bits in low limb */
-  bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */
-  cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB + 1;
-  MPFR_EXP(a) = MPFR_EXP(b);
+  an = (aq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of a */
+  sh = an*BITS_PER_MP_LIMB - aq;    /* non-significant bits in low limb */
+  bn = (bq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */
+  cn = (cq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of c */
 
-  if (MPFR_SIGN(a) * MPFR_SIGN(b) < 0) MPFR_CHANGE_SIGN(a);
+  MPFR_EXP(a) = MPFR_EXP(b);
+  MPFR_SET_SAME_SIGN(a, b);
 
-  /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit
-     through rounding */
-  dif = MPFR_PREC(a) - diff_exp;
+  if (aq <= diff_exp)
+  {
 
-#ifdef DEBUG
-  printf("diff_exp=%u dif=MPFR_PREC(a)-diff_exp=%d\n", diff_exp, dif);
-  printf("b= "); mpfr_print_raw(b); putchar('\n');
-  printf("c="); for (k=0;k<diff_exp;k++) putchar(' ');
-  if (MPFR_SIGN(c)>0) putchar(' '); mpfr_print_raw(c); putchar('\n');
-#endif
-  if (dif <= 0)
-    { 
     /* diff_exp>=MPFR_PREC(a): c does not overlap with a */
-    /* either MPFR_PREC(b)<=MPFR_PREC(a), and we can copy the mantissa of b directly 
-       into that of a, or MPFR_PREC(b)>MPFR_PREC(a) and we have to round b+c */
+    /* either MPFR_PREC(b)<=MPFR_PREC(a), and we can copy the mantissa of b
+       directly into that of a, or MPFR_PREC(b)>MPFR_PREC(a) and we have to
+       round b+c */
 
-    if (MPFR_PREC(b) <= MPFR_PREC(a)) {
+    if (bq <= aq)
+    {
 
-      MPN_COPY(ap+(an-bn), bp, bn);
-      /* fill low significant limbs with zero */
+      MPN_COPY(ap + (an - bn), bp, bn);
 
-      for (bp=ap;bn<an;bn++) *bp++=0;
+      /* fill low significant limbs with zero */
+      MPN_ZERO(ap, an - bn);
 
-      /* now take c into account */
-      if (rnd_mode==GMP_RNDN) { 
-	
+      /* now take c into account (c != 0) */
+      if (rnd_mode == GMP_RNDN)
+      {
 	/* to nearest */
 	/* if diff_exp > MPFR_PREC(a), no change */
 
-	if (diff_exp==MPFR_PREC(a)) {
+	if (aq == diff_exp)
+	{
+	  /* as c is normalized, we have to add one to the lsb of a
+             if c>1/2, or c=1/2 and lsb(a)=1 (round to even) */
+
+          mp_limb_t cc;
+          mp_limb_t *cp2 = cp + (cn-1);  /* highest limb */
+
+          cc = *cp2 - MP_LIMB_T_HIGHBIT;
+          while (cc == 0 && cp2 > cp) cc = *--cp2;
+
+          if (cc || (ap[0] & (ONE<<sh))) goto add_one_ulp;
+          /* mant(c) != 1/2 or mant(c) = 1/2: add 1 iff lsb(a)=1 */
+        } /* aq == diff_exp */
+      } /* rnd_mode == GMP_RNDN */
+      else if ((rnd_mode == GMP_RNDU && MPFR_ISNONNEG(b)) ||
+               (rnd_mode == GMP_RNDD && MPFR_ISNEG(b)))  /* round up */
+        goto add_one_ulp;
+      /* in the other cases (round to zero, or up/down with sign -/+),
+         nothing to do */
+    } /* bq <= aq */
+    else /* bq > aq */
+    {
+      mp_limb_t inv, bb, cc = 0;
+      int difs, r;
+      mp_exp_t difw;  /* mp_exp_t may be larger than mp_size_t */
 
-	  /* if c is not zero, then as it is normalized, we have to add
-	     one to the lsb of a if c>1/2, or c=1/2 and lsb(a)=1 (round to
-	     even) */
-	  
-	  if (MPFR_NOTZERO(c)) { 
+      /* MPFR_PREC(b)>MPFR_PREC(a) : we have to round b+c */
+      k = bn - an;
 
-	    /* c is not zero */
-	    /* check whether mant(c)=1/2 or not */
+      /* first copy the 'an' most significant limbs of b to a */
+      MPN_COPY(ap, bp + k, an);
 
-            mp_limb_t *cp2 = cp + (cn-1);  /* highest limb */
+      difs = diff_exp % BITS_PER_MP_LIMB;
+      difw = an - (mp_exp_t) (diff_exp / BITS_PER_MP_LIMB);
+      MPFR_ASSERTN(difw <= 1);
 
-	    cc = *cp2 - MP_LIMB_T_HIGHBIT;
-            while (cc == 0 && cp2 > cp) cc = *--cp2;
+      nulp = -1;
 
-	    if (cc || ((ap[0] >> sh) & ONE)) goto add_one_ulp;
-	    /* mant(c) != 1/2 or mant(c) = 1/2: add 1 iff lsb(a)=1 */
-	  }
-	}
+      inv = rnd_mode == GMP_RNDN ? (sh ? ONE<<(sh-1) : MP_LIMB_T_HIGHBIT) : 0;
+      if (sh)
+      {
+        bb = *ap & ((ONE<<sh)-1);
+        *ap &= ~bb; /* truncate last bits */
+        if (inv)
+        {
+          r = bb >> (sh-1);
+          bb ^= inv;
+          inv = 0;
+        }
+        if (difw > 0)
+        {
+          /* c overlaps with the lowest limb of a */
+          MPFR_ASSERTN(difs > 0);
+          cc = cp[--cn];
+          bb += cc >> (BITS_PER_MP_LIMB - difs);
+        }
+        if (bb >= ONE<<sh)
+        {
+          nulp = 1;
+          bb -= ONE<<sh;
+        }
+        else if (bb < (ONE<<sh)-1)
+        {
+          nulp = 0;
+        }
       }
-      else if ((MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) || 
-	       (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) {
-	
-	/* round up */
-	if (MPFR_NOTZERO(c)) goto add_one_ulp;
+
+      if (nulp < 0 || (nulp > 0 && bb == 0))
+        while (1) /* look for carry and/or sticky bit */
+        {
+          if (k == 0)
+          {
+            if (nulp < 0)
+              nulp = 0;
+            else /* nulp == 1 */
+            {
+              if (difw > 0)
+                bb = cc << difs;
+              while (bb == 0 && cn)
+                bb = cp[--cn];
+            }
+            break;
+          }
+          bb = bp[--k];
+          if (inv)
+          {
+            r = bb >> (BITS_PER_MP_LIMB - 1);
+            bb ^= inv;
+            inv = 0;
+          }
+          if (difw >= 0)
+          {
+            mp_limb_t b2;
+            if (difw > 0)
+            {
+              b2 = bb;
+              bb = b2 + (cc << difs);
+              if (bb < b2)
+                nulp = 1;
+            }
+            if (cn == 0)
+            {
+              if (nulp < 0)
+                nulp = 0;
+              else /* nulp == 1 */
+                while (bb == 0 && k)
+                  bb = bp[--k];
+              break;
+            }
+            cc = cp[--cn];
+            if (difs)
+            {
+              difw = 1;
+              bb = b2 + (cc >> (BITS_PER_MP_LIMB - difs));
+            }
+            else
+            {
+              bb = b2 + cc;
+            }
+            if (bb < b2)
+              nulp = 1;
+          }
+          else
+            difw++;
+          if (nulp > 0 && bb != 0)
+            break;
+          if (nulp < 0 && bb != (mp_limb_t) (-1))
+          {
+            nulp = 0;
+            break;
+          }
+        } /* while */
+
+      /* if nulp == 1, sticky bit = bb != 0 */
+
+      if (rnd_mode == GMP_RNDN)
+      { /* r: old rounding bit of b
+           ulp(a)/2 were added if r = 0
+           ulp(a)/2 were subtracted if r = 1
+             r    carry  stbit  -->  nulp to add
+             0      0      1              0
+             0      1      0         0 or 1 (to even)
+             0      1      1              1
+             1      0      1              1
+             1      1      0         1 or 2 (to even)
+             1      1      1              2 [*]
+           [*] or only 1 if adding the first ulp changes the exponent */
+
+        if (nulp == 0)
+          bb = 1;
+        nulp += r;
+        if (!bb && (ap[0] & (ONE<<sh)))
+          nulp++;  /* sticky bit is 0 --> round to even */
       }
-      /* in the other cases (round to zero, or up/down with sign -/+),
-         nothing to do */
-    }
-    else
+      else
       {
-      /* MPFR_PREC(b)>MPFR_PREC(a) : we have to round b+c */      
-      k = bn - an;
-
-      /* first copy the 'an' most significant limbs of b to a */
-      MPN_COPY(ap, bp + k, an);
-      { /* treat all rounding modes together */
-	mp_limb_t c2old; long int cout=0; int sign=0;
-	if (sh) {
-	  cc = *ap & ((ONE<<sh)-1);
-	  *ap &= ~cc; /* truncate last bits */
-	}
-	else cc=0;
-
-        dif += sh; 
-	if (dif > 0)
-	  {
-	    cn--;
-	    c2old = cp[cn]; /* last limb from c considered */
-	    cout += mpn_add_1(&cc, &cc, 1, c2old >> (BITS_PER_MP_LIMB-dif));
-	  }
-	else
-	  c2 = c2old = 0;
-	if (sh && rnd_mode == GMP_RNDN)
-	    cout -= mpn_sub_1 (&cc, &cc, 1, ONE << (sh - 1));
-	if (cout == 0)
-	  {
-	   dif += BITS_PER_MP_LIMB;
-	   while (cout==0 && (k || cn)) {
-	     cout = (cc>(mp_limb_t)1) ? 2 : cc;
-	     cc = (k) ? bp[--k] : 0;
-	     if (sh==0) {
-	       cout -= mpn_sub_1(&cc, &cc, 1, ONE << (BITS_PER_MP_LIMB-1));
-	       sh = 0;
-	     }
-	     /* next limb from c to consider is cp[cn-1], with lower part of
-		c2old */
-	     if (cn && (dif >= 0))
-	       {
-		 c2 = c2old << dif;
-		 cn--;
-		 c2old = cp[cn];
-		 c2 += c2old >> (BITS_PER_MP_LIMB - dif);
-	       }
-	     else
-	       dif += BITS_PER_MP_LIMB;
-	     cout += mpn_add_1(&cc, &cc, 1, c2);
-	   }
-	}
-	if (cout==0) cout=(cc!=0);
-        sign = (MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU)
-	  || (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD) || (rnd_mode==GMP_RNDN);
-	/* round towards infinity if dif=1, towards zero otherwise */
-	if ((sign==1) && (cout>0)) goto add_one_ulp;
-	else if (rnd_mode==GMP_RNDN && cout==0 && (*ap & (ONE<<sh)))
-	  goto add_one_ulp;
+        if ((nulp == 0 || bb != 0) &&
+               ((rnd_mode == GMP_RNDU && MPFR_ISNONNEG(b))
+             || (rnd_mode == GMP_RNDD && MPFR_ISNEG(b))))
+          nulp++;
       }
-    }
-  }
-  else { 
+      goto add_one_ulp;
+    } /* bq > aq */
+  } /* aq <= diff_exp */
+  else {
     /* diff_exp < MPFR_PREC(a) : c overlaps with a by dif bits */
     /* first copy upper part of c into a (after shift) */
     unsigned char overlap;
-    
+
+    dif = MPFR_PREC(a) - diff_exp;
     k = (dif-1)/BITS_PER_MP_LIMB + 1; /* only the highest k limbs from c
 					 have to be considered */
     cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB + 1;
@@ -212,7 +291,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
 #ifdef DEBUG
     printf("MPFR_PREC(c)=%d\n", MPFR_PREC(c));
 #endif
-    if (dif<=MPFR_PREC(c)) { 
+    if (dif<=MPFR_PREC(c)) {
       /* c has to be truncated */
       dif = dif % BITS_PER_MP_LIMB;
       dif = (dif) ? BITS_PER_MP_LIMB-dif-sh : -sh;
@@ -231,7 +310,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
       else MPN_COPY(ap, cp+(cn-k), k);
       overlap=1;
     }
-    else { 
+    else {
 
       /* c is not truncated, but we have to fill low limbs with 0 */
 
@@ -239,9 +318,9 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
       overlap = diff_exp%BITS_PER_MP_LIMB;
 
       /* warning: a shift of zero bit is not allowed */
-      MPN_ZERO(ap, an-k-cn); 
-      if (overlap) { 
-	cc=mpn_rshift(ap+(an-k-cn), cp, cn, overlap); 
+      MPN_ZERO(ap, an-k-cn);
+      if (overlap) {
+	cc=mpn_rshift(ap+(an-k-cn), cp, cn, overlap);
 	if (an-k-cn>0) ap[an-k-cn-1]=cc;
       }
       else MPN_COPY(ap+(an-k-cn), cp, cn);
@@ -258,18 +337,18 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
     }
     else
       /* MPFR_PREC(b) > MPFR_PREC(a): we have to truncate b */
-      cc = mpn_add_n(ap, ap, bp+(bn-an), an); 
+      cc = mpn_add_n(ap, ap, bp+(bn-an), an);
 
-    if (cc) { 
+    if (cc) {
 
       /* shift one bit to the right */
 
       c3 = (ap[0]&1) && (MPFR_PREC(a)%BITS_PER_MP_LIMB==0);
       mpn_rshift(ap, ap, an, 1);
-      ap[an-1] += ONE<<(BITS_PER_MP_LIMB-1);
+      ap[an-1] += MP_LIMB_T_HIGHBIT;
       MPFR_EXP(a)++;
     }
-    
+
     /* remains to do the rounding */
 
 #ifdef DEBUG
@@ -279,7 +358,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
 	|| (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) {
 
       int kc;
-      
+
       /* four cases: overlap =
          (0) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) <= MPFR_PREC(a)
          (1) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) > MPFR_PREC(a)
@@ -309,7 +388,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
 	  }
 	  if (kc==0 && dif) {
 	    /* it still remains cp[0]<<(BITS_PER_MP_LIMB-dif) */
-	    if (k!=0) cout += mpn_add_1(&cc, bp+(--k), 1, 
+	    if (k!=0) cout += mpn_add_1(&cc, bp+(--k), 1,
 				      cp[0]<<(BITS_PER_MP_LIMB-dif));
 	    else cc = cp[0]<<(BITS_PER_MP_LIMB-dif);
 	    if ((cout==0 && cc==0) || (~cout==0 && ~cc==0)) cout=cc;
@@ -333,21 +412,21 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
 	case 0: /* only b to round */
 	  k=bn-an; dif=0;
 	  goto to_nearest;
-        
+
 	  /* otherwise the result is exact: nothing to do */
 	}
     }
     /* else nothing to do: round towards zero, i.e. truncate last sh bits */
-    else 
+    else
       *ap &= ~((ONE<<sh)-1);
   }
   goto end_of_add;
-    
+
   to_nearest: /* 0 <= sh < BITS_PER_MP_LIMB : number of bits of a to truncate
                  bp[k] : last significant limb from b
 		 bn : number of limbs of b
 	      */
-        /* c3=1 whenever b+c gave a carry out in most significant limb 
+        /* c3=1 whenever b+c gave a carry out in most significant limb
 	   and the least significant bit (shifted right) was 1.
 	   This can occur only when BITS_PER_MP_LIMB divides MPFR_PREC(a),
 	   i.e. sh=0.
@@ -358,7 +437,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
 	  c2 = (rnd_mode==GMP_RNDN) ? ONE<<(sh-1) : 0;
 	}
 	else /* sh=0: no bit to truncate */
-	  { 
+	  {
 	    if (k) cc = bp[--k]; else cc = 0;
 	    c2 = (rnd_mode==GMP_RNDN) ? ONE<<(BITS_PER_MP_LIMB-1) : 0;
 	    if (c3 && (cc || c2==0)) cc=c2+1; /* will force adding one ulp */
@@ -376,20 +455,35 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp)
         goto end_of_add;
 
   add_one_ulp: /* add one unit in last place to a */
-    cc = mpn_add_1(ap, ap, an, ONE<<sh);
-    if (cc) {
-      ap[an-1] = (mp_limb_t)1 << (BITS_PER_MP_LIMB-1);
-      MPFR_EXP(a)++;
+    while (nulp--)
+    {
+      cc = mpn_add_1(ap, ap, an, ONE<<sh);
+      if (cc)
+      {
+        mp_exp_t exp = MPFR_EXP(a);
+        if (exp == __mpfr_emax)
+        {
+          mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a));
+          break;
+        }
+        else
+        {
+          MPFR_EXP(a)++;
+          ap[an-1] = MP_LIMB_T_HIGHBIT;
+          if (rnd_mode == GMP_RNDN)
+            break;  /* because ulp is doubled */
+        }
+      }
     }
 
  end_of_add:
-  TMP_FREE(marker); 
+  TMP_FREE(marker);
   return;
 }
 
 void
 #if __STDC__
-mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) 
+mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
 #else
 mpfr_add (a, b, c, rnd_mode)
      mpfr_ptr a;
@@ -398,62 +492,76 @@ mpfr_add (a, b, c, rnd_mode)
      mp_rnd_t rnd_mode;
 #endif
 {
-  int diff_exp;
-
   if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) {
     MPFR_SET_NAN(a); return;
   }
-  
+
   MPFR_CLEAR_NAN(a);
 
-  if (MPFR_IS_INF(b)) 
-    { 
-      if (MPFR_IS_INF(c)) 
+  if (MPFR_IS_INF(b))
+    {
+      if (MPFR_IS_INF(c))
 	{
-	  if (MPFR_SIGN(b) == MPFR_SIGN(c)) 
-	    { 
-	      MPFR_SET_INF(a); 
+	  if (MPFR_SIGN(b) == MPFR_SIGN(c))
+	    {
+	      MPFR_SET_INF(a);
 	      if (MPFR_SIGN(a) != MPFR_SIGN(b)) MPFR_CHANGE_SIGN(a);
 	    }
 	  else
-	    MPFR_SET_NAN(a); 
+	    MPFR_SET_NAN(a);
 	}
       else
 	{
-	  MPFR_SET_INF(a); 
+	  MPFR_SET_INF(a);
 	  if (MPFR_SIGN(b) != MPFR_SIGN(a)) MPFR_CHANGE_SIGN(a);
 	}
       return;
     }
-  else 
+  else
     if (MPFR_IS_INF(c))
       {
-	MPFR_SET_INF(a); 
+	MPFR_SET_INF(a);
 	if (MPFR_SIGN(c) != MPFR_SIGN(a)) MPFR_CHANGE_SIGN(a);
 	return;
       }
 
-  MPFR_CLEAR_INF(a); /* clear Inf flag */
+  if (MPFR_IS_ZERO(b)) { mpfr_set(a, c, rnd_mode); return; }
+  if (MPFR_IS_ZERO(c)) { mpfr_set(a, b, rnd_mode); return; }
 
-  if (!MPFR_NOTZERO(b)) { mpfr_set(a, c, rnd_mode); return; }
-  if (!MPFR_NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; }
+  MPFR_CLEAR_INF(a); /* clear Inf flag */
 
-  diff_exp = MPFR_EXP(b)-MPFR_EXP(c);
-  if (MPFR_SIGN(b) * MPFR_SIGN(c) < 0) { /* signs differ, it's a subtraction */
-    if (diff_exp<0) {
-      mpfr_sub1(a, c, b, rnd_mode, -diff_exp);
+  if (MPFR_SIGN(b) != MPFR_SIGN(c))
+  { /* signs differ, it's a subtraction */
+    if (MPFR_EXP(b) < MPFR_EXP(c))
+    {
+      mpfr_sub1(a, c, b, rnd_mode,
+                (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b));
     }
-    else if (diff_exp>0) mpfr_sub1(a, b, c, rnd_mode, diff_exp);
-    else { /* diff_exp=0 */
-      diff_exp = mpfr_cmp3(b,c,-1);
+    else if (MPFR_EXP(b) > MPFR_EXP(c))
+    {
+      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 diff_exp = mpfr_cmp3(b,c,-1);
       /* if b>0 and diff_exp>0 or b<0 and diff_exp<0: abs(b) > abs(c) */
       if (diff_exp==0) MPFR_SET_ZERO(a);
       else if (diff_exp * MPFR_SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0);
       else mpfr_sub1(a, c, b, rnd_mode, 0);
     }
   }
-  else /* signs are equal, it's an addition */
-    if (diff_exp<0) mpfr_add1(a, c, b, rnd_mode, -diff_exp);
-    else mpfr_add1(a, b, c, rnd_mode, diff_exp);
+  else
+  { /* signs are equal, it's an addition */
+    if (MPFR_EXP(b) < MPFR_EXP(c))
+    {
+      mpfr_add1(a, c, b, rnd_mode,
+                (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b));
+    }
+    else
+    {
+      mpfr_add1(a, b, c, rnd_mode,
+                (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c));
+    }
+  }
 }
-