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