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 16937-4bd9b4e) Lines: 2186 2336 93.6 %
Date: 2014-10-24 Functions: 201 216 93.1 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1430 1787 80.0 %

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

Generated by: LCOV version 1.9