path: root/patches/mpfr/3.0.0/160-mpfr_sub1.patch

diff -Naurd mpfr-3.0.0-a/PATCHES mpfr-3.0.0-b/PATCHES
--- mpfr-3.0.0-a/PATCHES	2010-10-21 20:59:32.000000000 +0000
+++ mpfr-3.0.0-b/PATCHES	2010-10-21 20:59:32.000000000 +0000
@@ -0,0 +1 @@
+mpfr_sub1
diff -Naurd mpfr-3.0.0-a/VERSION mpfr-3.0.0-b/VERSION
--- mpfr-3.0.0-a/VERSION	2010-10-21 20:28:38.000000000 +0000
+++ mpfr-3.0.0-b/VERSION	2010-10-21 20:59:32.000000000 +0000
@@ -1 +1 @@
-3.0.0-p5
+3.0.0-p6
diff -Naurd mpfr-3.0.0-a/mpfr.h mpfr-3.0.0-b/mpfr.h
--- mpfr-3.0.0-a/mpfr.h	2010-10-21 20:28:38.000000000 +0000
+++ mpfr-3.0.0-b/mpfr.h	2010-10-21 20:59:32.000000000 +0000
@@ -27,7 +27,7 @@
 #define MPFR_VERSION_MAJOR 3
 #define MPFR_VERSION_MINOR 0
 #define MPFR_VERSION_PATCHLEVEL 0
-#define MPFR_VERSION_STRING "3.0.0-p5"
+#define MPFR_VERSION_STRING "3.0.0-p6"
 
 /* Macros dealing with MPFR VERSION */
 #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c))
diff -Naurd mpfr-3.0.0-a/sub1.c mpfr-3.0.0-b/sub1.c
--- mpfr-3.0.0-a/sub1.c	2010-06-10 11:00:14.000000000 +0000
+++ mpfr-3.0.0-b/sub1.c	2010-10-21 20:59:32.000000000 +0000
@@ -37,7 +37,9 @@
   mp_size_t cancel2, an, bn, cn, cn0;
   mp_limb_t *ap, *bp, *cp;
   mp_limb_t carry, bb, cc, borrow = 0;
-  int inexact, shift_b, shift_c, is_exact = 1, down = 0, add_exp = 0;
+  int inexact, shift_b, shift_c, add_exp = 0;
+  int cmp_low = 0; /* used for rounding to nearest: 0 if low(b) = low(c),
+                      negative if low(b) < low(c), positive if low(b)>low(c) */
   int sh, k;
   MPFR_TMP_DECL(marker);
 
@@ -196,7 +198,8 @@
     }
 
 #ifdef DEBUG
-  printf ("shift_b=%d shift_c=%d diffexp=%lu\n", shift_b, shift_c,
+  printf ("rnd=%s shift_b=%d shift_c=%d diffexp=%lu\n",
+          mpfr_print_rnd_mode (rnd_mode), shift_b, shift_c,
           (unsigned long) diff_exp);
 #endif
 
@@ -307,17 +310,18 @@
     {
       if (MPFR_LIKELY(sh))
         {
-          is_exact = (carry == 0);
           /* can decide except when carry = 2^(sh-1) [middle]
              or carry = 0 [truncate, but cannot decide inexact flag] */
-          down = (carry < (MPFR_LIMB_ONE << (sh - 1)));
           if (carry > (MPFR_LIMB_ONE << (sh - 1)))
             goto add_one_ulp;
-          else if ((0 < carry) && down)
+          else if ((0 < carry) && (carry < (MPFR_LIMB_ONE << (sh - 1))))
             {
               inexact = -1; /* result if smaller than exact value */
               goto truncate;
             }
+          /* now carry = 2^(sh-1), in which case cmp_low=2,
+             or carry = 0, in which case cmp_low=0 */
+          cmp_low = (carry == 0) ? 0 : 2;
         }
     }
   else /* directed rounding: set rnd_mode to RNDZ iff toward zero */
@@ -344,12 +348,32 @@
   cn -= (long int) an + cancel2;
 
 #ifdef DEBUG
-  printf ("last %d bits from a are %lu, bn=%ld, cn=%ld\n",
+  printf ("last sh=%d bits from a are %lu, bn=%ld, cn=%ld\n",
           sh, (unsigned long) carry, (long) bn, (long) cn);
 #endif
 
+  /* for rounding to nearest, we couldn't conclude up to here in the following
+     cases:
+     1. sh = 0, then cmp_low=0: we can either truncate, subtract one ulp
+        or add one ulp: -1 ulp < low(b)-low(c) < 1 ulp
+     2. sh > 0 but the low sh bits from high(b)-high(c) equal 2^(sh-1):
+        -0.5 ulp <= -1/2^sh < low(b)-low(c)-0.5 < 1/2^sh <= 0.5 ulp
+        we can't decide the rounding, in that case cmp_low=2:
+        either we truncate and flag=-1, or we add one ulp and flag=1
+     3. the low sh>0 bits from high(b)-high(c) equal 0: we know we have to
+        truncate but we can't decide the ternary value, here cmp_low=0:
+        -0.5 ulp <= -1/2^sh < low(b)-low(c) < 1/2^sh <= 0.5 ulp
+        we always truncate and inexact can be any of -1,0,1
+  */
+
+  /* note: here cn might exceed cn0, in which case we consider a zero limb */
   for (k = 0; (bn > 0) || (cn > 0); k = 1)
     {
+      /* if cmp_low < 0, we know low(b) - low(c) < 0
+         if cmp_low > 0, we know low(b) - low(c) > 0
+            (more precisely if cmp_low = 2, low(b) - low(c) = 0.5 ulp so far)
+         if cmp_low = 0, so far low(b) - low(c) = 0 */
+
       /* get next limbs */
       bb = (bn > 0) ? bp[--bn] : 0;
       if ((cn > 0) && (cn-- <= cn0))
@@ -357,76 +381,115 @@
       else
         cc = 0;
 
-      /* down is set when low(b) < low(c) */
-      if (down == 0)
-        down = (bb < cc);
+      /* cmp_low compares low(b) and low(c) */
+      if (cmp_low == 0) /* case 1 or 3 */
+        cmp_low = (bb < cc) ? -2+k : (bb > cc) ? 1 : 0;
+
+      /* Case 1 for k=0 splits into 7 subcases:
+         1a: bb > cc + half
+         1b: bb = cc + half
+         1c: 0 < bb - cc < half
+         1d: bb = cc
+         1e: -half < bb - cc < 0
+         1f: bb - cc = -half
+         1g: bb - cc < -half
+
+         Case 2 splits into 3 subcases:
+         2a: bb > cc
+         2b: bb = cc
+         2c: bb < cc
+
+         Case 3 splits into 3 subcases:
+         3a: bb > cc
+         3b: bb = cc
+         3c: bb < cc
+      */
 
       /* the case rounding to nearest with sh=0 is special since one couldn't
          subtract above 1/2 ulp in the trailing limb of the result */
-      if ((rnd_mode == MPFR_RNDN) && sh == 0 && k == 0)
+      if (rnd_mode == MPFR_RNDN && sh == 0 && k == 0) /* case 1 for k=0 */
         {
           mp_limb_t half = MPFR_LIMB_HIGHBIT;
 
-          is_exact = (bb == cc);
-
           /* add one ulp if bb > cc + half
              truncate if cc - half < bb < cc + half
              sub one ulp if bb < cc - half
           */
 
-          if (down)
+          if (cmp_low < 0) /* bb < cc: -1 ulp < low(b) - low(c) < 0,
+                              cases 1e, 1f and 1g */
             {
               if (cc >= half)
                 cc -= half;
-              else
+              else /* since bb < cc < half, bb+half < 2*half */
                 bb += half;
+              /* now we have bb < cc + half:
+                 we have to subtract one ulp if bb < cc,
+                 and truncate if bb > cc */
             }
-          else /* bb >= cc */
+          else if (cmp_low >= 0) /* bb >= cc, cases 1a to 1d */
             {
               if (cc < half)
                 cc += half;
-              else
+              else /* since bb >= cc >= half, bb - half >= 0 */
                 bb -= half;
+              /* now we have bb > cc - half: we have to add one ulp if bb > cc,
+                 and truncate if bb < cc */
+              if (cmp_low > 0)
+                cmp_low = 2;
             }
         }
 
 #ifdef DEBUG
-      printf ("    bb=%lu cc=%lu down=%d is_exact=%d\n",
-              (unsigned long) bb, (unsigned long) cc, down, is_exact);
+      printf ("k=%u bb=%lu cc=%lu cmp_low=%d\n", k,
+              (unsigned long) bb, (unsigned long) cc, cmp_low);
 #endif
-      if (bb < cc)
+      if (cmp_low < 0) /* low(b) - low(c) < 0: either truncate or subtract
+                          one ulp */
         {
           if (rnd_mode == MPFR_RNDZ)
-            goto sub_one_ulp;
+            goto sub_one_ulp; /* set inexact=-1 */
           else if (rnd_mode != MPFR_RNDN) /* round away */
             {
               inexact = 1;
               goto truncate;
             }
-          else /* round to nearest: special case here since for sh=k=0
-                  bb = bb0 - MPFR_LIMB_HIGHBIT */
+          else /* round to nearest */
             {
-              if (is_exact && sh == 0)
-                {
-                  /* For k=0 we can't decide exactness since it may depend
-                     from low order bits.
-                     For k=1, the first low limbs matched: low(b)-low(c)<0. */
-                  if (k)
-                    {
-                      inexact = 1;
-                      goto truncate;
-                    }
-                }
-              else if (down && sh == 0)
-                goto sub_one_ulp;
-              else
-                {
-                  inexact = (is_exact) ? 1 : -1;
+              /* If cmp_low < 0 and bb > cc, then -0.5 ulp < low(b)-low(c) < 0,
+                 whatever the value of sh.
+                 If sh>0, then cmp_low < 0 implies that the initial neglected
+                 sh bits were 0 (otherwise cmp_low=2 initially), thus the
+                 weight of the new bits is less than 0.5 ulp too.
+                 If k > 0 (and sh=0) this means that either the first neglected
+                 limbs bb and cc were equal (thus cmp_low was 0 for k=0),
+                 or we had bb - cc = -0.5 ulp or 0.5 ulp.
+                 The last case is not possible here since we would have
+                 cmp_low > 0 which is sticky.
+                 In the first case (where we have cmp_low = -1), we truncate,
+                 whereas in the 2nd case we have cmp_low = -2 and we subtract
+                 one ulp.
+              */
+              if (bb > cc || sh > 0 || cmp_low == -1)
+                {  /* -0.5 ulp < low(b)-low(c) < 0,
+                      bb > cc corresponds to cases 1e and 1f1
+                      sh > 0 corresponds to cases 3c and 3b3
+                      cmp_low = -1 corresponds to case 1d3 (also 3b3) */
+                  inexact = 1;
                   goto truncate;
                 }
+              else if (bb < cc) /* here sh = 0 and low(b)-low(c) < -0.5 ulp,
+                                   this corresponds to cases 1g and 1f3 */
+                goto sub_one_ulp;
+              /* the only case where we can't conclude is sh=0 and bb=cc,
+                 i.e., we have low(b) - low(c) = -0.5 ulp (up to now), thus
+                 we don't know if we must truncate or subtract one ulp.
+                 Note: for sh=0 we can't have low(b) - low(c) = -0.5 ulp up to
+                 now, since low(b) - low(c) > 1/2^sh */
             }
         }
-      else if (bb > cc)
+      else if (cmp_low > 0) /* 0 < low(b) - low(c): either truncate or
+                               add one ulp */
         {
           if (rnd_mode == MPFR_RNDZ)
             {
@@ -437,34 +500,70 @@
             goto add_one_ulp;
           else /* round to nearest */
             {
-              if (is_exact)
+              if (bb > cc)
                 {
-                  inexact = -1;
-                  goto truncate;
+                  /* if sh=0, then bb>cc means that low(b)-low(c) > 0.5 ulp,
+                     and similarly when cmp_low=2 */
+                  if (cmp_low == 2) /* cases 1a, 1b1, 2a and 2b1 */
+                    goto add_one_ulp;
+                  /* sh > 0 and cmp_low > 0: this implies that the sh initial
+                     neglected bits were 0, and the remaining low(b)-low(c)>0,
+                     but its weight is less than 0.5 ulp */
+                  else /* 0 < low(b) - low(c) < 0.5 ulp, this corresponds to
+                          cases 3a, 1d1 and 3b1 */
+                    {
+                      inexact = -1;
+                      goto truncate;
+                    }
                 }
-              else if (down)
+              else if (bb < cc) /* 0 < low(b) - low(c) < 0.5 ulp, cases 1c,
+                                   1b3, 2b3 and 2c */
                 {
-                  inexact = 1;
+                  inexact = -1;
                   goto truncate;
                 }
-              else
-                goto add_one_ulp;
+              /* the only case where we can't conclude is bb=cc, i.e.,
+                 low(b) - low(c) = 0.5 ulp (up to now), thus we don't know
+                 if we must truncate or add one ulp. */
             }
         }
+      /* after k=0, we cannot conclude in the following cases, we split them
+         according to the values of bb and cc for k=1:
+         1b. sh=0 and cmp_low = 1 and bb-cc = half [around 0.5 ulp]
+             1b1. bb > cc: add one ulp, inex = 1
+             1b2: bb = cc: cannot conclude
+             1b3: bb < cc: truncate, inex = -1
+         1d. sh=0 and cmp_low = 0 and bb-cc = 0 [around 0]
+             1d1: bb > cc: truncate, inex = -1
+             1d2: bb = cc: cannot conclude
+             1d3: bb < cc: truncate, inex = +1
+         1f. sh=0 and cmp_low = -1 and bb-cc = -half [around -0.5 ulp]
+             1f1: bb > cc: truncate, inex = +1
+             1f2: bb = cc: cannot conclude
+             1f3: bb < cc: sub one ulp, inex = -1
+         2b. sh > 0 and cmp_low = 2 and bb=cc [around 0.5 ulp]
+             2b1. bb > cc: add one ulp, inex = 1
+             2b2: bb = cc: cannot conclude
+             2b3: bb < cc: truncate, inex = -1
+         3b. sh > 0 and cmp_low = 0 [around 0]
+             3b1. bb > cc: truncate, inex = -1
+             3b2: bb = cc: cannot conclude
+             3b3: bb < cc: truncate, inex = +1
+      */
     }
 
-  if ((rnd_mode == MPFR_RNDN) && !is_exact)
+  if ((rnd_mode == MPFR_RNDN) && cmp_low != 0)
     {
       /* even rounding rule */
       if ((ap[0] >> sh) & 1)
         {
-          if (down)
+          if (cmp_low < 0)
             goto sub_one_ulp;
           else
             goto add_one_ulp;
         }
       else
-        inexact = (down) ? 1 : -1;
+        inexact = (cmp_low > 0) ? -1 : 1;
     }
   else
     inexact = 0;
diff -Naurd mpfr-3.0.0-a/tests/tfma.c mpfr-3.0.0-b/tests/tfma.c
--- mpfr-3.0.0-a/tests/tfma.c	2010-06-10 11:00:13.000000000 +0000
+++ mpfr-3.0.0-b/tests/tfma.c	2010-10-21 20:59:32.000000000 +0000
@@ -337,6 +337,94 @@
   mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
 }
 
+static void
+bug20101018 (void)
+{
+  mpfr_t x, y, z, t, u;
+  int i;
+
+  mpfr_init2 (x, 64);
+  mpfr_init2 (y, 64);
+  mpfr_init2 (z, 64);
+  mpfr_init2 (t, 64);
+  mpfr_init2 (u, 64);
+
+  mpfr_set_str (x, "0xf.fffffffffffffffp-14766", 16, MPFR_RNDN);
+  mpfr_set_str (y, "-0xf.fffffffffffffffp+317", 16, MPFR_RNDN);
+  mpfr_set_str (z, "0x8.3ffffffffffe3ffp-14443", 16, MPFR_RNDN);
+  mpfr_set_str (t, "0x8.7ffffffffffc7ffp-14444", 16, MPFR_RNDN);
+  i = mpfr_fma (u, x, y, z, MPFR_RNDN);
+  if (mpfr_cmp (u, t) != 0)
+    {
+      printf ("Wrong result in bug20101018 (a)\n");
+      printf ("Expected ");
+      mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
+      printf ("\nGot      ");
+      mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
+      printf ("\n");
+      exit (1);
+    }
+  if (i <= 0)
+    {
+      printf ("Wrong ternary value in bug20101018 (a)\n");
+      printf ("Expected > 0\n");
+      printf ("Got      %d\n", i);
+      exit (1);
+    }
+
+  mpfr_set_str (x, "-0xf.fffffffffffffffp-11420", 16, MPFR_RNDN);
+  mpfr_set_str (y, "0xf.fffffffffffffffp+9863", 16, MPFR_RNDN);
+  mpfr_set_str (z, "0x8.fffff80ffffffffp-1551", 16, MPFR_RNDN);
+  mpfr_set_str (t, "0x9.fffff01ffffffffp-1552", 16, MPFR_RNDN);
+  i = mpfr_fma (u, x, y, z, MPFR_RNDN);
+  if (mpfr_cmp (u, t) != 0)
+    {
+      printf ("Wrong result in bug20101018 (b)\n");
+      printf ("Expected ");
+      mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
+      printf ("\nGot      ");
+      mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
+      printf ("\n");
+      exit (1);
+    }
+  if (i <= 0)
+    {
+      printf ("Wrong ternary value in bug20101018 (b)\n");
+      printf ("Expected > 0\n");
+      printf ("Got      %d\n", i);
+      exit (1);
+    }
+
+  mpfr_set_str (x, "0xf.fffffffffffffffp-2125", 16, MPFR_RNDN);
+  mpfr_set_str (y, "-0xf.fffffffffffffffp-6000", 16, MPFR_RNDN);
+  mpfr_set_str (z, "0x8p-8119", 16, MPFR_RNDN);
+  mpfr_set_str (t, "0x8.000000000000001p-8120", 16, MPFR_RNDN);
+  i = mpfr_fma (u, x, y, z, MPFR_RNDN);
+  if (mpfr_cmp (u, t) != 0)
+    {
+      printf ("Wrong result in bug20101018 (c)\n");
+      printf ("Expected ");
+      mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
+      printf ("\nGot      ");
+      mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
+      printf ("\n");
+      exit (1);
+    }
+  if (i <= 0)
+    {
+      printf ("Wrong ternary value in bug20101018 (c)\n");
+      printf ("Expected > 0\n");
+      printf ("Got      %d\n", i);
+      exit (1);
+    }
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+  mpfr_clear (t);
+  mpfr_clear (u);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -345,6 +433,8 @@
 
   tests_start_mpfr ();
 
+  bug20101018 ();
+
   mpfr_init (x);
   mpfr_init (s);
   mpfr_init (y);
diff -Naurd mpfr-3.0.0-a/tests/tsub.c mpfr-3.0.0-b/tests/tsub.c
--- mpfr-3.0.0-a/tests/tsub.c	2010-06-10 11:00:13.000000000 +0000
+++ mpfr-3.0.0-b/tests/tsub.c	2010-10-21 20:59:32.000000000 +0000
@@ -201,6 +201,8 @@
   if (mpfr_cmp (z, x))
     {
       printf ("Error in mpfr_sub (2)\n");
+      printf ("Expected "); mpfr_print_binary (x); puts ("");
+      printf ("Got      "); mpfr_print_binary (z); puts ("");
       exit (1);
     }
   mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
@@ -478,6 +480,156 @@
   mpfr_clear (u);
 }
 
+/* Bug found by Jakub Jelinek
+ * http://bugzilla.redhat.com/643657
+ * https://gforge.inria.fr/tracker/index.php?func=detail&aid=11301
+ * The consequence can be either an assertion failure (i = 2 in the
+ * testcase below, in debug mode) or an incorrectly rounded value.
+ */
+static void
+bug20101017 (void)
+{
+  mpfr_t a, b, c;
+  int inex;
+  int i;
+
+  mpfr_init2 (a, GMP_NUMB_BITS * 2);
+  mpfr_init2 (b, GMP_NUMB_BITS);
+  mpfr_init2 (c, GMP_NUMB_BITS);
+
+  /* a = 2^(2N) + k.2^(2N-1) + 2^N and b = 1
+     with N = GMP_NUMB_BITS and k = 0 or 1.
+     c = a - b should round to the same value as a. */
+
+  for (i = 2; i <= 3; i++)
+    {
+      mpfr_set_ui_2exp (a, i, GMP_NUMB_BITS - 1, MPFR_RNDN);
+      mpfr_add_ui (a, a, 1, MPFR_RNDN);
+      mpfr_mul_2ui (a, a, GMP_NUMB_BITS, MPFR_RNDN);
+      mpfr_set_ui (b, 1, MPFR_RNDN);
+      inex = mpfr_sub (c, a, b, MPFR_RNDN);
+      mpfr_set (b, a, MPFR_RNDN);
+      if (! mpfr_equal_p (c, b))
+        {
+          printf ("Error in bug20101017 for i = %d.\n", i);
+          printf ("Expected ");
+          mpfr_out_str (stdout, 16, 0, b, MPFR_RNDN);
+          putchar ('\n');
+          printf ("Got      ");
+          mpfr_out_str (stdout, 16, 0, c, MPFR_RNDN);
+          putchar ('\n');
+          exit (1);
+        }
+      if (inex >= 0)
+        {
+          printf ("Error in bug20101017 for i = %d: bad inex value.\n", i);
+          printf ("Expected negative, got %d.\n", inex);
+          exit (1);
+        }
+    }
+
+  mpfr_set_prec (a, 64);
+  mpfr_set_prec (b, 129);
+  mpfr_set_prec (c, 2);
+  mpfr_set_str_binary (b, "0.100000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000001E65");
+  mpfr_set_str_binary (c, "0.10E1");
+  inex = mpfr_sub (a, b, c, MPFR_RNDN);
+  if (mpfr_cmp_ui_2exp (a, 1, 64) != 0 || inex >= 0)
+    {
+      printf ("Error in mpfr_sub for b-c for b=2^64+1+2^(-64), c=1\n");
+      printf ("Expected result 2^64 with inex < 0\n");
+      printf ("Got "); mpfr_print_binary (a);
+      printf (" with inex=%d\n", inex);
+      exit (1);
+    }
+
+  mpfr_clears (a, b, c, (mpfr_ptr) 0);
+}
+
+/* hard test of rounding */
+static void
+check_rounding (void)
+{
+  mpfr_t a, b, c, res;
+  mpfr_prec_t p;
+  long k, l;
+  int i;
+
+#define MAXKL (2 * GMP_NUMB_BITS)
+  for (p = MPFR_PREC_MIN; p <= GMP_NUMB_BITS; p++)
+    {
+      mpfr_init2 (a, p);
+      mpfr_init2 (res, p);
+      mpfr_init2 (b, p + 1 + MAXKL);
+      mpfr_init2 (c, MPFR_PREC_MIN);
+
+      /* b = 2^p + 1 + 2^(-k), c = 2^(-l) */
+      for (k = 0; k <= MAXKL; k++)
+        for (l = 0; l <= MAXKL; l++)
+          {
+            mpfr_set_ui_2exp (b, 1, p, MPFR_RNDN);
+            mpfr_add_ui (b, b, 1, MPFR_RNDN);
+            mpfr_mul_2ui (b, b, k, MPFR_RNDN);
+            mpfr_add_ui (b, b, 1, MPFR_RNDN);
+            mpfr_div_2ui (b, b, k, MPFR_RNDN);
+            mpfr_set_ui_2exp (c, 1, -l, MPFR_RNDN);
+            i = mpfr_sub (a, b, c, MPFR_RNDN);
+            /* b - c = 2^p + 1 + 2^(-k) - 2^(-l), should be rounded to
+               2^p for l <= k, and 2^p+2 for l < k */
+            if (l <= k)
+              {
+                if (mpfr_cmp_ui_2exp (a, 1, p) != 0)
+                  {
+                    printf ("Wrong result in check_rounding\n");
+                    printf ("p=%lu k=%ld l=%ld\n", p, k, l);
+                    printf ("b="); mpfr_print_binary (b); puts ("");
+                    printf ("c="); mpfr_print_binary (c); puts ("");
+                    printf ("Expected 2^%lu\n", p);
+                    printf ("Got      "); mpfr_print_binary (a); puts ("");
+                    exit (1);
+                  }
+                if (i >= 0)
+                  {
+                    printf ("Wrong ternary value in check_rounding\n");
+                    printf ("p=%lu k=%ld l=%ld\n", p, k, l);
+                    printf ("b="); mpfr_print_binary (b); puts ("");
+                    printf ("c="); mpfr_print_binary (c); puts ("");
+                    printf ("a="); mpfr_print_binary (a); puts ("");
+                    printf ("Expected < 0, got %d\n", i);
+                    exit (1);
+                  }
+              }
+            else /* l < k */
+              {
+                mpfr_set_ui_2exp (res, 1, p, MPFR_RNDN);
+                mpfr_add_ui (res, res, 2, MPFR_RNDN);
+                if (mpfr_cmp (a, res) != 0)
+                  {
+                    printf ("Wrong result in check_rounding\n");
+                    printf ("b="); mpfr_print_binary (b); puts ("");
+                    printf ("c="); mpfr_print_binary (c); puts ("");
+                    printf ("Expected "); mpfr_print_binary (res); puts ("");
+                    printf ("Got      "); mpfr_print_binary (a); puts ("");
+                    exit (1);
+                  }
+                if (i <= 0)
+                  {
+                    printf ("Wrong ternary value in check_rounding\n");
+                    printf ("b="); mpfr_print_binary (b); puts ("");
+                    printf ("c="); mpfr_print_binary (c); puts ("");
+                    printf ("Expected > 0, got %d\n", i);
+                    exit (1);
+                  }
+              }
+          }
+
+      mpfr_clear (a);
+      mpfr_clear (res);
+      mpfr_clear (b);
+      mpfr_clear (c);
+    }
+}
+
 #define TEST_FUNCTION test_sub
 #define TWO_ARGS
 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
@@ -491,6 +643,8 @@
 
   tests_start_mpfr ();
 
+  bug20101017 ();
+  check_rounding ();
   check_diverse ();
   check_inexact ();
   bug_ddefour ();
diff -Naurd mpfr-3.0.0-a/version.c mpfr-3.0.0-b/version.c
--- mpfr-3.0.0-a/version.c	2010-10-21 20:28:38.000000000 +0000
+++ mpfr-3.0.0-b/version.c	2010-10-21 20:59:32.000000000 +0000
@@ -25,5 +25,5 @@
 const char *
 mpfr_get_version (void)
 {
-  return "3.0.0-p5";
+  return "3.0.0-p6";
 }