diff options
Diffstat (limited to 'add.c')
-rw-r--r-- | add.c | 470 |
1 files changed, 289 insertions, 181 deletions
@@ -25,184 +25,263 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -/* #define DEBUG */ - -extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, - mp_rnd_t, int)); -void mpfr_add1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t, int)); - #define ONE ((mp_limb_t) 1) +extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, + mp_rnd_t, mp_exp_unsigned_t)); +void mpfr_add1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, + mp_rnd_t, mp_exp_unsigned_t)); + /* signs of b and c are supposed equal, diff_exp is the difference between the exponents of b and c, which is supposed >= 0 */ -void +void #if __STDC__ -mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, - mp_rnd_t rnd_mode, int diff_exp) +mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, + mp_rnd_t rnd_mode, mp_exp_unsigned_t diff_exp) #else -mpfr_add1 (a, b, c, rnd_mode, diff_exp) +mpfr_add1 (a, b, c, rnd_mode, diff_exp) mpfr_ptr a; mpfr_srcptr b; mpfr_srcptr c; mp_rnd_t rnd_mode; - int diff_exp; + mp_exp_unsigned_t diff_exp; #endif { - mp_limb_t *ap, *bp, *cp, cc, c2=0, c3=0; - unsigned int an, bn, cn; - int sh, dif, k; - TMP_DECL(marker); - - TMP_MARK(marker); + mp_limb_t *ap, *bp, *cp; + mp_limb_t cc, c2 = 0, c3 = 0; /* should be removed later */ + mp_exp_t dif; /* should be removed later */ + mp_prec_t aq, bq, cq; + mp_size_t an, bn, cn, k; + int sh; + int nulp = 1; + TMP_DECL(marker); + + MPFR_ASSERTN(MPFR_NOTZERO(c)); + + TMP_MARK(marker); ap = MPFR_MANT(a); bp = MPFR_MANT(b); cp = MPFR_MANT(c); - if (ap == bp) { - bp = (mp_ptr) TMP_ALLOC(MPFR_ABSSIZE(b) * BYTES_PER_MP_LIMB); + + if (ap == bp) + { + bp = (mp_ptr) TMP_ALLOC(MPFR_ABSSIZE(b) * BYTES_PER_MP_LIMB); MPN_COPY (bp, ap, MPFR_ABSSIZE(b)); if (ap == cp) { cp = bp; } } else if (ap == cp) - { - cp = (mp_ptr) TMP_ALLOC (MPFR_ABSSIZE(c) * BYTES_PER_MP_LIMB); - MPN_COPY(cp, ap, MPFR_ABSSIZE(c)); - } + { + cp = (mp_ptr) TMP_ALLOC (MPFR_ABSSIZE(c) * BYTES_PER_MP_LIMB); + MPN_COPY(cp, ap, MPFR_ABSSIZE(c)); + } - an = (MPFR_PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */ + aq = MPFR_PREC(a); + bq = MPFR_PREC(b); + cq = MPFR_PREC(c); - sh = an*BITS_PER_MP_LIMB-MPFR_PREC(a); /* non-significant bits in low limb */ - bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */ - cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB + 1; - MPFR_EXP(a) = MPFR_EXP(b); + an = (aq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of a */ + sh = an*BITS_PER_MP_LIMB - aq; /* non-significant bits in low limb */ + bn = (bq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */ + cn = (cq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of c */ - if (MPFR_SIGN(a) * MPFR_SIGN(b) < 0) MPFR_CHANGE_SIGN(a); + MPFR_EXP(a) = MPFR_EXP(b); + MPFR_SET_SAME_SIGN(a, b); - /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit - through rounding */ - dif = MPFR_PREC(a) - diff_exp; + if (aq <= diff_exp) + { -#ifdef DEBUG - printf("diff_exp=%u dif=MPFR_PREC(a)-diff_exp=%d\n", diff_exp, dif); - printf("b= "); mpfr_print_raw(b); putchar('\n'); - printf("c="); for (k=0;k<diff_exp;k++) putchar(' '); - if (MPFR_SIGN(c)>0) putchar(' '); mpfr_print_raw(c); putchar('\n'); -#endif - if (dif <= 0) - { /* diff_exp>=MPFR_PREC(a): c does not overlap with a */ - /* either MPFR_PREC(b)<=MPFR_PREC(a), and we can copy the mantissa of b directly - into that of a, or MPFR_PREC(b)>MPFR_PREC(a) and we have to round b+c */ + /* either MPFR_PREC(b)<=MPFR_PREC(a), and we can copy the mantissa of b + directly into that of a, or MPFR_PREC(b)>MPFR_PREC(a) and we have to + round b+c */ - if (MPFR_PREC(b) <= MPFR_PREC(a)) { + if (bq <= aq) + { - MPN_COPY(ap+(an-bn), bp, bn); - /* fill low significant limbs with zero */ + MPN_COPY(ap + (an - bn), bp, bn); - for (bp=ap;bn<an;bn++) *bp++=0; + /* fill low significant limbs with zero */ + MPN_ZERO(ap, an - bn); - /* now take c into account */ - if (rnd_mode==GMP_RNDN) { - + /* now take c into account (c != 0) */ + if (rnd_mode == GMP_RNDN) + { /* to nearest */ /* if diff_exp > MPFR_PREC(a), no change */ - if (diff_exp==MPFR_PREC(a)) { + if (aq == diff_exp) + { + /* as c is normalized, we have to add one to the lsb of a + if c>1/2, or c=1/2 and lsb(a)=1 (round to even) */ + + mp_limb_t cc; + mp_limb_t *cp2 = cp + (cn-1); /* highest limb */ + + cc = *cp2 - MP_LIMB_T_HIGHBIT; + while (cc == 0 && cp2 > cp) cc = *--cp2; + + if (cc || (ap[0] & (ONE<<sh))) goto add_one_ulp; + /* mant(c) != 1/2 or mant(c) = 1/2: add 1 iff lsb(a)=1 */ + } /* aq == diff_exp */ + } /* rnd_mode == GMP_RNDN */ + else if ((rnd_mode == GMP_RNDU && MPFR_ISNONNEG(b)) || + (rnd_mode == GMP_RNDD && MPFR_ISNEG(b))) /* round up */ + goto add_one_ulp; + /* in the other cases (round to zero, or up/down with sign -/+), + nothing to do */ + } /* bq <= aq */ + else /* bq > aq */ + { + mp_limb_t inv, bb, cc = 0; + int difs, r; + mp_exp_t difw; /* mp_exp_t may be larger than mp_size_t */ - /* if c is not zero, then as it is normalized, we have to add - one to the lsb of a if c>1/2, or c=1/2 and lsb(a)=1 (round to - even) */ - - if (MPFR_NOTZERO(c)) { + /* MPFR_PREC(b)>MPFR_PREC(a) : we have to round b+c */ + k = bn - an; - /* c is not zero */ - /* check whether mant(c)=1/2 or not */ + /* first copy the 'an' most significant limbs of b to a */ + MPN_COPY(ap, bp + k, an); - mp_limb_t *cp2 = cp + (cn-1); /* highest limb */ + difs = diff_exp % BITS_PER_MP_LIMB; + difw = an - (mp_exp_t) (diff_exp / BITS_PER_MP_LIMB); + MPFR_ASSERTN(difw <= 1); - cc = *cp2 - MP_LIMB_T_HIGHBIT; - while (cc == 0 && cp2 > cp) cc = *--cp2; + nulp = -1; - if (cc || ((ap[0] >> sh) & ONE)) goto add_one_ulp; - /* mant(c) != 1/2 or mant(c) = 1/2: add 1 iff lsb(a)=1 */ - } - } + inv = rnd_mode == GMP_RNDN ? (sh ? ONE<<(sh-1) : MP_LIMB_T_HIGHBIT) : 0; + if (sh) + { + bb = *ap & ((ONE<<sh)-1); + *ap &= ~bb; /* truncate last bits */ + if (inv) + { + r = bb >> (sh-1); + bb ^= inv; + inv = 0; + } + if (difw > 0) + { + /* c overlaps with the lowest limb of a */ + MPFR_ASSERTN(difs > 0); + cc = cp[--cn]; + bb += cc >> (BITS_PER_MP_LIMB - difs); + } + if (bb >= ONE<<sh) + { + nulp = 1; + bb -= ONE<<sh; + } + else if (bb < (ONE<<sh)-1) + { + nulp = 0; + } } - else if ((MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) || - (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) { - - /* round up */ - if (MPFR_NOTZERO(c)) goto add_one_ulp; + + if (nulp < 0 || (nulp > 0 && bb == 0)) + while (1) /* look for carry and/or sticky bit */ + { + if (k == 0) + { + if (nulp < 0) + nulp = 0; + else /* nulp == 1 */ + { + if (difw > 0) + bb = cc << difs; + while (bb == 0 && cn) + bb = cp[--cn]; + } + break; + } + bb = bp[--k]; + if (inv) + { + r = bb >> (BITS_PER_MP_LIMB - 1); + bb ^= inv; + inv = 0; + } + if (difw >= 0) + { + mp_limb_t b2; + if (difw > 0) + { + b2 = bb; + bb = b2 + (cc << difs); + if (bb < b2) + nulp = 1; + } + if (cn == 0) + { + if (nulp < 0) + nulp = 0; + else /* nulp == 1 */ + while (bb == 0 && k) + bb = bp[--k]; + break; + } + cc = cp[--cn]; + if (difs) + { + difw = 1; + bb = b2 + (cc >> (BITS_PER_MP_LIMB - difs)); + } + else + { + bb = b2 + cc; + } + if (bb < b2) + nulp = 1; + } + else + difw++; + if (nulp > 0 && bb != 0) + break; + if (nulp < 0 && bb != (mp_limb_t) (-1)) + { + nulp = 0; + break; + } + } /* while */ + + /* if nulp == 1, sticky bit = bb != 0 */ + + if (rnd_mode == GMP_RNDN) + { /* r: old rounding bit of b + ulp(a)/2 were added if r = 0 + ulp(a)/2 were subtracted if r = 1 + r carry stbit --> nulp to add + 0 0 1 0 + 0 1 0 0 or 1 (to even) + 0 1 1 1 + 1 0 1 1 + 1 1 0 1 or 2 (to even) + 1 1 1 2 [*] + [*] or only 1 if adding the first ulp changes the exponent */ + + if (nulp == 0) + bb = 1; + nulp += r; + if (!bb && (ap[0] & (ONE<<sh))) + nulp++; /* sticky bit is 0 --> round to even */ } - /* in the other cases (round to zero, or up/down with sign -/+), - nothing to do */ - } - else + else { - /* MPFR_PREC(b)>MPFR_PREC(a) : we have to round b+c */ - k = bn - an; - - /* first copy the 'an' most significant limbs of b to a */ - MPN_COPY(ap, bp + k, an); - { /* treat all rounding modes together */ - mp_limb_t c2old; long int cout=0; int sign=0; - if (sh) { - cc = *ap & ((ONE<<sh)-1); - *ap &= ~cc; /* truncate last bits */ - } - else cc=0; - - dif += sh; - if (dif > 0) - { - cn--; - c2old = cp[cn]; /* last limb from c considered */ - cout += mpn_add_1(&cc, &cc, 1, c2old >> (BITS_PER_MP_LIMB-dif)); - } - else - c2 = c2old = 0; - if (sh && rnd_mode == GMP_RNDN) - cout -= mpn_sub_1 (&cc, &cc, 1, ONE << (sh - 1)); - if (cout == 0) - { - dif += BITS_PER_MP_LIMB; - while (cout==0 && (k || cn)) { - cout = (cc>(mp_limb_t)1) ? 2 : cc; - cc = (k) ? bp[--k] : 0; - if (sh==0) { - cout -= mpn_sub_1(&cc, &cc, 1, ONE << (BITS_PER_MP_LIMB-1)); - sh = 0; - } - /* next limb from c to consider is cp[cn-1], with lower part of - c2old */ - if (cn && (dif >= 0)) - { - c2 = c2old << dif; - cn--; - c2old = cp[cn]; - c2 += c2old >> (BITS_PER_MP_LIMB - dif); - } - else - dif += BITS_PER_MP_LIMB; - cout += mpn_add_1(&cc, &cc, 1, c2); - } - } - if (cout==0) cout=(cc!=0); - sign = (MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) - || (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD) || (rnd_mode==GMP_RNDN); - /* round towards infinity if dif=1, towards zero otherwise */ - if ((sign==1) && (cout>0)) goto add_one_ulp; - else if (rnd_mode==GMP_RNDN && cout==0 && (*ap & (ONE<<sh))) - goto add_one_ulp; + if ((nulp == 0 || bb != 0) && + ((rnd_mode == GMP_RNDU && MPFR_ISNONNEG(b)) + || (rnd_mode == GMP_RNDD && MPFR_ISNEG(b)))) + nulp++; } - } - } - else { + goto add_one_ulp; + } /* bq > aq */ + } /* aq <= diff_exp */ + else { /* diff_exp < MPFR_PREC(a) : c overlaps with a by dif bits */ /* first copy upper part of c into a (after shift) */ unsigned char overlap; - + + dif = MPFR_PREC(a) - diff_exp; k = (dif-1)/BITS_PER_MP_LIMB + 1; /* only the highest k limbs from c have to be considered */ cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB + 1; @@ -212,7 +291,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) #ifdef DEBUG printf("MPFR_PREC(c)=%d\n", MPFR_PREC(c)); #endif - if (dif<=MPFR_PREC(c)) { + if (dif<=MPFR_PREC(c)) { /* c has to be truncated */ dif = dif % BITS_PER_MP_LIMB; dif = (dif) ? BITS_PER_MP_LIMB-dif-sh : -sh; @@ -231,7 +310,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) else MPN_COPY(ap, cp+(cn-k), k); overlap=1; } - else { + else { /* c is not truncated, but we have to fill low limbs with 0 */ @@ -239,9 +318,9 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) overlap = diff_exp%BITS_PER_MP_LIMB; /* warning: a shift of zero bit is not allowed */ - MPN_ZERO(ap, an-k-cn); - if (overlap) { - cc=mpn_rshift(ap+(an-k-cn), cp, cn, overlap); + MPN_ZERO(ap, an-k-cn); + if (overlap) { + cc=mpn_rshift(ap+(an-k-cn), cp, cn, overlap); if (an-k-cn>0) ap[an-k-cn-1]=cc; } else MPN_COPY(ap+(an-k-cn), cp, cn); @@ -258,18 +337,18 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) } else /* MPFR_PREC(b) > MPFR_PREC(a): we have to truncate b */ - cc = mpn_add_n(ap, ap, bp+(bn-an), an); + cc = mpn_add_n(ap, ap, bp+(bn-an), an); - if (cc) { + if (cc) { /* shift one bit to the right */ c3 = (ap[0]&1) && (MPFR_PREC(a)%BITS_PER_MP_LIMB==0); mpn_rshift(ap, ap, an, 1); - ap[an-1] += ONE<<(BITS_PER_MP_LIMB-1); + ap[an-1] += MP_LIMB_T_HIGHBIT; MPFR_EXP(a)++; } - + /* remains to do the rounding */ #ifdef DEBUG @@ -279,7 +358,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) || (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) { int kc; - + /* four cases: overlap = (0) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) <= MPFR_PREC(a) (1) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) > MPFR_PREC(a) @@ -309,7 +388,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) } if (kc==0 && dif) { /* it still remains cp[0]<<(BITS_PER_MP_LIMB-dif) */ - if (k!=0) cout += mpn_add_1(&cc, bp+(--k), 1, + if (k!=0) cout += mpn_add_1(&cc, bp+(--k), 1, cp[0]<<(BITS_PER_MP_LIMB-dif)); else cc = cp[0]<<(BITS_PER_MP_LIMB-dif); if ((cout==0 && cc==0) || (~cout==0 && ~cc==0)) cout=cc; @@ -333,21 +412,21 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) case 0: /* only b to round */ k=bn-an; dif=0; goto to_nearest; - + /* otherwise the result is exact: nothing to do */ } } /* else nothing to do: round towards zero, i.e. truncate last sh bits */ - else + else *ap &= ~((ONE<<sh)-1); } goto end_of_add; - + to_nearest: /* 0 <= sh < BITS_PER_MP_LIMB : number of bits of a to truncate bp[k] : last significant limb from b bn : number of limbs of b */ - /* c3=1 whenever b+c gave a carry out in most significant limb + /* c3=1 whenever b+c gave a carry out in most significant limb and the least significant bit (shifted right) was 1. This can occur only when BITS_PER_MP_LIMB divides MPFR_PREC(a), i.e. sh=0. @@ -358,7 +437,7 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) c2 = (rnd_mode==GMP_RNDN) ? ONE<<(sh-1) : 0; } else /* sh=0: no bit to truncate */ - { + { if (k) cc = bp[--k]; else cc = 0; c2 = (rnd_mode==GMP_RNDN) ? ONE<<(BITS_PER_MP_LIMB-1) : 0; if (c3 && (cc || c2==0)) cc=c2+1; /* will force adding one ulp */ @@ -376,20 +455,35 @@ mpfr_add1 (a, b, c, rnd_mode, diff_exp) goto end_of_add; add_one_ulp: /* add one unit in last place to a */ - cc = mpn_add_1(ap, ap, an, ONE<<sh); - if (cc) { - ap[an-1] = (mp_limb_t)1 << (BITS_PER_MP_LIMB-1); - MPFR_EXP(a)++; + while (nulp--) + { + cc = mpn_add_1(ap, ap, an, ONE<<sh); + if (cc) + { + mp_exp_t exp = MPFR_EXP(a); + if (exp == __mpfr_emax) + { + mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); + break; + } + else + { + MPFR_EXP(a)++; + ap[an-1] = MP_LIMB_T_HIGHBIT; + if (rnd_mode == GMP_RNDN) + break; /* because ulp is doubled */ + } + } } end_of_add: - TMP_FREE(marker); + TMP_FREE(marker); return; } void #if __STDC__ -mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) +mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) #else mpfr_add (a, b, c, rnd_mode) mpfr_ptr a; @@ -398,62 +492,76 @@ mpfr_add (a, b, c, rnd_mode) mp_rnd_t rnd_mode; #endif { - int diff_exp; - if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { MPFR_SET_NAN(a); return; } - + MPFR_CLEAR_NAN(a); - if (MPFR_IS_INF(b)) - { - if (MPFR_IS_INF(c)) + if (MPFR_IS_INF(b)) + { + if (MPFR_IS_INF(c)) { - if (MPFR_SIGN(b) == MPFR_SIGN(c)) - { - MPFR_SET_INF(a); + if (MPFR_SIGN(b) == MPFR_SIGN(c)) + { + MPFR_SET_INF(a); if (MPFR_SIGN(a) != MPFR_SIGN(b)) MPFR_CHANGE_SIGN(a); } else - MPFR_SET_NAN(a); + MPFR_SET_NAN(a); } else { - MPFR_SET_INF(a); + MPFR_SET_INF(a); if (MPFR_SIGN(b) != MPFR_SIGN(a)) MPFR_CHANGE_SIGN(a); } return; } - else + else if (MPFR_IS_INF(c)) { - MPFR_SET_INF(a); + MPFR_SET_INF(a); if (MPFR_SIGN(c) != MPFR_SIGN(a)) MPFR_CHANGE_SIGN(a); return; } - MPFR_CLEAR_INF(a); /* clear Inf flag */ + if (MPFR_IS_ZERO(b)) { mpfr_set(a, c, rnd_mode); return; } + if (MPFR_IS_ZERO(c)) { mpfr_set(a, b, rnd_mode); return; } - if (!MPFR_NOTZERO(b)) { mpfr_set(a, c, rnd_mode); return; } - if (!MPFR_NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; } + MPFR_CLEAR_INF(a); /* clear Inf flag */ - diff_exp = MPFR_EXP(b)-MPFR_EXP(c); - if (MPFR_SIGN(b) * MPFR_SIGN(c) < 0) { /* signs differ, it's a subtraction */ - if (diff_exp<0) { - mpfr_sub1(a, c, b, rnd_mode, -diff_exp); + if (MPFR_SIGN(b) != MPFR_SIGN(c)) + { /* signs differ, it's a subtraction */ + if (MPFR_EXP(b) < MPFR_EXP(c)) + { + mpfr_sub1(a, c, b, rnd_mode, + (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b)); } - else if (diff_exp>0) mpfr_sub1(a, b, c, rnd_mode, diff_exp); - else { /* diff_exp=0 */ - diff_exp = mpfr_cmp3(b,c,-1); + else if (MPFR_EXP(b) > MPFR_EXP(c)) + { + mpfr_sub1(a, b, c, rnd_mode, + (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c)); + } + else + { /* MPFR_EXP(b) == MPFR_EXP(c) */ + int diff_exp = mpfr_cmp3(b,c,-1); /* if b>0 and diff_exp>0 or b<0 and diff_exp<0: abs(b) > abs(c) */ if (diff_exp==0) MPFR_SET_ZERO(a); else if (diff_exp * MPFR_SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0); else mpfr_sub1(a, c, b, rnd_mode, 0); } } - else /* signs are equal, it's an addition */ - if (diff_exp<0) mpfr_add1(a, c, b, rnd_mode, -diff_exp); - else mpfr_add1(a, b, c, rnd_mode, diff_exp); + else + { /* signs are equal, it's an addition */ + if (MPFR_EXP(b) < MPFR_EXP(c)) + { + mpfr_add1(a, c, b, rnd_mode, + (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b)); + } + else + { + mpfr_add1(a, b, c, rnd_mode, + (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c)); + } + } } - |