Annotation of OpenXM_contrib/gmp/mpfr/sub.c, Revision 1.1.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