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