diff -Naur mpfr-2.1.0-p5/sub1.c mpfr-2.1.0-p6/sub1.c
--- mpfr-2.1.0-p5/sub1.c	2004-03-26 13:27:16.000000000 +0000
+++ mpfr-2.1.0-p6/sub1.c	2005-01-23 23:18:41.392997864 +0000
@@ -1,6 +1,6 @@
 /* mpfr_sub1 -- internal function to perform a "real" subtraction
 
-Copyright 2001, 2002, 2003, 2004 Free Software Foundation.
+Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation.
 Contributed by the Spaces project, INRIA Lorraine.
 
 This file is part of the MPFR Library.
@@ -20,7 +20,6 @@
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
-
 #include "mpfr-impl.h"
 
 /* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
@@ -81,6 +80,60 @@
   else
     MPFR_SET_SAME_SIGN(a,b);
 
+  /* Check if c is too small.
+     A more precise test is to replace 2 by   
+      (rnd == GMP_RNDN) + mpfr_power2_raw (b) 
+      but it is more expensive and not very usefull */
+  if (MPFR_UNLIKELY (MPFR_GET_EXP (c) <= MPFR_GET_EXP (b) 
+		     - (mp_exp_t) MAX (MPFR_PREC (a), MPFR_PREC (b)) - 2))
+    {
+      /* Remember, we can't have an exact result! */
+      /*   A.AAAAAAAAAAAAAAAAA
+	 = B.BBBBBBBBBBBBBBB
+	  -                     C.CCCCCCCCCCCCC */
+      /* A = S*ABS(B) +/- ulp(a) */
+      inexact = mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (a)); 
+      if (inexact == 0)
+	{
+	  /* a = b (Exact)
+	     But we know it isn't (Since we have to remove `c')
+	     So if we round to Zero, we have to remove one ulp.
+	     Otherwise the result is correctly rounded. */
+	  if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a))) {
+	    mpfr_nexttoward (a, c);
+	    return -MPFR_INT_SIGN (a);
+	  }
+	  return MPFR_INT_SIGN (a);
+	}
+      else
+	{
+	  /*   A.AAAAAAAAAAAAAA
+             = B.BBBBBBBBBBBBBBB
+  	      -                   C.CCCCCCCCCCCCC */
+	  /* It isn't exact so Prec(b) > Prec(a) and the last
+	     Prec(b)-Prec(a) bits of `b' are not zeros.
+	     Which means that removing c from b can't generate a carry
+	     execpt in case of even rounding.
+	     In all other case the result and the inexact flag should be
+	     correct (We can't have an exact result).
+	     In case of EVEN rounding:
+	       1.BBBBBBBBBBBBBx10
+	     -                     1.CCCCCCCCCCCC 
+	     = 1.BBBBBBBBBBBBBx01  Rounded to Prec(b)
+	     = 1.BBBBBBBBBBBBBx    Nearest / Rounded to Prec(a)
+	     Set gives:
+	       1.BBBBBBBBBBBBB0   if inexact == EVEN_INEX  (x == 0)
+  	       1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1) 
+	     which means we get a wrong rounded result if x==1, 
+	     i.e. inexact= MPFR_EVEN_INEX */
+	  if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX*MPFR_INT_SIGN (a))) {
+	    mpfr_nexttoward (a, c);
+	    inexact = -MPFR_INT_SIGN (a);
+	  }
+	  return inexact;
+	}
+    }
+
   diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
 
   /* reserve a space to store b aligned with the result, i.e. shifted by
@@ -136,19 +189,20 @@
     }
 
 #ifdef DEBUG
-  printf("shift_b=%u shift_c=%u\n", shift_b, shift_c);
+  printf("shift_b=%u shift_c=%u diffexp=%lu\n", shift_b, shift_c, diff_exp);
 #endif
 
   MPFR_ASSERTD (ap != cp);
   MPFR_ASSERTD (bp != cp);
 
   /* here we have shift_c = (diff_exp - cancel) % BITS_PER_MP_LIMB,
+        0 <= shift_c < BITS_PER_MP_LIMB
      thus we want cancel2 = ceil((cancel - diff_exp) / BITS_PER_MP_LIMB) */
 
   cancel2 = (long int) (cancel - (diff_exp - shift_c)) / BITS_PER_MP_LIMB;
   /* the high cancel2 limbs from b should not be taken into account */
 #ifdef DEBUG
-  printf("cancel=%u cancel1=%u cancel2=%d\n", cancel, cancel1, cancel2);
+  printf("cancel=%lu cancel1=%lu cancel2=%ld\n", cancel, cancel1, cancel2);
 #endif
 
   /*               ap[an-1]        ap[0]
@@ -276,6 +330,7 @@
   bn -= an + cancel1;
   cn0 = cn;
   cn -= (long int) an + cancel2;
+
 #ifdef DEBUG
   printf("last %d bits from a are %lu, bn=%ld, cn=%ld\n", sh, carry, bn, cn);
 #endif