diff options
Diffstat (limited to 'sub.c')
| -rw-r--r-- | sub.c | 412 |
1 files changed, 412 insertions, 0 deletions
@@ -0,0 +1,412 @@ +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "mpfr.h" + +/* #define DEBUG2 */ + +#ifdef DEBUG2 +mp_limb_t *ap0, *ap1; +#define ck(x, dx, i) /* mp_limb_t *x; int dx; */ \ +{ \ + if ((x)<ap0 || (x)+(dx)-1>ap1) { \ + printf("error %d\n",(i)); \ + printf("%1.20e,%d,%1.20e,%d,%d,%d\n",mpfr_get_d(b),PREC(b),mpfr_get_d(c), \ + PREC(c),PREC(a),rnd_mode); \ + exit(1); } \ +} +#else +#define ck(x, dx, i) {} +#endif + +void mpfr_sub(a, b, c, rnd_mode) +mpfr_ptr a; mpfr_srcptr b, c; unsigned char rnd_mode; +{ + CHANGE_SIGN(c); + mpfr_add(a, b, c, rnd_mode); + CHANGE_SIGN(c); +} + +/* 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 mpn_sub_lshift_n (ap, bp, n, sh, an) mp_limb_t *ap, *bp; int n,sh,an; +{ + mp_limb_t c, bh; + + /* shift b to the left */ + if (sh) bh = mpn_lshift(bp, bp, n, sh); + c = (n<an) ? mpn_sub_n(ap, bp, ap, n) : mpn_sub_n(ap, bp+(n-an), ap, an); + /* shift back b to the right */ + if (sh) { + mpn_rshift(bp, bp, n, sh); + bp[n-1] += bh<<(mp_bits_per_limb-sh); + } + return c; +} + +/* signs of b and c differ, abs(b)>=abs(c), diff_exp>=0 */ +void mpfr_sub1(a, b, c, rnd_mode, diff_exp) +mpfr_t a, b, c; unsigned char rnd_mode; int diff_exp; +{ + mp_limb_t *ap, *bp, *cp, cc, c2; unsigned int an,bn,cn; + int sh,dif,k,cancel,cancel1,cancel2; + +#ifdef DEBUG2 + printf("b= "); if (SIGN(b)>=0) putchar(' '); + mpfr_print_raw(b); putchar('\n'); + printf("c= "); if (SIGN(c)>=0) putchar(' '); + for (k=0;k<diff_exp;k++) putchar(' '); mpfr_print_raw(c); + putchar('\n'); + printf("b=%1.20e c=%1.20e\n",mpfr_get_d(b),mpfr_get_d(c)); +#endif + cancel = mpfr_cmp2(b, c); + /* we have to take into account the first (PREC(a)+cancel) bits from b */ + cancel1 = cancel/mp_bits_per_limb; cancel2 = cancel%mp_bits_per_limb; + ap = MANT(a); + bp = MANT(b); + cp = MANT(c); + an = (PREC(a)-1)/mp_bits_per_limb+1; /* number of significant limbs of a */ +#ifdef DEBUG2 +ap0=ap; ap1=ap+an-1; +#endif + sh = an*mp_bits_per_limb-PREC(a); /* non-significant bits in low limb */ + bn = (PREC(b)-1)/mp_bits_per_limb+1; /* number of significant limbs of b */ + cn = (PREC(c)-1)/mp_bits_per_limb + 1; + EXP(a) = EXP(b)-cancel; + /* adjust sign to that of b */ + if (SIGN(a)*SIGN(b)<0) CHANGE_SIGN(a); + /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit + through rounding */ + dif = PREC(a)-diff_exp; +#ifdef DEBUG2 + printf("PREC(a)=%d an=%u PREC(b)=%d bn=%u PREC(c)=%d diff_exp=%u dif=%d\n", + PREC(a),an,PREC(b),bn,PREC(c),diff_exp,dif); +#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) { ck(ap+(an-bn+cancel1), bn-cancel1, 1); + mpn_lshift(ap+(an-bn+cancel1), bp, bn-cancel1, cancel2); } + else { ck(ap+(an-bn+cancel1), bn-cancel1, 2); + MPN_COPY(ap+(an-bn+cancel1), bp, bn-cancel1); } + /* fill low significant limbs with zero */ +ck(ap, an-bn+cancel1, 3); + 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 - (1<<(mp_bits_per_limb-1)); + if (cc==0) { + bp = cp+(PREC(c)-1)/mp_bits_per_limb; + while (cp<bp && cc==0) cc = *++cp; + } +ck(ap+an-1, 1, 4); + if (cc || (ap[an-1] & 1<<sh)) goto sub_one_ulp; + /* mant(c) > 1/2 or mant(c) = 1/2: subtract 1 iff lsb(a)=1 */ + } + } + } + else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || + (ISNEG(b) && rnd_mode==GMP_RNDD)) { + /* round up: nothing to do */ + } + 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 */ +ck(ap, an, 5); + 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) { +ck(ap, 1, 6); + cc = *ap & (((mp_limb_t)1<<sh)-1); + *ap &= ~cc; /* truncate last bits */ + cc -= (mp_limb_t)1<<(sh-1); + } + else /* no bit to truncate */ + cc = bp[--k] - (1<<(mp_bits_per_limb-1)); + if ((long)cc>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 */ +ck(ap, 1, 7); + if (NOTZERO(c) || (*ap & (1<<sh))) goto sub_one_ulp; + /* if trunc(b)-c is exactly 1/2*lsb(a) : round to even lsb */ + } + /* if cc<0 : trunc(b) < 1/2*lsb(a) -> round down, i.e. do nothing */ + } + else { /* round towards infinity or zero */ + if (sh) { +ck(ap, 1, 8); + cc = *ap & ((1<<sh)-1); + *ap &= ~cc; /* truncate last bits */ + } + else cc=0; + cn--; + c2 = (dif>-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 && cc<c2) goto sub_one_ulp; + } + } + } + else { /* case 2: diff_exp < PREC(a) : c overlaps with a by dif bits */ + /* first copy upper part of c into a (after shift) */ + int overlap; + dif += cancel; + k = (dif-1)/mp_bits_per_limb + 1; /* only the highest k limbs from c + have to be considered */ + if (k<an) { +ck(ap+k, an-k, 9); + MPN_ZERO(ap+k, an-k); /* do it now otherwise ap[k] may be + destroyed in case dif<0 */ + } +#ifdef DEBUG2 + printf("cancel=%d dif=%d k=%d cn=%d sh=%d\n",cancel,dif,k,cn,sh); +#endif + if (dif<=PREC(c)) { /* c has to be truncated */ + dif = dif % mp_bits_per_limb; + dif = (dif) ? mp_bits_per_limb-dif-sh : -sh; + /* we have to shift by dif bits to the right */ + if (dif>0) { +ck(ap, (k<=an) ? k : an, 10); + 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) { +ck(ap, k, 11); + cc = mpn_lshift(ap, cp+(cn-k), k, -dif); + if (k<an) { ck(ap+k, 1, 12); ap[k]=cc; } + /* put the non-significant bits in low limb for further rounding */ + ap[0] += cp[cn-k-1]>>(mp_bits_per_limb+dif); + } + else { ck(ap, k, 13); MPN_COPY(ap, cp+(cn-k), k); } + overlap=1; + } + else { /* c is not truncated, but we have to fill low limbs with 0 */ +ck(ap, k-cn, 14); + MPN_ZERO(ap, k-cn); + overlap = cancel-diff_exp; +#ifdef DEBUG2 + printf("0:a="); mpfr_print_raw(a); putchar('\n'); + printf("overlap=%d\n",overlap); +#endif + if (overlap>=0) { + cn -= overlap/mp_bits_per_limb; + /* warning: a shift of zero with mpn_lshift is not allowed */ + if (overlap) { + if (an<cn) { +ck(ap, an, 15); + mpn_lshift(ap, cp+(cn-an), an, overlap%mp_bits_per_limb); + ap[0] += cp[cn-an-1]>>(mp_bits_per_limb-overlap); + } + else { +ck(ap+(an-cn), cn, 15); + mpn_lshift(ap+(an-cn), cp, cn, overlap%mp_bits_per_limb); + } + } + else { ck(ap+(an-cn), cn, 15); MPN_COPY(ap+(an-cn), cp, cn); } + } + else { /* shift to the right by -overlap bits */ + overlap = -overlap; + k = overlap/mp_bits_per_limb; + cc=mpn_rshift(ap+(an-k-cn),cp,cn,overlap%mp_bits_per_limb); + if (an>k+cn) ap[an-k-cn-1]=cc; + } + overlap=0; + } +#ifdef DEBUG2 + printf("1:a="); mpfr_print_raw(a); putchar('\n'); +#endif + /* here overlap=1 iff ulp(c)<ulp(a) */ + /* then put high limbs to zero */ + /* now add 'an' upper limbs of b in place */ + if (PREC(b)<=PREC(a)+cancel) { int i, imax; + overlap += 2; + /* invert low limbs */ + imax = (int)an-(int)bn+cancel1; + for (i=0;i<imax;i++) { ck(ap+i, 1, 17); ap[i] = ~ap[i]; } +if (i) ck(ap, i, 18); + cc = (i) ? mpn_add_1(ap, ap, i, 1) : 1; +ck(ap+i, ((bn-cancel1)<an) ? (bn-cancel1) : an, 19); + mpn_sub_lshift_n(ap+i, bp, bn-cancel1, cancel2, an); +ck(ap+i, an-i, 20); + mpn_sub_1(ap+i, ap+i, an-i, 1-cc); + } + else { /* PREC(b) > PREC(a): we have to truncate b */ +ck(ap, an, 21); + mpn_sub_lshift_n(ap, bp+(bn-an-cancel1), an, cancel2, an); + } + /* remains to do the rounding */ +#ifdef DEBUG2 + 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 & ((1<<sh)-1); + *ap &= ~cc; /* truncate last bits */ + cc -= 1<<(sh-1); + while ((cc==0 || cc==-1) && k!=0 && kc!=0) { + kc--; + cc -= mpn_sub_1(&c2, bp+(--k), 1, (cp[kc]>>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) return; /* 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 & ((1<<sh)-1); + *ap &= ~cc; /* truncate last bits */ + if (cc) goto add_one_ulp; /* will happen most of the time */ + else { /* same four cases too */ + int kc = cn-k; /* remains kc limbs from c */ + switch (overlap) + { + case 1: /* both b and c to round */ + k = bn-an; /* remains k limbs from b */ + dif = diff_exp % mp_bits_per_limb; + while (cc==0 && k!=0 && kc!=0) { + kc--; + cc = bp[--k] - (cp[kc]>>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 & ((1<<sh)-1); + *ap &= ~cc; + if (cc==0) { /* otherwise neglected difference cannot be < 0 */ + /* take into account bp[0]..bp[bn-cancel1-1] shifted by cancel2 to left + and cp[0]..cp[cn-k-1] shifted by dif bits to right */ + switch (overlap) { + case 0: + case 2: + break; /* c is not truncated ==> 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 DEBUG2 +mpfr_print_raw(a); putchar('\n'); +#endif + if (sh) { + cc = *ap & ((1<<sh)-1); + *ap &= ~cc; /* truncate last bits */ + c2 = 1<<(sh-1); + } + else /* no bit to truncate */ + { cc = bp[--k]; c2 = 1<<(mp_bits_per_limb-1); } +#ifdef DEBUG2 + printf("cc=%lu c2=%lu k=%u\n",cc,c2,k); +#endif + if (cc>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 DEBUG2 + 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 & (1<<sh))) goto add_one_ulp; + } + else { + /* subtract: if cc>0, do nothing */ + if (cc==0 && (*ap & (1<<sh))) goto add_one_ulp; + } + } + goto end_of_sub; + + sub_one_ulp: + cc = mpn_sub_1(ap, ap, an, 1<<sh); + if (cc) { printf("carry(2) in mpfr_sub1\n"); exit(1); } + goto end_of_sub; + + add_one_ulp: /* add one unit in last place to a */ + cc = mpn_add_1(ap, ap, an, 1<<sh); + if (cc) { printf("carry(3) in mpfr_sub1\n"); exit(1); } + + end_of_sub: +#ifdef DEBUG2 +printf("b-c="); if (SIGN(a)>0) putchar(' '); mpfr_print_raw(a); putchar('\n'); +#endif + return; +} + |