![[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 / mpf|
Return to sub.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpf |
1.1 maekawa 1: /* mpf_sub -- Subtract two floats.
2:
3: Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
4:
5: This file is part of the GNU MP Library.
6:
7: The GNU MP 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 GNU MP 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 GNU MP 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 "gmp.h"
23: #include "gmp-impl.h"
24:
25: void
26: #if __STDC__
27: mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
28: #else
29: mpf_sub (r, u, v)
30: mpf_ptr r;
31: mpf_srcptr u;
32: mpf_srcptr v;
33: #endif
34: {
35: mp_srcptr up, vp;
36: mp_ptr rp, tp;
37: mp_size_t usize, vsize, rsize;
38: mp_size_t prec;
39: mp_exp_t exp;
40: mp_size_t ediff;
41: int negate;
42: TMP_DECL (marker);
43:
44: usize = u->_mp_size;
45: vsize = v->_mp_size;
46:
47: /* Handle special cases that don't work in generic code below. */
48: if (usize == 0)
49: {
50: mpf_neg (r, v);
51: return;
52: }
53: if (vsize == 0)
54: {
55: mpf_set (r, u);
56: return;
57: }
58:
59: /* If signs of U and V are different, perform addition. */
60: if ((usize ^ vsize) < 0)
61: {
62: __mpf_struct v_negated;
63: v_negated._mp_size = -vsize;
64: v_negated._mp_exp = v->_mp_exp;
65: v_negated._mp_d = v->_mp_d;
66: mpf_add (r, u, &v_negated);
67: return;
68: }
69:
70: TMP_MARK (marker);
71:
72: /* Signs are now known to be the same. */
73: negate = usize < 0;
74:
75: /* Make U be the operand with the largest exponent. */
76: if (u->_mp_exp < v->_mp_exp)
77: {
78: mpf_srcptr t;
79: t = u; u = v; v = t;
80: negate ^= 1;
81: usize = u->_mp_size;
82: vsize = v->_mp_size;
83: }
84:
85: usize = ABS (usize);
86: vsize = ABS (vsize);
87: up = u->_mp_d;
88: vp = v->_mp_d;
89: rp = r->_mp_d;
90: prec = r->_mp_prec + 1;
91: exp = u->_mp_exp;
92: ediff = u->_mp_exp - v->_mp_exp;
93:
94: /* If ediff is 0 or 1, we might have a situation where the operands are
95: extremely close. We need to scan the operands from the most significant
96: end ignore the initial parts that are equal. */
97: if (ediff <= 1)
98: {
99: if (ediff == 0)
100: {
101: /* Skip leading limbs in U and V that are equal. */
102: if (up[usize - 1] == vp[vsize - 1])
103: {
104: /* This loop normally exits immediately. Optimize for that. */
105: do
106: {
107: usize--;
108: vsize--;
109: exp--;
110:
111: if (usize == 0)
112: {
113: rsize = vsize;
114: tp = (mp_ptr) vp;
115: negate ^= 1;
116: goto normalize;
117: }
118: if (vsize == 0)
119: {
120: rsize = usize;
121: tp = (mp_ptr) up;
122: goto normalize;
123: }
124: }
125: while (up[usize - 1] == vp[vsize - 1]);
126: }
127:
128: if (up[usize - 1] < vp[vsize - 1])
129: {
130: /* For simplicity, swap U and V. Note that since the loop above
131: wouldn't have exited unless up[usize - 1] and vp[vsize - 1]
132: were non-equal, this if-statement catches all cases where U
133: is smaller than V. */
134: { mp_srcptr tp = up; up = vp; vp = tp; }
135: { mp_size_t tsize = usize; usize = vsize; vsize = tsize; }
136: negate ^= 1;
137: /* negating ediff not necessary since it is 0. */
138: }
139:
140: /* Check for
141: x+1 00000000 ...
142: x ffffffff ... */
143: if (up[usize - 1] != vp[vsize - 1] + 1)
144: goto general_case;
145: usize--;
146: vsize--;
147: exp--;
148: }
149: else /* ediff == 1 */
150: {
151: /* Check for
152: 1 00000000 ...
153: 0 ffffffff ... */
154:
155: if (up[usize - 1] != 1 || vp[vsize - 1] != ~(mp_limb_t) 0
156: || (usize >= 2 && up[usize - 2] != 0))
157: goto general_case;
158:
159: usize--;
160: exp--;
161: }
162:
163: /* Skip sequences of 00000000/ffffffff */
164: while (vsize != 0 && usize != 0 && up[usize - 1] == 0
165: && vp[vsize - 1] == ~(mp_limb_t) 0)
166: {
167: usize--;
168: vsize--;
169: exp--;
170: }
171:
172: if (usize == 0)
173: {
174: while (vsize != 0 && vp[vsize - 1] == ~(mp_limb_t) 0)
175: {
176: vsize--;
177: exp--;
178: }
179: }
180:
181: if (usize > prec - 1)
182: {
183: up += usize - (prec - 1);
184: usize = prec - 1;
185: }
186: if (vsize > prec - 1)
187: {
188: vp += vsize - (prec - 1);
189: vsize = prec - 1;
190: }
191:
192: tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
193: {
194: mp_limb_t cy_limb;
195: if (vsize == 0)
196: {
197: mp_size_t size, i;
198: size = usize;
199: for (i = 0; i < size; i++)
200: tp[i] = up[i];
201: tp[size] = 1;
202: rsize = size + 1;
203: exp++;
204: goto normalize;
205: }
206: if (usize == 0)
207: {
208: mp_size_t size, i;
209: size = vsize;
210: for (i = 0; i < size; i++)
211: tp[i] = ~vp[i];
212: cy_limb = 1 - mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
213: rsize = vsize;
214: if (cy_limb == 0)
215: {
216: tp[rsize] = 1;
217: rsize++;
218: exp++;
219: }
220: goto normalize;
221: }
222: if (usize >= vsize)
223: {
224: /* uuuu */
225: /* vv */
226: mp_size_t size;
227: size = usize - vsize;
228: MPN_COPY (tp, up, size);
229: cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize);
230: rsize = usize;
231: }
232: else /* (usize < vsize) */
233: {
234: /* uuuu */
235: /* vvvvvvv */
236: mp_size_t size, i;
237: size = vsize - usize;
238: for (i = 0; i < size; i++)
239: tp[i] = ~vp[i];
240: cy_limb = mpn_sub_n (tp + size, up, vp + size, usize);
241: cy_limb+= mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
242: cy_limb-= mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
243: rsize = vsize;
244: }
245: if (cy_limb == 0)
246: {
247: tp[rsize] = 1;
248: rsize++;
249: exp++;
250: }
251: goto normalize;
252: }
253: }
254:
255: general_case:
256: /* If U extends beyond PREC, ignore the part that does. */
257: if (usize > prec)
258: {
259: up += usize - prec;
260: usize = prec;
261: }
262:
263: /* If V extends beyond PREC, ignore the part that does.
264: Note that this may make vsize negative. */
265: if (vsize + ediff > prec)
266: {
267: vp += vsize + ediff - prec;
268: vsize = prec - ediff;
269: }
270:
271: /* Allocate temp space for the result. Allocate
272: just vsize + ediff later??? */
273: tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
274:
275: if (ediff >= prec)
276: {
277: /* V completely cancelled. */
278: if (tp != up)
279: MPN_COPY (rp, up, usize);
280: rsize = usize;
281: }
282: else
283: {
284: /* Locate the least significant non-zero limb in (the needed
285: parts of) U and V, to simplify the code below. */
286: for (;;)
287: {
288: if (vsize == 0)
289: {
290: MPN_COPY (rp, up, usize);
291: rsize = usize;
292: goto done;
293: }
294: if (vp[0] != 0)
295: break;
296: vp++, vsize--;
297: }
298: for (;;)
299: {
300: if (usize == 0)
301: {
302: MPN_COPY (rp, vp, vsize);
303: rsize = vsize;
304: negate ^= 1;
305: goto done;
306: }
307: if (up[0] != 0)
308: break;
309: up++, usize--;
310: }
311:
312: /* uuuu | uuuu | uuuu | uuuu | uuuu */
313: /* vvvvvvv | vv | vvvvv | v | vv */
314:
315: if (usize > ediff)
316: {
317: /* U and V partially overlaps. */
318: if (ediff == 0)
319: {
320: /* Have to compare the leading limbs of u and v
321: to determine whether to compute u - v or v - u. */
322: if (usize >= vsize)
323: {
324: /* uuuu */
325: /* vv */
326: mp_size_t size;
327: size = usize - vsize;
328: MPN_COPY (tp, up, size);
329: mpn_sub_n (tp + size, up + size, vp, vsize);
330: rsize = usize;
331: }
332: else /* (usize < vsize) */
333: {
334: /* uuuu */
335: /* vvvvvvv */
336: mp_size_t size, i;
337: size = vsize - usize;
338: tp[0] = -vp[0];
339: for (i = 1; i < size; i++)
340: tp[i] = ~vp[i];
341: mpn_sub_n (tp + size, up, vp + size, usize);
342: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
343: rsize = vsize;
344: }
345: }
346: else
347: {
348: if (vsize + ediff <= usize)
349: {
350: /* uuuu */
351: /* v */
352: mp_size_t size;
353: size = usize - ediff - vsize;
354: MPN_COPY (tp, up, size);
355: mpn_sub (tp + size, up + size, usize - size, vp, vsize);
356: rsize = usize;
357: }
358: else
359: {
360: /* uuuu */
361: /* vvvvv */
362: mp_size_t size, i;
363: size = vsize + ediff - usize;
364: tp[0] = -vp[0];
365: for (i = 1; i < size; i++)
366: tp[i] = ~vp[i];
367: mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
368: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
369: rsize = vsize + ediff;
370: }
371: }
372: }
373: else
374: {
375: /* uuuu */
376: /* vv */
377: mp_size_t size, i;
378: size = vsize + ediff - usize;
379: tp[0] = -vp[0];
380: for (i = 1; i < vsize; i++)
381: tp[i] = ~vp[i];
382: for (i = vsize; i < size; i++)
383: tp[i] = ~(mp_limb_t) 0;
384: mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
385: rsize = size + usize;
386: }
387:
388: normalize:
389: /* Full normalize. Optimize later. */
390: while (rsize != 0 && tp[rsize - 1] == 0)
391: {
392: rsize--;
393: exp--;
394: }
395: MPN_COPY (rp, tp, rsize);
396: }
397:
398: done:
399: r->_mp_size = negate ? -rsize : rsize;
400: r->_mp_exp = exp;
401: TMP_FREE (marker);
402: }