[doc/sum.txt] Added a new table concerning the correction term.

[src/sum.c] Updated rounding.

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

2 files changed, 116 insertions, 45 deletions

diff --git a/doc/sum.txt b/doc/sum.txt
index baca8d577..b51b95456 100644
--- a/doc/sum.txt
+++ b/doc/sum.txt
@@ -39,8 +39,39 @@ sum_raw (obtained from code refactoring for steps 3 to 6, and 8) should
 be introduced and described. Pre-rounding in (7) and rounding in (8)
 have been merged, and this [previously pre-rounding in (7)] is still
 done at the same time as sign handling with a single mpn operation.
-The table concerning the correction term is also incorrect (the initial
-truncation also has an influence, but it isn't taken into account yet).
+The table concerning the correction term and the ternary value in case
+of TMD should be changed to:
+
+     rnd   tmd    R    sst      correction   ternary
+     ───────────────────────────────────────────────
+      N     1     0     -            0          +
+      N     1     0     0            0          0
+      N     1     0     +            0          -
+      N     1     1     -           +1          +
+      N     1     1     0           +1          0
+      N     1     1     +           +1          -
+     ───────────────────────────────────────────────
+      N     2     0     -            0          -
+      N     2     0     0            ?          ?     (halfway case)
+      N     2     0     +           +1          +
+      N     2     1     -            0          -
+      N     2     1     0            ?          ?     (halfway case)
+      N     2     1     +           +1          +
+     ───────────────────────────────────────────────
+      D     1     0     -           -1          -
+      D     1     0     0            0          0
+      D     1     0     +            0          -
+      D     1     1     -            0          -
+      D     1     1     0           +1          0
+      D     1     1     +           +1          -
+     ───────────────────────────────────────────────
+      U     1     0     -            0          +
+      U     1     0     0            0          0
+      U     1     0     +           +1          +
+      U     1     1     -           +1          +
+      U     1     1     0           +1          0
+      U     1     1     +           +2          +
+     ───────────────────────────────────────────────
 
 The general ideas of the algorithm:
 
diff --git a/src/sum.c b/src/sum.c
index 22cdb01ce..6edc6eca4 100644
--- a/src/sum.c
+++ b/src/sum.c
@@ -541,6 +541,8 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
     sd = (mpfr_prec_t) sn * GMP_NUMB_BITS - sq;
     sh = cancel % GMP_NUMB_BITS;
 
+    MPFR_ASSERTD (sd >= 0 && sd < GMP_NUMB_BITS);
+
     if (MPFR_LIKELY (u > minexp))
       {
         mpfr_prec_t tq;
@@ -779,6 +781,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
         int sst;  /* sign of the secondary term */
 
         MPFR_ASSERTD (maxexp > MPFR_EXP_MIN);
+        MPFR_ASSERTD (tmd == 1 || tmd == 2);
 
         /* New accumulator size */
         ws = MPFR_PREC2LIMBS (wq - sq);
@@ -855,10 +858,9 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
             /* For halfway cases, let's virtually eliminate them
                by setting a sst equivalent to a non-halfway case,
                which depends on the last bit of the pre-rounded
-               result and the sign. */
+               result. */
             MPFR_ASSERTD (rnd == MPFR_RNDN && tmd == 2);
-            sst = (sump[0] & (MPFR_LIMB_ONE << sd)) ?
-              (pos ? 1 : -1) : (pos ? -1 : 1);
+            sst = (sump[0] & (MPFR_LIMB_ONE << sd)) ? 1 : -1;
           }
 
         MPFR_LOG_MSG (("[TMD] tmd=%d rbit=%d sst=%d\n",
@@ -866,25 +868,20 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
 
         MPFR_COV_SET (sum_tmd[(int) rnd][tmd-1][rbit][sst+1][pos]);
 
-        /* FIXME: Update the following. */
-
         inex =
           MPFR_IS_LIKE_RNDD (rnd, pos ? 1 : -1) ? (sst ? -1 : 0) :
           MPFR_IS_LIKE_RNDU (rnd, pos ? 1 : -1) ? (sst ?  1 : 0) :
           (MPFR_ASSERTD (rnd == MPFR_RNDN),
            tmd == 1 ? - sst : sst);
 
-        corr =
-          (tmd == 2 && rbit == 0 && sst > 0) ||
-          (rbit == 1 &&
-           ((MPFR_IS_LIKE_RNDD (rnd, pos ? 1 : -1) && sst >= 0) ||
-            (MPFR_IS_LIKE_RNDU (rnd, pos ? 1 : -1) && sst > 0))) ?
-          (pos ? 1 : -1) :
-          (tmd == 2 && rbit != 0 && sst < 0) ||
-          (rbit == 0 &&
-           ((MPFR_IS_LIKE_RNDU (rnd, pos ? 1 : -1) && sst <= 0) ||
-            (MPFR_IS_LIKE_RNDD (rnd, pos ? 1 : -1) && sst < 0))) ?
-          (pos ? -1 : 1) : 0;
+        if (tmd == 2 && sst == (rbit ? -1 : 1))
+          corr = 1 - (int) rbit;
+        else if (MPFR_IS_LIKE_RNDD (rnd, pos ? 1 : -1) && sst == -1)
+          corr = (int) rbit - 1;
+        else if (MPFR_IS_LIKE_RNDU (rnd, pos ? 1 : -1) && sst == +1)
+          corr = (int) rbit + 1;
+        else
+          corr = (int) rbit;
       }
 
     /* Sign handling (-> absolute value and sign), together with
@@ -892,9 +889,6 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
        as this is necessarily the case when the TMD did not occur. */
     MPFR_ASSERTD (corr >= -1 && corr <= 2);
 
-    /* FIXME: Support the case corr == 2. Fix the case corr == -1 on a
-       negative number with a change of the binade. */
-
     if (pos)
       {
         MPFR_ASSERTD (MPFR_LIMB_MSB (sump[sn-1]) != 0);
@@ -903,13 +897,24 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
 
         if (corr > 0)
           {
-            mp_limb_t carry_out;
+            mp_limb_t corr2, carry_out;
+
+            corr2 = (mp_limb_t) corr << sd;
+            /* If corr == 2 && sd == GMP_NUMB_BITS - 1, this overflows
+               to corr2 = 0. This case is taken into account below. */
+
+            carry_out = corr2 != 0 ? mpn_add_1 (sump, sump, sn, corr2) :
+              (MPFR_ASSERTD (sn > 1),
+               mpn_add_1 (sump + 1, sump + 1, sn - 1, MPFR_LIMB_ONE));
 
-            carry_out = mpn_add_1 (sump, sump, sn, MPFR_LIMB_ONE << sd);
             MPFR_ASSERTD (sump[sn-1] >> (GMP_NUMB_BITS - 1) == !carry_out);
+
             if (MPFR_UNLIKELY (carry_out))
               {
-                sump[sn-1] = MPFR_LIMB_HIGHBIT;
+                /* Note: The | is important when sump[sn-1] is not 0
+                   (this can occur with sn = 1 and corr = 2). TODO:
+                   Add something to make sure that this is tested. */
+                sump[sn-1] |= MPFR_LIMB_HIGHBIT;
                 e++;
               }
           }
@@ -917,6 +922,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
         if (corr < 0)
           {
             mpn_sub_1 (sump, sump, sn, MPFR_LIMB_ONE << sd);
+
             if (MPFR_UNLIKELY (MPFR_LIMB_MSB (sump[sn-1]) == 0))
               {
                 sump[sn-1] |= MPFR_LIMB_HIGHBIT;
@@ -926,40 +932,74 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd,
       }
     else
       {
-        mp_limb_t corr2;
+        MPFR_ASSERTD (MPFR_LIMB_MSB (sump[sn-1]) == 0);
+        MPFR_SET_NEG (sum);
 
         /* abs(x + corr) = - (x + corr) = com(x) + (1 - corr) */
 
         /* We want to avoid separate mpn_com (or mpn_neg) and mpn_add_1
            (or mpn_sub_1) operations, as they could yield two loops in
-           some particular cases involving a long sequence of 0's in the
-           low significant bits (except the least significant bit, which
-           doesn't matter). But we can just do the correction operation
-           on the least significant limb, then do either a mpn_com or a
-           mpn_neg on the remaining limbs, depending on the carry (BTW,
-           mpn_neg is just a mpn_com with an initial carry propagation:
-           after some point, mpn_neg does a complement). */
+           some particular cases involving a long sequence of 0's in
+           the low significant bits (except the least significant bit,
+           which doesn't matter). */
 
-        MPFR_ASSERTD (MPFR_LIMB_MSB (sump[sn-1]) == 0);
-        MPFR_SET_NEG (sum);
+        if (corr <= 1)
+          {
+            mp_limb_t corr2;
+
+            /* Here we can just do the correction operation on the
+               least significant limb, then do either a mpn_com or
+               a mpn_neg on the remaining limbs, depending on the
+               carry (BTW, mpn_neg is just a mpn_com with an initial
+               carry propagation: after some point, mpn_neg does a
+               complement). */
 
-        MPFR_ASSERTD (corr >= -1 && corr <= 1);
-        corr2 = (mp_limb_t) (1 - corr) << sd;
-        /* Note: If corr = -1, this can overflow to corr2 = 0.
-           This case is taken into account below. */
+            corr2 = (mp_limb_t) (1 - corr) << sd;
+            /* Note: If corr = -1, this can overflow to corr2 = 0.
+               This case is taken into account below. */
 
-        sump[0] = (~ (sump[0] | MPFR_LIMB_MASK (sd))) + corr2;
+            sump[0] = (~ (sump[0] | MPFR_LIMB_MASK (sd))) + corr2;
 
-        if (sump[0] < corr2 || (corr2 == 0 && corr < 0))
+            if (sump[0] < corr2 || (corr2 == 0 && corr < 0))
+              {
+                if (sn == 1 || ! mpn_neg (sump + 1, sump + 1, sn - 1))
+                  {
+                    /* Note: The | is important when sump[sn-1] is not 0
+                       (this can occur with sn = 1 and corr = -1). TODO:
+                       Add something to make sure that this is tested. */
+                    sump[sn-1] |= MPFR_LIMB_HIGHBIT;
+                    e++;
+                  }
+              }
+            else if (sn > 1)
+              mpn_com (sump + 1, sump + 1, sn - 1);
+          }
+        else  /* corr == 2 */
           {
-            if (sn == 1 || ! mpn_neg (sump + 1, sump + 1, sn - 1))
+            mp_limb_t corr2, c;
+            mp_size_t i = 1;
+
+            /* We want to compute com(x) - 1, but GMP doesn't have an
+               operation for that. The fact is that a sequence of low
+               significant bits 1 is invariant. Starting at the first
+               low significant bit 0, we can do the complement with
+               mpn_com. */
+
+            corr2 = MPFR_LIMB_ONE << sd;
+            c = ~ (sump[0] | MPFR_LIMB_MASK (sd));
+            sump[0] = c - corr2;
+
+            if (c == 0)
               {
-                e++;
-                sump[sn-1] = MPFR_LIMB_HIGHBIT;
+                while (MPFR_ASSERTD (i < sn), sump[i] == MPFR_LIMB_MAX)
+                  i++;
+                sump[i] = (~ sump[i]) - 1;
+                i++;
               }
+
+            if (i < sn)
+              mpn_com (sump + i, sump + i, sn - i);
           }
-        else if (sn > 1)
-          mpn_com (sump + 1, sump + 1, sn - 1);
       }
 
     MPFR_ASSERTD (MPFR_LIMB_MSB (sump[sn-1]) != 0);