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