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 19825-b77c7f8) Lines: 2249 2431 92.5 %
Date: 2016-12-04 05:49:01 Functions: 199 218 91.3 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11