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