+ Add function mpfr_print_mantissa_binary, for debugging reason.

+ Rename MPFR_ALLOC_SIZE in MPFR_MALLOC_SIZE.
+ Add conditionnal -DSMALL directive in mpfr-impl.h.
+ Add new function: sub1sp.
   Substraction in case of all the ops have the same prec.
+ Add its test (tsub1sp).
+ Modify a few the tests to avoid comparing mpfr results with double, for portability reason.


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

29 files changed, 1185 insertions, 54 deletions

diff --git a/Makefile.am b/Makefile.am
index 921281b8d..bef929333 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -7,7 +7,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h
 
 lib_LIBRARIES = libmpfr.a
 
-libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-test.h log_b2.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_one_ulp.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c  factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c
+libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-test.h log_b2.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_one_ulp.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c  factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c
 
 libmpfr_a_LIBADD = @LIBOBJS@
 
diff --git a/add.c b/add.c
index 42f5841ac..01e7455bf 100644
--- a/add.c
+++ b/add.c
@@ -79,9 +79,13 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
   MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));
   MPFR_CLEAR_FLAGS(a); /* clear flags */
 
-  if (MPFR_SIGN(b) != MPFR_SIGN(c))
+  if (MPFR_UNLIKELY(MPFR_SIGN(b) != MPFR_SIGN(c)))
     { /* signs differ, it's a subtraction */
-      return mpfr_sub1(a, b, c, rnd_mode);
+      if (MPFR_LIKELY(MPFR_PREC(a) == MPFR_PREC(b) 
+		      && MPFR_PREC(b) == MPFR_PREC(c)))
+	return mpfr_sub1sp(a,b,c,rnd_mode);
+      else
+	return mpfr_sub1(a, b, c, rnd_mode);
     }
   else
     { /* signs are equal, it's an addition */
diff --git a/clear.c b/clear.c
index 5db4f561e..df4b6ea6a 100644
--- a/clear.c
+++ b/clear.c
@@ -30,6 +30,6 @@ mpfr_clear (mpfr_ptr m)
 {
   /* be careful to always free an entire number of limbs */
   (*__gmp_free_func) 
-    (MPFR_GET_REAL_PTR(m), MPFR_ALLOC_SIZE(MPFR_GET_ALLOC_SIZE(m)));
+    (MPFR_GET_REAL_PTR(m), MPFR_MALLOC_SIZE(MPFR_GET_ALLOC_SIZE(m)));
   MPFR_MANT(m) = NULL;
 }
diff --git a/init2.c b/init2.c
index 9efcf6c47..dce3466a1 100644
--- a/init2.c
+++ b/init2.c
@@ -37,7 +37,7 @@ mpfr_init2 (mpfr_ptr x, mpfr_prec_t p)
   MPFR_ASSERTN(p >= MPFR_PREC_MIN && p <= MPFR_PREC_MAX);
 
   xsize = (mp_size_t) ((p - 1) / BITS_PER_MP_LIMB) + 1;
-  tmp   = (mp_ptr) (*__gmp_allocate_func)(MPFR_ALLOC_SIZE(xsize));
+  tmp   = (mp_ptr) (*__gmp_allocate_func)(MPFR_MALLOC_SIZE(xsize));
 
   MPFR_PREC(x) = p;
   MPFR_EXP (x) = MPFR_EXP_INVALID; /* make sure that the exp field has a
diff --git a/mpfr-impl.h b/mpfr-impl.h
index 03a5293ba..9d0828d12 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -105,6 +105,28 @@ typedef unsigned long int       mpfr_exp_unsigned_t;
 
 #define MPFR_INTPREC_MAX (ULONG_MAX & ~(unsigned long) (BITS_PER_MP_LIMB - 1))
 
+/* Redefine MPN_COPY if compilation for small vars
+   FIXME: Usefull ?*/
+#ifdef SMALL
+#undef MPN_COPY
+#define MPN_COPY(dest, src, n) \
+ do {                          \
+  if ((n) != 0)                \
+    {                          \
+      mp_size_t __n = ((n)-1); \
+      mp_ptr __dst = (dest);   \
+      mp_srcptr __src = (src); \
+      mp_limb_t __x = *__src++;\
+      if (__n)                 \
+	do {                   \
+	  *__dst++ = __x;      \
+	  __x      = *__src++; \
+	} while (--__n);       \
+      *__dst = __x;            \
+    }                          \
+ } while (0)                   
+#endif
+
 /* Assertions */
 
 /* Compile with -DWANT_ASSERT to check all assert statements */
@@ -309,6 +331,7 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO
  */
 
 #define MPFR_PREC(x) ((x)->_mpfr_prec)
+
 #define MPFR_EXP(x) ((x)->_mpfr_exp)
 #define MPFR_MANT(x) ((x)->_mpfr_d)
 
@@ -380,14 +403,16 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO
 #define MPFR_RET_NAN return (__gmpfr_flags |= MPFR_FLAGS_NAN), 0
 
 /* Heap Memory gestion */ /* Old ABSSIZE */
-#define MPFR_GET_ALLOC_SIZE(x) ( ((mp_size_t*) MPFR_MANT(x))[-1] + 0)
-#define MPFR_SET_ALLOC_SIZE(x, n) ( ((mp_size_t*) MPFR_MANT(x))[-1] = n)
-#define MPFR_ALLOC_SIZE(s) \
-  ((size_t) (sizeof(mp_size_t) + BYTES_PER_MP_LIMB*(s)))
+#define MPFR_GET_ALLOC_SIZE(x) \
+ ( ((mp_size_t*) MPFR_MANT(x))[-1] + 0)
+#define MPFR_SET_ALLOC_SIZE(x, n) \
+ ( ((mp_size_t*) MPFR_MANT(x))[-1] = n)
+#define MPFR_MALLOC_SIZE(s) \
+  ((size_t) (sizeof(mp_size_t)*2 + BYTES_PER_MP_LIMB*(s)))
 #define MPFR_SET_MANT_PTR(x,p) \
-   (MPFR_MANT(x) = (mp_limb_t*) ((mp_size_t*) p + 1))
+   (MPFR_MANT(x) = (mp_limb_t*) ((mp_size_t*) p + 2))
 #define MPFR_GET_REAL_PTR(x) \
-   ((mp_limb_t*) ((mp_size_t*) MPFR_MANT(x) - 1))
+   ((mp_limb_t*) ((mp_size_t*) MPFR_MANT(x) - 2))
 
 /* Temporary memory gestion */
 /* temporary allocate 1 limb at xp, and initialize mpfr variable x */
@@ -437,6 +462,7 @@ void mpfr_restore_emin_emax _MPFR_PROTO ((void));
 
 int mpfr_add1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_sub1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
+int mpfr_sub1sp _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_can_round_raw _MPFR_PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t,
 				     mp_rnd_t, mp_rnd_t, mp_prec_t));
 
@@ -460,6 +486,8 @@ long mpn_exp _MPFR_PROTO ((mp_limb_t *, mp_exp_t *, int,
 			   mp_exp_t, size_t));
 
 void mpfr_print_binary _MPFR_PROTO ((mpfr_srcptr));
+void mpfr_print_mant_binary _MPFR_PROTO ((const char *, const mp_limb_t *,
+					  mp_prec_t));
 void mpfr_set_str_binary _MPFR_PROTO ((mpfr_ptr, __gmp_const char *));
 
 int mpfr_round_raw _MPFR_PROTO ((mp_limb_t *, mp_limb_t *,
@@ -474,7 +502,8 @@ int mpfr_round_raw_4 _MPFR_PROTO ((mp_limb_t *, mp_limb_t *,
 #define mpfr_round_raw2(xp, xn, neg, r, prec) \
   mpfr_round_raw_2(0, (xp), (xn) * BITS_PER_MP_LIMB, (neg), (prec), (r) )
 
-int mpfr_check(mpfr_srcptr);
+int mpfr_check _MPFR_PROTO ((mpfr_srcptr));
+int mpfr_sub1sp _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
 
 #if defined (__cplusplus)
 }
diff --git a/mul.c b/mul.c
index 3aa8266d5..d7b180dc5 100644
--- a/mul.c
+++ b/mul.c
@@ -94,14 +94,14 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
 
   /* FIXME: Usefull since we do an exponent check after ?
    * It is usefull iff the precision is big, there is an overflow
-   * and we are doing further mults... Probable ? */
-  /*
+   * and we are doing further mults...*/
+#ifdef HUGE
   if (MPFR_UNLIKELY(bx + cx > __gmpfr_emax + 1))
     return mpfr_set_overflow (a, rnd_mode, sign_product);
   if (MPFR_UNLIKELY(bx + cx < __gmpfr_emin - 2))
   return mpfr_set_underflow (a, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode,
 			     sign_product);
-  */
+#endif
 
   aq = MPFR_PREC(a);
   bq = MPFR_PREC(b);
diff --git a/print_raw.c b/print_raw.c
index a7427ee58..09da6fcdc 100644
--- a/print_raw.c
+++ b/print_raw.c
@@ -23,9 +23,6 @@ MA 02111-1307, USA. */
 #include <stdio.h>
 #include <limits.h>
 
-#include "gmp.h"
-#include "gmp-impl.h"
-#include "mpfr.h"
 #include "mpfr-impl.h"
 
 void
@@ -76,3 +73,26 @@ mpfr_print_binary (mpfr_srcptr x)
         }
     }
 }
+
+void 
+mpfr_print_mant_binary(const char *str, const mp_limb_t *p, mp_prec_t r)
+{
+  int i;
+  mp_prec_t count = 0;
+  char c;
+  mp_size_t n = (r - 1) / BITS_PER_MP_LIMB + 1;
+
+  printf("%s ", str);
+  for(n-- ; n>=0 ; n--)
+    {
+      for(i = BITS_PER_MP_LIMB-1 ; i >=0 ; i--)
+        {
+          c = (p[n] & (1<<i)) ? '1' : '0';
+          putchar(c);
+          count++;
+          if (count == r)
+            putchar('[');
+        }
+    }
+  putchar('\n');
+}
diff --git a/round_prec.c b/round_prec.c
index a76594f4c..ca14a6eaf 100644
--- a/round_prec.c
+++ b/round_prec.c
@@ -63,7 +63,7 @@ mpfr_prec_round (mpfr_ptr x, mp_prec_t prec, mp_rnd_t rnd_mode)
     {
       /* Realloc mantissa */
       mp_ptr tmp = (mp_ptr) (*__gmp_reallocate_func)
-        (MPFR_GET_REAL_PTR(x),  MPFR_ALLOC_SIZE(ow), MPFR_ALLOC_SIZE(nw));
+        (MPFR_GET_REAL_PTR(x),  MPFR_MALLOC_SIZE(ow), MPFR_MALLOC_SIZE(nw));
       MPFR_SET_MANT_PTR(x, tmp); /* mant ptr must be set before alloc size */
       MPFR_SET_ALLOC_SIZE(x, nw); /* new number of allocated limbs */
     }
diff --git a/set_prec.c b/set_prec.c
index 6c9828931..113da143b 100644
--- a/set_prec.c
+++ b/set_prec.c
@@ -42,7 +42,7 @@ mpfr_set_prec (mpfr_ptr x, mpfr_prec_t p)
   if (xsize > xoldsize)
     {
       tmp = (mp_ptr) (*__gmp_reallocate_func)
-        (MPFR_GET_REAL_PTR(x), MPFR_ALLOC_SIZE(xoldsize), MPFR_ALLOC_SIZE(xsize));
+        (MPFR_GET_REAL_PTR(x), MPFR_MALLOC_SIZE(xoldsize), MPFR_MALLOC_SIZE(xsize));
       MPFR_SET_MANT_PTR(x, tmp);
       MPFR_SET_ALLOC_SIZE(x, xsize);
     }
diff --git a/sub.c b/sub.c
index 0c9128bce..409f87a6c 100644
--- a/sub.c
+++ b/sub.c
@@ -80,9 +80,13 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
   MPFR_CLEAR_FLAGS(a);
   MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));
 
-  if (MPFR_SIGN(b) == MPFR_SIGN(c))
+  if (MPFR_LIKELY(MPFR_SIGN(b) == MPFR_SIGN(c)))
     { /* signs are equal, it's a real subtraction */
-      return mpfr_sub1(a, b, c, rnd_mode);
+      if (MPFR_LIKELY(MPFR_PREC(a) == MPFR_PREC(b)
+                      && MPFR_PREC(b) == MPFR_PREC(c)))
+        return mpfr_sub1sp(a,b,c,rnd_mode);
+      else
+	return mpfr_sub1(a, b, c, rnd_mode);
     }
   else
     { /* signs differ, it's an addition */
diff --git a/sub1.c b/sub1.c
index 00c670466..c7d759bea 100644
--- a/sub1.c
+++ b/sub1.c
@@ -40,7 +40,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
   mp_size_t cancel2, an, bn, cn, cn0;
   mp_limb_t *ap, *bp, *cp;
   mp_limb_t carry, bb, cc, borrow = 0;
-  int inexact = 0, shift_b, shift_c, is_exact = 1, down = 0, add_exp = 0;
+  int inexact, shift_b, shift_c, is_exact = 1, down = 0, add_exp = 0;
   int sh, k;
   TMP_DECL(marker);
 
diff --git a/sub1sp.c b/sub1sp.c
new file mode 100644
index 000000000..91837bee2
--- /dev/null
+++ b/sub1sp.c
@@ -0,0 +1,711 @@
+/* mpfr_sub1sp -- internal function to perform a "real" substraction
+   All the op must have the same precision
+
+Copyright 2003 Free Software Foundation.
+Contributed by the Spaces project, INRIA Lorraine.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the MPFR Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <stdio.h>
+#include "mpfr-impl.h"
+#include "longlong.h"
+
+/* 
+   compute sgn(b)*(|b| - |c|) if |b|>|c| else -sgn(b)*(|c| -|b|)
+   Returns 0 iff result is exact,
+   a negative value when the result is less than the exact value,
+   a positive value otherwise.
+*/
+
+/*#define DEBUG*/
+
+#ifdef DEBUG
+# undef DEBUG
+# define DEBUG(x) (x)
+#else
+# define DEBUG(x) /**/
+#endif
+
+/* Rounding Sub */
+      
+/* A0...Ap-1
+ *          Cp Cp+1 ....
+ *             <- C'p+1 ->
+ * Cp = -1 if calculated from c mantissa
+ * Cp = 0  if 0 from a or c
+ * Cp = 1  if calculated from a.
+ * C'p+1 = Frist bit not null or 0 if there isn't one 
+ *
+ * Can't have Cp=-1 and C'p+1=1*/
+
+/* RND = GMP_RNDZ:
+ *  + if Cp=0 and C'p+1=0,1,  Truncate.
+ *  + if Cp=0 and C'p+1=-1,   SubOneUlp
+ *  + if Cp=-1,               SubOneUlp
+ *  + if Cp=1,                AddOneUlp
+ * RND = GMP_RNDA (Away)
+ *  + if Cp=0 and C'p+1=0,-1, Truncate
+ *  + if Cp=0 and C'p+1=1,    AddOneUlp
+ *  + if Cp=1,                AddOneUlp
+ *  + if Cp=-1,               Truncate
+ * RND = GMP_RNDN
+ *  + if Cp=0,                Truncate
+ *  + if Cp=1 and C'p+1=1,    AddOneUlp
+ *  + if Cp=1 and C'p+1=-1,   Truncate
+ *  + if Cp=1 and C'p+1=0,    Truncate if Ap-1=0, AddOneUlp else
+ *  + if Cp=-1 and C'p+1=-1,  SubOneUlp
+ *  + if Cp=-1 and C'p+1=0,   Truncate if Ap-1=0, SubOneUlp else
+ * 
+ * If AddOneUlp:
+ *   If carry, then it is 11111111111 + 1 = 10000000000000
+ *      ap[n-1]=MPFR_HIGHT_BIT
+ * If SubOneUlp:
+ *   If we lose one bit, it is 1000000000 - 1 = 0111111111111
+ *      Then shift, and put as last bit x which is calculated 
+ *              according Cp, Cp-1 and rnd_mode.
+ * If Truncate,
+ *    If it is a power of 2,
+ *       we may have to suboneulp in some special cases.
+ *
+ * To simplify, we don't use Cp = 1.
+ *
+ */
+
+int
+mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
+{
+  mp_exp_t bx,cx;
+  mp_exp_unsigned_t d;
+  mpfr_prec_t p, sh, cnt;
+  mp_size_t n;
+  mp_limb_t *ap, *bp, *cp;
+  mp_limb_t limb;
+  int inexact;
+  int bcp,bcp1; /* Cp and C'p+1 */
+  int bbcp, bbcp1; /* Cp+1 and C'p+2 */
+  TMP_DECL(marker);
+
+  TMP_MARK(marker);
+  
+  MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
+  MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));
+
+  /* Read prec and num of limbs */
+  p = MPFR_PREC(b);
+  n = (p-1)/BITS_PER_MP_LIMB+1;
+
+  /* Fast cmp of |b| and |c|*/
+  bx = MPFR_GET_EXP (b);
+  cx = MPFR_GET_EXP (c);
+  if (MPFR_UNLIKELY(bx == cx))
+    {
+      mp_size_t k = n - 1;
+      /* Check mantissa since exponent are equals */
+      bp = MPFR_MANT(b);
+      cp = MPFR_MANT(c);
+      while (k>=0 && bp[k] == cp[k])
+        k--;
+      if (MPFR_UNLIKELY(k < 0))
+        /* b == c ! */
+        {
+	  /* Return exact number 0 */
+	  if (rnd_mode == GMP_RNDD)
+	    MPFR_SET_NEG(a);
+	  else
+	    MPFR_SET_POS(a);
+	  MPFR_SET_ZERO(a);
+	  MPFR_RET(0);
+	}
+      else if (bp[k] > cp[k])
+        goto BGreater;
+      else
+        {
+	  MPFR_ASSERTD(bp[k]<cp[k]);
+	  goto CGreater;
+	}
+    }   
+  else if (MPFR_UNLIKELY(bx < cx))
+    {
+      /* Swap b and c and set sign */
+      mpfr_srcptr t;
+      mp_exp_t tx;
+    CGreater:
+      MPFR_SET_OPPOSITE_SIGN(a,b);
+      t  = b;  b  = c;  c  = t;
+      tx = bx; bx = cx; cx = tx;
+    }
+  else
+    {
+      /* b > c */
+    BGreater:
+      MPFR_SET_SAME_SIGN(a,b);
+    }
+
+  /* Now b > c */
+  MPFR_ASSERTD(bx >= cx);
+  d = (mp_exp_unsigned_t) bx - cx;
+  DEBUG( printf("New with diff=%lu\n", d) );
+
+  /* Transform RNDU and RNDD to RNDA or RNDZ */
+  if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_NEG(a)))
+    rnd_mode = GMP_RNDZ;
+
+  if (d == 0)
+    {
+      /* <-- b -->	
+	 <-- c --> : exact sub */
+      ap = MPFR_MANT(a);
+      mpn_sub_n (ap, MPFR_MANT(b), MPFR_MANT(c), n);
+      /* Normalize */
+    ExactNormalize:
+      limb = ap[n-1];
+      if (limb)
+	{
+	  /* First limb is not zero. */
+	  count_leading_zeros(cnt, limb);
+	  if (MPFR_LIKELY(cnt))
+	    {
+	      mpn_lshift(ap, ap, n, cnt); /* Normalize number */
+	      bx -= cnt; /* Update final expo */
+	    }
+	  sh = (-p % BITS_PER_MP_LIMB);
+	  if (MPFR_LIKELY(sh))
+	    ap[0] &= ~((MPFR_LIMB_ONE << sh) - MPFR_LIMB_ONE);
+	}
+      else
+	{
+	  /* First limb is zero */
+	  mp_size_t k = n-1, len;
+	  /* Find the first limb not equal to zero.
+	     FIXME: It is assume it exists (since |b| > |c| and same prec)*/
+	  do {
+	    MPFR_ASSERTD( k > 0 );
+	    limb = ap[--k];
+	  } while (limb == 0);
+	  MPFR_ASSERTD(limb != 0);
+          count_leading_zeros(cnt, limb);
+	  k++;
+	  len = n - k; /* Number of last limb */
+	  MPFR_ASSERTD(k >= 0);
+          mpn_lshift(ap+len, ap, k, cnt); /* Normalize the High Limb*/
+	  MPN_ZERO(ap, len); /* Zeroing the last limbs */
+	  bx -= cnt + len*BITS_PER_MP_LIMB; /* Update Expo */
+	  sh = (-p % BITS_PER_MP_LIMB);
+	  if (MPFR_LIKELY(sh))
+	    ap[len] &= ~((MPFR_LIMB_ONE << sh) - MPFR_LIMB_ONE);
+	}
+      /* Check expo overflow */
+      if (MPFR_UNLIKELY(bx < __gmpfr_emin))
+	{
+	  TMP_FREE(marker);
+	  /* inexact=0 */
+	  DEBUG( printf("(D==0 Underflow)\n") );
+	  if (rnd_mode == GMP_RNDN &&
+	      (bx < __gmpfr_emin - 1 ||
+	       (/*inexact >= 0 &&*/ mpfr_powerof2_raw (a))))
+	    rnd_mode = GMP_RNDZ;
+          return mpfr_set_underflow (a, rnd_mode, MPFR_SIGN(a));
+	}
+      MPFR_SET_EXP (a, bx);
+      /* No rounding is necessary since the result is exact */
+      MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
+      TMP_FREE(marker);
+      return 0;
+    }
+  else if (d == 1)
+    {
+      /* | <-- b -->
+	 |  <-- c --> */
+      mp_limb_t c0, mask;
+      mp_size_t k;
+      sh = -p % BITS_PER_MP_LIMB;
+      /* If we lose at least one bit, compute 2*b-c (Exact)
+       * else compute b-c/2 */
+      bp = MPFR_MANT(b);
+      cp = MPFR_MANT(c);
+      k = n-1;
+      limb = bp[k] - cp[k]/2;
+      if (limb > MPFR_LIMB_HIGHBIT)
+	{
+	  /* We can't lose precision: compute b-c/2 */
+	  /* Shift c in the allocated temporary block */
+	SubD1NoLose:
+          c0 = cp[0] & (MPFR_LIMB_ONE<<sh);
+	  cp = (mp_limb_t*) TMP_ALLOC(n * BYTES_PER_MP_LIMB);
+	  mpn_rshift(cp, MPFR_MANT(c), n, 1);
+          if (c0 == 0)
+            {
+	      /* Result is exact: no need of rounding! */
+	      ap = MPFR_MANT(a);
+	      mpn_sub_n (ap, bp, cp, n);
+	      MPFR_SET_EXP(a, bx); /* No expo overflow! */
+	      /* No truncate or normalize is needed */
+	      MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
+	      /* No rounding is necessary since the result is exact */
+	      TMP_FREE(marker);
+	      return 0;
+	    }
+	  ap = MPFR_MANT(a);
+	  mask = ~((MPFR_LIMB_ONE << sh) - MPFR_LIMB_ONE);
+	  cp[0] &= mask; /* Delete last bit of c */
+	  mpn_sub_n (ap, bp, cp, n);
+	  MPFR_SET_EXP(a, bx);                 /* No expo overflow! */
+	  MPFR_ASSERTD( !(ap[0] & ~mask) );    /* Check last bits */
+	  /* No normalize is needed */
+	  MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
+	  /* Rounding is necessary since c0 = 1*/	  
+	  /* Cp =-1 and C'p+1=0 */
+	  bcp = 1; bcp1 = 0;
+	  switch (rnd_mode)
+	    {
+	    case GMP_RNDN:
+	      /* Even Rule apply: Check Ap-1 */
+	      if ((ap[0] & (MPFR_LIMB_ONE<<sh)) == 0)
+		goto truncate;
+	    case GMP_RNDZ:
+	      goto sub_one_ulp;
+	    default:
+	      goto truncate;
+	    }
+	}
+      else if (limb < MPFR_LIMB_HIGHBIT)
+	{
+	  /* We lose at least one bit of prec */
+	  /* Calcul of 2*b-c (Exact) */
+          /* Shift b in the allocated temporary block */
+	SubD1Lose:
+          bp = (mp_limb_t*) TMP_ALLOC(n * BYTES_PER_MP_LIMB);
+          mpn_lshift(bp, MPFR_MANT(b), n, 1);
+          ap = MPFR_MANT(a);
+          mpn_sub_n (ap, bp, cp, n);
+	  bx--;
+	  goto ExactNormalize;
+	}
+      else
+	{
+	  /* Case: limb = 100000000000 */
+	  /* Check while b[k] == c'[k]/2 (C' is C shifted by 1) */
+	  /* If b[k]<c'[k] => We lose at least one bit*/
+	  /* If b[k]>c'[k] => We don't lose any bit */
+	  /* If k==-1 => We don't lose any bit 
+	     AND the result is 100000000000 0000000000 00000000000 */
+	  mp_limb_t carry;
+	  do {
+	    carry = cp[k]&MPFR_LIMB_ONE;
+	    k--;
+	  } while (k>=0 && 
+		   bp[k]==(carry=cp[k]/2+(carry<<(BITS_PER_MP_LIMB-1))));
+	  if (MPFR_UNLIKELY(k<0))
+	    {
+	      /* If carry then (sh==0 and Virtual c'[-1] > Virtual b[-1]) */
+	      if (MPFR_UNLIKELY(carry)) /* carry = cp[0]&MPFR_LIMB_ONE */
+		{
+		  MPFR_ASSERTD(sh == 0);
+		  goto SubD1Lose;
+		}
+	      /* Result is a power of 2 */
+	      ap = MPFR_MANT(a);
+	      MPN_ZERO(ap, n);
+              ap[n-1] = MPFR_LIMB_HIGHBIT;
+	      MPFR_SET_EXP(a, bx); /* No expo overflow! */
+	      /* No Normalize is needed*/
+	      /* No Rounding is needed */
+	      TMP_FREE(marker);
+	      return 0;
+	    }
+	  /* carry = cp[k]/2+(cp[k-1]&1)<<(BITS_PER_MP_LIMB-1) = c'[k]*/
+	  else if (bp[k]>carry)
+	    goto SubD1NoLose;
+	  else 
+	    {
+	      MPFR_ASSERTD(bp[k]<carry);
+	      goto SubD1Lose;
+	    }
+	}
+    }
+  else if (d >= p)
+    {
+      ap = MPFR_MANT(a);
+      sh = (-p) % BITS_PER_MP_LIMB;
+      /* We can't set A before since we use cp for rounding... */
+      /* Perform rounding: check if a=b or a=b-ulp(b) */
+      if (MPFR_UNLIKELY(d == p))
+	{
+	  /* cp == -1 and c'p+1 = ? */
+	  bcp  = 1;
+	  /* We need Cp+1 later for a very improbable case. */
+	  bbcp = (MPFR_MANT(c)[n-1] & (MPFR_LIMB_ONE<<(BITS_PER_MP_LIMB-2)));
+	  /* We need also C'p+1 for an even more unprobable case... */
+	  if (bbcp)
+	    bcp1 = 1;
+	  else
+	    {
+              cp = MPFR_MANT(c);
+              if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
+                {
+                  mp_size_t k = n-1;
+                  do {
+                    k--;
+                  } while (k>=0 && cp[k]==0);
+                  bcp1 = (k>=0);      
+                }
+	      else
+		bcp1 = 1;
+	    }
+	  DEBUG( printf("(D=P) Cp=-1 Cp+1=%d C'p+1=%d \n", bbcp, bcp1) );
+	  bp = MPFR_MANT(b);
+	  switch (rnd_mode)
+	    {
+	    case GMP_RNDN:
+	      if (bcp && bcp1==0)		
+		/* Cp=-1 and C'p+1=0: Even rule Apply! */
+		/* Check Ap-1 = Bp-1 */
+		if ((bp[0] & (MPFR_LIMB_ONE<<sh)) == 0)
+		  goto copy_truncate;
+	    case GMP_RNDZ:
+	      MPN_COPY(ap, bp, n);
+	      goto sub_one_ulp;
+	    default:
+	    copy_truncate:
+	      MPN_COPY(ap, bp, n);
+	      goto truncate;
+	    }
+	}
+      else
+	{
+	  /* Cp=0, Cp+1=-1 if d==p+1, C'p+1=-1 */
+          bcp = 0; bbcp = (d==p+1); bcp1 = 1;
+          DEBUG( printf("(D>P) Cp=%d Cp+1=%d C'p+1=%d\n", bcp, bbcp, bcp1) );
+	  /* Need to compute C'p+2 if d==p+1 and if rnd_mode=NEAREST
+	     (Because of a very improbable case) */
+	  if (MPFR_UNLIKELY(d==p+1 && rnd_mode==GMP_RNDN))
+	    {
+	      cp = MPFR_MANT(c);
+	      if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
+		{
+		  mp_size_t k = n-1;
+		  do {
+		    k--;
+		  } while (k>=0 && cp[k]==0);
+		  bbcp1 = (k>=0);
+		}
+	      else
+		bbcp1 = 1;
+	      DEBUG( printf("(D>P) C'p+2=%d\n", bbcp1) );
+	    }
+	  /* Copy mantissa B in A */
+	  MPN_COPY(ap, MPFR_MANT(b), n);
+	  /* Round */
+	  if (rnd_mode == GMP_RNDZ)
+	    goto sub_one_ulp;
+	  else /* rnd_mode = AWAY or NEAREST */
+	    goto truncate;
+	}
+    }
+  else
+    {
+      mpfr_prec_t dm;
+      mp_size_t m;
+      mp_limb_t mask;
+
+      /* General case: 2 <= d < p */
+      sh = (-p) % BITS_PER_MP_LIMB;
+      cp = (mp_limb_t*) TMP_ALLOC(n * BYTES_PER_MP_LIMB);
+
+      /* Shift c in temporary allocated place */
+      dm = d % BITS_PER_MP_LIMB;
+      m = d / BITS_PER_MP_LIMB;
+      if (MPFR_UNLIKELY(dm == 0))
+	{
+	  /* dm = 0 and m > 0: Just copy */
+	  MPFR_ASSERTD(m!=0);
+	  MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
+	  MPN_ZERO(cp+n-m, m);
+	}
+      else if (MPFR_LIKELY(m == 0))
+	{
+	  /* dm >=2 and m == 0: just shift */
+	  MPFR_ASSERTD(dm >= 2);
+	  mpn_rshift(cp, MPFR_MANT(c), n, dm);
+	}
+      else
+	{
+	  /* dm > 0 and m > 0: shift and zero  */
+	  mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
+	  MPN_ZERO(cp+n-m, m);
+	}
+
+      DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
+      DEBUG( mpfr_print_mant_binary("B=    ", MPFR_MANT(b), p) );
+      DEBUG( mpfr_print_mant_binary("After ", cp, p) );
+
+      /* Compute bcp=Cp and bcp1=C'p+1 */
+      if (MPFR_LIKELY(sh))
+	{
+	  /* Try to compute them from C' rather than C (FIXME: Faster?) */
+	  bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
+	  if (MPFR_LIKELY((cp[0]&((MPFR_LIMB_ONE<<(sh-1))-MPFR_LIMB_ONE))))
+	    bcp1 = 1;
+	  else
+	    {
+	      /* We can't compute C'p+1 from C'. Compute it from C */
+	      mp_limb_t *tp = MPFR_MANT(c);
+	      /* Start from bit x=p-d+sh in mantissa C
+		 (+sh since we have already looked sh bits in C'!) */
+	      mpfr_prec_t x = p-d+sh;
+	      if (MPFR_UNLIKELY(x>p))
+		/* We are already looked at all the bits of c, so C'p+1 = 0*/
+		bcp1 = 0;
+	      else
+		{
+		  mp_size_t kx = n-1 - (x / BITS_PER_MP_LIMB);
+		  mpfr_prec_t sx = (-x)%BITS_PER_MP_LIMB;
+		  /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
+		  if (tp[kx] & ((MPFR_LIMB_ONE<<sx)-MPFR_LIMB_ONE))
+		    bcp1 = -1;
+		  else
+		    {
+		      kx += (sx==0); /*If sx==0, tp[kx] hasn't been checked*/
+		      do { 
+			kx--;
+		      } while (kx>=0 && tp[kx]==0);
+		      bcp1 = (kx >= 0);
+		    }
+		}
+	    }
+	}
+      else
+	{
+          /* Compute Cp and C'p+1 from C with sh=0 */
+	  mp_limb_t *tp = MPFR_MANT(c);
+	  /* Start from bit x=p-d in mantissa C */
+	  mpfr_prec_t  x = p-d;
+	  mp_size_t   kx = n-1 - (x / BITS_PER_MP_LIMB);
+	  mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB);
+	  MPFR_ASSERTD(p >= d);
+	  bcp = (tp[kx] & (MPFR_LIMB_ONE<<sx));
+	  /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
+	  if (tp[kx] &((MPFR_LIMB_ONE<<sx)-MPFR_LIMB_ONE))
+	    bcp1 = 1;
+	  else
+	    {
+	      kx += (sx==0); /*If sx==0, tp[kx] hasn't been checked*/
+	      do {
+		kx--;
+	      } while (kx>=0 && tp[kx]==0);
+	      bcp1 = (kx>=0);
+	    }
+	}
+      DEBUG( printf("Cp=%d C'p+1=%d\n", bcp, bcp1) );
+
+      /* Check if we can lose a bit, and if so compute Cp+1 and C'p+2 */
+      bp = MPFR_MANT(b);
+      if (MPFR_UNLIKELY((bp[n-1]-cp[n-1]) <= MPFR_LIMB_HIGHBIT))
+	{
+	  /* We can lose a bit so we precompute Cp+1 and C'p+2 */
+	  /* Test for trivial case: since C'p+1=0, Cp+1=0 and C'p+2 =0 */
+	  if (MPFR_UNLIKELY(bcp1 == 0))
+	    {
+	      bbcp = 0;
+	      bbcp1 = 0;
+	    }
+ 	  else /* bcp1 != 0 */
+	    {
+	      /* We can lose a bit: 
+		 compute Cp+1 and C'p+2 from mantissa C */
+	      mp_limb_t *tp = MPFR_MANT(c);
+	      /* Start from bit x=(p+1)-d in mantissa C */
+	      mpfr_prec_t  x = p+1-d;
+	      mp_size_t   kx = n-1 - (x/BITS_PER_MP_LIMB);
+	      mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB);
+	      MPFR_ASSERTD(p > d);
+	      bbcp = (tp[kx] & (MPFR_LIMB_ONE<<sx)) ;
+	      /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
+	      /* If Cp+1=0, since C'p+1!=0, C'p+2=1 ! */
+	      if (!bbcp || tp[kx] &((MPFR_LIMB_ONE<<sx)-MPFR_LIMB_ONE))
+		bbcp1 = 1;
+	      else
+		{
+		  kx += (sx==0); /*If sx==0, tp[kx] hasn't been checked*/
+		  do {
+		    kx--;
+		  } while (kx>=0 && tp[kx]==0);
+		  bbcp1 = (kx>=0);
+		}
+	    } /*End of Bcp1 != 0*/
+	  DEBUG( printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp, bbcp1) );
+	} /* End of "can lose a bit" */
+
+      /* Clean shifted C' */
+      mask = ~((MPFR_LIMB_ONE << sh) - MPFR_LIMB_ONE);
+      cp[0] &= mask;
+
+      /* Substract the mantissa c from b in a */
+      ap = MPFR_MANT(a);
+      mpn_sub_n (ap, bp, cp, n);
+      DEBUG( mpfr_print_mant_binary("Sub=  ", ap, p) );
+
+      /* Normalize: we lose at max one bit*/
+      limb = ap[n-1];
+      if (MPFR_UNLIKELY(limb < ~limb))
+	{
+	  /* High bit is not set and we have to fix it! */
+	  /* Ap >= 010000xxx001 */
+	  mpn_lshift(ap, ap, n, 1);
+	  /* Ap >= 100000xxx010 */
+	  if (bcp) /* Check if Cp = -1 */
+	    /* Since Cp == -1, we have to substract one more */
+	    {
+	      mpn_sub_1(ap, ap, n, MPFR_LIMB_ONE<<sh);
+	      MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
+	    }
+	  /* Ap >= 10000xxx001 */
+	  /* Final exponent -1 since we have shifted the mantissa */
+	  bx--;
+	  /* Update bcp and bcp1 */
+	  bcp  = bbcp; 
+	  bcp1 = bbcp1;
+	  /* We dont't have anymore a valid Cp+1!
+	     But since Ap >= 100000xxx001, the final sub can't unnormalize!*/
+	}
+      MPFR_ASSERTD( !(ap[0] & ~mask) );
+
+      /* Rounding */
+      switch (rnd_mode)
+	{
+	case GMP_RNDN:
+	  if (bcp==0)
+	    goto truncate;
+	  else if ((bcp1) || ((ap[0] & (MPFR_LIMB_ONE<<sh)) != 0))
+	    goto sub_one_ulp;
+	  else
+	    goto truncate;
+	case GMP_RNDZ:
+	  if (MPFR_LIKELY(bcp || bcp1))
+	    goto sub_one_ulp;
+	default:
+	  goto truncate;
+	}      
+ }
+  MPFR_ASSERTN(0);
+
+  /* Sub one ulp to the result */
+ sub_one_ulp:
+  mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
+  /* Result should be smaller than exact value: inexact=-1 */
+  inexact = -1;
+  /* Check normalisation */
+  limb = ap[n-1];
+  if (MPFR_UNLIKELY(limb < ~limb))
+    {
+      /* ap was a power of 2, and we lose a bit */
+      /* Now it is 0111111111111111111[00000 */
+      mpn_lshift(ap, ap, n, 1);
+      bx--;
+      /* And the lost bit x depends on Cp+1, and Cp */
+      /* Compute Cp+1 if it isn't already compute (ie d==1) */
+      if (d == 1)
+	bbcp = 0;
+      DEBUG( printf("(SubOneUlp) Cp=%d, Cp+1=%d C'p+1=%d\n", bcp, bbcp, bcp1));
+      /* Compute the last bit (Since we have shifted the mantissa)
+	 we need one more bit!*/
+      if ((rnd_mode == GMP_RNDZ && bcp==0) 
+	  || (rnd_mode==GMP_RNDN && bbcp ==0))
+	{
+	  ap[0] |= MPFR_LIMB_ONE<<sh;
+	  if (rnd_mode == GMP_RNDN)
+	    inexact = 1;
+	  DEBUG( printf("(SubOneUlp) Last bit set\n") );
+	}
+      /* Result could be exact if C'p+1 = 0 and rnd == Zero 
+	 since we have had one more bit to the result */
+      /* Fixme: rnd_mode == GMP_RNDZ needed ? */
+      if (bcp1==0 && rnd_mode == GMP_RNDZ)
+	{
+	  DEBUG( printf("(SubOneUlp) Exact result\n") );
+	  inexact = 0;
+	}
+    }
+  
+  goto end_of_sub;
+
+ truncate:
+  /* Check if the result is an exact power of 2: 100000000000
+     in which cases, we could have to do sub_one_ulp due to some nasty reasons:
+     If Result is a Power of 2:
+      + If rnd = AWAY, 
+      |  If Cp=-1 and C'p+1 = 0, SubOneUlp and the result is EXACT.
+	 If Cp=-1 and C'p+1 =-1, SubOneUlp and the result is above.
+	 Otherwise truncate
+      + If rnd = NEAREST,
+         If Cp= 0 and Cp+1  =-1 and C'p+2=-1, SubOneUlp and the result is above
+  	 If cp=-1 and C'p+1 = 0, SubOneUlp and the result is exact.
+	 Otherwise truncate.
+      X bit should always be set if SubOneUlp*/
+  if (MPFR_UNLIKELY(ap[n-1] == MPFR_LIMB_HIGHBIT))
+    {
+      mp_size_t k = n-1;
+      do {
+	k--;
+      } while (k>=0 && ap[k]==0);
+      if (MPFR_UNLIKELY(k<0))
+	{
+	  /* It is a power of 2! */
+	  /* Compute Cp+1 if it isn't already compute (ie d==1) */
+	  if (d == 1)
+	    bbcp=0;
+	  DEBUG( printf("(Truncate) Cp=%d, Cp+1=%d C'p+1=%d C'p+2=%d\n", \
+		 bcp, bbcp, bcp1, bbcp1) );
+	  if (( ((rnd_mode == GMP_RNDN) || (rnd_mode != GMP_RNDZ)) && bcp) 
+	      ||
+	      ((rnd_mode == GMP_RNDN) && (bcp == 0) && (bbcp) && (bbcp1)))
+	    {
+	      DEBUG( printf("(Truncate) Do sub\n") );
+	      mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
+	      mpn_lshift(ap, ap, n, 1);
+	      ap[0] |= MPFR_LIMB_ONE<<sh;
+	      bx--;
+	      /* FIXME: Explain why it works (or why not)... */
+	      inexact = (bcp1 == 0) ? 0 : (rnd_mode==GMP_RNDN) ? -1 : 1;
+	      goto end_of_sub;
+	    }
+	}
+    }
+  
+  /* Calcul of Inexact flag.*/
+  inexact = MPFR_LIKELY(bcp || bcp1) ? 1 : 0;
+
+ end_of_sub:
+  /* Update Expo */
+  if (MPFR_UNLIKELY(bx < __gmpfr_emin))
+    {
+      TMP_FREE(marker);
+      DEBUG( printf("(Final Underflow)\n") );
+      /* Update rnd_mode */
+      if (rnd_mode == GMP_RNDN &&
+	  (bx < __gmpfr_emin - 1 ||
+	   (inexact >= 0 && mpfr_powerof2_raw (a))))
+	rnd_mode = GMP_RNDZ;
+      return mpfr_set_underflow (a, rnd_mode, MPFR_SIGN(a));
+    }
+  MPFR_SET_EXP (a, bx);
+
+  TMP_FREE(marker);
+  return inexact*MPFR_INT_SIGN(a);
+}
+
diff --git a/tests/Makefile.am b/tests/Makefile.am
index a67f8afd3..03a4800f1 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,6 +1,6 @@
 AUTOMAKE_OPTIONS = gnu
 
-check_PROGRAMS = tinits tsgn tcheck tisnan texceptions tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tadd tsub tmul tdiv tcmp tcmp2 tcan_round tround_prec tswap reuse tsqrt tnext tcomparisons teq tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tcmp_d tmul_2exp tfma tfrac tget_str tout_str tget_d tget_d_2exp tconst_log2 tconst_pi tconst_euler trandom ttrunc trint texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt mpf_compat mpfr_compat tzeta 
+check_PROGRAMS = tinits tsgn tcheck tisnan texceptions tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tadd tsub tmul tdiv tcmp tcmp2 tcan_round tround_prec tswap reuse tsqrt tnext tcomparisons teq tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tcmp_d tmul_2exp tfma tfrac tget_str tout_str tget_d tget_d_2exp tconst_log2 tconst_pi tconst_euler trandom ttrunc trint texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt mpf_compat mpfr_compat tzeta tsub1sp
 
 EXTRA_DIST = tgeneric.c mpf_compat.h
 
diff --git a/tests/tabs.c b/tests/tabs.c
index 46d4fcdb4..d26891bd8 100644
--- a/tests/tabs.c
+++ b/tests/tabs.c
@@ -96,7 +96,7 @@ main (int argc, char *argv[])
 
   mpfr_set_ui(x, 1, GMP_RNDN);
   mpfr_abs(x, x, GMP_RNDN);
-  if (mpfr_get_d1 (x) != 1.0)
+  if (mpfr_cmp_ui (x, 1))
     {
       printf ("Error in mpfr_abs(1.0)\n");
       exit (1);
@@ -104,7 +104,7 @@ main (int argc, char *argv[])
 
   mpfr_set_si(x, -1, GMP_RNDN);
   mpfr_abs(x, x, GMP_RNDN);
-  if (mpfr_get_d1 (x) != 1.0)
+  if (mpfr_cmp_ui (x, 1))
     {
       printf ("Error in mpfr_abs(-1.0)\n");
       exit (1);
diff --git a/tests/tadd.c b/tests/tadd.c
index 8faf2f62d..c910cc264 100644
--- a/tests/tadd.c
+++ b/tests/tadd.c
@@ -237,7 +237,7 @@ check64 (void)
   mpfr_set_str_binary(x, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
   mpfr_set_str_binary(t, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
   mpfr_sub(u, x, t, GMP_RNDU);
-  if (mpfr_get_d1 (u) != 9.4349060620538533806e167)
+  if (mpfr_cmp_ui_2exp(u, 1, 558))
     { /* 2^558 */
       printf ("Error (1) in mpfr_sub\n");
       exit (1);
@@ -427,7 +427,7 @@ check_same (void)
 
   mpfr_init(x); mpfr_set_ui(x, 1, GMP_RNDZ);
   mpfr_add(x, x, x, GMP_RNDZ);
-  if (mpfr_get_d1 (x) != 2.0)
+  if (mpfr_cmp_ui (x, 2)) 
     {
       printf ("Error when all 3 operands are equal\n");
       exit (1);
diff --git a/tests/tconst_log2.c b/tests/tconst_log2.c
index 450db9f80..71695d52d 100644
--- a/tests/tconst_log2.c
+++ b/tests/tconst_log2.c
@@ -88,7 +88,7 @@ main (int argc, char *argv[])
   /* check precision of 2 bits */
   mpfr_set_prec (x, 2);
   mpfr_const_log2 (x, GMP_RNDN);
-  if (mpfr_get_d1 (x) != 0.75)
+  if (mpfr_cmp_ui_2exp(x, 3, -2)) /* 3*2^-2 */
     {
       printf ("mpfr_const_log2 failed for prec=2, rnd=GMP_RNDN\n"
               "expected 0.75, got %f\n", mpfr_get_d1 (x));
diff --git a/tests/tdiv_ui.c b/tests/tdiv_ui.c
index 8e468bdfa..7796fb388 100644
--- a/tests/tdiv_ui.c
+++ b/tests/tdiv_ui.c
@@ -103,7 +103,7 @@ special (void)
 	{
 	  mpfr_set_prec (y, yprec);
 	  mpfr_div_ui (y, x, 1, GMP_RNDN);
-	  if (mpfr_get_d1 (x) != mpfr_get_d1 (y))
+	  if (mpfr_cmp(x,y))
 	    {
 	      printf ("division by 1.0 fails for xprec=%u, yprec=%u\n", xprec, yprec);
 	      printf ("expected "); mpfr_print_binary (x); puts ("");
diff --git a/tests/tmul.c b/tests/tmul.c
index dc6d05392..328326312 100644
--- a/tests/tmul.c
+++ b/tests/tmul.c
@@ -170,7 +170,7 @@ check_sign (void)
   mpfr_set_si (a, -1, GMP_RNDN);
   mpfr_set_ui (b, 2, GMP_RNDN);
   mpfr_mul(a, b, b, GMP_RNDN);
-  if (mpfr_get_d1 (a) != 4.0)
+  if (mpfr_cmp_ui (a, 4) )
     {
       printf ("2.0*2.0 gives %1.20e\n", mpfr_get_d1 (a));
       exit (1);
diff --git a/tests/tmul_ui.c b/tests/tmul_ui.c
index ad54b9af5..f8defca92 100644
--- a/tests/tmul_ui.c
+++ b/tests/tmul_ui.c
@@ -136,7 +136,7 @@ main (int argc, char *argv[])
   mpfr_set_str (x, /*1.0/3.0*/ 
 		"0.333333333333333333333333333333333", 10, GMP_RNDZ);
   mpfr_mul_ui (x, x, 3, GMP_RNDU);
-  if (mpfr_get_d1 (x) != 1.0)
+  if (mpfr_cmp_ui (x, 1))
     {
       printf ("Error in mpfr_mul_ui: U(Z(1/3)*3) does not give 1\n");
       exit (1);
@@ -217,7 +217,7 @@ main (int argc, char *argv[])
         {
           mpfr_set_prec (y, yprec);
           mpfr_mul_ui (y, x, 1, GMP_RNDN);
-          if (mpfr_get_d1 (x) != mpfr_get_d1 (y))
+          if (mpfr_cmp(x,y))
             {
               printf ("multiplication by 1.0 fails for xprec=%u, yprec=%u\n",
                       xprec, yprec);
diff --git a/tests/tpow.c b/tests/tpow.c
index 0adb8a4b2..beea9be20 100644
--- a/tests/tpow.c
+++ b/tests/tpow.c
@@ -165,7 +165,7 @@ special ()
   mpfr_set_str_binary (x, "0.111011000111100000111010000101010100110011010000011");
   mpfr_set_str_binary (y, "0.111110010100110000011101100011010111000010000100101");
   mpfr_set_str_binary (t, "0.1110110011110110001000110100100001001111010011111000010000011001");
-    ;
+  
   mpfr_pow (z, x, y, GMP_RNDN);
   if (mpfr_cmp (z, t))
     {
@@ -179,7 +179,7 @@ special ()
   mpfr_set_str (x, "5.68824667828621954868e-01", 10, GMP_RNDN);
   mpfr_set_str (y, "9.03327850535952658895e-01", 10, GMP_RNDN);
   mpfr_pow (z, x, y, GMP_RNDZ);
-  if (mpfr_get_d1 (z) != 0.60071044650456473235)
+  if (mpfr_cmp_d(z, 0.60071044650456473235))
     {
       printf ("Error in mpfr_pow for prec=53, rnd=GMP_RNDZ\n");
       exit (1);
@@ -365,7 +365,7 @@ main (void)
   for (p=2; p<100; p++)
     check_inexact (p);
 
-  /* underflows (); */
+  /*underflows ();*/
 
   overflows ();
 
diff --git a/tests/tset_d.c b/tests/tset_d.c
index 6d90326db..af44a8a6b 100644
--- a/tests/tset_d.c
+++ b/tests/tset_d.c
@@ -62,20 +62,20 @@ main (int argc, char *argv[])
    /* checks that rounds to nearest sets the last
      bit to zero in case of equal distance */
    mpfr_set_d (x, 5.0, GMP_RNDN);
-   if (mpfr_get_d1 (x) != 4.0)
+   if (mpfr_cmp_ui (x, 4))
      {
        printf ("Error in tset_d: got %g instead of 4.0\n", mpfr_get_d1 (x));
        exit (1);
      }
    mpfr_set_d (x, -5.0, GMP_RNDN);
-   if (mpfr_get_d1 (x) != -4.0)
+   if (mpfr_cmp_si (x, -4))
      {
        printf ("Error in tset_d: got %g instead of -4.0\n", mpfr_get_d1 (x));
        exit (1);
      }
 
    mpfr_set_d (x, 9.84891017624509146344e-01, GMP_RNDU);
-   if (mpfr_get_d1 (x) != 1.0)
+   if (mpfr_cmp_ui (x, 1))
      {
        printf ("Error in tset_d: got %g instead of 1.0\n", mpfr_get_d1 (x));
        exit (1);
@@ -83,7 +83,7 @@ main (int argc, char *argv[])
 
   mpfr_init2(z, 32);
   mpfr_set_d(z, 1.0, 0);
-  if (mpfr_get_d1 (z) != 1.0)
+  if (mpfr_cmp_ui (z, 1))
     {
       mpfr_print_binary (z); puts ("");
       printf ("Error: 1.0 != 1.0\n");
diff --git a/tests/tset_str.c b/tests/tset_str.c
index b9aa25b2b..13f9cfeb7 100644
--- a/tests/tset_str.c
+++ b/tests/tset_str.c
@@ -120,7 +120,7 @@ main (int argc, char *argv[])
   mpfr_set_str_binary (x, "+110101100.01010000101101000000100111001000101011101110E00");
 
   mpfr_set_str_binary (x, "1.0");
-  if (mpfr_get_d1 (x) != 1.0)
+  if (mpfr_cmp_ui (x, 1))
     {
       printf ("Error in mpfr_set_str_binary for s=1.0\n");
       mpfr_clear(x);
@@ -131,7 +131,7 @@ main (int argc, char *argv[])
   mpfr_set_str_binary (x, "+0000");
   mpfr_set_str_binary (x, "+0000E0");
   mpfr_set_str_binary (x, "0000E0");
-  if (mpfr_get_d1 (x) != 0.0)
+  if (mpfr_cmp_ui (x, 0))
     {
       printf ("Error in mpfr_set_str_binary for s=0.0\n");
       mpfr_clear (x);
diff --git a/tests/tsin.c b/tests/tsin.c
index ba5bc0409..33cd4d938 100644
--- a/tests/tsin.c
+++ b/tests/tsin.c
@@ -152,7 +152,7 @@ main (int argc, char *argv[])
 
   mpfr_set_str (x, "0.5", 10, GMP_RNDN);
   mpfr_sin (x, x, GMP_RNDD);
-  if (mpfr_get_d1 (x) != 0.375)
+  if (mpfr_cmp_ui_2exp (x, 3, -3)) /* x != 0.375 = 3/8 */
     {
       printf ("mpfr_sin(0.5, GMP_RNDD) failed with precision=2\n");
       exit (1);
diff --git a/tests/tsqrt.c b/tests/tsqrt.c
index 6d319ee2e..c100a37b4 100644
--- a/tests/tsqrt.c
+++ b/tests/tsqrt.c
@@ -172,7 +172,7 @@ special (void)
       mpfr_set_ui (z, 1, GMP_RNDN);
       mpfr_add_one_ulp (z, GMP_RNDN);
       mpfr_sqrt (x, z, GMP_RNDU);
-      if (mpfr_get_d1 (x) != 1.5)
+      if (mpfr_cmp_ui_2exp(x, 3, -1))
 	{
 	  printf ("Error: sqrt(1+ulp(1), up) should give 1.5 (prec=%u)\n",
                   (unsigned int) p);
diff --git a/tests/tsub.c b/tests/tsub.c
index 1d54f1c23..2126d6cb2 100644
--- a/tests/tsub.c
+++ b/tests/tsub.c
@@ -161,9 +161,9 @@ check_diverse (void)
   /* check in-place operations */
   mpfr_set_ui (x, 1, GMP_RNDN);
   mpfr_sub (x, x, x, GMP_RNDN);
-  if (mpfr_get_d1 (x) != 0.0)
+  if (mpfr_cmp_ui(x, 0))
     {
-      printf ("Error for mpfr_add (x, x, x, GMP_RNDN) with x=1.0\n");
+      printf ("Error for mpfr_sub (x, x, x, GMP_RNDN) with x=1.0\n");
       exit (1);
     }
 
diff --git a/tests/tsub1sp.c b/tests/tsub1sp.c
new file mode 100644
index 000000000..66ee6af64
--- /dev/null
+++ b/tests/tsub1sp.c
@@ -0,0 +1,365 @@
+/* Test file for mpfr_sub1sp.
+
+Copyright 2003 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the MPFR Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "mpfr-impl.h"
+#include "mpfr-test.h"
+
+void check_special(void);
+void check_random(mpfr_prec_t p);
+
+int main(void)
+{
+  mpfr_prec_t p;
+
+  tests_start_mpfr ();
+
+  check_special ();
+  for(p = 2 ; p < 200 ; p++)
+    check_random (p);
+
+  tests_end_mpfr ();
+  return 0;
+}
+
+#define STD_ERROR \
+            {\
+              printf("ERROR: for %s and p=%lu and i=%d:\nY=",\
+		      mpfr_print_rnd_mode(r), p, i);\
+	       mpfr_print_binary(y);\
+	       printf("\nZ="); mpfr_print_binary(z);\
+	       printf("\nReal: "); mpfr_print_binary(x2);\
+	       printf("\nGot : "); mpfr_print_binary(x);\
+	       putchar('\n');\
+               abort();\
+            }
+
+#define STD_ERROR2 \
+            {\
+              printf("ERROR: for %s and p=%lu and i=%d:\nY=",\
+                      mpfr_print_rnd_mode(r), p, i);\
+               mpfr_print_binary(y);\
+               printf("\nZ="); mpfr_print_binary(z);\
+               printf("\nR="); mpfr_print_binary(x);\
+               printf("\nWrong inexact flag. Real: %d. Got: %d", \
+		      inexact1, inexact2); \
+               abort();\
+            }
+
+void check_random(mpfr_prec_t p)
+{
+  mpfr_t x,y,z,x2;
+  mp_rnd_t r;
+  int i, inexact1, inexact2;
+  gmp_randstate_t state;
+
+  mpfr_inits2(p, x,y,z,x2,NULL);
+  gmp_randinit_lc_2exp_size (state, 128);
+  gmp_randseed_ui (state, 17422471);
+
+  for(i = 0 ; i < 1000 ; i++)
+    {
+      mpfr_urandomb(y, state);
+      mpfr_urandomb(z, state);
+      if (MPFR_IS_PURE_FP(y) && MPFR_IS_PURE_FP(z))
+	for(r = 0 ; r < 4 ; r++)
+	  {
+	    inexact1 = mpfr_sub1(x2, y, z, r);
+	    inexact2 = mpfr_sub1sp(x, y, z, r);
+	    if (mpfr_cmp(x, x2))
+	      STD_ERROR;
+	    if (inexact1 != inexact2)
+	      STD_ERROR2;
+	  }
+    }
+
+  mpfr_clears(x,y,z,x2,NULL);
+  gmp_randclear(state);
+}
+
+void check_special(void)
+{
+  mpfr_t x,y,z,x2;
+  mp_rnd_t r;
+  mpfr_prec_t p;
+  int i = -1, inexact1, inexact2;
+
+  mpfr_inits(x,y,z,x2,NULL);
+
+  for(r = 0 ; r < 4 ; r++)
+    {
+      p = 53;
+      mpfr_set_prec(x, 53);
+      mpfr_set_prec(x2, 53);
+      mpfr_set_prec(y, 53);
+      mpfr_set_prec(z, 53);
+
+      mpfr_set_str_binary (y, 
+       "0.10110111101101110010010010011011000001101101011011001E31");
+  
+      mpfr_sub1sp(x, y, y, r);
+      if (mpfr_cmp_ui(x, 0))
+	{
+	  printf("Error for x-x with p=%lu. Expected 0. Got:", p);
+	  mpfr_print_binary(x);
+	  abort();
+	}
+      
+      mpfr_set(z, y, r);
+      mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp_ui(x, 0))
+	{
+	  printf("Error for x-y with y=x and p=%lu. Expected 0. Got:", p);
+	  mpfr_print_binary(x);
+	  abort();
+	}
+      /* diff = 0 */
+      mpfr_set_str_binary (y,
+       "0.10110111101101110010010010011011001001101101011011001E31");
+      inexact1 = mpfr_sub1(x2, y, z, r);
+      inexact2 = mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp(x, x2))
+	STD_ERROR;
+      if (inexact1 != inexact2)
+	STD_ERROR2;
+      
+      /* Diff = 1 */
+      mpfr_set_str_binary (y,
+       "0.10110111101101110010010010011011000001101101011011001E30");
+      inexact1 = mpfr_sub1(x2, y, z, r);
+      inexact2 = mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp(x, x2))
+	STD_ERROR;
+      if (inexact1 != inexact2)
+	STD_ERROR2;
+      
+      /* Diff = 2 */
+      mpfr_set_str_binary (y,
+       "0.10110111101101110010010010011011000101101101011011001E32");
+      inexact1 = mpfr_sub1(x2, y, z, r);
+      inexact2 = mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp(x, x2))
+	STD_ERROR;
+      if (inexact1 != inexact2)
+	STD_ERROR2;
+      
+      /* Diff = 32 */
+      mpfr_set_str_binary (y,
+       "0.10110111101101110010010010011011000001101101011011001E63");
+      inexact1 = mpfr_sub1(x2, y, z, r);
+      inexact2 = mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp(x, x2))
+	STD_ERROR;
+      if (inexact1 != inexact2)
+	STD_ERROR2;
+      
+      /* Diff = 52 */
+      mpfr_set_str_binary (y,
+       "0.10110111101101110010010010011011010001101101011011001E83");
+      inexact1 = mpfr_sub1(x2, y, z, r);
+      inexact2 = mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp(x, x2))
+	STD_ERROR;
+      if (inexact1 != inexact2)
+	STD_ERROR2;
+      
+      /* Diff = 53 */
+      mpfr_set_str_binary (y,
+       "0.10110111101101110010010010011111000001101101011011001E31");
+      inexact1 = mpfr_sub1(x2, y, z, r);
+      inexact2 = mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp(x, x2))
+        STD_ERROR;
+      if (inexact1 != inexact2)
+	STD_ERROR2;
+
+      /* Diff > 200 */
+      mpfr_set_str_binary (y,
+       "0.10110111101101110010010010011011000001101101011011001E331");
+      inexact1 = mpfr_sub1(x2, y, z, r);
+      inexact2 = mpfr_sub1sp(x, y, z, r);
+      if (mpfr_cmp(x, x2))
+	STD_ERROR;
+      if (inexact1 != inexact2)
+	STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+       "0.10000000000000000000000000000000000000000000000000000E31");
+     mpfr_set_str_binary (z,
+       "0.11111111111111111111111111111111111111111111111111111E30");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+       "0.10000000000000000000000000000000000000000000000000000E31");
+     mpfr_set_str_binary (z,
+       "0.11111111111111111111111111111111111111111111111111111E29");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+       "0.10000000000000000000000000000000000000000000000000000E52");
+     mpfr_set_str_binary (z,
+       "0.10000000000010000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+	"0.11100000000000000000000000000000000000000000000000000E53");
+     mpfr_set_str_binary (z,
+	"0.10000000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(z, y, z, r);
+     mpfr_set(x, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+     
+     mpfr_set_str_binary (y,
+       "0.10000000000000000000000000000000000000000000000000000E53");
+     mpfr_set_str_binary (z,
+       "0.10100000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR; 
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+	"0.10000000000000000000000000000000000000000000000000000E54");
+     mpfr_set_str_binary (z,
+	"0.10100000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+     
+     p = 63;
+     mpfr_set_prec(x, p);
+     mpfr_set_prec(x2, p);
+     mpfr_set_prec(y, p);
+     mpfr_set_prec(z, p);
+     mpfr_set_str_binary (y,
+      "0.100000000000000000000000000000000000000000000000000000000000000E62");
+     mpfr_set_str_binary (z,
+      "0.110000000000000000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     p = 64;
+     mpfr_set_prec(x, 64);
+     mpfr_set_prec(x2, 64);
+     mpfr_set_prec(y, 64);
+     mpfr_set_prec(z, 64);
+
+     mpfr_set_str_binary (y,
+      "0.1100000000000000000000000000000000000000000000000000000000000000E31");
+     mpfr_set_str_binary (z,
+      "0.1111111111111111111111111110000000000000000000000000011111111111E29");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+      "0.1000000000000000000000000000000000000000000000000000000000000000E63");
+     mpfr_set_str_binary (z,
+      "0.1011000000000000000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+     
+     mpfr_set_str_binary (y,
+      "0.1000000000000000000000000000000000000000000000000000000000000000E63");
+     mpfr_set_str_binary (z,
+      "0.1110000000000000000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+        "0.10000000000000000000000000000000000000000000000000000000000000E63");
+     mpfr_set_str_binary (z,
+        "0.10000000000000000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+      "0.1000000000000000000000000000000000000000000000000000000000000000E64");
+     mpfr_set_str_binary (z,
+      "0.1010000000000000000000000000000000000000000000000000000000000000E00");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+     mpfr_set_str_binary (y,
+      "0.1000000000000000000000000000000000000000000000000000000000000000"
+			  "E-1073741823");
+     mpfr_set_str_binary (z,
+      "0.1010000000000000000000000000000000000000000000000000000000000000"
+			  "E-1073741823");
+     inexact1 = mpfr_sub1(x2, y, z, r);
+     inexact2 = mpfr_sub1sp(x, y, z, r);
+     if (mpfr_cmp(x, x2))
+       STD_ERROR;
+     if (inexact1 != inexact2)
+       STD_ERROR2;
+
+    }
+
+  mpfr_clears(x,y,z,x2,NULL);  
+}
diff --git a/tests/ttan.c b/tests/ttan.c
index 9b15ea9f6..8f435037a 100644
--- a/tests/ttan.c
+++ b/tests/ttan.c
@@ -82,7 +82,7 @@ main(int argc, char *argv[])
   mpfr_set_prec (x, 2);
   mpfr_set_str (x, "0.5", 10, GMP_RNDN);
   mpfr_tan (x, x, GMP_RNDD);
-  if (mpfr_get_d1 (x) != 0.5)
+  if (mpfr_cmp_ui_2exp(x, 1, -1))
     {
       printf ("mpfr_tan(0.5, GMP_RNDD) failed\n"
               "expected 0.5, got %f\n", mpfr_get_d1 (x));
diff --git a/tests/ttrunc.c b/tests/ttrunc.c
index 7d5010adc..347f4000c 100644
--- a/tests/ttrunc.c
+++ b/tests/ttrunc.c
@@ -45,7 +45,7 @@ main (void)
 
   mpfr_set_str (x, "0.5", 10, GMP_RNDN);
   mpfr_ceil(y, x);
-  if (mpfr_get_d1 (y) != 1.0)
+  if (mpfr_cmp_ui (y, 1))
     {
       printf ("Error in mpfr_ceil for x=0.5: expected 1.0, got %f\n",
               mpfr_get_d1 (y));
@@ -54,7 +54,7 @@ main (void)
 
   mpfr_set_ui (x, 0, GMP_RNDN);
   mpfr_ceil(y, x);
-  if (mpfr_get_d1 (y) != 0.0)
+  if (mpfr_cmp_ui(y,0))
     {
       printf ("Error in mpfr_ceil for x=0.0: expected 0.0, got %f\n",
               mpfr_get_d1 (y));
@@ -63,7 +63,7 @@ main (void)
 
   mpfr_set_ui (x, 1, GMP_RNDN);
   mpfr_ceil(y, x);
-  if (mpfr_get_d1 (y) != 1.0)
+  if (mpfr_cmp_ui(y,1))
     {
       printf ("Error in mpfr_ceil for x=1.0: expected 1.0, got %f\n",
               mpfr_get_d1 (y));
diff --git a/tests/tui_pow.c b/tests/tui_pow.c
index cc9af8fef..d072c2810 100644
--- a/tests/tui_pow.c
+++ b/tests/tui_pow.c
@@ -182,7 +182,7 @@ main (int argc, char *argv[])
   mpfr_set_prec (y, 2);
   mpfr_set_str (x, "-0.5", 10, GMP_RNDZ);
   mpfr_ui_pow (y, 4, x, GMP_RNDD);
-  if (mpfr_get_d1 (y) != 0.5)
+  if (mpfr_cmp_ui_2exp(y, 1, -1))
     {
       fprintf (stderr, "Error for 4^(-0.5), prec=2, GMP_RNDD\n");
       fprintf (stderr, "expected 0.5, got ");
@@ -197,7 +197,7 @@ main (int argc, char *argv[])
   mpfr_set_prec (y, 2);
   mpfr_set_str (x, "0.5", 10, GMP_RNDN);
   mpfr_ui_pow (y, 398441521, x, GMP_RNDN);
-  if (mpfr_get_d1 (y) != 16384.0)
+  if (mpfr_cmp_ui_2exp(y, 1, 14))
     {
       fprintf (stderr, "Error for 398441521^(0.5), prec=2, GMP_RNDN\n");
       fprintf (stderr, "expected 1.0e14, got ");
@@ -211,8 +211,6 @@ main (int argc, char *argv[])
 
   mpfr_set_prec (x, 2);
   mpfr_set_str (x, "0.5", 10, GMP_RNDN);
-  /*  mpfr_set_ui (x, 1, GMP_RNDN);
-      mpfr_div_2exp (x, x, 1, GMP_RNDN); /* x = 1/2 */
   check1 (x, 2, 398441521, GMP_RNDN);  /* 398441521 = 19961^2 */
 
   /* generic test */