Annotation of OpenXM_contrib/gmp/mpfr/sub.c, Revision 1.1
1.1 ! maekawa 1: /* mpfr_sub -- subtract 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: /* #define DEBUG */
! 28:
! 29: extern void mpfr_add1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
! 30: unsigned char, int));
! 31:
! 32: /* put in ap[0]..ap[an-1] the value of bp[0]..bp[n-1] shifted by sh bits
! 33: to the left minus ap[0]..ap[n-1], with 0 <= sh < mp_bits_per_limb, and
! 34: returns the borrow.
! 35: */
! 36: mp_limb_t
! 37: #if __STDC__
! 38: mpn_sub_lshift_n (mp_limb_t *ap, mp_limb_t *bp, int n, int sh, int an)
! 39: #else
! 40: mpn_sub_lshift_n (ap, bp, n, sh, an) mp_limb_t *ap, *bp; int n,sh,an;
! 41: #endif
! 42: {
! 43: mp_limb_t c, bh;
! 44:
! 45: /* shift b to the left */
! 46: if (sh) bh = mpn_lshift(bp, bp, n, sh);
! 47: c = (n<an) ? mpn_sub_n(ap, bp, ap, n) : mpn_sub_n(ap, bp+(n-an), ap, an);
! 48: /* shift back b to the right */
! 49: if (sh) {
! 50: mpn_rshift(bp, bp, n, sh);
! 51: bp[n-1] += bh<<(mp_bits_per_limb-sh);
! 52: }
! 53: return c;
! 54: }
! 55:
! 56: /* signs of b and c differ, abs(b)>=abs(c), diff_exp>=0 */
! 57: void
! 58: #if __STDC__
! 59: mpfr_sub1(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, unsigned char rnd_mode, int diff_exp)
! 60: #else
! 61: mpfr_sub1(a, b, c, rnd_mode, diff_exp)
! 62: mpfr_ptr a;
! 63: mpfr_srcptr b;
! 64: mpfr_srcptr c;
! 65: unsigned char rnd_mode;
! 66: int diff_exp;
! 67: #endif
! 68: {
! 69: mp_limb_t *ap, *bp, *cp, cc, c2; unsigned int an,bn,cn;
! 70: int sh,dif,k,cancel,cancel1,cancel2;
! 71: TMP_DECL(marker);
! 72:
! 73: #ifdef DEBUG
! 74: printf("b= "); if (SIGN(b)>=0) putchar(' ');
! 75: mpfr_print_raw(b); putchar('\n');
! 76: printf("c= "); if (SIGN(c)>=0) putchar(' ');
! 77: for (k=0;k<diff_exp;k++) putchar(' '); mpfr_print_raw(c);
! 78: putchar('\n');
! 79: printf("b=%1.20e c=%1.20e\n",mpfr_get_d(b),mpfr_get_d(c));
! 80: #endif
! 81: cancel = mpfr_cmp2(b, c);
! 82: /* we have to take into account the first (PREC(a)+cancel) bits from b */
! 83: cancel1 = cancel/mp_bits_per_limb; cancel2 = cancel%mp_bits_per_limb;
! 84: TMP_MARK(marker);
! 85: ap = MANT(a);
! 86: bp = MANT(b);
! 87: cp = MANT(c);
! 88: if (ap == bp) {
! 89: bp = (mp_ptr) TMP_ALLOC(ABSSIZE(b) * BYTES_PER_MP_LIMB);
! 90: MPN_COPY (bp, ap, ABSSIZE(b));
! 91: if (ap == cp) { cp = bp; }
! 92: }
! 93: else if (ap == cp)
! 94: {
! 95: cp = (mp_ptr) TMP_ALLOC (ABSSIZE(c) * BYTES_PER_MP_LIMB);
! 96: MPN_COPY(cp, ap, ABSSIZE(c));
! 97: }
! 98: an = (PREC(a)-1)/mp_bits_per_limb+1; /* number of significant limbs of a */
! 99: sh = an*mp_bits_per_limb-PREC(a); /* non-significant bits in low limb */
! 100: bn = (PREC(b)-1)/mp_bits_per_limb+1; /* number of significant limbs of b */
! 101: cn = (PREC(c)-1)/mp_bits_per_limb + 1;
! 102: EXP(a) = EXP(b)-cancel;
! 103: /* adjust sign to that of b */
! 104: if (SIGN(a)*SIGN(b)<0) CHANGE_SIGN(a);
! 105: /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit
! 106: through rounding */
! 107: dif = PREC(a)-diff_exp;
! 108: #ifdef DEBUG
! 109: printf("PREC(a)=%d an=%u PREC(b)=%d bn=%u PREC(c)=%d diff_exp=%u dif=%d cancel=%d\n",
! 110: PREC(a),an,PREC(b),bn,PREC(c),diff_exp,dif,cancel);
! 111: #endif
! 112: if (dif<=0) { /* diff_exp>=PREC(a): c does not overlap with a */
! 113: /* either PREC(b)<=PREC(a), and we can copy the mantissa of b directly
! 114: into that of a, or PREC(b)>PREC(a) and we have to round b-c */
! 115: if (PREC(b)<=PREC(a)+cancel) {
! 116: if (cancel2) mpn_lshift(ap+(an-bn+cancel1), bp, bn-cancel1, cancel2);
! 117: else MPN_COPY(ap+(an-bn+cancel1), bp, bn-cancel1);
! 118: /* fill low significant limbs with zero */
! 119: MPN_ZERO(ap, an-bn+cancel1);
! 120: /* now take c into account */
! 121: if (rnd_mode==GMP_RNDN) { /* to nearest */
! 122: /* if diff_exp > PREC(a), no change */
! 123: if (diff_exp==PREC(a)) {
! 124: /* if c is not zero, then as it is normalized, we have to subtract
! 125: one to the lsb of a if c>1/2, or c=1/2 and lsb(a)=1 (round to
! 126: even) */
! 127: if (NOTZERO(c)) { /* c is not zero */
! 128: /* check whether mant(c)=1/2 or not */
! 129: cc = *cp - ((mp_limb_t)1<<(mp_bits_per_limb-1));
! 130: if (cc==0) {
! 131: bp = cp+(PREC(c)-1)/mp_bits_per_limb;
! 132: while (cp<bp && cc==0) cc = *++cp;
! 133: }
! 134: if (cc || (ap[an-1] & (mp_limb_t)1<<sh)) goto sub_one_ulp;
! 135: /* mant(c) > 1/2 or mant(c) = 1/2: subtract 1 iff lsb(a)=1 */
! 136: }
! 137: }
! 138: else if (ap[an-1]==0) { /* case b=2^n */
! 139: ap[an-1] = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
! 140: EXP(a)++;
! 141: }
! 142: }
! 143: else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) ||
! 144: (ISNEG(b) && rnd_mode==GMP_RNDD)) {
! 145: /* round up: nothing to do */
! 146: if (ap[an-1]==0) { /* case b=2^n */
! 147: ap[an-1] = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
! 148: EXP(a)++;
! 149: }
! 150: }
! 151: else {
! 152: /* round down: subtract 1 ulp iff c<>0 */
! 153: if (NOTZERO(c)) goto sub_one_ulp;
! 154: }
! 155: }
! 156: else { /* PREC(b)>PREC(a) : we have to round b-c */
! 157: k=bn-an;
! 158: /* first copy the 'an' most significant limbs of b to a */
! 159: MPN_COPY(ap, bp+k, an);
! 160: if (rnd_mode==GMP_RNDN) { /* to nearest */
! 161: /* first check whether the truncated bits from b are 1/2*lsb(a) */
! 162: if (sh) {
! 163: cc = *ap & (((mp_limb_t)1<<sh)-1);
! 164: *ap &= ~cc; /* truncate last bits */
! 165: cc -= (mp_limb_t)1<<(sh-1);
! 166: }
! 167: else /* no bit to truncate */
! 168: cc = bp[--k] - ((mp_limb_t)1<<(mp_bits_per_limb-1));
! 169: if ((long)cc>0) { /* suppose sizeof(long)=sizeof(mp_limb_t) */
! 170: goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
! 171: }
! 172: else if (cc==0) {
! 173: while (k>1 && cc==0) cc=bp[--k];
! 174: /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
! 175: if (NOTZERO(c) || (*ap & ((mp_limb_t)1<<sh))) goto sub_one_ulp;
! 176: /* if trunc(b)-c is exactly 1/2*lsb(a) : round to even lsb */
! 177: }
! 178: /* if cc<0 : trunc(b) < 1/2*lsb(a) -> round down, i.e. do nothing */
! 179: }
! 180: else { /* round towards infinity or zero */
! 181: if (sh) {
! 182: cc = *ap & (((mp_limb_t)1<<sh)-1);
! 183: *ap &= ~cc; /* truncate last bits */
! 184: }
! 185: else cc=0;
! 186: cn--;
! 187: c2 = (dif>-sh) ? cp[cn]>>(mp_bits_per_limb-dif-sh) : 0;
! 188: while (cc==c2 && (k || cn)) {
! 189: if (k) cc = bp[--k];
! 190: if (cn) {
! 191: c2 = cp[cn]<<(dif+sh);
! 192: if (--cn) c2 += cp[cn]>>(mp_bits_per_limb-dif-sh);
! 193: }
! 194: }
! 195: dif = ((ISNONNEG(b) && rnd_mode==GMP_RNDU) ||
! 196: (ISNEG(b) && rnd_mode==GMP_RNDD));
! 197: /* round towards infinity if dif=1, towards zero otherwise */
! 198: if (dif && cc>c2) goto add_one_ulp;
! 199: else if (dif==0 && cc<c2) goto sub_one_ulp;
! 200: }
! 201: }
! 202: }
! 203: else { /* case 2: diff_exp < PREC(a) : c overlaps with a by dif bits */
! 204: /* first copy upper part of c into a (after shift) */
! 205: int overlap;
! 206: dif += cancel;
! 207: k = (dif-1)/mp_bits_per_limb + 1; /* only the highest k limbs from c
! 208: have to be considered */
! 209: if (k<an) {
! 210: MPN_ZERO(ap+k, an-k); /* do it now otherwise ap[k] may be
! 211: destroyed in case dif<0 */
! 212: }
! 213: #ifdef DEBUG
! 214: printf("cancel=%d dif=%d k=%d cn=%d sh=%d\n",cancel,dif,k,cn,sh);
! 215: #endif
! 216: if (dif<=PREC(c)) { /* c has to be truncated */
! 217: dif = dif % mp_bits_per_limb;
! 218: dif = (dif) ? mp_bits_per_limb-dif-sh : -sh;
! 219: /* we have to shift by dif bits to the right */
! 220: if (dif>0) {
! 221: mpn_rshift(ap, cp+(cn-k), (k<=an) ? k : an, dif);
! 222: if (k>an) ap[an-1] += cp[cn-k+an]<<(mp_bits_per_limb-dif);
! 223: }
! 224: else if (dif<0) {
! 225: cc = mpn_lshift(ap, cp+(cn-k), k, -dif);
! 226: if (k<an) ap[k]=cc;
! 227: /* put the non-significant bits in low limb for further rounding */
! 228: if (cn >= k+1)
! 229: ap[0] += cp[cn-k-1]>>(mp_bits_per_limb+dif);
! 230: }
! 231: else MPN_COPY(ap, cp+(cn-k), k);
! 232: overlap=1;
! 233: }
! 234: else { /* c is not truncated, but we have to fill low limbs with 0 */
! 235: MPN_ZERO(ap, k-cn);
! 236: overlap = cancel-diff_exp;
! 237: #ifdef DEBUG
! 238: printf("0:a="); mpfr_print_raw(a); putchar('\n');
! 239: printf("overlap=%d\n",overlap);
! 240: #endif
! 241: if (overlap>=0) {
! 242: cn -= overlap/mp_bits_per_limb;
! 243: overlap %= mp_bits_per_limb;
! 244: /* warning: a shift of zero with mpn_lshift is not allowed */
! 245: if (overlap) {
! 246: if (an<cn) {
! 247: mpn_lshift(ap, cp+(cn-an), an, overlap);
! 248: ap[0] += cp[cn-an-1]>>(mp_bits_per_limb-overlap);
! 249: }
! 250: else mpn_lshift(ap+(an-cn), cp, cn, overlap);
! 251: }
! 252: else MPN_COPY(ap+(an-cn), cp, cn);
! 253: }
! 254: else { /* shift to the right by -overlap bits */
! 255: overlap = -overlap;
! 256: k = overlap/mp_bits_per_limb;
! 257: overlap = overlap % mp_bits_per_limb;
! 258: if (overlap) cc = mpn_rshift(ap+(an-k-cn), cp, cn, overlap);
! 259: else {
! 260: MPN_COPY(ap+(an-k-cn), cp, cn);
! 261: cc = 0;
! 262: }
! 263: if (an>k+cn) ap[an-k-cn-1]=cc;
! 264: }
! 265: overlap=0;
! 266: }
! 267: #ifdef DEBUG
! 268: printf("1:a="); mpfr_print_raw(a); putchar('\n');
! 269: #endif
! 270: /* here overlap=1 iff ulp(c)<ulp(a) */
! 271: /* then put high limbs to zero */
! 272: /* now add 'an' upper limbs of b in place */
! 273: if (PREC(b)<=PREC(a)+cancel) { int i, imax;
! 274: overlap += 2;
! 275: /* invert low limbs */
! 276: imax = (int)an-(int)bn+cancel1;
! 277: for (i=0;i<imax;i++) ap[i] = ~ap[i];
! 278: cc = (i) ? mpn_add_1(ap, ap, i, 1) : 1;
! 279: mpn_sub_lshift_n(ap+i, bp, bn-cancel1, cancel2, an);
! 280: mpn_sub_1(ap+i, ap+i, an-i, (mp_limb_t)1-cc);
! 281: }
! 282: else /* PREC(b) > PREC(a): we have to truncate b */
! 283: mpn_sub_lshift_n(ap, bp+(bn-an-cancel1), an, cancel2, an);
! 284: /* remains to do the rounding */
! 285: #ifdef DEBUG
! 286: printf("2:a="); mpfr_print_raw(a); putchar('\n');
! 287: printf("overlap=%d\n",overlap);
! 288: #endif
! 289: if (rnd_mode==GMP_RNDN) { /* to nearest */
! 290: int kc;
! 291: /* four cases: overlap =
! 292: (0) PREC(b) > PREC(a) and diff_exp+PREC(c) <= PREC(a)
! 293: (1) PREC(b) > PREC(a) and diff_exp+PREC(c) > PREC(a)
! 294: (2) PREC(b) <= PREC(a) and diff_exp+PREC(c) <= PREC(a)
! 295: (3) PREC(b) <= PREC(a) and diff_exp+PREC(c) > PREC(a) */
! 296: switch (overlap)
! 297: {
! 298: case 1: /* both b and c to round */
! 299: kc = cn-k; /* remains kc limbs from c */
! 300: k = bn-an; /* remains k limbs from b */
! 301: /* truncate last bits and store the difference with 1/2*ulp in cc */
! 302: cc = *ap & (((mp_limb_t)1<<sh)-1);
! 303: *ap &= ~cc; /* truncate last bits */
! 304: cc -= (mp_limb_t)1<<(sh-1);
! 305: while ((cc==0 || cc==-1) && k!=0 && kc!=0) {
! 306: kc--;
! 307: cc -= mpn_sub_1(&c2, bp+(--k), 1, (cp[kc]>>dif) +
! 308: (cp[kc+1]<<(mp_bits_per_limb-dif)));
! 309: if (cc==0 || cc==-1) cc=c2;
! 310: }
! 311: if ((long)cc>0) goto add_one_ulp;
! 312: else if ((long)cc<-1) goto end_of_sub; /* the carry can be at most 1 */
! 313: else if (kc==0) goto round_b;
! 314: /* else round c: go through */
! 315: case 3: /* only c to round */
! 316: bp=cp; k=cn-k; goto to_nearest;
! 317: case 0: /* only b to round */
! 318: round_b:
! 319: k=bn-an; dif=0; goto to_nearest;
! 320: /* otherwise the result is exact: nothing to do */
! 321: }
! 322: }
! 323: else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) ||
! 324: (ISNEG(b) && rnd_mode==GMP_RNDD)) {
! 325: cc = *ap & (((mp_limb_t)1<<sh)-1);
! 326: *ap &= ~cc; /* truncate last bits */
! 327: if (cc) goto add_one_ulp; /* will happen most of the time */
! 328: else { /* same four cases too */
! 329: int kc = cn-k; /* remains kc limbs from c */
! 330: switch (overlap)
! 331: {
! 332: case 1: /* both b and c to round */
! 333: k = bn-an; /* remains k limbs from b */
! 334: dif = diff_exp % mp_bits_per_limb;
! 335: while (cc==0 && k!=0 && kc!=0) {
! 336: kc--;
! 337: cc = bp[--k] - (cp[kc]>>dif);
! 338: if (dif) cc -= (cp[kc+1]<<(mp_bits_per_limb-dif));
! 339: }
! 340: if (cc) goto add_one_ulp;
! 341: else if (kc==0) goto round_b2;
! 342: /* else round c: go through */
! 343: case 3: /* only c to round: nothing to do */
! 344: /* while (kc) if (cp[--kc]) goto add_one_ulp; */
! 345: /* if dif>0 : remains to check last dif bits from c */
! 346: /* if (dif>0 && (cp[0]<<(mp_bits_per_limb-dif))) goto add_one_ulp; */
! 347: break;
! 348: case 0: /* only b to round */
! 349: round_b2:
! 350: k=bn-an;
! 351: while (k) if (bp[--k]) goto add_one_ulp;
! 352: /* otherwise the result is exact: nothing to do */
! 353: }
! 354: }
! 355: }
! 356: /* else round to zero: remove 1 ulp if neglected bits from b-c are < 0 */
! 357: else {
! 358: cc = *ap & (((mp_limb_t)1<<sh)-1);
! 359: *ap &= ~cc;
! 360: if (cc==0) { /* otherwise neglected difference cannot be < 0 */
! 361: /* take into account bp[0]..bp[bn-cancel1-1] shifted by cancel2 to left
! 362: and cp[0]..cp[cn-k-1] shifted by dif bits to right */
! 363: switch (overlap) {
! 364: case 0:
! 365: case 2:
! 366: break; /* c is not truncated ==> no borrow */
! 367: case 1: /* both b and c are truncated */
! 368: break;
! 369: case 3: /* only c is truncated */
! 370: cn -= k; /* take into account cp[0]..cp[cn-1] shifted by dif bits
! 371: to the right */
! 372: cc = (dif>0) ? cp[cn]<<(mp_bits_per_limb-dif) : 0;
! 373: while (cc==0 && cn>0) cc = cp[--cn];
! 374: if (cc) goto sub_one_ulp;
! 375: break;
! 376: }
! 377: }
! 378: }
! 379: }
! 380: goto end_of_sub;
! 381:
! 382: to_nearest: /* 0 <= sh < mp_bits_per_limb : number of bits of a to truncate
! 383: bp[k] : last significant limb from b */
! 384: #ifdef DEBUG
! 385: mpfr_print_raw(a); putchar('\n');
! 386: #endif
! 387: if (sh) {
! 388: cc = *ap & (((mp_limb_t)1<<sh)-1);
! 389: *ap &= ~cc; /* truncate last bits */
! 390: c2 = (mp_limb_t)1<<(sh-1);
! 391: }
! 392: else /* no bit to truncate */
! 393: { if (k) cc = bp[--k]; else cc = 0; c2 = (mp_limb_t)1<<(mp_bits_per_limb-1); }
! 394: #ifdef DEBUG
! 395: printf("cc=%lu c2=%lu k=%u\n",cc,c2,k);
! 396: #endif
! 397: if (cc>c2) goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
! 398: else if (cc==c2) {
! 399: cc=0; while (k && cc==0) cc=bp[--k];
! 400: #ifdef DEBUG
! 401: printf("cc=%lu\n",cc);
! 402: #endif
! 403: /* special case of rouding c shifted to the right */
! 404: if (cc==0 && dif>0) cc=bp[0]<<(mp_bits_per_limb-dif);
! 405: /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
! 406: if (bp!=cp) {
! 407: if (cc || (*ap & ((mp_limb_t)1<<sh))) goto add_one_ulp;
! 408: }
! 409: else {
! 410: /* subtract: if cc>0, do nothing */
! 411: if (cc==0 && (*ap & ((mp_limb_t)1<<sh))) goto add_one_ulp;
! 412: }
! 413: }
! 414: goto end_of_sub;
! 415:
! 416: sub_one_ulp:
! 417: cc = mpn_sub_1(ap, ap, an, (mp_limb_t)1<<sh);
! 418: goto end_of_sub;
! 419:
! 420: add_one_ulp: /* add one unit in last place to a */
! 421: cc = mpn_add_1(ap, ap, an, (mp_limb_t)1<<sh);
! 422:
! 423: end_of_sub:
! 424: #ifdef DEBUG
! 425: printf("b-c="); if (SIGN(a)>0) putchar(' '); mpfr_print_raw(a); putchar('\n');
! 426: #endif
! 427: TMP_FREE(marker);
! 428: return;
! 429: }
! 430:
! 431: void
! 432: #if __STDC__
! 433: mpfr_sub(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, unsigned char rnd_mode)
! 434: #else
! 435: mpfr_sub(a, b, c, rnd_mode)
! 436: mpfr_ptr a;
! 437: mpfr_srcptr b;
! 438: mpfr_srcptr c;
! 439: unsigned char rnd_mode;
! 440: #endif
! 441: {
! 442: int diff_exp;
! 443:
! 444: if (FLAG_NAN(b) || FLAG_NAN(c)) { SET_NAN(a); return; }
! 445:
! 446: if (!NOTZERO(b)) { mpfr_neg(a, c, rnd_mode); return; }
! 447: if (!NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; }
! 448:
! 449: diff_exp = EXP(b)-EXP(c);
! 450: if (SIGN(b) == SIGN(c)) { /* signs are equal, it's a real subtraction */
! 451: if (diff_exp<0) {
! 452: /* exchange rounding modes towards +/- infinity */
! 453: if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD;
! 454: else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU;
! 455: mpfr_sub1(a, c, b, rnd_mode, -diff_exp);
! 456: CHANGE_SIGN(a);
! 457: }
! 458: else if (diff_exp>0) mpfr_sub1(a, b, c, rnd_mode, diff_exp);
! 459: else { /* diff_exp=0 */
! 460: diff_exp = mpfr_cmp3(b,c,1);
! 461: /* if b>0 and diff_exp>0 or b<0 and diff_exp<0: abs(b) > abs(c) */
! 462: if (diff_exp==0) SET_ZERO(a);
! 463: else if (diff_exp*SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0);
! 464: else {
! 465: /* exchange rounding modes towards +/- infinity */
! 466: if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD;
! 467: else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU;
! 468: mpfr_sub1(a, c, b, rnd_mode, 0);
! 469: CHANGE_SIGN(a);
! 470: }
! 471: }
! 472: }
! 473: else /* signs differ, it's an addition */
! 474: if (diff_exp<0) {
! 475: /* exchange rounding modes towards +/- infinity */
! 476: if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD;
! 477: else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU;
! 478: mpfr_add1(a, c, b, rnd_mode, -diff_exp);
! 479: CHANGE_SIGN(a);
! 480: }
! 481: else mpfr_add1(a, b, c, rnd_mode, diff_exp);
! 482: }
! 483:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sub.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr