"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/Flv.c" (4 Jan 2021, 32119 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 "Flv.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-2019  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 #include "pari.h"
   16 #include "paripriv.h"
   17 
   18 GEN
   19 Flv_to_ZV(GEN x)
   20 { pari_APPLY_type(t_VEC, utoi(x[i])) }
   21 
   22 GEN
   23 Flc_to_ZC(GEN x)
   24 { pari_APPLY_type(t_COL, utoi(x[i])) }
   25 
   26 GEN
   27 Flm_to_ZM(GEN x)
   28 { pari_APPLY_type(t_MAT, Flc_to_ZC(gel(x,i))) }
   29 
   30 GEN
   31 Flc_to_ZC_inplace(GEN z)
   32 {
   33   long i, l = lg(z);
   34   for (i=1; i<l; i++) gel(z,i) = utoi(z[i]);
   35   settyp(z, t_COL);
   36   return z;
   37 }
   38 
   39 GEN
   40 Flm_to_ZM_inplace(GEN z)
   41 {
   42   long i, l = lg(z);
   43   for (i=1; i<l; i++) Flc_to_ZC_inplace(gel(z, i));
   44   return z;
   45 }
   46 
   47 static GEN
   48 Flm_solve_upper_1(GEN U, GEN B, ulong p, ulong pi)
   49 { return Flm_Fl_mul_pre(B, Fl_inv(ucoeff(U, 1, 1), p), p, pi); }
   50 
   51 static GEN
   52 Flm_rsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
   53 {
   54   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
   55   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
   56   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
   57   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
   58   GEN B1 = rowslice(B, 1, 1);
   59   GEN B2 = rowslice(B, 2, 2);
   60   GEN X2 = Flm_Fl_mul_pre(B2, dinv, p, pi);
   61   GEN X1 = Flm_Fl_mul_pre(Flm_sub(B1, Flm_Fl_mul_pre(X2, b, p, pi), p),
   62                       ainv, p, pi);
   63   return vconcat(X1, X2);
   64 }
   65 
   66 /* solve U*X = B,  U upper triangular and invertible */
   67 static GEN
   68 Flm_rsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
   69 {
   70   long n = lg(U) - 1, n1;
   71   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
   72   pari_sp av = avma;
   73 
   74   if (n == 0) return B;
   75   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
   76   if (n == 2) return Flm_rsolve_upper_2(U, B, p, pi);
   77   n1 = (n + 1)/2;
   78   U2 = vecslice(U, n1 + 1, n);
   79   U22 = rowslice(U2, n1 + 1, n);
   80   B2 = rowslice(B, n1 + 1, n);
   81   X2 = Flm_rsolve_upper_pre(U22, B2, p, pi);
   82   U12 = rowslice(U2, 1, n1);
   83   B1 = rowslice(B, 1, n1);
   84   B1 = Flm_sub(B1, Flm_mul_pre(U12, X2, p, pi), p);
   85   if (gc_needed(av, 1)) gerepileall(av, 2, &B1, &X2);
   86   U11 = matslice(U, 1,n1, 1,n1);
   87   X1 = Flm_rsolve_upper_pre(U11, B1, p, pi);
   88   X = vconcat(X1, X2);
   89   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
   90   return X;
   91 }
   92 
   93 static GEN
   94 Flm_lsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
   95 {
   96   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
   97   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
   98   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
   99   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
  100   GEN B1 = vecslice(B, 1, 1);
  101   GEN B2 = vecslice(B, 2, 2);
  102   GEN X1 = Flm_Fl_mul_pre(B1, ainv, p, pi);
  103   GEN X2 = Flm_Fl_mul_pre(Flm_sub(B2, Flm_Fl_mul_pre(X1, b, p, pi), p),
  104                       dinv, p, pi);
  105   return shallowconcat(X1, X2);
  106 }
  107 
  108 /* solve X*U = B,  U upper triangular and invertible */
  109 static GEN
  110 Flm_lsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
  111 {
  112   long n = lg(U) - 1, n1;
  113   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
  114   pari_sp av = avma;
  115 
  116   if (n == 0) return B;
  117   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
  118   if (n == 2) return Flm_lsolve_upper_2(U, B, p, pi);
  119   n1 = (n + 1)/2;
  120   U2 = vecslice(U, n1 + 1, n);
  121   U11 = matslice(U, 1,n1, 1,n1);
  122   U12 = rowslice(U2, 1, n1);
  123   U22 = rowslice(U2, n1 + 1, n);
  124   B1 = vecslice(B, 1, n1);
  125   B2 = vecslice(B, n1 + 1, n);
  126   X1 = Flm_lsolve_upper_pre(U11, B1, p, pi);
  127   B2 = Flm_sub(B2, Flm_mul_pre(X1, U12, p, pi), p);
  128   if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
  129   X2 = Flm_lsolve_upper_pre(U22, B2, p, pi);
  130   X = shallowconcat(X1, X2);
  131   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
  132   return X;
  133 }
  134 
  135 static GEN
  136 Flm_rsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
  137 {
  138   GEN X1 = rowslice(A, 1, 1);
  139   GEN X2 = Flm_sub(rowslice(A, 2, 2),
  140                    Flm_Fl_mul_pre(X1, ucoeff(L, 2, 1), p, pi), p);
  141   return vconcat(X1, X2);
  142 }
  143 
  144 /* solve L*X = A,  L lower triangular with ones on the diagonal
  145 * (at least as many rows as columns) */
  146 static GEN
  147 Flm_rsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
  148 {
  149   long m = lg(L) - 1, m1, n;
  150   GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
  151   pari_sp av = avma;
  152 
  153   if (m == 0) return zero_Flm(0, lg(A) - 1);
  154   if (m == 1) return rowslice(A, 1, 1);
  155   if (m == 2) return Flm_rsolve_lower_unit_2(L, A, p, pi);
  156   m1 = (m + 1)/2;
  157   n = nbrows(L);
  158   L1 = vecslice(L, 1, m1);
  159   L11 = rowslice(L1, 1, m1);
  160   L21 = rowslice(L1, m1 + 1, n);
  161   A1 = rowslice(A, 1, m1);
  162   X1 = Flm_rsolve_lower_unit_pre(L11, A1, p, pi);
  163   A2 = rowslice(A, m1 + 1, n);
  164   A2 = Flm_sub(A2, Flm_mul_pre(L21, X1, p, pi), p);
  165   if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
  166   L22 = matslice(L, m1+1,n, m1+1,m);
  167   X2 = Flm_rsolve_lower_unit_pre(L22, A2, p, pi);
  168   X = vconcat(X1, X2);
  169   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
  170   return X;
  171 }
  172 
  173 static GEN
  174 Flm_lsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
  175 {
  176   GEN X2 = vecslice(A, 2, 2);
  177   GEN X1 = Flm_sub(vecslice(A, 1, 1),
  178                    Flm_Fl_mul_pre(X2, ucoeff(L, 2, 1), p, pi), p);
  179   return shallowconcat(X1, X2);
  180 }
  181 
  182 /* solve L*X = A,  L square lower triangular with ones on the diagonal */
  183 static GEN
  184 Flm_lsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
  185 {
  186   long m = lg(L) - 1, m1;
  187   GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
  188   pari_sp av = avma;
  189 
  190   if (m <= 1) return A;
  191   if (m == 2) return Flm_lsolve_lower_unit_2(L, A, p, pi);
  192   m1 = (m + 1)/2;
  193   L2 = vecslice(L, m1 + 1, m);
  194   L22 = rowslice(L2, m1 + 1, m);
  195   A2 = vecslice(A, m1 + 1, m);
  196   X2 = Flm_lsolve_lower_unit_pre(L22, A2, p, pi);
  197   if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
  198   L1 = vecslice(L, 1, m1);
  199   L21 = rowslice(L1, m1 + 1, m);
  200   A1 = vecslice(A, 1, m1);
  201   A1 = Flm_sub(A1, Flm_mul_pre(X2, L21, p, pi), p);
  202   L11 = rowslice(L1, 1, m1);
  203   if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
  204   X1 = Flm_lsolve_lower_unit_pre(L11, A1, p, pi);
  205   X = shallowconcat(X1, X2);
  206   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
  207   return X;
  208 }
  209 
  210 /* destroy A */
  211 static long
  212 Flm_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
  213 {
  214   long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
  215   ulong u, v;
  216 
  217   if (P) *P = identity_perm(n);
  218   *R = cgetg(m + 1, t_VECSMALL);
  219   for (j = 1, pr = 0; j <= n; j++)
  220   {
  221     for (pr++, pc = 0; pr <= m; pr++)
  222     {
  223       for (k = j; k <= n; k++) { v = ucoeff(A, pr, k); if (!pc && v) pc = k; }
  224       if (pc) break;
  225     }
  226     if (!pc) break;
  227     (*R)[j] = pr;
  228     if (pc != j)
  229     {
  230       swap(gel(A, j), gel(A, pc));
  231       if (P) lswap((*P)[j], (*P)[pc]);
  232     }
  233     u = Fl_inv(ucoeff(A, pr, j), p);
  234     for (i = pr + 1; i <= m; i++)
  235     {
  236       v = Fl_mul_pre(ucoeff(A, i, j), u, p, pi);
  237       ucoeff(A, i, j) = v;
  238       v = Fl_neg(v, p);
  239       for (k = j + 1; k <= n; k++)
  240         ucoeff(A, i, k) = Fl_addmul_pre(ucoeff(A, i, k),
  241                                         ucoeff(A, pr, k), v, p, pi);
  242     }
  243   }
  244   setlg(*R, j);
  245   *C = vecslice(A, 1, j - 1);
  246   if (U) *U = rowpermute(A, *R);
  247   return j - 1;
  248 }
  249 
  250 static const long Flm_CUP_LIMIT = 8;
  251 
  252 static long
  253 Flm_CUP_pre(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
  254 {
  255   long m = nbrows(A), m1, n = lg(A) - 1, i, r1, r2, r;
  256   GEN R1, C1, U1, P1, R2, C2, U2, P2;
  257   GEN A1, A2, B2, C21, U11, U12, T21, T22;
  258   pari_sp av = avma;
  259 
  260   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
  261     /* destroy A; not called at the outermost recursion level */
  262     return Flm_CUP_basecase(A, R, C, U, P, p, pi);
  263   m1 = (minss(m, n) + 1)/2;
  264   A1 = rowslice(A, 1, m1);
  265   A2 = rowslice(A, m1 + 1, m);
  266   r1 = Flm_CUP_pre(A1, &R1, &C1, &U1, &P1, p, pi);
  267   if (r1 == 0)
  268   {
  269     r2 = Flm_CUP_pre(A2, &R2, &C2, &U2, &P2, p, pi);
  270     *R = cgetg(r2 + 1, t_VECSMALL);
  271     for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
  272     *C = vconcat(zero_Flm(m1, r2), C2);
  273     *U = U2;
  274     *P = P2;
  275     r = r2;
  276   }
  277   else
  278   {
  279     U11 = vecslice(U1, 1, r1);
  280     U12 = vecslice(U1, r1 + 1, n);
  281     T21 = vecslicepermute(A2, P1, 1, r1);
  282     T22 = vecslicepermute(A2, P1, r1 + 1, n);
  283     C21 = Flm_lsolve_upper_pre(U11, T21, p, pi);
  284     if (gc_needed(av, 1))
  285       gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
  286     B2 = Flm_sub(T22, Flm_mul_pre(C21, U12, p, pi), p);
  287     r2 = Flm_CUP_pre(B2, &R2, &C2, &U2, &P2, p, pi);
  288     r = r1 + r2;
  289     *R = cgetg(r + 1, t_VECSMALL);
  290     for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
  291     for (;      i <= r; i++)  (*R)[i] = R2[i - r1] + m1;
  292     *C = shallowconcat(vconcat(C1, C21),
  293                        vconcat(zero_Flm(m1, r2), C2));
  294     *U = shallowconcat(vconcat(U11, zero_Flm(r2, r1)),
  295                        vconcat(vecpermute(U12, P2), U2));
  296     *P = cgetg(n + 1, t_VECSMALL);
  297     for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
  298     for (; i <= n; i++)       (*P)[i] = P1[P2[i - r1] + r1];
  299   }
  300   if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
  301   return r;
  302 }
  303 
  304 static long
  305 Flm_echelon_gauss(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
  306 { return Flm_CUP_basecase(A, R, C, NULL, NULL, p, pi); }
  307 
  308 /* complement of a strictly increasing subsequence of (1, 2, ..., n) */
  309 static GEN
  310 indexcompl(GEN v, long n)
  311 {
  312   long i, j, k, m = lg(v) - 1;
  313   GEN w = cgetg(n - m + 1, t_VECSMALL);
  314   for (i = j = k = 1; i <= n; i++)
  315     if (j <= m && v[j] == i) j++; else w[k++] = i;
  316   return w;
  317 }
  318 
  319 /* column echelon form */
  320 static long
  321 Flm_echelon_pre(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
  322 {
  323   long j, j1, j2, m = nbrows(A), n = lg(A) - 1, n1, r, r1, r2;
  324   GEN A1, A2, R1, R1c, C1, R2, C2;
  325   GEN A12, A22, B2, C11, C21, M12;
  326   pari_sp av = avma;
  327 
  328   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
  329     return Flm_echelon_gauss(Flm_copy(A), R, C, p, pi);
  330 
  331   n1 = (n + 1)/2;
  332   A1 = vecslice(A, 1, n1);
  333   A2 = vecslice(A, n1 + 1, n);
  334   r1 = Flm_echelon_pre(A1, &R1, &C1, p, pi);
  335   if (!r1) return Flm_echelon_pre(A2, R, C, p, pi);
  336   if (r1 == m) { *R = R1; *C = C1; return r1; }
  337 
  338   R1c = indexcompl(R1, m);
  339   C11 = rowpermute(C1, R1);
  340   C21 = rowpermute(C1, R1c);
  341   A12 = rowpermute(A2, R1);
  342   A22 = rowpermute(A2, R1c);
  343   M12 = Flm_rsolve_lower_unit_pre(C11, A12, p, pi);
  344   B2 = Flm_sub(A22, Flm_mul_pre(C21, M12, p, pi), p);
  345   r2 = Flm_echelon_pre(B2, &R2, &C2, p, pi);
  346   if (!r2) { *R = R1; *C = C1; r = r1; }
  347   else
  348   {
  349     R2 = perm_mul(R1c, R2);
  350     C2 = rowpermute(vconcat(zero_Flm(r1, r2), C2),
  351                     perm_inv(vecsmall_concat(R1, R1c)));
  352     r = r1 + r2;
  353     *R = cgetg(r + 1, t_VECSMALL);
  354     *C = cgetg(r + 1, t_MAT);
  355     for (j = j1 = j2 = 1; j <= r; j++)
  356       if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
  357       {
  358         gel(*C, j) = gel(C1, j1);
  359         (*R)[j] = R1[j1++];
  360       }
  361       else
  362       {
  363         gel(*C, j) = gel(C2, j2);
  364         (*R)[j] = R2[j2++];
  365       }
  366   }
  367   if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
  368   return r;
  369 }
  370 
  371 static void /* assume m < p */
  372 _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
  373 {
  374   uel(b,k) = Fl_addmul_pre(uel(b, k), m, uel(b, i), p, pi);
  375 }
  376 static void /* same m = 1 */
  377 _Fl_add(GEN b, long k, long i, ulong p)
  378 {
  379   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
  380 }
  381 static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
  382 _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
  383 {
  384   uel(b,k) += m * uel(b,i);
  385   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
  386 }
  387 static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
  388 _Fl_add_OK(GEN b, long k, long i, ulong p)
  389 {
  390   uel(b,k) += uel(b,i);
  391   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
  392 }
  393 
  394 /* assume 0 <= a[i,j] < p */
  395 static GEN
  396 Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
  397 {
  398   GEN u = cgetg(li+1,t_VECSMALL);
  399   ulong m = uel(b,li) % p;
  400   long i,j;
  401 
  402   uel(u,li) = (m * ucoeff(a,li,li)) % p;
  403   for (i = li-1; i > 0; i--)
  404   {
  405     m = p - uel(b,i)%p;
  406     for (j = i+1; j <= li; j++) {
  407       if (m & HIGHBIT) m %= p;
  408       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
  409     }
  410     m %= p;
  411     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
  412     uel(u,i) = m;
  413   }
  414   return u;
  415 }
  416 static GEN
  417 Fl_get_col(GEN a, GEN b, long li, ulong p)
  418 {
  419   GEN u = cgetg(li+1,t_VECSMALL);
  420   ulong m = uel(b,li) % p;
  421   long i,j;
  422 
  423   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
  424   for (i=li-1; i>0; i--)
  425   {
  426     m = b[i]%p;
  427     for (j = i+1; j <= li; j++)
  428       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
  429     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
  430     uel(u,i) = m;
  431   }
  432   return u;
  433 }
  434 
  435 static GEN
  436 Flm_ker_gauss_OK(GEN x, ulong p, long deplin)
  437 {
  438   GEN y, c, d;
  439   long i, j, k, r, t, m, n;
  440   ulong a;
  441 
  442   n = lg(x)-1;
  443   m=nbrows(x); r=0;
  444 
  445   c = zero_zv(m);
  446   d = cgetg(n+1, t_VECSMALL);
  447   a = 0; /* for gcc -Wall */
  448   for (k=1; k<=n; k++)
  449   {
  450     for (j=1; j<=m; j++)
  451       if (!c[j])
  452       {
  453         a = ucoeff(x,j,k) % p;
  454         if (a) break;
  455       }
  456     if (j > m)
  457     {
  458       if (deplin==1) {
  459         c = cgetg(n+1, t_VECSMALL);
  460         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
  461         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
  462         return c;
  463       }
  464       r++; d[k] = 0;
  465     }
  466     else
  467     {
  468       ulong piv = p - Fl_inv(a, p); /* -1/a */
  469       c[j] = k; d[k] = j;
  470       ucoeff(x,j,k) = p-1;
  471       if (piv != 1)
  472         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
  473       for (t=1; t<=m; t++)
  474       {
  475         if (t == j) continue;
  476 
  477         piv = ( ucoeff(x,t,k) %= p );
  478         if (!piv) continue;
  479         if (piv == 1)
  480           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
  481         else
  482           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
  483       }
  484     }
  485   }
  486   if (deplin==1) return NULL;
  487 
  488   y = cgetg(r+1, t_MAT);
  489   for (j=k=1; j<=r; j++,k++)
  490   {
  491     GEN C = cgetg(n+1, t_VECSMALL);
  492 
  493     gel(y,j) = C; while (d[k]) k++;
  494     for (i=1; i<k; i++)
  495       if (d[i])
  496         uel(C,i) = ucoeff(x,d[i],k) % p;
  497       else
  498         uel(C,i) = 0UL;
  499     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
  500   }
  501   if (deplin == 2) {
  502     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
  503     for (i = j = 1; j <= n; j++)
  504       if (d[j]) pc[i++] = j;
  505     return mkvec2(y, pc);
  506   }
  507   return y;
  508 }
  509 
  510 /* in place, destroy x */
  511 static GEN
  512 Flm_ker_gauss(GEN x, ulong p, long deplin)
  513 {
  514   GEN y, c, d;
  515   long i, j, k, r, t, m, n;
  516   ulong a, pi;
  517   n = lg(x)-1;
  518   if (!n) return cgetg(1,t_MAT);
  519   if (SMALL_ULONG(p)) return Flm_ker_gauss_OK(x, p, deplin);
  520   pi = get_Fl_red(p);
  521 
  522   m=nbrows(x); r=0;
  523 
  524   c = zero_zv(m);
  525   d = cgetg(n+1, t_VECSMALL);
  526   a = 0; /* for gcc -Wall */
  527   for (k=1; k<=n; k++)
  528   {
  529     for (j=1; j<=m; j++)
  530       if (!c[j])
  531       {
  532         a = ucoeff(x,j,k);
  533         if (a) break;
  534       }
  535     if (j > m)
  536     {
  537       if (deplin==1) {
  538         c = cgetg(n+1, t_VECSMALL);
  539         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k);
  540         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
  541         return c;
  542       }
  543       r++; d[k] = 0;
  544     }
  545     else
  546     {
  547       ulong piv = p - Fl_inv(a, p); /* -1/a */
  548       c[j] = k; d[k] = j;
  549       ucoeff(x,j,k) = p-1;
  550       if (piv != 1)
  551         for (i=k+1; i<=n; i++)
  552           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
  553       for (t=1; t<=m; t++)
  554       {
  555         if (t == j) continue;
  556 
  557         piv = ucoeff(x,t,k);
  558         if (!piv) continue;
  559         if (piv == 1)
  560           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
  561         else
  562           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
  563       }
  564     }
  565   }
  566   if (deplin==1) return NULL;
  567 
  568   y = cgetg(r+1, t_MAT);
  569   for (j=k=1; j<=r; j++,k++)
  570   {
  571     GEN C = cgetg(n+1, t_VECSMALL);
  572 
  573     gel(y,j) = C; while (d[k]) k++;
  574     for (i=1; i<k; i++)
  575       if (d[i])
  576         uel(C,i) = ucoeff(x,d[i],k);
  577       else
  578         uel(C,i) = 0UL;
  579     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
  580   }
  581   if (deplin == 2) {
  582     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
  583     for (i = j = 1; j <= n; j++)
  584       if (d[j]) pc[i++] = j;
  585     return mkvec2(y, pc);
  586   }
  587   return y;
  588 }
  589 
  590 GEN
  591 Flm_intersect(GEN x, GEN y, ulong p)
  592 {
  593   pari_sp av = avma;
  594   long j, lx = lg(x);
  595   GEN z;
  596 
  597   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
  598   z = Flm_ker(shallowconcat(x,y), p);
  599   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
  600   return gerepileupto(av, Flm_mul(x,z,p));
  601 }
  602 
  603 static GEN
  604 Flm_ker_echelon(GEN x, ulong p, long pivots) {
  605   pari_sp av = avma;
  606   GEN R, Rc, C, C1, C2, S, K;
  607   long n = lg(x) - 1, r;
  608   ulong pi = get_Fl_red(p);
  609   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
  610   Rc = indexcompl(R, n);
  611   C1 = rowpermute(C, R);
  612   C2 = rowpermute(C, Rc);
  613   S = Flm_lsolve_lower_unit_pre(C1, C2, p, pi);
  614   K = vecpermute(shallowconcat(Flm_neg(S, p), matid_Flm(n - r)),
  615                  perm_inv(vecsmall_concat(R, Rc)));
  616   K = Flm_transpose(K);
  617   if (pivots)
  618     return gerepilecopy(av, mkvec2(K, R));
  619   return gerepileupto(av, K);
  620 }
  621 
  622 static GEN
  623 Flm_deplin_echelon(GEN x, ulong p) {
  624   pari_sp av = avma;
  625   GEN R, Rc, C, C1, C2, s, v;
  626   long i, n = lg(x) - 1, r;
  627   ulong pi = get_Fl_red(p);
  628   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
  629   if (r == n) return gc_NULL(av);
  630   Rc = indexcompl(R, n);
  631   i = Rc[1];
  632   C1 = rowpermute(C, R);
  633   C2 = rowslice(C, i, i);
  634   s = Flm_row(Flm_lsolve_lower_unit_pre(C1, C2, p, pi), 1);
  635   v = vecsmallpermute(vecsmall_concat(Flv_neg(s, p), vecsmall_ei(n - r, 1)),
  636                  perm_inv(vecsmall_concat(R, Rc)));
  637   return gerepileuptoleaf(av, v);
  638 }
  639 
  640 static GEN
  641 Flm_ker_i(GEN x, ulong p, long deplin, long inplace) {
  642   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
  643     switch(deplin) {
  644     case 0: return Flm_ker_echelon(x, p, 0);
  645     case 1: return Flm_deplin_echelon(x, p);
  646     case 2: return Flm_ker_echelon(x, p, 1);
  647     }
  648   return Flm_ker_gauss(inplace? x: Flm_copy(x), p, deplin);
  649 }
  650 
  651 GEN
  652 Flm_ker_sp(GEN x, ulong p, long deplin) {
  653   return Flm_ker_i(x, p, deplin, 1);
  654 }
  655 
  656 GEN
  657 Flm_ker(GEN x, ulong p) {
  658   return Flm_ker_i(x, p, 0, 0);
  659 }
  660 
  661 GEN
  662 Flm_deplin(GEN x, ulong p) {
  663   return Flm_ker_i(x, p, 1, 0);
  664 }
  665 
  666 /* in place, destroy a, SMALL_ULONG(p) is TRUE */
  667 static ulong
  668 Flm_det_gauss_OK(GEN a, long nbco, ulong p)
  669 {
  670   long i,j,k, s = 1;
  671   ulong q, x = 1;
  672 
  673   for (i=1; i<nbco; i++)
  674   {
  675     for(k=i; k<=nbco; k++)
  676     {
  677       ulong c = ucoeff(a,k,i) % p;
  678       ucoeff(a,k,i) = c;
  679       if (c) break;
  680     }
  681     for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
  682     if (k > nbco) return ucoeff(a,i,i);
  683     if (k != i)
  684     { /* exchange the lines s.t. k = i */
  685       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
  686       s = -s;
  687     }
  688     q = ucoeff(a,i,i);
  689 
  690     if (x & HIGHMASK) x %= p;
  691     x *= q;
  692     q = Fl_inv(q,p);
  693     for (k=i+1; k<=nbco; k++)
  694     {
  695       ulong m = ucoeff(a,i,k) % p;
  696       if (!m) continue;
  697 
  698       m = p - ((m*q)%p);
  699       for (j=i+1; j<=nbco; j++)
  700       {
  701         ulong c = ucoeff(a,j,k);
  702         if (c & HIGHMASK) c %= p;
  703         ucoeff(a,j,k) = c  + m*ucoeff(a,j,i);
  704       }
  705     }
  706   }
  707   if (x & HIGHMASK) x %= p;
  708   q = ucoeff(a,nbco,nbco);
  709   if (q & HIGHMASK) q %= p;
  710   x = (x*q) % p;
  711   if (s < 0 && x) x = p - x;
  712   return x;
  713 }
  714 
  715 /* in place, destroy a */
  716 static ulong
  717 Flm_det_gauss(GEN a, ulong p)
  718 {
  719   long i,j,k, s = 1, nbco = lg(a)-1;
  720   ulong pi, q, x = 1;
  721 
  722   if (SMALL_ULONG(p)) return Flm_det_gauss_OK(a, nbco, p);
  723   pi = get_Fl_red(p);
  724   for (i=1; i<nbco; i++)
  725   {
  726     for(k=i; k<=nbco; k++)
  727       if (ucoeff(a,k,i)) break;
  728     if (k > nbco) return ucoeff(a,i,i);
  729     if (k != i)
  730     { /* exchange the lines s.t. k = i */
  731       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
  732       s = -s;
  733     }
  734     q = ucoeff(a,i,i);
  735 
  736     x = Fl_mul_pre(x, q, p, pi);
  737     q = Fl_inv(q,p);
  738     for (k=i+1; k<=nbco; k++)
  739     {
  740       ulong m = ucoeff(a,i,k);
  741       if (!m) continue;
  742 
  743       m = Fl_mul_pre(m, q, p, pi);
  744       for (j=i+1; j<=nbco; j++)
  745         ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
  746     }
  747   }
  748   if (s < 0) x = Fl_neg(x, p);
  749   return Fl_mul(x, ucoeff(a,nbco,nbco), p);
  750 }
  751 
  752 static ulong
  753 Flm_det_CUP(GEN a, ulong p) {
  754   GEN R, C, U, P;
  755   long i, n = lg(a) - 1, r;
  756   ulong d;
  757   ulong pi = get_Fl_red(p);
  758   r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi);
  759   if (r < n)
  760     d = 0;
  761   else {
  762     d = perm_sign(P) == 1? 1: p-1;
  763     for (i = 1; i <= n; i++)
  764       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
  765   }
  766   return d;
  767 }
  768 
  769 static ulong
  770 Flm_det_i(GEN x, ulong p, long inplace) {
  771   pari_sp av = avma;
  772   ulong d;
  773   if (lg(x) - 1 >= Flm_CUP_LIMIT)
  774     d = Flm_det_CUP(x, p);
  775   else
  776     d = Flm_det_gauss(inplace? x: Flm_copy(x), p);
  777   return gc_ulong(av, d);
  778 }
  779 
  780 ulong
  781 Flm_det_sp(GEN x, ulong p) { return Flm_det_i(x, p, 1); }
  782 ulong
  783 Flm_det(GEN x, ulong p) { return Flm_det_i(x, p, 0); }
  784 
  785 /* Destroy x */
  786 static GEN
  787 Flm_gauss_pivot(GEN x, ulong p, long *rr)
  788 {
  789   GEN c,d;
  790   long i,j,k,r,t,n,m;
  791 
  792   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
  793 
  794   m=nbrows(x); r=0;
  795   d=cgetg(n+1,t_VECSMALL);
  796   c = zero_zv(m);
  797   for (k=1; k<=n; k++)
  798   {
  799     for (j=1; j<=m; j++)
  800       if (!c[j])
  801       {
  802         ucoeff(x,j,k) %= p;
  803         if (ucoeff(x,j,k)) break;
  804       }
  805     if (j>m) { r++; d[k]=0; }
  806     else
  807     {
  808       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
  809       c[j]=k; d[k]=j;
  810       for (i=k+1; i<=n; i++)
  811         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
  812       for (t=1; t<=m; t++)
  813         if (!c[t]) /* no pivot on that line yet */
  814         {
  815           piv = ucoeff(x,t,k);
  816           if (piv)
  817           {
  818             ucoeff(x,t,k) = 0;
  819             for (i=k+1; i<=n; i++)
  820               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
  821                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
  822           }
  823         }
  824       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
  825     }
  826   }
  827   *rr = r; return gc_const((pari_sp)d, d);
  828 }
  829 
  830 static GEN
  831 Flm_pivots_CUP(GEN x, ulong p, long *rr)
  832 {
  833   long i, n = lg(x) - 1, r;
  834   GEN R, C, U, P, d = zero_zv(n);
  835   ulong pi = get_Fl_red(p);
  836   r = Flm_CUP_pre(x, &R, &C, &U, &P, p, pi);
  837   for(i = 1; i <= r; i++)
  838     d[P[i]] = R[i];
  839   *rr = n - r; return gc_const((pari_sp)d, d);
  840 }
  841 
  842 GEN
  843 Flm_pivots(GEN x, ulong p, long *rr, long inplace)
  844 {
  845   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
  846     return Flm_pivots_CUP(x, p, rr);
  847   return Flm_gauss_pivot(inplace? x: Flm_copy(x), p, rr);
  848 }
  849 
  850 long
  851 Flm_rank(GEN x, ulong p)
  852 {
  853   pari_sp av = avma;
  854   long r;
  855   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT) {
  856     GEN R, C;
  857     ulong pi = get_Fl_red(p);
  858     return gc_long(av, Flm_echelon_pre(x, &R, &C, p, pi));
  859   }
  860   (void) Flm_pivots(x, p, &r, 0);
  861   return gc_long(av, lg(x)-1 - r);
  862 }
  863 
  864 /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
  865  * reduced mod p */
  866 static GEN
  867 Flm_inv_upper_1_ind(GEN A, long index, ulong p)
  868 {
  869   long n = lg(A)-1, i = index, j;
  870   GEN u = const_vecsmall(n, 0);
  871   u[i] = 1;
  872   if (SMALL_ULONG(p))
  873     for (i--; i>0; i--)
  874     {
  875       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
  876       for (j=i+2; j<=n; j++)
  877       {
  878         if (m & HIGHMASK) m %= p;
  879         m += ucoeff(A,i,j) * uel(u,j);
  880       }
  881       u[i] = Fl_neg(m % p, p);
  882     }
  883   else
  884     for (i--; i>0; i--)
  885     {
  886       ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
  887       for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
  888       u[i] = Fl_neg(m, p);
  889     }
  890   return u;
  891 }
  892 static GEN
  893 Flm_inv_upper_1(GEN A, ulong p)
  894 {
  895   long i, l;
  896   GEN B = cgetg_copy(A, &l);
  897   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
  898   return B;
  899 }
  900 /* destroy a, b */
  901 static GEN
  902 Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
  903 {
  904   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
  905   ulong det = 1;
  906   GEN u;
  907 
  908   li = nbrows(a);
  909   bco = lg(b)-1;
  910   for (i=1; i<=aco; i++)
  911   {
  912     ulong invpiv;
  913     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
  914     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
  915     for (k = i; k <= li; k++)
  916     {
  917       ulong piv = ( ucoeff(a,k,i) %= p );
  918       if (piv)
  919       {
  920         ucoeff(a,k,i) = Fl_inv(piv, p);
  921         if (detp)
  922         {
  923           if (det & HIGHMASK) det %= p;
  924           det *= piv;
  925         }
  926         break;
  927       }
  928     }
  929     /* found a pivot on line k */
  930     if (k > li) return NULL;
  931     if (k != i)
  932     { /* swap lines so that k = i */
  933       s = -s;
  934       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
  935       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
  936     }
  937     if (i == aco) break;
  938 
  939     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
  940     for (k=i+1; k<=li; k++)
  941     {
  942       ulong m = ( ucoeff(a,k,i) %= p );
  943       if (!m) continue;
  944 
  945       m = Fl_mul(m, invpiv, p);
  946       if (m == 1) {
  947         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
  948         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
  949       } else {
  950         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
  951         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
  952       }
  953     }
  954   }
  955   if (detp)
  956   {
  957     det %=  p;
  958     if (s < 0 && det) det = p - det;
  959     *detp = det;
  960   }
  961   u = cgetg(bco+1,t_MAT);
  962   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
  963   return u;
  964 }
  965 
  966 /* destroy a, b */
  967 static GEN
  968 Flm_gauss_sp_i(GEN a, GEN b, ulong *detp, ulong p)
  969 {
  970   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
  971   ulong det = 1;
  972   GEN u;
  973   ulong pi;
  974   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
  975   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
  976   pi = get_Fl_red(p);
  977   li = nbrows(a);
  978   bco = lg(b)-1;
  979   for (i=1; i<=aco; i++)
  980   {
  981     ulong invpiv;
  982     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
  983     for (k = i; k <= li; k++)
  984     {
  985       ulong piv = ucoeff(a,k,i);
  986       if (piv)
  987       {
  988         ucoeff(a,k,i) = Fl_inv(piv, p);
  989         if (detp) det = Fl_mul_pre(det, piv, p, pi);
  990         break;
  991       }
  992     }
  993     /* found a pivot on line k */
  994     if (k > li) return NULL;
  995     if (k != i)
  996     { /* swap lines so that k = i */
  997       s = -s;
  998       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
  999       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
 1000     }
 1001     if (i == aco) break;
 1002 
 1003     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
 1004     for (k=i+1; k<=li; k++)
 1005     {
 1006       ulong m = ucoeff(a,k,i);
 1007       if (!m) continue;
 1008 
 1009       m = Fl_mul_pre(m, invpiv, p, pi);
 1010       if (m == 1) {
 1011         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
 1012         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
 1013       } else {
 1014         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
 1015         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
 1016       }
 1017     }
 1018   }
 1019   if (detp)
 1020   {
 1021     if (s < 0 && det) det = p - det;
 1022     *detp = det;
 1023   }
 1024   u = cgetg(bco+1,t_MAT);
 1025   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
 1026   return u;
 1027 }
 1028 
 1029 static GEN
 1030 Flm_gauss_from_CUP(GEN b, GEN R, GEN C, GEN U, GEN P, ulong p, ulong pi, ulong *detp)
 1031 {
 1032   GEN Y = Flm_rsolve_lower_unit_pre(rowpermute(C, R), rowpermute(b, R), p, pi);
 1033   GEN X = rowpermute(Flm_rsolve_upper_pre(U, Y, p, pi), perm_inv(P));
 1034   if (detp)
 1035   {
 1036     ulong d = perm_sign(P) == 1? 1: p-1;
 1037     long i, r = lg(R);
 1038     for (i = 1; i < r; i++)
 1039       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
 1040     *detp = d;
 1041   }
 1042   return X;
 1043 }
 1044 
 1045 static GEN
 1046 Flm_gauss_CUP(GEN a, GEN b, ulong *detp, ulong p) {
 1047   GEN R, C, U, P;
 1048   long n = lg(a) - 1, r;
 1049   ulong pi = get_Fl_red(p);
 1050   if (nbrows(a) < n || (r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi)) < n)
 1051     return NULL;
 1052   return Flm_gauss_from_CUP(b, R, C, U, P, p, pi, detp);
 1053 }
 1054 
 1055 GEN
 1056 Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p) {
 1057   pari_sp av = avma;
 1058   GEN x;
 1059   if (lg(a) - 1 >= Flm_CUP_LIMIT)
 1060     x = Flm_gauss_CUP(a, b, detp, p);
 1061   else
 1062     x = Flm_gauss_sp_i(a, b, detp, p);
 1063   if (!x) return gc_NULL(av);
 1064   return gerepileupto(av, x);
 1065 }
 1066 
 1067 GEN
 1068 Flm_gauss(GEN a, GEN b, ulong p) {
 1069   pari_sp av = avma;
 1070   GEN x;
 1071   if (lg(a) - 1 >= Flm_CUP_LIMIT)
 1072     x = Flm_gauss_CUP(a, b, NULL, p);
 1073   else {
 1074     a = RgM_shallowcopy(a);
 1075     b = RgM_shallowcopy(b);
 1076     x = Flm_gauss_sp(a, b, NULL, p);
 1077   }
 1078   if (!x) return gc_NULL(av);
 1079   return gerepileupto(av, x);
 1080 }
 1081 
 1082 static GEN
 1083 Flm_inv_i(GEN a, ulong *detp, ulong p, long inplace) {
 1084   pari_sp av = avma;
 1085   long n = lg(a) - 1;
 1086   GEN b, x;
 1087   if (n == 0) return cgetg(1, t_MAT);
 1088   b = matid_Flm(nbrows(a));
 1089   if (n >= Flm_CUP_LIMIT)
 1090     x = Flm_gauss_CUP(a, b, detp, p);
 1091   else {
 1092     if (!inplace)
 1093       a = RgM_shallowcopy(a);
 1094     x = Flm_gauss_sp(a, b, detp, p);
 1095   }
 1096   if (!x) return gc_NULL(av);
 1097   return gerepileupto(av, x);
 1098 }
 1099 
 1100 GEN
 1101 Flm_inv_sp(GEN a, ulong *detp, ulong p) {
 1102   return Flm_inv_i(a, detp, p, 1);
 1103 }
 1104 
 1105 GEN
 1106 Flm_inv(GEN a, ulong p) {
 1107   return Flm_inv_i(a, NULL, p, 0);
 1108 }
 1109 
 1110 GEN
 1111 Flm_Flc_gauss(GEN a, GEN b, ulong p) {
 1112   pari_sp av = avma;
 1113   GEN z = Flm_gauss(a, mkmat(b), p);
 1114   if (!z) return gc_NULL(av);
 1115   if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
 1116   return gerepileuptoleaf(av, gel(z,1));
 1117 }
 1118 
 1119 GEN
 1120 Flm_adjoint(GEN A, ulong p)
 1121 {
 1122   pari_sp av = avma;
 1123   GEN R, C, U, P, C1, U1, v, c, d;
 1124   long r, i, q, n = lg(A)-1, m;
 1125   ulong D;
 1126   ulong pi = get_Fl_red(p);
 1127   if (n == 0) return cgetg(1, t_MAT);
 1128   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
 1129   m = nbrows(A);
 1130   if (r == n)
 1131   {
 1132     GEN X = Flm_gauss_from_CUP(matid_Flm(m), R, C, U, P, p, pi, &D);
 1133     return gerepileupto(av, Flm_Fl_mul_pre(X, D, p, pi));
 1134   }
 1135   if (r < n-1) return zero_Flm(n, m);
 1136   for (q = n, i = 1; i < n; i++)
 1137     if (R[i] != i) { q = i; break; }
 1138   C1 = matslice(C, 1, q-1, 1, q-1);
 1139   c = vecslice(Flm_row(C, q), 1, q-1);
 1140   c = Flm_lsolve_lower_unit_pre(C1, Flm_transpose(mkmat(c)), p, pi);
 1141   d = cgetg(m+1, t_VECSMALL);
 1142   for (i=1; i<q; i++)    uel(d,i) = ucoeff(c,1,i);
 1143   uel(d,q) = p-1;
 1144   for (i=q+1; i<=m; i++) uel(d,i) = 0;
 1145   U1 = vecslice(U, 1, n-1);
 1146   v = gel(Flm_rsolve_upper_pre(U1, mkmat(gel(U,n)), p, pi),1);
 1147   v = vecsmall_append(v, p-1);
 1148   D = perm_sign(P) != (odd(q+n)?-1:1) ? p-1 : 1;
 1149   for (i = 1; i <= n-1; i++)
 1150     D = Fl_mul_pre(D, ucoeff(U1, i, i), p, pi);
 1151   d = Flv_Fl_mul(d, D, p);
 1152   return rowpermute(Flc_Flv_mul(v, d, p),perm_inv(P));
 1153 }
 1154 
 1155 static GEN
 1156 Flm_invimage_CUP(GEN A, GEN B, ulong p) {
 1157   pari_sp av = avma;
 1158   GEN R, Rc, C, U, P, B1, B2, C1, C2, X, Y, Z;
 1159   long r;
 1160   ulong pi = get_Fl_red(p);
 1161   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
 1162   Rc = indexcompl(R, nbrows(B));
 1163   C1 = rowpermute(C, R);
 1164   C2 = rowpermute(C, Rc);
 1165   B1 = rowpermute(B, R);
 1166   B2 = rowpermute(B, Rc);
 1167   Z = Flm_rsolve_lower_unit_pre(C1, B1, p, pi);
 1168   if (!gequal(Flm_mul_pre(C2, Z, p, pi), B2))
 1169     return NULL;
 1170   Y = vconcat(Flm_rsolve_upper_pre(vecslice(U, 1, r), Z, p, pi),
 1171               zero_Flm(lg(A) - 1 - r, lg(B) - 1));
 1172   X = rowpermute(Y, perm_inv(P));
 1173   return gerepileupto(av, X);
 1174 }
 1175 
 1176 GEN
 1177 Flm_invimage_i(GEN A, GEN B, ulong p)
 1178 {
 1179   GEN d, x, X, Y;
 1180   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
 1181 
 1182   if (!nB) return cgetg(1, t_MAT);
 1183   if (maxss(nA, nB) >= Flm_CUP_LIMIT && nbrows(B) >= Flm_CUP_LIMIT)
 1184     return Flm_invimage_CUP(A, B, p);
 1185 
 1186   x = Flm_ker(shallowconcat(Flm_neg(A,p), B), p);
 1187   /* AX = BY, Y in strict upper echelon form with pivots = 1.
 1188    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
 1189    * Y has at least nB columns and full rank */
 1190   nY = lg(x)-1;
 1191   if (nY < nB) return NULL;
 1192   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
 1193   d = cgetg(nB+1, t_VECSMALL);
 1194   for (i = nB, j = nY; i >= 1; i--, j--)
 1195   {
 1196     for (; j>=1; j--)
 1197       if (coeff(Y,i,j)) { d[i] = j; break; }
 1198     if (!j) return NULL;
 1199   }
 1200   /* reduce to the case Y square, upper triangular with 1s on diagonal */
 1201   Y = vecpermute(Y, d);
 1202   x = vecpermute(x, d);
 1203   X = rowslice(x, 1, nA);
 1204   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
 1205 }
 1206 GEN
 1207 Flm_invimage(GEN A, GEN B, ulong p)
 1208 {
 1209   pari_sp av = avma;
 1210   GEN X = Flm_invimage_i(A,B,p);
 1211   if (!X) return gc_NULL(av);
 1212   return gerepileupto(av, X);
 1213 }
 1214 
 1215 GEN
 1216 Flm_Flc_invimage(GEN A, GEN y, ulong p)
 1217 {
 1218   pari_sp av = avma;
 1219   long i, l = lg(A);
 1220   GEN M, x;
 1221   ulong t;
 1222 
 1223   if (l==1) return NULL;
 1224   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
 1225   M = cgetg(l+1,t_MAT);
 1226   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
 1227   gel(M,l) = y; M = Flm_ker(M,p);
 1228   i = lg(M)-1; if (!i) return gc_NULL(av);
 1229 
 1230   x = gel(M,i); t = x[l];
 1231   if (!t) return gc_NULL(av);
 1232 
 1233   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
 1234   if (t!=1) x = Flv_Fl_mul(x, t, p);
 1235   return gerepileuptoleaf(av, x);
 1236 }