Annotation of OpenXM_contrib/gmp/mpfr/sub1.c, Revision 1.1.1.1
1.1 ohara 1: /* mpfr_sub1 -- internal function to perform a "real" subtraction
2:
3: Copyright 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|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
29: Returns 0 iff result is exact,
30: a negative value when the result is less than the exact value,
31: a positive value otherwise.
32: */
33:
34: int
35: mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode,
36: int sub)
37: {
38: int sign;
39: mp_exp_unsigned_t diff_exp;
40: mp_prec_t cancel, cancel1;
41: mp_size_t cancel2, an, bn, cn, cn0;
42: mp_limb_t *ap, *bp, *cp;
43: mp_limb_t carry, bb, cc, borrow = 0;
44: int inexact = 0, shift_b, shift_c, is_exact = 1, down = 0, add_exp = 0;
45: int sh, k;
46: TMP_DECL(marker);
47:
48: TMP_MARK(marker);
49: ap = MPFR_MANT(a);
50: an = 1 + (MPFR_PREC(a) - 1) / BITS_PER_MP_LIMB;
51:
52: sign = mpfr_cmp2 (b, c, &cancel);
53: if (sign == 0)
54: {
55: if (rnd_mode == GMP_RNDD)
56: MPFR_SET_NEG(a);
57: else
58: MPFR_SET_POS(a);
59: MPFR_SET_ZERO(a);
60: MPFR_RET(0);
61: }
62:
63: /* If subtraction: sign(a) = sign * sign(b) */
64: if (sub && MPFR_SIGN(a) != sign * MPFR_SIGN(b))
65: MPFR_CHANGE_SIGN(a);
66:
67: if (sign < 0) /* swap b and c so that |b| > |c| */
68: {
69: mpfr_srcptr t;
70: t = b; b = c; c = t;
71: }
72:
73: /* If addition: sign(a) = sign of the larger argument in absolute value */
74: if (!sub)
75: MPFR_SET_SAME_SIGN(a, b);
76:
77: diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
78:
79: /* reserve a space to store b aligned with the result, i.e. shifted by
80: (-cancel) % BITS_PER_MP_LIMB to the right */
81: bn = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
82: shift_b = cancel % BITS_PER_MP_LIMB;
83: if (shift_b)
84: shift_b = BITS_PER_MP_LIMB - shift_b;
85: cancel1 = (cancel + shift_b) / BITS_PER_MP_LIMB;
86: /* the high cancel1 limbs from b should not be taken into account */
87: if (shift_b == 0)
88: bp = MPFR_MANT(b); /* no need of an extra space */
89: else
90: {
91: bp = TMP_ALLOC ((bn + 1) * BYTES_PER_MP_LIMB);
92: bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
93: }
94:
95: /* reserve a space to store c aligned with the result, i.e. shifted by
96: (diff_exp-cancel) % BITS_PER_MP_LIMB to the right */
97: cn = 1 + (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
98: shift_c = diff_exp - (cancel % BITS_PER_MP_LIMB);
99: shift_c = (shift_c + BITS_PER_MP_LIMB) % BITS_PER_MP_LIMB;
100: if (shift_c == 0)
101: cp = MPFR_MANT(c);
102: else
103: {
104: cp = TMP_ALLOC ((cn + 1) * BYTES_PER_MP_LIMB);
105: cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
106: }
107:
108: #ifdef DEBUG
109: printf("shift_b=%u shift_c=%u\n", shift_b, shift_c);
110: #endif
111:
112: /* ensure ap != bp and ap != cp */
113: if (ap == bp)
114: {
115: bp = (mp_ptr) TMP_ALLOC(bn * BYTES_PER_MP_LIMB);
116: MPN_COPY (bp, ap, bn);
117: /* ap == cp cannot occur since we would have b=c, which is detected
118: in mpfr_add or mpfr_sub */
119: }
120: else if (ap == cp)
121: {
122: cp = (mp_ptr) TMP_ALLOC (cn * BYTES_PER_MP_LIMB);
123: MPN_COPY(cp, ap, cn);
124: }
125:
126: /* here we have shift_c = (diff_exp - cancel) % BITS_PER_MP_LIMB,
127: thus we want cancel2 = ceil((cancel - diff_exp) / BITS_PER_MP_LIMB) */
128:
129: cancel2 = (long int) (cancel - (diff_exp - shift_c)) / BITS_PER_MP_LIMB;
130: /* the high cancel2 limbs from b should not be taken into account */
131: #ifdef DEBUG
132: printf("cancel=%u cancel1=%u cancel2=%d\n", cancel, cancel1, cancel2);
133: #endif
134:
135: /* ap[an-1] ap[0]
136: <----------------+-----------|---->
137: <----------PREC(a)----------><-sh->
138: cancel1
139: limbs bp[bn-cancel1-1]
140: <--...-----><----------------+-----------+----------->
141: cancel2
142: limbs cp[cn-cancel2-1] cancel2 >= 0
143: <--...--><----------------+----------------+---------------->
144: (-cancel2) cancel2 < 0
145: limbs <----------------+---------------->
146: */
147:
148: /* first part: put in ap[0..an-1] the value of high(b) - high(c),
149: where high(b) consists of the high an+cancel1 limbs of b,
150: and high(c) consists of the high an+cancel2 limbs of c.
151: */
152:
153: /* copy high(b) into a */
154: if (an + cancel1 <= bn) /* a: <----------------+-----------|---->
155: b: <-----------------------------------------> */
156: MPN_COPY (ap, bp + bn - (an + cancel1), an);
157: else /* a: <----------------+-----------|---->
158: b: <-------------------------> */
159: if (cancel1 < bn) /* otherwise b does not overlap with a */
160: {
161: MPN_ZERO (ap, an + cancel1 - bn);
162: MPN_COPY (ap + an + cancel1 - bn, bp, bn - cancel1);
163: }
164: else
165: MPN_ZERO (ap, an);
166:
167: #ifdef DEBUG
168: printf("after copying high(b), a="); mpfr_print_binary(a); putchar('\n');
169: #endif
170:
171: /* subtract high(c) */
172: if (an + cancel2 > 0) /* otherwise c does not overlap with a */
173: {
174: mp_limb_t *ap2;
175:
176: if (cancel2 >= 0)
177: {
178: if (an + cancel2 <= cn) /* a: <----------------------------->
179: c: <-----------------------------------------> */
180: mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
181: else /* a: <---------------------------->
182: c: <-------------------------> */
183: {
184: ap2 = ap + an + cancel2 - cn;
185: if (cn > cancel2)
186: mpn_sub_n (ap2, ap2, cp, cn - cancel2);
187: }
188: }
189: else /* cancel2 < 0 */
190: {
191: if (an + cancel2 <= cn) /* a: <----------------------------->
192: c: <-----------------------------> */
193: borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an + cancel2);
194: else /* a: <---------------------------->
195: c: <----------------> */
196: {
197: ap2 = ap + an + cancel2 - cn;
198: borrow = mpn_sub_n (ap2, ap2, cp, cn);
199: }
200: ap2 = ap + an + cancel2;
201: mpn_sub_1 (ap2, ap2, -cancel2, borrow);
202: }
203: }
204:
205: #ifdef DEBUG
206: printf("after subtracting high(c), a="); mpfr_print_binary(a); putchar('\n');
207: #endif
208:
209: /* now perform rounding */
210: sh = an * BITS_PER_MP_LIMB - MPFR_PREC(a); /* last unused bits from a */
211: carry = ap[0] & ((MP_LIMB_T_ONE << sh) - MP_LIMB_T_ONE);
212: ap[0] -= carry;
213:
214: if (rnd_mode == GMP_RNDN)
215: {
216: if (sh)
217: {
218: is_exact = (carry == 0);
219: /* can decide except when carry = 2^(sh-1) [middle]
220: or carry = 0 [truncate, but cannot decide inexact flag] */
221: down = (carry < (MP_LIMB_T_ONE << (sh - 1)));
222: if (carry > (MP_LIMB_T_ONE << (sh - 1)))
223: goto add_one_ulp;
224: else if ((0 < carry) && down)
225: {
226: inexact = -1; /* result if smaller than exact value */
227: goto truncate;
228: }
229: }
230: }
231: else /* directed rounding: set rnd_mode to RNDZ iff towards zero */
232: {
233: if (((rnd_mode == GMP_RNDD) && (MPFR_SIGN(a) > 0)) ||
234: ((rnd_mode == GMP_RNDU) && (MPFR_SIGN(a) < 0)))
235: rnd_mode = GMP_RNDZ;
236:
237: if (carry)
238: {
239: if (rnd_mode == GMP_RNDZ)
240: {
241: inexact = -1;
242: goto truncate;
243: }
244: else /* round away */
245: goto add_one_ulp;
246: }
247: }
248:
249: /* we have to consider the low (bn - (an+cancel1)) limbs from b,
250: and the (cn - (an+cancel2)) limbs from c. */
251: bn -= an + cancel1;
252: cn0 = cn;
253: cn -= (long int) an + cancel2;
254: #ifdef DEBUG
255: printf("last %d bits from a are %lu, bn=%ld, cn=%ld\n", sh, carry, bn, cn);
256: #endif
257:
258: for (k = 0; (bn > 0) || (cn > 0); k = 1)
259: {
260: bb = (bn > 0) ? bp[--bn] : 0;
261: if ((cn > 0) && (cn-- <= cn0))
262: cc = cp[cn];
263: else
264: cc = 0;
265:
266: if (down == 0)
267: down = (bb < cc);
268:
269: if ((rnd_mode == GMP_RNDN) && !k && sh == 0)
270: {
271: mp_limb_t half = GMP_LIMB_HIGHBIT;
272:
273: is_exact = (bb == cc);
274:
275: /* add one ulp if bb > cc + half
276: truncate if cc - half < bb < cc + half
277: sub one ulp if bb < cc - half
278: */
279:
280: if (down)
281: {
282: if (cc >= half)
283: cc -= half;
284: else
285: bb += half;
286: }
287: else /* bb >= cc */
288: {
289: if (cc < half)
290: cc += half;
291: else
292: bb -= half;
293: }
294: }
295:
296: #ifdef DEBUG
297: printf(" bb=%lu cc=%lu down=%d is_exact=%d\n", bb, cc, down, is_exact);
298: #endif
299: if (bb < cc)
300: {
301: if (rnd_mode == GMP_RNDZ)
302: goto sub_one_ulp;
303: else if (rnd_mode != GMP_RNDN) /* round away */
304: {
305: inexact = 1;
306: goto truncate;
307: }
308: else /* round to nearest */
309: {
310: if (is_exact && sh == 0)
311: {
312: inexact = 0;
313: goto truncate;
314: }
315: else if (down && sh == 0)
316: goto sub_one_ulp;
317: else
318: {
319: inexact = (is_exact) ? 1 : -1;
320: goto truncate;
321: }
322: }
323: }
324: else if (bb > cc)
325: {
326: if (rnd_mode == GMP_RNDZ)
327: {
328: inexact = -1;
329: goto truncate;
330: }
331: else if (rnd_mode != GMP_RNDN) /* round away */
332: goto add_one_ulp;
333: else /* round to nearest */
334: {
335: if (is_exact)
336: {
337: inexact = -1;
338: goto truncate;
339: }
340: else if (down)
341: {
342: inexact = 1;
343: goto truncate;
344: }
345: else
346: goto add_one_ulp;
347: }
348: }
349: }
350:
351: if ((rnd_mode == GMP_RNDN) && !is_exact)
352: {
353: /* even rounding rule */
354: if ((ap[0] >> sh) & 1)
355: {
356: if (down) goto sub_one_ulp;
357: else goto add_one_ulp;
358: }
359: else
360: inexact = (down) ? 1 : -1;
361: }
362: else
363: inexact = 0;
364: goto truncate;
365:
366: sub_one_ulp: /* add one unit in last place to a */
367: mpn_sub_1 (ap, ap, an, MP_LIMB_T_ONE << sh);
368: inexact = -1;
369: goto end_of_sub;
370:
371: add_one_ulp: /* add one unit in last place to a */
372: if (mpn_add_1 (ap, ap, an, MP_LIMB_T_ONE << sh)) /* result is a power of 2 */
373: {
374: ap[an-1] = GMP_LIMB_HIGHBIT;
375: add_exp = 1;
376: }
377: inexact = 1; /* result larger than exact value */
378:
379: truncate:
380: if ((ap[an-1] >> (BITS_PER_MP_LIMB - 1)) == 0) /* case 1 - epsilon */
381: {
382: ap[an-1] = GMP_LIMB_HIGHBIT;
383: add_exp = 1;
384: }
385:
386: end_of_sub:
387: /* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
388: care of underflows/overflows in that computation, and of the allowed
389: exponent range */
390: if (cancel)
391: {
392: mp_exp_t exp_b;
393:
394: cancel -= add_exp; /* still valid as unsigned long */
395: exp_b = MPFR_EXP(b); /* save it in case a equals b */
396: MPFR_EXP(a) = MPFR_EXP(b) - cancel;
397: if ((MPFR_EXP(a) > exp_b) /* underflow in type mp_exp_t */
398: || (MPFR_EXP(a) < __mpfr_emin))
399: {
400: TMP_FREE(marker);
401: return mpfr_set_underflow (a, rnd_mode, MPFR_SIGN(a));
402: }
403: }
404: else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
405: {
406: /* in case cancel = 0, add_exp can still be 1, in case b is just
407: below a power of two, c is very small, prec(a) < prec(b),
408: and rnd=away or nearest */
409: if (add_exp && MPFR_EXP(b) == __mpfr_emax)
410: {
411: TMP_FREE(marker);
412: return mpfr_set_overflow (a, rnd_mode, MPFR_SIGN(a));
413: }
414: MPFR_EXP(a) = MPFR_EXP(b) + add_exp;
415: }
416: TMP_FREE(marker);
417: #ifdef DEBUG
418: printf ("result is a="); mpfr_print_binary(a); putchar('\n');
419: #endif
420: /* check that result is msb-normalized */
421: MPFR_ASSERTN(ap[an-1] > ~ap[an-1]);
422: return inexact * MPFR_SIGN(a);
423: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sub1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr