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-bordeaux1.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 17932-6fdb233) Lines: 2212 2347 94.2 %
Date: 2015-08-03 Functions: 202 217 93.1 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1458 1797 81.1 %

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

Generated by: LCOV version 1.9