1/* mpfr_sub1sp -- internal function to perform a "real" substraction 2 All the op must have the same precision 3 4Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 5Contributed by the AriC and Caramel projects, INRIA. 6 7This file is part of the GNU MPFR Library. 8 9The GNU MPFR Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MPFR Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2251 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 31int mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode); 32int 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 134int 135mpfr_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 = MPFR_PREC2LIMBS (p); 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