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