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 } 811