Fixed bug in cache.c with directed rounding, and added corresponding

testcase in tests/tconst_pi.c [Ported changeset r6527 from the trunk]

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/2.4@6541 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 77 insertions, 6 deletions

diff --git a/cache.c b/cache.c
index 400bdc656..5fb247a15 100644
--- a/cache.c
+++ b/cache.c
@@ -62,21 +62,28 @@ mpfr_cache (mpfr_ptr dest, mpfr_cache_t cache, mp_rnd_t rnd)
 
       /* Update the cache. */
       pold = prec;
-      mpfr_prec_round (cache->x, pold, GMP_RNDN);
+      /* no need to keep the previous value */
+      mpfr_set_prec (cache->x, pold);
       cache->inexact = (*cache->func) (cache->x, GMP_RNDN);
     }
 
+  /* now pold >= prec is the precision of cache->x */
+
   /* First, check if the cache has the exact value (unlikely).
      Else the exact value is between (assuming x=cache->x > 0):
        x and x+ulp(x) if cache->inexact < 0,
        x-ulp(x) and x if cache->inexact > 0,
      and abs(x-exact) <= ulp(x)/2. */
+
+  /* we assume all cached constants are positive */
   MPFR_ASSERTN (MPFR_IS_POS (cache->x)); /* TODO... */
   sign = MPFR_SIGN (cache->x);
   MPFR_SET_EXP (dest, MPFR_GET_EXP (cache->x));
   MPFR_SET_SIGN (dest, sign);
+
+  /* round cache->x from precision pold down to precision prec */
   MPFR_RNDRAW_GEN (inexact, dest,
-                   MPFR_MANT (cache->x), MPFR_PREC (cache->x), rnd, sign,
+                   MPFR_MANT (cache->x), pold, rnd, sign,
                    if (MPFR_UNLIKELY (cache->inexact == 0))
                      {
                        if ((_sp[0] & _ulp) == 0)
@@ -89,7 +96,7 @@ mpfr_cache (mpfr_ptr dest, mpfr_cache_t cache, mp_rnd_t rnd)
                      }
                    else if (cache->inexact < 0)
                      goto addoneulp;
-                   else
+                   else /* cache->inexact > 0 */
                      {
                        inexact = -sign;
                        goto trunc_doit;
@@ -97,6 +104,7 @@ mpfr_cache (mpfr_ptr dest, mpfr_cache_t cache, mp_rnd_t rnd)
                    if (MPFR_UNLIKELY (++MPFR_EXP (dest) > __gmpfr_emax))
                      mpfr_overflow (dest, rnd, sign);
                   );
+
   if (MPFR_LIKELY (cache->inexact != 0))
     {
       switch (rnd)
@@ -107,7 +115,10 @@ mpfr_cache (mpfr_ptr dest, mpfr_cache_t cache, mp_rnd_t rnd)
             {
               inexact = cache->inexact;
               if (inexact > 0)
-                mpfr_nextbelow (dest);
+                {
+                  mpfr_nextbelow (dest);
+                  inexact = -inexact;
+                }
             }
           break;
         case GMP_RNDU:
@@ -115,7 +126,10 @@ mpfr_cache (mpfr_ptr dest, mpfr_cache_t cache, mp_rnd_t rnd)
             {
               inexact = cache->inexact;
               if (inexact < 0)
-                mpfr_nextabove (dest);
+                {
+                  mpfr_nextabove (dest);
+                  inexact = -inexact;
+                }
             }
           break;
         default: /* GMP_RNDN */
diff --git a/tests/tconst_pi.c b/tests/tconst_pi.c
index 1e18db1c8..bdad7614d 100644
--- a/tests/tconst_pi.c
+++ b/tests/tconst_pi.c
@@ -71,6 +71,60 @@ my_const_pi (mpfr_ptr x, mpfr_srcptr y, mp_rnd_t r)
 #define TEST_FUNCTION my_const_pi
 #include "tgeneric.c"
 
+static void
+bug20091030 (void)
+{
+  mpfr_t x, x_ref;
+  int inex, inex_ref;
+  mp_prec_t p;
+  mpfr_rnd_t r;
+
+  mpfr_free_cache ();
+  mpfr_init2 (x, MPFR_PREC_MIN);
+  for (p = MPFR_PREC_MIN; p <= 100; p++)
+    {
+      mpfr_set_prec (x, p);
+      inex = mpfr_const_pi (x, GMP_RNDU);
+      if (inex < 0)
+        {
+          printf ("Error, inex < 0 for RNDU (prec=%lu)\n", p);
+          exit (1);
+        }
+      inex = mpfr_const_pi (x, GMP_RNDD);
+      if (inex > 0)
+        {
+          printf ("Error, inex > 0 for RNDD (prec=%lu)\n", p);
+          exit (1);
+        }
+    }
+  mpfr_free_cache ();
+  mpfr_init2 (x_ref, MPFR_PREC_MIN);
+  for (p = MPFR_PREC_MIN; p <= 100; p++)
+    {
+      mpfr_set_prec (x, p + 10);
+      mpfr_const_pi (x, GMP_RNDN);
+      mpfr_set_prec (x, p);
+      mpfr_set_prec (x_ref, p);
+      for (r = 0; r < GMP_RND_MAX; r++)
+        {
+          inex = mpfr_const_pi (x, r);
+          inex_ref = mpfr_const_pi_internal (x_ref, r);
+          if (inex != inex_ref || mpfr_cmp (x, x_ref) != 0)
+            {
+              printf ("mpfr_const_pi and mpfr_const_pi_internal disagree\n");
+              printf ("mpfr_const_pi gives ");
+              mpfr_dump (x);
+              printf ("mpfr_const_pi_internal gives ");
+              mpfr_dump (x_ref);
+              printf ("inex=%d inex_ref=%d\n", inex, inex_ref);
+              exit (1);
+            }
+        }
+    }
+  mpfr_clear (x);
+  mpfr_clear (x_ref);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -118,10 +172,13 @@ main (int argc, char *argv[])
 
   mpfr_clear (x);
 
-  check_large();
+  bug20091030 ();
+
+  check_large ();
 
   test_generic (2, 200, 1);
 
   tests_end_mpfr ();
+
   return 0;
 }