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 20291-5fbfea9) Lines: 2398 2682 89.4 %
Date: 2017-02-25 05:49:34 Functions: 226 248 91.1 %
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          68 : gerepile_gauss(GEN x,long k,long t,pari_sp av, long j, GEN c)
      85             : {
      86          68 :   pari_sp tetpil = avma, A, bot;
      87          68 :   long u,i, n = lg(x)-1, m = n? nbrows(x): 0;
      88             :   size_t dec;
      89             : 
      90          68 :   if (DEBUGMEM > 1) pari_warn(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n);
      91        1751 :   for (u=t+1; u<=m; u++)
      92        1683 :     if (u==j || !c[u]) COPY(gcoeff(x,u,k));
      93        3544 :   for (u=1; u<=m; u++)
      94        3476 :     if (u==j || !c[u])
      95        2581 :       for (i=k+1; i<=n; i++) COPY(gcoeff(x,u,i));
      96             : 
      97          68 :   (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
      98          68 :   bot = pari_mainstack->bot;
      99        1751 :   for (u=t+1; u<=m; u++)
     100        1683 :     if (u==j || !c[u])
     101             :     {
     102        1683 :       A=(pari_sp)coeff(x,u,k);
     103        1683 :       if (A<av && A>=bot) coeff(x,u,k)+=dec;
     104             :     }
     105        3544 :   for (u=1; u<=m; u++)
     106        3476 :     if (u==j || !c[u])
     107      113660 :       for (i=k+1; i<=n; i++)
     108             :       {
     109      111079 :         A=(pari_sp)coeff(x,u,i);
     110      111079 :         if (A<av && A>=bot) coeff(x,u,i)+=dec;
     111             :       }
     112          68 : }
     113             : 
     114             : /*******************************************************************/
     115             : /*                                                                 */
     116             : /*                         GENERIC                                 */
     117             : /*                                                                 */
     118             : /*******************************************************************/
     119             : GEN
     120        2361 : gen_ker(GEN x, long deplin, void *E, const struct bb_field *ff)
     121             : {
     122        2361 :   pari_sp av0 = avma, av, tetpil;
     123             :   GEN y, c, d;
     124             :   long i, j, k, r, t, n, m;
     125             : 
     126        2361 :   n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
     127        2361 :   m=nbrows(x); r=0;
     128        2361 :   x = RgM_shallowcopy(x);
     129        2361 :   c = zero_zv(m);
     130        2361 :   d=new_chunk(n+1);
     131        2361 :   av=avma;
     132       30209 :   for (k=1; k<=n; k++)
     133             :   {
     134      745350 :     for (j=1; j<=m; j++)
     135      730449 :       if (!c[j])
     136             :       {
     137      478695 :         gcoeff(x,j,k) = ff->red(E, gcoeff(x,j,k));
     138      478695 :         if (!ff->equal0(gcoeff(x,j,k))) break;
     139             :       }
     140       27876 :     if (j>m)
     141             :     {
     142       14901 :       if (deplin)
     143             :       {
     144          28 :         GEN c = cgetg(n+1, t_COL), g0 = ff->s(E,0), g1=ff->s(E,1);
     145          28 :         for (i=1; i<k; i++) gel(c,i) = ff->red(E, gcoeff(x,d[i],k));
     146          28 :         gel(c,k) = g1; for (i=k+1; i<=n; i++) gel(c,i) = g0;
     147          28 :         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       12975 :       GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
     156       12975 :       c[j] = k; d[k] = j;
     157       12975 :       gcoeff(x,j,k) = ff->s(E,-1);
     158       12975 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
     159      530995 :       for (t=1; t<=m; t++)
     160             :       {
     161      518020 :         if (t==j) continue;
     162             : 
     163      505045 :         piv = ff->red(E,gcoeff(x,t,k));
     164      505045 :         if (ff->equal0(piv)) continue;
     165             : 
     166       85429 :         gcoeff(x,t,k) = ff->s(E,0);
     167     1491236 :         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       85429 :         if (gc_needed(av,1))
     170          10 :           gen_gerepile_gauss_ker(x,k,t,av,E,ff->red);
     171             :       }
     172             :     }
     173             :   }
     174        2333 :   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         423 : gen_Gauss_pivot(GEN x, long *rr, void *E, const struct bb_field *ff)
     197             : {
     198             :   pari_sp av;
     199             :   GEN c, d;
     200         423 :   long i, j, k, r, t, m, n = lg(x)-1;
     201             : 
     202         423 :   if (!n) { *rr = 0; return NULL; }
     203             : 
     204         423 :   m=nbrows(x); r=0;
     205         423 :   d = cgetg(n+1, t_VECSMALL);
     206         423 :   x = RgM_shallowcopy(x);
     207         423 :   c = zero_zv(m);
     208         423 :   av=avma;
     209        5489 :   for (k=1; k<=n; k++)
     210             :   {
     211       94903 :     for (j=1; j<=m; j++)
     212       93782 :       if (!c[j])
     213             :       {
     214       42298 :         gcoeff(x,j,k) = ff->red(E,gcoeff(x,j,k));
     215       42298 :         if (!ff->equal0(gcoeff(x,j,k))) break;
     216             :       }
     217        5066 :     if (j>m) { r++; d[k]=0; }
     218             :     else
     219             :     {
     220        3945 :       GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
     221        3945 :       GEN g0 = ff->s(E,0);
     222        3945 :       c[j] = k; d[k] = j;
     223        3945 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
     224      106195 :       for (t=1; t<=m; t++)
     225             :       {
     226      102250 :         if (c[t]) continue; /* already a pivot on that line */
     227             : 
     228       55377 :         piv = ff->red(E,gcoeff(x,t,k));
     229       55377 :         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          68 :           gerepile_gauss(x,k,t,av,j,c);
     235             :       }
     236        3945 :       for (i=k; i<=n; i++) gcoeff(x,j,i) = g0; /* dummy */
     237             :     }
     238             :   }
     239         423 :   *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       89807 : image_from_pivot(GEN x, GEN d, long r)
     515             : {
     516             :   GEN y;
     517             :   long j, k;
     518             : 
     519       89807 :   if (!d) return gcopy(x);
     520             :   /* d left on stack for efficiency */
     521       88379 :   r = lg(x)-1 - r; /* = dim Im(x) */
     522       88379 :   y = cgetg(r+1,t_MAT);
     523      739097 :   for (j=k=1; j<=r; k++)
     524      650718 :     if (d[k]) gel(y,j++) = gcopy(gel(x,k));
     525       88379 :   return y;
     526             : }
     527             : 
     528             : /*******************************************************************/
     529             : /*                                                                 */
     530             : /*                    LINEAR ALGEBRA MODULO P                      */
     531             : /*                                                                 */
     532             : /*******************************************************************/
     533             : 
     534             : static long
     535     1940528 : F2v_find_nonzero(GEN x0, GEN mask0, long l, long m)
     536             : {
     537     1940528 :   ulong *x = (ulong *)x0+2, *mask = (ulong *)mask0+2, e;
     538             :   long i, j;
     539     3528254 :   for (i = 0; i < l; i++)
     540             :   {
     541     3526751 :     e = *x++ & *mask++;
     542     3526751 :     if (e)
     543     1939025 :       for (j = 1; ; j++, e >>= 1) if (e & 1uL) return i*BITS_IN_LONG+j;
     544     8193040 :   }
     545        1503 :   return m+1;
     546             : }
     547             : 
     548             : /* in place, destroy x */
     549             : GEN
     550      288517 : F2m_ker_sp(GEN x, long deplin)
     551             : {
     552             :   GEN y, c, d;
     553             :   long i, j, k, l, r, m, n;
     554             : 
     555      288517 :   n = lg(x)-1;
     556      288517 :   m = mael(x,1,1); r=0;
     557             : 
     558      288517 :   d = cgetg(n+1, t_VECSMALL);
     559      288517 :   c = const_F2v(m); l = lg(c)-1;
     560     2080302 :   for (k=1; k<=n; k++)
     561             :   {
     562     1810594 :     GEN xk = gel(x,k);
     563     1810594 :     j = F2v_find_nonzero(xk, c, l, m);
     564     1810594 :     if (j>m)
     565             :     {
     566      395046 :       if (deplin) {
     567       18809 :         GEN c = zero_F2v(n);
     568       67233 :         for (i=1; i<k; i++)
     569       48424 :           if (F2v_coeff(xk, d[i]))
     570       24920 :             F2v_set(c, i);
     571       18809 :         F2v_set(c, k);
     572       18809 :         return c;
     573             :       }
     574      376237 :       r++; d[k] = 0;
     575             :     }
     576             :     else
     577             :     {
     578     1415548 :       F2v_clear(c,j); d[k] = j;
     579     1415548 :       F2v_clear(xk, j);
     580    76652094 :       for (i=k+1; i<=n; i++)
     581             :       {
     582    75236546 :         GEN xi = gel(x,i);
     583    75236546 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     584             :       }
     585     1415548 :       F2v_set(xk, j);
     586             :     }
     587             :   }
     588      269708 :   if (deplin) return NULL;
     589             : 
     590      269673 :   y = zero_F2m_copy(n,r);
     591      645910 :   for (j=k=1; j<=r; j++,k++)
     592             :   {
     593      376237 :     GEN C = gel(y,j); while (d[k]) k++;
     594    13862540 :     for (i=1; i<k; i++)
     595    13486303 :       if (d[i] && F2m_coeff(x,d[i],k))
     596     4967611 :         F2v_set(C,i);
     597      376237 :     F2v_set(C, k);
     598             :   }
     599      269673 :   return y;
     600             : }
     601             : 
     602             : static void /* assume m < p */
     603   143491313 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
     604             : {
     605   143491313 :   uel(b,k) = Fl_addmul_pre(m, uel(b,i), uel(b,k), p, pi);
     606   143491313 : }
     607             : static void /* same m = 1 */
     608     3169596 : _Fl_add(GEN b, long k, long i, ulong p)
     609             : {
     610     3169596 :   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
     611     3169596 : }
     612             : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     613   251474473 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
     614             : {
     615   251474473 :   uel(b,k) += m * uel(b,i);
     616   251474473 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     617   251474473 : }
     618             : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     619    78438125 : _Fl_add_OK(GEN b, long k, long i, ulong p)
     620             : {
     621    78438125 :   uel(b,k) += uel(b,i);
     622    78438125 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     623    78438125 : }
     624             : 
     625             : static GEN
     626      493945 : 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      493945 :   n = lg(x)-1;
     633      493945 :   m=nbrows(x); r=0;
     634             : 
     635      493945 :   c = zero_zv(m);
     636      493945 :   d = cgetg(n+1, t_VECSMALL);
     637      493945 :   a = 0; /* for gcc -Wall */
     638     3839348 :   for (k=1; k<=n; k++)
     639             :   {
     640    30630376 :     for (j=1; j<=m; j++)
     641    29422452 :       if (!c[j])
     642             :       {
     643    12034165 :         a = ucoeff(x,j,k) % p;
     644    12034165 :         if (a) break;
     645             :       }
     646     3385608 :     if (j > m)
     647             :     {
     648     1207924 :       if (deplin==1) {
     649       40205 :         c = cgetg(n+1, t_VECSMALL);
     650       40205 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
     651       40205 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     652       40205 :         return c;
     653             :       }
     654     1167719 :       r++; d[k] = 0;
     655             :     }
     656             :     else
     657             :     {
     658     2177684 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     659     2177684 :       c[j] = k; d[k] = j;
     660     2177684 :       ucoeff(x,j,k) = p-1;
     661     2177684 :       if (piv != 1)
     662     1747031 :         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
     663    61920334 :       for (t=1; t<=m; t++)
     664             :       {
     665    59742650 :         if (t == j) continue;
     666             : 
     667    57564966 :         piv = ( ucoeff(x,t,k) %= p );
     668    57564966 :         if (!piv) continue;
     669    20557825 :         if (piv == 1)
     670     4265155 :           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
     671             :         else
     672    16292670 :           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
     673             :       }
     674             :     }
     675             :   }
     676      453740 :   if (deplin==1) return NULL;
     677             : 
     678      453733 :   y = cgetg(r+1, t_MAT);
     679     1621452 :   for (j=k=1; j<=r; j++,k++)
     680             :   {
     681     1167719 :     GEN C = cgetg(n+1, t_VECSMALL);
     682             : 
     683     1167719 :     gel(y,j) = C; while (d[k]) k++;
     684     9940297 :     for (i=1; i<k; i++)
     685     8772578 :       if (d[i])
     686     6109305 :         uel(C,i) = ucoeff(x,d[i],k) % p;
     687             :       else
     688     2663273 :         uel(C,i) = 0UL;
     689     1167719 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     690             :   }
     691      453733 :   return deplin==2? mkvec2(y, d): y;
     692             : }
     693             : 
     694             : /* in place, destroy x */
     695             : GEN
     696      501871 : 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      501871 :   n = lg(x)-1;
     702      501871 :   if (!n) return cgetg(1,t_MAT);
     703      501871 :   if (SMALL_ULONG(p)) return Flm_ker_sp_OK(x, p, deplin);
     704        7926 :   pi = get_Fl_red(p);
     705             : 
     706        7926 :   m=nbrows(x); r=0;
     707             : 
     708        7926 :   c = zero_zv(m);
     709        7926 :   d = cgetg(n+1, t_VECSMALL);
     710        7926 :   a = 0; /* for gcc -Wall */
     711       87429 :   for (k=1; k<=n; k++)
     712             :   {
     713     2104383 :     for (j=1; j<=m; j++)
     714     2071863 :       if (!c[j])
     715             :       {
     716     1044915 :         a = ucoeff(x,j,k);
     717     1044915 :         if (a) break;
     718             :       }
     719       79510 :     if (j > m)
     720             :     {
     721       32520 :       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       32513 :       r++; d[k] = 0;
     728             :     }
     729             :     else
     730             :     {
     731       46990 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     732       46990 :       c[j] = k; d[k] = j;
     733       46990 :       ucoeff(x,j,k) = p-1;
     734       46990 :       if (piv != 1)
     735      759629 :         for (i=k+1; i<=n; i++)
     736      716242 :           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
     737     1895170 :       for (t=1; t<=m; t++)
     738             :       {
     739     1848180 :         if (t == j) continue;
     740             : 
     741     1801190 :         piv = ucoeff(x,t,k);
     742     1801190 :         if (!piv) continue;
     743      631930 :         if (piv == 1)
     744        8119 :           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
     745             :         else
     746      623811 :           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
     747             :       }
     748             :     }
     749             :   }
     750        7919 :   if (deplin==1) return NULL;
     751             : 
     752        7912 :   y = cgetg(r+1, t_MAT);
     753       40425 :   for (j=k=1; j<=r; j++,k++)
     754             :   {
     755       32513 :     GEN C = cgetg(n+1, t_VECSMALL);
     756             : 
     757       32513 :     gel(y,j) = C; while (d[k]) k++;
     758      966536 :     for (i=1; i<k; i++)
     759      934023 :       if (d[i])
     760      589777 :         uel(C,i) = ucoeff(x,d[i],k);
     761             :       else
     762      344246 :         uel(C,i) = 0UL;
     763       32513 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     764             :   }
     765        7912 :   return deplin==2? mkvec2(y, d): y;
     766             : }
     767             : 
     768             : GEN
     769       70287 : FpM_intersect(GEN x, GEN y, GEN p)
     770             : {
     771       70287 :   pari_sp av = avma;
     772       70287 :   long j, lx = lg(x);
     773             :   GEN z;
     774             : 
     775       70287 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     776       70287 :   z = FpM_ker(shallowconcat(x,y), p);
     777       70287 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     778       70287 :   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       55252 : 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      408748 : FpM_init(GEN a, GEN p, ulong *pp)
     909             : {
     910      408748 :   if (lgefint(p) == 3)
     911             :   {
     912      404214 :     *pp = uel(p,2);
     913      404214 :     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       23434 : F2m_gauss_pivot(GEN x, long *rr)
     953             : {
     954             :   GEN c, d;
     955             :   long i, j, k, l, r, m, n;
     956             : 
     957       23434 :   n = lg(x)-1; if (!n) { *rr=0; return NULL; }
     958       23434 :   m = mael(x,1,1); r=0;
     959             : 
     960       23434 :   d = cgetg(n+1, t_VECSMALL);
     961       23434 :   c = const_F2v(m); l = lg(c)-1;
     962      153368 :   for (k=1; k<=n; k++)
     963             :   {
     964      129934 :     GEN xk = gel(x,k);
     965      129934 :     j = F2v_find_nonzero(xk, c, l, m);
     966      129934 :     if (j>m) { r++; d[k] = 0; }
     967             :     else
     968             :     {
     969       83892 :       F2v_clear(c,j); d[k] = j;
     970     1427215 :       for (i=k+1; i<=n; i++)
     971             :       {
     972     1343323 :         GEN xi = gel(x,i);
     973     1343323 :         if (F2v_coeff(xi,j)) F2v_add_inplace(xi, xk);
     974             :       }
     975             :     }
     976             :   }
     977             : 
     978       23434 :   *rr = r; avma = (pari_sp)d; return d;
     979             : }
     980             : 
     981             : /* Destroy x */
     982             : static GEN
     983      210272 : 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      210272 :   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
     989             : 
     990      210272 :   m=nbrows(x); r=0;
     991      210272 :   d=cgetg(n+1,t_VECSMALL);
     992      210272 :   c = zero_zv(m);
     993     1848726 :   for (k=1; k<=n; k++)
     994             :   {
     995    17012454 :     for (j=1; j<=m; j++)
     996    16463184 :       if (!c[j])
     997             :       {
     998     4948453 :         ucoeff(x,j,k) %= p;
     999     4948453 :         if (ucoeff(x,j,k)) break;
    1000             :       }
    1001     1638454 :     if (j>m) { r++; d[k]=0; }
    1002             :     else
    1003             :     {
    1004     1089184 :       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
    1005     1089184 :       c[j]=k; d[k]=j;
    1006    13357082 :       for (i=k+1; i<=n; i++)
    1007    12267898 :         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
    1008    19628950 :       for (t=1; t<=m; t++)
    1009    18539766 :         if (!c[t]) /* no pivot on that line yet */
    1010             :         {
    1011    10956882 :           piv = ucoeff(x,t,k);
    1012    10956882 :           if (piv)
    1013             :           {
    1014     5355264 :             ucoeff(x,t,k) = 0;
    1015   223828095 :             for (i=k+1; i<=n; i++)
    1016   436945662 :               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
    1017   218472831 :                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
    1018             :           }
    1019             :         }
    1020     1089184 :       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
    1021             :     }
    1022             :   }
    1023      210272 :   *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      119645 : FpM_gauss_pivot(GEN x, GEN p, long *rr)
    1035             : {
    1036             :   ulong pp;
    1037      119645 :   if (lg(x)==1) { *rr = 0; return NULL; }
    1038      118252 :   x = FpM_init(x, p, &pp);
    1039      118252 :   switch(pp)
    1040             :   {
    1041          59 :   case 0: return FpM_gauss_pivot_gen(x, p, rr);
    1042       23350 :   case 2: return F2m_gauss_pivot(x, rr);
    1043       94843 :   default:return Flm_gauss_pivot(x, pp, rr);
    1044             :   }
    1045             : }
    1046             : 
    1047             : GEN
    1048       86860 : FpM_image(GEN x, GEN p)
    1049             : {
    1050             :   long r;
    1051       86860 :   GEN d = FpM_gauss_pivot(x,p,&r); /* d left on stack for efficiency */
    1052       86860 :   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        2184 : Flm_rank(GEN x, ulong p)
    1079             : {
    1080        2184 :   pari_sp av = avma;
    1081             :   long r;
    1082        2184 :   (void)Flm_gauss_pivot(Flm_copy(x),p,&r);
    1083        2184 :   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         258 : FlxqM_gauss_pivot(GEN x, GEN T, ulong p, long *rr)
    1097             : {
    1098             :   void *E;
    1099         258 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1100         258 :   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          78 : FqM_gauss_pivot_gen(GEN x, GEN T, GEN p, long *rr)
    1242             : {
    1243             :   void *E;
    1244          78 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1245          78 :   return gen_Gauss_pivot(x, rr, E, S);
    1246             : }
    1247             : static GEN
    1248         176 : FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr)
    1249             : {
    1250         176 :   if (lg(x)==1) { *rr = 0; return NULL; }
    1251         176 :   if (!T) return FpM_gauss_pivot(x, p, rr);
    1252         176 :   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          78 :   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      248192 : FpM_ker_i(GEN x, GEN p, long deplin)
    1325             : {
    1326      248192 :   pari_sp av = avma;
    1327             :   ulong pp;
    1328             :   GEN y;
    1329             : 
    1330      248192 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1331      248087 :   x = FpM_init(x, p, &pp);
    1332      248087 :   switch(pp)
    1333             :   {
    1334        1619 :   case 0: return FpM_ker_gen(x,p,deplin);
    1335             :   case 2:
    1336       56840 :     y = F2m_ker_sp(x, deplin);
    1337       56840 :     if (!y) return y;
    1338       56840 :     y = deplin? F2c_to_ZC(y): F2m_to_ZM(y);
    1339       56840 :     return gerepileupto(av, y);
    1340             :   default:
    1341      189628 :     y = Flm_ker_sp(x, pp, deplin);
    1342      189628 :     if (!y) return y;
    1343      189628 :     y = deplin? Flc_to_ZC(y): Flm_to_ZM(y);
    1344      189628 :     return gerepileupto(av, y);
    1345             :   }
    1346             : }
    1347             : 
    1348             : GEN
    1349      188191 : FpM_ker(GEN x, GEN p) { return FpM_ker_i(x,p,0); }
    1350             : 
    1351             : GEN
    1352       58993 : FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x,p,1); }
    1353             : 
    1354             : static GEN
    1355          29 : FqM_ker_gen(GEN x, GEN T, GEN p, long deplin)
    1356             : {
    1357             :   void *E;
    1358          29 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    1359          29 :   return gen_ker(x,deplin,E,S);
    1360             : }
    1361             : static GEN
    1362        1268 : FqM_ker_i(GEN x, GEN T, GEN p, long deplin)
    1363             : {
    1364        1268 :   if (!T) return FpM_ker_i(x,p,deplin);
    1365         260 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1366             : 
    1367         260 :   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          29 :   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          14 : FqM_deplin(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,1); }
    1384             : 
    1385             : static GEN
    1386         566 : FlxqM_ker_i(GEN x, GEN T, ulong p, long deplin)
    1387             : {
    1388             :   const struct bb_field *ff;
    1389             :   void *E;
    1390             : 
    1391         566 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1392         566 :   ff=get_Flxq_field(&E,T,p);
    1393         566 :   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             : GEN
    1403          14 : FlxqM_deplin(GEN x, GEN T, ulong p)
    1404             : {
    1405          14 :   return FlxqM_ker_i(x, T, p, 1);
    1406             : }
    1407             : 
    1408             : static GEN
    1409          35 : F2xqM_ker_i(GEN x, GEN T, long deplin)
    1410             : {
    1411             :   const struct bb_field *ff;
    1412             :   void *E;
    1413             : 
    1414          35 :   if (lg(x)==1) return cgetg(1,t_MAT);
    1415          35 :   ff = get_F2xq_field(&E,T);
    1416          35 :   return gen_ker(x,deplin, E, ff);
    1417             : }
    1418             : 
    1419             : GEN
    1420          21 : F2xqM_ker(GEN x, GEN T)
    1421             : {
    1422          21 :   return F2xqM_ker_i(x, T, 0);
    1423             : }
    1424             : 
    1425             : GEN
    1426          14 : F2xqM_deplin(GEN x, GEN T)
    1427             : {
    1428          14 :   return F2xqM_ker_i(x, T, 1);
    1429             : }
    1430             : 
    1431             : static GEN
    1432          28 : F2xqM_gauss_pivot(GEN x, GEN T, long *rr)
    1433             : {
    1434             :   void *E;
    1435          28 :   const struct bb_field *S = get_F2xq_field(&E,T);
    1436          28 :   return gen_Gauss_pivot(x, rr, E, S);
    1437             : }
    1438             : GEN
    1439           7 : F2xqM_image(GEN x, GEN T)
    1440             : {
    1441             :   long r;
    1442           7 :   GEN d = F2xqM_gauss_pivot(x,T,&r); /* d left on stack for efficiency */
    1443           7 :   return image_from_pivot(x,d,r);
    1444             : }
    1445             : long
    1446           7 : F2xqM_rank(GEN x, GEN T)
    1447             : {
    1448           7 :   pari_sp av = avma;
    1449             :   long r;
    1450           7 :   (void)F2xqM_gauss_pivot(x,T,&r);
    1451           7 :   avma = av; return lg(x)-1 - r;
    1452             : }
    1453             : /*******************************************************************/
    1454             : /*                                                                 */
    1455             : /*                       Solve A*X=B (Gauss pivot)                 */
    1456             : /*                                                                 */
    1457             : /*******************************************************************/
    1458             : /* x ~ 0 compared to reference y */
    1459             : int
    1460      653274 : approx_0(GEN x, GEN y)
    1461             : {
    1462      653274 :   long tx = typ(x);
    1463      653274 :   if (tx == t_COMPLEX)
    1464           0 :     return approx_0(gel(x,1), y) && approx_0(gel(x,2), y);
    1465      653414 :   return gequal0(x) ||
    1466      457139 :          (tx == t_REAL && gexpo(y) - gexpo(x) > bit_prec(x));
    1467             : }
    1468             : /* x a column, x0 same column in the original input matrix (for reference),
    1469             :  * c list of pivots so far */
    1470             : static long
    1471      657418 : gauss_get_pivot_max(GEN X, GEN X0, long ix, GEN c)
    1472             : {
    1473      657418 :   GEN p, r, x = gel(X,ix), x0 = gel(X0,ix);
    1474      657418 :   long i, k = 0, ex = - (long)HIGHEXPOBIT, lx = lg(x);
    1475      657418 :   if (c)
    1476             :   {
    1477       46550 :     for (i=1; i<lx; i++)
    1478       35805 :       if (!c[i])
    1479             :       {
    1480       18634 :         long e = gexpo(gel(x,i));
    1481       18634 :         if (e > ex) { ex = e; k = i; }
    1482             :       }
    1483             :   }
    1484             :   else
    1485             :   {
    1486     2228868 :     for (i=ix; i<lx; i++)
    1487             :     {
    1488     1582195 :       long e = gexpo(gel(x,i));
    1489     1582195 :       if (e > ex) { ex = e; k = i; }
    1490             :     }
    1491             :   }
    1492      657418 :   if (!k) return lx;
    1493      653274 :   p = gel(x,k);
    1494      653274 :   r = gel(x0,k); if (isrationalzero(r)) r = x0;
    1495      653274 :   return approx_0(p, r)? lx: k;
    1496             : }
    1497             : static long
    1498       63091 : gauss_get_pivot_padic(GEN X, GEN p, long ix, GEN c)
    1499             : {
    1500       63091 :   GEN x = gel(X, ix);
    1501       63091 :   long i, k = 0, ex = (long)HIGHVALPBIT, lx = lg(x);
    1502       63091 :   if (c)
    1503             :   {
    1504         504 :     for (i=1; i<lx; i++)
    1505         378 :       if (!c[i] && !gequal0(gel(x,i)))
    1506             :       {
    1507         245 :         long e = gvaluation(gel(x,i), p);
    1508         245 :         if (e < ex) { ex = e; k = i; }
    1509             :       }
    1510             :   }
    1511             :   else
    1512             :   {
    1513      443905 :     for (i=ix; i<lx; i++)
    1514      380940 :       if (!gequal0(gel(x,i)))
    1515             :       {
    1516      182434 :         long e = gvaluation(gel(x,i), p);
    1517      182434 :         if (e < ex) { ex = e; k = i; }
    1518             :       }
    1519             :   }
    1520       63091 :   return k? k: lx;
    1521             : }
    1522             : static long
    1523      294162 : gauss_get_pivot_NZ(GEN X, GEN x0/*unused*/, long ix, GEN c)
    1524             : {
    1525      294162 :   GEN x = gel(X, ix);
    1526      294162 :   long i, lx = lg(x);
    1527             :   (void)x0;
    1528      294162 :   if (c)
    1529             :   {
    1530     1544249 :     for (i=1; i<lx; i++)
    1531     1523214 :       if (!c[i] && !gequal0(gel(x,i))) return i;
    1532             :   }
    1533             :   else
    1534             :   {
    1535      290404 :     for (i=ix; i<lx; i++)
    1536      290397 :       if (!gequal0(gel(x,i))) return i;
    1537             :   }
    1538       21042 :   return lx;
    1539             : }
    1540             : 
    1541             : /* Return pivot seeking function appropriate for the domain of the RgM x
    1542             :  * (first non zero pivot, maximal pivot...)
    1543             :  * x0 is a reference point used when guessing whether x[i,j] ~ 0
    1544             :  * (iff x[i,j] << x0[i,j]); typical case: mateigen, Gauss pivot on x - vp.Id,
    1545             :  * but use original x when deciding whether a prospective pivot is non-0 */
    1546             : static pivot_fun
    1547      272176 : get_pivot_fun(GEN x, GEN x0, GEN *data)
    1548             : {
    1549      272176 :   long i, j, hx, lx = lg(x);
    1550      272176 :   int res = t_INT;
    1551      272176 :   GEN p = NULL;
    1552             : 
    1553      272176 :   *data = NULL;
    1554      272176 :   if (lx == 1) return &gauss_get_pivot_NZ;
    1555      272036 :   hx = lgcols(x);
    1556     1305215 :   for (j=1; j<lx; j++)
    1557             :   {
    1558     1033207 :     GEN xj = gel(x,j);
    1559     9205724 :     for (i=1; i<hx; i++)
    1560             :     {
    1561     8172545 :       GEN c = gel(xj,i);
    1562     8172545 :       switch(typ(c))
    1563             :       {
    1564             :         case t_REAL:
    1565     1820378 :           res = t_REAL;
    1566     1820378 :           break;
    1567             :         case t_COMPLEX:
    1568          14 :           if (typ(gel(c,1)) == t_REAL || typ(gel(c,2)) == t_REAL) res = t_REAL;
    1569          14 :           break;
    1570             :         case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT: case t_QUAD:
    1571             :         case t_POLMOD: /* exact types */
    1572     6192686 :           break;
    1573             :         case t_PADIC:
    1574      159439 :           p = gel(c,2);
    1575      159439 :           res = t_PADIC;
    1576      159439 :           break;
    1577          28 :         default: return &gauss_get_pivot_NZ;
    1578             :       }
    1579             :     }
    1580             :   }
    1581      272008 :   switch(res)
    1582             :   {
    1583      191454 :     case t_REAL: *data = x0; return &gauss_get_pivot_max;
    1584        8225 :     case t_PADIC: *data = p; return &gauss_get_pivot_padic;
    1585       72329 :     default: return &gauss_get_pivot_NZ;
    1586             :   }
    1587             : }
    1588             : 
    1589             : static GEN
    1590      488833 : get_col(GEN a, GEN b, GEN p, long li)
    1591             : {
    1592      488833 :   GEN u = cgetg(li+1,t_COL);
    1593             :   long i, j;
    1594             : 
    1595      488833 :   gel(u,li) = gdiv(gel(b,li), p);
    1596     2713989 :   for (i=li-1; i>0; i--)
    1597             :   {
    1598     2225156 :     pari_sp av = avma;
    1599     2225156 :     GEN m = gel(b,i);
    1600     2225156 :     for (j=i+1; j<=li; j++) m = gsub(m, gmul(gcoeff(a,i,j), gel(u,j)));
    1601     2225156 :     gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(a,i,i)));
    1602             :   }
    1603      488833 :   return u;
    1604             : }
    1605             : /* assume 0 <= a[i,j] < p */
    1606             : static GEN
    1607      427860 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
    1608             : {
    1609      427860 :   GEN u = cgetg(li+1,t_VECSMALL);
    1610      427860 :   ulong m = uel(b,li) % p;
    1611             :   long i,j;
    1612             : 
    1613      427860 :   uel(u,li) = (m * ucoeff(a,li,li)) % p;
    1614     6271869 :   for (i = li-1; i > 0; i--)
    1615             :   {
    1616     5844009 :     m = p - uel(b,i)%p;
    1617    67973071 :     for (j = i+1; j <= li; j++) {
    1618    62129062 :       if (m & HIGHBIT) m %= p;
    1619    62129062 :       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
    1620             :     }
    1621     5844009 :     m %= p;
    1622     5844009 :     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
    1623     5844009 :     uel(u,i) = m;
    1624             :   }
    1625      427860 :   return u;
    1626             : }
    1627             : static GEN
    1628     2991044 : Fl_get_col(GEN a, GEN b, long li, ulong p)
    1629             : {
    1630     2991044 :   GEN u = cgetg(li+1,t_VECSMALL);
    1631     2991044 :   ulong m = uel(b,li) % p;
    1632             :   long i,j;
    1633             : 
    1634     2991044 :   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
    1635    22499272 :   for (i=li-1; i>0; i--)
    1636             :   {
    1637    19508228 :     m = b[i]%p;
    1638   224163212 :     for (j = i+1; j <= li; j++)
    1639   204654984 :       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
    1640    19508228 :     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
    1641    19508228 :     uel(u,i) = m;
    1642             :   }
    1643     2991044 :   return u;
    1644             : }
    1645             : 
    1646             : /* bk -= m * bi */
    1647             : static void
    1648     4457481 : _submul(GEN b, long k, long i, GEN m)
    1649             : {
    1650     4457481 :   gel(b,k) = gsub(gel(b,k), gmul(m, gel(b,i)));
    1651     4457481 : }
    1652             : static int
    1653      699779 : init_gauss(GEN a, GEN *b, long *aco, long *li, int *iscol)
    1654             : {
    1655      699779 :   *iscol = *b ? (typ(*b) == t_COL): 0;
    1656      699779 :   *aco = lg(a) - 1;
    1657      699779 :   if (!*aco) /* a empty */
    1658             :   {
    1659         322 :     if (*b && lg(*b) != 1) pari_err_DIM("gauss");
    1660         322 :     *li = 0; return 0;
    1661             :   }
    1662      699457 :   *li = nbrows(a);
    1663      699457 :   if (*li < *aco) pari_err_INV("gauss [no left inverse]", a);
    1664      699457 :   if (*b)
    1665             :   {
    1666      686488 :     if (*li != *aco) pari_err_DIM("gauss");
    1667      686488 :     switch(typ(*b))
    1668             :     {
    1669             :       case t_MAT:
    1670       53843 :         if (lg(*b) == 1) return 0;
    1671       53843 :         *b = RgM_shallowcopy(*b);
    1672       53843 :         break;
    1673             :       case t_COL:
    1674      632645 :         *b = mkmat( leafcopy(*b) );
    1675      632645 :         break;
    1676           0 :       default: pari_err_TYPE("gauss",*b);
    1677             :     }
    1678      686488 :     if (nbrows(*b) != *li) pari_err_DIM("gauss");
    1679             :   }
    1680             :   else
    1681       12969 :     *b = matid(*li);
    1682      699457 :   return 1;
    1683             : }
    1684             : 
    1685             : static int
    1686      337694 : is_modular_solve(GEN a, GEN b, GEN *u)
    1687             : {
    1688      337694 :   GEN p = NULL;
    1689             :   ulong pp;
    1690      337694 :   if (!RgM_is_FpM(a, &p) || !p) return 0;
    1691         168 :   if (!b)
    1692             :   {
    1693          70 :     a = RgM_Fp_init(a, p, &pp);
    1694          70 :     switch(pp)
    1695             :     {
    1696             :     case 0:
    1697          14 :       a = FpM_inv(a,p);
    1698          14 :       if (a) a = FpM_to_mod(a, p);
    1699          14 :       break;
    1700             :     case 2:
    1701          35 :       a = F2m_inv(a);
    1702          35 :       if (a) a = F2m_to_mod(a);
    1703          35 :       break;
    1704             :     default:
    1705          21 :       a = Flm_inv(a,pp);
    1706          21 :       if (a) a = Flm_to_mod(a, pp);
    1707             :     }
    1708             :   }
    1709          98 :   else switch(typ(b))
    1710             :   {
    1711             :     case t_COL:
    1712          49 :       if (!RgV_is_FpV(b, &p)) return 0;
    1713          49 :       a = RgM_Fp_init(a, p, &pp);
    1714          49 :       switch(pp)
    1715             :       {
    1716             :       case 0:
    1717          14 :         b = RgC_to_FpC(b, p);
    1718          14 :         a = FpM_FpC_gauss(a,b,p);
    1719          14 :         if (a) a = FpC_to_mod(a, p);
    1720          14 :         break;
    1721             :       case 2:
    1722          14 :         b = RgV_to_F2v(b);
    1723          14 :         a = F2m_F2c_gauss(a,b);
    1724          14 :         if (a) a = F2c_to_mod(a);
    1725          14 :         break;
    1726             :       default:
    1727          21 :         b = RgV_to_Flv(b, pp);
    1728          21 :         a = Flm_Flc_gauss(a,b,pp);
    1729          21 :         if (a) a = Flc_to_mod(a, pp);
    1730          21 :         break;
    1731             :       }
    1732          49 :       break;
    1733             :     case t_MAT:
    1734          49 :       if (!RgM_is_FpM(b, &p)) return 0;
    1735          49 :       a = RgM_Fp_init(a, p, &pp);
    1736          49 :       switch(pp)
    1737             :       {
    1738             :       case 0:
    1739          14 :         b = RgM_to_FpM(b, p);
    1740          14 :         a = FpM_gauss(a,b,p);
    1741          14 :         if (a) a = FpM_to_mod(a, p);
    1742          14 :         break;
    1743             :       case 2:
    1744          14 :         b = RgM_to_F2m(b);
    1745          14 :         a = F2m_gauss(a,b);
    1746          14 :         if (a) a = F2m_to_mod(a);
    1747          14 :         break;
    1748             :       default:
    1749          21 :         b = RgM_to_Flm(b, pp);
    1750          21 :         a = Flm_gauss(a,b,pp);
    1751          21 :         if (a) a = Flm_to_mod(a, pp);
    1752          21 :         break;
    1753             :       }
    1754          49 :       break;
    1755           0 :     default: return 0;
    1756             :   }
    1757         168 :   *u = a; return 1;
    1758             : }
    1759             : /* Gaussan Elimination. If a is square, return a^(-1)*b;
    1760             :  * if a has more rows than columns and b is NULL, return c such that c a = Id.
    1761             :  * a is a (not necessarily square) matrix
    1762             :  * b is a matrix or column vector, NULL meaning: take the identity matrix,
    1763             :  *   effectively returning the inverse of a
    1764             :  * If a and b are empty, the result is the empty matrix.
    1765             :  *
    1766             :  * li: number of rows of a and b
    1767             :  * aco: number of columns of a
    1768             :  * bco: number of columns of b (if matrix)
    1769             :  */
    1770             : GEN
    1771      337694 : RgM_solve(GEN a, GEN b)
    1772             : {
    1773      337694 :   pari_sp av = avma;
    1774             :   long i, j, k, li, bco, aco;
    1775             :   int iscol;
    1776             :   pivot_fun pivot;
    1777      337694 :   GEN p, u, data, ff = NULL;
    1778             : 
    1779      337694 :   if (is_modular_solve(a,b,&u)) return gerepileupto(av, u);
    1780      337526 :   if (RgM_is_FFM(a, &ff)) {
    1781         189 :     if (!b)
    1782         105 :       return FFM_inv(a, ff);
    1783          84 :     if (typ(b) == t_COL && RgC_is_FFC(b, &ff))
    1784          56 :       return FFM_FFC_gauss(a, b, ff);
    1785          28 :     if (typ(b) == t_MAT && RgM_is_FFM(b, &ff))
    1786          28 :       return FFM_gauss(a, b, ff);
    1787             :   }
    1788      337337 :   avma = av;
    1789             : 
    1790      337337 :   if (lg(a)-1 == 2 && nbrows(a) == 2) {
    1791             :     /* 2x2 matrix, start by inverting a */
    1792       90060 :     GEN u = gcoeff(a,1,1), v = gcoeff(a,1,2);
    1793       90060 :     GEN w = gcoeff(a,2,1), x = gcoeff(a,2,2);
    1794       90060 :     GEN D = gsub(gmul(u,x), gmul(v,w)), ainv;
    1795       90060 :     if (gequal0(D)) return NULL;
    1796       90060 :     ainv = mkmat2(mkcol2(x, gneg(w)), mkcol2(gneg(v), u));
    1797       90060 :     ainv = gmul(ainv, ginv(D));
    1798       90060 :     if (b) ainv = gmul(ainv, b);
    1799       90060 :     return gerepileupto(av, ainv);
    1800             :   }
    1801             : 
    1802      247277 :   if (!init_gauss(a, &b, &aco, &li, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
    1803      247025 :   pivot = get_pivot_fun(a, a, &data);
    1804      247025 :   a = RgM_shallowcopy(a);
    1805      247025 :   bco = lg(b)-1;
    1806      247025 :   if(DEBUGLEVEL>4) err_printf("Entering gauss\n");
    1807             : 
    1808      247025 :   p = NULL; /* gcc -Wall */
    1809      880873 :   for (i=1; i<=aco; i++)
    1810             :   {
    1811             :     /* k is the line where we find the pivot */
    1812      880873 :     k = pivot(a, data, i, NULL);
    1813      880873 :     if (k > li) return NULL;
    1814      880866 :     if (k != i)
    1815             :     { /* exchange the lines s.t. k = i */
    1816      155612 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    1817      155612 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    1818             :     }
    1819      880866 :     p = gcoeff(a,i,i);
    1820      880866 :     if (i == aco) break;
    1821             : 
    1822     2411481 :     for (k=i+1; k<=li; k++)
    1823             :     {
    1824     1777633 :       GEN m = gcoeff(a,k,i);
    1825     1777633 :       if (!gequal0(m))
    1826             :       {
    1827      767657 :         m = gdiv(m,p);
    1828      767657 :         for (j=i+1; j<=aco; j++) _submul(gel(a,j),k,i,m);
    1829      767657 :         for (j=1;   j<=bco; j++) _submul(gel(b,j),k,i,m);
    1830             :       }
    1831             :     }
    1832      633848 :     if (gc_needed(av,1))
    1833             :     {
    1834           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gauss. i=%ld",i);
    1835           0 :       gerepileall(av,2, &a,&b);
    1836             :     }
    1837             :   }
    1838             : 
    1839      247018 :   if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
    1840      247018 :   u = cgetg(bco+1,t_MAT);
    1841      247018 :   for (j=1; j<=bco; j++) gel(u,j) = get_col(a,gel(b,j),p,aco);
    1842      247018 :   return gerepilecopy(av, iscol? gel(u,1): u);
    1843             : }
    1844             : 
    1845             : /* assume dim A >= 1, A invertible + upper triangular  */
    1846             : static GEN
    1847       16409 : RgM_inv_upper_ind(GEN A, long index)
    1848             : {
    1849       16409 :   long n = lg(A)-1, i = index, j;
    1850       16409 :   GEN u = zerocol(n);
    1851       16409 :   gel(u,i) = ginv(gcoeff(A,i,i));
    1852       45812 :   for (i--; i>0; i--)
    1853             :   {
    1854       29403 :     pari_sp av = avma;
    1855       29403 :     GEN m = gneg(gmul(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
    1856       29403 :     for (j=i+2; j<=n; j++) m = gsub(m, gmul(gcoeff(A,i,j),gel(u,j)));
    1857       29403 :     gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(A,i,i)));
    1858             :   }
    1859       16409 :   return u;
    1860             : }
    1861             : GEN
    1862        5922 : RgM_inv_upper(GEN A)
    1863             : {
    1864             :   long i, l;
    1865        5922 :   GEN B = cgetg_copy(A, &l);
    1866        5922 :   for (i = 1; i < l; i++) gel(B,i) = RgM_inv_upper_ind(A, i);
    1867        5922 :   return B;
    1868             : }
    1869             : 
    1870             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal  */
    1871             : static GEN
    1872          14 : FpM_inv_upper_1_ind(GEN A, long index, GEN p)
    1873             : {
    1874          14 :   long n = lg(A)-1, i = index, j;
    1875          14 :   GEN u = zerocol(n);
    1876          14 :   gel(u,i) = gen_1;
    1877          21 :   for (i--; i>0; i--)
    1878             :   {
    1879           7 :     pari_sp av = avma;
    1880           7 :     GEN m = negi(mulii(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
    1881           7 :     for (j=i+2; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    1882           7 :     gel(u,i) = gerepileuptoint(av, modii(m,p));
    1883             :   }
    1884          14 :   return u;
    1885             : }
    1886             : static GEN
    1887           7 : FpM_inv_upper_1(GEN A, GEN p)
    1888             : {
    1889             :   long i, l;
    1890           7 :   GEN B = cgetg_copy(A, &l);
    1891           7 :   for (i = 1; i < l; i++) gel(B,i) = FpM_inv_upper_1_ind(A, i, p);
    1892           7 :   return B;
    1893             : }
    1894             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
    1895             :  * reduced mod p */
    1896             : static GEN
    1897          28 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
    1898             : {
    1899          28 :   long n = lg(A)-1, i = index, j;
    1900          28 :   GEN u = const_vecsmall(n, 0);
    1901          28 :   u[i] = 1;
    1902          28 :   if (SMALL_ULONG(p))
    1903          21 :     for (i--; i>0; i--)
    1904             :     {
    1905           7 :       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
    1906           7 :       for (j=i+2; j<=n; j++)
    1907             :       {
    1908           0 :         if (m & HIGHMASK) m %= p;
    1909           0 :         m += ucoeff(A,i,j) * uel(u,j);
    1910             :       }
    1911           7 :       u[i] = Fl_neg(m % p, p);
    1912             :     }
    1913             :   else
    1914          21 :     for (i--; i>0; i--)
    1915             :     {
    1916           7 :       ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
    1917           7 :       for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
    1918           7 :       u[i] = Fl_neg(m, p);
    1919             :     }
    1920          28 :   return u;
    1921             : }
    1922             : static GEN
    1923          14 : F2m_inv_upper_1_ind(GEN A, long index)
    1924             : {
    1925          14 :   pari_sp av = avma;
    1926          14 :   long n = lg(A)-1, i = index, j;
    1927          14 :   GEN u = const_vecsmall(n, 0);
    1928          14 :   u[i] = 1;
    1929          21 :   for (i--; i>0; i--)
    1930             :   {
    1931           7 :     ulong m = F2m_coeff(A,i,i+1) & uel(u,i+1); /* j = i+1 */
    1932           7 :     for (j=i+2; j<=n; j++) m ^= F2m_coeff(A,i,j) & uel(u,j);
    1933           7 :     u[i] = m & 1;
    1934             :   }
    1935          14 :   return gerepileuptoleaf(av, Flv_to_F2v(u));
    1936             : }
    1937             : static GEN
    1938          14 : Flm_inv_upper_1(GEN A, ulong p)
    1939             : {
    1940             :   long i, l;
    1941          14 :   GEN B = cgetg_copy(A, &l);
    1942          14 :   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
    1943          14 :   return B;
    1944             : }
    1945             : static GEN
    1946           7 : F2m_inv_upper_1(GEN A)
    1947             : {
    1948             :   long i, l;
    1949           7 :   GEN B = cgetg_copy(A, &l);
    1950           7 :   for (i = 1; i < l; i++) gel(B,i) = F2m_inv_upper_1_ind(A, i);
    1951           7 :   return B;
    1952             : }
    1953             : 
    1954             : static GEN
    1955      982978 : split_realimag_col(GEN z, long r1, long r2)
    1956             : {
    1957      982978 :   long i, ru = r1+r2;
    1958      982978 :   GEN x = cgetg(ru+r2+1,t_COL), y = x + r2;
    1959     2903068 :   for (i=1; i<=r1; i++) {
    1960     1920090 :     GEN a = gel(z,i);
    1961     1920090 :     if (typ(a) == t_COMPLEX) a = gel(a,1); /* paranoia: a should be real */
    1962     1920090 :     gel(x,i) = a;
    1963             :   }
    1964     1740861 :   for (   ; i<=ru; i++) {
    1965      757883 :     GEN b, a = gel(z,i);
    1966      757883 :     if (typ(a) == t_COMPLEX) { b = gel(a,2); a = gel(a,1); } else b = gen_0;
    1967      757883 :     gel(x,i) = a;
    1968      757883 :     gel(y,i) = b;
    1969             :   }
    1970      982978 :   return x;
    1971             : }
    1972             : GEN
    1973      486769 : split_realimag(GEN x, long r1, long r2)
    1974             : {
    1975             :   long i,l; GEN y;
    1976      486769 :   if (typ(x) == t_COL) return split_realimag_col(x,r1,r2);
    1977      239668 :   y = cgetg_copy(x, &l);
    1978      239668 :   for (i=1; i<l; i++) gel(y,i) = split_realimag_col(gel(x,i), r1, r2);
    1979      239668 :   return y;
    1980             : }
    1981             : 
    1982             : /* assume M = (r1+r2) x (r1+2r2) matrix and y compatible vector or matrix
    1983             :  * r1 first lines of M,y are real. Solve the system obtained by splitting
    1984             :  * real and imaginary parts. */
    1985             : GEN
    1986      237704 : RgM_solve_realimag(GEN M, GEN y)
    1987             : {
    1988      237704 :   long l = lg(M), r2 = l - lgcols(M), r1 = l-1 - 2*r2;
    1989      237704 :   return RgM_solve(split_realimag(M, r1,r2),
    1990             :                    split_realimag(y, r1,r2));
    1991             : }
    1992             : 
    1993             : GEN
    1994         294 : gauss(GEN a, GEN b)
    1995             : {
    1996             :   GEN z;
    1997         294 :   if (typ(a)!=t_MAT) pari_err_TYPE("gauss",a);
    1998         378 :   if (RgM_is_ZM(a) && b &&
    1999         133 :       ((typ(b) == t_COL && RgV_is_ZV(b)) || (typ(b) == t_MAT && RgM_is_ZM(b))))
    2000          84 :     z = ZM_gauss(a,b);
    2001             :   else
    2002         210 :     z = RgM_solve(a,b);
    2003         294 :   if (!z) pari_err_INV("gauss",a);
    2004         231 :   return z;
    2005             : }
    2006             : 
    2007             : static GEN
    2008       52448 : F2_get_col(GEN b, GEN d, long li, long aco)
    2009             : {
    2010       52448 :   long i, l = nbits2lg(aco);
    2011       52448 :   GEN u = cgetg(l, t_VECSMALL);
    2012       52448 :   u[1] = aco;
    2013      992944 :   for (i = 1; i <= li; i++)
    2014      940496 :     if (d[i]) /* d[i] can still be 0 if li > aco */
    2015             :     {
    2016      940461 :       if (F2v_coeff(b, i))
    2017      328748 :         F2v_set(u, d[i]);
    2018             :       else
    2019      611713 :         F2v_clear(u, d[i]);
    2020             :     }
    2021       52448 :   return u;
    2022             : }
    2023             : 
    2024             : /* destroy a, b */
    2025             : static GEN
    2026        9429 : F2m_gauss_sp(GEN a, GEN b)
    2027             : {
    2028        9429 :   long i, j, k, l, li, bco, aco = lg(a)-1;
    2029             :   GEN u, d;
    2030             : 
    2031        9429 :   if (!aco) return cgetg(1,t_MAT);
    2032        9429 :   li = gel(a,1)[1];
    2033        9429 :   d = zero_Flv(li);
    2034        9429 :   bco = lg(b)-1;
    2035       61891 :   for (i=1; i<=aco; i++)
    2036             :   {
    2037       52476 :     GEN ai = vecsmall_copy(gel(a,i));
    2038       52476 :     if (!d[i] && F2v_coeff(ai, i))
    2039       28363 :       k = i;
    2040             :     else
    2041       24113 :       for (k = 1; k <= li; k++) if (!d[k] && F2v_coeff(ai,k)) break;
    2042             :     /* found a pivot on row k */
    2043       52476 :     if (k > li) return NULL;
    2044       52462 :     d[k] = i;
    2045             : 
    2046             :     /* Clear k-th row but column-wise instead of line-wise */
    2047             :     /*  a_ij -= a_i1*a1j/a_11
    2048             :        line-wise grouping:  L_j -= a_1j/a_11*L_1
    2049             :        column-wise:         C_i -= a_i1/a_11*C_1
    2050             :     */
    2051       52462 :     F2v_clear(ai,k);
    2052      992965 :     for (l=1; l<=aco; l++)
    2053             :     {
    2054      940503 :       GEN al = gel(a,l);
    2055      940503 :       if (F2v_coeff(al,k)) F2v_add_inplace(al,ai);
    2056             :     }
    2057      992930 :     for (l=1; l<=bco; l++)
    2058             :     {
    2059      940468 :       GEN bl = gel(b,l);
    2060      940468 :       if (F2v_coeff(bl,k)) F2v_add_inplace(bl,ai);
    2061             :     }
    2062             :   }
    2063        9415 :   u = cgetg(bco+1,t_MAT);
    2064        9415 :   for (j = 1; j <= bco; j++) gel(u,j) = F2_get_col(gel(b,j), d, li, aco);
    2065        9415 :   return u;
    2066             : }
    2067             : 
    2068             : GEN
    2069          28 : F2m_gauss(GEN a, GEN b)
    2070             : {
    2071          28 :   pari_sp av = avma;
    2072          28 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    2073          28 :   return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), F2m_copy(b)));
    2074             : }
    2075             : GEN
    2076          14 : F2m_F2c_gauss(GEN a, GEN b)
    2077             : {
    2078          14 :   pari_sp av = avma;
    2079          14 :   GEN z = F2m_gauss(a, mkmat(b));
    2080          14 :   if (!z) { avma = av; return NULL; }
    2081           7 :   if (lg(z) == 1) { avma = av; return cgetg(1,t_VECSMALL); }
    2082           7 :   return gerepileuptoleaf(av, gel(z,1));
    2083             : }
    2084             : 
    2085             : GEN
    2086          35 : F2m_inv(GEN a)
    2087             : {
    2088          35 :   pari_sp av = avma;
    2089          35 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    2090          35 :   return gerepileupto(av, F2m_gauss_sp(F2m_copy(a), matid_F2m(gel(a,1)[1])));
    2091             : }
    2092             : 
    2093             : /* destroy a, b */
    2094             : static GEN
    2095       44400 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
    2096             : {
    2097       44400 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
    2098       44400 :   ulong det = 1;
    2099             :   GEN u;
    2100             : 
    2101       44400 :   li = nbrows(a);
    2102       44400 :   bco = lg(b)-1;
    2103      427881 :   for (i=1; i<=aco; i++)
    2104             :   {
    2105             :     ulong invpiv;
    2106             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
    2107      427881 :     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
    2108      823102 :     for (k = i; k <= li; k++)
    2109             :     {
    2110      823095 :       ulong piv = ( ucoeff(a,k,i) %= p );
    2111      823095 :       if (piv)
    2112             :       {
    2113      427874 :         ucoeff(a,k,i) = Fl_inv(piv, p);
    2114      427874 :         if (detp)
    2115             :         {
    2116           0 :           if (det & HIGHMASK) det %= p;
    2117           0 :           det *= piv;
    2118             :         }
    2119      427874 :         break;
    2120             :       }
    2121             :     }
    2122             :     /* found a pivot on line k */
    2123      427881 :     if (k > li) return NULL;
    2124      427874 :     if (k != i)
    2125             :     { /* swap lines so that k = i */
    2126      136963 :       s = -s;
    2127      136963 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    2128      136963 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    2129             :     }
    2130      427874 :     if (i == aco) break;
    2131             : 
    2132      383481 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    2133     3305517 :     for (k=i+1; k<=li; k++)
    2134             :     {
    2135     2922036 :       ulong m = ( ucoeff(a,k,i) %= p );
    2136     2922036 :       if (!m) continue;
    2137             : 
    2138      727105 :       m = Fl_mul(m, invpiv, p);
    2139      727105 :       if (m == 1) {
    2140      182410 :         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
    2141      182410 :         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
    2142             :       } else {
    2143      544695 :         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
    2144      544695 :         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
    2145             :       }
    2146             :     }
    2147             :   }
    2148       44393 :   if (detp)
    2149             :   {
    2150           0 :     det %=  p;
    2151           0 :     if (s < 0 && det) det = p - det;
    2152           0 :     *detp = det;
    2153             :   }
    2154       44393 :   u = cgetg(bco+1,t_MAT);
    2155       44393 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
    2156       44393 :   return u;
    2157             : }
    2158             : 
    2159             : /* destroy a, b */
    2160             : static GEN
    2161      713531 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p)
    2162             : {
    2163      713531 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
    2164      713531 :   ulong det = 1;
    2165             :   GEN u;
    2166             :   ulong pi;
    2167      713531 :   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
    2168      713531 :   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
    2169      669131 :   pi = get_Fl_red(p);
    2170      669131 :   li = nbrows(a);
    2171      669131 :   bco = lg(b)-1;
    2172     2991058 :   for (i=1; i<=aco; i++)
    2173             :   {
    2174             :     ulong invpiv;
    2175             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
    2176     4230716 :     for (k = i; k <= li; k++)
    2177             :     {
    2178     4230716 :       ulong piv = ucoeff(a,k,i);
    2179     4230716 :       if (piv)
    2180             :       {
    2181     2991058 :         ucoeff(a,k,i) = Fl_inv(piv, p);
    2182     2991058 :         if (detp) det = Fl_mul_pre(det, piv, p, pi);
    2183     2991058 :         break;
    2184             :       }
    2185             :     }
    2186             :     /* found a pivot on line k */
    2187     2991058 :     if (k > li) return NULL;
    2188     2991058 :     if (k != i)
    2189             :     { /* swap lines so that k = i */
    2190      437332 :       s = -s;
    2191      437332 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    2192      437332 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    2193             :     }
    2194     2991058 :     if (i == aco) break;
    2195             : 
    2196     2321927 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    2197    12076055 :     for (k=i+1; k<=li; k++)
    2198             :     {
    2199     9754128 :       ulong m = ucoeff(a,k,i);
    2200     9754128 :       if (!m) continue;
    2201             : 
    2202     3045190 :       m = Fl_mul_pre(m, invpiv, p, pi);
    2203     3045190 :       if (m == 1) {
    2204      111432 :         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
    2205      111432 :         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
    2206             :       } else {
    2207     2933758 :         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
    2208     2933758 :         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
    2209             :       }
    2210             :     }
    2211             :   }
    2212      669131 :   if (detp)
    2213             :   {
    2214           0 :     if (s < 0 && det) det = p - det;
    2215           0 :     *detp = det;
    2216             :   }
    2217      669131 :   u = cgetg(bco+1,t_MAT);
    2218      669131 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
    2219      669131 :   return u;
    2220             : }
    2221             : 
    2222             : GEN
    2223          42 : Flm_gauss(GEN a, GEN b, ulong p) {
    2224          42 :   return Flm_gauss_sp(RgM_shallowcopy(a), RgM_shallowcopy(b), NULL, p);
    2225             : }
    2226             : static GEN
    2227      683323 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
    2228      683323 :   if (lg(a) == 1) return cgetg(1,t_MAT);
    2229      683323 :   return Flm_gauss_sp(a, matid_Flm(nbrows(a)), detp, p);
    2230             : }
    2231             : GEN
    2232      466656 : Flm_inv(GEN a, ulong p) {
    2233      466656 :   return Flm_inv_sp(RgM_shallowcopy(a), NULL, p);
    2234             : }
    2235             : GEN
    2236          21 : Flm_Flc_gauss(GEN a, GEN b, ulong p) {
    2237          21 :   pari_sp av = avma;
    2238          21 :   GEN z = Flm_gauss(a, mkmat(b), p);
    2239          21 :   if (!z) { avma = av; return NULL; }
    2240          14 :   if (lg(z) == 1) { avma = av; return cgetg(1,t_VECSMALL); }
    2241          14 :   return gerepileuptoleaf(av, gel(z,1));
    2242             : }
    2243             : 
    2244             : static GEN
    2245        2807 : FpM_gauss_gen(GEN a, GEN b, GEN p)
    2246             : {
    2247             :   void *E;
    2248        2807 :   const struct bb_field *S = get_Fp_field(&E,p);
    2249        2807 :   return gen_Gauss(a,b, E, S);
    2250             : }
    2251             : /* a an FpM, lg(a)>1; b an FpM or NULL (replace by identity) */
    2252             : static GEN
    2253       42339 : FpM_gauss_i(GEN a, GEN b, GEN p, ulong *pp)
    2254             : {
    2255       42339 :   long n = nbrows(a);
    2256       42339 :   a = FpM_init(a,p,pp);
    2257       42339 :   switch(*pp)
    2258             :   {
    2259             :   case 0:
    2260        2807 :     if (!b) b = matid(n);
    2261        2807 :     return FpM_gauss_gen(a,b,p);
    2262             :   case 2:
    2263        9366 :     if (b) b = ZM_to_F2m(b); else b = matid_F2m(n);
    2264        9366 :     return F2m_gauss_sp(a,b);
    2265             :   default:
    2266       30166 :     if (b) b = ZM_to_Flm(b, *pp); else b = matid_Flm(n);
    2267       30166 :     return Flm_gauss_sp(a,b, NULL, *pp);
    2268             :   }
    2269             : }
    2270             : GEN
    2271          14 : FpM_gauss(GEN a, GEN b, GEN p)
    2272             : {
    2273          14 :   pari_sp av = avma;
    2274             :   ulong pp;
    2275             :   GEN u;
    2276          14 :   if (lg(a) == 1 || lg(b)==1) return cgetg(1, t_MAT);
    2277          14 :   u = FpM_gauss_i(a, b, p, &pp);
    2278          14 :   if (!u) { avma = av; return NULL; }
    2279          14 :   switch(pp)
    2280             :   {
    2281          14 :   case 0: return gerepilecopy(av, u);
    2282           0 :   case 2:  u = F2m_to_ZM(u); break;
    2283           0 :   default: u = Flm_to_ZM(u); break;
    2284             :   }
    2285           0 :   return gerepileupto(av, u);
    2286             : }
    2287             : GEN
    2288       42311 : FpM_inv(GEN a, GEN p)
    2289             : {
    2290       42311 :   pari_sp av = avma;
    2291             :   ulong pp;
    2292             :   GEN u;
    2293       42311 :   if (lg(a) == 1) return cgetg(1, t_MAT);
    2294       42311 :   u = FpM_gauss_i(a, NULL, p, &pp);
    2295       42311 :   if (!u) { avma = av; return NULL; }
    2296       42311 :   switch(pp)
    2297             :   {
    2298        2779 :   case 0: return gerepilecopy(av, u);
    2299        9366 :   case 2:  u = F2m_to_ZM(u); break;
    2300       30166 :   default: u = Flm_to_ZM(u); break;
    2301             :   }
    2302       39532 :   return gerepileupto(av, u);
    2303             : }
    2304             : 
    2305             : GEN
    2306          14 : FpM_FpC_gauss(GEN a, GEN b, GEN p)
    2307             : {
    2308          14 :   pari_sp av = avma;
    2309             :   ulong pp;
    2310             :   GEN u;
    2311          14 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2312          14 :   u = FpM_gauss_i(a, mkmat(b), p, &pp);
    2313          14 :   if (!u) { avma = av; return NULL; }
    2314          14 :   switch(pp)
    2315             :   {
    2316          14 :   case 0: return gerepilecopy(av, gel(u,1));
    2317           0 :   case 2:  u = F2c_to_ZC(gel(u,1)); break;
    2318           0 :   default: u = Flc_to_ZC(gel(u,1)); break;
    2319             :   }
    2320           0 :   return gerepileupto(av, u);
    2321             : }
    2322             : 
    2323             : static GEN
    2324          56 : FlxqM_gauss_gen(GEN a, GEN b, GEN T, ulong p)
    2325             : {
    2326             :   void *E;
    2327          56 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    2328          56 :   return gen_Gauss(a, b, E, S);
    2329             : }
    2330             : GEN
    2331           7 : FlxqM_gauss(GEN a, GEN b, GEN T, ulong p)
    2332             : {
    2333           7 :   pari_sp av = avma;
    2334           7 :   long n = lg(a)-1;
    2335             :   GEN u;
    2336           7 :   if (!n || lg(b)==1) { avma = av; return cgetg(1, t_MAT); }
    2337           7 :   u = FlxqM_gauss_gen(a, b, T, p);
    2338           7 :   if (!u) { avma = av; return NULL; }
    2339           7 :   return gerepilecopy(av, u);
    2340             : }
    2341             : GEN
    2342          35 : FlxqM_inv(GEN a, GEN T, ulong p)
    2343             : {
    2344          35 :   pari_sp av = avma;
    2345             :   GEN u;
    2346          35 :   if (lg(a) == 1) { avma = av; return cgetg(1, t_MAT); }
    2347          35 :   u = FlxqM_gauss_gen(a, matid_FlxqM(nbrows(a),T,p), T,p);
    2348          35 :   if (!u) { avma = av; return NULL; }
    2349          28 :   return gerepilecopy(av, u);
    2350             : }
    2351             : GEN
    2352          14 : FlxqM_FlxqC_gauss(GEN a, GEN b, GEN T, ulong p)
    2353             : {
    2354          14 :   pari_sp av = avma;
    2355             :   GEN u;
    2356          14 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2357          14 :   u = FlxqM_gauss_gen(a, mkmat(b), T, p);
    2358          14 :   if (!u) { avma = av; return NULL; }
    2359           7 :   return gerepilecopy(av, gel(u,1));
    2360             : }
    2361             : 
    2362             : static GEN
    2363          56 : FqM_gauss_gen(GEN a, GEN b, GEN T, GEN p)
    2364             : {
    2365             :   void *E;
    2366          56 :   const struct bb_field *S = get_Fq_field(&E,T,p);
    2367          56 :   return gen_Gauss(a,b,E,S);
    2368             : }
    2369             : GEN
    2370           7 : FqM_gauss(GEN a, GEN b, GEN T, GEN p)
    2371             : {
    2372           7 :   pari_sp av = avma;
    2373             :   GEN u;
    2374             :   long n;
    2375           7 :   if (!T) return FpM_gauss(a,b,p);
    2376           7 :   n = lg(a)-1; if (!n || lg(b)==1) return cgetg(1, t_MAT);
    2377           7 :   u = FqM_gauss_gen(a,b,T,p);
    2378           7 :   if (!u) { avma = av; return NULL; }
    2379           7 :   return gerepilecopy(av, u);
    2380             : }
    2381             : GEN
    2382          35 : FqM_inv(GEN a, GEN T, GEN p)
    2383             : {
    2384          35 :   pari_sp av = avma;
    2385             :   GEN u;
    2386          35 :   if (!T) return FpM_inv(a,p);
    2387          35 :   if (lg(a) == 1) return cgetg(1, t_MAT);
    2388          35 :   u = FqM_gauss_gen(a,matid(nbrows(a)),T,p);
    2389          35 :   if (!u) { avma = av; return NULL; }
    2390          28 :   return gerepilecopy(av, u);
    2391             : }
    2392             : GEN
    2393          14 : FqM_FqC_gauss(GEN a, GEN b, GEN T, GEN p)
    2394             : {
    2395          14 :   pari_sp av = avma;
    2396             :   GEN u;
    2397          14 :   if (!T) return FpM_FpC_gauss(a,b,p);
    2398          14 :   if (lg(a) == 1) return cgetg(1, t_COL);
    2399          14 :   u = FqM_gauss_gen(a,mkmat(b),T,p);
    2400          14 :   if (!u) { avma = av; return NULL; }
    2401           7 :   return gerepilecopy(av, gel(u,1));
    2402             : }
    2403             : 
    2404             : /* Dixon p-adic lifting algorithm.
    2405             :  * Numer. Math. 40, 137-141 (1982), DOI: 10.1007/BF01459082 */
    2406             : GEN
    2407      452502 : ZM_gauss(GEN a, GEN b0)
    2408             : {
    2409      452502 :   pari_sp av = avma, av2;
    2410             :   int iscol;
    2411             :   long n, ncol, i, m, elim;
    2412             :   ulong p;
    2413      452502 :   GEN N, C, delta, xb, nb, nmin, res, b = b0;
    2414             :   forprime_t S;
    2415             : 
    2416      452502 :   if (!init_gauss(a, &b, &n, &ncol, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
    2417      452432 :   nb = gen_0; ncol = lg(b);
    2418      914908 :   for (i = 1; i < ncol; i++)
    2419             :   {
    2420      462476 :     GEN ni = gnorml2(gel(b, i));
    2421      462476 :     if (cmpii(nb, ni) < 0) nb = ni;
    2422             :   }
    2423      452432 :   if (!signe(nb)) { avma = av; return gcopy(b0); }
    2424      452432 :   delta = gen_1; nmin = nb;
    2425     1846285 :   for (i = 1; i <= n; i++)
    2426             :   {
    2427     1393853 :     GEN ni = gnorml2(gel(a, i));
    2428     1393853 :     if (cmpii(ni, nmin) < 0)
    2429             :     {
    2430        1894 :       delta = mulii(delta, nmin); nmin = ni;
    2431             :     }
    2432             :     else
    2433     1391959 :       delta = mulii(delta, ni);
    2434             :   }
    2435      452432 :   if (!signe(nmin)) return NULL;
    2436      452418 :   elim = expi(delta)+1;
    2437      452418 :   av2 = avma;
    2438      452418 :   init_modular_big(&S);
    2439             :   for(;;)
    2440             :   {
    2441      452418 :     p = u_forprime_next(&S);
    2442      452418 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    2443      452418 :     if (C) break;
    2444           0 :     elim -= expu(p);
    2445           0 :     if (elim < 0) return NULL;
    2446           0 :     avma = av2;
    2447           0 :   }
    2448             :   /* N.B. Our delta/lambda are SQUARES of those in the paper
    2449             :    * log(delta lambda) / log p, where lambda is 3+sqrt(5) / 2,
    2450             :    * whose log is < 1, hence + 1 (to cater for rounding errors) */
    2451      904836 :   m = (long)ceil((rtodbl(logr_abs(itor(delta,LOWDEFAULTPREC))) + 1)
    2452      452418 :                  / log((double)p));
    2453      452418 :   xb = ZlM_gauss(a, b, p, m, C);
    2454      452418 :   N = powuu(p, m);
    2455      452418 :   delta = sqrti(delta);
    2456      452418 :   if (iscol)
    2457      449691 :     res = FpC_ratlift(gel(xb,1), N, delta,delta, NULL);
    2458             :   else
    2459        2727 :     res = FpM_ratlift(xb, N, delta,delta, NULL);
    2460      452418 :   return gerepileupto(av, res);
    2461             : }
    2462             : 
    2463             : /* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */
    2464             : GEN
    2465      102576 : ZM_inv(GEN M, GEN dM)
    2466             : {
    2467      102576 :   pari_sp av2, av = avma;
    2468             :   GEN Hp,q,H;
    2469             :   ulong p;
    2470      102576 :   long lM = lg(M), stable = 0;
    2471      102576 :   int negate = 0;
    2472             :   forprime_t S;
    2473             :   pari_timer ti;
    2474             : 
    2475      102576 :   if (lM == 1) return cgetg(1,t_MAT);
    2476             : 
    2477             :   /* HACK: include dM = -1 ! */
    2478      101603 :   if (dM && is_pm1(dM))
    2479             :   {
    2480             :     /* modular algorithm computes M^{-1}, NOT multiplied by det(M) = -1.
    2481             :      * We will correct (negate) at the end. */
    2482       90849 :     if (signe(dM) < 0) negate = 1;
    2483       90849 :     dM = gen_1;
    2484             :   }
    2485      101603 :   init_modular_big(&S);
    2486      101603 :   av2 = avma;
    2487      101603 :   H = NULL;
    2488      101603 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2489      311486 :   while ((p = u_forprime_next(&S)))
    2490             :   {
    2491             :     ulong dMp;
    2492             :     GEN Mp;
    2493      209883 :     Mp = ZM_to_Flm(M,p);
    2494      209883 :     if (dM == gen_1)
    2495      184281 :       Hp = Flm_inv_sp(Mp, NULL, p);
    2496             :     else
    2497             :     {
    2498       25602 :       if (dM) {
    2499       25602 :         dMp = umodiu(dM,p); if (!dMp) continue;
    2500       25602 :         Hp = Flm_inv_sp(Mp, NULL, p);
    2501       25602 :         if (!Hp) pari_err_INV("ZM_inv", Mp);
    2502             :       } else {
    2503           0 :         Hp = Flm_inv_sp(Mp, &dMp, p);
    2504           0 :         if (!Hp) continue;
    2505             :       }
    2506       25602 :       if (dMp != 1) Flm_Fl_mul_inplace(Hp, dMp, p);
    2507             :     }
    2508             : 
    2509      209883 :     if (!H)
    2510             :     {
    2511      101603 :       H = ZM_init_CRT(Hp, p);
    2512      101603 :       q = utoipos(p);
    2513             :     }
    2514             :     else
    2515      108280 :       stable = ZM_incremental_CRT(&H, Hp, &q, p);
    2516      209883 :     if (DEBUGLEVEL>5) timer_printf(&ti, "ZM_inv mod %lu (stable=%ld)", p,stable);
    2517      209883 :     if (stable) {/* DONE ? */
    2518      101603 :       if (dM != gen_1)
    2519      112357 :       { if (ZM_isscalar(ZM_mul(M, H), dM)) break; }
    2520             :       else
    2521       90849 :       { if (ZM_isidentity(ZM_mul(M, H))) break; }
    2522             :     }
    2523             : 
    2524      108280 :     if (gc_needed(av,2))
    2525             :     {
    2526           9 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv");
    2527           9 :       gerepileall(av2, 2, &H, &q);
    2528             :     }
    2529             :   }
    2530      101603 :   if (!p) pari_err_OVERFLOW("ZM_inv [ran out of primes]");
    2531      101603 :   if (DEBUGLEVEL>5) err_printf("ZM_inv done\n");
    2532      101603 :   if (negate)
    2533           0 :     return gerepileupto(av, ZM_neg(H));
    2534             :   else
    2535      101603 :     return gerepilecopy(av, H);
    2536             : }
    2537             : 
    2538             : /* to be used when denom(M^(-1)) << det(M) and a sharp multiple is
    2539             :  * not available. Return H primitive such that M*H = den*Id */
    2540             : GEN
    2541        6713 : ZM_inv_ratlift(GEN M, GEN *pden)
    2542             : {
    2543        6713 :   pari_sp av2, av = avma;
    2544             :   GEN Hp, q, H;
    2545             :   ulong p;
    2546        6713 :   long lM = lg(M);
    2547             :   forprime_t S;
    2548        6713 :   if (lM == 1) { *pden = gen_1; return cgetg(1,t_MAT); }
    2549             : 
    2550        6020 :   init_modular_big(&S);
    2551        6020 :   av2 = avma;
    2552        6020 :   H = NULL;
    2553       12804 :   while ((p = u_forprime_next(&S)))
    2554             :   {
    2555             :     GEN Mp, B, Hr;
    2556        6784 :     Mp = ZM_to_Flm(M,p);
    2557        6784 :     Hp = Flm_inv_sp(Mp, NULL, p);
    2558        6784 :     if (!Hp) continue;
    2559        6784 :     if (!H)
    2560             :     {
    2561        6020 :       H = ZM_init_CRT(Hp, p);
    2562        6020 :       q = utoipos(p);
    2563             :     }
    2564             :     else
    2565         764 :       ZM_incremental_CRT(&H, Hp, &q, p);
    2566        6784 :     B = sqrti(shifti(q,-1));
    2567        6784 :     Hr = FpM_ratlift(H,q,B,B,NULL);
    2568        6784 :     if (DEBUGLEVEL>5) err_printf("ZM_inv mod %lu (ratlift=%ld)\n", p,!!Hr);
    2569        6784 :     if (Hr) {/* DONE ? */
    2570        6083 :       GEN Hl = Q_remove_denom(Hr, pden);
    2571        6083 :       if (*pden)
    2572        4676 :       { if (ZM_isscalar(ZM_mul(M, Hl), *pden)) { H = Hl; break; }}
    2573             :       else
    2574        1407 :       { if (ZM_isidentity(ZM_mul(M, Hl))) { H = Hl; *pden = gen_1; break; } }
    2575             :     }
    2576             : 
    2577         764 :     if (gc_needed(av,2))
    2578             :     {
    2579           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv_ratlift");
    2580           0 :       gerepileall(av2, 2, &H, &q);
    2581             :     }
    2582             :   }
    2583        6020 :   gerepileall(av, 2, &H, pden);
    2584        6020 :   return H;
    2585             : }
    2586             : 
    2587             : /* same as above, M rational */
    2588             : GEN
    2589        6447 : QM_inv(GEN M, GEN dM)
    2590             : {
    2591        6447 :   pari_sp av = avma;
    2592        6447 :   GEN cM, pM = Q_primitive_part(M, &cM);
    2593        6447 :   if (!cM) return ZM_inv(pM,dM);
    2594        2646 :   return gerepileupto(av, ZM_inv(pM, gdiv(dM,cM)));
    2595             : }
    2596             : 
    2597             : GEN
    2598        1877 : ZM_ker(GEN M)
    2599             : {
    2600        1877 :   pari_sp av2, av = avma;
    2601             :   GEN q, H, D;
    2602             :   forprime_t S;
    2603        1877 :   av2 = avma;
    2604        1877 :   H = NULL; D = NULL;
    2605        1877 :   if (lg(M)==1) return cgetg(1, t_MAT);
    2606        1870 :   init_modular_big(&S);
    2607             :   for(;;)
    2608             :   {
    2609             :     GEN Kp, Hp, Dp, Mp, Hr, B;
    2610        2107 :     ulong p = u_forprime_next(&S);
    2611        2107 :     Mp = ZM_to_Flm(M, p);
    2612        2107 :     Kp = Flm_ker_sp(Mp, p, 2);
    2613        2107 :     Hp = gel(Kp,1); Dp = gel(Kp,2);
    2614        2107 :     if (H && (lg(Hp)>lg(H) || (lg(Hp)==lg(H) && vecsmall_lexcmp(Dp,D)>0))) continue;
    2615        2107 :     if (!H || (lg(Hp)<lg(H) || vecsmall_lexcmp(Dp,D)<0))
    2616             :     {
    2617        1870 :       H = ZM_init_CRT(Hp, p); D = Dp;
    2618        1870 :       q = utoipos(p);
    2619             :     }
    2620             :     else
    2621         237 :       ZM_incremental_CRT(&H, Hp, &q, p);
    2622        2107 :     B = sqrti(shifti(q,-1));
    2623        2107 :     Hr = FpM_ratlift(H, q, B, B, NULL);
    2624        2107 :     if (DEBUGLEVEL>5) err_printf("ZM_ker mod %lu (ratlift=%ld)\n", p,!!Hr);
    2625        2107 :     if (Hr) {/* DONE ? */
    2626        1889 :       GEN Hl = vec_Q_primpart(Q_remove_denom(Hr, NULL));
    2627        1889 :       GEN MH = ZM_mul(M, Hl);
    2628        1889 :       if (gequal0(MH)) { H = Hl;  break; }
    2629             :     }
    2630         237 :     if (gc_needed(av,2))
    2631             :     {
    2632           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_ker");
    2633           0 :       gerepileall(av2, 3, &H, &D, &q);
    2634             :     }
    2635         237 :   }
    2636        1870 :   return gerepilecopy(av, H);
    2637             : }
    2638             : 
    2639             : /* x a ZM. Return a multiple of the determinant of the lattice generated by
    2640             :  * the columns of x. From Algorithm 2.2.6 in GTM138 */
    2641             : GEN
    2642          42 : detint(GEN A)
    2643             : {
    2644          42 :   if (typ(A) != t_MAT) pari_err_TYPE("detint",A);
    2645          42 :   RgM_check_ZM(A, "detint");
    2646          42 :   return ZM_detmult(A);
    2647             : }
    2648             : GEN
    2649        7623 : ZM_detmult(GEN A)
    2650             : {
    2651        7623 :   pari_sp av1, av = avma;
    2652             :   GEN B, c, v, piv;
    2653        7623 :   long rg, i, j, k, m, n = lg(A) - 1;
    2654             : 
    2655        7623 :   if (!n) return gen_1;
    2656        7623 :   m = nbrows(A);
    2657        7623 :   if (n < m) return gen_0;
    2658        7602 :   c = zero_zv(m);
    2659        7602 :   av1 = avma;
    2660        7602 :   B = zeromatcopy(m,m);
    2661        7602 :   v = cgetg(m+1, t_COL);
    2662        7602 :   piv = gen_1; rg = 0;
    2663       52927 :   for (k=1; k<=n; k++)
    2664             :   {
    2665       52920 :     GEN pivprec = piv;
    2666       52920 :     long t = 0;
    2667      461895 :     for (i=1; i<=m; i++)
    2668             :     {
    2669      408975 :       pari_sp av2 = avma;
    2670             :       GEN vi;
    2671      408975 :       if (c[i]) continue;
    2672             : 
    2673      231000 :       vi = mulii(piv, gcoeff(A,i,k));
    2674     2141580 :       for (j=1; j<=m; j++)
    2675     1910580 :         if (c[j]) vi = addii(vi, mulii(gcoeff(B,j,i),gcoeff(A,j,k)));
    2676      231000 :       if (!t && signe(vi)) t = i;
    2677      231000 :       gel(v,i) = gerepileuptoint(av2, vi);
    2678             :     }
    2679       52920 :     if (!t) continue;
    2680             :     /* at this point c[t] = 0 */
    2681             : 
    2682       52892 :     if (++rg >= m) { /* full rank; mostly done */
    2683        7595 :       GEN det = gel(v,t); /* last on stack */
    2684        7595 :       if (++k > n)
    2685        7539 :         det = absi(det);
    2686             :       else
    2687             :       {
    2688             :         /* improve further; at this point c[i] is set for all i != t */
    2689          56 :         gcoeff(B,t,t) = piv; v = centermod(gel(B,t), det);
    2690         140 :         for ( ; k<=n; k++)
    2691          84 :           det = gcdii(det, ZV_dotproduct(v, gel(A,k)));
    2692             :       }
    2693        7595 :       return gerepileuptoint(av, det);
    2694             :     }
    2695             : 
    2696       45297 :     piv = gel(v,t);
    2697      401345 :     for (i=1; i<=m; i++)
    2698             :     {
    2699             :       GEN mvi;
    2700      356048 :       if (c[i] || i == t) continue;
    2701             : 
    2702      178024 :       gcoeff(B,t,i) = mvi = negi(gel(v,i));
    2703     1679223 :       for (j=1; j<=m; j++)
    2704     1501199 :         if (c[j]) /* implies j != t */
    2705             :         {
    2706      381717 :           pari_sp av2 = avma;
    2707      381717 :           GEN z = addii(mulii(gcoeff(B,j,i), piv), mulii(gcoeff(B,j,t), mvi));
    2708      381717 :           if (rg > 1) z = diviiexact(z, pivprec);
    2709      381717 :           gcoeff(B,j,i) = gerepileuptoint(av2, z);
    2710             :         }
    2711             :     }
    2712       45297 :     c[t] = k;
    2713       45297 :     if (gc_needed(av,1))
    2714             :     {
    2715           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"detint. k=%ld",k);
    2716           0 :       gerepileall(av1, 2, &piv,&B); v = zerovec(m);
    2717             :     }
    2718             :   }
    2719           7 :   avma = av; return gen_0;
    2720             : }
    2721             : 
    2722             : /* Reduce x modulo (invertible) y */
    2723             : GEN
    2724        6398 : closemodinvertible(GEN x, GEN y)
    2725             : {
    2726        6398 :   return gmul(y, ground(RgM_solve(y,x)));
    2727             : }
    2728             : GEN
    2729           7 : reducemodinvertible(GEN x, GEN y)
    2730             : {
    2731           7 :   return gsub(x, closemodinvertible(x,y));
    2732             : }
    2733             : GEN
    2734           0 : reducemodlll(GEN x,GEN y)
    2735             : {
    2736           0 :   return reducemodinvertible(x, ZM_lll(y, 0.75, LLL_INPLACE));
    2737             : }
    2738             : 
    2739             : /*******************************************************************/
    2740             : /*                                                                 */
    2741             : /*                    KERNEL of an m x n matrix                    */
    2742             : /*          return n - rk(x) linearly independent vectors          */
    2743             : /*                                                                 */
    2744             : /*******************************************************************/
    2745             : static GEN
    2746         126 : deplin_aux(GEN x0)
    2747             : {
    2748         126 :   pari_sp av = avma;
    2749         126 :   long i, j, k, nl, nc = lg(x0)-1;
    2750             :   GEN D, x, y, c, l, d, ck;
    2751             : 
    2752         126 :   if (!nc) { avma=av; return cgetg(1,t_COL); }
    2753          91 :   x = RgM_shallowcopy(x0);
    2754          91 :   nl = nbrows(x);
    2755          91 :   d = const_vec(nl, gen_1); /* pivot list */
    2756          91 :   c = zero_zv(nl);
    2757          91 :   l = cgetg(nc+1, t_VECSMALL); /* not initialized */
    2758          91 :   ck = NULL; /* gcc -Wall */
    2759         399 :   for (k=1; k<=nc; k++)
    2760             :   {
    2761         378 :     ck = gel(x,k);
    2762        1176 :     for (j=1; j<k; j++)
    2763             :     {
    2764         798 :       GEN cj = gel(x,j), piv = gel(d,j), q = gel(ck,l[j]);
    2765        9828 :       for (i=1; i<=nl; i++)
    2766        9030 :         if (i!=l[j]) gel(ck,i) = gsub(gmul(piv, gel(ck,i)), gmul(q, gel(cj,i)));
    2767             :     }
    2768             : 
    2769         378 :     i = gauss_get_pivot_NZ(x, NULL, k, c);
    2770         378 :     if (i > nl) break;
    2771             : 
    2772         308 :     gel(d,k) = gel(ck,i);
    2773         308 :     c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */
    2774             :   }
    2775          91 :   if (k > nc) { avma = av; return cgetg(1,t_COL); }
    2776          70 :   if (k == 1) { avma = av; return scalarcol_shallow(gen_1,nc); }
    2777          70 :   y = cgetg(nc+1,t_COL);
    2778          70 :   gel(y,1) = gcopy(gel(ck, l[1]));
    2779         203 :   for (D=gel(d,1),j=2; j<k; j++)
    2780             :   {
    2781         133 :     gel(y,j) = gmul(gel(ck, l[j]), D);
    2782         133 :     D = gmul(D, gel(d,j));
    2783             :   }
    2784          70 :   gel(y,j) = gneg(D);
    2785          70 :   for (j++; j<=nc; j++) gel(y,j) = gen_0;
    2786          70 :   y = primitive_part(y, &c);
    2787          70 :   return c? gerepileupto(av, y): gerepilecopy(av, y);
    2788             : }
    2789             : static GEN
    2790           0 : RgV_deplin(GEN v)
    2791             : {
    2792           0 :   pari_sp av = avma;
    2793           0 :   long n = lg(v)-1;
    2794           0 :   GEN y, p = NULL;
    2795           0 :   if (n <= 1)
    2796             :   {
    2797           0 :     if (n == 1 && gequal0(gel(v,1))) return mkcol(gen_1);
    2798           0 :     return cgetg(1, t_COL);
    2799             :   }
    2800           0 :   if (gequal0(gel(v,1))) return scalarcol_shallow(gen_1, n);
    2801           0 :   v = primpart(mkvec2(gel(v,1),gel(v,2)));
    2802           0 :   if (RgV_is_FpV(v, &p) && p) v = centerlift(v);
    2803           0 :   y = zerocol(n);
    2804           0 :   gel(y,1) = gneg(gel(v,2));
    2805           0 :   gel(y,2) = gcopy(gel(v,1));
    2806           0 :   return gerepileupto(av, y);
    2807             : 
    2808             : }
    2809             : GEN
    2810         231 : deplin(GEN x)
    2811             : {
    2812         231 :   GEN p = NULL, ff = NULL;
    2813         231 :   switch(typ(x))
    2814             :   {
    2815         231 :     case t_MAT: break;
    2816           0 :     case t_VEC: return RgV_deplin(x);
    2817           0 :     default: pari_err_TYPE("deplin",x);
    2818             :   }
    2819         231 :   if (RgM_is_FpM(x, &p) && p)
    2820             :   {
    2821          63 :     pari_sp av = avma;
    2822             :     ulong pp;
    2823          63 :     x = RgM_Fp_init(x, p, &pp);
    2824          63 :     switch(pp)
    2825             :     {
    2826             :     case 0:
    2827          14 :       x = FpM_ker_gen(x,p,1);
    2828          14 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2829           7 :       x = FpC_center(x,p,shifti(p,-1));
    2830           7 :       break;
    2831             :     case 2:
    2832          14 :       x = F2m_ker_sp(x,1);
    2833          14 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2834           7 :       x = F2c_to_ZC(x); break;
    2835             :     default:
    2836          35 :       x = Flm_ker_sp(x,pp,1);
    2837          35 :       if (!x) { avma = av; return cgetg(1,t_COL); }
    2838          21 :       x = Flv_center(x, pp, pp>>1);
    2839          21 :       x = zc_to_ZC(x);
    2840          21 :       break;
    2841             :     }
    2842          35 :     return gerepileupto(av, x);
    2843             :   }
    2844         168 :   if (RgM_is_FFM(x, &ff))
    2845             :   {
    2846          42 :     x = FFM_deplin(x, ff);
    2847          42 :     if (!x) return cgetg(1, t_COL);
    2848          21 :     return x;
    2849             :   }
    2850         126 :   return deplin_aux(x);
    2851             : }
    2852             : 
    2853             : /*******************************************************************/
    2854             : /*                                                                 */
    2855             : /*         GAUSS REDUCTION OF MATRICES  (m lines x n cols)         */
    2856             : /*           (kernel, image, complementary image, rank)            */
    2857             : /*                                                                 */
    2858             : /*******************************************************************/
    2859             : /* return the transform of x under a standard Gauss pivot.
    2860             :  * x0 is a reference point when guessing whether x[i,j] ~ 0
    2861             :  * (iff x[i,j] << x0[i,j])
    2862             :  * Set r = dim ker(x). d[k] contains the index of the first non-zero pivot
    2863             :  * in column k */
    2864             : static GEN
    2865        8050 : gauss_pivot_ker(GEN x, GEN x0, GEN *dd, long *rr)
    2866             : {
    2867             :   GEN c, d, p, data;
    2868             :   pari_sp av;
    2869             :   long i, j, k, r, t, n, m;
    2870             :   pivot_fun pivot;
    2871             : 
    2872        8050 :   n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
    2873        8015 :   m=nbrows(x); r=0;
    2874        8015 :   pivot = get_pivot_fun(x, x0, &data);
    2875        8015 :   x = RgM_shallowcopy(x);
    2876        8015 :   c = zero_zv(m);
    2877        8015 :   d = cgetg(n+1,t_VECSMALL);
    2878        8015 :   av=avma;
    2879       48216 :   for (k=1; k<=n; k++)
    2880             :   {
    2881       40201 :     j = pivot(x, data, k, c);
    2882       40201 :     if (j > m)
    2883             :     {
    2884       13573 :       r++; d[k]=0;
    2885       74081 :       for(j=1; j<k; j++)
    2886       60508 :         if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
    2887             :     }
    2888             :     else
    2889             :     { /* pivot for column k on row j */
    2890       26628 :       c[j]=k; d[k]=j; p = gdiv(gen_m1,gcoeff(x,j,k));
    2891       26628 :       gcoeff(x,j,k) = gen_m1;
    2892             :       /* x[j,] /= - x[j,k] */
    2893       26628 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
    2894     1817053 :       for (t=1; t<=m; t++)
    2895     1790425 :         if (t!=j)
    2896             :         { /* x[t,] -= 1 / x[j,k] x[j,] */
    2897     1763797 :           p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
    2898    23485175 :           for (i=k+1; i<=n; i++)
    2899    21721378 :             gcoeff(x,t,i) = gadd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
    2900     1763797 :           if (gc_needed(av,1)) gerepile_gauss_ker(x,k,t,av);
    2901             :         }
    2902             :     }
    2903             :   }
    2904        8015 :   *dd=d; *rr=r; return x;
    2905             : }
    2906             : 
    2907             : static GEN FpM_gauss_pivot(GEN x, GEN p, long *rr);
    2908             : static GEN FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr);
    2909             : static GEN F2m_gauss_pivot(GEN x, long *rr);
    2910             : static GEN Flm_gauss_pivot(GEN x, ulong p, long *rry);
    2911             : 
    2912             : /* r = dim ker(x).
    2913             :  * Returns d:
    2914             :  *   d[k] != 0 contains the index of a non-zero pivot in column k
    2915             :  *   d[k] == 0 if column k is a linear combination of the (k-1) first ones */
    2916             : GEN
    2917       15490 : RgM_pivots(GEN x0, GEN data, long *rr, pivot_fun pivot)
    2918             : {
    2919             :   GEN x, c, d, p;
    2920       15490 :   long i, j, k, r, t, m, n = lg(x0)-1;
    2921             :   pari_sp av;
    2922             : 
    2923       15490 :   if (RgM_is_ZM(x0)) return ZM_pivots(x0, rr);
    2924       13516 :   if (!n) { *rr = 0; return NULL; }
    2925             : 
    2926       13516 :   d = cgetg(n+1, t_VECSMALL);
    2927       13516 :   x = RgM_shallowcopy(x0);
    2928       13516 :   m = nbrows(x); r = 0;
    2929       13516 :   c = zero_zv(m);
    2930       13516 :   av = avma;
    2931      885935 :   for (k=1; k<=n; k++)
    2932             :   {
    2933      872419 :     j = pivot(x, data, k, c);
    2934      872419 :     if (j > m) { r++; d[k] = 0; }
    2935             :     else
    2936             :     {
    2937       35042 :       c[j] = k; d[k] = j; p = gdiv(gen_m1, gcoeff(x,j,k));
    2938       35042 :       for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
    2939             : 
    2940      282422 :       for (t=1; t<=m; t++)
    2941      247380 :         if (!c[t]) /* no pivot on that line yet */
    2942             :         {
    2943      130128 :           p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
    2944     6812212 :           for (i=k+1; i<=n; i++)
    2945     6682084 :             gcoeff(x,t,i) = gadd(gcoeff(x,t,i), gmul(p, gcoeff(x,j,i)));
    2946      130128 :           if (gc_needed(av,1)) gerepile_gauss(x,k,t,av,j,c);
    2947             :         }
    2948       35042 :       for (i=k; i<=n; i++) gcoeff(x,j,i) = gen_0; /* dummy */
    2949             :     }
    2950             :   }
    2951       13516 :   *rr = r; avma = (pari_sp)d; return d;
    2952             : }
    2953             : 
    2954             : static long
    2955       99535 : ZM_count_0_cols(GEN M)
    2956             : {
    2957       99535 :   long i, l = lg(M), n = 0;
    2958      504648 :   for (i = 1; i < l; i++)
    2959      405113 :     if (ZV_equal0(gel(M,i))) n++;
    2960       99535 :   return n;
    2961             : }
    2962             : 
    2963             : static void indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol);
    2964             : /* As RgM_pivots, integer entries. Set *rr = dim Ker M0 */
    2965             : GEN
    2966      102760 : ZM_pivots(GEN M0, long *rr)
    2967             : {
    2968      102760 :   GEN d, dbest = NULL;
    2969             :   long m, n, i, imax, rmin, rbest, zc;
    2970      102760 :   int beenthere = 0;
    2971      102760 :   pari_sp av, av0 = avma;
    2972             :   forprime_t S;
    2973             : 
    2974      102760 :   rbest = n = lg(M0)-1;
    2975      102760 :   if (n == 0) { *rr = 0; return NULL; }
    2976       99535 :   zc = ZM_count_0_cols(M0);
    2977       99535 :   if (n == zc) { *rr = zc; return zero_zv(n); }
    2978             : 
    2979       99234 :   m = nbrows(M0);
    2980       99234 :   rmin = maxss(zc, n-m);
    2981       99234 :   init_modular(&S);
    2982       99234 :   imax = (n < (1<<4))? 1: (n>>3); /* heuristic */
    2983             : 
    2984             :   for(;;)
    2985             :   {
    2986             :     GEN row, col, M, KM, IM, RHS, X, cX;
    2987             :     long rk;
    2988      103102 :     for (av = avma, i = 0;; avma = av, i++)
    2989             :     {
    2990      103102 :       ulong p = u_forprime_next(&S);
    2991             :       long rp;
    2992      103102 :       if (!p) pari_err_OVERFLOW("ZM_pivots [ran out of primes]");
    2993      103102 :       d = Flm_gauss_pivot(ZM_to_Flm(M0, p), p, &rp);
    2994      103102 :       if (rp == rmin) { rbest = rp; goto END; } /* maximal rank, return */
    2995        6595 :       if (rp < rbest) { /* save best r so far */
    2996        2762 :         rbest = rp;
    2997        2762 :         if (dbest) gunclone(dbest);
    2998        2762 :         dbest = gclone(d);
    2999        5489 :         if (beenthere) break;
    3000             :       }
    3001        6595 :       if (!beenthere && i >= imax) break;
    3002        3868 :     }
    3003        2727 :     beenthere = 1;
    3004             :     /* Dubious case: there is (probably) a non trivial kernel */
    3005        2727 :     indexrank_all(m,n, rbest, dbest, &row, &col);
    3006        2727 :     M = rowpermute(vecpermute(M0, col), row);
    3007        2727 :     rk = n - rbest; /* (probable) dimension of image */
    3008        2727 :     IM = vecslice(M,1,rk);
    3009        2727 :     KM = vecslice(M,rk+1, n);
    3010        2727 :     M = rowslice(IM, 1,rk); /* square maximal rank */
    3011        2727 :     X = ZM_gauss(M, rowslice(KM, 1,rk));
    3012        2727 :     X = Q_remove_denom(X, &cX);
    3013        2727 :     RHS = rowslice(KM,rk+1,m);
    3014        2727 :     if (cX) RHS = ZM_Z_mul(RHS, cX);
    3015        2727 :     if (ZM_equal(ZM_mul(rowslice(IM,rk+1,m), X), RHS))
    3016             :     {
    3017        2727 :       d = vecsmall_copy(dbest);
    3018        2727 :       goto END;
    3019             :     }
    3020           0 :     avma = av;
    3021           0 :   }
    3022             : END:
    3023       99234 :   *rr = rbest; if (dbest) gunclone(dbest);
    3024       99234 :   return gerepileuptoleaf(av0, d);
    3025             : }
    3026             : 
    3027             : /* set *pr = dim Ker x */
    3028             : static GEN
    3029        7175 : gauss_pivot(GEN x, long *pr) {
    3030             :   GEN data;
    3031        7175 :   pivot_fun pivot = get_pivot_fun(x, x, &data);
    3032        7175 :   return RgM_pivots(x, data, pr, pivot);
    3033             : }
    3034             : 
    3035             : /* compute ker(x), x0 is a reference point when guessing whether x[i,j] ~ 0
    3036             :  * (iff x[i,j] << x0[i,j]) */
    3037             : static GEN
    3038        8050 : ker_aux(GEN x, GEN x0)
    3039             : {
    3040        8050 :   pari_sp av = avma;
    3041             :   GEN d,y;
    3042             :   long i,j,k,r,n;
    3043             : 
    3044        8050 :   x = gauss_pivot_ker(x,x0,&d,&r);
    3045        8050 :   if (!r) { avma=av; return cgetg(1,t_MAT); }
    3046        7693 :   n = lg(x)-1; y=cgetg(r+1,t_MAT);
    3047       21266 :   for (j=k=1; j<=r; j++,k++)
    3048             :   {
    3049       13573 :     GEN p = cgetg(n+1,t_COL);
    3050             : 
    3051       13573 :     gel(y,j) = p; while (d[k]) k++;
    3052       74081 :     for (i=1; i<k; i++)
    3053       60508 :       if (d[i])
    3054             :       {
    3055       31073 :         GEN p1=gcoeff(x,d[i],k);
    3056       31073 :         gel(p,i) = gcopy(p1); gunclone(p1);
    3057             :       }
    3058             :       else
    3059       29435 :         gel(p,i) = gen_0;
    3060       13573 :     gel(p,k) = gen_1; for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
    3061             :   }
    3062        7693 :   return gerepileupto(av,y);
    3063             : }
    3064             : GEN
    3065        7966 : ker(GEN x)
    3066             : {
    3067        7966 :   pari_sp av = avma;
    3068        7966 :   GEN p = NULL, ff = NULL;
    3069        7966 :   if (RgM_is_FpM(x, &p) && p)
    3070             :   {
    3071             :     ulong pp;
    3072          35 :     x = RgM_Fp_init(x, p, &pp);
    3073          35 :     switch(pp)
    3074             :     {
    3075          14 :     case 0: x = FpM_to_mod(FpM_ker_gen(x,p,0),p); break;
    3076           7 :     case 2: x = F2m_to_mod(F2m_ker_sp(x,0)); break;
    3077          14 :     default:x = Flm_to_mod(Flm_ker_sp(x,pp,0), pp); break;
    3078             :     }
    3079          35 :     return gerepileupto(av, x);
    3080             :   }
    3081        7931 :   if (RgM_is_FFM(x, &ff)) return FFM_ker(x, ff);
    3082        7896 :   return ker_aux(x,x);
    3083             : }
    3084             : GEN
    3085         147 : matker0(GEN x,long flag)
    3086             : {
    3087         147 :   if (typ(x)!=t_MAT) pari_err_TYPE("matker",x);
    3088         147 :   if (!flag) return ker(x);
    3089           7 :   RgM_check_ZM(x, "matker");
    3090           7 :   return ZM_ker(x);
    3091             : }
    3092             : 
    3093             : GEN
    3094        2926 : image(GEN x)
    3095             : {
    3096        2926 :   pari_sp av = avma;
    3097        2926 :   GEN d, ff = NULL, p = NULL;
    3098             :   long r;
    3099             : 
    3100        2926 :   if (typ(x)!=t_MAT) pari_err_TYPE("matimage",x);
    3101        2926 :   if (RgM_is_FpM(x, &p) && p)
    3102             :   {
    3103             :     ulong pp;
    3104          35 :     x = RgM_Fp_init(x, p, &pp);
    3105          35 :     switch(pp)
    3106             :     {
    3107          14 :     case 0: x = FpM_to_mod(FpM_image(x,p), p); break;
    3108           7 :     case 2: x = F2m_to_mod(F2m_image(x)); break;
    3109          14 :     default:x = Flm_to_mod(Flm_image(x,pp), pp);
    3110             :     }
    3111          35 :     return gerepileupto(av, x);
    3112             :   }
    3113        2891 :   if (RgM_is_FFM(x, &ff)) return FFM_image(x, ff);
    3114        2870 :   d = gauss_pivot(x,&r); /* d left on stack for efficiency */
    3115        2870 :   return image_from_pivot(x,d,r);
    3116             : }
    3117             : 
    3118             : static GEN
    3119          84 : imagecompl_aux(GEN x, GEN(*PIVOT)(GEN,long*))
    3120             : {
    3121          84 :   pari_sp av = avma;
    3122             :   GEN d,y;
    3123             :   long j,i,r;
    3124             : 
    3125          84 :   if (typ(x)!=t_MAT) pari_err_TYPE("imagecompl",x);
    3126          84 :   (void)new_chunk(lg(x) * 4 + 1); /* HACK */
    3127          84 :   d = PIVOT(x,&r); /* if (!d) then r = 0 */
    3128          84 :   avma = av; y = cgetg(r+1,t_VECSMALL);
    3129         126 :   for (i=j=1; j<=r; i++)
    3130          42 :     if (!d[i]) y[j++] = i;
    3131          84 :   return y;
    3132             : }
    3133             : GEN
    3134          84 : imagecompl(GEN x) { return imagecompl_aux(x, &gauss_pivot); }
    3135             : GEN
    3136           0 : ZM_imagecompl(GEN x) { return imagecompl_aux(x, &ZM_pivots); }
    3137             : 
    3138             : GEN
    3139        1456 : RgM_RgC_invimage(GEN A, GEN y)
    3140             : {
    3141        1456 :   pari_sp av = avma;
    3142        1456 :   long i, l = lg(A);
    3143        1456 :   GEN M, x, t, p = NULL;
    3144             : 
    3145        1456 :   if (RgM_is_FpM(A, &p) && RgV_is_FpV(y, &p) && p)
    3146             :   {
    3147             :     ulong pp;
    3148          28 :     A = RgM_Fp_init(A,p,&pp);
    3149          28 :     switch(pp)
    3150             :     {
    3151             :     case 0:
    3152           7 :       y = RgC_to_FpC(y,p);
    3153           7 :       x = FpM_FpC_invimage(A, y, p);
    3154           7 :       if (x) x = FpC_to_mod(x,p);
    3155           7 :       break;
    3156             :     case 2:
    3157           7 :       y = RgV_to_F2v(y);
    3158           7 :       x = F2m_F2c_invimage(A, y);
    3159           7 :       if (x) x = F2c_to_mod(x);
    3160           7 :       break;
    3161             :     default:
    3162          14 :       y = RgV_to_Flv(y,pp);
    3163          14 :       x = Flm_Flc_invimage(A, y, pp);
    3164          14 :       if (x) x = Flc_to_mod(x,pp);
    3165             :     }
    3166          28 :     if (!x) { avma = av; return NULL; }
    3167          28 :     return gerepileupto(av, x);
    3168             :   }
    3169             : 
    3170        1428 :   if (l==1) return NULL;
    3171        1428 :   if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
    3172             :   {
    3173        1428 :     GEN ff = NULL;
    3174        1428 :     if (RgM_is_FFM(A, &ff) && RgC_is_FFC(y, &ff))
    3175          63 :       return FFM_FFC_invimage(A, y, ff);
    3176             :   }
    3177        1365 :   M = ker(shallowconcat(A, y));
    3178        1365 :   i = lg(M)-1;
    3179        1365 :   if (!i) { avma = av; return NULL; }
    3180             : 
    3181        1365 :   x = gel(M,i); t = gel(x,l);
    3182        1365 :   if (gequal0(t)) { avma = av; return NULL; }
    3183             : 
    3184        1337 :   t = gneg_i(t); setlg(x,l);
    3185        1337 :   return gerepileupto(av, RgC_Rg_div(x, t));
    3186             : }
    3187             : GEN
    3188       39859 : FpM_FpC_invimage(GEN A, GEN y, GEN p)
    3189             : {
    3190       39859 :   pari_sp av = avma;
    3191       39859 :   long i, l = lg(A);
    3192             :   GEN M, x, t;
    3193             : 
    3194       39859 :   if (lgefint(p) == 3)
    3195             :   {
    3196       39852 :     ulong pp = p[2];
    3197       39852 :     A = ZM_to_Flm(A, pp);
    3198       39852 :     y = ZV_to_Flv(y, pp);
    3199       39852 :     x = Flm_Flc_invimage(A,y,pp);
    3200       39852 :     if (!x) { avma = av; return NULL; }
    3201       39852 :     return gerepileupto(av, Flc_to_ZC(x));
    3202             :   }
    3203           7 :   if (l==1) return NULL;
    3204           7 :   if (lg(y) != lgcols(A)) pari_err_DIM("FpM_FpC_invimage");
    3205           7 :   M = FpM_ker(shallowconcat(A,y),p);
    3206           7 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3207             : 
    3208           7 :   x = gel(M,i); t = gel(x,l);
    3209           7 :   if (!signe(t)) { avma = av; return NULL; }
    3210             : 
    3211           7 :   setlg(x,l); t = Fp_inv(negi(t),p);
    3212           7 :   if (is_pm1(t)) return gerepilecopy(av, x);
    3213           7 :   return gerepileupto(av, FpC_Fp_mul(x, t, p));
    3214             : }
    3215             : GEN
    3216       44444 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
    3217             : {
    3218       44444 :   pari_sp av = avma;
    3219       44444 :   long i, l = lg(A);
    3220             :   GEN M, x;
    3221             :   ulong t;
    3222             : 
    3223       44444 :   if (l==1) return NULL;
    3224       44444 :   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
    3225       44444 :   M = cgetg(l+1,t_MAT);
    3226       44444 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    3227       44444 :   gel(M,l) = y; M = Flm_ker(M,p);
    3228       44444 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3229             : 
    3230       44444 :   x = gel(M,i); t = x[l];
    3231       44444 :   if (!t) { avma = av; return NULL; }
    3232             : 
    3233       44444 :   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
    3234       44444 :   if (t!=1) x = Flv_Fl_mul(x, t, p);
    3235       44444 :   return gerepileuptoleaf(av, x);
    3236             : }
    3237             : GEN
    3238          21 : F2m_F2c_invimage(GEN A, GEN y)
    3239             : {
    3240          21 :   pari_sp av = avma;
    3241          21 :   long i, l = lg(A);
    3242             :   GEN M, x;
    3243             : 
    3244          21 :   if (l==1) return NULL;
    3245          21 :   if (lg(y) != lgcols(A)) pari_err_DIM("F2m_F2c_invimage");
    3246          21 :   M = cgetg(l+1,t_MAT);
    3247          21 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    3248          21 :   gel(M,l) = y; M = F2m_ker(M);
    3249          21 :   i = lg(M)-1; if (!i) { avma = av; return NULL; }
    3250             : 
    3251          21 :   x = gel(M,i);
    3252          21 :   if (!F2v_coeff(x,l)) { avma = av; return NULL; }
    3253          21 :   F2v_clear(x, x[1]); x[1]--; /* remove last coord */
    3254          21 :   return gerepileuptoleaf(av, x);
    3255             : }
    3256             : 
    3257             : /* Return X such that m X = v (t_COL or t_MAT), resp. an empty t_COL / t_MAT
    3258             :  * if no solution exist */
    3259             : GEN
    3260        1568 : inverseimage(GEN m, GEN v)
    3261             : {
    3262             :   GEN y;
    3263        1568 :   if (typ(m)!=t_MAT) pari_err_TYPE("inverseimage",m);
    3264        1568 :   switch(typ(v))
    3265             :   {
    3266             :     case t_COL:
    3267        1456 :       y = RgM_RgC_invimage(m,v);
    3268        1456 :       return y? y: cgetg(1,t_COL);
    3269             :     case t_MAT:
    3270         112 :       y = RgM_invimage(m, v);
    3271         112 :       return y? y: cgetg(1,t_MAT);
    3272             :   }
    3273           0 :   pari_err_TYPE("inverseimage",v);
    3274             :   return NULL;/*LCOV_EXCL_LINE*/
    3275             : }
    3276             : 
    3277             : static GEN
    3278          21 : Flm_invimage_i(GEN A, GEN B, ulong p)
    3279             : {
    3280             :   GEN d, x, X, Y;
    3281          21 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3282          21 :   x = Flm_ker_sp(shallowconcat(Flm_neg(A,p), B), p, 0);
    3283             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3284             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3285             :    * Y has at least nB columns and full rank */
    3286          21 :   nY = lg(x)-1;
    3287          21 :   if (nY < nB) return NULL;
    3288          21 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3289          21 :   d = cgetg(nB+1, t_VECSMALL);
    3290          56 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3291             :   {
    3292          49 :     for (; j>=1; j--)
    3293          42 :       if (coeff(Y,i,j)) { d[i] = j; break; }
    3294          42 :     if (!j) return NULL;
    3295             :   }
    3296             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3297          14 :   Y = vecpermute(Y, d);
    3298          14 :   x = vecpermute(x, d);
    3299          14 :   X = rowslice(x, 1, nA);
    3300          14 :   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
    3301             : }
    3302             : 
    3303             : static GEN
    3304           7 : F2m_invimage_i(GEN A, GEN B)
    3305             : {
    3306             :   GEN d, x, X, Y;
    3307           7 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3308           7 :   x = F2m_ker_sp(shallowconcat(A, B), 0);
    3309             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3310             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3311             :    * Y has at least nB columns and full rank */
    3312           7 :   nY = lg(x)-1;
    3313           7 :   if (nY < nB) return NULL;
    3314             : 
    3315             :   /* implicitly: Y = rowslice(x, nA+1, nA+nB), nB rows */
    3316           7 :   d = cgetg(nB+1, t_VECSMALL);
    3317          21 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3318             :   {
    3319          14 :     for (; j>=1; j--)
    3320          14 :       if (F2m_coeff(x,nA+i,j)) { d[i] = j; break; } /* Y[i,j] */
    3321          14 :     if (!j) return NULL;
    3322             :   }
    3323           7 :   x = vecpermute(x, d);
    3324             : 
    3325           7 :   X = F2m_rowslice(x, 1, nA);
    3326           7 :   Y = F2m_rowslice(x, nA+1, nA+nB);
    3327           7 :   return F2m_mul(X, F2m_inv_upper_1(Y));
    3328             : }
    3329             : GEN
    3330           0 : Flm_invimage(GEN A, GEN B, ulong p)
    3331             : {
    3332           0 :   pari_sp av = avma;
    3333           0 :   GEN X = Flm_invimage_i(A,B,p);
    3334           0 :   if (!X) { avma = av; return NULL; }
    3335           0 :   return gerepileupto(av, X);
    3336             : }
    3337             : GEN
    3338           0 : F2m_invimage(GEN A, GEN B)
    3339             : {
    3340           0 :   pari_sp av = avma;
    3341           0 :   GEN X = F2m_invimage_i(A,B);
    3342           0 :   if (!X) { avma = av; return NULL; }
    3343           0 :   return gerepileupto(av, X);
    3344             : }
    3345             : static GEN
    3346          14 : FpM_invimage_i(GEN A, GEN B, GEN p)
    3347             : {
    3348             :   GEN d, x, X, Y;
    3349          14 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3350          14 :   if (lgefint(p) == 3)
    3351             :   {
    3352           0 :     ulong pp = p[2];
    3353           0 :     A = ZM_to_Flm(A, pp);
    3354           0 :     B = ZM_to_Flm(B, pp);
    3355           0 :     x = Flm_invimage_i(A, B, pp);
    3356           0 :     return x? Flm_to_ZM(x): NULL;
    3357             :   }
    3358          14 :   x = FpM_ker(shallowconcat(ZM_neg(A), B), p);
    3359             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3360             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3361             :    * Y has at least nB columns and full rank */
    3362          14 :   nY = lg(x)-1;
    3363          14 :   if (nY < nB) return NULL;
    3364          14 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3365          14 :   d = cgetg(nB+1, t_VECSMALL);
    3366          35 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3367             :   {
    3368          35 :     for (; j>=1; j--)
    3369          28 :       if (signe(gcoeff(Y,i,j))) { d[i] = j; break; }
    3370          28 :     if (!j) return NULL;
    3371             :   }
    3372             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3373           7 :   Y = vecpermute(Y, d);
    3374           7 :   x = vecpermute(x, d);
    3375           7 :   X = rowslice(x, 1, nA);
    3376           7 :   return FpM_mul(X, FpM_inv_upper_1(Y,p), p);
    3377             : }
    3378             : GEN
    3379           0 : FpM_invimage(GEN A, GEN B, GEN p)
    3380             : {
    3381           0 :   pari_sp av = avma;
    3382           0 :   GEN X = FpM_invimage_i(A,B,p);
    3383           0 :   if (!X) { avma = av; return NULL; }
    3384           0 :   return gerepileupto(av, X);
    3385             : }
    3386             : 
    3387             : /* find Z such that A Z = B. Return NULL if no solution */
    3388             : GEN
    3389         126 : RgM_invimage(GEN A, GEN B)
    3390             : {
    3391         126 :   pari_sp av = avma;
    3392             :   GEN d, x, X, Y;
    3393         126 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    3394         126 :   GEN p = NULL;
    3395         126 :   if (RgM_is_FpM(A, &p) && RgM_is_FpM(B, &p) && p)
    3396             :   {
    3397             :     ulong pp;
    3398          42 :     A = RgM_Fp_init(A,p,&pp);
    3399          42 :     switch(pp)
    3400             :     {
    3401             :     case 0:
    3402          14 :       B = RgM_to_FpM(B,p);
    3403          14 :       x = FpM_invimage_i(A, B, p);
    3404          14 :       if (x) x = FpM_to_mod(x, p);
    3405          14 :     break;
    3406             :     case 2:
    3407           7 :       B = RgM_to_F2m(B);
    3408           7 :       x = F2m_invimage_i(A, B);
    3409           7 :       if (x) x = F2m_to_mod(x);
    3410           7 :       break;
    3411             :     default:
    3412          21 :       B = RgM_to_Flm(B,pp);
    3413          21 :       x = Flm_invimage_i(A, B, pp);
    3414          21 :       if (x) x = Flm_to_mod(x,pp);
    3415          21 :       break;
    3416             :     }
    3417          42 :     if (!x) { avma = av; return NULL; }
    3418          28 :     return gerepileupto(av, x);
    3419             :   }
    3420             :   {
    3421          84 :     GEN ff = NULL;
    3422          84 :     if (RgM_is_FFM(A, &ff) && RgM_is_FFM(B, &ff))
    3423          21 :       return FFM_invimage(A, B, ff);
    3424             :   }
    3425          63 :   x = ker(shallowconcat(RgM_neg(A), B));
    3426             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    3427             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    3428             :    * Y has at least nB columns and full rank */
    3429          63 :   nY = lg(x)-1;
    3430          63 :   if (nY < nB) { avma = av; return NULL; }
    3431          63 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    3432          63 :   d = cgetg(nB+1, t_VECSMALL);
    3433         140 :   for (i = nB, j = nY; i >= 1; i--, j--)
    3434             :   {
    3435         182 :     for (; j>=1; j--)
    3436         140 :       if (!gequal0(gcoeff(Y,i,j))) { d[i] = j; break; }
    3437         119 :     if (!j) { avma = av; return NULL; }
    3438             :   }
    3439             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    3440          21 :   Y = vecpermute(Y, d);
    3441          21 :   x = vecpermute(x, d);
    3442          21 :   X = rowslice(x, 1, nA);
    3443          21 :   return gerepileupto(av, RgM_mul(X, RgM_inv_upper(Y)));
    3444             : }
    3445             : 
    3446             : /* r = dim Ker x, n = nbrows(x) */
    3447             : static GEN
    3448       32960 : get_suppl(GEN x, GEN d, long n, long r, GEN(*ei)(long,long))
    3449             : {
    3450             :   pari_sp av;
    3451             :   GEN y, c;
    3452       32960 :   long j, k, rx = lg(x)-1; /* != 0 due to init_suppl() */
    3453             : 
    3454       32960 :   if (rx == n && r == 0) return gcopy(x);
    3455       26186 :   y = cgetg(n+1, t_MAT);
    3456       26186 :   av = avma; c = zero_zv(n);
    3457             :   /* c = lines containing pivots (could get it from gauss_pivot, but cheap)
    3458             :    * In theory r = 0 and d[j] > 0 for all j, but why take chances? */
    3459      208654 :   for (k = j = 1; j<=rx; j++)
    3460      182468 :     if (d[j]) { c[ d[j] ] = 1; gel(y,k++) = gel(x,j); }
    3461      286727 :   for (j=1; j<=n; j++)
    3462      260541 :     if (!c[j]) gel(y,k++) = (GEN)j; /* HACK */
    3463       26186 :   avma = av;
    3464             : 
    3465       26186 :   rx -= r;
    3466       26186 :   for (j=1; j<=rx; j++) gel(y,j) = gcopy(gel(y,j));
    3467       26186 :   for (   ; j<=n; j++)  gel(y,j) = ei(n, y[j]);
    3468       26186 :   return y;
    3469             : }
    3470             : 
    3471             : static void
    3472       32960 : init_suppl(GEN x)
    3473             : {
    3474       32960 :   if (lg(x) == 1) pari_err_IMPL("suppl [empty matrix]");
    3475             :   /* HACK: avoid overwriting d from gauss_pivot() after avma=av */
    3476       32960 :   (void)new_chunk(lgcols(x) * 2);
    3477       32960 : }
    3478             : 
    3479             : /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
    3480             :  * whose first k columns are given by x. If rank(x) < k, undefined result. */
    3481             : GEN
    3482          91 : suppl(GEN x)
    3483             : {
    3484          91 :   pari_sp av = avma;
    3485          91 :   GEN d, X = x, p = NULL, ff = NULL;
    3486             :   long r;
    3487             : 
    3488          91 :   if (typ(x)!=t_MAT) pari_err_TYPE("suppl",x);
    3489          91 :   if (RgM_is_FpM(x, &p) && p)
    3490             :   {
    3491             :     ulong pp;
    3492          28 :     x = RgM_Fp_init(x, p, &pp);
    3493          28 :     switch(pp)
    3494             :     {
    3495           7 :     case 0: x = FpM_to_mod(FpM_suppl(x,p), p); break;
    3496           7 :     case 2: x = F2m_to_mod(F2m_suppl(x)); break;
    3497          14 :     default:x = Flm_to_mod(Flm_suppl(x,pp), pp); break;
    3498             :     }
    3499          28 :     return gerepileupto(av, x);
    3500             :   }
    3501          63 :   if (RgM_is_FFM(x, &ff)) return FFM_suppl(x, ff);
    3502          42 :   avma = av; init_suppl(x);
    3503          42 :   d = gauss_pivot(X,&r);
    3504          42 :   avma = av; return get_suppl(X,d,nbrows(X),r,&col_ei);
    3505             : }
    3506             : GEN
    3507       32750 : FpM_suppl(GEN x, GEN p)
    3508             : {
    3509       32750 :   pari_sp av = avma;
    3510             :   GEN d;
    3511             :   long r;
    3512       32750 :   init_suppl(x); d = FpM_gauss_pivot(x,p, &r);
    3513       32750 :   avma = av; return get_suppl(x,d,nbrows(x),r,&col_ei);
    3514             : }
    3515             : GEN
    3516          42 : Flm_suppl(GEN x, ulong p)
    3517             : {
    3518          42 :   pari_sp av = avma;
    3519             :   GEN d;
    3520             :   long r;
    3521          42 :   init_suppl(x); d = Flm_gauss_pivot(Flm_copy(x),p, &r);
    3522          42 :   avma = av; return get_suppl(x,d,nbrows(x),r,&vecsmall_ei);
    3523             : 
    3524             : }
    3525             : GEN
    3526           7 : F2m_suppl(GEN x)
    3527             : {
    3528           7 :   pari_sp av = avma;
    3529             :   GEN d;
    3530             :   long r;
    3531           7 :   init_suppl(x); d = F2m_gauss_pivot(F2m_copy(x), &r);
    3532           7 :   avma = av; return get_suppl(x,d,mael(x,1,1),r,&F2v_ei);
    3533             : }
    3534             : 
    3535             : /* variable number to be filled in later */
    3536             : static GEN
    3537          14 : _FlxC_ei(long n, long i)
    3538             : {
    3539          14 :   GEN x = cgetg(n + 1, t_COL);
    3540             :   long j;
    3541          42 :   for (j = 1; j <= n; j++)
    3542          28 :     gel(x, j) = (j == i)? pol1_Flx(0): pol0_Flx(0);
    3543          14 :   return x;
    3544             : }
    3545             : 
    3546             : GEN
    3547           7 : F2xqM_suppl(GEN x, GEN T)
    3548             : {
    3549           7 :   pari_sp av = avma;
    3550             :   GEN d, y;
    3551           7 :   long n = nbrows(x), r, sv = get_Flx_var(T);
    3552             : 
    3553           7 :   init_suppl(x);
    3554           7 :   d = F2xqM_gauss_pivot(x, T, &r);
    3555           7 :   avma = av;
    3556           7 :   y = get_suppl(x, d, n, r, &_FlxC_ei);
    3557           7 :   if (sv) {
    3558             :     long i, j;
    3559          21 :     for (j = r + 1; j <= n; j++) {
    3560          42 :       for (i = 1; i <= n; i++)
    3561          28 :         gcoeff(y, i, j)[1] = sv;
    3562             :     }
    3563             :   }
    3564           7 :   return y;
    3565             : }
    3566             : 
    3567             : GEN
    3568           7 : FlxqM_suppl(GEN x, GEN T, ulong p)
    3569             : {
    3570           7 :   pari_sp av = avma;
    3571             :   GEN d, y;
    3572           7 :   long n = nbrows(x), r, sv = get_Flx_var(T);
    3573             : 
    3574           7 :   init_suppl(x);
    3575           7 :   d = FlxqM_gauss_pivot(x, T, p, &r);
    3576           7 :   avma = av;
    3577           7 :   y = get_suppl(x, d, n, r, &_FlxC_ei);
    3578           7 :   if (sv) {
    3579             :     long i, j;
    3580          21 :     for (j = r + 1; j <= n; j++) {
    3581          42 :       for (i = 1; i <= n; i++)
    3582          28 :         gcoeff(y, i, j)[1] = sv;
    3583             :     }
    3584             :   }
    3585           7 :   return y;
    3586             : }
    3587             : 
    3588             : GEN
    3589         581 : FqM_suppl(GEN x, GEN T, GEN p)
    3590             : {
    3591         581 :   pari_sp av = avma;
    3592             :   GEN d;
    3593             :   long r;
    3594             : 
    3595         581 :   if (!T) return FpM_suppl(x,p);
    3596         105 :   init_suppl(x);
    3597         105 :   d = FqM_gauss_pivot(x,T,p,&r);
    3598         105 :   avma = av; return get_suppl(x,d,nbrows(x),r,&col_ei);
    3599             : }
    3600             : 
    3601             : GEN
    3602           7 : image2(GEN x)
    3603             : {
    3604           7 :   pari_sp av = avma;
    3605             :   long k, n, i;
    3606             :   GEN A, B;
    3607             : 
    3608           7 :   if (typ(x)!=t_MAT) pari_err_TYPE("image2",x);
    3609           7 :   if (lg(x) == 1) return cgetg(1,t_MAT);
    3610           7 :   A = ker(x); k = lg(A)-1;
    3611           7 :   if (!k) { avma = av; return gcopy(x); }
    3612           7 :   A = suppl(A); n = lg(A)-1;
    3613           7 :   B = cgetg(n-k+1, t_MAT);
    3614           7 :   for (i = k+1; i <= n; i++) gel(B,i-k) = RgM_RgC_mul(x, gel(A,i));
    3615           7 :   return gerepileupto(av, B);
    3616             : }
    3617             : 
    3618             : GEN
    3619         119 : matimage0(GEN x,long flag)
    3620             : {
    3621         119 :   switch(flag)
    3622             :   {
    3623         112 :     case 0: return image(x);
    3624           7 :     case 1: return image2(x);
    3625           0 :     default: pari_err_FLAG("matimage");
    3626             :   }
    3627             :   return NULL; /* LCOV_EXCL_LINE */
    3628             : }
    3629             : 
    3630             : long
    3631         154 : rank(GEN x)
    3632             : {
    3633         154 :   pari_sp av = avma;
    3634             :   long r;
    3635         154 :   GEN ff = NULL, p = NULL;
    3636             : 
    3637         154 :   if (typ(x)!=t_MAT) pari_err_TYPE("rank",x);
    3638         154 :   if (RgM_is_FpM(x, &p) && p)
    3639             :   {
    3640             :     ulong pp;
    3641          84 :     x = RgM_Fp_init(x,p,&pp);
    3642          84 :     switch(pp)
    3643             :     {
    3644           7 :     case 0: r = FpM_rank(x,p); break;
    3645          63 :     case 2: r = F2m_rank(x); break;
    3646          14 :     default:r = Flm_rank(x,pp); break;
    3647             :     }
    3648          84 :     avma = av; return r;
    3649             :   }
    3650          70 :   if (RgM_is_FFM(x, &ff)) return FFM_rank(x, ff);
    3651          49 :   (void)gauss_pivot(x, &r);
    3652          49 :   avma = av; return lg(x)-1 - r;
    3653             : }
    3654             : 
    3655             : /* d a t_VECSMALL of integers in 1..n. Return the vector of the d[i]
    3656             :  * followed by the missing indices */
    3657             : static GEN
    3658        5454 : perm_complete(GEN d, long n)
    3659             : {
    3660        5454 :   GEN y = cgetg(n+1, t_VECSMALL);
    3661        5454 :   long i, j = 1, k = n, l = lg(d);
    3662        5454 :   pari_sp av = avma;
    3663        5454 :   char *T = stack_calloc(n+1);
    3664        5454 :   for (i = 1; i < l; i++) T[d[i]] = 1;
    3665       68290 :   for (i = 1; i <= n; i++)
    3666       62836 :     if (T[i]) y[j++] = i; else y[k--] = i;
    3667        5454 :   avma = av; return y;
    3668             : }
    3669             : 
    3670             : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
    3671             : static GEN
    3672       20619 : indexrank0(long n, long r, GEN d)
    3673             : {
    3674       20619 :   GEN p1, p2, res = cgetg(3,t_VEC);
    3675             :   long i, j;
    3676             : 
    3677       20619 :   r = n - r; /* now r = dim Im(x) */
    3678       20619 :   p1 = cgetg(r+1,t_VECSMALL); gel(res,1) = p1;
    3679       20619 :   p2 = cgetg(r+1,t_VECSMALL); gel(res,2) = p2;
    3680       20619 :   if (d)
    3681             :   {
    3682      117457 :     for (i=0,j=1; j<=n; j++)
    3683       97566 :       if (d[j]) { i++; p1[i] = d[j]; p2[i] = j; }
    3684       19891 :     vecsmall_sort(p1);
    3685             :   }
    3686       20619 :   return res;
    3687             : }
    3688             : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
    3689             : static GEN
    3690         840 : indeximage0(long n, long r, GEN d)
    3691             : {
    3692             :   long i, j;
    3693             :   GEN v;
    3694             : 
    3695         840 :   r = n - r; /* now r = dim Im(x) */
    3696         840 :   v = cgetg(r+1,t_VECSMALL);
    3697        7651 :   if (d) for (i=j=1; j<=n; j++)
    3698        6811 :     if (d[j]) v[i++] = j;
    3699         840 :   return v;
    3700             : }
    3701             : /* x an m x n t_MAT, n > 0, r = dim Ker(x), d from gauss_pivot */
    3702             : static void
    3703        2727 : indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol)
    3704             : {
    3705        2727 :   GEN IR = indexrank0(n, r, d);
    3706        2727 :   *prow = perm_complete(gel(IR,1), m);
    3707        2727 :   *pcol = perm_complete(gel(IR,2), n);
    3708        2727 : }
    3709             : static void
    3710       18732 : init_indexrank(GEN x) {
    3711       18732 :   (void)new_chunk(3 + 2*lg(x)); /* HACK */
    3712       18732 : }
    3713             : 
    3714             : GEN
    3715        4179 : indexrank(GEN x) {
    3716             :   pari_sp av;
    3717             :   long r;
    3718        4179 :   GEN d, p = NULL, ff = NULL;
    3719        4179 :   if (typ(x)!=t_MAT) pari_err_TYPE("indexrank",x);
    3720        4179 :   if (RgM_is_FpM(x, &p) && p)
    3721             :   {
    3722             :     ulong pp;
    3723          28 :     x = RgM_Fp_init(x,p,&pp);
    3724          28 :     switch(pp)
    3725             :     {
    3726           7 :     case 0:  return FpM_indexrank(x, p);
    3727           7 :     case 2:  return F2m_indexrank(x);
    3728          14 :     default: return Flm_indexrank(x, pp);
    3729             :     }
    3730             :   }
    3731        4151 :   if (RgM_is_FFM(x, &ff)) return FFM_indexrank(x, ff);
    3732        4130 :   av = avma;
    3733        4130 :   init_indexrank(x);
    3734        4130 :   d = gauss_pivot(x, &r);
    3735        4130 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3736             : }
    3737             : 
    3738             : GEN
    3739          28 : FpM_indexrank(GEN x, GEN p) {
    3740          28 :   pari_sp av = avma;
    3741             :   long r;
    3742             :   GEN d;
    3743          28 :   init_indexrank(x);
    3744          28 :   d = FpM_gauss_pivot(x,p,&r);
    3745          28 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3746             : }
    3747             : GEN
    3748       10052 : Flm_indexrank(GEN x, ulong p) {
    3749       10052 :   pari_sp av = avma;
    3750             :   long r;
    3751             :   GEN d;
    3752       10052 :   init_indexrank(x);
    3753       10052 :   d = Flm_gauss_pivot(Flm_copy(x),p,&r);
    3754       10052 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3755             : }
    3756             : GEN
    3757           7 : F2m_indexrank(GEN x) {
    3758           7 :   pari_sp av = avma;
    3759             :   long r;
    3760             :   GEN d;
    3761           7 :   init_indexrank(x);
    3762           7 :   d = F2m_gauss_pivot(F2m_copy(x),&r);
    3763           7 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3764             : }
    3765             : 
    3766             : GEN
    3767           7 : F2xqM_indexrank(GEN x, GEN T) {
    3768           7 :   pari_sp av = avma;
    3769             :   long r;
    3770             :   GEN d;
    3771           7 :   init_indexrank(x);
    3772           7 :   d = F2xqM_gauss_pivot(x, T, &r);
    3773           7 :   avma = av; return indexrank0(lg(x) - 1, r, d);
    3774             : }
    3775             : 
    3776             : GEN
    3777           7 : FlxqM_indexrank(GEN x, GEN T, ulong p) {
    3778           7 :   pari_sp av = avma;
    3779             :   long r;
    3780             :   GEN d;
    3781           7 :   init_indexrank(x);
    3782           7 :   d = FlxqM_gauss_pivot(x, T, p, &r);
    3783           7 :   avma = av; return indexrank0(lg(x) - 1, r, d);
    3784             : }
    3785             : 
    3786             : GEN
    3787           7 : FqM_indexrank(GEN x, GEN T, GEN p) {
    3788           7 :   pari_sp av = avma;
    3789             :   long r;
    3790             :   GEN d;
    3791           7 :   init_indexrank(x);
    3792           7 :   d = FqM_gauss_pivot(x, T, p, &r);
    3793           7 :   avma = av; return indexrank0(lg(x) - 1, r, d);
    3794             : }
    3795             : 
    3796             : GEN
    3797         840 : ZM_indeximage(GEN x) {
    3798         840 :   pari_sp av = avma;
    3799             :   long r;
    3800             :   GEN d;
    3801         840 :   init_indexrank(x);
    3802         840 :   d = ZM_pivots(x,&r);
    3803         840 :   avma = av; return indeximage0(lg(x)-1, r, d);
    3804             : }
    3805             : long
    3806       34046 : ZM_rank(GEN x) {
    3807       34046 :   pari_sp av = avma;
    3808             :   long r;
    3809       34046 :   (void)ZM_pivots(x,&r);
    3810       34046 :   avma = av; return lg(x)-1-r;
    3811             : }
    3812             : GEN
    3813        3654 : ZM_indexrank(GEN x) {
    3814        3654 :   pari_sp av = avma;
    3815             :   long r;
    3816             :   GEN d;
    3817        3654 :   init_indexrank(x);
    3818        3654 :   d = ZM_pivots(x,&r);
    3819        3654 :   avma = av; return indexrank0(lg(x)-1, r, d);
    3820             : }
    3821             : 
    3822             : /*******************************************************************/
    3823             : /*                                                                 */
    3824             : /*                             ZabM                                */
    3825             : /*                                                                 */
    3826             : /*******************************************************************/
    3827             : 
    3828             : static GEN
    3829           0 : FpVM_ratlift(GEN a, GEN q, long vs)
    3830             : {
    3831             :   GEN B, y;
    3832           0 :   long i, j, l = lg(a), n;
    3833           0 :   B = sqrti(shifti(q,-1));
    3834           0 :   y = cgetg(l, t_MAT);
    3835           0 :   if (l==1) return y;
    3836           0 :   n = lgcols(a);
    3837           0 :   for (i=1; i<l; i++)
    3838             :   {
    3839           0 :     GEN yi = cgetg(n, t_COL);
    3840           0 :     for (j=1; j<n; j++)
    3841             :     {
    3842           0 :       GEN v = FpC_ratlift(gmael(a,i,j), q, B, B, NULL);
    3843           0 :       if (!v) return NULL;
    3844           0 :       gel(yi, j) = RgV_to_RgX(v, vs);
    3845             :     }
    3846           0 :     gel(y,i) = yi;
    3847             :   }
    3848           0 :   return y;
    3849             : }
    3850             : 
    3851             : static GEN
    3852           0 : FlmV_recover(GEN a, GEN M, ulong p)
    3853             : {
    3854           0 :   GEN a1 = gel(a,1);
    3855           0 :   long i, j, k, l = lg(a1), n, lM = lg(M);
    3856           0 :   GEN y = cgetg(l, t_MAT);
    3857           0 :   if (l==1) return y;
    3858           0 :   n = lgcols(a1);
    3859           0 :   for (i=1; i<l; i++)
    3860             :   {
    3861           0 :     GEN yi = cgetg(n, t_COL);
    3862           0 :     for (j=1; j<n; j++)
    3863             :     {
    3864           0 :       GEN v = cgetg(lM, t_COL);
    3865           0 :       for (k=1; k<lM; k++) gel(v,k) = gmael(gel(a,k),i,j);
    3866           0 :       gel(yi, j) = Flm_Flc_mul(M, v, p);
    3867             :     }
    3868           0 :     gel(y,i) = yi;
    3869             :   }
    3870           0 :   return y;
    3871             : }
    3872             : 
    3873             : static GEN
    3874           0 : FlkM_inv(GEN M, GEN P, ulong p)
    3875             : {
    3876           0 :   ulong pi = get_Fl_red(p);
    3877           0 :   GEN R = Flx_roots(P, p);
    3878           0 :   long l = lg(R), i;
    3879           0 :   GEN W = Flv_invVandermonde(R, 1UL, p);
    3880           0 :   GEN V = cgetg(l, t_VEC);
    3881           0 :   for(i=1; i<l; i++)
    3882             :   {
    3883           0 :     gel(V, i) = Flm_inv(FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,i), degpol(P), p, pi), p, pi), p);
    3884           0 :     if (!gel(V, i)) return NULL;
    3885             :   }
    3886           0 :   return FlmV_recover(V,W,p);
    3887             : }
    3888             : 
    3889             : GEN
    3890           0 : ZabM_inv(GEN M, GEN P, long n, GEN *pden)
    3891             : {
    3892           0 :   pari_sp av2, av = avma;
    3893             :   GEN q, H;
    3894           0 :   ulong m = LONG_MAX>>1;
    3895           0 :   ulong p= 1 + m - (m % n);
    3896           0 :   long lM = lg(M);
    3897           0 :   if (lM == 1) { *pden = gen_1; return cgetg(1,t_MAT); }
    3898             : 
    3899           0 :   av2 = avma;
    3900           0 :   H = NULL;
    3901             :   for(;;)
    3902             :   {
    3903             :     GEN Hp, Pp, Mp, Hr;
    3904           0 :     do p += n; while(!uisprime(p));
    3905           0 :     Pp = ZX_to_Flx(P, p);
    3906           0 :     Mp = FqM_to_FlxM(M, P, utoi(p));
    3907           0 :     Hp = FlkM_inv(Mp, Pp, p);
    3908           0 :     if (!Hp) continue;
    3909           0 :     if (!H)
    3910             :     {
    3911           0 :       H = ZVM_init_CRT(Hp, p);
    3912           0 :       q = utoipos(p);
    3913             :     }
    3914             :     else
    3915           0 :       ZVM_incremental_CRT(&H, Hp, &q, p);
    3916           0 :     Hr = FpVM_ratlift(H, q, varn(P));
    3917           0 :     if (DEBUGLEVEL>5) err_printf("ZabM_inv mod %ld (ratlift=%ld)\n", p,!!Hr);
    3918           0 :     if (Hr) {/* DONE ? */
    3919           0 :       GEN Hl = Q_remove_denom(Hr, pden);
    3920           0 :       GEN MH = RgXQM_mul(M, Hl,P);
    3921           0 :       if (*pden)
    3922           0 :       { if (RgM_isscalar(MH, *pden)) { H = Hl; break; }}
    3923             :       else
    3924           0 :       { if (RgM_isidentity(MH)) { H = Hl; *pden = gen_1; break; } }
    3925             :     }
    3926             : 
    3927           0 :     if (gc_needed(av,2))
    3928             :     {
    3929           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_inv");
    3930           0 :       gerepileall(av2, 2, &H, &q);
    3931             :     }
    3932           0 :   }
    3933           0 :   gerepileall(av, 2, &H, pden);
    3934           0 :   return H;
    3935             : }
    3936             : 
    3937             : static GEN
    3938           0 : FlkM_ker(GEN M, GEN P, ulong p)
    3939             : {
    3940           0 :   ulong pi = get_Fl_red(p);
    3941           0 :   GEN R = Flx_roots(P, p);
    3942           0 :   long l = lg(R), i, dP = degpol(P), r;
    3943             :   GEN M1, K, D;
    3944           0 :   GEN W = Flv_invVandermonde(R, 1UL, p);
    3945           0 :   GEN V = cgetg(l, t_VEC);
    3946           0 :   M1 = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,1), dP, p, pi), p, pi);
    3947           0 :   K = Flm_ker_sp(M1, p, 2);
    3948           0 :   r = lg(gel(K,1)); D = gel(K,2);
    3949           0 :   gel(V, 1) = gel(K,1);
    3950           0 :   for(i=2; i<l; i++)
    3951             :   {
    3952           0 :     GEN Mi = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,i), dP, p, pi), p, pi);
    3953           0 :     GEN K = Flm_ker_sp(Mi, p, 2);
    3954           0 :     if (lg(gel(K,1)) != r || !zv_equal(D, gel(K,2))) return NULL;
    3955           0 :     gel(V, i) = gel(K,1);
    3956             :   }
    3957           0 :   return mkvec2(FlmV_recover(V,W,p), D);
    3958             : }
    3959             : 
    3960             : GEN
    3961           0 : ZabM_ker(GEN M, GEN P, long n)
    3962             : {
    3963           0 :   pari_sp av2, av = avma;
    3964             :   GEN q, H, D;
    3965           0 :   ulong m = LONG_MAX>>1;
    3966           0 :   ulong p= 1 + m - (m % n);
    3967           0 :   av2 = avma;
    3968           0 :   H = NULL; D = NULL;
    3969             :   for(;;)
    3970             :   {
    3971             :     GEN Kp, Hp, Dp, Pp, Mp, Hr;
    3972           0 :     do p += n; while(!uisprime(p));
    3973           0 :     Pp = ZX_to_Flx(P, p);
    3974           0 :     Mp = FqM_to_FlxM(M, P, utoi(p));
    3975           0 :     Kp = FlkM_ker(Mp, Pp, p);
    3976           0 :     if (!Kp) continue;
    3977           0 :     Hp = gel(Kp,1); Dp = gel(Kp,2);
    3978           0 :     if (H && (lg(Hp)>lg(H) || (lg(Hp)==lg(H) && vecsmall_lexcmp(Dp,D)>0))) continue;
    3979           0 :     if (!H || (lg(Hp)<lg(H) || vecsmall_lexcmp(Dp,D)<0))
    3980             :     {
    3981           0 :       H = ZVM_init_CRT(Hp, p); D = Dp;
    3982           0 :       q = utoipos(p);
    3983             :     }
    3984             :     else
    3985           0 :       ZVM_incremental_CRT(&H, Hp, &q, p);
    3986           0 :     Hr = FpVM_ratlift(H, q, varn(P));
    3987           0 :     if (DEBUGLEVEL>5) err_printf("ZabM_ker mod %ld (ratlift=%ld)\n", p,!!Hr);
    3988           0 :     if (Hr) {/* DONE ? */
    3989           0 :       GEN Hl = Q_remove_denom(Hr, NULL);
    3990           0 :       GEN MH = RgXQM_mul(M, Hl,P);
    3991           0 :       if (gequal0(MH)) { H = Hl;  break; }
    3992             :     }
    3993             : 
    3994           0 :     if (gc_needed(av,2))
    3995             :     {
    3996           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_ker");
    3997           0 :       gerepileall(av2, 3, &H, &D, &q);
    3998             :     }
    3999           0 :   }
    4000           0 :   return gerepilecopy(av, H);
    4001             : }
    4002             : 
    4003             : GEN
    4004           0 : ZabM_indexrank(GEN M, GEN P, long n)
    4005             : {
    4006           0 :   pari_sp av = avma;
    4007           0 :   ulong m = LONG_MAX>>1;
    4008           0 :   ulong p = 1+m-(m%n), D = degpol(P);
    4009             :   GEN v;
    4010             :   do
    4011             :   {
    4012             :     ulong pi;
    4013             :     GEN R, Mp, K;
    4014           0 :     do p += n; while (!uisprime(p));
    4015           0 :     pi = get_Fl_red(p);
    4016           0 :     R = Flx_roots(ZX_to_Flx(P, p), p);
    4017           0 :     Mp = FqM_to_FlxM(M, P, utoipos(p));
    4018           0 :     K = FlxM_eval_powers_pre(Mp, Fl_powers_pre(uel(R,1), D,p,pi), p,pi);
    4019           0 :     v = Flm_indexrank(K, p);
    4020           0 :   } while (lg(gel(v,2)) < lg(M));
    4021           0 :   return gerepileupto(av, v);
    4022             : }
    4023             : 
    4024             : #if 0
    4025             : GEN
    4026             : ZabM_gauss(GEN M, GEN P, long n, GEN *den)
    4027             : {
    4028             :   pari_sp av = avma;
    4029             :   GEN v, S, W;
    4030             :   v = ZabM_indexrank(M, P, n);
    4031             :   S = shallowmatextract(M,gel(v,1),gel(v,2));
    4032             :   W = ZabM_inv(S, P, n, den);
    4033             :   gerepileall(av,2,&W,den);
    4034             :   return W;
    4035             : }
    4036             : #endif
    4037             : 
    4038             : GEN
    4039           0 : ZabM_pseudoinv(GEN M, GEN P, long n, GEN *den)
    4040             : {
    4041           0 :   pari_sp av = avma;
    4042             :   GEN v, S, W, z, v1;
    4043             :   long l, i;
    4044           0 :   v = ZabM_indexrank(M, P, n);
    4045           0 :   S = shallowmatextract(M,gel(v,1),gel(v,2));
    4046           0 :   W = ZabM_inv(S, P, n, den);
    4047           0 :   z = zeromatcopy(lg(M)-1,lgcols(M)-1);
    4048           0 :   v1 = gel(v,1); l = lg(v1);
    4049           0 :   for(i=1; i<l; i++) gel(z, v1[i]) = gel(W, i);
    4050           0 :   gerepileall(av,2,&z,den);
    4051           0 :   return z;
    4052             : }
    4053             : 
    4054             : GEN
    4055           0 : ZM_pseudoinv(GEN M, GEN *den)
    4056             : {
    4057           0 :   pari_sp av = avma;
    4058             :   GEN v, S, W, z, v1;
    4059             :   long l, i;
    4060           0 :   v = ZM_indexrank(M);
    4061           0 :   S = shallowmatextract(M,gel(v,1),gel(v,2));
    4062           0 :   W = ZM_inv_ratlift(S, den);
    4063           0 :   z = zeromatcopy(lg(M)-1,lgcols(M)-1);
    4064           0 :   v1 = gel(v,1); l = lg(v1);
    4065           0 :   for(i=1; i<l; i++) gel(z, v1[i]) = gel(W, i);
    4066           0 :   gerepileall(av,2,&z,den);
    4067           0 :   return z;
    4068             : }
    4069             : 
    4070             : /*******************************************************************/
    4071             : /*                                                                 */
    4072             : /*                   Structured Elimination                        */
    4073             : /*                                                                 */
    4074             : /*******************************************************************/
    4075             : 
    4076             : static void
    4077       94143 : rem_col(GEN c, long i, GEN iscol, GEN Wrow, long *rcol, long *rrow)
    4078             : {
    4079       94143 :   long lc = lg(c), k;
    4080       94143 :   iscol[i] = 0; (*rcol)--;
    4081      832818 :   for (k = 1; k < lc; ++k)
    4082             :   {
    4083      738675 :     Wrow[c[k]]--;
    4084      738675 :     if (Wrow[c[k]]==0) (*rrow)--;
    4085             :   }
    4086       94143 : }
    4087             : 
    4088             : static void
    4089        5632 : rem_singleton(GEN M, GEN iscol, GEN Wrow, long *rcol, long *rrow)
    4090             : {
    4091             :   long i, j;
    4092        5632 :   long nbcol = lg(iscol)-1, last;
    4093             :   do
    4094             :   {
    4095        7535 :     last = 0;
    4096    16320651 :     for (i = 1; i <= nbcol; ++i)
    4097    16313116 :       if (iscol[i])
    4098             :       {
    4099     8812408 :         GEN c = gmael(M, i, 1);
    4100     8812408 :         long lc = lg(c);
    4101    81915276 :         for (j = 1; j < lc; ++j)
    4102    73115559 :           if (Wrow[c[j]] == 1)
    4103             :           {
    4104       12691 :             rem_col(c, i, iscol, Wrow, rcol, rrow);
    4105       12691 :             last=1; break;
    4106             :           }
    4107             :       }
    4108        7535 :   } while (last);
    4109        5632 : }
    4110             : 
    4111             : static GEN
    4112        5510 : fill_wcol(GEN M, GEN iscol, GEN Wrow, long *w, GEN wcol)
    4113             : {
    4114        5510 :   long nbcol = lg(iscol)-1;
    4115             :   long i, j, m, last;
    4116             :   GEN per;
    4117       13320 :   for (m = 2, last=0; !last ; m++)
    4118             :   {
    4119    17483777 :     for (i = 1; i <= nbcol; ++i)
    4120             :     {
    4121    17475967 :       wcol[i] = 0;
    4122    17475967 :       if (iscol[i])
    4123             :       {
    4124     9388007 :         GEN c = gmael(M, i, 1);
    4125     9388007 :         long lc = lg(c);
    4126    87122850 :         for (j = 1; j < lc; ++j)
    4127    77734843 :           if (Wrow[c[j]] == m) {  wcol[i]++; last = 1; }
    4128             :       }
    4129             :     }
    4130             :   }
    4131        5510 :   per = vecsmall_indexsort(wcol);
    4132        5510 :   *w = wcol[per[nbcol]];
    4133        5510 :   return per;
    4134             : }
    4135             : 
    4136             : /* M is a RgMs with nbrow rows, A a list of row indices.
    4137             :    Eliminate rows of M with a single entry that do not belong to A,
    4138             :    and the corresponding columns. Also eliminate columns until #colums=#rows.
    4139             :    Return pcol and prow:
    4140             :    pcol is a map from the new columns indices to the old one.
    4141             :    prow is a map from the old rows indices to the new one (0 if removed).
    4142             : */
    4143             : 
    4144             : void
    4145         122 : RgMs_structelim_col(GEN M, long nbcol, long nbrow, GEN A, GEN *p_col, GEN *p_row)
    4146             : {
    4147             :   long i,j,k;
    4148         122 :   long lA = lg(A);
    4149         122 :   GEN prow = cgetg(nbrow+1, t_VECSMALL);
    4150         122 :   GEN pcol = zero_zv(nbcol);
    4151         122 :   pari_sp av = avma;
    4152         122 :   long rcol = nbcol, rrow = 0, imin = nbcol - usqrt(nbcol);
    4153         122 :   GEN iscol = const_vecsmall(nbcol, 1);
    4154         122 :   GEN Wrow  = zero_zv(nbrow);
    4155         122 :   GEN wcol = cgetg(nbcol+1, t_VECSMALL);
    4156         122 :   pari_sp av2=avma;
    4157      120085 :   for (i = 1; i <= nbcol; ++i)
    4158             :   {
    4159      119963 :     GEN F = gmael(M, i, 1);
    4160      119963 :     long l = lg(F)-1;
    4161     1059728 :     for (j = 1; j <= l; ++j)
    4162      939765 :       Wrow[F[j]]++;
    4163             :   }
    4164         122 :   for (j = 1; j < lA; ++j)
    4165             :   {
    4166         122 :     if (Wrow[A[j]] == 0) { *p_col=NULL; return; }
    4167           0 :     Wrow[A[j]] = -1;
    4168             :   }
    4169      221992 :   for (i = 1; i <= nbrow; ++i)
    4170      221870 :     if (Wrow[i])
    4171       65679 :       rrow++;
    4172         122 :   rem_singleton(M, iscol, Wrow, &rcol, &rrow);
    4173         122 :   if (rcol<rrow) pari_err_BUG("RgMs_structelim, rcol<rrow");
    4174        5754 :   for (; rcol>rrow;)
    4175             :   {
    4176             :     long w;
    4177        5510 :     GEN per = fill_wcol(M, iscol, Wrow, &w, wcol);
    4178       86962 :     for (i = nbcol; i>=imin && wcol[per[i]]>=w && rcol>rrow; i--)
    4179       81452 :       rem_col(gmael(M, per[i], 1), per[i], iscol, Wrow, &rcol, &rrow);
    4180        5510 :     rem_singleton(M, iscol, Wrow, &rcol, &rrow);
    4181        5510 :     avma = av2;
    4182             :   }
    4183      120085 :   for (j = 1, i = 1; i <= nbcol; ++i)
    4184      119963 :     if (iscol[i])
    4185       25820 :       pcol[j++] = i;
    4186         122 :   setlg(pcol,j);
    4187      221992 :   for (k = 1, i = 1; i <= nbrow; ++i)
    4188      221870 :     prow[i] = Wrow[i] ? k++: 0;
    4189         122 :   avma = av;
    4190         122 :   *p_col = pcol; *p_row = prow;
    4191             : }
    4192             : 
    4193             : void
    4194           0 : RgMs_structelim(GEN M, long nbrow, GEN A, GEN *p_col, GEN *p_row)
    4195             : {
    4196           0 :   RgMs_structelim_col(M, lg(M)-1, nbrow, A, p_col, p_row);
    4197           0 : }
    4198             : 
    4199             : /*******************************************************************/
    4200             : /*                                                                 */
    4201             : /*                        EIGENVECTORS                             */
    4202             : /*   (independent eigenvectors, sorted by increasing eigenvalue)   */
    4203             : /*                                                                 */
    4204             : /*******************************************************************/
    4205             : 
    4206             : GEN
    4207          63 : mateigen(GEN x, long flag, long prec)
    4208             : {
    4209             :   GEN y, R, T;
    4210          63 :   long k, l, ex, n = lg(x);
    4211          63 :   pari_sp av = avma;
    4212             : 
    4213          63 :   if (typ(x)!=t_MAT) pari_err_TYPE("eigen",x);
    4214          63 :   if (n != 1 && n != lgcols(x)) pari_err_DIM("eigen");
    4215          63 :   if (flag < 0 || flag > 1) pari_err_FLAG("mateigen");
    4216          63 :   if (n == 1)
    4217             :   {
    4218          14 :     if (flag) retmkvec2(cgetg(1,t_VEC), cgetg(1,t_MAT));
    4219           7 :     return cgetg(1,t_VEC);
    4220             :   }
    4221          49 :   if (n == 2)
    4222             :   {
    4223          14 :     if (flag) retmkvec2(mkveccopy(gcoeff(x,1,1)), matid(1));
    4224           7 :     return matid(1);
    4225             :   }
    4226             : 
    4227          35 :   ex = 16 - prec2nbits(prec);
    4228          35 :   T = charpoly(x,0);
    4229          35 :   if (RgX_is_QX(T))
    4230             :   {
    4231          28 :     T = ZX_radical( Q_primpart(T) );
    4232          28 :     R = nfrootsQ(T);
    4233          28 :     if (lg(R)-1 < degpol(T))
    4234             :     { /* add missing complex roots */
    4235          14 :       GEN r = cleanroots(RgX_div(T, roots_to_pol(R, 0)), prec);
    4236          14 :       settyp(r, t_VEC);
    4237          14 :       R = shallowconcat(R, r);
    4238             :     }
    4239             :   }
    4240             :   else
    4241             :   {
    4242           7 :     GEN r1, v = vectrunc_init(lg(T));
    4243             :     long e;
    4244           7 :     R = cleanroots(T,prec);
    4245           7 :     r1 = NULL;
    4246          21 :     for (k = 1; k < lg(R); k++)
    4247             :     {
    4248          14 :       GEN r2 = gel(R,k), r = grndtoi(r2, &e);
    4249          14 :       if (e < ex) r2 = r;
    4250          14 :       if (r1)
    4251             :       {
    4252           7 :         r = gsub(r1,r2);
    4253           7 :         if (gequal0(r) || gexpo(r) < ex) continue;
    4254             :       }
    4255          14 :       vectrunc_append(v, r2);
    4256          14 :       r1 = r2;
    4257             :     }
    4258           7 :     R = v;
    4259             :   }
    4260             :   /* R = distinct complex roots of charpoly(x) */
    4261          35 :   l = lg(R); y = cgetg(l, t_VEC);
    4262         189 :   for (k = 1; k < l; k++)
    4263             :   {
    4264         154 :     GEN F = ker_aux(RgM_Rg_sub_shallow(x, gel(R,k)), x);
    4265         154 :     long d = lg(F)-1;
    4266         154 :     if (!d) pari_err_PREC("mateigen");
    4267         154 :     gel(y,k) = F;
    4268         154 :     if (flag) gel(R,k) = const_vec(d, gel(R,k));
    4269             :   }
    4270          35 :   y = shallowconcat1(y);
    4271          35 :   if (lg(y) > n) pari_err_PREC("mateigen");
    4272             :   /* lg(y) < n if x is not diagonalizable */
    4273          35 :   if (flag) y = mkvec2(shallowconcat1(R), y);
    4274          35 :   return gerepilecopy(av,y);
    4275             : }
    4276             : GEN
    4277           0 : eigen(GEN x, long prec) { return mateigen(x, 0, prec); }
    4278             : 
    4279             : /*******************************************************************/
    4280             : /*                                                                 */
    4281             : /*                           DETERMINANT                           */
    4282             : /*                                                                 */
    4283             : /*******************************************************************/
    4284             : 
    4285             : GEN
    4286        1617 : det0(GEN a,long flag)
    4287             : {
    4288        1617 :   switch(flag)
    4289             :   {
    4290        1603 :     case 0: return det(a);
    4291          14 :     case 1: return det2(a);
    4292           0 :     default: pari_err_FLAG("matdet");
    4293             :   }
    4294             :   return NULL; /* LCOV_EXCL_LINE */
    4295             : }
    4296             : 
    4297             : /* M a 2x2 matrix, returns det(M) */
    4298             : static GEN
    4299        2486 : RgM_det2(GEN M)
    4300             : {
    4301        2486 :   pari_sp av = avma;
    4302        2486 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    4303        2486 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    4304        2486 :   return gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
    4305             : }
    4306             : /* M a 2x2 ZM, returns det(M) */
    4307             : static GEN
    4308        7448 : ZM_det2(GEN M)
    4309             : {
    4310        7448 :   pari_sp av = avma;
    4311        7448 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    4312        7448 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    4313        7448 :   return gerepileuptoint(av, subii(mulii(a,d), mulii(b, c)));
    4314             : }
    4315             : /* M a 3x3 ZM, return det(M) */
    4316             : static GEN
    4317        2933 : ZM_det3(GEN M)
    4318             : {
    4319        2933 :   pari_sp av = avma;
    4320        2933 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2), c = gcoeff(M,1,3);
    4321        2933 :   GEN d = gcoeff(M,2,1), e = gcoeff(M,2,2), f = gcoeff(M,2,3);
    4322        2933 :   GEN g = gcoeff(M,3,1), h = gcoeff(M,3,2), i = gcoeff(M,3,3);
    4323        2933 :   GEN t, D = signe(i)? mulii(subii(mulii(a,e), mulii(b,d)), i): gen_0;
    4324        2933 :   if (signe(g))
    4325             :   {
    4326        2730 :     t = mulii(subii(mulii(b,f), mulii(c,e)), g);
    4327        2730 :     D = addii(D, t);
    4328             :   }
    4329        2933 :   if (signe(h))
    4330             :   {
    4331        2751 :     t = mulii(subii(mulii(c,d), mulii(a,f)), h);
    4332        2751 :     D = addii(D, t);
    4333             :   }
    4334        2933 :   return gerepileuptoint(av, D);
    4335             : }
    4336             : 
    4337             : static GEN
    4338        9863 : det_simple_gauss(GEN a, GEN data, pivot_fun pivot)
    4339             : {
    4340        9863 :   pari_sp av = avma;
    4341        9863 :   long i,j,k, s = 1, nbco = lg(a)-1;
    4342        9863 :   GEN p, x = gen_1;
    4343             : 
    4344        9863 :   a = RgM_shallowcopy(a);
    4345       78526 :   for (i=1; i<nbco; i++)
    4346             :   {
    4347       68670 :     k = pivot(a, data, i, NULL);
    4348       68670 :     if (k > nbco) return gerepilecopy(av, gcoeff(a,i,i));
    4349       68663 :     if (k != i)
    4350             :     { /* exchange the lines s.t. k = i */
    4351        4522 :       for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    4352        4522 :       s = -s;
    4353             :     }
    4354       68663 :     p = gcoeff(a,i,i);
    4355             : 
    4356       68663 :     x = gmul(x,p);
    4357      401633 :     for (k=i+1; k<=nbco; k++)
    4358             :     {
    4359      332970 :       GEN m = gcoeff(a,i,k);
    4360      332970 :       if (gequal0(m)) continue;
    4361             : 
    4362      116173 :       m = gdiv(m,p);
    4363      772623 :       for (j=i+1; j<=nbco; j++)
    4364      656450 :         gcoeff(a,j,k) = gsub(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
    4365             :     }
    4366       68663 :     if (gc_needed(av,2))
    4367             :     {
    4368           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
    4369           0 :       gerepileall(av,2, &a,&x);
    4370             :     }
    4371             :   }
    4372        9856 :   if (s < 0) x = gneg_i(x);
    4373        9856 :   return gerepileupto(av, gmul(x, gcoeff(a,nbco,nbco)));
    4374             : }
    4375             : 
    4376             : GEN
    4377        3280 : det2(GEN a)
    4378             : {
    4379             :   GEN data;
    4380             :   pivot_fun pivot;
    4381        3280 :   long n = lg(a)-1;
    4382        3280 :   if (typ(a)!=t_MAT) pari_err_TYPE("det2",a);
    4383        3280 :   if (!n) return gen_1;
    4384        3280 :   if (n != nbrows(a)) pari_err_DIM("det2");
    4385        3280 :   if (n == 1) return gcopy(gcoeff(a,1,1));
    4386        3280 :   if (n == 2) return RgM_det2(a);
    4387        1690 :   pivot = get_pivot_fun(a, a, &data);
    4388        1690 :   return det_simple_gauss(a, data, pivot);
    4389             : }
    4390             : 
    4391             : static GEN
    4392         672 : mydiv(GEN x, GEN y)
    4393             : {
    4394         672 :   long tx = typ(x), ty = typ(y);
    4395         672 :   if (tx == ty && tx == t_POL && varn(x) == varn(y)) return RgX_div(x,y);
    4396         672 :   return gdiv(x,y);
    4397             : }
    4398             : 
    4399             : /* Assumes a a square t_MAT of dimension n > 0. Returns det(a) using
    4400             :  * Gauss-Bareiss. */
    4401             : static GEN
    4402         196 : det_bareiss(GEN a)
    4403             : {
    4404         196 :   pari_sp av = avma;
    4405         196 :   long nbco = lg(a)-1,i,j,k,s = 1;
    4406             :   GEN p, pprec;
    4407             : 
    4408         196 :   a = RgM_shallowcopy(a);
    4409         770 :   for (pprec=gen_1,i=1; i<nbco; i++,pprec=p)
    4410             :   {
    4411             :     GEN ci;
    4412         574 :     int diveuc = (gequal1(pprec)==0);
    4413             : 
    4414         574 :     p = gcoeff(a,i,i);
    4415         574 :     if (gequal0(p))
    4416             :     {
    4417           0 :       k=i+1; while (k<=nbco && gequal0(gcoeff(a,i,k))) k++;
    4418           0 :       if (k>nbco) return gerepilecopy(av, p);
    4419           0 :       swap(gel(a,k), gel(a,i)); s = -s;
    4420           0 :       p = gcoeff(a,i,i);
    4421             :     }
    4422         574 :     ci = gel(a,i);
    4423        1708 :     for (k=i+1; k<=nbco; k++)
    4424             :     {
    4425        1134 :       GEN ck = gel(a,k), m = gel(ck,i);
    4426        1134 :       if (gequal0(m))
    4427             :       {
    4428           0 :         if (gequal1(p))
    4429             :         {
    4430           0 :           if (diveuc)
    4431           0 :             gel(a,k) = mydiv(gel(a,k), pprec);
    4432             :         }
    4433             :         else
    4434           0 :           for (j=i+1; j<=nbco; j++)
    4435             :           {
    4436           0 :             GEN p1 = gmul(p, gel(ck,j));
    4437           0 :             if (diveuc) p1 = mydiv(p1,pprec);
    4438           0 :             gel(ck,j) = p1;
    4439             :           }
    4440             :       }
    4441             :       else
    4442        3752 :         for (j=i+1; j<=nbco; j++)
    4443             :         {
    4444        2618 :           pari_sp av2 = avma;
    4445        2618 :           GEN p1 = gsub(gmul(p,gel(ck,j)), gmul(m,gel(ci,j)));
    4446        2618 :           if (diveuc) p1 = mydiv(p1,pprec);
    4447        2618 :           gel(ck,j) = gerepileupto(av2, p1);
    4448             :         }
    4449        1134 :       if (gc_needed(av,2))
    4450             :       {
    4451           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
    4452           0 :         gerepileall(av,2, &a,&pprec);
    4453           0 :         ci = gel(a,i);
    4454           0 :         p = gcoeff(a,i,i);
    4455             :       }
    4456             :     }
    4457             :   }
    4458         196 :   p = gcoeff(a,nbco,nbco);
    4459         196 :   p = (s < 0)? gneg(p): gcopy(p);
    4460         196 :   return gerepileupto(av, p);
    4461             : }
    4462             : 
    4463             : /* count non-zero entries in col j, at most 'max' of them.
    4464             :  * Return their indices */
    4465             : static GEN
    4466         609 : col_count_non_zero(GEN a, long j, long max)
    4467             : {
    4468         609 :   GEN v = cgetg(max+1, t_VECSMALL);
    4469         609 :   GEN c = gel(a,j);
    4470         609 :   long i, l = lg(a), k = 1;
    4471        2695 :   for (i = 1; i < l; i++)
    4472        2632 :     if (!gequal0(gel(c,i)))
    4473             :     {
    4474        2324 :       if (k > max) return NULL; /* fail */
    4475        1778 :       v[k++] = i;
    4476             :     }
    4477          63 :   setlg(v, k); return v;
    4478             : }
    4479             : /* count non-zero entries in row i, at most 'max' of them.
    4480             :  * Return their indices */
    4481             : static GEN
    4482         560 : row_count_non_zero(GEN a, long i, long max)
    4483             : {
    4484         560 :   GEN v = cgetg(max+1, t_VECSMALL);
    4485         560 :   long j, l = lg(a), k = 1;
    4486        2422 :   for (j = 1; j < l; j++)
    4487        2366 :     if (!gequal0(gcoeff(a,i,j)))
    4488             :     {
    4489        2212 :       if (k > max) return NULL; /* fail */
    4490        1708 :       v[k++] = j;
    4491             :     }
    4492          56 :   setlg(v, k); return v;
    4493             : }
    4494             : 
    4495             : static GEN det_develop(GEN a, long max, double bound);
    4496             : /* (-1)^(i+j) a[i,j] * det RgM_minor(a,i,j) */
    4497             : static GEN
    4498         203 : coeff_det(GEN a, long i, long j, long max, double bound)
    4499             : {
    4500         203 :   GEN c = gcoeff(a, i, j);
    4501         203 :   c = gmul(c, det_develop(RgM_minor(a, i,j), max, bound));
    4502         203 :   if (odd(i+j)) c = gneg(c);
    4503         203 :   return c;
    4504             : }
    4505             : /* a square t_MAT, 'bound' a rough upper bound for the number of
    4506             :  * multiplications we are willing to pay while developing rows/columns before
    4507             :  * switching to Gaussian elimination */
    4508             : static GEN
    4509         301 : det_develop(GEN M, long max, double bound)
    4510             : {
    4511         301 :   pari_sp av = avma;
    4512         301 :   long i,j, n = lg(M)-1, lbest = max+2, best_col = 0, best_row = 0;
    4513         301 :   GEN best = NULL;
    4514             : 
    4515         301 :   if (bound < 1.) return det_bareiss(M); /* too costly now */
    4516             : 
    4517         175 :   switch(n)
    4518             :   {
    4519           0 :     case 0: return gen_1;
    4520           0 :     case 1: return gcopy(gcoeff(M,1,1));
    4521          14 :     case 2: return RgM_det2(M);
    4522             :   }
    4523         161 :   if (max > ((n+2)>>1)) max = (n+2)>>1;
    4524         735 :   for (j = 1; j <= n; j++)
    4525             :   {
    4526         609 :     pari_sp av2 = avma;
    4527         609 :     GEN v = col_count_non_zero(M, j, max);
    4528             :     long lv;
    4529         609 :     if (!v || (lv = lg(v)) >= lbest) { avma = av2; continue; }
    4530          63 :     if (lv == 1) { avma = av; return gen_0; }
    4531          63 :     if (lv == 2) {
    4532          35 :       avma = av;
    4533          35 :       return gerepileupto(av, coeff_det(M,v[1],j,max,bound));
    4534             :     }
    4535          28 :     best = v; lbest = lv; best_col = j;
    4536             :   }
    4537         686 :   for (i = 1; i <= n; i++)
    4538             :   {
    4539         560 :     pari_sp av2 = avma;
    4540         560 :     GEN v = row_count_non_zero(M, i, max);
    4541             :     long lv;
    4542         560 :     if (!v || (lv = lg(v)) >= lbest) { avma = av2; continue; }
    4543          42 :     if (lv == 1) { avma = av; return gen_0; }
    4544          42 :     if (lv == 2) {
    4545           0 :       avma = av;
    4546           0 :       return gerepileupto(av, coeff_det(M,i,v[1],max,bound));
    4547             :     }
    4548          42 :     best = v; lbest = lv; best_row = i;
    4549             :   }
    4550         126 :   if (best_row)
    4551             :   {
    4552          42 :     double d = lbest-1;
    4553          42 :     GEN s = NULL;
    4554             :     long k;
    4555          42 :     bound /= d*d*d;
    4556         168 :     for (k = 1; k < lbest; k++)
    4557             :     {
    4558         126 :       GEN c = coeff_det(M, best_row, best[k], max, bound);
    4559         126 :       s = s? gadd(s, c): c;
    4560             :     }
    4561          42 :     return gerepileupto(av, s);
    4562             :   }
    4563          84 :   if (best_col)
    4564             :   {
    4565          14 :     double d = lbest-1;
    4566          14 :     GEN s = NULL;
    4567             :     long k;
    4568          14 :     bound /= d*d*d;
    4569          56 :     for (k = 1; k < lbest; k++)
    4570             :     {
    4571          42 :       GEN c = coeff_det(M, best[k], best_col, max, bound);
    4572          42 :       s = s? gadd(s, c): c;
    4573             :     }
    4574          14 :     return gerepileupto(av, s);
    4575             :   }
    4576          70 :   return det_bareiss(M);
    4577             : }
    4578             : 
    4579             : /* area of parallelogram bounded by (v1,v2) */
    4580             : static GEN
    4581       51464 : parallelogramarea(GEN v1, GEN v2)
    4582       51464 : { return gsub(gmul(gnorml2(v1), gnorml2(v2)), gsqr(RgV_dotproduct(v1, v2))); }
    4583             : 
    4584             : /* Square of Hadamard bound for det(a), a square matrix.
    4585             :  * Slightly improvement: instead of using the column norms, use the area of
    4586             :  * the parallelogram formed by pairs of consecutive vectors */
    4587             : GEN
    4588       16716 : RgM_Hadamard(GEN a)
    4589             : {
    4590       16716 :   pari_sp av = avma;
    4591       16716 :   long n = lg(a)-1, i;
    4592             :   GEN B;
    4593       16716 :   if (n == 0) return gen_1;
    4594       16716 :   if (n == 1) return gsqr(gcoeff(a,1,1));
    4595       16716 :   a = RgM_gtofp(a, LOWDEFAULTPREC);
    4596       16716 :   B = gen_1;
    4597       68180 :   for (i = 1; i <= n/2; i++)
    4598       51464 :     B = gmul(B, parallelogramarea(gel(a,2*i-1), gel(a,2*i)));
    4599       16716 :   if (odd(n)) B = gmul(B, gnorml2(gel(a, n)));
    4600       16716 :   return gerepileuptoint(av, ceil_safe(B));
    4601             : }
    4602             : 
    4603             : /* assume dim(a) = n > 0 */
    4604             : static GEN
    4605       27916 : ZM_det_i(GEN M, long n)
    4606             : {
    4607       27916 :   const long DIXON_THRESHOLD = 40;
    4608       27916 :   pari_sp av = avma, av2;
    4609             :   long i;
    4610       27916 :   ulong p, compp, Dp = 1;
    4611             :   forprime_t S;
    4612             :   GEN D, h, q, v, comp;
    4613       27916 :   if (n == 1) return icopy(gcoeff(M,1,1));
    4614       27097 :   if (n == 2) return ZM_det2(M);
    4615       19649 :   if (n == 3) return ZM_det3(M);
    4616       16716 :   h = RgM_Hadamard(M);
    4617       16716 :   if (!signe(h)) { avma = av; return gen_0; }
    4618       16716 :   h = sqrti(h); q = gen_1;
    4619       16716 :   init_modular_big(&S);
    4620       16716 :   p = 0; /* -Wall */
    4621       33432 :   while( cmpii(q, h) <= 0 && (p = u_forprime_next(&S)) )
    4622             :   {
    4623       16716 :     av2 = avma; Dp = Flm_det(ZM_to_Flm(M, p), p);
    4624       16716 :     avma = av2;
    4625       16716 :     if (Dp) break;
    4626           0 :     q = muliu(q, p);
    4627             :   }
    4628       16716 :   if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
    4629       16716 :   if (!Dp) { avma = av; return gen_0; }
    4630       16716 :   if (n <= DIXON_THRESHOLD)
    4631       16716 :     D = q;
    4632             :   else
    4633             :   {
    4634           0 :     av2 = avma;
    4635           0 :     v = cgetg(n+1, t_COL);
    4636           0 :     gel(v, 1) = gen_1; /* ensure content(v) = 1 */
    4637           0 :     for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
    4638           0 :     D = Q_denom(ZM_gauss(M, v));
    4639           0 :     if (expi(D) < expi(h) >> 1)
    4640             :     { /* First try unlucky, try once more */
    4641           0 :       for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
    4642           0 :       D = lcmii(D, Q_denom(ZM_gauss(M, v)));
    4643             :     }
    4644           0 :     D = gerepileuptoint(av2, D);
    4645           0 :     if (q != gen_1) D = lcmii(D, q);
    4646             :   }
    4647             :   /* determinant is a multiple of D */
    4648       16716 :   h = shifti(divii(h, D), 1);
    4649             : 
    4650       16716 :   compp = Fl_div(Dp, umodiu(D,p), p);
    4651       16716 :   comp = Z_init_CRT(compp, p);
    4652       16716 :   q = utoipos(p);
    4653       85128 :   while (cmpii(q, h) <= 0)
    4654             :   {
    4655       51696 :     p = u_forprime_next(&S);
    4656       51696 :     if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
    4657       51696 :     Dp = umodiu(D, p);
    4658       51696 :     if (!Dp) continue;
    4659       51696 :     av2 = avma;
    4660       51696 :     compp = Fl_div(Flm_det(ZM_to_Flm(M, p), p), Dp, p);
    4661       51696 :     avma = av2;
    4662       51696 :     (void) Z_incremental_CRT(&comp, compp, &q, p);
    4663             :   }
    4664       16716 :   return gerepileuptoint(av, mulii(comp, D));
    4665             : }
    4666             : 
    4667             : static long
    4668          98 : det_init_max(long n)
    4669             : {
    4670          98 :   if (n > 100) return 0;
    4671          98 :   if (n > 50) return 1;
    4672          98 :   if (n > 30) return 4;
    4673          98 :   return 7;
    4674             : }
    4675             : 
    4676             : GEN
    4677       11125 : det(GEN a)
    4678             : {
    4679       11125 :   long n = lg(a)-1;
    4680             :   double B;
    4681       11125 :   GEN data, ff = NULL, p = NULL;
    4682             :   pivot_fun pivot;
    4683             : 
    4684       11125 :   if (typ(a)!=t_MAT) pari_err_TYPE("det",a);
    4685       11125 :   if (!n) return gen_1;
    4686       11083 :   if (n != nbrows(a)) pari_err_DIM("det");
    4687       11076 :   if (n == 1) return gcopy(gcoeff(a,1,1));
    4688       10882 :   if (n == 2) return RgM_det2(a);
    4689       10000 :   if (RgM_is_FpM(a, &p))
    4690             :   {
    4691        1708 :     pari_sp av = avma;
    4692             :     ulong pp, d;
    4693        1708 :     if (!p) return ZM_det_i(a, n); /* ZM */
    4694             :     /* FpM */
    4695        1477 :     a = RgM_Fp_init(a,p,&pp);
    4696        1477 :     switch(pp)
    4697             :     {
    4698          49 :     case 0: return gerepileupto(av, Fp_to_mod(FpM_det(a,p),p)); break;
    4699          14 :     case 2: d = F2m_det(a); break;
    4700        1414 :     default:d = Flm_det(a,pp); break;
    4701             :     }
    4702        1428 :     avma = av; return mkintmodu(d, pp);
    4703             :   }
    4704        8292 :   if (RgM_is_FFM(a, &ff)) return FFM_det(a, ff);
    4705        8271 :   pivot = get_pivot_fun(a, a, &data);
    4706        8271 :   if (pivot != gauss_get_pivot_NZ) return det_simple_gauss(a, data, pivot);
    4707          98 :   B = (double)n;
    4708          98 :   return det_develop(a, det_init_max(n), B*B*B);
    4709             : }
    4710             : 
    4711             : GEN
    4712       27685 : ZM_det(GEN a)
    4713             : {
    4714       27685 :   long n = lg(a)-1;
    4715       27685 :   if (!n) return gen_1;
    4716       27685 :   return ZM_det_i(a, n);
    4717             : }
    4718             : 
    4719             : /* return a solution of congruence system sum M_{ i,j } X_j = Y_i mod D_i
    4720             :  * If ptu1 != NULL, put in *ptu1 a Z-basis of the homogeneous system */
    4721             : static GEN
    4722          49 : gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1)
    4723             : {
    4724          49 :   pari_sp av = avma;
    4725             :   long n, m, j, l, lM;
    4726             :   GEN delta, H, U, u1, u2, x;
    4727             : 
    4728          49 :   if (typ(M)!=t_MAT) pari_err_TYPE("gaussmodulo",M);
    4729          49 :   lM = lg(M);
    4730          49 :   if (lM == 1)
    4731             :   {
    4732          28 :     switch(typ(Y))
    4733             :     {
    4734           7 :       case t_INT: break;
    4735          14 :       case t_COL: if (lg(Y) != 1) pari_err_DIM("gaussmodulo");
    4736          14 :                   break;
    4737           7 :       default: pari_err_TYPE("gaussmodulo",Y);
    4738             :     }
    4739          21 :     switch(typ(D))
    4740             :     {
    4741           7 :       case t_INT: break;
    4742           7 :       case t_COL: if (lg(D) != 1) pari_err_DIM("gaussmodulo");
    4743           7 :                   break;
    4744           7 :       default: pari_err_TYPE("gaussmodulo",D);
    4745             :     }
    4746          14 :     if (ptu1) *ptu1 = cgetg(1, t_MAT);
    4747          14 :     return gen_0;
    4748             :   }
    4749          21 :   n = nbrows(M);
    4750          21 :   switch(typ(D))
    4751             :   {
    4752             :     case t_COL:
    4753          14 :       if (lg(D)-1!=n) pari_err_DIM("gaussmodulo");
    4754          14 :       delta = diagonal_shallow(D); break;
    4755           7 :     case t_INT: delta = scalarmat_shallow(D,n); break;
    4756           0 :     default: pari_err_TYPE("gaussmodulo",D);
    4757             :       return NULL; /* LCOV_EXCL_LINE */
    4758             :   }
    4759          21 :   switch(typ(Y))
    4760             :   {
    4761           7 :     case t_INT: Y = const_col(n, Y); break;
    4762             :     case t_COL:
    4763          14 :       if (lg(Y)-1!=n) pari_err_DIM("gaussmodulo");
    4764          14 :       break;
    4765           0 :     default: pari_err_TYPE("gaussmodulo",Y);
    4766             :       return NULL; /* LCOV_EXCL_LINE */
    4767             :   }
    4768          21 :   H = ZM_hnfall_i(shallowconcat(M,delta), &U, 1);
    4769          21 :   Y = hnf_solve(H,Y); if (!Y) return gen_0;
    4770          21 :   l = lg(H); /* may be smaller than lM if some moduli are 0 */
    4771          21 :   n = l-1;
    4772          21 :   m = lg(U)-1 - n;
    4773          21 :   u1 = cgetg(m+1,t_MAT);
    4774          21 :   u2 = cgetg(n+1,t_MAT);
    4775          21 :   for (j=1; j<=m; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u1,j) = c; }
    4776          21 :   U += m;
    4777          21 :   for (j=1; j<=n; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u2,j) = c; }
    4778             :   /*       (u1 u2)
    4779             :    * (M D) (*  * ) = (0 H) */
    4780          21 :   u1 = ZM_lll(u1, 0.75, LLL_INPLACE);
    4781          21 :   Y = ZM_ZC_mul(u2,Y);
    4782          21 :   x = ZC_reducemodmatrix(Y, u1);
    4783          21 :   if (!ptu1) x = gerepileupto(av, x);
    4784             :   else
    4785             :   {
    4786          14 :     gerepileall(av, 2, &x, &u1);
    4787          14 :     *ptu1 = u1;
    4788             :   }
    4789          21 :   return x;
    4790             : }
    4791             : 
    4792             : GEN
    4793          49 : matsolvemod0(GEN M, GEN D, GEN Y, long flag)
    4794             : {
    4795             :   pari_sp av;
    4796             :   GEN p1,y;
    4797             : 
    4798          49 :   if (!flag) return gaussmoduloall(M,D,Y,NULL);
    4799             : 
    4800          42 :   av=avma; y = cgetg(3,t_VEC);
    4801          42 :   p1 = gaussmoduloall(M,D,Y, (GEN*)y+2);
    4802          28 :   if (p1==gen_0) { avma=av; return gen_0; }
    4803          14 :   gel(y,1) = p1; return y;
    4804             : }
    4805             : 
    4806             : GEN
    4807           0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod0(M,D,Y,1); }
    4808             : 
    4809             : GEN
    4810           0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod0(M,D,Y,0); }

Generated by: LCOV version 1.11