1/* mpfr_sub1 -- internal function to perform a "real" subtraction 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-impl.h" 24 25/* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c) 26 Returns 0 iff result is exact, 27 a negative value when the result is less than the exact value, 28 a positive value otherwise. 29*/ 30 31int 32mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 33{ 34 int sign; 35 mpfr_uexp_t diff_exp; 36 mpfr_prec_t cancel, cancel1; 37 mp_size_t cancel2, an, bn, cn, cn0; 38 mp_limb_t *ap, *bp, *cp; 39 mp_limb_t carry, bb, cc; 40 int inexact, shift_b, shift_c, add_exp = 0; 41 int cmp_low = 0; /* used for rounding to nearest: 0 if low(b) = low(c), 42 negative if low(b) < low(c), positive if low(b)>low(c) */ 43 int sh, k; 44 MPFR_TMP_DECL(marker); 45 46 MPFR_TMP_MARK(marker); 47 ap = MPFR_MANT(a); 48 an = MPFR_LIMB_SIZE(a); 49 50 sign = mpfr_cmp2 (b, c, &cancel); 51 if (MPFR_UNLIKELY(sign == 0)) 52 { 53 if (rnd_mode == MPFR_RNDD) 54 MPFR_SET_NEG (a); 55 else 56 MPFR_SET_POS (a); 57 MPFR_SET_ZERO (a); 58 MPFR_RET (0); 59 } 60 61 /* 62 * If subtraction: sign(a) = sign * sign(b) 63 * If addition: sign(a) = sign of the larger argument in absolute value. 64 * 65 * Both cases can be simplidied in: 66 * if (sign>0) 67 * if addition: sign(a) = sign * sign(b) = sign(b) 68 * if subtraction, b is greater, so sign(a) = sign(b) 69 * else 70 * if subtraction, sign(a) = - sign(b) 71 * if addition, sign(a) = sign(c) (since c is greater) 72 * But if it is an addition, sign(b) and sign(c) are opposed! 73 * So sign(a) = - sign(b) 74 */ 75 76 if (sign < 0) /* swap b and c so that |b| > |c| */ 77 { 78 mpfr_srcptr t; 79 MPFR_SET_OPPOSITE_SIGN (a,b); 80 t = b; b = c; c = t; 81 } 82 else 83 MPFR_SET_SAME_SIGN (a,b); 84 85 /* Check if c is too small. 86 A more precise test is to replace 2 by 87 (rnd == MPFR_RNDN) + mpfr_power2_raw (b) 88 but it is more expensive and not very useful */ 89 if (MPFR_UNLIKELY (MPFR_GET_EXP (c) <= MPFR_GET_EXP (b) 90 - (mpfr_exp_t) MAX (MPFR_PREC (a), MPFR_PREC (b)) - 2)) 91 { 92 /* Remember, we can't have an exact result! */ 93 /* A.AAAAAAAAAAAAAAAAA 94 = B.BBBBBBBBBBBBBBB 95 - C.CCCCCCCCCCCCC */ 96 /* A = S*ABS(B) +/- ulp(a) */ 97 MPFR_SET_EXP (a, MPFR_GET_EXP (b)); 98 MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), MPFR_PREC (b), 99 rnd_mode, MPFR_SIGN (a), 100 if (MPFR_UNLIKELY ( ++MPFR_EXP (a) > __gmpfr_emax)) 101 inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a))); 102 /* inexact = mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (a)); */ 103 if (inexact == 0) 104 { 105 /* a = b (Exact) 106 But we know it isn't (Since we have to remove `c') 107 So if we round to Zero, we have to remove one ulp. 108 Otherwise the result is correctly rounded. */ 109 if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a))) 110 { 111 mpfr_nexttozero (a); 112 MPFR_RET (- MPFR_INT_SIGN (a)); 113 } 114 MPFR_RET (MPFR_INT_SIGN (a)); 115 } 116 else 117 { 118 /* A.AAAAAAAAAAAAAA 119 = B.BBBBBBBBBBBBBBB 120 - C.CCCCCCCCCCCCC */ 121 /* It isn't exact so Prec(b) > Prec(a) and the last 122 Prec(b)-Prec(a) bits of `b' are not zeros. 123 Which means that removing c from b can't generate a carry 124 execpt in case of even rounding. 125 In all other case the result and the inexact flag should be 126 correct (We can't have an exact result). 127 In case of EVEN rounding: 128 1.BBBBBBBBBBBBBx10 129 - 1.CCCCCCCCCCCC 130 = 1.BBBBBBBBBBBBBx01 Rounded to Prec(b) 131 = 1.BBBBBBBBBBBBBx Nearest / Rounded to Prec(a) 132 Set gives: 133 1.BBBBBBBBBBBBB0 if inexact == EVEN_INEX (x == 0) 134 1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1) 135 which means we get a wrong rounded result if x==1, 136 i.e. inexact= MPFR_EVEN_INEX */ 137 if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX*MPFR_INT_SIGN (a))) 138 { 139 mpfr_nexttozero (a); 140 inexact = -MPFR_INT_SIGN (a); 141 } 142 MPFR_RET (inexact); 143 } 144 } 145 146 diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c); 147 148 /* reserve a space to store b aligned with the result, i.e. shifted by 149 (-cancel) % GMP_NUMB_BITS to the right */ 150 bn = MPFR_LIMB_SIZE (b); 151 MPFR_UNSIGNED_MINUS_MODULO (shift_b, cancel); 152 cancel1 = (cancel + shift_b) / GMP_NUMB_BITS; 153 154 /* the high cancel1 limbs from b should not be taken into account */ 155 if (MPFR_UNLIKELY (shift_b == 0)) 156 { 157 bp = MPFR_MANT(b); /* no need of an extra space */ 158 /* Ensure ap != bp */ 159 if (MPFR_UNLIKELY (ap == bp)) 160 { 161 bp = MPFR_TMP_LIMBS_ALLOC (bn); 162 MPN_COPY (bp, ap, bn); 163 } 164 } 165 else 166 { 167 bp = MPFR_TMP_LIMBS_ALLOC (bn + 1); 168 bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b); 169 } 170 171 /* reserve a space to store c aligned with the result, i.e. shifted by 172 (diff_exp-cancel) % GMP_NUMB_BITS to the right */ 173 cn = MPFR_LIMB_SIZE(c); 174 if ((UINT_MAX % GMP_NUMB_BITS) == (GMP_NUMB_BITS-1) 175 && ((-(unsigned) 1)%GMP_NUMB_BITS > 0)) 176 shift_c = ((mpfr_uexp_t) diff_exp - cancel) % GMP_NUMB_BITS; 177 else 178 { 179 shift_c = diff_exp - (cancel % GMP_NUMB_BITS); 180 shift_c = (shift_c + GMP_NUMB_BITS) % GMP_NUMB_BITS; 181 } 182 MPFR_ASSERTD( shift_c >= 0 && shift_c < GMP_NUMB_BITS); 183 184 if (MPFR_UNLIKELY(shift_c == 0)) 185 { 186 cp = MPFR_MANT(c); 187 /* Ensure ap != cp */ 188 if (ap == cp) 189 { 190 cp = MPFR_TMP_LIMBS_ALLOC (cn); 191 MPN_COPY(cp, ap, cn); 192 } 193 } 194 else 195 { 196 cp = MPFR_TMP_LIMBS_ALLOC (cn + 1); 197 cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c); 198 } 199 200#ifdef DEBUG 201 printf ("rnd=%s shift_b=%d shift_c=%d diffexp=%lu\n", 202 mpfr_print_rnd_mode (rnd_mode), shift_b, shift_c, 203 (unsigned long) diff_exp); 204#endif 205 206 MPFR_ASSERTD (ap != cp); 207 MPFR_ASSERTD (bp != cp); 208 209 /* here we have shift_c = (diff_exp - cancel) % GMP_NUMB_BITS, 210 0 <= shift_c < GMP_NUMB_BITS 211 thus we want cancel2 = ceil((cancel - diff_exp) / GMP_NUMB_BITS) */ 212 213 /* Possible optimization with a C99 compiler (i.e. well-defined 214 integer division): if MPFR_PREC_MAX is reduced to 215 ((mpfr_prec_t)((mpfr_uprec_t)(~(mpfr_uprec_t)0)>>1) - GMP_NUMB_BITS + 1) 216 and diff_exp is of type mpfr_exp_t (no need for mpfr_uexp_t, since 217 the sum or difference of 2 exponents must be representable, as used 218 by the multiplication code), then the computation of cancel2 could 219 be simplified to 220 cancel2 = (cancel - (diff_exp - shift_c)) / GMP_NUMB_BITS; 221 because cancel, diff_exp and shift_c are all non-negative and 222 these variables are signed. */ 223 224 MPFR_ASSERTD (cancel >= 0); 225 if (cancel >= diff_exp) 226 /* Note that cancel is signed and will be converted to mpfr_uexp_t 227 (type of diff_exp) in the expression below, so that this will 228 work even if cancel is very large and diff_exp = 0. */ 229 cancel2 = (cancel - diff_exp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS; 230 else 231 cancel2 = - (mp_size_t) ((diff_exp - cancel) / GMP_NUMB_BITS); 232 /* the high cancel2 limbs from b should not be taken into account */ 233#ifdef DEBUG 234 printf ("cancel=%lu cancel1=%lu cancel2=%ld\n", 235 (unsigned long) cancel, (unsigned long) cancel1, (long) cancel2); 236#endif 237 238 /* ap[an-1] ap[0] 239 <----------------+-----------|----> 240 <----------PREC(a)----------><-sh-> 241 cancel1 242 limbs bp[bn-cancel1-1] 243 <--...-----><----------------+-----------+-----------> 244 cancel2 245 limbs cp[cn-cancel2-1] cancel2 >= 0 246 <--...--><----------------+----------------+----------------> 247 (-cancel2) cancel2 < 0 248 limbs <----------------+----------------> 249 */ 250 251 /* first part: put in ap[0..an-1] the value of high(b) - high(c), 252 where high(b) consists of the high an+cancel1 limbs of b, 253 and high(c) consists of the high an+cancel2 limbs of c. 254 */ 255 256 /* copy high(b) into a */ 257 if (MPFR_LIKELY(an + (mp_size_t) cancel1 <= bn)) 258 /* a: <----------------+-----------|----> 259 b: <-----------------------------------------> */ 260 MPN_COPY (ap, bp + bn - (an + cancel1), an); 261 else 262 /* a: <----------------+-----------|----> 263 b: <-------------------------> */ 264 if ((mp_size_t) cancel1 < bn) /* otherwise b does not overlap with a */ 265 { 266 MPN_ZERO (ap, an + cancel1 - bn); 267 MPN_COPY (ap + (an + cancel1 - bn), bp, bn - cancel1); 268 } 269 else 270 MPN_ZERO (ap, an); 271 272#ifdef DEBUG 273 printf("after copying high(b), a="); mpfr_print_binary(a); putchar('\n'); 274#endif 275 276 /* subtract high(c) */ 277 if (MPFR_LIKELY(an + cancel2 > 0)) /* otherwise c does not overlap with a */ 278 { 279 mp_limb_t *ap2; 280 281 if (cancel2 >= 0) 282 { 283 if (an + cancel2 <= cn) 284 /* a: <-----------------------------> 285 c: <-----------------------------------------> */ 286 mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an); 287 else 288 /* a: <----------------------------> 289 c: <-------------------------> */ 290 { 291 ap2 = ap + an + (cancel2 - cn); 292 if (cn > cancel2) 293 mpn_sub_n (ap2, ap2, cp, cn - cancel2); 294 } 295 } 296 else /* cancel2 < 0 */ 297 { 298 mp_limb_t borrow; 299 300 if (an + cancel2 <= cn) 301 /* a: <-----------------------------> 302 c: <-----------------------------> */ 303 borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2), 304 an + cancel2); 305 else 306 /* a: <----------------------------> 307 c: <----------------> */ 308 { 309 ap2 = ap + an + cancel2 - cn; 310 borrow = mpn_sub_n (ap2, ap2, cp, cn); 311 } 312 ap2 = ap + an + cancel2; 313 mpn_sub_1 (ap2, ap2, -cancel2, borrow); 314 } 315 } 316 317#ifdef DEBUG 318 printf("after subtracting high(c), a="); 319 mpfr_print_binary(a); 320 putchar('\n'); 321#endif 322 323 /* now perform rounding */ 324 sh = (mpfr_prec_t) an * GMP_NUMB_BITS - MPFR_PREC(a); 325 /* last unused bits from a */ 326 carry = ap[0] & MPFR_LIMB_MASK (sh); 327 ap[0] -= carry; 328 329 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN)) 330 { 331 if (MPFR_LIKELY(sh)) 332 { 333 /* can decide except when carry = 2^(sh-1) [middle] 334 or carry = 0 [truncate, but cannot decide inexact flag] */ 335 if (carry > (MPFR_LIMB_ONE << (sh - 1))) 336 goto add_one_ulp; 337 else if ((0 < carry) && (carry < (MPFR_LIMB_ONE << (sh - 1)))) 338 { 339 inexact = -1; /* result if smaller than exact value */ 340 goto truncate; 341 } 342 /* now carry = 2^(sh-1), in which case cmp_low=2, 343 or carry = 0, in which case cmp_low=0 */ 344 cmp_low = (carry == 0) ? 0 : 2; 345 } 346 } 347 else /* directed rounding: set rnd_mode to RNDZ iff toward zero */ 348 { 349 if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_NEG(a))) 350 rnd_mode = MPFR_RNDZ; 351 352 if (carry) 353 { 354 if (rnd_mode == MPFR_RNDZ) 355 { 356 inexact = -1; 357 goto truncate; 358 } 359 else /* round away */ 360 goto add_one_ulp; 361 } 362 } 363 364 /* we have to consider the low (bn - (an+cancel1)) limbs from b, 365 and the (cn - (an+cancel2)) limbs from c. */ 366 bn -= an + cancel1; 367 cn0 = cn; 368 cn -= an + cancel2; 369 370#ifdef DEBUG 371 printf ("last sh=%d bits from a are %lu, bn=%ld, cn=%ld\n", 372 sh, (unsigned long) carry, (long) bn, (long) cn); 373#endif 374 375 /* for rounding to nearest, we couldn't conclude up to here in the following 376 cases: 377 1. sh = 0, then cmp_low=0: we can either truncate, subtract one ulp 378 or add one ulp: -1 ulp < low(b)-low(c) < 1 ulp 379 2. sh > 0 but the low sh bits from high(b)-high(c) equal 2^(sh-1): 380 -0.5 ulp <= -1/2^sh < low(b)-low(c)-0.5 < 1/2^sh <= 0.5 ulp 381 we can't decide the rounding, in that case cmp_low=2: 382 either we truncate and flag=-1, or we add one ulp and flag=1 383 3. the low sh>0 bits from high(b)-high(c) equal 0: we know we have to 384 truncate but we can't decide the ternary value, here cmp_low=0: 385 -0.5 ulp <= -1/2^sh < low(b)-low(c) < 1/2^sh <= 0.5 ulp 386 we always truncate and inexact can be any of -1,0,1 387 */ 388 389 /* note: here cn might exceed cn0, in which case we consider a zero limb */ 390 for (k = 0; (bn > 0) || (cn > 0); k = 1) 391 { 392 /* if cmp_low < 0, we know low(b) - low(c) < 0 393 if cmp_low > 0, we know low(b) - low(c) > 0 394 (more precisely if cmp_low = 2, low(b) - low(c) = 0.5 ulp so far) 395 if cmp_low = 0, so far low(b) - low(c) = 0 */ 396 397 /* get next limbs */ 398 bb = (bn > 0) ? bp[--bn] : 0; 399 if ((cn > 0) && (cn-- <= cn0)) 400 cc = cp[cn]; 401 else 402 cc = 0; 403 404 /* cmp_low compares low(b) and low(c) */ 405 if (cmp_low == 0) /* case 1 or 3 */ 406 cmp_low = (bb < cc) ? -2+k : (bb > cc) ? 1 : 0; 407 408 /* Case 1 for k=0 splits into 7 subcases: 409 1a: bb > cc + half 410 1b: bb = cc + half 411 1c: 0 < bb - cc < half 412 1d: bb = cc 413 1e: -half < bb - cc < 0 414 1f: bb - cc = -half 415 1g: bb - cc < -half 416 417 Case 2 splits into 3 subcases: 418 2a: bb > cc 419 2b: bb = cc 420 2c: bb < cc 421 422 Case 3 splits into 3 subcases: 423 3a: bb > cc 424 3b: bb = cc 425 3c: bb < cc 426 */ 427 428 /* the case rounding to nearest with sh=0 is special since one couldn't 429 subtract above 1/2 ulp in the trailing limb of the result */ 430 if (rnd_mode == MPFR_RNDN && sh == 0 && k == 0) /* case 1 for k=0 */ 431 { 432 mp_limb_t half = MPFR_LIMB_HIGHBIT; 433 434 /* add one ulp if bb > cc + half 435 truncate if cc - half < bb < cc + half 436 sub one ulp if bb < cc - half 437 */ 438 439 if (cmp_low < 0) /* bb < cc: -1 ulp < low(b) - low(c) < 0, 440 cases 1e, 1f and 1g */ 441 { 442 if (cc >= half) 443 cc -= half; 444 else /* since bb < cc < half, bb+half < 2*half */ 445 bb += half; 446 /* now we have bb < cc + half: 447 we have to subtract one ulp if bb < cc, 448 and truncate if bb > cc */ 449 } 450 else if (cmp_low >= 0) /* bb >= cc, cases 1a to 1d */ 451 { 452 if (cc < half) 453 cc += half; 454 else /* since bb >= cc >= half, bb - half >= 0 */ 455 bb -= half; 456 /* now we have bb > cc - half: we have to add one ulp if bb > cc, 457 and truncate if bb < cc */ 458 if (cmp_low > 0) 459 cmp_low = 2; 460 } 461 } 462 463#ifdef DEBUG 464 printf ("k=%u bb=%lu cc=%lu cmp_low=%d\n", k, 465 (unsigned long) bb, (unsigned long) cc, cmp_low); 466#endif 467 if (cmp_low < 0) /* low(b) - low(c) < 0: either truncate or subtract 468 one ulp */ 469 { 470 if (rnd_mode == MPFR_RNDZ) 471 goto sub_one_ulp; /* set inexact=-1 */ 472 else if (rnd_mode != MPFR_RNDN) /* round away */ 473 { 474 inexact = 1; 475 goto truncate; 476 } 477 else /* round to nearest */ 478 { 479 /* If cmp_low < 0 and bb > cc, then -0.5 ulp < low(b)-low(c) < 0, 480 whatever the value of sh. 481 If sh>0, then cmp_low < 0 implies that the initial neglected 482 sh bits were 0 (otherwise cmp_low=2 initially), thus the 483 weight of the new bits is less than 0.5 ulp too. 484 If k > 0 (and sh=0) this means that either the first neglected 485 limbs bb and cc were equal (thus cmp_low was 0 for k=0), 486 or we had bb - cc = -0.5 ulp or 0.5 ulp. 487 The last case is not possible here since we would have 488 cmp_low > 0 which is sticky. 489 In the first case (where we have cmp_low = -1), we truncate, 490 whereas in the 2nd case we have cmp_low = -2 and we subtract 491 one ulp. 492 */ 493 if (bb > cc || sh > 0 || cmp_low == -1) 494 { /* -0.5 ulp < low(b)-low(c) < 0, 495 bb > cc corresponds to cases 1e and 1f1 496 sh > 0 corresponds to cases 3c and 3b3 497 cmp_low = -1 corresponds to case 1d3 (also 3b3) */ 498 inexact = 1; 499 goto truncate; 500 } 501 else if (bb < cc) /* here sh = 0 and low(b)-low(c) < -0.5 ulp, 502 this corresponds to cases 1g and 1f3 */ 503 goto sub_one_ulp; 504 /* the only case where we can't conclude is sh=0 and bb=cc, 505 i.e., we have low(b) - low(c) = -0.5 ulp (up to now), thus 506 we don't know if we must truncate or subtract one ulp. 507 Note: for sh=0 we can't have low(b) - low(c) = -0.5 ulp up to 508 now, since low(b) - low(c) > 1/2^sh */ 509 } 510 } 511 else if (cmp_low > 0) /* 0 < low(b) - low(c): either truncate or 512 add one ulp */ 513 { 514 if (rnd_mode == MPFR_RNDZ) 515 { 516 inexact = -1; 517 goto truncate; 518 } 519 else if (rnd_mode != MPFR_RNDN) /* round away */ 520 goto add_one_ulp; 521 else /* round to nearest */ 522 { 523 if (bb > cc) 524 { 525 /* if sh=0, then bb>cc means that low(b)-low(c) > 0.5 ulp, 526 and similarly when cmp_low=2 */ 527 if (cmp_low == 2) /* cases 1a, 1b1, 2a and 2b1 */ 528 goto add_one_ulp; 529 /* sh > 0 and cmp_low > 0: this implies that the sh initial 530 neglected bits were 0, and the remaining low(b)-low(c)>0, 531 but its weight is less than 0.5 ulp */ 532 else /* 0 < low(b) - low(c) < 0.5 ulp, this corresponds to 533 cases 3a, 1d1 and 3b1 */ 534 { 535 inexact = -1; 536 goto truncate; 537 } 538 } 539 else if (bb < cc) /* 0 < low(b) - low(c) < 0.5 ulp, cases 1c, 540 1b3, 2b3 and 2c */ 541 { 542 inexact = -1; 543 goto truncate; 544 } 545 /* the only case where we can't conclude is bb=cc, i.e., 546 low(b) - low(c) = 0.5 ulp (up to now), thus we don't know 547 if we must truncate or add one ulp. */ 548 } 549 } 550 /* after k=0, we cannot conclude in the following cases, we split them 551 according to the values of bb and cc for k=1: 552 1b. sh=0 and cmp_low = 1 and bb-cc = half [around 0.5 ulp] 553 1b1. bb > cc: add one ulp, inex = 1 554 1b2: bb = cc: cannot conclude 555 1b3: bb < cc: truncate, inex = -1 556 1d. sh=0 and cmp_low = 0 and bb-cc = 0 [around 0] 557 1d1: bb > cc: truncate, inex = -1 558 1d2: bb = cc: cannot conclude 559 1d3: bb < cc: truncate, inex = +1 560 1f. sh=0 and cmp_low = -1 and bb-cc = -half [around -0.5 ulp] 561 1f1: bb > cc: truncate, inex = +1 562 1f2: bb = cc: cannot conclude 563 1f3: bb < cc: sub one ulp, inex = -1 564 2b. sh > 0 and cmp_low = 2 and bb=cc [around 0.5 ulp] 565 2b1. bb > cc: add one ulp, inex = 1 566 2b2: bb = cc: cannot conclude 567 2b3: bb < cc: truncate, inex = -1 568 3b. sh > 0 and cmp_low = 0 [around 0] 569 3b1. bb > cc: truncate, inex = -1 570 3b2: bb = cc: cannot conclude 571 3b3: bb < cc: truncate, inex = +1 572 */ 573 } 574 575 if ((rnd_mode == MPFR_RNDN) && cmp_low != 0) 576 { 577 /* even rounding rule */ 578 if ((ap[0] >> sh) & 1) 579 { 580 if (cmp_low < 0) 581 goto sub_one_ulp; 582 else 583 goto add_one_ulp; 584 } 585 else 586 inexact = (cmp_low > 0) ? -1 : 1; 587 } 588 else 589 inexact = 0; 590 goto truncate; 591 592 sub_one_ulp: /* sub one unit in last place to a */ 593 mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh); 594 inexact = -1; 595 goto end_of_sub; 596 597 add_one_ulp: /* add one unit in last place to a */ 598 if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh))) 599 /* result is a power of 2: 11111111111111 + 1 = 1000000000000000 */ 600 { 601 ap[an-1] = MPFR_LIMB_HIGHBIT; 602 add_exp = 1; 603 } 604 inexact = 1; /* result larger than exact value */ 605 606 truncate: 607 if (MPFR_UNLIKELY((ap[an-1] >> (GMP_NUMB_BITS - 1)) == 0)) 608 /* case 1 - epsilon */ 609 { 610 ap[an-1] = MPFR_LIMB_HIGHBIT; 611 add_exp = 1; 612 } 613 614 end_of_sub: 615 /* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking 616 care of underflows/overflows in that computation, and of the allowed 617 exponent range */ 618 if (MPFR_LIKELY(cancel)) 619 { 620 mpfr_exp_t exp_a; 621 622 cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */ 623 exp_a = MPFR_GET_EXP (b) - cancel; 624 if (MPFR_UNLIKELY(exp_a < __gmpfr_emin)) 625 { 626 MPFR_TMP_FREE(marker); 627 if (rnd_mode == MPFR_RNDN && 628 (exp_a < __gmpfr_emin - 1 || 629 (inexact >= 0 && mpfr_powerof2_raw (a)))) 630 rnd_mode = MPFR_RNDZ; 631 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); 632 } 633 MPFR_SET_EXP (a, exp_a); 634 } 635 else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */ 636 { 637 /* in case cancel = 0, add_exp can still be 1, in case b is just 638 below a power of two, c is very small, prec(a) < prec(b), 639 and rnd=away or nearest */ 640 mpfr_exp_t exp_b; 641 642 exp_b = MPFR_GET_EXP (b); 643 if (MPFR_UNLIKELY(add_exp && exp_b == __gmpfr_emax)) 644 { 645 MPFR_TMP_FREE(marker); 646 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 647 } 648 MPFR_SET_EXP (a, exp_b + add_exp); 649 } 650 MPFR_TMP_FREE(marker); 651#ifdef DEBUG 652 printf ("result is a="); mpfr_print_binary(a); putchar('\n'); 653#endif 654 /* check that result is msb-normalized */ 655 MPFR_ASSERTD(ap[an-1] > ~ap[an-1]); 656 MPFR_RET (inexact * MPFR_INT_SIGN (a)); 657} 658