sub1sp.c

   1 /* mpfr_sub1sp -- internal function to perform a "real" substraction
   2    All the op must have the same precision
   3
   4 Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao projects, INRIA.
   6
   7 This file is part of the GNU MPFR Library.
   8
   9 The GNU MPFR Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 3 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  23
  24 #define MPFR_NEED_LONGLONG_H
  25 #include "mpfr-impl.h"
  26
  27 /* Check if we have to check the result of mpfr_sub1sp with mpfr_sub1 */
  28 #ifdef WANT_ASSERT
  29 # if WANT_ASSERT >= 2
  30
  31 int mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode);
  32 int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
  33 {
  34   mpfr_t tmpa, tmpb, tmpc;
  35   int inexb, inexc, inexact, inexact2;
  36
  37   mpfr_init2 (tmpa, MPFR_PREC (a));
  38   mpfr_init2 (tmpb, MPFR_PREC (b));
  39   mpfr_init2 (tmpc, MPFR_PREC (c));
  40
  41   inexb = mpfr_set (tmpb, b, MPFR_RNDN);
  42   MPFR_ASSERTN (inexb == 0);
  43
  44   inexc = mpfr_set (tmpc, c, MPFR_RNDN);
  45   MPFR_ASSERTN (inexc == 0);
  46
  47   inexact2 = mpfr_sub1 (tmpa, tmpb, tmpc, rnd_mode);
  48   inexact  = mpfr_sub1sp2(a, b, c, rnd_mode);
  49
  50   if (mpfr_cmp (tmpa, a) || inexact != inexact2)
  51     {
  52       fprintf (stderr, "sub1 & sub1sp return different values for %s\n"
  53                "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
  54                mpfr_print_rnd_mode (rnd_mode), (unsigned long) MPFR_PREC (a),
  55                (unsigned long) MPFR_PREC (b), (unsigned long) MPFR_PREC (c));
  56       mpfr_fprint_binary (stderr, tmpb);
  57       fprintf (stderr, "\nC = ");
  58       mpfr_fprint_binary (stderr, tmpc);
  59       fprintf (stderr, "\nSub1  : ");
  60       mpfr_fprint_binary (stderr, tmpa);
  61       fprintf (stderr, "\nSub1sp: ");
  62       mpfr_fprint_binary (stderr, a);
  63       fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
  64                inexact, inexact2);
  65       MPFR_ASSERTN (0);
  66     }
  67   mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
  68   return inexact;
  69 }
  70 #  define mpfr_sub1sp mpfr_sub1sp2
  71 # endif
  72 #endif
  73
  74 /* Debugging support */
  75 #ifdef DEBUG
  76 # undef DEBUG
  77 # define DEBUG(x) (x)
  78 #else
  79 # define DEBUG(x) /**/
  80 #endif
  81
  82 /* Rounding Sub */
  83
  84 /*
  85    compute sgn(b)*(|b| - |c|) if |b|>|c| else -sgn(b)*(|c| -|b|)
  86    Returns 0 iff result is exact,
  87    a negative value when the result is less than the exact value,
  88    a positive value otherwise.
  89 */
  90
  91 /* A0...Ap-1
  92  *          Cp Cp+1 ....
  93  *             <- C'p+1 ->
  94  * Cp = -1 if calculated from c mantissa
  95  * Cp = 0  if 0 from a or c
  96  * Cp = 1  if calculated from a.
  97  * C'p+1 = First bit not null or 0 if there isn't one
  98  *
  99  * Can't have Cp=-1 and C'p+1=1*/
 100
 101 /* RND = MPFR_RNDZ:
 102  *  + if Cp=0 and C'p+1=0,1,  Truncate.
 103  *  + if Cp=0 and C'p+1=-1,   SubOneUlp
 104  *  + if Cp=-1,               SubOneUlp
 105  *  + if Cp=1,                AddOneUlp
 106  * RND = MPFR_RNDA (Away)
 107  *  + if Cp=0 and C'p+1=0,-1, Truncate
 108  *  + if Cp=0 and C'p+1=1,    AddOneUlp
 109  *  + if Cp=1,                AddOneUlp
 110  *  + if Cp=-1,               Truncate
 111  * RND = MPFR_RNDN
 112  *  + if Cp=0,                Truncate
 113  *  + if Cp=1 and C'p+1=1,    AddOneUlp
 114  *  + if Cp=1 and C'p+1=-1,   Truncate
 115  *  + if Cp=1 and C'p+1=0,    Truncate if Ap-1=0, AddOneUlp else
 116  *  + if Cp=-1 and C'p+1=-1,  SubOneUlp
 117  *  + if Cp=-1 and C'p+1=0,   Truncate if Ap-1=0, SubOneUlp else
 118  *
 119  * If AddOneUlp:
 120  *   If carry, then it is 11111111111 + 1 = 10000000000000
 121  *      ap[n-1]=MPFR_HIGHT_BIT
 122  * If SubOneUlp:
 123  *   If we lose one bit, it is 1000000000 - 1 = 0111111111111
 124  *      Then shift, and put as last bit x which is calculated
 125  *              according Cp, Cp-1 and rnd_mode.
 126  * If Truncate,
 127  *    If it is a power of 2,
 128  *       we may have to suboneulp in some special cases.
 129  *
 130  * To simplify, we don't use Cp = 1.
 131  *
 132  */
 133
 134 int
 135 mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
 136 {
 137   mpfr_exp_t bx,cx;
 138   mpfr_uexp_t d;
 139   mpfr_prec_t p, sh, cnt;
 140   mp_size_t n;
 141   mp_limb_t *ap, *bp, *cp;
 142   mp_limb_t limb;
 143   int inexact;
 144   mp_limb_t bcp,bcp1; /* Cp and C'p+1 */
 145   mp_limb_t bbcp = (mp_limb_t) -1, bbcp1 = (mp_limb_t) -1; /* Cp+1 and C'p+2,
 146     gcc claims that they might be used uninitialized. We fill them with invalid
 147     values, which should produce a failure if so. See README.dev file. */
 148
 149   MPFR_TMP_DECL(marker);
 150
 151   MPFR_TMP_MARK(marker);
 152
 153   MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
 154   MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));
 155
 156   /* Read prec and num of limbs */
 157   p = MPFR_PREC(b);
 158   n = (p-1)/GMP_NUMB_BITS+1;
 159
 160   /* Fast cmp of |b| and |c|*/
 161   bx = MPFR_GET_EXP (b);
 162   cx = MPFR_GET_EXP (c);
 163   if (MPFR_UNLIKELY(bx == cx))
 164     {
 165       mp_size_t k = n - 1;
 166       /* Check mantissa since exponent are equals */
 167       bp = MPFR_MANT(b);
 168       cp = MPFR_MANT(c);
 169       while (k>=0 && MPFR_UNLIKELY(bp[k] == cp[k]))
 170         k--;
 171       if (MPFR_UNLIKELY(k < 0))
 172         /* b == c ! */
 173         {
 174           /* Return exact number 0 */
 175           if (rnd_mode == MPFR_RNDD)
 176             MPFR_SET_NEG(a);
 177           else
 178             MPFR_SET_POS(a);
 179           MPFR_SET_ZERO(a);
 180           MPFR_RET(0);
 181         }
 182       else if (bp[k] > cp[k])
 183         goto BGreater;
 184       else
 185         {
 186           MPFR_ASSERTD(bp[k]<cp[k]);
 187           goto CGreater;
 188         }
 189     }
 190   else if (MPFR_UNLIKELY(bx < cx))
 191     {
 192       /* Swap b and c and set sign */
 193       mpfr_srcptr t;
 194       mpfr_exp_t tx;
 195     CGreater:
 196       MPFR_SET_OPPOSITE_SIGN(a,b);
 197       t  = b;  b  = c;  c  = t;
 198       tx = bx; bx = cx; cx = tx;
 199     }
 200   else
 201     {
 202       /* b > c */
 203     BGreater:
 204       MPFR_SET_SAME_SIGN(a,b);
 205     }
 206
 207   /* Now b > c */
 208   MPFR_ASSERTD(bx >= cx);
 209   d = (mpfr_uexp_t) bx - cx;
 210   DEBUG (printf ("New with diff=%lu\n", (unsigned long) d));
 211
 212   if (MPFR_UNLIKELY(d <= 1))
 213     {
 214       if (MPFR_LIKELY(d < 1))
 215         {
 216           /* <-- b -->
 217              <-- c --> : exact sub */
 218           ap = MPFR_MANT(a);
 219           mpn_sub_n (ap, MPFR_MANT(b), MPFR_MANT(c), n);
 220           /* Normalize */
 221         ExactNormalize:
 222           limb = ap[n-1];
 223           if (MPFR_LIKELY(limb))
 224             {
 225               /* First limb is not zero. */
 226               count_leading_zeros(cnt, limb);
 227               /* cnt could be == 0 <= SubD1Lose */
 228               if (MPFR_LIKELY(cnt))
 229                 {
 230                   mpn_lshift(ap, ap, n, cnt); /* Normalize number */
 231                   bx -= cnt; /* Update final expo */
 232                 }
 233               /* Last limb should be ok */
 234               MPFR_ASSERTD(!(ap[0] & MPFR_LIMB_MASK((unsigned int) (-p)
 235                                                     % GMP_NUMB_BITS)));
 236             }
 237           else
 238             {
 239               /* First limb is zero */
 240               mp_size_t k = n-1, len;
 241               /* Find the first limb not equal to zero.
 242                  FIXME:It is assume it exists (since |b| > |c| and same prec)*/
 243               do
 244                 {
 245                   MPFR_ASSERTD( k > 0 );
 246                   limb = ap[--k];
 247                 }
 248               while (limb == 0);
 249               MPFR_ASSERTD(limb != 0);
 250               count_leading_zeros(cnt, limb);
 251               k++;
 252               len = n - k; /* Number of last limb */
 253               MPFR_ASSERTD(k >= 0);
 254               if (MPFR_LIKELY(cnt))
 255                 mpn_lshift(ap+len, ap, k, cnt); /* Normalize the High Limb*/
 256               else
 257                 {
 258                   /* Must use DECR since src and dest may overlap & dest>=src*/
 259                   MPN_COPY_DECR(ap+len, ap, k);
 260                 }
 261               MPN_ZERO(ap, len); /* Zeroing the last limbs */
 262               bx -= cnt + len*GMP_NUMB_BITS; /* Update Expo */
 263               /* Last limb should be ok */
 264               MPFR_ASSERTD(!(ap[len]&MPFR_LIMB_MASK((unsigned int) (-p)
 265                                                     % GMP_NUMB_BITS)));
 266             }
 267           /* Check expo underflow */
 268           if (MPFR_UNLIKELY(bx < __gmpfr_emin))
 269             {
 270               MPFR_TMP_FREE(marker);
 271               /* inexact=0 */
 272               DEBUG( printf("(D==0 Underflow)\n") );
 273               if (rnd_mode == MPFR_RNDN &&
 274                   (bx < __gmpfr_emin - 1 ||
 275                    (/*inexact >= 0 &&*/ mpfr_powerof2_raw (a))))
 276                 rnd_mode = MPFR_RNDZ;
 277               return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
 278             }
 279           MPFR_SET_EXP (a, bx);
 280           /* No rounding is necessary since the result is exact */
 281           MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
 282           MPFR_TMP_FREE(marker);
 283           return 0;
 284         }
 285       else /* if (d == 1) */
 286         {
 287           /* | <-- b -->
 288              |  <-- c --> */
 289           mp_limb_t c0, mask;
 290           mp_size_t k;
 291           MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 292           /* If we lose at least one bit, compute 2*b-c (Exact)
 293            * else compute b-c/2 */
 294           bp = MPFR_MANT(b);
 295           cp = MPFR_MANT(c);
 296           k = n-1;
 297           limb = bp[k] - cp[k]/2;
 298           if (limb > MPFR_LIMB_HIGHBIT)
 299             {
 300               /* We can't lose precision: compute b-c/2 */
 301               /* Shift c in the allocated temporary block */
 302             SubD1NoLose:
 303               c0 = cp[0] & (MPFR_LIMB_ONE<<sh);
 304               cp = (mp_limb_t*) MPFR_TMP_ALLOC(n * BYTES_PER_MP_LIMB);
 305               mpn_rshift(cp, MPFR_MANT(c), n, 1);
 306               if (MPFR_LIKELY(c0 == 0))
 307                 {
 308                   /* Result is exact: no need of rounding! */
 309                   ap = MPFR_MANT(a);
 310                   mpn_sub_n (ap, bp, cp, n);
 311                   MPFR_SET_EXP(a, bx); /* No expo overflow! */
 312                   /* No truncate or normalize is needed */
 313                   MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
 314                   /* No rounding is necessary since the result is exact */
 315                   MPFR_TMP_FREE(marker);
 316                   return 0;
 317                 }
 318               ap = MPFR_MANT(a);
 319               mask = ~MPFR_LIMB_MASK(sh);
 320               cp[0] &= mask; /* Delete last bit of c */
 321               mpn_sub_n (ap, bp, cp, n);
 322               MPFR_SET_EXP(a, bx);                 /* No expo overflow! */
 323               MPFR_ASSERTD( !(ap[0] & ~mask) );    /* Check last bits */
 324               /* No normalize is needed */
 325               MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
 326               /* Rounding is necessary since c0 = 1*/
 327               /* Cp =-1 and C'p+1=0 */
 328               bcp = 1; bcp1 = 0;
 329               if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 330                 {
 331                   /* Even Rule apply: Check Ap-1 */
 332                   if (MPFR_LIKELY( (ap[0] & (MPFR_LIMB_ONE<<sh)) == 0) )
 333                     goto truncate;
 334                   else
 335                     goto sub_one_ulp;
 336                 }
 337               MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 338               if (rnd_mode == MPFR_RNDZ)
 339                 goto sub_one_ulp;
 340               else
 341                 goto truncate;
 342             }
 343           else if (MPFR_LIKELY(limb < MPFR_LIMB_HIGHBIT))
 344             {
 345               /* We lose at least one bit of prec */
 346               /* Calcul of 2*b-c (Exact) */
 347               /* Shift b in the allocated temporary block */
 348             SubD1Lose:
 349               bp = (mp_limb_t*) MPFR_TMP_ALLOC (n * BYTES_PER_MP_LIMB);
 350               mpn_lshift (bp, MPFR_MANT(b), n, 1);
 351               ap = MPFR_MANT(a);
 352               mpn_sub_n (ap, bp, cp, n);
 353               bx--;
 354               goto ExactNormalize;
 355             }
 356           else
 357             {
 358               /* Case: limb = 100000000000 */
 359               /* Check while b[k] == c'[k] (C' is C shifted by 1) */
 360               /* If b[k]<c'[k] => We lose at least one bit*/
 361               /* If b[k]>c'[k] => We don't lose any bit */
 362               /* If k==-1 => We don't lose any bit
 363                  AND the result is 100000000000 0000000000 00000000000 */
 364               mp_limb_t carry;
 365               do {
 366                 carry = cp[k]&MPFR_LIMB_ONE;
 367                 k--;
 368               } while (k>=0 &&
 369                        bp[k]==(carry=cp[k]/2+(carry<<(GMP_NUMB_BITS-1))));
 370               if (MPFR_UNLIKELY(k<0))
 371                 {
 372                   /*If carry then (sh==0 and Virtual c'[-1] > Virtual b[-1]) */
 373                   if (MPFR_UNLIKELY(carry)) /* carry = cp[0]&MPFR_LIMB_ONE */
 374                     {
 375                       /* FIXME: Can be faster? */
 376                       MPFR_ASSERTD(sh == 0);
 377                       goto SubD1Lose;
 378                     }
 379                   /* Result is a power of 2 */
 380                   ap = MPFR_MANT (a);
 381                   MPN_ZERO (ap, n);
 382                   ap[n-1] = MPFR_LIMB_HIGHBIT;
 383                   MPFR_SET_EXP (a, bx); /* No expo overflow! */
 384                   /* No Normalize is needed*/
 385                   /* No Rounding is needed */
 386                   MPFR_TMP_FREE (marker);
 387                   return 0;
 388                 }
 389               /* carry = cp[k]/2+(cp[k-1]&1)<<(GMP_NUMB_BITS-1) = c'[k]*/
 390               else if (bp[k] > carry)
 391                 goto SubD1NoLose;
 392               else
 393                 {
 394                   MPFR_ASSERTD(bp[k]<carry);
 395                   goto SubD1Lose;
 396                 }
 397             }
 398         }
 399     }
 400   else if (MPFR_UNLIKELY(d >= p))
 401     {
 402       ap = MPFR_MANT(a);
 403       MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 404       /* We can't set A before since we use cp for rounding... */
 405       /* Perform rounding: check if a=b or a=b-ulp(b) */
 406       if (MPFR_UNLIKELY(d == p))
 407         {
 408           /* cp == -1 and c'p+1 = ? */
 409           bcp  = 1;
 410           /* We need Cp+1 later for a very improbable case. */
 411           bbcp = (MPFR_MANT(c)[n-1] & (MPFR_LIMB_ONE<<(GMP_NUMB_BITS-2)));
 412           /* We need also C'p+1 for an even more unprobable case... */
 413           if (MPFR_LIKELY( bbcp ))
 414             bcp1 = 1;
 415           else
 416             {
 417               cp = MPFR_MANT(c);
 418               if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
 419                 {
 420                   mp_size_t k = n-1;
 421                   do {
 422                     k--;
 423                   } while (k>=0 && cp[k]==0);
 424                   bcp1 = (k>=0);
 425                 }
 426               else
 427                 bcp1 = 1;
 428             }
 429           DEBUG( printf("(D=P) Cp=-1 Cp+1=%d C'p+1=%d \n", bbcp!=0, bcp1!=0) );
 430           bp = MPFR_MANT (b);
 431
 432           /* Even if src and dest overlap, it is ok using MPN_COPY */
 433           if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 434             {
 435               if (MPFR_UNLIKELY( bcp && bcp1==0 ))
 436                 /* Cp=-1 and C'p+1=0: Even rule Apply! */
 437                 /* Check Ap-1 = Bp-1 */
 438                 if ((bp[0] & (MPFR_LIMB_ONE<<sh)) == 0)
 439                   {
 440                     MPN_COPY(ap, bp, n);
 441                     goto truncate;
 442                   }
 443               MPN_COPY(ap, bp, n);
 444               goto sub_one_ulp;
 445             }
 446           MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 447           if (rnd_mode == MPFR_RNDZ)
 448             {
 449               MPN_COPY(ap, bp, n);
 450               goto sub_one_ulp;
 451             }
 452           else
 453             {
 454               MPN_COPY(ap, bp, n);
 455               goto truncate;
 456             }
 457         }
 458       else
 459         {
 460           /* Cp=0, Cp+1=-1 if d==p+1, C'p+1=-1 */
 461           bcp = 0; bbcp = (d==p+1); bcp1 = 1;
 462           DEBUG( printf("(D>P) Cp=%d Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0) );
 463           /* Need to compute C'p+2 if d==p+1 and if rnd_mode=NEAREST
 464              (Because of a very improbable case) */
 465           if (MPFR_UNLIKELY(d==p+1 && rnd_mode==MPFR_RNDN))
 466             {
 467               cp = MPFR_MANT(c);
 468               if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
 469                 {
 470                   mp_size_t k = n-1;
 471                   do {
 472                     k--;
 473                   } while (k>=0 && cp[k]==0);
 474                   bbcp1 = (k>=0);
 475                 }
 476               else
 477                 bbcp1 = 1;
 478               DEBUG( printf("(D>P) C'p+2=%d\n", bbcp1!=0) );
 479             }
 480           /* Copy mantissa B in A */
 481           MPN_COPY(ap, MPFR_MANT(b), n);
 482           /* Round */
 483           if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 484             goto truncate;
 485           MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 486           if (rnd_mode == MPFR_RNDZ)
 487             goto sub_one_ulp;
 488           else /* rnd_mode = AWAY */
 489             goto truncate;
 490         }
 491     }
 492   else
 493     {
 494       mpfr_uexp_t dm;
 495       mp_size_t m;
 496       mp_limb_t mask;
 497
 498       /* General case: 2 <= d < p */
 499       MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 500       cp = (mp_limb_t*) MPFR_TMP_ALLOC(n * BYTES_PER_MP_LIMB);
 501
 502       /* Shift c in temporary allocated place */
 503       dm = d % GMP_NUMB_BITS;
 504       m = d / GMP_NUMB_BITS;
 505       if (MPFR_UNLIKELY(dm == 0))
 506         {
 507           /* dm = 0 and m > 0: Just copy */
 508           MPFR_ASSERTD(m!=0);
 509           MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
 510           MPN_ZERO(cp+n-m, m);
 511         }
 512       else if (MPFR_LIKELY(m == 0))
 513         {
 514           /* dm >=2 and m == 0: just shift */
 515           MPFR_ASSERTD(dm >= 2);
 516           mpn_rshift(cp, MPFR_MANT(c), n, dm);
 517         }
 518       else
 519         {
 520           /* dm > 0 and m > 0: shift and zero  */
 521           mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
 522           MPN_ZERO(cp+n-m, m);
 523         }
 524
 525       DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
 526       DEBUG( mpfr_print_mant_binary("B=    ", MPFR_MANT(b), p) );
 527       DEBUG( mpfr_print_mant_binary("After ", cp, p) );
 528
 529       /* Compute bcp=Cp and bcp1=C'p+1 */
 530       if (MPFR_LIKELY(sh))
 531         {
 532           /* Try to compute them from C' rather than C (FIXME: Faster?) */
 533           bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
 534           if (MPFR_LIKELY( cp[0] & MPFR_LIMB_MASK(sh-1) ))
 535             bcp1 = 1;
 536           else
 537             {
 538               /* We can't compute C'p+1 from C'. Compute it from C */
 539               /* Start from bit x=p-d+sh in mantissa C
 540                  (+sh since we have already looked sh bits in C'!) */
 541               mpfr_prec_t x = p-d+sh-1;
 542               if (MPFR_LIKELY(x>p))
 543                 /* We are already looked at all the bits of c, so C'p+1 = 0*/
 544                 bcp1 = 0;
 545               else
 546                 {
 547                   mp_limb_t *tp = MPFR_MANT(c);
 548                   mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
 549                   mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 550                   DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
 551                                  (unsigned long) x, (long) kx,
 552                                  (unsigned long) sx));
 553                   /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
 554                   if (tp[kx] & MPFR_LIMB_MASK(sx))
 555                     bcp1 = 1;
 556                   else
 557                     {
 558                       /*kx += (sx==0);*/
 559                       /*If sx==0, tp[kx] hasn't been checked*/
 560                       do {
 561                         kx--;
 562                       } while (kx>=0 && tp[kx]==0);
 563                       bcp1 = (kx >= 0);
 564                     }
 565                 }
 566             }
 567         }
 568       else
 569         {
 570           /* Compute Cp and C'p+1 from C with sh=0 */
 571           mp_limb_t *tp = MPFR_MANT(c);
 572           /* Start from bit x=p-d in mantissa C */
 573           mpfr_prec_t  x = p-d;
 574           mp_size_t   kx = n-1 - (x / GMP_NUMB_BITS);
 575           mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 576           MPFR_ASSERTD(p >= d);
 577           bcp = (tp[kx] & (MPFR_LIMB_ONE<<sx));
 578           /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
 579           if (tp[kx] & MPFR_LIMB_MASK(sx))
 580             bcp1 = 1;
 581           else
 582             {
 583               /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
 584               do {
 585                 kx--;
 586               } while (kx>=0 && tp[kx]==0);
 587               bcp1 = (kx>=0);
 588             }
 589         }
 590       DEBUG( printf("sh=%lu Cp=%d C'p+1=%d\n", sh, bcp!=0, bcp1!=0) );
 591
 592       /* Check if we can lose a bit, and if so compute Cp+1 and C'p+2 */
 593       bp = MPFR_MANT(b);
 594       if (MPFR_UNLIKELY((bp[n-1]-cp[n-1]) <= MPFR_LIMB_HIGHBIT))
 595         {
 596           /* We can lose a bit so we precompute Cp+1 and C'p+2 */
 597           /* Test for trivial case: since C'p+1=0, Cp+1=0 and C'p+2 =0 */
 598           if (MPFR_LIKELY(bcp1 == 0))
 599             {
 600               bbcp = 0;
 601               bbcp1 = 0;
 602             }
 603           else /* bcp1 != 0 */
 604             {
 605               /* We can lose a bit:
 606                  compute Cp+1 and C'p+2 from mantissa C */
 607               mp_limb_t *tp = MPFR_MANT(c);
 608               /* Start from bit x=(p+1)-d in mantissa C */
 609               mpfr_prec_t x  = p+1-d;
 610               mp_size_t kx = n-1 - (x/GMP_NUMB_BITS);
 611               mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 612               MPFR_ASSERTD(p > d);
 613               DEBUG (printf ("(pre) x=%lu Kx=%ld Sx=%lu\n",
 614                              (unsigned long) x, (long) kx,
 615                              (unsigned long) sx));
 616               bbcp = (tp[kx] & (MPFR_LIMB_ONE<<sx)) ;
 617               /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
 618               /* If Cp+1=0, since C'p+1!=0, C'p+2=1 ! */
 619               if (MPFR_LIKELY(bbcp==0 || (tp[kx]&MPFR_LIMB_MASK(sx))))
 620                 bbcp1 = 1;
 621               else
 622                 {
 623                   /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
 624                   do {
 625                     kx--;
 626                   } while (kx>=0 && tp[kx]==0);
 627                   bbcp1 = (kx>=0);
 628                   DEBUG (printf ("(Pre) Scan done for %ld\n", (long) kx));
 629                 }
 630             } /*End of Bcp1 != 0*/
 631           DEBUG( printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp!=0, bbcp1!=0) );
 632         } /* End of "can lose a bit" */
 633
 634       /* Clean shifted C' */
 635       mask = ~MPFR_LIMB_MASK (sh);
 636       cp[0] &= mask;
 637
 638       /* Substract the mantissa c from b in a */
 639       ap = MPFR_MANT(a);
 640       mpn_sub_n (ap, bp, cp, n);
 641       DEBUG( mpfr_print_mant_binary("Sub=  ", ap, p) );
 642
 643      /* Normalize: we lose at max one bit*/
 644       if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
 645         {
 646           /* High bit is not set and we have to fix it! */
 647           /* Ap >= 010000xxx001 */
 648           mpn_lshift(ap, ap, n, 1);
 649           /* Ap >= 100000xxx010 */
 650           if (MPFR_UNLIKELY(bcp!=0)) /* Check if Cp = -1 */
 651             /* Since Cp == -1, we have to substract one more */
 652             {
 653               mpn_sub_1(ap, ap, n, MPFR_LIMB_ONE<<sh);
 654               MPFR_ASSERTD(MPFR_LIMB_MSB(ap[n-1]) != 0);
 655             }
 656           /* Ap >= 10000xxx001 */
 657           /* Final exponent -1 since we have shifted the mantissa */
 658           bx--;
 659           /* Update bcp and bcp1 */
 660           MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
 661           MPFR_ASSERTN(bbcp1 != (mp_limb_t) -1);
 662           bcp  = bbcp;
 663           bcp1 = bbcp1;
 664           /* We dont't have anymore a valid Cp+1!
 665              But since Ap >= 100000xxx001, the final sub can't unnormalize!*/
 666         }
 667       MPFR_ASSERTD( !(ap[0] & ~mask) );
 668
 669       /* Rounding */
 670       if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 671         {
 672           if (MPFR_LIKELY(bcp==0))
 673             goto truncate;
 674           else if ((bcp1) || ((ap[0] & (MPFR_LIMB_ONE<<sh)) != 0))
 675             goto sub_one_ulp;
 676           else
 677             goto truncate;
 678         }
 679
 680       /* Update rounding mode */
 681       MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 682       if (rnd_mode == MPFR_RNDZ && (MPFR_LIKELY(bcp || bcp1)))
 683         goto sub_one_ulp;
 684       goto truncate;
 685     }
 686   MPFR_RET_NEVER_GO_HERE ();
 687
 688   /* Sub one ulp to the result */
 689  sub_one_ulp:
 690   mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
 691   /* Result should be smaller than exact value: inexact=-1 */
 692   inexact = -1;
 693   /* Check normalisation */
 694   if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
 695     {
 696       /* ap was a power of 2, and we lose a bit */
 697       /* Now it is 0111111111111111111[00000 */
 698       mpn_lshift(ap, ap, n, 1);
 699       bx--;
 700       /* And the lost bit x depends on Cp+1, and Cp */
 701       /* Compute Cp+1 if it isn't already compute (ie d==1) */
 702       /* FIXME: Is this case possible? */
 703       if (MPFR_UNLIKELY(d == 1))
 704         bbcp = 0;
 705       DEBUG( printf("(SubOneUlp)Cp=%d, Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0));
 706       /* Compute the last bit (Since we have shifted the mantissa)
 707          we need one more bit!*/
 708       MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
 709       if ( (rnd_mode == MPFR_RNDZ && bcp==0)
 710            || (rnd_mode==MPFR_RNDN && bbcp==0)
 711            || (bcp && bcp1==0) ) /*Exact result*/
 712         {
 713           ap[0] |= MPFR_LIMB_ONE<<sh;
 714           if (rnd_mode == MPFR_RNDN)
 715             inexact = 1;
 716           DEBUG( printf("(SubOneUlp) Last bit set\n") );
 717         }
 718       /* Result could be exact if C'p+1 = 0 and rnd == Zero
 719          since we have had one more bit to the result */
 720       /* Fixme: rnd_mode == MPFR_RNDZ needed ? */
 721       if (bcp1==0 && rnd_mode==MPFR_RNDZ)
 722         {
 723           DEBUG( printf("(SubOneUlp) Exact result\n") );
 724           inexact = 0;
 725         }
 726     }
 727
 728   goto end_of_sub;
 729
 730  truncate:
 731   /* Check if the result is an exact power of 2: 100000000000
 732      in which cases, we could have to do sub_one_ulp due to some nasty reasons:
 733      If Result is a Power of 2:
 734       + If rnd = AWAY,
 735       |  If Cp=-1 and C'p+1 = 0, SubOneUlp and the result is EXACT.
 736          If Cp=-1 and C'p+1 =-1, SubOneUlp and the result is above.
 737          Otherwise truncate
 738       + If rnd = NEAREST,
 739          If Cp= 0 and Cp+1  =-1 and C'p+2=-1, SubOneUlp and the result is above
 740          If cp=-1 and C'p+1 = 0, SubOneUlp and the result is exact.
 741          Otherwise truncate.
 742       X bit should always be set if SubOneUlp*/
 743   if (MPFR_UNLIKELY(ap[n-1] == MPFR_LIMB_HIGHBIT))
 744     {
 745       mp_size_t k = n-1;
 746       do {
 747         k--;
 748       } while (k>=0 && ap[k]==0);
 749       if (MPFR_UNLIKELY(k<0))
 750         {
 751           /* It is a power of 2! */
 752           /* Compute Cp+1 if it isn't already compute (ie d==1) */
 753           /* FIXME: Is this case possible? */
 754           if (d == 1)
 755             bbcp=0;
 756           DEBUG( printf("(Truncate) Cp=%d, Cp+1=%d C'p+1=%d C'p+2=%d\n", \
 757                  bcp!=0, bbcp!=0, bcp1!=0, bbcp1!=0) );
 758           MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
 759           MPFR_ASSERTN((rnd_mode != MPFR_RNDN) || (bcp != 0) || (bbcp == 0) || (bbcp1 != (mp_limb_t) -1));
 760           if (((rnd_mode != MPFR_RNDZ) && bcp)
 761               ||
 762               ((rnd_mode == MPFR_RNDN) && (bcp == 0) && (bbcp) && (bbcp1)))
 763             {
 764               DEBUG( printf("(Truncate) Do sub\n") );
 765               mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
 766               mpn_lshift(ap, ap, n, 1);
 767               ap[0] |= MPFR_LIMB_ONE<<sh;
 768               bx--;
 769               /* FIXME: Explain why it works (or why not)... */
 770               inexact = (bcp1 == 0) ? 0 : (rnd_mode==MPFR_RNDN) ? -1 : 1;
 771               goto end_of_sub;
 772             }
 773         }
 774     }
 775
 776   /* Calcul of Inexact flag.*/
 777   inexact = MPFR_LIKELY(bcp || bcp1) ? 1 : 0;
 778
 779  end_of_sub:
 780   /* Update Expo */
 781   /* FIXME: Is this test really useful?
 782       If d==0      : Exact case. This is never called.
 783       if 1 < d < p : bx=MPFR_EXP(b) or MPFR_EXP(b)-1 > MPFR_EXP(c) > emin
 784       if d == 1    : bx=MPFR_EXP(b). If we could lose any bits, the exact
 785                      normalisation is called.
 786       if d >=  p   : bx=MPFR_EXP(b) >= MPFR_EXP(c) + p > emin
 787      After SubOneUlp, we could have one bit less.
 788       if 1 < d < p : bx >= MPFR_EXP(b)-2 >= MPFR_EXP(c) > emin
 789       if d == 1    : bx >= MPFR_EXP(b)-1 = MPFR_EXP(c) > emin.
 790       if d >=  p   : bx >= MPFR_EXP(b)-1 > emin since p>=2.
 791   */
 792   MPFR_ASSERTD( bx >= __gmpfr_emin);
 793   /*
 794     if (MPFR_UNLIKELY(bx < __gmpfr_emin))
 795     {
 796       DEBUG( printf("(Final Underflow)\n") );
 797       if (rnd_mode == MPFR_RNDN &&
 798           (bx < __gmpfr_emin - 1 ||
 799            (inexact >= 0 && mpfr_powerof2_raw (a))))
 800         rnd_mode = MPFR_RNDZ;
 801       MPFR_TMP_FREE(marker);
 802       return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
 803     }
 804   */
 805   MPFR_SET_EXP (a, bx);
 806
 807   MPFR_TMP_FREE(marker);
 808   MPFR_RET (inexact * MPFR_INT_SIGN (a));
 809 }