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 3953004 : F2c_to_ZC(GEN x)
44 : {
45 3953004 : long l = x[1]+1, lx = lg(x);
46 3953004 : GEN z = cgetg(l, t_COL);
47 : long i, j, k;
48 7930282 : for (i=2, k=1; i<lx; i++)
49 32089561 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
50 28112283 : gel(z,k) = (x[i]&(1UL<<j))? gen_1: gen_0;
51 3953004 : return z;
52 : }
53 : GEN
54 4690 : F2c_to_mod(GEN x)
55 : {
56 4690 : long l = x[1]+1, lx = lg(x);
57 4690 : GEN z = cgetg(l, t_COL);
58 4690 : GEN _0 = mkintmod(gen_0,gen_2);
59 4690 : GEN _1 = mkintmod(gen_1,gen_2);
60 : long i, j, k;
61 19030 : for (i=2, k=1; i<lx; i++)
62 604825 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
63 590485 : gel(z,k) = (x[i]&(1UL<<j))? _1: _0;
64 4690 : 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 4818887 : F2m_to_ZM(GEN x) { pari_APPLY_same(F2c_to_ZC(gel(x,i))) }
100 : GEN
101 4963 : 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 10786994 : ZV_to_F2v(GEN x)
121 : {
122 10786994 : long i, j, k, l = lg(x)-1;
123 10786994 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
124 10786994 : z[1] = l;
125 99757366 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
126 : {
127 88970372 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
128 88970372 : if (mpodd(gel(x,i))) z[k] |= 1UL<<j;
129 : }
130 10786994 : return z;
131 : }
132 :
133 : GEN
134 9079 : RgV_to_F2v(GEN x)
135 : {
136 9079 : long l = lg(x)-1;
137 9079 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
138 : long i,j,k;
139 9079 : z[1] = l;
140 1441062 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
141 : {
142 1431983 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
143 1431983 : if (Rg_to_F2(gel(x,i))) z[k] |= 1UL<<j;
144 : }
145 9079 : return z;
146 : }
147 :
148 : GEN
149 6489 : Flv_to_F2v(GEN x)
150 : {
151 6489 : long l = lg(x)-1;
152 6489 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
153 : long i,j,k;
154 6489 : z[1] = l;
155 10235295 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
156 : {
157 10228806 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
158 10228806 : if (x[i]&1L) z[k] |= 1UL<<j;
159 : }
160 6489 : return z;
161 : }
162 :
163 : GEN
164 13534764 : ZM_to_F2m(GEN x) { pari_APPLY_same(ZV_to_F2v(gel(x,i))) }
165 : GEN
166 9548 : 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 1666556 : const_F2v(long m)
172 : {
173 1666556 : long i, l = nbits2lg(m);
174 1666556 : GEN c = cgetg(l, t_VECSMALL);
175 1666556 : c[1] = m;
176 3346904 : for (i = 2; i < l; i++) uel(c,i) = -1UL;
177 1666556 : if (remsBIL(m)) uel(c,l-1) = (1UL<<remsBIL(m))-1UL;
178 1666556 : return c;
179 : }
180 :
181 : /* Allow lg(y)<lg(x) */
182 : void
183 27364870 : F2v_add_inplace(GEN x, GEN y)
184 : {
185 27364870 : long n = lg(y);
186 27364870 : long r = (n-2)&7L, q = n-r, i;
187 37379828 : for (i = 2; i < q; i += 8)
188 : {
189 10014958 : x[ i] ^= y[ i]; x[1+i] ^= y[1+i]; x[2+i] ^= y[2+i]; x[3+i] ^= y[3+i];
190 10014958 : 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 27364870 : switch (r)
193 : {
194 1786441 : case 7: x[i] ^= y[i]; i++; case 6: x[i] ^= y[i]; i++;
195 3495114 : case 5: x[i] ^= y[i]; i++; case 4: x[i] ^= y[i]; i++;
196 10134846 : case 3: x[i] ^= y[i]; i++; case 2: x[i] ^= y[i]; i++;
197 24591589 : case 1: x[i] ^= y[i]; i++;
198 : }
199 27364870 : }
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 1771 : F2v_negimply_inplace(GEN x, GEN y)
244 : {
245 1771 : long n = lg(y);
246 1771 : long r = (n-2)&7L, q = n-r, i;
247 92531 : 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 1771 : 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 301 : case 3: x[i] &= ~y[i]; i++; case 2: x[i] &= ~y[i]; i++;
257 1771 : case 1: x[i] &= ~y[i]; i++;
258 : }
259 1771 : }
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 0 : return thuemorseu(c);
270 : }
271 :
272 : ulong
273 2681 : F2v_hamming(GEN H)
274 : {
275 2681 : long i, n=0, l=lg(H);
276 1095157 : for (i=2; i<l; i++) n += hammingu(uel(H,i));
277 2681 : return n;
278 : }
279 :
280 : int
281 9450 : F2v_subset(GEN x, GEN y)
282 : {
283 9450 : long i, n = lg(y);
284 11233 : for (i = 2; i < n; i ++)
285 10064 : if ((x[i] & y[i]) != x[i]) return 0;
286 1169 : return 1;
287 : }
288 :
289 : GEN
290 230605 : matid_F2m(long n)
291 : {
292 230605 : GEN y = cgetg(n+1,t_MAT);
293 : long i;
294 230605 : if (n < 0) pari_err_DOMAIN("matid_F2m", "dimension","<",gen_0,stoi(n));
295 990018 : for (i=1; i<=n; i++) { gel(y,i) = zero_F2v(n); F2v_set(gel(y,i),i); }
296 230605 : return y;
297 : }
298 :
299 : GEN
300 0 : F2m_row(GEN x, long j)
301 : {
302 0 : long i, l = lg(x);
303 0 : GEN V = zero_F2v(l-1);
304 0 : for(i = 1; i < l; i++)
305 0 : if (F2m_coeff(x,j,i))
306 0 : F2v_set(V,i);
307 0 : return V;
308 : }
309 :
310 : GEN
311 0 : F2m_transpose(GEN x)
312 : {
313 0 : long i, dx, lx = lg(x);
314 : GEN y;
315 0 : if (lx == 1) return cgetg(1,t_MAT);
316 0 : dx = coeff(x,1,1); y = cgetg(dx+1,t_MAT);
317 0 : for (i=1; i<=dx; i++) gel(y,i) = F2m_row(x,i);
318 0 : return y;
319 : }
320 :
321 : INLINE GEN
322 2195609 : F2m_F2c_mul_i(GEN x, GEN y, long lx, long l)
323 : {
324 : long j;
325 2195609 : GEN z = NULL;
326 :
327 23711664 : for (j=1; j<lx; j++)
328 : {
329 21516055 : if (!F2v_coeff(y,j)) continue;
330 5200747 : if (!z) z = vecsmall_copy(gel(x,j));
331 3273505 : else F2v_add_inplace(z,gel(x,j));
332 : }
333 2195609 : if (!z) z = zero_F2v(l);
334 2195609 : return z;
335 : }
336 :
337 : GEN
338 723926 : F2m_mul(GEN x, GEN y)
339 : {
340 723926 : long i,j,l,lx=lg(x), ly=lg(y);
341 : GEN z;
342 723926 : if (ly==1) return cgetg(1,t_MAT);
343 723926 : z = cgetg(ly,t_MAT);
344 723926 : if (lx==1)
345 : {
346 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
347 0 : return z;
348 : }
349 723926 : l = coeff(x,1,1);
350 2919535 : for (j=1; j<ly; j++) gel(z,j) = F2m_F2c_mul_i(x, gel(y,j), lx, l);
351 723926 : return z;
352 : }
353 :
354 : GEN
355 0 : F2m_F2c_mul(GEN x, GEN y)
356 : {
357 0 : long l, lx = lg(x);
358 0 : if (lx==1) return cgetg(1,t_VECSMALL);
359 0 : l = coeff(x,1,1);
360 0 : return F2m_F2c_mul_i(x, y, lx, l);
361 : }
362 :
363 : static GEN
364 0 : _F2m_mul(void *data, GEN x, GEN y) { (void) data; return F2m_mul(x,y); }
365 : static GEN
366 0 : _F2m_sqr(void *data, GEN x) { (void) data; return F2m_mul(x,x); }
367 : GEN
368 0 : F2m_powu(GEN x, ulong n)
369 : {
370 0 : if (!n) return matid(lg(x)-1);
371 0 : return gen_powu(x, n,NULL, &_F2m_sqr, &_F2m_mul);
372 : }
373 :
374 : static long
375 6349946 : F2v_find_nonzero(GEN x0, GEN mask0, long m)
376 : {
377 6349946 : ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
378 6349946 : long i, l = lg(x0)-2;
379 8967118 : for (i = 0; i < l; i++)
380 : {
381 7013960 : e = *x++ & *mask++;
382 7013960 : if (e) return i*BITS_IN_LONG+vals(e)+1;
383 : }
384 1953158 : return m+1;
385 : }
386 :
387 : /* in place, destroy x */
388 : GEN
389 1150946 : F2m_ker_sp(GEN x, long deplin)
390 : {
391 : GEN y, c, d;
392 : long i, j, k, r, m, n;
393 :
394 1150946 : n = lg(x)-1; if (n==0) return x;
395 1150946 : m = mael(x,1,1); r=0;
396 :
397 1150946 : d = cgetg(n+1, t_VECSMALL);
398 1150946 : c = const_F2v(m);
399 5827511 : for (k=1; k<=n; k++)
400 : {
401 4978997 : GEN xk = gel(x,k);
402 4978997 : j = F2v_find_nonzero(xk, c, m);
403 4978997 : if (j>m)
404 : {
405 1478472 : if (deplin) {
406 302432 : GEN v = zero_F2v(n);
407 921381 : for (i=1; i<k; i++)
408 618949 : if (F2v_coeff(xk, d[i])) F2v_set(v, i);
409 302432 : F2v_set(v, k); return v;
410 : }
411 1176040 : r++; d[k] = 0;
412 : }
413 : else
414 : {
415 3500525 : F2v_clear(c,j); d[k] = j;
416 3500525 : F2v_clear(xk, j);
417 50038601 : for (i=k+1; i<=n; i++)
418 : {
419 46538076 : GEN xi = gel(x,i);
420 46538076 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
421 : }
422 3500525 : F2v_set(xk, j);
423 : }
424 : }
425 848514 : if (deplin) return NULL;
426 :
427 846946 : y = zero_F2m_copy(n,r);
428 2022986 : for (j=k=1; j<=r; j++,k++)
429 : {
430 3407516 : GEN C = gel(y,j); while (d[k]) k++;
431 9649773 : for (i=1; i<k; i++)
432 8473733 : if (d[i] && F2m_coeff(x,d[i],k)) F2v_set(C,i);
433 1176040 : F2v_set(C, k);
434 : }
435 846946 : return y;
436 : }
437 :
438 : /* not memory clean */
439 : GEN
440 201003 : F2m_ker(GEN x) { return F2m_ker_sp(F2m_copy(x), 0); }
441 : GEN
442 0 : F2m_deplin(GEN x) { return F2m_ker_sp(F2m_copy(x), 1); }
443 :
444 : ulong
445 1764 : F2m_det_sp(GEN x) { return !F2m_ker_sp(x, 1); }
446 :
447 : ulong
448 0 : F2m_det(GEN x)
449 : {
450 0 : pari_sp av = avma;
451 0 : ulong d = F2m_det_sp(F2m_copy(x));
452 0 : return gc_ulong(av, d);
453 : }
454 :
455 : /* Destroy x */
456 : GEN
457 356448 : F2m_gauss_pivot(GEN x, long *rr)
458 : {
459 : GEN c, d;
460 : long i, j, k, r, m, n;
461 :
462 356448 : n = lg(x)-1; if (!n) { *rr=0; return NULL; }
463 356448 : m = mael(x,1,1); r=0;
464 :
465 356448 : d = cgetg(n+1, t_VECSMALL);
466 356448 : c = const_F2v(m);
467 1727397 : for (k=1; k<=n; k++)
468 : {
469 1370949 : GEN xk = gel(x,k);
470 1370949 : j = F2v_find_nonzero(xk, c, m);
471 1370949 : if (j>m) { r++; d[k] = 0; }
472 : else
473 : {
474 896263 : F2v_clear(c,j); d[k] = j;
475 6362263 : for (i=k+1; i<=n; i++)
476 : {
477 5466000 : GEN xi = gel(x,i);
478 5466000 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
479 : }
480 : }
481 : }
482 :
483 356448 : *rr = r; return gc_const((pari_sp)d, d);
484 : }
485 :
486 : long
487 63 : F2m_rank(GEN x)
488 : {
489 63 : pari_sp av = avma;
490 : long r;
491 63 : (void)F2m_gauss_pivot(F2m_copy(x),&r);
492 63 : return gc_long(av, lg(x)-1 - r);
493 : }
494 :
495 : static GEN
496 14 : F2m_inv_upper_1_ind(GEN A, long index)
497 : {
498 14 : pari_sp av = avma;
499 14 : long n = lg(A)-1, i = index, j;
500 14 : GEN u = const_vecsmall(n, 0);
501 14 : u[i] = 1;
502 21 : for (i--; i>0; i--)
503 : {
504 7 : ulong m = F2m_coeff(A,i,i+1) & uel(u,i+1); /* j = i+1 */
505 7 : for (j=i+2; j<=n; j++) m ^= F2m_coeff(A,i,j) & uel(u,j);
506 7 : u[i] = m & 1;
507 : }
508 14 : return gc_leaf(av, Flv_to_F2v(u));
509 : }
510 : static GEN
511 7 : F2m_inv_upper_1(GEN A)
512 : {
513 : long i, l;
514 7 : GEN B = cgetg_copy(A, &l);
515 21 : for (i = 1; i < l; i++) gel(B,i) = F2m_inv_upper_1_ind(A, i);
516 7 : return B;
517 : }
518 :
519 : static GEN
520 759392 : F2_get_col(GEN b, GEN d, long li, long aco)
521 : {
522 759392 : long i, l = nbits2lg(aco);
523 759392 : GEN u = cgetg(l, t_VECSMALL);
524 759392 : u[1] = aco;
525 5649264 : for (i = 1; i <= li; i++)
526 4889872 : if (d[i]) /* d[i] can still be 0 if li > aco */
527 : {
528 4736473 : if (F2v_coeff(b, i))
529 1564875 : F2v_set(u, d[i]);
530 : else
531 3171598 : F2v_clear(u, d[i]);
532 : }
533 759392 : return u;
534 : }
535 :
536 : /* destroy a, b */
537 : GEN
538 230640 : F2m_gauss_sp(GEN a, GEN b)
539 : {
540 230640 : long i, j, k, l, li, bco, aco = lg(a)-1;
541 : GEN u, d;
542 :
543 230640 : if (!aco) return cgetg(1,t_MAT);
544 230640 : li = gel(a,1)[1];
545 230640 : d = zero_Flv(li);
546 230640 : bco = lg(b)-1;
547 978431 : for (i=1; i<=aco; i++)
548 : {
549 747812 : GEN ai = vecsmall_copy(gel(a,i));
550 747812 : if (!d[i] && F2v_coeff(ai, i))
551 501342 : k = i;
552 : else
553 1187093 : for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
554 : /* found a pivot on row k */
555 747812 : if (k > li) return NULL;
556 747791 : d[k] = i;
557 :
558 : /* Clear k-th row but column-wise instead of line-wise */
559 : /* a_ij -= a_i1*a1j/a_11
560 : line-wise grouping: L_j -= a_1j/a_11*L_1
561 : column-wise: C_i -= a_i1/a_11*C_1
562 : */
563 747791 : F2v_clear(ai,k);
564 5456449 : for (l=1; l<=aco; l++)
565 : {
566 4708658 : GEN al = gel(a,l);
567 4708658 : if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
568 : }
569 5484523 : for (l=1; l<=bco; l++)
570 : {
571 4736732 : GEN bl = gel(b,l);
572 4736732 : if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
573 : }
574 : }
575 230619 : u = cgetg(bco+1,t_MAT);
576 990011 : for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
577 230619 : return u;
578 : }
579 :
580 : GEN
581 35 : F2m_gauss(GEN a, GEN b)
582 : {
583 35 : pari_sp av = avma;
584 35 : if (lg(a) == 1) return cgetg(1,t_MAT);
585 35 : return gc_upto(av, F2m_gauss_sp(F2m_copy(a), F2m_copy(b)));
586 : }
587 : GEN
588 14 : F2m_F2c_gauss(GEN a, GEN b)
589 : {
590 14 : pari_sp av = avma;
591 14 : GEN z = F2m_gauss(a, mkmat(b));
592 14 : if (!z) return gc_NULL(av);
593 7 : if (lg(z) == 1) retgc_const(av, cgetg(1, t_VECSMALL));
594 7 : return gc_leaf(av, gel(z,1));
595 : }
596 :
597 : GEN
598 189 : F2m_inv(GEN a)
599 : {
600 189 : pari_sp av = avma;
601 189 : if (lg(a) == 1) return cgetg(1,t_MAT);
602 189 : return gc_upto(av, F2m_gauss_sp(F2m_copy(a), matid_F2m(gel(a,1)[1])));
603 : }
604 :
605 : GEN
606 7 : F2m_invimage_i(GEN A, GEN B)
607 : {
608 : GEN d, x, X, Y;
609 7 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
610 7 : x = F2m_ker_sp(shallowconcat(A, B), 0);
611 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
612 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
613 : * Y has at least nB columns and full rank */
614 7 : nY = lg(x)-1;
615 7 : if (nY < nB) return NULL;
616 :
617 : /* implicitly: Y = rowslice(x, nA+1, nA+nB), nB rows */
618 7 : d = cgetg(nB+1, t_VECSMALL);
619 21 : for (i = nB, j = nY; i >= 1; i--, j--)
620 : {
621 14 : for (; j>=1; j--)
622 14 : if (F2m_coeff(x,nA+i,j)) { d[i] = j; break; } /* Y[i,j] */
623 14 : if (!j) return NULL;
624 : }
625 7 : x = vecpermute(x, d);
626 :
627 7 : X = F2m_rowslice(x, 1, nA);
628 7 : Y = F2m_rowslice(x, nA+1, nA+nB);
629 7 : return F2m_mul(X, F2m_inv_upper_1(Y));
630 : }
631 : GEN
632 0 : F2m_invimage(GEN A, GEN B)
633 : {
634 0 : pari_sp av = avma;
635 0 : GEN X = F2m_invimage_i(A,B);
636 0 : if (!X) return gc_NULL(av);
637 0 : return gc_upto(av, X);
638 : }
639 :
640 : GEN
641 196681 : F2m_F2c_invimage(GEN A, GEN y)
642 : {
643 196681 : pari_sp av = avma;
644 196681 : long i, l = lg(A);
645 : GEN M, x;
646 :
647 196681 : if (l==1) return NULL;
648 196681 : if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
649 196681 : M = cgetg(l+1,t_MAT);
650 1022686 : for (i=1; i<l; i++) gel(M,i) = gel(A,i);
651 196681 : gel(M,l) = y; M = F2m_ker(M);
652 196681 : i = lg(M)-1; if (!i) return gc_NULL(av);
653 :
654 196681 : x = gel(M,i);
655 196681 : if (!F2v_coeff(x,l)) return gc_NULL(av);
656 196681 : F2v_clear(x, x[1]); x[1]--; /* remove last coord */
657 196681 : return gc_leaf(av, x);
658 : }
659 :
660 : /* Block Lanczos algorithm for kernel of sparse matrix (F2Ms)
661 : Based on lanczos.cpp by Jason Papadopoulos
662 : <https://github.com/sagemath/FlintQS/blob/master/src/lanczos.cpp>
663 : Copyright Jason Papadopoulos 2006
664 : Released under the GNU General Public License v2 or later version.
665 : */
666 :
667 : /* F2Ms are vector of vecsmall which represents nonzero entries of columns
668 : * F2w are vecsmall whoses entries are columns of a n x BIL F2m
669 : * F2wB are F2w in the special case where dim = BIL.
670 : */
671 :
672 : #define BIL BITS_IN_LONG
673 :
674 : static GEN
675 232 : F2w_transpose_F2m(GEN x)
676 : {
677 232 : long i, j, l = lg(x)-1;
678 232 : GEN z = cgetg(BIL+1, t_MAT);
679 14568 : for (j = 1; j <= BIL; j++)
680 14336 : gel(z,j) = zero_F2v(l);
681 398714 : for (i = 1; i <= l; i++)
682 : {
683 398482 : ulong xi = uel(x,i);
684 25378386 : for(j = 1; j <= BIL; j++)
685 24979904 : if (xi&(1UL<<(j-1)))
686 7141106 : F2v_set(gel(z, j), i);
687 : }
688 232 : return z;
689 : }
690 :
691 : static GEN
692 9030 : F2wB_mul(GEN a, GEN b)
693 : {
694 : long i, j;
695 9030 : GEN c = cgetg(BIL+1, t_VECSMALL);
696 564294 : for (i = 1; i <= BIL; i++)
697 : {
698 555264 : ulong s = 0, ai = a[i];
699 29236175 : for (j = 1; ai; j++, ai>>=1)
700 28680911 : if (ai & 1)
701 14451765 : s ^= b[j];
702 555264 : c[i] = s;
703 : }
704 9030 : return c;
705 : }
706 :
707 : static void
708 6020 : precompute_F2w_F2wB(GEN x, GEN c)
709 : {
710 : ulong z, xk;
711 : ulong i, j, k, index;
712 6020 : x++; c++;
713 52292 : for (j = 0; j < (BIL>>3); j++)
714 : {
715 11891904 : for (i = 0; i < 256; i++)
716 : {
717 11845632 : k = 0;
718 11845632 : index = i;
719 11845632 : z = 0;
720 94811328 : while (index) {
721 82965696 : xk = x[k];
722 82965696 : if (index & 1)
723 47382528 : z ^= xk;
724 82965696 : index >>= 1;
725 82965696 : k++;
726 : }
727 11845632 : c[i] = z;
728 : }
729 46272 : x += 8; c += 256;
730 : }
731 6020 : }
732 :
733 : static void
734 6020 : F2w_F2wB_mul_add_inplace(GEN v, GEN x, GEN y)
735 : {
736 6020 : long i, n = lg(y)-1;
737 : ulong word;
738 6020 : GEN c = cgetg(1+(BIL<<5), t_VECSMALL);
739 6020 : precompute_F2w_F2wB(x, c);
740 6020 : c++;
741 11095216 : for (i = 1; i <= n; i++)
742 : {
743 11089196 : word = v[i];
744 11089196 : y[i] ^= c[ 0*256 + ((word>> 0) & 0xff) ]
745 11089196 : ^ c[ 1*256 + ((word>> 8) & 0xff) ]
746 11089196 : ^ c[ 2*256 + ((word>>16) & 0xff) ]
747 11089196 : ^ c[ 3*256 + ((word>>24) & 0xff) ]
748 : #ifdef LONG_IS_64BIT
749 10616620 : ^ c[ 4*256 + ((word>>32) & 0xff) ]
750 10616620 : ^ c[ 5*256 + ((word>>40) & 0xff) ]
751 10616620 : ^ c[ 6*256 + ((word>>48) & 0xff) ]
752 10616620 : ^ c[ 7*256 + ((word>>56) ) ]
753 : #endif
754 : ;
755 : }
756 6020 : }
757 :
758 : /* Return x*y~, which is a F2wB */
759 : static GEN
760 4575 : F2w_transmul(GEN x, GEN y)
761 : {
762 4575 : long i, j, n = lg(x)-1;
763 4575 : GEN z = zero_zv(BIL);
764 4575 : pari_sp av = avma;
765 4575 : GEN c = zero_zv(BIL<<5) + 1;
766 4575 : GEN xy = z + 1;
767 :
768 8417271 : for (i = 1; i <= n; i++)
769 : {
770 8412696 : ulong xi = x[i];
771 8412696 : ulong yi = y[i];
772 8412696 : c[ 0*256 + ( xi & 0xff) ] ^= yi;
773 8412696 : c[ 1*256 + ((xi >> 8) & 0xff) ] ^= yi;
774 8412696 : c[ 2*256 + ((xi >> 16) & 0xff) ] ^= yi;
775 8412696 : c[ 3*256 + ((xi >> 24) & 0xff) ] ^= yi;
776 : #ifdef LONG_IS_64BIT
777 8054399 : c[ 4*256 + ((xi >> 32) & 0xff) ] ^= yi;
778 8054399 : c[ 5*256 + ((xi >> 40) & 0xff) ] ^= yi;
779 8054399 : c[ 6*256 + ((xi >> 48) & 0xff) ] ^= yi;
780 8054399 : c[ 7*256 + ((xi >> 56) ) ] ^= yi;
781 : #endif
782 : }
783 41175 : for(i = 0; i < 8; i++)
784 : {
785 36600 : ulong a0 = 0, a1 = 0, a2 = 0, a3 = 0;
786 : #ifdef LONG_IS_64BIT
787 33736 : ulong a4 = 0, a5 = 0, a6 = 0, a7 = 0;
788 : #endif
789 9406200 : for (j = 0; j < 256; j++) {
790 9369600 : if ((j >> i) & 1) {
791 4684800 : a0 ^= c[0*256 + j];
792 4684800 : a1 ^= c[1*256 + j];
793 4684800 : a2 ^= c[2*256 + j];
794 4684800 : a3 ^= c[3*256 + j];
795 : #ifdef LONG_IS_64BIT
796 4318208 : a4 ^= c[4*256 + j];
797 4318208 : a5 ^= c[5*256 + j];
798 4318208 : a6 ^= c[6*256 + j];
799 4318208 : a7 ^= c[7*256 + j];
800 : #endif
801 : }
802 : }
803 36600 : xy[ 0] = a0; xy[ 8] = a1; xy[16] = a2; xy[24] = a3;
804 : #ifdef LONG_IS_64BIT
805 33736 : xy[32] = a4; xy[40] = a5; xy[48] = a6; xy[56] = a7;
806 : #endif
807 36600 : xy++;
808 : }
809 4575 : return gc_const(av, z);
810 : }
811 :
812 : static GEN
813 1506 : identity_F2wB(void)
814 : {
815 : long i;
816 1506 : GEN M = cgetg(BIL+1, t_VECSMALL);
817 94114 : for (i = 1; i <= BIL; i++)
818 92608 : uel(M,i) = 1UL<<(i-1);
819 1506 : return M;
820 : }
821 :
822 : static GEN
823 1506 : find_nonsingular_sub(GEN t, GEN last_s, GEN *pt_s)
824 : {
825 1506 : long i, j, dim = 0;
826 : ulong mask, row_i, row_j;
827 1506 : long last_dim = lg(last_s)-1;
828 1506 : GEN s = cgetg(BIL+1, t_VECSMALL);
829 1506 : GEN M1 = identity_F2wB();
830 1506 : pari_sp av = avma;
831 1506 : GEN cols = cgetg(BIL+1, t_VECSMALL);
832 1506 : GEN M0 = zv_copy(t);
833 :
834 1506 : mask = 0;
835 93092 : for (i = 1; i <= last_dim; i++)
836 : {
837 91586 : cols[BIL + 1 - i] = last_s[i];
838 91586 : mask |= 1UL<<(last_s[i]-1);
839 : }
840 94114 : for (i = j = 1; i <= BIL; i++)
841 92608 : if (!(mask & (1UL<<(i-1))))
842 1022 : cols[j++] = i;
843 :
844 : /* compute the inverse of t[][] */
845 :
846 94114 : for (i = 1; i <= BIL; i++)
847 : {
848 92608 : mask = 1UL<<(cols[i]-1);
849 92608 : row_i = cols[i];
850 187213 : for (j = i; j <= BIL; j++)
851 : {
852 184720 : row_j = cols[j];
853 184720 : if (uel(M0,row_j) & mask)
854 : {
855 90115 : swap(gel(M0, row_j), gel(M0, row_i));
856 90115 : swap(gel(M1, row_j), gel(M1, row_i));
857 90115 : break;
858 : }
859 : }
860 92608 : if (j <= BIL)
861 : {
862 5740675 : for (j = 1; j <= BIL; j++)
863 : {
864 5650560 : row_j = cols[j];
865 5650560 : if (row_i != row_j && (M0[row_j] & mask))
866 : {
867 2736838 : uel(M0,row_j) ^= uel(M0,row_i);
868 2736838 : uel(M1,row_j) ^= uel(M1,row_i);
869 : }
870 : }
871 90115 : s[++dim] = cols[i];
872 90115 : continue;
873 : }
874 2493 : for (j = i; j <= BIL; j++)
875 : {
876 2493 : row_j = cols[j];
877 2493 : if (uel(M1,row_j) & mask)
878 : {
879 2493 : swap(gel(M0, row_j), gel(M0, row_i));
880 2493 : swap(gel(M1, row_j), gel(M1, row_i));
881 2493 : break;
882 : }
883 : }
884 2493 : if (j > BIL) return NULL;
885 158013 : for (j = 1; j <= BIL; j++)
886 : {
887 155520 : row_j = cols[j];
888 155520 : if (row_i != row_j && (M1[row_j] & mask))
889 : {
890 0 : uel(M0,row_j) ^= uel(M0,row_i);
891 0 : uel(M1,row_j) ^= uel(M1,row_i);
892 : }
893 : }
894 2493 : M0[row_i] = M1[row_i] = 0;
895 : }
896 1506 : mask = 0;
897 91621 : for (i = 1; i <= dim; i++)
898 90115 : mask |= 1UL<<(s[i]-1);
899 93092 : for (i = 1; i <= last_dim; i++)
900 91586 : mask |= 1UL<<(last_s[i]-1);
901 1506 : if (mask != (ulong)(-1))
902 1 : return NULL; /* Failure */
903 1505 : setlg(s, dim+1);
904 1505 : set_avma(av);
905 1505 : *pt_s = s;
906 1505 : return M1;
907 : }
908 :
909 : /* Compute x * A~ */
910 : static GEN
911 1739 : F2w_F2Ms_transmul(GEN x, GEN A, long nbrow)
912 : {
913 1739 : long i, j, l = lg(A);
914 1739 : GEN b = zero_zv(nbrow);
915 3151972 : for (i = 1; i < l; i++)
916 : {
917 3150233 : GEN c = gel(A,i);
918 3150233 : long lc = lg(c);
919 3150233 : ulong xi = x[i];
920 59364045 : for (j = 1; j < lc; j++)
921 56213812 : b[c[j]] ^= xi;
922 : }
923 1739 : return b;
924 : }
925 :
926 : /* Compute x * A */
927 : static GEN
928 1623 : F2w_F2Ms_mul(GEN x, GEN A)
929 : {
930 1623 : long i, j, l = lg(A);
931 1623 : GEN b = cgetg(l, t_VECSMALL);
932 2963766 : for (i = 1; i < l; i++)
933 : {
934 2962143 : GEN c = gel(A,i);
935 2962143 : long lc = lg(c);
936 2962143 : ulong s = 0;
937 55857469 : for (j = 1; j < lc; j++)
938 52895326 : s ^= x[c[j]];
939 2962143 : b[i] = s;
940 : }
941 1623 : return b;
942 : }
943 :
944 : static void
945 3010 : F2wB_addid_inplace(GEN f)
946 : {
947 : long i;
948 188098 : for (i = 1; i <= BIL; i++)
949 185088 : uel(f,i) ^= 1UL<<(i-1);
950 3010 : }
951 :
952 : static void
953 3010 : F2w_mask_inplace(GEN f, ulong m)
954 : {
955 3010 : long i, l = lg(f);
956 2867853 : for (i = 1; i < l; i++)
957 2864843 : uel(f,i) &= m;
958 3010 : }
959 :
960 : static GEN
961 59 : block_lanczos(GEN B, ulong nbrow)
962 : {
963 59 : pari_sp av = avma, av2;
964 : GEN v0, v1, v2, vnext, x, w;
965 : GEN winv0, winv1, winv2;
966 : GEN vt_a_v0, vt_a_v1, vt_a2_v0, vt_a2_v1;
967 : GEN d, e, f, f2, s0;
968 : long i, iter;
969 59 : long n = lg(B)-1;
970 : long dim0;
971 : ulong mask0, mask1;
972 59 : v1 = zero_zv(n);
973 59 : v2 = zero_zv(n);
974 59 : vt_a_v1 = zero_zv(BIL);
975 59 : vt_a2_v1 = zero_zv(BIL);
976 59 : winv1 = zero_zv(BIL);
977 59 : winv2 = zero_zv(BIL);
978 59 : s0 = identity_zv(BIL);
979 59 : mask1 = (ulong)(-1);
980 :
981 59 : x = random_zv(n);
982 59 : w = F2w_F2Ms_mul(F2w_F2Ms_transmul(x, B, nbrow), B);
983 59 : v0 = w;
984 59 : av2 = avma;
985 59 : for (iter=1;;iter++)
986 : {
987 1564 : vnext = F2w_F2Ms_mul(F2w_F2Ms_transmul(v0, B, nbrow), B);
988 1564 : vt_a_v0 = F2w_transmul(v0, vnext);
989 1564 : if (zv_equal0(vt_a_v0)) break; /* success */
990 1506 : vt_a2_v0 = F2w_transmul(vnext, vnext);
991 1506 : winv0 = find_nonsingular_sub(vt_a_v0, s0, &s0);
992 1506 : if (!winv0) return gc_NULL(av); /* failure */
993 1505 : dim0 = lg(s0)-1;
994 1505 : mask0 = 0;
995 91618 : for (i = 1; i <= dim0; i++)
996 90113 : mask0 |= 1UL<<(s0[i]-1);
997 1505 : d = cgetg(BIL+1, t_VECSMALL);
998 94049 : for (i = 1; i <= BIL; i++)
999 92544 : d[i] = (vt_a2_v0[i] & mask0) ^ vt_a_v0[i];
1000 :
1001 1505 : d = F2wB_mul(winv0, d);
1002 1505 : F2wB_addid_inplace(d);
1003 1505 : e = F2wB_mul(winv1, vt_a_v0);
1004 1505 : F2w_mask_inplace(e, mask0);
1005 1505 : f = F2wB_mul(vt_a_v1, winv1);
1006 1505 : F2wB_addid_inplace(f);
1007 1505 : f = F2wB_mul(winv2, f);
1008 1505 : f2 = cgetg(BIL+1, t_VECSMALL);
1009 94049 : for (i = 1; i <= BIL; i++)
1010 92544 : f2[i] = ((vt_a2_v1[i] & mask1) ^ vt_a_v1[i]) & mask0;
1011 :
1012 1505 : f = F2wB_mul(f, f2);
1013 1505 : F2w_mask_inplace(vnext, mask0);
1014 1505 : F2w_F2wB_mul_add_inplace(v0, d, vnext);
1015 1505 : F2w_F2wB_mul_add_inplace(v1, e, vnext);
1016 1505 : F2w_F2wB_mul_add_inplace(v2, f, vnext);
1017 1505 : d = F2wB_mul(winv0, F2w_transmul(v0, w));
1018 1505 : F2w_F2wB_mul_add_inplace(v0, d, x);
1019 1505 : v2 = v1; v1 = v0; v0 = vnext;
1020 1505 : winv2 = winv1; winv1 = winv0;
1021 1505 : vt_a_v1 = vt_a_v0;
1022 1505 : vt_a2_v1 = vt_a2_v0;
1023 1505 : mask1 = mask0;
1024 1505 : (void)gc_all(av2, 9, &x, &s0, &v0, &v1, &v2,
1025 : &winv1, &winv2, &vt_a_v1, &vt_a2_v1);
1026 : }
1027 58 : if (DEBUGLEVEL >= 5)
1028 0 : err_printf("Lanczos halted after %ld iterations\n", iter);
1029 58 : v1 = F2w_F2Ms_transmul(x, B, nbrow);
1030 58 : v2 = F2w_F2Ms_transmul(v0, B, nbrow);
1031 58 : x = shallowconcat(F2w_transpose_F2m(x), F2w_transpose_F2m(v0));
1032 58 : v1 = shallowconcat(F2w_transpose_F2m(v1), F2w_transpose_F2m(v2));
1033 58 : s0 = gel(F2m_indexrank(x), 2);
1034 58 : x = shallowextract(x, s0);
1035 58 : v1 = shallowextract(v1, s0);
1036 58 : return F2m_mul(x, F2m_ker(v1));
1037 : }
1038 :
1039 : static GEN
1040 3250 : F2v_inflate(GEN v, GEN p, long n)
1041 : {
1042 3250 : long i, l = lg(p)-1;
1043 3250 : GEN w = zero_F2v(n);
1044 5508992 : for (i=1; i<=l; i++)
1045 5505742 : if (F2v_coeff(v,i))
1046 2751001 : F2v_set(w, p[i]);
1047 3250 : return w;
1048 : }
1049 :
1050 : static GEN
1051 58 : F2m_inflate(GEN x, GEN p, long n)
1052 3308 : { pari_APPLY_same(F2v_inflate(gel(x,i), p, n)) }
1053 :
1054 : GEN
1055 3394 : F2Ms_ker(GEN M, long nbrow)
1056 : {
1057 3394 : pari_sp av = avma;
1058 3394 : long nbcol = lg(M)-1;
1059 : GEN Mp, R, Rp, p;
1060 3394 : if (nbrow <= 640)
1061 3336 : return gc_upto(av, F2m_ker(F2Ms_to_F2m(M, nbrow)));
1062 58 : p = F2Ms_colelim(M, nbrow);
1063 58 : Mp = vecpermute(M, p);
1064 : do
1065 : {
1066 59 : R = block_lanczos(Mp, nbrow);
1067 59 : } while(!R);
1068 58 : Rp = F2m_inflate(R, p, nbcol);
1069 58 : return gc_GEN(av, Rp);
1070 : }
1071 :
1072 : GEN
1073 0 : F2m_to_F2Ms(GEN M)
1074 : {
1075 0 : long ncol = lg(M)-1;
1076 0 : GEN B = cgetg(ncol+1, t_MAT);
1077 : long i, j, k;
1078 0 : for(i = 1; i <= ncol; i++)
1079 : {
1080 0 : GEN D, V = gel(M,i);
1081 0 : long n = F2v_hamming(V), l = V[1];
1082 0 : D = cgetg(n+1, t_VECSMALL);
1083 0 : for (j=1, k=1; j<=l; j++)
1084 0 : if( F2v_coeff(V,j))
1085 0 : D[k++] = j;
1086 0 : gel(B, i) = D;
1087 : }
1088 0 : return B;
1089 : }
1090 :
1091 : GEN
1092 3336 : F2Ms_to_F2m(GEN M, long nrow)
1093 : {
1094 3336 : long i, j, l = lg(M);
1095 3336 : GEN B = cgetg(l, t_MAT);
1096 312643 : for(i = 1; i < l; i++)
1097 : {
1098 309307 : GEN Bi = zero_F2v(nrow), Mi = gel(M,i);
1099 309307 : long l = lg(Mi);
1100 3696308 : for (j = 1; j < l; j++)
1101 3387001 : F2v_set(Bi, Mi[j]);
1102 309307 : gel(B, i) = Bi;
1103 : }
1104 3336 : return B;
1105 : }
|
