Annotation of OpenXM_contrib/gmp/mpfr/add.c, Revision 1.1
1.1 ! maekawa 1: /* mpfr_add -- add two floating-point numbers
! 2:
! 3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
! 4:
! 5: This file is part of the MPFR Library.
! 6:
! 7: The MPFR Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Library General Public License as published by
! 9: the Free Software Foundation; either version 2 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The MPFR Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Library General Public License
! 18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
! 19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: #include <stdio.h>
! 23: #include "gmp.h"
! 24: #include "gmp-impl.h"
! 25: #include "mpfr.h"
! 26:
! 27: extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
! 28: unsigned char, int));
! 29:
! 30: #define ONE ((mp_limb_t) 1)
! 31:
! 32: /* signs of b and c are supposed equal,
! 33: diff_exp is the difference between the exponents of b and c,
! 34: which is supposed >= 0 */
! 35:
! 36: void
! 37: #if __STDC__
! 38: mpfr_add1(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
! 39: unsigned char rnd_mode, int diff_exp)
! 40: #else
! 41: mpfr_add1(a, b, c, rnd_mode, diff_exp)
! 42: mpfr_ptr a;
! 43: mpfr_srcptr b;
! 44: mpfr_srcptr c;
! 45: unsigned char rnd_mode;
! 46: int diff_exp;
! 47: #endif
! 48: {
! 49: mp_limb_t *ap, *bp, *cp, cc, c2, c3=0; unsigned int an,bn,cn; int sh,dif,k;
! 50: TMP_DECL(marker);
! 51:
! 52: TMP_MARK(marker);
! 53: ap = MANT(a);
! 54: bp = MANT(b);
! 55: cp = MANT(c);
! 56: if (ap == bp) {
! 57: bp = (mp_ptr) TMP_ALLOC(ABSSIZE(b) * BYTES_PER_MP_LIMB);
! 58: MPN_COPY (bp, ap, ABSSIZE(b));
! 59: if (ap == cp) { cp = bp; }
! 60: }
! 61: else if (ap == cp)
! 62: {
! 63: cp = (mp_ptr) TMP_ALLOC (ABSSIZE(c) * BYTES_PER_MP_LIMB);
! 64: MPN_COPY(cp, ap, ABSSIZE(c));
! 65: }
! 66:
! 67: an = (PREC(a)-1)/mp_bits_per_limb+1; /* number of significant limbs of a */
! 68:
! 69: sh = an*mp_bits_per_limb-PREC(a); /* non-significant bits in low limb */
! 70: bn = (PREC(b)-1)/mp_bits_per_limb+1; /* number of significant limbs of b */
! 71: EXP(a) = EXP(b);
! 72:
! 73: if (SIGN(a)!=SIGN(b)) CHANGE_SIGN(a);
! 74:
! 75: /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit
! 76: through rounding */
! 77: dif = PREC(a)-diff_exp;
! 78:
! 79: if (dif<=0) {
! 80:
! 81: /* diff_exp>=PREC(a): c does not overlap with a */
! 82: /* either PREC(b)<=PREC(a), and we can copy the mantissa of b directly
! 83: into that of a, or PREC(b)>PREC(a) and we have to round b+c */
! 84:
! 85: if (PREC(b)<=PREC(a)) {
! 86:
! 87: MPN_COPY(ap+(an-bn), bp, bn);
! 88: /* fill low significant limbs with zero */
! 89:
! 90: for (bp=ap;bn<an;bn++) *bp++=0;
! 91:
! 92: /* now take c into account */
! 93: if (rnd_mode==GMP_RNDN) {
! 94:
! 95: /* to nearest */
! 96: /* if diff_exp > PREC(a), no change */
! 97:
! 98: if (diff_exp==PREC(a)) {
! 99:
! 100: /* if c is not zero, then as it is normalized, we have to add
! 101: one to the lsb of a if c>1/2, or c=1/2 and lsb(a)=1 (round to
! 102: even) */
! 103:
! 104: if (NOTZERO(c)) {
! 105:
! 106: /* c is not zero */
! 107: /* check whether mant(c)=1/2 or not */
! 108:
! 109: cc = *cp - (ONE<<(mp_bits_per_limb-1));
! 110: if (cc==0) {
! 111: bp = cp+(PREC(c)-1)/mp_bits_per_limb;
! 112: while (cp<bp && cc==0) cc = *++cp;
! 113: }
! 114:
! 115: if (cc || (ap[an-1] & (ONE<<sh))) goto add_one_ulp;
! 116: /* mant(c) != 1/2 or mant(c) = 1/2: add 1 iff lsb(a)=1 */
! 117: }
! 118: }
! 119: }
! 120: else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) ||
! 121: (ISNEG(b) && rnd_mode==GMP_RNDD)) {
! 122:
! 123: /* round up */
! 124: if (NOTZERO(c)) goto add_one_ulp;
! 125: }
! 126: /* in the other cases (round to zero, or up/down with sign -/+),
! 127: nothing to do */
! 128: }
! 129: else {
! 130:
! 131: /* PREC(b)>PREC(a) : we have to round b+c */
! 132: k=bn-an;
! 133:
! 134: /* first copy the 'an' most significant limbs of b to a */
! 135: MPN_COPY(ap, bp+k, an);
! 136: if (rnd_mode==GMP_RNDN) {
! 137:
! 138: /* to nearest */
! 139: /* first check whether the truncated bits from b are 1/2*lsb(a) */
! 140:
! 141: if (sh) {
! 142: cc = *ap & ((ONE<<sh)-1);
! 143: *ap &= ~cc; /* truncate last bits */
! 144: cc -= ONE<<(sh-1);
! 145: }
! 146: else /* no bit to truncate */
! 147: cc = bp[--k] - (ONE<<(mp_bits_per_limb-1));
! 148:
! 149: if ((long)cc>0) goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
! 150: else if (cc==0) {
! 151:
! 152: while (k>1 && cc==0) cc=bp[--k];
! 153:
! 154: /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
! 155: if (NOTZERO(c) || (*ap & (ONE<<sh))) goto add_one_ulp;
! 156: /* if trunc(b)+c is exactly 1/2*lsb(a) : round to even lsb */
! 157: }
! 158:
! 159: /* if cc<0 : trunc(b) < 1/2*lsb(a) -> round down, i.e. do nothing */
! 160: }
! 161: else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) ||
! 162: (ISNEG(b) && rnd_mode==GMP_RNDD)) {
! 163:
! 164: /* first check whether trunc(b)+c is zero or not */
! 165: if (sh) {
! 166: cc = *ap & ((ONE<<sh)-1); *ap &= ~cc; /* truncate last bits */
! 167: }
! 168: else cc = bp[--k] - (ONE<<(mp_bits_per_limb-1));
! 169: while (cc==0 && k>1) cc=bp[--k];
! 170: if (cc || NOTZERO(c)) goto add_one_ulp;
! 171: }
! 172:
! 173: /* in the other cases (round to zero, or up/down with sign -/+),
! 174: nothing to do, since b and c don't overlap, there can't be any
! 175: carry */
! 176:
! 177: }
! 178: }
! 179: else {
! 180: /* diff_exp < PREC(a) : c overlaps with a by dif bits */
! 181: /* first copy upper part of c into a (after shift) */
! 182: unsigned char overlap;
! 183:
! 184: k = (dif-1)/mp_bits_per_limb + 1; /* only the highest k limbs from c
! 185: have to be considered */
! 186: cn = (PREC(c)-1)/mp_bits_per_limb + 1;
! 187: MPN_ZERO(ap+k, an-k); /* do it now otherwise ap[k] may be destroyed
! 188: in case dif<0 */
! 189:
! 190: if (dif<=PREC(c)) {
! 191: /* c has to be truncated */
! 192: dif = dif % mp_bits_per_limb;
! 193: dif = (dif) ? mp_bits_per_limb-dif-sh : -sh;
! 194:
! 195: /* we have to shift by dif bits to the right */
! 196:
! 197: if (dif>0) mpn_rshift(ap, cp+(cn-k), k, dif);
! 198: else if (dif<0) {
! 199: ap[k] = mpn_lshift(ap, cp+(cn-k), k, -dif);
! 200:
! 201: /* put the non-significant bits in low limb for further rounding */
! 202:
! 203: if (cn >= k+1)
! 204: ap[0] += cp[cn-k-1]>>(mp_bits_per_limb+dif);
! 205: }
! 206: else MPN_COPY(ap, cp+(cn-k), k);
! 207: overlap=1;
! 208: }
! 209: else {
! 210:
! 211: /* c is not truncated, but we have to fill low limbs with 0 */
! 212:
! 213: k = diff_exp/mp_bits_per_limb;
! 214: overlap = diff_exp%mp_bits_per_limb;
! 215:
! 216: /* warning: a shift of zero bit is not allowed */
! 217: MPN_ZERO(ap, an-k-cn);
! 218: if (overlap) {
! 219: cc=mpn_rshift(ap+(an-k-cn), cp, cn, overlap);
! 220: if (an-k-cn>0) ap[an-k-cn-1]=cc;
! 221: }
! 222: else MPN_COPY(ap+(an-k-cn), cp, cn);
! 223: overlap=0;
! 224: }
! 225:
! 226: /* here overlap=1 iff ulp(c)<ulp(a) */
! 227: /* then put high limbs to zero */
! 228: /* now add 'an' upper limbs of b in place */
! 229:
! 230: if (PREC(b)<=PREC(a)) {
! 231: overlap += 2;
! 232: cc = mpn_add_n(ap+(an-bn), ap+(an-bn), bp, bn);
! 233: }
! 234: else
! 235: /* PREC(b) > PREC(a): we have to truncate b */
! 236: cc = mpn_add_n(ap, ap, bp+(bn-an), an);
! 237:
! 238: if (cc) {
! 239:
! 240: /* shift one bit to the right */
! 241:
! 242: c3 = (ap[0]&1) && (PREC(a)%mp_bits_per_limb==0);
! 243: mpn_rshift(ap, ap, an, 1);
! 244: ap[an-1] += ONE<<(mp_bits_per_limb-1);
! 245: EXP(a)++;
! 246: }
! 247:
! 248: /* remains to do the rounding */
! 249:
! 250: if (rnd_mode==GMP_RNDN) {
! 251:
! 252: /* to nearest */
! 253:
! 254: int kc;
! 255:
! 256: /* four cases: overlap =
! 257: (0) PREC(b) > PREC(a) and diff_exp+PREC(c) <= PREC(a)
! 258: (1) PREC(b) > PREC(a) and diff_exp+PREC(c) > PREC(a)
! 259: (2) PREC(b) <= PREC(a) and diff_exp+PREC(c) <= PREC(a)
! 260: (3) PREC(b) <= PREC(a) and diff_exp+PREC(c) > PREC(a) */
! 261:
! 262: switch (overlap)
! 263: {
! 264: case 1: /* both b and c to round */
! 265: kc = cn-k; /* remains kc limbs from c */
! 266: k = bn-an; /* remains k limbs from b */
! 267:
! 268: /* truncate last bits and store the difference with 1/2*ulp in cc */
! 269:
! 270: cc = *ap & ((ONE<<sh)-1);
! 271: *ap &= ~cc; /* truncate last bits */
! 272: cc -= ONE<<(sh-1);
! 273: while ((cc==0 || cc==-1) && k!=0 && kc!=0) {
! 274: kc--;
! 275: cc += mpn_add_1(&c2, bp+(--k), 1,(cp[kc+1]<<(mp_bits_per_limb-dif))
! 276: +(cp[kc]>>dif));
! 277: if (cc==0 || cc==-1) cc=c2;
! 278: }
! 279: if ((long)cc>0) goto add_one_ulp;
! 280: else if ((long)cc<-1)
! 281: { TMP_FREE(marker); return; /* the carry can be at most 1 */ }
! 282: else if (kc==0) goto round_b;
! 283:
! 284: /* else round c: go through */
! 285:
! 286: case 3: /* only c to round */
! 287: bp=cp; k=cn-k; goto to_nearest;
! 288:
! 289: case 0: /* only b to round */
! 290: round_b:
! 291: k=bn-an; dif=0; goto to_nearest;
! 292:
! 293: /* otherwise the result is exact: nothing to do */
! 294: }
! 295: }
! 296: else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) ||
! 297: (ISNEG(b) && rnd_mode==GMP_RNDD)) {
! 298: cc = *ap & ((ONE<<sh)-1);
! 299: *ap &= ~cc; /* truncate last bits */
! 300: if (cc || c3) goto add_one_ulp; /* will happen most of the time */
! 301: else {
! 302:
! 303: /* same four cases too */
! 304:
! 305: int kc = cn-k; /* remains kc limbs from c */
! 306: switch (overlap)
! 307: {
! 308: case 1: /* both b and c to round */
! 309: k = bn-an; /* remains k limbs from b */
! 310: while (cc==0 && k!=0 && kc!=0) {
! 311: kc--;
! 312: cc = mpn_add_1(&c2, bp+(--k), 1,(cp[kc+1]<<(mp_bits_per_limb-dif))
! 313: + (cp[kc]>>dif));
! 314: }
! 315: if (cc) goto add_one_ulp;
! 316: else if (kc==0) goto round_b2;
! 317: /* else round c: go through */
! 318: case 3: /* only c to round */
! 319: while (kc) if (cp[--kc]) goto add_one_ulp;
! 320: /* if dif>0 : remains to check last dif bits from c */
! 321: if (dif>0 && (cp[0]<<(mp_bits_per_limb-dif))) goto add_one_ulp;
! 322: break;
! 323: case 0: /* only b to round */
! 324: round_b2:
! 325: k=bn-an;
! 326: while (k) if (bp[--k]) goto add_one_ulp;
! 327: /* otherwise the result is exact: nothing to do */
! 328: }
! 329: }
! 330: }
! 331: /* else nothing to do: round towards zero, i.e. truncate last sh bits */
! 332: else
! 333: *ap &= ~((ONE<<sh)-1);
! 334: }
! 335: goto end_of_add;
! 336:
! 337: to_nearest: /* 0 <= sh < mp_bits_per_limb : number of bits of a to truncate
! 338: bp[k] : last significant limb from b */
! 339: if (sh) {
! 340: cc = *ap & ((ONE<<sh)-1);
! 341: *ap &= ~cc; /* truncate last bits */
! 342: c2 = ONE<<(sh-1);
! 343: }
! 344: else /* no bit to truncate */
! 345: { if (k) cc = bp[--k]; else cc = 0; c2 = ONE<<(mp_bits_per_limb-1); }
! 346: if (cc>c2) goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
! 347: else if (cc==c2) {
! 348: cc=0; while (k && cc==0) cc=bp[--k];
! 349: /* special case of rouding c shifted to the right */
! 350: if (cc==0 && dif>0) cc=cp[0]<<(mp_bits_per_limb-dif);
! 351: /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
! 352: if (cc || (*ap & (ONE<<sh))) goto add_one_ulp;
! 353: }
! 354: goto end_of_add;
! 355:
! 356: add_one_ulp: /* add one unit in last place to a */
! 357: cc = mpn_add_1(ap, ap, an, ONE<<sh);
! 358: if (cc) { fprintf(stderr, "carry(3) in mpfr_add\n"); exit(1); }
! 359:
! 360: end_of_add:
! 361: TMP_FREE(marker);
! 362: return;
! 363: }
! 364:
! 365: void
! 366: #if __STDC__
! 367: mpfr_add(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
! 368: unsigned char rnd_mode)
! 369: #else
! 370: mpfr_add(a, b, c, rnd_mode)
! 371: mpfr_ptr a;
! 372: mpfr_srcptr b;
! 373: mpfr_srcptr c;
! 374: unsigned char rnd_mode;
! 375: #endif
! 376: {
! 377: int diff_exp;
! 378:
! 379: if (FLAG_NAN(b) || FLAG_NAN(c)) {
! 380: SET_NAN(a); return;
! 381: }
! 382:
! 383: if (!NOTZERO(b)) { mpfr_set(a, c, rnd_mode); return; }
! 384: if (!NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; }
! 385:
! 386: diff_exp = EXP(b)-EXP(c);
! 387: if (SIGN(b) != SIGN(c)) { /* signs differ, it's a subtraction */
! 388: if (diff_exp<0) {
! 389: mpfr_sub1(a, c, b, rnd_mode, -diff_exp);
! 390: }
! 391: else if (diff_exp>0) mpfr_sub1(a, b, c, rnd_mode, diff_exp);
! 392: else { /* diff_exp=0 */
! 393: diff_exp = mpfr_cmp3(b,c,-1);
! 394: /* if b>0 and diff_exp>0 or b<0 and diff_exp<0: abs(b) > abs(c) */
! 395: if (diff_exp==0) SET_ZERO(a);
! 396: else if (diff_exp*SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0);
! 397: else mpfr_sub1(a, c, b, rnd_mode, 0);
! 398: }
! 399: }
! 400: else /* signs are equal, it's an addition */
! 401: if (diff_exp<0) mpfr_add1(a, c, b, rnd_mode, -diff_exp);
! 402: else mpfr_add1(a, b, c, rnd_mode, diff_exp);
! 403: }
! 404:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to add.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr