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.16.2 lcov report (development 29115-f22e516b23) Lines: 555 634 87.5 %
Date: 2024-03-29 08:06:26 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             : 
      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     3485604 : F2c_to_ZC(GEN x)
      46             : {
      47     3485604 :   long l = x[1]+1, lx = lg(x);
      48     3485604 :   GEN  z = cgetg(l, t_COL);
      49             :   long i, j, k;
      50     6971199 :   for (i=2, k=1; i<lx; i++)
      51    20068205 :     for (j=0; j<BITS_IN_LONG && k<l; j++,k++)
      52    16582567 :       gel(z,k) = (x[i]&(1UL<<j))? gen_1: gen_0;
      53     3485561 :   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     4298534 : 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     9616698 : ZV_to_F2v(GEN x)
     123             : {
     124     9616698 :   long i, j, k, l = lg(x)-1;
     125     9616698 :   GEN z = cgetg(nbits2lg(l), t_VECSMALL);
     126     9616551 :   z[1] = l;
     127    66432118 :   for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j++)
     128             :   {
     129    56815359 :     if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
     130    56815359 :     if (mpodd(gel(x,i))) z[k] |= 1UL<<j;
     131             :   }
     132     9616759 :   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    12238057 : 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     1633742 : const_F2v(long m)
     174             : {
     175     1633742 :   long i, l = nbits2lg(m);
     176     1633729 :   GEN c = cgetg(l, t_VECSMALL);
     177     1633688 :   c[1] = m;
     178     3276345 :   for (i = 2; i < l; i++) uel(c,i) = -1UL;
     179     1633688 :   if (remsBIL(m)) uel(c,l-1) = (1UL<<remsBIL(m))-1UL;
     180     1633682 :   return c;
     181             : }
     182             : 
     183             : /* Allow lg(y)<lg(x) */
     184             : void
     185    22063614 : F2v_add_inplace(GEN x, GEN y)
     186             : {
     187    22063614 :   long n = lg(y);
     188    22063614 :   long r = (n-2)&7L, q = n-r, i;
     189    31534791 :   for (i = 2; i < q; i += 8)
     190             :   {
     191     9471177 :     x[  i] ^= y[  i]; x[1+i] ^= y[1+i]; x[2+i] ^= y[2+i]; x[3+i] ^= y[3+i];
     192     9471177 :     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    22063614 :   switch (r)
     195             :   {
     196     1768843 :     case 7: x[i] ^= y[i]; i++; case 6: x[i] ^= y[i]; i++;
     197     3455634 :     case 5: x[i] ^= y[i]; i++; case 4: x[i] ^= y[i]; i++;
     198     7969843 :     case 3: x[i] ^= y[i]; i++; case 2: x[i] ^= y[i]; i++;
     199    19285288 :     case 1: x[i] ^= y[i]; i++;
     200             :   }
     201    22063614 : }
     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         602 : F2v_negimply_inplace(GEN x, GEN y)
     246             : {
     247         602 :   long n = lg(y);
     248         602 :   long r = (n-2)&7L, q = n-r, i;
     249       91362 :   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         602 :   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         134 :     case 3: x[i] &= ~y[i]; i++; case 2: x[i] &= ~y[i]; i++;
     259         602 :     case 1: x[i] &= ~y[i]; i++;
     260             :   }
     261         602 : }
     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        1708 : F2v_hamming(GEN H)
     284             : {
     285        1708 :   long i, n=0, l=lg(H);
     286     1093086 :   for (i=2; i<l; i++) n += hammingl(uel(H,i));
     287        1708 :   return n;
     288             : }
     289             : 
     290             : int
     291        3171 : F2v_subset(GEN x, GEN y)
     292             : {
     293        3171 :   long i, n = lg(y);
     294        3738 :   for (i = 2; i < n; i ++)
     295        3374 :     if ((x[i] & y[i]) != x[i]) return 0;
     296         364 :   return 1;
     297             : }
     298             : 
     299             : GEN
     300      224763 : matid_F2m(long n)
     301             : {
     302      224763 :   GEN y = cgetg(n+1,t_MAT);
     303             :   long i;
     304      224763 :   if (n < 0) pari_err_DOMAIN("matid_F2m", "dimension","<",gen_0,stoi(n));
     305      943081 :   for (i=1; i<=n; i++) { gel(y,i) = zero_F2v(n); F2v_set(gel(y,i),i); }
     306      224769 :   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     1782903 : F2m_F2c_mul_i(GEN x, GEN y, long lx, long l)
     333             : {
     334             :   long j;
     335     1782903 :   GEN z = NULL;
     336             : 
     337    12105429 :   for (j=1; j<lx; j++)
     338             :   {
     339    10322528 :     if (!F2v_coeff(y,j)) continue;
     340     2854116 :     if (!z) z = vecsmall_copy(gel(x,j));
     341     1333531 :     else F2v_add_inplace(z,gel(x,j));
     342             :   }
     343     1782901 :   if (!z) z = zero_F2v(l);
     344     1782903 :   return z;
     345             : }
     346             : 
     347             : GEN
     348      674828 : F2m_mul(GEN x, GEN y)
     349             : {
     350      674828 :   long i,j,l,lx=lg(x), ly=lg(y);
     351             :   GEN z;
     352      674828 :   if (ly==1) return cgetg(1,t_MAT);
     353      674828 :   z = cgetg(ly,t_MAT);
     354      674820 :   if (lx==1)
     355             :   {
     356           0 :     for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
     357           0 :     return z;
     358             :   }
     359      674820 :   l = coeff(x,1,1);
     360     2457721 :   for (j=1; j<ly; j++) gel(z,j) = F2m_F2c_mul_i(x, gel(y,j), lx, l);
     361      674821 :   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     6059690 : F2v_find_nonzero(GEN x0, GEN mask0, long m)
     386             : {
     387     6059690 :   ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
     388     6059690 :   long i, l = lg(x0)-2;
     389     8523191 :   for (i = 0; i < l; i++)
     390             :   {
     391     6634832 :     e = *x++ & *mask++;
     392     6634832 :     if (e) return i*BITS_IN_LONG+vals(e)+1;
     393             :   }
     394     1888359 :   return m+1;
     395             : }
     396             : 
     397             : /* in place, destroy x */
     398             : GEN
     399     1133998 : F2m_ker_sp(GEN x, long deplin)
     400             : {
     401             :   GEN y, c, d;
     402             :   long i, j, k, r, m, n;
     403             : 
     404     1133998 :   n = lg(x)-1;
     405     1133998 :   m = mael(x,1,1); r=0;
     406             : 
     407     1133998 :   d = cgetg(n+1, t_VECSMALL);
     408     1133987 :   c = const_F2v(m);
     409     5582238 :   for (k=1; k<=n; k++)
     410             :   {
     411     4745959 :     GEN xk = gel(x,k);
     412     4745959 :     j = F2v_find_nonzero(xk, c, m);
     413     4746005 :     if (j>m)
     414             :     {
     415     1434556 :       if (deplin) {
     416      297730 :         GEN v = zero_F2v(n);
     417      902740 :         for (i=1; i<k; i++)
     418      605013 :           if (F2v_coeff(xk, d[i])) F2v_set(v, i);
     419      297727 :         F2v_set(v, k); return v;
     420             :       }
     421     1136826 :       r++; d[k] = 0;
     422             :     }
     423             :     else
     424             :     {
     425     3311449 :       F2v_clear(c,j); d[k] = j;
     426     3311543 :       F2v_clear(xk, j);
     427    43097703 :       for (i=k+1; i<=n; i++)
     428             :       {
     429    39786170 :         GEN xi = gel(x,i);
     430    39786170 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     431             :       }
     432     3311533 :       F2v_set(xk, j);
     433             :     }
     434             :   }
     435      836279 :   if (deplin) return NULL;
     436             : 
     437      834844 :   y = zero_F2m_copy(n,r);
     438     1971675 :   for (j=k=1; j<=r; j++,k++)
     439             :   {
     440     3216598 :     GEN C = gel(y,j); while (d[k]) k++;
     441     7991415 :     for (i=1; i<k; i++)
     442     6854593 :       if (d[i] && F2m_coeff(x,d[i],k)) F2v_set(C,i);
     443     1136822 :     F2v_set(C, k);
     444             :   }
     445      834846 :   return y;
     446             : }
     447             : 
     448             : /* not memory clean */
     449             : GEN
     450      194948 : 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      345450 : F2m_gauss_pivot(GEN x, long *rr)
     468             : {
     469             :   GEN c, d;
     470             :   long i, j, k, r, m, n;
     471             : 
     472      345450 :   n = lg(x)-1; if (!n) { *rr=0; return NULL; }
     473      345450 :   m = mael(x,1,1); r=0;
     474             : 
     475      345450 :   d = cgetg(n+1, t_VECSMALL);
     476      345451 :   c = const_F2v(m);
     477     1659256 :   for (k=1; k<=n; k++)
     478             :   {
     479     1313809 :     GEN xk = gel(x,k);
     480     1313809 :     j = F2v_find_nonzero(xk, c, m);
     481     1313818 :     if (j>m) { r++; d[k] = 0; }
     482             :     else
     483             :     {
     484      859802 :       F2v_clear(c,j); d[k] = j;
     485     6003127 :       for (i=k+1; i<=n; i++)
     486             :       {
     487     5143326 :         GEN xi = gel(x,i);
     488     5143326 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     489             :       }
     490             :     }
     491             :   }
     492             : 
     493      345447 :   *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      718289 : F2_get_col(GEN b, GEN d, long li, long aco)
     531             : {
     532      718289 :   long i, l = nbits2lg(aco);
     533      718293 :   GEN u = cgetg(l, t_VECSMALL);
     534      718297 :   u[1] = aco;
     535     4961528 :   for (i = 1; i <= li; i++)
     536     4243252 :     if (d[i]) /* d[i] can still be 0 if li > aco */
     537             :     {
     538     4174358 :       if (F2v_coeff(b, i))
     539     1422830 :         F2v_set(u, d[i]);
     540             :       else
     541     2751527 :         F2v_clear(u, d[i]);
     542             :     }
     543      718276 :   return u;
     544             : }
     545             : 
     546             : /* destroy a, b */
     547             : GEN
     548      224804 : F2m_gauss_sp(GEN a, GEN b)
     549             : {
     550      224804 :   long i, j, k, l, li, bco, aco = lg(a)-1;
     551             :   GEN u, d;
     552             : 
     553      224804 :   if (!aco) return cgetg(1,t_MAT);
     554      224804 :   li = gel(a,1)[1];
     555      224804 :   d = zero_Flv(li);
     556      224803 :   bco = lg(b)-1;
     557      934415 :   for (i=1; i<=aco; i++)
     558             :   {
     559      709658 :     GEN ai = vecsmall_copy(gel(a,i));
     560      709658 :     if (!d[i] && F2v_coeff(ai, i))
     561      478161 :       k = i;
     562             :     else
     563     1059924 :       for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
     564             :     /* found a pivot on row k */
     565      709660 :     if (k > li) return NULL;
     566      709639 :     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      709639 :     F2v_clear(ai,k);
     574     4866851 :     for (l=1; l<=aco; l++)
     575             :     {
     576     4157233 :       GEN al = gel(a,l);
     577     4157233 :       if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
     578             :     }
     579     4884238 :     for (l=1; l<=bco; l++)
     580             :     {
     581     4174626 :       GEN bl = gel(b,l);
     582     4174626 :       if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
     583             :     }
     584             :   }
     585      224757 :   u = cgetg(bco+1,t_MAT);
     586      943054 :   for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
     587      224768 :   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      193096 : F2m_F2c_invimage(GEN A, GEN y)
     652             : {
     653      193096 :   pari_sp av = avma;
     654      193096 :   long i, l = lg(A);
     655             :   GEN M, x;
     656             : 
     657      193096 :   if (l==1) return NULL;
     658      193096 :   if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
     659      193096 :   M = cgetg(l+1,t_MAT);
     660      999448 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
     661      193096 :   gel(M,l) = y; M = F2m_ker(M);
     662      193099 :   i = lg(M)-1; if (!i) return gc_NULL(av);
     663             : 
     664      193099 :   x = gel(M,i);
     665      193099 :   if (!F2v_coeff(x,l)) return gc_NULL(av);
     666      193100 :   F2v_clear(x, x[1]); x[1]--; /* remove last coord */
     667      193100 :   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         212 : F2w_transpose_F2m(GEN x)
     686             : {
     687         212 :   long i, j, l = lg(x)-1;
     688         212 :   GEN z = cgetg(BIL+1, t_MAT);
     689       13140 :   for (j = 1; j <= BIL; j++)
     690       12928 :     gel(z,j) = zero_F2v(l);
     691      322968 :   for (i = 1; i <= l; i++)
     692             :   {
     693      322756 :     ulong xi = uel(x,i);
     694    20312388 :     for(j = 1; j <= BIL; j++)
     695    19989632 :       if (xi&(1UL<<(j-1)))
     696     5974601 :         F2v_set(gel(z, j), i);
     697             :   }
     698         212 :   return z;
     699             : }
     700             : 
     701             : static GEN
     702        7548 : F2wB_mul(GEN a, GEN b)
     703             : {
     704             :   long i, j;
     705        7548 :   GEN c = cgetg(BIL+1, t_VECSMALL);
     706      456444 :   for (i = 1; i <= BIL; i++)
     707             :   {
     708      448896 :     ulong s = 0, ai = a[i];
     709    23322104 :     for (j = 1; ai; j++, ai>>=1)
     710    22873208 :       if (ai & 1)
     711    11518466 :         s ^= b[j];
     712      448896 :     c[i] = s;
     713             :   }
     714        7548 :   return c;
     715             : }
     716             : 
     717             : static void
     718        5032 : precompute_F2w_F2wB(GEN x, GEN c)
     719             : {
     720             :   ulong z, xk;
     721             :   ulong i, j, k, index;
     722        5032 :   x++; c++;
     723       42440 :   for (j = 0; j < (BIL>>3); j++)
     724             :   {
     725     9613856 :     for (i = 0; i < 256; i++)
     726             :     {
     727     9576448 :       k = 0;
     728     9576448 :       index = i;
     729     9576448 :       z = 0;
     730    76648992 :       while (index) {
     731    67072544 :         xk = x[k];
     732    67072544 :         if (index & 1)
     733    38305792 :           z ^= xk;
     734    67072544 :         index >>= 1;
     735    67072544 :         k++;
     736             :       }
     737     9576448 :       c[i] = z;
     738             :     }
     739       37408 :     x += 8; c += 256;
     740             :   }
     741        5032 : }
     742             : 
     743             : static void
     744        5032 : F2w_F2wB_mul_add_inplace(GEN v, GEN x, GEN y)
     745             : {
     746        5032 :   long i, n = lg(y)-1;
     747             :   ulong word;
     748        5032 :   GEN c = cgetg(1+(BIL<<5), t_VECSMALL);
     749        5032 :   precompute_F2w_F2wB(x, c);
     750        5032 :   c++;
     751     8364096 :   for (i = 1; i <= n; i++)
     752             :   {
     753     8359064 :     word = v[i];
     754     8359064 :     y[i] ^=  c[ 0*256 + ((word>> 0) & 0xff) ]
     755     8359064 :            ^ c[ 1*256 + ((word>> 8) & 0xff) ]
     756     8359064 :            ^ c[ 2*256 + ((word>>16) & 0xff) ]
     757     8359064 :            ^ c[ 3*256 + ((word>>24) & 0xff) ]
     758             : #ifdef LONG_IS_64BIT
     759     7649592 :            ^ c[ 4*256 + ((word>>32) & 0xff) ]
     760     7649592 :            ^ c[ 5*256 + ((word>>40) & 0xff) ]
     761     7649592 :            ^ c[ 6*256 + ((word>>48) & 0xff) ]
     762     7649592 :            ^ c[ 7*256 + ((word>>56)       ) ]
     763             : #endif
     764             :            ;
     765             :   }
     766        5032 : }
     767             : 
     768             : /* Return x*y~, which is a F2wB */
     769             : static GEN
     770        3829 : F2w_transmul(GEN x, GEN y)
     771             : {
     772        3829 :   long i, j, n = lg(x)-1;
     773        3829 :   GEN z = zero_zv(BIL);
     774        3829 :   pari_sp av = avma;
     775        3829 :   GEN c = zero_zv(BIL<<5) + 1;
     776        3829 :   GEN xy = z + 1;
     777             : 
     778     6350950 :   for (i = 1; i <= n; i++)
     779             :   {
     780     6347121 :     ulong xi = x[i];
     781     6347121 :     ulong yi = y[i];
     782     6347121 :     c[ 0*256 + ( xi        & 0xff) ] ^= yi;
     783     6347121 :     c[ 1*256 + ((xi >>  8) & 0xff) ] ^= yi;
     784     6347121 :     c[ 2*256 + ((xi >> 16) & 0xff) ] ^= yi;
     785     6347121 :     c[ 3*256 + ((xi >> 24) & 0xff) ] ^= yi;
     786             : #ifdef LONG_IS_64BIT
     787     5808264 :     c[ 4*256 + ((xi >> 32) & 0xff) ] ^= yi;
     788     5808264 :     c[ 5*256 + ((xi >> 40) & 0xff) ] ^= yi;
     789     5808264 :     c[ 6*256 + ((xi >> 48) & 0xff) ] ^= yi;
     790     5808264 :     c[ 7*256 + ((xi >> 56)       ) ] ^= yi;
     791             : #endif
     792             :   }
     793       34461 :   for(i = 0; i < 8; i++)
     794             :   {
     795       30632 :     ulong a0 = 0, a1 = 0, a2 = 0, a3 = 0;
     796             : #ifdef LONG_IS_64BIT
     797       26304 :     ulong a4 = 0, a5 = 0, a6 = 0, a7 = 0;
     798             : #endif
     799     7872424 :     for (j = 0; j < 256; j++) {
     800     7841792 :       if ((j >> i) & 1) {
     801     3920896 :         a0 ^= c[0*256 + j];
     802     3920896 :         a1 ^= c[1*256 + j];
     803     3920896 :         a2 ^= c[2*256 + j];
     804     3920896 :         a3 ^= c[3*256 + j];
     805             : #ifdef LONG_IS_64BIT
     806     3366912 :         a4 ^= c[4*256 + j];
     807     3366912 :         a5 ^= c[5*256 + j];
     808     3366912 :         a6 ^= c[6*256 + j];
     809     3366912 :         a7 ^= c[7*256 + j];
     810             : #endif
     811             :       }
     812             :     }
     813       30632 :     xy[ 0] = a0; xy[ 8] = a1; xy[16] = a2; xy[24] = a3;
     814             : #ifdef LONG_IS_64BIT
     815       26304 :     xy[32] = a4; xy[40] = a5; xy[48] = a6; xy[56] = a7;
     816             : #endif
     817       30632 :     xy++;
     818             :   }
     819        3829 :   return gc_const(av, z);
     820             : }
     821             : 
     822             : static GEN
     823        1259 : identity_F2wB(void)
     824             : {
     825             :   long i;
     826        1259 :   GEN M = cgetg(BIL+1, t_VECSMALL);
     827       76107 :   for (i = 1; i <= BIL; i++)
     828       74848 :     uel(M,i) = 1UL<<(i-1);
     829        1259 :   return M;
     830             : }
     831             : 
     832             : static GEN
     833        1259 : find_nonsingular_sub(GEN t, GEN last_s, GEN *pt_s)
     834             : {
     835        1259 :   long i, j, dim = 0;
     836             :   ulong mask, row_i, row_j;
     837        1259 :   long last_dim = lg(last_s)-1;
     838        1259 :   GEN s = cgetg(BIL+1, t_VECSMALL);
     839        1259 :   GEN M1 = identity_F2wB();
     840        1259 :   pari_sp av = avma;
     841        1259 :   GEN cols = cgetg(BIL+1, t_VECSMALL);
     842        1259 :   GEN M0 = zv_copy(t);
     843             : 
     844        1259 :   mask = 0;
     845       75102 :   for (i = 1; i <= last_dim; i++)
     846             :   {
     847       73843 :     cols[BIL + 1 - i] = last_s[i];
     848       73843 :     mask |= 1UL<<(last_s[i]-1);
     849             :   }
     850       76107 :   for (i = j = 1; i <= BIL; i++)
     851       74848 :     if (!(mask & (1UL<<(i-1))))
     852        1005 :       cols[j++] = i;
     853             : 
     854             :   /* compute the inverse of t[][] */
     855             : 
     856       76107 :   for (i = 1; i <= BIL; i++)
     857             :   {
     858       74848 :     mask = 1UL<<(cols[i]-1);
     859       74848 :     row_i = cols[i];
     860      147699 :     for (j = i; j <= BIL; j++)
     861             :     {
     862      145617 :       row_j = cols[j];
     863      145617 :       if (uel(M0,row_j) & mask)
     864             :       {
     865       72766 :         swap(gel(M0, row_j), gel(M0, row_i));
     866       72766 :         swap(gel(M1, row_j), gel(M1, row_i));
     867       72766 :         break;
     868             :       }
     869             :     }
     870       74848 :     if (j <= BIL)
     871             :     {
     872     4553598 :       for (j = 1; j <= BIL; j++)
     873             :       {
     874     4480832 :         row_j = cols[j];
     875     4480832 :         if (row_i != row_j && (M0[row_j] & mask))
     876             :         {
     877     2168618 :           uel(M0,row_j) ^= uel(M0,row_i);
     878     2168618 :           uel(M1,row_j) ^= uel(M1,row_i);
     879             :         }
     880             :       }
     881       72766 :       s[++dim] = cols[i];
     882       72766 :       continue;
     883             :     }
     884        2082 :     for (j = i; j <= BIL; j++)
     885             :     {
     886        2082 :       row_j = cols[j];
     887        2082 :       if (uel(M1,row_j) & mask)
     888             :       {
     889        2082 :         swap(gel(M0, row_j), gel(M0, row_i));
     890        2082 :         swap(gel(M1, row_j), gel(M1, row_i));
     891        2082 :         break;
     892             :       }
     893             :     }
     894        2082 :     if (j > BIL) return NULL;
     895      128226 :     for (j = 1; j <= BIL; j++)
     896             :     {
     897      126144 :       row_j = cols[j];
     898      126144 :       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        2082 :     M0[row_i] = M1[row_i] = 0;
     905             :   }
     906        1259 :   mask = 0;
     907       74025 :   for (i = 1; i <= dim; i++)
     908       72766 :     mask |= 1UL<<(s[i]-1);
     909       75102 :   for (i = 1; i <= last_dim; i++)
     910       73843 :     mask |= 1UL<<(last_s[i]-1);
     911        1259 :   if (mask != (ulong)(-1))
     912           1 :     return NULL; /* Failure */
     913        1258 :   setlg(s, dim+1);
     914        1258 :   set_avma(av);
     915        1258 :   *pt_s = s;
     916        1258 :   return M1;
     917             : }
     918             : 
     919             : /* Compute x * A~ */
     920             : static GEN
     921        1472 : F2w_F2Ms_transmul(GEN x, GEN A, long nbrow)
     922             : {
     923        1472 :   long i, j, l = lg(A);
     924        1472 :   GEN b = zero_zv(nbrow);
     925     2396998 :   for (i = 1; i < l; i++)
     926             :   {
     927     2395526 :     GEN c = gel(A,i);
     928     2395526 :     long lc = lg(c);
     929     2395526 :     ulong xi = x[i];
     930    44876922 :     for (j = 1; j < lc; j++)
     931    42481396 :       b[c[j]] ^= xi;
     932             :   }
     933        1472 :   return b;
     934             : }
     935             : 
     936             : /* Compute x * A */
     937             : static GEN
     938        1366 : F2w_F2Ms_mul(GEN x, GEN A)
     939             : {
     940        1366 :   long i, j, l = lg(A);
     941        1366 :   GEN b = cgetg(l, t_VECSMALL);
     942     2244934 :   for (i = 1; i < l; i++)
     943             :   {
     944     2243568 :     GEN c = gel(A,i);
     945     2243568 :     long lc = lg(c);
     946     2243568 :     ulong s = 0;
     947    42065986 :     for (j = 1; j < lc; j++)
     948    39822418 :       s ^= x[c[j]];
     949     2243568 :     b[i] = s;
     950             :   }
     951        1366 :   return b;
     952             : }
     953             : 
     954             : static void
     955        2516 : F2wB_addid_inplace(GEN f)
     956             : {
     957             :   long i;
     958      152148 :   for (i = 1; i <= BIL; i++)
     959      149632 :     uel(f,i) ^= 1UL<<(i-1);
     960        2516 : }
     961             : 
     962             : static void
     963        2516 : F2w_mask_inplace(GEN f, ulong m)
     964             : {
     965        2516 :   long i, l = lg(f);
     966     2167098 :   for (i = 1; i < l; i++)
     967     2164582 :     uel(f,i) &= m;
     968        2516 : }
     969             : 
     970             : static GEN
     971          54 : block_lanczos(GEN B, ulong nbrow)
     972             : {
     973          54 :   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          54 :   long n = lg(B)-1;
     980             :   long dim0;
     981             :   ulong mask0, mask1;
     982          54 :   v1 = zero_zv(n);
     983          54 :   v2 = zero_zv(n);
     984          54 :   vt_a_v1 = zero_zv(BIL);
     985          54 :   vt_a2_v1 = zero_zv(BIL);
     986          54 :   winv1 = zero_zv(BIL);
     987          54 :   winv2 = zero_zv(BIL);
     988          54 :   s0 = identity_zv(BIL);
     989          54 :   mask1 = (ulong)(-1);
     990             : 
     991          54 :   x = random_zv(n);
     992          54 :   w = F2w_F2Ms_mul(F2w_F2Ms_transmul(x, B, nbrow), B);
     993          54 :   v0 = w;
     994          54 :   av2 = avma;
     995          54 :   for (iter=1;;iter++)
     996             :   {
     997        1312 :     vnext = F2w_F2Ms_mul(F2w_F2Ms_transmul(v0, B, nbrow), B);
     998        1312 :     vt_a_v0  = F2w_transmul(v0, vnext);
     999        1312 :     if (zv_equal0(vt_a_v0)) break; /* success */
    1000        1259 :     vt_a2_v0 = F2w_transmul(vnext, vnext);
    1001        1259 :     winv0 = find_nonsingular_sub(vt_a_v0, s0, &s0);
    1002        1259 :     if (!winv0) return gc_NULL(av); /* failure */
    1003        1258 :     dim0 = lg(s0)-1;
    1004        1258 :     mask0 = 0;
    1005       74020 :     for (i = 1; i <= dim0; i++)
    1006       72762 :       mask0 |= 1UL<<(s0[i]-1);
    1007        1258 :     d = cgetg(BIL+1, t_VECSMALL);
    1008       76074 :     for (i = 1; i <= BIL; i++)
    1009       74816 :       d[i] = (vt_a2_v0[i] & mask0) ^ vt_a_v0[i];
    1010             : 
    1011        1258 :     d = F2wB_mul(winv0, d);
    1012        1258 :     F2wB_addid_inplace(d);
    1013        1258 :     e = F2wB_mul(winv1, vt_a_v0);
    1014        1258 :     F2w_mask_inplace(e, mask0);
    1015        1258 :     f = F2wB_mul(vt_a_v1, winv1);
    1016        1258 :     F2wB_addid_inplace(f);
    1017        1258 :     f = F2wB_mul(winv2, f);
    1018        1258 :     f2 = cgetg(BIL+1, t_VECSMALL);
    1019       76074 :     for (i = 1; i <= BIL; i++)
    1020       74816 :       f2[i] = ((vt_a2_v1[i] & mask1) ^ vt_a_v1[i]) & mask0;
    1021             : 
    1022        1258 :     f = F2wB_mul(f, f2);
    1023        1258 :     F2w_mask_inplace(vnext, mask0);
    1024        1258 :     F2w_F2wB_mul_add_inplace(v0, d, vnext);
    1025        1258 :     F2w_F2wB_mul_add_inplace(v1, e, vnext);
    1026        1258 :     F2w_F2wB_mul_add_inplace(v2, f, vnext);
    1027        1258 :     d = F2wB_mul(winv0, F2w_transmul(v0, w));
    1028        1258 :     F2w_F2wB_mul_add_inplace(v0, d, x);
    1029        1258 :     v2 = v1; v1 = v0; v0 = vnext;
    1030        1258 :     winv2 = winv1; winv1 = winv0;
    1031        1258 :     vt_a_v1 = vt_a_v0;
    1032        1258 :     vt_a2_v1 = vt_a2_v0;
    1033        1258 :     mask1 = mask0;
    1034        1258 :     gerepileall(av2, 9, &x, &s0, &v0, &v1, &v2,
    1035             :                         &winv1, &winv2, &vt_a_v1, &vt_a2_v1);
    1036             :   }
    1037          53 :   if (DEBUGLEVEL >= 5)
    1038           0 :     err_printf("Lanczos halted after %ld iterations\n", iter);
    1039          53 :   v1 = F2w_F2Ms_transmul(x, B, nbrow);
    1040          53 :   v2 = F2w_F2Ms_transmul(v0, B, nbrow);
    1041          53 :   x  = shallowconcat(F2w_transpose_F2m(x), F2w_transpose_F2m(v0));
    1042          53 :   v1 = shallowconcat(F2w_transpose_F2m(v1), F2w_transpose_F2m(v2));
    1043          53 :   s0 = gel(F2m_indexrank(x), 2);
    1044          53 :   x = shallowextract(x, s0);
    1045          53 :   v1 = shallowextract(v1, s0);
    1046          53 :   return F2m_mul(x, F2m_ker(v1));
    1047             : }
    1048             : 
    1049             : static GEN
    1050        2892 : F2v_inflate(GEN v, GEN p, long n)
    1051             : {
    1052        2892 :   long i, l = lg(p)-1;
    1053        2892 :   GEN w = zero_F2v(n);
    1054     4341298 :   for (i=1; i<=l; i++)
    1055     4338406 :     if (F2v_coeff(v,i))
    1056     2167619 :       F2v_set(w, p[i]);
    1057        2892 :   return w;
    1058             : }
    1059             : 
    1060             : static GEN
    1061          53 : F2m_inflate(GEN x, GEN p, long n)
    1062        2945 : { pari_APPLY_same(F2v_inflate(gel(x,i), p, n)) }
    1063             : 
    1064             : GEN
    1065         963 : F2Ms_ker(GEN M, long nbrow)
    1066             : {
    1067         963 :   pari_sp av = avma;
    1068         963 :   long nbcol = lg(M)-1;
    1069             :   GEN Mp, R, Rp, p;
    1070         963 :   if (nbrow <= 640)
    1071         910 :     return gerepileupto(av, F2m_ker(F2Ms_to_F2m(M, nbrow)));
    1072          53 :   p = F2Ms_colelim(M, nbrow);
    1073          53 :   Mp = vecpermute(M, p);
    1074             :   do
    1075             :   {
    1076          54 :     R = block_lanczos(Mp, nbrow);
    1077          54 :   } while(!R);
    1078          53 :   Rp = F2m_inflate(R, p, nbcol);
    1079          53 :   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         910 : F2Ms_to_F2m(GEN M, long nrow)
    1103             : {
    1104         910 :   long i, j, l = lg(M);
    1105         910 :   GEN B = cgetg(l, t_MAT);
    1106      151466 :   for(i = 1; i < l; i++)
    1107             :   {
    1108      150556 :     GEN Bi = zero_F2v(nrow), Mi = gel(M,i);
    1109      150556 :     long l = lg(Mi);
    1110     1997205 :     for (j = 1; j < l; j++)
    1111     1846649 :       F2v_set(Bi, Mi[j]);
    1112      150556 :     gel(B, i) = Bi;
    1113             :   }
    1114         910 :   return B;
    1115             : }

Generated by: LCOV version 1.14