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