![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ui_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 ui_sub.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpf |
1.1 maekawa 1: /* mpf_ui_sub -- Subtract a float from an unsigned long int.
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_ui_sub (mpf_ptr r, unsigned long int u, mpf_srcptr v)
28: #else
29: mpf_ui_sub (r, u, v)
30: mpf_ptr r;
31: unsigned long int 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 uexp;
40: mp_size_t ediff;
41: int negate;
42: mp_limb_t ulimb;
43: TMP_DECL (marker);
44:
45: vsize = v->_mp_size;
46:
47: /* Handle special cases that don't work in generic code below. */
48: if (u == 0)
49: {
50: mpf_neg (r, v);
51: return;
52: }
53: if (vsize == 0)
54: {
55: mpf_set_ui (r, u);
56: return;
57: }
58:
59: /* If signs of U and V are different, perform addition. */
60: if (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_ui (r, &v_negated, u);
67: return;
68: }
69:
70: TMP_MARK (marker);
71:
72: /* Signs are now known to be the same. */
73:
74: ulimb = u;
75: /* Make U be the operand with the largest exponent. */
76: if (1 < v->_mp_exp)
77: {
78: negate = 1;
79: usize = ABS (vsize);
80: vsize = 1;
81: up = v->_mp_d;
82: vp = &ulimb;
83: rp = r->_mp_d;
84: prec = r->_mp_prec + 1;
85: uexp = v->_mp_exp;
86: ediff = uexp - 1;
87: }
88: else
89: {
90: negate = 0;
91: usize = 1;
92: vsize = ABS (vsize);
93: up = &ulimb;
94: vp = v->_mp_d;
95: rp = r->_mp_d;
96: prec = r->_mp_prec;
97: uexp = 1;
98: ediff = 1 - v->_mp_exp;
99: }
100:
101: /* Ignore leading limbs in U and V that are equal. Doing
102: this helps increase the precision of the result. */
103: if (ediff == 0)
104: {
105: /* This loop normally exits immediately. Optimize for that. */
106: for (;;)
107: {
108: usize--;
109: vsize--;
110: if (up[usize] != vp[vsize])
111: break;
112: uexp--;
113: if (usize == 0)
114: goto Lu0;
115: if (vsize == 0)
116: goto Lv0;
117: }
118: usize++;
119: vsize++;
120: /* Note that either operand (but not both operands) might now have
121: leading zero limbs. It matters only that U is unnormalized if
122: vsize is now zero, and vice versa. And it is only in that case
123: that we have to adjust uexp. */
124: if (vsize == 0)
125: Lv0:
126: while (usize != 0 && up[usize - 1] == 0)
127: usize--, uexp--;
128: if (usize == 0)
129: Lu0:
130: while (vsize != 0 && vp[vsize - 1] == 0)
131: vsize--, uexp--;
132: }
133:
134: /* If U extends beyond PREC, ignore the part that does. */
135: if (usize > prec)
136: {
137: up += usize - prec;
138: usize = prec;
139: }
140:
141: /* If V extends beyond PREC, ignore the part that does.
142: Note that this may make vsize negative. */
143: if (vsize + ediff > prec)
144: {
145: vp += vsize + ediff - prec;
146: vsize = prec - ediff;
147: }
148:
149: /* Allocate temp space for the result. Allocate
150: just vsize + ediff later??? */
151: tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
152:
153: if (ediff >= prec)
154: {
155: /* V completely cancelled. */
156: if (tp != up)
157: MPN_COPY (rp, up, usize);
158: rsize = usize;
159: }
160: else
161: {
162: /* Locate the least significant non-zero limb in (the needed
163: parts of) U and V, to simplify the code below. */
164: for (;;)
165: {
166: if (vsize == 0)
167: {
168: MPN_COPY (rp, up, usize);
169: rsize = usize;
170: goto done;
171: }
172: if (vp[0] != 0)
173: break;
174: vp++, vsize--;
175: }
176: for (;;)
177: {
178: if (usize == 0)
179: {
180: MPN_COPY (rp, vp, vsize);
181: rsize = vsize;
182: negate ^= 1;
183: goto done;
184: }
185: if (up[0] != 0)
186: break;
187: up++, usize--;
188: }
189:
190: /* uuuu | uuuu | uuuu | uuuu | uuuu */
191: /* vvvvvvv | vv | vvvvv | v | vv */
192:
193: if (usize > ediff)
194: {
195: /* U and V partially overlaps. */
196: if (ediff == 0)
197: {
198: /* Have to compare the leading limbs of u and v
199: to determine whether to compute u - v or v - u. */
200: if (usize > vsize)
201: {
202: /* uuuu */
203: /* vv */
204: int cmp;
205: cmp = mpn_cmp (up + usize - vsize, vp, vsize);
206: if (cmp >= 0)
207: {
208: mp_size_t size;
209: size = usize - vsize;
210: MPN_COPY (tp, up, size);
211: mpn_sub_n (tp + size, up + size, vp, vsize);
212: rsize = usize;
213: }
214: else
215: {
216: /* vv */ /* Swap U and V. */
217: /* uuuu */
218: mp_size_t size, i;
219: size = usize - vsize;
220: tp[0] = -up[0];
221: for (i = 1; i < size; i++)
222: tp[i] = ~up[i];
223: mpn_sub_n (tp + size, vp, up + size, vsize);
224: mpn_sub_1 (tp + size, tp + size, vsize, (mp_limb_t) 1);
225: negate ^= 1;
226: rsize = usize;
227: }
228: }
229: else if (usize < vsize)
230: {
231: /* uuuu */
232: /* vvvvvvv */
233: int cmp;
234: cmp = mpn_cmp (up, vp + vsize - usize, usize);
235: if (cmp > 0)
236: {
237: mp_size_t size, i;
238: size = vsize - usize;
239: tp[0] = -vp[0];
240: for (i = 1; i < size; i++)
241: tp[i] = ~vp[i];
242: mpn_sub_n (tp + size, up, vp + size, usize);
243: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
244: rsize = vsize;
245: }
246: else
247: {
248: /* vvvvvvv */ /* Swap U and V. */
249: /* uuuu */
250: /* This is the only place we can get 0.0. */
251: mp_size_t size;
252: size = vsize - usize;
253: MPN_COPY (tp, vp, size);
254: mpn_sub_n (tp + size, vp + size, up, usize);
255: negate ^= 1;
256: rsize = vsize;
257: }
258: }
259: else
260: {
261: /* uuuu */
262: /* vvvv */
263: int cmp;
264: cmp = mpn_cmp (up, vp + vsize - usize, usize);
265: if (cmp > 0)
266: {
267: mpn_sub_n (tp, up, vp, usize);
268: rsize = usize;
269: }
270: else
271: {
272: mpn_sub_n (tp, vp, up, usize);
273: negate ^= 1;
274: rsize = usize;
275: /* can give zero */
276: }
277: }
278: }
279: else
280: {
281: if (vsize + ediff <= usize)
282: {
283: /* uuuu */
284: /* v */
285: mp_size_t size;
286: size = usize - ediff - vsize;
287: MPN_COPY (tp, up, size);
288: mpn_sub (tp + size, up + size, usize - size, vp, vsize);
289: rsize = usize;
290: }
291: else
292: {
293: /* uuuu */
294: /* vvvvv */
295: mp_size_t size, i;
296: size = vsize + ediff - usize;
297: tp[0] = -vp[0];
298: for (i = 1; i < size; i++)
299: tp[i] = ~vp[i];
300: mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
301: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
302: rsize = vsize + ediff;
303: }
304: }
305: }
306: else
307: {
308: /* uuuu */
309: /* vv */
310: mp_size_t size, i;
311: size = vsize + ediff - usize;
312: tp[0] = -vp[0];
313: for (i = 1; i < vsize; i++)
314: tp[i] = ~vp[i];
315: for (i = vsize; i < size; i++)
316: tp[i] = ~(mp_limb_t) 0;
317: mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
318: rsize = size + usize;
319: }
320:
321: /* Full normalize. Optimize later. */
322: while (rsize != 0 && tp[rsize - 1] == 0)
323: {
324: rsize--;
325: uexp--;
326: }
327: MPN_COPY (rp, tp, rsize);
328: }
329:
330: done:
331: r->_mp_size = negate ? -rsize : rsize;
332: r->_mp_exp = uexp;
333: TMP_FREE (marker);
334: }