Fix bug of add_q and sub_q with special values (NAN, INF and ZERO).

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

1 files changed, 53 insertions, 4 deletions

diff --git a/gmp_op.c b/gmp_op.c
index 00534badc..0b6fc79bc 100644
--- a/gmp_op.c
+++ b/gmp_op.c
@@ -108,10 +108,32 @@ mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
   mpfr_t t,q;
   mp_prec_t p = MPFR_PREC(y)+10;
   int res;
+
+  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
+    {
+      if (MPFR_IS_NAN(x))
+	{
+	  MPFR_SET_NAN(y);
+	  MPFR_RET_NAN;
+	}
+      else if (MPFR_IS_INF(x))
+	{
+	  MPFR_SET_INF (y);
+	  MPFR_SET_SAME_SIGN (y, x);
+	  MPFR_RET (0);
+	}
+      else
+	{
+	  MPFR_ASSERTD( MPFR_IS_ZERO(x));
+	  return mpfr_set_q (y, z, rnd_mode);
+	}
+    }
+
   mpfr_inits2(p, t, q, NULL);
   do {
-    mpfr_set_q(q, z, GMP_RNDN);  /* Error <= 1/2ulp */
-    mpfr_add(t, x, q, GMP_RNDN); /* Error <= 1 ulp  */
+    mpfr_set_q(q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
+    mpfr_add(t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */
+    /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) */
     res = mpfr_can_round(t, p-1, GMP_RNDN, GMP_RNDZ, 
 			 MPFR_PREC(y) + (rnd_mode == GMP_RNDN) );
     if (!res)
@@ -132,10 +154,37 @@ mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode)
   mpfr_t t,q;
   mp_prec_t p = MPFR_PREC(y)+10;
   int res;
+
+  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
+    {
+      if (MPFR_IS_NAN(x))
+        {
+          MPFR_SET_NAN (y);
+          MPFR_RET_NAN;
+        }
+      else if (MPFR_IS_INF(x))
+        {
+          MPFR_SET_INF (y);
+          MPFR_SET_SAME_SIGN (y, x);
+          MPFR_RET (0);
+        }
+      else
+        {
+          MPFR_ASSERTD( MPFR_IS_ZERO(x));
+	  if (rnd_mode == GMP_RNDU) 
+	    rnd_mode = GMP_RNDD;
+	  else if (rnd_mode == GMP_RNDD)
+	    rnd_mode = GMP_RNDU;
+          res =  mpfr_set_q (y, z, rnd_mode);
+	  MPFR_CHANGE_SIGN (y);
+	  return -res;
+        }
+    }
+
   mpfr_inits2(p, t, q, NULL);
   do {
-    mpfr_set_q(q, z, GMP_RNDN);  /* Error <= 1/2ulp */
-    mpfr_sub(t, x, q, GMP_RNDN); /* Error <= 1 ulp  */
+    mpfr_set_q(q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
+    mpfr_sub(t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */
     res = mpfr_can_round(t, p-1, GMP_RNDN, GMP_RNDZ,
                          MPFR_PREC(y) + (rnd_mode == GMP_RNDN) );
     if (!res)