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 19047-48ce0fc) Lines: 2262 2418 93.5 %
Date: 2016-06-28 Functions: 203 220 92.3 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1506 1869 80.6 %

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

Generated by: LCOV version 1.9