"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/alglin1.c" (12 Dec 2020, 140426 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "alglin1.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

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