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 3478517 : F2c_to_ZC(GEN x)
46 : {
47 3478517 : long l = x[1]+1, lx = lg(x);
48 3478517 : GEN z = cgetg(l, t_COL);
49 : long i, j, k;
50 6956881 : for (i=2, k=1; i<lx; i++)
51 20001753 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
52 16523274 : gel(z,k) = (x[i]&(1UL<<j))? gen_1: gen_0;
53 3478402 : 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 3654 : F2v_to_Flv(GEN x)
90 : {
91 3654 : long l = x[1]+1, lx = lg(x);
92 3654 : GEN z = cgetg(l, t_VECSMALL);
93 : long i,j,k;
94 7332 : for (i=2, k=1; i<lx; i++)
95 66685 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
96 63007 : z[k] = (x[i]>>j)&1UL;
97 3654 : return z;
98 : }
99 :
100 : GEN
101 4290300 : 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 1568 : F2m_to_Flm(GEN x) { pari_APPLY_same(F2v_to_Flv(gel(x,i))) }
106 :
107 : GEN
108 1218 : RgV_F2v_extract_shallow(GEN V, GEN x)
109 : {
110 1218 : long n = F2v_hamming(x), m = 1;
111 1218 : long l = x[1]+1, lx = lg(x);
112 1218 : GEN W = cgetg(n+1, t_VEC);
113 : long i,j,k;
114 2436 : for (i=2, k=1; i<lx; i++)
115 15764 : for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
116 14546 : if ((x[i]>>j)&1UL)
117 2499 : gel(W, m++) = gel(V,k);
118 1218 : return W;
119 : }
120 :
121 : GEN
122 9592866 : ZV_to_F2v(GEN x)
123 : {
124 9592866 : long i, j, k, l = lg(x)-1;
125 9592866 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
126 9592691 : z[1] = l;
127 66193508 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
128 : {
129 56600746 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
130 56600746 : if (mpodd(gel(x,i))) z[k] |= 1UL<<j;
131 : }
132 9592762 : return z;
133 : }
134 :
135 : GEN
136 7798 : RgV_to_F2v(GEN x)
137 : {
138 7798 : long l = lg(x)-1;
139 7798 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
140 : long i,j,k;
141 7798 : z[1] = l;
142 1410206 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
143 : {
144 1402408 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
145 1402408 : if (Rg_to_F2(gel(x,i))) z[k] |= 1UL<<j;
146 : }
147 7798 : return z;
148 : }
149 :
150 : GEN
151 6027 : Flv_to_F2v(GEN x)
152 : {
153 6027 : long l = lg(x)-1;
154 6027 : GEN z = cgetg(nbits2lg(l), t_VECSMALL);
155 : long i,j,k;
156 6027 : z[1] = l;
157 10230689 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
158 : {
159 10224662 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
160 10224662 : if (x[i]&1L) z[k] |= 1UL<<j;
161 : }
162 6027 : return z;
163 : }
164 :
165 : GEN
166 12212034 : ZM_to_F2m(GEN x) { pari_APPLY_same(ZV_to_F2v(gel(x,i))) }
167 : GEN
168 8113 : RgM_to_F2m(GEN x) { pari_APPLY_same(RgV_to_F2v(gel(x,i))) }
169 : GEN
170 0 : Flm_to_F2m(GEN x) { pari_APPLY_same(Flv_to_F2v(gel(x,i))) }
171 :
172 : GEN
173 1634497 : const_F2v(long m)
174 : {
175 1634497 : long i, l = nbits2lg(m);
176 1634490 : GEN c = cgetg(l, t_VECSMALL);
177 1634438 : c[1] = m;
178 3278355 : for (i = 2; i < l; i++) uel(c,i) = -1UL;
179 1634438 : if (remsBIL(m)) uel(c,l-1) = (1UL<<remsBIL(m))-1UL;
180 1634420 : return c;
181 : }
182 :
183 : /* Allow lg(y)<lg(x) */
184 : void
185 23804937 : F2v_add_inplace(GEN x, GEN y)
186 : {
187 23804937 : long n = lg(y);
188 23804937 : long r = (n-2)&7L, q = n-r, i;
189 33242468 : for (i = 2; i < q; i += 8)
190 : {
191 9437531 : x[ i] ^= y[ i]; x[1+i] ^= y[1+i]; x[2+i] ^= y[2+i]; x[3+i] ^= y[3+i];
192 9437531 : x[4+i] ^= y[4+i]; x[5+i] ^= y[5+i]; x[6+i] ^= y[6+i]; x[7+i] ^= y[7+i];
193 : }
194 23804937 : switch (r)
195 : {
196 1782356 : case 7: x[i] ^= y[i]; i++; case 6: x[i] ^= y[i]; i++;
197 3541588 : case 5: x[i] ^= y[i]; i++; case 4: x[i] ^= y[i]; i++;
198 9171932 : case 3: x[i] ^= y[i]; i++; case 2: x[i] ^= y[i]; i++;
199 21026636 : case 1: x[i] ^= y[i]; i++;
200 : }
201 23804937 : }
202 :
203 : /* Allow lg(y)<lg(x) */
204 : void
205 14448 : F2v_and_inplace(GEN x, GEN y)
206 : {
207 14448 : long n = lg(y);
208 14448 : long r = (n-2)&7L, q = n-r, i;
209 14448 : for (i = 2; i < q; i += 8)
210 : {
211 0 : x[ i] &= y[ i]; x[1+i] &= y[1+i]; x[2+i] &= y[2+i]; x[3+i] &= y[3+i];
212 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];
213 : }
214 14448 : switch (r)
215 : {
216 0 : case 7: x[i] &= y[i]; i++; case 6: x[i] &= y[i]; i++;
217 0 : case 5: x[i] &= y[i]; i++; case 4: x[i] &= y[i]; i++;
218 2064 : case 3: x[i] &= y[i]; i++; case 2: x[i] &= y[i]; i++;
219 14448 : case 1: x[i] &= y[i]; i++;
220 : }
221 14448 : }
222 :
223 : /* Allow lg(y)<lg(x) */
224 : void
225 0 : F2v_or_inplace(GEN x, GEN y)
226 : {
227 0 : long n = lg(y);
228 0 : long r = (n-2)&7L, q = n-r, i;
229 0 : for (i = 2; i < q; i += 8)
230 : {
231 0 : x[ i] |= y[ i]; x[1+i] |= y[1+i]; x[2+i] |= y[2+i]; x[3+i] |= y[3+i];
232 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];
233 : }
234 0 : switch (r)
235 : {
236 0 : case 7: x[i] |= y[i]; i++; case 6: x[i] |= y[i]; i++;
237 0 : case 5: x[i] |= y[i]; i++; case 4: x[i] |= y[i]; i++;
238 0 : case 3: x[i] |= y[i]; i++; case 2: x[i] |= y[i]; i++;
239 0 : case 1: x[i] |= y[i]; i++;
240 : }
241 0 : }
242 :
243 : /* Allow lg(y)<lg(x) */
244 : void
245 1148 : F2v_negimply_inplace(GEN x, GEN y)
246 : {
247 1148 : long n = lg(y);
248 1148 : long r = (n-2)&7L, q = n-r, i;
249 91908 : for (i = 2; i < q; i += 8)
250 : {
251 90760 : x[ i] &= ~y[ i]; x[1+i] &= ~y[1+i]; x[2+i] &= ~y[2+i]; x[3+i] &= ~y[3+i];
252 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];
253 : }
254 1148 : switch (r)
255 : {
256 56 : case 7: x[i] &= ~y[i]; i++; case 6: x[i] &= ~y[i]; i++;
257 56 : case 5: x[i] &= ~y[i]; i++; case 4: x[i] &= ~y[i]; i++;
258 212 : case 3: x[i] &= ~y[i]; i++; case 2: x[i] &= ~y[i]; i++;
259 1148 : case 1: x[i] &= ~y[i]; i++;
260 : }
261 1148 : }
262 :
263 : ulong
264 0 : F2v_dotproduct(GEN x, GEN y)
265 : {
266 0 : long i, lx = lg(x);
267 : ulong c;
268 0 : if (lx <= 2) return 0;
269 0 : c = uel(x,2) & uel(y,2);
270 0 : for (i=3; i<lx; i++) c ^= uel(x,i) & uel(y,i);
271 : #ifdef LONG_IS_64BIT
272 0 : c ^= c >> 32;
273 : #endif
274 0 : c ^= c >> 16;
275 0 : c ^= c >> 8;
276 0 : c ^= c >> 4;
277 0 : c ^= c >> 2;
278 0 : c ^= c >> 1;
279 0 : return c & 1;
280 : }
281 :
282 : ulong
283 2156 : F2v_hamming(GEN H)
284 : {
285 2156 : long i, n=0, l=lg(H);
286 1094046 : for (i=2; i<l; i++) n += hammingl(uel(H,i));
287 2156 : return n;
288 : }
289 :
290 : int
291 6342 : F2v_subset(GEN x, GEN y)
292 : {
293 6342 : long i, n = lg(y);
294 7476 : for (i = 2; i < n; i ++)
295 6748 : if ((x[i] & y[i]) != x[i]) return 0;
296 728 : return 1;
297 : }
298 :
299 : GEN
300 224269 : matid_F2m(long n)
301 : {
302 224269 : GEN y = cgetg(n+1,t_MAT);
303 : long i;
304 224269 : if (n < 0) pari_err_DOMAIN("matid_F2m", "dimension","<",gen_0,stoi(n));
305 940416 : for (i=1; i<=n; i++) { gel(y,i) = zero_F2v(n); F2v_set(gel(y,i),i); }
306 224269 : return y;
307 : }
308 :
309 : GEN
310 0 : F2m_row(GEN x, long j)
311 : {
312 0 : long i, l = lg(x);
313 0 : GEN V = zero_F2v(l-1);
314 0 : for(i = 1; i < l; i++)
315 0 : if (F2m_coeff(x,j,i))
316 0 : F2v_set(V,i);
317 0 : return V;
318 : }
319 :
320 : GEN
321 0 : F2m_transpose(GEN x)
322 : {
323 0 : long i, dx, lx = lg(x);
324 : GEN y;
325 0 : if (lx == 1) return cgetg(1,t_MAT);
326 0 : dx = coeff(x,1,1); y = cgetg(dx+1,t_MAT);
327 0 : for (i=1; i<=dx; i++) gel(y,i) = F2m_row(x,i);
328 0 : return y;
329 : }
330 :
331 : INLINE GEN
332 1778345 : F2m_F2c_mul_i(GEN x, GEN y, long lx, long l)
333 : {
334 : long j;
335 1778345 : GEN z = NULL;
336 :
337 12035900 : for (j=1; j<lx; j++)
338 : {
339 10257568 : if (!F2v_coeff(y,j)) continue;
340 2840896 : if (!z) z = vecsmall_copy(gel(x,j));
341 1324544 : else F2v_add_inplace(z,gel(x,j));
342 : }
343 1778332 : if (!z) z = zero_F2v(l);
344 1778336 : return z;
345 : }
346 :
347 : GEN
348 673830 : F2m_mul(GEN x, GEN y)
349 : {
350 673830 : long i,j,l,lx=lg(x), ly=lg(y);
351 : GEN z;
352 673830 : if (ly==1) return cgetg(1,t_MAT);
353 673830 : z = cgetg(ly,t_MAT);
354 673828 : if (lx==1)
355 : {
356 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
357 0 : return z;
358 : }
359 673828 : l = coeff(x,1,1);
360 2452164 : for (j=1; j<ly; j++) gel(z,j) = F2m_F2c_mul_i(x, gel(y,j), lx, l);
361 673821 : return z;
362 : }
363 :
364 : GEN
365 0 : F2m_F2c_mul(GEN x, GEN y)
366 : {
367 0 : long l, lx = lg(x);
368 0 : if (lx==1) return cgetg(1,t_VECSMALL);
369 0 : l = coeff(x,1,1);
370 0 : return F2m_F2c_mul_i(x, y, lx, l);
371 : }
372 :
373 : static GEN
374 0 : _F2m_mul(void *data, GEN x, GEN y) { (void) data; return F2m_mul(x,y); }
375 : static GEN
376 0 : _F2m_sqr(void *data, GEN x) { (void) data; return F2m_mul(x,x); }
377 : GEN
378 0 : F2m_powu(GEN x, ulong n)
379 : {
380 0 : if (!n) return matid(lg(x)-1);
381 0 : return gen_powu(x, n,NULL, &_F2m_sqr, &_F2m_mul);
382 : }
383 :
384 : static long
385 6154252 : F2v_find_nonzero(GEN x0, GEN mask0, long m)
386 : {
387 6154252 : ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
388 6154252 : long i, l = lg(x0)-2;
389 8641041 : for (i = 0; i < l; i++)
390 : {
391 6742199 : e = *x++ & *mask++;
392 6742199 : if (e) return i*BITS_IN_LONG+vals(e)+1;
393 : }
394 1898842 : return m+1;
395 : }
396 :
397 : /* in place, destroy x */
398 : GEN
399 1134935 : F2m_ker_sp(GEN x, long deplin)
400 : {
401 : GEN y, c, d;
402 : long i, j, k, r, m, n;
403 :
404 1134935 : n = lg(x)-1;
405 1134935 : m = mael(x,1,1); r=0;
406 :
407 1134935 : d = cgetg(n+1, t_VECSMALL);
408 1134919 : c = const_F2v(m);
409 5683735 : for (k=1; k<=n; k++)
410 : {
411 4846155 : GEN xk = gel(x,k);
412 4846155 : j = F2v_find_nonzero(xk, c, m);
413 4846196 : if (j>m)
414 : {
415 1447415 : if (deplin) {
416 297368 : GEN v = zero_F2v(n);
417 900737 : for (i=1; i<k; i++)
418 603375 : if (F2v_coeff(xk, d[i])) F2v_set(v, i);
419 297362 : F2v_set(v, k); return v;
420 : }
421 1150047 : r++; d[k] = 0;
422 : }
423 : else
424 : {
425 3398781 : F2v_clear(c,j); d[k] = j;
426 3398899 : F2v_clear(xk, j);
427 47128977 : for (i=k+1; i<=n; i++)
428 : {
429 43730080 : GEN xi = gel(x,i);
430 43730080 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
431 : }
432 3398897 : F2v_set(xk, j);
433 : }
434 : }
435 837580 : if (deplin) return NULL;
436 :
437 836145 : y = zero_F2m_copy(n,r);
438 1986193 : for (j=k=1; j<=r; j++,k++)
439 : {
440 3318490 : GEN C = gel(y,j); while (d[k]) k++;
441 8936313 : for (i=1; i<k; i++)
442 7786269 : if (d[i] && F2m_coeff(x,d[i],k)) F2v_set(C,i);
443 1150044 : F2v_set(C, k);
444 : }
445 836146 : return y;
446 : }
447 :
448 : /* not memory clean */
449 : GEN
450 196362 : F2m_ker(GEN x) { return F2m_ker_sp(F2m_copy(x), 0); }
451 : GEN
452 0 : F2m_deplin(GEN x) { return F2m_ker_sp(F2m_copy(x), 1); }
453 :
454 : ulong
455 1631 : F2m_det_sp(GEN x) { return !F2m_ker_sp(x, 1); }
456 :
457 : ulong
458 0 : F2m_det(GEN x)
459 : {
460 0 : pari_sp av = avma;
461 0 : ulong d = F2m_det_sp(F2m_copy(x));
462 0 : return gc_ulong(av, d);
463 : }
464 :
465 : /* Destroy x */
466 : GEN
467 345455 : F2m_gauss_pivot(GEN x, long *rr)
468 : {
469 : GEN c, d;
470 : long i, j, k, r, m, n;
471 :
472 345455 : n = lg(x)-1; if (!n) { *rr=0; return NULL; }
473 345455 : m = mael(x,1,1); r=0;
474 :
475 345455 : d = cgetg(n+1, t_VECSMALL);
476 345454 : c = const_F2v(m);
477 1653671 : for (k=1; k<=n; k++)
478 : {
479 1308216 : GEN xk = gel(x,k);
480 1308216 : j = F2v_find_nonzero(xk, c, m);
481 1308223 : if (j>m) { r++; d[k] = 0; }
482 : else
483 : {
484 856562 : F2v_clear(c,j); d[k] = j;
485 5933862 : for (i=k+1; i<=n; i++)
486 : {
487 5077299 : GEN xi = gel(x,i);
488 5077299 : if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
489 : }
490 : }
491 : }
492 :
493 345455 : *rr = r; return gc_const((pari_sp)d, d);
494 : }
495 :
496 : long
497 63 : F2m_rank(GEN x)
498 : {
499 63 : pari_sp av = avma;
500 : long r;
501 63 : (void)F2m_gauss_pivot(F2m_copy(x),&r);
502 63 : return gc_long(av, lg(x)-1 - r);
503 : }
504 :
505 : static GEN
506 14 : F2m_inv_upper_1_ind(GEN A, long index)
507 : {
508 14 : pari_sp av = avma;
509 14 : long n = lg(A)-1, i = index, j;
510 14 : GEN u = const_vecsmall(n, 0);
511 14 : u[i] = 1;
512 21 : for (i--; i>0; i--)
513 : {
514 7 : ulong m = F2m_coeff(A,i,i+1) & uel(u,i+1); /* j = i+1 */
515 7 : for (j=i+2; j<=n; j++) m ^= F2m_coeff(A,i,j) & uel(u,j);
516 7 : u[i] = m & 1;
517 : }
518 14 : return gerepileuptoleaf(av, Flv_to_F2v(u));
519 : }
520 : static GEN
521 7 : F2m_inv_upper_1(GEN A)
522 : {
523 : long i, l;
524 7 : GEN B = cgetg_copy(A, &l);
525 21 : for (i = 1; i < l; i++) gel(B,i) = F2m_inv_upper_1_ind(A, i);
526 7 : return B;
527 : }
528 :
529 : static GEN
530 716121 : F2_get_col(GEN b, GEN d, long li, long aco)
531 : {
532 716121 : long i, l = nbits2lg(aco);
533 716122 : GEN u = cgetg(l, t_VECSMALL);
534 716134 : u[1] = aco;
535 4942006 : for (i = 1; i <= li; i++)
536 4225893 : if (d[i]) /* d[i] can still be 0 if li > aco */
537 : {
538 4156999 : if (F2v_coeff(b, i))
539 1417937 : F2v_set(u, d[i]);
540 : else
541 2739053 : F2v_clear(u, d[i]);
542 : }
543 716113 : return u;
544 : }
545 :
546 : /* destroy a, b */
547 : GEN
548 224304 : F2m_gauss_sp(GEN a, GEN b)
549 : {
550 224304 : long i, j, k, l, li, bco, aco = lg(a)-1;
551 : GEN u, d;
552 :
553 224304 : if (!aco) return cgetg(1,t_MAT);
554 224304 : li = gel(a,1)[1];
555 224304 : d = zero_Flv(li);
556 224302 : bco = lg(b)-1;
557 931761 : for (i=1; i<=aco; i++)
558 : {
559 707478 : GEN ai = vecsmall_copy(gel(a,i));
560 707478 : if (!d[i] && F2v_coeff(ai, i))
561 476799 : k = i;
562 : else
563 1055968 : for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
564 : /* found a pivot on row k */
565 707484 : if (k > li) return NULL;
566 707463 : d[k] = i;
567 :
568 : /* Clear k-th row but column-wise instead of line-wise */
569 : /* a_ij -= a_i1*a1j/a_11
570 : line-wise grouping: L_j -= a_1j/a_11*L_1
571 : column-wise: C_i -= a_i1/a_11*C_1
572 : */
573 707463 : F2v_clear(ai,k);
574 4847293 : for (l=1; l<=aco; l++)
575 : {
576 4139849 : GEN al = gel(a,l);
577 4139849 : if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
578 : }
579 4864691 : for (l=1; l<=bco; l++)
580 : {
581 4157232 : GEN bl = gel(b,l);
582 4157232 : if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
583 : }
584 : }
585 224283 : u = cgetg(bco+1,t_MAT);
586 940394 : for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
587 224275 : return u;
588 : }
589 :
590 : GEN
591 35 : F2m_gauss(GEN a, GEN b)
592 : {
593 35 : pari_sp av = avma;
594 35 : if (lg(a) == 1) return cgetg(1,t_MAT);
595 35 : return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), F2m_copy(b)));
596 : }
597 : GEN
598 14 : F2m_F2c_gauss(GEN a, GEN b)
599 : {
600 14 : pari_sp av = avma;
601 14 : GEN z = F2m_gauss(a, mkmat(b));
602 14 : if (!z) return gc_NULL(av);
603 7 : if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
604 7 : return gerepileuptoleaf(av, gel(z,1));
605 : }
606 :
607 : GEN
608 35 : F2m_inv(GEN a)
609 : {
610 35 : pari_sp av = avma;
611 35 : if (lg(a) == 1) return cgetg(1,t_MAT);
612 35 : return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), matid_F2m(gel(a,1)[1])));
613 : }
614 :
615 : GEN
616 7 : F2m_invimage_i(GEN A, GEN B)
617 : {
618 : GEN d, x, X, Y;
619 7 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
620 7 : x = F2m_ker_sp(shallowconcat(A, B), 0);
621 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
622 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
623 : * Y has at least nB columns and full rank */
624 7 : nY = lg(x)-1;
625 7 : if (nY < nB) return NULL;
626 :
627 : /* implicitly: Y = rowslice(x, nA+1, nA+nB), nB rows */
628 7 : d = cgetg(nB+1, t_VECSMALL);
629 21 : for (i = nB, j = nY; i >= 1; i--, j--)
630 : {
631 14 : for (; j>=1; j--)
632 14 : if (F2m_coeff(x,nA+i,j)) { d[i] = j; break; } /* Y[i,j] */
633 14 : if (!j) return NULL;
634 : }
635 7 : x = vecpermute(x, d);
636 :
637 7 : X = F2m_rowslice(x, 1, nA);
638 7 : Y = F2m_rowslice(x, nA+1, nA+nB);
639 7 : return F2m_mul(X, F2m_inv_upper_1(Y));
640 : }
641 : GEN
642 0 : F2m_invimage(GEN A, GEN B)
643 : {
644 0 : pari_sp av = avma;
645 0 : GEN X = F2m_invimage_i(A,B);
646 0 : if (!X) return gc_NULL(av);
647 0 : return gerepileupto(av, X);
648 : }
649 :
650 : GEN
651 192964 : F2m_F2c_invimage(GEN A, GEN y)
652 : {
653 192964 : pari_sp av = avma;
654 192964 : long i, l = lg(A);
655 : GEN M, x;
656 :
657 192964 : if (l==1) return NULL;
658 192964 : if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
659 192964 : M = cgetg(l+1,t_MAT);
660 997287 : for (i=1; i<l; i++) gel(M,i) = gel(A,i);
661 192964 : gel(M,l) = y; M = F2m_ker(M);
662 192962 : i = lg(M)-1; if (!i) return gc_NULL(av);
663 :
664 192962 : x = gel(M,i);
665 192962 : if (!F2v_coeff(x,l)) return gc_NULL(av);
666 192962 : F2v_clear(x, x[1]); x[1]--; /* remove last coord */
667 192960 : return gerepileuptoleaf(av, x);
668 : }
669 :
670 : /* Block Lanczos algorithm for kernel of sparse matrix (F2Ms)
671 : Based on lanczos.cpp by Jason Papadopoulos
672 : <https://github.com/sagemath/FlintQS/blob/master/src/lanczos.cpp>
673 : Copyright Jason Papadopoulos 2006
674 : Released under the GNU General Public License v2 or later version.
675 : */
676 :
677 : /* F2Ms are vector of vecsmall which represents nonzero entries of columns
678 : * F2w are vecsmall whoses entries are columns of a n x BIL F2m
679 : * F2wB are F2w in the special case where dim = BIL.
680 : */
681 :
682 : #define BIL BITS_IN_LONG
683 :
684 : static GEN
685 184 : F2w_transpose_F2m(GEN x)
686 : {
687 184 : long i, j, l = lg(x)-1;
688 184 : GEN z = cgetg(BIL+1, t_MAT);
689 11448 : for (j = 1; j <= BIL; j++)
690 11264 : gel(z,j) = zero_F2v(l);
691 291482 : for (i = 1; i <= l; i++)
692 : {
693 291298 : ulong xi = uel(x,i);
694 18411426 : for(j = 1; j <= BIL; j++)
695 18120128 : if (xi&(1UL<<(j-1)))
696 5326957 : F2v_set(gel(z, j), i);
697 : }
698 184 : return z;
699 : }
700 :
701 : static GEN
702 6612 : F2wB_mul(GEN a, GEN b)
703 : {
704 : long i, j;
705 6612 : GEN c = cgetg(BIL+1, t_VECSMALL);
706 407124 : for (i = 1; i <= BIL; i++)
707 : {
708 400512 : ulong s = 0, ai = a[i];
709 21101212 : for (j = 1; ai; j++, ai>>=1)
710 20700700 : if (ai & 1)
711 10421264 : s ^= b[j];
712 400512 : c[i] = s;
713 : }
714 6612 : return c;
715 : }
716 :
717 : static void
718 4408 : precompute_F2w_F2wB(GEN x, GEN c)
719 : {
720 : ulong z, xk;
721 : ulong i, j, k, index;
722 4408 : x++; c++;
723 37784 : for (j = 0; j < (BIL>>3); j++)
724 : {
725 8577632 : for (i = 0; i < 256; i++)
726 : {
727 8544256 : k = 0;
728 8544256 : index = i;
729 8544256 : z = 0;
730 68387424 : while (index) {
731 59843168 : xk = x[k];
732 59843168 : if (index & 1)
733 34177024 : z ^= xk;
734 59843168 : index >>= 1;
735 59843168 : k++;
736 : }
737 8544256 : c[i] = z;
738 : }
739 33376 : x += 8; c += 256;
740 : }
741 4408 : }
742 :
743 : static void
744 4408 : F2w_F2wB_mul_add_inplace(GEN v, GEN x, GEN y)
745 : {
746 4408 : long i, n = lg(y)-1;
747 : ulong word;
748 4408 : GEN c = cgetg(1+(BIL<<5), t_VECSMALL);
749 4408 : precompute_F2w_F2wB(x, c);
750 4408 : c++;
751 7725680 : for (i = 1; i <= n; i++)
752 : {
753 7721272 : word = v[i];
754 7721272 : y[i] ^= c[ 0*256 + ((word>> 0) & 0xff) ]
755 7721272 : ^ c[ 1*256 + ((word>> 8) & 0xff) ]
756 7721272 : ^ c[ 2*256 + ((word>>16) & 0xff) ]
757 7721272 : ^ c[ 3*256 + ((word>>24) & 0xff) ]
758 : #ifdef LONG_IS_64BIT
759 7248696 : ^ c[ 4*256 + ((word>>32) & 0xff) ]
760 7248696 : ^ c[ 5*256 + ((word>>40) & 0xff) ]
761 7248696 : ^ c[ 6*256 + ((word>>48) & 0xff) ]
762 7248696 : ^ c[ 7*256 + ((word>>56) ) ]
763 : #endif
764 : ;
765 : }
766 4408 : }
767 :
768 : /* Return x*y~, which is a F2wB */
769 : static GEN
770 3352 : F2w_transmul(GEN x, GEN y)
771 : {
772 3352 : long i, j, n = lg(x)-1;
773 3352 : GEN z = zero_zv(BIL);
774 3352 : pari_sp av = avma;
775 3352 : GEN c = zero_zv(BIL<<5) + 1;
776 3352 : GEN xy = z + 1;
777 :
778 5862977 : for (i = 1; i <= n; i++)
779 : {
780 5859625 : ulong xi = x[i];
781 5859625 : ulong yi = y[i];
782 5859625 : c[ 0*256 + ( xi & 0xff) ] ^= yi;
783 5859625 : c[ 1*256 + ((xi >> 8) & 0xff) ] ^= yi;
784 5859625 : c[ 2*256 + ((xi >> 16) & 0xff) ] ^= yi;
785 5859625 : c[ 3*256 + ((xi >> 24) & 0xff) ] ^= yi;
786 : #ifdef LONG_IS_64BIT
787 5501328 : c[ 4*256 + ((xi >> 32) & 0xff) ] ^= yi;
788 5501328 : c[ 5*256 + ((xi >> 40) & 0xff) ] ^= yi;
789 5501328 : c[ 6*256 + ((xi >> 48) & 0xff) ] ^= yi;
790 5501328 : c[ 7*256 + ((xi >> 56) ) ] ^= yi;
791 : #endif
792 : }
793 30168 : for(i = 0; i < 8; i++)
794 : {
795 26816 : ulong a0 = 0, a1 = 0, a2 = 0, a3 = 0;
796 : #ifdef LONG_IS_64BIT
797 23952 : ulong a4 = 0, a5 = 0, a6 = 0, a7 = 0;
798 : #endif
799 6891712 : for (j = 0; j < 256; j++) {
800 6864896 : if ((j >> i) & 1) {
801 3432448 : a0 ^= c[0*256 + j];
802 3432448 : a1 ^= c[1*256 + j];
803 3432448 : a2 ^= c[2*256 + j];
804 3432448 : a3 ^= c[3*256 + j];
805 : #ifdef LONG_IS_64BIT
806 3065856 : a4 ^= c[4*256 + j];
807 3065856 : a5 ^= c[5*256 + j];
808 3065856 : a6 ^= c[6*256 + j];
809 3065856 : a7 ^= c[7*256 + j];
810 : #endif
811 : }
812 : }
813 26816 : xy[ 0] = a0; xy[ 8] = a1; xy[16] = a2; xy[24] = a3;
814 : #ifdef LONG_IS_64BIT
815 23952 : xy[32] = a4; xy[40] = a5; xy[48] = a6; xy[56] = a7;
816 : #endif
817 26816 : xy++;
818 : }
819 3352 : return gc_const(av, z);
820 : }
821 :
822 : static GEN
823 1102 : identity_F2wB(void)
824 : {
825 : long i;
826 1102 : GEN M = cgetg(BIL+1, t_VECSMALL);
827 67854 : for (i = 1; i <= BIL; i++)
828 66752 : uel(M,i) = 1UL<<(i-1);
829 1102 : return M;
830 : }
831 :
832 : static GEN
833 1102 : find_nonsingular_sub(GEN t, GEN last_s, GEN *pt_s)
834 : {
835 1102 : long i, j, dim = 0;
836 : ulong mask, row_i, row_j;
837 1102 : long last_dim = lg(last_s)-1;
838 1102 : GEN s = cgetg(BIL+1, t_VECSMALL);
839 1102 : GEN M1 = identity_F2wB();
840 1102 : pari_sp av = avma;
841 1102 : GEN cols = cgetg(BIL+1, t_VECSMALL);
842 1102 : GEN M0 = zv_copy(t);
843 :
844 1102 : mask = 0;
845 66944 : for (i = 1; i <= last_dim; i++)
846 : {
847 65842 : cols[BIL + 1 - i] = last_s[i];
848 65842 : mask |= 1UL<<(last_s[i]-1);
849 : }
850 67854 : for (i = j = 1; i <= BIL; i++)
851 66752 : if (!(mask & (1UL<<(i-1))))
852 910 : cols[j++] = i;
853 :
854 : /* compute the inverse of t[][] */
855 :
856 67854 : for (i = 1; i <= BIL; i++)
857 : {
858 66752 : mask = 1UL<<(cols[i]-1);
859 66752 : row_i = cols[i];
860 130411 : for (j = i; j <= BIL; j++)
861 : {
862 128665 : row_j = cols[j];
863 128665 : if (uel(M0,row_j) & mask)
864 : {
865 65006 : swap(gel(M0, row_j), gel(M0, row_i));
866 65006 : swap(gel(M1, row_j), gel(M1, row_i));
867 65006 : break;
868 : }
869 : }
870 66752 : if (j <= BIL)
871 : {
872 4108590 : for (j = 1; j <= BIL; j++)
873 : {
874 4043584 : row_j = cols[j];
875 4043584 : if (row_i != row_j && (M0[row_j] & mask))
876 : {
877 1957597 : uel(M0,row_j) ^= uel(M0,row_i);
878 1957597 : uel(M1,row_j) ^= uel(M1,row_i);
879 : }
880 : }
881 65006 : s[++dim] = cols[i];
882 65006 : continue;
883 : }
884 1746 : for (j = i; j <= BIL; j++)
885 : {
886 1746 : row_j = cols[j];
887 1746 : if (uel(M1,row_j) & mask)
888 : {
889 1746 : swap(gel(M0, row_j), gel(M0, row_i));
890 1746 : swap(gel(M1, row_j), gel(M1, row_i));
891 1746 : break;
892 : }
893 : }
894 1746 : if (j > BIL) return NULL;
895 109458 : for (j = 1; j <= BIL; j++)
896 : {
897 107712 : row_j = cols[j];
898 107712 : if (row_i != row_j && (M1[row_j] & mask))
899 : {
900 0 : uel(M0,row_j) ^= uel(M0,row_i);
901 0 : uel(M1,row_j) ^= uel(M1,row_i);
902 : }
903 : }
904 1746 : M0[row_i] = M1[row_i] = 0;
905 : }
906 1102 : mask = 0;
907 66108 : for (i = 1; i <= dim; i++)
908 65006 : mask |= 1UL<<(s[i]-1);
909 66944 : for (i = 1; i <= last_dim; i++)
910 65842 : mask |= 1UL<<(last_s[i]-1);
911 1102 : if (mask != (ulong)(-1))
912 0 : return NULL; /* Failure */
913 1102 : setlg(s, dim+1);
914 1102 : set_avma(av);
915 1102 : *pt_s = s;
916 1102 : return M1;
917 : }
918 :
919 : /* Compute x * A~ */
920 : static GEN
921 1286 : F2w_F2Ms_transmul(GEN x, GEN A, long nbrow)
922 : {
923 1286 : long i, j, l = lg(A);
924 1286 : GEN b = zero_zv(nbrow);
925 2206288 : for (i = 1; i < l; i++)
926 : {
927 2205002 : GEN c = gel(A,i);
928 2205002 : long lc = lg(c);
929 2205002 : ulong xi = x[i];
930 41407938 : for (j = 1; j < lc; j++)
931 39202936 : b[c[j]] ^= xi;
932 : }
933 1286 : return b;
934 : }
935 :
936 : /* Compute x * A */
937 : static GEN
938 1194 : F2w_F2Ms_mul(GEN x, GEN A)
939 : {
940 1194 : long i, j, l = lg(A);
941 1194 : GEN b = cgetg(l, t_VECSMALL);
942 2068854 : for (i = 1; i < l; i++)
943 : {
944 2067660 : GEN c = gel(A,i);
945 2067660 : long lc = lg(c);
946 2067660 : ulong s = 0;
947 38865718 : for (j = 1; j < lc; j++)
948 36798058 : s ^= x[c[j]];
949 2067660 : b[i] = s;
950 : }
951 1194 : return b;
952 : }
953 :
954 : static void
955 2204 : F2wB_addid_inplace(GEN f)
956 : {
957 : long i;
958 135708 : for (i = 1; i <= BIL; i++)
959 133504 : uel(f,i) ^= 1UL<<(i-1);
960 2204 : }
961 :
962 : static void
963 2204 : F2w_mask_inplace(GEN f, ulong m)
964 : {
965 2204 : long i, l = lg(f);
966 1999274 : for (i = 1; i < l; i++)
967 1997070 : uel(f,i) &= m;
968 2204 : }
969 :
970 : static GEN
971 46 : block_lanczos(GEN B, ulong nbrow)
972 : {
973 46 : pari_sp av = avma, av2;
974 : GEN v0, v1, v2, vnext, x, w;
975 : GEN winv0, winv1, winv2;
976 : GEN vt_a_v0, vt_a_v1, vt_a2_v0, vt_a2_v1;
977 : GEN d, e, f, f2, s0;
978 : long i, iter;
979 46 : long n = lg(B)-1;
980 : long dim0;
981 : ulong mask0, mask1;
982 46 : v1 = zero_zv(n);
983 46 : v2 = zero_zv(n);
984 46 : vt_a_v1 = zero_zv(BIL);
985 46 : vt_a2_v1 = zero_zv(BIL);
986 46 : winv1 = zero_zv(BIL);
987 46 : winv2 = zero_zv(BIL);
988 46 : s0 = identity_zv(BIL);
989 46 : mask1 = (ulong)(-1);
990 :
991 46 : x = random_zv(n);
992 46 : w = F2w_F2Ms_mul(F2w_F2Ms_transmul(x, B, nbrow), B);
993 46 : v0 = w;
994 46 : av2 = avma;
995 46 : for (iter=1;;iter++)
996 : {
997 1148 : vnext = F2w_F2Ms_mul(F2w_F2Ms_transmul(v0, B, nbrow), B);
998 1148 : vt_a_v0 = F2w_transmul(v0, vnext);
999 1148 : if (zv_equal0(vt_a_v0)) break; /* success */
1000 1102 : vt_a2_v0 = F2w_transmul(vnext, vnext);
1001 1102 : winv0 = find_nonsingular_sub(vt_a_v0, s0, &s0);
1002 1102 : if (!winv0) return gc_NULL(av); /* failure */
1003 1102 : dim0 = lg(s0)-1;
1004 1102 : mask0 = 0;
1005 66108 : for (i = 1; i <= dim0; i++)
1006 65006 : mask0 |= 1UL<<(s0[i]-1);
1007 1102 : d = cgetg(BIL+1, t_VECSMALL);
1008 67854 : for (i = 1; i <= BIL; i++)
1009 66752 : d[i] = (vt_a2_v0[i] & mask0) ^ vt_a_v0[i];
1010 :
1011 1102 : d = F2wB_mul(winv0, d);
1012 1102 : F2wB_addid_inplace(d);
1013 1102 : e = F2wB_mul(winv1, vt_a_v0);
1014 1102 : F2w_mask_inplace(e, mask0);
1015 1102 : f = F2wB_mul(vt_a_v1, winv1);
1016 1102 : F2wB_addid_inplace(f);
1017 1102 : f = F2wB_mul(winv2, f);
1018 1102 : f2 = cgetg(BIL+1, t_VECSMALL);
1019 67854 : for (i = 1; i <= BIL; i++)
1020 66752 : f2[i] = ((vt_a2_v1[i] & mask1) ^ vt_a_v1[i]) & mask0;
1021 :
1022 1102 : f = F2wB_mul(f, f2);
1023 1102 : F2w_mask_inplace(vnext, mask0);
1024 1102 : F2w_F2wB_mul_add_inplace(v0, d, vnext);
1025 1102 : F2w_F2wB_mul_add_inplace(v1, e, vnext);
1026 1102 : F2w_F2wB_mul_add_inplace(v2, f, vnext);
1027 1102 : d = F2wB_mul(winv0, F2w_transmul(v0, w));
1028 1102 : F2w_F2wB_mul_add_inplace(v0, d, x);
1029 1102 : v2 = v1; v1 = v0; v0 = vnext;
1030 1102 : winv2 = winv1; winv1 = winv0;
1031 1102 : vt_a_v1 = vt_a_v0;
1032 1102 : vt_a2_v1 = vt_a2_v0;
1033 1102 : mask1 = mask0;
1034 1102 : gerepileall(av2, 9, &x, &s0, &v0, &v1, &v2,
1035 : &winv1, &winv2, &vt_a_v1, &vt_a2_v1);
1036 : }
1037 46 : if (DEBUGLEVEL >= 5)
1038 0 : err_printf("Lanczos halted after %ld iterations\n", iter);
1039 46 : v1 = F2w_F2Ms_transmul(x, B, nbrow);
1040 46 : v2 = F2w_F2Ms_transmul(v0, B, nbrow);
1041 46 : x = shallowconcat(F2w_transpose_F2m(x), F2w_transpose_F2m(v0));
1042 46 : v1 = shallowconcat(F2w_transpose_F2m(v1), F2w_transpose_F2m(v2));
1043 46 : s0 = gel(F2m_indexrank(x), 2);
1044 46 : x = shallowextract(x, s0);
1045 46 : v1 = shallowextract(v1, s0);
1046 46 : return F2m_mul(x, F2m_ker(v1));
1047 : }
1048 :
1049 : static GEN
1050 2512 : F2v_inflate(GEN v, GEN p, long n)
1051 : {
1052 2512 : long i, l = lg(p)-1;
1053 2512 : GEN w = zero_F2v(n);
1054 3944198 : for (i=1; i<=l; i++)
1055 3941686 : if (F2v_coeff(v,i))
1056 1968382 : F2v_set(w, p[i]);
1057 2512 : return w;
1058 : }
1059 :
1060 : static GEN
1061 46 : F2m_inflate(GEN x, GEN p, long n)
1062 2558 : { pari_APPLY_same(F2v_inflate(gel(x,i), p, n)) }
1063 :
1064 : GEN
1065 2495 : F2Ms_ker(GEN M, long nbrow)
1066 : {
1067 2495 : pari_sp av = avma;
1068 2495 : long nbcol = lg(M)-1;
1069 : GEN Mp, R, Rp, p;
1070 2495 : if (nbrow <= 640)
1071 2449 : return gerepileupto(av, F2m_ker(F2Ms_to_F2m(M, nbrow)));
1072 46 : p = F2Ms_colelim(M, nbrow);
1073 46 : Mp = vecpermute(M, p);
1074 : do
1075 : {
1076 46 : R = block_lanczos(Mp, nbrow);
1077 46 : } while(!R);
1078 46 : Rp = F2m_inflate(R, p, nbcol);
1079 46 : return gerepilecopy(av, Rp);
1080 : }
1081 :
1082 : GEN
1083 0 : F2m_to_F2Ms(GEN M)
1084 : {
1085 0 : long ncol = lg(M)-1;
1086 0 : GEN B = cgetg(ncol+1, t_MAT);
1087 : long i, j, k;
1088 0 : for(i = 1; i <= ncol; i++)
1089 : {
1090 0 : GEN D, V = gel(M,i);
1091 0 : long n = F2v_hamming(V), l = V[1];
1092 0 : D = cgetg(n+1, t_VECSMALL);
1093 0 : for (j=1, k=1; j<=l; j++)
1094 0 : if( F2v_coeff(V,j))
1095 0 : D[k++] = j;
1096 0 : gel(B, i) = D;
1097 : }
1098 0 : return B;
1099 : }
1100 :
1101 : GEN
1102 2449 : F2Ms_to_F2m(GEN M, long nrow)
1103 : {
1104 2449 : long i, j, l = lg(M);
1105 2449 : GEN B = cgetg(l, t_MAT);
1106 257475 : for(i = 1; i < l; i++)
1107 : {
1108 255026 : GEN Bi = zero_F2v(nrow), Mi = gel(M,i);
1109 255026 : long l = lg(Mi);
1110 3122670 : for (j = 1; j < l; j++)
1111 2867644 : F2v_set(Bi, Mi[j]);
1112 255026 : gel(B, i) = Bi;
1113 : }
1114 2449 : return B;
1115 : }
|
