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