Line data Source code
1 : /* Copyright (C) 2012-2019 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 : #define DEBUGLEVEL DEBUGLEVEL_mat
19 :
20 : /***********************************************************************/
21 : /** **/
22 : /** F2v **/
23 : /** **/
24 : /***********************************************************************/
25 : /* F2v objects are defined as follows:
26 : * An F2v is a t_VECSMALL:
27 : * v[0] = codeword
28 : * v[1] = number of components
29 : * x[2] = a_0...a_31 x[3] = a_32..a_63, etc. on 32bit
30 : * x[2] = a_0...a_63 x[3] = a_64..a_127, etc. on 64bit
31 : * where the a_i are bits. */
32 :
33 : int
34 4816 : F2v_equal0(GEN V)
35 : {
36 4816 : long l = lg(V);
37 5982 : while (--l > 1)
38 5422 : if (V[l]) return 0;
39 560 : return 1;
40 : }
41 :
42 : GEN
43 3950785 : F2c_to_ZC(GEN x)
44 : {
45 3950785 : long l = x[1]+1, lx = lg(x);
46 3950785 : GEN z = cgetg(l, t_COL);
47 : long i, j, k;
48 7924899 : for (i=2, k=1; i<lx; i++)
49 31783670 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
50 27809468 : gel(z,k) = (x[i]&(1UL<<j))? gen_1: gen_0;
51 3950697 : return z;
52 : }
53 : GEN
54 4704 : F2c_to_mod(GEN x)
55 : {
56 4704 : long l = x[1]+1, lx = lg(x);
57 4704 : GEN z = cgetg(l, t_COL);
58 4704 : GEN _0 = mkintmod(gen_0,gen_2);
59 4704 : GEN _1 = mkintmod(gen_1,gen_2);
60 : long i, j, k;
61 19058 : for (i=2, k=1; i<lx; i++)
62 604867 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
63 590513 : gel(z,k) = (x[i]&(1UL<<j))? _1: _0;
64 4704 : return z;
65 : }
66 :
67 : /* x[a..b], a <= b */
68 : GEN
69 28 : F2v_slice(GEN x, long a, long b)
70 : {
71 28 : long i,j,k, l = b-a+1;
72 28 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
73 28 : z[1] = l;
74 98 : for(i=a,k=1,j=BITS_IN_LONG; i<=b; i++,j++)
75 : {
76 70 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
77 70 : if (F2v_coeff(x,i)) z[k] |= 1UL<<j;
78 : }
79 28 : return z;
80 : }
81 : /* x[a..b,], a <= b */
82 : GEN
83 14 : F2m_rowslice(GEN x, long a, long b)
84 42 : { pari_APPLY_same(F2v_slice(gel(x,i),a,b)) }
85 :
86 : GEN
87 3927 : F2v_to_Flv(GEN x)
88 : {
89 3927 : long l = x[1]+1, lx = lg(x);
90 3927 : GEN z = cgetg(l, t_VECSMALL);
91 : long i,j,k;
92 7911 : for (i=2, k=1; i<lx; i++)
93 75034 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
94 71050 : z[k] = (x[i]>>j)&1UL;
95 3927 : return z;
96 : }
97 :
98 : GEN
99 4817197 : F2m_to_ZM(GEN x) { pari_APPLY_same(F2c_to_ZC(gel(x,i))) }
100 : GEN
101 4984 : F2m_to_mod(GEN x) { pari_APPLY_same(F2c_to_mod(gel(x,i))) }
102 : GEN
103 1596 : F2m_to_Flm(GEN x) { pari_APPLY_same(F2v_to_Flv(gel(x,i))) }
104 :
105 : GEN
106 1316 : RgV_F2v_extract_shallow(GEN V, GEN x)
107 : {
108 1316 : long n = F2v_hamming(x), m = 1;
109 1316 : long l = x[1]+1, lx = lg(x);
110 1316 : GEN W = cgetg(n+1, t_VEC);
111 : long i,j,k;
112 2632 : for (i=2, k=1; i<lx; i++)
113 17836 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
114 16520 : if ((x[i]>>j)&1UL)
115 2597 : gel(W, m++) = gel(V,k);
116 1316 : return W;
117 : }
118 :
119 : GEN
120 10764775 : ZV_to_F2v(GEN x)
121 : {
122 10764775 : long i, j, k, l = lg(x)-1;
123 10764775 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
124 10764418 : z[1] = l;
125 98197254 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
126 : {
127 87432550 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
128 87432550 : if (mpodd(gel(x,i))) z[k] |= 1UL<<j;
129 : }
130 10764704 : return z;
131 : }
132 :
133 : GEN
134 9177 : RgV_to_F2v(GEN x)
135 : {
136 9177 : long l = lg(x)-1;
137 9177 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
138 : long i,j,k;
139 9177 : z[1] = l;
140 1441356 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
141 : {
142 1432179 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
143 1432179 : if (Rg_to_F2(gel(x,i))) z[k] |= 1UL<<j;
144 : }
145 9177 : return z;
146 : }
147 :
148 : GEN
149 6482 : Flv_to_F2v(GEN x)
150 : {
151 6482 : long l = lg(x)-1;
152 6482 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
153 : long i,j,k;
154 6482 : z[1] = l;
155 10235078 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
156 : {
157 10228596 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
158 10228596 : if (x[i]&1L) z[k] |= 1UL<<j;
159 : }
160 6482 : return z;
161 : }
162 :
163 : GEN
164 13513738 : ZM_to_F2m(GEN x) { pari_APPLY_same(ZV_to_F2v(gel(x,i))) }
165 : GEN
166 9674 : RgM_to_F2m(GEN x) { pari_APPLY_same(RgV_to_F2v(gel(x,i))) }
167 : GEN
168 0 : Flm_to_F2m(GEN x) { pari_APPLY_same(Flv_to_F2v(gel(x,i))) }
169 :
170 : GEN
171 1665006 : const_F2v(long m)
172 : {
173 1665006 : long i, l = nbits2lg(m);
174 1665005 : GEN c = cgetg(l, t_VECSMALL);
175 1664940 : c[1] = m;
176 3342762 : for (i = 2; i < l; i++) uel(c,i) = -1UL;
177 1664940 : if (remsBIL(m)) uel(c,l-1) = (1UL<<remsBIL(m))-1UL;
178 1664927 : return c;
179 : }
180 :
181 : /* Allow lg(y)<lg(x) */
182 : void
183 26522252 : F2v_add_inplace(GEN x, GEN y)
184 : {
185 26522252 : long n = lg(y);
186 26522252 : long r = (n-2)&7L, q = n-r, i;
187 36006969 : for (i = 2; i < q; i += 8)
188 : {
189 9484717 : x[ i] ^= y[ i]; x[1+i] ^= y[1+i]; x[2+i] ^= y[2+i]; x[3+i] ^= y[3+i];
190 9484717 : x[4+i] ^= y[4+i]; x[5+i] ^= y[5+i]; x[6+i] ^= y[6+i]; x[7+i] ^= y[7+i];
191 : }
192 26522252 : switch (r)
193 : {
194 1777803 : case 7: x[i] ^= y[i]; i++; case 6: x[i] ^= y[i]; i++;
195 3476908 : case 5: x[i] ^= y[i]; i++; case 4: x[i] ^= y[i]; i++;
196 9590453 : case 3: x[i] ^= y[i]; i++; case 2: x[i] ^= y[i]; i++;
197 23749774 : case 1: x[i] ^= y[i]; i++;
198 : }
199 26522252 : }
200 :
201 : /* Allow lg(y)<lg(x) */
202 : void
203 14448 : F2v_and_inplace(GEN x, GEN y)
204 : {
205 14448 : long n = lg(y);
206 14448 : long r = (n-2)&7L, q = n-r, i;
207 14448 : for (i = 2; i < q; i += 8)
208 : {
209 0 : x[ i] &= y[ i]; x[1+i] &= y[1+i]; x[2+i] &= y[2+i]; x[3+i] &= y[3+i];
210 0 : x[4+i] &= y[4+i]; x[5+i] &= y[5+i]; x[6+i] &= y[6+i]; x[7+i] &= y[7+i];
211 : }
212 14448 : switch (r)
213 : {
214 0 : case 7: x[i] &= y[i]; i++; case 6: x[i] &= y[i]; i++;
215 0 : case 5: x[i] &= y[i]; i++; case 4: x[i] &= y[i]; i++;
216 2064 : case 3: x[i] &= y[i]; i++; case 2: x[i] &= y[i]; i++;
217 14448 : case 1: x[i] &= y[i]; i++;
218 : }
219 14448 : }
220 :
221 : /* Allow lg(y)<lg(x) */
222 : void
223 0 : F2v_or_inplace(GEN x, GEN y)
224 : {
225 0 : long n = lg(y);
226 0 : long r = (n-2)&7L, q = n-r, i;
227 0 : for (i = 2; i < q; i += 8)
228 : {
229 0 : x[ i] |= y[ i]; x[1+i] |= y[1+i]; x[2+i] |= y[2+i]; x[3+i] |= y[3+i];
230 0 : x[4+i] |= y[4+i]; x[5+i] |= y[5+i]; x[6+i] |= y[6+i]; x[7+i] |= y[7+i];
231 : }
232 0 : switch (r)
233 : {
234 0 : case 7: x[i] |= y[i]; i++; case 6: x[i] |= y[i]; i++;
235 0 : case 5: x[i] |= y[i]; i++; case 4: x[i] |= y[i]; i++;
236 0 : case 3: x[i] |= y[i]; i++; case 2: x[i] |= y[i]; i++;
237 0 : case 1: x[i] |= y[i]; i++;
238 : }
239 0 : }
240 :
241 : /* Allow lg(y)<lg(x) */
242 : void
243 1148 : F2v_negimply_inplace(GEN x, GEN y)
244 : {
245 1148 : long n = lg(y);
246 1148 : long r = (n-2)&7L, q = n-r, i;
247 91908 : for (i = 2; i < q; i += 8)
248 : {
249 90760 : x[ i] &= ~y[ i]; x[1+i] &= ~y[1+i]; x[2+i] &= ~y[2+i]; x[3+i] &= ~y[3+i];
250 90760 : x[4+i] &= ~y[4+i]; x[5+i] &= ~y[5+i]; x[6+i] &= ~y[6+i]; x[7+i] &= ~y[7+i];
251 : }
252 1148 : switch (r)
253 : {
254 56 : case 7: x[i] &= ~y[i]; i++; case 6: x[i] &= ~y[i]; i++;
255 56 : case 5: x[i] &= ~y[i]; i++; case 4: x[i] &= ~y[i]; i++;
256 212 : case 3: x[i] &= ~y[i]; i++; case 2: x[i] &= ~y[i]; i++;
257 1148 : case 1: x[i] &= ~y[i]; i++;
258 : }
259 1148 : }
260 :
261 : ulong
262 0 : F2v_dotproduct(GEN x, GEN y)
263 : {
264 0 : long i, lx = lg(x);
265 : ulong c;
266 0 : if (lx <= 2) return 0;
267 0 : c = uel(x,2) & uel(y,2);
268 0 : for (i=3; i<lx; i++) c ^= uel(x,i) & uel(y,i);
269 : #ifdef LONG_IS_64BIT
270 0 : c ^= c >> 32;
271 : #endif
272 0 : c ^= c >> 16;
273 0 : c ^= c >> 8;
274 0 : c ^= c >> 4;
275 0 : c ^= c >> 2;
276 0 : c ^= c >> 1;
277 0 : return c & 1;
278 : }
279 :
280 : ulong
281 2254 : F2v_hamming(GEN H)
282 : {
283 2254 : long i, n=0, l=lg(H);
284 1094242 : for (i=2; i<l; i++) n += hammingu(uel(H,i));
285 2254 : return n;
286 : }
287 :
288 : int
289 6342 : F2v_subset(GEN x, GEN y)
290 : {
291 6342 : long i, n = lg(y);
292 7476 : for (i = 2; i < n; i ++)
293 6748 : if ((x[i] & y[i]) != x[i]) return 0;
294 728 : return 1;
295 : }
296 :
297 : GEN
298 230895 : matid_F2m(long n)
299 : {
300 230895 : GEN y = cgetg(n+1,t_MAT);
301 : long i;
302 230893 : if (n < 0) pari_err_DOMAIN("matid_F2m", "dimension","<",gen_0,stoi(n));
303 989104 : for (i=1; i<=n; i++) { gel(y,i) = zero_F2v(n); F2v_set(gel(y,i),i); }
304 230896 : return y;
305 : }
306 :
307 : GEN
308 0 : F2m_row(GEN x, long j)
309 : {
310 0 : long i, l = lg(x);
311 0 : GEN V = zero_F2v(l-1);
312 0 : for(i = 1; i < l; i++)
313 0 : if (F2m_coeff(x,j,i))
314 0 : F2v_set(V,i);
315 0 : return V;
316 : }
317 :
318 : GEN
319 0 : F2m_transpose(GEN x)
320 : {
321 0 : long i, dx, lx = lg(x);
322 : GEN y;
323 0 : if (lx == 1) return cgetg(1,t_MAT);
324 0 : dx = coeff(x,1,1); y = cgetg(dx+1,t_MAT);
325 0 : for (i=1; i<=dx; i++) gel(y,i) = F2m_row(x,i);
326 0 : return y;
327 : }
328 :
329 : INLINE GEN
330 2193495 : F2m_F2c_mul_i(GEN x, GEN y, long lx, long l)
331 : {
332 : long j;
333 2193495 : GEN z = NULL;
334 :
335 23480323 : for (j=1; j<lx; j++)
336 : {
337 21286837 : if (!F2v_coeff(y,j)) continue;
338 5143262 : if (!z) z = vecsmall_copy(gel(x,j));
339 3218215 : else F2v_add_inplace(z,gel(x,j));
340 : }
341 2193486 : if (!z) z = zero_F2v(l);
342 2193495 : return z;
343 : }
344 :
345 : GEN
346 724309 : F2m_mul(GEN x, GEN y)
347 : {
348 724309 : long i,j,l,lx=lg(x), ly=lg(y);
349 : GEN z;
350 724309 : if (ly==1) return cgetg(1,t_MAT);
351 724309 : z = cgetg(ly,t_MAT);
352 724297 : if (lx==1)
353 : {
354 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
355 0 : return z;
356 : }
357 724297 : l = coeff(x,1,1);
358 2917790 : for (j=1; j<ly; j++) gel(z,j) = F2m_F2c_mul_i(x, gel(y,j), lx, l);
359 724300 : return z;
360 : }
361 :
362 : GEN
363 0 : F2m_F2c_mul(GEN x, GEN y)
364 : {
365 0 : long l, lx = lg(x);
366 0 : if (lx==1) return cgetg(1,t_VECSMALL);
367 0 : l = coeff(x,1,1);
368 0 : return F2m_F2c_mul_i(x, y, lx, l);
369 : }
370 :
371 : static GEN
372 0 : _F2m_mul(void *data, GEN x, GEN y) { (void) data; return F2m_mul(x,y); }
373 : static GEN
374 0 : _F2m_sqr(void *data, GEN x) { (void) data; return F2m_mul(x,x); }
375 : GEN
376 0 : F2m_powu(GEN x, ulong n)
377 : {
378 0 : if (!n) return matid(lg(x)-1);
379 0 : return gen_powu(x, n,NULL, &_F2m_sqr, &_F2m_mul);
380 : }
381 :
382 : static long
383 6322939 : F2v_find_nonzero(GEN x0, GEN mask0, long m)
384 : {
385 6322939 : ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
386 6322939 : long i, l = lg(x0)-2;
387 8861553 : for (i = 0; i < l; i++)
388 : {
389 6917822 : e = *x++ & *mask++;
390 6917822 : if (e) return i*BITS_IN_LONG+vals(e)+1;
391 : }
392 1943731 : return m+1;
393 : }
394 :
395 : /* in place, destroy x */
396 : GEN
397 1149042 : F2m_ker_sp(GEN x, long deplin)
398 : {
399 : GEN y, c, d;
400 : long i, j, k, r, m, n;
401 :
402 1149042 : n = lg(x)-1;
403 1149042 : m = mael(x,1,1); r=0;
404 :
405 1149042 : d = cgetg(n+1, t_VECSMALL);
406 1149035 : c = const_F2v(m);
407 5807363 : for (k=1; k<=n; k++)
408 : {
409 4961015 : GEN xk = gel(x,k);
410 4961015 : j = F2v_find_nonzero(xk, c, m);
411 4961039 : if (j>m)
412 : {
413 1473006 : if (deplin) {
414 302707 : GEN v = zero_F2v(n);
415 920292 : for (i=1; i<k; i++)
416 617588 : if (F2v_coeff(xk, d[i])) F2v_set(v, i);
417 302704 : F2v_set(v, k); return v;
418 : }
419 1170299 : r++; d[k] = 0;
420 : }
421 : else
422 : {
423 3488033 : F2v_clear(c,j); d[k] = j;
424 3488109 : F2v_clear(xk, j);
425 48662559 : for (i=k+1; i<=n; i++)
426 : {
427 45174432 : GEN xi = gel(x,i);
428 45174432 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
429 : }
430 3488127 : F2v_set(xk, j);
431 : }
432 : }
433 846348 : if (deplin) return NULL;
434 :
435 844780 : y = zero_F2m_copy(n,r);
436 2015076 : for (j=k=1; j<=r; j++,k++)
437 : {
438 3388797 : GEN C = gel(y,j); while (d[k]) k++;
439 9360754 : for (i=1; i<k; i++)
440 8190463 : if (d[i] && F2m_coeff(x,d[i],k)) F2v_set(C,i);
441 1170291 : F2v_set(C, k);
442 : }
443 844783 : return y;
444 : }
445 :
446 : /* not memory clean */
447 : GEN
448 200054 : F2m_ker(GEN x) { return F2m_ker_sp(F2m_copy(x), 0); }
449 : GEN
450 0 : F2m_deplin(GEN x) { return F2m_ker_sp(F2m_copy(x), 1); }
451 :
452 : ulong
453 1764 : F2m_det_sp(GEN x) { return !F2m_ker_sp(x, 1); }
454 :
455 : ulong
456 0 : F2m_det(GEN x)
457 : {
458 0 : pari_sp av = avma;
459 0 : ulong d = F2m_det_sp(F2m_copy(x));
460 0 : return gc_ulong(av, d);
461 : }
462 :
463 : /* Destroy x */
464 : GEN
465 356540 : F2m_gauss_pivot(GEN x, long *rr)
466 : {
467 : GEN c, d;
468 : long i, j, k, r, m, n;
469 :
470 356540 : n = lg(x)-1; if (!n) { *rr=0; return NULL; }
471 356540 : m = mael(x,1,1); r=0;
472 :
473 356540 : d = cgetg(n+1, t_VECSMALL);
474 356536 : c = const_F2v(m);
475 1718576 : for (k=1; k<=n; k++)
476 : {
477 1362037 : GEN xk = gel(x,k);
478 1362037 : j = F2v_find_nonzero(xk, c, m);
479 1362043 : if (j>m) { r++; d[k] = 0; }
480 : else
481 : {
482 891132 : F2v_clear(c,j); d[k] = j;
483 5979258 : for (i=k+1; i<=n; i++)
484 : {
485 5088148 : GEN xi = gel(x,i);
486 5088148 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
487 : }
488 : }
489 : }
490 :
491 356539 : *rr = r; return gc_const((pari_sp)d, d);
492 : }
493 :
494 : long
495 63 : F2m_rank(GEN x)
496 : {
497 63 : pari_sp av = avma;
498 : long r;
499 63 : (void)F2m_gauss_pivot(F2m_copy(x),&r);
500 63 : return gc_long(av, lg(x)-1 - r);
501 : }
502 :
503 : static GEN
504 14 : F2m_inv_upper_1_ind(GEN A, long index)
505 : {
506 14 : pari_sp av = avma;
507 14 : long n = lg(A)-1, i = index, j;
508 14 : GEN u = const_vecsmall(n, 0);
509 14 : u[i] = 1;
510 21 : for (i--; i>0; i--)
511 : {
512 7 : ulong m = F2m_coeff(A,i,i+1) & uel(u,i+1); /* j = i+1 */
513 7 : for (j=i+2; j<=n; j++) m ^= F2m_coeff(A,i,j) & uel(u,j);
514 7 : u[i] = m & 1;
515 : }
516 14 : return gc_leaf(av, Flv_to_F2v(u));
517 : }
518 : static GEN
519 7 : F2m_inv_upper_1(GEN A)
520 : {
521 : long i, l;
522 7 : GEN B = cgetg_copy(A, &l);
523 21 : for (i = 1; i < l; i++) gel(B,i) = F2m_inv_upper_1_ind(A, i);
524 7 : return B;
525 : }
526 :
527 : static GEN
528 758201 : F2_get_col(GEN b, GEN d, long li, long aco)
529 : {
530 758201 : long i, l = nbits2lg(aco);
531 758205 : GEN u = cgetg(l, t_VECSMALL);
532 758214 : u[1] = aco;
533 5525086 : for (i = 1; i <= li; i++)
534 4766893 : if (d[i]) /* d[i] can still be 0 if li > aco */
535 : {
536 4613494 : if (F2v_coeff(b, i))
537 1516064 : F2v_set(u, d[i]);
538 : else
539 3097416 : F2v_clear(u, d[i]);
540 : }
541 758193 : return u;
542 : }
543 :
544 : /* destroy a, b */
545 : GEN
546 230952 : F2m_gauss_sp(GEN a, GEN b)
547 : {
548 230952 : long i, j, k, l, li, bco, aco = lg(a)-1;
549 : GEN u, d;
550 :
551 230952 : if (!aco) return cgetg(1,t_MAT);
552 230952 : li = gel(a,1)[1];
553 230952 : d = zero_Flv(li);
554 230950 : bco = lg(b)-1;
555 977555 : for (i=1; i<=aco; i++)
556 : {
557 746658 : GEN ai = vecsmall_copy(gel(a,i));
558 746657 : if (!d[i] && F2v_coeff(ai, i))
559 501268 : k = i;
560 : else
561 1145876 : for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
562 : /* found a pivot on row k */
563 746670 : if (k > li) return NULL;
564 746635 : d[k] = i;
565 :
566 : /* Clear k-th row but column-wise instead of line-wise */
567 : /* a_ij -= a_i1*a1j/a_11
568 : line-wise grouping: L_j -= a_1j/a_11*L_1
569 : column-wise: C_i -= a_i1/a_11*C_1
570 : */
571 746635 : F2v_clear(ai,k);
572 5332317 : for (l=1; l<=aco; l++)
573 : {
574 4585711 : GEN al = gel(a,l);
575 4585711 : if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
576 : }
577 5360353 : for (l=1; l<=bco; l++)
578 : {
579 4613748 : GEN bl = gel(b,l);
580 4613748 : if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
581 : }
582 : }
583 230897 : u = cgetg(bco+1,t_MAT);
584 989109 : for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
585 230909 : return u;
586 : }
587 :
588 : GEN
589 56 : F2m_gauss(GEN a, GEN b)
590 : {
591 56 : pari_sp av = avma;
592 56 : if (lg(a) == 1) return cgetg(1,t_MAT);
593 56 : return gc_upto(av, F2m_gauss_sp(F2m_copy(a), F2m_copy(b)));
594 : }
595 : GEN
596 28 : F2m_F2c_gauss(GEN a, GEN b)
597 : {
598 28 : pari_sp av = avma;
599 28 : GEN z = F2m_gauss(a, mkmat(b));
600 28 : if (!z) return gc_NULL(av);
601 14 : if (lg(z) == 1) retgc_const(av, cgetg(1, t_VECSMALL));
602 14 : return gc_leaf(av, gel(z,1));
603 : }
604 :
605 : GEN
606 189 : F2m_inv(GEN a)
607 : {
608 189 : pari_sp av = avma;
609 189 : if (lg(a) == 1) return cgetg(1,t_MAT);
610 189 : return gc_upto(av, F2m_gauss_sp(F2m_copy(a), matid_F2m(gel(a,1)[1])));
611 : }
612 :
613 : GEN
614 7 : F2m_invimage_i(GEN A, GEN B)
615 : {
616 : GEN d, x, X, Y;
617 7 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
618 7 : x = F2m_ker_sp(shallowconcat(A, B), 0);
619 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
620 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
621 : * Y has at least nB columns and full rank */
622 7 : nY = lg(x)-1;
623 7 : if (nY < nB) return NULL;
624 :
625 : /* implicitly: Y = rowslice(x, nA+1, nA+nB), nB rows */
626 7 : d = cgetg(nB+1, t_VECSMALL);
627 21 : for (i = nB, j = nY; i >= 1; i--, j--)
628 : {
629 14 : for (; j>=1; j--)
630 14 : if (F2m_coeff(x,nA+i,j)) { d[i] = j; break; } /* Y[i,j] */
631 14 : if (!j) return NULL;
632 : }
633 7 : x = vecpermute(x, d);
634 :
635 7 : X = F2m_rowslice(x, 1, nA);
636 7 : Y = F2m_rowslice(x, nA+1, nA+nB);
637 7 : return F2m_mul(X, F2m_inv_upper_1(Y));
638 : }
639 : GEN
640 0 : F2m_invimage(GEN A, GEN B)
641 : {
642 0 : pari_sp av = avma;
643 0 : GEN X = F2m_invimage_i(A,B);
644 0 : if (!X) return gc_NULL(av);
645 0 : return gc_upto(av, X);
646 : }
647 :
648 : GEN
649 195759 : F2m_F2c_invimage(GEN A, GEN y)
650 : {
651 195759 : pari_sp av = avma;
652 195759 : long i, l = lg(A);
653 : GEN M, x;
654 :
655 195759 : if (l==1) return NULL;
656 195759 : if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
657 195759 : M = cgetg(l+1,t_MAT);
658 1016855 : for (i=1; i<l; i++) gel(M,i) = gel(A,i);
659 195757 : gel(M,l) = y; M = F2m_ker(M);
660 195758 : i = lg(M)-1; if (!i) return gc_NULL(av);
661 :
662 195758 : x = gel(M,i);
663 195758 : if (!F2v_coeff(x,l)) return gc_NULL(av);
664 195755 : F2v_clear(x, x[1]); x[1]--; /* remove last coord */
665 195754 : return gc_leaf(av, x);
666 : }
667 :
668 : /* Block Lanczos algorithm for kernel of sparse matrix (F2Ms)
669 : Based on lanczos.cpp by Jason Papadopoulos
670 : <https://github.com/sagemath/FlintQS/blob/master/src/lanczos.cpp>
671 : Copyright Jason Papadopoulos 2006
672 : Released under the GNU General Public License v2 or later version.
673 : */
674 :
675 : /* F2Ms are vector of vecsmall which represents nonzero entries of columns
676 : * F2w are vecsmall whoses entries are columns of a n x BIL F2m
677 : * F2wB are F2w in the special case where dim = BIL.
678 : */
679 :
680 : #define BIL BITS_IN_LONG
681 :
682 : static GEN
683 184 : F2w_transpose_F2m(GEN x)
684 : {
685 184 : long i, j, l = lg(x)-1;
686 184 : GEN z = cgetg(BIL+1, t_MAT);
687 11448 : for (j = 1; j <= BIL; j++)
688 11264 : gel(z,j) = zero_F2v(l);
689 291482 : for (i = 1; i <= l; i++)
690 : {
691 291298 : ulong xi = uel(x,i);
692 18411426 : for(j = 1; j <= BIL; j++)
693 18120128 : if (xi&(1UL<<(j-1)))
694 5323988 : F2v_set(gel(z, j), i);
695 : }
696 184 : return z;
697 : }
698 :
699 : static GEN
700 6612 : F2wB_mul(GEN a, GEN b)
701 : {
702 : long i, j;
703 6612 : GEN c = cgetg(BIL+1, t_VECSMALL);
704 407124 : for (i = 1; i <= BIL; i++)
705 : {
706 400512 : ulong s = 0, ai = a[i];
707 21073964 : for (j = 1; ai; j++, ai>>=1)
708 20673452 : if (ai & 1)
709 10408460 : s ^= b[j];
710 400512 : c[i] = s;
711 : }
712 6612 : return c;
713 : }
714 :
715 : static void
716 4408 : precompute_F2w_F2wB(GEN x, GEN c)
717 : {
718 : ulong z, xk;
719 : ulong i, j, k, index;
720 4408 : x++; c++;
721 37784 : for (j = 0; j < (BIL>>3); j++)
722 : {
723 8577632 : for (i = 0; i < 256; i++)
724 : {
725 8544256 : k = 0;
726 8544256 : index = i;
727 8544256 : z = 0;
728 68387424 : while (index) {
729 59843168 : xk = x[k];
730 59843168 : if (index & 1)
731 34177024 : z ^= xk;
732 59843168 : index >>= 1;
733 59843168 : k++;
734 : }
735 8544256 : c[i] = z;
736 : }
737 33376 : x += 8; c += 256;
738 : }
739 4408 : }
740 :
741 : static void
742 4408 : F2w_F2wB_mul_add_inplace(GEN v, GEN x, GEN y)
743 : {
744 4408 : long i, n = lg(y)-1;
745 : ulong word;
746 4408 : GEN c = cgetg(1+(BIL<<5), t_VECSMALL);
747 4408 : precompute_F2w_F2wB(x, c);
748 4408 : c++;
749 7725680 : for (i = 1; i <= n; i++)
750 : {
751 7721272 : word = v[i];
752 7721272 : y[i] ^= c[ 0*256 + ((word>> 0) & 0xff) ]
753 7721272 : ^ c[ 1*256 + ((word>> 8) & 0xff) ]
754 7721272 : ^ c[ 2*256 + ((word>>16) & 0xff) ]
755 7721272 : ^ c[ 3*256 + ((word>>24) & 0xff) ]
756 : #ifdef LONG_IS_64BIT
757 7248696 : ^ c[ 4*256 + ((word>>32) & 0xff) ]
758 7248696 : ^ c[ 5*256 + ((word>>40) & 0xff) ]
759 7248696 : ^ c[ 6*256 + ((word>>48) & 0xff) ]
760 7248696 : ^ c[ 7*256 + ((word>>56) ) ]
761 : #endif
762 : ;
763 : }
764 4408 : }
765 :
766 : /* Return x*y~, which is a F2wB */
767 : static GEN
768 3352 : F2w_transmul(GEN x, GEN y)
769 : {
770 3352 : long i, j, n = lg(x)-1;
771 3352 : GEN z = zero_zv(BIL);
772 3352 : pari_sp av = avma;
773 3352 : GEN c = zero_zv(BIL<<5) + 1;
774 3352 : GEN xy = z + 1;
775 :
776 5862977 : for (i = 1; i <= n; i++)
777 : {
778 5859625 : ulong xi = x[i];
779 5859625 : ulong yi = y[i];
780 5859625 : c[ 0*256 + ( xi & 0xff) ] ^= yi;
781 5859625 : c[ 1*256 + ((xi >> 8) & 0xff) ] ^= yi;
782 5859625 : c[ 2*256 + ((xi >> 16) & 0xff) ] ^= yi;
783 5859625 : c[ 3*256 + ((xi >> 24) & 0xff) ] ^= yi;
784 : #ifdef LONG_IS_64BIT
785 5501328 : c[ 4*256 + ((xi >> 32) & 0xff) ] ^= yi;
786 5501328 : c[ 5*256 + ((xi >> 40) & 0xff) ] ^= yi;
787 5501328 : c[ 6*256 + ((xi >> 48) & 0xff) ] ^= yi;
788 5501328 : c[ 7*256 + ((xi >> 56) ) ] ^= yi;
789 : #endif
790 : }
791 30168 : for(i = 0; i < 8; i++)
792 : {
793 26816 : ulong a0 = 0, a1 = 0, a2 = 0, a3 = 0;
794 : #ifdef LONG_IS_64BIT
795 23952 : ulong a4 = 0, a5 = 0, a6 = 0, a7 = 0;
796 : #endif
797 6891712 : for (j = 0; j < 256; j++) {
798 6864896 : if ((j >> i) & 1) {
799 3432448 : a0 ^= c[0*256 + j];
800 3432448 : a1 ^= c[1*256 + j];
801 3432448 : a2 ^= c[2*256 + j];
802 3432448 : a3 ^= c[3*256 + j];
803 : #ifdef LONG_IS_64BIT
804 3065856 : a4 ^= c[4*256 + j];
805 3065856 : a5 ^= c[5*256 + j];
806 3065856 : a6 ^= c[6*256 + j];
807 3065856 : a7 ^= c[7*256 + j];
808 : #endif
809 : }
810 : }
811 26816 : xy[ 0] = a0; xy[ 8] = a1; xy[16] = a2; xy[24] = a3;
812 : #ifdef LONG_IS_64BIT
813 23952 : xy[32] = a4; xy[40] = a5; xy[48] = a6; xy[56] = a7;
814 : #endif
815 26816 : xy++;
816 : }
817 3352 : return gc_const(av, z);
818 : }
819 :
820 : static GEN
821 1102 : identity_F2wB(void)
822 : {
823 : long i;
824 1102 : GEN M = cgetg(BIL+1, t_VECSMALL);
825 67854 : for (i = 1; i <= BIL; i++)
826 66752 : uel(M,i) = 1UL<<(i-1);
827 1102 : return M;
828 : }
829 :
830 : static GEN
831 1102 : find_nonsingular_sub(GEN t, GEN last_s, GEN *pt_s)
832 : {
833 1102 : long i, j, dim = 0;
834 : ulong mask, row_i, row_j;
835 1102 : long last_dim = lg(last_s)-1;
836 1102 : GEN s = cgetg(BIL+1, t_VECSMALL);
837 1102 : GEN M1 = identity_F2wB();
838 1102 : pari_sp av = avma;
839 1102 : GEN cols = cgetg(BIL+1, t_VECSMALL);
840 1102 : GEN M0 = zv_copy(t);
841 :
842 1102 : mask = 0;
843 66975 : for (i = 1; i <= last_dim; i++)
844 : {
845 65873 : cols[BIL + 1 - i] = last_s[i];
846 65873 : mask |= 1UL<<(last_s[i]-1);
847 : }
848 67854 : for (i = j = 1; i <= BIL; i++)
849 66752 : if (!(mask & (1UL<<(i-1))))
850 879 : cols[j++] = i;
851 :
852 : /* compute the inverse of t[][] */
853 :
854 67854 : for (i = 1; i <= BIL; i++)
855 : {
856 66752 : mask = 1UL<<(cols[i]-1);
857 66752 : row_i = cols[i];
858 130042 : for (j = i; j <= BIL; j++)
859 : {
860 128296 : row_j = cols[j];
861 128296 : if (uel(M0,row_j) & mask)
862 : {
863 65006 : swap(gel(M0, row_j), gel(M0, row_i));
864 65006 : swap(gel(M1, row_j), gel(M1, row_i));
865 65006 : break;
866 : }
867 : }
868 66752 : if (j <= BIL)
869 : {
870 4108590 : for (j = 1; j <= BIL; j++)
871 : {
872 4043584 : row_j = cols[j];
873 4043584 : if (row_i != row_j && (M0[row_j] & mask))
874 : {
875 1958783 : uel(M0,row_j) ^= uel(M0,row_i);
876 1958783 : uel(M1,row_j) ^= uel(M1,row_i);
877 : }
878 : }
879 65006 : s[++dim] = cols[i];
880 65006 : continue;
881 : }
882 1746 : for (j = i; j <= BIL; j++)
883 : {
884 1746 : row_j = cols[j];
885 1746 : if (uel(M1,row_j) & mask)
886 : {
887 1746 : swap(gel(M0, row_j), gel(M0, row_i));
888 1746 : swap(gel(M1, row_j), gel(M1, row_i));
889 1746 : break;
890 : }
891 : }
892 1746 : if (j > BIL) return NULL;
893 109458 : for (j = 1; j <= BIL; j++)
894 : {
895 107712 : row_j = cols[j];
896 107712 : if (row_i != row_j && (M1[row_j] & mask))
897 : {
898 0 : uel(M0,row_j) ^= uel(M0,row_i);
899 0 : uel(M1,row_j) ^= uel(M1,row_i);
900 : }
901 : }
902 1746 : M0[row_i] = M1[row_i] = 0;
903 : }
904 1102 : mask = 0;
905 66108 : for (i = 1; i <= dim; i++)
906 65006 : mask |= 1UL<<(s[i]-1);
907 66975 : for (i = 1; i <= last_dim; i++)
908 65873 : mask |= 1UL<<(last_s[i]-1);
909 1102 : if (mask != (ulong)(-1))
910 0 : return NULL; /* Failure */
911 1102 : setlg(s, dim+1);
912 1102 : set_avma(av);
913 1102 : *pt_s = s;
914 1102 : return M1;
915 : }
916 :
917 : /* Compute x * A~ */
918 : static GEN
919 1286 : F2w_F2Ms_transmul(GEN x, GEN A, long nbrow)
920 : {
921 1286 : long i, j, l = lg(A);
922 1286 : GEN b = zero_zv(nbrow);
923 2206288 : for (i = 1; i < l; i++)
924 : {
925 2205002 : GEN c = gel(A,i);
926 2205002 : long lc = lg(c);
927 2205002 : ulong xi = x[i];
928 41407938 : for (j = 1; j < lc; j++)
929 39202936 : b[c[j]] ^= xi;
930 : }
931 1286 : return b;
932 : }
933 :
934 : /* Compute x * A */
935 : static GEN
936 1194 : F2w_F2Ms_mul(GEN x, GEN A)
937 : {
938 1194 : long i, j, l = lg(A);
939 1194 : GEN b = cgetg(l, t_VECSMALL);
940 2068854 : for (i = 1; i < l; i++)
941 : {
942 2067660 : GEN c = gel(A,i);
943 2067660 : long lc = lg(c);
944 2067660 : ulong s = 0;
945 38865718 : for (j = 1; j < lc; j++)
946 36798058 : s ^= x[c[j]];
947 2067660 : b[i] = s;
948 : }
949 1194 : return b;
950 : }
951 :
952 : static void
953 2204 : F2wB_addid_inplace(GEN f)
954 : {
955 : long i;
956 135708 : for (i = 1; i <= BIL; i++)
957 133504 : uel(f,i) ^= 1UL<<(i-1);
958 2204 : }
959 :
960 : static void
961 2204 : F2w_mask_inplace(GEN f, ulong m)
962 : {
963 2204 : long i, l = lg(f);
964 1999274 : for (i = 1; i < l; i++)
965 1997070 : uel(f,i) &= m;
966 2204 : }
967 :
968 : static GEN
969 46 : block_lanczos(GEN B, ulong nbrow)
970 : {
971 46 : pari_sp av = avma, av2;
972 : GEN v0, v1, v2, vnext, x, w;
973 : GEN winv0, winv1, winv2;
974 : GEN vt_a_v0, vt_a_v1, vt_a2_v0, vt_a2_v1;
975 : GEN d, e, f, f2, s0;
976 : long i, iter;
977 46 : long n = lg(B)-1;
978 : long dim0;
979 : ulong mask0, mask1;
980 46 : v1 = zero_zv(n);
981 46 : v2 = zero_zv(n);
982 46 : vt_a_v1 = zero_zv(BIL);
983 46 : vt_a2_v1 = zero_zv(BIL);
984 46 : winv1 = zero_zv(BIL);
985 46 : winv2 = zero_zv(BIL);
986 46 : s0 = identity_zv(BIL);
987 46 : mask1 = (ulong)(-1);
988 :
989 46 : x = random_zv(n);
990 46 : w = F2w_F2Ms_mul(F2w_F2Ms_transmul(x, B, nbrow), B);
991 46 : v0 = w;
992 46 : av2 = avma;
993 46 : for (iter=1;;iter++)
994 : {
995 1148 : vnext = F2w_F2Ms_mul(F2w_F2Ms_transmul(v0, B, nbrow), B);
996 1148 : vt_a_v0 = F2w_transmul(v0, vnext);
997 1148 : if (zv_equal0(vt_a_v0)) break; /* success */
998 1102 : vt_a2_v0 = F2w_transmul(vnext, vnext);
999 1102 : winv0 = find_nonsingular_sub(vt_a_v0, s0, &s0);
1000 1102 : if (!winv0) return gc_NULL(av); /* failure */
1001 1102 : dim0 = lg(s0)-1;
1002 1102 : mask0 = 0;
1003 66108 : for (i = 1; i <= dim0; i++)
1004 65006 : mask0 |= 1UL<<(s0[i]-1);
1005 1102 : d = cgetg(BIL+1, t_VECSMALL);
1006 67854 : for (i = 1; i <= BIL; i++)
1007 66752 : d[i] = (vt_a2_v0[i] & mask0) ^ vt_a_v0[i];
1008 :
1009 1102 : d = F2wB_mul(winv0, d);
1010 1102 : F2wB_addid_inplace(d);
1011 1102 : e = F2wB_mul(winv1, vt_a_v0);
1012 1102 : F2w_mask_inplace(e, mask0);
1013 1102 : f = F2wB_mul(vt_a_v1, winv1);
1014 1102 : F2wB_addid_inplace(f);
1015 1102 : f = F2wB_mul(winv2, f);
1016 1102 : f2 = cgetg(BIL+1, t_VECSMALL);
1017 67854 : for (i = 1; i <= BIL; i++)
1018 66752 : f2[i] = ((vt_a2_v1[i] & mask1) ^ vt_a_v1[i]) & mask0;
1019 :
1020 1102 : f = F2wB_mul(f, f2);
1021 1102 : F2w_mask_inplace(vnext, mask0);
1022 1102 : F2w_F2wB_mul_add_inplace(v0, d, vnext);
1023 1102 : F2w_F2wB_mul_add_inplace(v1, e, vnext);
1024 1102 : F2w_F2wB_mul_add_inplace(v2, f, vnext);
1025 1102 : d = F2wB_mul(winv0, F2w_transmul(v0, w));
1026 1102 : F2w_F2wB_mul_add_inplace(v0, d, x);
1027 1102 : v2 = v1; v1 = v0; v0 = vnext;
1028 1102 : winv2 = winv1; winv1 = winv0;
1029 1102 : vt_a_v1 = vt_a_v0;
1030 1102 : vt_a2_v1 = vt_a2_v0;
1031 1102 : mask1 = mask0;
1032 1102 : (void)gc_all(av2, 9, &x, &s0, &v0, &v1, &v2,
1033 : &winv1, &winv2, &vt_a_v1, &vt_a2_v1);
1034 : }
1035 46 : if (DEBUGLEVEL >= 5)
1036 0 : err_printf("Lanczos halted after %ld iterations\n", iter);
1037 46 : v1 = F2w_F2Ms_transmul(x, B, nbrow);
1038 46 : v2 = F2w_F2Ms_transmul(v0, B, nbrow);
1039 46 : x = shallowconcat(F2w_transpose_F2m(x), F2w_transpose_F2m(v0));
1040 46 : v1 = shallowconcat(F2w_transpose_F2m(v1), F2w_transpose_F2m(v2));
1041 46 : s0 = gel(F2m_indexrank(x), 2);
1042 46 : x = shallowextract(x, s0);
1043 46 : v1 = shallowextract(v1, s0);
1044 46 : return F2m_mul(x, F2m_ker(v1));
1045 : }
1046 :
1047 : static GEN
1048 2512 : F2v_inflate(GEN v, GEN p, long n)
1049 : {
1050 2512 : long i, l = lg(p)-1;
1051 2512 : GEN w = zero_F2v(n);
1052 3944198 : for (i=1; i<=l; i++)
1053 3941686 : if (F2v_coeff(v,i))
1054 1967731 : F2v_set(w, p[i]);
1055 2512 : return w;
1056 : }
1057 :
1058 : static GEN
1059 46 : F2m_inflate(GEN x, GEN p, long n)
1060 2558 : { pari_APPLY_same(F2v_inflate(gel(x,i), p, n)) }
1061 :
1062 : GEN
1063 3369 : F2Ms_ker(GEN M, long nbrow)
1064 : {
1065 3369 : pari_sp av = avma;
1066 3369 : long nbcol = lg(M)-1;
1067 : GEN Mp, R, Rp, p;
1068 3369 : if (nbrow <= 640)
1069 3323 : return gc_upto(av, F2m_ker(F2Ms_to_F2m(M, nbrow)));
1070 46 : p = F2Ms_colelim(M, nbrow);
1071 46 : Mp = vecpermute(M, p);
1072 : do
1073 : {
1074 46 : R = block_lanczos(Mp, nbrow);
1075 46 : } while(!R);
1076 46 : Rp = F2m_inflate(R, p, nbcol);
1077 46 : return gc_GEN(av, Rp);
1078 : }
1079 :
1080 : GEN
1081 0 : F2m_to_F2Ms(GEN M)
1082 : {
1083 0 : long ncol = lg(M)-1;
1084 0 : GEN B = cgetg(ncol+1, t_MAT);
1085 : long i, j, k;
1086 0 : for(i = 1; i <= ncol; i++)
1087 : {
1088 0 : GEN D, V = gel(M,i);
1089 0 : long n = F2v_hamming(V), l = V[1];
1090 0 : D = cgetg(n+1, t_VECSMALL);
1091 0 : for (j=1, k=1; j<=l; j++)
1092 0 : if( F2v_coeff(V,j))
1093 0 : D[k++] = j;
1094 0 : gel(B, i) = D;
1095 : }
1096 0 : return B;
1097 : }
1098 :
1099 : GEN
1100 3323 : F2Ms_to_F2m(GEN M, long nrow)
1101 : {
1102 3323 : long i, j, l = lg(M);
1103 3323 : GEN B = cgetg(l, t_MAT);
1104 308467 : for(i = 1; i < l; i++)
1105 : {
1106 305144 : GEN Bi = zero_F2v(nrow), Mi = gel(M,i);
1107 305144 : long l = lg(Mi);
1108 3630379 : for (j = 1; j < l; j++)
1109 3325235 : F2v_set(Bi, Mi[j]);
1110 305144 : gel(B, i) = Bi;
1111 : }
1112 3323 : return B;
1113 : }
|
