Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - F2v.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.0 lcov report (development 29789-eb252ba62d) Lines: 554 634 87.4 %
Date: 2024-12-03 09:08:48 Functions: 55 68 80.9 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.16