=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpfr/Attic/sub.c,v retrieving revision 1.1 retrieving revision 1.1.1.2 diff -u -p -r1.1 -r1.1.1.2 --- OpenXM_contrib/gmp/mpfr/Attic/sub.c 2000/09/09 14:12:19 1.1 +++ OpenXM_contrib/gmp/mpfr/Attic/sub.c 2003/08/25 16:06:07 1.1.1.2 @@ -1,483 +1,111 @@ /* mpfr_sub -- subtract two floating-point numbers -Copyright (C) 1999 PolKA project, Inria Lorraine and Loria +Copyright 2001 Free Software Foundation. +Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. The MPFR Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. -You should have received a copy of the GNU Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" +#include "mpfr-impl.h" -/* #define DEBUG */ - -extern void mpfr_add1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, - unsigned char, int)); - -/* put in ap[0]..ap[an-1] the value of bp[0]..bp[n-1] shifted by sh bits - to the left minus ap[0]..ap[n-1], with 0 <= sh < mp_bits_per_limb, and - returns the borrow. -*/ -mp_limb_t -#if __STDC__ -mpn_sub_lshift_n (mp_limb_t *ap, mp_limb_t *bp, int n, int sh, int an) -#else -mpn_sub_lshift_n (ap, bp, n, sh, an) mp_limb_t *ap, *bp; int n,sh,an; -#endif +int +mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) { - mp_limb_t c, bh; + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } - /* shift b to the left */ - if (sh) bh = mpn_lshift(bp, bp, n, sh); - c = (n=abs(c), diff_exp>=0 */ -void -#if __STDC__ -mpfr_sub1(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, unsigned char rnd_mode, int diff_exp) -#else -mpfr_sub1(a, b, c, rnd_mode, diff_exp) - mpfr_ptr a; - mpfr_srcptr b; - mpfr_srcptr c; - unsigned char rnd_mode; - int diff_exp; -#endif -{ - mp_limb_t *ap, *bp, *cp, cc, c2; unsigned int an,bn,cn; - int sh,dif,k,cancel,cancel1,cancel2; - TMP_DECL(marker); - -#ifdef DEBUG - printf("b= "); if (SIGN(b)>=0) putchar(' '); - mpfr_print_raw(b); putchar('\n'); - printf("c= "); if (SIGN(c)>=0) putchar(' '); - for (k=0;k=prec(a), i.e. c only affects the last bit - through rounding */ - dif = PREC(a)-diff_exp; -#ifdef DEBUG - printf("PREC(a)=%d an=%u PREC(b)=%d bn=%u PREC(c)=%d diff_exp=%u dif=%d cancel=%d\n", - PREC(a),an,PREC(b),bn,PREC(c),diff_exp,dif,cancel); -#endif - if (dif<=0) { /* diff_exp>=PREC(a): c does not overlap with a */ - /* either PREC(b)<=PREC(a), and we can copy the mantissa of b directly - into that of a, or PREC(b)>PREC(a) and we have to round b-c */ - if (PREC(b)<=PREC(a)+cancel) { - if (cancel2) mpn_lshift(ap+(an-bn+cancel1), bp, bn-cancel1, cancel2); - else MPN_COPY(ap+(an-bn+cancel1), bp, bn-cancel1); - /* fill low significant limbs with zero */ - MPN_ZERO(ap, an-bn+cancel1); - /* now take c into account */ - if (rnd_mode==GMP_RNDN) { /* to nearest */ - /* if diff_exp > PREC(a), no change */ - if (diff_exp==PREC(a)) { - /* if c is not zero, then as it is normalized, we have to subtract - one to the lsb of a if c>1/2, or c=1/2 and lsb(a)=1 (round to - even) */ - if (NOTZERO(c)) { /* c is not zero */ - /* check whether mant(c)=1/2 or not */ - cc = *cp - ((mp_limb_t)1<<(mp_bits_per_limb-1)); - if (cc==0) { - bp = cp+(PREC(c)-1)/mp_bits_per_limb; - while (cp 1/2 or mant(c) = 1/2: subtract 1 iff lsb(a)=1 */ - } - } - else if (ap[an-1]==0) { /* case b=2^n */ - ap[an-1] = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1); - EXP(a)++; - } + else + if (MPFR_IS_INF(c)) + { + MPFR_SET_INF(a); + if (MPFR_SIGN(c) == MPFR_SIGN(a)) + MPFR_CHANGE_SIGN(a); + MPFR_RET(0); /* exact */ } - else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || - (ISNEG(b) && rnd_mode==GMP_RNDD)) { - /* round up: nothing to do */ - if (ap[an-1]==0) { /* case b=2^n */ - ap[an-1] = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1); - EXP(a)++; - } - } - else { - /* round down: subtract 1 ulp iff c<>0 */ - if (NOTZERO(c)) goto sub_one_ulp; - } - } - else { /* PREC(b)>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); - if (rnd_mode==GMP_RNDN) { /* to nearest */ - /* first check whether the truncated bits from b are 1/2*lsb(a) */ - if (sh) { - cc = *ap & (((mp_limb_t)1<0) { /* suppose sizeof(long)=sizeof(mp_limb_t) */ - goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */ - } - else if (cc==0) { - while (k>1 && cc==0) cc=bp[--k]; - /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */ - if (NOTZERO(c) || (*ap & ((mp_limb_t)1< round down, i.e. do nothing */ - } - else { /* round towards infinity or zero */ - if (sh) { - cc = *ap & (((mp_limb_t)1<-sh) ? cp[cn]>>(mp_bits_per_limb-dif-sh) : 0; - while (cc==c2 && (k || cn)) { - if (k) cc = bp[--k]; - if (cn) { - c2 = cp[cn]<<(dif+sh); - if (--cn) c2 += cp[cn]>>(mp_bits_per_limb-dif-sh); - } - } - dif = ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || - (ISNEG(b) && rnd_mode==GMP_RNDD)); - /* round towards infinity if dif=1, towards zero otherwise */ - if (dif && cc>c2) goto add_one_ulp; - else if (dif==0 && cc0) { - mpn_rshift(ap, cp+(cn-k), (k<=an) ? k : an, dif); - if (k>an) ap[an-1] += cp[cn-k+an]<<(mp_bits_per_limb-dif); - } - else if (dif<0) { - cc = mpn_lshift(ap, cp+(cn-k), k, -dif); - if (k= k+1) - ap[0] += cp[cn-k-1]>>(mp_bits_per_limb+dif); - } - else MPN_COPY(ap, cp+(cn-k), k); - overlap=1; - } - else { /* c is not truncated, but we have to fill low limbs with 0 */ - MPN_ZERO(ap, k-cn); - overlap = cancel-diff_exp; -#ifdef DEBUG - printf("0:a="); mpfr_print_raw(a); putchar('\n'); - printf("overlap=%d\n",overlap); -#endif - if (overlap>=0) { - cn -= overlap/mp_bits_per_limb; - overlap %= mp_bits_per_limb; - /* warning: a shift of zero with mpn_lshift is not allowed */ - if (overlap) { - if (an>(mp_bits_per_limb-overlap); - } - else mpn_lshift(ap+(an-cn), cp, cn, overlap); - } - else MPN_COPY(ap+(an-cn), cp, cn); - } - else { /* shift to the right by -overlap bits */ - overlap = -overlap; - k = overlap/mp_bits_per_limb; - overlap = overlap % mp_bits_per_limb; - if (overlap) cc = mpn_rshift(ap+(an-k-cn), cp, cn, overlap); - else { - MPN_COPY(ap+(an-k-cn), cp, cn); - cc = 0; - } - if (an>k+cn) ap[an-k-cn-1]=cc; - } - overlap=0; - } -#ifdef DEBUG - printf("1:a="); mpfr_print_raw(a); putchar('\n'); -#endif - /* here overlap=1 iff ulp(c) PREC(a): we have to truncate b */ - mpn_sub_lshift_n(ap, bp+(bn-an-cancel1), an, cancel2, an); - /* remains to do the rounding */ -#ifdef DEBUG - printf("2:a="); mpfr_print_raw(a); putchar('\n'); - printf("overlap=%d\n",overlap); -#endif - if (rnd_mode==GMP_RNDN) { /* to nearest */ - int kc; - /* four cases: overlap = - (0) PREC(b) > PREC(a) and diff_exp+PREC(c) <= PREC(a) - (1) PREC(b) > PREC(a) and diff_exp+PREC(c) > PREC(a) - (2) PREC(b) <= PREC(a) and diff_exp+PREC(c) <= PREC(a) - (3) PREC(b) <= PREC(a) and diff_exp+PREC(c) > PREC(a) */ - switch (overlap) - { - case 1: /* both b and c to round */ - kc = cn-k; /* remains kc limbs from c */ - k = bn-an; /* remains k limbs from b */ - /* truncate last bits and store the difference with 1/2*ulp in cc */ - cc = *ap & (((mp_limb_t)1<>dif) + - (cp[kc+1]<<(mp_bits_per_limb-dif))); - if (cc==0 || cc==-1) cc=c2; - } - if ((long)cc>0) goto add_one_ulp; - else if ((long)cc<-1) goto end_of_sub; /* the carry can be at most 1 */ - else if (kc==0) goto round_b; - /* else round c: go through */ - case 3: /* only c to round */ - bp=cp; k=cn-k; goto to_nearest; - case 0: /* only b to round */ - round_b: - k=bn-an; dif=0; goto to_nearest; - /* otherwise the result is exact: nothing to do */ - } - } - else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || - (ISNEG(b) && rnd_mode==GMP_RNDD)) { - cc = *ap & (((mp_limb_t)1<>dif); - if (dif) cc -= (cp[kc+1]<<(mp_bits_per_limb-dif)); - } - if (cc) goto add_one_ulp; - else if (kc==0) goto round_b2; - /* else round c: go through */ - case 3: /* only c to round: nothing to do */ - /* while (kc) if (cp[--kc]) goto add_one_ulp; */ - /* if dif>0 : remains to check last dif bits from c */ - /* if (dif>0 && (cp[0]<<(mp_bits_per_limb-dif))) goto add_one_ulp; */ - break; - case 0: /* only b to round */ - round_b2: - k=bn-an; - while (k) if (bp[--k]) goto add_one_ulp; - /* otherwise the result is exact: nothing to do */ - } - } - } - /* else round to zero: remove 1 ulp if neglected bits from b-c are < 0 */ - else { - cc = *ap & (((mp_limb_t)1< no borrow */ - case 1: /* both b and c are truncated */ - break; - case 3: /* only c is truncated */ - cn -= k; /* take into account cp[0]..cp[cn-1] shifted by dif bits - to the right */ - cc = (dif>0) ? cp[cn]<<(mp_bits_per_limb-dif) : 0; - while (cc==0 && cn>0) cc = cp[--cn]; - if (cc) goto sub_one_ulp; - break; - } - } - } - } - goto end_of_sub; - - to_nearest: /* 0 <= sh < mp_bits_per_limb : number of bits of a to truncate - bp[k] : last significant limb from b */ -#ifdef DEBUG -mpfr_print_raw(a); putchar('\n'); -#endif - if (sh) { - cc = *ap & (((mp_limb_t)1<c2) goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */ - else if (cc==c2) { - cc=0; while (k && cc==0) cc=bp[--k]; -#ifdef DEBUG - printf("cc=%lu\n",cc); -#endif - /* special case of rouding c shifted to the right */ - if (cc==0 && dif>0) cc=bp[0]<<(mp_bits_per_limb-dif); - /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */ - if (bp!=cp) { - if (cc || (*ap & ((mp_limb_t)1<0, do nothing */ - if (cc==0 && (*ap & ((mp_limb_t)1< 0) ? -1 : 1) : + ((MPFR_SIGN(b) > 0 && MPFR_SIGN(c) < 0) ? 1 : -1))) + MPFR_CHANGE_SIGN(a); + MPFR_CLEAR_INF(a); + MPFR_SET_ZERO(a); + MPFR_RET(0); /* 0 - 0 is exact */ + } + return mpfr_neg (a, c, rnd_mode); + } - end_of_sub: -#ifdef DEBUG -printf("b-c="); if (SIGN(a)>0) putchar(' '); mpfr_print_raw(a); putchar('\n'); -#endif - TMP_FREE(marker); - return; -} + if (MPFR_IS_ZERO(c)) + { + return mpfr_set (a, b, rnd_mode); + } -void -#if __STDC__ -mpfr_sub(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, unsigned char rnd_mode) -#else -mpfr_sub(a, b, c, rnd_mode) - mpfr_ptr a; - mpfr_srcptr b; - mpfr_srcptr c; - unsigned char rnd_mode; -#endif -{ - int diff_exp; + MPFR_CLEAR_INF(a); - if (FLAG_NAN(b) || FLAG_NAN(c)) { SET_NAN(a); return; } - - if (!NOTZERO(b)) { mpfr_neg(a, c, rnd_mode); return; } - if (!NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; } - - diff_exp = EXP(b)-EXP(c); - if (SIGN(b) == SIGN(c)) { /* signs are equal, it's a real subtraction */ - if (diff_exp<0) { - /* exchange rounding modes towards +/- infinity */ - if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD; - else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU; - mpfr_sub1(a, c, b, rnd_mode, -diff_exp); - CHANGE_SIGN(a); + if (MPFR_SIGN(b) == MPFR_SIGN(c)) + { /* signs are equal, it's a real subtraction */ + return mpfr_sub1(a, b, c, rnd_mode, 1); } - 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); - /* if b>0 and diff_exp>0 or b<0 and diff_exp<0: abs(b) > abs(c) */ - if (diff_exp==0) SET_ZERO(a); - else if (diff_exp*SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0); - else { - /* exchange rounding modes towards +/- infinity */ - if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD; - else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU; - mpfr_sub1(a, c, b, rnd_mode, 0); - CHANGE_SIGN(a); - } + else + { /* signs differ, it's an addition */ + if (MPFR_EXP(b) < MPFR_EXP(c)) + { /* exchange rounding modes towards +/- infinity */ + int inexact; + if (rnd_mode == GMP_RNDU) + rnd_mode = GMP_RNDD; + else if (rnd_mode == GMP_RNDD) + rnd_mode = GMP_RNDU; + inexact = mpfr_add1(a, c, b, rnd_mode, + (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b)); + MPFR_CHANGE_SIGN(a); + return -inexact; + } + else + { + return mpfr_add1(a, b, c, rnd_mode, + (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c)); + } } - } - else /* signs differ, it's an addition */ - if (diff_exp<0) { - /* exchange rounding modes towards +/- infinity */ - if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD; - else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU; - mpfr_add1(a, c, b, rnd_mode, -diff_exp); - CHANGE_SIGN(a); - } - else mpfr_add1(a, b, c, rnd_mode, diff_exp); } -