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

Generated by: LCOV version 1.11