Branch data Line data Source code
1 : : #line 2 "../src/kernel/none/add.c"
2 : : /* Copyright (C) 2002-2003 The PARI group.
3 : :
4 : : This file is part of the PARI/GP package.
5 : :
6 : : PARI/GP is free software; you can redistribute it and/or modify it under the
7 : : terms of the GNU General Public License as published by the Free Software
8 : : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
9 : : ANY WARRANTY WHATSOEVER.
10 : :
11 : : Check the License for details. You should have received a copy of it, along
12 : : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 : :
15 : : INLINE GEN
16 : 1093505999 : icopy_sign(GEN x, long sx)
17 : : {
18 : 1093505999 : GEN y=icopy(x);
19 : 1094238649 : setsigne(y,sx);
20 : 1094238649 : return y;
21 : : }
22 : :
23 : : GEN
24 : 81440981 : addsi_sign(long x, GEN y, long sy)
25 : : {
26 : : long sx,ly;
27 : : GEN z;
28 : :
29 [ + + ]: 81440981 : if (!x) return icopy_sign(y, sy);
30 [ + + ]: 80724428 : if (!sy) return stoi(x);
31 [ + + ]: 77275178 : if (x<0) { sx=-1; x=-x; } else sx=1;
32 [ + + ]: 77275178 : if (sx==sy)
33 : : {
34 : 53975059 : z = adduispec(x,y+2, lgefint(y)-2);
35 : 53975059 : setsigne(z,sy); return z;
36 : : }
37 : 23300119 : ly=lgefint(y);
38 [ + + ]: 23300119 : if (ly==3)
39 : : {
40 : 3457826 : const long d = y[2] - x;
41 [ + + ]: 3457826 : if (!d) return gen_0;
42 : 3148079 : z=cgeti(3);
43 [ + + ][ + + ]: 3148079 : if (y[2] < 0 || d > 0) {
44 : 2709105 : z[1] = evalsigne(sy) | evallgefint(3);
45 : 2709105 : z[2] = d;
46 : : }
47 : : else {
48 : 438974 : z[1] = evalsigne(-sy) | evallgefint(3);
49 : 438974 : z[2] =-d;
50 : : }
51 : 3148079 : return z;
52 : : }
53 : 19842293 : z = subiuspec(y+2,x, ly-2);
54 : 81439186 : setsigne(z,sy); return z;
55 : : }
56 : : GEN
57 : 85873289 : addui_sign(ulong x, GEN y, long sy)
58 : : {
59 : : long ly;
60 : : GEN z;
61 : :
62 [ + + ]: 85873289 : if (!x) return icopy_sign(y, sy);
63 [ + + ]: 84747388 : if (!sy) return utoipos(x);
64 [ + + ]: 41819755 : if (sy == 1) return adduispec(x,y+2, lgefint(y)-2);
65 : 24031699 : ly=lgefint(y);
66 [ + + ]: 24031699 : if (ly==3)
67 : : {
68 : 24022341 : const ulong t = y[2];
69 [ + + ]: 24022341 : if (x == t) return gen_0;
70 : 3392935 : z=cgeti(3);
71 [ + + ]: 3392935 : if (x < t) {
72 : 3392627 : z[1] = evalsigne(-1) | evallgefint(3);
73 : 3392627 : z[2] = t - x;
74 : : }
75 : : else {
76 : 308 : z[1] = evalsigne(1) | evallgefint(3);
77 : 308 : z[2] = x - t;
78 : : }
79 : 3392935 : return z;
80 : : }
81 : 9358 : z = subiuspec(y+2,x, ly-2);
82 : 85873289 : setsigne(z,-1); return z;
83 : : }
84 : :
85 : : /* return gen_0 when the sign is 0 */
86 : : GEN
87 : 6554822531 : addii_sign(GEN x, long sx, GEN y, long sy)
88 : : {
89 : : long lx,ly;
90 : : GEN z;
91 : :
92 [ + + ][ + + ]: 6554822531 : if (!sx) return sy? icopy_sign(y, sy): gen_0;
93 [ + + ]: 5776035998 : if (!sy) return icopy_sign(x, sx);
94 : 4953862307 : lx = lgefint(x);
95 : 4953862307 : ly = lgefint(y);
96 [ + + ]: 4953862307 : if (sx==sy)
97 : 2362533165 : z = addiispec(x+2,y+2,lx-2,ly-2);
98 : : else
99 : : { /* sx != sy */
100 : 2591329142 : long i = cmpiispec(x+2,y+2,lx-2,ly-2);
101 [ + + ]: 2592326310 : if (!i) return gen_0;
102 : : /* we must ensure |x| > |y| for subiispec */
103 [ + + ]: 2499757185 : if (i < 0) {
104 : 960819266 : sx = sy;
105 : 960819266 : z = subiispec(y+2,x+2,ly-2,lx-2);
106 : : }
107 : : else
108 : 1538937919 : z = subiispec(x+2,y+2,lx-2,ly-2);
109 : : }
110 : 6552394531 : setsigne(z,sx); return z;
111 : : }
112 : :
113 : : INLINE GEN
114 : 63800461 : rcopy_sign(GEN x, long sx) { GEN y = rcopy(x); setsigne(y,sx); return y; }
115 : :
116 : : GEN
117 : 192393888 : addir_sign(GEN x, long sx, GEN y, long sy)
118 : : {
119 : : long e, l, ly;
120 : : GEN z;
121 : :
122 [ + + ]: 192393888 : if (!sx) return rcopy_sign(y, sy);
123 : 141094730 : e = expo(y) - expi(x);
124 [ + + ]: 141094730 : if (!sy)
125 : : {
126 [ + + ]: 6360005 : if (e >= 0) return rcopy_sign(y, sy);
127 : 6359917 : z = itor(x, nbits2prec(-e));
128 : 6359917 : setsigne(z, sx); return z;
129 : : }
130 : :
131 : 134734725 : ly = lg(y);
132 [ + + ]: 134734725 : if (e > 0)
133 : : {
134 : 14990659 : l = ly - divsBIL(e);
135 [ + + ]: 14990659 : if (l < 3) return rcopy_sign(y, sy);
136 : : }
137 : 119744066 : else l = ly + nbits2extraprec(-e);
138 : 134544495 : z = (GEN)avma;
139 : 134544495 : y = addrr_sign(itor(x,l), sx, y, sy);
140 [ + + ]: 746130791 : ly = lg(y); while (ly--) *--z = y[ly];
141 : 192393888 : avma = (pari_sp)z; return z;
142 : : }
143 : :
144 : : static GEN
145 : 169738023 : addsr_sign(long x, GEN y, long sy)
146 : : {
147 : : long e, l, ly, sx;
148 : : GEN z;
149 : :
150 [ + + ]: 169738023 : if (!x) return rcopy_sign(y, sy);
151 [ + + ]: 169737729 : if (x < 0) { sx = -1; x = -x; } else sx = 1;
152 : 169737729 : e = expo(y) - expu(x);
153 [ + + ]: 169737729 : if (!sy)
154 : : {
155 [ - + ]: 346502 : if (e >= 0) return rcopy_sign(y, sy);
156 [ + + ]: 346502 : if (sx == -1) x = -x;
157 : 346502 : return stor(x, nbits2prec(-e));
158 : : }
159 : :
160 : 169391227 : ly = lg(y);
161 [ + + ]: 169391227 : if (e > 0)
162 : : {
163 : 563612 : l = ly - divsBIL(e);
164 [ + + ]: 563612 : if (l < 3) return rcopy_sign(y, sy);
165 : : }
166 : 168827615 : else l = ly + nbits2extraprec(-e);
167 : 169366399 : z = (GEN)avma;
168 : 169366399 : y = addrr_sign(stor(x,l), sx, y, sy);
169 [ + + ]: 1376555821 : ly = lg(y); while (ly--) *--z = y[ly];
170 : 169738023 : avma = (pari_sp)z; return z;
171 : : }
172 : :
173 : : GEN
174 : 166372963 : addsr(long x, GEN y) { return addsr_sign(x, y, signe(y)); }
175 : :
176 : : GEN
177 : 3365060 : subsr(long x, GEN y) { return addsr_sign(x, y, -signe(y)); }
178 : :
179 : : /* return x + 1, assuming x > 0 is a normalized t_REAL of exponent 0 */
180 : : GEN
181 : 0 : addrex01(GEN x)
182 : : {
183 : 0 : long l = lg(x);
184 : 0 : GEN y = cgetr(l);
185 : 0 : y[1] = evalsigne(1) | _evalexpo(1);
186 : 0 : y[2] = HIGHBIT | ((uel(x,2) & ~HIGHBIT) >> 1);
187 : 0 : shift_right(y, x, 3,l, x[2], 1);
188 : 0 : return y;
189 : : }
190 : : /* return subrs(x,1) to the min of (prec(x), prec(x-1) + 1),
191 : : * assuming x > 1 is a normalized t_REAL of exponent 0
192 : : * [ goal: avoid the loss of significant bits form subrs ]*/
193 : : GEN
194 : 0 : subrex01(GEN x)
195 : : {
196 : 0 : long i, sh, k, ly, lx = lg(x);
197 : : ulong u;
198 : : GEN y;
199 : 0 : k = 2;
200 : 0 : u = uel(x,2) & (~HIGHBIT);
201 [ # # ]: 0 : while (!u) u = x[++k]; /* terminates: x not a power of 2 */
202 [ # # ]: 0 : ly = (k == 2)? lx: lx - k+3; /* NB: +3, not +2: 1 extra word */
203 : 0 : y = cgetr(ly);
204 [ # # ][ # # ]: 0 : sh = bfffo(u);
[ # # ][ # # ]
205 [ # # ][ # # ]: 0 : if (sh)
206 : 0 : shift_left(y+2, x+k, 0, lx-k-1, 0, sh);
207 : : else
208 [ # # ]: 0 : { for (i = 2; i < lx-k+2; i++) y[i] = x[k-2 + i]; }
209 [ # # ]: 0 : for (i = lx-k+2; i < ly; i++) y[i] = 0;
210 : 0 : y[1] = evalsigne(1) | evalexpo(- (bit_accuracy(k) + sh));
211 : 0 : return y;
212 : : }
213 : :
214 : : GEN
215 : 2627077567 : addrr_sign(GEN x, long sx, GEN y, long sy)
216 : : {
217 : 2627077567 : long lx, ex = expo(x);
218 : 2627077567 : long ly, ey = expo(y), e = ey - ex;
219 : : long i, j, lz, ez, m;
220 : : int extend, f2;
221 : : GEN z;
222 : : LOCAL_OVERFLOW;
223 : :
224 [ + + ]: 2627077567 : if (!sy)
225 : : {
226 [ + + ]: 41670430 : if (!sx)
227 : : {
228 [ + + ]: 10051276 : if (e > 0) ex = ey;
229 : 10051276 : return real_0_bit(ex);
230 : : }
231 [ + + ]: 31619154 : if (e >= 0) return real_0_bit(ey);
232 : 31523344 : lz = nbits2prec(-e);
233 [ + + ]: 31523344 : lx = lg(x); if (lz > lx) lz = lx;
234 [ + + ]: 132159135 : z = cgetr(lz); while(--lz) z[lz] = x[lz];
235 : 31523344 : setsigne(z,sx); return z;
236 : : }
237 [ + + ]: 2585407137 : if (!sx)
238 : : {
239 [ + + ]: 17097828 : if (e <= 0) return real_0_bit(ex);
240 : 16904343 : lz = nbits2prec(e);
241 [ + + ]: 16904343 : ly = lg(y); if (lz > ly) lz = ly;
242 [ + + ]: 129130467 : z = cgetr(lz); while (--lz) z[lz] = y[lz];
243 : 16904343 : setsigne(z,sy); return z;
244 : : }
245 : :
246 [ + + ]: 2568309309 : if (e < 0) { swap(x,y); lswap(sx,sy); ey=ex; e=-e; }
247 : : /* now ey >= ex */
248 : 2568309309 : lx = lg(x);
249 : 2568309309 : ly = lg(y);
250 : : /* If exponents differ, need to shift one argument, here x. If
251 : : * extend = 1: extension of x,z by m < BIL bits (round to 1 word) */
252 : : /* in this case, lz = lx + d + 1, otherwise lx + d */
253 : 2568309309 : extend = 0;
254 [ + + ]: 2568309309 : if (e)
255 : : {
256 : 2264863880 : long d = dvmdsBIL(e, &m), l = ly-d;
257 [ + + ]: 2264863880 : if (l <= 2) return rcopy_sign(y, sy);
258 [ + + ]: 2252578017 : if (l > lx) { lz = lx + d + 1; extend = 1; }
259 : 1795537249 : else { lz = ly; lx = l; }
260 [ + + ]: 2252578017 : if (m)
261 : : { /* shift x right m bits */
262 : 2251297019 : const pari_sp av = avma;
263 : 2251297019 : const ulong sh = BITS_IN_LONG-m;
264 : 2251297019 : GEN p1 = x; x = new_chunk(lx + lz + 1);
265 : 2251297019 : shift_right(x,p1,2,lx, 0,m);
266 [ + + ]: 2251297019 : if (extend) x[lx] = p1[lx-1] << sh;
267 : 2251297019 : avma = av; /* HACK: cgetr(lz) will not overwrite x */
268 : : }
269 : : }
270 : : else
271 : : { /* d = 0 */
272 : 303445429 : m = 0;
273 [ + + ]: 303445429 : if (lx > ly) lx = ly;
274 : 303445429 : lz = lx;
275 : : }
276 : :
277 [ + + ]: 2556023446 : if (sx == sy)
278 : : { /* addition */
279 : 1071295070 : i = lz-1;
280 : 1071295070 : j = lx-1;
281 [ + + ]: 1071295070 : if (extend) {
282 : 209383007 : ulong garde = addll(x[lx], y[i]);
283 [ + + ][ + + ]: 209383007 : if (m < 4) /* don't extend for few correct bits */
284 : 38973674 : z = cgetr(--lz);
285 : : else
286 : : {
287 : 170409333 : z = cgetr(lz);
288 : 170409333 : z[i] = garde;
289 : : }
290 : : }
291 : : else
292 : : {
293 : 861912063 : z = cgetr(lz);
294 : 861912063 : z[i] = addll(x[j], y[i]); j--;
295 : : }
296 : 1071295070 : i--;
297 [ + + + ]: 3033419455 : for (; j>=2; i--,j--) z[i] = addllx(x[j],y[i]);
298 [ + + ]: 1071295070 : if (overflow)
299 : : {
300 : 309624206 : z[1] = 1; /* stops since z[1] != 0 */
301 [ + + ]: 310746028 : for (;;) { z[i] = uel(y,i)+1; if (z[i--]) break; }
302 [ + + ]: 309624206 : if (i <= 0)
303 : : {
304 : 308807688 : shift_right(z,z, 2,lz, 1,1);
305 : 308807688 : z[1] = evalsigne(sx) | evalexpo(ey+1); return z;
306 : : }
307 : : }
308 [ + + ]: 853204791 : for (; i>=2; i--) z[i] = y[i];
309 : 762487382 : z[1] = evalsigne(sx) | evalexpo(ey); return z;
310 : : }
311 : :
312 : : /* subtraction */
313 [ + + ]: 1484728376 : if (e) f2 = 1;
314 : : else
315 : : {
316 [ + + ][ + + ]: 214394736 : i = 2; while (i < lx && x[i] == y[i]) i++;
317 [ + + ]: 183866241 : if (i==lx) return real_0_bit(ey+1 - bit_accuracy(lx));
318 : 180379249 : f2 = (uel(y,i) > uel(x,i));
319 : : }
320 : : /* result is non-zero. f2 = (y > x) */
321 : 1481241384 : i = lz-1; z = cgetr(lz);
322 [ + + ]: 1481241384 : if (f2)
323 : : {
324 : 1396998354 : j = lx-1;
325 [ + + ]: 1396998354 : if (extend) z[i] = subll(y[i], x[lx]);
326 : 1149340593 : else z[i] = subll(y[i], x[j--]);
327 [ + + + ]: 3272098200 : for (i--; j>=2; i--) z[i] = subllx(y[i], x[j--]);
328 [ + + ]: 1396998354 : if (overflow) /* stops since y[1] != 0 */
329 [ + + ]: 8878395 : for (;;) { z[i] = uel(y,i)-1; if (y[i--]) break; }
330 [ + + ]: 1438831875 : for (; i>=2; i--) z[i] = y[i];
331 : 1396998354 : sx = sy;
332 : : }
333 : : else
334 : : {
335 [ - + ]: 84243030 : if (extend) z[i] = subll(x[lx], y[i]);
336 : 84243030 : else z[i] = subll(x[i], y[i]);
337 [ + + + ]: 180352080 : for (i--; i>=2; i--) z[i] = subllx(x[i], y[i]);
338 : : }
339 : :
340 [ + + ]: 1500399707 : x = z+2; i = 0; while (!x[i]) i++;
341 : 1481241384 : lz -= i; z += i;
342 [ + + ][ + + ]: 1481241384 : j = bfffo(z[2]); /* need to shift left by j bits to normalize mantissa */
[ + + ][ + + ]
343 : 1481241384 : ez = ey - (j | (i * BITS_IN_LONG));
344 [ + + ][ + + ]: 1481241384 : if (extend)
345 : : { /* z was extended by d+1 words [should be e bits = d words + m bits] */
346 : : /* not worth keeping extra word if less than 5 significant bits in there */
347 [ + + ][ + + ]: 247657761 : if (m - j < 5 && lz > 3)
348 : 70828818 : { /* shorten z */
349 : 70828818 : ulong last = (ulong)z[--lz]; /* cancelled word */
350 : :
351 : : /* if we need to shift anyway, shorten from left
352 : : * If not, shorten from right, neutralizing last word of z */
353 [ + + ]: 70828818 : if (j == 0)
354 : : /* stackdummy((pari_sp)(z + lz+1), (pari_sp)(z + lz)); */
355 : 35292309 : z[lz] = evaltyp(t_VECSMALL) | _evallg(1);
356 : : else
357 : : {
358 : 35536509 : GEN t = z;
359 : 35536509 : z++; shift_left(z,t,2,lz-1, last,j);
360 : : }
361 [ + + ]: 70828818 : if ((last<<j) & HIGHBIT)
362 : : { /* round up */
363 : 31892691 : i = lz-1;
364 [ + + ][ + - ]: 31957973 : while (++((ulong*)z)[i] == 0 && i > 1) i--;
365 [ + + ]: 31892691 : if (i == 1) { ez++; z[2] = (long)HIGHBIT; }
366 : : }
367 : : }
368 [ + + ]: 247657761 : else if (j) shift_left(z,z,2,lz-1, 0,j);
369 : : }
370 [ + + ]: 1233583623 : else if (j) shift_left(z,z,2,lz-1, 0,j);
371 : 1481241384 : z[1] = evalsigne(sx) | evalexpo(ez);
372 : 1481241384 : z[0] = evaltyp(t_REAL) | evallg(lz);
373 : 2627077567 : avma = (pari_sp)z; return z;
374 : : }
|
