Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - Flv.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.13.0 lcov report (development 25825-edc109b529) Lines: 730 759 96.2 %
Date: 2020-09-21 06:08:33 Functions: 57 61 93.4 %
Legend: Lines: hit not hit

          Line data    Source code
       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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : GEN
      18      257067 : Flv_to_ZV(GEN x)
      19      594724 : { pari_APPLY_type(t_VEC, utoi(x[i])) }
      20             : 
      21             : GEN
      22    23210952 : Flc_to_ZC(GEN x)
      23   194321752 : { pari_APPLY_type(t_COL, utoi(x[i])) }
      24             : 
      25             : GEN
      26     8744509 : Flm_to_ZM(GEN x)
      27    31565990 : { pari_APPLY_type(t_MAT, Flc_to_ZC(gel(x,i))) }
      28             : 
      29             : GEN
      30      123378 : Flc_to_ZC_inplace(GEN z)
      31             : {
      32      123378 :   long i, l = lg(z);
      33      768413 :   for (i=1; i<l; i++) gel(z,i) = utoi(z[i]);
      34      123378 :   settyp(z, t_COL);
      35      123378 :   return z;
      36             : }
      37             : 
      38             : GEN
      39       56166 : Flm_to_ZM_inplace(GEN z)
      40             : {
      41       56166 :   long i, l = lg(z);
      42      179530 :   for (i=1; i<l; i++) Flc_to_ZC_inplace(gel(z, i));
      43       56166 :   return z;
      44             : }
      45             : 
      46             : static GEN
      47      845382 : Flm_solve_upper_1(GEN U, GEN B, ulong p, ulong pi)
      48      845382 : { return Flm_Fl_mul_pre(B, Fl_inv(ucoeff(U, 1, 1), p), p, pi); }
      49             : 
      50             : static GEN
      51     1080463 : Flm_rsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      52             : {
      53     1080463 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      54     1080463 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      55     1080492 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      56     1080493 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
      57     1080488 :   GEN B1 = rowslice(B, 1, 1);
      58     1080437 :   GEN B2 = rowslice(B, 2, 2);
      59     1080440 :   GEN X2 = Flm_Fl_mul_pre(B2, dinv, p, pi);
      60     1080475 :   GEN X1 = Flm_Fl_mul_pre(Flm_sub(B1, Flm_Fl_mul_pre(X2, b, p, pi), p),
      61             :                       ainv, p, pi);
      62     1080463 :   return vconcat(X1, X2);
      63             : }
      64             : 
      65             : /* solve U*X = B,  U upper triangular and invertible */
      66             : static GEN
      67     2712728 : Flm_rsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
      68             : {
      69     2712728 :   long n = lg(U) - 1, n1;
      70             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
      71     2712728 :   pari_sp av = avma;
      72             : 
      73     2712728 :   if (n == 0) return B;
      74     2712721 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
      75     2324823 :   if (n == 2) return Flm_rsolve_upper_2(U, B, p, pi);
      76     1244364 :   n1 = (n + 1)/2;
      77     1244364 :   U2 = vecslice(U, n1 + 1, n);
      78     1244609 :   U22 = rowslice(U2, n1 + 1, n);
      79     1244572 :   B2 = rowslice(B, n1 + 1, n);
      80     1244556 :   X2 = Flm_rsolve_upper_pre(U22, B2, p, pi);
      81     1244646 :   U12 = rowslice(U2, 1, n1);
      82     1244572 :   B1 = rowslice(B, 1, n1);
      83     1244576 :   B1 = Flm_sub(B1, Flm_mul_pre(U12, X2, p, pi), p);
      84     1244553 :   if (gc_needed(av, 1)) gerepileall(av, 2, &B1, &X2);
      85     1244553 :   U11 = matslice(U, 1,n1, 1,n1);
      86     1244549 :   X1 = Flm_rsolve_upper_pre(U11, B1, p, pi);
      87     1244657 :   X = vconcat(X1, X2);
      88     1244673 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
      89     1244673 :   return X;
      90             : }
      91             : 
      92             : static GEN
      93     1292048 : Flm_lsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      94             : {
      95     1292048 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      96     1292048 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      97     1292087 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      98     1292072 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
      99     1292054 :   GEN B1 = vecslice(B, 1, 1);
     100     1292039 :   GEN B2 = vecslice(B, 2, 2);
     101     1291985 :   GEN X1 = Flm_Fl_mul_pre(B1, ainv, p, pi);
     102     1292023 :   GEN X2 = Flm_Fl_mul_pre(Flm_sub(B2, Flm_Fl_mul_pre(X1, b, p, pi), p),
     103             :                       dinv, p, pi);
     104     1292008 :   return shallowconcat(X1, X2);
     105             : }
     106             : 
     107             : /* solve X*U = B,  U upper triangular and invertible */
     108             : static GEN
     109     3020673 : Flm_lsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
     110             : {
     111     3020673 :   long n = lg(U) - 1, n1;
     112             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
     113     3020673 :   pari_sp av = avma;
     114             : 
     115     3020673 :   if (n == 0) return B;
     116     3020673 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
     117     2563180 :   if (n == 2) return Flm_lsolve_upper_2(U, B, p, pi);
     118     1271126 :   n1 = (n + 1)/2;
     119     1271126 :   U2 = vecslice(U, n1 + 1, n);
     120     1271371 :   U11 = matslice(U, 1,n1, 1,n1);
     121     1271340 :   U12 = rowslice(U2, 1, n1);
     122     1271346 :   U22 = rowslice(U2, n1 + 1, n);
     123     1271340 :   B1 = vecslice(B, 1, n1);
     124     1271332 :   B2 = vecslice(B, n1 + 1, n);
     125     1271313 :   X1 = Flm_lsolve_upper_pre(U11, B1, p, pi);
     126     1271326 :   B2 = Flm_sub(B2, Flm_mul_pre(X1, U12, p, pi), p);
     127     1271320 :   if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
     128     1271320 :   X2 = Flm_lsolve_upper_pre(U22, B2, p, pi);
     129     1271345 :   X = shallowconcat(X1, X2);
     130     1271346 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     131     1271347 :   return X;
     132             : }
     133             : 
     134             : static GEN
     135     2231527 : Flm_rsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     136             : {
     137     2231527 :   GEN X1 = rowslice(A, 1, 1);
     138     2231540 :   GEN X2 = Flm_sub(rowslice(A, 2, 2),
     139     2231513 :                    Flm_Fl_mul_pre(X1, ucoeff(L, 2, 1), p, pi), p);
     140     2231478 :   return vconcat(X1, X2);
     141             : }
     142             : 
     143             : /* solve L*X = A,  L lower triangular with ones on the diagonal
     144             : * (at least as many rows as columns) */
     145             : static GEN
     146     5242109 : Flm_rsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     147             : {
     148     5242109 :   long m = lg(L) - 1, m1, n;
     149             :   GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
     150     5242109 :   pari_sp av = avma;
     151             : 
     152     5242109 :   if (m == 0) return zero_Flm(0, lg(A) - 1);
     153     5242109 :   if (m == 1) return rowslice(A, 1, 1);
     154     4541530 :   if (m == 2) return Flm_rsolve_lower_unit_2(L, A, p, pi);
     155     2310005 :   m1 = (m + 1)/2;
     156     2310005 :   n = nbrows(L);
     157     2310213 :   L1 = vecslice(L, 1, m1);
     158     2310218 :   L11 = rowslice(L1, 1, m1);
     159     2310196 :   L21 = rowslice(L1, m1 + 1, n);
     160     2310188 :   A1 = rowslice(A, 1, m1);
     161     2310172 :   X1 = Flm_rsolve_lower_unit_pre(L11, A1, p, pi);
     162     2310249 :   A2 = rowslice(A, m1 + 1, n);
     163     2310166 :   A2 = Flm_sub(A2, Flm_mul_pre(L21, X1, p, pi), p);
     164     2310124 :   if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
     165     2310124 :   L22 = matslice(L, m1+1,n, m1+1,m);
     166     2310168 :   X2 = Flm_rsolve_lower_unit_pre(L22, A2, p, pi);
     167     2310208 :   X = vconcat(X1, X2);
     168     2310261 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     169     2310261 :   return X;
     170             : }
     171             : 
     172             : static GEN
     173      984296 : Flm_lsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     174             : {
     175      984296 :   GEN X2 = vecslice(A, 2, 2);
     176      984296 :   GEN X1 = Flm_sub(vecslice(A, 1, 1),
     177      984296 :                    Flm_Fl_mul_pre(X2, ucoeff(L, 2, 1), p, pi), p);
     178      984296 :   return shallowconcat(X1, X2);
     179             : }
     180             : 
     181             : /* solve L*X = A,  L square lower triangular with ones on the diagonal */
     182             : static GEN
     183     2431833 : Flm_lsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     184             : {
     185     2431833 :   long m = lg(L) - 1, m1;
     186             :   GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
     187     2431833 :   pari_sp av = avma;
     188             : 
     189     2431833 :   if (m <= 1) return A;
     190     2113244 :   if (m == 2) return Flm_lsolve_lower_unit_2(L, A, p, pi);
     191     1128948 :   m1 = (m + 1)/2;
     192     1128948 :   L2 = vecslice(L, m1 + 1, m);
     193     1128948 :   L22 = rowslice(L2, m1 + 1, m);
     194     1128948 :   A2 = vecslice(A, m1 + 1, m);
     195     1128948 :   X2 = Flm_lsolve_lower_unit_pre(L22, A2, p, pi);
     196     1128948 :   if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
     197     1128948 :   L1 = vecslice(L, 1, m1);
     198     1128948 :   L21 = rowslice(L1, m1 + 1, m);
     199     1128948 :   A1 = vecslice(A, 1, m1);
     200     1128948 :   A1 = Flm_sub(A1, Flm_mul_pre(X2, L21, p, pi), p);
     201     1128948 :   L11 = rowslice(L1, 1, m1);
     202     1128948 :   if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
     203     1128948 :   X1 = Flm_lsolve_lower_unit_pre(L11, A1, p, pi);
     204     1128948 :   X = shallowconcat(X1, X2);
     205     1128948 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     206     1128948 :   return X;
     207             : }
     208             : 
     209             : /* destroy A */
     210             : static long
     211     1459208 : Flm_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     212             : {
     213     1459208 :   long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
     214             :   ulong u, v;
     215             : 
     216     1459210 :   if (P) *P = identity_perm(n);
     217     1459204 :   *R = cgetg(m + 1, t_VECSMALL);
     218     7556298 :   for (j = 1, pr = 0; j <= n; j++)
     219             :   {
     220     9766055 :     for (pr++, pc = 0; pr <= m; pr++)
     221             :     {
     222    75406774 :       for (k = j; k <= n; k++) { v = ucoeff(A, pr, k); if (!pc && v) pc = k; }
     223     9035555 :       if (pc) break;
     224             :     }
     225     6827746 :     if (!pc) break;
     226     6097216 :     (*R)[j] = pr;
     227     6097216 :     if (pc != j)
     228             :     {
     229     1028092 :       swap(gel(A, j), gel(A, pc));
     230     1028092 :       if (P) lswap((*P)[j], (*P)[pc]);
     231             :     }
     232     6097216 :     u = Fl_inv(ucoeff(A, pr, j), p);
     233    46118935 :     for (i = pr + 1; i <= m; i++)
     234             :     {
     235    40021830 :       v = Fl_mul_pre(ucoeff(A, i, j), u, p, pi);
     236    40021475 :       ucoeff(A, i, j) = v;
     237    40021475 :       v = Fl_neg(v, p);
     238   204285534 :       for (k = j + 1; k <= n; k++)
     239   164266658 :         ucoeff(A, i, k) = Fl_addmul_pre(ucoeff(A, i, k),
     240   164263680 :                                         ucoeff(A, pr, k), v, p, pi);
     241             :     }
     242             :   }
     243     1459082 :   setlg(*R, j);
     244     1459208 :   *C = vecslice(A, 1, j - 1);
     245     1459180 :   if (U) *U = rowpermute(A, *R);
     246     1459189 :   return j - 1;
     247             : }
     248             : 
     249             : static const long Flm_CUP_LIMIT = 8;
     250             : 
     251             : static long
     252     1308163 : Flm_CUP_pre(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     253             : {
     254     1308163 :   long m = nbrows(A), m1, n = lg(A) - 1, i, r1, r2, r;
     255             :   GEN R1, C1, U1, P1, R2, C2, U2, P2;
     256             :   GEN A1, A2, B2, C21, U11, U12, T21, T22;
     257     1308161 :   pari_sp av = avma;
     258             : 
     259     1308161 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     260             :     /* destroy A; not called at the outermost recursion level */
     261      820875 :     return Flm_CUP_basecase(A, R, C, U, P, p, pi);
     262      487286 :   m1 = (minss(m, n) + 1)/2;
     263      487284 :   A1 = rowslice(A, 1, m1);
     264      487278 :   A2 = rowslice(A, m1 + 1, m);
     265      487281 :   r1 = Flm_CUP_pre(A1, &R1, &C1, &U1, &P1, p, pi);
     266      487283 :   if (r1 == 0)
     267             :   {
     268        9072 :     r2 = Flm_CUP_pre(A2, &R2, &C2, &U2, &P2, p, pi);
     269        9072 :     *R = cgetg(r2 + 1, t_VECSMALL);
     270       14763 :     for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
     271        9072 :     *C = vconcat(zero_Flm(m1, r2), C2);
     272        9072 :     *U = U2;
     273        9072 :     *P = P2;
     274        9072 :     r = r2;
     275             :   }
     276             :   else
     277             :   {
     278      478211 :     U11 = vecslice(U1, 1, r1);
     279      478212 :     U12 = vecslice(U1, r1 + 1, n);
     280      478212 :     T21 = vecslicepermute(A2, P1, 1, r1);
     281      478214 :     T22 = vecslicepermute(A2, P1, r1 + 1, n);
     282      478211 :     C21 = Flm_lsolve_upper_pre(U11, T21, p, pi);
     283      478208 :     if (gc_needed(av, 1))
     284           0 :       gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
     285      478208 :     B2 = Flm_sub(T22, Flm_mul_pre(C21, U12, p, pi), p);
     286      478207 :     r2 = Flm_CUP_pre(B2, &R2, &C2, &U2, &P2, p, pi);
     287      478207 :     r = r1 + r2;
     288      478207 :     *R = cgetg(r + 1, t_VECSMALL);
     289     3519931 :     for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
     290     3257415 :     for (;      i <= r; i++)  (*R)[i] = R2[i - r1] + m1;
     291      478206 :     *C = shallowconcat(vconcat(C1, C21),
     292             :                        vconcat(zero_Flm(m1, r2), C2));
     293      478212 :     *U = shallowconcat(vconcat(U11, zero_Flm(r2, r1)),
     294             :                        vconcat(vecpermute(U12, P2), U2));
     295      478217 :     *P = cgetg(n + 1, t_VECSMALL);
     296     3519981 :     for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
     297     7460148 :     for (; i <= n; i++)       (*P)[i] = P1[P2[i - r1] + r1];
     298             :   }
     299      487285 :   if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
     300      487286 :   return r;
     301             : }
     302             : 
     303             : static long
     304      638299 : Flm_echelon_gauss(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     305      638299 : { return Flm_CUP_basecase(A, R, C, NULL, NULL, p, pi); }
     306             : 
     307             : /* complement of a strictly increasing subsequence of (1, 2, ..., n) */
     308             : static GEN
     309      572085 : indexcompl(GEN v, long n)
     310             : {
     311      572085 :   long i, j, k, m = lg(v) - 1;
     312      572085 :   GEN w = cgetg(n - m + 1, t_VECSMALL);
     313    12109926 :   for (i = j = k = 1; i <= n; i++)
     314    11537841 :     if (j <= m && v[j] == i) j++; else w[k++] = i;
     315      572085 :   return w;
     316             : }
     317             : 
     318             : /* column echelon form */
     319             : static long
     320     1099801 : Flm_echelon_pre(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     321             : {
     322     1099801 :   long j, j1, j2, m = nbrows(A), n = lg(A) - 1, n1, r, r1, r2;
     323             :   GEN A1, A2, R1, R1c, C1, R2, C2;
     324             :   GEN A12, A22, B2, C11, C21, M12;
     325     1099801 :   pari_sp av = avma;
     326             : 
     327     1099801 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     328      638299 :     return Flm_echelon_gauss(Flm_copy(A), R, C, p, pi);
     329             : 
     330      461502 :   n1 = (n + 1)/2;
     331      461502 :   A1 = vecslice(A, 1, n1);
     332      461502 :   A2 = vecslice(A, n1 + 1, n);
     333      461502 :   r1 = Flm_echelon_pre(A1, &R1, &C1, p, pi);
     334      461502 :   if (!r1) return Flm_echelon_pre(A2, R, C, p, pi);
     335      402981 :   if (r1 == m) { *R = R1; *C = C1; return r1; }
     336             : 
     337      398141 :   R1c = indexcompl(R1, m);
     338      398141 :   C11 = rowpermute(C1, R1);
     339      398141 :   C21 = rowpermute(C1, R1c);
     340      398141 :   A12 = rowpermute(A2, R1);
     341      398141 :   A22 = rowpermute(A2, R1c);
     342      398141 :   M12 = Flm_rsolve_lower_unit_pre(C11, A12, p, pi);
     343      398141 :   B2 = Flm_sub(A22, Flm_mul_pre(C21, M12, p, pi), p);
     344      398141 :   r2 = Flm_echelon_pre(B2, &R2, &C2, p, pi);
     345      398141 :   if (!r2) { *R = R1; *C = C1; r = r1; }
     346             :   else
     347             :   {
     348      361294 :     R2 = perm_mul(R1c, R2);
     349      361294 :     C2 = rowpermute(vconcat(zero_Flm(r1, r2), C2),
     350             :                     perm_inv(vecsmall_concat(R1, R1c)));
     351      361294 :     r = r1 + r2;
     352      361294 :     *R = cgetg(r + 1, t_VECSMALL);
     353      361294 :     *C = cgetg(r + 1, t_MAT);
     354     4579442 :     for (j = j1 = j2 = 1; j <= r; j++)
     355     4218148 :       if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
     356             :       {
     357     2268043 :         gel(*C, j) = gel(C1, j1);
     358     2268043 :         (*R)[j] = R1[j1++];
     359             :       }
     360             :       else
     361             :       {
     362     1950105 :         gel(*C, j) = gel(C2, j2);
     363     1950105 :         (*R)[j] = R2[j2++];
     364             :       }
     365             :   }
     366      398141 :   if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
     367      398141 :   return r;
     368             : }
     369             : 
     370             : static void /* assume m < p */
     371     5140775 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
     372             : {
     373     5140775 :   uel(b,k) = Fl_addmul_pre(uel(b, k), m, uel(b, i), p, pi);
     374     5140783 : }
     375             : static void /* same m = 1 */
     376      177553 : _Fl_add(GEN b, long k, long i, ulong p)
     377             : {
     378      177553 :   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
     379      177553 : }
     380             : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     381     2035710 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
     382             : {
     383     2035710 :   uel(b,k) += m * uel(b,i);
     384     2035710 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     385     2035710 : }
     386             : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     387      559794 : _Fl_add_OK(GEN b, long k, long i, ulong p)
     388             : {
     389      559794 :   uel(b,k) += uel(b,i);
     390      559794 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     391      559794 : }
     392             : 
     393             : /* assume 0 <= a[i,j] < p */
     394             : static GEN
     395      134013 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
     396             : {
     397      134013 :   GEN u = cgetg(li+1,t_VECSMALL);
     398      134013 :   ulong m = uel(b,li) % p;
     399             :   long i,j;
     400             : 
     401      134013 :   uel(u,li) = (m * ucoeff(a,li,li)) % p;
     402      518378 :   for (i = li-1; i > 0; i--)
     403             :   {
     404      384365 :     m = p - uel(b,i)%p;
     405     1295029 :     for (j = i+1; j <= li; j++) {
     406      910664 :       if (m & HIGHBIT) m %= p;
     407      910664 :       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
     408             :     }
     409      384365 :     m %= p;
     410      384365 :     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
     411      384365 :     uel(u,i) = m;
     412             :   }
     413      134013 :   return u;
     414             : }
     415             : static GEN
     416     1084676 : Fl_get_col(GEN a, GEN b, long li, ulong p)
     417             : {
     418     1084676 :   GEN u = cgetg(li+1,t_VECSMALL);
     419     1084676 :   ulong m = uel(b,li) % p;
     420             :   long i,j;
     421             : 
     422     1084676 :   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
     423     3249827 :   for (i=li-1; i>0; i--)
     424             :   {
     425     2165151 :     m = b[i]%p;
     426     5992246 :     for (j = i+1; j <= li; j++)
     427     3827095 :       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
     428     2165151 :     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
     429     2165151 :     uel(u,i) = m;
     430             :   }
     431     1084676 :   return u;
     432             : }
     433             : 
     434             : static GEN
     435      129353 : Flm_ker_gauss_OK(GEN x, ulong p, long deplin)
     436             : {
     437             :   GEN y, c, d;
     438             :   long i, j, k, r, t, m, n;
     439             :   ulong a;
     440             : 
     441      129353 :   n = lg(x)-1;
     442      129353 :   m=nbrows(x); r=0;
     443             : 
     444      129353 :   c = zero_zv(m);
     445      129353 :   d = cgetg(n+1, t_VECSMALL);
     446      129353 :   a = 0; /* for gcc -Wall */
     447      628288 :   for (k=1; k<=n; k++)
     448             :   {
     449     2144964 :     for (j=1; j<=m; j++)
     450     1923219 :       if (!c[j])
     451             :       {
     452     1125540 :         a = ucoeff(x,j,k) % p;
     453     1125540 :         if (a) break;
     454             :       }
     455      522717 :     if (j > m)
     456             :     {
     457      221745 :       if (deplin==1) {
     458       23782 :         c = cgetg(n+1, t_VECSMALL);
     459       83001 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
     460       54876 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     461       23782 :         return c;
     462             :       }
     463      197963 :       r++; d[k] = 0;
     464             :     }
     465             :     else
     466             :     {
     467      300972 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     468      300972 :       c[j] = k; d[k] = j;
     469      300972 :       ucoeff(x,j,k) = p-1;
     470      300972 :       if (piv != 1)
     471      984055 :         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
     472     1962217 :       for (t=1; t<=m; t++)
     473             :       {
     474     1661245 :         if (t == j) continue;
     475             : 
     476     1360273 :         piv = ( ucoeff(x,t,k) %= p );
     477     1360273 :         if (!piv) continue;
     478      679822 :         if (piv == 1)
     479      613675 :           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
     480             :         else
     481     2114429 :           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
     482             :       }
     483             :     }
     484             :   }
     485      105571 :   if (deplin==1) return NULL;
     486             : 
     487      105564 :   y = cgetg(r+1, t_MAT);
     488      303527 :   for (j=k=1; j<=r; j++,k++)
     489             :   {
     490      197963 :     GEN C = cgetg(n+1, t_VECSMALL);
     491             : 
     492      397789 :     gel(y,j) = C; while (d[k]) k++;
     493      730306 :     for (i=1; i<k; i++)
     494      532343 :       if (d[i])
     495      351809 :         uel(C,i) = ucoeff(x,d[i],k) % p;
     496             :       else
     497      180534 :         uel(C,i) = 0UL;
     498      444218 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     499             :   }
     500      105564 :   if (deplin == 2) {
     501           0 :     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
     502           0 :     for (i = j = 1; j <= n; j++)
     503           0 :       if (d[j]) pc[i++] = j;
     504           0 :     return mkvec2(y, pc);
     505             :   }
     506      105564 :   return y;
     507             : }
     508             : 
     509             : /* in place, destroy x */
     510             : static GEN
     511      272938 : Flm_ker_gauss(GEN x, ulong p, long deplin)
     512             : {
     513             :   GEN y, c, d;
     514             :   long i, j, k, r, t, m, n;
     515             :   ulong a, pi;
     516      272938 :   n = lg(x)-1;
     517      272938 :   if (!n) return cgetg(1,t_MAT);
     518      272938 :   if (SMALL_ULONG(p)) return Flm_ker_gauss_OK(x, p, deplin);
     519      143585 :   pi = get_Fl_red(p);
     520             : 
     521      143585 :   m=nbrows(x); r=0;
     522             : 
     523      143583 :   c = zero_zv(m);
     524      143578 :   d = cgetg(n+1, t_VECSMALL);
     525      143556 :   a = 0; /* for gcc -Wall */
     526      461310 :   for (k=1; k<=n; k++)
     527             :   {
     528      914072 :     for (j=1; j<=m; j++)
     529      812914 :       if (!c[j])
     530             :       {
     531      635285 :         a = ucoeff(x,j,k);
     532      635285 :         if (a) break;
     533             :       }
     534      317727 :     if (j > m)
     535             :     {
     536      101177 :       if (deplin==1) {
     537           7 :         c = cgetg(n+1, t_VECSMALL);
     538          21 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k);
     539           7 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     540           7 :         return c;
     541             :       }
     542      101170 :       r++; d[k] = 0;
     543             :     }
     544             :     else
     545             :     {
     546      216550 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     547      216575 :       c[j] = k; d[k] = j;
     548      216575 :       ucoeff(x,j,k) = p-1;
     549      216575 :       if (piv != 1)
     550      396524 :         for (i=k+1; i<=n; i++)
     551      182950 :           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
     552      938466 :       for (t=1; t<=m; t++)
     553             :       {
     554      721882 :         if (t == j) continue;
     555             : 
     556      505305 :         piv = ucoeff(x,t,k);
     557      505305 :         if (!piv) continue;
     558      287224 :         if (piv == 1)
     559       30505 :           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
     560             :         else
     561      720389 :           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
     562             :       }
     563             :     }
     564             :   }
     565      143583 :   if (deplin==1) return NULL;
     566             : 
     567      143576 :   y = cgetg(r+1, t_MAT);
     568      244733 :   for (j=k=1; j<=r; j++,k++)
     569             :   {
     570      101167 :     GEN C = cgetg(n+1, t_VECSMALL);
     571             : 
     572      161891 :     gel(y,j) = C; while (d[k]) k++;
     573      211630 :     for (i=1; i<k; i++)
     574      110471 :       if (d[i])
     575       75618 :         uel(C,i) = ucoeff(x,d[i],k);
     576             :       else
     577       34853 :         uel(C,i) = 0UL;
     578      162793 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     579             :   }
     580      143566 :   if (deplin == 2) {
     581      139912 :     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
     582      434465 :     for (i = j = 1; j <= n; j++)
     583      294557 :       if (d[j]) pc[i++] = j;
     584      139908 :     return mkvec2(y, pc);
     585             :   }
     586        3654 :   return y;
     587             : }
     588             : 
     589             : GEN
     590           0 : Flm_intersect(GEN x, GEN y, ulong p)
     591             : {
     592           0 :   pari_sp av = avma;
     593           0 :   long j, lx = lg(x);
     594             :   GEN z;
     595             : 
     596           0 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     597           0 :   z = Flm_ker(shallowconcat(x,y), p);
     598           0 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     599           0 :   return gerepileupto(av, Flm_mul(x,z,p));
     600             : }
     601             : 
     602             : static GEN
     603      150684 : Flm_ker_echelon(GEN x, ulong p, long pivots) {
     604      150684 :   pari_sp av = avma;
     605             :   GEN R, Rc, C, C1, C2, S, K;
     606      150684 :   long n = lg(x) - 1, r;
     607      150684 :   ulong pi = get_Fl_red(p);
     608      150684 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     609      150684 :   Rc = indexcompl(R, n);
     610      150684 :   C1 = rowpermute(C, R);
     611      150684 :   C2 = rowpermute(C, Rc);
     612      150684 :   S = Flm_lsolve_lower_unit_pre(C1, C2, p, pi);
     613      150684 :   K = vecpermute(shallowconcat(Flm_neg(S, p), matid_Flm(n - r)),
     614             :                  perm_inv(vecsmall_concat(R, Rc)));
     615      150684 :   K = Flm_transpose(K);
     616      150684 :   if (pivots)
     617       31426 :     return gerepilecopy(av, mkvec2(K, R));
     618      119258 :   return gerepileupto(av, K);
     619             : }
     620             : 
     621             : static GEN
     622       23246 : Flm_deplin_echelon(GEN x, ulong p) {
     623       23246 :   pari_sp av = avma;
     624             :   GEN R, Rc, C, C1, C2, s, v;
     625       23246 :   long i, n = lg(x) - 1, r;
     626       23246 :   ulong pi = get_Fl_red(p);
     627       23246 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     628       23246 :   if (r == n) return gc_NULL(av);
     629       23239 :   Rc = indexcompl(R, n);
     630       23239 :   i = Rc[1];
     631       23239 :   C1 = rowpermute(C, R);
     632       23239 :   C2 = rowslice(C, i, i);
     633       23239 :   s = Flm_row(Flm_lsolve_lower_unit_pre(C1, C2, p, pi), 1);
     634       23239 :   v = vecsmallpermute(vecsmall_concat(Flv_neg(s, p), vecsmall_ei(n - r, 1)),
     635             :                  perm_inv(vecsmall_concat(R, Rc)));
     636       23239 :   return gerepileuptoleaf(av, v);
     637             : }
     638             : 
     639             : static GEN
     640      446868 : Flm_ker_i(GEN x, ulong p, long deplin, long inplace) {
     641      446868 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     642      173930 :     switch(deplin) {
     643      119258 :     case 0: return Flm_ker_echelon(x, p, 0);
     644       23246 :     case 1: return Flm_deplin_echelon(x, p);
     645       31426 :     case 2: return Flm_ker_echelon(x, p, 1);
     646             :     }
     647      272938 :   return Flm_ker_gauss(inplace? x: Flm_copy(x), p, deplin);
     648             : }
     649             : 
     650             : GEN
     651      386520 : Flm_ker_sp(GEN x, ulong p, long deplin) {
     652      386520 :   return Flm_ker_i(x, p, deplin, 1);
     653             : }
     654             : 
     655             : GEN
     656       60349 : Flm_ker(GEN x, ulong p) {
     657       60349 :   return Flm_ker_i(x, p, 0, 0);
     658             : }
     659             : 
     660             : GEN
     661           0 : Flm_deplin(GEN x, ulong p) {
     662           0 :   return Flm_ker_i(x, p, 1, 0);
     663             : }
     664             : 
     665             : /* in place, destroy a, SMALL_ULONG(p) is TRUE */
     666             : static ulong
     667        1743 : Flm_det_gauss_OK(GEN a, long nbco, ulong p)
     668             : {
     669        1743 :   long i,j,k, s = 1;
     670        1743 :   ulong q, x = 1;
     671             : 
     672        7266 :   for (i=1; i<nbco; i++)
     673             :   {
     674        7140 :     for(k=i; k<=nbco; k++)
     675             :     {
     676        6958 :       ulong c = ucoeff(a,k,i) % p;
     677        6958 :       ucoeff(a,k,i) = c;
     678        6958 :       if (c) break;
     679             :     }
     680       17528 :     for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
     681        5705 :     if (k > nbco) return ucoeff(a,i,i);
     682        5523 :     if (k != i)
     683             :     { /* exchange the lines s.t. k = i */
     684        3122 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     685         784 :       s = -s;
     686             :     }
     687        5523 :     q = ucoeff(a,i,i);
     688             : 
     689        5523 :     if (x & HIGHMASK) x %= p;
     690        5523 :     x *= q;
     691        5523 :     q = Fl_inv(q,p);
     692       18417 :     for (k=i+1; k<=nbco; k++)
     693             :     {
     694       12894 :       ulong m = ucoeff(a,i,k) % p;
     695       12894 :       if (!m) continue;
     696             : 
     697        8708 :       m = p - ((m*q)%p);
     698       34489 :       for (j=i+1; j<=nbco; j++)
     699             :       {
     700       25781 :         ulong c = ucoeff(a,j,k);
     701       25781 :         if (c & HIGHMASK) c %= p;
     702       25781 :         ucoeff(a,j,k) = c  + m*ucoeff(a,j,i);
     703             :       }
     704             :     }
     705             :   }
     706        1561 :   if (x & HIGHMASK) x %= p;
     707        1561 :   q = ucoeff(a,nbco,nbco);
     708        1561 :   if (q & HIGHMASK) q %= p;
     709        1561 :   x = (x*q) % p;
     710        1561 :   if (s < 0 && x) x = p - x;
     711        1561 :   return x;
     712             : }
     713             : 
     714             : /* in place, destroy a */
     715             : static ulong
     716       50567 : Flm_det_gauss(GEN a, ulong p)
     717             : {
     718       50567 :   long i,j,k, s = 1, nbco = lg(a)-1;
     719       50567 :   ulong pi, q, x = 1;
     720             : 
     721       50567 :   if (SMALL_ULONG(p)) return Flm_det_gauss_OK(a, nbco, p);
     722       48824 :   pi = get_Fl_red(p);
     723      314782 :   for (i=1; i<nbco; i++)
     724             :   {
     725      266360 :     for(k=i; k<=nbco; k++)
     726      266361 :       if (ucoeff(a,k,i)) break;
     727      265957 :     if (k > nbco) return ucoeff(a,i,i);
     728      265957 :     if (k != i)
     729             :     { /* exchange the lines s.t. k = i */
     730        1090 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     731         212 :       s = -s;
     732             :     }
     733      265957 :     q = ucoeff(a,i,i);
     734             : 
     735      265957 :     x = Fl_mul_pre(x, q, p, pi);
     736      265958 :     q = Fl_inv(q,p);
     737     1145276 :     for (k=i+1; k<=nbco; k++)
     738             :     {
     739      879318 :       ulong m = ucoeff(a,i,k);
     740      879318 :       if (!m) continue;
     741             : 
     742      839544 :       m = Fl_mul_pre(m, q, p, pi);
     743     4316392 :       for (j=i+1; j<=nbco; j++)
     744     3476852 :         ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
     745             :     }
     746             :   }
     747       48825 :   if (s < 0) x = Fl_neg(x, p);
     748       48825 :   return Fl_mul(x, ucoeff(a,nbco,nbco), p);
     749             : }
     750             : 
     751             : static ulong
     752       35577 : Flm_det_CUP(GEN a, ulong p) {
     753             :   GEN R, C, U, P;
     754       35577 :   long i, n = lg(a) - 1, r;
     755             :   ulong d;
     756       35577 :   ulong pi = get_Fl_red(p);
     757       35577 :   r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi);
     758       35579 :   if (r < n)
     759          42 :     d = 0;
     760             :   else {
     761       35537 :     d = perm_sign(P) == 1? 1: p-1;
     762      376406 :     for (i = 1; i <= n; i++)
     763      340881 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
     764             :   }
     765       35567 :   return d;
     766             : }
     767             : 
     768             : static ulong
     769       86145 : Flm_det_i(GEN x, ulong p, long inplace) {
     770       86145 :   pari_sp av = avma;
     771             :   ulong d;
     772       86145 :   if (lg(x) - 1 >= Flm_CUP_LIMIT)
     773       35578 :     d = Flm_det_CUP(x, p);
     774             :   else
     775       50567 :     d = Flm_det_gauss(inplace? x: Flm_copy(x), p);
     776       86138 :   return gc_ulong(av, d);
     777             : }
     778             : 
     779             : ulong
     780       86145 : Flm_det_sp(GEN x, ulong p) { return Flm_det_i(x, p, 1); }
     781             : ulong
     782           0 : Flm_det(GEN x, ulong p) { return Flm_det_i(x, p, 0); }
     783             : 
     784             : /* Destroy x */
     785             : static GEN
     786      199357 : Flm_gauss_pivot(GEN x, ulong p, long *rr)
     787             : {
     788             :   GEN c,d;
     789             :   long i,j,k,r,t,n,m;
     790             : 
     791      199357 :   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
     792             : 
     793      199245 :   m=nbrows(x); r=0;
     794      199245 :   d=cgetg(n+1,t_VECSMALL);
     795      199245 :   c = zero_zv(m);
     796      923044 :   for (k=1; k<=n; k++)
     797             :   {
     798     2584821 :     for (j=1; j<=m; j++)
     799     2384728 :       if (!c[j])
     800             :       {
     801     1315517 :         ucoeff(x,j,k) %= p;
     802     1315517 :         if (ucoeff(x,j,k)) break;
     803             :       }
     804      723799 :     if (j>m) { r++; d[k]=0; }
     805             :     else
     806             :     {
     807      523706 :       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
     808      523706 :       c[j]=k; d[k]=j;
     809     1656659 :       for (i=k+1; i<=n; i++)
     810     1132953 :         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
     811     4335662 :       for (t=1; t<=m; t++)
     812     3811956 :         if (!c[t]) /* no pivot on that line yet */
     813             :         {
     814     2607287 :           piv = ucoeff(x,t,k);
     815     2607287 :           if (piv)
     816             :           {
     817     1450027 :             ucoeff(x,t,k) = 0;
     818     4517449 :             for (i=k+1; i<=n; i++)
     819     3067422 :               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
     820     3067422 :                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
     821             :           }
     822             :         }
     823     2180365 :       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
     824             :     }
     825             :   }
     826      199245 :   *rr = r; return gc_const((pari_sp)d, d);
     827             : }
     828             : 
     829             : static GEN
     830       74255 : Flm_pivots_CUP(GEN x, ulong p, long *rr)
     831             : {
     832       74255 :   long i, n = lg(x) - 1, r;
     833       74255 :   GEN R, C, U, P, d = zero_zv(n);
     834       74255 :   ulong pi = get_Fl_red(p);
     835       74255 :   r = Flm_CUP_pre(x, &R, &C, &U, &P, p, pi);
     836      888150 :   for(i = 1; i <= r; i++)
     837      813895 :     d[P[i]] = R[i];
     838       74255 :   *rr = n - r; return gc_const((pari_sp)d, d);
     839             : }
     840             : 
     841             : GEN
     842      273612 : Flm_pivots(GEN x, ulong p, long *rr, long inplace)
     843             : {
     844      273612 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     845       74255 :     return Flm_pivots_CUP(x, p, rr);
     846      199357 :   return Flm_gauss_pivot(inplace? x: Flm_copy(x), p, rr);
     847             : }
     848             : 
     849             : long
     850       15183 : Flm_rank(GEN x, ulong p)
     851             : {
     852       15183 :   pari_sp av = avma;
     853             :   long r;
     854       15183 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT) {
     855             :     GEN R, C;
     856        7707 :     ulong pi = get_Fl_red(p);
     857        7707 :     return gc_long(av, Flm_echelon_pre(x, &R, &C, p, pi));
     858             :   }
     859        7476 :   (void) Flm_pivots(x, p, &r, 0);
     860        7476 :   return gc_long(av, lg(x)-1 - r);
     861             : }
     862             : 
     863             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
     864             :  * reduced mod p */
     865             : static GEN
     866          28 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
     867             : {
     868          28 :   long n = lg(A)-1, i = index, j;
     869          28 :   GEN u = const_vecsmall(n, 0);
     870          28 :   u[i] = 1;
     871          28 :   if (SMALL_ULONG(p))
     872          21 :     for (i--; i>0; i--)
     873             :     {
     874           7 :       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
     875           7 :       for (j=i+2; j<=n; j++)
     876             :       {
     877           0 :         if (m & HIGHMASK) m %= p;
     878           0 :         m += ucoeff(A,i,j) * uel(u,j);
     879             :       }
     880           7 :       u[i] = Fl_neg(m % p, p);
     881             :     }
     882             :   else
     883          21 :     for (i--; i>0; i--)
     884             :     {
     885           7 :       ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
     886           7 :       for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
     887           7 :       u[i] = Fl_neg(m, p);
     888             :     }
     889          28 :   return u;
     890             : }
     891             : static GEN
     892          14 : Flm_inv_upper_1(GEN A, ulong p)
     893             : {
     894             :   long i, l;
     895          14 :   GEN B = cgetg_copy(A, &l);
     896          42 :   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
     897          14 :   return B;
     898             : }
     899             : /* destroy a, b */
     900             : static GEN
     901       35926 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
     902             : {
     903       35926 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     904       35926 :   ulong det = 1;
     905             :   GEN u;
     906             : 
     907       35926 :   li = nbrows(a);
     908       35926 :   bco = lg(b)-1;
     909      120599 :   for (i=1; i<=aco; i++)
     910             :   {
     911             :     ulong invpiv;
     912             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
     913      301876 :     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
     914      139216 :     for (k = i; k <= li; k++)
     915             :     {
     916      139209 :       ulong piv = ( ucoeff(a,k,i) %= p );
     917      139209 :       if (piv)
     918             :       {
     919      120592 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     920      120592 :         if (detp)
     921             :         {
     922           0 :           if (det & HIGHMASK) det %= p;
     923           0 :           det *= piv;
     924             :         }
     925      120592 :         break;
     926             :       }
     927             :     }
     928             :     /* found a pivot on line k */
     929      120599 :     if (k > li) return NULL;
     930      120592 :     if (k != i)
     931             :     { /* swap lines so that k = i */
     932       12507 :       s = -s;
     933       55042 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     934       78906 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     935             :     }
     936      120592 :     if (i == aco) break;
     937             : 
     938       84673 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
     939      287817 :     for (k=i+1; k<=li; k++)
     940             :     {
     941      203144 :       ulong m = ( ucoeff(a,k,i) %= p );
     942      203144 :       if (!m) continue;
     943             : 
     944       64520 :       m = Fl_mul(m, invpiv, p);
     945       64520 :       if (m == 1) {
     946       45405 :         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
     947       71441 :         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
     948             :       } else {
     949      203747 :         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
     950      355669 :         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
     951             :       }
     952             :     }
     953             :   }
     954       35919 :   if (detp)
     955             :   {
     956           0 :     det %=  p;
     957           0 :     if (s < 0 && det) det = p - det;
     958           0 :     *detp = det;
     959             :   }
     960       35919 :   u = cgetg(bco+1,t_MAT);
     961      169932 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
     962       35919 :   return u;
     963             : }
     964             : 
     965             : /* destroy a, b */
     966             : static GEN
     967      436747 : Flm_gauss_sp_i(GEN a, GEN b, ulong *detp, ulong p)
     968             : {
     969      436747 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     970      436747 :   ulong det = 1;
     971             :   GEN u;
     972             :   ulong pi;
     973      436747 :   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
     974      436747 :   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
     975      400821 :   pi = get_Fl_red(p);
     976      400821 :   li = nbrows(a);
     977      400821 :   bco = lg(b)-1;
     978     1084671 :   for (i=1; i<=aco; i++)
     979             :   {
     980             :     ulong invpiv;
     981             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
     982     1135369 :     for (k = i; k <= li; k++)
     983             :     {
     984     1135362 :       ulong piv = ucoeff(a,k,i);
     985     1135362 :       if (piv)
     986             :       {
     987     1084664 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     988     1084664 :         if (detp) det = Fl_mul_pre(det, piv, p, pi);
     989     1084664 :         break;
     990             :       }
     991             :     }
     992             :     /* found a pivot on line k */
     993     1084671 :     if (k > li) return NULL;
     994     1084664 :     if (k != i)
     995             :     { /* swap lines so that k = i */
     996       45323 :       s = -s;
     997      168703 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     998      195013 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     999             :     }
    1000     1084664 :     if (i == aco) break;
    1001             : 
    1002      683850 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    1003     1766517 :     for (k=i+1; k<=li; k++)
    1004             :     {
    1005     1082667 :       ulong m = ucoeff(a,k,i);
    1006     1082667 :       if (!m) continue;
    1007             : 
    1008      899238 :       m = Fl_mul_pre(m, invpiv, p, pi);
    1009      899238 :       if (m == 1) {
    1010       88893 :         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
    1011      131057 :         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
    1012             :       } else {
    1013     2554947 :         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
    1014     3878237 :         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
    1015             :       }
    1016             :     }
    1017             :   }
    1018      400814 :   if (detp)
    1019             :   {
    1020           0 :     if (s < 0 && det) det = p - det;
    1021           0 :     *detp = det;
    1022             :   }
    1023      400814 :   u = cgetg(bco+1,t_MAT);
    1024     1485490 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
    1025      400814 :   return u;
    1026             : }
    1027             : 
    1028             : static GEN
    1029      223746 : Flm_gauss_from_CUP(GEN b, GEN R, GEN C, GEN U, GEN P, ulong p, ulong pi, ulong *detp)
    1030             : {
    1031      223746 :   GEN Y = Flm_rsolve_lower_unit_pre(rowpermute(C, R), rowpermute(b, R), p, pi);
    1032      223747 :   GEN X = rowpermute(Flm_rsolve_upper_pre(U, Y, p, pi), perm_inv(P));
    1033      223751 :   if (detp)
    1034             :   {
    1035      175238 :     ulong d = perm_sign(P) == 1? 1: p-1;
    1036      175237 :     long i, r = lg(R);
    1037     2112070 :     for (i = 1; i < r; i++)
    1038     1936839 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
    1039      175231 :     *detp = d;
    1040             :   }
    1041      223744 :   return X;
    1042             : }
    1043             : 
    1044             : static GEN
    1045       48522 : Flm_gauss_CUP(GEN a, GEN b, ulong *detp, ulong p) {
    1046             :   GEN R, C, U, P;
    1047       48522 :   long n = lg(a) - 1, r;
    1048       48522 :   ulong pi = get_Fl_red(p);
    1049       48522 :   if (nbrows(a) < n || (r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi)) < n)
    1050          14 :     return NULL;
    1051       48508 :   return Flm_gauss_from_CUP(b, R, C, U, P, p, pi, detp);
    1052             : }
    1053             : 
    1054             : GEN
    1055      457710 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p) {
    1056      457710 :   pari_sp av = avma;
    1057             :   GEN x;
    1058      457710 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1059       20963 :     x = Flm_gauss_CUP(a, b, detp, p);
    1060             :   else
    1061      436747 :     x = Flm_gauss_sp_i(a, b, detp, p);
    1062      457710 :   if (!x) return gc_NULL(av);
    1063      457696 :   return gerepileupto(av, x);
    1064             : }
    1065             : 
    1066             : GEN
    1067          63 : Flm_gauss(GEN a, GEN b, ulong p) {
    1068          63 :   pari_sp av = avma;
    1069             :   GEN x;
    1070          63 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1071          14 :     x = Flm_gauss_CUP(a, b, NULL, p);
    1072             :   else {
    1073          49 :     a = RgM_shallowcopy(a);
    1074          49 :     b = RgM_shallowcopy(b);
    1075          49 :     x = Flm_gauss_sp(a, b, NULL, p);
    1076             :   }
    1077          63 :   if (!x) return gc_NULL(av);
    1078          49 :   return gerepileupto(av, x);
    1079             : }
    1080             : 
    1081             : static GEN
    1082      439841 : Flm_inv_i(GEN a, ulong *detp, ulong p, long inplace) {
    1083      439841 :   pari_sp av = avma;
    1084      439841 :   long n = lg(a) - 1;
    1085             :   GEN b, x;
    1086      439841 :   if (n == 0) return cgetg(1, t_MAT);
    1087      439841 :   b = matid_Flm(nbrows(a));
    1088      439841 :   if (n >= Flm_CUP_LIMIT)
    1089       27545 :     x = Flm_gauss_CUP(a, b, detp, p);
    1090             :   else {
    1091      412296 :     if (!inplace)
    1092       11494 :       a = RgM_shallowcopy(a);
    1093      412296 :     x = Flm_gauss_sp(a, b, detp, p);
    1094             :   }
    1095      439841 :   if (!x) return gc_NULL(av);
    1096      439827 :   return gerepileupto(av, x);
    1097             : }
    1098             : 
    1099             : GEN
    1100      416825 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
    1101      416825 :   return Flm_inv_i(a, detp, p, 1);
    1102             : }
    1103             : 
    1104             : GEN
    1105       23016 : Flm_inv(GEN a, ulong p) {
    1106       23016 :   return Flm_inv_i(a, NULL, p, 0);
    1107             : }
    1108             : 
    1109             : GEN
    1110          21 : Flm_Flc_gauss(GEN a, GEN b, ulong p) {
    1111          21 :   pari_sp av = avma;
    1112          21 :   GEN z = Flm_gauss(a, mkmat(b), p);
    1113          21 :   if (!z) return gc_NULL(av);
    1114          14 :   if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
    1115          14 :   return gerepileuptoleaf(av, gel(z,1));
    1116             : }
    1117             : 
    1118             : GEN
    1119      175253 : Flm_adjoint(GEN A, ulong p)
    1120             : {
    1121      175253 :   pari_sp av = avma;
    1122             :   GEN R, C, U, P, C1, U1, v, c, d;
    1123      175253 :   long r, i, q, n = lg(A)-1, m;
    1124             :   ulong D;
    1125      175253 :   ulong pi = get_Fl_red(p);
    1126      175253 :   if (n == 0) return cgetg(1, t_MAT);
    1127      175253 :   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
    1128      175249 :   m = nbrows(A);
    1129      175249 :   if (r == n)
    1130             :   {
    1131      175235 :     GEN X = Flm_gauss_from_CUP(matid_Flm(m), R, C, U, P, p, pi, &D);
    1132      175236 :     return gerepileupto(av, Flm_Fl_mul_pre(X, D, p, pi));
    1133             :   }
    1134          14 :   if (r < n-1) return zero_Flm(n, m);
    1135          28 :   for (q = n, i = 1; i < n; i++)
    1136          14 :     if (R[i] != i) { q = i; break; }
    1137          14 :   C1 = matslice(C, 1, q-1, 1, q-1);
    1138          14 :   c = vecslice(Flm_row(C, q), 1, q-1);
    1139          14 :   c = Flm_lsolve_lower_unit_pre(C1, Flm_transpose(mkmat(c)), p, pi);
    1140          14 :   d = cgetg(m+1, t_VECSMALL);
    1141          28 :   for (i=1; i<q; i++)    uel(d,i) = ucoeff(c,1,i);
    1142          14 :   uel(d,q) = p-1;
    1143          21 :   for (i=q+1; i<=m; i++) uel(d,i) = 0;
    1144          14 :   U1 = vecslice(U, 1, n-1);
    1145          14 :   v = gel(Flm_rsolve_upper_pre(U1, mkmat(gel(U,n)), p, pi),1);
    1146          14 :   v = vecsmall_append(v, p-1);
    1147          14 :   D = perm_sign(P) != (odd(q+n)?-1:1) ? p-1 : 1;
    1148          28 :   for (i = 1; i <= n-1; i++)
    1149          14 :     D = Fl_mul_pre(D, ucoeff(U1, i, i), p, pi);
    1150          14 :   d = Flv_Fl_mul(d, D, p);
    1151          14 :   return rowpermute(Flc_Flv_mul(v, d, p),perm_inv(P));
    1152             : }
    1153             : 
    1154             : static GEN
    1155          21 : Flm_invimage_CUP(GEN A, GEN B, ulong p) {
    1156          21 :   pari_sp av = avma;
    1157             :   GEN R, Rc, C, U, P, B1, B2, C1, C2, X, Y, Z;
    1158             :   long r;
    1159          21 :   ulong pi = get_Fl_red(p);
    1160          21 :   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
    1161          21 :   Rc = indexcompl(R, nbrows(B));
    1162          21 :   C1 = rowpermute(C, R);
    1163          21 :   C2 = rowpermute(C, Rc);
    1164          21 :   B1 = rowpermute(B, R);
    1165          21 :   B2 = rowpermute(B, Rc);
    1166          21 :   Z = Flm_rsolve_lower_unit_pre(C1, B1, p, pi);
    1167          21 :   if (!gequal(Flm_mul_pre(C2, Z, p, pi), B2))
    1168          14 :     return NULL;
    1169           7 :   Y = vconcat(Flm_rsolve_upper_pre(vecslice(U, 1, r), Z, p, pi),
    1170           7 :               zero_Flm(lg(A) - 1 - r, lg(B) - 1));
    1171           7 :   X = rowpermute(Y, perm_inv(P));
    1172           7 :   return gerepileupto(av, X);
    1173             : }
    1174             : 
    1175             : GEN
    1176          42 : Flm_invimage_i(GEN A, GEN B, ulong p)
    1177             : {
    1178             :   GEN d, x, X, Y;
    1179          42 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    1180             : 
    1181          42 :   if (!nB) return cgetg(1, t_MAT);
    1182          42 :   if (nA + nB >= Flm_CUP_LIMIT && nbrows(B) >= Flm_CUP_LIMIT)
    1183          21 :     return Flm_invimage_CUP(A, B, p);
    1184             : 
    1185          21 :   x = Flm_ker_sp(shallowconcat(Flm_neg(A,p), B), p, 0);
    1186             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    1187             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    1188             :    * Y has at least nB columns and full rank */
    1189          21 :   nY = lg(x)-1;
    1190          21 :   if (nY < nB) return NULL;
    1191          21 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    1192          21 :   d = cgetg(nB+1, t_VECSMALL);
    1193          56 :   for (i = nB, j = nY; i >= 1; i--, j--)
    1194             :   {
    1195          49 :     for (; j>=1; j--)
    1196          42 :       if (coeff(Y,i,j)) { d[i] = j; break; }
    1197          42 :     if (!j) return NULL;
    1198             :   }
    1199             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    1200          14 :   Y = vecpermute(Y, d);
    1201          14 :   x = vecpermute(x, d);
    1202          14 :   X = rowslice(x, 1, nA);
    1203          14 :   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
    1204             : }
    1205             : GEN
    1206           0 : Flm_invimage(GEN A, GEN B, ulong p)
    1207             : {
    1208           0 :   pari_sp av = avma;
    1209           0 :   GEN X = Flm_invimage_i(A,B,p);
    1210           0 :   if (!X) return gc_NULL(av);
    1211           0 :   return gerepileupto(av, X);
    1212             : }
    1213             : 
    1214             : GEN
    1215       40292 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
    1216             : {
    1217       40292 :   pari_sp av = avma;
    1218       40292 :   long i, l = lg(A);
    1219             :   GEN M, x;
    1220             :   ulong t;
    1221             : 
    1222       40292 :   if (l==1) return NULL;
    1223       40292 :   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
    1224       40292 :   M = cgetg(l+1,t_MAT);
    1225      455099 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    1226       40292 :   gel(M,l) = y; M = Flm_ker(M,p);
    1227       40292 :   i = lg(M)-1; if (!i) return gc_NULL(av);
    1228             : 
    1229       40292 :   x = gel(M,i); t = x[l];
    1230       40292 :   if (!t) return gc_NULL(av);
    1231             : 
    1232       40285 :   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
    1233       40285 :   if (t!=1) x = Flv_Fl_mul(x, t, p);
    1234       40285 :   return gerepileuptoleaf(av, x);
    1235             : }

Generated by: LCOV version 1.13