Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - alglin1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19383-e934d0b) Lines: 2244 2420 92.7 %
Date: 2016-09-01 06:11:39 Functions: 201 220 91.4 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11