Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : #include "pari.h"
15 : #include "paripriv.h"
16 : /*********************************************************************/
17 : /** **/
18 : /** BINARY DECOMPOSITION **/
19 : /** **/
20 : /*********************************************************************/
21 :
22 : INLINE GEN
23 630 : inegate(GEN z) { return subsi(-1,z); }
24 :
25 : GEN
26 42 : binary_zv(GEN x)
27 : {
28 : GEN xp, z;
29 : long i, k, lx;
30 42 : if (!signe(x)) return cgetg(1,t_VECSMALL);
31 28 : xp = int_LSW(x);
32 28 : lx = lgefint(x);
33 28 : k = expi(x)+2;
34 28 : z = cgetg(k, t_VECSMALL);
35 28 : k--;
36 174 : for(i = 2; i < lx; i++)
37 : {
38 174 : ulong u = *xp;
39 : long j;
40 174 : for (j=0; j<BITS_IN_LONG && k; j++) z[k--] = (u>>j)&1UL;
41 174 : if (!k) break;
42 146 : xp = int_nextW(xp);
43 : }
44 28 : return z;
45 : }
46 : static GEN
47 42 : F2v_to_ZV_inplace(GEN v)
48 : {
49 42 : long i, l = lg(v);
50 42 : v[0] = evaltyp(t_VEC) | _evallg(l);
51 42 : for (i = 1; i < l; i++) gel(v,i) = v[i]? gen_1: gen_0;
52 42 : return v;
53 : }
54 : /* "vector" of l bits (possibly no code word) to non-negative t_INT */
55 : GEN
56 0 : bits_to_int(GEN x, long l)
57 : {
58 : long i, j, lz;
59 : GEN z, zp;
60 :
61 0 : if (!l) return gen_0;
62 0 : lz = nbits2lg(l);
63 0 : z = cgetg(lz, t_INT);
64 0 : z[1] = evalsigne(1) | evallgefint(lz);
65 0 : zp = int_LSW(z); *zp = 0;
66 0 : for(i=l,j=0; i; i--,j++)
67 : {
68 0 : if (j==BITS_IN_LONG) { j=0; zp = int_nextW(zp); *zp = 0; }
69 0 : if (x[i]) *zp |= 1UL<<j;
70 : }
71 0 : return int_normalize(z, 0);
72 : }
73 : /* "vector" of l < BITS_IN_LONG bits (possibly no code word) to non-negative
74 : * ulong */
75 : ulong
76 0 : bits_to_u(GEN v, long l)
77 : {
78 0 : ulong u = 0;
79 : long i;
80 0 : for (i = 1; i <= l; i++) u = (u <<1) | v[i];
81 0 : return u;
82 : }
83 :
84 : /* set BITS_IN_LONG bits starting at word *w plus *r bits,
85 : * clearing subsequent bits in the last word touched */
86 : INLINE void
87 0 : int_set_ulong(ulong a, GEN *w, long *r)
88 : {
89 0 : if (*r) {
90 0 : **w |= (a << *r);
91 0 : *w = int_nextW(*w);
92 0 : **w = (a >> (BITS_IN_LONG - *r));
93 : } else {
94 0 : **w = a;
95 0 : *w = int_nextW(*w);
96 : }
97 0 : }
98 :
99 : /* set k bits starting at word *w plus *r bits,
100 : * clearing subsequent bits in the last word touched */
101 : INLINE void
102 123537 : int_set_bits(ulong a, long k, GEN *w, long *r)
103 : {
104 123537 : if (*r) {
105 99079 : **w |= a << *r;
106 99079 : a >>= BITS_IN_LONG - *r;
107 : } else {
108 24458 : **w = a;
109 24458 : a = 0;
110 : }
111 123537 : *r += k;
112 123537 : if (*r >= BITS_IN_LONG) {
113 81927 : *w = int_nextW(*w);
114 81927 : *r -= BITS_IN_LONG;
115 113427 : for (; *r >= BITS_IN_LONG; *r -= BITS_IN_LONG) {
116 31500 : **w = a;
117 31500 : a = 0;
118 31500 : *w = int_nextW(*w);
119 : }
120 81927 : if (*r)
121 81902 : **w = a;
122 : }
123 123537 : }
124 :
125 : /* set k bits from z (t_INT) starting at word *w plus *r bits,
126 : * clearing subsequent bits in the last word touched */
127 : INLINE void
128 1631 : int_set_int(GEN z, long k, GEN *w, long *r)
129 : {
130 1631 : long l = lgefint(z) - 2;
131 : GEN y;
132 1631 : if (!l) {
133 896 : int_set_bits(0, k, w, r);
134 2527 : return;
135 : }
136 735 : y = int_LSW(z);
137 735 : for (; l > 1; l--) {
138 0 : int_set_ulong((ulong) *y, w, r);
139 0 : y = int_nextW(y);
140 0 : k -= BITS_IN_LONG;
141 : }
142 735 : if (k)
143 735 : int_set_bits((ulong) *y, k, w, r);
144 : }
145 :
146 : GEN
147 24426 : nv_fromdigits_2k(GEN x, long k)
148 : {
149 24426 : long l = lg(x) - 1, r;
150 : GEN w, z;
151 24426 : if (k == 1) return bits_to_int(x, l);
152 24426 : if (!l) return gen_0;
153 24426 : z = cgetipos(nbits2lg(k * l));
154 24426 : w = int_LSW(z);
155 24426 : r = 0;
156 146332 : for (; l; l--)
157 121906 : int_set_bits(uel(x, l), k, &w, &r);
158 24426 : return int_normalize(z, 0);
159 : }
160 :
161 : GEN
162 7 : fromdigits_2k(GEN x, long k)
163 : {
164 : long l, m;
165 : GEN w, y, z;
166 7 : for (l = lg(x) - 1; l && !signe(gel(x, 1)); x++, l--);
167 7 : if (!l) return gen_0;
168 7 : m = expi(gel(x, 1)) + 1;
169 7 : z = cgetipos(nbits2lg(k * (l - 1) + m));
170 7 : w = int_LSW(z);
171 7 : if (!(k & (BITS_IN_LONG - 1))) {
172 0 : long i, j, t = k >> TWOPOTBITS_IN_LONG;
173 0 : for (; l; l--) {
174 0 : j = lgefint(gel(x, l)) - 2;
175 0 : y = int_LSW(gel(x, l));
176 0 : for (i = 0; i < j; i++) {
177 0 : *w = *y;
178 0 : y = int_nextW(y);
179 0 : w = int_nextW(w);
180 : }
181 0 : for (; i < t; i++) {
182 0 : *w = 0;
183 0 : w = int_nextW(w);
184 : }
185 : }
186 : } else {
187 7 : long r = 0;
188 1631 : for (; l > 1; l--)
189 1624 : int_set_int(gel(x, l), k, &w, &r);
190 7 : int_set_int(gel(x, 1), m, &w, &r);
191 : }
192 7 : return int_normalize(z, 0);
193 : }
194 :
195 : GEN
196 77 : binaire(GEN x)
197 : {
198 : ulong m,u;
199 77 : long i,lx,ex,ly,tx=typ(x);
200 : GEN y,p1,p2;
201 :
202 77 : switch(tx)
203 : {
204 : case t_INT:
205 42 : return F2v_to_ZV_inplace( binary_zv(x) );
206 : case t_REAL:
207 21 : ex = expo(x);
208 21 : if (!signe(x)) return const_vec(maxss(-ex,0), gen_0);
209 :
210 14 : lx=lg(x); y=cgetg(3,t_VEC);
211 14 : if (ex > bit_prec(x)) pari_err_PREC("binary");
212 14 : p1 = cgetg(maxss(ex,0)+2,t_VEC);
213 14 : p2 = cgetg(bit_prec(x)-ex,t_VEC);
214 14 : gel(y,1) = p1;
215 14 : gel(y,2) = p2;
216 14 : ly = -ex; ex++; m = HIGHBIT;
217 14 : if (ex<=0)
218 : {
219 7 : gel(p1,1) = gen_0; for (i=1; i <= -ex; i++) gel(p2,i) = gen_0;
220 7 : i=2;
221 : }
222 : else
223 : {
224 7 : ly=1;
225 14 : for (i=2; i<lx && ly<=ex; i++)
226 : {
227 7 : m=HIGHBIT; u=x[i];
228 : do
229 7 : { gel(p1,ly) = (m & u) ? gen_1 : gen_0; ly++; }
230 7 : while ((m>>=1) && ly<=ex);
231 : }
232 7 : ly=1;
233 7 : if (m) i--; else m=HIGHBIT;
234 : }
235 46 : for (; i<lx; i++)
236 : {
237 32 : u=x[i];
238 1785 : do { gel(p2,ly) = m & u ? gen_1 : gen_0; ly++; } while (m>>=1);
239 32 : m=HIGHBIT;
240 : }
241 14 : break;
242 :
243 : case t_VEC: case t_COL: case t_MAT:
244 7 : y = cgetg_copy(x, &lx);
245 7 : for (i=1; i<lx; i++) gel(y,i) = binaire(gel(x,i));
246 7 : break;
247 7 : default: pari_err_TYPE("binary",x);
248 0 : return NULL; /* not reached */
249 : }
250 21 : return y;
251 : }
252 :
253 : /* extract k bits (as ulong) starting at word *w plus *r bits */
254 : INLINE ulong
255 117223 : int_get_bits(long k, GEN *w, long *r)
256 : {
257 117223 : ulong mask = (1UL << k) - 1;
258 117223 : ulong a = (((ulong) **w) >> *r) & mask;
259 117223 : *r += k;
260 117223 : if (*r >= BITS_IN_LONG) {
261 25455 : *r -= BITS_IN_LONG;
262 25455 : *w = int_nextW(*w);
263 25455 : if (*r)
264 21800 : a |= ((ulong)**w << (k - *r)) & mask;
265 : }
266 117223 : return a;
267 : }
268 :
269 : /* extract BITS_IN_LONG bits starting at word *w plus *r bits */
270 : INLINE ulong
271 92446 : int_get_ulong(GEN *w, long *r)
272 : {
273 92446 : ulong a = ((ulong) **w) >> *r;
274 92446 : *w = int_nextW(*w);
275 92446 : if (*r)
276 80742 : a |= ((ulong)**w << (BITS_IN_LONG - *r));
277 92446 : return a;
278 : }
279 :
280 : /* extract k bits (as t_INT) starting at word *w plus *r bits */
281 : INLINE GEN
282 63000 : int_get_int(long k, GEN *w, long *r)
283 : {
284 63000 : GEN z = cgetipos(nbits2lg(k));
285 63000 : GEN y = int_LSW(z);
286 155446 : for (; k >= BITS_IN_LONG; k -= BITS_IN_LONG) {
287 92446 : *y = int_get_ulong(w, r);
288 92446 : y = int_nextW(y);
289 : }
290 63000 : if (k)
291 58121 : *y = int_get_bits(k, w, r);
292 63000 : return int_normalize(z, 0);
293 : }
294 :
295 : /* assume k < BITS_IN_LONG */
296 : GEN
297 6577 : binary_2k_nv(GEN x, long k)
298 : {
299 : long l, n, r;
300 : GEN v, w;
301 6577 : if (k == 1) return binary_zv(x);
302 6577 : if (!signe(x)) return cgetg(1, t_VECSMALL);
303 6577 : n = expi(x) + 1;
304 6577 : l = (n + k - 1) / k;
305 6577 : v = cgetg(l + 1, t_VECSMALL);
306 6577 : w = int_LSW(x);
307 6577 : r = 0;
308 59102 : for (; l > 1; l--) {
309 52525 : uel(v, l) = int_get_bits(k, &w, &r);
310 52525 : n -= k;
311 : }
312 6577 : uel(v, 1) = int_get_bits(n, &w, &r);
313 6577 : return v;
314 : }
315 :
316 : GEN
317 7014 : binary_2k(GEN x, long k)
318 : {
319 : long l, n;
320 : GEN v, w, y, z;
321 7014 : if (k == 1) return binaire(x);
322 7014 : if (!signe(x)) return cgetg(1, t_VEC);
323 7014 : n = expi(x) + 1;
324 7014 : l = (n + k - 1) / k;
325 7014 : v = cgetg(l + 1, t_VEC);
326 7014 : w = int_LSW(x);
327 7014 : if (!(k & (BITS_IN_LONG - 1))) {
328 14 : long m, t = k >> TWOPOTBITS_IN_LONG, u = lgefint(x) - 2;
329 56 : for (; l; l--) {
330 42 : m = minss(t, u);
331 42 : z = cgetipos(m + 2);
332 42 : y = int_LSW(z);
333 88 : for (; m; m--) {
334 46 : *y = *w;
335 46 : y = int_nextW(y);
336 46 : w = int_nextW(w);
337 : }
338 42 : gel(v, l) = int_normalize(z, 0);
339 42 : u -= t;
340 : }
341 : } else {
342 7000 : long r = 0;
343 63000 : for (; l > 1; l--, n -= k)
344 56000 : gel(v, l) = int_get_int(k, &w, &r);
345 7000 : gel(v, 1) = int_get_int(n, &w, &r);
346 : }
347 7014 : return v;
348 : }
349 :
350 : /* return 1 if bit n of x is set, 0 otherwise */
351 : long
352 91 : bittest(GEN x, long n)
353 : {
354 91 : if (typ(x) != t_INT) pari_err_TYPE("bittest",x);
355 91 : if (!signe(x) || n < 0) return 0;
356 91 : if (signe(x) < 0)
357 : {
358 7 : pari_sp ltop=avma;
359 7 : long b = !int_bit(inegate(x),n);
360 7 : avma=ltop;
361 7 : return b;
362 : }
363 84 : return int_bit(x, n);
364 : }
365 :
366 : GEN
367 91 : gbittest(GEN x, long n) { return map_proto_lGL(bittest,x,n); }
368 :
369 : /***********************************************************************/
370 : /** **/
371 : /** BITMAP OPS **/
372 : /** x & y (and), x | y (or), x ^ y (xor), ~x (neg), x & ~y (negimply) **/
373 : /** **/
374 : /***********************************************************************/
375 : /* Truncate a non-negative integer to a number of bits. */
376 : static GEN
377 35 : ibittrunc(GEN x, long bits)
378 : {
379 35 : long lowbits, known_zero_words, xl = lgefint(x) - 2;
380 35 : long len_out = nbits2nlong(bits);
381 :
382 35 : if (xl < len_out)
383 8 : return x;
384 : /* Check whether mask is trivial */
385 27 : lowbits = bits & (BITS_IN_LONG-1);
386 27 : if (!lowbits) {
387 6 : if (xl == len_out)
388 6 : return x;
389 21 : } else if (len_out <= xl) {
390 21 : GEN xi = int_W(x, len_out-1);
391 : /* Non-trival mask is given by a formula, if x is not
392 : normalized, this works even in the exceptional case */
393 21 : *xi &= (1L << lowbits) - 1;
394 21 : if (*xi && xl == len_out) return x;
395 : }
396 : /* Normalize */
397 21 : known_zero_words = xl - len_out;
398 21 : if (known_zero_words < 0) known_zero_words = 0;
399 21 : return int_normalize(x, known_zero_words);
400 : }
401 :
402 : GEN
403 112 : gbitneg(GEN x, long bits)
404 : {
405 112 : const ulong uzero = 0;
406 : long lowbits, xl, len_out, i;
407 :
408 112 : if (typ(x) != t_INT) pari_err_TYPE("bitwise negation",x);
409 105 : if (bits < -1)
410 7 : pari_err_DOMAIN("bitwise negation","exponent","<",gen_m1,stoi(bits));
411 98 : if (bits == -1) return inegate(x);
412 56 : if (bits == 0) return gen_0;
413 56 : if (signe(x) < 0) { /* Consider as if mod big power of 2 */
414 21 : pari_sp ltop = avma;
415 21 : return gerepileuptoint(ltop, ibittrunc(inegate(x), bits));
416 : }
417 35 : xl = lgefint(x);
418 35 : len_out = nbits2lg(bits);
419 35 : lowbits = bits & (BITS_IN_LONG-1);
420 35 : if (len_out > xl) /* Need to grow */
421 : {
422 21 : GEN out, outp, xp = int_MSW(x);
423 21 : out = cgetipos(len_out);
424 21 : outp = int_MSW(out);
425 21 : if (!lowbits)
426 7 : *outp = ~uzero;
427 : else
428 14 : *outp = (1L << lowbits) - 1;
429 32 : for (i = 3; i < len_out - xl + 2; i++)
430 : {
431 11 : outp = int_precW(outp); *outp = ~uzero;
432 : }
433 35 : for ( ; i < len_out; i++)
434 : {
435 14 : outp = int_precW(outp); *outp = ~*xp;
436 14 : xp = int_precW(xp);
437 : }
438 21 : return out;
439 : }
440 14 : x = icopy(x);
441 14 : for (i = 2; i < xl; i++) x[i] = ~x[i];
442 14 : return ibittrunc(int_normalize(x,0), bits);
443 : }
444 :
445 : /* bitwise 'and' of two positive integers (any integers, but we ignore sign).
446 : * Inputs are not necessary normalized. */
447 : GEN
448 35849240 : ibitand(GEN x, GEN y)
449 : {
450 : long lx, ly, lout;
451 : long *xp, *yp, *outp;
452 : GEN out;
453 : long i;
454 :
455 35849240 : if (!signe(x) || !signe(y)) return gen_0;
456 35849198 : lx=lgefint(x); ly=lgefint(y);
457 35849198 : lout = minss(lx,ly); /* > 2 */
458 35849198 : xp = int_LSW(x);
459 35849198 : yp = int_LSW(y);
460 35849198 : out = cgetipos(lout);
461 35849198 : outp = int_LSW(out);
462 74192628 : for (i=2; i<lout; i++)
463 : {
464 38343430 : *outp = (*xp) & (*yp);
465 38343430 : outp = int_nextW(outp);
466 38343430 : xp = int_nextW(xp);
467 38343430 : yp = int_nextW(yp);
468 : }
469 35849198 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
470 35849198 : return out;
471 : }
472 :
473 : /* bitwise 'or' of absolute values of two integers */
474 : GEN
475 105 : ibitor(GEN x, GEN y)
476 : {
477 : long lx, ly;
478 : long *xp, *yp, *outp;
479 : GEN out;
480 : long i;
481 105 : if (!signe(x)) return absi(y);
482 77 : if (!signe(y)) return absi(x);
483 :
484 77 : lx = lgefint(x); xp = int_LSW(x);
485 77 : ly = lgefint(y); yp = int_LSW(y);
486 77 : if (lx < ly) swapspec(xp,yp,lx,ly);
487 : /* lx > 2 */
488 77 : out = cgetipos(lx);
489 77 : outp = int_LSW(out);
490 202 : for (i=2;i<ly;i++)
491 : {
492 125 : *outp = (*xp) | (*yp);
493 125 : outp = int_nextW(outp);
494 125 : xp = int_nextW(xp);
495 125 : yp = int_nextW(yp);
496 : }
497 149 : for ( ;i<lx;i++)
498 : {
499 72 : *outp = *xp;
500 72 : outp = int_nextW(outp);
501 72 : xp = int_nextW(xp);
502 : }
503 : /* If input is normalized, this is not needed */
504 77 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
505 77 : return out;
506 : }
507 :
508 : /* bitwise 'xor' of absolute values of two integers */
509 : GEN
510 147 : ibitxor(GEN x, GEN y)
511 : {
512 : long lx, ly;
513 : long *xp, *yp, *outp;
514 : GEN out;
515 : long i;
516 147 : if (!signe(x)) return absi(y);
517 105 : if (!signe(y)) return absi(x);
518 :
519 105 : lx = lgefint(x); xp = int_LSW(x);
520 105 : ly = lgefint(y); yp = int_LSW(y);
521 105 : if (lx < ly) swapspec(xp,yp,lx,ly);
522 : /* lx > 2 */
523 105 : out = cgetipos(lx);
524 105 : outp = int_LSW(out);
525 282 : for (i=2;i<ly;i++)
526 : {
527 177 : *outp = (*xp) ^ (*yp);
528 177 : outp = int_nextW(outp);
529 177 : xp = int_nextW(xp);
530 177 : yp = int_nextW(yp);
531 : }
532 201 : for ( ;i<lx;i++)
533 : {
534 96 : *outp = *xp;
535 96 : outp = int_nextW(outp);
536 96 : xp = int_nextW(xp);
537 : }
538 105 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
539 105 : return out;
540 : }
541 :
542 : /* bitwise 'negimply' of absolute values of two integers */
543 : /* "negimply(x,y)" is ~(x => y) == ~(~x | y) == x & ~y */
544 : GEN
545 203 : ibitnegimply(GEN x, GEN y)
546 : {
547 : long lx, ly, lin;
548 : long *xp, *yp, *outp;
549 : GEN out;
550 : long i;
551 203 : if (!signe(x)) return gen_0;
552 161 : if (!signe(y)) return absi(x);
553 :
554 147 : lx = lgefint(x); xp = int_LSW(x);
555 147 : ly = lgefint(y); yp = int_LSW(y);
556 147 : lin = minss(lx,ly);
557 147 : out = cgetipos(lx);
558 147 : outp = int_LSW(out);
559 390 : for (i=2; i<lin; i++)
560 : {
561 243 : *outp = (*xp) & ~(*yp);
562 243 : outp = int_nextW(outp);
563 243 : xp = int_nextW(xp);
564 243 : yp = int_nextW(yp);
565 : }
566 211 : for ( ;i<lx;i++)
567 : {
568 64 : *outp = *xp;
569 64 : outp = int_nextW(outp);
570 64 : xp = int_nextW(xp);
571 : }
572 147 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
573 147 : return out;
574 : }
575 :
576 : static int
577 35849695 : signs(GEN x, GEN y) { return (((signe(x) >= 0) << 1) | (signe(y) >= 0)); }
578 : static void
579 35849891 : checkint2(const char *f,GEN x, GEN y)
580 35849891 : { if (typ(x)!=t_INT || typ(y)!=t_INT) pari_err_TYPE2(f,x,y); }
581 :
582 : GEN
583 196 : gbitor(GEN x, GEN y)
584 : {
585 196 : pari_sp ltop = avma;
586 : GEN z;
587 :
588 196 : checkint2("bitwise or",x,y);
589 147 : switch (signs(x, y))
590 : {
591 : case 3: /*1,1*/
592 70 : return ibitor(x,y);
593 : case 2: /*1,-1*/
594 42 : z = ibitnegimply(inegate(y),x);
595 42 : break;
596 : case 1: /*-1,1*/
597 14 : z = ibitnegimply(inegate(x),y);
598 14 : break;
599 : default: /*-1,-1*/
600 21 : z = ibitand(inegate(x),inegate(y));
601 21 : break;
602 : }
603 77 : return gerepileuptoint(ltop, inegate(z));
604 : }
605 :
606 : GEN
607 35849303 : gbitand(GEN x, GEN y)
608 : {
609 35849303 : pari_sp ltop = avma;
610 : GEN z;
611 :
612 35849303 : checkint2("bitwise and",x,y);
613 35849254 : switch (signs(x, y))
614 : {
615 : case 3: /*1,1*/
616 35849177 : return ibitand(x,y);
617 : case 2: /*1,-1*/
618 42 : z = ibitnegimply(x,inegate(y));
619 42 : break;
620 : case 1: /*-1,1*/
621 14 : z = ibitnegimply(y,inegate(x));
622 14 : break;
623 : default: /*-1,-1*/
624 21 : z = inegate(ibitor(inegate(x),inegate(y)));
625 21 : break;
626 : }
627 77 : return gerepileuptoint(ltop, z);
628 : }
629 :
630 : GEN
631 196 : gbitxor(GEN x, GEN y)
632 : {
633 196 : pari_sp ltop = avma;
634 : GEN z;
635 :
636 196 : checkint2("bitwise xor",x,y);
637 147 : switch (signs(x, y))
638 : {
639 : case 3: /*1,1*/
640 70 : return ibitxor(x,y);
641 : case 2: /*1,-1*/
642 42 : z = inegate(ibitxor(x,inegate(y)));
643 42 : break;
644 : case 1: /*-1,1*/
645 14 : z = inegate(ibitxor(inegate(x),y));
646 14 : break;
647 : default: /*-1,-1*/
648 21 : z = ibitxor(inegate(x),inegate(y));
649 21 : break;
650 : }
651 77 : return gerepileuptoint(ltop,z);
652 : }
653 :
654 : /* x & ~y */
655 : GEN
656 196 : gbitnegimply(GEN x, GEN y)
657 : {
658 196 : pari_sp ltop = avma;
659 : GEN z;
660 :
661 196 : checkint2("bitwise negated imply",x,y);
662 147 : switch (signs(x, y))
663 : {
664 : case 3: /*1,1*/
665 70 : return ibitnegimply(x,y);
666 : case 2: /*1,-1*/
667 42 : z = ibitand(x,inegate(y));
668 42 : break;
669 : case 1: /*-1,1*/
670 14 : z = inegate(ibitor(y,inegate(x)));
671 14 : break;
672 : default: /*-1,-1*/
673 21 : z = ibitnegimply(inegate(y),inegate(x));
674 21 : break;
675 : }
676 77 : return gerepileuptoint(ltop,z);
677 : }
678 :
679 : INLINE long
680 14811 : hamming_word(ulong w)
681 : {
682 : #if 0
683 : return __builtin_popcountl(w);
684 : #endif
685 : static long byte_weight[] = {
686 : 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
687 : 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
688 : 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
689 : 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
690 : 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
691 : 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
692 : 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
693 : 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
694 : };
695 14811 : long sum = 0;
696 14811 : while (w) { sum += byte_weight[w & 255]; w >>= 8; }
697 14811 : return sum;
698 : }
699 :
700 : /* number of non-zero entries among x[a], ..., x[b] */
701 : static long
702 714 : hamming_slice(GEN x, long a, long b)
703 : {
704 714 : long i, nb = 0;
705 71442 : for (i = a; i <= b; i++)
706 70728 : if (!gequal0(gel(x,i))) nb++;
707 714 : return nb;
708 : }
709 : static long
710 7 : hamming_mat(GEN x)
711 : {
712 7 : long i, lx = lg(x), nb = 0;
713 7 : for (i = 1; i < lx; i++) nb += hammingweight(gel(x,i));
714 7 : return nb;
715 : }
716 : static long
717 7 : hamming_vecsmall(GEN x)
718 : {
719 7 : long i, lx = lg(x), nb = 0;
720 35 : for (i = 1; i < lx; i++)
721 28 : if (x[i]) nb++;
722 7 : return nb;
723 : }
724 : static long
725 14665 : hamming_int(GEN n)
726 : {
727 14665 : long lx = lgefint(n), i, sum;
728 14665 : if (lx == 2) return 0;
729 14665 : sum = hamming_word(n[2]);
730 14665 : for (i = 3; i < lx; i++) sum += hamming_word(n[i]);
731 14665 : return sum;
732 : }
733 :
734 : long
735 15400 : hammingweight(GEN n)
736 : {
737 15400 : switch(typ(n))
738 : {
739 14665 : case t_INT: return hamming_int(n);
740 : case t_VEC:
741 707 : case t_COL: return hamming_slice(n, 1, lg(n)-1);
742 7 : case t_POL: return hamming_slice(n, 2, lg(n)-1);
743 7 : case t_VECSMALL: return hamming_vecsmall(n);
744 7 : case t_MAT: return hamming_mat(n);
745 : }
746 7 : pari_err_TYPE("hammingweight", n);
747 0 : return 0;/*notreached*/
748 : }
|
