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