![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to add1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr|
Return to add1.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpfr |
1.1 ohara 1: /* mpfr_add1 -- internal function to perform a "real" addition
2:
3: Copyright 1999, 2000, 2001 Free Software Foundation.
4: Contributed by the Spaces project, INRIA Lorraine.
5:
6: This file is part of the MPFR Library.
7:
8: The MPFR Library is free software; you can redistribute it and/or modify
9: it under the terms of the GNU Lesser General Public License as published by
10: the Free Software Foundation; either version 2.1 of the License, or (at your
11: option) any later version.
12:
13: The MPFR Library is distributed in the hope that it will be useful, but
14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16: License for more details.
17:
18: You should have received a copy of the GNU Lesser General Public License
19: along with the MPFR Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: #include "gmp.h"
24: #include "gmp-impl.h"
25: #include "mpfr.h"
26: #include "mpfr-impl.h"
27:
28: /* compute sign(b) * (|b| + |c|)
29: diff_exp is the difference between the exponents of b and c,
30: which is >= 0 */
31:
32: int
33: mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
34: mp_rnd_t rnd_mode, mp_exp_unsigned_t diff_exp)
35: {
36: mp_limb_t *ap, *bp, *cp;
37: mp_prec_t aq, bq, cq, aq2;
38: mp_size_t an, bn, cn;
39: mp_exp_t difw;
40: int sh, rb, fb, inex;
41: TMP_DECL(marker);
42:
43: MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_NOTZERO(b));
44: MPFR_ASSERTN(MPFR_IS_FP(c) && MPFR_NOTZERO(c));
45:
46: TMP_MARK(marker);
47: ap = MPFR_MANT(a);
48: bp = MPFR_MANT(b);
49: cp = MPFR_MANT(c);
50:
51: if (ap == bp)
52: {
53: bp = (mp_ptr) TMP_ALLOC((size_t) MPFR_ABSSIZE(b) * BYTES_PER_MP_LIMB);
54: MPN_COPY (bp, ap, MPFR_ABSSIZE(b));
55: if (ap == cp)
56: { cp = bp; }
57: }
58: else if (ap == cp)
59: {
60: cp = (mp_ptr) TMP_ALLOC ((size_t) MPFR_ABSSIZE(c) * BYTES_PER_MP_LIMB);
61: MPN_COPY(cp, ap, MPFR_ABSSIZE(c));
62: }
63:
64: aq = MPFR_PREC(a);
65: bq = MPFR_PREC(b);
66: cq = MPFR_PREC(c);
67:
68: an = (aq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of a */
69: aq2 = an * BITS_PER_MP_LIMB;
70: sh = aq2 - aq; /* non-significant bits in low limb */
71: bn = (bq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */
72: cn = (cq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of c */
73:
74: MPFR_EXP(a) = MPFR_EXP(b);
75: MPFR_SET_SAME_SIGN(a, b);
76:
77: /*
78: * 1. Compute the significant part A', the non-significant bits of A
79: * are taken into account.
80: *
81: * 2. Perform the rounding. At each iteration, we remember:
82: * _ r = rounding bit
83: * _ f = following bits (same value)
84: * where the result has the form: [number A]rfff...fff + a remaining
85: * value in the interval [0,2) ulp. We consider the most significant
86: * bits of the remaining value to update the result; a possible carry
87: * is immediately taken into account and A is updated accordingly. As
88: * soon as the bits f don't have the same value, A can be rounded.
89: * Variables:
90: * _ rb = rounding bit (0 or 1).
91: * _ fb = following bits (0 or 1), then sticky bit.
92: * If fb == 0, the only thing that can change is the sticky bit.
93: */
94:
95: rb = fb = -1; /* means: not initialized */
96:
97: if (aq2 <= diff_exp)
98: { /* c does not overlap with a' */
99: if (an > bn)
100: { /* a has more limbs than b */
101: /* copy b to the most significant limbs of a */
102: MPN_COPY(ap + (an - bn), bp, bn);
103: /* zero the least significant limbs of a */
104: MPN_ZERO(ap, an - bn);
105: }
106: else /* an <= bn */
107: {
108: /* copy the most significant limbs of b to a */
109: MPN_COPY(ap, bp + (bn - an), an);
110: }
111: }
112: else /* aq2 > diff_exp */
113: { /* c overlaps with a' */
114: mp_limb_t *a2p;
115: mp_limb_t cc;
116: mp_prec_t dif;
117: mp_size_t difn, k;
118: int shift;
119:
120: /* copy c (shifted) into a */
121:
122: dif = aq2 - diff_exp;
123: /* dif is the number of bits of c which overlap with a' */
124:
125: difn = (dif-1)/BITS_PER_MP_LIMB + 1;
126: /* only the highest difn limbs from c have to be considered */
127: if (difn > cn)
128: {
129: /* c doesn't have enough limbs; take into account the virtual
130: zero limbs now by zeroing the least significant limbs of a' */
131: MPFR_ASSERTN(difn - cn <= an);
132: MPN_ZERO(ap, difn - cn);
133: difn = cn;
134: }
135:
136: k = diff_exp / BITS_PER_MP_LIMB;
137:
138: /* zero the most significant k limbs of a */
139: a2p = ap + (an - k);
140: MPN_ZERO(a2p, k);
141:
142: shift = diff_exp % BITS_PER_MP_LIMB;
143:
144: if (shift)
145: {
146: MPFR_ASSERTN(a2p - difn >= ap);
147: cc = mpn_rshift(a2p - difn, cp + (cn - difn), difn, shift);
148: if (a2p - difn > ap)
149: *(a2p - difn - 1) = cc;
150: }
151: else
152: MPN_COPY(a2p - difn, cp + (cn - difn), difn);
153:
154: /* add b to a */
155:
156: cc = an > bn
157: ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn)
158: : mpn_add_n(ap, ap, bp + (bn - an), an);
159:
160: if (cc) /* carry */
161: {
162: mp_exp_t exp = MPFR_EXP(a);
163: if (exp == __mpfr_emax)
164: {
165: inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a));
166: goto end_of_add;
167: }
168: MPFR_EXP(a)++;
169: rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */
170: if (sh)
171: {
172: mp_limb_t mask, bb;
173:
174: mask = (MP_LIMB_T_ONE << sh) - 1;
175: bb = ap[0] & mask;
176: ap[0] &= (~mask) << 1;
177: if (bb == 0)
178: fb = 0;
179: else if (bb == mask)
180: fb = 1;
181: }
182: mpn_rshift(ap, ap, an, 1);
183: ap[an-1] += GMP_LIMB_HIGHBIT;
184: if (sh && fb < 0)
185: goto rounding;
186: } /* cc */
187: } /* aq2 > diff_exp */
188:
189: /* non-significant bits of a */
190: if (rb < 0 && sh)
191: {
192: mp_limb_t mask, bb;
193:
194: mask = (MP_LIMB_T_ONE << sh) - 1;
195: bb = ap[0] & mask;
196: ap[0] &= ~mask;
197: rb = bb >> (sh - 1);
198: if (sh > 1)
199: {
200: mask >>= 1;
201: bb &= mask;
202: if (bb == 0)
203: fb = 0;
204: else if (bb == mask)
205: fb = 1;
206: else
207: goto rounding;
208: }
209: }
210:
211: /* determine rounding and sticky bits (and possible carry) */
212:
213: difw = (mp_exp_t) an - (mp_exp_t) (diff_exp / BITS_PER_MP_LIMB);
214: /* difw is the number of limbs from b (regarded as having an infinite
215: precision) that have already been combined with c; -n if the next
216: n limbs from b won't be combined with c. */
217:
218: if (bn > an)
219: { /* there are still limbs from b that haven't been taken into account */
220: mp_size_t bk;
221:
222: if (fb == 0 && difw <= 0)
223: {
224: fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */
225: goto rounding;
226: }
227:
228: bk = bn - an; /* index of lowest considered limb from b, > 0 */
229: while (difw < 0)
230: { /* ulp(next limb from b) > msb(c) */
231: mp_limb_t bb;
232:
233: bb = bp[--bk];
234:
235: MPFR_ASSERTD(fb != 0);
236: if (fb > 0)
237: {
238: if (bb != MP_LIMB_T_MAX)
239: {
240: fb = 1; /* c hasn't been taken into account
241: ==> sticky bit != 0 */
242: goto rounding;
243: }
244: }
245: else /* fb not initialized yet */
246: {
247: if (rb < 0) /* rb not initialized yet */
248: {
249: rb = bb >> (BITS_PER_MP_LIMB - 1);
250: bb |= GMP_LIMB_HIGHBIT;
251: }
252: fb = 1;
253: if (bb != MP_LIMB_T_MAX)
254: goto rounding;
255: }
256:
257: if (bk == 0)
258: { /* b has entirely been read */
259: fb = 1; /* c hasn't been taken into account
260: ==> sticky bit != 0 */
261: goto rounding;
262: }
263:
264: difw++;
265: } /* while */
266: MPFR_ASSERTN(bk > 0 && difw >= 0);
267:
268: if (difw <= cn)
269: {
270: mp_size_t ck;
271: mp_limb_t cprev;
272: int difs;
273:
274: ck = cn - difw;
275: difs = diff_exp % BITS_PER_MP_LIMB;
276:
277: if (difs == 0 && ck == 0)
278: goto c_read;
279:
280: cprev = ck == cn ? 0 : cp[ck];
281:
282: if (fb < 0)
283: {
284: mp_limb_t bb, cc;
285:
286: if (difs)
287: {
288: cc = cprev << (BITS_PER_MP_LIMB - difs);
289: if (--ck >= 0)
290: {
291: cprev = cp[ck];
292: cc += cprev >> difs;
293: }
294: }
295: else
296: cc = cp[--ck];
297:
298: bb = bp[--bk] + cc;
299:
300: if (bb < cc /* carry */
301: && (rb < 0 || (rb ^= 1) == 0)
302: && mpn_add_1(ap, ap, an, MP_LIMB_T_ONE << sh))
303: {
304: mp_exp_t exp = MPFR_EXP(a);
305: if (exp == __mpfr_emax)
306: {
307: inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a));
308: goto end_of_add;
309: }
310: MPFR_EXP(a)++;
311: ap[an-1] = GMP_LIMB_HIGHBIT;
312: rb = 0;
313: }
314:
315: if (rb < 0) /* rb not initialized yet */
316: {
317: rb = bb >> (BITS_PER_MP_LIMB - 1);
318: bb <<= 1;
319: bb |= bb >> (BITS_PER_MP_LIMB - 1);
320: }
321:
322: fb = bb != 0;
323: if (fb && bb != MP_LIMB_T_MAX)
324: goto rounding;
325: } /* fb < 0 */
326:
327: while (bk > 0)
328: {
329: mp_limb_t bb, cc;
330:
331: if (difs)
332: {
333: if (ck < 0)
334: goto c_read;
335: cc = cprev << (BITS_PER_MP_LIMB - difs);
336: if (--ck >= 0)
337: {
338: cprev = cp[ck];
339: cc += cprev >> difs;
340: }
341: }
342: else
343: {
344: if (ck == 0)
345: goto c_read;
346: cc = cp[--ck];
347: }
348:
349: bb = bp[--bk] + cc;
350: if (bb < cc) /* carry */
351: {
352: fb ^= 1;
353: if (fb)
354: goto rounding;
355: rb ^= 1;
356: if (rb == 0 && mpn_add_1(ap, ap, an, MP_LIMB_T_ONE << sh))
357: {
358: mp_exp_t exp = MPFR_EXP(a);
359: if (exp == __mpfr_emax)
360: {
361: inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a));
362: goto end_of_add;
363: }
364: MPFR_EXP(a)++;
365: ap[an-1] = GMP_LIMB_HIGHBIT;
366: }
367: } /* bb < cc */
368:
369: if (!fb && bb != 0)
370: {
371: fb = 1;
372: goto rounding;
373: }
374: if (fb && bb != MP_LIMB_T_MAX)
375: goto rounding;
376: } /* while */
377:
378: /* b has entirely been read */
379:
380: if (fb || ck < 0)
381: goto rounding;
382: if (difs && cprev << (BITS_PER_MP_LIMB - difs))
383: {
384: fb = 1;
385: goto rounding;
386: }
387: while (ck)
388: {
389: if (cp[--ck])
390: {
391: fb = 1;
392: goto rounding;
393: }
394: } /* while */
395: } /* difw <= cn */
396: else
397: { /* c has entirely been read */
398: c_read:
399: if (fb < 0) /* fb not initialized yet */
400: {
401: mp_limb_t bb;
402:
403: MPFR_ASSERTN(bk > 0);
404: bb = bp[--bk];
405: if (rb < 0) /* rb not initialized yet */
406: {
407: rb = bb >> (BITS_PER_MP_LIMB - 1);
408: bb &= ~GMP_LIMB_HIGHBIT;
409: }
410: fb = bb != 0;
411: } /* fb < 0 */
412: if (fb)
413: goto rounding;
414: while (bk)
415: {
416: if (bp[--bk])
417: {
418: fb = 1;
419: goto rounding;
420: }
421: } /* while */
422: } /* difw > cn */
423: } /* bn > an */
424: else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
425: { /* b has entirely been read */
426: if (difw > cn)
427: { /* c has entirely been read */
428: if (rb < 0)
429: rb = 0;
430: fb = 0;
431: }
432: else if (diff_exp > aq2)
433: { /* b is followed by at least a zero bit, then by c */
434: if (rb < 0)
435: rb = 0;
436: fb = 1;
437: }
438: else
439: {
440: mp_size_t ck;
441: int difs;
442:
443: MPFR_ASSERTN(difw >= 0 && cn >= difw);
444: ck = cn - difw;
445: difs = diff_exp % BITS_PER_MP_LIMB;
446:
447: if (difs == 0 && ck == 0)
448: { /* c has entirely been read */
449: if (rb < 0)
450: rb = 0;
451: fb = 0;
452: }
453: else
454: {
455: mp_limb_t cc;
456:
457: cc = difs ? (MPFR_ASSERTN(ck < cn),
458: cp[ck] << (BITS_PER_MP_LIMB - difs)) : cp[--ck];
459: if (rb < 0)
460: {
461: rb = cc >> (BITS_PER_MP_LIMB - 1);
462: cc &= ~GMP_LIMB_HIGHBIT;
463: }
464: while (cc == 0)
465: {
466: if (ck == 0)
467: {
468: fb = 0;
469: goto rounding;
470: }
471: cc = cp[--ck];
472: } /* while */
473: fb = 1;
474: }
475: }
476: } /* fb != 1 */
477:
478: rounding:
479: switch (rnd_mode)
480: {
481: case GMP_RNDN:
482: if (fb == 0)
483: {
484: if (rb == 0)
485: {
486: inex = 0;
487: goto end_of_add;
488: }
489: /* round to even */
490: if (ap[0] & (MP_LIMB_T_ONE << sh))
491: goto rndn_away;
492: else
493: goto rndn_zero;
494: }
495: if (rb == 0)
496: {
497: rndn_zero:
498: inex = MPFR_ISNEG(a) ? 1 : -1;
499: goto end_of_add;
500: }
501: else
502: {
503: rndn_away:
504: inex = MPFR_ISNONNEG(a) ? 1 : -1;
505: goto add_one_ulp;
506: }
507:
508: case GMP_RNDZ:
509: inex = rb || fb ? (MPFR_ISNEG(a) ? 1 : -1) : 0;
510: goto end_of_add;
511:
512: case GMP_RNDU:
513: inex = rb || fb;
514: if (inex && MPFR_ISNONNEG(a))
515: goto add_one_ulp;
516: else
517: goto end_of_add;
518:
519: case GMP_RNDD:
520: inex = - (rb || fb);
521: if (inex && MPFR_ISNEG(a))
522: goto add_one_ulp;
523: else
524: goto end_of_add;
525:
526: default:
527: MPFR_ASSERTN(0);
528: inex = 0;
529: goto end_of_add;
530: }
531:
532: add_one_ulp: /* add one unit in last place to a */
533: if (mpn_add_1(ap, ap, an, MP_LIMB_T_ONE << sh))
534: {
535: mp_exp_t exp = MPFR_EXP(a);
536: if (exp == __mpfr_emax)
537: inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a));
538: else
539: {
540: MPFR_EXP(a)++;
541: ap[an-1] = GMP_LIMB_HIGHBIT;
542: }
543: }
544:
545: end_of_add:
546: TMP_FREE(marker);
547: return inex;
548: }