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