1/* mpfr_add1 -- internal function to perform a "real" addition 2 3Copyright 1999, 2000, 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|), assuming b and c have same sign, 26 and are not NaN, Inf, nor zero. */ 27int 28mpfr_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