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