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 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 - alglin1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19230-c71492b) Lines: 2255 2419 93.2 %
Date: 2016-07-30 07:10:28 Functions: 203 220 92.3 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000, 2012  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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                         LINEAR ALGEBRA                         **/
      17             : /**                          (first part)                          **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /*******************************************************************/
      24             : /*                                                                 */
      25             : /*                         GEREPILE                                */
      26             : /*                                                                 */
      27             : /*******************************************************************/
      28             : 
      29             : static void
      30          56 : gerepile_mat(pari_sp av, pari_sp tetpil, GEN x, long k, long m, long n, long t)
      31             : {
      32          56 :   pari_sp A, bot = pari_mainstack->bot;
      33             :   long u, i;
      34             :   size_t dec;
      35             : 
      36          56 :   (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
      37             : 
      38        3889 :   for (u=t+1; u<=m; u++)
      39             :   {
      40        3833 :     A = (pari_sp)coeff(x,u,k);
      41        3833 :     if (A < av && A >= bot) coeff(x,u,k) += dec;
      42             :   }
      43        3545 :   for (i=k+1; i<=n; i++)
      44      553993 :     for (u=1; u<=m; u++)
      45             :     {
      46      550504 :       A = (pari_sp)coeff(x,u,i);
      47      550504 :       if (A < av && A >= bot) coeff(x,u,i) += dec;
      48             :     }
      49          56 : }
      50             : 
      51             : static void
      52          56 : gen_gerepile_gauss_ker(GEN x, long k, long t, pari_sp av, void *E, GEN (*copy)(void*, GEN))
      53             : {
      54          56 :   pari_sp tetpil = avma;
      55          56 :   long u,i, n = lg(x)-1, m = n? nbrows(x): 0;
      56             : 
      57          56 :   if (DEBUGMEM > 1) pari_warn(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
      58          56 :   for (u=t+1; u<=m; u++) gcoeff(x,u,k) = copy(E,gcoeff(x,u,k));
      59        3545 :   for (i=k+1; i<=n; i++)
      60        3489 :     for (u=1; u<=m; u++) gcoeff(x,u,i) = copy(E,gcoeff(x,u,i));
      61          56 :   gerepile_mat(av,tetpil,x,k,m,n,t);
      62          56 : }
      63             : 
      64             : /* special gerepile for huge matrices */
      65             : 
      66             : #define COPY(x) {\
      67             :   GEN _t = (x); if (!is_universal_constant(_t)) x = gcopy(_t); \
      68             : }
      69             : 
      70             : INLINE GEN
      71      506368 : _copy(void *E, GEN x)
      72             : {
      73      506368 :   (void) E; COPY(x);
      74      506368 :   return x;
      75             : }
      76             : 
      77             : static void
      78          46 : gerepile_gauss_ker(GEN x, long k, long t, pari_sp av)
      79             : {
      80          46 :   gen_gerepile_gauss_ker(x, k, t, av, NULL, &_copy);
      81          46 : }
      82             : 
      83             : static void
      84          89 : gerepile_gauss(GEN x,long k,long t,pari_sp av, long j, GEN c)
      85             : {
      86          89 :   pari_sp tetpil = avma, A, bot;
      87          89 :   long u,i, n = lg(x)-1, m = n? nbrows(x): 0;
      88             :   size_t dec;
      89             : 
      90          89 :   if (DEBUGMEM > 1) pari_warn(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n);
      91        1599 :   for (u=t+1; u<=m; u++)
      92        1510 :     if (u==j || !c[u]) COPY(gcoeff(x,u,k));
      93        4633 :   for (u=1; u<=m; u++)
      94        4544 :     if (u==j || !c[u])
      95        3446 :       for (i=k+1; i<=n; i++) COPY(gcoeff(x,u,i));
      96             : 
      97          89 :   (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
      98          89 :   bot = pari_mainstack->bot;
      99        1599 :   for (u=t+1; u<=m; u++)
     100        1510 :     if (u==j || !c[u])
     101             :     {
     102        1510 :       A=(pari_sp)coeff(x,u,k);
     103        1510 :       if (A<av && A>=bot) coeff(x,u,k)+=dec;
     104             :     }
     105        4633 :   for (u=1; u<=m; u++)
     106        4544 :     if (u==j || !c[u])
     107      149194 :       for (i=k+1; i<=n; i++)
     108             :       {
     109      145748 :         A=(pari_sp)coeff(x,u,i);
     110      145748 :         if (A<av && A>=bot) coeff(x,u,i)+=dec;
     111             :       }
     112          89 : }
     113             : 
     114             : /*******************************************************************/
     115             : /*                                                                 */
     116             : /*                         GENERIC                                 */
     117             : /*                                                                 */
     118             : /*******************************************************************/
     119             : GEN
     120        1992 : gen_ker(GEN x, long deplin, void *E, const struct bb_field *ff)
     121             : {
     122        1992 :   pari_sp av0 = avma, av, tetpil;
     123             :   GEN y, c, d;
     124             :   long i, j, k, r, t, n, m;
     125             : 
     126        1992 :   n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
     127        1992 :   m=nbrows(x); r=0;
     128        1992 :   x = RgM_shallowcopy(x);
     129        1992 :   c = zero_zv(m);
     130        1992 :   d=new_chunk(n+1);
     131        1992 :   av=avma;
     132       19633 :   for (k=1; k<=n; k++)
     133             :   {
     134      297934 :     for (j=1; j<=m; j++)
     135      291336 :       if (!c[j])
     136             :       {
     137       98816 :         gcoeff(x,j,k) = ff->red(E, gcoeff(x,j,k));
     138       98816 :         if (!ff->equal0(gcoeff(x,j,k))) break;
     139             :       }
     140       17648 :     if (j>m)
     141             :     {
     142        6598 :       if (deplin)
     143             :       {
     144           7 :         GEN c = cgetg(n+1, t_COL), g0 = ff->s(E,0), g1=ff->s(E,1);
     145           7 :         for (i=1; i<k; i++) gel(c,i) = ff->red(E, gcoeff(x,d[i],k));
     146           7 :         gel(c,k) = g1; for (i=k+1; i<=n; i++) gel(c,i) = g0;
     147           7 :         return gerepileupto(av0, c);
     148             :       }
     149        6591 :       r++; d[k]=0;
     150       52934 :       for(j=1; j<k; j++)
     151       46343 :         if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
     152             :     }
     153             :     else
     154             :     {
     155       11050 :       GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
     156       11050 :       c[j] = k; d[k] = j;
     157       11050 :       gcoeff(x,j,k) = ff->s(E,-1);
     158       11050 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
     159      427675 :       for (t=1; t<=m; t++)
     160             :       {
     161      416625 :         if (t==j) continue;
     162             : 
     163      405575 :         piv = ff->red(E,gcoeff(x,t,k));
     164      405575 :         if (ff->equal0(piv)) continue;
     165             : 
     166       64996 :         gcoeff(x,t,k) = ff->s(E,0);
     167      783032 :         for (i=k+1; i<=n; i++)
     168      718036 :            gcoeff(x,t,i) = ff->add(E, gcoeff(x,t,i), ff->mul(E,piv,gcoeff(x,j,i)));
     169       64996 :         if (gc_needed(av,1))
     170          10 :           gen_gerepile_gauss_ker(x,k,t,av,E,ff->red);
     171             :       }
     172             :     }
     173             :   }
     174        1985 :   if (deplin) { avma = av0; return NULL; }
     175             : 
     176        1978 :   tetpil=avma; y=cgetg(r+1,t_MAT);
     177        8569 :   for (j=k=1; j<=r; j++,k++)
     178             :   {
     179        6591 :     GEN C = cgetg(n+1,t_COL);
     180        6591 :     GEN g0 = ff->s(E,0), g1 = ff->s(E,1);
     181        6591 :     gel(y,j) = C; while (d[k]) k++;
     182       52934 :     for (i=1; i<k; i++)
     183       46343 :       if (d[i])
     184             :       {
     185       17736 :         GEN p1=gcoeff(x,d[i],k);
     186       17736 :         gel(C,i) = ff->red(E,p1); gunclone(p1);
     187             :       }
     188             :       else
     189       28607 :         gel(C,i) = g0;
     190        6591 :     gel(C,k) = g1; for (i=k+1; i<=n; i++) gel(C,i) = g0;
     191             :   }
     192        1978 :   return gerepile(av0,tetpil,y);
     193             : }
     194             : 
     195             : GEN
     196         402 : gen_Gauss_pivot(GEN x, long *rr, void *E, const struct bb_field *ff)
     197             : {
     198             :   pari_sp av;
     199             :   GEN c, d;
     200         402 :   long i, j, k, r, t, m, n = lg(x)-1;
     201             : 
     202         402 :   if (!n) { *rr = 0; return NULL; }
     203             : 
     204         402 :   m=nbrows(x); r=0;
     205         402 :   d = cgetg(n+1, t_VECSMALL);
     206         402 :   x = RgM_shallowcopy(x);
     207         402 :   c = zero_zv(m);
     208         402 :   av=avma;
     209        5489 :   for (k=1; k<=n; k++)
     210             :   {
     211       95092 :     for (j=1; j<=m; j++)
     212       93992 :       if (!c[j])
     213             :       {
     214       42382 :         gcoeff(x,j,k) = ff->red(E,gcoeff(x,j,k));
     215       42382 :         if (!ff->equal0(gcoeff(x,j,k))) break;
     216             :       }
     217        5087 :     if (j>m) { r++; d[k]=0; }
     218             :     else
     219             :     {
     220        3987 :       GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
     221        3987 :       GEN g0 = ff->s(E,0);
     222        3987 :       c[j] = k; d[k] = j;
     223        3987 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
     224      106573 :       for (t=1; t<=m; t++)
     225             :       {
     226      102586 :         if (c[t]) continue; /* already a pivot on that line */
     227             : 
     228       55545 :         piv = ff->red(E,gcoeff(x,t,k));
     229       55545 :         if (ff->equal0(piv)) continue;
     230       53165 :         gcoeff(x,t,k) = g0;
     231     1310708 :         for (i=k+1; i<=n; i++)
     232     1257543 :           gcoeff(x,t,i) = ff->add(E,gcoeff(x,t,i), ff->mul(E,piv,gcoeff(x,j,i)));
     233       53165 :         if (gc_needed(av,1))
     234          89 :           gerepile_gauss(x,k,t,av,j,c);
     235             :       }
     236        3987 :       for (i=k; i<=n; i++) gcoeff(x,j,i) = g0; /* dummy */
     237             :     }
     238             :   }
     239         402 :   *rr = r; avma = (pari_sp)d; return d;
     240             : }
     241             : 
     242             : GEN
     243          70 : gen_det(GEN a, void *E, const struct bb_field *ff)
     244             : {
     245          70 :   pari_sp av = avma;
     246          70 :   long i,j,k, s = 1, nbco = lg(a)-1;
     247          70 :   GEN q, x = ff->s(E,1);
     248          70 :   a = RgM_shallowcopy(a);
     249         385 :   for (i=1; i<nbco; i++)
     250             :   {
     251         322 :     for(k=i; k<=nbco; k++)
     252             :     {
     253         322 :       gcoeff(a,k,i) = ff->red(E,gcoeff(a,k,i));
     254         322 :       if (!ff->equal0(gcoeff(a,k,i))) break;
     255             :     }
     256         315 :     if (k > nbco) return gerepileupto(av, gcoeff(a,i,i));
     257         315 :     if (k != i)
     258             :     { /* exchange the lines s.t. k = i */
     259           7 :       for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     260           7 :       s = -s;
     261             :     }
     262         315 :     q = gcoeff(a,i,i);
     263             : 
     264         315 :     x = ff->red(E,ff->mul(E,x,q));
     265         315 :     q = ff->inv(E,q);
     266        1400 :     for (k=i+1; k<=nbco; k++)
     267             :     {
     268        1085 :       GEN m = ff->red(E,gcoeff(a,i,k));
     269        1085 :       if (ff->equal0(m)) continue;
     270             : 
     271         959 :       m = ff->neg(E, ff->mul(E,m, q));
     272        5264 :       for (j=i+1; j<=nbco; j++)
     273             :       {
     274        4305 :         gcoeff(a,j,k) = ff->add(E, gcoeff(a,j,k), ff->mul(E,m,gcoeff(a,j,i)));
     275        4305 :         if (gc_needed(av,1))
     276             :         {
     277           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
     278           0 :           gerepileall(av,4, &a,&x,&q,&m);
     279             :         }
     280             :       }
     281             :     }
     282             :   }
     283          70 :   if (s < 0) x = ff->neg(E,x);
     284          70 :   return gerepileupto(av, ff->red(E,ff->mul(E, x, gcoeff(a,nbco,nbco))));
     285             : }
     286             : 
     287             : INLINE void
     288       89332 : _gen_addmul(GEN b, long k, long i, GEN m, void *E, const struct bb_field *ff)
     289             : {
     290       89332 :   gel(b,i) = ff->red(E,gel(b,i));
     291       89332 :   gel(b,k) = ff->add(E,gel(b,k), ff->mul(E,m, gel(b,i)));
     292       89332 : }
     293             : 
     294             : static GEN
     295        7710 : _gen_get_col(GEN a, GEN b, long li, void *E, const struct bb_field *ff)
     296             : {
     297        7710 :   GEN u = cgetg(li+1,t_COL);
     298        7710 :   pari_sp av = avma;
     299             :   long i, j;
     300             : 
     301        7710 :   gel(u,li) = gerepileupto(av, ff->red(E,ff->mul(E,gel(b,li), gcoeff(a,li,li))));
     302       28770 :   for (i=li-1; i>0; i--)
     303             :   {
     304       21060 :     pari_sp av = avma;
     305       21060 :     GEN m = gel(b,i);
     306       21060 :     for (j=i+1; j<=li; j++) m = ff->add(E,m, ff->neg(E,ff->mul(E,gcoeff(a,i,j), gel(u,j))));
     307       21060 :     m = ff->red(E, m);
     308       21060 :     gel(u,i) = gerepileupto(av, ff->red(E,ff->mul(E,m, gcoeff(a,i,i))));
     309             :   }
     310        7710 :   return u;
     311             : }
     312             : 
     313             : GEN
     314        2849 : gen_Gauss(GEN a, GEN b, void *E, const struct bb_field *ff)
     315             : {
     316             :   long i, j, k, li, bco, aco;
     317        2849 :   GEN u, g0 = ff->s(E,0);
     318        2849 :   pari_sp av = avma;
     319        2849 :   a = RgM_shallowcopy(a);
     320        2849 :   b = RgM_shallowcopy(b);
     321        2849 :   aco = lg(a)-1; bco = lg(b)-1; li = nbrows(a);
     322        7710 :   for (i=1; i<=aco; i++)
     323             :   {
     324             :     GEN invpiv;
     325        7800 :     for (k = i; k <= li; k++)
     326             :     {
     327        7800 :       GEN piv = ff->red(E,gcoeff(a,k,i));
     328        7800 :       if (!ff->equal0(piv)) { gcoeff(a,k,i) = ff->inv(E,piv); break; }
     329          90 :       gcoeff(a,k,i) = g0;
     330             :     }
     331             :     /* found a pivot on line k */
     332        7710 :     if (k > li) return NULL;
     333        7710 :     if (k != i)
     334             :     { /* swap lines so that k = i */
     335          54 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     336          54 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     337             :     }
     338        7710 :     if (i == aco) break;
     339             : 
     340        4861 :     invpiv = gcoeff(a,i,i); /* 1/piv mod p */
     341       15440 :     for (k=i+1; k<=li; k++)
     342             :     {
     343       10579 :       GEN m = ff->red(E,gcoeff(a,k,i)); gcoeff(a,k,i) = g0;
     344       10579 :       if (ff->equal0(m)) continue;
     345             : 
     346        9031 :       m = ff->red(E,ff->neg(E,ff->mul(E,m, invpiv)));
     347        9031 :       for (j=i+1; j<=aco; j++) _gen_addmul(gel(a,j),k,i,m,E,ff);
     348        9031 :       for (j=1  ; j<=bco; j++) _gen_addmul(gel(b,j),k,i,m,E,ff);
     349             :     }
     350        4861 :     if (gc_needed(av,1))
     351             :     {
     352           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_Gauss. i=%ld",i);
     353           0 :       gerepileall(av,2, &a,&b);
     354             :     }
     355             :   }
     356             : 
     357        2849 :   if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
     358        2849 :   u = cgetg(bco+1,t_MAT);
     359        2849 :   for (j=1; j<=bco; j++) gel(u,j) = _gen_get_col(a, gel(b,j), aco, E, ff);
     360        2849 :   return u;
     361             : }
     362             : 
     363             : /* compatible t_MAT * t_COL, lgA = lg(A) = lg(B) > 1, l = lgcols(A) */
     364             : static GEN
     365         967 : gen_matcolmul_i(GEN A, GEN B, ulong lgA, ulong l,
     366             :                 void *E, const struct bb_field *ff)
     367             : {
     368         967 :   GEN C = cgetg(l, t_COL);
     369             :   ulong i;
     370        4484 :   for (i = 1; i < l; i++) {
     371        3517 :     pari_sp av = avma;
     372        3517 :     GEN e = ff->mul(E, gcoeff(A, i, 1), gel(B, 1));
     373             :     ulong k;
     374       16543 :     for(k = 2; k < lgA; k++)
     375       13026 :       e = ff->add(E, e, ff->mul(E, gcoeff(A, i, k), gel(B, k)));
     376        3517 :     gel(C, i) = gerepileupto(av, ff->red(E, e));
     377             :   }
     378         967 :   return C;
     379             : }
     380             : 
     381             : GEN
     382          21 : gen_matcolmul(GEN A, GEN B, void *E, const struct bb_field *ff)
     383             : {
     384          21 :   ulong lgA = lg(A);
     385          21 :   if (lgA != (ulong)lg(B))
     386           0 :     pari_err_OP("operation 'gen_matcolmul'", A, B);
     387          21 :   if (lgA == 1)
     388           0 :     return cgetg(1, t_COL);
     389          21 :   return gen_matcolmul_i(A, B, lgA, lgcols(A), E, ff);
     390             : }
     391             : 
     392             : GEN
     393         327 : gen_matmul(GEN A, GEN B, void *E, const struct bb_field *ff)
     394             : {
     395         327 :   ulong j, l, lgA, lgB = lg(B);
     396             :   GEN C;
     397         327 :   if (lgB == 1)
     398           0 :     return cgetg(1, t_MAT);
     399         327 :   lgA = lg(A);
     400         327 :   if (lgA != (ulong)lgcols(B))
     401           0 :     pari_err_OP("operation 'gen_matmul'", A, B);
     402         327 :   if (lgA == 1)
     403           0 :     return zeromat(0, lgB - 1);
     404         327 :   l = lgcols(A);
     405         327 :   C = cgetg(lgB, t_MAT);
     406        1273 :   for(j = 1; j < lgB; j++)
     407         946 :     gel(C, j) = gen_matcolmul_i(A, gel(B, j), lgA, l, E, ff);
     408         327 :   return C;
     409             : }
     410             : 
     411             : static GEN
     412       50780 : image_from_pivot(GEN x, GEN d, long r)
     413             : {
     414             :   GEN y;
     415             :   long j, k;
     416             : 
     417       50780 :   if (!d) return gcopy(x);
     418             :   /* d left on stack for efficiency */
     419       49380 :   r = lg(x)-1 - r; /* = dim Im(x) */
     420       49380 :   y = cgetg(r+1,t_MAT);
     421      462789 :   for (j=k=1; j<=r; k++)
     422      413409 :     if (d[k]) gel(y,j++) = gcopy(gel(x,k));
     423       49380 :   return y;
     424             : }
     425             : 
     426             : /*******************************************************************/
     427             : /*                                                                 */
     428             : /*                    LINEAR ALGEBRA MODULO P                      */
     429             : /*                                                                 */
     430             : /*******************************************************************/
     431             : 
     432             : static long
     433     1076611 : F2v_find_nonzero(GEN x0, GEN mask0, long l, long m)
     434             : {
     435     1076611 :   ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
     436             :   long i, j;
     437     2499346 :   for (i = 0; i < l; i++)
     438             :   {
     439     2497715 :     e = *x++ & *mask++;
     440     2497715 :     if (e)
     441     1074980 :       for (j = 1; ; j++, e >>= 1) if (e & 1uL) return i*BITS_IN_LONG+j;
     442     6636687 :   }
     443        1631 :   return m+1;
     444             : }
     445             : 
     446             : /* in place, destroy x */
     447             : GEN
     448      130319 : F2m_ker_sp(GEN x, long deplin)
     449             : {
     450             :   GEN y, c, d;
     451             :   long i, j, k, l, r, m, n;
     452             : 
     453      130319 :   n = lg(x)-1;
     454      130319 :   m = mael(x,1,1); r=0;
     455             : 
     456      130319 :   d = cgetg(n+1, t_VECSMALL);
     457      130319 :   c = const_F2v(m); l = lg(c)-1;
     458     1030000 :   for (k=1; k<=n; k++)
     459             :   {
     460      921505 :     GEN xk = gel(x,k);
     461      921505 :     j = F2v_find_nonzero(xk, c, l, m);
     462      921505 :     if (j>m)
     463             :     {
     464      245679 :       if (deplin) {
     465       21824 :         GEN c = zero_F2v(n);
     466       79479 :         for (i=1; i<k; i++)
     467       57655 :           if (F2v_coeff(xk, d[i]))
     468       31582 :             F2v_set(c, i);
     469       21824 :         F2v_set(c, k);
     470       21824 :         return c;
     471             :       }
     472      223855 :       r++; d[k] = 0;
     473             :     }
     474             :     else
     475             :     {
     476      675826 :       F2v_clear(c,j); d[k] = j;
     477      675826 :       F2v_clear(xk, j);
     478    72567913 :       for (i=k+1; i<=n; i++)
     479             :       {
     480    71892087 :         GEN xi = gel(x,i);
     481    71892087 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     482             :       }
     483      675826 :       F2v_set(xk, j);
     484             :     }
     485             :   }
     486      108495 :   if (deplin) return NULL;
     487             : 
     488      108460 :   y = zero_F2m_copy(n,r);
     489      332315 :   for (j=k=1; j<=r; j++,k++)
     490             :   {
     491      223855 :     GEN C = gel(y,j); while (d[k]) k++;
     492    12757902 :     for (i=1; i<k; i++)
     493    12534047 :       if (d[i] && F2m_coeff(x,d[i],k))
     494     4446124 :         F2v_set(C,i);
     495      223855 :     F2v_set(C, k);
     496             :   }
     497      108460 :   return y;
     498             : }
     499             : 
     500             : static void /* assume m < p */
     501   119074147 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
     502             : {
     503   119074147 :   uel(b,k) = Fl_addmul_pre(m, uel(b,i), uel(b,k), p, pi);
     504   119074147 : }
     505             : static void /* same m = 1 */
     506     2119030 : _Fl_add(GEN b, long k, long i, ulong p)
     507             : {
     508     2119030 :   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
     509     2119030 : }
     510             : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     511   228803907 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
     512             : {
     513   228803907 :   uel(b,k) += m * uel(b,i);
     514   228803907 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     515   228803907 : }
     516             : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     517    73938890 : _Fl_add_OK(GEN b, long k, long i, ulong p)
     518             : {
     519    73938890 :   uel(b,k) += uel(b,i);
     520    73938890 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     521    73938890 : }
     522             : 
     523             : static GEN
     524      418135 : Flm_ker_sp_OK(GEN x, ulong p, long deplin)
     525             : {
     526             :   GEN y, c, d;
     527             :   long i, j, k, r, t, m, n;
     528             :   ulong a;
     529             : 
     530      418135 :   n = lg(x)-1;
     531      418135 :   m=nbrows(x); r=0;
     532             : 
     533      418135 :   c = zero_zv(m);
     534      418135 :   d = new_chunk(n+1);
     535      418135 :   a = 0; /* for gcc -Wall */
     536     3070823 :   for (k=1; k<=n; k++)
     537             :   {
     538    24836093 :     for (j=1; j<=m; j++)
     539    23808239 :       if (!c[j])
     540             :       {
     541     9340962 :         a = ucoeff(x,j,k) % p;
     542     9340962 :         if (a) break;
     543             :       }
     544     2682022 :     if (j > m)
     545             :     {
     546     1027854 :       if (deplin) {
     547       29334 :         c = cgetg(n+1, t_VECSMALL);
     548       29334 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
     549       29334 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     550       29334 :         return c;
     551             :       }
     552      998520 :       r++; d[k] = 0;
     553             :     }
     554             :     else
     555             :     {
     556     1654168 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     557     1654168 :       c[j] = k; d[k] = j;
     558     1654168 :       ucoeff(x,j,k) = p-1;
     559     1654168 :       if (piv != 1)
     560     1275723 :         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
     561    55520622 :       for (t=1; t<=m; t++)
     562             :       {
     563    53866454 :         if (t == j) continue;
     564             : 
     565    52212286 :         piv = ( ucoeff(x,t,k) %= p );
     566    52212286 :         if (!piv) continue;
     567    18023882 :         if (piv == 1)
     568     3840932 :           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
     569             :         else
     570    14182950 :           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
     571             :       }
     572             :     }
     573             :   }
     574      388801 :   if (deplin) return NULL;
     575             : 
     576      387030 :   y = cgetg(r+1, t_MAT);
     577     1385550 :   for (j=k=1; j<=r; j++,k++)
     578             :   {
     579      998520 :     GEN C = cgetg(n+1, t_VECSMALL);
     580             : 
     581      998520 :     gel(y,j) = C; while (d[k]) k++;
     582     8758314 :     for (i=1; i<k; i++)
     583     7759794 :       if (d[i])
     584     5362955 :         uel(C,i) = ucoeff(x,d[i],k) % p;
     585             :       else
     586     2396839 :         uel(C,i) = 0UL;
     587      998520 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     588             :   }
     589      387030 :   return y;
     590             : }
     591             : 
     592             : /* in place, destroy x */
     593             : GEN
     594      423970 : Flm_ker_sp(GEN x, ulong p, long deplin)
     595             : {
     596             :   GEN y, c, d;
     597             :   long i, j, k, r, t, m, n;
     598             :   ulong a, pi;
     599      423970 :   if (SMALL_ULONG(p)) return Flm_ker_sp_OK(x, p, deplin);
     600        5835 :   pi = get_Fl_red(p);
     601             : 
     602        5835 :   n = lg(x)-1;
     603        5835 :   m=nbrows(x); r=0;
     604             : 
     605        5835 :   c = zero_zv(m);
     606        5835 :   d = new_chunk(n+1);
     607        5835 :   a = 0; /* for gcc -Wall */
     608       48849 :   for (k=1; k<=n; k++)
     609             :   {
     610      257130 :     for (j=1; j<=m; j++)
     611      240676 :       if (!c[j])
     612             :       {
     613      134303 :         a = ucoeff(x,j,k);
     614      134303 :         if (a) break;
     615             :       }
     616       43021 :     if (j > m)
     617             :     {
     618       16454 :       if (deplin) {
     619           7 :         c = cgetg(n+1, t_VECSMALL);
     620           7 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k);
     621           7 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     622           7 :         return c;
     623             :       }
     624       16447 :       r++; d[k] = 0;
     625             :     }
     626             :     else
     627             :     {
     628       26567 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     629       26567 :       c[j] = k; d[k] = j;
     630       26567 :       ucoeff(x,j,k) = p-1;
     631       26567 :       if (piv != 1)
     632      139231 :         for (i=k+1; i<=n; i++)
     633      113362 :           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
     634      295546 :       for (t=1; t<=m; t++)
     635             :       {
     636      268979 :         if (t == j) continue;
     637             : 
     638      242412 :         piv = ucoeff(x,t,k);
     639      242412 :         if (!piv) continue;
     640      142664 :         if (piv == 1)
     641         521 :           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
     642             :         else
     643      142143 :           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
     644             :       }
     645             :     }
     646             :   }
     647        5828 :   if (deplin) return NULL;
     648             : 
     649        5821 :   y = cgetg(r+1, t_MAT);
     650       22268 :   for (j=k=1; j<=r; j++,k++)
     651             :   {
     652       16447 :     GEN C = cgetg(n+1, t_VECSMALL);
     653             : 
     654       16447 :     gel(y,j) = C; while (d[k]) k++;
     655       86123 :     for (i=1; i<k; i++)
     656       69676 :       if (d[i])
     657       41893 :         uel(C,i) = ucoeff(x,d[i],k);
     658             :       else
     659       27783 :         uel(C,i) = 0UL;
     660       16447 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     661             :   }
     662        5821 :   return y;
     663             : }
     664             : 
     665             : GEN
     666       47670 : FpM_intersect(GEN x, GEN y, GEN p)
     667             : {
     668       47670 :   pari_sp av = avma;
     669       47670 :   long j, lx = lg(x);
     670             :   GEN z;
     671             : 
     672       47670 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     673       47670 :   z = FpM_ker(shallowconcat(x,y), p);
     674       47670 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     675       47670 :   return gerepileupto(av, FpM_mul(x,z,p));
     676             : }
     677             : GEN
     678           0 : Flm_intersect(GEN x, GEN y, ulong p)
     679             : {
     680           0 :   pari_sp av = avma;
     681           0 :   long j, lx = lg(x);
     682             :   GEN z;
     683             : 
     684           0 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     685           0 :   z = Flm_ker(shallowconcat(x,y), p);
     686           0 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     687           0 :   return gerepileupto(av, Flm_mul(x,z,p));
     688             : }
     689             : 
     690             : /* not memory clean */
     691             : GEN
     692        1449 : F2m_ker(GEN x) { return F2m_ker_sp(F2m_copy(x), 0); }
     693             : GEN
     694           0 : F2m_deplin(GEN x) { return F2m_ker_sp(F2m_copy(x), 1); }
     695             : GEN
     696       45341 : Flm_ker(GEN x, ulong p) { return Flm_ker_sp(Flm_copy(x), p, 0); }
     697             : GEN
     698        3192 : Flm_deplin(GEN x, ulong p) { return Flm_ker_sp(Flm_copy(x), p, 1); }
     699             : 
     700             : ulong
     701          28 : F2m_det_sp(GEN x) { return !F2m_ker_sp(x, 1); }
     702             : 
     703             : ulong
     704          14 : F2m_det(GEN x)
     705             : {
     706          14 :   pari_sp av = avma;
     707          14 :   ulong d = F2m_det_sp(F2m_copy(x));
     708          14 :   avma = av; return d;
     709             : }
     710             : 
     711             : /* in place, destroy a, SMALL_ULONG(p) is TRUE */
     712             : static ulong
     713        1407 : Flm_det_sp_OK(GEN a, long nbco, ulong p)
     714             : {
     715        1407 :   long i,j,k, s = 1;
     716        1407 :   ulong q, x = 1;
     717             : 
     718       11123 :   for (i=1; i<nbco; i++)
     719             :   {
     720        9898 :     for(k=i; k<=nbco; k++)
     721             :     {
     722        9877 :       ulong c = ucoeff(a,k,i) % p;
     723        9877 :       ucoeff(a,k,i) = c;
     724        9877 :       if (c) break;
     725             :     }
     726        9737 :     for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
     727        9737 :     if (k > nbco) return ucoeff(a,i,i);
     728        9716 :     if (k != i)
     729             :     { /* exchange the lines s.t. k = i */
     730          77 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     731          77 :       s = -s;
     732             :     }
     733        9716 :     q = ucoeff(a,i,i);
     734             : 
     735        9716 :     if (x & HIGHMASK) x %= p;
     736        9716 :     x *= q;
     737        9716 :     q = Fl_inv(q,p);
     738       48587 :     for (k=i+1; k<=nbco; k++)
     739             :     {
     740       38871 :       ulong m = ucoeff(a,i,k) % p;
     741       38871 :       if (!m) continue;
     742             : 
     743       34265 :       m = p - ((m*q)%p);
     744      197379 :       for (j=i+1; j<=nbco; j++)
     745             :       {
     746      163114 :         ulong c = ucoeff(a,j,k);
     747      163114 :         if (c & HIGHMASK) c %= p;
     748      163114 :         ucoeff(a,j,k) = c  + m*ucoeff(a,j,i);
     749             :       }
     750             :     }
     751             :   }
     752        1386 :   if (x & HIGHMASK) x %= p;
     753        1386 :   q = ucoeff(a,nbco,nbco);
     754        1386 :   if (q & HIGHMASK) q %= p;
     755        1386 :   x = (x*q) % p;
     756        1386 :   if (s < 0 && x) x = p - x;
     757        1386 :   return x;
     758             : }
     759             : /* in place, destroy a */
     760             : ulong
     761       70577 : Flm_det_sp(GEN a, ulong p)
     762             : {
     763       70577 :   long i,j,k, s = 1, nbco = lg(a)-1;
     764       70577 :   ulong pi, q, x = 1;
     765             : 
     766       70577 :   if (SMALL_ULONG(p)) return Flm_det_sp_OK(a, nbco, p);
     767       69170 :   pi = get_Fl_red(p);
     768      520518 :   for (i=1; i<nbco; i++)
     769             :   {
     770      453218 :     for(k=i; k<=nbco; k++)
     771      453218 :       if (ucoeff(a,k,i)) break;
     772      451348 :     if (k > nbco) return ucoeff(a,i,i);
     773      451348 :     if (k != i)
     774             :     { /* exchange the lines s.t. k = i */
     775         668 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     776         668 :       s = -s;
     777             :     }
     778      451348 :     q = ucoeff(a,i,i);
     779             : 
     780      451348 :     x = Fl_mul_pre(x, q, p, pi);
     781      451348 :     q = Fl_inv(q,p);
     782     2235383 :     for (k=i+1; k<=nbco; k++)
     783             :     {
     784     1784035 :       ulong m = ucoeff(a,i,k);
     785     1784035 :       if (!m) continue;
     786             : 
     787     1653697 :       m = Fl_mul_pre(m, q, p, pi);
     788     9942449 :       for (j=i+1; j<=nbco; j++)
     789     8288752 :         ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
     790             :     }
     791             :   }
     792       69170 :   if (s < 0) x = Fl_neg(x, p);
     793       69170 :   return Fl_mul(x, ucoeff(a,nbco,nbco), p);
     794             : }
     795             : 
     796             : ulong
     797       70570 : Flm_det(GEN x, ulong p)
     798             : {
     799       70570 :   pari_sp av = avma;
     800       70570 :   ulong d = Flm_det_sp(Flm_copy(x), p);
     801       70570 :   avma = av; return d;
     802             : }
     803             : 
     804             : static GEN
     805      302004 : FpM_init(GEN a, GEN p, ulong *pp)
     806             : {
     807      302004 :   if (lgefint(p) == 3)
     808             :   {
     809      297475 :     *pp = uel(p,2);
     810      297475 :     return (*pp==2)? ZM_to_F2m(a): ZM_to_Flm(a, *pp);
     811             :   }
     812        4529 :   *pp = 0; return a;
     813             : }
     814             : GEN
     815        2114 : RgM_Fp_init(GEN a, GEN p, ulong *pp)
     816             : {
     817        2114 :   if (lgefint(p) == 3)
     818             :   {
     819        1916 :     *pp = uel(p,2);
     820        1916 :     return (*pp==2)? RgM_to_F2m(a): RgM_to_Flm(a, *pp);
     821             :   }
     822         198 :   *pp = 0; return RgM_to_FpM(a,p);
     823             : }
     824             : 
     825             : static GEN
     826          49 : FpM_det_gen(GEN a, GEN p)
     827             : {
     828             :   void *E;
     829          49 :   const struct bb_field *S = get_Fp_field(&E,p);
     830          49 :   return gen_det(a, E, S);
     831             : }
     832             : GEN
     833          70 : FpM_det(GEN a, GEN p)
     834             : {
     835          70 :   pari_sp av = avma;
     836             :   ulong pp, d;
     837          70 :   a = FpM_init(a, p, &pp);
     838          70 :   switch(pp)
     839             :   {
     840          49 :   case 0: return FpM_det_gen(a, p);
     841          14 :   case 2: d = F2m_det_sp(a); break;
     842           7 :   default:d = Flm_det_sp(a,pp); break;
     843             :   }
     844          21 :   avma = av; return utoi(d);
     845             : }
     846             : 
     847             : /* Destroy x */
     848             : static GEN
     849       25273 : F2m_gauss_pivot(GEN x, long *rr)
     850             : {
     851             :   GEN c, d;
     852             :   long i, j, k, l, r, m, n;
     853             : 
     854       25273 :   n = lg(x)-1; if (!n) { *rr=0; return NULL; }
     855       25273 :   m = mael(x,1,1); r=0;
     856             : 
     857       25273 :   d = cgetg(n+1, t_VECSMALL);
     858       25273 :   c = const_F2v(m); l = lg(c)-1;
     859      180379 :   for (k=1; k<=n; k++)
     860             :   {
     861      155106 :     GEN xk = gel(x,k);
     862      155106 :     j = F2v_find_nonzero(xk, c, l, m);
     863      155106 :     if (j>m) { r++; d[k] = 0; }
     864             :     else
     865             :     {
     866       99929 :       F2v_clear(c,j); d[k] = j;
     867     1599317 :       for (i=k+1; i<=n; i++)
     868             :       {
     869     1499388 :         GEN xi = gel(x,i);
     870     1499388 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     871             :       }
     872             :     }
     873             :   }
     874             : 
     875       25273 :   *rr = r; avma = (pari_sp)d; return d;
     876             : }
     877             : 
     878             : /* Destroy x */
     879             : static GEN
     880      157157 : Flm_gauss_pivot(GEN x, ulong p, long *rr)
     881             : {
     882             :   GEN c,d;
     883             :   long i,j,k,r,t,n,m;
     884             : 
     885      157157 :   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
     886             : 
     887      157157 :   m=nbrows(x); r=0;
     888      157157 :   d=cgetg(n+1,t_VECSMALL);
     889      157157 :   c = zero_zv(m);
     890     1387910 :   for (k=1; k<=n; k++)
     891             :   {
     892    14334846 :     for (j=1; j<=m; j++)
     893    13945702 :       if (!c[j])
     894             :       {
     895     4062235 :         ucoeff(x,j,k) %= p;
     896     4062235 :         if (ucoeff(x,j,k)) break;
     897             :       }
     898     1230753 :     if (j>m) { r++; d[k]=0; }
     899             :     else
     900             :     {
     901      841609 :       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
     902      841609 :       c[j]=k; d[k]=j;
     903    11351663 :       for (i=k+1; i<=n; i++)
     904    10510054 :         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
     905    17177558 :       for (t=1; t<=m; t++)
     906    16335949 :         if (!c[t]) /* no pivot on that line yet */
     907             :         {
     908     9741785 :           piv = ucoeff(x,t,k);
     909     9741785 :           if (piv)
     910             :           {
     911     5011010 :             ucoeff(x,t,k) = 0;
     912   220587801 :             for (i=k+1; i<=n; i++)
     913   431153582 :               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
     914   215576791 :                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
     915             :           }
     916             :         }
     917      841609 :       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
     918             :     }
     919             :   }
     920      157157 :   *rr = r; avma = (pari_sp)d; return d;
     921             : }
     922             : 
     923             : static GEN
     924          59 : FpM_gauss_pivot_gen(GEN x, GEN p, long *rr)
     925             : {
     926             :   void *E;
     927          59 :   const struct bb_field *S = get_Fp_field(&E,p);
     928          59 :   return gen_Gauss_pivot(x, rr, E, S);
     929             : }
     930             : static GEN
     931       74316 : FpM_gauss_pivot(GEN x, GEN p, long *rr)
     932             : {
     933             :   ulong pp;
     934       74316 :   if (lg(x)==1) { *rr = 0; return NULL; }
     935       72951 :   x = FpM_init(x, p, &pp);
     936       72951 :   switch(pp)
     937             :   {
     938          59 :   case 0: return FpM_gauss_pivot_gen(x, p, rr);
     939       25189 :   case 2: return F2m_gauss_pivot(x, rr);
     940       47703 :   default:return Flm_gauss_pivot(x, pp, rr);
     941             :   }
     942             : }
     943             : 
     944             : GEN
     945       50402 : FpM_image(GEN x, GEN p)
     946             : {
     947             :   long r;
     948       50402 :   GEN d = FpM_gauss_pivot(x,p,&r); /* d left on stack for efficiency */
     949       50402 :   return image_from_pivot(x,d,r);
     950             : }
     951             : GEN
     952          28 : Flm_image(GEN x, ulong p)
     953             : {
     954             :   long r;
     955          28 :   GEN d = Flm_gauss_pivot(Flm_copy(x),p,&r); /* d left on stack for efficiency */
     956          28 :   return image_from_pivot(x,d,r);
     957             : }
     958             : GEN
     959           7 : F2m_image(GEN x)
     960             : {
     961             :   long r;
     962           7 :   GEN d = F2m_gauss_pivot(F2m_copy(x),&r); /* d left on stack for efficiency */
     963           7 :   return image_from_pivot(x,d,r);
     964             : }
     965             : 
     966             : long
     967           7 : FpM_rank(GEN x, GEN p)
     968             : {
     969           7 :   pari_sp av = avma;
     970             :   long r;
     971           7 :   (void)FpM_gauss_pivot(x,p,&r);
     972           7 :   avma = av; return lg(x)-1 - r;
     973             : }
     974             : long
     975          14 : Flm_rank(GEN x, ulong p)
     976             : {
     977          14 :   pari_sp av = avma;
     978             :   long r;
     979          14 :   (void)Flm_gauss_pivot(Flm_copy(x),p,&r);
     980          14 :   avma = av; return lg(x)-1 - r;
     981             : }
     982             : long
     983          63 : F2m_rank(GEN x)
     984             : {
     985          63 :   pari_sp av = avma;
     986             :   long r;
     987          63 :   (void)F2m_gauss_pivot(F2m_copy(x),&r);
     988          63 :   avma = av; return lg(x)-1 - r;
     989             : }
     990             : 
     991             : 
     992             : static GEN
     993         265 : FlxqM_gauss_pivot(GEN x, GEN T, ulong p, long *rr)
     994             : {
     995             :   void *E;
     996         265 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
     997         265 :   return gen_Gauss_pivot(x, rr, E, S);
     998             : }
     999             : 
    1000             : GEN
    1001           7 : FlxqM_image(GEN x, GEN T, ulong p)
    1002             : {
    1003             :   long r;
    1004           7 :   GEN d = FlxqM_gauss_pivot(x,T,p,&r); /* d left on stack for efficiency */
    1005           7 :   return image_from_pivot(x,d,r);
    1006             : }
    1007             : 
    1008             : long
    1009         139 : FlxqM_rank(GEN x, GEN T, ulong p)
    1010             : {
    1011         139 :   pari_sp av = avma;
    1012             :   long r;
    1013         139 :   (void)FlxqM_gauss_pivot(x,T,p,&r);
    1014         139 :   avma = av; return lg(x)-1 - r;
    1015             : }
    1016             : GEN
    1017           7 : FlxqM_det(GEN a, GEN T, ulong p)
    1018             : {
    1019             :   void *E;
    1020           7 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1021           7 :   return gen_det(a, E, S);
    1022             : }
    1023             : 
    1024             : GEN
    1025           7 : FlxqM_FlxqC_mul(GEN A, GEN B, GEN T, ulong p) {
    1026             :   void *E;
    1027           7 :   const struct bb_field *ff = get_Flxq_field(&E, T, p);
    1028           7 :   return gen_matcolmul(A, B, E, ff);
    1029             : }
    1030             : 
    1031             : GEN
    1032         142 : FlxqM_mul(GEN A, GEN B, GEN T, ulong p) {
    1033             :   void *E;
    1034             :   const struct bb_field *ff;
    1035         142 :   long n = lg(A) - 1;
    1036             : 
    1037         142 :   if (n == 0)
    1038           0 :     return cgetg(1, t_MAT);
    1039         142 :   if (n > 1)
    1040         142 :     return FlxqM_mul_Kronecker(A, B, T, p);
    1041           0 :   ff = get_Flxq_field(&E, T, p);
    1042           0 :   return gen_matmul(A, B, E, ff);
    1043             : }
    1044             : 
    1045             : GEN
    1046           7 : F2xqM_det(GEN a, GEN T)
    1047             : {
    1048             :   void *E;
    1049           7 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1050           7 :   return gen_det(a, E, S);
    1051             : }
    1052             : 
    1053             : static GEN
    1054          14 : F2xqM_gauss_gen(GEN a, GEN b, GEN T)
    1055             : {
    1056             :   void *E;
    1057          14 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1058          14 :   return gen_Gauss(a, b, E, S);
    1059             : }
    1060             : GEN
    1061          14 : F2xqM_inv(GEN a, GEN T)
    1062             : {
    1063          14 :   pari_sp av = avma;
    1064             :   GEN u;
    1065          14 :   if (lg(a) == 1) { avma = av; return cgetg(1, t_MAT); }
    1066          14 :   u = F2xqM_gauss_gen(a, matid_F2xqM(nbrows(a),T), T);
    1067          14 :   if (!u) { avma = av; return NULL; }
    1068          14 :   return gerepilecopy(av, u);
    1069             : }
    1070             : 
    1071             : GEN
    1072           7 : F2xqM_F2xqC_mul(GEN A, GEN B, GEN T) {
    1073             :   void *E;
    1074           7 :   const struct bb_field *ff = get_F2xq_field(&E, T);
    1075           7 :   return gen_matcolmul(A, B, E, ff);
    1076             : }
    1077             : 
    1078             : GEN
    1079           7 : F2xqM_mul(GEN A, GEN B, GEN T) {
    1080             :   void *E;
    1081           7 :   const struct bb_field *ff = get_F2xq_field(&E, T);
    1082           7 :   return gen_matmul(A, B, E, ff);
    1083             : }
    1084             : 
    1085             : static GEN
    1086          64 : FqM_gauss_pivot_gen(GEN x, GEN T, GEN p, long *rr)
    1087             : {
    1088             :   void *E;
    1089          64 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1090          64 :   return gen_Gauss_pivot(x, rr, E, S);
    1091             : }
    1092             : static GEN
    1093         183 : FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr)
    1094             : {
    1095         183 :   if (lg(x)==1) { *rr = 0; return NULL; }
    1096         183 :   if (!T) return FpM_gauss_pivot(x, p, rr);
    1097         183 :   if (lgefint(p) == 3)
    1098             :   {
    1099         119 :     pari_sp av = avma;
    1100         119 :     ulong pp = uel(p,2);
    1101         119 :     GEN Tp = ZXT_to_FlxT(T, pp);
    1102         119 :     GEN d = FlxqM_gauss_pivot(FqM_to_FlxM(x, T, p), Tp, pp, rr);
    1103         119 :     return d ? gerepileuptoleaf(av, d): d;
    1104             :   }
    1105          64 :   return FqM_gauss_pivot_gen(x, T, p, rr);
    1106             : }
    1107             : 
    1108             : GEN
    1109           7 : FqM_image(GEN x, GEN T, GEN p)
    1110             : {
    1111             :   long r;
    1112           7 :   GEN d = FqM_gauss_pivot(x,T,p,&r); /* d left on stack for efficiency */
    1113           7 :   return image_from_pivot(x,d,r);
    1114             : }
    1115             : 
    1116             : long
    1117          57 : FqM_rank(GEN x, GEN T, GEN p)
    1118             : {
    1119          57 :   pari_sp av = avma;
    1120             :   long r;
    1121          57 :   (void)FqM_gauss_pivot(x,T,p,&r);
    1122          57 :   avma = av; return lg(x)-1 - r;
    1123             : }
    1124             : 
    1125             : GEN
    1126           7 : FqM_det(GEN x, GEN T, GEN p)
    1127             : {
    1128             :   void *E;
    1129           7 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1130           7 :   return gen_det(x, E, S);
    1131             : }
    1132             : 
    1133             : GEN
    1134           7 : FqM_FqC_mul(GEN A, GEN B, GEN T, GEN p) {
    1135             :   void *E;
    1136           7 :   const struct bb_field *ff = get_Fq_field(&E, T, p);
    1137           7 :   return gen_matcolmul(A, B, E, ff);
    1138             : }
    1139             : 
    1140             : GEN
    1141         320 : FqM_mul(GEN A, GEN B, GEN T, GEN p) {
    1142             :   void *E;
    1143         320 :   const struct bb_field *ff = get_Fq_field(&E, T, p);
    1144         320 :   return gen_matmul(A, B, E, ff);
    1145             : }
    1146             : 
    1147             : static GEN
    1148        1642 : FpM_ker_gen(GEN x, GEN p, long deplin)
    1149             : {
    1150             :   void *E;
    1151        1642 :   const struct bb_field *S = get_Fp_field(&E,p);
    1152        1642 :   return gen_ker(x, deplin, E, S);
    1153             : }
    1154             : static GEN
    1155      195816 : FpM_ker_i(GEN x, GEN p, long deplin)
    1156             : {
    1157      195816 :   pari_sp av = avma;
    1158             :   ulong pp;
    1159             :   GEN y;
    1160             : 
    1161      195816 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1162      195711 :   x = FpM_init(x, p, &pp);
    1163      195711 :   switch(pp)
    1164             :   {
    1165        1614 :   case 0: return FpM_ker_gen(x,p,deplin);
    1166             :   case 2:
    1167       65990 :     y = F2m_ker_sp(x, deplin);
    1168       65990 :     if (!y) return y;
    1169       65990 :     y = deplin? F2c_to_ZC(y): F2m_to_ZM(y);
    1170       65990 :     return gerepileupto(av, y);
    1171             :   default:
    1172      128107 :     y = Flm_ker_sp(x, pp, deplin);
    1173      128107 :     if (!y) return y;
    1174      128107 :     y = deplin? Flc_to_ZC(y): Flm_to_ZM(y);
    1175      128107 :     return gerepileupto(av, y);
    1176             :   }
    1177             : }
    1178             : 
    1179             : GEN
    1180      145120 : FpM_ker(GEN x, GEN p) { return FpM_ker_i(x,p,0); }
    1181             : 
    1182             : GEN
    1183       49709 : FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x,p,1); }
    1184             : 
    1185             : static GEN
    1186          15 : FqM_ker_gen(GEN x, GEN T, GEN p, long deplin)
    1187             : {
    1188             :   void *E;
    1189          15 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1190          15 :   return gen_ker(x,deplin,E,S);
    1191             : }
    1192             : static GEN
    1193        1275 : FqM_ker_i(GEN x, GEN T, GEN p, long deplin)
    1194             : {
    1195        1275 :   if (!T) return FpM_ker_i(x,p,deplin);
    1196         288 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1197             : 
    1198         288 :   if (lgefint(p)==3)
    1199             :   {
    1200         273 :     pari_sp ltop=avma;
    1201         273 :     ulong l= p[2];
    1202         273 :     GEN Ml = FqM_to_FlxM(x, T, p);
    1203         273 :     GEN Tl = ZXT_to_FlxT(T,l);
    1204         273 :     GEN p1 = FlxM_to_ZXM(FlxqM_ker(Ml,Tl,l));
    1205         273 :     return gerepileupto(ltop,p1);
    1206             :   }
    1207          15 :   return FqM_ker_gen(x, T, p, deplin);
    1208             : }
    1209             : 
    1210             : GEN
    1211        1275 : FqM_ker(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,0); }
    1212             : 
    1213             : GEN
    1214           0 : FqM_deplin(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,1); }
    1215             : 
    1216             : static GEN
    1217         328 : FlxqM_ker_i(GEN x, GEN T, ulong p, long deplin)
    1218             : {
    1219             :   const struct bb_field *ff;
    1220             :   void *E;
    1221             : 
    1222         328 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1223         328 :   ff=get_Flxq_field(&E,T,p);
    1224         328 :   return gen_ker(x,deplin, E, ff);
    1225             : }
    1226             : 
    1227             : GEN
    1228         328 : FlxqM_ker(GEN x, GEN T, ulong p)
    1229             : {
    1230         328 :   return FlxqM_ker_i(x, T, p, 0);
    1231             : }
    1232             : 
    1233             : static GEN
    1234           7 : F2xqM_ker_i(GEN x, GEN T, long deplin)
    1235             : {
    1236             :   const struct bb_field *ff;
    1237             :   void *E;
    1238             : 
    1239           7 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1240           7 :   ff = get_F2xq_field(&E,T);
    1241           7 :   return gen_ker(x,deplin, E, ff);
    1242             : }
    1243             : 
    1244             : GEN
    1245           7 : F2xqM_ker(GEN x, GEN T)
    1246             : {
    1247           7 :   return F2xqM_ker_i(x, T, 0);
    1248             : }
    1249             : static GEN
    1250          14 : F2xqM_gauss_pivot(GEN x, GEN T, long *rr)
    1251             : {
    1252             :   void *E;
    1253          14 :   const struct bb_field *S = get_F2xq_field(&E,T);
    1254          14 :   return gen_Gauss_pivot(x, rr, E, S);
    1255             : }
    1256             : GEN
    1257           7 : F2xqM_image(GEN x, GEN T)
    1258             : {
    1259             :   long r;
    1260           7 :   GEN d = F2xqM_gauss_pivot(x,T,&r); /* d left on stack for efficiency */
    1261           7 :   return image_from_pivot(x,d,r);
    1262             : }
    1263             : long
    1264           7 : F2xqM_rank(GEN x, GEN T)
    1265             : {
    1266           7 :   pari_sp av = avma;
    1267             :   long r;
    1268           7 :   (void)F2xqM_gauss_pivot(x,T,&r);
    1269           7 :   avma = av; return lg(x)-1 - r;
    1270             : }
    1271             : /*******************************************************************/
    1272             : /*                                                                 */
    1273             : /*                       Solve A*X=B (Gauss pivot)                 */
    1274             : /*                                                                 */
    1275             : /*******************************************************************/
    1276             : /* x ~ 0 compared to reference y */
    1277             : int
    1278      533270 : approx_0(GEN x, GEN y)
    1279             : {
    1280      533270 :   long tx = typ(x);
    1281      533270 :   if (tx == t_COMPLEX)
    1282           0 :     return approx_0(gel(x,1), y) && approx_0(gel(x,2), y);
    1283      533410 :   return gequal0(x) ||
    1284      365983 :          (tx == t_REAL && gexpo(y) - gexpo(x) > bit_prec(x));
    1285             : }
    1286             : /* x a column, x0 same column in the original input matrix (for reference),
    1287             :  * c list of pivots so far */
    1288             : static long
    1289      533508 : gauss_get_pivot_max(GEN X, GEN X0, long ix, GEN c)
    1290             : {
    1291      533508 :   GEN p, r, x = gel(X,ix), x0 = gel(X0,ix);
    1292      533508 :   long i, k = 0, ex = - (long)HIGHEXPOBIT, lx = lg(x);
    1293      533508 :   if (c)
    1294             :   {
    1295       27699 :     for (i=1; i<lx; i++)
    1296       25480 :       if (!c[i])
    1297             :       {
    1298       13580 :         long e = gexpo(gel(x,i));
    1299       13580 :         if (e > ex) { ex = e; k = i; }
    1300             :       }
    1301             :   }
    1302             :   else
    1303             :   {
    1304     1816622 :     for (i=ix; i<lx; i++)
    1305             :     {
    1306     1285333 :       long e = gexpo(gel(x,i));
    1307     1285333 :       if (e > ex) { ex = e; k = i; }
    1308             :     }
    1309             :   }
    1310      533508 :   if (!k) return lx;
    1311      533270 :   p = gel(x,k);
    1312      533270 :   r = gel(x0,k); if (isrationalzero(r)) r = x0;
    1313      533270 :   return approx_0(p, r)? lx: k;
    1314             : }
    1315             : static long
    1316         287 : gauss_get_pivot_padic(GEN X, GEN p, long ix, GEN c)
    1317             : {
    1318         287 :   GEN x = gel(X, ix);
    1319         287 :   long i, k = 0, ex = (long)HIGHVALPBIT, lx = lg(x);
    1320         287 :   if (c)
    1321             :   {
    1322         504 :     for (i=1; i<lx; i++)
    1323         378 :       if (!c[i] && !gequal0(gel(x,i)))
    1324             :       {
    1325         245 :         long e = gvaluation(gel(x,i), p);
    1326         245 :         if (e < ex) { ex = e; k = i; }
    1327             :       }
    1328             :   }
    1329             :   else
    1330             :   {
    1331         588 :     for (i=ix; i<lx; i++)
    1332         427 :       if (!gequal0(gel(x,i)))
    1333             :       {
    1334         399 :         long e = gvaluation(gel(x,i), p);
    1335         399 :         if (e < ex) { ex = e; k = i; }
    1336             :       }
    1337             :   }
    1338         287 :   return k? k: lx;
    1339             : }
    1340             : static long
    1341      310780 : gauss_get_pivot_NZ(GEN X, GEN x0/*unused*/, long ix, GEN c)
    1342             : {
    1343      310780 :   GEN x = gel(X, ix);
    1344      310780 :   long i, lx = lg(x);
    1345             :   (void)x0;
    1346      310780 :   if (c)
    1347             :   {
    1348      197911 :     for (i=1; i<lx; i++)
    1349      187992 :       if (!c[i] && !gequal0(gel(x,i))) return i;
    1350             :   }
    1351             :   else
    1352             :   {
    1353      326006 :     for (i=ix; i<lx; i++)
    1354      325999 :       if (!gequal0(gel(x,i))) return i;
    1355             :   }
    1356        9926 :   return lx;
    1357             : }
    1358             : 
    1359             : /* Return pivot seeking function appropriate for the domain of the RgM x
    1360             :  * (first non zero pivot, maximal pivot...)
    1361             :  * x0 is a reference point used when guessing whether x[i,j] ~ 0
    1362             :  * (iff x[i,j] << x0[i,j]); typical case: mateigen, Gauss pivot on x - vp.Id,
    1363             :  * but use original x when deciding whether a prospective pivot is non-0 */
    1364             : static pivot_fun
    1365      236494 : get_pivot_fun(GEN x, GEN x0, GEN *data)
    1366             : {
    1367      236494 :   long i, j, hx, lx = lg(x);
    1368      236494 :   int res = t_INT;
    1369      236494 :   GEN p = NULL;
    1370             : 
    1371      236494 :   *data = NULL;
    1372      236494 :   if (lx == 1) return &gauss_get_pivot_NZ;
    1373      236354 :   hx = lgcols(x);
    1374     1083467 :   for (j=1; j<lx; j++)
    1375             :   {
    1376      847232 :     GEN xj = gel(x,j);
    1377     5170062 :     for (i=1; i<hx; i++)
    1378             :     {
    1379     4322949 :       GEN c = gel(xj,i);
    1380     4322949 :       switch(typ(c))
    1381             :       {
    1382             :         case t_REAL:
    1383     1422737 :           res = t_REAL;
    1384     1422737 :           break;
    1385             :         case t_COMPLEX:
    1386          14 :           if (typ(gel(c,1)) == t_REAL || typ(gel(c,2)) == t_REAL) res = t_REAL;
    1387          14 :           break;
    1388             :         case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT: case t_QUAD:
    1389             :         case t_POLMOD: /* exact types */
    1390     2899085 :           break;
    1391             :         case t_PADIC:
    1392         994 :           p = gel(c,2);
    1393         994 :           res = t_PADIC;
    1394         994 :           break;
    1395         119 :         default: return &gauss_get_pivot_NZ;
    1396             :       }
    1397             :     }
    1398             :   }
    1399      236235 :   switch(res)
    1400             :   {
    1401      159628 :     case t_REAL: *data = x0; return &gauss_get_pivot_max;
    1402          98 :     case t_PADIC: *data = p; return &gauss_get_pivot_padic;
    1403       76509 :     default: return &gauss_get_pivot_NZ;
    1404             :   }
    1405             : }
    1406             : 
    1407             : static GEN
    1408      497457 : get_col(GEN a, GEN b, GEN p, long li)
    1409             : {
    1410      497457 :   GEN u = cgetg(li+1,t_COL);
    1411             :   long i, j;
    1412             : 
    1413      497457 :   gel(u,li) = gdiv(gel(b,li), p);
    1414     2686836 :   for (i=li-1; i>0; i--)
    1415             :   {
    1416     2189379 :     pari_sp av = avma;
    1417     2189379 :     GEN m = gel(b,i);
    1418     2189379 :     for (j=i+1; j<=li; j++) m = gsub(m, gmul(gcoeff(a,i,j), gel(u,j)));
    1419     2189379 :     gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(a,i,i)));
    1420             :   }
    1421      497457 :   return u;
    1422             : }
    1423             : /* assume 0 <= a[i,j] < p */
    1424             : static GEN
    1425      342117 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
    1426             : {
    1427      342117 :   GEN u = cgetg(li+1,t_VECSMALL);
    1428      342117 :   ulong m = uel(b,li) % p;
    1429             :   long i,j;
    1430             : 
    1431      342117 :   uel(u,li) = (m * ucoeff(a,li,li)) % p;
    1432     5491068 :   for (i = li-1; i > 0; i--)
    1433             :   {
    1434     5148951 :     m = p - uel(b,i)%p;
    1435    63893254 :     for (j = i+1; j <= li; j++) {
    1436    58744303 :       if (m & HIGHBIT) m %= p;
    1437    58744303 :       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
    1438             :     }
    1439     5148951 :     m %= p;
    1440     5148951 :     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
    1441     5148951 :     uel(u,i) = m;
    1442             :   }
    1443      342117 :   return u;
    1444             : }
    1445             : static GEN
    1446     1642643 : Fl_get_col(GEN a, GEN b, long li, ulong p)
    1447             : {
    1448     1642643 :   GEN u = cgetg(li+1,t_VECSMALL);
    1449     1642643 :   ulong m = uel(b,li) % p;
    1450             :   long i,j;
    1451             : 
    1452     1642643 :   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
    1453    17814783 :   for (i=li-1; i>0; i--)
    1454             :   {
    1455    16172140 :     m = b[i]%p;
    1456   217035879 :     for (j = i+1; j <= li; j++)
    1457   200863739 :       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
    1458    16172140 :     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
    1459    16172140 :     uel(u,i) = m;
    1460             :   }
    1461     1642643 :   return u;
    1462             : }
    1463             : 
    1464             : /* bk -= m * bi */
    1465             : static void
    1466     2510726 : _submul(GEN b, long k, long i, GEN m)
    1467             : {
    1468     2510726 :   gel(b,k) = gsub(gel(b,k), gmul(m, gel(b,i)));
    1469     2510726 : }
    1470             : static int
    1471      230831 : init_gauss(GEN a, GEN *b, long *aco, long *li, int *iscol)
    1472             : {
    1473      230831 :   *iscol = *b ? (typ(*b) == t_COL): 0;
    1474      230831 :   *aco = lg(a) - 1;
    1475      230831 :   if (!*aco) /* a empty */
    1476             :   {
    1477         826 :     if (*b && lg(*b) != 1) pari_err_DIM("gauss");
    1478         826 :     *li = 0; return 0;
    1479             :   }
    1480      230005 :   *li = nbrows(a);
    1481      230005 :   if (*li < *aco) pari_err_INV("gauss [no left inverse]", a);
    1482      230005 :   if (*b)
    1483             :   {
    1484      209906 :     if (*li != *aco) pari_err_DIM("gauss");
    1485      209906 :     switch(typ(*b))
    1486             :     {
    1487             :       case t_MAT:
    1488       54057 :         if (lg(*b) == 1) return 0;
    1489       54057 :         *b = RgM_shallowcopy(*b);
    1490       54057 :         break;
    1491             :       case t_COL:
    1492      155849 :         *b = mkmat( leafcopy(*b) );
    1493      155849 :         break;
    1494           0 :       default: pari_err_TYPE("gauss",*b);
    1495             :     }
    1496      209906 :     if (nbrows(*b) != *li) pari_err_DIM("gauss");
    1497             :   }
    1498             :   else
    1499       20099 :     *b = matid(*li);
    1500      230005 :   return 1;
    1501             : }
    1502             : 
    1503             : static int
    1504      321842 : is_modular_solve(GEN a, GEN b, GEN *u)
    1505             : {
    1506      321842 :   GEN p = NULL;
    1507             :   ulong pp;
    1508      321842 :   if (!RgM_is_FpM(a, &p) || !p) return 0;
    1509         168 :   if (!b)
    1510             :   {
    1511          70 :     a = RgM_Fp_init(a, p, &pp);
    1512          70 :     switch(pp)
    1513             :     {
    1514             :     case 0:
    1515          14 :       a = FpM_inv(a,p);
    1516          14 :       if (a) a = FpM_to_mod(a, p);
    1517          14 :       break;
    1518             :     case 2:
    1519          35 :       a = F2m_inv(a);
    1520          35 :       if (a) a = F2m_to_mod(a);
    1521          35 :       break;
    1522             :     default:
    1523          21 :       a = Flm_inv(a,pp);
    1524          21 :       if (a) a = Flm_to_mod(a, pp);
    1525             :     }
    1526             :   }
    1527          98 :   else switch(typ(b))
    1528             :   {
    1529             :     case t_COL:
    1530          49 :       if (!RgV_is_FpV(b, &p)) return 0;
    1531          49 :       a = RgM_Fp_init(a, p, &pp);
    1532          49 :       switch(pp)
    1533             :       {
    1534             :       case 0:
    1535          14 :         b = RgC_to_FpC(b, p);
    1536          14 :         a = FpM_FpC_gauss(a,b,p);
    1537          14 :         if (a) a = FpC_to_mod(a, p);
    1538          14 :         break;
    1539             :       case 2:
    1540          14 :         b = RgV_to_F2v(b);
    1541          14 :         a = F2m_F2c_gauss(a,b);
    1542          14 :         if (a) a = F2c_to_mod(a);
    1543          14 :         break;
    1544             :       default:
    1545          21 :         b = RgV_to_Flv(b, pp);
    1546          21 :         a = Flm_Flc_gauss(a,b,pp);
    1547          21 :         if (a) a = Flc_to_mod(a, pp);
    1548          21 :         break;
    1549             :       }
    1550          49 :       break;
    1551             :     case t_MAT:
    1552          49 :       if (!RgM_is_FpM(b, &p)) return 0;
    1553          49 :       a = RgM_Fp_init(a, p, &pp);
    1554          49 :       switch(pp)
    1555             :       {
    1556             :       case 0:
    1557          14 :         b = RgM_to_FpM(b, p);
    1558          14 :         a = FpM_gauss(a,b,p);
    1559          14 :         if (a) a = FpM_to_mod(a, p);
    1560          14 :         break;
    1561             :       case 2:
    1562          14 :         b = RgM_to_F2m(b);
    1563          14 :         a = F2m_gauss(a,b);
    1564          14 :         if (a) a = F2m_to_mod(a);
    1565          14 :         break;
    1566             :       default:
    1567          21 :         b = RgM_to_Flm(b, pp);
    1568          21 :         a = Flm_gauss(a,b,pp);
    1569          21 :         if (a) a = Flm_to_mod(a, pp);
    1570          21 :         break;
    1571             :       }
    1572          49 :       break;
    1573           0 :     default: return 0;
    1574             :   }
    1575         168 :   *u = a; return 1;
    1576             : }
    1577             : /* Gaussan Elimination. If a is square, return a^(-1)*b;
    1578             :  * if a has more rows than columns and b is NULL, return c such that c a = Id.
    1579             :  * a is a (not necessarily square) matrix
    1580             :  * b is a matrix or column vector, NULL meaning: take the identity matrix,
    1581             :  *   effectively returning the inverse of a
    1582             :  * If a and b are empty, the result is the empty matrix.
    1583             :  *
    1584             :  * li: number of rows of a and b
    1585             :  * aco: number of columns of a
    1586             :  * bco: number of columns of b (if matrix)
    1587             :  */
    1588             : GEN
    1589      321842 : RgM_solve(GEN a, GEN b)
    1590             : {
    1591      321842 :   pari_sp av = avma;
    1592             :   long i, j, k, li, bco, aco;
    1593             :   int iscol;
    1594             :   pivot_fun pivot;
    1595      321842 :   GEN p, u, data, ff = NULL;
    1596             : 
    1597      321842 :   if (is_modular_solve(a,b,&u)) return gerepileupto(av, u);
    1598      321674 :   if (!b && RgM_is_FFM(a, &ff)) return FFM_inv(a, ff);
    1599      321632 :   avma = av;
    1600             : 
    1601      321632 :   if (lg(a)-1 == 2 && nbrows(a) == 2) {
    1602             :     /* 2x2 matrix, start by inverting a */
    1603       92982 :     GEN u = gcoeff(a,1,1), v = gcoeff(a,1,2);
    1604       92982 :     GEN w = gcoeff(a,2,1), x = gcoeff(a,2,2);
    1605       92982 :     GEN D = gsub(gmul(u,x), gmul(v,w)), ainv;
    1606       92982 :     if (gequal0(D)) return NULL;
    1607       92982 :     ainv = mkmat2(mkcol2(x, gneg(w)), mkcol2(gneg(v), u));
    1608       92982 :     ainv = gmul(ainv, ginv(D));
    1609       92982 :     if (b) ainv = gmul(ainv, b);
    1610       92982 :     return gerepileupto(av, ainv);
    1611             :   }
    1612             : 
    1613      228650 :   if (!init_gauss(a, &b, &aco, &li, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
    1614      227894 :   pivot = get_pivot_fun(a, a, &data);
    1615      227894 :   a = RgM_shallowcopy(a);
    1616      227894 :   bco = lg(b)-1;
    1617      227894 :   if(DEBUGLEVEL>4) err_printf("Entering gauss\n");
    1618             : 
    1619      227894 :   p = NULL; /* gcc -Wall */
    1620      807797 :   for (i=1; i<=aco; i++)
    1621             :   {
    1622             :     /* k is the line where we find the pivot */
    1623      807797 :     k = pivot(a, data, i, NULL);
    1624      807797 :     if (k > li) return NULL;
    1625      807790 :     if (k != i)
    1626             :     { /* exchange the lines s.t. k = i */
    1627      122308 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    1628      122308 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    1629             :     }
    1630      807790 :     p = gcoeff(a,i,i);
    1631      807790 :     if (i == aco) break;
    1632             : 
    1633     2167564 :     for (k=i+1; k<=li; k++)
    1634             :     {
    1635     1587661 :       GEN m = gcoeff(a,k,i);
    1636     1587661 :       if (!gequal0(m))
    1637             :       {
    1638      586027 :         m = gdiv(m,p);
    1639      586027 :         for (j=i+1; j<=aco; j++) _submul(gel(a,j),k,i,m);
    1640      586027 :         for (j=1;   j<=bco; j++) _submul(gel(b,j),k,i,m);
    1641             :       }
    1642             :     }
    1643      579903 :     if (gc_needed(av,1))
    1644             :     {
    1645           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gauss. i=%ld",i);
    1646           0 :       gerepileall(av,2, &a,&b);
    1647             :     }
    1648             :   }
    1649             : 
    1650      227887 :   if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
    1651      227887 :   u = cgetg(bco+1,t_MAT);
    1652      227887 :   for (j=1; j<=bco; j++) gel(u,j) = get_col(a,gel(b,j),p,aco);
    1653      227887 :   return gerepilecopy(av, iscol? gel(u,1): u);
    1654             : }
    1655             : 
    1656             : /* assume dim A >= 1, A invertible + upper triangular  */
    1657             : static GEN
    1658       17756 : RgM_inv_upper_ind(GEN A, long index)
    1659             : {
    1660       17756 :   long n = lg(A)-1, i = index, j;
    1661       17756 :   GEN u = zerocol(n);
    1662       17756 :   gel(u,i) = ginv(gcoeff(A,i,i));
    1663       46434 :   for (i--; i>0; i--)
    1664             :   {
    1665       28678 :     pari_sp av = avma;
    1666       28678 :     GEN m = gneg(gmul(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
    1667       28678 :     for (j=i+2; j<=n; j++) m = gsub(m, gmul(gcoeff(A,i,j),gel(u,j)));
    1668       28678 :     gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(A,i,i)));
    1669             :   }
    1670       17756 :   return u;
    1671             : }
    1672             : GEN
    1673        6924 : RgM_inv_upper(GEN A)
    1674             : {
    1675             :   long i, l;
    1676        6924 :   GEN B = cgetg_copy(A, &l);
    1677        6924 :   for (i = 1; i < l; i++) gel(B,i) = RgM_inv_upper_ind(A, i);
    1678        6924 :   return B;
    1679             : }
    1680             : 
    1681             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal  */
    1682             : static GEN
    1683          14 : FpM_inv_upper_1_ind(GEN A, long index, GEN p)
    1684             : {
    1685          14 :   long n = lg(A)-1, i = index, j;
    1686          14 :   GEN u = zerocol(n);
    1687          14 :   gel(u,i) = gen_1;
    1688          21 :   for (i--; i>0; i--)
    1689             :   {
    1690           7 :     pari_sp av = avma;
    1691           7 :     GEN m = negi(mulii(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
    1692           7 :     for (j=i+2; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    1693           7 :     gel(u,i) = gerepileuptoint(av, modii(m,p));
    1694             :   }
    1695          14 :   return u;
    1696             : }
    1697             : static GEN
    1698           7 : FpM_inv_upper_1(GEN A, GEN p)
    1699             : {
    1700             :   long i, l;
    1701           7 :   GEN B = cgetg_copy(A, &l);
    1702           7 :   for (i = 1; i < l; i++) gel(B,i) = FpM_inv_upper_1_ind(A, i, p);
    1703           7 :   return B;
    1704             : }
    1705             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
    1706             :  * reduced mod p */
    1707             : static GEN
    1708          28 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
    1709             : {
    1710          28 :   long n = lg(A)-1, i = index, j;
    1711          28 :   GEN u = const_vecsmall(n, 0);
    1712          28 :   u[i] = 1;
    1713          28 :   if (SMALL_ULONG(p))
    1714          21 :     for (i--; i>0; i--)
    1715             :     {
    1716           7 :       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
    1717           7 :       for (j=i+2; j<=n; j++)
    1718             :       {
    1719           0 :         if (m & HIGHMASK) m %= p;
    1720           0 :         m += ucoeff(A,i,j) * uel(u,j);
    1721             :       }
    1722           7 :       u[i] = Fl_neg(m % p, p);
    1723             :     }
    1724             :   else
    1725          21 :     for (i--; i>0; i--)
    1726             :     {
    1727           7 :       ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
    1728           7 :       for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
    1729           7 :       u[i] = Fl_neg(m, p);
    1730             :     }
    1731          28 :   return u;
    1732             : }
    1733             : static GEN
    1734          14 : F2m_inv_upper_1_ind(GEN A, long index)
    1735             : {
    1736          14 :   pari_sp av = avma;
    1737          14 :   long n = lg(A)-1, i = index, j;
    1738          14 :   GEN u = const_vecsmall(n, 0);
    1739          14 :   u[i] = 1;
    1740          21 :   for (i--; i>0; i--)
    1741             :   {
    1742           7 :     ulong m = F2m_coeff(A,i,i+1) & uel(u,i+1); /* j = i+1 */
    1743           7 :     for (j=i+2; j<=n; j++) m ^= F2m_coeff(A,i,j) & uel(u,j);
    1744           7 :     u[i] = m & 1;
    1745             :   }
    1746          14 :   return gerepileuptoleaf(av, Flv_to_F2v(u));
    1747             : }
    1748             : static GEN
    1749          14 : Flm_inv_upper_1(GEN A, ulong p)
    1750             : {
    1751             :   long i, l;
    1752          14 :   GEN B = cgetg_copy(A, &l);
    1753          14 :   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
    1754          14 :   return B;
    1755             : }
    1756             : static GEN
    1757           7 : F2m_inv_upper_1(GEN A)
    1758             : {
    1759             :   long i, l;
    1760           7 :   GEN B = cgetg_copy(A, &l);
    1761           7 :   for (i = 1; i < l; i++) gel(B,i) = F2m_inv_upper_1_ind(A, i);
    1762           7 :   return B;
    1763             : }
    1764             : 
    1765             : static GEN
    1766      848355 : split_realimag_col(GEN z, long r1, long r2)
    1767             : {
    1768      848355 :   long i, ru = r1+r2;
    1769      848355 :   GEN x = cgetg(ru+r2+1,t_COL), y = x + r2;
    1770     2725001 :   for (i=1; i<=r1; i++) {
    1771     1876646 :     GEN a = gel(z,i);
    1772     1876646 :     if (typ(a) == t_COMPLEX) a = gel(a,1); /* paranoia: a should be real */
    1773     1876646 :     gel(x,i) = a;
    1774             :   }
    1775     1316821 :   for (   ; i<=ru; i++) {
    1776      468466 :     GEN b, a = gel(z,i);
    1777      468466 :     if (typ(a) == t_COMPLEX) { b = gel(a,2); a = gel(a,1); } else b = gen_0;
    1778      468466 :     gel(x,i) = a;
    1779      468466 :     gel(y,i) = b;
    1780             :   }
    1781      848355 :   return x;
    1782             : }
    1783             : GEN
    1784      437011 : split_realimag(GEN x, long r1, long r2)
    1785             : {
    1786             :   long i,l; GEN y;
    1787      437011 :   if (typ(x) == t_COL) return split_realimag_col(x,r1,r2);
    1788      214601 :   y = cgetg_copy(x, &l);
    1789      214601 :   for (i=1; i<l; i++) gel(y,i) = split_realimag_col(gel(x,i), r1, r2);
    1790      214601 :   return y;
    1791             : }
    1792             : 
    1793             : /* assume M = (r1+r2) x (r1+2r2) matrix and y compatible vector or matrix
    1794             :  * r1 first lines of M,y are real. Solve the system obtained by splitting
    1795             :  * real and imaginary parts. */
    1796             : GEN
    1797      212825 : RgM_solve_realimag(GEN M, GEN y)
    1798             : {
    1799      212825 :   long l = lg(M), r2 = l - lgcols(M), r1 = l-1 - 2*r2;
    1800      212825 :   return RgM_solve(split_realimag(M, r1,r2),
    1801             :                    split_realimag(y, r1,r2));
    1802             : }
    1803             : 
    1804             : GEN
    1805         210 : gauss(GEN a, GEN b)
    1806             : {
    1807             :   GEN z;
    1808         210 :   if (typ(a)!=t_MAT) pari_err_TYPE("gauss",a);
    1809         294 :   if (RgM_is_ZM(a) && b &&
    1810         133 :       ((typ(b) == t_COL && RgV_is_ZV(b)) || (typ(b) == t_MAT && RgM_is_ZM(b))))
    1811          84 :     z = ZM_gauss(a,b);
    1812             :   else
    1813         126 :     z = RgM_solve(a,b);
    1814         210 :   if (!z) pari_err_INV("gauss",a);
    1815         182 :   return z;
    1816             : }
    1817             : 
    1818             : static GEN
    1819       66074 : F2_get_col(GEN b, GEN d, long li, long aco)
    1820             : {
    1821       66074 :   long i, l = nbits2lg(aco);
    1822       66074 :   GEN u = cgetg(l, t_VECSMALL);
    1823       66074 :   u[1] = aco;
    1824     1109189 :   for (i = 1; i <= li; i++)
    1825     1043115 :     if (d[i]) /* d[i] can still be 0 if li > aco */
    1826             :     {
    1827     1043080 :       if (F2v_coeff(b, i))
    1828      357125 :         F2v_set(u, d[i]);
    1829             :       else
    1830      685955 :         F2v_clear(u, d[i]);
    1831             :     }
    1832       66074 :   return u;
    1833             : }
    1834             : 
    1835             : /* destroy a, b */
    1836             : static GEN
    1837       12913 : F2m_gauss_sp(GEN a, GEN b)
    1838             : {
    1839       12913 :   long i, j, k, l, li, bco, aco = lg(a)-1;
    1840             :   GEN u, d;
    1841             : 
    1842       12913 :   if (!aco) return cgetg(1,t_MAT);
    1843       12913 :   li = gel(a,1)[1];
    1844       12913 :   d = zero_Flv(li);
    1845       12913 :   bco = lg(b)-1;
    1846       79022 :   for (i=1; i<=aco; i++)
    1847             :   {
    1848       66137 :     GEN ai = vecsmall_copy(gel(a,i));
    1849       66137 :     if (!d[i] && F2v_coeff(ai, i))
    1850       35663 :       k = i;
    1851             :     else
    1852       30474 :       for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
    1853             :     /* found a pivot on row k */
    1854       66137 :     if (k > li) return NULL;
    1855       66109 :     d[k] = i;
    1856             : 
    1857             :     /* Clear k-th row but column-wise instead of line-wise */
    1858             :     /*  a_ij -= a_i1*a1j/a_11
    1859             :        line-wise grouping:  L_j -= a_1j/a_11*L_1
    1860             :        column-wise:         C_i -= a_i1/a_11*C_1
    1861             :     */
    1862       66109 :     F2v_clear(ai,k);
    1863     1109273 :     for (l=1; l<=aco; l++)
    1864             :     {
    1865     1043164 :       GEN al = gel(a,l);
    1866     1043164 :       if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
    1867             :     }
    1868     1109217 :     for (l=1; l<=bco; l++)
    1869             :     {
    1870     1043108 :       GEN bl = gel(b,l);
    1871     1043108 :       if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
    1872             :     }
    1873             :   }
    1874       12885 :   u = cgetg(bco+1,t_MAT);
    1875       12885 :   for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
    1876       12885 :   return u;
    1877             : }
    1878             : 
    1879             : GEN
    1880          28 : F2m_gauss(GEN a, GEN b)
    1881             : {
    1882          28 :   pari_sp av = avma;
    1883          28 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    1884          28 :   return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), F2m_copy(b)));
    1885             : }
    1886             : GEN
    1887          14 : F2m_F2c_gauss(GEN a, GEN b)
    1888             : {
    1889          14 :   pari_sp av = avma;
    1890          14 :   GEN z = F2m_gauss(a, mkmat(b));
    1891          14 :   if (!z) { avma = av; return NULL; }
    1892           7 :   if (lg(z) == 1) { avma = av; return cgetg(1,t_VECSMALL); }
    1893           7 :   return gerepileuptoleaf(av, gel(z,1));
    1894             : }
    1895             : 
    1896             : GEN
    1897          35 : F2m_inv(GEN a)
    1898             : {
    1899          35 :   pari_sp av = avma;
    1900          35 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    1901          35 :   return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), matid_F2m(gel(a,1)[1])));
    1902             : }
    1903             : 
    1904             : /* destroy a, b */
    1905             : static GEN
    1906       32843 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
    1907             : {
    1908       32843 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
    1909       32843 :   ulong det = 1;
    1910             :   GEN u;
    1911             : 
    1912       32843 :   li = nbrows(a);
    1913       32843 :   bco = lg(b)-1;
    1914      342138 :   for (i=1; i<=aco; i++)
    1915             :   {
    1916             :     ulong invpiv;
    1917             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
    1918      342138 :     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
    1919      677925 :     for (k = i; k <= li; k++)
    1920             :     {
    1921      677918 :       ulong piv = ( ucoeff(a,k,i) %= p );
    1922      677918 :       if (piv)
    1923             :       {
    1924      342131 :         ucoeff(a,k,i) = Fl_inv(piv, p);
    1925      342131 :         if (detp)
    1926             :         {
    1927           0 :           if (det & HIGHMASK) det %= p;
    1928           0 :           det *= piv;
    1929             :         }
    1930      342131 :         break;
    1931             :       }
    1932             :     }
    1933             :     /* found a pivot on line k */
    1934      342138 :     if (k > li) return NULL;
    1935      342131 :     if (k != i)
    1936             :     { /* swap lines so that k = i */
    1937      113893 :       s = -s;
    1938      113893 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    1939      113893 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    1940             :     }
    1941      342131 :     if (i == aco) break;
    1942             : 
    1943      309295 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    1944     2883802 :     for (k=i+1; k<=li; k++)
    1945             :     {
    1946     2574507 :       ulong m = ( ucoeff(a,k,i) %= p );
    1947     2574507 :       if (!m) continue;
    1948             : 
    1949      629913 :       m = Fl_mul(m, invpiv, p);
    1950      629913 :       if (m == 1) {
    1951      168650 :         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
    1952      168650 :         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
    1953             :       } else {
    1954      461263 :         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
    1955      461263 :         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
    1956             :       }
    1957             :     }
    1958             :   }
    1959       32836 :   if (detp)
    1960             :   {
    1961           0 :     det %=  p;
    1962           0 :     if (s < 0 && det) det = p - det;
    1963           0 :     *detp = det;
    1964             :   }
    1965       32836 :   u = cgetg(bco+1,t_MAT);
    1966       32836 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
    1967       32836 :   return u;
    1968             : }
    1969             : 
    1970             : /* destroy a, b */
    1971             : static GEN
    1972      250128 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p)
    1973             : {
    1974      250128 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
    1975      250128 :   ulong det = 1;
    1976             :   GEN u;
    1977             :   ulong pi;
    1978      250128 :   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
    1979      250128 :   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
    1980      217285 :   pi = get_Fl_red(p);
    1981      217285 :   li = nbrows(a);
    1982      217285 :   bco = lg(b)-1;
    1983     1642657 :   for (i=1; i<=aco; i++)
    1984             :   {
    1985             :     ulong invpiv;
    1986             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
    1987     2814460 :     for (k = i; k <= li; k++)
    1988             :     {
    1989     2814460 :       ulong piv = ucoeff(a,k,i);
    1990     2814460 :       if (piv)
    1991             :       {
    1992     1642657 :         ucoeff(a,k,i) = Fl_inv(piv, p);
    1993     1642657 :         if (detp) det = Fl_mul_pre(det, piv, p, pi);
    1994     1642657 :         break;
    1995             :       }
    1996             :     }
    1997             :     /* found a pivot on line k */
    1998     1642657 :     if (k > li) return NULL;
    1999     1642657 :     if (k != i)
    2000             :     { /* swap lines so that k = i */
    2001      380869 :       s = -s;
    2002      380869 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    2003      380869 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    2004             :     }
    2005     1642657 :     if (i == aco) break;
    2006             : 
    2007     1425372 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    2008     9511456 :     for (k=i+1; k<=li; k++)
    2009             :     {
    2010     8086084 :       ulong m = ucoeff(a,k,i);
    2011     8086084 :       if (!m) continue;
    2012             : 
    2013     1701251 :       m = Fl_mul_pre(m, invpiv, p, pi);
    2014     1701251 :       if (m == 1) {
    2015       76176 :         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
    2016       76176 :         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
    2017             :       } else {
    2018     1625075 :         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
    2019     1625075 :         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
    2020             :       }
    2021             :     }
    2022             :   }
    2023      217285 :   if (detp)
    2024             :   {
    2025           0 :     if (s < 0 && det) det = p - det;
    2026           0 :     *detp = det;
    2027             :   }
    2028      217285 :   u = cgetg(bco+1,t_MAT);
    2029      217285 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
    2030      217285 :   return u;
    2031             : }
    2032             : 
    2033             : GEN
    2034          42 : Flm_gauss(GEN a, GEN b, ulong p) {
    2035          42 :   return Flm_gauss_sp(RgM_shallowcopy(a), RgM_shallowcopy(b), NULL, p);
    2036             : }
    2037             : static GEN
    2038      232471 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
    2039      232471 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    2040      232471 :   return Flm_gauss_sp(a, matid_Flm(nbrows(a)), detp, p);
    2041             : }
    2042             : GEN
    2043       17329 : Flm_inv(GEN a, ulong p) {
    2044       17329 :   return Flm_inv_sp(RgM_shallowcopy(a), NULL, p);
    2045             : }
    2046             : GEN
    2047          21 : Flm_Flc_gauss(GEN a, GEN b, ulong p) {
    2048          21 :   pari_sp av = avma;
    2049          21 :   GEN z = Flm_gauss(a, mkmat(b), p);
    2050          21 :   if (!z) { avma = av; return NULL; }
    2051          14 :   if (lg(z) == 1) { avma = av; return cgetg(1,t_VECSMALL); }
    2052          14 :   return gerepileuptoleaf(av, gel(z,1));
    2053             : }
    2054             : 
    2055             : static GEN
    2056        2807 : FpM_gauss_gen(GEN a, GEN b, GEN p)
    2057             : {
    2058             :   void *E;
    2059        2807 :   const struct bb_field *S = get_Fp_field(&E,p);
    2060        2807 :   return gen_Gauss(a,b, E, S);
    2061             : }
    2062             : /* a an FpM, lg(a)>1; b an FpM or NULL (replace by identity) */
    2063             : static GEN
    2064       33272 : FpM_gauss_i(GEN a, GEN b, GEN p, ulong *pp)
    2065             : {
    2066       33272 :   long n = nbrows(a);
    2067       33272 :   a = FpM_init(a,p,pp);
    2068       33272 :   switch(*pp)
    2069             :   {
    2070             :   case 0:
    2071        2807 :     if (!b) b = matid(n);
    2072        2807 :     return FpM_gauss_gen(a,b,p);
    2073             :   case 2:
    2074       12850 :     if (b) b = ZM_to_F2m(b); else b = matid_F2m(n);
    2075       12850 :     return F2m_gauss_sp(a,b);
    2076             :   default:
    2077       17615 :     if (b) b = ZM_to_Flm(b, *pp); else b = matid_Flm(n);
    2078       17615 :     return Flm_gauss_sp(a,b, NULL, *pp);
    2079             :   }
    2080             : }
    2081             : GEN
    2082          21 : FpM_gauss(GEN a, GEN b, GEN p)
    2083             : {
    2084          21 :   pari_sp av = avma;
    2085             :   ulong pp;
    2086             :   GEN u;
    2087          21 :   if (lg(a) == 1 || lg(b)==1) return cgetg(1, t_MAT);
    2088          21 :   u = FpM_gauss_i(a, b, p, &pp);
    2089          21 :   if (!u) { avma = av; return NULL; }
    2090          14 :   switch(pp)
    2091             :   {
    2092          14 :   case 0: return gerepilecopy(av, u);
    2093           0 :   case 2:  u = F2m_to_ZM(u); break;
    2094           0 :   default: u = Flm_to_ZM(u); break;
    2095             :   }
    2096           0 :   return gerepileupto(av, u);
    2097             : }
    2098             : GEN
    2099       33223 : FpM_inv(GEN a, GEN p)
    2100             : {
    2101       33223 :   pari_sp av = avma;
    2102             :   ulong pp;
    2103             :   GEN u;
    2104       33223 :   if (lg(a) == 1) return cgetg(1, t_MAT);
    2105       33223 :   u = FpM_gauss_i(a, NULL, p, &pp);
    2106       33223 :   if (!u) { avma = av; return NULL; }
    2107       33223 :   switch(pp)
    2108             :   {
    2109        2779 :   case 0: return gerepilecopy(av, u);
    2110       12829 :   case 2:  u = F2m_to_ZM(u); break;
    2111       17615 :   default: u = Flm_to_ZM(u); break;
    2112             :   }
    2113       30444 :   return gerepileupto(av, u);
    2114             : }
    2115             : 
    2116             : GEN
    2117          28 : FpM_FpC_gauss(GEN a, GEN b, GEN p)
    2118             : {
    2119          28 :   pari_sp av = avma;
    2120             :   ulong pp;
    2121             :   GEN u;
    2122          28 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2123          28 :   u = FpM_gauss_i(a, mkmat(b), p, &pp);
    2124          28 :   if (!u) { avma = av; return NULL; }
    2125          21 :   switch(pp)
    2126             :   {
    2127          14 :   case 0: return gerepilecopy(av, gel(u,1));
    2128           7 :   case 2:  u = F2c_to_ZC(gel(u,1)); break;
    2129           0 :   default: u = Flc_to_ZC(gel(u,1)); break;
    2130             :   }
    2131           7 :   return gerepileupto(av, u);
    2132             : }
    2133             : 
    2134             : static GEN
    2135          14 : FlxqM_gauss_gen(GEN a, GEN b, GEN T, ulong p)
    2136             : {
    2137             :   void *E;
    2138          14 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    2139          14 :   return gen_Gauss(a, b, E, S);
    2140             : }
    2141             : GEN
    2142           0 : FlxqM_gauss(GEN a, GEN b, GEN T, ulong p)
    2143             : {
    2144           0 :   pari_sp av = avma;
    2145           0 :   long n = lg(a)-1;
    2146             :   GEN u;
    2147           0 :   if (!n || lg(b)==1) { avma = av; return cgetg(1, t_MAT); }
    2148           0 :   u = FlxqM_gauss_gen(a, b, T, p);
    2149           0 :   if (!u) { avma = av; return NULL; }
    2150           0 :   return gerepilecopy(av, u);
    2151             : }
    2152             : GEN
    2153          14 : FlxqM_inv(GEN a, GEN T, ulong p)
    2154             : {
    2155          14 :   pari_sp av = avma;
    2156             :   GEN u;
    2157          14 :   if (lg(a) == 1) { avma = av; return cgetg(1, t_MAT); }
    2158          14 :   u = FlxqM_gauss_gen(a, matid_FlxqM(nbrows(a),T,p), T,p);
    2159          14 :   if (!u) { avma = av; return NULL; }
    2160          14 :   return gerepilecopy(av, u);
    2161             : }
    2162             : GEN
    2163           0 : FlxqM_FlxqC_gauss(GEN a, GEN b, GEN T, ulong p)
    2164             : {
    2165           0 :   pari_sp av = avma;
    2166             :   GEN u;
    2167           0 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2168           0 :   u = FlxqM_gauss_gen(a, mkmat(b), T, p);
    2169           0 :   if (!u) { avma = av; return NULL; }
    2170           0 :   return gerepilecopy(av, gel(u,1));
    2171             : }
    2172             : 
    2173             : static GEN
    2174          14 : FqM_gauss_gen(GEN a, GEN b, GEN T, GEN p)
    2175             : {
    2176             :   void *E;
    2177          14 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    2178          14 :   return gen_Gauss(a,b,E,S);
    2179             : }
    2180             : GEN
    2181           7 : FqM_gauss(GEN a, GEN b, GEN T, GEN p)
    2182             : {
    2183           7 :   pari_sp av = avma;
    2184             :   GEN u;
    2185             :   long n;
    2186           7 :   if (!T) return FpM_gauss(a,b,p);
    2187           0 :   n = lg(a)-1; if (!n || lg(b)==1) return cgetg(1, t_MAT);
    2188           0 :   u = FqM_gauss_gen(a,b,T,p);
    2189           0 :   if (!u) { avma = av; return NULL; }
    2190           0 :   return gerepilecopy(av, u);
    2191             : }
    2192             : GEN
    2193          14 : FqM_inv(GEN a, GEN T, GEN p)
    2194             : {
    2195          14 :   pari_sp av = avma;
    2196             :   GEN u;
    2197          14 :   if (!T) return FpM_inv(a,p);
    2198          14 :   if (lg(a) == 1) return cgetg(1, t_MAT);
    2199          14 :   u = FqM_gauss_gen(a,matid(nbrows(a)),T,p);
    2200          14 :   if (!u) { avma = av; return NULL; }
    2201          14 :   return gerepilecopy(av, u);
    2202             : }
    2203             : GEN
    2204          14 : FqM_FqC_gauss(GEN a, GEN b, GEN T, GEN p)
    2205             : {
    2206          14 :   pari_sp av = avma;
    2207             :   GEN u;
    2208          14 :   if (!T) return FpM_FpC_gauss(a,b,p);
    2209           0 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2210           0 :   u = FqM_gauss_gen(a,mkmat(b),T,p);
    2211           0 :   if (!u) { avma = av; return NULL; }
    2212           0 :   return gerepilecopy(av, gel(u,1));
    2213             : }
    2214             : 
    2215             : /* Dixon p-adic lifting algorithm.
    2216             :  * Numer. Math. 40, 137-141 (1982), DOI: 10.1007/BF01459082 */
    2217             : GEN
    2218        2181 : ZM_gauss(GEN a, GEN b0)
    2219             : {
    2220        2181 :   pari_sp av = avma, av2;
    2221             :   int iscol;
    2222             :   long n, ncol, i, m, elim;
    2223             :   ulong p;
    2224        2181 :   GEN N, C, delta, xb, nb, nmin, res, b = b0;
    2225             :   forprime_t S;
    2226             : 
    2227        2181 :   if (!init_gauss(a, &b, &n, &ncol, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
    2228        2111 :   nb = gen_0; ncol = lg(b);
    2229       13048 :   for (i = 1; i < ncol; i++)
    2230             :   {
    2231       10937 :     GEN ni = gnorml2(gel(b, i));
    2232       10937 :     if (cmpii(nb, ni) < 0) nb = ni;
    2233             :   }
    2234        2111 :   if (!signe(nb)) { avma = av; return gcopy(b0); }
    2235        2111 :   delta = gen_1; nmin = nb;
    2236       14064 :   for (i = 1; i <= n; i++)
    2237             :   {
    2238       11953 :     GEN ni = gnorml2(gel(a, i));
    2239       11953 :     if (cmpii(ni, nmin) < 0)
    2240             :     {
    2241        1453 :       delta = mulii(delta, nmin); nmin = ni;
    2242             :     }
    2243             :     else
    2244       10500 :       delta = mulii(delta, ni);
    2245             :   }
    2246        2111 :   if (!signe(nmin)) return NULL;
    2247        2097 :   elim = expi(delta)+1;
    2248        2097 :   av2 = avma;
    2249        2097 :   init_modular_big(&S);
    2250             :   for(;;)
    2251             :   {
    2252        2097 :     p = u_forprime_next(&S);
    2253        2097 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    2254        2097 :     if (C) break;
    2255           0 :     elim -= expu(p);
    2256           0 :     if (elim < 0) return NULL;
    2257           0 :     avma = av2;
    2258           0 :   }
    2259             :   /* N.B. Our delta/lambda are SQUARES of those in the paper
    2260             :    * log(delta lambda) / log p, where lambda is 3+sqrt(5) / 2,
    2261             :    * whose log is < 1, hence + 1 (to cater for rounding errors) */
    2262        4194 :   m = (long)ceil((rtodbl(logr_abs(itor(delta,LOWDEFAULTPREC))) + 1)
    2263        2097 :                  / log((double)p));
    2264        2097 :   xb = ZlM_gauss(a, b, p, m, C);
    2265        2097 :   N = powuu(p, m);
    2266        2097 :   delta = sqrti(delta);
    2267        2097 :   if (iscol)
    2268           0 :     res = FpC_ratlift(gel(xb,1), N, delta,delta, NULL);
    2269             :   else
    2270        2097 :     res = FpM_ratlift(xb, N, delta,delta, NULL);
    2271        2097 :   return gerepileupto(av, res);
    2272             : }
    2273             : 
    2274             : /* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */
    2275             : GEN
    2276       99227 : ZM_inv(GEN M, GEN dM)
    2277             : {
    2278       99227 :   pari_sp av2, av = avma;
    2279             :   GEN Hp,q,H;
    2280             :   ulong p;
    2281       99227 :   long lM = lg(M), stable = 0;
    2282       99227 :   int negate = 0;
    2283             :   forprime_t S;
    2284             :   pari_timer ti;
    2285             : 
    2286       99227 :   if (lM == 1) return cgetg(1,t_MAT);
    2287             : 
    2288             :   /* HACK: include dM = -1 ! */
    2289       99122 :   if (dM && is_pm1(dM))
    2290             :   {
    2291             :     /* modular algorithm computes M^{-1}, NOT multiplied by det(M) = -1.
    2292             :      * We will correct (negate) at the end. */
    2293       89007 :     if (signe(dM) < 0) negate = 1;
    2294       89007 :     dM = gen_1;
    2295             :   }
    2296       99122 :   init_modular_big(&S);
    2297       99122 :   av2 = avma;
    2298       99122 :   H = NULL;
    2299       99122 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2300      307536 :   while ((p = u_forprime_next(&S)))
    2301             :   {
    2302             :     ulong dMp;
    2303             :     GEN Mp;
    2304      208414 :     Mp = ZM_to_Flm(M,p);
    2305      208414 :     if (dM == gen_1)
    2306      180595 :       Hp = Flm_inv_sp(Mp, NULL, p);
    2307             :     else
    2308             :     {
    2309       27819 :       if (dM) {
    2310       27819 :         dMp = umodiu(dM,p); if (!dMp) continue;
    2311       27819 :         Hp = Flm_inv_sp(Mp, NULL, p);
    2312       27819 :         if (!Hp) pari_err_INV("ZM_inv", Mp);
    2313             :       } else {
    2314           0 :         Hp = Flm_inv_sp(Mp, &dMp, p);
    2315           0 :         if (!Hp) continue;
    2316             :       }
    2317       27819 :       if (dMp != 1) Flm_Fl_mul_inplace(Hp, dMp, p);
    2318             :     }
    2319             : 
    2320      208414 :     if (!H)
    2321             :     {
    2322       99122 :       H = ZM_init_CRT(Hp, p);
    2323       99122 :       q = utoipos(p);
    2324             :     }
    2325             :     else
    2326      109292 :       stable = ZM_incremental_CRT(&H, Hp, &q, p);
    2327      208414 :     if (DEBUGLEVEL>5) timer_printf(&ti, "ZM_inv mod %lu (stable=%ld)", p,stable);
    2328      208414 :     if (stable) {/* DONE ? */
    2329       99122 :       if (dM != gen_1)
    2330      109237 :       { if (ZM_isscalar(ZM_mul(M, H), dM)) break; }
    2331             :       else
    2332       89007 :       { if (ZM_isidentity(ZM_mul(M, H))) break; }
    2333             :     }
    2334             : 
    2335      109292 :     if (gc_needed(av,2))
    2336             :     {
    2337          30 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv");
    2338          30 :       gerepileall(av2, 2, &H, &q);
    2339             :     }
    2340             :   }
    2341       99122 :   if (!p) pari_err_OVERFLOW("ZM_inv [ran out of primes]");
    2342       99122 :   if (DEBUGLEVEL>5) err_printf("ZM_inv done\n");
    2343       99122 :   if (negate)
    2344           0 :     return gerepileupto(av, ZM_neg(H));
    2345             :   else
    2346       99122 :     return gerepilecopy(av, H);
    2347             : }
    2348             : 
    2349             : /* to be used when denom(M^(-1)) << det(M) and a sharp multiple is
    2350             :  * not available. Return H primitive such that M*H = den*Id */
    2351             : GEN
    2352        6643 : ZM_inv_ratlift(GEN M, GEN *pden)
    2353             : {
    2354        6643 :   pari_sp av2, av = avma;
    2355             :   GEN Hp, q, H;
    2356             :   ulong p;
    2357        6643 :   long lM = lg(M);
    2358             :   forprime_t S;
    2359        6643 :   if (lM == 1) { *pden = gen_1; return cgetg(1,t_MAT); }
    2360             : 
    2361        5964 :   init_modular_big(&S);
    2362        5964 :   av2 = avma;
    2363        5964 :   H = NULL;
    2364       12692 :   while ((p = u_forprime_next(&S)))
    2365             :   {
    2366             :     GEN Mp, B, Hr;
    2367        6728 :     Mp = ZM_to_Flm(M,p);
    2368        6728 :     Hp = Flm_inv_sp(Mp, NULL, p);
    2369        6728 :     if (!Hp) continue;
    2370        6728 :     if (!H)
    2371             :     {
    2372        5964 :       H = ZM_init_CRT(Hp, p);
    2373        5964 :       q = utoipos(p);
    2374             :     }
    2375             :     else
    2376         764 :       ZM_incremental_CRT(&H, Hp, &q, p);
    2377        6728 :     B = sqrti(shifti(q,-1));
    2378        6728 :     Hr = FpM_ratlift(H,q,B,B,NULL);
    2379        6728 :     if (DEBUGLEVEL>5) err_printf("ZM_inv mod %lu (ratlift=%ld)\n", p,!!Hr);
    2380        6728 :     if (Hr) {/* DONE ? */
    2381        6027 :       GEN Hl = Q_remove_denom(Hr, pden);
    2382        6027 :       if (*pden)
    2383        4676 :       { if (ZM_isscalar(ZM_mul(M, Hl), *pden)) { H = Hl; break; }}
    2384             :       else
    2385        1351 :       { if (ZM_isidentity(ZM_mul(M, Hl))) { H = Hl; *pden = gen_1; break; } }
    2386             :     }
    2387             : 
    2388         764 :     if (gc_needed(av,2))
    2389             :     {
    2390           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv_ratlift");
    2391           0 :       gerepileall(av2, 2, &H, &q);
    2392             :     }
    2393             :   }
    2394        5964 :   gerepileall(av, 2, &H, pden);
    2395        5964 :   return H;
    2396             : }
    2397             : 
    2398             : /* same as above, M rational */
    2399             : GEN
    2400        5971 : QM_inv(GEN M, GEN dM)
    2401             : {
    2402        5971 :   pari_sp av = avma;
    2403        5971 :   GEN cM, pM = Q_primitive_part(M, &cM);
    2404        5971 :   if (!cM) return ZM_inv(pM,dM);
    2405        2429 :   return gerepileupto(av, ZM_inv(pM, gdiv(dM,cM)));
    2406             : }
    2407             : 
    2408             : /* x a ZM. Return a multiple of the determinant of the lattice generated by
    2409             :  * the columns of x. From Algorithm 2.2.6 in GTM138 */
    2410             : GEN
    2411          42 : detint(GEN A)
    2412             : {
    2413          42 :   if (typ(A) != t_MAT) pari_err_TYPE("detint",A);
    2414          42 :   RgM_check_ZM(A, "detint");
    2415          42 :   return ZM_detmult(A);
    2416             : }
    2417             : GEN
    2418      226936 : ZM_detmult(GEN A)
    2419             : {
    2420      226936 :   pari_sp av1, av = avma;
    2421             :   GEN B, c, v, piv;
    2422      226936 :   long rg, i, j, k, m, n = lg(A) - 1;
    2423             : 
    2424      226936 :   if (!n) return gen_1;
    2425      226936 :   m = nbrows(A);
    2426      226936 :   if (n < m) return gen_0;
    2427      226915 :   c = zero_zv(m);
    2428      226915 :   av1 = avma;
    2429      226915 :   B = zeromatcopy(m,m);
    2430      226915 :   v = cgetg(m+1, t_COL);
    2431      226915 :   piv = gen_1; rg = 0;
    2432      668994 :   for (k=1; k<=n; k++)
    2433             :   {
    2434      668987 :     GEN pivprec = piv;
    2435      668987 :     long t = 0;
    2436     3028325 :     for (i=1; i<=m; i++)
    2437             :     {
    2438     2359338 :       pari_sp av2 = avma;
    2439             :       GEN vi;
    2440     2359338 :       if (c[i]) continue;
    2441             : 
    2442     1514215 :       vi = mulii(piv, gcoeff(A,i,k));
    2443     8390876 :       for (j=1; j<=m; j++)
    2444     6876661 :         if (c[j]) vi = addii(vi, mulii(gcoeff(B,j,i),gcoeff(A,j,k)));
    2445     1514215 :       if (!t && signe(vi)) t = i;
    2446     1514215 :       gel(v,i) = gerepileuptoint(av2, vi);
    2447             :     }
    2448      668987 :     if (!t) continue;
    2449             :     /* at this point c[t] = 0 */
    2450             : 
    2451      668959 :     if (++rg >= m) { /* full rank; mostly done */
    2452      226908 :       GEN det = gel(v,t); /* last on stack */
    2453      226908 :       if (++k > n)
    2454      226838 :         det = absi(det);
    2455             :       else
    2456             :       {
    2457             :         /* improve further; at this point c[i] is set for all i != t */
    2458          70 :         gcoeff(B,t,t) = piv; v = centermod(gel(B,t), det);
    2459         182 :         for ( ; k<=n; k++)
    2460         112 :           det = gcdii(det, ZV_dotproduct(v, gel(A,k)));
    2461             :       }
    2462      226908 :       return gerepileuptoint(av, det);
    2463             :     }
    2464             : 
    2465      442051 :     piv = gel(v,t);
    2466     2132395 :     for (i=1; i<=m; i++)
    2467             :     {
    2468             :       GEN mvi;
    2469     1690344 :       if (c[i] || i == t) continue;
    2470             : 
    2471      845172 :       gcoeff(B,t,i) = mvi = negi(gel(v,i));
    2472     5362089 :       for (j=1; j<=m; j++)
    2473     4516917 :         if (c[j]) /* implies j != t */
    2474             :         {
    2475      942191 :           pari_sp av2 = avma;
    2476      942191 :           GEN z = addii(mulii(gcoeff(B,j,i), piv), mulii(gcoeff(B,j,t), mvi));
    2477      942191 :           if (rg > 1) z = diviiexact(z, pivprec);
    2478      942191 :           gcoeff(B,j,i) = gerepileuptoint(av2, z);
    2479             :         }
    2480             :     }
    2481      442051 :     c[t] = k;
    2482      442051 :     if (gc_needed(av,1))
    2483             :     {
    2484           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"detint. k=%ld",k);
    2485           0 :       gerepileall(av1, 2, &piv,&B); v = zerovec(m);
    2486             :     }
    2487             :   }
    2488           7 :   avma = av; return gen_0;
    2489             : }
    2490             : 
    2491             : /* Reduce x modulo (invertible) y */
    2492             : GEN
    2493        6111 : closemodinvertible(GEN x, GEN y)
    2494             : {
    2495        6111 :   return gmul(y, ground(RgM_solve(y,x)));
    2496             : }
    2497             : GEN
    2498           7 : reducemodinvertible(GEN x, GEN y)
    2499             : {
    2500           7 :   return gsub(x, closemodinvertible(x,y));
    2501             : }
    2502             : GEN
    2503           0 : reducemodlll(GEN x,GEN y)
    2504             : {
    2505           0 :   return reducemodinvertible(x, ZM_lll(y, 0.75, LLL_INPLACE));
    2506             : }
    2507             : 
    2508             : /*******************************************************************/
    2509             : /*                                                                 */
    2510             : /*                    KERNEL of an m x n matrix                    */
    2511             : /*          return n - rk(x) linearly independent vectors          */
    2512             : /*                                                                 */
    2513             : /*******************************************************************/
    2514             : /* x has INTEGER coefficients. Gauss-Bareiss */
    2515             : GEN
    2516        2395 : keri(GEN x)
    2517             : {
    2518             :   pari_sp av, av0;
    2519             :   GEN c, l, y, p, pp;
    2520             :   long i, j, k, r, t, n, m;
    2521             : 
    2522        2395 :   n = lg(x)-1; if (!n) return cgetg(1,t_MAT);
    2523        2395 :   av0 = avma; m = nbrows(x);
    2524        2395 :   pp = cgetg(n+1,t_COL);
    2525        2395 :   x = RgM_shallowcopy(x);
    2526        2395 :   c = zero_zv(m);
    2527        2395 :   l = cgetg(n+1, t_VECSMALL);
    2528        2395 :   av = avma;
    2529       37356 :   for (r=0, p=gen_1, k=1; k<=n; k++)
    2530             :   {
    2531       34961 :     j = 1;
    2532       34961 :     while ( j <= m && (c[j] || !signe(gcoeff(x,j,k))) ) j++;
    2533       34961 :     if (j > m)
    2534             :     {
    2535       15072 :       r++; l[k] = 0;
    2536      852813 :       for(j=1; j<k; j++)
    2537      837741 :         if (l[j]) gcoeff(x,l[j],k) = gclone(gcoeff(x,l[j],k));
    2538       15072 :       gel(pp,k) = gclone(p);
    2539             :     }
    2540             :     else
    2541             :     {
    2542       19889 :       GEN p0 = p;
    2543       19889 :       c[j] = k; l[k] = j; p = gcoeff(x,j,k);
    2544     1652918 :       for (t=1; t<=m; t++)
    2545     1633029 :         if (t!=j)
    2546             :         {
    2547     1613140 :           GEN q = gcoeff(x,t,k);
    2548   167353595 :           for (i=k+1; i<=n; i++)
    2549             :           {
    2550   165740455 :             pari_sp av1 = avma;
    2551   165740455 :             GEN p1 = subii(mulii(p,gcoeff(x,t,i)), mulii(q,gcoeff(x,j,i)));
    2552   165740455 :             gcoeff(x,t,i) = gerepileuptoint(av1, diviiexact(p1,p0));
    2553             :           }
    2554     1613140 :           if (gc_needed(av,1))
    2555             :           {
    2556          40 :             GEN _p0 = gclone(p0);
    2557          40 :             GEN _p  = gclone(p);
    2558          40 :             gerepile_gauss_ker(x,k,t,av);
    2559          40 :             p = icopy(_p);  gunclone(_p);
    2560          40 :             p0= icopy(_p0); gunclone(_p0);
    2561             :           }
    2562             :         }
    2563             :     }
    2564             :   }
    2565        2395 :   if (!r) { avma = av0; y = cgetg(1,t_MAT); return y; }
    2566             : 
    2567             :   /* non trivial kernel */
    2568        2388 :   y = cgetg(r+1,t_MAT);
    2569       17460 :   for (j=k=1; j<=r; j++,k++)
    2570             :   {
    2571       15072 :     p = cgetg(n+1, t_COL);
    2572       15072 :     gel(y,j) = p; while (l[k]) k++;
    2573      852813 :     for (i=1; i<k; i++)
    2574      837741 :       if (l[i])
    2575             :       {
    2576      535243 :         c = gcoeff(x,l[i],k);
    2577      535243 :         gel(p,i) = icopy(c); gunclone(c);
    2578             :       }
    2579             :       else
    2580      302498 :         gel(p,i) = gen_0;
    2581       15072 :     gel(p,k) = negi(gel(pp,k)); gunclone(gel(pp,k));
    2582       15072 :     for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
    2583             :   }
    2584        2388 :   return gerepileupto(av0, y);
    2585             : }
    2586             : 
    2587             : static GEN
    2588          98 : deplin_aux(GEN x0)
    2589             : {
    2590          98 :   pari_sp av = avma;
    2591          98 :   long i, j, k, nl, nc = lg(x0)-1;
    2592             :   GEN D, x, y, c, l, d, ck;
    2593             : 
    2594          98 :   if (!nc) { avma=av; return cgetg(1,t_COL); }
    2595          63 :   x = RgM_shallowcopy(x0);
    2596          63 :   nl = nbrows(x);
    2597          63 :   d = const_vec(nl, gen_1); /* pivot list */
    2598          63 :   c = zero_zv(nl);
    2599          63 :   l = cgetg(nc+1, t_VECSMALL); /* not initialized */
    2600          63 :   ck = NULL; /* gcc -Wall */
    2601         308 :   for (k=1; k<=nc; k++)
    2602             :   {
    2603         294 :     ck = gel(x,k);
    2604        1008 :     for (j=1; j<k; j++)
    2605             :     {
    2606         714 :       GEN cj = gel(x,j), piv = gel(d,j), q = gel(ck,l[j]);
    2607        9492 :       for (i=1; i<=nl; i++)
    2608        8778 :         if (i!=l[j]) gel(ck,i) = gsub(gmul(piv, gel(ck,i)), gmul(q, gel(cj,i)));
    2609             :     }
    2610             : 
    2611         294 :     i = gauss_get_pivot_NZ(x, NULL, k, c);
    2612         294 :     if (i > nl) break;
    2613             : 
    2614         245 :     gel(d,k) = gel(ck,i);
    2615         245 :     c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */
    2616             :   }
    2617          63 :   if (k > nc) { avma = av; return cgetg(1,t_COL); }
    2618          49 :   if (k == 1) { avma = av; return scalarcol_shallow(gen_1,nc); }
    2619          49 :   y = cgetg(nc+1,t_COL);
    2620          49 :   gel(y,1) = gcopy(gel(ck, l[1]));
    2621         161 :   for (D=gel(d,1),j=2; j<k; j++)
    2622             :   {
    2623         112 :     gel(y,j) = gmul(gel(ck, l[j]), D);
    2624         112 :     D = gmul(D, gel(d,j));
    2625             :   }
    2626          49 :   gel(y,j) = gneg(D);
    2627          49 :   for (j++; j<=nc; j++) gel(y,j) = gen_0;
    2628          49 :   y = primitive_part(y, &c);
    2629          49 :   return c? gerepileupto(av, y): gerepilecopy(av, y);
    2630             : }
    2631             : static GEN
    2632           0 : RgV_deplin(GEN v)
    2633             : {
    2634           0 :   pari_sp av = avma;
    2635           0 :   long n = lg(v)-1;
    2636           0 :   GEN y, p = NULL;
    2637           0 :   if (n <= 1)
    2638             :   {
    2639           0 :     if (n == 1 && gequal0(gel(v,1))) return mkcol(gen_1);
    2640           0 :     return cgetg(1, t_COL);
    2641             :   }
    2642           0 :   if (gequal0(gel(v,1))) return scalarcol_shallow(gen_1, n);
    2643           0 :   v = primpart(mkvec2(gel(v,1),gel(v,2)));
    2644           0 :   if (RgV_is_FpV(v, &p) && p) v = centerlift(v);
    2645           0 :   y = zerocol(n);
    2646           0 :   gel(y,1) = gneg(gel(v,2));
    2647           0 :   gel(y,2) = gcopy(gel(v,1));
    2648           0 :   return gerepileupto(av, y);
    2649             : 
    2650             : }
    2651             : GEN
    2652         161 : deplin(GEN x)
    2653             : {
    2654         161 :   GEN p = NULL;
    2655         161 :   switch(typ(x))
    2656             :   {
    2657         161 :     case t_MAT: break;
    2658           0 :     case t_VEC: return RgV_deplin(x);
    2659           0 :     default: pari_err_TYPE("deplin",x);
    2660             :   }
    2661         161 :   if (RgM_is_FpM(x, &p) && p)
    2662             :   {
    2663          63 :     pari_sp av = avma;
    2664             :     ulong pp;
    2665          63 :     x = RgM_Fp_init(x, p, &pp);
    2666          63 :     switch(pp)
    2667             :     {
    2668             :     case 0:
    2669          14 :       x = FpM_ker_gen(x,p,1);
    2670          14 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2671           7 :       x = FpC_center(x,p,shifti(p,-1));
    2672           7 :       break;
    2673             :     case 2:
    2674          14 :       x = F2m_ker_sp(x,1);
    2675          14 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2676           7 :       x = F2c_to_ZC(x); break;
    2677             :     default:
    2678          35 :       x = Flm_ker_sp(x,pp,1);
    2679          35 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2680          21 :       x = Flv_center(x, pp, pp>>1);
    2681          21 :       x = zc_to_ZC(x);
    2682          21 :       break;
    2683             :     }
    2684          35 :     return gerepileupto(av, x);
    2685             :   }
    2686          98 :   return deplin_aux(x);
    2687             : }
    2688             : /*******************************************************************/
    2689             : /*                                                                 */
    2690             : /*         GAUSS REDUCTION OF MATRICES  (m lines x n cols)         */
    2691             : /*           (kernel, image, complementary image, rank)            */
    2692             : /*                                                                 */
    2693             : /*******************************************************************/
    2694             : /* return the transform of x under a standard Gauss pivot.
    2695             :  * x0 is a reference point when guessing whether x[i,j] ~ 0
    2696             :  * (iff x[i,j] << x0[i,j])
    2697             :  * Set r = dim ker(x). d[k] contains the index of the first non-zero pivot
    2698             :  * in column k */
    2699             : static GEN
    2700        6433 : gauss_pivot_ker(GEN x, GEN x0, GEN *dd, long *rr)
    2701             : {
    2702             :   GEN c, d, p, data;
    2703             :   pari_sp av;
    2704             :   long i, j, k, r, t, n, m;
    2705             :   pivot_fun pivot;
    2706             : 
    2707        6433 :   n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
    2708        6398 :   m=nbrows(x); r=0;
    2709        6398 :   pivot = get_pivot_fun(x, x0, &data);
    2710        6398 :   x = RgM_shallowcopy(x);
    2711        6398 :   c = zero_zv(m);
    2712        6398 :   d = cgetg(n+1,t_VECSMALL);
    2713        6398 :   av=avma;
    2714       37758 :   for (k=1; k<=n; k++)
    2715             :   {
    2716       31360 :     j = pivot(x, data, k, c);
    2717       31360 :     if (j > m)
    2718             :     {
    2719       10234 :       r++; d[k]=0;
    2720       55594 :       for(j=1; j<k; j++)
    2721       45360 :         if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
    2722             :     }
    2723             :     else
    2724             :     { /* pivot for column k on row j */
    2725       21126 :       c[j]=k; d[k]=j; p = gdiv(gen_m1,gcoeff(x,j,k));
    2726       21126 :       gcoeff(x,j,k) = gen_m1;
    2727             :       /* x[j,] /= - x[j,k] */
    2728       21126 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
    2729      274456 :       for (t=1; t<=m; t++)
    2730      253330 :         if (t!=j)
    2731             :         { /* x[t,] -= 1 / x[j,k] x[j,] */
    2732      232204 :           p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
    2733     2968077 :           for (i=k+1; i<=n; i++)
    2734     2735873 :             gcoeff(x,t,i) = gadd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
    2735      232204 :           if (gc_needed(av,1)) gerepile_gauss_ker(x,k,t,av);
    2736             :         }
    2737             :     }
    2738             :   }
    2739        6398 :   *dd=d; *rr=r; return x;
    2740             : }
    2741             : 
    2742             : static GEN FpM_gauss_pivot(GEN x, GEN p, long *rr);
    2743             : static GEN FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr);
    2744             : static GEN F2m_gauss_pivot(GEN x, long *rr);
    2745             : static GEN Flm_gauss_pivot(GEN x, ulong p, long *rry);
    2746             : 
    2747             : /* r = dim ker(x).
    2748             :  * Returns d:
    2749             :  *   d[k] != 0 contains the index of a non-zero pivot in column k
    2750             :  *   d[k] == 0 if column k is a linear combination of the (k-1) first ones */
    2751             : GEN
    2752        7635 : RgM_pivots(GEN x0, GEN data, long *rr, pivot_fun pivot)
    2753             : {
    2754             :   GEN x, c, d, p;
    2755        7635 :   long i, j, k, r, t, m, n = lg(x0)-1;
    2756             :   pari_sp av;
    2757             : 
    2758        7635 :   if (RgM_is_ZM(x0)) return ZM_pivots(x0, rr);
    2759        7159 :   if (!n) { *rr = 0; return NULL; }
    2760             : 
    2761        7159 :   d = cgetg(n+1, t_VECSMALL);
    2762        7159 :   x = RgM_shallowcopy(x0);
    2763        7159 :   m = nbrows(x); r = 0;
    2764        7159 :   c = zero_zv(m);
    2765        7159 :   av = avma;
    2766      815787 :   for (k=1; k<=n; k++)
    2767             :   {
    2768      808628 :     j = pivot(x, data, k, c);
    2769      808628 :     if (j > m) { r++; d[k] = 0; }
    2770             :     else
    2771             :     {
    2772       19335 :       c[j] = k; d[k] = j; p = gdiv(gen_m1, gcoeff(x,j,k));
    2773       19335 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
    2774             : 
    2775      105151 :       for (t=1; t<=m; t++)
    2776       85816 :         if (!c[t]) /* no pivot on that line yet */
    2777             :         {
    2778       39624 :           p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
    2779     5104693 :           for (i=k+1; i<=n; i++)
    2780     5065069 :             gcoeff(x,t,i) = gadd(gcoeff(x,t,i), gmul(p, gcoeff(x,j,i)));
    2781       39624 :           if (gc_needed(av,1)) gerepile_gauss(x,k,t,av,j,c);
    2782             :         }
    2783       19335 :       for (i=k; i<=n; i++) gcoeff(x,j,i) = gen_0; /* dummy */
    2784             :     }
    2785             :   }
    2786        7159 :   *rr = r; avma = (pari_sp)d; return d;
    2787             : }
    2788             : 
    2789             : static long
    2790       96927 : ZM_count_0_cols(GEN M)
    2791             : {
    2792       96927 :   long i, l = lg(M), n = 0;
    2793      483115 :   for (i = 1; i < l; i++)
    2794      386188 :     if (ZV_equal0(gel(M,i))) n++;
    2795       96927 :   return n;
    2796             : }
    2797             : 
    2798             : static void indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol);
    2799             : /* As RgM_pivots, integer entries. Set *rr = dim Ker M0 */
    2800             : GEN
    2801      100112 : ZM_pivots(GEN M0, long *rr)
    2802             : {
    2803      100112 :   GEN d, dbest = NULL;
    2804             :   long m, n, i, imax, rmin, rbest, zc;
    2805      100112 :   int beenthere = 0;
    2806      100112 :   pari_sp av, av0 = avma;
    2807             :   forprime_t S;
    2808             : 
    2809      100112 :   rbest = n = lg(M0)-1;
    2810      100112 :   if (n == 0) { *rr = 0; return NULL; }
    2811       96927 :   zc = ZM_count_0_cols(M0);
    2812       96927 :   if (n == zc) { *rr = zc; return zero_zv(n); }
    2813             : 
    2814       96843 :   m = nbrows(M0);
    2815       96843 :   rmin = maxss(zc, n-m);
    2816       96843 :   init_modular(&S);
    2817       96843 :   imax = (n < (1<<4))? 1: (n>>3); /* heuristic */
    2818             : 
    2819             :   for(;;)
    2820             :   {
    2821             :     GEN row, col, M, KM, IM, RHS, X, cX;
    2822             :     long rk;
    2823      100053 :     for (av = avma, i = 0;; avma = av, i++)
    2824             :     {
    2825      100053 :       ulong p = u_forprime_next(&S);
    2826             :       long rp;
    2827      100053 :       if (!p) pari_err_OVERFLOW("ZM_pivots [ran out of primes]");
    2828      100053 :       d = Flm_gauss_pivot(ZM_to_Flm(M0, p), p, &rp);
    2829      100053 :       if (rp == rmin) { rbest = rp; goto END; } /* maximal rank, return */
    2830        5307 :       if (rp < rbest) { /* save best r so far */
    2831        2132 :         rbest = rp;
    2832        2132 :         if (dbest) gunclone(dbest);
    2833        2132 :         dbest = gclone(d);
    2834        4229 :         if (beenthere) break;
    2835             :       }
    2836        5307 :       if (!beenthere && i >= imax) break;
    2837        3210 :     }
    2838        2097 :     beenthere = 1;
    2839             :     /* Dubious case: there is (probably) a non trivial kernel */
    2840        2097 :     indexrank_all(m,n, rbest, dbest, &row, &col);
    2841        2097 :     M = rowpermute(vecpermute(M0, col), row);
    2842        2097 :     rk = n - rbest; /* (probable) dimension of image */
    2843        2097 :     IM = vecslice(M,1,rk);
    2844        2097 :     KM = vecslice(M,rk+1, n);
    2845        2097 :     M = rowslice(IM, 1,rk); /* square maximal rank */
    2846        2097 :     X = ZM_gauss(M, rowslice(KM, 1,rk));
    2847        2097 :     X = Q_remove_denom(X, &cX);
    2848        2097 :     RHS = rowslice(KM,rk+1,m);
    2849        2097 :     if (cX) RHS = ZM_Z_mul(RHS, cX);
    2850        2097 :     if (ZM_equal(ZM_mul(rowslice(IM,rk+1,m), X), RHS))
    2851             :     {
    2852        2097 :       d = vecsmall_copy(dbest);
    2853        2097 :       goto END;
    2854             :     }
    2855           0 :     avma = av;
    2856           0 :   }
    2857             : END:
    2858       96843 :   *rr = rbest; if (dbest) gunclone(dbest);
    2859       96843 :   return gerepileuptoleaf(av0, d);
    2860             : }
    2861             : 
    2862             : /* set *pr = dim Ker x */
    2863             : static GEN
    2864         546 : gauss_pivot(GEN x, long *pr) {
    2865             :   GEN data;
    2866         546 :   pivot_fun pivot = get_pivot_fun(x, x, &data);
    2867         546 :   return RgM_pivots(x, data, pr, pivot);
    2868             : }
    2869             : 
    2870             : /* compute ker(x), x0 is a reference point when guessing whether x[i,j] ~ 0
    2871             :  * (iff x[i,j] << x0[i,j]) */
    2872             : static GEN
    2873        6433 : ker_aux(GEN x, GEN x0)
    2874             : {
    2875        6433 :   pari_sp av = avma;
    2876             :   GEN d,y;
    2877             :   long i,j,k,r,n;
    2878             : 
    2879        6433 :   x = gauss_pivot_ker(x,x0,&d,&r);
    2880        6433 :   if (!r) { avma=av; return cgetg(1,t_MAT); }
    2881        6216 :   n = lg(x)-1; y=cgetg(r+1,t_MAT);
    2882       16450 :   for (j=k=1; j<=r; j++,k++)
    2883             :   {
    2884       10234 :     GEN p = cgetg(n+1,t_COL);
    2885             : 
    2886       10234 :     gel(y,j) = p; while (d[k]) k++;
    2887       55594 :     for (i=1; i<k; i++)
    2888       45360 :       if (d[i])
    2889             :       {
    2890       24668 :         GEN p1=gcoeff(x,d[i],k);
    2891       24668 :         gel(p,i) = gcopy(p1); gunclone(p1);
    2892             :       }
    2893             :       else
    2894       20692 :         gel(p,i) = gen_0;
    2895       10234 :     gel(p,k) = gen_1; for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
    2896             :   }
    2897        6216 :   return gerepileupto(av,y);
    2898             : }
    2899             : GEN
    2900        6335 : ker(GEN x)
    2901             : {
    2902        6335 :   pari_sp av = avma;
    2903        6335 :   GEN p = NULL, ff = NULL;
    2904        6335 :   if (RgM_is_FpM(x, &p) && p)
    2905             :   {
    2906             :     ulong pp;
    2907          35 :     x = RgM_Fp_init(x, p, &pp);
    2908          35 :     switch(pp)
    2909             :     {
    2910          14 :     case 0: x = FpM_to_mod(FpM_ker_gen(x,p,0),p); break;
    2911           7 :     case 2: x = F2m_to_mod(F2m_ker_sp(x,0)); break;
    2912          14 :     default:x = Flm_to_mod(Flm_ker_sp(x,pp,0), pp); break;
    2913             :     }
    2914          35 :     return gerepileupto(av, x);
    2915             :   }
    2916        6300 :   if (RgM_is_FFM(x, &ff)) return FFM_ker(x, ff);
    2917        6279 :   return ker_aux(x,x);
    2918             : }
    2919             : GEN
    2920         133 : matker0(GEN x,long flag)
    2921             : {
    2922         133 :   if (typ(x)!=t_MAT) pari_err_TYPE("matker",x);
    2923         133 :   if (!flag) return ker(x);
    2924           7 :   RgM_check_ZM(x, "keri");
    2925           7 :   return keri(x);
    2926             : }
    2927             : 
    2928             : GEN
    2929         378 : image(GEN x)
    2930             : {
    2931         378 :   pari_sp av = avma;
    2932         378 :   GEN d, ff = NULL, p = NULL;
    2933             :   long r;
    2934             : 
    2935         378 :   if (typ(x)!=t_MAT) pari_err_TYPE("matimage",x);
    2936         378 :   if (RgM_is_FpM(x, &p) && p)
    2937             :   {
    2938             :     ulong pp;
    2939          35 :     x = RgM_Fp_init(x, p, &pp);
    2940          35 :     switch(pp)
    2941             :     {
    2942          14 :     case 0: x = FpM_to_mod(FpM_image(x,p), p); break;
    2943           7 :     case 2: x = F2m_to_mod(F2m_image(x)); break;
    2944          14 :     default:x = Flm_to_mod(Flm_image(x,pp), pp);
    2945             :     }
    2946          35 :     return gerepileupto(av, x);
    2947             :   }
    2948         343 :   if (RgM_is_FFM(x, &ff)) return FFM_image(x, ff);
    2949         322 :   d = gauss_pivot(x,&r); /* d left on stack for efficiency */
    2950         322 :   return image_from_pivot(x,d,r);
    2951             : }
    2952             : 
    2953             : static GEN
    2954          84 : imagecompl_aux(GEN x, GEN(*PIVOT)(GEN,long*))
    2955             : {
    2956          84 :   pari_sp av = avma;
    2957             :   GEN d,y;
    2958             :   long j,i,r;
    2959             : 
    2960          84 :   if (typ(x)!=t_MAT) pari_err_TYPE("imagecompl",x);
    2961          84 :   (void)new_chunk(lg(x) * 4 + 1); /* HACK */
    2962          84 :   d = PIVOT(x,&r); /* if (!d) then r = 0 */
    2963          84 :   avma = av; y = cgetg(r+1,t_VECSMALL);
    2964         126 :   for (i=j=1; j<=r; i++)
    2965          42 :     if (!d[i]) y[j++] = i;
    2966          84 :   return y;
    2967             : }
    2968             : GEN
    2969          84 : imagecompl(GEN x) { return imagecompl_aux(x, &gauss_pivot); }
    2970             : GEN
    2971           0 : ZM_imagecompl(GEN x) { return imagecompl_aux(x, &ZM_pivots); }
    2972             : 
    2973             : /* permutation giving imagecompl(x') | image(x'), x' = transpose of x */
    2974             : static GEN
    2975       59214 : imagecomplspec_aux(GEN x, long *nlze, GEN(*PIVOT)(GEN,long*))
    2976             : {
    2977       59214 :   pari_sp av = avma;
    2978             :   GEN d,y;
    2979             :   long i,j,k,l,r;
    2980             : 
    2981       59214 :   if (typ(x)!=t_MAT) pari_err_TYPE("imagecomplspec",x);
    2982       59214 :   x = shallowtrans(x); l = lg(x);
    2983       59214 :   d = PIVOT(x,&r);
    2984       59214 :   *nlze = r;
    2985       59214 :   avma = av; /* HACK: shallowtrans(x) big enough to avoid overwriting d */
    2986       59214 :   if (!d) return identity_perm(l-1);
    2987       56848 :   y = cgetg(l,t_VECSMALL);
    2988      307555 :   for (i=j=1, k=r+1; i<l; i++)
    2989      250707 :     if (d[i]) y[k++]=i; else y[j++]=i;
    2990       56848 :   return y;
    2991             : }
    2992             : GEN
    2993           0 : imagecomplspec(GEN x, long *nlze)
    2994           0 : { return imagecomplspec_aux(x,nlze,&gauss_pivot); }
    2995             : GEN
    2996       59214 : ZM_imagecomplspec(GEN x, long *nlze)
    2997       59214 : { return imagecomplspec_aux(x,nlze,&ZM_pivots); }
    2998             : 
    2999             : GEN
    3000        1666 : RgM_RgC_invimage(GEN A, GEN y)
    3001             : {
    3002        1666 :   pari_sp av = avma;
    3003        1666 :   long i, l = lg(A);
    3004        1666 :   GEN M, x, t, p = NULL;
    3005             : 
    3006        1666 :   if (RgM_is_FpM(A, &p) && RgV_is_FpV(y, &p) && p)
    3007             :   {
    3008             :     ulong pp;
    3009          28 :     A = RgM_Fp_init(A,p,&pp);
    3010          28 :     switch(pp)
    3011             :     {
    3012             :     case 0:
    3013           7 :       y = RgC_to_FpC(y,p);
    3014           7 :       x = FpM_FpC_invimage(A, y, p);
    3015           7 :       if (x) x = FpC_to_mod(x,p);
    3016           7 :       break;
    3017             :     case 2:
    3018           7 :       y = RgV_to_F2v(y);
    3019           7 :       x = F2m_F2c_invimage(A, y);
    3020           7 :       if (x) x = F2c_to_mod(x);
    3021           7 :       break;
    3022             :     default:
    3023          14 :       y = RgV_to_Flv(y,pp);
    3024          14 :       x = Flm_Flc_invimage(A, y, pp);
    3025          14 :       if (x) x = Flc_to_mod(x,pp);
    3026             :     }
    3027          28 :     if (!x) { avma = av; return NULL; }
    3028          28 :     return gerepileupto(av, x);
    3029             :   }
    3030             : 
    3031        1638 :   if (l==1) return NULL;
    3032        1638 :   if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
    3033        1638 :   M = ker(shallowconcat(A, y));
    3034        1638 :   i = lg(M)-1;
    3035        1638 :   if (!i) { avma = av; return NULL; }
    3036             : 
    3037        1638 :   x = gel(M,i); t = gel(x,l);
    3038        1638 :   if (gequal0(t)) { avma = av; return NULL; }
    3039             : 
    3040        1610 :   t = gneg_i(t); setlg(x,l);
    3041        1610 :   return gerepileupto(av, RgC_Rg_div(x, t));
    3042             : }
    3043             : GEN
    3044       31852 : FpM_FpC_invimage(GEN A, GEN y, GEN p)
    3045             : {
    3046       31852 :   pari_sp av = avma;
    3047       31852 :   long i, l = lg(A);
    3048             :   GEN M, x, t;
    3049             : 
    3050       31852 :   if (lgefint(p) == 3)
    3051             :   {
    3052       31845 :     ulong pp = p[2];
    3053       31845 :     A = ZM_to_Flm(A, pp);
    3054       31845 :     y = ZV_to_Flv(y, pp);
    3055       31845 :     x = Flm_Flc_invimage(A,y,pp);
    3056       31845 :     if (!x) { avma = av; return NULL; }
    3057       31845 :     return gerepileupto(av, Flc_to_ZC(x));
    3058             :   }
    3059           7 :   if (l==1) return NULL;
    3060           7 :   if (lg(y) != lgcols(A)) pari_err_DIM("FpM_FpC_invimage");
    3061           7 :   M = FpM_ker(shallowconcat(A,y),p);
    3062           7 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3063             : 
    3064           7 :   x = gel(M,i); t = gel(x,l);
    3065           7 :   if (!signe(t)) { avma = av; return NULL; }
    3066             : 
    3067           7 :   setlg(x,l); t = Fp_inv(negi(t),p);
    3068           7 :   if (is_pm1(t)) return gerepilecopy(av, x);
    3069           7 :   return gerepileupto(av, FpC_Fp_mul(x, t, p));
    3070             : }
    3071             : GEN
    3072       36430 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
    3073             : {
    3074       36430 :   pari_sp av = avma;
    3075       36430 :   long i, l = lg(A);
    3076             :   GEN M, x;
    3077             :   ulong t;
    3078             : 
    3079       36430 :   if (l==1) return NULL;
    3080       36430 :   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
    3081       36430 :   M = cgetg(l+1,t_MAT);
    3082       36430 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    3083       36430 :   gel(M,l) = y; M = Flm_ker(M,p);
    3084       36430 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3085             : 
    3086       36430 :   x = gel(M,i); t = x[l];
    3087       36430 :   if (!t) { avma = av; return NULL; }
    3088             : 
    3089       36430 :   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
    3090       36430 :   if (t!=1) x = Flv_Fl_mul(x, t, p);
    3091       36430 :   return gerepileuptoleaf(av, x);
    3092             : }
    3093             : GEN
    3094          21 : F2m_F2c_invimage(GEN A, GEN y)
    3095             : {
    3096          21 :   pari_sp av = avma;
    3097          21 :   long i, l = lg(A);
    3098             :   GEN M, x;
    3099             : 
    3100          21 :   if (l==1) return NULL;
    3101          21 :   if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
    3102          21 :   M = cgetg(l+1,t_MAT);
    3103          21 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    3104          21 :   gel(M,l) = y; M = F2m_ker(M);
    3105          21 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3106             : 
    3107          21 :   x = gel(M,i);
    3108          21 :   if (!F2v_coeff(x,l)) { avma = av; return NULL; }
    3109          21 :   F2v_clear(x, x[1]); x[1]--; /* remove last coord */
    3110          21 :   return gerepileuptoleaf(av, x);
    3111             : }
    3112             : 
    3113             : /* Return X such that m X = v (t_COL or t_MAT), resp. an empty t_COL / t_MAT
    3114             :  * if no solution exist */
    3115             : GEN
    3116        1428 : inverseimage(GEN m, GEN v)
    3117             : {
    3118             :   GEN y;
    3119        1428 :   if (typ(m)!=t_MAT) pari_err_TYPE("inverseimage",m);
    3120        1428 :   switch(typ(v))
    3121             :   {
    3122             :     case t_COL:
    3123        1393 :       y = RgM_RgC_invimage(m,v);
    3124        1393 :       return y? y: cgetg(1,t_COL);
    3125             :     case t_MAT:
    3126          35 :       y = RgM_invimage(m, v);
    3127          35 :       return y? y: cgetg(1,t_MAT);
    3128             :   }
    3129           0 :   pari_err_TYPE("inverseimage",v);
    3130           0 :   return NULL;/*not reached*/
    3131             : }
    3132             : 
    3133             : static GEN
    3134          14 : Flm_invimage_i(GEN A, GEN B, ulong p)
    3135             : {
    3136             :   GEN d, x, X, Y;
    3137          14 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3138          14 :   x = Flm_ker_sp(shallowconcat(Flm_neg(A,p), B), p, 0);
    3139             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3140             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3141             :    * Y has at least nB columns and full rank */
    3142          14 :   nY = lg(x)-1;
    3143          14 :   if (nY < nB) return NULL;
    3144          14 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3145          14 :   d = cgetg(nB+1, t_VECSMALL);
    3146          42 :   for (i = nB, j = nY; i >= 1; i--)
    3147             :   {
    3148          42 :     for (; j>=1; j--)
    3149          42 :       if (coeff(Y,i,j)) { d[i] = j; break; }
    3150          28 :     if (!j) return NULL;
    3151             :   }
    3152             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3153          14 :   Y = vecpermute(Y, d);
    3154          14 :   x = vecpermute(x, d);
    3155          14 :   X = rowslice(x, 1, nA);
    3156          14 :   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
    3157             : }
    3158             : 
    3159             : static GEN
    3160           7 : F2m_invimage_i(GEN A, GEN B)
    3161             : {
    3162             :   GEN d, x, X, Y;
    3163           7 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3164           7 :   x = F2m_ker_sp(shallowconcat(A, B), 0);
    3165             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3166             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3167             :    * Y has at least nB columns and full rank */
    3168           7 :   nY = lg(x)-1;
    3169           7 :   if (nY < nB) return NULL;
    3170             : 
    3171             :   /* implicitly: Y = rowslice(x, nA+1, nA+nB), nB rows */
    3172           7 :   d = cgetg(nB+1, t_VECSMALL);
    3173          21 :   for (i = nB, j = nY; i >= 1; i--)
    3174             :   {
    3175          21 :     for (; j>=1; j--)
    3176          21 :       if (F2m_coeff(x,nA+i,j)) { d[i] = j; break; } /* Y[i,j] */
    3177          14 :     if (!j) return NULL;
    3178             :   }
    3179           7 :   x = vecpermute(x, d);
    3180             : 
    3181           7 :   X = F2m_rowslice(x, 1, nA);
    3182           7 :   Y = F2m_rowslice(x, nA+1, nA+nB);
    3183           7 :   return F2m_mul(X, F2m_inv_upper_1(Y));
    3184             : }
    3185             : GEN
    3186           0 : Flm_invimage(GEN A, GEN B, ulong p)
    3187             : {
    3188           0 :   pari_sp av = avma;
    3189           0 :   GEN X = Flm_invimage_i(A,B,p);
    3190           0 :   if (!X) { avma = av; return NULL; }
    3191           0 :   return gerepileupto(av, X);
    3192             : }
    3193             : GEN
    3194           0 : F2m_invimage(GEN A, GEN B)
    3195             : {
    3196           0 :   pari_sp av = avma;
    3197           0 :   GEN X = F2m_invimage_i(A,B);
    3198           0 :   if (!X) { avma = av; return NULL; }
    3199           0 :   return gerepileupto(av, X);
    3200             : }
    3201             : static GEN
    3202           7 : FpM_invimage_i(GEN A, GEN B, GEN p)
    3203             : {
    3204             :   GEN d, x, X, Y;
    3205           7 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3206           7 :   if (lgefint(p) == 3)
    3207             :   {
    3208           0 :     ulong pp = p[2];
    3209           0 :     A = ZM_to_Flm(A, pp);
    3210           0 :     B = ZM_to_Flm(B, pp);
    3211           0 :     x = Flm_invimage_i(A, B, pp);
    3212           0 :     return x? Flm_to_ZM(x): NULL;
    3213             :   }
    3214           7 :   x = FpM_ker(shallowconcat(ZM_neg(A), B), p);
    3215             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3216             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3217             :    * Y has at least nB columns and full rank */
    3218           7 :   nY = lg(x)-1;
    3219           7 :   if (nY < nB) return NULL;
    3220           7 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3221           7 :   d = cgetg(nB+1, t_VECSMALL);
    3222          21 :   for (i = nB, j = nY; i >= 1; i--)
    3223             :   {
    3224          21 :     for (; j>=1; j--)
    3225          21 :       if (signe(gcoeff(Y,i,j))) { d[i] = j; break; }
    3226          14 :     if (!j) return NULL;
    3227             :   }
    3228             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3229           7 :   Y = vecpermute(Y, d);
    3230           7 :   x = vecpermute(x, d);
    3231           7 :   X = rowslice(x, 1, nA);
    3232           7 :   return FpM_mul(X, FpM_inv_upper_1(Y,p), p);
    3233             : }
    3234             : GEN
    3235           0 : FpM_invimage(GEN A, GEN B, GEN p)
    3236             : {
    3237           0 :   pari_sp av = avma;
    3238           0 :   GEN X = FpM_invimage_i(A,B,p);
    3239           0 :   if (!X) { avma = av; return NULL; }
    3240           0 :   return gerepileupto(av, X);
    3241             : }
    3242             : 
    3243             : /* find Z such that A Z = B. Return NULL if no solution */
    3244             : GEN
    3245          35 : RgM_invimage(GEN A, GEN B)
    3246             : {
    3247          35 :   pari_sp av = avma;
    3248             :   GEN d, x, X, Y;
    3249          35 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3250          35 :   GEN p = NULL;
    3251          35 :   if (RgM_is_FpM(A, &p) && RgM_is_FpM(B, &p) && p)
    3252             :   {
    3253             :     ulong pp;
    3254          28 :     A = RgM_Fp_init(A,p,&pp);
    3255          28 :     switch(pp)
    3256             :     {
    3257             :     case 0:
    3258           7 :       B = RgM_to_FpM(B,p);
    3259           7 :       x = FpM_invimage_i(A, B, p);
    3260           7 :       if (x) x = FpM_to_mod(x, p);
    3261           7 :     break;
    3262             :     case 2:
    3263           7 :       B = RgM_to_F2m(B);
    3264           7 :       x = F2m_invimage_i(A, B);
    3265           7 :       if (x) x = F2m_to_mod(x);
    3266           7 :       break;
    3267             :     default:
    3268          14 :       B = RgM_to_Flm(B,pp);
    3269          14 :       x = Flm_invimage_i(A, B, pp);
    3270          14 :       if (x) x = Flm_to_mod(x,pp);
    3271          14 :       break;
    3272             :     }
    3273          28 :     if (!x) { avma = av; return NULL; }
    3274          28 :     return gerepileupto(av, x);
    3275             :   }
    3276           7 :   x = ker(shallowconcat(RgM_neg(A), B));
    3277             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3278             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3279             :    * Y has at least nB columns and full rank */
    3280           7 :   nY = lg(x)-1;
    3281           7 :   if (nY < nB) { avma = av; return NULL; }
    3282           7 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3283           7 :   d = cgetg(nB+1, t_VECSMALL);
    3284          21 :   for (i = nB, j = nY; i >= 1; i--)
    3285             :   {
    3286          21 :     for (; j>=1; j--)
    3287          21 :       if (!gequal0(gcoeff(Y,i,j))) { d[i] = j; break; }
    3288          14 :     if (!j) { avma = av; return NULL; }
    3289             :   }
    3290             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3291           7 :   Y = vecpermute(Y, d);
    3292           7 :   x = vecpermute(x, d);
    3293           7 :   X = rowslice(x, 1, nA);
    3294           7 :   return gerepileupto(av, RgM_mul(X, RgM_inv_upper(Y)));
    3295             : }
    3296             : 
    3297             : /* r = dim Ker x, n = nbrows(x) */
    3298             : static GEN
    3299       24061 : get_suppl(GEN x, GEN d, long n, long r, GEN(*ei)(long,long))
    3300             : {
    3301             :   pari_sp av;
    3302             :   GEN y, c;
    3303       24061 :   long j, k, rx = lg(x)-1; /* != 0 due to init_suppl() */
    3304             : 
    3305       24061 :   if (rx == n && r == 0) return gcopy(x);
    3306       21277 :   y = cgetg(n+1, t_MAT);
    3307       21277 :   av = avma; c = zero_zv(n);
    3308             :   /* c = lines containing pivots (could get it from gauss_pivot, but cheap)
    3309             :    * In theory r = 0 and d[j] > 0 for all j, but why take chances? */
    3310      173598 :   for (k = j = 1; j<=rx; j++)
    3311      152321 :     if (d[j]) { c[ d[j] ] = 1; gel(y,k++) = gel(x,j); }
    3312      230022 :   for (j=1; j<=n; j++)
    3313      208745 :     if (!c[j]) gel(y,k++) = (GEN)j; /* HACK */
    3314       21277 :   avma = av;
    3315             : 
    3316       21277 :   rx -= r;
    3317       21277 :   for (j=1; j<=rx; j++) gel(y,j) = gcopy(gel(y,j));
    3318       21277 :   for (   ; j<=n; j++)  gel(y,j) = ei(n, y[j]);
    3319       21277 :   return y;
    3320             : }
    3321             : 
    3322             : static void
    3323       24061 : init_suppl(GEN x)
    3324             : {
    3325       24061 :   if (lg(x) == 1) pari_err_IMPL("suppl [empty matrix]");
    3326             :   /* HACK: avoid overwriting d from gauss_pivot() after avma=av */
    3327       24061 :   (void)new_chunk(lgcols(x) * 2);
    3328       24061 : }
    3329             : 
    3330             : /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
    3331             :  * whose first k columns are given by x. If rank(x) < k, undefined result. */
    3332             : GEN
    3333          70 : suppl(GEN x)
    3334             : {
    3335          70 :   pari_sp av = avma;
    3336          70 :   GEN d, X = x, p = NULL;
    3337             :   long r;
    3338             : 
    3339          70 :   if (typ(x)!=t_MAT) pari_err_TYPE("suppl",x);
    3340          70 :   if (RgM_is_FpM(x, &p) && p)
    3341             :   {
    3342             :     ulong pp;
    3343          28 :     x = RgM_Fp_init(x, p, &pp);
    3344          28 :     switch(pp)
    3345             :     {
    3346           7 :     case 0: x = FpM_to_mod(FpM_suppl(x,p), p); break;
    3347           7 :     case 2: x = F2m_to_mod(F2m_suppl(x)); break;
    3348          14 :     default:x = Flm_to_mod(Flm_suppl(x,pp), pp); break;
    3349             :     }
    3350          28 :     return gerepileupto(av, x);
    3351             :   }
    3352          42 :   avma = av; init_suppl(x);
    3353          42 :   d = gauss_pivot(X,&r);
    3354          42 :   avma = av; return get_suppl(X,d,nbrows(X),r,&col_ei);
    3355             : }
    3356             : GEN
    3357       23879 : FpM_suppl(GEN x, GEN p)
    3358             : {
    3359       23879 :   pari_sp av = avma;
    3360             :   GEN d;
    3361             :   long r;
    3362       23879 :   init_suppl(x); d = FpM_gauss_pivot(x,p, &r);
    3363       23879 :   avma = av; return get_suppl(x,d,nbrows(x),r,&col_ei);
    3364             : }
    3365             : GEN
    3366          14 : Flm_suppl(GEN x, ulong p)
    3367             : {
    3368          14 :   pari_sp av = avma;
    3369             :   GEN d;
    3370             :   long r;
    3371          14 :   init_suppl(x); d = Flm_gauss_pivot(Flm_copy(x),p, &r);
    3372          14 :   avma = av; return get_suppl(x,d,nbrows(x),r,&vecsmall_ei);
    3373             : 
    3374             : }
    3375             : GEN
    3376           7 : F2m_suppl(GEN x)
    3377             : {
    3378           7 :   pari_sp av = avma;
    3379             :   GEN d;
    3380             :   long r;
    3381           7 :   init_suppl(x); d = F2m_gauss_pivot(F2m_copy(x), &r);
    3382           7 :   avma = av; return get_suppl(x,d,mael(x,1,1),r,&F2v_ei);
    3383             : }
    3384             : 
    3385             : GEN
    3386         581 : FqM_suppl(GEN x, GEN T, GEN p)
    3387             : {
    3388         581 :   pari_sp av = avma;
    3389             :   GEN d;
    3390             :   long r;
    3391             : 
    3392         581 :   if (!T) return FpM_suppl(x,p);
    3393         119 :   init_suppl(x);
    3394         119 :   d = FqM_gauss_pivot(x,T,p,&r);
    3395         119 :   avma = av; return get_suppl(x,d,nbrows(x),r,&col_ei);
    3396             : }
    3397             : 
    3398             : GEN
    3399           7 : image2(GEN x)
    3400             : {
    3401           7 :   pari_sp av = avma;
    3402             :   long k, n, i;
    3403             :   GEN A, B;
    3404             : 
    3405           7 :   if (typ(x)!=t_MAT) pari_err_TYPE("image2",x);
    3406           7 :   if (lg(x) == 1) return cgetg(1,t_MAT);
    3407           7 :   A = ker(x); k = lg(A)-1;
    3408           7 :   if (!k) { avma = av; return gcopy(x); }
    3409           7 :   A = suppl(A); n = lg(A)-1;
    3410           7 :   B = cgetg(n-k+1, t_MAT);
    3411           7 :   for (i = k+1; i <= n; i++) gel(B,i-k) = RgM_RgC_mul(x, gel(A,i));
    3412           7 :   return gerepileupto(av, B);
    3413             : }
    3414             : 
    3415             : GEN
    3416         119 : matimage0(GEN x,long flag)
    3417             : {
    3418         119 :   switch(flag)
    3419             :   {
    3420         112 :     case 0: return image(x);
    3421           7 :     case 1: return image2(x);
    3422           0 :     default: pari_err_FLAG("matimage");
    3423             :   }
    3424           0 :   return NULL; /* not reached */
    3425             : }
    3426             : 
    3427             : long
    3428         154 : rank(GEN x)
    3429             : {
    3430         154 :   pari_sp av = avma;
    3431             :   long r;
    3432         154 :   GEN ff = NULL, p = NULL;
    3433             : 
    3434         154 :   if (typ(x)!=t_MAT) pari_err_TYPE("rank",x);
    3435         154 :   if (RgM_is_FpM(x, &p) && p)
    3436             :   {
    3437             :     ulong pp;
    3438          84 :     x = RgM_Fp_init(x,p,&pp);
    3439          84 :     switch(pp)
    3440             :     {
    3441           7 :     case 0: r = FpM_rank(x,p); break;
    3442          63 :     case 2: r = F2m_rank(x); break;
    3443          14 :     default:r = Flm_rank(x,pp); break;
    3444             :     }
    3445          84 :     avma = av; return r;
    3446             :   }
    3447          70 :   if (RgM_is_FFM(x, &ff)) return FFM_rank(x, ff);
    3448          49 :   (void)gauss_pivot(x, &r);
    3449          49 :   avma = av; return lg(x)-1 - r;
    3450             : }
    3451             : 
    3452             : /* d a t_VECSMALL of integers in 1..n. Return the vector of the d[i]
    3453             :  * followed by the missing indices */
    3454             : static GEN
    3455        4194 : perm_complete(GEN d, long n)
    3456             : {
    3457        4194 :   GEN y = cgetg(n+1, t_VECSMALL);
    3458        4194 :   long i, j = 1, k = n, l = lg(d);
    3459        4194 :   pari_sp av = avma;
    3460        4194 :   char *T = stack_calloc(n+1);
    3461        4194 :   for (i = 1; i < l; i++) T[d[i]] = 1;
    3462       56635 :   for (i = 1; i <= n; i++)
    3463       52441 :     if (T[i]) y[j++] = i; else y[k--] = i;
    3464        4194 :   avma = av; return y;
    3465             : }
    3466             : 
    3467             : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
    3468             : static GEN
    3469       15110 : indexrank0(long n, long r, GEN d)
    3470             : {
    3471       15110 :   GEN p1, p2, res = cgetg(3,t_VEC);
    3472             :   long i, j;
    3473             : 
    3474       15110 :   r = n - r; /* now r = dim Im(x) */
    3475       15110 :   p1 = cgetg(r+1,t_VECSMALL); gel(res,1) = p1;
    3476       15110 :   p2 = cgetg(r+1,t_VECSMALL); gel(res,2) = p2;
    3477       15110 :   if (d)
    3478             :   {
    3479       96261 :     for (i=0,j=1; j<=n; j++)
    3480       81865 :       if (d[j]) { i++; p1[i] = d[j]; p2[i] = j; }
    3481       14396 :     vecsmall_sort(p1);
    3482             :   }
    3483       15110 :   return res;
    3484             : }
    3485             : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
    3486             : static GEN
    3487         819 : indeximage0(long n, long r, GEN d)
    3488             : {
    3489             :   long i, j;
    3490             :   GEN v;
    3491             : 
    3492         819 :   r = n - r; /* now r = dim Im(x) */
    3493         819 :   v = cgetg(r+1,t_VECSMALL);
    3494        7525 :   if (d) for (i=j=1; j<=n; j++)
    3495        6706 :     if (d[j]) v[i++] = j;
    3496         819 :   return v;
    3497             : }
    3498             : /* x an m x n t_MAT, n > 0, r = dim Ker(x), d from gauss_pivot */
    3499             : static void
    3500        2097 : indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol)
    3501             : {
    3502        2097 :   GEN IR = indexrank0(n, r, d);
    3503        2097 :   *prow = perm_complete(gel(IR,1), m);
    3504        2097 :   *pcol = perm_complete(gel(IR,2), n);
    3505        2097 : }
    3506             : static void
    3507       13832 : init_indexrank(GEN x) {
    3508       13832 :   (void)new_chunk(3 + 2*lg(x)); /* HACK */
    3509       13832 : }
    3510             : 
    3511             : GEN
    3512          77 : indexrank(GEN x) {
    3513          77 :   pari_sp av = avma;
    3514             :   long r;
    3515          77 :   GEN d, p = NULL;
    3516          77 :   if (typ(x)!=t_MAT) pari_err_TYPE("indexrank",x);
    3517          77 :   init_indexrank(x);
    3518          77 :   if (RgM_is_FpM(x, &p) && p)
    3519          28 :   {
    3520             :     ulong pp;
    3521          28 :     x = RgM_Fp_init(x,p,&pp);
    3522          28 :     switch(pp)
    3523             :     {
    3524           7 :     case 0: d = FpM_gauss_pivot(x,p,&r); break;
    3525           7 :     case 2: d = F2m_gauss_pivot(x,&r); break;
    3526          14 :     default:d = Flm_gauss_pivot(x,pp,&r); break;
    3527             :     }
    3528             :   }
    3529             :   else
    3530          49 :     d = gauss_pivot(x,&r);
    3531          77 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3532             : }
    3533             : 
    3534             : GEN
    3535          21 : FpM_indexrank(GEN x, GEN p) {
    3536          21 :   pari_sp av = avma;
    3537             :   long r;
    3538             :   GEN d;
    3539          21 :   init_indexrank(x);
    3540          21 :   d = FpM_gauss_pivot(x,p,&r);
    3541          21 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3542             : }
    3543             : GEN
    3544        9331 : Flm_indexrank(GEN x, ulong p) {
    3545        9331 :   pari_sp av = avma;
    3546             :   long r;
    3547             :   GEN d;
    3548        9331 :   init_indexrank(x);
    3549        9331 :   d = Flm_gauss_pivot(Flm_copy(x),p,&r);
    3550        9331 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3551             : }
    3552             : GEN
    3553           0 : F2m_indexrank(GEN x) {
    3554           0 :   pari_sp av = avma;
    3555             :   long r;
    3556             :   GEN d;
    3557           0 :   init_indexrank(x);
    3558           0 :   d = F2m_gauss_pivot(F2m_copy(x),&r);
    3559           0 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3560             : }
    3561             : 
    3562             : GEN
    3563         819 : ZM_indeximage(GEN x) {
    3564         819 :   pari_sp av = avma;
    3565             :   long r;
    3566             :   GEN d;
    3567         819 :   init_indexrank(x);
    3568         819 :   d = ZM_pivots(x,&r);
    3569         819 :   avma = av; return indeximage0(lg(x)-1, r, d);
    3570             : }
    3571             : long
    3572       34592 : ZM_rank(GEN x) {
    3573       34592 :   pari_sp av = avma;
    3574             :   long r;
    3575       34592 :   (void)ZM_pivots(x,&r);
    3576       34592 :   avma = av; return lg(x)-1-r;
    3577             : }
    3578             : GEN
    3579        3584 : ZM_indexrank(GEN x) {
    3580        3584 :   pari_sp av = avma;
    3581             :   long r;
    3582             :   GEN d;
    3583        3584 :   init_indexrank(x);
    3584        3584 :   d = ZM_pivots(x,&r);
    3585        3584 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3586             : }
    3587             : 
    3588             : /*******************************************************************/
    3589             : /*                                                                 */
    3590             : /*                   Structured Elimination                        */
    3591             : /*                                                                 */
    3592             : /*******************************************************************/
    3593             : 
    3594             : static void
    3595       90853 : rem_col(GEN c, long i, GEN iscol, GEN Wrow, long *rcol, long *rrow)
    3596             : {
    3597       90853 :   long lc = lg(c), k;
    3598       90853 :   iscol[i] = 0; (*rcol)--;
    3599      805784 :   for (k = 1; k < lc; ++k)
    3600             :   {
    3601      714931 :     Wrow[c[k]]--;
    3602      714931 :     if (Wrow[c[k]]==0) (*rrow)--;
    3603             :   }
    3604       90853 : }
    3605             : 
    3606             : static void
    3607        5432 : rem_singleton(GEN M, GEN iscol, GEN Wrow, long *rcol, long *rrow)
    3608             : {
    3609             :   long i, j;
    3610        5432 :   long nbcol = lg(iscol)-1, last;
    3611             :   do
    3612             :   {
    3613        7084 :     last = 0;
    3614    16086497 :     for (i = 1; i <= nbcol; ++i)
    3615    16079413 :       if (iscol[i])
    3616             :       {
    3617     8670011 :         GEN c = gmael(M, i, 1);
    3618     8670011 :         long lc = lg(c);
    3619    80693291 :         for (j = 1; j < lc; ++j)
    3620    72034011 :           if (Wrow[c[j]] == 1)
    3621             :           {
    3622       10731 :             rem_col(c, i, iscol, Wrow, rcol, rrow);
    3623       10731 :             last=1; break;
    3624             :           }
    3625             :       }
    3626        7084 :   } while (last);
    3627        5432 : }
    3628             : 
    3629             : static GEN
    3630        5320 : fill_wcol(GEN M, GEN iscol, GEN Wrow, long *w, GEN wcol)
    3631             : {
    3632        5320 :   long nbcol = lg(iscol)-1;
    3633             :   long i, j, m, last;
    3634             :   GEN per;
    3635       12950 :   for (m = 2, last=0; !last ; m++)
    3636             :   {
    3637    17345083 :     for (i = 1; i <= nbcol; ++i)
    3638             :     {
    3639    17337453 :       wcol[i] = 0;
    3640    17337453 :       if (iscol[i])
    3641             :       {
    3642     9300690 :         GEN c = gmael(M, i, 1);
    3643     9300690 :         long lc = lg(c);
    3644    86334969 :         for (j = 1; j < lc; ++j)
    3645    77034279 :           if (Wrow[c[j]] == m) {  wcol[i]++; last = 1; }
    3646             :       }
    3647             :     }
    3648             :   }
    3649        5320 :   per = vecsmall_indexsort(wcol);
    3650        5320 :   *w = wcol[per[nbcol]];
    3651        5320 :   return per;
    3652             : }
    3653             : 
    3654             : /* M is a RgMs with nbrow rows, A a list of row indices.
    3655             :    Eliminate rows of M with a single entry that do not belong to A,
    3656             :    and the corresponding columns. Also eliminate columns until #colums=#rows.
    3657             :    Return pcol and prow:
    3658             :    pcol is a map from the new columns indices to the old one.
    3659             :    prow is a map from the old rows indices to the new one (0 if removed).
    3660             : */
    3661             : 
    3662             : void
    3663         112 : RgMs_structelim_col(GEN M, long nbcol, long nbrow, GEN A, GEN *p_col, GEN *p_row)
    3664             : {
    3665             :   long i,j,k;
    3666         112 :   long lA = lg(A);
    3667         112 :   GEN prow = cgetg(nbrow+1, t_VECSMALL);
    3668         112 :   GEN pcol = zero_zv(nbcol);
    3669         112 :   pari_sp av = avma;
    3670         112 :   long rcol = nbcol, rrow = 0, imin = nbcol - usqrt(nbcol);
    3671         112 :   GEN iscol = const_vecsmall(nbcol, 1);
    3672         112 :   GEN Wrow  = zero_zv(nbrow);
    3673         112 :   GEN wcol = cgetg(nbcol+1, t_VECSMALL);
    3674         112 :   pari_sp av2=avma;
    3675      115822 :   for (i = 1; i <= nbcol; ++i)
    3676             :   {
    3677      115710 :     GEN F = gmael(M, i, 1);
    3678      115710 :     long l = lg(F)-1;
    3679     1024709 :     for (j = 1; j <= l; ++j)
    3680      908999 :       Wrow[F[j]]++;
    3681             :   }
    3682         112 :   for (j = 1; j < lA; ++j)
    3683             :   {
    3684         112 :     if (Wrow[A[j]] == 0) { *p_col=NULL; return; }
    3685           0 :     Wrow[A[j]] = -1;
    3686             :   }
    3687      192031 :   for (i = 1; i <= nbrow; ++i)
    3688      191919 :     if (Wrow[i])
    3689       62489 :       rrow++;
    3690         112 :   rem_singleton(M, iscol, Wrow, &rcol, &rrow);
    3691         112 :   if (rcol<rrow) pari_err_BUG("RgMs_structelim, rcol<rrow");
    3692        5544 :   for (; rcol>rrow;)
    3693             :   {
    3694             :     long w;
    3695        5320 :     GEN per = fill_wcol(M, iscol, Wrow, &w, wcol);
    3696       85442 :     for (i = nbcol; i>=imin && wcol[per[i]]>=w && rcol>rrow; i--)
    3697       80122 :       rem_col(gmael(M, per[i], 1), per[i], iscol, Wrow, &rcol, &rrow);
    3698        5320 :     rem_singleton(M, iscol, Wrow, &rcol, &rrow);
    3699        5320 :     avma = av2;
    3700             :   }
    3701      115822 :   for (j = 1, i = 1; i <= nbcol; ++i)
    3702      115710 :     if (iscol[i])
    3703       24857 :       pcol[j++] = i;
    3704         112 :   setlg(pcol,j);
    3705      192031 :   for (k = 1, i = 1; i <= nbrow; ++i)
    3706      191919 :     prow[i] = Wrow[i] ? k++: 0;
    3707         112 :   avma = av;
    3708         112 :   *p_col = pcol; *p_row = prow;
    3709             : }
    3710             : 
    3711             : void
    3712           0 : RgMs_structelim(GEN M, long nbrow, GEN A, GEN *p_col, GEN *p_row)
    3713             : {
    3714           0 :   RgMs_structelim_col(M, lg(M)-1, nbrow, A, p_col, p_row);
    3715           0 : }
    3716             : 
    3717             : /*******************************************************************/
    3718             : /*                                                                 */
    3719             : /*                        EIGENVECTORS                             */
    3720             : /*   (independent eigenvectors, sorted by increasing eigenvalue)   */
    3721             : /*                                                                 */
    3722             : /*******************************************************************/
    3723             : 
    3724             : GEN
    3725          63 : mateigen(GEN x, long flag, long prec)
    3726             : {
    3727             :   GEN y, R, T;
    3728          63 :   long k, l, ex, n = lg(x);
    3729          63 :   pari_sp av = avma;
    3730             : 
    3731          63 :   if (typ(x)!=t_MAT) pari_err_TYPE("eigen",x);
    3732          63 :   if (n != 1 && n != lgcols(x)) pari_err_DIM("eigen");
    3733          63 :   if (flag < 0 || flag > 1) pari_err_FLAG("mateigen");
    3734          63 :   if (n == 1)
    3735             :   {
    3736          14 :     if (flag) retmkvec2(cgetg(1,t_VEC), cgetg(1,t_MAT));
    3737           7 :     return cgetg(1,t_VEC);
    3738             :   }
    3739          49 :   if (n == 2)
    3740             :   {
    3741          14 :     if (flag) retmkvec2(mkveccopy(gcoeff(x,1,1)), matid(1));
    3742           7 :     return matid(1);
    3743             :   }
    3744             : 
    3745          35 :   ex = 16 - prec2nbits(prec);
    3746          35 :   T = charpoly(x,0);
    3747          35 :   if (RgX_is_QX(T))
    3748             :   {
    3749          28 :     T = Q_primpart(T);
    3750          28 :     (void)ZX_gcd_all(T, ZX_deriv(T),  &T);
    3751          28 :     R = nfrootsQ(T);
    3752          28 :     if (lg(R)-1 < degpol(T))
    3753             :     { /* add missing complex roots */
    3754          14 :       GEN r = cleanroots(RgX_div(T, roots_to_pol(R, 0)), prec);
    3755          14 :       settyp(r, t_VEC);
    3756          14 :       R = shallowconcat(R, r);
    3757             :     }
    3758             :   }
    3759             :   else
    3760             :   {
    3761           7 :     GEN r1, v = vectrunc_init(lg(T));
    3762             :     long e;
    3763           7 :     R = cleanroots(T,prec);
    3764           7 :     r1 = NULL;
    3765          21 :     for (k = 1; k < lg(R); k++)
    3766             :     {
    3767          14 :       GEN r2 = gel(R,k), r = grndtoi(r2, &e);
    3768          14 :       if (e < ex) r2 = r;
    3769          14 :       if (r1)
    3770             :       {
    3771           7 :         r = gsub(r1,r2);
    3772           7 :         if (gequal0(r) || gexpo(r) < ex) continue;
    3773             :       }
    3774          14 :       vectrunc_append(v, r2);
    3775          14 :       r1 = r2;
    3776             :     }
    3777           7 :     R = v;
    3778             :   }
    3779             :   /* R = distinct complex roots of charpoly(x) */
    3780          35 :   l = lg(R); y = cgetg(l, t_VEC);
    3781         189 :   for (k = 1; k < l; k++)
    3782             :   {
    3783         154 :     GEN F = ker_aux(RgM_Rg_sub_shallow(x, gel(R,k)), x);
    3784         154 :     long d = lg(F)-1;
    3785         154 :     if (!d) pari_err_PREC("mateigen");
    3786         154 :     gel(y,k) = F;
    3787         154 :     if (flag) gel(R,k) = const_vec(d, gel(R,k));
    3788             :   }
    3789          35 :   y = shallowconcat1(y);
    3790          35 :   if (lg(y) > n) pari_err_PREC("mateigen");
    3791             :   /* lg(y) < n if x is not diagonalizable */
    3792          35 :   if (flag) y = mkvec2(shallowconcat1(R), y);
    3793          35 :   return gerepilecopy(av,y);
    3794             : }
    3795             : GEN
    3796           0 : eigen(GEN x, long prec) { return mateigen(x, 0, prec); }
    3797             : 
    3798             : /*******************************************************************/
    3799             : /*                                                                 */
    3800             : /*                           DETERMINANT                           */
    3801             : /*                                                                 */
    3802             : /*******************************************************************/
    3803             : 
    3804             : GEN
    3805        1610 : det0(GEN a,long flag)
    3806             : {
    3807        1610 :   switch(flag)
    3808             :   {
    3809        1596 :     case 0: return det(a);
    3810          14 :     case 1: return det2(a);
    3811           0 :     default: pari_err_FLAG("matdet");
    3812             :   }
    3813           0 :   return NULL; /* not reached */
    3814             : }
    3815             : 
    3816             : /* M a 2x2 matrix, returns det(M) */
    3817             : static GEN
    3818        2312 : RgM_det2(GEN M)
    3819             : {
    3820        2312 :   pari_sp av = avma;
    3821        2312 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    3822        2312 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    3823        2312 :   return gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
    3824             : }
    3825             : /* M a 2x2 ZM, returns det(M) */
    3826             : static GEN
    3827        7434 : ZM_det2(GEN M)
    3828             : {
    3829        7434 :   pari_sp av = avma;
    3830        7434 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    3831        7434 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    3832        7434 :   return gerepileuptoint(av, subii(mulii(a,d), mulii(b, c)));
    3833             : }
    3834             : /* M a 3x3 ZM, return det(M) */
    3835             : static GEN
    3836        3003 : ZM_det3(GEN M)
    3837             : {
    3838        3003 :   pari_sp av = avma;
    3839        3003 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2), c = gcoeff(M,1,3);
    3840        3003 :   GEN d = gcoeff(M,2,1), e = gcoeff(M,2,2), f = gcoeff(M,2,3);
    3841        3003 :   GEN g = gcoeff(M,3,1), h = gcoeff(M,3,2), i = gcoeff(M,3,3);
    3842        3003 :   GEN t, D = signe(i)? mulii(subii(mulii(a,e), mulii(b,d)), i): gen_0;
    3843        3003 :   if (signe(g))
    3844             :   {
    3845        2821 :     t = mulii(subii(mulii(b,f), mulii(c,e)), g);
    3846        2821 :     D = addii(D, t);
    3847             :   }
    3848        3003 :   if (signe(h))
    3849             :   {
    3850        2856 :     t = mulii(subii(mulii(c,d), mulii(a,f)), h);
    3851        2856 :     D = addii(D, t);
    3852             :   }
    3853        3003 :   return gerepileuptoint(av, D);
    3854             : }
    3855             : 
    3856             : static GEN
    3857        1467 : det_simple_gauss(GEN a, GEN data, pivot_fun pivot)
    3858             : {
    3859        1467 :   pari_sp av = avma;
    3860        1467 :   long i,j,k, s = 1, nbco = lg(a)-1;
    3861        1467 :   GEN p, x = gen_1;
    3862             : 
    3863        1467 :   a = RgM_shallowcopy(a);
    3864        6353 :   for (i=1; i<nbco; i++)
    3865             :   {
    3866        4893 :     k = pivot(a, data, i, NULL);
    3867        4893 :     if (k > nbco) return gerepilecopy(av, gcoeff(a,i,i));
    3868        4886 :     if (k != i)
    3869             :     { /* exchange the lines s.t. k = i */
    3870        1935 :       for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    3871        1935 :       s = -s;
    3872             :     }
    3873        4886 :     p = gcoeff(a,i,i);
    3874             : 
    3875        4886 :     x = gmul(x,p);
    3876       17621 :     for (k=i+1; k<=nbco; k++)
    3877             :     {
    3878       12735 :       GEN m = gcoeff(a,i,k);
    3879       12735 :       if (gequal0(m)) continue;
    3880             : 
    3881       12574 :       m = gdiv(m,p);
    3882       61725 :       for (j=i+1; j<=nbco; j++)
    3883       49151 :         gcoeff(a,j,k) = gsub(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
    3884             :     }
    3885        4886 :     if (gc_needed(av,2))
    3886             :     {
    3887           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
    3888           0 :       gerepileall(av,2, &a,&x);
    3889             :     }
    3890             :   }
    3891        1460 :   if (s < 0) x = gneg_i(x);
    3892        1460 :   return gerepileupto(av, gmul(x, gcoeff(a,nbco,nbco)));
    3893             : }
    3894             : 
    3895             : GEN
    3896        2722 : det2(GEN a)
    3897             : {
    3898             :   GEN data;
    3899             :   pivot_fun pivot;
    3900        2722 :   long n = lg(a)-1;
    3901        2722 :   if (typ(a)!=t_MAT) pari_err_TYPE("det2",a);
    3902        2722 :   if (!n) return gen_1;
    3903        2722 :   if (n != nbrows(a)) pari_err_DIM("det2");
    3904        2722 :   if (n == 1) return gcopy(gcoeff(a,1,1));
    3905        2722 :   if (n == 2) return RgM_det2(a);
    3906        1425 :   pivot = get_pivot_fun(a, a, &data);
    3907        1425 :   return det_simple_gauss(a, data, pivot);
    3908             : }
    3909             : 
    3910             : static GEN
    3911        2268 : mydiv(GEN x, GEN y)
    3912             : {
    3913        2268 :   long tx = typ(x), ty = typ(y);
    3914        2268 :   if (tx == ty && tx == t_POL && varn(x) == varn(y)) return RgX_div(x,y);
    3915        1715 :   return gdiv(x,y);
    3916             : }
    3917             : 
    3918             : /* Assumes a a square t_MAT of dimension n > 0. Returns det(a) using
    3919             :  * Gauss-Bareiss. */
    3920             : static GEN
    3921         294 : det_bareiss(GEN a)
    3922             : {
    3923         294 :   pari_sp av = avma;
    3924         294 :   long nbco = lg(a)-1,i,j,k,s = 1;
    3925             :   GEN p, pprec;
    3926             : 
    3927         294 :   a = RgM_shallowcopy(a);
    3928        1113 :   for (pprec=gen_1,i=1; i<nbco; i++,pprec=p)
    3929             :   {
    3930             :     GEN ci;
    3931         819 :     int diveuc = (gequal1(pprec)==0);
    3932             : 
    3933         819 :     p = gcoeff(a,i,i);
    3934         819 :     if (gequal0(p))
    3935             :     {
    3936           0 :       k=i+1; while (k<=nbco && gequal0(gcoeff(a,i,k))) k++;
    3937           0 :       if (k>nbco) return gerepilecopy(av, p);
    3938           0 :       swap(gel(a,k), gel(a,i)); s = -s;
    3939           0 :       p = gcoeff(a,i,i);
    3940             :     }
    3941         819 :     ci = gel(a,i);
    3942        2618 :     for (k=i+1; k<=nbco; k++)
    3943             :     {
    3944        1799 :       GEN ck = gel(a,k), m = gel(ck,i);
    3945        1799 :       if (gequal0(m))
    3946             :       {
    3947           7 :         if (gequal1(p))
    3948             :         {
    3949           0 :           if (diveuc)
    3950           0 :             gel(a,k) = mydiv(gel(a,k), pprec);
    3951             :         }
    3952             :         else
    3953          42 :           for (j=i+1; j<=nbco; j++)
    3954             :           {
    3955          35 :             GEN p1 = gmul(p, gel(ck,j));
    3956          35 :             if (diveuc) p1 = mydiv(p1,pprec);
    3957          35 :             gel(ck,j) = p1;
    3958             :           }
    3959             :       }
    3960             :       else
    3961        7056 :         for (j=i+1; j<=nbco; j++)
    3962             :         {
    3963        5264 :           pari_sp av2 = avma;
    3964        5264 :           GEN p1 = gsub(gmul(p,gel(ck,j)), gmul(m,gel(ci,j)));
    3965        5264 :           if (diveuc) p1 = mydiv(p1,pprec);
    3966        5264 :           gel(ck,j) = gerepileupto(av2, p1);
    3967             :         }
    3968        1799 :       if (gc_needed(av,2))
    3969             :       {
    3970           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
    3971           0 :         gerepileall(av,2, &a,&pprec);
    3972           0 :         ci = gel(a,i);
    3973           0 :         p = gcoeff(a,i,i);
    3974             :       }
    3975             :     }
    3976             :   }
    3977         294 :   p = gcoeff(a,nbco,nbco);
    3978         294 :   p = (s < 0)? gneg(p): gcopy(p);
    3979         294 :   return gerepileupto(av, p);
    3980             : }
    3981             : 
    3982             : /* count non-zero entries in col j, at most 'max' of them.
    3983             :  * Return their indices */
    3984             : static GEN
    3985        1232 : col_count_non_zero(GEN a, long j, long max)
    3986             : {
    3987        1232 :   GEN v = cgetg(max+1, t_VECSMALL);
    3988        1232 :   GEN c = gel(a,j);
    3989        1232 :   long i, l = lg(a), k = 1;
    3990        5663 :   for (i = 1; i < l; i++)
    3991        5348 :     if (!gequal0(gel(c,i)))
    3992             :     {
    3993        4396 :       if (k > max) return NULL; /* fail */
    3994        3479 :       v[k++] = i;
    3995             :     }
    3996         315 :   setlg(v, k); return v;
    3997             : }
    3998             : /* count non-zero entries in row i, at most 'max' of them.
    3999             :  * Return their indices */
    4000             : static GEN
    4001        1057 : row_count_non_zero(GEN a, long i, long max)
    4002             : {
    4003        1057 :   GEN v = cgetg(max+1, t_VECSMALL);
    4004        1057 :   long j, l = lg(a), k = 1;
    4005        4760 :   for (j = 1; j < l; j++)
    4006        4550 :     if (!gequal0(gcoeff(a,i,j)))
    4007             :     {
    4008        4144 :       if (k > max) return NULL; /* fail */
    4009        3297 :       v[k++] = j;
    4010             :     }
    4011         210 :   setlg(v, k); return v;
    4012             : }
    4013             : 
    4014             : static GEN det_develop(GEN a, long max, double bound);
    4015             : /* (-1)^(i+j) a[i,j] * det RgM_minor(a,i,j) */
    4016             : static GEN
    4017         455 : coeff_det(GEN a, long i, long j, long max, double bound)
    4018             : {
    4019         455 :   GEN c = gcoeff(a, i, j);
    4020         455 :   c = gmul(c, det_develop(RgM_minor(a, i,j), max, bound));
    4021         455 :   if (odd(i+j)) c = gneg(c);
    4022         455 :   return c;
    4023             : }
    4024             : /* a square t_MAT, 'bound' a rough upper bound for the number of
    4025             :  * multiplications we are willing to pay while developing rows/columns before
    4026             :  * switching to Gaussian elimination */
    4027             : static GEN
    4028         644 : det_develop(GEN M, long max, double bound)
    4029             : {
    4030         644 :   pari_sp av = avma;
    4031         644 :   long i,j, n = lg(M)-1, lbest = max+2, best_col = 0, best_row = 0;
    4032         644 :   GEN best = NULL;
    4033             : 
    4034         644 :   if (bound < 1.) return det_bareiss(M); /* too costly now */
    4035             : 
    4036         462 :   switch(n)
    4037             :   {
    4038           0 :     case 0: return gen_1;
    4039           0 :     case 1: return gcopy(gcoeff(M,1,1));
    4040          77 :     case 2: return RgM_det2(M);
    4041             :   }
    4042         385 :   if (max > ((n+2)>>1)) max = (n+2)>>1;
    4043        1470 :   for (j = 1; j <= n; j++)
    4044             :   {
    4045        1232 :     pari_sp av2 = avma;
    4046        1232 :     GEN v = col_count_non_zero(M, j, max);
    4047             :     long lv;
    4048        1232 :     if (!v || (lv = lg(v)) >= lbest) { avma = av2; continue; }
    4049         273 :     if (lv == 1) { avma = av; return gen_0; }
    4050         273 :     if (lv == 2) {
    4051         147 :       avma = av;
    4052         147 :       return gerepileupto(av, coeff_det(M,v[1],j,max,bound));
    4053             :     }
    4054         126 :     best = v; lbest = lv; best_col = j;
    4055             :   }
    4056        1295 :   for (i = 1; i <= n; i++)
    4057             :   {
    4058        1057 :     pari_sp av2 = avma;
    4059        1057 :     GEN v = row_count_non_zero(M, i, max);
    4060             :     long lv;
    4061        1057 :     if (!v || (lv = lg(v)) >= lbest) { avma = av2; continue; }
    4062          42 :     if (lv == 1) { avma = av; return gen_0; }
    4063          42 :     if (lv == 2) {
    4064           0 :       avma = av;
    4065           0 :       return gerepileupto(av, coeff_det(M,i,v[1],max,bound));
    4066             :     }
    4067          42 :     best = v; lbest = lv; best_row = i;
    4068             :   }
    4069         238 :   if (best_row)
    4070             :   {
    4071          42 :     double d = lbest-1;
    4072          42 :     GEN s = NULL;
    4073             :     long k;
    4074          42 :     bound /= d*d*d;
    4075         168 :     for (k = 1; k < lbest; k++)
    4076             :     {
    4077         126 :       GEN c = coeff_det(M, best_row, best[k], max, bound);
    4078         126 :       s = s? gadd(s, c): c;
    4079             :     }
    4080          42 :     return gerepileupto(av, s);
    4081             :   }
    4082         196 :   if (best_col)
    4083             :   {
    4084          84 :     double d = lbest-1;
    4085          84 :     GEN s = NULL;
    4086             :     long k;
    4087          84 :     bound /= d*d*d;
    4088         266 :     for (k = 1; k < lbest; k++)
    4089             :     {
    4090         182 :       GEN c = coeff_det(M, best[k], best_col, max, bound);
    4091         182 :       s = s? gadd(s, c): c;
    4092             :     }
    4093          84 :     return gerepileupto(av, s);
    4094             :   }
    4095         112 :   return det_bareiss(M);
    4096             : }
    4097             : 
    4098             : /* area of parallelogram bounded by (v1,v2) */
    4099             : static GEN
    4100       52507 : parallelogramarea(GEN v1, GEN v2)
    4101       52507 : { return gsub(gmul(gnorml2(v1), gnorml2(v2)), gsqr(RgV_dotproduct(v1, v2))); }
    4102             : 
    4103             : /* Square of Hadamard bound for det(a), a square matrix.
    4104             :  * Slightly improvement: instead of using the column norms, use the area of
    4105             :  * the parallelogram formed by pairs of consecutive vectors */
    4106             : GEN
    4107       17101 : RgM_Hadamard(GEN a)
    4108             : {
    4109       17101 :   pari_sp av = avma;
    4110       17101 :   long n = lg(a)-1, i;
    4111             :   GEN B;
    4112       17101 :   if (n == 0) return gen_1;
    4113       17101 :   if (n == 1) return gsqr(gcoeff(a,1,1));
    4114       17101 :   a = RgM_gtofp(a, LOWDEFAULTPREC);
    4115       17101 :   B = gen_1;
    4116       69608 :   for (i = 1; i <= n/2; i++)
    4117       52507 :     B = gmul(B, parallelogramarea(gel(a,2*i-1), gel(a,2*i)));
    4118       17101 :   if (odd(n)) B = gmul(B, gnorml2(gel(a, n)));
    4119       17101 :   return gerepileuptoint(av, ceil_safe(B));
    4120             : }
    4121             : 
    4122             : /* assume dim(a) = n > 0 */
    4123             : static GEN
    4124       28357 : ZM_det_i(GEN M, long n)
    4125             : {
    4126       28357 :   const long DIXON_THRESHOLD = 40;
    4127       28357 :   pari_sp av = avma, av2;
    4128             :   long i;
    4129       28357 :   ulong p, compp, Dp = 1;
    4130             :   forprime_t S;
    4131             :   GEN D, h, q, v, comp;
    4132       28357 :   if (n == 1) return icopy(gcoeff(M,1,1));
    4133       27538 :   if (n == 2) return ZM_det2(M);
    4134       20104 :   if (n == 3) return ZM_det3(M);
    4135       17101 :   h = RgM_Hadamard(M);
    4136       17101 :   if (!signe(h)) { avma = av; return gen_0; }
    4137       17101 :   h = sqrti(h); q = gen_1;
    4138       17101 :   init_modular_big(&S);
    4139       17101 :   p = 0; /* -Wall */
    4140       34202 :   while( cmpii(q, h) <= 0 && (p = u_forprime_next(&S)) )
    4141             :   {
    4142       17101 :     av2 = avma; Dp = Flm_det(ZM_to_Flm(M, p), p);
    4143       17101 :     avma = av2;
    4144       17101 :     if (Dp) break;
    4145           0 :     q = muliu(q, p);
    4146             :   }
    4147       17101 :   if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
    4148       17101 :   if (!Dp) { avma = av; return gen_0; }
    4149       17101 :   if (n <= DIXON_THRESHOLD)
    4150       17101 :     D = q;
    4151             :   else
    4152             :   {
    4153           0 :     av2 = avma;
    4154           0 :     v = cgetg(n+1, t_COL);
    4155           0 :     gel(v, 1) = gen_1; /* ensure content(v) = 1 */
    4156           0 :     for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
    4157           0 :     D = Q_denom(ZM_gauss(M, v));
    4158           0 :     if (expi(D) < expi(h) >> 1)
    4159             :     { /* First try unlucky, try once more */
    4160           0 :       for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
    4161           0 :       D = lcmii(D, Q_denom(ZM_gauss(M, v)));
    4162             :     }
    4163           0 :     D = gerepileuptoint(av2, D);
    4164           0 :     if (q != gen_1) D = lcmii(D, q);
    4165             :   }
    4166             :   /* determinant is M multiple of D */
    4167       17101 :   h = shifti(divii(h, D), 1);
    4168             : 
    4169       17101 :   compp = Fl_div(Dp, umodiu(D,p), p);
    4170       17101 :   comp = Z_init_CRT(compp, p);
    4171       17101 :   q = utoipos(p);
    4172       86257 :   while (cmpii(q, h) <= 0)
    4173             :   {
    4174       52055 :     p = u_forprime_next(&S);
    4175       52055 :     if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
    4176       52055 :     Dp = umodiu(D, p);
    4177       52055 :     if (!Dp) continue;
    4178       52055 :     av2 = avma;
    4179       52055 :     compp = Fl_div(Flm_det(ZM_to_Flm(M, p), p), Dp, p);
    4180       52055 :     avma = av2;
    4181       52055 :     (void) Z_incremental_CRT(&comp, compp, &q, p);
    4182             :   }
    4183       17101 :   return gerepileuptoint(av, mulii(comp, D));
    4184             : }
    4185             : 
    4186             : static long
    4187         189 : det_init_max(long n)
    4188             : {
    4189         189 :   if (n > 100) return 0;
    4190         189 :   if (n > 50) return 1;
    4191         189 :   if (n > 30) return 4;
    4192         189 :   return 7;
    4193             : }
    4194             : 
    4195             : GEN
    4196        3464 : det(GEN a)
    4197             : {
    4198        3464 :   long n = lg(a)-1;
    4199             :   double B;
    4200        3464 :   GEN data, ff = NULL, p = NULL;
    4201             :   pivot_fun pivot;
    4202             : 
    4203        3464 :   if (typ(a)!=t_MAT) pari_err_TYPE("det",a);
    4204        3464 :   if (!n) return gen_1;
    4205        3415 :   if (n != nbrows(a)) pari_err_DIM("det");
    4206        3408 :   if (n == 1) return gcopy(gcoeff(a,1,1));
    4207        3129 :   if (n == 2) return RgM_det2(a);
    4208        2191 :   if (RgM_is_FpM(a, &p))
    4209             :   {
    4210        1939 :     pari_sp av = avma;
    4211             :     ulong pp, d;
    4212        1939 :     if (!p) return ZM_det_i(a, n); /* ZM */
    4213             :     /* FpM */
    4214        1477 :     a = RgM_Fp_init(a,p,&pp);
    4215        1477 :     switch(pp)
    4216             :     {
    4217          49 :     case 0: return gerepileupto(av, Fp_to_mod(FpM_det(a,p),p)); break;
    4218          14 :     case 2: d = F2m_det(a); break;
    4219        1414 :     default:d = Flm_det(a,pp); break;
    4220             :     }
    4221        1428 :     avma = av; return mkintmodu(d, pp);
    4222             :   }
    4223         252 :   if (RgM_is_FFM(a, &ff)) return FFM_det(a, ff);
    4224         231 :   pivot = get_pivot_fun(a, a, &data);
    4225         231 :   if (pivot != gauss_get_pivot_NZ) return det_simple_gauss(a, data, pivot);
    4226         189 :   B = (double)n;
    4227         189 :   return det_develop(a, det_init_max(n), B*B*B);
    4228             : }
    4229             : 
    4230             : GEN
    4231       27895 : ZM_det(GEN a)
    4232             : {
    4233       27895 :   long n = lg(a)-1;
    4234       27895 :   if (!n) return gen_1;
    4235       27895 :   return ZM_det_i(a, n);
    4236             : }
    4237             : 
    4238             : /* return a solution of congruence system sum M_{ i,j } X_j = Y_i mod D_i
    4239             :  * If ptu1 != NULL, put in *ptu1 a Z-basis of the homogeneous system */
    4240             : static GEN
    4241          49 : gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1)
    4242             : {
    4243          49 :   pari_sp av = avma;
    4244             :   long n, m, j, l, lM;
    4245             :   GEN delta, H, U, u1, u2, x;
    4246             : 
    4247          49 :   if (typ(M)!=t_MAT) pari_err_TYPE("gaussmodulo",M);
    4248          49 :   lM = lg(M);
    4249          49 :   if (lM == 1)
    4250             :   {
    4251          28 :     switch(typ(Y))
    4252             :     {
    4253           7 :       case t_INT: break;
    4254          14 :       case t_COL: if (lg(Y) != 1) pari_err_DIM("gaussmodulo");
    4255          14 :                   break;
    4256           7 :       default: pari_err_TYPE("gaussmodulo",Y);
    4257             :     }
    4258          21 :     switch(typ(D))
    4259             :     {
    4260           7 :       case t_INT: break;
    4261           7 :       case t_COL: if (lg(D) != 1) pari_err_DIM("gaussmodulo");
    4262           7 :                   break;
    4263           7 :       default: pari_err_TYPE("gaussmodulo",D);
    4264             :     }
    4265          14 :     if (ptu1) *ptu1 = cgetg(1, t_MAT);
    4266          14 :     return gen_0;
    4267             :   }
    4268          21 :   n = nbrows(M);
    4269          21 :   switch(typ(D))
    4270             :   {
    4271             :     case t_COL:
    4272          14 :       if (lg(D)-1!=n) pari_err_DIM("gaussmodulo");
    4273          14 :       delta = diagonal_shallow(D); break;
    4274           7 :     case t_INT: delta = scalarmat_shallow(D,n); break;
    4275           0 :     default: pari_err_TYPE("gaussmodulo",D);
    4276           0 :       return NULL; /* not reached */
    4277             :   }
    4278          21 :   switch(typ(Y))
    4279             :   {
    4280           7 :     case t_INT: Y = const_col(n, Y); break;
    4281             :     case t_COL:
    4282          14 :       if (lg(Y)-1!=n) pari_err_DIM("gaussmodulo");
    4283          14 :       break;
    4284           0 :     default: pari_err_TYPE("gaussmodulo",Y);
    4285           0 :       return NULL; /* not reached */
    4286             :   }
    4287          21 :   H = ZM_hnfall(shallowconcat(M,delta), &U, 1);
    4288          21 :   Y = hnf_solve(H,Y); if (!Y) return gen_0;
    4289          21 :   l = lg(H); /* may be smaller than lM if some moduli are 0 */
    4290          21 :   n = l-1;
    4291          21 :   m = lg(U)-1 - n;
    4292          21 :   u1 = cgetg(m+1,t_MAT);
    4293          21 :   u2 = cgetg(n+1,t_MAT);
    4294          21 :   for (j=1; j<=m; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u1,j) = c; }
    4295          21 :   U += m;
    4296          21 :   for (j=1; j<=n; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u2,j) = c; }
    4297             :   /*       (u1 u2)
    4298             :    * (M D) (*  * ) = (0 H) */
    4299          21 :   u1 = ZM_lll(u1, 0.75, LLL_INPLACE);
    4300          21 :   Y = ZM_ZC_mul(u2,Y);
    4301          21 :   x = ZC_reducemodmatrix(Y, u1);
    4302          21 :   if (!ptu1) x = gerepileupto(av, x);
    4303             :   else
    4304             :   {
    4305          14 :     gerepileall(av, 2, &x, &u1);
    4306          14 :     *ptu1 = u1;
    4307             :   }
    4308          21 :   return x;
    4309             : }
    4310             : 
    4311             : GEN
    4312          49 : matsolvemod0(GEN M, GEN D, GEN Y, long flag)
    4313             : {
    4314             :   pari_sp av;
    4315             :   GEN p1,y;
    4316             : 
    4317          49 :   if (!flag) return gaussmoduloall(M,D,Y,NULL);
    4318             : 
    4319          42 :   av=avma; y = cgetg(3,t_VEC);
    4320          42 :   p1 = gaussmoduloall(M,D,Y, (GEN*)y+2);
    4321          28 :   if (p1==gen_0) { avma=av; return gen_0; }
    4322          14 :   gel(y,1) = p1; return y;
    4323             : }
    4324             : 
    4325             : GEN
    4326           0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod0(M,D,Y,1); }
    4327             : 
    4328             : GEN
    4329           0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod0(M,D,Y,0); }

Generated by: LCOV version 1.11