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 :
32 : where the a_i are bits.
33 : */
34 :
35 : int
36 4816 : F2v_equal0(GEN V)
37 : {
38 4816 : long l = lg(V);
39 5982 : while (--l > 1)
40 5422 : if (V[l]) return 0;
41 560 : return 1;
42 : }
43 :
44 : GEN
45 3665649 : F2c_to_ZC(GEN x)
46 : {
47 3665649 : long l = x[1]+1, lx = lg(x);
48 3665649 : GEN z = cgetg(l, t_COL);
49 : long i, j, k;
50 7331360 : for (i=2, k=1; i<lx; i++)
51 22206135 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
52 18540410 : gel(z,k) = (x[i]&(1UL<<j))? gen_1: gen_0;
53 3665635 : return z;
54 : }
55 : GEN
56 3115 : F2c_to_mod(GEN x)
57 : {
58 3115 : long l = x[1]+1, lx = lg(x);
59 3115 : GEN z = cgetg(l, t_COL);
60 3115 : GEN _0 = mkintmod(gen_0,gen_2);
61 3115 : GEN _1 = mkintmod(gen_1,gen_2);
62 : long i, j, k;
63 15830 : for (i=2, k=1; i<lx; i++)
64 573625 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
65 560910 : gel(z,k) = (x[i]&(1UL<<j))? _1: _0;
66 3115 : return z;
67 : }
68 :
69 : /* x[a..b], a <= b */
70 : GEN
71 28 : F2v_slice(GEN x, long a, long b)
72 : {
73 28 : long i,j,k, l = b-a+1;
74 28 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
75 28 : z[1] = l;
76 98 : for(i=a,k=1,j=BITS_IN_LONG; i<=b; i++,j++)
77 : {
78 70 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
79 70 : if (F2v_coeff(x,i)) z[k] |= 1UL<<j;
80 : }
81 28 : return z;
82 : }
83 : /* x[a..b,], a <= b */
84 : GEN
85 14 : F2m_rowslice(GEN x, long a, long b)
86 42 : { pari_APPLY_same(F2v_slice(gel(x,i),a,b)) }
87 :
88 : GEN
89 3556 : F2v_to_Flv(GEN x)
90 : {
91 3556 : long l = x[1]+1, lx = lg(x);
92 3556 : GEN z = cgetg(l, t_VECSMALL);
93 : long i,j,k;
94 7136 : for (i=2, k=1; i<lx; i++)
95 65887 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
96 62307 : z[k] = (x[i]>>j)&1UL;
97 3556 : return z;
98 : }
99 :
100 : GEN
101 4486503 : F2m_to_ZM(GEN x) { pari_APPLY_same(F2c_to_ZC(gel(x,i))) }
102 : GEN
103 3234 : F2m_to_mod(GEN x) { pari_APPLY_same(F2c_to_mod(gel(x,i))) }
104 : GEN
105 1512 : F2m_to_Flm(GEN x) { pari_APPLY_same(F2v_to_Flv(gel(x,i))) }
106 :
107 : GEN
108 10236999 : ZV_to_F2v(GEN x)
109 : {
110 10236999 : long l = lg(x)-1;
111 10236999 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
112 : long i,j,k;
113 10236637 : z[1] = l;
114 73247050 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
115 : {
116 63010202 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
117 63010202 : if (mpodd(gel(x,i))) z[k] |= 1UL<<j;
118 : }
119 10236848 : return z;
120 : }
121 :
122 : GEN
123 7798 : RgV_to_F2v(GEN x)
124 : {
125 7798 : long l = lg(x)-1;
126 7798 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
127 : long i,j,k;
128 7798 : z[1] = l;
129 1410206 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
130 : {
131 1402408 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
132 1402408 : if (Rg_to_F2(gel(x,i))) z[k] |= 1UL<<j;
133 : }
134 7798 : return z;
135 : }
136 :
137 : GEN
138 5929 : Flv_to_F2v(GEN x)
139 : {
140 5929 : long l = lg(x)-1;
141 5929 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
142 : long i,j,k;
143 5929 : z[1] = l;
144 60235 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
145 : {
146 54306 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
147 54306 : if (x[i]&1L) z[k] |= 1UL<<j;
148 : }
149 5929 : return z;
150 : }
151 :
152 : GEN
153 12909979 : ZM_to_F2m(GEN x) { pari_APPLY_same(ZV_to_F2v(gel(x,i))) }
154 : GEN
155 8113 : RgM_to_F2m(GEN x) { pari_APPLY_same(RgV_to_F2v(gel(x,i))) }
156 : GEN
157 0 : Flm_to_F2m(GEN x) { pari_APPLY_same(Flv_to_F2v(gel(x,i))) }
158 :
159 : GEN
160 1665031 : const_F2v(long m)
161 : {
162 1665031 : long i, l = nbits2lg(m);
163 1665026 : GEN c = cgetg(l, t_VECSMALL);
164 1665003 : c[1] = m;
165 3339006 : for (i = 2; i < l; i++) uel(c,i) = -1UL;
166 1665003 : if (remsBIL(m)) uel(c,l-1) = (1UL<<remsBIL(m))-1UL;
167 1664996 : return c;
168 : }
169 :
170 : /* Allow lg(y)<lg(x) */
171 : void
172 23453150 : F2v_add_inplace(GEN x, GEN y)
173 : {
174 23453150 : long n = lg(y);
175 23453150 : long r = (n-2)&7L, q = n-r, i;
176 32928934 : for (i = 2; i < q; i += 8)
177 : {
178 9475784 : x[ i] ^= y[ i]; x[1+i] ^= y[1+i]; x[2+i] ^= y[2+i]; x[3+i] ^= y[3+i];
179 9475784 : x[4+i] ^= y[4+i]; x[5+i] ^= y[5+i]; x[6+i] ^= y[6+i]; x[7+i] ^= y[7+i];
180 : }
181 23453150 : switch (r)
182 : {
183 1768814 : case 7: x[i] ^= y[i]; i++; case 6: x[i] ^= y[i]; i++;
184 3533358 : case 5: x[i] ^= y[i]; i++; case 4: x[i] ^= y[i]; i++;
185 8515511 : case 3: x[i] ^= y[i]; i++; case 2: x[i] ^= y[i]; i++;
186 20673632 : case 1: x[i] ^= y[i]; i++;
187 : }
188 23453150 : }
189 :
190 : /* Allow lg(y)<lg(x) */
191 : void
192 14448 : F2v_and_inplace(GEN x, GEN y)
193 : {
194 14448 : long n = lg(y);
195 14448 : long r = (n-2)&7L, q = n-r, i;
196 14448 : for (i = 2; i < q; i += 8)
197 : {
198 0 : x[ i] &= y[ i]; x[1+i] &= y[1+i]; x[2+i] &= y[2+i]; x[3+i] &= y[3+i];
199 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];
200 : }
201 14448 : switch (r)
202 : {
203 0 : case 7: x[i] &= y[i]; i++; case 6: x[i] &= y[i]; i++;
204 0 : case 5: x[i] &= y[i]; i++; case 4: x[i] &= y[i]; i++;
205 2064 : case 3: x[i] &= y[i]; i++; case 2: x[i] &= y[i]; i++;
206 14448 : case 1: x[i] &= y[i]; i++;
207 : }
208 14448 : }
209 :
210 : /* Allow lg(y)<lg(x) */
211 : void
212 0 : F2v_or_inplace(GEN x, GEN y)
213 : {
214 0 : long n = lg(y);
215 0 : long r = (n-2)&7L, q = n-r, i;
216 0 : for (i = 2; i < q; i += 8)
217 : {
218 0 : x[ i] |= y[ i]; x[1+i] |= y[1+i]; x[2+i] |= y[2+i]; x[3+i] |= y[3+i];
219 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];
220 : }
221 0 : switch (r)
222 : {
223 0 : case 7: x[i] |= y[i]; i++; case 6: x[i] |= y[i]; i++;
224 0 : case 5: x[i] |= y[i]; i++; case 4: x[i] |= y[i]; i++;
225 0 : case 3: x[i] |= y[i]; i++; case 2: x[i] |= y[i]; i++;
226 0 : case 1: x[i] |= y[i]; i++;
227 : }
228 0 : }
229 :
230 : /* Allow lg(y)<lg(x) */
231 : void
232 546 : F2v_negimply_inplace(GEN x, GEN y)
233 : {
234 546 : long n = lg(y);
235 546 : long r = (n-2)&7L, q = n-r, i;
236 546 : for (i = 2; i < q; i += 8)
237 : {
238 0 : x[ i] &= ~y[ i]; x[1+i] &= ~y[1+i]; x[2+i] &= ~y[2+i]; x[3+i] &= ~y[3+i];
239 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];
240 : }
241 546 : switch (r)
242 : {
243 0 : case 7: x[i] &= ~y[i]; i++; case 6: x[i] &= ~y[i]; i++;
244 0 : case 5: x[i] &= ~y[i]; i++; case 4: x[i] &= ~y[i]; i++;
245 78 : case 3: x[i] &= ~y[i]; i++; case 2: x[i] &= ~y[i]; i++;
246 546 : case 1: x[i] &= ~y[i]; i++;
247 : }
248 546 : }
249 :
250 : ulong
251 0 : F2v_dotproduct(GEN x, GEN y)
252 : {
253 0 : long i, lx = lg(x);
254 : ulong c;
255 0 : if (lx <= 2) return 0;
256 0 : c = uel(x,2) & uel(y,2);
257 0 : for (i=3; i<lx; i++) c ^= uel(x,i) & uel(y,i);
258 : #ifdef LONG_IS_64BIT
259 0 : c ^= c >> 32;
260 : #endif
261 0 : c ^= c >> 16;
262 0 : c ^= c >> 8;
263 0 : c ^= c >> 4;
264 0 : c ^= c >> 2;
265 0 : c ^= c >> 1;
266 0 : return c & 1;
267 : }
268 :
269 : ulong
270 406 : F2v_hamming(GEN H)
271 : {
272 406 : long i, n=0, l=lg(H);
273 870 : for (i=2; i<l; i++) n += hammingl(uel(H,i));
274 406 : return n;
275 : }
276 :
277 : int
278 3171 : F2v_subset(GEN x, GEN y)
279 : {
280 3171 : long i, n = lg(y);
281 3738 : for (i = 2; i < n; i ++)
282 3374 : if ((x[i] & y[i]) != x[i]) return 0;
283 364 : return 1;
284 : }
285 :
286 : GEN
287 227341 : matid_F2m(long n)
288 : {
289 227341 : GEN y = cgetg(n+1,t_MAT);
290 : long i;
291 227340 : if (n < 0) pari_err_DOMAIN("matid_F2m", "dimension","<",gen_0,stoi(n));
292 985009 : for (i=1; i<=n; i++) { gel(y,i) = zero_F2v(n); F2v_set(gel(y,i),i); }
293 227345 : return y;
294 : }
295 :
296 : GEN
297 0 : F2m_row(GEN x, long j)
298 : {
299 0 : long i, l = lg(x);
300 0 : GEN V = zero_F2v(l-1);
301 0 : for(i = 1; i < l; i++)
302 0 : if (F2m_coeff(x,j,i))
303 0 : F2v_set(V,i);
304 0 : return V;
305 : }
306 :
307 : GEN
308 0 : F2m_transpose(GEN x)
309 : {
310 0 : long i, dx, lx = lg(x);
311 : GEN y;
312 0 : if (lx == 1) return cgetg(1,t_MAT);
313 0 : dx = coeff(x,1,1); y = cgetg(dx+1,t_MAT);
314 0 : for (i=1; i<=dx; i++) gel(y,i) = F2m_row(x,i);
315 0 : return y;
316 : }
317 :
318 : INLINE GEN
319 1899790 : F2m_F2c_mul_i(GEN x, GEN y, long lx, long l)
320 : {
321 : long j;
322 1899790 : GEN z = NULL;
323 :
324 13613093 : for (j=1; j<lx; j++)
325 : {
326 11713311 : if (!F2v_coeff(y,j)) continue;
327 3139410 : if (!z) z = vecsmall_copy(gel(x,j));
328 1508874 : else F2v_add_inplace(z,gel(x,j));
329 : }
330 1899782 : if (!z) z = zero_F2v(l);
331 1899784 : return z;
332 : }
333 :
334 : GEN
335 689083 : F2m_mul(GEN x, GEN y)
336 : {
337 689083 : long i,j,l,lx=lg(x), ly=lg(y);
338 : GEN z;
339 689083 : if (ly==1) return cgetg(1,t_MAT);
340 689083 : z = cgetg(ly,t_MAT);
341 689081 : if (lx==1)
342 : {
343 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
344 0 : return z;
345 : }
346 689081 : l = coeff(x,1,1);
347 2588865 : for (j=1; j<ly; j++) gel(z,j) = F2m_F2c_mul_i(x, gel(y,j), lx, l);
348 689078 : return z;
349 : }
350 :
351 : GEN
352 0 : F2m_F2c_mul(GEN x, GEN y)
353 : {
354 0 : long l, lx = lg(x);
355 0 : if (lx==1) return cgetg(1,t_VECSMALL);
356 0 : l = coeff(x,1,1);
357 0 : return F2m_F2c_mul_i(x, y, lx, l);
358 : }
359 :
360 : static GEN
361 0 : _F2m_mul(void *data, GEN x, GEN y) { (void) data; return F2m_mul(x,y); }
362 : static GEN
363 0 : _F2m_sqr(void *data, GEN x) { (void) data; return F2m_mul(x,x); }
364 : GEN
365 0 : F2m_powu(GEN x, ulong n)
366 : {
367 0 : if (!n) return matid(lg(x)-1);
368 0 : return gen_powu(x, n,NULL, &_F2m_sqr, &_F2m_mul);
369 : }
370 :
371 : static long
372 6352404 : F2v_find_nonzero(GEN x0, GEN mask0, long m)
373 : {
374 6352404 : ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
375 6352404 : long i, l = lg(x0)-2;
376 8914316 : for (i = 0; i < l; i++)
377 : {
378 6941322 : e = *x++ & *mask++;
379 6941322 : if (e) return i*BITS_IN_LONG+vals(e)+1;
380 : }
381 1972994 : return m+1;
382 : }
383 :
384 : /* in place, destroy x */
385 : GEN
386 1148889 : F2m_ker_sp(GEN x, long deplin)
387 : {
388 : GEN y, c, d;
389 : long i, j, k, r, m, n;
390 :
391 1148889 : n = lg(x)-1;
392 1148889 : m = mael(x,1,1); r=0;
393 :
394 1148889 : d = cgetg(n+1, t_VECSMALL);
395 1148883 : c = const_F2v(m);
396 5760417 : for (k=1; k<=n; k++)
397 : {
398 4916299 : GEN xk = gel(x,k);
399 4916299 : j = F2v_find_nonzero(xk, c, m);
400 4916319 : if (j>m)
401 : {
402 1467857 : if (deplin) {
403 304789 : GEN v = zero_F2v(n);
404 953653 : for (i=1; i<k; i++)
405 648865 : if (F2v_coeff(xk, d[i])) F2v_set(v, i);
406 304788 : F2v_set(v, k); return v;
407 : }
408 1163068 : r++; d[k] = 0;
409 : }
410 : else
411 : {
412 3448462 : F2v_clear(c,j); d[k] = j;
413 3448516 : F2v_clear(xk, j);
414 45280425 : for (i=k+1; i<=n; i++)
415 : {
416 41831884 : GEN xi = gel(x,i);
417 41831884 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
418 : }
419 3448541 : F2v_set(xk, j);
420 : }
421 : }
422 844118 : if (deplin) return NULL;
423 :
424 842683 : y = zero_F2m_copy(n,r);
425 2005763 : for (j=k=1; j<=r; j++,k++)
426 : {
427 3322255 : GEN C = gel(y,j); while (d[k]) k++;
428 8465826 : for (i=1; i<k; i++)
429 7302761 : if (d[i] && F2m_coeff(x,d[i],k)) F2v_set(C,i);
430 1163065 : F2v_set(C, k);
431 : }
432 842692 : return y;
433 : }
434 :
435 : /* not memory clean */
436 : GEN
437 199036 : F2m_ker(GEN x) { return F2m_ker_sp(F2m_copy(x), 0); }
438 : GEN
439 0 : F2m_deplin(GEN x) { return F2m_ker_sp(F2m_copy(x), 1); }
440 :
441 : ulong
442 1631 : F2m_det_sp(GEN x) { return !F2m_ker_sp(x, 1); }
443 :
444 : ulong
445 0 : F2m_det(GEN x)
446 : {
447 0 : pari_sp av = avma;
448 0 : ulong d = F2m_det_sp(F2m_copy(x));
449 0 : return gc_ulong(av, d);
450 : }
451 :
452 : /* Destroy x */
453 : GEN
454 354835 : F2m_gauss_pivot(GEN x, long *rr)
455 : {
456 : GEN c, d;
457 : long i, j, k, r, m, n;
458 :
459 354835 : n = lg(x)-1; if (!n) { *rr=0; return NULL; }
460 354835 : m = mael(x,1,1); r=0;
461 :
462 354835 : d = cgetg(n+1, t_VECSMALL);
463 354835 : c = const_F2v(m);
464 1791028 : for (k=1; k<=n; k++)
465 : {
466 1436193 : GEN xk = gel(x,k);
467 1436193 : j = F2v_find_nonzero(xk, c, m);
468 1436195 : if (j>m) { r++; d[k] = 0; }
469 : else
470 : {
471 930899 : F2v_clear(c,j); d[k] = j;
472 6625238 : for (i=k+1; i<=n; i++)
473 : {
474 5694341 : GEN xi = gel(x,i);
475 5694341 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
476 : }
477 : }
478 : }
479 :
480 354835 : *rr = r; return gc_const((pari_sp)d, d);
481 : }
482 :
483 : long
484 63 : F2m_rank(GEN x)
485 : {
486 63 : pari_sp av = avma;
487 : long r;
488 63 : (void)F2m_gauss_pivot(F2m_copy(x),&r);
489 63 : return gc_long(av, lg(x)-1 - r);
490 : }
491 :
492 : static GEN
493 14 : F2m_inv_upper_1_ind(GEN A, long index)
494 : {
495 14 : pari_sp av = avma;
496 14 : long n = lg(A)-1, i = index, j;
497 14 : GEN u = const_vecsmall(n, 0);
498 14 : u[i] = 1;
499 21 : for (i--; i>0; i--)
500 : {
501 7 : ulong m = F2m_coeff(A,i,i+1) & uel(u,i+1); /* j = i+1 */
502 7 : for (j=i+2; j<=n; j++) m ^= F2m_coeff(A,i,j) & uel(u,j);
503 7 : u[i] = m & 1;
504 : }
505 14 : return gerepileuptoleaf(av, Flv_to_F2v(u));
506 : }
507 : static GEN
508 7 : F2m_inv_upper_1(GEN A)
509 : {
510 : long i, l;
511 7 : GEN B = cgetg_copy(A, &l);
512 21 : for (i = 1; i < l; i++) gel(B,i) = F2m_inv_upper_1_ind(A, i);
513 7 : return B;
514 : }
515 :
516 : static GEN
517 757638 : F2_get_col(GEN b, GEN d, long li, long aco)
518 : {
519 757638 : long i, l = nbits2lg(aco);
520 757638 : GEN u = cgetg(l, t_VECSMALL);
521 757651 : u[1] = aco;
522 5451763 : for (i = 1; i <= li; i++)
523 4694115 : if (d[i]) /* d[i] can still be 0 if li > aco */
524 : {
525 4625278 : if (F2v_coeff(b, i))
526 1527462 : F2v_set(u, d[i]);
527 : else
528 3097816 : F2v_clear(u, d[i]);
529 : }
530 757648 : return u;
531 : }
532 :
533 : /* destroy a, b */
534 : GEN
535 227376 : F2m_gauss_sp(GEN a, GEN b)
536 : {
537 227376 : long i, j, k, l, li, bco, aco = lg(a)-1;
538 : GEN u, d;
539 :
540 227376 : if (!aco) return cgetg(1,t_MAT);
541 227376 : li = gel(a,1)[1];
542 227376 : d = zero_Flv(li);
543 227377 : bco = lg(b)-1;
544 976370 : for (i=1; i<=aco; i++)
545 : {
546 749017 : GEN ai = vecsmall_copy(gel(a,i));
547 749013 : if (!d[i] && F2v_coeff(ai, i))
548 499252 : k = i;
549 : else
550 1132496 : for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
551 : /* found a pivot on row k */
552 749018 : if (k > li) return NULL;
553 748997 : d[k] = i;
554 :
555 : /* Clear k-th row but column-wise instead of line-wise */
556 : /* a_ij -= a_i1*a1j/a_11
557 : line-wise grouping: L_j -= a_1j/a_11*L_1
558 : column-wise: C_i -= a_i1/a_11*C_1
559 : */
560 748997 : F2v_clear(ai,k);
561 5357161 : for (l=1; l<=aco; l++)
562 : {
563 4608167 : GEN al = gel(a,l);
564 4608167 : if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
565 : }
566 5374533 : for (l=1; l<=bco; l++)
567 : {
568 4625540 : GEN bl = gel(b,l);
569 4625540 : if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
570 : }
571 : }
572 227353 : u = cgetg(bco+1,t_MAT);
573 984994 : for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
574 227359 : return u;
575 : }
576 :
577 : GEN
578 35 : F2m_gauss(GEN a, GEN b)
579 : {
580 35 : pari_sp av = avma;
581 35 : if (lg(a) == 1) return cgetg(1,t_MAT);
582 35 : return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), F2m_copy(b)));
583 : }
584 : GEN
585 14 : F2m_F2c_gauss(GEN a, GEN b)
586 : {
587 14 : pari_sp av = avma;
588 14 : GEN z = F2m_gauss(a, mkmat(b));
589 14 : if (!z) return gc_NULL(av);
590 7 : if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
591 7 : return gerepileuptoleaf(av, gel(z,1));
592 : }
593 :
594 : GEN
595 35 : F2m_inv(GEN a)
596 : {
597 35 : pari_sp av = avma;
598 35 : if (lg(a) == 1) return cgetg(1,t_MAT);
599 35 : return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), matid_F2m(gel(a,1)[1])));
600 : }
601 :
602 : GEN
603 7 : F2m_invimage_i(GEN A, GEN B)
604 : {
605 : GEN d, x, X, Y;
606 7 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
607 7 : x = F2m_ker_sp(shallowconcat(A, B), 0);
608 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
609 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
610 : * Y has at least nB columns and full rank */
611 7 : nY = lg(x)-1;
612 7 : if (nY < nB) return NULL;
613 :
614 : /* implicitly: Y = rowslice(x, nA+1, nA+nB), nB rows */
615 7 : d = cgetg(nB+1, t_VECSMALL);
616 21 : for (i = nB, j = nY; i >= 1; i--, j--)
617 : {
618 14 : for (; j>=1; j--)
619 14 : if (F2m_coeff(x,nA+i,j)) { d[i] = j; break; } /* Y[i,j] */
620 14 : if (!j) return NULL;
621 : }
622 7 : x = vecpermute(x, d);
623 :
624 7 : X = F2m_rowslice(x, 1, nA);
625 7 : Y = F2m_rowslice(x, nA+1, nA+nB);
626 7 : return F2m_mul(X, F2m_inv_upper_1(Y));
627 : }
628 : GEN
629 0 : F2m_invimage(GEN A, GEN B)
630 : {
631 0 : pari_sp av = avma;
632 0 : GEN X = F2m_invimage_i(A,B);
633 0 : if (!X) return gc_NULL(av);
634 0 : return gerepileupto(av, X);
635 : }
636 :
637 : GEN
638 197080 : F2m_F2c_invimage(GEN A, GEN y)
639 : {
640 197080 : pari_sp av = avma;
641 197080 : long i, l = lg(A);
642 : GEN M, x;
643 :
644 197080 : if (l==1) return NULL;
645 197080 : if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
646 197079 : M = cgetg(l+1,t_MAT);
647 1060687 : for (i=1; i<l; i++) gel(M,i) = gel(A,i);
648 197079 : gel(M,l) = y; M = F2m_ker(M);
649 197085 : i = lg(M)-1; if (!i) return gc_NULL(av);
650 :
651 197085 : x = gel(M,i);
652 197085 : if (!F2v_coeff(x,l)) return gc_NULL(av);
653 197082 : F2v_clear(x, x[1]); x[1]--; /* remove last coord */
654 197083 : return gerepileuptoleaf(av, x);
655 : }
656 :
657 : /* Block Lanczos algorithm for kernel of sparse matrix (F2Ms)
658 : Based on lanczos.cpp by Jason Papadopoulos
659 : <https://github.com/sagemath/FlintQS/blob/master/src/lanczos.cpp>
660 : Copyright Jason Papadopoulos 2006
661 : Released under the GNU General Public License v2 or later version.
662 : */
663 :
664 : /* F2Ms are vector of vecsmall which represents nonzero entries of columns
665 : * F2w are vecsmall whoses entries are columns of a n x BIL F2m
666 : * F2wB are F2w in the special case where dim = BIL.
667 : */
668 :
669 : #define BIL BITS_IN_LONG
670 :
671 : static GEN
672 212 : F2w_transpose_F2m(GEN x)
673 : {
674 212 : long i, j, l = lg(x)-1;
675 212 : GEN z = cgetg(BIL+1, t_MAT);
676 13140 : for (j = 1; j <= BIL; j++)
677 12928 : gel(z,j) = zero_F2v(l);
678 322968 : for (i = 1; i <= l; i++)
679 : {
680 322756 : ulong xi = uel(x,i);
681 20312388 : for(j = 1; j <= BIL; j++)
682 19989632 : if (xi&(1UL<<(j-1)))
683 5941281 : F2v_set(gel(z, j), i);
684 : }
685 212 : return z;
686 : }
687 :
688 : static GEN
689 7704 : F2wB_mul(GEN a, GEN b)
690 : {
691 : long i, j;
692 7704 : GEN c = cgetg(BIL+1, t_VECSMALL);
693 466584 : for (i = 1; i <= BIL; i++)
694 : {
695 458880 : ulong s = 0, ai = a[i];
696 23872948 : for (j = 1; ai; j++, ai>>=1)
697 23414068 : if (ai & 1)
698 11794101 : s ^= b[j];
699 458880 : c[i] = s;
700 : }
701 7704 : return c;
702 : }
703 :
704 : static void
705 5136 : precompute_F2w_F2wB(GEN x, GEN c)
706 : {
707 : ulong z, xk;
708 : ulong i, j, k, index;
709 5136 : x++; c++;
710 43376 : for (j = 0; j < (BIL>>3); j++)
711 : {
712 9827680 : for (i = 0; i < 256; i++)
713 : {
714 9789440 : k = 0;
715 9789440 : index = i;
716 9789440 : z = 0;
717 78353760 : while (index) {
718 68564320 : xk = x[k];
719 68564320 : if (index & 1)
720 39157760 : z ^= xk;
721 68564320 : index >>= 1;
722 68564320 : k++;
723 : }
724 9789440 : c[i] = z;
725 : }
726 38240 : x += 8; c += 256;
727 : }
728 5136 : }
729 :
730 : static void
731 5136 : F2w_F2wB_mul_add_inplace(GEN v, GEN x, GEN y)
732 : {
733 5136 : long i, n = lg(y)-1;
734 : ulong word;
735 5136 : GEN c = cgetg(1+(BIL<<5), t_VECSMALL);
736 5136 : precompute_F2w_F2wB(x, c);
737 5136 : c++;
738 8455408 : for (i = 1; i <= n; i++)
739 : {
740 8450272 : word = v[i];
741 16900544 : y[i] ^= c[ 0*256 + ((word>> 0) & 0xff) ]
742 8450272 : ^ c[ 1*256 + ((word>> 8) & 0xff) ]
743 8450272 : ^ c[ 2*256 + ((word>>16) & 0xff) ]
744 8450272 : ^ c[ 3*256 + ((word>>24) & 0xff) ]
745 : #ifdef LONG_IS_64BIT
746 7740800 : ^ c[ 4*256 + ((word>>32) & 0xff) ]
747 7740800 : ^ c[ 5*256 + ((word>>40) & 0xff) ]
748 7740800 : ^ c[ 6*256 + ((word>>48) & 0xff) ]
749 7740800 : ^ c[ 7*256 + ((word>>56) ) ]
750 : #endif
751 : ;
752 : }
753 5136 : }
754 :
755 : /* Return x*y~, which is a F2wB */
756 : static GEN
757 3911 : F2w_transmul(GEN x, GEN y)
758 : {
759 3911 : long i, j, n = lg(x)-1;
760 3911 : GEN z = zero_zv(BIL);
761 3911 : pari_sp av = avma;
762 3911 : GEN c = zero_zv(BIL<<5) + 1;
763 3911 : GEN xy = z + 1;
764 :
765 6422946 : for (i = 1; i <= n; i++)
766 : {
767 6419035 : ulong xi = x[i];
768 6419035 : ulong yi = y[i];
769 6419035 : c[ 0*256 + ( xi & 0xff) ] ^= yi;
770 6419035 : c[ 1*256 + ((xi >> 8) & 0xff) ] ^= yi;
771 6419035 : c[ 2*256 + ((xi >> 16) & 0xff) ] ^= yi;
772 6419035 : c[ 3*256 + ((xi >> 24) & 0xff) ] ^= yi;
773 : #ifdef LONG_IS_64BIT
774 5880178 : c[ 4*256 + ((xi >> 32) & 0xff) ] ^= yi;
775 5880178 : c[ 5*256 + ((xi >> 40) & 0xff) ] ^= yi;
776 5880178 : c[ 6*256 + ((xi >> 48) & 0xff) ] ^= yi;
777 5880178 : c[ 7*256 + ((xi >> 56) ) ] ^= yi;
778 : #endif
779 : }
780 35199 : for(i = 0; i < 8; i++)
781 : {
782 31288 : ulong a0 = 0, a1 = 0, a2 = 0, a3 = 0;
783 : #ifdef LONG_IS_64BIT
784 26960 : ulong a4 = 0, a5 = 0, a6 = 0, a7 = 0;
785 : #endif
786 8041016 : for (j = 0; j < 256; j++) {
787 8009728 : if ((j >> i) & 1) {
788 4004864 : a0 ^= c[0*256 + j];
789 4004864 : a1 ^= c[1*256 + j];
790 4004864 : a2 ^= c[2*256 + j];
791 4004864 : a3 ^= c[3*256 + j];
792 : #ifdef LONG_IS_64BIT
793 3450880 : a4 ^= c[4*256 + j];
794 3450880 : a5 ^= c[5*256 + j];
795 3450880 : a6 ^= c[6*256 + j];
796 3450880 : a7 ^= c[7*256 + j];
797 : #endif
798 : }
799 : }
800 31288 : xy[ 0] = a0; xy[ 8] = a1; xy[16] = a2; xy[24] = a3;
801 : #ifdef LONG_IS_64BIT
802 26960 : xy[32] = a4; xy[40] = a5; xy[48] = a6; xy[56] = a7;
803 : #endif
804 31288 : xy++;
805 : }
806 3911 : return gc_const(av, z);
807 : }
808 :
809 : static GEN
810 1287 : identity_F2wB(void)
811 : {
812 : long i;
813 1287 : GEN M = cgetg(BIL+1, t_VECSMALL);
814 77927 : for (i = 1; i <= BIL; i++)
815 76640 : uel(M,i) = 1UL<<(i-1);
816 1287 : return M;
817 : }
818 :
819 : static GEN
820 1287 : find_nonsingular_sub(GEN t, GEN last_s, GEN *pt_s)
821 : {
822 1287 : long i, j, dim = 0;
823 : ulong mask, row_i, row_j;
824 1287 : long last_dim = lg(last_s)-1;
825 1287 : GEN s = cgetg(BIL+1, t_VECSMALL);
826 1287 : GEN M1 = identity_F2wB();
827 1287 : pari_sp av = avma;
828 1287 : GEN cols = cgetg(BIL+1, t_VECSMALL);
829 1287 : GEN M0 = zv_copy(t);
830 :
831 1287 : mask = 0;
832 76891 : for (i = 1; i <= last_dim; i++)
833 : {
834 75604 : cols[BIL + 1 - i] = last_s[i];
835 75604 : mask |= 1UL<<(last_s[i]-1);
836 : }
837 77927 : for (i = j = 1; i <= BIL; i++)
838 76640 : if (!(mask & (1UL<<(i-1))))
839 1036 : cols[j++] = i;
840 :
841 : /* compute the inverse of t[][] */
842 :
843 77927 : for (i = 1; i <= BIL; i++)
844 : {
845 76640 : mask = 1UL<<(cols[i]-1);
846 76640 : row_i = cols[i];
847 153840 : for (j = i; j <= BIL; j++)
848 : {
849 151608 : row_j = cols[j];
850 151608 : if (uel(M0,row_j) & mask)
851 : {
852 74408 : swap(gel(M0, row_j), gel(M0, row_i));
853 74408 : swap(gel(M1, row_j), gel(M1, row_i));
854 74408 : break;
855 : }
856 : }
857 76640 : if (j <= BIL)
858 : {
859 4660328 : for (j = 1; j <= BIL; j++)
860 : {
861 4585920 : row_j = cols[j];
862 4585920 : if (row_i != row_j && (M0[row_j] & mask))
863 : {
864 2220765 : uel(M0,row_j) ^= uel(M0,row_i);
865 2220765 : uel(M1,row_j) ^= uel(M1,row_i);
866 : }
867 : }
868 74408 : s[++dim] = cols[i];
869 74408 : continue;
870 : }
871 2232 : for (j = i; j <= BIL; j++)
872 : {
873 2232 : row_j = cols[j];
874 2232 : if (uel(M1,row_j) & mask)
875 : {
876 2232 : swap(gel(M0, row_j), gel(M0, row_i));
877 2232 : swap(gel(M1, row_j), gel(M1, row_i));
878 2232 : break;
879 : }
880 : }
881 2232 : if (j > BIL) return NULL;
882 137976 : for (j = 1; j <= BIL; j++)
883 : {
884 135744 : row_j = cols[j];
885 135744 : if (row_i != row_j && (M1[row_j] & mask))
886 : {
887 0 : uel(M0,row_j) ^= uel(M0,row_i);
888 0 : uel(M1,row_j) ^= uel(M1,row_i);
889 : }
890 : }
891 2232 : M0[row_i] = M1[row_i] = 0;
892 : }
893 1287 : mask = 0;
894 75695 : for (i = 1; i <= dim; i++)
895 74408 : mask |= 1UL<<(s[i]-1);
896 76891 : for (i = 1; i <= last_dim; i++)
897 75604 : mask |= 1UL<<(last_s[i]-1);
898 1287 : if (mask != (ulong)(-1))
899 3 : return NULL; /* Failure */
900 1284 : setlg(s, dim+1);
901 1284 : set_avma(av);
902 1284 : *pt_s = s;
903 1284 : return M1;
904 : }
905 :
906 : /* Compute x * A~ */
907 : static GEN
908 1502 : F2w_F2Ms_transmul(GEN x, GEN A, long nbrow)
909 : {
910 1502 : long i, j, l = lg(A);
911 1502 : GEN b = zero_zv(nbrow);
912 2423338 : for (i = 1; i < l; i++)
913 : {
914 2421836 : GEN c = gel(A,i);
915 2421836 : long lc = lg(c);
916 2421836 : ulong xi = x[i];
917 45327612 : for (j = 1; j < lc; j++)
918 42905776 : b[c[j]] ^= xi;
919 : }
920 1502 : return b;
921 : }
922 :
923 : /* Compute x * A */
924 : static GEN
925 1396 : F2w_F2Ms_mul(GEN x, GEN A)
926 : {
927 1396 : long i, j, l = lg(A);
928 1396 : GEN b = cgetg(l, t_VECSMALL);
929 2271274 : for (i = 1; i < l; i++)
930 : {
931 2269878 : GEN c = gel(A,i);
932 2269878 : long lc = lg(c);
933 2269878 : ulong s = 0;
934 42516676 : for (j = 1; j < lc; j++)
935 40246798 : s ^= x[c[j]];
936 2269878 : b[i] = s;
937 : }
938 1396 : return b;
939 : }
940 :
941 : static void
942 2568 : F2wB_addid_inplace(GEN f)
943 : {
944 : long i;
945 155528 : for (i = 1; i <= BIL; i++)
946 152960 : uel(f,i) ^= 1UL<<(i-1);
947 2568 : }
948 :
949 : static void
950 2568 : F2w_mask_inplace(GEN f, ulong m)
951 : {
952 2568 : long i, l = lg(f);
953 2191616 : for (i = 1; i < l; i++)
954 2189048 : uel(f,i) &= m;
955 2568 : }
956 :
957 : static GEN
958 56 : block_lanczos(GEN B, ulong nbrow)
959 : {
960 56 : pari_sp av = avma, av2;
961 : GEN v0, v1, v2, vnext, x, w;
962 : GEN winv0, winv1, winv2;
963 : GEN vt_a_v0, vt_a_v1, vt_a2_v0, vt_a2_v1;
964 : GEN d, e, f, f2, s0;
965 : long i, iter;
966 56 : long n = lg(B)-1;
967 : long dim0;
968 : ulong mask0, mask1;
969 56 : v1 = zero_zv(n);
970 56 : v2 = zero_zv(n);
971 56 : vt_a_v1 = zero_zv(BIL);
972 56 : vt_a2_v1 = zero_zv(BIL);
973 56 : winv1 = zero_zv(BIL);
974 56 : winv2 = zero_zv(BIL);
975 56 : s0 = identity_zv(BIL);
976 56 : mask1 = (ulong)(-1);
977 :
978 56 : x = random_zv(n);
979 56 : w = F2w_F2Ms_mul(F2w_F2Ms_transmul(x, B, nbrow), B);
980 56 : v0 = w;
981 56 : av2 = avma;
982 56 : for (iter=1;;iter++)
983 : {
984 1340 : vnext = F2w_F2Ms_mul(F2w_F2Ms_transmul(v0, B, nbrow), B);
985 1340 : vt_a_v0 = F2w_transmul(v0, vnext);
986 1340 : if (zv_equal0(vt_a_v0)) break; /* success */
987 1287 : vt_a2_v0 = F2w_transmul(vnext, vnext);
988 1287 : winv0 = find_nonsingular_sub(vt_a_v0, s0, &s0);
989 1287 : if (!winv0) return gc_NULL(av); /* failure */
990 1284 : dim0 = lg(s0)-1;
991 1284 : mask0 = 0;
992 75688 : for (i = 1; i <= dim0; i++)
993 74404 : mask0 |= 1UL<<(s0[i]-1);
994 1284 : d = cgetg(BIL+1, t_VECSMALL);
995 77764 : for (i = 1; i <= BIL; i++)
996 76480 : d[i] = (vt_a2_v0[i] & mask0) ^ vt_a_v0[i];
997 :
998 1284 : d = F2wB_mul(winv0, d);
999 1284 : F2wB_addid_inplace(d);
1000 1284 : e = F2wB_mul(winv1, vt_a_v0);
1001 1284 : F2w_mask_inplace(e, mask0);
1002 1284 : f = F2wB_mul(vt_a_v1, winv1);
1003 1284 : F2wB_addid_inplace(f);
1004 1284 : f = F2wB_mul(winv2, f);
1005 1284 : f2 = cgetg(BIL+1, t_VECSMALL);
1006 77764 : for (i = 1; i <= BIL; i++)
1007 76480 : f2[i] = ((vt_a2_v1[i] & mask1) ^ vt_a_v1[i]) & mask0;
1008 :
1009 1284 : f = F2wB_mul(f, f2);
1010 1284 : F2w_mask_inplace(vnext, mask0);
1011 1284 : F2w_F2wB_mul_add_inplace(v0, d, vnext);
1012 1284 : F2w_F2wB_mul_add_inplace(v1, e, vnext);
1013 1284 : F2w_F2wB_mul_add_inplace(v2, f, vnext);
1014 1284 : d = F2wB_mul(winv0, F2w_transmul(v0, w));
1015 1284 : F2w_F2wB_mul_add_inplace(v0, d, x);
1016 1284 : v2 = v1; v1 = v0; v0 = vnext;
1017 1284 : winv2 = winv1; winv1 = winv0;
1018 1284 : vt_a_v1 = vt_a_v0;
1019 1284 : vt_a2_v1 = vt_a2_v0;
1020 1284 : mask1 = mask0;
1021 1284 : gerepileall(av2, 9, &x, &s0, &v0, &v1, &v2,
1022 : &winv1, &winv2, &vt_a_v1, &vt_a2_v1);
1023 : }
1024 53 : if (DEBUGLEVEL >= 5)
1025 0 : err_printf("Lanczos halted after %ld iterations\n", iter);
1026 53 : v1 = F2w_F2Ms_transmul(x, B, nbrow);
1027 53 : v2 = F2w_F2Ms_transmul(v0, B, nbrow);
1028 53 : x = shallowconcat(F2w_transpose_F2m(x), F2w_transpose_F2m(v0));
1029 53 : v1 = shallowconcat(F2w_transpose_F2m(v1), F2w_transpose_F2m(v2));
1030 53 : s0 = gel(F2m_indexrank(x), 2);
1031 53 : x = shallowextract(x, s0);
1032 53 : v1 = shallowextract(v1, s0);
1033 53 : return F2m_mul(x, F2m_ker(v1));
1034 : }
1035 :
1036 : static GEN
1037 2892 : F2v_inflate(GEN v, GEN p, long n)
1038 : {
1039 2892 : long i, l = lg(p)-1;
1040 2892 : GEN w = zero_F2v(n);
1041 4341298 : for (i=1; i<=l; i++)
1042 4338406 : if (F2v_coeff(v,i))
1043 2166436 : F2v_set(w, p[i]);
1044 2892 : return w;
1045 : }
1046 :
1047 : static GEN
1048 53 : F2m_inflate(GEN x, GEN p, long n)
1049 2945 : { pari_APPLY_same(F2v_inflate(gel(x,i), p, n)) }
1050 :
1051 : GEN
1052 1094 : F2Ms_ker(GEN M, long nbrow)
1053 : {
1054 1094 : pari_sp av = avma;
1055 1094 : long nbcol = lg(M)-1;
1056 : GEN Mp, R, Rp, p;
1057 1094 : if (nbrow <= 640)
1058 1041 : return gerepileupto(av, F2m_ker(F2Ms_to_F2m(M, nbrow)));
1059 53 : p = F2Ms_colelim(M, nbrow);
1060 53 : Mp = vecpermute(M, p);
1061 : do
1062 : {
1063 56 : R = block_lanczos(Mp, nbrow);
1064 56 : } while(!R);
1065 53 : Rp = F2m_inflate(R, p, nbcol);
1066 53 : return gerepilecopy(av, Rp);
1067 : }
1068 :
1069 : GEN
1070 0 : F2m_to_F2Ms(GEN M)
1071 : {
1072 0 : long ncol = lg(M)-1;
1073 0 : GEN B = cgetg(ncol+1, t_MAT);
1074 : long i, j, k;
1075 0 : for(i = 1; i <= ncol; i++)
1076 : {
1077 0 : GEN D, V = gel(M,i);
1078 0 : long n = F2v_hamming(V), l = V[1];
1079 0 : D = cgetg(n+1, t_VECSMALL);
1080 0 : for (j=1, k=1; j<=l; j++)
1081 0 : if( F2v_coeff(V,j))
1082 0 : D[k++] = j;
1083 0 : gel(B, i) = D;
1084 : }
1085 0 : return B;
1086 : }
1087 :
1088 : GEN
1089 1041 : F2Ms_to_F2m(GEN M, long nrow)
1090 : {
1091 1041 : long i, j, l = lg(M);
1092 1041 : GEN B = cgetg(l, t_MAT);
1093 169790 : for(i = 1; i < l; i++)
1094 : {
1095 168749 : GEN Bi = zero_F2v(nrow), Mi = gel(M,i);
1096 168749 : long l = lg(Mi);
1097 2226857 : for (j = 1; j < l; j++)
1098 2058108 : F2v_set(Bi, Mi[j]);
1099 168749 : gel(B, i) = Bi;
1100 : }
1101 1041 : return B;
1102 : }
|
