Branch data 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 : 360 : inegate(GEN z) { return subsi(-1,z); }
24 : :
25 : : GEN
26 : 44 : binaire(GEN x)
27 : : {
28 : : ulong m,u;
29 : 44 : long i,lx,ex,ly,tx=typ(x);
30 : : GEN y,p1,p2;
31 : :
32 [ + + + + ]: 44 : switch(tx)
33 : : {
34 : : case t_INT:
35 [ + + ]: 24 : if (!signe(x))
36 : 8 : return cgetg(1,t_VEC);
37 : : else
38 : : {
39 : 16 : GEN xp=int_MSW(x);
40 : 16 : lx=lgefint(x);
41 : 16 : y = cgetg(2 + expi(x), t_VEC); ly=1;
42 [ + + ][ + + ]: 16 : m=HIGHBIT>>bfffo(*xp); u=*xp;
[ + + ][ - + ]
43 [ + + ][ + + ]: 272 : do { gel(y,ly) = m & u ? gen_1 : gen_0; ly++; } while (m>>=1);
44 [ + + ]: 32 : for (i=3; i<lx; i++)
45 : : {
46 : 16 : m=HIGHBIT; xp=int_precW(xp); u=*xp;
47 [ + + ][ + + ]: 800 : do { gel(y,ly) = m & u ? gen_1 : gen_0; ly++; } while (m>>=1);
48 : : }
49 : 16 : break;
50 : : }
51 : : case t_REAL:
52 : 12 : ex=expo(x);
53 [ + + ]: 12 : if (!signe(x)) return const_vec(maxss(-ex,0), gen_0);
54 : :
55 : 8 : lx=lg(x); y=cgetg(3,t_VEC);
56 [ - + ]: 8 : if (ex > bit_prec(x)) pari_err_PREC("binary");
57 : 8 : p1 = cgetg(maxss(ex,0)+2,t_VEC);
58 : 8 : p2 = cgetg(bit_prec(x)-ex,t_VEC);
59 : 8 : gel(y,1) = p1;
60 : 8 : gel(y,2) = p2;
61 : 8 : ly = -ex; ex++; m = HIGHBIT;
62 [ + + ]: 8 : if (ex<=0)
63 : : {
64 [ + + ]: 32 : gel(p1,1) = gen_0; for (i=1; i <= -ex; i++) gel(p2,i) = gen_0;
65 : 4 : i=2;
66 : : }
67 : : else
68 : : {
69 : 4 : ly=1;
70 [ + - ][ + + ]: 8 : for (i=2; i<lx && ly<=ex; i++)
71 : : {
72 : 4 : m=HIGHBIT; u=x[i];
73 : : do
74 [ + - ]: 4 : { gel(p1,ly) = (m & u) ? gen_1 : gen_0; ly++; }
75 [ + - ][ - + ]: 4 : while ((m>>=1) && ly<=ex);
76 : : }
77 : 4 : ly=1;
78 [ + - ]: 4 : if (m) i--; else m=HIGHBIT;
79 : : }
80 [ + + ]: 28 : for (; i<lx; i++)
81 : : {
82 : 20 : u=x[i];
83 [ + + ][ + + ]: 1020 : do { gel(p2,ly) = m & u ? gen_1 : gen_0; ly++; } while (m>>=1);
84 : 20 : m=HIGHBIT;
85 : : }
86 : 8 : break;
87 : :
88 : : case t_VEC: case t_COL: case t_MAT:
89 : 4 : y = cgetg_copy(x, &lx);
90 [ + + ]: 12 : for (i=1; i<lx; i++) gel(y,i) = binaire(gel(x,i));
91 : 4 : break;
92 : 4 : default: pari_err_TYPE("binary",x);
93 : 0 : return NULL; /* not reached */
94 : : }
95 : 40 : return y;
96 : : }
97 : :
98 : : /* return 1 if bit n of x is set, 0 otherwise */
99 : : long
100 : 56 : bittest(GEN x, long n)
101 : : {
102 [ - + ]: 56 : if (typ(x) != t_INT) pari_err_TYPE("bittest",x);
103 [ + - ][ - + ]: 56 : if (!signe(x) || n < 0) return 0;
104 [ + + ]: 56 : if (signe(x) < 0)
105 : : {
106 : 4 : pari_sp ltop=avma;
107 : 4 : long b = !int_bit(inegate(x),n);
108 : 4 : avma=ltop;
109 : 4 : return b;
110 : : }
111 : 56 : return int_bit(x, n);
112 : : }
113 : :
114 : : GEN
115 : 56 : gbittest(GEN x, long n) { return map_proto_lGL(bittest,x,n); }
116 : :
117 : : /***********************************************************************/
118 : : /** **/
119 : : /** BITMAP OPS **/
120 : : /** x & y (and), x | y (or), x ^ y (xor), ~x (neg), x & ~y (negimply) **/
121 : : /** **/
122 : : /***********************************************************************/
123 : : /* Truncate a non-negative integer to a number of bits. */
124 : : static GEN
125 : 20 : ibittrunc(GEN x, long bits)
126 : : {
127 : 20 : long lowbits, known_zero_words, xl = lgefint(x) - 2;
128 : 20 : long len_out = nbits2nlong(bits);
129 : :
130 [ + + ]: 20 : if (xl < len_out)
131 : 5 : return x;
132 : : /* Check whether mask is trivial */
133 : 15 : lowbits = bits & (BITS_IN_LONG-1);
134 [ + + ]: 15 : if (!lowbits) {
135 [ + - ]: 3 : if (xl == len_out)
136 : 3 : return x;
137 [ + - ]: 12 : } else if (len_out <= xl) {
138 : 12 : GEN xi = int_W(x, len_out-1);
139 : : /* Non-trival mask is given by a formula, if x is not
140 : : normalized, this works even in the exceptional case */
141 : 12 : *xi &= (1L << lowbits) - 1;
142 [ + + ][ - + ]: 12 : if (*xi && xl == len_out) return x;
143 : : }
144 : : /* Normalize */
145 : 12 : known_zero_words = xl - len_out;
146 [ - + ]: 12 : if (known_zero_words < 0) known_zero_words = 0;
147 : 20 : return int_normalize(x, known_zero_words);
148 : : }
149 : :
150 : : GEN
151 : 64 : gbitneg(GEN x, long bits)
152 : : {
153 : 64 : const ulong uzero = 0;
154 : : long lowbits, xl, len_out, i;
155 : :
156 [ + + ]: 64 : if (typ(x) != t_INT) pari_err_TYPE("bitwise negation",x);
157 [ + + ]: 60 : if (bits < -1)
158 : 4 : pari_err_DOMAIN("bitwise negation","exponent","<",gen_m1,stoi(bits));
159 [ + + ]: 56 : if (bits == -1) return inegate(x);
160 [ - + ]: 32 : if (bits == 0) return gen_0;
161 [ + + ]: 32 : if (signe(x) < 0) { /* Consider as if mod big power of 2 */
162 : 12 : pari_sp ltop = avma;
163 : 12 : return gerepileuptoint(ltop, ibittrunc(inegate(x), bits));
164 : : }
165 : 20 : xl = lgefint(x);
166 : 20 : len_out = nbits2lg(bits);
167 : 20 : lowbits = bits & (BITS_IN_LONG-1);
168 [ + + ]: 20 : if (len_out > xl) /* Need to grow */
169 : : {
170 : 12 : GEN out, outp, xp = int_MSW(x);
171 : 12 : out = cgetipos(len_out);
172 : 12 : outp = int_MSW(out);
173 [ + + ]: 12 : if (!lowbits)
174 : 4 : *outp = ~uzero;
175 : : else
176 : 8 : *outp = (1L << lowbits) - 1;
177 [ + + ]: 20 : for (i = 3; i < len_out - xl + 2; i++)
178 : : {
179 : 8 : outp = int_precW(outp); *outp = ~uzero;
180 : : }
181 [ + + ]: 20 : for ( ; i < len_out; i++)
182 : : {
183 : 8 : outp = int_precW(outp); *outp = ~*xp;
184 : 8 : xp = int_precW(xp);
185 : : }
186 : 12 : return out;
187 : : }
188 : 8 : x = icopy(x);
189 [ + + ]: 31 : for (i = 2; i < xl; i++) x[i] = ~x[i];
190 : 56 : return ibittrunc(int_normalize(x,0), bits);
191 : : }
192 : :
193 : : /* bitwise 'and' of two positive integers (any integers, but we ignore sign).
194 : : * Inputs are not necessary normalized. */
195 : : GEN
196 : 76 : ibitand(GEN x, GEN y)
197 : : {
198 : : long lx, ly, lout;
199 : : long *xp, *yp, *outp;
200 : : GEN out;
201 : : long i;
202 : :
203 [ + + ][ - + ]: 76 : if (!signe(x) || !signe(y)) return gen_0;
204 : 52 : lx=lgefint(x); ly=lgefint(y);
205 : 52 : lout = minss(lx,ly); /* > 2 */
206 : 52 : xp = int_LSW(x);
207 : 52 : yp = int_LSW(y);
208 : 52 : out = cgetipos(lout);
209 : 52 : outp = int_LSW(out);
210 [ + + ]: 149 : for (i=2; i<lout; i++)
211 : : {
212 : 97 : *outp = (*xp) & (*yp);
213 : 97 : outp = int_nextW(outp);
214 : 97 : xp = int_nextW(xp);
215 : 97 : yp = int_nextW(yp);
216 : : }
217 [ + + ]: 52 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
218 : 76 : return out;
219 : : }
220 : :
221 : : /* bitwise 'or' of absolute values of two integers */
222 : : GEN
223 : 60 : ibitor(GEN x, GEN y)
224 : : {
225 : : long lx, ly;
226 : : long *xp, *yp, *outp;
227 : : GEN out;
228 : : long i;
229 [ + + ]: 60 : if (!signe(x)) return absi(y);
230 [ - + ]: 44 : if (!signe(y)) return absi(x);
231 : :
232 : 44 : lx = lgefint(x); xp = int_LSW(x);
233 : 44 : ly = lgefint(y); yp = int_LSW(y);
234 [ + + ]: 44 : if (lx < ly) swapspec(xp,yp,lx,ly);
235 : : /* lx > 2 */
236 : 44 : out = cgetipos(lx);
237 : 44 : outp = int_LSW(out);
238 [ + + ]: 118 : for (i=2;i<ly;i++)
239 : : {
240 : 74 : *outp = (*xp) | (*yp);
241 : 74 : outp = int_nextW(outp);
242 : 74 : xp = int_nextW(xp);
243 : 74 : yp = int_nextW(yp);
244 : : }
245 [ + + ]: 89 : for ( ;i<lx;i++)
246 : : {
247 : 45 : *outp = *xp;
248 : 45 : outp = int_nextW(outp);
249 : 45 : xp = int_nextW(xp);
250 : : }
251 : : /* If input is normalized, this is not needed */
252 [ - + ]: 44 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
253 : 60 : return out;
254 : : }
255 : :
256 : : /* bitwise 'xor' of absolute values of two integers */
257 : : GEN
258 : 84 : ibitxor(GEN x, GEN y)
259 : : {
260 : : long lx, ly;
261 : : long *xp, *yp, *outp;
262 : : GEN out;
263 : : long i;
264 [ + + ]: 84 : if (!signe(x)) return absi(y);
265 [ - + ]: 60 : if (!signe(y)) return absi(x);
266 : :
267 : 60 : lx = lgefint(x); xp = int_LSW(x);
268 : 60 : ly = lgefint(y); yp = int_LSW(y);
269 [ + + ]: 60 : if (lx < ly) swapspec(xp,yp,lx,ly);
270 : : /* lx > 2 */
271 : 60 : out = cgetipos(lx);
272 : 60 : outp = int_LSW(out);
273 [ + + ]: 165 : for (i=2;i<ly;i++)
274 : : {
275 : 105 : *outp = (*xp) ^ (*yp);
276 : 105 : outp = int_nextW(outp);
277 : 105 : xp = int_nextW(xp);
278 : 105 : yp = int_nextW(yp);
279 : : }
280 [ + + ]: 120 : for ( ;i<lx;i++)
281 : : {
282 : 60 : *outp = *xp;
283 : 60 : outp = int_nextW(outp);
284 : 60 : xp = int_nextW(xp);
285 : : }
286 [ + + ]: 60 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
287 : 84 : return out;
288 : : }
289 : :
290 : : /* bitwise 'negimply' of absolute values of two integers */
291 : : /* "negimply(x,y)" is ~(x => y) == ~(~x | y) == x & ~y */
292 : : GEN
293 : 116 : ibitnegimply(GEN x, GEN y)
294 : : {
295 : : long lx, ly, lin;
296 : : long *xp, *yp, *outp;
297 : : GEN out;
298 : : long i;
299 [ + + ]: 116 : if (!signe(x)) return gen_0;
300 [ + + ]: 92 : if (!signe(y)) return absi(x);
301 : :
302 : 84 : lx = lgefint(x); xp = int_LSW(x);
303 : 84 : ly = lgefint(y); yp = int_LSW(y);
304 : 84 : lin = minss(lx,ly);
305 : 84 : out = cgetipos(lx);
306 : 84 : outp = int_LSW(out);
307 [ + + ]: 228 : for (i=2; i<lin; i++)
308 : : {
309 : 144 : *outp = (*xp) & ~(*yp);
310 : 144 : outp = int_nextW(outp);
311 : 144 : xp = int_nextW(xp);
312 : 144 : yp = int_nextW(yp);
313 : : }
314 [ + + ]: 124 : for ( ;i<lx;i++)
315 : : {
316 : 40 : *outp = *xp;
317 : 40 : outp = int_nextW(outp);
318 : 40 : xp = int_nextW(xp);
319 : : }
320 [ + + ]: 84 : if ( !*int_MSW(out) ) out = int_normalize(out, 1);
321 : 116 : return out;
322 : : }
323 : :
324 : : static int
325 [ + + ]: 336 : signs(GEN x, GEN y) { return (((signe(x) >= 0) << 1) | (signe(y) >= 0)); }
326 : : static void
327 : 448 : checkint2(const char *f,GEN x, GEN y)
328 [ + + ][ + + ]: 448 : { if (typ(x)!=t_INT || typ(y)!=t_INT) pari_err_TYPE2(f,x,y); }
329 : :
330 : : GEN
331 : 112 : gbitor(GEN x, GEN y)
332 : : {
333 : 112 : pari_sp ltop = avma;
334 : : GEN z;
335 : :
336 : 112 : checkint2("bitwise or",x,y);
337 [ + + + + ]: 84 : switch (signs(x, y))
338 : : {
339 : : case 3: /*1,1*/
340 : 40 : return ibitor(x,y);
341 : : case 2: /*1,-1*/
342 : 24 : z = ibitnegimply(inegate(y),x);
343 : 24 : break;
344 : : case 1: /*-1,1*/
345 : 8 : z = ibitnegimply(inegate(x),y);
346 : 8 : break;
347 : : default: /*-1,-1*/
348 : 12 : z = ibitand(inegate(x),inegate(y));
349 : 12 : break;
350 : : }
351 : 84 : return gerepileuptoint(ltop, inegate(z));
352 : : }
353 : :
354 : : GEN
355 : 112 : gbitand(GEN x, GEN y)
356 : : {
357 : 112 : pari_sp ltop = avma;
358 : : GEN z;
359 : :
360 : 112 : checkint2("bitwise and",x,y);
361 [ + + + + ]: 84 : switch (signs(x, y))
362 : : {
363 : : case 3: /*1,1*/
364 : 40 : return ibitand(x,y);
365 : : case 2: /*1,-1*/
366 : 24 : z = ibitnegimply(x,inegate(y));
367 : 24 : break;
368 : : case 1: /*-1,1*/
369 : 8 : z = ibitnegimply(y,inegate(x));
370 : 8 : break;
371 : : default: /*-1,-1*/
372 : 12 : z = inegate(ibitor(inegate(x),inegate(y)));
373 : 12 : break;
374 : : }
375 : 84 : return gerepileuptoint(ltop, z);
376 : : }
377 : :
378 : : GEN
379 : 112 : gbitxor(GEN x, GEN y)
380 : : {
381 : 112 : pari_sp ltop = avma;
382 : : GEN z;
383 : :
384 : 112 : checkint2("bitwise xor",x,y);
385 [ + + + + ]: 84 : switch (signs(x, y))
386 : : {
387 : : case 3: /*1,1*/
388 : 40 : return ibitxor(x,y);
389 : : case 2: /*1,-1*/
390 : 24 : z = inegate(ibitxor(x,inegate(y)));
391 : 24 : break;
392 : : case 1: /*-1,1*/
393 : 8 : z = inegate(ibitxor(inegate(x),y));
394 : 8 : break;
395 : : default: /*-1,-1*/
396 : 12 : z = ibitxor(inegate(x),inegate(y));
397 : 12 : break;
398 : : }
399 : 84 : return gerepileuptoint(ltop,z);
400 : : }
401 : :
402 : : /* x & ~y */
403 : : GEN
404 : 112 : gbitnegimply(GEN x, GEN y)
405 : : {
406 : 112 : pari_sp ltop = avma;
407 : : GEN z;
408 : :
409 : 112 : checkint2("bitwise negated imply",x,y);
410 [ + + + + ]: 84 : switch (signs(x, y))
411 : : {
412 : : case 3: /*1,1*/
413 : 40 : return ibitnegimply(x,y);
414 : : case 2: /*1,-1*/
415 : 24 : z = ibitand(x,inegate(y));
416 : 24 : break;
417 : : case 1: /*-1,1*/
418 : 8 : z = inegate(ibitor(y,inegate(x)));
419 : 8 : break;
420 : : default: /*-1,-1*/
421 : 12 : z = ibitnegimply(inegate(y),inegate(x));
422 : 12 : break;
423 : : }
424 : 84 : return gerepileuptoint(ltop,z);
425 : : }
426 : :
427 : : INLINE long
428 : 4 : hamming_word(ulong w)
429 : : {
430 : : #if 0
431 : : return __builtin_popcountl(w);
432 : : #endif
433 : : static long byte_weight[] = {
434 : : 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,
435 : : 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,
436 : : 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,
437 : : 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,
438 : : 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,
439 : : 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,
440 : : 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,
441 : : 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
442 : : };
443 : 4 : long sum = 0;
444 [ + + ]: 8 : while (w) { sum += byte_weight[w & 255]; w >>= 8; }
445 : 4 : return sum;
446 : : }
447 : :
448 : : /* number of non-zero entries among x[a], ..., x[b] */
449 : : static long
450 : 408 : hamming_slice(GEN x, long a, long b)
451 : : {
452 : 408 : long i, nb = 0;
453 [ + + ]: 40824 : for (i = a; i <= b; i++)
454 [ + + ]: 40416 : if (!gequal0(gel(x,i))) nb++;
455 : 408 : return nb;
456 : : }
457 : : static long
458 : 4 : hamming_mat(GEN x)
459 : : {
460 : 4 : long i, lx = lg(x), nb = 0;
461 [ + + ]: 404 : for (i = 1; i < lx; i++) nb += hammingweight(gel(x,i));
462 : 4 : return nb;
463 : : }
464 : : static long
465 : 4 : hamming_vecsmall(GEN x)
466 : : {
467 : 4 : long i, lx = lg(x), nb = 0;
468 [ + + ]: 20 : for (i = 1; i < lx; i++)
469 [ + + ]: 16 : if (x[i]) nb++;
470 : 4 : return nb;
471 : : }
472 : : static long
473 : 4 : hamming_int(GEN n)
474 : : {
475 : 4 : long lx = lgefint(n), i, sum;
476 [ - + ]: 4 : if (lx == 2) return 0;
477 : 4 : sum = hamming_word(n[2]);
478 [ - + ]: 4 : for (i = 3; i < lx; i++) sum += hamming_word(n[i]);
479 : 4 : return sum;
480 : : }
481 : :
482 : : long
483 : 424 : hammingweight(GEN n)
484 : : {
485 [ + + + + : 424 : switch(typ(n))
+ + ]
486 : : {
487 : 4 : case t_INT: return hamming_int(n);
488 : : case t_VEC:
489 : 404 : case t_COL: return hamming_slice(n, 1, lg(n)-1);
490 : 4 : case t_POL: return hamming_slice(n, 2, lg(n)-1);
491 : 4 : case t_VECSMALL: return hamming_vecsmall(n);
492 : 4 : case t_MAT: return hamming_mat(n);
493 : : }
494 : 4 : pari_err_TYPE("hammingweight", n);
495 : 420 : return 0;/*notreached*/
496 : : }
|
