contrib/mpfr/src/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, 2011 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Caramel 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));
 155   MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
 156
 157   /* Read prec and num of limbs */
 158   p = MPFR_PREC(b);
 159   n = (p-1)/GMP_NUMB_BITS+1;
 160
 161   /* Fast cmp of |b| and |c|*/
 162   bx = MPFR_GET_EXP (b);
 163   cx = MPFR_GET_EXP (c);
 164   if (MPFR_UNLIKELY(bx == cx))
 165     {
 166       mp_size_t k = n - 1;
 167       /* Check mantissa since exponent are equals */
 168       bp = MPFR_MANT(b);
 169       cp = MPFR_MANT(c);
 170       while (k>=0 && MPFR_UNLIKELY(bp[k] == cp[k]))
 171         k--;
 172       if (MPFR_UNLIKELY(k < 0))
 173         /* b == c ! */
 174         {
 175           /* Return exact number 0 */
 176           if (rnd_mode == MPFR_RNDD)
 177             MPFR_SET_NEG(a);
 178           else
 179             MPFR_SET_POS(a);
 180           MPFR_SET_ZERO(a);
 181           MPFR_RET(0);
 182         }
 183       else if (bp[k] > cp[k])
 184         goto BGreater;
 185       else
 186         {
 187           MPFR_ASSERTD(bp[k]<cp[k]);
 188           goto CGreater;
 189         }
 190     }
 191   else if (MPFR_UNLIKELY(bx < cx))
 192     {
 193       /* Swap b and c and set sign */
 194       mpfr_srcptr t;
 195       mpfr_exp_t tx;
 196     CGreater:
 197       MPFR_SET_OPPOSITE_SIGN(a,b);
 198       t  = b;  b  = c;  c  = t;
 199       tx = bx; bx = cx; cx = tx;
 200     }
 201   else
 202     {
 203       /* b > c */
 204     BGreater:
 205       MPFR_SET_SAME_SIGN(a,b);
 206     }
 207
 208   /* Now b > c */
 209   MPFR_ASSERTD(bx >= cx);
 210   d = (mpfr_uexp_t) bx - cx;
 211   DEBUG (printf ("New with diff=%lu\n", (unsigned long) d));
 212
 213   if (MPFR_UNLIKELY(d <= 1))
 214     {
 215       if (MPFR_LIKELY(d < 1))
 216         {
 217           /* <-- b -->
 218              <-- c --> : exact sub */
 219           ap = MPFR_MANT(a);
 220           mpn_sub_n (ap, MPFR_MANT(b), MPFR_MANT(c), n);
 221           /* Normalize */
 222         ExactNormalize:
 223           limb = ap[n-1];
 224           if (MPFR_LIKELY(limb))
 225             {
 226               /* First limb is not zero. */
 227               count_leading_zeros(cnt, limb);
 228               /* cnt could be == 0 <= SubD1Lose */
 229               if (MPFR_LIKELY(cnt))
 230                 {
 231                   mpn_lshift(ap, ap, n, cnt); /* Normalize number */
 232                   bx -= cnt; /* Update final expo */
 233                 }
 234               /* Last limb should be ok */
 235               MPFR_ASSERTD(!(ap[0] & MPFR_LIMB_MASK((unsigned int) (-p)
 236                                                     % GMP_NUMB_BITS)));
 237             }
 238           else
 239             {
 240               /* First limb is zero */
 241               mp_size_t k = n-1, len;
 242               /* Find the first limb not equal to zero.
 243                  FIXME:It is assume it exists (since |b| > |c| and same prec)*/
 244               do
 245                 {
 246                   MPFR_ASSERTD( k > 0 );
 247                   limb = ap[--k];
 248                 }
 249               while (limb == 0);
 250               MPFR_ASSERTD(limb != 0);
 251               count_leading_zeros(cnt, limb);
 252               k++;
 253               len = n - k; /* Number of last limb */
 254               MPFR_ASSERTD(k >= 0);
 255               if (MPFR_LIKELY(cnt))
 256                 mpn_lshift(ap+len, ap, k, cnt); /* Normalize the High Limb*/
 257               else
 258                 {
 259                   /* Must use DECR since src and dest may overlap & dest>=src*/
 260                   MPN_COPY_DECR(ap+len, ap, k);
 261                 }
 262               MPN_ZERO(ap, len); /* Zeroing the last limbs */
 263               bx -= cnt + len*GMP_NUMB_BITS; /* Update Expo */
 264               /* Last limb should be ok */
 265               MPFR_ASSERTD(!(ap[len]&MPFR_LIMB_MASK((unsigned int) (-p)
 266                                                     % GMP_NUMB_BITS)));
 267             }
 268           /* Check expo underflow */
 269           if (MPFR_UNLIKELY(bx < __gmpfr_emin))
 270             {
 271               MPFR_TMP_FREE(marker);
 272               /* inexact=0 */
 273               DEBUG( printf("(D==0 Underflow)\n") );
 274               if (rnd_mode == MPFR_RNDN &&
 275                   (bx < __gmpfr_emin - 1 ||
 276                    (/*inexact >= 0 &&*/ mpfr_powerof2_raw (a))))
 277                 rnd_mode = MPFR_RNDZ;
 278               return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
 279             }
 280           MPFR_SET_EXP (a, bx);
 281           /* No rounding is necessary since the result is exact */
 282           MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
 283           MPFR_TMP_FREE(marker);
 284           return 0;
 285         }
 286       else /* if (d == 1) */
 287         {
 288           /* | <-- b -->
 289              |  <-- c --> */
 290           mp_limb_t c0, mask;
 291           mp_size_t k;
 292           MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 293           /* If we lose at least one bit, compute 2*b-c (Exact)
 294            * else compute b-c/2 */
 295           bp = MPFR_MANT(b);
 296           cp = MPFR_MANT(c);
 297           k = n-1;
 298           limb = bp[k] - cp[k]/2;
 299           if (limb > MPFR_LIMB_HIGHBIT)
 300             {
 301               /* We can't lose precision: compute b-c/2 */
 302               /* Shift c in the allocated temporary block */
 303             SubD1NoLose:
 304               c0 = cp[0] & (MPFR_LIMB_ONE<<sh);
 305               cp = MPFR_TMP_LIMBS_ALLOC (n);
 306               mpn_rshift(cp, MPFR_MANT(c), n, 1);
 307               if (MPFR_LIKELY(c0 == 0))
 308                 {
 309                   /* Result is exact: no need of rounding! */
 310                   ap = MPFR_MANT(a);
 311                   mpn_sub_n (ap, bp, cp, n);
 312                   MPFR_SET_EXP(a, bx); /* No expo overflow! */
 313                   /* No truncate or normalize is needed */
 314                   MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
 315                   /* No rounding is necessary since the result is exact */
 316                   MPFR_TMP_FREE(marker);
 317                   return 0;
 318                 }
 319               ap = MPFR_MANT(a);
 320               mask = ~MPFR_LIMB_MASK(sh);
 321               cp[0] &= mask; /* Delete last bit of c */
 322               mpn_sub_n (ap, bp, cp, n);
 323               MPFR_SET_EXP(a, bx);                 /* No expo overflow! */
 324               MPFR_ASSERTD( !(ap[0] & ~mask) );    /* Check last bits */
 325               /* No normalize is needed */
 326               MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
 327               /* Rounding is necessary since c0 = 1*/
 328               /* Cp =-1 and C'p+1=0 */
 329               bcp = 1; bcp1 = 0;
 330               if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 331                 {
 332                   /* Even Rule apply: Check Ap-1 */
 333                   if (MPFR_LIKELY( (ap[0] & (MPFR_LIMB_ONE<<sh)) == 0) )
 334                     goto truncate;
 335                   else
 336                     goto sub_one_ulp;
 337                 }
 338               MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 339               if (rnd_mode == MPFR_RNDZ)
 340                 goto sub_one_ulp;
 341               else
 342                 goto truncate;
 343             }
 344           else if (MPFR_LIKELY(limb < MPFR_LIMB_HIGHBIT))
 345             {
 346               /* We lose at least one bit of prec */
 347               /* Calcul of 2*b-c (Exact) */
 348               /* Shift b in the allocated temporary block */
 349             SubD1Lose:
 350               bp = MPFR_TMP_LIMBS_ALLOC (n);
 351               mpn_lshift (bp, MPFR_MANT(b), n, 1);
 352               ap = MPFR_MANT(a);
 353               mpn_sub_n (ap, bp, cp, n);
 354               bx--;
 355               goto ExactNormalize;
 356             }
 357           else
 358             {
 359               /* Case: limb = 100000000000 */
 360               /* Check while b[k] == c'[k] (C' is C shifted by 1) */
 361               /* If b[k]<c'[k] => We lose at least one bit*/
 362               /* If b[k]>c'[k] => We don't lose any bit */
 363               /* If k==-1 => We don't lose any bit
 364                  AND the result is 100000000000 0000000000 00000000000 */
 365               mp_limb_t carry;
 366               do {
 367                 carry = cp[k]&MPFR_LIMB_ONE;
 368                 k--;
 369               } while (k>=0 &&
 370                        bp[k]==(carry=cp[k]/2+(carry<<(GMP_NUMB_BITS-1))));
 371               if (MPFR_UNLIKELY(k<0))
 372                 {
 373                   /*If carry then (sh==0 and Virtual c'[-1] > Virtual b[-1]) */
 374                   if (MPFR_UNLIKELY(carry)) /* carry = cp[0]&MPFR_LIMB_ONE */
 375                     {
 376                       /* FIXME: Can be faster? */
 377                       MPFR_ASSERTD(sh == 0);
 378                       goto SubD1Lose;
 379                     }
 380                   /* Result is a power of 2 */
 381                   ap = MPFR_MANT (a);
 382                   MPN_ZERO (ap, n);
 383                   ap[n-1] = MPFR_LIMB_HIGHBIT;
 384                   MPFR_SET_EXP (a, bx); /* No expo overflow! */
 385                   /* No Normalize is needed*/
 386                   /* No Rounding is needed */
 387                   MPFR_TMP_FREE (marker);
 388                   return 0;
 389                 }
 390               /* carry = cp[k]/2+(cp[k-1]&1)<<(GMP_NUMB_BITS-1) = c'[k]*/
 391               else if (bp[k] > carry)
 392                 goto SubD1NoLose;
 393               else
 394                 {
 395                   MPFR_ASSERTD(bp[k]<carry);
 396                   goto SubD1Lose;
 397                 }
 398             }
 399         }
 400     }
 401   else if (MPFR_UNLIKELY(d >= p))
 402     {
 403       ap = MPFR_MANT(a);
 404       MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 405       /* We can't set A before since we use cp for rounding... */
 406       /* Perform rounding: check if a=b or a=b-ulp(b) */
 407       if (MPFR_UNLIKELY(d == p))
 408         {
 409           /* cp == -1 and c'p+1 = ? */
 410           bcp  = 1;
 411           /* We need Cp+1 later for a very improbable case. */
 412           bbcp = (MPFR_MANT(c)[n-1] & (MPFR_LIMB_ONE<<(GMP_NUMB_BITS-2)));
 413           /* We need also C'p+1 for an even more unprobable case... */
 414           if (MPFR_LIKELY( bbcp ))
 415             bcp1 = 1;
 416           else
 417             {
 418               cp = MPFR_MANT(c);
 419               if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
 420                 {
 421                   mp_size_t k = n-1;
 422                   do {
 423                     k--;
 424                   } while (k>=0 && cp[k]==0);
 425                   bcp1 = (k>=0);
 426                 }
 427               else
 428                 bcp1 = 1;
 429             }
 430           DEBUG( printf("(D=P) Cp=-1 Cp+1=%d C'p+1=%d \n", bbcp!=0, bcp1!=0) );
 431           bp = MPFR_MANT (b);
 432
 433           /* Even if src and dest overlap, it is ok using MPN_COPY */
 434           if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 435             {
 436               if (MPFR_UNLIKELY( bcp && bcp1==0 ))
 437                 /* Cp=-1 and C'p+1=0: Even rule Apply! */
 438                 /* Check Ap-1 = Bp-1 */
 439                 if ((bp[0] & (MPFR_LIMB_ONE<<sh)) == 0)
 440                   {
 441                     MPN_COPY(ap, bp, n);
 442                     goto truncate;
 443                   }
 444               MPN_COPY(ap, bp, n);
 445               goto sub_one_ulp;
 446             }
 447           MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 448           if (rnd_mode == MPFR_RNDZ)
 449             {
 450               MPN_COPY(ap, bp, n);
 451               goto sub_one_ulp;
 452             }
 453           else
 454             {
 455               MPN_COPY(ap, bp, n);
 456               goto truncate;
 457             }
 458         }
 459       else
 460         {
 461           /* Cp=0, Cp+1=-1 if d==p+1, C'p+1=-1 */
 462           bcp = 0; bbcp = (d==p+1); bcp1 = 1;
 463           DEBUG( printf("(D>P) Cp=%d Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0) );
 464           /* Need to compute C'p+2 if d==p+1 and if rnd_mode=NEAREST
 465              (Because of a very improbable case) */
 466           if (MPFR_UNLIKELY(d==p+1 && rnd_mode==MPFR_RNDN))
 467             {
 468               cp = MPFR_MANT(c);
 469               if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
 470                 {
 471                   mp_size_t k = n-1;
 472                   do {
 473                     k--;
 474                   } while (k>=0 && cp[k]==0);
 475                   bbcp1 = (k>=0);
 476                 }
 477               else
 478                 bbcp1 = 1;
 479               DEBUG( printf("(D>P) C'p+2=%d\n", bbcp1!=0) );
 480             }
 481           /* Copy mantissa B in A */
 482           MPN_COPY(ap, MPFR_MANT(b), n);
 483           /* Round */
 484           if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 485             goto truncate;
 486           MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 487           if (rnd_mode == MPFR_RNDZ)
 488             goto sub_one_ulp;
 489           else /* rnd_mode = AWAY */
 490             goto truncate;
 491         }
 492     }
 493   else
 494     {
 495       mpfr_uexp_t dm;
 496       mp_size_t m;
 497       mp_limb_t mask;
 498
 499       /* General case: 2 <= d < p */
 500       MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 501       cp = MPFR_TMP_LIMBS_ALLOC (n);
 502
 503       /* Shift c in temporary allocated place */
 504       dm = d % GMP_NUMB_BITS;
 505       m = d / GMP_NUMB_BITS;
 506       if (MPFR_UNLIKELY(dm == 0))
 507         {
 508           /* dm = 0 and m > 0: Just copy */
 509           MPFR_ASSERTD(m!=0);
 510           MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
 511           MPN_ZERO(cp+n-m, m);
 512         }
 513       else if (MPFR_LIKELY(m == 0))
 514         {
 515           /* dm >=2 and m == 0: just shift */
 516           MPFR_ASSERTD(dm >= 2);
 517           mpn_rshift(cp, MPFR_MANT(c), n, dm);
 518         }
 519       else
 520         {
 521           /* dm > 0 and m > 0: shift and zero  */
 522           mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
 523           MPN_ZERO(cp+n-m, m);
 524         }
 525
 526       DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
 527       DEBUG( mpfr_print_mant_binary("B=    ", MPFR_MANT(b), p) );
 528       DEBUG( mpfr_print_mant_binary("After ", cp, p) );
 529
 530       /* Compute bcp=Cp and bcp1=C'p+1 */
 531       if (MPFR_LIKELY(sh))
 532         {
 533           /* Try to compute them from C' rather than C (FIXME: Faster?) */
 534           bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
 535           if (MPFR_LIKELY( cp[0] & MPFR_LIMB_MASK(sh-1) ))
 536             bcp1 = 1;
 537           else
 538             {
 539               /* We can't compute C'p+1 from C'. Compute it from C */
 540               /* Start from bit x=p-d+sh in mantissa C
 541                  (+sh since we have already looked sh bits in C'!) */
 542               mpfr_prec_t x = p-d+sh-1;
 543               if (MPFR_LIKELY(x>p))
 544                 /* We are already looked at all the bits of c, so C'p+1 = 0*/
 545                 bcp1 = 0;
 546               else
 547                 {
 548                   mp_limb_t *tp = MPFR_MANT(c);
 549                   mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
 550                   mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 551                   DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
 552                                  (unsigned long) x, (long) kx,
 553                                  (unsigned long) sx));
 554                   /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
 555                   if (tp[kx] & MPFR_LIMB_MASK(sx))
 556                     bcp1 = 1;
 557                   else
 558                     {
 559                       /*kx += (sx==0);*/
 560                       /*If sx==0, tp[kx] hasn't been checked*/
 561                       do {
 562                         kx--;
 563                       } while (kx>=0 && tp[kx]==0);
 564                       bcp1 = (kx >= 0);
 565                     }
 566                 }
 567             }
 568         }
 569       else
 570         {
 571           /* Compute Cp and C'p+1 from C with sh=0 */
 572           mp_limb_t *tp = MPFR_MANT(c);
 573           /* Start from bit x=p-d in mantissa C */
 574           mpfr_prec_t  x = p-d;
 575           mp_size_t   kx = n-1 - (x / GMP_NUMB_BITS);
 576           mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 577           MPFR_ASSERTD(p >= d);
 578           bcp = (tp[kx] & (MPFR_LIMB_ONE<<sx));
 579           /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
 580           if (tp[kx] & MPFR_LIMB_MASK(sx))
 581             bcp1 = 1;
 582           else
 583             {
 584               /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
 585               do {
 586                 kx--;
 587               } while (kx>=0 && tp[kx]==0);
 588               bcp1 = (kx>=0);
 589             }
 590         }
 591       DEBUG( printf("sh=%lu Cp=%d C'p+1=%d\n", sh, bcp!=0, bcp1!=0) );
 592
 593       /* Check if we can lose a bit, and if so compute Cp+1 and C'p+2 */
 594       bp = MPFR_MANT(b);
 595       if (MPFR_UNLIKELY((bp[n-1]-cp[n-1]) <= MPFR_LIMB_HIGHBIT))
 596         {
 597           /* We can lose a bit so we precompute Cp+1 and C'p+2 */
 598           /* Test for trivial case: since C'p+1=0, Cp+1=0 and C'p+2 =0 */
 599           if (MPFR_LIKELY(bcp1 == 0))
 600             {
 601               bbcp = 0;
 602               bbcp1 = 0;
 603             }
 604           else /* bcp1 != 0 */
 605             {
 606               /* We can lose a bit:
 607                  compute Cp+1 and C'p+2 from mantissa C */
 608               mp_limb_t *tp = MPFR_MANT(c);
 609               /* Start from bit x=(p+1)-d in mantissa C */
 610               mpfr_prec_t x  = p+1-d;
 611               mp_size_t kx = n-1 - (x/GMP_NUMB_BITS);
 612               mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 613               MPFR_ASSERTD(p > d);
 614               DEBUG (printf ("(pre) x=%lu Kx=%ld Sx=%lu\n",
 615                              (unsigned long) x, (long) kx,
 616                              (unsigned long) sx));
 617               bbcp = (tp[kx] & (MPFR_LIMB_ONE<<sx)) ;
 618               /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
 619               /* If Cp+1=0, since C'p+1!=0, C'p+2=1 ! */
 620               if (MPFR_LIKELY(bbcp==0 || (tp[kx]&MPFR_LIMB_MASK(sx))))
 621                 bbcp1 = 1;
 622               else
 623                 {
 624                   /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
 625                   do {
 626                     kx--;
 627                   } while (kx>=0 && tp[kx]==0);
 628                   bbcp1 = (kx>=0);
 629                   DEBUG (printf ("(Pre) Scan done for %ld\n", (long) kx));
 630                 }
 631             } /*End of Bcp1 != 0*/
 632           DEBUG( printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp!=0, bbcp1!=0) );
 633         } /* End of "can lose a bit" */
 634
 635       /* Clean shifted C' */
 636       mask = ~MPFR_LIMB_MASK (sh);
 637       cp[0] &= mask;
 638
 639       /* Substract the mantissa c from b in a */
 640       ap = MPFR_MANT(a);
 641       mpn_sub_n (ap, bp, cp, n);
 642       DEBUG( mpfr_print_mant_binary("Sub=  ", ap, p) );
 643
 644      /* Normalize: we lose at max one bit*/
 645       if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
 646         {
 647           /* High bit is not set and we have to fix it! */
 648           /* Ap >= 010000xxx001 */
 649           mpn_lshift(ap, ap, n, 1);
 650           /* Ap >= 100000xxx010 */
 651           if (MPFR_UNLIKELY(bcp!=0)) /* Check if Cp = -1 */
 652             /* Since Cp == -1, we have to substract one more */
 653             {
 654               mpn_sub_1(ap, ap, n, MPFR_LIMB_ONE<<sh);
 655               MPFR_ASSERTD(MPFR_LIMB_MSB(ap[n-1]) != 0);
 656             }
 657           /* Ap >= 10000xxx001 */
 658           /* Final exponent -1 since we have shifted the mantissa */
 659           bx--;
 660           /* Update bcp and bcp1 */
 661           MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
 662           MPFR_ASSERTN(bbcp1 != (mp_limb_t) -1);
 663           bcp  = bbcp;
 664           bcp1 = bbcp1;
 665           /* We dont't have anymore a valid Cp+1!
 666              But since Ap >= 100000xxx001, the final sub can't unnormalize!*/
 667         }
 668       MPFR_ASSERTD( !(ap[0] & ~mask) );
 669
 670       /* Rounding */
 671       if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 672         {
 673           if (MPFR_LIKELY(bcp==0))
 674             goto truncate;
 675           else if ((bcp1) || ((ap[0] & (MPFR_LIMB_ONE<<sh)) != 0))
 676             goto sub_one_ulp;
 677           else
 678             goto truncate;
 679         }
 680
 681       /* Update rounding mode */
 682       MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
 683       if (rnd_mode == MPFR_RNDZ && (MPFR_LIKELY(bcp || bcp1)))
 684         goto sub_one_ulp;
 685       goto truncate;
 686     }
 687   MPFR_RET_NEVER_GO_HERE ();
 688
 689   /* Sub one ulp to the result */
 690  sub_one_ulp:
 691   mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
 692   /* Result should be smaller than exact value: inexact=-1 */
 693   inexact = -1;
 694   /* Check normalisation */
 695   if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
 696     {
 697       /* ap was a power of 2, and we lose a bit */
 698       /* Now it is 0111111111111111111[00000 */
 699       mpn_lshift(ap, ap, n, 1);
 700       bx--;
 701       /* And the lost bit x depends on Cp+1, and Cp */
 702       /* Compute Cp+1 if it isn't already compute (ie d==1) */
 703       /* FIXME: Is this case possible? */
 704       if (MPFR_UNLIKELY(d == 1))
 705         bbcp = 0;
 706       DEBUG( printf("(SubOneUlp)Cp=%d, Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0));
 707       /* Compute the last bit (Since we have shifted the mantissa)
 708          we need one more bit!*/
 709       MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
 710       if ( (rnd_mode == MPFR_RNDZ && bcp==0)
 711            || (rnd_mode==MPFR_RNDN && bbcp==0)
 712            || (bcp && bcp1==0) ) /*Exact result*/
 713         {
 714           ap[0] |= MPFR_LIMB_ONE<<sh;
 715           if (rnd_mode == MPFR_RNDN)
 716             inexact = 1;
 717           DEBUG( printf("(SubOneUlp) Last bit set\n") );
 718         }
 719       /* Result could be exact if C'p+1 = 0 and rnd == Zero
 720          since we have had one more bit to the result */
 721       /* Fixme: rnd_mode == MPFR_RNDZ needed ? */
 722       if (bcp1==0 && rnd_mode==MPFR_RNDZ)
 723         {
 724           DEBUG( printf("(SubOneUlp) Exact result\n") );
 725           inexact = 0;
 726         }
 727     }
 728
 729   goto end_of_sub;
 730
 731  truncate:
 732   /* Check if the result is an exact power of 2: 100000000000
 733      in which cases, we could have to do sub_one_ulp due to some nasty reasons:
 734      If Result is a Power of 2:
 735       + If rnd = AWAY,
 736       |  If Cp=-1 and C'p+1 = 0, SubOneUlp and the result is EXACT.
 737          If Cp=-1 and C'p+1 =-1, SubOneUlp and the result is above.
 738          Otherwise truncate
 739       + If rnd = NEAREST,
 740          If Cp= 0 and Cp+1  =-1 and C'p+2=-1, SubOneUlp and the result is above
 741          If cp=-1 and C'p+1 = 0, SubOneUlp and the result is exact.
 742          Otherwise truncate.
 743       X bit should always be set if SubOneUlp*/
 744   if (MPFR_UNLIKELY(ap[n-1] == MPFR_LIMB_HIGHBIT))
 745     {
 746       mp_size_t k = n-1;
 747       do {
 748         k--;
 749       } while (k>=0 && ap[k]==0);
 750       if (MPFR_UNLIKELY(k<0))
 751         {
 752           /* It is a power of 2! */
 753           /* Compute Cp+1 if it isn't already compute (ie d==1) */
 754           /* FIXME: Is this case possible? */
 755           if (d == 1)
 756             bbcp=0;
 757           DEBUG( printf("(Truncate) Cp=%d, Cp+1=%d C'p+1=%d C'p+2=%d\n", \
 758                  bcp!=0, bbcp!=0, bcp1!=0, bbcp1!=0) );
 759           MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
 760           MPFR_ASSERTN((rnd_mode != MPFR_RNDN) || (bcp != 0) || (bbcp == 0) || (bbcp1 != (mp_limb_t) -1));
 761           if (((rnd_mode != MPFR_RNDZ) && bcp)
 762               ||
 763               ((rnd_mode == MPFR_RNDN) && (bcp == 0) && (bbcp) && (bbcp1)))
 764             {
 765               DEBUG( printf("(Truncate) Do sub\n") );
 766               mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
 767               mpn_lshift(ap, ap, n, 1);
 768               ap[0] |= MPFR_LIMB_ONE<<sh;
 769               bx--;
 770               /* FIXME: Explain why it works (or why not)... */
 771               inexact = (bcp1 == 0) ? 0 : (rnd_mode==MPFR_RNDN) ? -1 : 1;
 772               goto end_of_sub;
 773             }
 774         }
 775     }
 776
 777   /* Calcul of Inexact flag.*/
 778   inexact = MPFR_LIKELY(bcp || bcp1) ? 1 : 0;
 779
 780  end_of_sub:
 781   /* Update Expo */
 782   /* FIXME: Is this test really useful?
 783       If d==0      : Exact case. This is never called.
 784       if 1 < d < p : bx=MPFR_EXP(b) or MPFR_EXP(b)-1 > MPFR_EXP(c) > emin
 785       if d == 1    : bx=MPFR_EXP(b). If we could lose any bits, the exact
 786                      normalisation is called.
 787       if d >=  p   : bx=MPFR_EXP(b) >= MPFR_EXP(c) + p > emin
 788      After SubOneUlp, we could have one bit less.
 789       if 1 < d < p : bx >= MPFR_EXP(b)-2 >= MPFR_EXP(c) > emin
 790       if d == 1    : bx >= MPFR_EXP(b)-1 = MPFR_EXP(c) > emin.
 791       if d >=  p   : bx >= MPFR_EXP(b)-1 > emin since p>=2.
 792   */
 793   MPFR_ASSERTD( bx >= __gmpfr_emin);
 794   /*
 795     if (MPFR_UNLIKELY(bx < __gmpfr_emin))
 796     {
 797       DEBUG( printf("(Final Underflow)\n") );
 798       if (rnd_mode == MPFR_RNDN &&
 799           (bx < __gmpfr_emin - 1 ||
 800            (inexact >= 0 && mpfr_powerof2_raw (a))))
 801         rnd_mode = MPFR_RNDZ;
 802       MPFR_TMP_FREE(marker);
 803       return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
 804     }
 805   */
 806   MPFR_SET_EXP (a, bx);
 807
 808   MPFR_TMP_FREE(marker);
 809   MPFR_RET (inexact * MPFR_INT_SIGN (a));
 810 }