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.10.0 lcov report (development 20084-3366ab4) Lines: 2378 2528 94.1 %
Date: 2017-01-19 05:50:13 Functions: 217 232 93.5 %
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          16 : gerepile_mat(pari_sp av, pari_sp tetpil, GEN x, long k, long m, long n, long t)
      31             : {
      32          16 :   pari_sp A, bot = pari_mainstack->bot;
      33             :   long u, i;
      34             :   size_t dec;
      35             : 
      36          16 :   (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
      37             : 
      38         487 :   for (u=t+1; u<=m; u++)
      39             :   {
      40         471 :     A = (pari_sp)coeff(x,u,k);
      41         471 :     if (A < av && A >= bot) coeff(x,u,k) += dec;
      42             :   }
      43        1104 :   for (i=k+1; i<=n; i++)
      44       77248 :     for (u=1; u<=m; u++)
      45             :     {
      46       76160 :       A = (pari_sp)coeff(x,u,i);
      47       76160 :       if (A < av && A >= bot) coeff(x,u,i) += dec;
      48             :     }
      49          16 : }
      50             : 
      51             : static void
      52          16 : gen_gerepile_gauss_ker(GEN x, long k, long t, pari_sp av, void *E, GEN (*copy)(void*, GEN))
      53             : {
      54          16 :   pari_sp tetpil = avma;
      55          16 :   long u,i, n = lg(x)-1, m = n? nbrows(x): 0;
      56             : 
      57          16 :   if (DEBUGMEM > 1) pari_warn(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
      58          16 :   for (u=t+1; u<=m; u++) gcoeff(x,u,k) = copy(E,gcoeff(x,u,k));
      59        1104 :   for (i=k+1; i<=n; i++)
      60        1088 :     for (u=1; u<=m; u++) gcoeff(x,u,i) = copy(E,gcoeff(x,u,i));
      61          16 :   gerepile_mat(av,tetpil,x,k,m,n,t);
      62          16 : }
      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       28662 : _copy(void *E, GEN x)
      72             : {
      73       28662 :   (void) E; COPY(x);
      74       28662 :   return x;
      75             : }
      76             : 
      77             : static void
      78           6 : gerepile_gauss_ker(GEN x, long k, long t, pari_sp av)
      79             : {
      80           6 :   gen_gerepile_gauss_ker(x, k, t, av, NULL, &_copy);
      81           6 : }
      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        2319 : gen_ker(GEN x, long deplin, void *E, const struct bb_field *ff)
     121             : {
     122        2319 :   pari_sp av0 = avma, av, tetpil;
     123             :   GEN y, c, d;
     124             :   long i, j, k, r, t, n, m;
     125             : 
     126        2319 :   n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
     127        2319 :   m=nbrows(x); r=0;
     128        2319 :   x = RgM_shallowcopy(x);
     129        2319 :   c = zero_zv(m);
     130        2319 :   d=new_chunk(n+1);
     131        2319 :   av=avma;
     132       30125 :   for (k=1; k<=n; k++)
     133             :   {
     134      745224 :     for (j=1; j<=m; j++)
     135      730344 :       if (!c[j])
     136             :       {
     137      478611 :         gcoeff(x,j,k) = ff->red(E, gcoeff(x,j,k));
     138      478611 :         if (!ff->equal0(gcoeff(x,j,k))) break;
     139             :       }
     140       27813 :     if (j>m)
     141             :     {
     142       14880 :       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       14873 :       r++; d[k]=0;
     150      268697 :       for(j=1; j<k; j++)
     151      253824 :         if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
     152             :     }
     153             :     else
     154             :     {
     155       12933 :       GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
     156       12933 :       c[j] = k; d[k] = j;
     157       12933 :       gcoeff(x,j,k) = ff->s(E,-1);
     158       12933 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
     159      530869 :       for (t=1; t<=m; t++)
     160             :       {
     161      517936 :         if (t==j) continue;
     162             : 
     163      505003 :         piv = ff->red(E,gcoeff(x,t,k));
     164      505003 :         if (ff->equal0(piv)) continue;
     165             : 
     166       85408 :         gcoeff(x,t,k) = ff->s(E,0);
     167     1491215 :         for (i=k+1; i<=n; i++)
     168     1405807 :            gcoeff(x,t,i) = ff->add(E, gcoeff(x,t,i), ff->mul(E,piv,gcoeff(x,j,i)));
     169       85408 :         if (gc_needed(av,1))
     170          10 :           gen_gerepile_gauss_ker(x,k,t,av,E,ff->red);
     171             :       }
     172             :     }
     173             :   }
     174        2312 :   if (deplin) { avma = av0; return NULL; }
     175             : 
     176        2305 :   tetpil=avma; y=cgetg(r+1,t_MAT);
     177       17178 :   for (j=k=1; j<=r; j++,k++)
     178             :   {
     179       14873 :     GEN C = cgetg(n+1,t_COL);
     180       14873 :     GEN g0 = ff->s(E,0), g1 = ff->s(E,1);
     181       14873 :     gel(y,j) = C; while (d[k]) k++;
     182      268697 :     for (i=1; i<k; i++)
     183      253824 :       if (d[i])
     184             :       {
     185       52428 :         GEN p1=gcoeff(x,d[i],k);
     186       52428 :         gel(C,i) = ff->red(E,p1); gunclone(p1);
     187             :       }
     188             :       else
     189      201396 :         gel(C,i) = g0;
     190       14873 :     gel(C,k) = g1; for (i=k+1; i<=n; i++) gel(C,i) = g0;
     191             :   }
     192        2305 :   return gerepile(av0,tetpil,y);
     193             : }
     194             : 
     195             : GEN
     196         381 : gen_Gauss_pivot(GEN x, long *rr, void *E, const struct bb_field *ff)
     197             : {
     198             :   pari_sp av;
     199             :   GEN c, d;
     200         381 :   long i, j, k, r, t, m, n = lg(x)-1;
     201             : 
     202         381 :   if (!n) { *rr = 0; return NULL; }
     203             : 
     204         381 :   m=nbrows(x); r=0;
     205         381 :   d = cgetg(n+1, t_VECSMALL);
     206         381 :   x = RgM_shallowcopy(x);
     207         381 :   c = zero_zv(m);
     208         381 :   av=avma;
     209        5384 :   for (k=1; k<=n; k++)
     210             :   {
     211       94777 :     for (j=1; j<=m; j++)
     212       93677 :       if (!c[j])
     213             :       {
     214       42193 :         gcoeff(x,j,k) = ff->red(E,gcoeff(x,j,k));
     215       42193 :         if (!ff->equal0(gcoeff(x,j,k))) break;
     216             :       }
     217        5003 :     if (j>m) { r++; d[k]=0; }
     218             :     else
     219             :     {
     220        3903 :       GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
     221        3903 :       GEN g0 = ff->s(E,0);
     222        3903 :       c[j] = k; d[k] = j;
     223        3903 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
     224      106069 :       for (t=1; t<=m; t++)
     225             :       {
     226      102166 :         if (c[t]) continue; /* already a pivot on that line */
     227             : 
     228       55335 :         piv = ff->red(E,gcoeff(x,t,k));
     229       55335 :         if (ff->equal0(piv)) continue;
     230       53186 :         gcoeff(x,t,k) = g0;
     231     1310736 :         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       53186 :         if (gc_needed(av,1))
     234          89 :           gerepile_gauss(x,k,t,av,j,c);
     235             :       }
     236        3903 :       for (i=k; i<=n; i++) gcoeff(x,j,i) = g0; /* dummy */
     237             :     }
     238             :   }
     239         381 :   *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      107693 : _gen_addmul(GEN b, long k, long i, GEN m, void *E, const struct bb_field *ff)
     290             : {
     291      107693 :   gel(b,i) = ff->red(E,gel(b,i));
     292      107693 :   gel(b,k) = ff->add(E,gel(b,k), ff->mul(E,m, gel(b,i)));
     293      107693 : }
     294             : 
     295             : static GEN
     296       13492 : _gen_get_col(GEN a, GEN b, long li, void *E, const struct bb_field *ff)
     297             : {
     298       13492 :   GEN u = cgetg(li+1,t_COL);
     299       13492 :   pari_sp av = avma;
     300             :   long i, j;
     301             : 
     302       13492 :   gel(u,li) = gerepileupto(av, ff->red(E,ff->mul(E,gel(b,li), gcoeff(a,li,li))));
     303       60221 :   for (i=li-1; i>0; i--)
     304             :   {
     305       46729 :     pari_sp av = avma;
     306       46729 :     GEN m = gel(b,i);
     307       46729 :     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       46729 :     m = ff->red(E, m);
     309       46729 :     gel(u,i) = gerepileupto(av, ff->red(E,ff->mul(E,m, gcoeff(a,i,i))));
     310             :   }
     311       13492 :   return u;
     312             : }
     313             : 
     314             : GEN
     315        4235 : gen_Gauss(GEN a, GEN b, void *E, const struct bb_field *ff)
     316             : {
     317             :   long i, j, k, li, bco, aco;
     318        4235 :   GEN u, g0 = ff->s(E,0);
     319        4235 :   pari_sp av = avma;
     320        4235 :   a = RgM_shallowcopy(a);
     321        4235 :   b = RgM_shallowcopy(b);
     322        4235 :   aco = lg(a)-1; bco = lg(b)-1; li = nbrows(a);
     323       13611 :   for (i=1; i<=aco; i++)
     324             :   {
     325             :     GEN invpiv;
     326       14534 :     for (k = i; k <= li; k++)
     327             :     {
     328       14478 :       GEN piv = ff->red(E,gcoeff(a,k,i));
     329       14478 :       if (!ff->equal0(piv)) { gcoeff(a,k,i) = ff->inv(E,piv); break; }
     330         923 :       gcoeff(a,k,i) = g0;
     331             :     }
     332             :     /* found a pivot on line k */
     333       13611 :     if (k > li) return NULL;
     334       13555 :     if (k != i)
     335             :     { /* swap lines so that k = i */
     336         740 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     337         740 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     338             :     }
     339       13555 :     if (i == aco) break;
     340             : 
     341        9376 :     invpiv = gcoeff(a,i,i); /* 1/piv mod p */
     342       32849 :     for (k=i+1; k<=li; k++)
     343             :     {
     344       23473 :       GEN m = ff->red(E,gcoeff(a,k,i)); gcoeff(a,k,i) = g0;
     345       23473 :       if (ff->equal0(m)) continue;
     346             : 
     347       10907 :       m = ff->red(E,ff->neg(E,ff->mul(E,m, invpiv)));
     348       10907 :       for (j=i+1; j<=aco; j++) _gen_addmul(gel(a,j),k,i,m,E,ff);
     349       10907 :       for (j=1  ; j<=bco; j++) _gen_addmul(gel(b,j),k,i,m,E,ff);
     350             :     }
     351        9376 :     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        4179 :   if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
     359        4179 :   u = cgetg(bco+1,t_MAT);
     360        4179 :   for (j=1; j<=bco; j++) gel(u,j) = _gen_get_col(a, gel(b,j), aco, E, ff);
     361        4179 :   return u;
     362             : }
     363             : 
     364             : /* compatible t_MAT * t_COL, lgA = lg(A) = lg(B) > 1, l = lgcols(A) */
     365             : static GEN
     366       38599 : gen_matcolmul_i(GEN A, GEN B, ulong lgA, ulong l,
     367             :                 void *E, const struct bb_field *ff)
     368             : {
     369       38599 :   GEN C = cgetg(l, t_COL);
     370             :   ulong i;
     371      279626 :   for (i = 1; i < l; i++) {
     372      241027 :     pari_sp av = avma;
     373      241027 :     GEN e = ff->mul(E, gcoeff(A, i, 1), gel(B, 1));
     374             :     ulong k;
     375     1653367 :     for(k = 2; k < lgA; k++)
     376     1412340 :       e = ff->add(E, e, ff->mul(E, gcoeff(A, i, k), gel(B, k)));
     377      241027 :     gel(C, i) = gerepileupto(av, ff->red(E, e));
     378             :   }
     379       38599 :   return C;
     380             : }
     381             : 
     382             : GEN
     383       20986 : gen_matcolmul(GEN A, GEN B, void *E, const struct bb_field *ff)
     384             : {
     385       20986 :   ulong lgA = lg(A);
     386       20986 :   if (lgA != (ulong)lg(B))
     387           0 :     pari_err_OP("operation 'gen_matcolmul'", A, B);
     388       20986 :   if (lgA == 1)
     389           0 :     return cgetg(1, t_COL);
     390       20986 :   return gen_matcolmul_i(A, B, lgA, lgcols(A), E, ff);
     391             : }
     392             : 
     393             : GEN
     394        1496 : gen_matmul(GEN A, GEN B, void *E, const struct bb_field *ff)
     395             : {
     396        1496 :   ulong j, l, lgA, lgB = lg(B);
     397             :   GEN C;
     398        1496 :   if (lgB == 1)
     399           0 :     return cgetg(1, t_MAT);
     400        1496 :   lgA = lg(A);
     401        1496 :   if (lgA != (ulong)lgcols(B))
     402           0 :     pari_err_OP("operation 'gen_matmul'", A, B);
     403        1496 :   if (lgA == 1)
     404           0 :     return zeromat(0, lgB - 1);
     405        1496 :   l = lgcols(A);
     406        1496 :   C = cgetg(lgB, t_MAT);
     407       19109 :   for(j = 1; j < lgB; j++)
     408       17613 :     gel(C, j) = gen_matcolmul_i(A, gel(B, j), lgA, l, E, ff);
     409        1496 :   return C;
     410             : }
     411             : 
     412             : static GEN
     413          63 : gen_colneg(GEN A, void *E, const struct bb_field *ff)
     414             : {
     415             :   long i, l;
     416          63 :   GEN B = cgetg_copy(A, &l);
     417         252 :   for (i = 1; i < l; i++)
     418         189 :     gel(B, i) = ff->neg(E, gel(A, i));
     419          63 :   return B;
     420             : }
     421             : 
     422             : static GEN
     423          21 : gen_matneg(GEN A, void *E, const struct bb_field *ff)
     424             : {
     425             :   long i, l;
     426          21 :   GEN B = cgetg_copy(A, &l);
     427          84 :   for (i = 1; i < l; i++)
     428          63 :     gel(B, i) = gen_colneg(gel(A, i), E, ff);
     429          21 :   return B;
     430             : }
     431             : 
     432             : /* assume dim A >= 1, A invertible + upper triangular  */
     433             : static GEN
     434          63 : gen_matinv_upper_ind(GEN A, long index, void *E, const struct bb_field *ff)
     435             : {
     436          63 :   long n = lg(A) - 1, i, j;
     437          63 :   GEN u = cgetg(n + 1, t_COL);
     438         126 :   for (i = n; i > index; i--)
     439          63 :     gel(u, i) = ff->s(E, 0);
     440          63 :   gel(u, i) = ff->inv(E, gcoeff(A, i, i));
     441         126 :   for (i--; i > 0; i--) {
     442          63 :     pari_sp av = avma;
     443          63 :     GEN m = ff->neg(E, ff->mul(E, gcoeff(A, i, i + 1), gel(u, i + 1)));
     444         105 :     for (j = i + 2; j <= n; j++)
     445          42 :       m = ff->add(E, m, ff->neg(E, ff->mul(E, gcoeff(A, i, j), gel(u, j))));
     446          63 :     gel(u, i) = gerepileupto(av, ff->red(E, ff->mul(E, m, ff->inv(E, gcoeff(A, i, i)))));
     447             :   }
     448          63 :   return u;
     449             : }
     450             : 
     451             : static GEN
     452          21 : gen_matinv_upper(GEN A, void *E, const struct bb_field *ff)
     453             : {
     454             :   long i, l;
     455          21 :   GEN B = cgetg_copy(A, &l);
     456          84 :   for (i = 1; i < l; i++)
     457          63 :     gel(B,i) = gen_matinv_upper_ind(A, i, E, ff);
     458          21 :   return B;
     459             : }
     460             : 
     461             : /* find z such that A z = y. Return NULL if no solution */
     462             : GEN
     463          63 : gen_matcolinvimage(GEN A, GEN y, void *E, const struct bb_field *ff)
     464             : {
     465          63 :   pari_sp av = avma;
     466          63 :   long i, l = lg(A);
     467             :   GEN M, x, t;
     468             : 
     469          63 :   M = gen_ker(shallowconcat(A, y), 0, E, ff);
     470          63 :   i = lg(M) - 1;
     471          63 :   if (!i) { avma = av; return NULL; }
     472             : 
     473          63 :   x = gel(M, i);
     474          63 :   t = gel(x, l);
     475          63 :   if (ff->equal0(t)) { avma = av; return NULL; }
     476             : 
     477          42 :   t = ff->neg(E, ff->inv(E, t));
     478          42 :   setlg(x, l);
     479         147 :   for (i = 1; i < l; i++)
     480         105 :     gel(x, i) = ff->red(E, ff->mul(E, t, gel(x, i)));
     481          42 :   return gerepilecopy(av, x);
     482             : }
     483             : 
     484             : /* find Z such that A Z = B. Return NULL if no solution */
     485             : GEN
     486          21 : gen_matinvimage(GEN A, GEN B, void *E, const struct bb_field *ff)
     487             : {
     488          21 :   pari_sp av = avma;
     489             :   GEN d, x, X, Y;
     490             :   long i, j, nY, nA, nB;
     491          21 :   x = gen_ker(shallowconcat(gen_matneg(A, E, ff), B), 0, E, ff);
     492             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
     493             :    * We must find T such that Y T = Id_nB then X T = Z. This exists
     494             :    * iff Y has at least nB columns and full rank. */
     495          21 :   nY = lg(x) - 1;
     496          21 :   nB = lg(B) - 1;
     497          21 :   if (nY < nB) { avma = av; return NULL; }
     498          21 :   nA = lg(A) - 1;
     499          21 :   Y = rowslice(x, nA + 1, nA + nB); /* nB rows */
     500          21 :   d = cgetg(nB + 1, t_VECSMALL);
     501          84 :   for (i = nB, j = nY; i >= 1; i--, j--) {
     502          63 :     for (; j >= 1; j--)
     503          63 :       if (!ff->equal0(gcoeff(Y, i, j))) { d[i] = j; break; }
     504          63 :     if (!j) { avma = av; return NULL; }
     505             :   }
     506             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
     507          21 :   Y = vecpermute(Y, d);
     508          21 :   x = vecpermute(x, d);
     509          21 :   X = rowslice(x, 1, nA);
     510          21 :   return gerepileupto(av, gen_matmul(X, gen_matinv_upper(Y, E, ff), E, ff));
     511             : }
     512             : 
     513             : static GEN
     514       80679 : image_from_pivot(GEN x, GEN d, long r)
     515             : {
     516             :   GEN y;
     517             :   long j, k;
     518             : 
     519       80679 :   if (!d) return gcopy(x);
     520             :   /* d left on stack for efficiency */
     521       79251 :   r = lg(x)-1 - r; /* = dim Im(x) */
     522       79251 :   y = cgetg(r+1,t_MAT);
     523      695557 :   for (j=k=1; j<=r; k++)
     524      616306 :     if (d[k]) gel(y,j++) = gcopy(gel(x,k));
     525       79251 :   return y;
     526             : }
     527             : 
     528             : /*******************************************************************/
     529             : /*                                                                 */
     530             : /*                    LINEAR ALGEBRA MODULO P                      */
     531             : /*                                                                 */
     532             : /*******************************************************************/
     533             : 
     534             : static long
     535     1928994 : F2v_find_nonzero(GEN x0, GEN mask0, long l, long m)
     536             : {
     537     1928994 :   ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
     538             :   long i, j;
     539     3489724 :   for (i = 0; i < l; i++)
     540             :   {
     541     3488160 :     e = *x++ & *mask++;
     542     3488160 :     if (e)
     543     1927430 :       for (j = 1; ; j++, e >>= 1) if (e & 1uL) return i*BITS_IN_LONG+j;
     544     8049104 :   }
     545        1564 :   return m+1;
     546             : }
     547             : 
     548             : /* in place, destroy x */
     549             : GEN
     550      287803 : F2m_ker_sp(GEN x, long deplin)
     551             : {
     552             :   GEN y, c, d;
     553             :   long i, j, k, l, r, m, n;
     554             : 
     555      287803 :   n = lg(x)-1;
     556      287803 :   m = mael(x,1,1); r=0;
     557             : 
     558      287803 :   d = cgetg(n+1, t_VECSMALL);
     559      287803 :   c = const_F2v(m); l = lg(c)-1;
     560     2070518 :   for (k=1; k<=n; k++)
     561             :   {
     562     1801132 :     GEN xk = gel(x,k);
     563     1801132 :     j = F2v_find_nonzero(xk, c, l, m);
     564     1801132 :     if (j>m)
     565             :     {
     566      393745 :       if (deplin) {
     567       18417 :         GEN c = zero_F2v(n);
     568       65175 :         for (i=1; i<k; i++)
     569       46758 :           if (F2v_coeff(xk, d[i]))
     570       23842 :             F2v_set(c, i);
     571       18417 :         F2v_set(c, k);
     572       18417 :         return c;
     573             :       }
     574      375328 :       r++; d[k] = 0;
     575             :     }
     576             :     else
     577             :     {
     578     1407387 :       F2v_clear(c,j); d[k] = j;
     579     1407387 :       F2v_clear(xk, j);
     580    75123868 :       for (i=k+1; i<=n; i++)
     581             :       {
     582    73716481 :         GEN xi = gel(x,i);
     583    73716481 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     584             :       }
     585     1407387 :       F2v_set(xk, j);
     586             :     }
     587             :   }
     588      269386 :   if (deplin) return NULL;
     589             : 
     590      269351 :   y = zero_F2m_copy(n,r);
     591      644679 :   for (j=k=1; j<=r; j++,k++)
     592             :   {
     593      375328 :     GEN C = gel(y,j); while (d[k]) k++;
     594    13376204 :     for (i=1; i<k; i++)
     595    13000876 :       if (d[i] && F2m_coeff(x,d[i],k))
     596     4798805 :         F2v_set(C,i);
     597      375328 :     F2v_set(C, k);
     598             :   }
     599      269351 :   return y;
     600             : }
     601             : 
     602             : static void /* assume m < p */
     603   143588176 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
     604             : {
     605   143588176 :   uel(b,k) = Fl_addmul_pre(m, uel(b,i), uel(b,k), p, pi);
     606   143588176 : }
     607             : static void /* same m = 1 */
     608     3168853 : _Fl_add(GEN b, long k, long i, ulong p)
     609             : {
     610     3168853 :   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
     611     3168853 : }
     612             : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     613   250837302 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
     614             : {
     615   250837302 :   uel(b,k) += m * uel(b,i);
     616   250837302 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     617   250837302 : }
     618             : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     619    78842284 : _Fl_add_OK(GEN b, long k, long i, ulong p)
     620             : {
     621    78842284 :   uel(b,k) += uel(b,i);
     622    78842284 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     623    78842284 : }
     624             : 
     625             : static GEN
     626      487870 : Flm_ker_sp_OK(GEN x, ulong p, long deplin)
     627             : {
     628             :   GEN y, c, d;
     629             :   long i, j, k, r, t, m, n;
     630             :   ulong a;
     631             : 
     632      487870 :   n = lg(x)-1;
     633      487870 :   m=nbrows(x); r=0;
     634             : 
     635      487870 :   c = zero_zv(m);
     636      487870 :   d = cgetg(n+1, t_VECSMALL);
     637      487870 :   a = 0; /* for gcc -Wall */
     638     3780893 :   for (k=1; k<=n; k++)
     639             :   {
     640    29477870 :     for (j=1; j<=m; j++)
     641    28286368 :       if (!c[j])
     642             :       {
     643    11129192 :         a = ucoeff(x,j,k) % p;
     644    11129192 :         if (a) break;
     645             :       }
     646     3332661 :     if (j > m)
     647             :     {
     648     1191502 :       if (deplin==1) {
     649       39638 :         c = cgetg(n+1, t_VECSMALL);
     650       39638 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
     651       39638 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     652       39638 :         return c;
     653             :       }
     654     1151864 :       r++; d[k] = 0;
     655             :     }
     656             :     else
     657             :     {
     658     2141159 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     659     2141159 :       c[j] = k; d[k] = j;
     660     2141159 :       ucoeff(x,j,k) = p-1;
     661     2141159 :       if (piv != 1)
     662     1718116 :         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
     663    61469582 :       for (t=1; t<=m; t++)
     664             :       {
     665    59328423 :         if (t == j) continue;
     666             : 
     667    57187264 :         piv = ( ucoeff(x,t,k) %= p );
     668    57187264 :         if (!piv) continue;
     669    20296645 :         if (piv == 1)
     670     4242279 :           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
     671             :         else
     672    16054366 :           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
     673             :       }
     674             :     }
     675             :   }
     676      448232 :   if (deplin==1) return NULL;
     677             : 
     678      448225 :   y = cgetg(r+1, t_MAT);
     679     1600089 :   for (j=k=1; j<=r; j++,k++)
     680             :   {
     681     1151864 :     GEN C = cgetg(n+1, t_VECSMALL);
     682             : 
     683     1151864 :     gel(y,j) = C; while (d[k]) k++;
     684     9801869 :     for (i=1; i<k; i++)
     685     8650005 :       if (d[i])
     686     6032117 :         uel(C,i) = ucoeff(x,d[i],k) % p;
     687             :       else
     688     2617888 :         uel(C,i) = 0UL;
     689     1151864 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     690             :   }
     691      448225 :   return deplin==2? mkvec2(y, d): y;
     692             : }
     693             : 
     694             : /* in place, destroy x */
     695             : GEN
     696      495550 : Flm_ker_sp(GEN x, ulong p, long deplin)
     697             : {
     698             :   GEN y, c, d;
     699             :   long i, j, k, r, t, m, n;
     700             :   ulong a, pi;
     701      495550 :   n = lg(x)-1;
     702      495550 :   if (!n) return cgetg(1,t_MAT);
     703      495550 :   if (SMALL_ULONG(p)) return Flm_ker_sp_OK(x, p, deplin);
     704        7680 :   pi = get_Fl_red(p);
     705             : 
     706        7680 :   m=nbrows(x); r=0;
     707             : 
     708        7680 :   c = zero_zv(m);
     709        7680 :   d = cgetg(n+1, t_VECSMALL);
     710        7680 :   a = 0; /* for gcc -Wall */
     711       86158 :   for (k=1; k<=n; k++)
     712             :   {
     713     2101142 :     for (j=1; j<=m; j++)
     714     2068875 :       if (!c[j])
     715             :       {
     716     1043876 :         a = ucoeff(x,j,k);
     717     1043876 :         if (a) break;
     718             :       }
     719       78485 :     if (j > m)
     720             :     {
     721       32267 :       if (deplin==1) {
     722           7 :         c = cgetg(n+1, t_VECSMALL);
     723           7 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k);
     724           7 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     725           7 :         return c;
     726             :       }
     727       32260 :       r++; d[k] = 0;
     728             :     }
     729             :     else
     730             :     {
     731       46218 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     732       46218 :       c[j] = k; d[k] = j;
     733       46218 :       ucoeff(x,j,k) = p-1;
     734       46218 :       if (piv != 1)
     735      757020 :         for (i=k+1; i<=n; i++)
     736      714363 :           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
     737     1890514 :       for (t=1; t<=m; t++)
     738             :       {
     739     1844296 :         if (t == j) continue;
     740             : 
     741     1798078 :         piv = ucoeff(x,t,k);
     742     1798078 :         if (!piv) continue;
     743      628994 :         if (piv == 1)
     744        8033 :           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
     745             :         else
     746      620961 :           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
     747             :       }
     748             :     }
     749             :   }
     750        7673 :   if (deplin==1) return NULL;
     751             : 
     752        7666 :   y = cgetg(r+1, t_MAT);
     753       39926 :   for (j=k=1; j<=r; j++,k++)
     754             :   {
     755       32260 :     GEN C = cgetg(n+1, t_VECSMALL);
     756             : 
     757       32260 :     gel(y,j) = C; while (d[k]) k++;
     758      965490 :     for (i=1; i<k; i++)
     759      933230 :       if (d[i])
     760      588991 :         uel(C,i) = ucoeff(x,d[i],k);
     761             :       else
     762      344239 :         uel(C,i) = 0UL;
     763       32260 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     764             :   }
     765        7666 :   return deplin==2? mkvec2(y, d): y;
     766             : }
     767             : 
     768             : GEN
     769       68586 : FpM_intersect(GEN x, GEN y, GEN p)
     770             : {
     771       68586 :   pari_sp av = avma;
     772       68586 :   long j, lx = lg(x);
     773             :   GEN z;
     774             : 
     775       68586 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     776       68586 :   z = FpM_ker(shallowconcat(x,y), p);
     777       68586 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     778       68586 :   return gerepileupto(av, FpM_mul(x,z,p));
     779             : }
     780             : GEN
     781           0 : Flm_intersect(GEN x, GEN y, ulong p)
     782             : {
     783           0 :   pari_sp av = avma;
     784           0 :   long j, lx = lg(x);
     785             :   GEN z;
     786             : 
     787           0 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     788           0 :   z = Flm_ker(shallowconcat(x,y), p);
     789           0 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     790           0 :   return gerepileupto(av, Flm_mul(x,z,p));
     791             : }
     792             : 
     793             : /* not memory clean */
     794             : GEN
     795        1551 : F2m_ker(GEN x) { return F2m_ker_sp(F2m_copy(x), 0); }
     796             : GEN
     797           0 : F2m_deplin(GEN x) { return F2m_ker_sp(F2m_copy(x), 1); }
     798             : GEN
     799       54677 : Flm_ker(GEN x, ulong p) { return Flm_ker_sp(Flm_copy(x), p, 0); }
     800             : GEN
     801           0 : Flm_deplin(GEN x, ulong p) { return Flm_ker_sp(Flm_copy(x), p, 1); }
     802             : 
     803             : ulong
     804          28 : F2m_det_sp(GEN x) { return !F2m_ker_sp(x, 1); }
     805             : 
     806             : ulong
     807          14 : F2m_det(GEN x)
     808             : {
     809          14 :   pari_sp av = avma;
     810          14 :   ulong d = F2m_det_sp(F2m_copy(x));
     811          14 :   avma = av; return d;
     812             : }
     813             : 
     814             : /* in place, destroy a, SMALL_ULONG(p) is TRUE */
     815             : static ulong
     816        1407 : Flm_det_sp_OK(GEN a, long nbco, ulong p)
     817             : {
     818        1407 :   long i,j,k, s = 1;
     819        1407 :   ulong q, x = 1;
     820             : 
     821       11123 :   for (i=1; i<nbco; i++)
     822             :   {
     823        9898 :     for(k=i; k<=nbco; k++)
     824             :     {
     825        9877 :       ulong c = ucoeff(a,k,i) % p;
     826        9877 :       ucoeff(a,k,i) = c;
     827        9877 :       if (c) break;
     828             :     }
     829        9737 :     for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
     830        9737 :     if (k > nbco) return ucoeff(a,i,i);
     831        9716 :     if (k != i)
     832             :     { /* exchange the lines s.t. k = i */
     833          77 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     834          77 :       s = -s;
     835             :     }
     836        9716 :     q = ucoeff(a,i,i);
     837             : 
     838        9716 :     if (x & HIGHMASK) x %= p;
     839        9716 :     x *= q;
     840        9716 :     q = Fl_inv(q,p);
     841       48587 :     for (k=i+1; k<=nbco; k++)
     842             :     {
     843       38871 :       ulong m = ucoeff(a,i,k) % p;
     844       38871 :       if (!m) continue;
     845             : 
     846       34265 :       m = p - ((m*q)%p);
     847      197379 :       for (j=i+1; j<=nbco; j++)
     848             :       {
     849      163114 :         ulong c = ucoeff(a,j,k);
     850      163114 :         if (c & HIGHMASK) c %= p;
     851      163114 :         ucoeff(a,j,k) = c  + m*ucoeff(a,j,i);
     852             :       }
     853             :     }
     854             :   }
     855        1386 :   if (x & HIGHMASK) x %= p;
     856        1386 :   q = ucoeff(a,nbco,nbco);
     857        1386 :   if (q & HIGHMASK) q %= p;
     858        1386 :   x = (x*q) % p;
     859        1386 :   if (s < 0 && x) x = p - x;
     860        1386 :   return x;
     861             : }
     862             : /* in place, destroy a */
     863             : ulong
     864       69833 : Flm_det_sp(GEN a, ulong p)
     865             : {
     866       69833 :   long i,j,k, s = 1, nbco = lg(a)-1;
     867       69833 :   ulong pi, q, x = 1;
     868             : 
     869       69833 :   if (SMALL_ULONG(p)) return Flm_det_sp_OK(a, nbco, p);
     870       68426 :   pi = get_Fl_red(p);
     871      515029 :   for (i=1; i<nbco; i++)
     872             :   {
     873      447969 :     for(k=i; k<=nbco; k++)
     874      447969 :       if (ucoeff(a,k,i)) break;
     875      446603 :     if (k > nbco) return ucoeff(a,i,i);
     876      446603 :     if (k != i)
     877             :     { /* exchange the lines s.t. k = i */
     878         295 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     879         295 :       s = -s;
     880             :     }
     881      446603 :     q = ucoeff(a,i,i);
     882             : 
     883      446603 :     x = Fl_mul_pre(x, q, p, pi);
     884      446603 :     q = Fl_inv(q,p);
     885     2210399 :     for (k=i+1; k<=nbco; k++)
     886             :     {
     887     1763796 :       ulong m = ucoeff(a,i,k);
     888     1763796 :       if (!m) continue;
     889             : 
     890     1634297 :       m = Fl_mul_pre(m, q, p, pi);
     891     9807539 :       for (j=i+1; j<=nbco; j++)
     892     8173242 :         ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
     893             :     }
     894             :   }
     895       68426 :   if (s < 0) x = Fl_neg(x, p);
     896       68426 :   return Fl_mul(x, ucoeff(a,nbco,nbco), p);
     897             : }
     898             : 
     899             : ulong
     900       69826 : Flm_det(GEN x, ulong p)
     901             : {
     902       69826 :   pari_sp av = avma;
     903       69826 :   ulong d = Flm_det_sp(Flm_copy(x), p);
     904       69826 :   avma = av; return d;
     905             : }
     906             : 
     907             : static GEN
     908      384143 : FpM_init(GEN a, GEN p, ulong *pp)
     909             : {
     910      384143 :   if (lgefint(p) == 3)
     911             :   {
     912      379609 :     *pp = uel(p,2);
     913      379609 :     return (*pp==2)? ZM_to_F2m(a): ZM_to_Flm(a, *pp);
     914             :   }
     915        4534 :   *pp = 0; return a;
     916             : }
     917             : GEN
     918        2128 : RgM_Fp_init(GEN a, GEN p, ulong *pp)
     919             : {
     920        2128 :   if (lgefint(p) == 3)
     921             :   {
     922        1923 :     *pp = uel(p,2);
     923        1923 :     return (*pp==2)? RgM_to_F2m(a): RgM_to_Flm(a, *pp);
     924             :   }
     925         205 :   *pp = 0; return RgM_to_FpM(a,p);
     926             : }
     927             : 
     928             : static GEN
     929          49 : FpM_det_gen(GEN a, GEN p)
     930             : {
     931             :   void *E;
     932          49 :   const struct bb_field *S = get_Fp_field(&E,p);
     933          49 :   return gen_det(a, E, S);
     934             : }
     935             : GEN
     936          70 : FpM_det(GEN a, GEN p)
     937             : {
     938          70 :   pari_sp av = avma;
     939             :   ulong pp, d;
     940          70 :   a = FpM_init(a, p, &pp);
     941          70 :   switch(pp)
     942             :   {
     943          49 :   case 0: return FpM_det_gen(a, p);
     944          14 :   case 2: d = F2m_det_sp(a); break;
     945           7 :   default:d = Flm_det_sp(a,pp); break;
     946             :   }
     947          21 :   avma = av; return utoi(d);
     948             : }
     949             : 
     950             : /* Destroy x */
     951             : static GEN
     952       23126 : F2m_gauss_pivot(GEN x, long *rr)
     953             : {
     954             :   GEN c, d;
     955             :   long i, j, k, l, r, m, n;
     956             : 
     957       23126 :   n = lg(x)-1; if (!n) { *rr=0; return NULL; }
     958       23126 :   m = mael(x,1,1); r=0;
     959             : 
     960       23126 :   d = cgetg(n+1, t_VECSMALL);
     961       23126 :   c = const_F2v(m); l = lg(c)-1;
     962      150988 :   for (k=1; k<=n; k++)
     963             :   {
     964      127862 :     GEN xk = gel(x,k);
     965      127862 :     j = F2v_find_nonzero(xk, c, l, m);
     966      127862 :     if (j>m) { r++; d[k] = 0; }
     967             :     else
     968             :     {
     969       82800 :       F2v_clear(c,j); d[k] = j;
     970     1419263 :       for (i=k+1; i<=n; i++)
     971             :       {
     972     1336463 :         GEN xi = gel(x,i);
     973     1336463 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     974             :       }
     975             :     }
     976             :   }
     977             : 
     978       23126 :   *rr = r; avma = (pari_sp)d; return d;
     979             : }
     980             : 
     981             : /* Destroy x */
     982             : static GEN
     983      199237 : Flm_gauss_pivot(GEN x, ulong p, long *rr)
     984             : {
     985             :   GEN c,d;
     986             :   long i,j,k,r,t,n,m;
     987             : 
     988      199237 :   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
     989             : 
     990      199237 :   m=nbrows(x); r=0;
     991      199237 :   d=cgetg(n+1,t_VECSMALL);
     992      199237 :   c = zero_zv(m);
     993     1784519 :   for (k=1; k<=n; k++)
     994             :   {
     995    16744928 :     for (j=1; j<=m; j++)
     996    16215380 :       if (!c[j])
     997             :       {
     998     4834242 :         ucoeff(x,j,k) %= p;
     999     4834242 :         if (ucoeff(x,j,k)) break;
    1000             :       }
    1001     1585282 :     if (j>m) { r++; d[k]=0; }
    1002             :     else
    1003             :     {
    1004     1055734 :       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
    1005     1055734 :       c[j]=k; d[k]=j;
    1006    13188941 :       for (i=k+1; i<=n; i++)
    1007    12133207 :         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
    1008    19359172 :       for (t=1; t<=m; t++)
    1009    18303438 :         if (!c[t]) /* no pivot on that line yet */
    1010             :         {
    1011    10832607 :           piv = ucoeff(x,t,k);
    1012    10832607 :           if (piv)
    1013             :           {
    1014     5317755 :             ucoeff(x,t,k) = 0;
    1015   223698549 :             for (i=k+1; i<=n; i++)
    1016   436761588 :               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
    1017   218380794 :                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
    1018             :           }
    1019             :         }
    1020     1055734 :       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
    1021             :     }
    1022             :   }
    1023      199237 :   *rr = r; avma = (pari_sp)d; return d;
    1024             : }
    1025             : 
    1026             : static GEN
    1027          59 : FpM_gauss_pivot_gen(GEN x, GEN p, long *rr)
    1028             : {
    1029             :   void *E;
    1030          59 :   const struct bb_field *S = get_Fp_field(&E,p);
    1031          59 :   return gen_Gauss_pivot(x, rr, E, S);
    1032             : }
    1033             : static GEN
    1034      108543 : FpM_gauss_pivot(GEN x, GEN p, long *rr)
    1035             : {
    1036             :   ulong pp;
    1037      108543 :   if (lg(x)==1) { *rr = 0; return NULL; }
    1038      107150 :   x = FpM_init(x, p, &pp);
    1039      107150 :   switch(pp)
    1040             :   {
    1041          59 :   case 0: return FpM_gauss_pivot_gen(x, p, rr);
    1042       23042 :   case 2: return F2m_gauss_pivot(x, rr);
    1043       84049 :   default:return Flm_gauss_pivot(x, pp, rr);
    1044             :   }
    1045             : }
    1046             : 
    1047             : GEN
    1048       77732 : FpM_image(GEN x, GEN p)
    1049             : {
    1050             :   long r;
    1051       77732 :   GEN d = FpM_gauss_pivot(x,p,&r); /* d left on stack for efficiency */
    1052       77732 :   return image_from_pivot(x,d,r);
    1053             : }
    1054             : GEN
    1055          49 : Flm_image(GEN x, ulong p)
    1056             : {
    1057             :   long r;
    1058          49 :   GEN d = Flm_gauss_pivot(Flm_copy(x),p,&r); /* d left on stack for efficiency */
    1059          49 :   return image_from_pivot(x,d,r);
    1060             : }
    1061             : GEN
    1062           7 : F2m_image(GEN x)
    1063             : {
    1064             :   long r;
    1065           7 :   GEN d = F2m_gauss_pivot(F2m_copy(x),&r); /* d left on stack for efficiency */
    1066           7 :   return image_from_pivot(x,d,r);
    1067             : }
    1068             : 
    1069             : long
    1070           7 : FpM_rank(GEN x, GEN p)
    1071             : {
    1072           7 :   pari_sp av = avma;
    1073             :   long r;
    1074           7 :   (void)FpM_gauss_pivot(x,p,&r);
    1075           7 :   avma = av; return lg(x)-1 - r;
    1076             : }
    1077             : long
    1078        2188 : Flm_rank(GEN x, ulong p)
    1079             : {
    1080        2188 :   pari_sp av = avma;
    1081             :   long r;
    1082        2188 :   (void)Flm_gauss_pivot(Flm_copy(x),p,&r);
    1083        2188 :   avma = av; return lg(x)-1 - r;
    1084             : }
    1085             : long
    1086          63 : F2m_rank(GEN x)
    1087             : {
    1088          63 :   pari_sp av = avma;
    1089             :   long r;
    1090          63 :   (void)F2m_gauss_pivot(F2m_copy(x),&r);
    1091          63 :   avma = av; return lg(x)-1 - r;
    1092             : }
    1093             : 
    1094             : 
    1095             : static GEN
    1096         244 : FlxqM_gauss_pivot(GEN x, GEN T, ulong p, long *rr)
    1097             : {
    1098             :   void *E;
    1099         244 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1100         244 :   return gen_Gauss_pivot(x, rr, E, S);
    1101             : }
    1102             : 
    1103             : GEN
    1104           7 : FlxqM_image(GEN x, GEN T, ulong p)
    1105             : {
    1106             :   long r;
    1107           7 :   GEN d = FlxqM_gauss_pivot(x,T,p,&r); /* d left on stack for efficiency */
    1108           7 :   return image_from_pivot(x,d,r);
    1109             : }
    1110             : 
    1111             : long
    1112         139 : FlxqM_rank(GEN x, GEN T, ulong p)
    1113             : {
    1114         139 :   pari_sp av = avma;
    1115             :   long r;
    1116         139 :   (void)FlxqM_gauss_pivot(x,T,p,&r);
    1117         139 :   avma = av; return lg(x)-1 - r;
    1118             : }
    1119             : GEN
    1120           7 : FlxqM_det(GEN a, GEN T, ulong p)
    1121             : {
    1122             :   void *E;
    1123           7 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1124           7 :   return gen_det(a, E, S);
    1125             : }
    1126             : 
    1127             : GEN
    1128          21 : FlxqM_FlxqC_invimage(GEN A, GEN B, GEN T, ulong p) {
    1129             :   void *E;
    1130          21 :   const struct bb_field *ff = get_Flxq_field(&E, T, p);
    1131          21 :   return gen_matcolinvimage(A, B, E, ff);
    1132             : }
    1133             : 
    1134             : GEN
    1135          21 : FlxqM_FlxqC_mul(GEN A, GEN B, GEN T, ulong p) {
    1136             :   void *E;
    1137          21 :   const struct bb_field *ff = get_Flxq_field(&E, T, p);
    1138          21 :   return gen_matcolmul(A, B, E, ff);
    1139             : }
    1140             : 
    1141             : GEN
    1142         142 : FlxqM_mul(GEN A, GEN B, GEN T, ulong p) {
    1143             :   void *E;
    1144             :   const struct bb_field *ff;
    1145         142 :   long n = lg(A) - 1;
    1146             : 
    1147         142 :   if (n == 0)
    1148           0 :     return cgetg(1, t_MAT);
    1149         142 :   if (n > 1)
    1150         142 :     return FlxqM_mul_Kronecker(A, B, T, p);
    1151           0 :   ff = get_Flxq_field(&E, T, p);
    1152           0 :   return gen_matmul(A, B, E, ff);
    1153             : }
    1154             : 
    1155             : GEN
    1156           7 : FlxqM_invimage(GEN A, GEN B, GEN T, ulong p) {
    1157             :   void *E;
    1158           7 :   const struct bb_field *ff = get_Flxq_field(&E, T, p);
    1159           7 :   return gen_matinvimage(A, B, E, ff);
    1160             : }
    1161             : 
    1162             : GEN
    1163           7 : F2xqM_det(GEN a, GEN T)
    1164             : {
    1165             :   void *E;
    1166           7 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1167           7 :   return gen_det(a, E, S);
    1168             : }
    1169             : 
    1170             : static GEN
    1171          77 : F2xqM_gauss_gen(GEN a, GEN b, GEN T)
    1172             : {
    1173             :   void *E;
    1174          77 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1175          77 :   return gen_Gauss(a, b, E, S);
    1176             : }
    1177             : 
    1178             : GEN
    1179          14 : F2xqM_gauss(GEN a, GEN b, GEN T)
    1180             : {
    1181          14 :   pari_sp av = avma;
    1182          14 :   long n = lg(a)-1;
    1183             :   GEN u;
    1184          14 :   if (!n || lg(b)==1) { avma = av; return cgetg(1, t_MAT); }
    1185          14 :   u = F2xqM_gauss_gen(a, b, T);
    1186          14 :   if (!u) { avma = av; return NULL; }
    1187           7 :   return gerepilecopy(av, u);
    1188             : }
    1189             : 
    1190             : GEN
    1191          35 : F2xqM_inv(GEN a, GEN T)
    1192             : {
    1193          35 :   pari_sp av = avma;
    1194             :   GEN u;
    1195          35 :   if (lg(a) == 1) { avma = av; return cgetg(1, t_MAT); }
    1196          35 :   u = F2xqM_gauss_gen(a, matid_F2xqM(nbrows(a),T), T);
    1197          35 :   if (!u) { avma = av; return NULL; }
    1198          28 :   return gerepilecopy(av, u);
    1199             : }
    1200             : 
    1201             : GEN
    1202          28 : F2xqM_F2xqC_gauss(GEN a, GEN b, GEN T)
    1203             : {
    1204          28 :   pari_sp av = avma;
    1205             :   GEN u;
    1206          28 :   if (lg(a) == 1) return cgetg(1, t_COL);
    1207          28 :   u = F2xqM_gauss_gen(a, mkmat(b), T);
    1208          28 :   if (!u) { avma = av; return NULL; }
    1209          14 :   return gerepilecopy(av, gel(u,1));
    1210             : }
    1211             : 
    1212             : GEN
    1213          21 : F2xqM_F2xqC_invimage(GEN A, GEN B, GEN T) {
    1214             :   void *E;
    1215          21 :   const struct bb_field *ff = get_F2xq_field(&E, T);
    1216          21 :   return gen_matcolinvimage(A, B, E, ff);
    1217             : }
    1218             : 
    1219             : GEN
    1220          21 : F2xqM_F2xqC_mul(GEN A, GEN B, GEN T) {
    1221             :   void *E;
    1222          21 :   const struct bb_field *ff = get_F2xq_field(&E, T);
    1223          21 :   return gen_matcolmul(A, B, E, ff);
    1224             : }
    1225             : 
    1226             : GEN
    1227           7 : F2xqM_mul(GEN A, GEN B, GEN T) {
    1228             :   void *E;
    1229           7 :   const struct bb_field *ff = get_F2xq_field(&E, T);
    1230           7 :   return gen_matmul(A, B, E, ff);
    1231             : }
    1232             : 
    1233             : GEN
    1234           7 : F2xqM_invimage(GEN A, GEN B, GEN T) {
    1235             :   void *E;
    1236           7 :   const struct bb_field *ff = get_F2xq_field(&E, T);
    1237           7 :   return gen_matinvimage(A, B, E, ff);
    1238             : }
    1239             : 
    1240             : static GEN
    1241          64 : FqM_gauss_pivot_gen(GEN x, GEN T, GEN p, long *rr)
    1242             : {
    1243             :   void *E;
    1244          64 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1245          64 :   return gen_Gauss_pivot(x, rr, E, S);
    1246             : }
    1247             : static GEN
    1248         162 : FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr)
    1249             : {
    1250         162 :   if (lg(x)==1) { *rr = 0; return NULL; }
    1251         162 :   if (!T) return FpM_gauss_pivot(x, p, rr);
    1252         162 :   if (lgefint(p) == 3)
    1253             :   {
    1254          98 :     pari_sp av = avma;
    1255          98 :     ulong pp = uel(p,2);
    1256          98 :     GEN Tp = ZXT_to_FlxT(T, pp);
    1257          98 :     GEN d = FlxqM_gauss_pivot(FqM_to_FlxM(x, T, p), Tp, pp, rr);
    1258          98 :     return d ? gerepileuptoleaf(av, d): d;
    1259             :   }
    1260          64 :   return FqM_gauss_pivot_gen(x, T, p, rr);
    1261             : }
    1262             : 
    1263             : GEN
    1264           7 : FqM_image(GEN x, GEN T, GEN p)
    1265             : {
    1266             :   long r;
    1267           7 :   GEN d = FqM_gauss_pivot(x,T,p,&r); /* d left on stack for efficiency */
    1268           7 :   return image_from_pivot(x,d,r);
    1269             : }
    1270             : 
    1271             : long
    1272          57 : FqM_rank(GEN x, GEN T, GEN p)
    1273             : {
    1274          57 :   pari_sp av = avma;
    1275             :   long r;
    1276          57 :   (void)FqM_gauss_pivot(x,T,p,&r);
    1277          57 :   avma = av; return lg(x)-1 - r;
    1278             : }
    1279             : 
    1280             : GEN
    1281           7 : FqM_det(GEN x, GEN T, GEN p)
    1282             : {
    1283             :   void *E;
    1284           7 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1285           7 :   return gen_det(x, E, S);
    1286             : }
    1287             : 
    1288             : GEN
    1289          21 : FqM_FqC_invimage(GEN A, GEN B, GEN T, GEN p) {
    1290             :   void *E;
    1291          21 :   const struct bb_field *ff = get_Fq_field(&E, T, p);
    1292          21 :   return gen_matcolinvimage(A, B, E, ff);
    1293             : }
    1294             : 
    1295             : GEN
    1296          21 : FqM_FqC_mul(GEN A, GEN B, GEN T, GEN p) {
    1297             :   void *E;
    1298          21 :   const struct bb_field *ff = get_Fq_field(&E, T, p);
    1299          21 :   return gen_matcolmul(A, B, E, ff);
    1300             : }
    1301             : 
    1302             : GEN
    1303         320 : FqM_mul(GEN A, GEN B, GEN T, GEN p) {
    1304             :   void *E;
    1305         320 :   const struct bb_field *ff = get_Fq_field(&E, T, p);
    1306         320 :   return gen_matmul(A, B, E, ff);
    1307             : }
    1308             : 
    1309             : GEN
    1310           7 : FqM_invimage(GEN A, GEN B, GEN T, GEN p) {
    1311             :   void *E;
    1312           7 :   const struct bb_field *ff = get_Fq_field(&E, T, p);
    1313           7 :   return gen_matinvimage(A, B, E, ff);
    1314             : }
    1315             : 
    1316             : static GEN
    1317        1647 : FpM_ker_gen(GEN x, GEN p, long deplin)
    1318             : {
    1319             :   void *E;
    1320        1647 :   const struct bb_field *S = get_Fp_field(&E,p);
    1321        1647 :   return gen_ker(x, deplin, E, S);
    1322             : }
    1323             : static GEN
    1324      238434 : FpM_ker_i(GEN x, GEN p, long deplin)
    1325             : {
    1326      238434 :   pari_sp av = avma;
    1327             :   ulong pp;
    1328             :   GEN y;
    1329             : 
    1330      238434 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1331      238329 :   x = FpM_init(x, p, &pp);
    1332      238329 :   switch(pp)
    1333             :   {
    1334        1619 :   case 0: return FpM_ker_gen(x,p,deplin);
    1335             :   case 2:
    1336       56224 :     y = F2m_ker_sp(x, deplin);
    1337       56224 :     if (!y) return y;
    1338       56224 :     y = deplin? F2c_to_ZC(y): F2m_to_ZM(y);
    1339       56224 :     return gerepileupto(av, y);
    1340             :   default:
    1341      180486 :     y = Flm_ker_sp(x, pp, deplin);
    1342      180486 :     if (!y) return y;
    1343      180486 :     y = deplin? Flc_to_ZC(y): Flm_to_ZM(y);
    1344      180486 :     return gerepileupto(av, y);
    1345             :   }
    1346             : }
    1347             : 
    1348             : GEN
    1349      179392 : FpM_ker(GEN x, GEN p) { return FpM_ker_i(x,p,0); }
    1350             : 
    1351             : GEN
    1352       58034 : FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x,p,1); }
    1353             : 
    1354             : static GEN
    1355          15 : FqM_ker_gen(GEN x, GEN T, GEN p, long deplin)
    1356             : {
    1357             :   void *E;
    1358          15 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1359          15 :   return gen_ker(x,deplin,E,S);
    1360             : }
    1361             : static GEN
    1362        1254 : FqM_ker_i(GEN x, GEN T, GEN p, long deplin)
    1363             : {
    1364        1254 :   if (!T) return FpM_ker_i(x,p,deplin);
    1365         246 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1366             : 
    1367         246 :   if (lgefint(p)==3)
    1368             :   {
    1369         231 :     pari_sp ltop=avma;
    1370         231 :     ulong l= p[2];
    1371         231 :     GEN Ml = FqM_to_FlxM(x, T, p);
    1372         231 :     GEN Tl = ZXT_to_FlxT(T,l);
    1373         231 :     GEN p1 = FlxM_to_ZXM(FlxqM_ker(Ml,Tl,l));
    1374         231 :     return gerepileupto(ltop,p1);
    1375             :   }
    1376          15 :   return FqM_ker_gen(x, T, p, deplin);
    1377             : }
    1378             : 
    1379             : GEN
    1380        1254 : FqM_ker(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,0); }
    1381             : 
    1382             : GEN
    1383           0 : FqM_deplin(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,1); }
    1384             : 
    1385             : static GEN
    1386         552 : FlxqM_ker_i(GEN x, GEN T, ulong p, long deplin)
    1387             : {
    1388             :   const struct bb_field *ff;
    1389             :   void *E;
    1390             : 
    1391         552 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1392         552 :   ff=get_Flxq_field(&E,T,p);
    1393         552 :   return gen_ker(x,deplin, E, ff);
    1394             : }
    1395             : 
    1396             : GEN
    1397         552 : FlxqM_ker(GEN x, GEN T, ulong p)
    1398             : {
    1399         552 :   return FlxqM_ker_i(x, T, p, 0);
    1400             : }
    1401             : 
    1402             : static GEN
    1403          21 : F2xqM_ker_i(GEN x, GEN T, long deplin)
    1404             : {
    1405             :   const struct bb_field *ff;
    1406             :   void *E;
    1407             : 
    1408          21 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1409          21 :   ff = get_F2xq_field(&E,T);
    1410          21 :   return gen_ker(x,deplin, E, ff);
    1411             : }
    1412             : 
    1413             : GEN
    1414          21 : F2xqM_ker(GEN x, GEN T)
    1415             : {
    1416          21 :   return F2xqM_ker_i(x, T, 0);
    1417             : }
    1418             : static GEN
    1419          14 : F2xqM_gauss_pivot(GEN x, GEN T, long *rr)
    1420             : {
    1421             :   void *E;
    1422          14 :   const struct bb_field *S = get_F2xq_field(&E,T);
    1423          14 :   return gen_Gauss_pivot(x, rr, E, S);
    1424             : }
    1425             : GEN
    1426           7 : F2xqM_image(GEN x, GEN T)
    1427             : {
    1428             :   long r;
    1429           7 :   GEN d = F2xqM_gauss_pivot(x,T,&r); /* d left on stack for efficiency */
    1430           7 :   return image_from_pivot(x,d,r);
    1431             : }
    1432             : long
    1433           7 : F2xqM_rank(GEN x, GEN T)
    1434             : {
    1435           7 :   pari_sp av = avma;
    1436             :   long r;
    1437           7 :   (void)F2xqM_gauss_pivot(x,T,&r);
    1438           7 :   avma = av; return lg(x)-1 - r;
    1439             : }
    1440             : /*******************************************************************/
    1441             : /*                                                                 */
    1442             : /*                       Solve A*X=B (Gauss pivot)                 */
    1443             : /*                                                                 */
    1444             : /*******************************************************************/
    1445             : /* x ~ 0 compared to reference y */
    1446             : int
    1447      657593 : approx_0(GEN x, GEN y)
    1448             : {
    1449      657593 :   long tx = typ(x);
    1450      657593 :   if (tx == t_COMPLEX)
    1451           0 :     return approx_0(gel(x,1), y) && approx_0(gel(x,2), y);
    1452      657733 :   return gequal0(x) ||
    1453      460422 :          (tx == t_REAL && gexpo(y) - gexpo(x) > bit_prec(x));
    1454             : }
    1455             : /* x a column, x0 same column in the original input matrix (for reference),
    1456             :  * c list of pivots so far */
    1457             : static long
    1458      661737 : gauss_get_pivot_max(GEN X, GEN X0, long ix, GEN c)
    1459             : {
    1460      661737 :   GEN p, r, x = gel(X,ix), x0 = gel(X0,ix);
    1461      661737 :   long i, k = 0, ex = - (long)HIGHEXPOBIT, lx = lg(x);
    1462      661737 :   if (c)
    1463             :   {
    1464       46550 :     for (i=1; i<lx; i++)
    1465       35805 :       if (!c[i])
    1466             :       {
    1467       18634 :         long e = gexpo(gel(x,i));
    1468       18634 :         if (e > ex) { ex = e; k = i; }
    1469             :       }
    1470             :   }
    1471             :   else
    1472             :   {
    1473     2243911 :     for (i=ix; i<lx; i++)
    1474             :     {
    1475     1592919 :       long e = gexpo(gel(x,i));
    1476     1592919 :       if (e > ex) { ex = e; k = i; }
    1477             :     }
    1478             :   }
    1479      661737 :   if (!k) return lx;
    1480      657593 :   p = gel(x,k);
    1481      657593 :   r = gel(x0,k); if (isrationalzero(r)) r = x0;
    1482      657593 :   return approx_0(p, r)? lx: k;
    1483             : }
    1484             : static long
    1485       63091 : gauss_get_pivot_padic(GEN X, GEN p, long ix, GEN c)
    1486             : {
    1487       63091 :   GEN x = gel(X, ix);
    1488       63091 :   long i, k = 0, ex = (long)HIGHVALPBIT, lx = lg(x);
    1489       63091 :   if (c)
    1490             :   {
    1491         504 :     for (i=1; i<lx; i++)
    1492         378 :       if (!c[i] && !gequal0(gel(x,i)))
    1493             :       {
    1494         245 :         long e = gvaluation(gel(x,i), p);
    1495         245 :         if (e < ex) { ex = e; k = i; }
    1496             :       }
    1497             :   }
    1498             :   else
    1499             :   {
    1500      443905 :     for (i=ix; i<lx; i++)
    1501      380940 :       if (!gequal0(gel(x,i)))
    1502             :       {
    1503      182434 :         long e = gvaluation(gel(x,i), p);
    1504      182434 :         if (e < ex) { ex = e; k = i; }
    1505             :       }
    1506             :   }
    1507       63091 :   return k? k: lx;
    1508             : }
    1509             : static long
    1510      292181 : gauss_get_pivot_NZ(GEN X, GEN x0/*unused*/, long ix, GEN c)
    1511             : {
    1512      292181 :   GEN x = gel(X, ix);
    1513      292181 :   long i, lx = lg(x);
    1514             :   (void)x0;
    1515      292181 :   if (c)
    1516             :   {
    1517     1278123 :     for (i=1; i<lx; i++)
    1518     1257634 :       if (!c[i] && !gequal0(gel(x,i))) return i;
    1519             :   }
    1520             :   else
    1521             :   {
    1522      292280 :     for (i=ix; i<lx; i++)
    1523      292273 :       if (!gequal0(gel(x,i))) return i;
    1524             :   }
    1525       20496 :   return lx;
    1526             : }
    1527             : 
    1528             : /* Return pivot seeking function appropriate for the domain of the RgM x
    1529             :  * (first non zero pivot, maximal pivot...)
    1530             :  * x0 is a reference point used when guessing whether x[i,j] ~ 0
    1531             :  * (iff x[i,j] << x0[i,j]); typical case: mateigen, Gauss pivot on x - vp.Id,
    1532             :  * but use original x when deciding whether a prospective pivot is non-0 */
    1533             : static pivot_fun
    1534      273219 : get_pivot_fun(GEN x, GEN x0, GEN *data)
    1535             : {
    1536      273219 :   long i, j, hx, lx = lg(x);
    1537      273219 :   int res = t_INT;
    1538      273219 :   GEN p = NULL;
    1539             : 
    1540      273219 :   *data = NULL;
    1541      273219 :   if (lx == 1) return &gauss_get_pivot_NZ;
    1542      273079 :   hx = lgcols(x);
    1543     1308071 :   for (j=1; j<lx; j++)
    1544             :   {
    1545     1035020 :     GEN xj = gel(x,j);
    1546     9169940 :     for (i=1; i<hx; i++)
    1547             :     {
    1548     8134948 :       GEN c = gel(xj,i);
    1549     8134948 :       switch(typ(c))
    1550             :       {
    1551             :         case t_REAL:
    1552     1833272 :           res = t_REAL;
    1553     1833272 :           break;
    1554             :         case t_COMPLEX:
    1555          14 :           if (typ(gel(c,1)) == t_REAL || typ(gel(c,2)) == t_REAL) res = t_REAL;
    1556          14 :           break;
    1557             :         case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT: case t_QUAD:
    1558             :         case t_POLMOD: /* exact types */
    1559     6142195 :           break;
    1560             :         case t_PADIC:
    1561      159439 :           p = gel(c,2);
    1562      159439 :           res = t_PADIC;
    1563      159439 :           break;
    1564          28 :         default: return &gauss_get_pivot_NZ;
    1565             :       }
    1566             :     }
    1567             :   }
    1568      273051 :   switch(res)
    1569             :   {
    1570      192546 :     case t_REAL: *data = x0; return &gauss_get_pivot_max;
    1571        8225 :     case t_PADIC: *data = p; return &gauss_get_pivot_padic;
    1572       72280 :     default: return &gauss_get_pivot_NZ;
    1573             :   }
    1574             : }
    1575             : 
    1576             : static GEN
    1577      488616 : get_col(GEN a, GEN b, GEN p, long li)
    1578             : {
    1579      488616 :   GEN u = cgetg(li+1,t_COL);
    1580             :   long i, j;
    1581             : 
    1582      488616 :   gel(u,li) = gdiv(gel(b,li), p);
    1583     2697777 :   for (i=li-1; i>0; i--)
    1584             :   {
    1585     2209161 :     pari_sp av = avma;
    1586     2209161 :     GEN m = gel(b,i);
    1587     2209161 :     for (j=i+1; j<=li; j++) m = gsub(m, gmul(gcoeff(a,i,j), gel(u,j)));
    1588     2209161 :     gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(a,i,i)));
    1589             :   }
    1590      488616 :   return u;
    1591             : }
    1592             : /* assume 0 <= a[i,j] < p */
    1593             : static GEN
    1594      411634 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
    1595             : {
    1596      411634 :   GEN u = cgetg(li+1,t_VECSMALL);
    1597      411634 :   ulong m = uel(b,li) % p;
    1598             :   long i,j;
    1599             : 
    1600      411634 :   uel(u,li) = (m * ucoeff(a,li,li)) % p;
    1601     6133577 :   for (i = li-1; i > 0; i--)
    1602             :   {
    1603     5721943 :     m = p - uel(b,i)%p;
    1604    66799395 :     for (j = i+1; j <= li; j++) {
    1605    61077452 :       if (m & HIGHBIT) m %= p;
    1606    61077452 :       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
    1607             :     }
    1608     5721943 :     m %= p;
    1609     5721943 :     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
    1610     5721943 :     uel(u,i) = m;
    1611             :   }
    1612      411634 :   return u;
    1613             : }
    1614             : static GEN
    1615     3031112 : Fl_get_col(GEN a, GEN b, long li, ulong p)
    1616             : {
    1617     3031112 :   GEN u = cgetg(li+1,t_VECSMALL);
    1618     3031112 :   ulong m = uel(b,li) % p;
    1619             :   long i,j;
    1620             : 
    1621     3031112 :   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
    1622    22591770 :   for (i=li-1; i>0; i--)
    1623             :   {
    1624    19560658 :     m = b[i]%p;
    1625   224283717 :     for (j = i+1; j <= li; j++)
    1626   204723059 :       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
    1627    19560658 :     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
    1628    19560658 :     uel(u,i) = m;
    1629             :   }
    1630     3031112 :   return u;
    1631             : }
    1632             : 
    1633             : /* bk -= m * bi */
    1634             : static void
    1635     4378983 : _submul(GEN b, long k, long i, GEN m)
    1636             : {
    1637     4378983 :   gel(b,k) = gsub(gel(b,k), gmul(m, gel(b,i)));
    1638     4378983 : }
    1639             : static int
    1640      719190 : init_gauss(GEN a, GEN *b, long *aco, long *li, int *iscol)
    1641             : {
    1642      719190 :   *iscol = *b ? (typ(*b) == t_COL): 0;
    1643      719190 :   *aco = lg(a) - 1;
    1644      719190 :   if (!*aco) /* a empty */
    1645             :   {
    1646         322 :     if (*b && lg(*b) != 1) pari_err_DIM("gauss");
    1647         322 :     *li = 0; return 0;
    1648             :   }
    1649      718868 :   *li = nbrows(a);
    1650      718868 :   if (*li < *aco) pari_err_INV("gauss [no left inverse]", a);
    1651      718868 :   if (*b)
    1652             :   {
    1653      705969 :     if (*li != *aco) pari_err_DIM("gauss");
    1654      705969 :     switch(typ(*b))
    1655             :     {
    1656             :       case t_MAT:
    1657       53661 :         if (lg(*b) == 1) return 0;
    1658       53661 :         *b = RgM_shallowcopy(*b);
    1659       53661 :         break;
    1660             :       case t_COL:
    1661      652308 :         *b = mkmat( leafcopy(*b) );
    1662      652308 :         break;
    1663           0 :       default: pari_err_TYPE("gauss",*b);
    1664             :     }
    1665      705969 :     if (nbrows(*b) != *li) pari_err_DIM("gauss");
    1666             :   }
    1667             :   else
    1668       12899 :     *b = matid(*li);
    1669      718868 :   return 1;
    1670             : }
    1671             : 
    1672             : static int
    1673      338695 : is_modular_solve(GEN a, GEN b, GEN *u)
    1674             : {
    1675      338695 :   GEN p = NULL;
    1676             :   ulong pp;
    1677      338695 :   if (!RgM_is_FpM(a, &p) || !p) return 0;
    1678         168 :   if (!b)
    1679             :   {
    1680          70 :     a = RgM_Fp_init(a, p, &pp);
    1681          70 :     switch(pp)
    1682             :     {
    1683             :     case 0:
    1684          14 :       a = FpM_inv(a,p);
    1685          14 :       if (a) a = FpM_to_mod(a, p);
    1686          14 :       break;
    1687             :     case 2:
    1688          35 :       a = F2m_inv(a);
    1689          35 :       if (a) a = F2m_to_mod(a);
    1690          35 :       break;
    1691             :     default:
    1692          21 :       a = Flm_inv(a,pp);
    1693          21 :       if (a) a = Flm_to_mod(a, pp);
    1694             :     }
    1695             :   }
    1696          98 :   else switch(typ(b))
    1697             :   {
    1698             :     case t_COL:
    1699          49 :       if (!RgV_is_FpV(b, &p)) return 0;
    1700          49 :       a = RgM_Fp_init(a, p, &pp);
    1701          49 :       switch(pp)
    1702             :       {
    1703             :       case 0:
    1704          14 :         b = RgC_to_FpC(b, p);
    1705          14 :         a = FpM_FpC_gauss(a,b,p);
    1706          14 :         if (a) a = FpC_to_mod(a, p);
    1707          14 :         break;
    1708             :       case 2:
    1709          14 :         b = RgV_to_F2v(b);
    1710          14 :         a = F2m_F2c_gauss(a,b);
    1711          14 :         if (a) a = F2c_to_mod(a);
    1712          14 :         break;
    1713             :       default:
    1714          21 :         b = RgV_to_Flv(b, pp);
    1715          21 :         a = Flm_Flc_gauss(a,b,pp);
    1716          21 :         if (a) a = Flc_to_mod(a, pp);
    1717          21 :         break;
    1718             :       }
    1719          49 :       break;
    1720             :     case t_MAT:
    1721          49 :       if (!RgM_is_FpM(b, &p)) return 0;
    1722          49 :       a = RgM_Fp_init(a, p, &pp);
    1723          49 :       switch(pp)
    1724             :       {
    1725             :       case 0:
    1726          14 :         b = RgM_to_FpM(b, p);
    1727          14 :         a = FpM_gauss(a,b,p);
    1728          14 :         if (a) a = FpM_to_mod(a, p);
    1729          14 :         break;
    1730             :       case 2:
    1731          14 :         b = RgM_to_F2m(b);
    1732          14 :         a = F2m_gauss(a,b);
    1733          14 :         if (a) a = F2m_to_mod(a);
    1734          14 :         break;
    1735             :       default:
    1736          21 :         b = RgM_to_Flm(b, pp);
    1737          21 :         a = Flm_gauss(a,b,pp);
    1738          21 :         if (a) a = Flm_to_mod(a, pp);
    1739          21 :         break;
    1740             :       }
    1741          49 :       break;
    1742           0 :     default: return 0;
    1743             :   }
    1744         168 :   *u = a; return 1;
    1745             : }
    1746             : /* Gaussan Elimination. If a is square, return a^(-1)*b;
    1747             :  * if a has more rows than columns and b is NULL, return c such that c a = Id.
    1748             :  * a is a (not necessarily square) matrix
    1749             :  * b is a matrix or column vector, NULL meaning: take the identity matrix,
    1750             :  *   effectively returning the inverse of a
    1751             :  * If a and b are empty, the result is the empty matrix.
    1752             :  *
    1753             :  * li: number of rows of a and b
    1754             :  * aco: number of columns of a
    1755             :  * bco: number of columns of b (if matrix)
    1756             :  */
    1757             : GEN
    1758      338695 : RgM_solve(GEN a, GEN b)
    1759             : {
    1760      338695 :   pari_sp av = avma;
    1761             :   long i, j, k, li, bco, aco;
    1762             :   int iscol;
    1763             :   pivot_fun pivot;
    1764      338695 :   GEN p, u, data, ff = NULL;
    1765             : 
    1766      338695 :   if (is_modular_solve(a,b,&u)) return gerepileupto(av, u);
    1767      338527 :   if (RgM_is_FFM(a, &ff)) {
    1768         189 :     if (!b)
    1769         105 :       return FFM_inv(a, ff);
    1770          84 :     if (typ(b) == t_COL && RgC_is_FFC(b, &ff))
    1771          56 :       return FFM_FFC_gauss(a, b, ff);
    1772          28 :     if (typ(b) == t_MAT && RgM_is_FFM(b, &ff))
    1773          28 :       return FFM_gauss(a, b, ff);
    1774             :   }
    1775      338338 :   avma = av;
    1776             : 
    1777      338338 :   if (lg(a)-1 == 2 && nbrows(a) == 2) {
    1778             :     /* 2x2 matrix, start by inverting a */
    1779       90018 :     GEN u = gcoeff(a,1,1), v = gcoeff(a,1,2);
    1780       90018 :     GEN w = gcoeff(a,2,1), x = gcoeff(a,2,2);
    1781       90018 :     GEN D = gsub(gmul(u,x), gmul(v,w)), ainv;
    1782       90018 :     if (gequal0(D)) return NULL;
    1783       90018 :     ainv = mkmat2(mkcol2(x, gneg(w)), mkcol2(gneg(v), u));
    1784       90018 :     ainv = gmul(ainv, ginv(D));
    1785       90018 :     if (b) ainv = gmul(ainv, b);
    1786       90018 :     return gerepileupto(av, ainv);
    1787             :   }
    1788             : 
    1789      248320 :   if (!init_gauss(a, &b, &aco, &li, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
    1790      248068 :   pivot = get_pivot_fun(a, a, &data);
    1791      248068 :   a = RgM_shallowcopy(a);
    1792      248068 :   bco = lg(b)-1;
    1793      248068 :   if(DEBUGLEVEL>4) err_printf("Entering gauss\n");
    1794             : 
    1795      248068 :   p = NULL; /* gcc -Wall */
    1796      884282 :   for (i=1; i<=aco; i++)
    1797             :   {
    1798             :     /* k is the line where we find the pivot */
    1799      884282 :     k = pivot(a, data, i, NULL);
    1800      884282 :     if (k > li) return NULL;
    1801      884275 :     if (k != i)
    1802             :     { /* exchange the lines s.t. k = i */
    1803      156529 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    1804      156529 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    1805             :     }
    1806      884275 :     p = gcoeff(a,i,i);
    1807      884275 :     if (i == aco) break;
    1808             : 
    1809     2409500 :     for (k=i+1; k<=li; k++)
    1810             :     {
    1811     1773286 :       GEN m = gcoeff(a,k,i);
    1812     1773286 :       if (!gequal0(m))
    1813             :       {
    1814      770100 :         m = gdiv(m,p);
    1815      770100 :         for (j=i+1; j<=aco; j++) _submul(gel(a,j),k,i,m);
    1816      770100 :         for (j=1;   j<=bco; j++) _submul(gel(b,j),k,i,m);
    1817             :       }
    1818             :     }
    1819      636214 :     if (gc_needed(av,1))
    1820             :     {
    1821           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gauss. i=%ld",i);
    1822           0 :       gerepileall(av,2, &a,&b);
    1823             :     }
    1824             :   }
    1825             : 
    1826      248061 :   if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
    1827      248061 :   u = cgetg(bco+1,t_MAT);
    1828      248061 :   for (j=1; j<=bco; j++) gel(u,j) = get_col(a,gel(b,j),p,aco);
    1829      248061 :   return gerepilecopy(av, iscol? gel(u,1): u);
    1830             : }
    1831             : 
    1832             : /* assume dim A >= 1, A invertible + upper triangular  */
    1833             : static GEN
    1834       16689 : RgM_inv_upper_ind(GEN A, long index)
    1835             : {
    1836       16689 :   long n = lg(A)-1, i = index, j;
    1837       16689 :   GEN u = zerocol(n);
    1838       16689 :   gel(u,i) = ginv(gcoeff(A,i,i));
    1839       46386 :   for (i--; i>0; i--)
    1840             :   {
    1841       29697 :     pari_sp av = avma;
    1842       29697 :     GEN m = gneg(gmul(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
    1843       29697 :     for (j=i+2; j<=n; j++) m = gsub(m, gmul(gcoeff(A,i,j),gel(u,j)));
    1844       29697 :     gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(A,i,i)));
    1845             :   }
    1846       16689 :   return u;
    1847             : }
    1848             : GEN
    1849        6013 : RgM_inv_upper(GEN A)
    1850             : {
    1851             :   long i, l;
    1852        6013 :   GEN B = cgetg_copy(A, &l);
    1853        6013 :   for (i = 1; i < l; i++) gel(B,i) = RgM_inv_upper_ind(A, i);
    1854        6013 :   return B;
    1855             : }
    1856             : 
    1857             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal  */
    1858             : static GEN
    1859          14 : FpM_inv_upper_1_ind(GEN A, long index, GEN p)
    1860             : {
    1861          14 :   long n = lg(A)-1, i = index, j;
    1862          14 :   GEN u = zerocol(n);
    1863          14 :   gel(u,i) = gen_1;
    1864          21 :   for (i--; i>0; i--)
    1865             :   {
    1866           7 :     pari_sp av = avma;
    1867           7 :     GEN m = negi(mulii(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
    1868           7 :     for (j=i+2; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    1869           7 :     gel(u,i) = gerepileuptoint(av, modii(m,p));
    1870             :   }
    1871          14 :   return u;
    1872             : }
    1873             : static GEN
    1874           7 : FpM_inv_upper_1(GEN A, GEN p)
    1875             : {
    1876             :   long i, l;
    1877           7 :   GEN B = cgetg_copy(A, &l);
    1878           7 :   for (i = 1; i < l; i++) gel(B,i) = FpM_inv_upper_1_ind(A, i, p);
    1879           7 :   return B;
    1880             : }
    1881             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
    1882             :  * reduced mod p */
    1883             : static GEN
    1884          28 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
    1885             : {
    1886          28 :   long n = lg(A)-1, i = index, j;
    1887          28 :   GEN u = const_vecsmall(n, 0);
    1888          28 :   u[i] = 1;
    1889          28 :   if (SMALL_ULONG(p))
    1890          21 :     for (i--; i>0; i--)
    1891             :     {
    1892           7 :       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
    1893           7 :       for (j=i+2; j<=n; j++)
    1894             :       {
    1895           0 :         if (m & HIGHMASK) m %= p;
    1896           0 :         m += ucoeff(A,i,j) * uel(u,j);
    1897             :       }
    1898           7 :       u[i] = Fl_neg(m % p, p);
    1899             :     }
    1900             :   else
    1901          21 :     for (i--; i>0; i--)
    1902             :     {
    1903           7 :       ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
    1904           7 :       for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
    1905           7 :       u[i] = Fl_neg(m, p);
    1906             :     }
    1907          28 :   return u;
    1908             : }
    1909             : static GEN
    1910          14 : F2m_inv_upper_1_ind(GEN A, long index)
    1911             : {
    1912          14 :   pari_sp av = avma;
    1913          14 :   long n = lg(A)-1, i = index, j;
    1914          14 :   GEN u = const_vecsmall(n, 0);
    1915          14 :   u[i] = 1;
    1916          21 :   for (i--; i>0; i--)
    1917             :   {
    1918           7 :     ulong m = F2m_coeff(A,i,i+1) & uel(u,i+1); /* j = i+1 */
    1919           7 :     for (j=i+2; j<=n; j++) m ^= F2m_coeff(A,i,j) & uel(u,j);
    1920           7 :     u[i] = m & 1;
    1921             :   }
    1922          14 :   return gerepileuptoleaf(av, Flv_to_F2v(u));
    1923             : }
    1924             : static GEN
    1925          14 : Flm_inv_upper_1(GEN A, ulong p)
    1926             : {
    1927             :   long i, l;
    1928          14 :   GEN B = cgetg_copy(A, &l);
    1929          14 :   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
    1930          14 :   return B;
    1931             : }
    1932             : static GEN
    1933           7 : F2m_inv_upper_1(GEN A)
    1934             : {
    1935             :   long i, l;
    1936           7 :   GEN B = cgetg_copy(A, &l);
    1937           7 :   for (i = 1; i < l; i++) gel(B,i) = F2m_inv_upper_1_ind(A, i);
    1938           7 :   return B;
    1939             : }
    1940             : 
    1941             : static GEN
    1942      988438 : split_realimag_col(GEN z, long r1, long r2)
    1943             : {
    1944      988438 :   long i, ru = r1+r2;
    1945      988438 :   GEN x = cgetg(ru+r2+1,t_COL), y = x + r2;
    1946     2908836 :   for (i=1; i<=r1; i++) {
    1947     1920398 :     GEN a = gel(z,i);
    1948     1920398 :     if (typ(a) == t_COMPLEX) a = gel(a,1); /* paranoia: a should be real */
    1949     1920398 :     gel(x,i) = a;
    1950             :   }
    1951     1756996 :   for (   ; i<=ru; i++) {
    1952      768558 :     GEN b, a = gel(z,i);
    1953      768558 :     if (typ(a) == t_COMPLEX) { b = gel(a,2); a = gel(a,1); } else b = gen_0;
    1954      768558 :     gel(x,i) = a;
    1955      768558 :     gel(y,i) = b;
    1956             :   }
    1957      988438 :   return x;
    1958             : }
    1959             : GEN
    1960      488967 : split_realimag(GEN x, long r1, long r2)
    1961             : {
    1962             :   long i,l; GEN y;
    1963      488967 :   if (typ(x) == t_COL) return split_realimag_col(x,r1,r2);
    1964      240767 :   y = cgetg_copy(x, &l);
    1965      240767 :   for (i=1; i<l; i++) gel(y,i) = split_realimag_col(gel(x,i), r1, r2);
    1966      240767 :   return y;
    1967             : }
    1968             : 
    1969             : /* assume M = (r1+r2) x (r1+2r2) matrix and y compatible vector or matrix
    1970             :  * r1 first lines of M,y are real. Solve the system obtained by splitting
    1971             :  * real and imaginary parts. */
    1972             : GEN
    1973      238803 : RgM_solve_realimag(GEN M, GEN y)
    1974             : {
    1975      238803 :   long l = lg(M), r2 = l - lgcols(M), r1 = l-1 - 2*r2;
    1976      238803 :   return RgM_solve(split_realimag(M, r1,r2),
    1977             :                    split_realimag(y, r1,r2));
    1978             : }
    1979             : 
    1980             : GEN
    1981         294 : gauss(GEN a, GEN b)
    1982             : {
    1983             :   GEN z;
    1984         294 :   if (typ(a)!=t_MAT) pari_err_TYPE("gauss",a);
    1985         378 :   if (RgM_is_ZM(a) && b &&
    1986         133 :       ((typ(b) == t_COL && RgV_is_ZV(b)) || (typ(b) == t_MAT && RgM_is_ZM(b))))
    1987          84 :     z = ZM_gauss(a,b);
    1988             :   else
    1989         210 :     z = RgM_solve(a,b);
    1990         294 :   if (!z) pari_err_INV("gauss",a);
    1991         231 :   return z;
    1992             : }
    1993             : 
    1994             : static GEN
    1995       50950 : F2_get_col(GEN b, GEN d, long li, long aco)
    1996             : {
    1997       50950 :   long i, l = nbits2lg(aco);
    1998       50950 :   GEN u = cgetg(l, t_VECSMALL);
    1999       50950 :   u[1] = aco;
    2000      979140 :   for (i = 1; i <= li; i++)
    2001      928190 :     if (d[i]) /* d[i] can still be 0 if li > aco */
    2002             :     {
    2003      928155 :       if (F2v_coeff(b, i))
    2004      326956 :         F2v_set(u, d[i]);
    2005             :       else
    2006      601199 :         F2v_clear(u, d[i]);
    2007             :     }
    2008       50950 :   return u;
    2009             : }
    2010             : 
    2011             : /* destroy a, b */
    2012             : static GEN
    2013        9219 : F2m_gauss_sp(GEN a, GEN b)
    2014             : {
    2015        9219 :   long i, j, k, l, li, bco, aco = lg(a)-1;
    2016             :   GEN u, d;
    2017             : 
    2018        9219 :   if (!aco) return cgetg(1,t_MAT);
    2019        9219 :   li = gel(a,1)[1];
    2020        9219 :   d = zero_Flv(li);
    2021        9219 :   bco = lg(b)-1;
    2022       60183 :   for (i=1; i<=aco; i++)
    2023             :   {
    2024       50978 :     GEN ai = vecsmall_copy(gel(a,i));
    2025       50978 :     if (!d[i] && F2v_coeff(ai, i))
    2026       26865 :       k = i;
    2027             :     else
    2028       24113 :       for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
    2029             :     /* found a pivot on row k */
    2030       50978 :     if (k > li) return NULL;
    2031       50964 :     d[k] = i;
    2032             : 
    2033             :     /* Clear k-th row but column-wise instead of line-wise */
    2034             :     /*  a_ij -= a_i1*a1j/a_11
    2035             :        line-wise grouping:  L_j -= a_1j/a_11*L_1
    2036             :        column-wise:         C_i -= a_i1/a_11*C_1
    2037             :     */
    2038       50964 :     F2v_clear(ai,k);
    2039      979161 :     for (l=1; l<=aco; l++)
    2040             :     {
    2041      928197 :       GEN al = gel(a,l);
    2042      928197 :       if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
    2043             :     }
    2044      979126 :     for (l=1; l<=bco; l++)
    2045             :     {
    2046      928162 :       GEN bl = gel(b,l);
    2047      928162 :       if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
    2048             :     }
    2049             :   }
    2050        9205 :   u = cgetg(bco+1,t_MAT);
    2051        9205 :   for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
    2052        9205 :   return u;
    2053             : }
    2054             : 
    2055             : GEN
    2056          28 : F2m_gauss(GEN a, GEN b)
    2057             : {
    2058          28 :   pari_sp av = avma;
    2059          28 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    2060          28 :   return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), F2m_copy(b)));
    2061             : }
    2062             : GEN
    2063          14 : F2m_F2c_gauss(GEN a, GEN b)
    2064             : {
    2065          14 :   pari_sp av = avma;
    2066          14 :   GEN z = F2m_gauss(a, mkmat(b));
    2067          14 :   if (!z) { avma = av; return NULL; }
    2068           7 :   if (lg(z) == 1) { avma = av; return cgetg(1,t_VECSMALL); }
    2069           7 :   return gerepileuptoleaf(av, gel(z,1));
    2070             : }
    2071             : 
    2072             : GEN
    2073          35 : F2m_inv(GEN a)
    2074             : {
    2075          35 :   pari_sp av = avma;
    2076          35 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    2077          35 :   return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), matid_F2m(gel(a,1)[1])));
    2078             : }
    2079             : 
    2080             : /* destroy a, b */
    2081             : static GEN
    2082       40865 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
    2083             : {
    2084       40865 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
    2085       40865 :   ulong det = 1;
    2086             :   GEN u;
    2087             : 
    2088       40865 :   li = nbrows(a);
    2089       40865 :   bco = lg(b)-1;
    2090      411655 :   for (i=1; i<=aco; i++)
    2091             :   {
    2092             :     ulong invpiv;
    2093             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
    2094      411655 :     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
    2095      794563 :     for (k = i; k <= li; k++)
    2096             :     {
    2097      794556 :       ulong piv = ( ucoeff(a,k,i) %= p );
    2098      794556 :       if (piv)
    2099             :       {
    2100      411648 :         ucoeff(a,k,i) = Fl_inv(piv, p);
    2101      411648 :         if (detp)
    2102             :         {
    2103           0 :           if (det & HIGHMASK) det %= p;
    2104           0 :           det *= piv;
    2105             :         }
    2106      411648 :         break;
    2107             :       }
    2108             :     }
    2109             :     /* found a pivot on line k */
    2110      411655 :     if (k > li) return NULL;
    2111      411648 :     if (k != i)
    2112             :     { /* swap lines so that k = i */
    2113      135402 :       s = -s;
    2114      135402 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    2115      135402 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    2116             :     }
    2117      411648 :     if (i == aco) break;
    2118             : 
    2119      370790 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    2120     3231793 :     for (k=i+1; k<=li; k++)
    2121             :     {
    2122     2861003 :       ulong m = ( ucoeff(a,k,i) %= p );
    2123     2861003 :       if (!m) continue;
    2124             : 
    2125      713693 :       m = Fl_mul(m, invpiv, p);
    2126      713693 :       if (m == 1) {
    2127      181003 :         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
    2128      181003 :         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
    2129             :       } else {
    2130      532690 :         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
    2131      532690 :         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
    2132             :       }
    2133             :     }
    2134             :   }
    2135       40858 :   if (detp)
    2136             :   {
    2137           0 :     det %=  p;
    2138           0 :     if (s < 0 && det) det = p - det;
    2139           0 :     *detp = det;
    2140             :   }
    2141       40858 :   u = cgetg(bco+1,t_MAT);
    2142       40858 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
    2143       40858 :   return u;
    2144             : }
    2145             : 
    2146             : /* destroy a, b */
    2147             : static GEN
    2148      728329 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p)
    2149             : {
    2150      728329 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
    2151      728329 :   ulong det = 1;
    2152             :   GEN u;
    2153             :   ulong pi;
    2154      728329 :   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
    2155      728329 :   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
    2156      687464 :   pi = get_Fl_red(p);
    2157      687464 :   li = nbrows(a);
    2158      687464 :   bco = lg(b)-1;
    2159     3031126 :   for (i=1; i<=aco; i++)
    2160             :   {
    2161             :     ulong invpiv;
    2162             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
    2163     4271442 :     for (k = i; k <= li; k++)
    2164             :     {
    2165     4271442 :       ulong piv = ucoeff(a,k,i);
    2166     4271442 :       if (piv)
    2167             :       {
    2168     3031126 :         ucoeff(a,k,i) = Fl_inv(piv, p);
    2169     3031126 :         if (detp) det = Fl_mul_pre(det, piv, p, pi);
    2170     3031126 :         break;
    2171             :       }
    2172             :     }
    2173             :     /* found a pivot on line k */
    2174     3031126 :     if (k > li) return NULL;
    2175     3031126 :     if (k != i)
    2176             :     { /* swap lines so that k = i */
    2177      437899 :       s = -s;
    2178      437899 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    2179      437899 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    2180             :     }
    2181     3031126 :     if (i == aco) break;
    2182             : 
    2183     2343662 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    2184    12124005 :     for (k=i+1; k<=li; k++)
    2185             :     {
    2186     9780343 :       ulong m = ucoeff(a,k,i);
    2187     9780343 :       if (!m) continue;
    2188             : 
    2189     3070516 :       m = Fl_mul_pre(m, invpiv, p, pi);
    2190     3070516 :       if (m == 1) {
    2191      111453 :         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
    2192      111453 :         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
    2193             :       } else {
    2194     2959063 :         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
    2195     2959063 :         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
    2196             :       }
    2197             :     }
    2198             :   }
    2199      687464 :   if (detp)
    2200             :   {
    2201           0 :     if (s < 0 && det) det = p - det;
    2202           0 :     *detp = det;
    2203             :   }
    2204      687464 :   u = cgetg(bco+1,t_MAT);
    2205      687464 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
    2206      687464 :   return u;
    2207             : }
    2208             : 
    2209             : GEN
    2210          42 : Flm_gauss(GEN a, GEN b, ulong p) {
    2211          42 :   return Flm_gauss_sp(RgM_shallowcopy(a), RgM_shallowcopy(b), NULL, p);
    2212             : }
    2213             : static GEN
    2214      701656 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
    2215      701656 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    2216      701656 :   return Flm_gauss_sp(a, matid_Flm(nbrows(a)), detp, p);
    2217             : }
    2218             : GEN
    2219      485024 : Flm_inv(GEN a, ulong p) {
    2220      485024 :   return Flm_inv_sp(RgM_shallowcopy(a), NULL, p);
    2221             : }
    2222             : GEN
    2223          21 : Flm_Flc_gauss(GEN a, GEN b, ulong p) {
    2224          21 :   pari_sp av = avma;
    2225          21 :   GEN z = Flm_gauss(a, mkmat(b), p);
    2226          21 :   if (!z) { avma = av; return NULL; }
    2227          14 :   if (lg(z) == 1) { avma = av; return cgetg(1,t_VECSMALL); }
    2228          14 :   return gerepileuptoleaf(av, gel(z,1));
    2229             : }
    2230             : 
    2231             : static GEN
    2232        2807 : FpM_gauss_gen(GEN a, GEN b, GEN p)
    2233             : {
    2234             :   void *E;
    2235        2807 :   const struct bb_field *S = get_Fp_field(&E,p);
    2236        2807 :   return gen_Gauss(a,b, E, S);
    2237             : }
    2238             : /* a an FpM, lg(a)>1; b an FpM or NULL (replace by identity) */
    2239             : static GEN
    2240       38594 : FpM_gauss_i(GEN a, GEN b, GEN p, ulong *pp)
    2241             : {
    2242       38594 :   long n = nbrows(a);
    2243       38594 :   a = FpM_init(a,p,pp);
    2244       38594 :   switch(*pp)
    2245             :   {
    2246             :   case 0:
    2247        2807 :     if (!b) b = matid(n);
    2248        2807 :     return FpM_gauss_gen(a,b,p);
    2249             :   case 2:
    2250        9156 :     if (b) b = ZM_to_F2m(b); else b = matid_F2m(n);
    2251        9156 :     return F2m_gauss_sp(a,b);
    2252             :   default:
    2253       26631 :     if (b) b = ZM_to_Flm(b, *pp); else b = matid_Flm(n);
    2254       26631 :     return Flm_gauss_sp(a,b, NULL, *pp);
    2255             :   }
    2256             : }
    2257             : GEN
    2258          14 : FpM_gauss(GEN a, GEN b, GEN p)
    2259             : {
    2260          14 :   pari_sp av = avma;
    2261             :   ulong pp;
    2262             :   GEN u;
    2263          14 :   if (lg(a) == 1 || lg(b)==1) return cgetg(1, t_MAT);
    2264          14 :   u = FpM_gauss_i(a, b, p, &pp);
    2265          14 :   if (!u) { avma = av; return NULL; }
    2266          14 :   switch(pp)
    2267             :   {
    2268          14 :   case 0: return gerepilecopy(av, u);
    2269           0 :   case 2:  u = F2m_to_ZM(u); break;
    2270           0 :   default: u = Flm_to_ZM(u); break;
    2271             :   }
    2272           0 :   return gerepileupto(av, u);
    2273             : }
    2274             : GEN
    2275       38566 : FpM_inv(GEN a, GEN p)
    2276             : {
    2277       38566 :   pari_sp av = avma;
    2278             :   ulong pp;
    2279             :   GEN u;
    2280       38566 :   if (lg(a) == 1) return cgetg(1, t_MAT);
    2281       38566 :   u = FpM_gauss_i(a, NULL, p, &pp);
    2282       38566 :   if (!u) { avma = av; return NULL; }
    2283       38566 :   switch(pp)
    2284             :   {
    2285        2779 :   case 0: return gerepilecopy(av, u);
    2286        9156 :   case 2:  u = F2m_to_ZM(u); break;
    2287       26631 :   default: u = Flm_to_ZM(u); break;
    2288             :   }
    2289       35787 :   return gerepileupto(av, u);
    2290             : }
    2291             : 
    2292             : GEN
    2293          14 : FpM_FpC_gauss(GEN a, GEN b, GEN p)
    2294             : {
    2295          14 :   pari_sp av = avma;
    2296             :   ulong pp;
    2297             :   GEN u;
    2298          14 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2299          14 :   u = FpM_gauss_i(a, mkmat(b), p, &pp);
    2300          14 :   if (!u) { avma = av; return NULL; }
    2301          14 :   switch(pp)
    2302             :   {
    2303          14 :   case 0: return gerepilecopy(av, gel(u,1));
    2304           0 :   case 2:  u = F2c_to_ZC(gel(u,1)); break;
    2305           0 :   default: u = Flc_to_ZC(gel(u,1)); break;
    2306             :   }
    2307           0 :   return gerepileupto(av, u);
    2308             : }
    2309             : 
    2310             : static GEN
    2311          56 : FlxqM_gauss_gen(GEN a, GEN b, GEN T, ulong p)
    2312             : {
    2313             :   void *E;
    2314          56 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    2315          56 :   return gen_Gauss(a, b, E, S);
    2316             : }
    2317             : GEN
    2318           7 : FlxqM_gauss(GEN a, GEN b, GEN T, ulong p)
    2319             : {
    2320           7 :   pari_sp av = avma;
    2321           7 :   long n = lg(a)-1;
    2322             :   GEN u;
    2323           7 :   if (!n || lg(b)==1) { avma = av; return cgetg(1, t_MAT); }
    2324           7 :   u = FlxqM_gauss_gen(a, b, T, p);
    2325           7 :   if (!u) { avma = av; return NULL; }
    2326           7 :   return gerepilecopy(av, u);
    2327             : }
    2328             : GEN
    2329          35 : FlxqM_inv(GEN a, GEN T, ulong p)
    2330             : {
    2331          35 :   pari_sp av = avma;
    2332             :   GEN u;
    2333          35 :   if (lg(a) == 1) { avma = av; return cgetg(1, t_MAT); }
    2334          35 :   u = FlxqM_gauss_gen(a, matid_FlxqM(nbrows(a),T,p), T,p);
    2335          35 :   if (!u) { avma = av; return NULL; }
    2336          28 :   return gerepilecopy(av, u);
    2337             : }
    2338             : GEN
    2339          14 : FlxqM_FlxqC_gauss(GEN a, GEN b, GEN T, ulong p)
    2340             : {
    2341          14 :   pari_sp av = avma;
    2342             :   GEN u;
    2343          14 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2344          14 :   u = FlxqM_gauss_gen(a, mkmat(b), T, p);
    2345          14 :   if (!u) { avma = av; return NULL; }
    2346           7 :   return gerepilecopy(av, gel(u,1));
    2347             : }
    2348             : 
    2349             : static GEN
    2350          56 : FqM_gauss_gen(GEN a, GEN b, GEN T, GEN p)
    2351             : {
    2352             :   void *E;
    2353          56 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    2354          56 :   return gen_Gauss(a,b,E,S);
    2355             : }
    2356             : GEN
    2357           7 : FqM_gauss(GEN a, GEN b, GEN T, GEN p)
    2358             : {
    2359           7 :   pari_sp av = avma;
    2360             :   GEN u;
    2361             :   long n;
    2362           7 :   if (!T) return FpM_gauss(a,b,p);
    2363           7 :   n = lg(a)-1; if (!n || lg(b)==1) return cgetg(1, t_MAT);
    2364           7 :   u = FqM_gauss_gen(a,b,T,p);
    2365           7 :   if (!u) { avma = av; return NULL; }
    2366           7 :   return gerepilecopy(av, u);
    2367             : }
    2368             : GEN
    2369          35 : FqM_inv(GEN a, GEN T, GEN p)
    2370             : {
    2371          35 :   pari_sp av = avma;
    2372             :   GEN u;
    2373          35 :   if (!T) return FpM_inv(a,p);
    2374          35 :   if (lg(a) == 1) return cgetg(1, t_MAT);
    2375          35 :   u = FqM_gauss_gen(a,matid(nbrows(a)),T,p);
    2376          35 :   if (!u) { avma = av; return NULL; }
    2377          28 :   return gerepilecopy(av, u);
    2378             : }
    2379             : GEN
    2380          14 : FqM_FqC_gauss(GEN a, GEN b, GEN T, GEN p)
    2381             : {
    2382          14 :   pari_sp av = avma;
    2383             :   GEN u;
    2384          14 :   if (!T) return FpM_FpC_gauss(a,b,p);
    2385          14 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2386          14 :   u = FqM_gauss_gen(a,mkmat(b),T,p);
    2387          14 :   if (!u) { avma = av; return NULL; }
    2388           7 :   return gerepilecopy(av, gel(u,1));
    2389             : }
    2390             : 
    2391             : /* Dixon p-adic lifting algorithm.
    2392             :  * Numer. Math. 40, 137-141 (1982), DOI: 10.1007/BF01459082 */
    2393             : GEN
    2394      470870 : ZM_gauss(GEN a, GEN b0)
    2395             : {
    2396      470870 :   pari_sp av = avma, av2;
    2397             :   int iscol;
    2398             :   long n, ncol, i, m, elim;
    2399             :   ulong p;
    2400      470870 :   GEN N, C, delta, xb, nb, nmin, res, b = b0;
    2401             :   forprime_t S;
    2402             : 
    2403      470870 :   if (!init_gauss(a, &b, &n, &ncol, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
    2404      470800 :   nb = gen_0; ncol = lg(b);
    2405      951623 :   for (i = 1; i < ncol; i++)
    2406             :   {
    2407      480823 :     GEN ni = gnorml2(gel(b, i));
    2408      480823 :     if (cmpii(nb, ni) < 0) nb = ni;
    2409             :   }
    2410      470800 :   if (!signe(nb)) { avma = av; return gcopy(b0); }
    2411      470800 :   delta = gen_1; nmin = nb;
    2412     1904371 :   for (i = 1; i <= n; i++)
    2413             :   {
    2414     1433571 :     GEN ni = gnorml2(gel(a, i));
    2415     1433571 :     if (cmpii(ni, nmin) < 0)
    2416             :     {
    2417        1803 :       delta = mulii(delta, nmin); nmin = ni;
    2418             :     }
    2419             :     else
    2420     1431768 :       delta = mulii(delta, ni);
    2421             :   }
    2422      470800 :   if (!signe(nmin)) return NULL;
    2423      470786 :   elim = expi(delta)+1;
    2424      470786 :   av2 = avma;
    2425      470786 :   init_modular_big(&S);
    2426             :   for(;;)
    2427             :   {
    2428      470786 :     p = u_forprime_next(&S);
    2429      470786 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    2430      470786 :     if (C) break;
    2431           0 :     elim -= expu(p);
    2432           0 :     if (elim < 0) return NULL;
    2433           0 :     avma = av2;
    2434           0 :   }
    2435             :   /* N.B. Our delta/lambda are SQUARES of those in the paper
    2436             :    * log(delta lambda) / log p, where lambda is 3+sqrt(5) / 2,
    2437             :    * whose log is < 1, hence + 1 (to cater for rounding errors) */
    2438      941572 :   m = (long)ceil((rtodbl(logr_abs(itor(delta,LOWDEFAULTPREC))) + 1)
    2439      470786 :                  / log((double)p));
    2440      470786 :   xb = ZlM_gauss(a, b, p, m, C);
    2441      470786 :   N = powuu(p, m);
    2442      470786 :   delta = sqrti(delta);
    2443      470786 :   if (iscol)
    2444      468234 :     res = FpC_ratlift(gel(xb,1), N, delta,delta, NULL);
    2445             :   else
    2446        2552 :     res = FpM_ratlift(xb, N, delta,delta, NULL);
    2447      470786 :   return gerepileupto(av, res);
    2448             : }
    2449             : 
    2450             : /* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */
    2451             : GEN
    2452      102520 : ZM_inv(GEN M, GEN dM)
    2453             : {
    2454      102520 :   pari_sp av2, av = avma;
    2455             :   GEN Hp,q,H;
    2456             :   ulong p;
    2457      102520 :   long lM = lg(M), stable = 0;
    2458      102520 :   int negate = 0;
    2459             :   forprime_t S;
    2460             :   pari_timer ti;
    2461             : 
    2462      102520 :   if (lM == 1) return cgetg(1,t_MAT);
    2463             : 
    2464             :   /* HACK: include dM = -1 ! */
    2465      101554 :   if (dM && is_pm1(dM))
    2466             :   {
    2467             :     /* modular algorithm computes M^{-1}, NOT multiplied by det(M) = -1.
    2468             :      * We will correct (negate) at the end. */
    2469       90919 :     if (signe(dM) < 0) negate = 1;
    2470       90919 :     dM = gen_1;
    2471             :   }
    2472      101554 :   init_modular_big(&S);
    2473      101554 :   av2 = avma;
    2474      101554 :   H = NULL;
    2475      101554 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2476      311402 :   while ((p = u_forprime_next(&S)))
    2477             :   {
    2478             :     ulong dMp;
    2479             :     GEN Mp;
    2480      209848 :     Mp = ZM_to_Flm(M,p);
    2481      209848 :     if (dM == gen_1)
    2482      184421 :       Hp = Flm_inv_sp(Mp, NULL, p);
    2483             :     else
    2484             :     {
    2485       25427 :       if (dM) {
    2486       25427 :         dMp = umodiu(dM,p); if (!dMp) continue;
    2487       25427 :         Hp = Flm_inv_sp(Mp, NULL, p);
    2488       25427 :         if (!Hp) pari_err_INV("ZM_inv", Mp);
    2489             :       } else {
    2490           0 :         Hp = Flm_inv_sp(Mp, &dMp, p);
    2491           0 :         if (!Hp) continue;
    2492             :       }
    2493       25427 :       if (dMp != 1) Flm_Fl_mul_inplace(Hp, dMp, p);
    2494             :     }
    2495             : 
    2496      209848 :     if (!H)
    2497             :     {
    2498      101554 :       H = ZM_init_CRT(Hp, p);
    2499      101554 :       q = utoipos(p);
    2500             :     }
    2501             :     else
    2502      108294 :       stable = ZM_incremental_CRT(&H, Hp, &q, p);
    2503      209848 :     if (DEBUGLEVEL>5) timer_printf(&ti, "ZM_inv mod %lu (stable=%ld)", p,stable);
    2504      209848 :     if (stable) {/* DONE ? */
    2505      101554 :       if (dM != gen_1)
    2506      112189 :       { if (ZM_isscalar(ZM_mul(M, H), dM)) break; }
    2507             :       else
    2508       90919 :       { if (ZM_isidentity(ZM_mul(M, H))) break; }
    2509             :     }
    2510             : 
    2511      108294 :     if (gc_needed(av,2))
    2512             :     {
    2513           9 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv");
    2514           9 :       gerepileall(av2, 2, &H, &q);
    2515             :     }
    2516             :   }
    2517      101554 :   if (!p) pari_err_OVERFLOW("ZM_inv [ran out of primes]");
    2518      101554 :   if (DEBUGLEVEL>5) err_printf("ZM_inv done\n");
    2519      101554 :   if (negate)
    2520           0 :     return gerepileupto(av, ZM_neg(H));
    2521             :   else
    2522      101554 :     return gerepilecopy(av, H);
    2523             : }
    2524             : 
    2525             : /* to be used when denom(M^(-1)) << det(M) and a sharp multiple is
    2526             :  * not available. Return H primitive such that M*H = den*Id */
    2527             : GEN
    2528        6713 : ZM_inv_ratlift(GEN M, GEN *pden)
    2529             : {
    2530        6713 :   pari_sp av2, av = avma;
    2531             :   GEN Hp, q, H;
    2532             :   ulong p;
    2533        6713 :   long lM = lg(M);
    2534             :   forprime_t S;
    2535        6713 :   if (lM == 1) { *pden = gen_1; return cgetg(1,t_MAT); }
    2536             : 
    2537        6020 :   init_modular_big(&S);
    2538        6020 :   av2 = avma;
    2539        6020 :   H = NULL;
    2540       12804 :   while ((p = u_forprime_next(&S)))
    2541             :   {
    2542             :     GEN Mp, B, Hr;
    2543        6784 :     Mp = ZM_to_Flm(M,p);
    2544        6784 :     Hp = Flm_inv_sp(Mp, NULL, p);
    2545        6784 :     if (!Hp) continue;
    2546        6784 :     if (!H)
    2547             :     {
    2548        6020 :       H = ZM_init_CRT(Hp, p);
    2549        6020 :       q = utoipos(p);
    2550             :     }
    2551             :     else
    2552         764 :       ZM_incremental_CRT(&H, Hp, &q, p);
    2553        6784 :     B = sqrti(shifti(q,-1));
    2554        6784 :     Hr = FpM_ratlift(H,q,B,B,NULL);
    2555        6784 :     if (DEBUGLEVEL>5) err_printf("ZM_inv mod %lu (ratlift=%ld)\n", p,!!Hr);
    2556        6784 :     if (Hr) {/* DONE ? */
    2557        6083 :       GEN Hl = Q_remove_denom(Hr, pden);
    2558        6083 :       if (*pden)
    2559        4676 :       { if (ZM_isscalar(ZM_mul(M, Hl), *pden)) { H = Hl; break; }}
    2560             :       else
    2561        1407 :       { if (ZM_isidentity(ZM_mul(M, Hl))) { H = Hl; *pden = gen_1; break; } }
    2562             :     }
    2563             : 
    2564         764 :     if (gc_needed(av,2))
    2565             :     {
    2566           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv_ratlift");
    2567           0 :       gerepileall(av2, 2, &H, &q);
    2568             :     }
    2569             :   }
    2570        6020 :   gerepileall(av, 2, &H, pden);
    2571        6020 :   return H;
    2572             : }
    2573             : 
    2574             : /* same as above, M rational */
    2575             : GEN
    2576        6426 : QM_inv(GEN M, GEN dM)
    2577             : {
    2578        6426 :   pari_sp av = avma;
    2579        6426 :   GEN cM, pM = Q_primitive_part(M, &cM);
    2580        6426 :   if (!cM) return ZM_inv(pM,dM);
    2581        2632 :   return gerepileupto(av, ZM_inv(pM, gdiv(dM,cM)));
    2582             : }
    2583             : 
    2584             : GEN
    2585        1631 : ZM_ker_ratlift(GEN M)
    2586             : {
    2587        1631 :   pari_sp av2, av = avma;
    2588             :   GEN q, H, D;
    2589             :   forprime_t S;
    2590        1631 :   av2 = avma;
    2591        1631 :   H = NULL; D = NULL;
    2592        1631 :   if (lg(M)==1) return cgetg(1, t_MAT);
    2593        1624 :   init_modular_big(&S);
    2594             :   for(;;)
    2595             :   {
    2596             :     GEN Kp, Hp, Dp, Mp, Hr, B;
    2597        1861 :     ulong p = u_forprime_next(&S);
    2598        1861 :     Mp = ZM_to_Flm(M, p);
    2599        1861 :     Kp = Flm_ker_sp(Mp, p, 2);
    2600        1861 :     Hp = gel(Kp,1); Dp = gel(Kp,2);
    2601        1861 :     if (H && (lg(Hp)>lg(H) || (lg(Hp)==lg(H) && vecsmall_lexcmp(Dp,D)>0))) continue;
    2602        1861 :     if (!H || (lg(Hp)<lg(H) || vecsmall_lexcmp(Dp,D)<0))
    2603             :     {
    2604        1624 :       H = ZM_init_CRT(Hp, p); D = Dp;
    2605        1624 :       q = utoipos(p);
    2606             :     }
    2607             :     else
    2608         237 :       ZM_incremental_CRT(&H, Hp, &q, p);
    2609        1861 :     B = sqrti(shifti(q,-1));
    2610        1861 :     Hr = FpM_ratlift(H, q, B, B, NULL);
    2611        1861 :     if (DEBUGLEVEL>5) err_printf("ZM_ker mod %lu (ratlift=%ld)\n", p,!!Hr);
    2612        1861 :     if (Hr) {/* DONE ? */
    2613        1643 :       GEN Hl = vec_Q_primpart(Q_remove_denom(Hr, NULL));
    2614        1643 :       GEN MH = ZM_mul(M, Hl);
    2615        1643 :       if (gequal0(MH)) { H = Hl;  break; }
    2616             :     }
    2617         237 :     if (gc_needed(av,2))
    2618             :     {
    2619           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_ker_ratlift");
    2620           0 :       gerepileall(av2, 3, &H, &D, &q);
    2621             :     }
    2622         237 :   }
    2623        1624 :   return gerepilecopy(av, H);
    2624             : }
    2625             : 
    2626             : /* x a ZM. Return a multiple of the determinant of the lattice generated by
    2627             :  * the columns of x. From Algorithm 2.2.6 in GTM138 */
    2628             : GEN
    2629          42 : detint(GEN A)
    2630             : {
    2631          42 :   if (typ(A) != t_MAT) pari_err_TYPE("detint",A);
    2632          42 :   RgM_check_ZM(A, "detint");
    2633          42 :   return ZM_detmult(A);
    2634             : }
    2635             : GEN
    2636        7623 : ZM_detmult(GEN A)
    2637             : {
    2638        7623 :   pari_sp av1, av = avma;
    2639             :   GEN B, c, v, piv;
    2640        7623 :   long rg, i, j, k, m, n = lg(A) - 1;
    2641             : 
    2642        7623 :   if (!n) return gen_1;
    2643        7623 :   m = nbrows(A);
    2644        7623 :   if (n < m) return gen_0;
    2645        7602 :   c = zero_zv(m);
    2646        7602 :   av1 = avma;
    2647        7602 :   B = zeromatcopy(m,m);
    2648        7602 :   v = cgetg(m+1, t_COL);
    2649        7602 :   piv = gen_1; rg = 0;
    2650       52927 :   for (k=1; k<=n; k++)
    2651             :   {
    2652       52920 :     GEN pivprec = piv;
    2653       52920 :     long t = 0;
    2654      461895 :     for (i=1; i<=m; i++)
    2655             :     {
    2656      408975 :       pari_sp av2 = avma;
    2657             :       GEN vi;
    2658      408975 :       if (c[i]) continue;
    2659             : 
    2660      231000 :       vi = mulii(piv, gcoeff(A,i,k));
    2661     2141580 :       for (j=1; j<=m; j++)
    2662     1910580 :         if (c[j]) vi = addii(vi, mulii(gcoeff(B,j,i),gcoeff(A,j,k)));
    2663      231000 :       if (!t && signe(vi)) t = i;
    2664      231000 :       gel(v,i) = gerepileuptoint(av2, vi);
    2665             :     }
    2666       52920 :     if (!t) continue;
    2667             :     /* at this point c[t] = 0 */
    2668             : 
    2669       52892 :     if (++rg >= m) { /* full rank; mostly done */
    2670        7595 :       GEN det = gel(v,t); /* last on stack */
    2671        7595 :       if (++k > n)
    2672        7539 :         det = absi(det);
    2673             :       else
    2674             :       {
    2675             :         /* improve further; at this point c[i] is set for all i != t */
    2676          56 :         gcoeff(B,t,t) = piv; v = centermod(gel(B,t), det);
    2677         140 :         for ( ; k<=n; k++)
    2678          84 :           det = gcdii(det, ZV_dotproduct(v, gel(A,k)));
    2679             :       }
    2680        7595 :       return gerepileuptoint(av, det);
    2681             :     }
    2682             : 
    2683       45297 :     piv = gel(v,t);
    2684      401345 :     for (i=1; i<=m; i++)
    2685             :     {
    2686             :       GEN mvi;
    2687      356048 :       if (c[i] || i == t) continue;
    2688             : 
    2689      178024 :       gcoeff(B,t,i) = mvi = negi(gel(v,i));
    2690     1679223 :       for (j=1; j<=m; j++)
    2691     1501199 :         if (c[j]) /* implies j != t */
    2692             :         {
    2693      381717 :           pari_sp av2 = avma;
    2694      381717 :           GEN z = addii(mulii(gcoeff(B,j,i), piv), mulii(gcoeff(B,j,t), mvi));
    2695      381717 :           if (rg > 1) z = diviiexact(z, pivprec);
    2696      381717 :           gcoeff(B,j,i) = gerepileuptoint(av2, z);
    2697             :         }
    2698             :     }
    2699       45297 :     c[t] = k;
    2700       45297 :     if (gc_needed(av,1))
    2701             :     {
    2702           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"detint. k=%ld",k);
    2703           0 :       gerepileall(av1, 2, &piv,&B); v = zerovec(m);
    2704             :     }
    2705             :   }
    2706           7 :   avma = av; return gen_0;
    2707             : }
    2708             : 
    2709             : /* Reduce x modulo (invertible) y */
    2710             : GEN
    2711        6398 : closemodinvertible(GEN x, GEN y)
    2712             : {
    2713        6398 :   return gmul(y, ground(RgM_solve(y,x)));
    2714             : }
    2715             : GEN
    2716           7 : reducemodinvertible(GEN x, GEN y)
    2717             : {
    2718           7 :   return gsub(x, closemodinvertible(x,y));
    2719             : }
    2720             : GEN
    2721           0 : reducemodlll(GEN x,GEN y)
    2722             : {
    2723           0 :   return reducemodinvertible(x, ZM_lll(y, 0.75, LLL_INPLACE));
    2724             : }
    2725             : 
    2726             : /*******************************************************************/
    2727             : /*                                                                 */
    2728             : /*                    KERNEL of an m x n matrix                    */
    2729             : /*          return n - rk(x) linearly independent vectors          */
    2730             : /*                                                                 */
    2731             : /*******************************************************************/
    2732             : /* x has INTEGER coefficients. Gauss-Bareiss */
    2733             : GEN
    2734         246 : keri(GEN x)
    2735             : {
    2736             :   pari_sp av, av0;
    2737             :   GEN c, l, y, p, pp;
    2738             :   long i, j, k, r, t, n, m;
    2739             : 
    2740         246 :   n = lg(x)-1; if (!n) return cgetg(1,t_MAT);
    2741         246 :   av0 = avma; m = nbrows(x);
    2742         246 :   pp = cgetg(n+1,t_COL);
    2743         246 :   x = RgM_shallowcopy(x);
    2744         246 :   c = zero_zv(m);
    2745         246 :   l = cgetg(n+1, t_VECSMALL);
    2746         246 :   av = avma;
    2747        1271 :   for (r=0, p=gen_1, k=1; k<=n; k++)
    2748             :   {
    2749        1025 :     j = 1;
    2750        1025 :     while ( j <= m && (c[j] || !signe(gcoeff(x,j,k))) ) j++;
    2751        1025 :     if (j > m)
    2752             :     {
    2753         253 :       r++; l[k] = 0;
    2754        1046 :       for(j=1; j<k; j++)
    2755         793 :         if (l[j]) gcoeff(x,l[j],k) = gclone(gcoeff(x,l[j],k));
    2756         253 :       gel(pp,k) = gclone(p);
    2757             :     }
    2758             :     else
    2759             :     {
    2760         772 :       GEN p0 = p;
    2761         772 :       c[j] = k; l[k] = j; p = gcoeff(x,j,k);
    2762        4656 :       for (t=1; t<=m; t++)
    2763        3884 :         if (t!=j)
    2764             :         {
    2765        3112 :           GEN q = gcoeff(x,t,k);
    2766       12337 :           for (i=k+1; i<=n; i++)
    2767             :           {
    2768        9225 :             pari_sp av1 = avma;
    2769        9225 :             GEN p1 = subii(mulii(p,gcoeff(x,t,i)), mulii(q,gcoeff(x,j,i)));
    2770        9225 :             gcoeff(x,t,i) = gerepileuptoint(av1, diviiexact(p1,p0));
    2771             :           }
    2772        3112 :           if (gc_needed(av,1))
    2773             :           {
    2774           0 :             GEN _p0 = gclone(p0);
    2775           0 :             GEN _p  = gclone(p);
    2776           0 :             gerepile_gauss_ker(x,k,t,av);
    2777           0 :             p = icopy(_p);  gunclone(_p);
    2778           0 :             p0= icopy(_p0); gunclone(_p0);
    2779             :           }
    2780             :         }
    2781             :     }
    2782             :   }
    2783         246 :   if (!r) { avma = av0; y = cgetg(1,t_MAT); return y; }
    2784             : 
    2785             :   /* non trivial kernel */
    2786         246 :   y = cgetg(r+1,t_MAT);
    2787         499 :   for (j=k=1; j<=r; j++,k++)
    2788             :   {
    2789         253 :     p = cgetg(n+1, t_COL);
    2790         253 :     gel(y,j) = p; while (l[k]) k++;
    2791        1046 :     for (i=1; i<k; i++)
    2792         793 :       if (l[i])
    2793             :       {
    2794         786 :         c = gcoeff(x,l[i],k);
    2795         786 :         gel(p,i) = icopy(c); gunclone(c);
    2796             :       }
    2797             :       else
    2798           7 :         gel(p,i) = gen_0;
    2799         253 :     gel(p,k) = negi(gel(pp,k)); gunclone(gel(pp,k));
    2800         253 :     for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
    2801             :   }
    2802         246 :   return gerepileupto(av0, y);
    2803             : }
    2804             : 
    2805             : static GEN
    2806         126 : deplin_aux(GEN x0)
    2807             : {
    2808         126 :   pari_sp av = avma;
    2809         126 :   long i, j, k, nl, nc = lg(x0)-1;
    2810             :   GEN D, x, y, c, l, d, ck;
    2811             : 
    2812         126 :   if (!nc) { avma=av; return cgetg(1,t_COL); }
    2813          91 :   x = RgM_shallowcopy(x0);
    2814          91 :   nl = nbrows(x);
    2815          91 :   d = const_vec(nl, gen_1); /* pivot list */
    2816          91 :   c = zero_zv(nl);
    2817          91 :   l = cgetg(nc+1, t_VECSMALL); /* not initialized */
    2818          91 :   ck = NULL; /* gcc -Wall */
    2819         399 :   for (k=1; k<=nc; k++)
    2820             :   {
    2821         378 :     ck = gel(x,k);
    2822        1176 :     for (j=1; j<k; j++)
    2823             :     {
    2824         798 :       GEN cj = gel(x,j), piv = gel(d,j), q = gel(ck,l[j]);
    2825        9828 :       for (i=1; i<=nl; i++)
    2826        9030 :         if (i!=l[j]) gel(ck,i) = gsub(gmul(piv, gel(ck,i)), gmul(q, gel(cj,i)));
    2827             :     }
    2828             : 
    2829         378 :     i = gauss_get_pivot_NZ(x, NULL, k, c);
    2830         378 :     if (i > nl) break;
    2831             : 
    2832         308 :     gel(d,k) = gel(ck,i);
    2833         308 :     c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */
    2834             :   }
    2835          91 :   if (k > nc) { avma = av; return cgetg(1,t_COL); }
    2836          70 :   if (k == 1) { avma = av; return scalarcol_shallow(gen_1,nc); }
    2837          70 :   y = cgetg(nc+1,t_COL);
    2838          70 :   gel(y,1) = gcopy(gel(ck, l[1]));
    2839         203 :   for (D=gel(d,1),j=2; j<k; j++)
    2840             :   {
    2841         133 :     gel(y,j) = gmul(gel(ck, l[j]), D);
    2842         133 :     D = gmul(D, gel(d,j));
    2843             :   }
    2844          70 :   gel(y,j) = gneg(D);
    2845          70 :   for (j++; j<=nc; j++) gel(y,j) = gen_0;
    2846          70 :   y = primitive_part(y, &c);
    2847          70 :   return c? gerepileupto(av, y): gerepilecopy(av, y);
    2848             : }
    2849             : static GEN
    2850           0 : RgV_deplin(GEN v)
    2851             : {
    2852           0 :   pari_sp av = avma;
    2853           0 :   long n = lg(v)-1;
    2854           0 :   GEN y, p = NULL;
    2855           0 :   if (n <= 1)
    2856             :   {
    2857           0 :     if (n == 1 && gequal0(gel(v,1))) return mkcol(gen_1);
    2858           0 :     return cgetg(1, t_COL);
    2859             :   }
    2860           0 :   if (gequal0(gel(v,1))) return scalarcol_shallow(gen_1, n);
    2861           0 :   v = primpart(mkvec2(gel(v,1),gel(v,2)));
    2862           0 :   if (RgV_is_FpV(v, &p) && p) v = centerlift(v);
    2863           0 :   y = zerocol(n);
    2864           0 :   gel(y,1) = gneg(gel(v,2));
    2865           0 :   gel(y,2) = gcopy(gel(v,1));
    2866           0 :   return gerepileupto(av, y);
    2867             : 
    2868             : }
    2869             : GEN
    2870         189 : deplin(GEN x)
    2871             : {
    2872         189 :   GEN p = NULL;
    2873         189 :   switch(typ(x))
    2874             :   {
    2875         189 :     case t_MAT: break;
    2876           0 :     case t_VEC: return RgV_deplin(x);
    2877           0 :     default: pari_err_TYPE("deplin",x);
    2878             :   }
    2879         189 :   if (RgM_is_FpM(x, &p) && p)
    2880             :   {
    2881          63 :     pari_sp av = avma;
    2882             :     ulong pp;
    2883          63 :     x = RgM_Fp_init(x, p, &pp);
    2884          63 :     switch(pp)
    2885             :     {
    2886             :     case 0:
    2887          14 :       x = FpM_ker_gen(x,p,1);
    2888          14 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2889           7 :       x = FpC_center(x,p,shifti(p,-1));
    2890           7 :       break;
    2891             :     case 2:
    2892          14 :       x = F2m_ker_sp(x,1);
    2893          14 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2894           7 :       x = F2c_to_ZC(x); break;
    2895             :     default:
    2896          35 :       x = Flm_ker_sp(x,pp,1);
    2897          35 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2898          21 :       x = Flv_center(x, pp, pp>>1);
    2899          21 :       x = zc_to_ZC(x);
    2900          21 :       break;
    2901             :     }
    2902          35 :     return gerepileupto(av, x);
    2903             :   }
    2904         126 :   return deplin_aux(x);
    2905             : }
    2906             : /*******************************************************************/
    2907             : /*                                                                 */
    2908             : /*         GAUSS REDUCTION OF MATRICES  (m lines x n cols)         */
    2909             : /*           (kernel, image, complementary image, rank)            */
    2910             : /*                                                                 */
    2911             : /*******************************************************************/
    2912             : /* return the transform of x under a standard Gauss pivot.
    2913             :  * x0 is a reference point when guessing whether x[i,j] ~ 0
    2914             :  * (iff x[i,j] << x0[i,j])
    2915             :  * Set r = dim ker(x). d[k] contains the index of the first non-zero pivot
    2916             :  * in column k */
    2917             : static GEN
    2918        8050 : gauss_pivot_ker(GEN x, GEN x0, GEN *dd, long *rr)
    2919             : {
    2920             :   GEN c, d, p, data;
    2921             :   pari_sp av;
    2922             :   long i, j, k, r, t, n, m;
    2923             :   pivot_fun pivot;
    2924             : 
    2925        8050 :   n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
    2926        8015 :   m=nbrows(x); r=0;
    2927        8015 :   pivot = get_pivot_fun(x, x0, &data);
    2928        8015 :   x = RgM_shallowcopy(x);
    2929        8015 :   c = zero_zv(m);
    2930        8015 :   d = cgetg(n+1,t_VECSMALL);
    2931        8015 :   av=avma;
    2932       48230 :   for (k=1; k<=n; k++)
    2933             :   {
    2934       40215 :     j = pivot(x, data, k, c);
    2935       40215 :     if (j > m)
    2936             :     {
    2937       13573 :       r++; d[k]=0;
    2938       82054 :       for(j=1; j<k; j++)
    2939       68481 :         if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
    2940             :     }
    2941             :     else
    2942             :     { /* pivot for column k on row j */
    2943       26642 :       c[j]=k; d[k]=j; p = gdiv(gen_m1,gcoeff(x,j,k));
    2944       26642 :       gcoeff(x,j,k) = gen_m1;
    2945             :       /* x[j,] /= - x[j,k] */
    2946       26642 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
    2947     1817095 :       for (t=1; t<=m; t++)
    2948     1790453 :         if (t!=j)
    2949             :         { /* x[t,] -= 1 / x[j,k] x[j,] */
    2950     1763811 :           p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
    2951    28946568 :           for (i=k+1; i<=n; i++)
    2952    27182757 :             gcoeff(x,t,i) = gadd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
    2953     1763811 :           if (gc_needed(av,1)) gerepile_gauss_ker(x,k,t,av);
    2954             :         }
    2955             :     }
    2956             :   }
    2957        8015 :   *dd=d; *rr=r; return x;
    2958             : }
    2959             : 
    2960             : static GEN FpM_gauss_pivot(GEN x, GEN p, long *rr);
    2961             : static GEN FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr);
    2962             : static GEN F2m_gauss_pivot(GEN x, long *rr);
    2963             : static GEN Flm_gauss_pivot(GEN x, ulong p, long *rry);
    2964             : 
    2965             : /* r = dim ker(x).
    2966             :  * Returns d:
    2967             :  *   d[k] != 0 contains the index of a non-zero pivot in column k
    2968             :  *   d[k] == 0 if column k is a linear combination of the (k-1) first ones */
    2969             : GEN
    2970       15490 : RgM_pivots(GEN x0, GEN data, long *rr, pivot_fun pivot)
    2971             : {
    2972             :   GEN x, c, d, p;
    2973       15490 :   long i, j, k, r, t, m, n = lg(x0)-1;
    2974             :   pari_sp av;
    2975             : 
    2976       15490 :   if (RgM_is_ZM(x0)) return ZM_pivots(x0, rr);
    2977       13530 :   if (!n) { *rr = 0; return NULL; }
    2978             : 
    2979       13530 :   d = cgetg(n+1, t_VECSMALL);
    2980       13530 :   x = RgM_shallowcopy(x0);
    2981       13530 :   m = nbrows(x); r = 0;
    2982       13530 :   c = zero_zv(m);
    2983       13530 :   av = avma;
    2984      884864 :   for (k=1; k<=n; k++)
    2985             :   {
    2986      871334 :     j = pivot(x, data, k, c);
    2987      871334 :     if (j > m) { r++; d[k] = 0; }
    2988             :     else
    2989             :     {
    2990       34503 :       c[j] = k; d[k] = j; p = gdiv(gen_m1, gcoeff(x,j,k));
    2991       34503 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
    2992             : 
    2993      267834 :       for (t=1; t<=m; t++)
    2994      233331 :         if (!c[t]) /* no pivot on that line yet */
    2995             :         {
    2996      123401 :           p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
    2997     6777226 :           for (i=k+1; i<=n; i++)
    2998     6653825 :             gcoeff(x,t,i) = gadd(gcoeff(x,t,i), gmul(p, gcoeff(x,j,i)));
    2999      123401 :           if (gc_needed(av,1)) gerepile_gauss(x,k,t,av,j,c);
    3000             :         }
    3001       34503 :       for (i=k; i<=n; i++) gcoeff(x,j,i) = gen_0; /* dummy */
    3002             :     }
    3003             :   }
    3004       13530 :   *rr = r; avma = (pari_sp)d; return d;
    3005             : }
    3006             : 
    3007             : static long
    3008       99465 : ZM_count_0_cols(GEN M)
    3009             : {
    3010       99465 :   long i, l = lg(M), n = 0;
    3011      503941 :   for (i = 1; i < l; i++)
    3012      404476 :     if (ZV_equal0(gel(M,i))) n++;
    3013       99465 :   return n;
    3014             : }
    3015             : 
    3016             : static void indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol);
    3017             : /* As RgM_pivots, integer entries. Set *rr = dim Ker M0 */
    3018             : GEN
    3019      102690 : ZM_pivots(GEN M0, long *rr)
    3020             : {
    3021      102690 :   GEN d, dbest = NULL;
    3022             :   long m, n, i, imax, rmin, rbest, zc;
    3023      102690 :   int beenthere = 0;
    3024      102690 :   pari_sp av, av0 = avma;
    3025             :   forprime_t S;
    3026             : 
    3027      102690 :   rbest = n = lg(M0)-1;
    3028      102690 :   if (n == 0) { *rr = 0; return NULL; }
    3029       99465 :   zc = ZM_count_0_cols(M0);
    3030       99465 :   if (n == zc) { *rr = zc; return zero_zv(n); }
    3031             : 
    3032       99164 :   m = nbrows(M0);
    3033       99164 :   rmin = maxss(zc, n-m);
    3034       99164 :   init_modular(&S);
    3035       99164 :   imax = (n < (1<<4))? 1: (n>>3); /* heuristic */
    3036             : 
    3037             :   for(;;)
    3038             :   {
    3039             :     GEN row, col, M, KM, IM, RHS, X, cX;
    3040             :     long rk;
    3041      102857 :     for (av = avma, i = 0;; avma = av, i++)
    3042             :     {
    3043      102857 :       ulong p = u_forprime_next(&S);
    3044             :       long rp;
    3045      102857 :       if (!p) pari_err_OVERFLOW("ZM_pivots [ran out of primes]");
    3046      102857 :       d = Flm_gauss_pivot(ZM_to_Flm(M0, p), p, &rp);
    3047      102857 :       if (rp == rmin) { rbest = rp; goto END; } /* maximal rank, return */
    3048        6245 :       if (rp < rbest) { /* save best r so far */
    3049        2587 :         rbest = rp;
    3050        2587 :         if (dbest) gunclone(dbest);
    3051        2587 :         dbest = gclone(d);
    3052        5139 :         if (beenthere) break;
    3053             :       }
    3054        6245 :       if (!beenthere && i >= imax) break;
    3055        3693 :     }
    3056        2552 :     beenthere = 1;
    3057             :     /* Dubious case: there is (probably) a non trivial kernel */
    3058        2552 :     indexrank_all(m,n, rbest, dbest, &row, &col);
    3059        2552 :     M = rowpermute(vecpermute(M0, col), row);
    3060        2552 :     rk = n - rbest; /* (probable) dimension of image */
    3061        2552 :     IM = vecslice(M,1,rk);
    3062        2552 :     KM = vecslice(M,rk+1, n);
    3063        2552 :     M = rowslice(IM, 1,rk); /* square maximal rank */
    3064        2552 :     X = ZM_gauss(M, rowslice(KM, 1,rk));
    3065        2552 :     X = Q_remove_denom(X, &cX);
    3066        2552 :     RHS = rowslice(KM,rk+1,m);
    3067        2552 :     if (cX) RHS = ZM_Z_mul(RHS, cX);
    3068        2552 :     if (ZM_equal(ZM_mul(rowslice(IM,rk+1,m), X), RHS))
    3069             :     {
    3070        2552 :       d = vecsmall_copy(dbest);
    3071        2552 :       goto END;
    3072             :     }
    3073           0 :     avma = av;
    3074           0 :   }
    3075             : END:
    3076       99164 :   *rr = rbest; if (dbest) gunclone(dbest);
    3077       99164 :   return gerepileuptoleaf(av0, d);
    3078             : }
    3079             : 
    3080             : /* set *pr = dim Ker x */
    3081             : static GEN
    3082        7175 : gauss_pivot(GEN x, long *pr) {
    3083             :   GEN data;
    3084        7175 :   pivot_fun pivot = get_pivot_fun(x, x, &data);
    3085        7175 :   return RgM_pivots(x, data, pr, pivot);
    3086             : }
    3087             : 
    3088             : /* compute ker(x), x0 is a reference point when guessing whether x[i,j] ~ 0
    3089             :  * (iff x[i,j] << x0[i,j]) */
    3090             : static GEN
    3091        8050 : ker_aux(GEN x, GEN x0)
    3092             : {
    3093        8050 :   pari_sp av = avma;
    3094             :   GEN d,y;
    3095             :   long i,j,k,r,n;
    3096             : 
    3097        8050 :   x = gauss_pivot_ker(x,x0,&d,&r);
    3098        8050 :   if (!r) { avma=av; return cgetg(1,t_MAT); }
    3099        7693 :   n = lg(x)-1; y=cgetg(r+1,t_MAT);
    3100       21266 :   for (j=k=1; j<=r; j++,k++)
    3101             :   {
    3102       13573 :     GEN p = cgetg(n+1,t_COL);
    3103             : 
    3104       13573 :     gel(y,j) = p; while (d[k]) k++;
    3105       82054 :     for (i=1; i<k; i++)
    3106       68481 :       if (d[i])
    3107             :       {
    3108       39046 :         GEN p1=gcoeff(x,d[i],k);
    3109       39046 :         gel(p,i) = gcopy(p1); gunclone(p1);
    3110             :       }
    3111             :       else
    3112       29435 :         gel(p,i) = gen_0;
    3113       13573 :     gel(p,k) = gen_1; for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
    3114             :   }
    3115        7693 :   return gerepileupto(av,y);
    3116             : }
    3117             : GEN
    3118        7966 : ker(GEN x)
    3119             : {
    3120        7966 :   pari_sp av = avma;
    3121        7966 :   GEN p = NULL, ff = NULL;
    3122        7966 :   if (RgM_is_FpM(x, &p) && p)
    3123             :   {
    3124             :     ulong pp;
    3125          35 :     x = RgM_Fp_init(x, p, &pp);
    3126          35 :     switch(pp)
    3127             :     {
    3128          14 :     case 0: x = FpM_to_mod(FpM_ker_gen(x,p,0),p); break;
    3129           7 :     case 2: x = F2m_to_mod(F2m_ker_sp(x,0)); break;
    3130          14 :     default:x = Flm_to_mod(Flm_ker_sp(x,pp,0), pp); break;
    3131             :     }
    3132          35 :     return gerepileupto(av, x);
    3133             :   }
    3134        7931 :   if (RgM_is_FFM(x, &ff)) return FFM_ker(x, ff);
    3135        7896 :   return ker_aux(x,x);
    3136             : }
    3137             : GEN
    3138         147 : matker0(GEN x,long flag)
    3139             : {
    3140         147 :   if (typ(x)!=t_MAT) pari_err_TYPE("matker",x);
    3141         147 :   if (!flag) return ker(x);
    3142           7 :   RgM_check_ZM(x, "keri");
    3143           7 :   return keri(x);
    3144             : }
    3145             : 
    3146             : GEN
    3147        2926 : image(GEN x)
    3148             : {
    3149        2926 :   pari_sp av = avma;
    3150        2926 :   GEN d, ff = NULL, p = NULL;
    3151             :   long r;
    3152             : 
    3153        2926 :   if (typ(x)!=t_MAT) pari_err_TYPE("matimage",x);
    3154        2926 :   if (RgM_is_FpM(x, &p) && p)
    3155             :   {
    3156             :     ulong pp;
    3157          35 :     x = RgM_Fp_init(x, p, &pp);
    3158          35 :     switch(pp)
    3159             :     {
    3160          14 :     case 0: x = FpM_to_mod(FpM_image(x,p), p); break;
    3161           7 :     case 2: x = F2m_to_mod(F2m_image(x)); break;
    3162          14 :     default:x = Flm_to_mod(Flm_image(x,pp), pp);
    3163             :     }
    3164          35 :     return gerepileupto(av, x);
    3165             :   }
    3166        2891 :   if (RgM_is_FFM(x, &ff)) return FFM_image(x, ff);
    3167        2870 :   d = gauss_pivot(x,&r); /* d left on stack for efficiency */
    3168        2870 :   return image_from_pivot(x,d,r);
    3169             : }
    3170             : 
    3171             : static GEN
    3172          84 : imagecompl_aux(GEN x, GEN(*PIVOT)(GEN,long*))
    3173             : {
    3174          84 :   pari_sp av = avma;
    3175             :   GEN d,y;
    3176             :   long j,i,r;
    3177             : 
    3178          84 :   if (typ(x)!=t_MAT) pari_err_TYPE("imagecompl",x);
    3179          84 :   (void)new_chunk(lg(x) * 4 + 1); /* HACK */
    3180          84 :   d = PIVOT(x,&r); /* if (!d) then r = 0 */
    3181          84 :   avma = av; y = cgetg(r+1,t_VECSMALL);
    3182         126 :   for (i=j=1; j<=r; i++)
    3183          42 :     if (!d[i]) y[j++] = i;
    3184          84 :   return y;
    3185             : }
    3186             : GEN
    3187          84 : imagecompl(GEN x) { return imagecompl_aux(x, &gauss_pivot); }
    3188             : GEN
    3189           0 : ZM_imagecompl(GEN x) { return imagecompl_aux(x, &ZM_pivots); }
    3190             : 
    3191             : GEN
    3192        1456 : RgM_RgC_invimage(GEN A, GEN y)
    3193             : {
    3194        1456 :   pari_sp av = avma;
    3195        1456 :   long i, l = lg(A);
    3196        1456 :   GEN M, x, t, p = NULL;
    3197             : 
    3198        1456 :   if (RgM_is_FpM(A, &p) && RgV_is_FpV(y, &p) && p)
    3199             :   {
    3200             :     ulong pp;
    3201          28 :     A = RgM_Fp_init(A,p,&pp);
    3202          28 :     switch(pp)
    3203             :     {
    3204             :     case 0:
    3205           7 :       y = RgC_to_FpC(y,p);
    3206           7 :       x = FpM_FpC_invimage(A, y, p);
    3207           7 :       if (x) x = FpC_to_mod(x,p);
    3208           7 :       break;
    3209             :     case 2:
    3210           7 :       y = RgV_to_F2v(y);
    3211           7 :       x = F2m_F2c_invimage(A, y);
    3212           7 :       if (x) x = F2c_to_mod(x);
    3213           7 :       break;
    3214             :     default:
    3215          14 :       y = RgV_to_Flv(y,pp);
    3216          14 :       x = Flm_Flc_invimage(A, y, pp);
    3217          14 :       if (x) x = Flc_to_mod(x,pp);
    3218             :     }
    3219          28 :     if (!x) { avma = av; return NULL; }
    3220          28 :     return gerepileupto(av, x);
    3221             :   }
    3222             : 
    3223        1428 :   if (l==1) return NULL;
    3224        1428 :   if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
    3225             :   {
    3226        1428 :     GEN ff = NULL;
    3227        1428 :     if (RgM_is_FFM(A, &ff) && RgC_is_FFC(y, &ff))
    3228          63 :       return FFM_FFC_invimage(A, y, ff);
    3229             :   }
    3230        1365 :   M = ker(shallowconcat(A, y));
    3231        1365 :   i = lg(M)-1;
    3232        1365 :   if (!i) { avma = av; return NULL; }
    3233             : 
    3234        1365 :   x = gel(M,i); t = gel(x,l);
    3235        1365 :   if (gequal0(t)) { avma = av; return NULL; }
    3236             : 
    3237        1337 :   t = gneg_i(t); setlg(x,l);
    3238        1337 :   return gerepileupto(av, RgC_Rg_div(x, t));
    3239             : }
    3240             : GEN
    3241       39285 : FpM_FpC_invimage(GEN A, GEN y, GEN p)
    3242             : {
    3243       39285 :   pari_sp av = avma;
    3244       39285 :   long i, l = lg(A);
    3245             :   GEN M, x, t;
    3246             : 
    3247       39285 :   if (lgefint(p) == 3)
    3248             :   {
    3249       39278 :     ulong pp = p[2];
    3250       39278 :     A = ZM_to_Flm(A, pp);
    3251       39278 :     y = ZV_to_Flv(y, pp);
    3252       39278 :     x = Flm_Flc_invimage(A,y,pp);
    3253       39278 :     if (!x) { avma = av; return NULL; }
    3254       39278 :     return gerepileupto(av, Flc_to_ZC(x));
    3255             :   }
    3256           7 :   if (l==1) return NULL;
    3257           7 :   if (lg(y) != lgcols(A)) pari_err_DIM("FpM_FpC_invimage");
    3258           7 :   M = FpM_ker(shallowconcat(A,y),p);
    3259           7 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3260             : 
    3261           7 :   x = gel(M,i); t = gel(x,l);
    3262           7 :   if (!signe(t)) { avma = av; return NULL; }
    3263             : 
    3264           7 :   setlg(x,l); t = Fp_inv(negi(t),p);
    3265           7 :   if (is_pm1(t)) return gerepilecopy(av, x);
    3266           7 :   return gerepileupto(av, FpC_Fp_mul(x, t, p));
    3267             : }
    3268             : GEN
    3269       43870 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
    3270             : {
    3271       43870 :   pari_sp av = avma;
    3272       43870 :   long i, l = lg(A);
    3273             :   GEN M, x;
    3274             :   ulong t;
    3275             : 
    3276       43870 :   if (l==1) return NULL;
    3277       43870 :   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
    3278       43870 :   M = cgetg(l+1,t_MAT);
    3279       43870 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    3280       43870 :   gel(M,l) = y; M = Flm_ker(M,p);
    3281       43870 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3282             : 
    3283       43870 :   x = gel(M,i); t = x[l];
    3284       43870 :   if (!t) { avma = av; return NULL; }
    3285             : 
    3286       43870 :   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
    3287       43870 :   if (t!=1) x = Flv_Fl_mul(x, t, p);
    3288       43870 :   return gerepileuptoleaf(av, x);
    3289             : }
    3290             : GEN
    3291          21 : F2m_F2c_invimage(GEN A, GEN y)
    3292             : {
    3293          21 :   pari_sp av = avma;
    3294          21 :   long i, l = lg(A);
    3295             :   GEN M, x;
    3296             : 
    3297          21 :   if (l==1) return NULL;
    3298          21 :   if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
    3299          21 :   M = cgetg(l+1,t_MAT);
    3300          21 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    3301          21 :   gel(M,l) = y; M = F2m_ker(M);
    3302          21 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3303             : 
    3304          21 :   x = gel(M,i);
    3305          21 :   if (!F2v_coeff(x,l)) { avma = av; return NULL; }
    3306          21 :   F2v_clear(x, x[1]); x[1]--; /* remove last coord */
    3307          21 :   return gerepileuptoleaf(av, x);
    3308             : }
    3309             : 
    3310             : /* Return X such that m X = v (t_COL or t_MAT), resp. an empty t_COL / t_MAT
    3311             :  * if no solution exist */
    3312             : GEN
    3313        1568 : inverseimage(GEN m, GEN v)
    3314             : {
    3315             :   GEN y;
    3316        1568 :   if (typ(m)!=t_MAT) pari_err_TYPE("inverseimage",m);
    3317        1568 :   switch(typ(v))
    3318             :   {
    3319             :     case t_COL:
    3320        1456 :       y = RgM_RgC_invimage(m,v);
    3321        1456 :       return y? y: cgetg(1,t_COL);
    3322             :     case t_MAT:
    3323         112 :       y = RgM_invimage(m, v);
    3324         112 :       return y? y: cgetg(1,t_MAT);
    3325             :   }
    3326           0 :   pari_err_TYPE("inverseimage",v);
    3327             :   return NULL;/*LCOV_EXCL_LINE*/
    3328             : }
    3329             : 
    3330             : static GEN
    3331          21 : Flm_invimage_i(GEN A, GEN B, ulong p)
    3332             : {
    3333             :   GEN d, x, X, Y;
    3334          21 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3335          21 :   x = Flm_ker_sp(shallowconcat(Flm_neg(A,p), B), p, 0);
    3336             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3337             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3338             :    * Y has at least nB columns and full rank */
    3339          21 :   nY = lg(x)-1;
    3340          21 :   if (nY < nB) return NULL;
    3341          21 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3342          21 :   d = cgetg(nB+1, t_VECSMALL);
    3343          56 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3344             :   {
    3345          49 :     for (; j>=1; j--)
    3346          42 :       if (coeff(Y,i,j)) { d[i] = j; break; }
    3347          42 :     if (!j) return NULL;
    3348             :   }
    3349             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3350          14 :   Y = vecpermute(Y, d);
    3351          14 :   x = vecpermute(x, d);
    3352          14 :   X = rowslice(x, 1, nA);
    3353          14 :   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
    3354             : }
    3355             : 
    3356             : static GEN
    3357           7 : F2m_invimage_i(GEN A, GEN B)
    3358             : {
    3359             :   GEN d, x, X, Y;
    3360           7 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3361           7 :   x = F2m_ker_sp(shallowconcat(A, B), 0);
    3362             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3363             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3364             :    * Y has at least nB columns and full rank */
    3365           7 :   nY = lg(x)-1;
    3366           7 :   if (nY < nB) return NULL;
    3367             : 
    3368             :   /* implicitly: Y = rowslice(x, nA+1, nA+nB), nB rows */
    3369           7 :   d = cgetg(nB+1, t_VECSMALL);
    3370          21 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3371             :   {
    3372          14 :     for (; j>=1; j--)
    3373          14 :       if (F2m_coeff(x,nA+i,j)) { d[i] = j; break; } /* Y[i,j] */
    3374          14 :     if (!j) return NULL;
    3375             :   }
    3376           7 :   x = vecpermute(x, d);
    3377             : 
    3378           7 :   X = F2m_rowslice(x, 1, nA);
    3379           7 :   Y = F2m_rowslice(x, nA+1, nA+nB);
    3380           7 :   return F2m_mul(X, F2m_inv_upper_1(Y));
    3381             : }
    3382             : GEN
    3383           0 : Flm_invimage(GEN A, GEN B, ulong p)
    3384             : {
    3385           0 :   pari_sp av = avma;
    3386           0 :   GEN X = Flm_invimage_i(A,B,p);
    3387           0 :   if (!X) { avma = av; return NULL; }
    3388           0 :   return gerepileupto(av, X);
    3389             : }
    3390             : GEN
    3391           0 : F2m_invimage(GEN A, GEN B)
    3392             : {
    3393           0 :   pari_sp av = avma;
    3394           0 :   GEN X = F2m_invimage_i(A,B);
    3395           0 :   if (!X) { avma = av; return NULL; }
    3396           0 :   return gerepileupto(av, X);
    3397             : }
    3398             : static GEN
    3399          14 : FpM_invimage_i(GEN A, GEN B, GEN p)
    3400             : {
    3401             :   GEN d, x, X, Y;
    3402          14 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3403          14 :   if (lgefint(p) == 3)
    3404             :   {
    3405           0 :     ulong pp = p[2];
    3406           0 :     A = ZM_to_Flm(A, pp);
    3407           0 :     B = ZM_to_Flm(B, pp);
    3408           0 :     x = Flm_invimage_i(A, B, pp);
    3409           0 :     return x? Flm_to_ZM(x): NULL;
    3410             :   }
    3411          14 :   x = FpM_ker(shallowconcat(ZM_neg(A), B), p);
    3412             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3413             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3414             :    * Y has at least nB columns and full rank */
    3415          14 :   nY = lg(x)-1;
    3416          14 :   if (nY < nB) return NULL;
    3417          14 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3418          14 :   d = cgetg(nB+1, t_VECSMALL);
    3419          35 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3420             :   {
    3421          35 :     for (; j>=1; j--)
    3422          28 :       if (signe(gcoeff(Y,i,j))) { d[i] = j; break; }
    3423          28 :     if (!j) return NULL;
    3424             :   }
    3425             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3426           7 :   Y = vecpermute(Y, d);
    3427           7 :   x = vecpermute(x, d);
    3428           7 :   X = rowslice(x, 1, nA);
    3429           7 :   return FpM_mul(X, FpM_inv_upper_1(Y,p), p);
    3430             : }
    3431             : GEN
    3432           0 : FpM_invimage(GEN A, GEN B, GEN p)
    3433             : {
    3434           0 :   pari_sp av = avma;
    3435           0 :   GEN X = FpM_invimage_i(A,B,p);
    3436           0 :   if (!X) { avma = av; return NULL; }
    3437           0 :   return gerepileupto(av, X);
    3438             : }
    3439             : 
    3440             : /* find Z such that A Z = B. Return NULL if no solution */
    3441             : GEN
    3442         126 : RgM_invimage(GEN A, GEN B)
    3443             : {
    3444         126 :   pari_sp av = avma;
    3445             :   GEN d, x, X, Y;
    3446         126 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3447         126 :   GEN p = NULL;
    3448         126 :   if (RgM_is_FpM(A, &p) && RgM_is_FpM(B, &p) && p)
    3449             :   {
    3450             :     ulong pp;
    3451          42 :     A = RgM_Fp_init(A,p,&pp);
    3452          42 :     switch(pp)
    3453             :     {
    3454             :     case 0:
    3455          14 :       B = RgM_to_FpM(B,p);
    3456          14 :       x = FpM_invimage_i(A, B, p);
    3457          14 :       if (x) x = FpM_to_mod(x, p);
    3458          14 :     break;
    3459             :     case 2:
    3460           7 :       B = RgM_to_F2m(B);
    3461           7 :       x = F2m_invimage_i(A, B);
    3462           7 :       if (x) x = F2m_to_mod(x);
    3463           7 :       break;
    3464             :     default:
    3465          21 :       B = RgM_to_Flm(B,pp);
    3466          21 :       x = Flm_invimage_i(A, B, pp);
    3467          21 :       if (x) x = Flm_to_mod(x,pp);
    3468          21 :       break;
    3469             :     }
    3470          42 :     if (!x) { avma = av; return NULL; }
    3471          28 :     return gerepileupto(av, x);
    3472             :   }
    3473             :   {
    3474          84 :     GEN ff = NULL;
    3475          84 :     if (RgM_is_FFM(A, &ff) && RgM_is_FFM(B, &ff))
    3476          21 :       return FFM_invimage(A, B, ff);
    3477             :   }
    3478          63 :   x = ker(shallowconcat(RgM_neg(A), B));
    3479             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3480             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3481             :    * Y has at least nB columns and full rank */
    3482          63 :   nY = lg(x)-1;
    3483          63 :   if (nY < nB) { avma = av; return NULL; }
    3484          63 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3485          63 :   d = cgetg(nB+1, t_VECSMALL);
    3486         140 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3487             :   {
    3488         182 :     for (; j>=1; j--)
    3489         140 :       if (!gequal0(gcoeff(Y,i,j))) { d[i] = j; break; }
    3490         119 :     if (!j) { avma = av; return NULL; }
    3491             :   }
    3492             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3493          21 :   Y = vecpermute(Y, d);
    3494          21 :   x = vecpermute(x, d);
    3495          21 :   X = rowslice(x, 1, nA);
    3496          21 :   return gerepileupto(av, RgM_mul(X, RgM_inv_upper(Y)));
    3497             : }
    3498             : 
    3499             : /* r = dim Ker x, n = nbrows(x) */
    3500             : static GEN
    3501       30965 : get_suppl(GEN x, GEN d, long n, long r, GEN(*ei)(long,long))
    3502             : {
    3503             :   pari_sp av;
    3504             :   GEN y, c;
    3505       30965 :   long j, k, rx = lg(x)-1; /* != 0 due to init_suppl() */
    3506             : 
    3507       30965 :   if (rx == n && r == 0) return gcopy(x);
    3508       25878 :   y = cgetg(n+1, t_MAT);
    3509       25878 :   av = avma; c = zero_zv(n);
    3510             :   /* c = lines containing pivots (could get it from gauss_pivot, but cheap)
    3511             :    * In theory r = 0 and d[j] > 0 for all j, but why take chances? */
    3512      206981 :   for (k = j = 1; j<=rx; j++)
    3513      181103 :     if (d[j]) { c[ d[j] ] = 1; gel(y,k++) = gel(x,j); }
    3514      283045 :   for (j=1; j<=n; j++)
    3515      257167 :     if (!c[j]) gel(y,k++) = (GEN)j; /* HACK */
    3516       25878 :   avma = av;
    3517             : 
    3518       25878 :   rx -= r;
    3519       25878 :   for (j=1; j<=rx; j++) gel(y,j) = gcopy(gel(y,j));
    3520       25878 :   for (   ; j<=n; j++)  gel(y,j) = ei(n, y[j]);
    3521       25878 :   return y;
    3522             : }
    3523             : 
    3524             : static void
    3525       30965 : init_suppl(GEN x)
    3526             : {
    3527       30965 :   if (lg(x) == 1) pari_err_IMPL("suppl [empty matrix]");
    3528             :   /* HACK: avoid overwriting d from gauss_pivot() after avma=av */
    3529       30965 :   (void)new_chunk(lgcols(x) * 2);
    3530       30965 : }
    3531             : 
    3532             : /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
    3533             :  * whose first k columns are given by x. If rank(x) < k, undefined result. */
    3534             : GEN
    3535          70 : suppl(GEN x)
    3536             : {
    3537          70 :   pari_sp av = avma;
    3538          70 :   GEN d, X = x, p = NULL;
    3539             :   long r;
    3540             : 
    3541          70 :   if (typ(x)!=t_MAT) pari_err_TYPE("suppl",x);
    3542          70 :   if (RgM_is_FpM(x, &p) && p)
    3543             :   {
    3544             :     ulong pp;
    3545          28 :     x = RgM_Fp_init(x, p, &pp);
    3546          28 :     switch(pp)
    3547             :     {
    3548           7 :     case 0: x = FpM_to_mod(FpM_suppl(x,p), p); break;
    3549           7 :     case 2: x = F2m_to_mod(F2m_suppl(x)); break;
    3550          14 :     default:x = Flm_to_mod(Flm_suppl(x,pp), pp); break;
    3551             :     }
    3552          28 :     return gerepileupto(av, x);
    3553             :   }
    3554          42 :   avma = av; init_suppl(x);
    3555          42 :   d = gauss_pivot(X,&r);
    3556          42 :   avma = av; return get_suppl(X,d,nbrows(X),r,&col_ei);
    3557             : }
    3558             : GEN
    3559       30776 : FpM_suppl(GEN x, GEN p)
    3560             : {
    3561       30776 :   pari_sp av = avma;
    3562             :   GEN d;
    3563             :   long r;
    3564       30776 :   init_suppl(x); d = FpM_gauss_pivot(x,p, &r);
    3565       30776 :   avma = av; return get_suppl(x,d,nbrows(x),r,&col_ei);
    3566             : }
    3567             : GEN
    3568          42 : Flm_suppl(GEN x, ulong p)
    3569             : {
    3570          42 :   pari_sp av = avma;
    3571             :   GEN d;
    3572             :   long r;
    3573          42 :   init_suppl(x); d = Flm_gauss_pivot(Flm_copy(x),p, &r);
    3574          42 :   avma = av; return get_suppl(x,d,nbrows(x),r,&vecsmall_ei);
    3575             : 
    3576             : }
    3577             : GEN
    3578           7 : F2m_suppl(GEN x)
    3579             : {
    3580           7 :   pari_sp av = avma;
    3581             :   GEN d;
    3582             :   long r;
    3583           7 :   init_suppl(x); d = F2m_gauss_pivot(F2m_copy(x), &r);
    3584           7 :   avma = av; return get_suppl(x,d,mael(x,1,1),r,&F2v_ei);
    3585             : }
    3586             : 
    3587             : GEN
    3588         574 : FqM_suppl(GEN x, GEN T, GEN p)
    3589             : {
    3590         574 :   pari_sp av = avma;
    3591             :   GEN d;
    3592             :   long r;
    3593             : 
    3594         574 :   if (!T) return FpM_suppl(x,p);
    3595          98 :   init_suppl(x);
    3596          98 :   d = FqM_gauss_pivot(x,T,p,&r);
    3597          98 :   avma = av; return get_suppl(x,d,nbrows(x),r,&col_ei);
    3598             : }
    3599             : 
    3600             : GEN
    3601           7 : image2(GEN x)
    3602             : {
    3603           7 :   pari_sp av = avma;
    3604             :   long k, n, i;
    3605             :   GEN A, B;
    3606             : 
    3607           7 :   if (typ(x)!=t_MAT) pari_err_TYPE("image2",x);
    3608           7 :   if (lg(x) == 1) return cgetg(1,t_MAT);
    3609           7 :   A = ker(x); k = lg(A)-1;
    3610           7 :   if (!k) { avma = av; return gcopy(x); }
    3611           7 :   A = suppl(A); n = lg(A)-1;
    3612           7 :   B = cgetg(n-k+1, t_MAT);
    3613           7 :   for (i = k+1; i <= n; i++) gel(B,i-k) = RgM_RgC_mul(x, gel(A,i));
    3614           7 :   return gerepileupto(av, B);
    3615             : }
    3616             : 
    3617             : GEN
    3618         119 : matimage0(GEN x,long flag)
    3619             : {
    3620         119 :   switch(flag)
    3621             :   {
    3622         112 :     case 0: return image(x);
    3623           7 :     case 1: return image2(x);
    3624           0 :     default: pari_err_FLAG("matimage");
    3625             :   }
    3626             :   return NULL; /* LCOV_EXCL_LINE */
    3627             : }
    3628             : 
    3629             : long
    3630         154 : rank(GEN x)
    3631             : {
    3632         154 :   pari_sp av = avma;
    3633             :   long r;
    3634         154 :   GEN ff = NULL, p = NULL;
    3635             : 
    3636         154 :   if (typ(x)!=t_MAT) pari_err_TYPE("rank",x);
    3637         154 :   if (RgM_is_FpM(x, &p) && p)
    3638             :   {
    3639             :     ulong pp;
    3640          84 :     x = RgM_Fp_init(x,p,&pp);
    3641          84 :     switch(pp)
    3642             :     {
    3643           7 :     case 0: r = FpM_rank(x,p); break;
    3644          63 :     case 2: r = F2m_rank(x); break;
    3645          14 :     default:r = Flm_rank(x,pp); break;
    3646             :     }
    3647          84 :     avma = av; return r;
    3648             :   }
    3649          70 :   if (RgM_is_FFM(x, &ff)) return FFM_rank(x, ff);
    3650          49 :   (void)gauss_pivot(x, &r);
    3651          49 :   avma = av; return lg(x)-1 - r;
    3652             : }
    3653             : 
    3654             : /* d a t_VECSMALL of integers in 1..n. Return the vector of the d[i]
    3655             :  * followed by the missing indices */
    3656             : static GEN
    3657        5104 : perm_complete(GEN d, long n)
    3658             : {
    3659        5104 :   GEN y = cgetg(n+1, t_VECSMALL);
    3660        5104 :   long i, j = 1, k = n, l = lg(d);
    3661        5104 :   pari_sp av = avma;
    3662        5104 :   char *T = stack_calloc(n+1);
    3663        5104 :   for (i = 1; i < l; i++) T[d[i]] = 1;
    3664       65532 :   for (i = 1; i <= n; i++)
    3665       60428 :     if (T[i]) y[j++] = i; else y[k--] = i;
    3666        5104 :   avma = av; return y;
    3667             : }
    3668             : 
    3669             : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
    3670             : static GEN
    3671       20423 : indexrank0(long n, long r, GEN d)
    3672             : {
    3673       20423 :   GEN p1, p2, res = cgetg(3,t_VEC);
    3674             :   long i, j;
    3675             : 
    3676       20423 :   r = n - r; /* now r = dim Im(x) */
    3677       20423 :   p1 = cgetg(r+1,t_VECSMALL); gel(res,1) = p1;
    3678       20423 :   p2 = cgetg(r+1,t_VECSMALL); gel(res,2) = p2;
    3679       20423 :   if (d)
    3680             :   {
    3681      116519 :     for (i=0,j=1; j<=n; j++)
    3682       96824 :       if (d[j]) { i++; p1[i] = d[j]; p2[i] = j; }
    3683       19695 :     vecsmall_sort(p1);
    3684             :   }
    3685       20423 :   return res;
    3686             : }
    3687             : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
    3688             : static GEN
    3689         840 : indeximage0(long n, long r, GEN d)
    3690             : {
    3691             :   long i, j;
    3692             :   GEN v;
    3693             : 
    3694         840 :   r = n - r; /* now r = dim Im(x) */
    3695         840 :   v = cgetg(r+1,t_VECSMALL);
    3696        7651 :   if (d) for (i=j=1; j<=n; j++)
    3697        6811 :     if (d[j]) v[i++] = j;
    3698         840 :   return v;
    3699             : }
    3700             : /* x an m x n t_MAT, n > 0, r = dim Ker(x), d from gauss_pivot */
    3701             : static void
    3702        2552 : indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol)
    3703             : {
    3704        2552 :   GEN IR = indexrank0(n, r, d);
    3705        2552 :   *prow = perm_complete(gel(IR,1), m);
    3706        2552 :   *pcol = perm_complete(gel(IR,2), n);
    3707        2552 : }
    3708             : static void
    3709       18711 : init_indexrank(GEN x) {
    3710       18711 :   (void)new_chunk(3 + 2*lg(x)); /* HACK */
    3711       18711 : }
    3712             : 
    3713             : GEN
    3714        4158 : indexrank(GEN x) {
    3715        4158 :   pari_sp av = avma;
    3716             :   long r;
    3717        4158 :   GEN d, p = NULL;
    3718        4158 :   if (typ(x)!=t_MAT) pari_err_TYPE("indexrank",x);
    3719        4158 :   init_indexrank(x);
    3720        4158 :   if (RgM_is_FpM(x, &p) && p)
    3721          28 :   {
    3722             :     ulong pp;
    3723          28 :     x = RgM_Fp_init(x,p,&pp);
    3724          28 :     switch(pp)
    3725             :     {
    3726           7 :     case 0: d = FpM_gauss_pivot(x,p,&r); break;
    3727           7 :     case 2: d = F2m_gauss_pivot(x,&r); break;
    3728          14 :     default:d = Flm_gauss_pivot(x,pp,&r); break;
    3729             :     }
    3730             :   }
    3731             :   else
    3732        4130 :     d = gauss_pivot(x,&r);
    3733        4158 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3734             : }
    3735             : 
    3736             : GEN
    3737          21 : FpM_indexrank(GEN x, GEN p) {
    3738          21 :   pari_sp av = avma;
    3739             :   long r;
    3740             :   GEN d;
    3741          21 :   init_indexrank(x);
    3742          21 :   d = FpM_gauss_pivot(x,p,&r);
    3743          21 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3744             : }
    3745             : GEN
    3746       10038 : Flm_indexrank(GEN x, ulong p) {
    3747       10038 :   pari_sp av = avma;
    3748             :   long r;
    3749             :   GEN d;
    3750       10038 :   init_indexrank(x);
    3751       10038 :   d = Flm_gauss_pivot(Flm_copy(x),p,&r);
    3752       10038 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3753             : }
    3754             : GEN
    3755           0 : F2m_indexrank(GEN x) {
    3756           0 :   pari_sp av = avma;
    3757             :   long r;
    3758             :   GEN d;
    3759           0 :   init_indexrank(x);
    3760           0 :   d = F2m_gauss_pivot(F2m_copy(x),&r);
    3761           0 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3762             : }
    3763             : 
    3764             : GEN
    3765         840 : ZM_indeximage(GEN x) {
    3766         840 :   pari_sp av = avma;
    3767             :   long r;
    3768             :   GEN d;
    3769         840 :   init_indexrank(x);
    3770         840 :   d = ZM_pivots(x,&r);
    3771         840 :   avma = av; return indeximage0(lg(x)-1, r, d);
    3772             : }
    3773             : long
    3774       34018 : ZM_rank(GEN x) {
    3775       34018 :   pari_sp av = avma;
    3776             :   long r;
    3777       34018 :   (void)ZM_pivots(x,&r);
    3778       34018 :   avma = av; return lg(x)-1-r;
    3779             : }
    3780             : GEN
    3781        3654 : ZM_indexrank(GEN x) {
    3782        3654 :   pari_sp av = avma;
    3783             :   long r;
    3784             :   GEN d;
    3785        3654 :   init_indexrank(x);
    3786        3654 :   d = ZM_pivots(x,&r);
    3787        3654 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3788             : }
    3789             : 
    3790             : /*******************************************************************/
    3791             : /*                                                                 */
    3792             : /*                   Structured Elimination                        */
    3793             : /*                                                                 */
    3794             : /*******************************************************************/
    3795             : 
    3796             : static void
    3797       94025 : rem_col(GEN c, long i, GEN iscol, GEN Wrow, long *rcol, long *rrow)
    3798             : {
    3799       94025 :   long lc = lg(c), k;
    3800       94025 :   iscol[i] = 0; (*rcol)--;
    3801      832065 :   for (k = 1; k < lc; ++k)
    3802             :   {
    3803      738040 :     Wrow[c[k]]--;
    3804      738040 :     if (Wrow[c[k]]==0) (*rrow)--;
    3805             :   }
    3806       94025 : }
    3807             : 
    3808             : static void
    3809        5629 : rem_singleton(GEN M, GEN iscol, GEN Wrow, long *rcol, long *rrow)
    3810             : {
    3811             :   long i, j;
    3812        5629 :   long nbcol = lg(iscol)-1, last;
    3813             :   do
    3814             :   {
    3815        7513 :     last = 0;
    3816    16428764 :     for (i = 1; i <= nbcol; ++i)
    3817    16421251 :       if (iscol[i])
    3818             :       {
    3819     8822943 :         GEN c = gmael(M, i, 1);
    3820     8822943 :         long lc = lg(c);
    3821    82051253 :         for (j = 1; j < lc; ++j)
    3822    73240899 :           if (Wrow[c[j]] == 1)
    3823             :           {
    3824       12589 :             rem_col(c, i, iscol, Wrow, rcol, rrow);
    3825       12589 :             last=1; break;
    3826             :           }
    3827             :       }
    3828        7513 :   } while (last);
    3829        5629 : }
    3830             : 
    3831             : static GEN
    3832        5508 : fill_wcol(GEN M, GEN iscol, GEN Wrow, long *w, GEN wcol)
    3833             : {
    3834        5508 :   long nbcol = lg(iscol)-1;
    3835             :   long i, j, m, last;
    3836             :   GEN per;
    3837       13314 :   for (m = 2, last=0; !last ; m++)
    3838             :   {
    3839    17657317 :     for (i = 1; i <= nbcol; ++i)
    3840             :     {
    3841    17649511 :       wcol[i] = 0;
    3842    17649511 :       if (iscol[i])
    3843             :       {
    3844     9416700 :         GEN c = gmael(M, i, 1);
    3845     9416700 :         long lc = lg(c);
    3846    87430408 :         for (j = 1; j < lc; ++j)
    3847    78013708 :           if (Wrow[c[j]] == m) {  wcol[i]++; last = 1; }
    3848             :       }
    3849             :     }
    3850             :   }
    3851        5508 :   per = vecsmall_indexsort(wcol);
    3852        5508 :   *w = wcol[per[nbcol]];
    3853        5508 :   return per;
    3854             : }
    3855             : 
    3856             : /* M is a RgMs with nbrow rows, A a list of row indices.
    3857             :    Eliminate rows of M with a single entry that do not belong to A,
    3858             :    and the corresponding columns. Also eliminate columns until #colums=#rows.
    3859             :    Return pcol and prow:
    3860             :    pcol is a map from the new columns indices to the old one.
    3861             :    prow is a map from the old rows indices to the new one (0 if removed).
    3862             : */
    3863             : 
    3864             : void
    3865         121 : RgMs_structelim_col(GEN M, long nbcol, long nbrow, GEN A, GEN *p_col, GEN *p_row)
    3866             : {
    3867             :   long i,j,k;
    3868         121 :   long lA = lg(A);
    3869         121 :   GEN prow = cgetg(nbrow+1, t_VECSMALL);
    3870         121 :   GEN pcol = zero_zv(nbcol);
    3871         121 :   pari_sp av = avma;
    3872         121 :   long rcol = nbcol, rrow = 0, imin = nbcol - usqrt(nbcol);
    3873         121 :   GEN iscol = const_vecsmall(nbcol, 1);
    3874         121 :   GEN Wrow  = zero_zv(nbrow);
    3875         121 :   GEN wcol = cgetg(nbcol+1, t_VECSMALL);
    3876         121 :   pari_sp av2=avma;
    3877      119875 :   for (i = 1; i <= nbcol; ++i)
    3878             :   {
    3879      119754 :     GEN F = gmael(M, i, 1);
    3880      119754 :     long l = lg(F)-1;
    3881     1058234 :     for (j = 1; j <= l; ++j)
    3882      938480 :       Wrow[F[j]]++;
    3883             :   }
    3884         121 :   for (j = 1; j < lA; ++j)
    3885             :   {
    3886         121 :     if (Wrow[A[j]] == 0) { *p_col=NULL; return; }
    3887           0 :     Wrow[A[j]] = -1;
    3888             :   }
    3889      221716 :   for (i = 1; i <= nbrow; ++i)
    3890      221595 :     if (Wrow[i])
    3891       65538 :       rrow++;
    3892         121 :   rem_singleton(M, iscol, Wrow, &rcol, &rrow);
    3893         121 :   if (rcol<rrow) pari_err_BUG("RgMs_structelim, rcol<rrow");
    3894        5750 :   for (; rcol>rrow;)
    3895             :   {
    3896             :     long w;
    3897        5508 :     GEN per = fill_wcol(M, iscol, Wrow, &w, wcol);
    3898       86944 :     for (i = nbcol; i>=imin && wcol[per[i]]>=w && rcol>rrow; i--)
    3899       81436 :       rem_col(gmael(M, per[i], 1), per[i], iscol, Wrow, &rcol, &rrow);
    3900        5508 :     rem_singleton(M, iscol, Wrow, &rcol, &rrow);
    3901        5508 :     avma = av2;
    3902             :   }
    3903      119875 :   for (j = 1, i = 1; i <= nbcol; ++i)
    3904      119754 :     if (iscol[i])
    3905       25729 :       pcol[j++] = i;
    3906         121 :   setlg(pcol,j);
    3907      221716 :   for (k = 1, i = 1; i <= nbrow; ++i)
    3908      221595 :     prow[i] = Wrow[i] ? k++: 0;
    3909         121 :   avma = av;
    3910         121 :   *p_col = pcol; *p_row = prow;
    3911             : }
    3912             : 
    3913             : void
    3914           0 : RgMs_structelim(GEN M, long nbrow, GEN A, GEN *p_col, GEN *p_row)
    3915             : {
    3916           0 :   RgMs_structelim_col(M, lg(M)-1, nbrow, A, p_col, p_row);
    3917           0 : }
    3918             : 
    3919             : /*******************************************************************/
    3920             : /*                                                                 */
    3921             : /*                        EIGENVECTORS                             */
    3922             : /*   (independent eigenvectors, sorted by increasing eigenvalue)   */
    3923             : /*                                                                 */
    3924             : /*******************************************************************/
    3925             : 
    3926             : GEN
    3927          63 : mateigen(GEN x, long flag, long prec)
    3928             : {
    3929             :   GEN y, R, T;
    3930          63 :   long k, l, ex, n = lg(x);
    3931          63 :   pari_sp av = avma;
    3932             : 
    3933          63 :   if (typ(x)!=t_MAT) pari_err_TYPE("eigen",x);
    3934          63 :   if (n != 1 && n != lgcols(x)) pari_err_DIM("eigen");
    3935          63 :   if (flag < 0 || flag > 1) pari_err_FLAG("mateigen");
    3936          63 :   if (n == 1)
    3937             :   {
    3938          14 :     if (flag) retmkvec2(cgetg(1,t_VEC), cgetg(1,t_MAT));
    3939           7 :     return cgetg(1,t_VEC);
    3940             :   }
    3941          49 :   if (n == 2)
    3942             :   {
    3943          14 :     if (flag) retmkvec2(mkveccopy(gcoeff(x,1,1)), matid(1));
    3944           7 :     return matid(1);
    3945             :   }
    3946             : 
    3947          35 :   ex = 16 - prec2nbits(prec);
    3948          35 :   T = charpoly(x,0);
    3949          35 :   if (RgX_is_QX(T))
    3950             :   {
    3951          28 :     T = ZX_radical( Q_primpart(T) );
    3952          28 :     R = nfrootsQ(T);
    3953          28 :     if (lg(R)-1 < degpol(T))
    3954             :     { /* add missing complex roots */
    3955          14 :       GEN r = cleanroots(RgX_div(T, roots_to_pol(R, 0)), prec);
    3956          14 :       settyp(r, t_VEC);
    3957          14 :       R = shallowconcat(R, r);
    3958             :     }
    3959             :   }
    3960             :   else
    3961             :   {
    3962           7 :     GEN r1, v = vectrunc_init(lg(T));
    3963             :     long e;
    3964           7 :     R = cleanroots(T,prec);
    3965           7 :     r1 = NULL;
    3966          21 :     for (k = 1; k < lg(R); k++)
    3967             :     {
    3968          14 :       GEN r2 = gel(R,k), r = grndtoi(r2, &e);
    3969          14 :       if (e < ex) r2 = r;
    3970          14 :       if (r1)
    3971             :       {
    3972           7 :         r = gsub(r1,r2);
    3973           7 :         if (gequal0(r) || gexpo(r) < ex) continue;
    3974             :       }
    3975          14 :       vectrunc_append(v, r2);
    3976          14 :       r1 = r2;
    3977             :     }
    3978           7 :     R = v;
    3979             :   }
    3980             :   /* R = distinct complex roots of charpoly(x) */
    3981          35 :   l = lg(R); y = cgetg(l, t_VEC);
    3982         189 :   for (k = 1; k < l; k++)
    3983             :   {
    3984         154 :     GEN F = ker_aux(RgM_Rg_sub_shallow(x, gel(R,k)), x);
    3985         154 :     long d = lg(F)-1;
    3986         154 :     if (!d) pari_err_PREC("mateigen");
    3987         154 :     gel(y,k) = F;
    3988         154 :     if (flag) gel(R,k) = const_vec(d, gel(R,k));
    3989             :   }
    3990          35 :   y = shallowconcat1(y);
    3991          35 :   if (lg(y) > n) pari_err_PREC("mateigen");
    3992             :   /* lg(y) < n if x is not diagonalizable */
    3993          35 :   if (flag) y = mkvec2(shallowconcat1(R), y);
    3994          35 :   return gerepilecopy(av,y);
    3995             : }
    3996             : GEN
    3997           0 : eigen(GEN x, long prec) { return mateigen(x, 0, prec); }
    3998             : 
    3999             : /*******************************************************************/
    4000             : /*                                                                 */
    4001             : /*                           DETERMINANT                           */
    4002             : /*                                                                 */
    4003             : /*******************************************************************/
    4004             : 
    4005             : GEN
    4006        1617 : det0(GEN a,long flag)
    4007             : {
    4008        1617 :   switch(flag)
    4009             :   {
    4010        1603 :     case 0: return det(a);
    4011          14 :     case 1: return det2(a);
    4012           0 :     default: pari_err_FLAG("matdet");
    4013             :   }
    4014             :   return NULL; /* LCOV_EXCL_LINE */
    4015             : }
    4016             : 
    4017             : /* M a 2x2 matrix, returns det(M) */
    4018             : static GEN
    4019        2486 : RgM_det2(GEN M)
    4020             : {
    4021        2486 :   pari_sp av = avma;
    4022        2486 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    4023        2486 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    4024        2486 :   return gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
    4025             : }
    4026             : /* M a 2x2 ZM, returns det(M) */
    4027             : static GEN
    4028        7448 : ZM_det2(GEN M)
    4029             : {
    4030        7448 :   pari_sp av = avma;
    4031        7448 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    4032        7448 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    4033        7448 :   return gerepileuptoint(av, subii(mulii(a,d), mulii(b, c)));
    4034             : }
    4035             : /* M a 3x3 ZM, return det(M) */
    4036             : static GEN
    4037        2933 : ZM_det3(GEN M)
    4038             : {
    4039        2933 :   pari_sp av = avma;
    4040        2933 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2), c = gcoeff(M,1,3);
    4041        2933 :   GEN d = gcoeff(M,2,1), e = gcoeff(M,2,2), f = gcoeff(M,2,3);
    4042        2933 :   GEN g = gcoeff(M,3,1), h = gcoeff(M,3,2), i = gcoeff(M,3,3);
    4043        2933 :   GEN t, D = signe(i)? mulii(subii(mulii(a,e), mulii(b,d)), i): gen_0;
    4044        2933 :   if (signe(g))
    4045             :   {
    4046        2730 :     t = mulii(subii(mulii(b,f), mulii(c,e)), g);
    4047        2730 :     D = addii(D, t);
    4048             :   }
    4049        2933 :   if (signe(h))
    4050             :   {
    4051        2751 :     t = mulii(subii(mulii(c,d), mulii(a,f)), h);
    4052        2751 :     D = addii(D, t);
    4053             :   }
    4054        2933 :   return gerepileuptoint(av, D);
    4055             : }
    4056             : 
    4057             : static GEN
    4058        9863 : det_simple_gauss(GEN a, GEN data, pivot_fun pivot)
    4059             : {
    4060        9863 :   pari_sp av = avma;
    4061        9863 :   long i,j,k, s = 1, nbco = lg(a)-1;
    4062        9863 :   GEN p, x = gen_1;
    4063             : 
    4064        9863 :   a = RgM_shallowcopy(a);
    4065       78526 :   for (i=1; i<nbco; i++)
    4066             :   {
    4067       68670 :     k = pivot(a, data, i, NULL);
    4068       68670 :     if (k > nbco) return gerepilecopy(av, gcoeff(a,i,i));
    4069       68663 :     if (k != i)
    4070             :     { /* exchange the lines s.t. k = i */
    4071        4522 :       for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    4072        4522 :       s = -s;
    4073             :     }
    4074       68663 :     p = gcoeff(a,i,i);
    4075             : 
    4076       68663 :     x = gmul(x,p);
    4077      401633 :     for (k=i+1; k<=nbco; k++)
    4078             :     {
    4079      332970 :       GEN m = gcoeff(a,i,k);
    4080      332970 :       if (gequal0(m)) continue;
    4081             : 
    4082      116173 :       m = gdiv(m,p);
    4083      772623 :       for (j=i+1; j<=nbco; j++)
    4084      656450 :         gcoeff(a,j,k) = gsub(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
    4085             :     }
    4086       68663 :     if (gc_needed(av,2))
    4087             :     {
    4088           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
    4089           0 :       gerepileall(av,2, &a,&x);
    4090             :     }
    4091             :   }
    4092        9856 :   if (s < 0) x = gneg_i(x);
    4093        9856 :   return gerepileupto(av, gmul(x, gcoeff(a,nbco,nbco)));
    4094             : }
    4095             : 
    4096             : GEN
    4097        3280 : det2(GEN a)
    4098             : {
    4099             :   GEN data;
    4100             :   pivot_fun pivot;
    4101        3280 :   long n = lg(a)-1;
    4102        3280 :   if (typ(a)!=t_MAT) pari_err_TYPE("det2",a);
    4103        3280 :   if (!n) return gen_1;
    4104        3280 :   if (n != nbrows(a)) pari_err_DIM("det2");
    4105        3280 :   if (n == 1) return gcopy(gcoeff(a,1,1));
    4106        3280 :   if (n == 2) return RgM_det2(a);
    4107        1690 :   pivot = get_pivot_fun(a, a, &data);
    4108        1690 :   return det_simple_gauss(a, data, pivot);
    4109             : }
    4110             : 
    4111             : static GEN
    4112         672 : mydiv(GEN x, GEN y)
    4113             : {
    4114         672 :   long tx = typ(x), ty = typ(y);
    4115         672 :   if (tx == ty && tx == t_POL && varn(x) == varn(y)) return RgX_div(x,y);
    4116         672 :   return gdiv(x,y);
    4117             : }
    4118             : 
    4119             : /* Assumes a a square t_MAT of dimension n > 0. Returns det(a) using
    4120             :  * Gauss-Bareiss. */
    4121             : static GEN
    4122         196 : det_bareiss(GEN a)
    4123             : {
    4124         196 :   pari_sp av = avma;
    4125         196 :   long nbco = lg(a)-1,i,j,k,s = 1;
    4126             :   GEN p, pprec;
    4127             : 
    4128         196 :   a = RgM_shallowcopy(a);
    4129         770 :   for (pprec=gen_1,i=1; i<nbco; i++,pprec=p)
    4130             :   {
    4131             :     GEN ci;
    4132         574 :     int diveuc = (gequal1(pprec)==0);
    4133             : 
    4134         574 :     p = gcoeff(a,i,i);
    4135         574 :     if (gequal0(p))
    4136             :     {
    4137           0 :       k=i+1; while (k<=nbco && gequal0(gcoeff(a,i,k))) k++;
    4138           0 :       if (k>nbco) return gerepilecopy(av, p);
    4139           0 :       swap(gel(a,k), gel(a,i)); s = -s;
    4140           0 :       p = gcoeff(a,i,i);
    4141             :     }
    4142         574 :     ci = gel(a,i);
    4143        1708 :     for (k=i+1; k<=nbco; k++)
    4144             :     {
    4145        1134 :       GEN ck = gel(a,k), m = gel(ck,i);
    4146        1134 :       if (gequal0(m))
    4147             :       {
    4148           0 :         if (gequal1(p))
    4149             :         {
    4150           0 :           if (diveuc)
    4151           0 :             gel(a,k) = mydiv(gel(a,k), pprec);
    4152             :         }
    4153             :         else
    4154           0 :           for (j=i+1; j<=nbco; j++)
    4155             :           {
    4156           0 :             GEN p1 = gmul(p, gel(ck,j));
    4157           0 :             if (diveuc) p1 = mydiv(p1,pprec);
    4158           0 :             gel(ck,j) = p1;
    4159             :           }
    4160             :       }
    4161             :       else
    4162        3752 :         for (j=i+1; j<=nbco; j++)
    4163             :         {
    4164        2618 :           pari_sp av2 = avma;
    4165        2618 :           GEN p1 = gsub(gmul(p,gel(ck,j)), gmul(m,gel(ci,j)));
    4166        2618 :           if (diveuc) p1 = mydiv(p1,pprec);
    4167        2618 :           gel(ck,j) = gerepileupto(av2, p1);
    4168             :         }
    4169        1134 :       if (gc_needed(av,2))
    4170             :       {
    4171           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
    4172           0 :         gerepileall(av,2, &a,&pprec);
    4173           0 :         ci = gel(a,i);
    4174           0 :         p = gcoeff(a,i,i);
    4175             :       }
    4176             :     }
    4177             :   }
    4178         196 :   p = gcoeff(a,nbco,nbco);
    4179         196 :   p = (s < 0)? gneg(p): gcopy(p);
    4180         196 :   return gerepileupto(av, p);
    4181             : }
    4182             : 
    4183             : /* count non-zero entries in col j, at most 'max' of them.
    4184             :  * Return their indices */
    4185             : static GEN
    4186         609 : col_count_non_zero(GEN a, long j, long max)
    4187             : {
    4188         609 :   GEN v = cgetg(max+1, t_VECSMALL);
    4189         609 :   GEN c = gel(a,j);
    4190         609 :   long i, l = lg(a), k = 1;
    4191        2695 :   for (i = 1; i < l; i++)
    4192        2632 :     if (!gequal0(gel(c,i)))
    4193             :     {
    4194        2324 :       if (k > max) return NULL; /* fail */
    4195        1778 :       v[k++] = i;
    4196             :     }
    4197          63 :   setlg(v, k); return v;
    4198             : }
    4199             : /* count non-zero entries in row i, at most 'max' of them.
    4200             :  * Return their indices */
    4201             : static GEN
    4202         560 : row_count_non_zero(GEN a, long i, long max)
    4203             : {
    4204         560 :   GEN v = cgetg(max+1, t_VECSMALL);
    4205         560 :   long j, l = lg(a), k = 1;
    4206        2422 :   for (j = 1; j < l; j++)
    4207        2366 :     if (!gequal0(gcoeff(a,i,j)))
    4208             :     {
    4209        2212 :       if (k > max) return NULL; /* fail */
    4210        1708 :       v[k++] = j;
    4211             :     }
    4212          56 :   setlg(v, k); return v;
    4213             : }
    4214             : 
    4215             : static GEN det_develop(GEN a, long max, double bound);
    4216             : /* (-1)^(i+j) a[i,j] * det RgM_minor(a,i,j) */
    4217             : static GEN
    4218         203 : coeff_det(GEN a, long i, long j, long max, double bound)
    4219             : {
    4220         203 :   GEN c = gcoeff(a, i, j);
    4221         203 :   c = gmul(c, det_develop(RgM_minor(a, i,j), max, bound));
    4222         203 :   if (odd(i+j)) c = gneg(c);
    4223         203 :   return c;
    4224             : }
    4225             : /* a square t_MAT, 'bound' a rough upper bound for the number of
    4226             :  * multiplications we are willing to pay while developing rows/columns before
    4227             :  * switching to Gaussian elimination */
    4228             : static GEN
    4229         301 : det_develop(GEN M, long max, double bound)
    4230             : {
    4231         301 :   pari_sp av = avma;
    4232         301 :   long i,j, n = lg(M)-1, lbest = max+2, best_col = 0, best_row = 0;
    4233         301 :   GEN best = NULL;
    4234             : 
    4235         301 :   if (bound < 1.) return det_bareiss(M); /* too costly now */
    4236             : 
    4237         175 :   switch(n)
    4238             :   {
    4239           0 :     case 0: return gen_1;
    4240           0 :     case 1: return gcopy(gcoeff(M,1,1));
    4241          14 :     case 2: return RgM_det2(M);
    4242             :   }
    4243         161 :   if (max > ((n+2)>>1)) max = (n+2)>>1;
    4244         735 :   for (j = 1; j <= n; j++)
    4245             :   {
    4246         609 :     pari_sp av2 = avma;
    4247         609 :     GEN v = col_count_non_zero(M, j, max);
    4248             :     long lv;
    4249         609 :     if (!v || (lv = lg(v)) >= lbest) { avma = av2; continue; }
    4250          63 :     if (lv == 1) { avma = av; return gen_0; }
    4251          63 :     if (lv == 2) {
    4252          35 :       avma = av;
    4253          35 :       return gerepileupto(av, coeff_det(M,v[1],j,max,bound));
    4254             :     }
    4255          28 :     best = v; lbest = lv; best_col = j;
    4256             :   }
    4257         686 :   for (i = 1; i <= n; i++)
    4258             :   {
    4259         560 :     pari_sp av2 = avma;
    4260         560 :     GEN v = row_count_non_zero(M, i, max);
    4261             :     long lv;
    4262         560 :     if (!v || (lv = lg(v)) >= lbest) { avma = av2; continue; }
    4263          42 :     if (lv == 1) { avma = av; return gen_0; }
    4264          42 :     if (lv == 2) {
    4265           0 :       avma = av;
    4266           0 :       return gerepileupto(av, coeff_det(M,i,v[1],max,bound));
    4267             :     }
    4268          42 :     best = v; lbest = lv; best_row = i;
    4269             :   }
    4270         126 :   if (best_row)
    4271             :   {
    4272          42 :     double d = lbest-1;
    4273          42 :     GEN s = NULL;
    4274             :     long k;
    4275          42 :     bound /= d*d*d;
    4276         168 :     for (k = 1; k < lbest; k++)
    4277             :     {
    4278         126 :       GEN c = coeff_det(M, best_row, best[k], max, bound);
    4279         126 :       s = s? gadd(s, c): c;
    4280             :     }
    4281          42 :     return gerepileupto(av, s);
    4282             :   }
    4283          84 :   if (best_col)
    4284             :   {
    4285          14 :     double d = lbest-1;
    4286          14 :     GEN s = NULL;
    4287             :     long k;
    4288          14 :     bound /= d*d*d;
    4289          56 :     for (k = 1; k < lbest; k++)
    4290             :     {
    4291          42 :       GEN c = coeff_det(M, best[k], best_col, max, bound);
    4292          42 :       s = s? gadd(s, c): c;
    4293             :     }
    4294          14 :     return gerepileupto(av, s);
    4295             :   }
    4296          70 :   return det_bareiss(M);
    4297             : }
    4298             : 
    4299             : /* area of parallelogram bounded by (v1,v2) */
    4300             : static GEN
    4301       51464 : parallelogramarea(GEN v1, GEN v2)
    4302       51464 : { return gsub(gmul(gnorml2(v1), gnorml2(v2)), gsqr(RgV_dotproduct(v1, v2))); }
    4303             : 
    4304             : /* Square of Hadamard bound for det(a), a square matrix.
    4305             :  * Slightly improvement: instead of using the column norms, use the area of
    4306             :  * the parallelogram formed by pairs of consecutive vectors */
    4307             : GEN
    4308       16716 : RgM_Hadamard(GEN a)
    4309             : {
    4310       16716 :   pari_sp av = avma;
    4311       16716 :   long n = lg(a)-1, i;
    4312             :   GEN B;
    4313       16716 :   if (n == 0) return gen_1;
    4314       16716 :   if (n == 1) return gsqr(gcoeff(a,1,1));
    4315       16716 :   a = RgM_gtofp(a, LOWDEFAULTPREC);
    4316       16716 :   B = gen_1;
    4317       68180 :   for (i = 1; i <= n/2; i++)
    4318       51464 :     B = gmul(B, parallelogramarea(gel(a,2*i-1), gel(a,2*i)));
    4319       16716 :   if (odd(n)) B = gmul(B, gnorml2(gel(a, n)));
    4320       16716 :   return gerepileuptoint(av, ceil_safe(B));
    4321             : }
    4322             : 
    4323             : /* assume dim(a) = n > 0 */
    4324             : static GEN
    4325       27916 : ZM_det_i(GEN M, long n)
    4326             : {
    4327       27916 :   const long DIXON_THRESHOLD = 40;
    4328       27916 :   pari_sp av = avma, av2;
    4329             :   long i;
    4330       27916 :   ulong p, compp, Dp = 1;
    4331             :   forprime_t S;
    4332             :   GEN D, h, q, v, comp;
    4333       27916 :   if (n == 1) return icopy(gcoeff(M,1,1));
    4334       27097 :   if (n == 2) return ZM_det2(M);
    4335       19649 :   if (n == 3) return ZM_det3(M);
    4336       16716 :   h = RgM_Hadamard(M);
    4337       16716 :   if (!signe(h)) { avma = av; return gen_0; }
    4338       16716 :   h = sqrti(h); q = gen_1;
    4339       16716 :   init_modular_big(&S);
    4340       16716 :   p = 0; /* -Wall */
    4341       33432 :   while( cmpii(q, h) <= 0 && (p = u_forprime_next(&S)) )
    4342             :   {
    4343       16716 :     av2 = avma; Dp = Flm_det(ZM_to_Flm(M, p), p);
    4344       16716 :     avma = av2;
    4345       16716 :     if (Dp) break;
    4346           0 :     q = muliu(q, p);
    4347             :   }
    4348       16716 :   if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
    4349       16716 :   if (!Dp) { avma = av; return gen_0; }
    4350       16716 :   if (n <= DIXON_THRESHOLD)
    4351       16716 :     D = q;
    4352             :   else
    4353             :   {
    4354           0 :     av2 = avma;
    4355           0 :     v = cgetg(n+1, t_COL);
    4356           0 :     gel(v, 1) = gen_1; /* ensure content(v) = 1 */
    4357           0 :     for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
    4358           0 :     D = Q_denom(ZM_gauss(M, v));
    4359           0 :     if (expi(D) < expi(h) >> 1)
    4360             :     { /* First try unlucky, try once more */
    4361           0 :       for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
    4362           0 :       D = lcmii(D, Q_denom(ZM_gauss(M, v)));
    4363             :     }
    4364           0 :     D = gerepileuptoint(av2, D);
    4365           0 :     if (q != gen_1) D = lcmii(D, q);
    4366             :   }
    4367             :   /* determinant is a multiple of D */
    4368       16716 :   h = shifti(divii(h, D), 1);
    4369             : 
    4370       16716 :   compp = Fl_div(Dp, umodiu(D,p), p);
    4371       16716 :   comp = Z_init_CRT(compp, p);
    4372       16716 :   q = utoipos(p);
    4373       85128 :   while (cmpii(q, h) <= 0)
    4374             :   {
    4375       51696 :     p = u_forprime_next(&S);
    4376       51696 :     if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
    4377       51696 :     Dp = umodiu(D, p);
    4378       51696 :     if (!Dp) continue;
    4379       51696 :     av2 = avma;
    4380       51696 :     compp = Fl_div(Flm_det(ZM_to_Flm(M, p), p), Dp, p);
    4381       51696 :     avma = av2;
    4382       51696 :     (void) Z_incremental_CRT(&comp, compp, &q, p);
    4383             :   }
    4384       16716 :   return gerepileuptoint(av, mulii(comp, D));
    4385             : }
    4386             : 
    4387             : static long
    4388          98 : det_init_max(long n)
    4389             : {
    4390          98 :   if (n > 100) return 0;
    4391          98 :   if (n > 50) return 1;
    4392          98 :   if (n > 30) return 4;
    4393          98 :   return 7;
    4394             : }
    4395             : 
    4396             : GEN
    4397       11125 : det(GEN a)
    4398             : {
    4399       11125 :   long n = lg(a)-1;
    4400             :   double B;
    4401       11125 :   GEN data, ff = NULL, p = NULL;
    4402             :   pivot_fun pivot;
    4403             : 
    4404       11125 :   if (typ(a)!=t_MAT) pari_err_TYPE("det",a);
    4405       11125 :   if (!n) return gen_1;
    4406       11083 :   if (n != nbrows(a)) pari_err_DIM("det");
    4407       11076 :   if (n == 1) return gcopy(gcoeff(a,1,1));
    4408       10882 :   if (n == 2) return RgM_det2(a);
    4409       10000 :   if (RgM_is_FpM(a, &p))
    4410             :   {
    4411        1708 :     pari_sp av = avma;
    4412             :     ulong pp, d;
    4413        1708 :     if (!p) return ZM_det_i(a, n); /* ZM */
    4414             :     /* FpM */
    4415        1477 :     a = RgM_Fp_init(a,p,&pp);
    4416        1477 :     switch(pp)
    4417             :     {
    4418          49 :     case 0: return gerepileupto(av, Fp_to_mod(FpM_det(a,p),p)); break;
    4419          14 :     case 2: d = F2m_det(a); break;
    4420        1414 :     default:d = Flm_det(a,pp); break;
    4421             :     }
    4422        1428 :     avma = av; return mkintmodu(d, pp);
    4423             :   }
    4424        8292 :   if (RgM_is_FFM(a, &ff)) return FFM_det(a, ff);
    4425        8271 :   pivot = get_pivot_fun(a, a, &data);
    4426        8271 :   if (pivot != gauss_get_pivot_NZ) return det_simple_gauss(a, data, pivot);
    4427          98 :   B = (double)n;
    4428          98 :   return det_develop(a, det_init_max(n), B*B*B);
    4429             : }
    4430             : 
    4431             : GEN
    4432       27685 : ZM_det(GEN a)
    4433             : {
    4434       27685 :   long n = lg(a)-1;
    4435       27685 :   if (!n) return gen_1;
    4436       27685 :   return ZM_det_i(a, n);
    4437             : }
    4438             : 
    4439             : /* return a solution of congruence system sum M_{ i,j } X_j = Y_i mod D_i
    4440             :  * If ptu1 != NULL, put in *ptu1 a Z-basis of the homogeneous system */
    4441             : static GEN
    4442          49 : gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1)
    4443             : {
    4444          49 :   pari_sp av = avma;
    4445             :   long n, m, j, l, lM;
    4446             :   GEN delta, H, U, u1, u2, x;
    4447             : 
    4448          49 :   if (typ(M)!=t_MAT) pari_err_TYPE("gaussmodulo",M);
    4449          49 :   lM = lg(M);
    4450          49 :   if (lM == 1)
    4451             :   {
    4452          28 :     switch(typ(Y))
    4453             :     {
    4454           7 :       case t_INT: break;
    4455          14 :       case t_COL: if (lg(Y) != 1) pari_err_DIM("gaussmodulo");
    4456          14 :                   break;
    4457           7 :       default: pari_err_TYPE("gaussmodulo",Y);
    4458             :     }
    4459          21 :     switch(typ(D))
    4460             :     {
    4461           7 :       case t_INT: break;
    4462           7 :       case t_COL: if (lg(D) != 1) pari_err_DIM("gaussmodulo");
    4463           7 :                   break;
    4464           7 :       default: pari_err_TYPE("gaussmodulo",D);
    4465             :     }
    4466          14 :     if (ptu1) *ptu1 = cgetg(1, t_MAT);
    4467          14 :     return gen_0;
    4468             :   }
    4469          21 :   n = nbrows(M);
    4470          21 :   switch(typ(D))
    4471             :   {
    4472             :     case t_COL:
    4473          14 :       if (lg(D)-1!=n) pari_err_DIM("gaussmodulo");
    4474          14 :       delta = diagonal_shallow(D); break;
    4475           7 :     case t_INT: delta = scalarmat_shallow(D,n); break;
    4476           0 :     default: pari_err_TYPE("gaussmodulo",D);
    4477             :       return NULL; /* LCOV_EXCL_LINE */
    4478             :   }
    4479          21 :   switch(typ(Y))
    4480             :   {
    4481           7 :     case t_INT: Y = const_col(n, Y); break;
    4482             :     case t_COL:
    4483          14 :       if (lg(Y)-1!=n) pari_err_DIM("gaussmodulo");
    4484          14 :       break;
    4485           0 :     default: pari_err_TYPE("gaussmodulo",Y);
    4486             :       return NULL; /* LCOV_EXCL_LINE */
    4487             :   }
    4488          21 :   H = ZM_hnfall_i(shallowconcat(M,delta), &U, 1);
    4489          21 :   Y = hnf_solve(H,Y); if (!Y) return gen_0;
    4490          21 :   l = lg(H); /* may be smaller than lM if some moduli are 0 */
    4491          21 :   n = l-1;
    4492          21 :   m = lg(U)-1 - n;
    4493          21 :   u1 = cgetg(m+1,t_MAT);
    4494          21 :   u2 = cgetg(n+1,t_MAT);
    4495          21 :   for (j=1; j<=m; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u1,j) = c; }
    4496          21 :   U += m;
    4497          21 :   for (j=1; j<=n; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u2,j) = c; }
    4498             :   /*       (u1 u2)
    4499             :    * (M D) (*  * ) = (0 H) */
    4500          21 :   u1 = ZM_lll(u1, 0.75, LLL_INPLACE);
    4501          21 :   Y = ZM_ZC_mul(u2,Y);
    4502          21 :   x = ZC_reducemodmatrix(Y, u1);
    4503          21 :   if (!ptu1) x = gerepileupto(av, x);
    4504             :   else
    4505             :   {
    4506          14 :     gerepileall(av, 2, &x, &u1);
    4507          14 :     *ptu1 = u1;
    4508             :   }
    4509          21 :   return x;
    4510             : }
    4511             : 
    4512             : GEN
    4513          49 : matsolvemod0(GEN M, GEN D, GEN Y, long flag)
    4514             : {
    4515             :   pari_sp av;
    4516             :   GEN p1,y;
    4517             : 
    4518          49 :   if (!flag) return gaussmoduloall(M,D,Y,NULL);
    4519             : 
    4520          42 :   av=avma; y = cgetg(3,t_VEC);
    4521          42 :   p1 = gaussmoduloall(M,D,Y, (GEN*)y+2);
    4522          28 :   if (p1==gen_0) { avma=av; return gen_0; }
    4523          14 :   gel(y,1) = p1; return y;
    4524             : }
    4525             : 
    4526             : GEN
    4527           0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod0(M,D,Y,1); }
    4528             : 
    4529             : GEN
    4530           0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod0(M,D,Y,0); }

Generated by: LCOV version 1.11