more changes for interface with NTL

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

10 files changed, 278 insertions, 171 deletions

diff --git a/algorithms.bib b/algorithms.bib
index c3f0e2e98..cab724934 100644
--- a/algorithms.bib
+++ b/algorithms.bib
@@ -72,3 +72,13 @@
   year = 	 1994
 }
 
+@Misc{GoSe04,
+  author =	 {Xavier Gourdon and Pascal Sebah},
+  title =	 {The logarithmic Constant: $\log 2$}
+  month =	 jan,
+  year =	 2004,
+  url =         {http://numbers.computation.free.fr/Constants/constants.html}
+}
+
+
+
diff --git a/algorithms.tex b/algorithms.tex
index dc84b3071..86a2cda44 100644
--- a/algorithms.tex
+++ b/algorithms.tex
@@ -3063,7 +3063,7 @@ Hence we have
 This constant is used in the exponential function,
 and in the base $2$ exponential and logarithm.
 
-We use the following formula:
+We use the following formula (formula (30) from \cite{GoSe04}):
 \[ \log 2 = \frac{3}{4} \sum_{n \geq 0} (-1)^n \frac{{n!}^2}{2^n (2n+1)!}. \]
 Let $w$ be the working precision. We take $N = \lfloor w/3 \rfloor + 1$,
 and evaluate exactly using binary spitting:
diff --git a/mpfr-test.h b/mpfr-test.h
index 29e148133..bcf8e78a8 100644
--- a/mpfr-test.h
+++ b/mpfr-test.h
@@ -90,4 +90,66 @@ int mpfr_cmp_str _MPFR_PROTO ((mpfr_srcptr x, const char *, int, mp_rnd_t));
 }
 #endif
 
+/* define CHECK_EXTERNAL if you want to check mpfr against another library
+   with correct rounding. You'll probably have to modify mpfr_print_raw()
+   and/or test_add() below:
+   * mpfr_print_raw() prints each number as "p m e" where p is the precision,
+     m the mantissa (as a binary integer with sign), and e the exponent.
+     The corresponding number is m*2^e. Example: "2 10 -6" represents
+     2*2^(-6) with a precision of 2 bits.
+   * test_add() outputs "b c a" on one line, for each addition a <- b + c.
+     Currently it only prints such a line for rounding to nearest, when
+     the inputs b and c are not NaN and/or Inf.
+*/
+#ifdef CHECK_EXTERNAL
+static void
+mpfr_print_raw (mpfr_srcptr x)
+{
+  printf ("%lu ", MPFR_PREC (x));
+  if (MPFR_IS_NAN (x))
+    {
+      printf ("@NaN@");
+      return;
+    }
+
+  if (MPFR_SIGN (x) < 0)
+    printf ("-");
+
+  if (MPFR_IS_INF (x))
+    printf ("@Inf@");
+  else if (MPFR_IS_ZERO (x))
+    printf ("0 0");
+  else
+    {
+      mp_limb_t *mx;
+      mp_prec_t px;
+      mp_size_t n;
+
+      mx = MPFR_MANT (x);
+      px = MPFR_PREC (x);
+
+      for (n = (px - 1) / BITS_PER_MP_LIMB; ; n--)
+        {
+          mp_limb_t wd, t;
+
+          MPFR_ASSERTN (n >= 0);
+          wd = mx[n];
+          for (t = MPFR_LIMB_HIGHBIT; t != 0; t >>= 1)
+            {
+              printf ((wd & t) == 0 ? "0" : "1");
+              if (--px == 0)
+                {
+                  mp_exp_t ex;
+
+                  ex = MPFR_GET_EXP (x);
+                  MPFR_ASSERTN (ex >= LONG_MIN && ex <= LONG_MAX);
+                  printf (" %ld", (long) ex - (long) MPFR_PREC (x));
+                  return;
+                }
+            }
+        }
+    }
+}
+#endif
+
 #endif
diff --git a/tests/RRTest.c b/tests/RRTest.c
index ad519e783..884b5f46d 100644
--- a/tests/RRTest.c
+++ b/tests/RRTest.c
@@ -68,17 +68,16 @@ int main()
       if (++line % 1000 == 0)
         cout << "line " << line << endl;
       ReadRR (b);
-      ReadRR (c);
+      //      ReadRR (c);
       p = ReadRR (a);
-      AddPrec (d, b, c, p);
+      SqrRootPrec (d, b, p);
       if (d != a)
         {
           cerr << "error at line " << line << " for b="; Output(b);
-          cerr << " c="; Output(c);
+          //          cerr << " c="; Output(c);
           cerr << "expected "; Output(a);
           cerr << "got      "; Output(d);
-          cerr << "prec(d)=" << d.precision() << endl;
-          exit (1);
+          cerr << "prec(d)=" << p << endl;
         }
     }
 }
diff --git a/tests/tadd.c b/tests/tadd.c
index 4f8b9caed..149810424 100644
--- a/tests/tadd.c
+++ b/tests/tadd.c
@@ -32,68 +32,6 @@ MA 02111-1307, USA. */
 
 static int usesp;
 
-/* define CHECK_EXTERNAL if you want to check mpfr against another library
-   with correct rounding. You'll probably have to modify mpfr_print_raw()
-   and/or test_add() below:
-   * mpfr_print_raw() prints each number as "p m e" where p is the precision,
-     m the mantissa (as a binary integer with sign), and e the exponent.
-     The corresponding number is m*2^e. Example: "2 10 -6" represents
-     2*2^(-6) with a precision of 2 bits.
-   * test_add() outputs "b c a" on one line, for each addition a <- b + c.
-     Currently it only prints such a line for rounding to nearest, when
-     the inputs b and c are not NaN and/or Inf.
-*/
-#ifdef CHECK_EXTERNAL
-static void
-mpfr_print_raw (mpfr_srcptr x)
-{
-  printf ("%lu ", MPFR_PREC (x));
-  if (MPFR_IS_NAN (x))
-    {
-      printf ("@NaN@");
-      return;
-    }
-
-  if (MPFR_SIGN (x) < 0)
-    printf ("-");
-
-  if (MPFR_IS_INF (x))
-    printf ("@Inf@");
-  else if (MPFR_IS_ZERO (x))
-    printf ("0 0");
-  else
-    {
-      mp_limb_t *mx;
-      mp_prec_t px;
-      mp_size_t n;
-
-      mx = MPFR_MANT (x);
-      px = MPFR_PREC (x);
-
-      for (n = (px - 1) / BITS_PER_MP_LIMB; ; n--)
-        {
-          mp_limb_t wd, t;
-
-          MPFR_ASSERTN (n >= 0);
-          wd = mx[n];
-          for (t = MPFR_LIMB_HIGHBIT; t != 0; t >>= 1)
-            {
-              printf ((wd & t) == 0 ? "0" : "1");
-              if (--px == 0)
-                {
-                  mp_exp_t ex;
-
-                  ex = MPFR_GET_EXP (x);
-                  MPFR_ASSERTN (ex >= LONG_MIN && ex <= LONG_MAX);
-                  printf (" %ld", (long) ex - (long) MPFR_PREC (x));
-                  return;
-                }
-            }
-        }
-    }
-}
-#endif
-
 static int
 test_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
 {
diff --git a/tests/tconst_log2.c b/tests/tconst_log2.c
index 203331593..202d00381 100644
--- a/tests/tconst_log2.c
+++ b/tests/tconst_log2.c
@@ -76,7 +76,7 @@ check_large (void)
   mpfr_prec_round (y, 25000, GMP_RNDN);
   if (mpfr_cmp (x, y))
     {
-      printf ("const_pi: error for large prec\n");
+      printf ("const_log2: error for large prec\n");
       exit (1);
     }
 
diff --git a/tests/tdiv.c b/tests/tdiv.c
index 75e363b85..1d4c78191 100644
--- a/tests/tdiv.c
+++ b/tests/tdiv.c
@@ -25,6 +25,31 @@ MA 02111-1307, USA. */
 
 #include "mpfr-test.h"
 
+#ifdef CHECK_EXTERNAL
+static int
+test_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
+{
+  int res;
+  int ok = rnd_mode == GMP_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
+  if (ok)
+    {
+      mpfr_print_raw (b);
+      printf (" ");
+      mpfr_print_raw (c);
+    }
+  res = mpfr_div (a, b, c, rnd_mode);
+  if (ok)
+    {
+      printf (" ");
+      mpfr_print_raw (a);
+      printf ("\n");
+    }
+  return res;
+}
+#else
+#define test_div mpfr_div
+#endif
+
 #define check53(n, d, rnd, res) check4(n, d, rnd, 53, res)
 
 /* return 0 iff a and b are of the same sign */
@@ -48,7 +73,7 @@ check4 (const char *Ns, const char *Ds, mp_rnd_t rnd_mode, int p,
   mpfr_inits2 (p, q, n, d, NULL);
   mpfr_set_str1 (n, Ns);
   mpfr_set_str1 (d, Ds);
-  mpfr_div(q, n, d, rnd_mode);
+  test_div(q, n, d, rnd_mode);
   if (mpfr_cmp_str (q, Qs, ((p==53) ? 10 : 2), GMP_RNDN) )
     {
       printf ("mpfr_div failed for n=%s, d=%s, p=%d, rnd_mode=%s\n",
@@ -71,7 +96,7 @@ check24 (const char *Ns, const char *Ds, mp_rnd_t rnd_mode, const char *Qs)
 
   mpfr_set_str1 (n, Ns);
   mpfr_set_str1 (d, Ds);
-  mpfr_div(q, n, d, rnd_mode);
+  test_div(q, n, d, rnd_mode);
   if (mpfr_cmp_str1 (q, Qs) )
     {
       printf ("mpfr_div failed for n=%s, d=%s, prec=24, rnd_mode=%s\n",
@@ -174,7 +199,7 @@ check_64(void)
 
   mpfr_set_str_binary(x, "1.00100100110110101001010010101111000001011100100101010000000000E54");
   mpfr_set_str_binary(y, "1.00000000000000000000000000000000000000000000000000000000000000E584");
-  mpfr_div(z, x, y, GMP_RNDU);
+  test_div(z, x, y, GMP_RNDU);
   if (mpfr_cmp_str (z, "0.1001001001101101010010100101011110000010111001001010100000000000E-529", 2, GMP_RNDN))
     {
       printf("Error for tdiv for GMP_RNDU and p=64\nx=");
@@ -199,13 +224,13 @@ check_convergence (void)
   mpfr_set_str_binary(x, "0.1011111101011010101000001010011111101000011100011101010011111011000011001010000000111100100111110011001010110100100001001000111001E6944");
   mpfr_init2(y, 130);
   mpfr_set_ui(y, 5, GMP_RNDN);
-  mpfr_div(x, x, y, GMP_RNDD); /* exact division */
+  test_div(x, x, y, GMP_RNDD); /* exact division */
 
   mpfr_set_prec(x, 64);
   mpfr_set_prec(y, 64);
   mpfr_set_str_binary(x, "0.10010010011011010100101001010111100000101110010010101E55");
   mpfr_set_str_binary(y, "0.1E585");
-  mpfr_div(x, x, y, GMP_RNDN);
+  test_div(x, x, y, GMP_RNDN);
   mpfr_set_str_binary(y, "0.10010010011011010100101001010111100000101110010010101E-529");
   if (mpfr_cmp (x, y))
     {
@@ -223,7 +248,7 @@ check_convergence (void)
       for (j = 0;j < GMP_RND_MAX; j++)
         {
           mpfr_set_ui (y, 1, GMP_RNDN);
-          mpfr_div (y, x, y, (mp_rnd_t) j);
+          test_div (y, x, y, (mp_rnd_t) j);
           if (mpfr_cmp_ui (y, 1))
             {
               printf ("mpfr_div failed for x=1.0, y=1.0, prec=%d rnd=%s\n",
@@ -300,7 +325,7 @@ check_hard (void)
 		{
 		  for (rnd = 0; rnd < GMP_RND_MAX; rnd++)
 		    {
-		      inex = mpfr_div (q, u, v, (mp_rnd_t) rnd);
+		      inex = test_div (q, u, v, (mp_rnd_t) rnd);
 		      inex2 = get_inexact (q, u, v);
 		      if (inex_cmp (inex, inex2))
 			{
@@ -355,7 +380,7 @@ check_lowr (void)
         }
       while (mpfr_cmp_ui (tmp, 0) == 0);
       mpfr_mul (x, z, tmp, GMP_RNDN); /* exact */
-      c = mpfr_div (z2, x, tmp, GMP_RNDN);
+      c = test_div (z2, x, tmp, GMP_RNDN);
 
       if (c || mpfr_cmp (z2, z))
         {
@@ -378,7 +403,7 @@ check_lowr (void)
         }
       while (mpfr_cmp_ui (tmp, 0) == 0);
       mpfr_mul (x, z, tmp, GMP_RNDN); /* exact */
-      c = mpfr_div (z2, x, tmp, GMP_RNDN);
+      c = test_div (z2, x, tmp, GMP_RNDN);
       /* since z2 has one less bit that z, either the division is exact
 	 if z is representable on 9 bits, or we have an even round case */
 
@@ -442,8 +467,8 @@ check_lowr (void)
       mpfr_set(y, tmp, GMP_RNDD);
       mpfr_add_one_ulp(x, GMP_RNDN);
 
-      c = mpfr_div(z2, x, y, GMP_RNDD);
-      mpfr_div(z3, x, y, GMP_RNDD);
+      c = test_div(z2, x, y, GMP_RNDD);
+      test_div(z3, x, y, GMP_RNDD);
       mpfr_set(z, z3, GMP_RNDD);
 
       if (c != -1 || mpfr_cmp(z2, z))
@@ -456,9 +481,9 @@ check_lowr (void)
         }
 
       mpfr_set (y, tmp, GMP_RNDU);
-      mpfr_div (z3, x, y, GMP_RNDU);
+      test_div (z3, x, y, GMP_RNDU);
       mpfr_set (z, z3, GMP_RNDU);
-      c = mpfr_div (z2, x, y, GMP_RNDU);
+      c = test_div (z2, x, y, GMP_RNDU);
       if (c != 1 || mpfr_cmp (z2, z))
         {
           printf ("Error in mpfr_div rnd=GMP_RNDU\n");
@@ -499,7 +524,7 @@ check_inexact (void)
   mpfr_set_prec (u, 127);
   mpfr_set_str_binary (u, "0.1000001100110110110101110110101101111000110000001111111110000000011111001010110100110010111111111101000001011011101011101101000E-2");
   mpfr_set_prec (y, 95);
-  inexact = mpfr_div (y, x, u, GMP_RNDN);
+  inexact = test_div (y, x, u, GMP_RNDN);
   if (inexact != (cmp = get_inexact (y, x, u)))
     {
       printf ("Wrong inexact flag (0): expected %d, got %d\n", cmp, inexact);
@@ -514,7 +539,7 @@ check_inexact (void)
   mpfr_set_prec (u, 2);
   mpfr_set_str_binary (u, "0.1E0");
   mpfr_set_prec (y, 28);
-  if ((inexact = mpfr_div (y, x, u, GMP_RNDN) >= 0))
+  if ((inexact = test_div (y, x, u, GMP_RNDN) >= 0))
     {
       printf ("Wrong inexact flag (1): expected -1, got %d\n",
               inexact);
@@ -526,7 +551,7 @@ check_inexact (void)
   mpfr_set_prec (u, 15);
   mpfr_set_str_binary (u, "0.101101000001100E-1");
   mpfr_set_prec (y, 92);
-  if ((inexact = mpfr_div (y, x, u, GMP_RNDN)) <= 0)
+  if ((inexact = test_div (y, x, u, GMP_RNDN)) <= 0)
     {
       printf ("Wrong inexact flag for rnd=GMP_RNDN(1): expected 1, got %d\n",
               inexact);
@@ -550,7 +575,7 @@ check_inexact (void)
 	      mpfr_set_prec (z, py + pu);
 		{
                   rnd = (mp_rnd_t) RND_RAND ();
-		  inexact = mpfr_div (y, x, u, rnd);
+		  inexact = test_div (y, x, u, rnd);
 		  if (mpfr_mul (z, y, u, rnd))
 		    {
 		      printf ("z <- y * u should be exact\n");
@@ -594,13 +619,13 @@ check_nan (void)
   /* 1/nan == nan */
   mpfr_set_ui (a, 1L, GMP_RNDN);
   MPFR_SET_NAN (d);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_nan_p (q));
 
   /* nan/1 == nan */
   MPFR_SET_NAN (a);
   mpfr_set_ui (d, 1L, GMP_RNDN);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_nan_p (q));
 
   /* +inf/1 == +inf */
@@ -608,7 +633,7 @@ check_nan (void)
   MPFR_SET_INF (a);
   MPFR_SET_POS (a);
   mpfr_set_ui (d, 1L, GMP_RNDN);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_inf_p (q));
   MPFR_ASSERTN (mpfr_sgn (q) > 0);
 
@@ -617,14 +642,14 @@ check_nan (void)
   MPFR_CLEAR_FLAGS (d);
   MPFR_SET_INF (d);
   MPFR_SET_POS (d);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_number_p (q));
   MPFR_ASSERTN (mpfr_sgn (q) == 0);
 
   /* 0/0 == nan */
   mpfr_set_ui (a, 0L, GMP_RNDN);
   mpfr_set_ui (d, 0L, GMP_RNDN);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_nan_p (q));
 
   /* +inf/+inf == nan */
@@ -634,33 +659,33 @@ check_nan (void)
   MPFR_CLEAR_FLAGS (d);
   MPFR_SET_INF (d);
   MPFR_SET_POS (d);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_nan_p (q));
 
   /* 1/+0 = +Inf */
   mpfr_set_ui (a, 1, GMP_RNDZ);
   mpfr_set_ui (d, 0, GMP_RNDZ);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
 
   /* 1/-0 = -Inf */
   mpfr_set_ui (a, 1, GMP_RNDZ);
   mpfr_set_ui (d, 0, GMP_RNDZ);
   mpfr_neg (d, d, GMP_RNDZ);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
 
   /* -1/+0 = -Inf */
   mpfr_set_si (a, -1, GMP_RNDZ);
   mpfr_set_ui (d, 0, GMP_RNDZ);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
 
   /* -1/-0 = +Inf */
   mpfr_set_si (a, -1, GMP_RNDZ);
   mpfr_set_ui (d, 0, GMP_RNDZ);
   mpfr_neg (d, d, GMP_RNDZ);
-  MPFR_ASSERTN (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_div (q, a, d, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
 
   /* check overflow */
@@ -669,7 +694,7 @@ check_nan (void)
   mpfr_set_ui (a, 1, GMP_RNDZ);
   mpfr_set_ui (d, 1, GMP_RNDZ);
   mpfr_div_2exp (d, d, 1, GMP_RNDZ);
-  mpfr_div (q, a, d, GMP_RNDU); /* 1 / 0.5 = 2 -> overflow */
+  test_div (q, a, d, GMP_RNDU); /* 1 / 0.5 = 2 -> overflow */
   MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
   set_emax (emax);
 
@@ -679,9 +704,9 @@ check_nan (void)
   mpfr_set_ui (a, 1, GMP_RNDZ);
   mpfr_div_2exp (a, a, 2, GMP_RNDZ);
   mpfr_set_ui (d, 2, GMP_RNDZ);
-  mpfr_div (q, a, d, GMP_RNDZ); /* 0.5*2^(-2) -> underflow */
+  test_div (q, a, d, GMP_RNDZ); /* 0.5*2^(-2) -> underflow */
   MPFR_ASSERTN (mpfr_cmp_ui (q, 0) == 0 && MPFR_IS_POS (q));
-  mpfr_div (q, a, d, GMP_RNDN); /* 0.5*2^(-2) -> underflow */
+  test_div (q, a, d, GMP_RNDN); /* 0.5*2^(-2) -> underflow */
   MPFR_ASSERTN (mpfr_cmp_ui (q, 0) == 0 && MPFR_IS_POS (q));
   set_emin (emin);
 
diff --git a/tests/tmul.c b/tests/tmul.c
index 55502c594..c1a82102e 100644
--- a/tests/tmul.c
+++ b/tests/tmul.c
@@ -24,6 +24,31 @@ MA 02111-1307, USA. */
 
 #include "mpfr-test.h"
 
+#ifdef CHECK_EXTERNAL
+static int
+test_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
+{
+  int res;
+  int ok = rnd_mode == GMP_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
+  if (ok)
+    {
+      mpfr_print_raw (b);
+      printf (" ");
+      mpfr_print_raw (c);
+    }
+  res = mpfr_mul (a, b, c, rnd_mode);
+  if (ok)
+    {
+      printf (" ");
+      mpfr_print_raw (a);
+      printf ("\n");
+    }
+  return res;
+}
+#else
+#define test_mul mpfr_mul
+#endif
+
 /* Workaround for sparc gcc 2.95.x bug, see notes in tadd.c. */
 #define check(x,y,rnd_mode,px,py,pz,res)  pcheck(x,y,res,rnd_mode,px,py,pz)
 
@@ -39,7 +64,7 @@ pcheck (const char *xs, const char *ys, const char *res, mp_rnd_t rnd_mode,
   mpfr_init2 (zz, pz);
   mpfr_set_str1 (xx, xs);
   mpfr_set_str1 (yy, ys);
-  mpfr_mul(zz, xx, yy, rnd_mode);
+  test_mul(zz, xx, yy, rnd_mode);
   if (mpfr_cmp_str1 (zz, res) )
     {
       printf ("(1)mpfr_mul failed for x=%s y=%s with rnd=%s\n",
@@ -70,7 +95,7 @@ check53 (const char *xs, const char *ys, mp_rnd_t rnd_mode, const char *zs)
   mpfr_inits2 (53, xx, yy, zz, NULL);
   mpfr_set_str1 (xx, xs);
   mpfr_set_str1 (yy, ys);
-  mpfr_mul (zz, xx, yy, rnd_mode);
+  test_mul (zz, xx, yy, rnd_mode);
   if (mpfr_cmp_str1 (zz, zs) )
     {
       printf ("(2) mpfr_mul failed for x=%s y=%s with rnd=%s\n",
@@ -102,7 +127,7 @@ check24 (const char *xs, const char *ys, mp_rnd_t rnd_mode, const char *zs)
   mpfr_inits2 (24, xx, yy, zz, NULL);
   mpfr_set_str1 (xx, xs);
   mpfr_set_str1 (yy, ys);
-  mpfr_mul (zz, xx, yy, rnd_mode);
+  test_mul (zz, xx, yy, rnd_mode);
   if (mpfr_cmp_str1 (zz, zs) )
     {
       printf ("(3) mpfr_mul failed for x=%s y=%s with "
@@ -171,7 +196,7 @@ check_sign (void)
   mpfr_init2 (b, 53);
   mpfr_set_si (a, -1, GMP_RNDN);
   mpfr_set_ui (b, 2, GMP_RNDN);
-  mpfr_mul(a, b, b, GMP_RNDN);
+  test_mul(a, b, b, GMP_RNDN);
   if (mpfr_cmp_ui (a, 4) )
     {
       printf ("2.0*2.0 gives \n");
@@ -201,7 +226,7 @@ check_exact (void)
   mpfr_set_prec (c, 32);
   mpfr_set_str_binary (a, "1.1000111011000100e-1");
   mpfr_set_str_binary (b, "1.0010001111100111e-1");
-  if (mpfr_mul (c, a, b, GMP_RNDZ))
+  if (test_mul (c, a, b, GMP_RNDZ))
     {
       printf ("wrong return value (1)\n");
       exit (1);
@@ -218,8 +243,8 @@ check_exact (void)
           mpfr_random (a);
           mpfr_random (b);
           rnd = (mp_rnd_t) RND_RAND ();
-          inexact = mpfr_mul (c, a, b, rnd);
-          if (mpfr_mul (d, a, b, rnd)) /* should be always exact */
+          inexact = test_mul (c, a, b, rnd);
+          if (test_mul (d, a, b, rnd)) /* should be always exact */
             {
               printf ("unexpected inexact return value\n");
               exit (1);
@@ -267,7 +292,7 @@ check_max(void)
   mpfr_set_str1 (yy, "0.68750");
   mpfr_mul_2si(yy, yy, MPFR_EMAX_DEFAULT - MPFR_EMAX_DEFAULT/2 + 1, GMP_RNDN);
   mpfr_clear_flags();
-  mpfr_mul(zz, xx, yy, GMP_RNDU);
+  test_mul(zz, xx, yy, GMP_RNDU);
   if (!(mpfr_overflow_p() && MPFR_IS_INF(zz)))
     {
       printf("check_max failed (should be an overflow)\n");
@@ -275,7 +300,7 @@ check_max(void)
     }
 
   mpfr_clear_flags();
-  mpfr_mul(zz, xx, yy, GMP_RNDD);
+  test_mul(zz, xx, yy, GMP_RNDD);
   if (mpfr_overflow_p() || MPFR_IS_INF(zz))
     {
       printf("check_max failed (should NOT be an overflow)\n");
@@ -303,7 +328,7 @@ check_max(void)
   set_emin (0);
   mpfr_set_str_binary (xx, "0.1E0");
   mpfr_set_str_binary (yy, "0.1E0");
-  mpfr_mul (zz, xx, yy, GMP_RNDN);
+  test_mul (zz, xx, yy, GMP_RNDN);
   /* exact result is 0.1E-1, which should round to 0 */
   MPFR_ASSERTN(mpfr_cmp_ui (zz, 0) == 0 && MPFR_IS_POS(zz));
   set_emin (emin);
@@ -315,7 +340,7 @@ check_max(void)
   mpfr_set_str_binary (xx, "0.1E0");
   mpfr_nextabove (xx);
   mpfr_set_str_binary (yy, "0.1E0");
-  mpfr_mul (zz, xx, yy, GMP_RNDN);
+  test_mul (zz, xx, yy, GMP_RNDN);
   /* exact result is just above 0.1E-1, which should round to minfloat */
   MPFR_ASSERTN(mpfr_cmp (zz, yy) == 0);
   set_emin (emin);
@@ -337,7 +362,7 @@ check_min(void)
   mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT/2, GMP_RNDN);
   mpfr_set_str1(yy, "0.9375");
   mpfr_mul_2si(yy, yy, MPFR_EMIN_DEFAULT - MPFR_EMIN_DEFAULT/2 - 1, GMP_RNDN);
-  mpfr_mul(zz, xx, yy, GMP_RNDD);
+  test_mul(zz, xx, yy, GMP_RNDD);
   if (mpfr_sgn(zz) != 0)
     {
       printf("check_min failed: got ");
@@ -346,7 +371,7 @@ check_min(void)
       exit(1);
     }
 
-  mpfr_mul(zz, xx, yy, GMP_RNDU);
+  test_mul(zz, xx, yy, GMP_RNDU);
   mpfr_set_str1 (xx, "0.5");
   mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT, GMP_RNDN);
   if (mpfr_sgn(xx) <= 0)
@@ -381,32 +406,32 @@ check_nans (void)
   /* nan * 0 == nan */
   mpfr_set_nan (x);
   mpfr_set_ui (y, 0L, GMP_RNDN);
-  mpfr_mul (p, x, y, GMP_RNDN);
+  test_mul (p, x, y, GMP_RNDN);
   MPFR_ASSERTN (mpfr_nan_p (p));
 
   /* 1 * nan == nan */
   mpfr_set_ui (x, 1L, GMP_RNDN);
   mpfr_set_nan (y);
-  mpfr_mul (p, x, y, GMP_RNDN);
+  test_mul (p, x, y, GMP_RNDN);
   MPFR_ASSERTN (mpfr_nan_p (p));
 
   /* 0 * +inf == nan */
   mpfr_set_ui (x, 0L, GMP_RNDN);
   mpfr_set_nan (y);
-  mpfr_mul (p, x, y, GMP_RNDN);
+  test_mul (p, x, y, GMP_RNDN);
   MPFR_ASSERTN (mpfr_nan_p (p));
 
   /* +1 * +inf == +inf */
   mpfr_set_ui (x, 1L, GMP_RNDN);
   mpfr_set_inf (y, 1);
-  mpfr_mul (p, x, y, GMP_RNDN);
+  test_mul (p, x, y, GMP_RNDN);
   MPFR_ASSERTN (mpfr_inf_p (p));
   MPFR_ASSERTN (mpfr_sgn (p) > 0);
 
   /* -1 * +inf == -inf */
   mpfr_set_si (x, -1L, GMP_RNDN);
   mpfr_set_inf (y, 1);
-  mpfr_mul (p, x, y, GMP_RNDN);
+  test_mul (p, x, y, GMP_RNDN);
   MPFR_ASSERTN (mpfr_inf_p (p));
   MPFR_ASSERTN (mpfr_sgn (p) < 0);
 
diff --git a/tests/tsqrt.c b/tests/tsqrt.c
index 859059794..d0a1d7fbd 100644
--- a/tests/tsqrt.c
+++ b/tests/tsqrt.c
@@ -24,6 +24,29 @@ MA 02111-1307, USA. */
 
 #include "mpfr-test.h"
 
+#ifdef CHECK_EXTERNAL
+static int
+test_sqrt (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode)
+{
+  int res;
+  int ok = rnd_mode == GMP_RNDN && mpfr_number_p (b);
+  if (ok)
+    {
+      mpfr_print_raw (b);
+    }
+  res = mpfr_sqrt (a, b, rnd_mode);
+  if (ok)
+    {
+      printf (" ");
+      mpfr_print_raw (a);
+      printf ("\n");
+    }
+  return res;
+}
+#else
+#define test_sqrt mpfr_sqrt
+#endif
+
 static void
 check3 (const char *as, mp_rnd_t rnd_mode, const char *qs)
 {
@@ -31,7 +54,7 @@ check3 (const char *as, mp_rnd_t rnd_mode, const char *qs)
 
   mpfr_init2 (q, 53);
   mpfr_set_str1 (q, as);
-  mpfr_sqrt (q, q, rnd_mode);
+  test_sqrt (q, q, rnd_mode);
   if (mpfr_cmp_str1 (q, qs) )
     {
       printf ("mpfr_sqrt failed for a=%s, rnd_mode=%s\n",
@@ -51,7 +74,7 @@ check4 (const char *as, mp_rnd_t rnd_mode, const char *Qs)
 
   mpfr_init2(q, 53);
   mpfr_set_str1 (q, as);
-  mpfr_sqrt(q, q, rnd_mode);
+  test_sqrt(q, q, rnd_mode);
   if (mpfr_cmp_str (q, Qs, 16, GMP_RNDN))
     {
       printf("mpfr_sqrt failed for a=%s, rnd_mode=%s\n",
@@ -72,7 +95,7 @@ check24 (const char *as, mp_rnd_t rnd_mode, const char *qs)
 
   mpfr_init2(q, 24);
   mpfr_set_str1(q, as);
-  mpfr_sqrt(q, q, rnd_mode);
+  test_sqrt(q, q, rnd_mode);
   if (mpfr_cmp_str1 (q, qs))
     {
       printf("mpfr_sqrt failed for a=%s, prec=24, rnd_mode=%s\n",
@@ -148,7 +171,7 @@ special (void)
   mpfr_set_prec (x, 64);
   mpfr_set_str_binary (x, "1010000010100011011001010101010010001100001101011101110001011001E-1");
   mpfr_set_prec (y, 32);
-  mpfr_sqrt (y, x, GMP_RNDN);
+  test_sqrt (y, x, GMP_RNDN);
   if (mpfr_cmp_ui (y, 2405743844UL))
     {
       printf ("Error for n^2+n+1/2 with n=2405743843\n");
@@ -158,7 +181,7 @@ special (void)
   mpfr_set_prec (x, 65);
   mpfr_set_str_binary (x, "10100000101000110110010101010100100011000011010111011100010110001E-2");
   mpfr_set_prec (y, 32);
-  mpfr_sqrt (y, x, GMP_RNDN);
+  test_sqrt (y, x, GMP_RNDN);
   if (mpfr_cmp_ui (y, 2405743844UL))
     {
       printf ("Error for n^2+n+1/4 with n=2405743843\n");
@@ -169,7 +192,7 @@ special (void)
   mpfr_set_prec (x, 66);
   mpfr_set_str_binary (x, "101000001010001101100101010101001000110000110101110111000101100011E-3");
   mpfr_set_prec (y, 32);
-  mpfr_sqrt (y, x, GMP_RNDN);
+  test_sqrt (y, x, GMP_RNDN);
   if (mpfr_cmp_ui (y, 2405743844UL))
     {
       printf ("Error for n^2+n+1/4+1/8 with n=2405743843\n");
@@ -180,7 +203,7 @@ special (void)
   mpfr_set_prec (x, 66);
   mpfr_set_str_binary (x, "101000001010001101100101010101001000110000110101110111000101100001E-3");
   mpfr_set_prec (y, 32);
-  mpfr_sqrt (y, x, GMP_RNDN);
+  test_sqrt (y, x, GMP_RNDN);
   if (mpfr_cmp_ui (y, 2405743843UL))
     {
       printf ("Error for n^2+n+1/8 with n=2405743843\n");
@@ -190,7 +213,7 @@ special (void)
 
   mpfr_set_prec (x, 27);
   mpfr_set_str_binary (x, "0.110100111010101000010001011");
-  if ((inexact = mpfr_sqrt (x, x, GMP_RNDZ)) >= 0)
+  if ((inexact = test_sqrt (x, x, GMP_RNDZ)) >= 0)
     {
       printf ("Wrong inexact flag: expected -1, got %d\n", inexact);
       exit (1);
@@ -202,7 +225,7 @@ special (void)
       mpfr_set_prec (z, p);
       mpfr_set_ui (z, 1, GMP_RNDN);
       mpfr_add_one_ulp (z, GMP_RNDN);
-      mpfr_sqrt (x, z, GMP_RNDU);
+      test_sqrt (x, z, GMP_RNDU);
       if (mpfr_cmp_ui_2exp(x, 3, -1))
 	{
 	  printf ("Error: sqrt(1+ulp(1), up) should give 1.5 (prec=%u)\n",
@@ -215,7 +238,7 @@ special (void)
   /* check inexact flag */
   mpfr_set_prec (x, 5);
   mpfr_set_str_binary (x, "1.1001E-2");
-  if ((inexact = mpfr_sqrt (x, x, GMP_RNDN)))
+  if ((inexact = test_sqrt (x, x, GMP_RNDN)))
     {
       printf ("Wrong inexact flag: expected 0, got %d\n", inexact);
       exit (1);
@@ -227,7 +250,7 @@ special (void)
   /* checks the sign is correctly set */
   mpfr_set_si (x, 1,  GMP_RNDN);
   mpfr_set_si (z, -1, GMP_RNDN);
-  mpfr_sqrt (z, x, GMP_RNDN);
+  test_sqrt (z, x, GMP_RNDN);
   if (mpfr_cmp_ui (z, 0) < 0)
     {
       printf ("Error: square root of 1 gives ");
@@ -240,13 +263,13 @@ special (void)
   mpfr_set_prec (z, 160);
   mpfr_set_str_binary (z, "0.1011010100000100100100100110011001011100100100000011000111011001011101101101110000110100001000100001100001011000E1");
   mpfr_set_prec (x, 160);
-  mpfr_sqrt(x, z, GMP_RNDN);
-  mpfr_sqrt(z, x, GMP_RNDN);
+  test_sqrt(x, z, GMP_RNDN);
+  test_sqrt(z, x, GMP_RNDN);
 
   mpfr_set_prec (x, 53);
   mpfr_set_str (x, "8093416094703476.0", 10, GMP_RNDN);
   mpfr_div_2exp (x, x, 1075, GMP_RNDN);
-  mpfr_sqrt (x, x, GMP_RNDN);
+  test_sqrt (x, x, GMP_RNDN);
   mpfr_set_str (z, "1e55596835b5ef@-141", 16, GMP_RNDN);
   if (mpfr_cmp (x, z))
     {
@@ -259,7 +282,7 @@ special (void)
   mpfr_set_prec (x, 33);
   mpfr_set_str_binary (x, "0.111011011011110001100111111001000e-10");
   mpfr_set_prec (z, 157);
-  inexact = mpfr_sqrt (z, x, GMP_RNDN);
+  inexact = test_sqrt (z, x, GMP_RNDN);
   mpfr_set_prec (x, 157);
   mpfr_set_str_binary (x, "0.11110110101100101111001011100011100011100001101010111011010000100111011000111110100001001011110011111100101110010110010110011001011011010110010000011001101E-5");
   if (mpfr_cmp (x, z))
@@ -281,13 +304,13 @@ special (void)
       mpfr_set_ui (x, 1, GMP_RNDN);
       mpfr_nextabove (x);
       /* 1.0 < x <= 1.5 thus 1 < sqrt(x) <= 1.23 */
-      inexact = mpfr_sqrt (z, x, GMP_RNDN);
+      inexact = test_sqrt (z, x, GMP_RNDN);
       MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0);
-      inexact = mpfr_sqrt (z, x, GMP_RNDZ);
+      inexact = test_sqrt (z, x, GMP_RNDZ);
       MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0);
-      inexact = mpfr_sqrt (z, x, GMP_RNDU);
+      inexact = test_sqrt (z, x, GMP_RNDU);
       MPFR_ASSERTN(inexact > 0 && mpfr_cmp_ui_2exp (z, 3, -1) == 0);
-      inexact = mpfr_sqrt (z, x, GMP_RNDD);
+      inexact = test_sqrt (z, x, GMP_RNDD);
       MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0);
     }
 
@@ -305,19 +328,19 @@ special (void)
       mpfr_add_ui (x, x, 1, GMP_RNDN);
       /* now x = 2^(mp_bits_per_limb - 1) + 1 (mp_bits_per_limb bits) */
       MPFR_ASSERTN(mpfr_mul (x, x, x, GMP_RNDN) == 0); /* exact */
-      inexact = mpfr_sqrt (z, x, GMP_RNDN);
+      inexact = test_sqrt (z, x, GMP_RNDN);
       /* even rule: z should be 2^(mp_bits_per_limb - 1) */
       MPFR_ASSERTN(inexact < 0 &&
                    mpfr_cmp_ui_2exp (z, 1, mp_bits_per_limb - 1) == 0);
       mpfr_nextbelow (x);
       /* now x is just below [2^(mp_bits_per_limb - 1) + 1]^2 */
-      inexact = mpfr_sqrt (z, x, GMP_RNDN);
+      inexact = test_sqrt (z, x, GMP_RNDN);
       MPFR_ASSERTN(inexact < 0 &&
                    mpfr_cmp_ui_2exp (z, 1, mp_bits_per_limb - 1) == 0);
       mpfr_nextabove (x);
       mpfr_nextabove (x);
       /* now x is just above [2^(mp_bits_per_limb - 1) + 1]^2 */
-      inexact = mpfr_sqrt (z, x, GMP_RNDN);
+      inexact = test_sqrt (z, x, GMP_RNDN);
       if (mpfr_cmp (z, y))
         {
           printf ("Error for sqrt(x) in rounding to nearest\n");
@@ -336,7 +359,7 @@ special (void)
   mpfr_set_prec (x, 1000);
   mpfr_set_ui (x, 9, GMP_RNDN);
   mpfr_set_prec (y, 10);
-  inexact = mpfr_sqrt (y, x, GMP_RNDN);
+  inexact = test_sqrt (y, x, GMP_RNDN);
   if (mpfr_cmp_ui (y, 3) || inexact != 0)
     {
       printf ("Error in sqrt(9:1000) for prec=10\n");
@@ -344,7 +367,7 @@ special (void)
     }
   mpfr_set_prec (y, mp_bits_per_limb);
   mpfr_nextabove (x);
-  inexact = mpfr_sqrt (y, x, GMP_RNDN);
+  inexact = test_sqrt (y, x, GMP_RNDN);
   if (mpfr_cmp_ui (y, 3) || inexact >= 0)
     {
       printf ("Error in sqrt(9:1000) for prec=%u\n", mp_bits_per_limb);
@@ -355,7 +378,7 @@ special (void)
   mpfr_set_ui (y, 1, GMP_RNDN);
   mpfr_nextabove (y);
   mpfr_set (x, y, GMP_RNDN);
-  inexact = mpfr_sqrt (y, x, GMP_RNDN);
+  inexact = test_sqrt (y, x, GMP_RNDN);
   if (mpfr_cmp_ui (y, 1) || inexact >= 0)
     {
       printf ("Error in sqrt(1) for prec=%u\n", mp_bits_per_limb);
@@ -380,7 +403,7 @@ check_inexact (mp_prec_t p)
   mpfr_init2 (z, 2*p);
   mpfr_random (x);
   rnd = (mp_rnd_t) RND_RAND();
-  inexact = mpfr_sqrt (y, x, rnd);
+  inexact = test_sqrt (y, x, rnd);
   if (mpfr_mul (z, y, y, rnd)) /* exact since prec(z) = 2*prec(y) */
     {
       printf ("Error: multiplication should be exact\n");
@@ -416,32 +439,32 @@ check_nan (void)
   /* sqrt(NaN) == NaN */
   MPFR_CLEAR_FLAGS (x);
   MPFR_SET_NAN (x);
-  MPFR_ASSERTN (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_nan_p (got));
 
   /* sqrt(-1) == NaN */
   mpfr_set_si (x, -1L, GMP_RNDZ);
-  MPFR_ASSERTN (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_nan_p (got));
 
   /* sqrt(+inf) == +inf */
   MPFR_CLEAR_FLAGS (x);
   MPFR_SET_INF (x);
   MPFR_SET_POS (x);
-  MPFR_ASSERTN (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_inf_p (got));
 
   /* sqrt(-inf) == NaN */
   MPFR_CLEAR_FLAGS (x);
   MPFR_SET_INF (x);
   MPFR_SET_NEG (x);
-  MPFR_ASSERTN (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_nan_p (got));
 
   /* sqrt(-0) == 0 */
   mpfr_set_si (x, 0L, GMP_RNDZ);
   MPFR_SET_NEG (x);
-  MPFR_ASSERTN (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
+  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
   MPFR_ASSERTN (mpfr_number_p (got));
   MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0);
 
diff --git a/tests/tsub.c b/tests/tsub.c
index bdf14e00a..e0cb81fad 100644
--- a/tests/tsub.c
+++ b/tests/tsub.c
@@ -24,6 +24,31 @@ MA 02111-1307, USA. */
 
 #include "mpfr-test.h"
 
+#ifdef CHECK_EXTERNAL
+static int
+test_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
+{
+  int res;
+  int ok = rnd_mode == GMP_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
+  if (ok)
+    {
+      mpfr_print_raw (b);
+      printf (" ");
+      mpfr_print_raw (c);
+    }
+  res = mpfr_sub (a, b, c, rnd_mode);
+  if (ok)
+    {
+      printf (" ");
+      mpfr_print_raw (a);
+      printf ("\n");
+    }
+  return res;
+}
+#else
+#define test_sub mpfr_sub
+#endif
+
 static void
 check_diverse (void)
 {
@@ -40,7 +65,7 @@ check_diverse (void)
   mpfr_set_prec (z, 2);
   mpfr_setmax (y, __gmpfr_emax);
   mpfr_set_str_binary (z, "0.1E-10"); /* tiny */
-  mpfr_sub (x, y, z, GMP_RNDN); /* should round to 2^emax, i.e. overflow */
+  test_sub (x, y, z, GMP_RNDN); /* should round to 2^emax, i.e. overflow */
   if (!mpfr_inf_p (x) || mpfr_sgn (x) < 0)
     {
       printf ("Error in mpfr_sub(a,b,c,RNDN) for b=maxfloat, prec(a)<prec(b), c tiny\n");
@@ -53,7 +78,7 @@ check_diverse (void)
   mpfr_set_prec (z, 2);
   mpfr_set_ui (y, 1, GMP_RNDN);
   mpfr_set_si (z, -2, GMP_RNDN);
-  mpfr_sub (x, y, z, GMP_RNDD);
+  test_sub (x, y, z, GMP_RNDD);
   if (mpfr_cmp_ui (x, 3))
     {
       printf ("Error in mpfr_sub(1,-1,RNDD)\n");
@@ -65,7 +90,7 @@ check_diverse (void)
   mpfr_set_prec (z, 288);
   mpfr_set_str_binary (y, "0.111000110011000001000111101010111011110011101001101111111110000011100101000001001010110010101010011001010100000001110011110001010101101010001011101110100100001011110100110000101101100011010001001011011010101010000010001101001000110010010111111011110001111101001000101101001100101100101000E80");
   mpfr_set_str_binary (z, "0.100001111111101001011010001100110010100111001110000110011101001011010100001000000100111011010110110010000000000010101101011000010000110001110010100001100101011100100100001011000100011110000001010101000100011101001000010111100000111000111011001000100100011000100000010010111000000100100111E-258");
-  inexact = mpfr_sub (x, y, z, GMP_RNDN);
+  inexact = test_sub (x, y, z, GMP_RNDN);
   if (inexact <= 0)
     {
       printf ("Wrong inexact flag for prec=288\n");
@@ -77,7 +102,7 @@ check_diverse (void)
   mpfr_set_prec (z, 63);
   mpfr_set_str_binary (x, "0.101101111011011100100100100111E31");
   mpfr_set_str_binary (y, "0.111110010010100100110101101010001001100101110001000101110111111E-1");
-  mpfr_sub (z, x, y, GMP_RNDN);
+  test_sub (z, x, y, GMP_RNDN);
   mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
   if (mpfr_cmp (z, y))
     {
@@ -117,7 +142,7 @@ check_diverse (void)
   mpfr_set_prec (y, 32);
   mpfr_set_str_binary (x, "0.10110111101001110100100101111000E0");
   mpfr_set_str_binary (y, "0.10001100100101000100110111000100E0");
-  if ((inexact = mpfr_sub (x, x, y, GMP_RNDN)))
+  if ((inexact = test_sub (x, x, y, GMP_RNDN)))
     {
       printf ("Wrong inexact flag (2): got %d instead of 0\n", inexact);
       exit (1);
@@ -127,7 +152,7 @@ check_diverse (void)
   mpfr_set_prec (y, 32);
   mpfr_set_str_binary (x, "0.11111000110111011000100111011010E0");
   mpfr_set_str_binary (y, "0.10011111101111000100001000000000E-3");
-  if ((inexact = mpfr_sub (x, x, y, GMP_RNDN)))
+  if ((inexact = test_sub (x, x, y, GMP_RNDN)))
     {
       printf ("Wrong inexact flag (1): got %d instead of 0\n", inexact);
       exit (1);
@@ -164,13 +189,13 @@ check_diverse (void)
   mpfr_set_prec (z, 64);
   mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
   mpfr_set_str_binary (y, "0.1110111011110001110111011111111111101000011001011100101100101100");
-  mpfr_sub (z, x, y, GMP_RNDZ);
+  test_sub (z, x, y, GMP_RNDZ);
   if (mpfr_cmp_ui (z, 1))
     {
       printf ("Error in mpfr_sub (1)\n");
       exit (1);
     }
-  mpfr_sub (z, x, y, GMP_RNDU);
+  test_sub (z, x, y, GMP_RNDU);
   mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001");
   if (mpfr_cmp (z, x))
     {
@@ -178,7 +203,7 @@ check_diverse (void)
       exit (1);
     }
   mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
-  mpfr_sub (z, x, y, GMP_RNDN);
+  test_sub (z, x, y, GMP_RNDN);
   if (mpfr_cmp_ui (z, 1))
     {
       printf ("Error in mpfr_sub (3)\n");
@@ -186,7 +211,7 @@ check_diverse (void)
     }
   mpfr_set_prec (x, 66);
   mpfr_set_str_binary (x, "1.11101110111100011101110111111111111010000110010111001011001010111");
-  mpfr_sub (z, x, y, GMP_RNDN);
+  test_sub (z, x, y, GMP_RNDN);
   if (mpfr_cmp_ui (z, 1))
     {
       printf ("Error in mpfr_sub (4)\n");
@@ -195,7 +220,7 @@ check_diverse (void)
 
   /* check in-place operations */
   mpfr_set_ui (x, 1, GMP_RNDN);
-  mpfr_sub (x, x, x, GMP_RNDN);
+  test_sub (x, x, x, GMP_RNDN);
   if (mpfr_cmp_ui(x, 0))
     {
       printf ("Error for mpfr_sub (x, x, x, GMP_RNDN) with x=1.0\n");
@@ -207,7 +232,7 @@ check_diverse (void)
   mpfr_set_prec (z, 53);
   mpfr_set_str1 (x, "1.229318102e+09");
   mpfr_set_str1 (y, "2.32221184180698677665e+05");
-  mpfr_sub (z, x, y, GMP_RNDN);
+  test_sub (z, x, y, GMP_RNDN);
   if (mpfr_cmp_str1 (z, "1229085880.815819263458251953125"))
     {
       printf ("Error in mpfr_sub (1.22e9 - 2.32e5)\n");
@@ -222,7 +247,7 @@ check_diverse (void)
   mpfr_set_prec (z, 54);
   mpfr_set_str_binary (x, "0.11111100100000000011000011100000101101010001000111E-401");
   mpfr_set_str_binary (y, "0.10110000100100000101101100011111111011101000111000101E-464");
-  mpfr_sub (z, x, y, GMP_RNDN);
+  test_sub (z, x, y, GMP_RNDN);
   if (mpfr_cmp (z, x)) {
     printf ("mpfr_sub(z, x, y) failed for prec(x)=112, prec(y)=98\n");
     printf ("expected "); mpfr_print_binary (x); puts ("");
@@ -239,7 +264,7 @@ check_diverse (void)
   mpfr_set_prec (y, 5);
   mpfr_set_str_binary (x, "1e-12");
   mpfr_set_ui (y, 1, GMP_RNDN);
-  mpfr_sub (x, y, x, GMP_RNDD);
+  test_sub (x, y, x, GMP_RNDD);
   mpfr_set_str_binary (y, "0.11111");
   if (mpfr_cmp (x, y))
     {
@@ -264,12 +289,12 @@ check_diverse (void)
   mpfr_set_prec (z, 10);
   mpfr_set_ui (y, 0, GMP_RNDN);
   mpfr_set_str_binary (z, "0.100000000000000000000100E15");
-  if (mpfr_sub (x, y, z, GMP_RNDN) <= 0)
+  if (test_sub (x, y, z, GMP_RNDN) <= 0)
     {
       printf ("Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n");
       exit (1);
     }
-  if (mpfr_sub (x, z, y, GMP_RNDN) >= 0)
+  if (test_sub (x, z, y, GMP_RNDN) >= 0)
     {
       printf ("Wrong inexact flag in x=mpfr_sub(z,0) for prec(z)>prec(x)\n");
       exit (1);
@@ -296,8 +321,8 @@ bug_ddefour(void)
     mpfr_mul_2exp( ex, ex, 906, GMP_RNDN);
     mpfr_log( tot, ex, GMP_RNDN);
     mpfr_set( ex1, tot, GMP_RNDN); /* ex1 = high(tot) */
-    mpfr_sub( ex2, tot, ex1, GMP_RNDN); /* ex2 = high(tot - ex1) */
-    mpfr_sub( tot1, tot, ex1, GMP_RNDN); /* tot1 = tot - ex1 */
+    test_sub( ex2, tot, ex1, GMP_RNDN); /* ex2 = high(tot - ex1) */
+    test_sub( tot1, tot, ex1, GMP_RNDN); /* tot1 = tot - ex1 */
     mpfr_set( ex3, tot1, GMP_RNDN); /* ex3 = high(tot - ex1) */
 
     if (mpfr_cmp(ex2, ex3))
@@ -334,8 +359,8 @@ check_two_sum (mp_prec_t p)
   if (mpfr_cmpabs (x, y) < 0)
     mpfr_swap (x, y);
   rnd = GMP_RNDN;
-  inexact = mpfr_sub (u, x, y, rnd);
-  mpfr_sub (v, u, x, rnd);
+  inexact = test_sub (u, x, y, rnd);
+  test_sub (v, u, x, rnd);
   mpfr_add (w, v, y, rnd);
   /* as u = (x-y) - w, we should have inexact and w of opposite signs */
   if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
@@ -380,7 +405,7 @@ check_inexact (void)
   mpfr_set_prec (y, 2);
   mpfr_set_si (y, -1, GMP_RNDN);
   mpfr_div_2exp (y, y, 4, GMP_RNDN); /* y = -1/16 */
-  inexact = mpfr_sub (y, y, x, GMP_RNDN); /* y = round(-7/16) = -1/2 */
+  inexact = test_sub (y, y, x, GMP_RNDN); /* y = round(-7/16) = -1/2 */
   if (inexact >= 0)
     {
       printf ("Error: wrong inexact flag for -1/16 - (6/16)\n");
@@ -412,14 +437,14 @@ check_inexact (void)
               pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u));
               mpfr_set_prec (z, pz);
               rnd = (mp_rnd_t) RND_RAND();
-              if (mpfr_sub (z, x, u, rnd))
+              if (test_sub (z, x, u, rnd))
                 {
                   printf ("z <- x - u should be exact\n");
                   exit (1);
                 }
                 {
                   rnd = (mp_rnd_t) RND_RAND ();
-                  inexact = mpfr_sub (y, x, u, rnd);
+                  inexact = test_sub (y, x, u, rnd);
                   cmp = mpfr_cmp (y, z);
                   if (((inexact == 0) && (cmp != 0)) ||
                       ((inexact > 0) && (cmp <= 0)) ||