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 18579-f9e35ad) Lines: 2266 2421 93.6 %
Date: 2016-02-06 Functions: 203 220 92.3 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1504 1871 80.4 %

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

Generated by: LCOV version 1.9