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