1/* mpfr_add1 -- internal function to perform a "real" addition 2 3Copyright 1999-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba 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 20https://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|), assuming that b and c 26 are not NaN, Inf, nor zero. Assumes EXP(b) >= EXP(c). 27*/ 28MPFR_HOT_FUNCTION_ATTR int 29mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 30{ 31 mp_limb_t *ap, *bp, *cp; 32 mpfr_prec_t aq, bq, cq, aq2; 33 mp_size_t an, bn, cn; 34 mpfr_exp_t difw, exp, diff_exp; 35 int sh, rb, fb, inex; 36 MPFR_TMP_DECL(marker); 37 38 MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); 39 MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); 40 MPFR_ASSERTD (! MPFR_UBF_EXP_LESS_P (b, c)); 41 42 if (MPFR_UNLIKELY (MPFR_IS_UBF (b))) 43 { 44 exp = MPFR_UBF_GET_EXP (b); 45 if (exp > __gmpfr_emax) 46 return mpfr_overflow (a, rnd_mode, MPFR_SIGN (b));; 47 } 48 else 49 exp = MPFR_GET_EXP (b); 50 51 MPFR_ASSERTD (exp <= __gmpfr_emax); 52 53 MPFR_TMP_MARK(marker); 54 55 aq = MPFR_GET_PREC (a); 56 bq = MPFR_GET_PREC (b); 57 cq = MPFR_GET_PREC (c); 58 59 an = MPFR_PREC2LIMBS (aq); /* number of limbs of a */ 60 aq2 = (mpfr_prec_t) an * GMP_NUMB_BITS; 61 sh = aq2 - aq; /* non-significant bits in low limb */ 62 63 bn = MPFR_PREC2LIMBS (bq); /* number of limbs of b */ 64 cn = MPFR_PREC2LIMBS (cq); /* number of limbs of c */ 65 66 ap = MPFR_MANT(a); 67 bp = MPFR_MANT(b); 68 cp = MPFR_MANT(c); 69 70 if (MPFR_UNLIKELY(ap == bp)) 71 { 72 bp = MPFR_TMP_LIMBS_ALLOC (bn); 73 MPN_COPY (bp, ap, bn); 74 if (ap == cp) 75 { cp = bp; } 76 } 77 else if (ap == cp) 78 { 79 cp = MPFR_TMP_LIMBS_ALLOC (cn); 80 MPN_COPY(cp, ap, cn); 81 } 82 83 MPFR_SET_SAME_SIGN(a, b); 84 MPFR_UPDATE2_RND_MODE (rnd_mode, MPFR_SIGN (b)); 85 /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ, MPFR_RNDA or MPFR_RNDF. */ 86 if (MPFR_UNLIKELY (MPFR_IS_UBF (c))) 87 { 88 MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX); 89 diff_exp = mpfr_ubf_diff_exp (b, c); 90 } 91 else 92 diff_exp = exp - MPFR_GET_EXP (c); 93 94 MPFR_ASSERTD (diff_exp >= 0); 95 96 /* 97 * 1. Compute the significant part A', the non-significant bits of A 98 * are taken into account. 99 * 100 * 2. Perform the rounding. At each iteration, we remember: 101 * _ r = rounding bit 102 * _ f = following bits (same value) 103 * where the result has the form: [number A]rfff...fff + a remaining 104 * value in the interval [0,2) ulp. We consider the most significant 105 * bits of the remaining value to update the result; a possible carry 106 * is immediately taken into account and A is updated accordingly. As 107 * soon as the bits f don't have the same value, A can be rounded. 108 * Variables: 109 * _ rb = rounding bit (0 or 1). 110 * _ fb = following bits (0 or 1), then sticky bit. 111 * If fb == 0, the only thing that can change is the sticky bit. 112 */ 113 114 rb = fb = -1; /* means: not initialized */ 115 116 if (MPFR_UNLIKELY (MPFR_UEXP (aq2) <= diff_exp)) 117 { /* c does not overlap with a' */ 118 if (MPFR_UNLIKELY(an > bn)) 119 { /* a has more limbs than b */ 120 /* copy b to the most significant limbs of a */ 121 MPN_COPY(ap + (an - bn), bp, bn); 122 /* zero the least significant limbs of a */ 123 MPN_ZERO(ap, an - bn); 124 } 125 else /* an <= bn */ 126 { 127 /* copy the most significant limbs of b to a */ 128 MPN_COPY(ap, bp + (bn - an), an); 129 } 130 } 131 else /* aq2 > diff_exp */ 132 { /* c overlaps with a' */ 133 mp_limb_t *a2p; 134 mp_limb_t cc; 135 mpfr_prec_t dif; 136 mp_size_t difn, k; 137 int shift; 138 139 /* copy c (shifted) into a */ 140 141 dif = aq2 - diff_exp; 142 /* dif is the number of bits of c which overlap with a' */ 143 144 difn = MPFR_PREC2LIMBS (dif); 145 /* only the highest difn limbs from c have to be considered */ 146 if (MPFR_UNLIKELY(difn > cn)) 147 { 148 /* c doesn't have enough limbs; take into account the virtual 149 zero limbs now by zeroing the least significant limbs of a' */ 150 MPFR_ASSERTD(difn - cn <= an); 151 MPN_ZERO(ap, difn - cn); 152 difn = cn; 153 } 154 k = diff_exp / GMP_NUMB_BITS; 155 156 /* zero the most significant k limbs of a */ 157 a2p = ap + (an - k); 158 MPN_ZERO(a2p, k); 159 160 shift = diff_exp % GMP_NUMB_BITS; 161 162 if (MPFR_LIKELY(shift)) 163 { 164 MPFR_ASSERTD(a2p - difn >= ap); 165 cc = mpn_rshift(a2p - difn, cp + (cn - difn), difn, shift); 166 if (MPFR_UNLIKELY(a2p - difn > ap)) 167 *(a2p - difn - 1) = cc; 168 } 169 else 170 MPN_COPY(a2p - difn, cp + (cn - difn), difn); 171 172 /* add b to a */ 173 cc = an > bn 174 ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn) 175 : mpn_add_n(ap, ap, bp + (bn - an), an); 176 177 if (MPFR_UNLIKELY(cc)) /* carry */ 178 { 179 if (MPFR_UNLIKELY(exp == __gmpfr_emax)) 180 { 181 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 182 goto end_of_add; 183 } 184 exp++; 185 rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */ 186 if (MPFR_LIKELY(sh)) 187 { 188 mp_limb_t mask, bb; 189 190 mask = MPFR_LIMB_MASK (sh); 191 bb = ap[0] & mask; 192 ap[0] &= MPFR_LIMB_LSHIFT (~mask, 1); 193 if (bb == 0) 194 fb = 0; 195 else if (bb == mask) 196 fb = 1; 197 } 198 mpn_rshift(ap, ap, an, 1); 199 ap[an-1] += MPFR_LIMB_HIGHBIT; 200 if (sh && fb < 0) 201 goto rounding; 202 } /* cc */ 203 } /* aq2 > diff_exp */ 204 205 /* zero the non-significant bits of a */ 206 if (MPFR_LIKELY(rb < 0 && sh)) 207 { 208 mp_limb_t mask, bb; 209 210 mask = MPFR_LIMB_MASK (sh); 211 bb = ap[0] & mask; 212 ap[0] &= ~mask; 213 rb = bb >> (sh - 1); 214 if (MPFR_LIKELY(sh > 1)) 215 { 216 mask >>= 1; 217 bb &= mask; 218 if (bb == 0) 219 fb = 0; 220 else if (bb == mask) 221 fb = 1; 222 else 223 goto rounding; 224 } 225 } 226 227 /* Determine rounding and sticky bits (and possible carry). 228 In faithful rounding, we may stop two bits after ulp(a): 229 the approximation is regarded as the number formed by a, 230 the rounding bit rb and an additional bit fb; and the 231 corresponding error is < 1/2 ulp of the unrounded result. */ 232 233 difw = (mpfr_exp_t) an - (mpfr_exp_t) (diff_exp / GMP_NUMB_BITS); 234 /* difw is the number of limbs from b (regarded as having an infinite 235 precision) that have already been combined with c; -n if the next 236 n limbs from b won't be combined with c. */ 237 238 if (MPFR_UNLIKELY(bn > an)) 239 { /* there are still limbs from b that haven't been taken into account */ 240 mp_size_t bk; 241 242 if (fb == 0 && difw <= 0) 243 { 244 fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */ 245 goto rounding; 246 } 247 248 bk = bn - an; /* index of lowest considered limb from b, > 0 */ 249 while (difw < 0) 250 { /* ulp(next limb from b) > msb(c) */ 251 mp_limb_t bb; 252 253 bb = bp[--bk]; 254 255 MPFR_ASSERTD(fb != 0); 256 if (fb > 0) 257 { 258 /* Note: Here, we can round to nearest, but the loop may still 259 be necessary to determine whether there is a carry from c, 260 which will have an effect on the ternary value. However, in 261 faithful rounding, we do not have to determine the ternary 262 value, so that we can end the loop here. */ 263 if (bb != MPFR_LIMB_MAX || rnd_mode == MPFR_RNDF) 264 goto rounding; 265 } 266 else /* fb not initialized yet */ 267 { 268 if (rb < 0) /* rb not initialized yet */ 269 { 270 rb = bb >> (GMP_NUMB_BITS - 1); 271 bb |= MPFR_LIMB_HIGHBIT; 272 } 273 fb = 1; 274 if (bb != MPFR_LIMB_MAX) 275 goto rounding; 276 } 277 278 if (bk == 0) 279 { /* b has entirely been read */ 280 fb = 1; /* c hasn't been taken into account 281 ==> sticky bit != 0 */ 282 goto rounding; 283 } 284 285 difw++; 286 } /* while */ 287 MPFR_ASSERTD(bk > 0 && difw >= 0); 288 289 if (difw <= cn) 290 { 291 mp_size_t ck; 292 mp_limb_t cprev; 293 int difs; 294 295 ck = cn - difw; 296 difs = diff_exp % GMP_NUMB_BITS; 297 298 if (difs == 0 && ck == 0) 299 goto c_read; 300 301 cprev = ck == cn ? 0 : cp[ck]; 302 303 if (fb < 0) 304 { 305 mp_limb_t bb, cc; 306 307 if (difs) 308 { 309 cc = cprev << (GMP_NUMB_BITS - difs); 310 if (--ck >= 0) 311 { 312 cprev = cp[ck]; 313 cc += cprev >> difs; 314 } 315 } 316 else 317 cc = cp[--ck]; 318 319 bb = bp[--bk] + cc; 320 321 if (bb < cc /* carry */ 322 && (rb < 0 || (rb ^= 1) == 0) 323 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh)) 324 { 325 if (exp == __gmpfr_emax) 326 { 327 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 328 goto end_of_add; 329 } 330 exp++; 331 ap[an-1] = MPFR_LIMB_HIGHBIT; 332 rb = 0; 333 } 334 335 if (rb < 0) /* rb not initialized yet */ 336 { 337 rb = bb >> (GMP_NUMB_BITS - 1); 338 bb <<= 1; 339 bb |= bb >> (GMP_NUMB_BITS - 1); 340 } 341 342 fb = bb != 0; 343 if (fb && bb != MPFR_LIMB_MAX) 344 goto rounding; 345 } /* fb < 0 */ 346 347 /* At least two bits after ulp(a) have been read, which is 348 sufficient for faithful rounding, as we do not need to 349 determine on which side of a breakpoint the result is. */ 350 if (rnd_mode == MPFR_RNDF) 351 goto rounding; 352 353 while (bk > 0) 354 { 355 mp_limb_t bb, cc; 356 357 if (difs) 358 { 359 if (ck < 0) 360 goto c_read; 361 cc = cprev << (GMP_NUMB_BITS - difs); 362 if (--ck >= 0) 363 { 364 cprev = cp[ck]; 365 cc += cprev >> difs; 366 } 367 } 368 else 369 { 370 if (ck == 0) 371 goto c_read; 372 cc = cp[--ck]; 373 } 374 375 bb = bp[--bk] + cc; 376 if (bb < cc) /* carry */ 377 { 378 fb ^= 1; 379 if (fb) 380 goto rounding; 381 rb ^= 1; 382 if (rb == 0 && mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)) 383 { 384 if (MPFR_UNLIKELY(exp == __gmpfr_emax)) 385 { 386 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 387 goto end_of_add; 388 } 389 exp++; 390 ap[an-1] = MPFR_LIMB_HIGHBIT; 391 } 392 } /* bb < cc */ 393 394 if (!fb && bb != 0) 395 { 396 fb = 1; 397 goto rounding; 398 } 399 if (fb && bb != MPFR_LIMB_MAX) 400 goto rounding; 401 } /* while */ 402 403 /* b has entirely been read */ 404 405 if (fb || ck < 0) 406 goto rounding; 407 if (difs && MPFR_LIMB_LSHIFT(cprev, GMP_NUMB_BITS - difs) != 0) 408 { 409 fb = 1; 410 goto rounding; 411 } 412 while (ck) 413 { 414 if (cp[--ck]) 415 { 416 fb = 1; 417 goto rounding; 418 } 419 } /* while */ 420 } /* difw <= cn */ 421 else 422 { /* c has entirely been read */ 423 c_read: 424 if (fb < 0) /* fb not initialized yet */ 425 { 426 mp_limb_t bb; 427 428 MPFR_ASSERTD(bk > 0); 429 bb = bp[--bk]; 430 if (rb < 0) /* rb not initialized yet */ 431 { 432 rb = bb >> (GMP_NUMB_BITS - 1); 433 bb &= ~MPFR_LIMB_HIGHBIT; 434 } 435 fb = bb != 0; 436 } /* fb < 0 */ 437 if (fb || rnd_mode == MPFR_RNDF) 438 goto rounding; 439 while (bk) 440 { 441 if (bp[--bk]) 442 { 443 fb = 1; 444 goto rounding; 445 } 446 } /* while */ 447 } /* difw > cn */ 448 } /* bn > an */ 449 else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */ 450 { /* b has entirely been read */ 451 if (difw > cn) 452 { /* c has entirely been read */ 453 if (rb < 0) 454 rb = 0; 455 fb = 0; 456 } 457 else if (diff_exp > MPFR_UEXP (aq2)) 458 { /* b is followed by at least a zero bit, then by c */ 459 if (rb < 0) 460 rb = 0; 461 fb = 1; 462 } 463 else 464 { 465 mp_size_t ck; 466 int difs; 467 468 MPFR_ASSERTD(difw >= 0 && cn >= difw); 469 ck = cn - difw; 470 difs = diff_exp % GMP_NUMB_BITS; 471 472 if (difs == 0 && ck == 0) 473 { /* c has entirely been read */ 474 if (rb < 0) 475 rb = 0; 476 fb = 0; 477 } 478 else 479 { 480 mp_limb_t cc; 481 482 cc = difs ? (MPFR_ASSERTD(ck < cn), 483 cp[ck] << (GMP_NUMB_BITS - difs)) : cp[--ck]; 484 if (rb < 0) 485 { 486 rb = cc >> (GMP_NUMB_BITS - 1); 487 cc &= ~MPFR_LIMB_HIGHBIT; 488 } 489 if (cc == 0 && rnd_mode == MPFR_RNDF) 490 { 491 fb = 0; 492 goto rounding; 493 } 494 while (cc == 0) 495 { 496 if (ck == 0) 497 { 498 fb = 0; 499 goto rounding; 500 } 501 cc = cp[--ck]; 502 } /* while */ 503 fb = 1; 504 } 505 } 506 } /* fb != 1 */ 507 508 rounding: 509 /* rnd_mode should be one of MPFR_RNDN, MPFR_RNDF, MPFR_RNDZ or MPFR_RNDA */ 510 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN || rnd_mode == MPFR_RNDF)) 511 { 512 if (fb == 0) 513 { 514 if (rb == 0) 515 { 516 inex = 0; 517 goto set_exponent; 518 } 519 /* round to even */ 520 if (ap[0] & (MPFR_LIMB_ONE << sh)) 521 goto rndn_away; 522 else 523 goto rndn_zero; 524 } 525 if (rb == 0) 526 { 527 rndn_zero: 528 inex = MPFR_IS_NEG(a) ? 1 : -1; 529 goto set_exponent; 530 } 531 else 532 { 533 rndn_away: 534 inex = MPFR_IS_POS(a) ? 1 : -1; 535 goto add_one_ulp; 536 } 537 } 538 else if (rnd_mode == MPFR_RNDZ) 539 { 540 inex = rb || fb ? (MPFR_IS_NEG(a) ? 1 : -1) : 0; 541 goto set_exponent; 542 } 543 else 544 { 545 MPFR_ASSERTN (rnd_mode == MPFR_RNDA); 546 inex = rb || fb ? (MPFR_IS_POS(a) ? 1 : -1) : 0; 547 if (inex) 548 goto add_one_ulp; 549 else 550 goto set_exponent; 551 } 552 553 add_one_ulp: /* add one unit in last place to a */ 554 if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh))) 555 { 556 if (MPFR_UNLIKELY(exp == __gmpfr_emax)) 557 { 558 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 559 goto end_of_add; 560 } 561 exp++; 562 ap[an-1] = MPFR_LIMB_HIGHBIT; 563 } 564 565 set_exponent: 566 if (MPFR_UNLIKELY (exp < __gmpfr_emin)) /* possible if b and c are UBF's */ 567 { 568 if (rnd_mode == MPFR_RNDN && 569 (exp < __gmpfr_emin - 1 || 570 (inex >= 0 && mpfr_powerof2_raw (a)))) 571 rnd_mode = MPFR_RNDZ; 572 inex = mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); 573 goto end_of_add; 574 } 575 MPFR_SET_EXP (a, exp); 576 577 end_of_add: 578 MPFR_TMP_FREE(marker); 579 MPFR_RET (inex); 580} 581