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.12.1 lcov report (development 24988-2584e74448) Lines: 733 762 96.2 %
Date: 2020-01-26 05:57:03 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      257081 : Flv_to_ZV(GEN x)
      19      257081 : { pari_APPLY_type(t_VEC, utoi(x[i])) }
      20             : 
      21             : GEN
      22    23535794 : Flc_to_ZC(GEN x)
      23    23535794 : { pari_APPLY_type(t_COL, utoi(x[i])) }
      24             : 
      25             : GEN
      26     8955387 : Flm_to_ZM(GEN x)
      27     8955387 : { 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      123378 :   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       56166 :   for (i=1; i<l; i++) Flc_to_ZC_inplace(gel(z, i));
      43       56166 :   return z;
      44             : }
      45             : 
      46             : static GEN
      47      733920 : Flm_solve_upper_1(GEN U, GEN B, ulong p, ulong pi)
      48      733920 : { return Flm_Fl_mul_pre(B, Fl_inv(ucoeff(U, 1, 1), p), p, pi); }
      49             : 
      50             : static GEN
      51      982102 : Flm_rsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      52             : {
      53      982102 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      54      982102 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      55      982113 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      56      982112 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
      57      982112 :   GEN B1 = rowslice(B, 1, 1);
      58      982099 :   GEN B2 = rowslice(B, 2, 2);
      59      982103 :   GEN X2 = Flm_Fl_mul_pre(B2, dinv, p, pi);
      60      982067 :   GEN X1 = Flm_Fl_mul_pre(Flm_sub(B1, Flm_Fl_mul_pre(X2, b, p, pi), p),
      61             :                       ainv, p, pi);
      62      982078 :   return vconcat(X1, X2);
      63             : }
      64             : 
      65             : /* solve U*X = B,  U upper triangular and invertible */
      66             : static GEN
      67     2427006 : Flm_rsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
      68             : {
      69     2427006 :   long n = lg(U) - 1, n1;
      70             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
      71     2427006 :   pari_sp av = avma;
      72             : 
      73     2427006 :   if (n == 0) return B;
      74     2426999 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
      75     2089538 :   if (n == 2) return Flm_rsolve_upper_2(U, B, p, pi);
      76     1107438 :   n1 = (n + 1)/2;
      77     1107438 :   U2 = vecslice(U, n1 + 1, n);
      78     1107444 :   U22 = rowslice(U2, n1 + 1, n);
      79     1107433 :   B2 = rowslice(B, n1 + 1, n);
      80     1107427 :   X2 = Flm_rsolve_upper_pre(U22, B2, p, pi);
      81     1107432 :   U12 = rowslice(U2, 1, n1);
      82     1107415 :   B1 = rowslice(B, 1, n1);
      83     1107420 :   B1 = Flm_sub(B1, Flm_mul_pre(U12, X2, p, pi), p);
      84     1107360 :   if (gc_needed(av, 1)) gerepileall(av, 2, &B1, &X2);
      85     1107360 :   U11 = matslice(U, 1,n1, 1,n1);
      86     1107414 :   X1 = Flm_rsolve_upper_pre(U11, B1, p, pi);
      87     1107453 :   X = vconcat(X1, X2);
      88     1107463 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
      89     1107467 :   return X;
      90             : }
      91             : 
      92             : static GEN
      93     1120059 : Flm_lsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      94             : {
      95     1120059 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      96     1120059 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      97     1120083 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      98     1120079 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
      99     1120068 :   GEN B1 = vecslice(B, 1, 1);
     100     1120043 :   GEN B2 = vecslice(B, 2, 2);
     101     1120034 :   GEN X1 = Flm_Fl_mul_pre(B1, ainv, p, pi);
     102     1119986 :   GEN X2 = Flm_Fl_mul_pre(Flm_sub(B2, Flm_Fl_mul_pre(X1, b, p, pi), p),
     103             :                       dinv, p, pi);
     104     1120005 :   return shallowconcat(X1, X2);
     105             : }
     106             : 
     107             : /* solve X*U = B,  U upper triangular and invertible */
     108             : static GEN
     109     2617139 : Flm_lsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
     110             : {
     111     2617139 :   long n = lg(U) - 1, n1;
     112             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
     113     2617139 :   pari_sp av = avma;
     114             : 
     115     2617139 :   if (n == 0) return B;
     116     2617139 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
     117     2220676 :   if (n == 2) return Flm_lsolve_upper_2(U, B, p, pi);
     118     1100613 :   n1 = (n + 1)/2;
     119     1100613 :   U2 = vecslice(U, n1 + 1, n);
     120     1100611 :   U11 = matslice(U, 1,n1, 1,n1);
     121     1100604 :   U12 = rowslice(U2, 1, n1);
     122     1100597 :   U22 = rowslice(U2, n1 + 1, n);
     123     1100591 :   B1 = vecslice(B, 1, n1);
     124     1100588 :   B2 = vecslice(B, n1 + 1, n);
     125     1100582 :   X1 = Flm_lsolve_upper_pre(U11, B1, p, pi);
     126     1100559 :   B2 = Flm_sub(B2, Flm_mul_pre(X1, U12, p, pi), p);
     127     1100550 :   if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
     128     1100550 :   X2 = Flm_lsolve_upper_pre(U22, B2, p, pi);
     129     1100562 :   X = shallowconcat(X1, X2);
     130     1100603 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     131     1100605 :   return X;
     132             : }
     133             : 
     134             : static GEN
     135     1941091 : Flm_rsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     136             : {
     137     1941091 :   GEN X1 = rowslice(A, 1, 1);
     138     1941088 :   GEN X2 = Flm_sub(rowslice(A, 2, 2),
     139     1941088 :                    Flm_Fl_mul_pre(X1, ucoeff(L, 2, 1), p, pi), p);
     140     1941067 :   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     4572575 : Flm_rsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     147             : {
     148     4572575 :   long m = lg(L) - 1, m1, n;
     149             :   GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
     150     4572575 :   pari_sp av = avma;
     151             : 
     152     4572575 :   if (m == 0) return zero_Flm(0, lg(A) - 1);
     153     4572575 :   if (m == 1) return rowslice(A, 1, 1);
     154     3957420 :   if (m == 2) return Flm_rsolve_lower_unit_2(L, A, p, pi);
     155     2016334 :   m1 = (m + 1)/2;
     156     2016334 :   n = nbrows(L);
     157     2016334 :   L1 = vecslice(L, 1, m1);
     158     2016336 :   L11 = rowslice(L1, 1, m1);
     159     2016326 :   L21 = rowslice(L1, m1 + 1, n);
     160     2016314 :   A1 = rowslice(A, 1, m1);
     161     2016320 :   X1 = Flm_rsolve_lower_unit_pre(L11, A1, p, pi);
     162     2016341 :   A2 = rowslice(A, m1 + 1, n);
     163     2016311 :   A2 = Flm_sub(A2, Flm_mul_pre(L21, X1, p, pi), p);
     164     2016285 :   if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
     165     2016285 :   L22 = matslice(L, m1+1,n, m1+1,m);
     166     2016285 :   X2 = Flm_rsolve_lower_unit_pre(L22, A2, p, pi);
     167     2016331 :   X = vconcat(X1, X2);
     168     2016347 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     169     2016348 :   return X;
     170             : }
     171             : 
     172             : static GEN
     173      785937 : Flm_lsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     174             : {
     175      785937 :   GEN X2 = vecslice(A, 2, 2);
     176      785937 :   GEN X1 = Flm_sub(vecslice(A, 1, 1),
     177      785937 :                    Flm_Fl_mul_pre(X2, ucoeff(L, 2, 1), p, pi), p);
     178      785937 :   return shallowconcat(X1, X2);
     179             : }
     180             : 
     181             : /* solve L*X = A,  L square lower triangular with ones on the diagonal */
     182             : static GEN
     183     1982745 : Flm_lsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     184             : {
     185     1982745 :   long m = lg(L) - 1, m1;
     186             :   GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
     187     1982745 :   pari_sp av = avma;
     188             : 
     189     1982745 :   if (m <= 1) return A;
     190     1708104 :   if (m == 2) return Flm_lsolve_lower_unit_2(L, A, p, pi);
     191      922167 :   m1 = (m + 1)/2;
     192      922167 :   L2 = vecslice(L, m1 + 1, m);
     193      922167 :   L22 = rowslice(L2, m1 + 1, m);
     194      922167 :   A2 = vecslice(A, m1 + 1, m);
     195      922167 :   X2 = Flm_lsolve_lower_unit_pre(L22, A2, p, pi);
     196      922167 :   if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
     197      922167 :   L1 = vecslice(L, 1, m1);
     198      922167 :   L21 = rowslice(L1, m1 + 1, m);
     199      922167 :   A1 = vecslice(A, 1, m1);
     200      922167 :   A1 = Flm_sub(A1, Flm_mul_pre(X2, L21, p, pi), p);
     201      922167 :   L11 = rowslice(L1, 1, m1);
     202      922166 :   if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
     203      922166 :   X1 = Flm_lsolve_lower_unit_pre(L11, A1, p, pi);
     204      922167 :   X = shallowconcat(X1, X2);
     205      922167 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     206      922167 :   return X;
     207             : }
     208             : 
     209             : /* destroy A */
     210             : static long
     211     1260556 : Flm_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     212             : {
     213     1260556 :   long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
     214             :   ulong u, v;
     215             : 
     216     1260557 :   if (P) *P = identity_perm(n);
     217     1260556 :   *R = cgetg(m + 1, t_VECSMALL);
     218     6509136 :   for (j = 1, pr = 0; j <= n; j++)
     219             :   {
     220     8559208 :     for (pr++, pc = 0; pr <= m; pr++)
     221             :     {
     222     7919789 :       for (k = j; k <= n; k++) { v = ucoeff(A, pr, k); if (!pc && v) pc = k; }
     223     7919789 :       if (pc) break;
     224             :     }
     225     5887970 :     if (!pc) break;
     226     5248553 :     (*R)[j] = pr;
     227     5248553 :     if (pc != j)
     228             :     {
     229      899999 :       swap(gel(A, j), gel(A, pc));
     230      899999 :       if (P) lswap((*P)[j], (*P)[pc]);
     231             :     }
     232     5248553 :     u = Fl_inv(ucoeff(A, pr, j), p);
     233    39414043 :     for (i = pr + 1; i <= m; i++)
     234             :     {
     235    34165462 :       v = Fl_mul_pre(ucoeff(A, i, j), u, p, pi);
     236    34165362 :       ucoeff(A, i, j) = v;
     237    34165362 :       v = Fl_neg(v, p);
     238   175607278 :       for (k = j + 1; k <= n; k++)
     239   282883652 :         ucoeff(A, i, k) = Fl_addmul_pre(ucoeff(A, i, k),
     240   141441826 :                                         ucoeff(A, pr, k), v, p, pi);
     241             :     }
     242             :   }
     243     1260583 :   setlg(*R, j);
     244     1260553 :   *C = vecslice(A, 1, j - 1);
     245     1260550 :   if (U) *U = rowpermute(A, *R);
     246     1260559 :   return j - 1;
     247             : }
     248             : 
     249             : static const long Flm_CUP_LIMIT = 8;
     250             : 
     251             : static long
     252     1153243 : Flm_CUP_pre(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     253             : {
     254     1153243 :   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     1153241 :   pari_sp av = avma;
     258             : 
     259     1153241 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     260             :     /* destroy A; not called at the outermost recursion level */
     261      730573 :     return Flm_CUP_basecase(A, R, C, U, P, p, pi);
     262      422668 :   m1 = (minss(m, n) + 1)/2;
     263      422667 :   A1 = rowslice(A, 1, m1);
     264      422669 :   A2 = rowslice(A, m1 + 1, m);
     265      422667 :   r1 = Flm_CUP_pre(A1, &R1, &C1, &U1, &P1, p, pi);
     266      422672 :   if (r1 == 0)
     267             :   {
     268        6728 :     r2 = Flm_CUP_pre(A2, &R2, &C2, &U2, &P2, p, pi);
     269        6728 :     *R = cgetg(r2 + 1, t_VECSMALL);
     270        6728 :     for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
     271        6728 :     *C = vconcat(zero_Flm(m1, r2), C2);
     272        6728 :     *U = U2;
     273        6728 :     *P = P2;
     274        6728 :     r = r2;
     275             :   }
     276             :   else
     277             :   {
     278      415944 :     U11 = vecslice(U1, 1, r1);
     279      415944 :     U12 = vecslice(U1, r1 + 1, n);
     280      415943 :     T21 = vecslicepermute(A2, P1, 1, r1);
     281      415943 :     T22 = vecslicepermute(A2, P1, r1 + 1, n);
     282      415942 :     C21 = Flm_lsolve_upper_pre(U11, T21, p, pi);
     283      415941 :     if (gc_needed(av, 1))
     284           0 :       gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
     285      415941 :     B2 = Flm_sub(T22, Flm_mul_pre(C21, U12, p, pi), p);
     286      415938 :     r2 = Flm_CUP_pre(B2, &R2, &C2, &U2, &P2, p, pi);
     287      415944 :     r = r1 + r2;
     288      415944 :     *R = cgetg(r + 1, t_VECSMALL);
     289      415945 :     for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
     290      415945 :     for (;      i <= r; i++)  (*R)[i] = R2[i - r1] + m1;
     291      415945 :     *C = shallowconcat(vconcat(C1, C21),
     292             :                        vconcat(zero_Flm(m1, r2), C2));
     293      415945 :     *U = shallowconcat(vconcat(U11, zero_Flm(r2, r1)),
     294             :                        vconcat(vecpermute(U12, P2), U2));
     295      415944 :     *P = cgetg(n + 1, t_VECSMALL);
     296      415944 :     for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
     297      415944 :     for (; i <= n; i++)       (*P)[i] = P1[P2[i - r1] + r1];
     298             :   }
     299      422672 :   if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
     300      422672 :   return r;
     301             : }
     302             : 
     303             : static long
     304      529974 : Flm_echelon_gauss(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     305      529974 : { 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      466223 : indexcompl(GEN v, long n)
     310             : {
     311      466223 :   long i, j, k, m = lg(v) - 1;
     312      466223 :   GEN w = cgetg(n - m + 1, t_VECSMALL);
     313    10127487 :   for (i = j = k = 1; i <= n; i++)
     314     9661264 :     if (j <= m && v[j] == i) j++; else w[k++] = i;
     315      466223 :   return w;
     316             : }
     317             : 
     318             : /* column echelon form */
     319             : static long
     320      918052 : Flm_echelon_pre(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     321             : {
     322      918052 :   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      918052 :   pari_sp av = avma;
     326             : 
     327      918052 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     328      529974 :     return Flm_echelon_gauss(Flm_copy(A), R, C, p, pi);
     329             : 
     330      388078 :   n1 = (n + 1)/2;
     331      388078 :   A1 = vecslice(A, 1, n1);
     332      388078 :   A2 = vecslice(A, n1 + 1, n);
     333      388078 :   r1 = Flm_echelon_pre(A1, &R1, &C1, p, pi);
     334      388078 :   if (!r1) return Flm_echelon_pre(A2, R, C, p, pi);
     335      331649 :   if (r1 == m) { *R = R1; *C = C1; return r1; }
     336             : 
     337      327805 :   R1c = indexcompl(R1, m);
     338      327805 :   C11 = rowpermute(C1, R1);
     339      327805 :   C21 = rowpermute(C1, R1c);
     340      327805 :   A12 = rowpermute(A2, R1);
     341      327805 :   A22 = rowpermute(A2, R1c);
     342      327805 :   M12 = Flm_rsolve_lower_unit_pre(C11, A12, p, pi);
     343      327805 :   B2 = Flm_sub(A22, Flm_mul_pre(C21, M12, p, pi), p);
     344      327805 :   r2 = Flm_echelon_pre(B2, &R2, &C2, p, pi);
     345      327805 :   if (!r2) { *R = R1; *C = C1; r = r1; }
     346             :   else
     347             :   {
     348      294762 :     R2 = perm_mul(R1c, R2);
     349      294762 :     C2 = rowpermute(vconcat(zero_Flm(r1, r2), C2),
     350             :                     perm_inv(vecsmall_concat(R1, R1c)));
     351      294762 :     r = r1 + r2;
     352      294762 :     *R = cgetg(r + 1, t_VECSMALL);
     353      294762 :     *C = cgetg(r + 1, t_MAT);
     354     3792561 :     for (j = j1 = j2 = 1; j <= r; j++)
     355     3497799 :       if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
     356             :       {
     357     1892427 :         gel(*C, j) = gel(C1, j1);
     358     1892427 :         (*R)[j] = R1[j1++];
     359             :       }
     360             :       else
     361             :       {
     362     1605372 :         gel(*C, j) = gel(C2, j2);
     363     1605372 :         (*R)[j] = R2[j2++];
     364             :       }
     365             :   }
     366      327805 :   if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
     367      327805 :   return r;
     368             : }
     369             : 
     370             : static void /* assume m < p */
     371     6483309 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
     372             : {
     373     6483309 :   uel(b,k) = Fl_addmul_pre(uel(b, k), m, uel(b, i), p, pi);
     374     6483316 : }
     375             : static void /* same m = 1 */
     376      244427 : _Fl_add(GEN b, long k, long i, ulong p)
     377             : {
     378      244427 :   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
     379      244427 : }
     380             : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     381     1815925 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
     382             : {
     383     1815925 :   uel(b,k) += m * uel(b,i);
     384     1815925 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     385     1815925 : }
     386             : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     387      585904 : _Fl_add_OK(GEN b, long k, long i, ulong p)
     388             : {
     389      585904 :   uel(b,k) += uel(b,i);
     390      585904 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     391      585904 : }
     392             : 
     393             : /* assume 0 <= a[i,j] < p */
     394             : static GEN
     395      126507 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
     396             : {
     397      126507 :   GEN u = cgetg(li+1,t_VECSMALL);
     398      126507 :   ulong m = uel(b,li) % p;
     399             :   long i,j;
     400             : 
     401      126507 :   uel(u,li) = (m * ucoeff(a,li,li)) % p;
     402      482574 :   for (i = li-1; i > 0; i--)
     403             :   {
     404      356067 :     m = p - uel(b,i)%p;
     405     1196025 :     for (j = i+1; j <= li; j++) {
     406      839958 :       if (m & HIGHBIT) m %= p;
     407      839958 :       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
     408             :     }
     409      356067 :     m %= p;
     410      356067 :     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
     411      356067 :     uel(u,i) = m;
     412             :   }
     413      126507 :   return u;
     414             : }
     415             : static GEN
     416     1369429 : Fl_get_col(GEN a, GEN b, long li, ulong p)
     417             : {
     418     1369429 :   GEN u = cgetg(li+1,t_VECSMALL);
     419     1369429 :   ulong m = uel(b,li) % p;
     420             :   long i,j;
     421             : 
     422     1369429 :   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
     423     4141804 :   for (i=li-1; i>0; i--)
     424             :   {
     425     2772375 :     m = b[i]%p;
     426     7703122 :     for (j = i+1; j <= li; j++)
     427     4930747 :       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
     428     2772375 :     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
     429     2772375 :     uel(u,i) = m;
     430             :   }
     431     1369429 :   return u;
     432             : }
     433             : 
     434             : static GEN
     435      113047 : 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      113047 :   n = lg(x)-1;
     442      113047 :   m=nbrows(x); r=0;
     443             : 
     444      113047 :   c = zero_zv(m);
     445      113047 :   d = cgetg(n+1, t_VECSMALL);
     446      113047 :   a = 0; /* for gcc -Wall */
     447      555654 :   for (k=1; k<=n; k++)
     448             :   {
     449     1930083 :     for (j=1; j<=m; j++)
     450     1737129 :       if (!c[j])
     451             :       {
     452     1033719 :         a = ucoeff(x,j,k) % p;
     453     1033719 :         if (a) break;
     454             :       }
     455      461341 :     if (j > m)
     456             :     {
     457      192954 :       if (deplin==1) {
     458       18734 :         c = cgetg(n+1, t_VECSMALL);
     459       18734 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
     460       18734 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     461       18734 :         return c;
     462             :       }
     463      174220 :       r++; d[k] = 0;
     464             :     }
     465             :     else
     466             :     {
     467      268387 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     468      268387 :       c[j] = k; d[k] = j;
     469      268387 :       ucoeff(x,j,k) = p-1;
     470      268387 :       if (piv != 1)
     471      200225 :         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
     472     1797383 :       for (t=1; t<=m; t++)
     473             :       {
     474     1528996 :         if (t == j) continue;
     475             : 
     476     1260609 :         piv = ( ucoeff(x,t,k) %= p );
     477     1260609 :         if (!piv) continue;
     478      624590 :         if (piv == 1)
     479      149070 :           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
     480             :         else
     481      475520 :           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
     482             :       }
     483             :     }
     484             :   }
     485       94313 :   if (deplin==1) return NULL;
     486             : 
     487       94306 :   y = cgetg(r+1, t_MAT);
     488      268526 :   for (j=k=1; j<=r; j++,k++)
     489             :   {
     490      174220 :     GEN C = cgetg(n+1, t_VECSMALL);
     491             : 
     492      174220 :     gel(y,j) = C; while (d[k]) k++;
     493      632797 :     for (i=1; i<k; i++)
     494      458577 :       if (d[i])
     495      303316 :         uel(C,i) = ucoeff(x,d[i],k) % p;
     496             :       else
     497      155261 :         uel(C,i) = 0UL;
     498      174220 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     499             :   }
     500       94306 :   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       94306 :   return y;
     507             : }
     508             : 
     509             : /* in place, destroy x */
     510             : static GEN
     511      256626 : 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      256626 :   n = lg(x)-1;
     517      256626 :   if (!n) return cgetg(1,t_MAT);
     518      256409 :   if (SMALL_ULONG(p)) return Flm_ker_gauss_OK(x, p, deplin);
     519      143362 :   pi = get_Fl_red(p);
     520             : 
     521      143359 :   m=nbrows(x); r=0;
     522             : 
     523      143359 :   c = zero_zv(m);
     524      143357 :   d = cgetg(n+1, t_VECSMALL);
     525      143357 :   a = 0; /* for gcc -Wall */
     526      459140 :   for (k=1; k<=n; k++)
     527             :   {
     528      870663 :     for (j=1; j<=m; j++)
     529      769951 :       if (!c[j])
     530             :       {
     531      597735 :         a = ucoeff(x,j,k);
     532      597735 :         if (a) break;
     533             :       }
     534      315783 :     if (j > m)
     535             :     {
     536      100708 :       if (deplin==1) {
     537           7 :         c = cgetg(n+1, t_VECSMALL);
     538           7 :         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      100701 :       r++; d[k] = 0;
     543             :     }
     544             :     else
     545             :     {
     546      215075 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     547      215077 :       c[j] = k; d[k] = j;
     548      215077 :       ucoeff(x,j,k) = p-1;
     549      215077 :       if (piv != 1)
     550      391717 :         for (i=k+1; i<=n; i++)
     551      179037 :           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
     552      890161 :       for (t=1; t<=m; t++)
     553             :       {
     554      675079 :         if (t == j) continue;
     555             : 
     556      460004 :         piv = ucoeff(x,t,k);
     557      460004 :         if (!piv) continue;
     558      272929 :         if (piv == 1)
     559       12482 :           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
     560             :         else
     561      260447 :           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
     562             :       }
     563             :     }
     564             :   }
     565      143357 :   if (deplin==1) return NULL;
     566             : 
     567      143350 :   y = cgetg(r+1, t_MAT);
     568      244042 :   for (j=k=1; j<=r; j++,k++)
     569             :   {
     570      100697 :     GEN C = cgetg(n+1, t_VECSMALL);
     571             : 
     572      100701 :     gel(y,j) = C; while (d[k]) k++;
     573      209009 :     for (i=1; i<k; i++)
     574      108308 :       if (d[i])
     575       73842 :         uel(C,i) = ucoeff(x,d[i],k);
     576             :       else
     577       34466 :         uel(C,i) = 0UL;
     578      100701 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     579             :   }
     580      143345 :   if (deplin == 2) {
     581      139547 :     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
     582      431610 :     for (i = j = 1; j <= n; j++)
     583      292063 :       if (d[j]) pc[i++] = j;
     584      139547 :     return mkvec2(y, pc);
     585             :   }
     586        3798 :   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      121845 : Flm_ker_echelon(GEN x, ulong p, long pivots) {
     604      121845 :   pari_sp av = avma;
     605             :   GEN R, Rc, C, C1, C2, S, K;
     606      121845 :   long n = lg(x) - 1, r;
     607      121845 :   ulong pi = get_Fl_red(p);
     608      121845 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     609      121845 :   Rc = indexcompl(R, n);
     610      121845 :   C1 = rowpermute(C, R);
     611      121845 :   C2 = rowpermute(C, Rc);
     612      121845 :   S = Flm_lsolve_lower_unit_pre(C1, C2, p, pi);
     613      121845 :   K = vecpermute(shallowconcat(Flm_neg(S, p), matid_Flm(n - r)),
     614             :                  perm_inv(vecsmall_concat(R, Rc)));
     615      121845 :   K = Flm_transpose(K);
     616      121845 :   if (pivots)
     617       28975 :     return gerepilecopy(av, mkvec2(K, R));
     618       92870 :   return gerepileupto(av, K);
     619             : }
     620             : 
     621             : static GEN
     622       16559 : Flm_deplin_echelon(GEN x, ulong p) {
     623       16559 :   pari_sp av = avma;
     624             :   GEN R, Rc, C, C1, C2, s, v;
     625       16559 :   long i, n = lg(x) - 1, r;
     626       16559 :   ulong pi = get_Fl_red(p);
     627       16559 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     628       16559 :   if (r == n) return gc_NULL(av);
     629       16552 :   Rc = indexcompl(R, n);
     630       16552 :   i = Rc[1];
     631       16552 :   C1 = rowpermute(C, R);
     632       16552 :   C2 = rowslice(C, i, i);
     633       16552 :   s = Flm_row(Flm_lsolve_lower_unit_pre(C1, C2, p, pi), 1);
     634       16552 :   v = vecsmallpermute(vecsmall_concat(Flv_neg(s, p), vecsmall_ei(n - r, 1)),
     635             :                  perm_inv(vecsmall_concat(R, Rc)));
     636       16552 :   return gerepileuptoleaf(av, v);
     637             : }
     638             : 
     639             : static GEN
     640      395029 : Flm_ker_i(GEN x, ulong p, long deplin, long inplace) {
     641      395029 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     642      138404 :     switch(deplin) {
     643       92870 :     case 0: return Flm_ker_echelon(x, p, 0);
     644       16559 :     case 1: return Flm_deplin_echelon(x, p);
     645       28975 :     case 2: return Flm_ker_echelon(x, p, 1);
     646             :     }
     647      256625 :   return Flm_ker_gauss(inplace? x: Flm_copy(x), p, deplin);
     648             : }
     649             : 
     650             : GEN
     651      344434 : Flm_ker_sp(GEN x, ulong p, long deplin) {
     652      344434 :   return Flm_ker_i(x, p, deplin, 1);
     653             : }
     654             : 
     655             : GEN
     656       50595 : Flm_ker(GEN x, ulong p) {
     657       50595 :   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        5705 :     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         784 :       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       50625 : Flm_det_gauss(GEN a, ulong p)
     717             : {
     718       50625 :   long i,j,k, s = 1, nbco = lg(a)-1;
     719       50625 :   ulong pi, q, x = 1;
     720             : 
     721       50625 :   if (SMALL_ULONG(p)) return Flm_det_gauss_OK(a, nbco, p);
     722       48882 :   pi = get_Fl_red(p);
     723      315032 :   for (i=1; i<nbco; i++)
     724             :   {
     725      266553 :     for(k=i; k<=nbco; k++)
     726      266552 :       if (ucoeff(a,k,i)) break;
     727      266150 :     if (k > nbco) return ucoeff(a,i,i);
     728      266150 :     if (k != i)
     729             :     { /* exchange the lines s.t. k = i */
     730         212 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     731         212 :       s = -s;
     732             :     }
     733      266150 :     q = ucoeff(a,i,i);
     734             : 
     735      266150 :     x = Fl_mul_pre(x, q, p, pi);
     736      266151 :     q = Fl_inv(q,p);
     737     1146132 :     for (k=i+1; k<=nbco; k++)
     738             :     {
     739      879980 :       ulong m = ucoeff(a,i,k);
     740      879980 :       if (!m) continue;
     741             : 
     742      839916 :       m = Fl_mul_pre(m, q, p, pi);
     743     4321576 :       for (j=i+1; j<=nbco; j++)
     744     3481662 :         ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
     745             :     }
     746             :   }
     747       48882 :   if (s < 0) x = Fl_neg(x, p);
     748       48882 :   return Fl_mul(x, ucoeff(a,nbco,nbco), p);
     749             : }
     750             : 
     751             : static ulong
     752       35591 : Flm_det_CUP(GEN a, ulong p) {
     753             :   GEN R, C, U, P;
     754       35591 :   long i, n = lg(a) - 1, r;
     755             :   ulong d;
     756       35591 :   ulong pi = get_Fl_red(p);
     757       35591 :   r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi);
     758       35591 :   if (r < n)
     759          42 :     d = 0;
     760             :   else {
     761       35549 :     d = perm_sign(P) == 1? 1: p-1;
     762      376794 :     for (i = 1; i <= n; i++)
     763      341245 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
     764             :   }
     765       35591 :   return d;
     766             : }
     767             : 
     768             : static ulong
     769       86216 : Flm_det_i(GEN x, ulong p, long inplace) {
     770       86216 :   pari_sp av = avma;
     771             :   ulong d;
     772       86216 :   if (lg(x) - 1 >= Flm_CUP_LIMIT)
     773       35591 :     d = Flm_det_CUP(x, p);
     774             :   else
     775       50625 :     d = Flm_det_gauss(inplace? x: Flm_copy(x), p);
     776       86216 :   return gc_ulong(av, d);
     777             : }
     778             : 
     779             : ulong
     780       86216 : 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      206210 : 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      206210 :   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
     792             : 
     793      206098 :   m=nbrows(x); r=0;
     794      206098 :   d=cgetg(n+1,t_VECSMALL);
     795      206098 :   c = zero_zv(m);
     796      932364 :   for (k=1; k<=n; k++)
     797             :   {
     798     2492398 :     for (j=1; j<=m; j++)
     799     2309039 :       if (!c[j])
     800             :       {
     801     1289297 :         ucoeff(x,j,k) %= p;
     802     1289297 :         if (ucoeff(x,j,k)) break;
     803             :       }
     804      726266 :     if (j>m) { r++; d[k]=0; }
     805             :     else
     806             :     {
     807      542907 :       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
     808      542907 :       c[j]=k; d[k]=j;
     809     1615354 :       for (i=k+1; i<=n; i++)
     810     1072447 :         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
     811     4306419 :       for (t=1; t<=m; t++)
     812     3763512 :         if (!c[t]) /* no pivot on that line yet */
     813             :         {
     814     2529363 :           piv = ucoeff(x,t,k);
     815     2529363 :           if (piv)
     816             :           {
     817     1438403 :             ucoeff(x,t,k) = 0;
     818     4457567 :             for (i=k+1; i<=n; i++)
     819     6038328 :               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
     820     3019164 :                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
     821             :           }
     822             :         }
     823      542907 :       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
     824             :     }
     825             :   }
     826      206098 :   *rr = r; set_avma((pari_sp)d); return d;
     827             : }
     828             : 
     829             : static GEN
     830       60167 : Flm_pivots_CUP(GEN x, ulong p, long *rr) {
     831             :   pari_sp av;
     832       60167 :   long i, n = lg(x) - 1, r;
     833       60167 :   GEN R, C, U, P, d = zero_zv(n);
     834       60167 :   ulong pi = get_Fl_red(p);
     835       60167 :   av = avma;
     836       60167 :   r = Flm_CUP_pre(x, &R, &C, &U, &P, p, pi);
     837      716321 :   for(i = 1; i <= r; i++)
     838      656154 :     d[P[i]] = R[i];
     839       60167 :   set_avma(av);
     840       60167 :   *rr = n - r;
     841       60167 :   return d;
     842             : }
     843             : 
     844             : GEN
     845      266377 : Flm_pivots(GEN x, ulong p, long *rr, long inplace) {
     846      266377 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     847       60167 :     return Flm_pivots_CUP(x, p, rr);
     848      206210 :   return Flm_gauss_pivot(inplace? x: Flm_copy(x), p, rr);
     849             : }
     850             : 
     851             : long
     852       14546 : Flm_rank(GEN x, ulong p)
     853             : {
     854       14546 :   pari_sp av = avma;
     855             :   long r;
     856       14546 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT) {
     857             :     GEN R, C;
     858        7336 :     ulong pi = get_Fl_red(p);
     859        7336 :     return gc_long(av, Flm_echelon_pre(x, &R, &C, p, pi));
     860             :   }
     861        7210 :   (void) Flm_pivots(x, p, &r, 0);
     862        7210 :   return gc_long(av, lg(x)-1 - r);
     863             : }
     864             : 
     865             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
     866             :  * reduced mod p */
     867             : static GEN
     868          28 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
     869             : {
     870          28 :   long n = lg(A)-1, i = index, j;
     871          28 :   GEN u = const_vecsmall(n, 0);
     872          28 :   u[i] = 1;
     873          28 :   if (SMALL_ULONG(p))
     874          21 :     for (i--; i>0; i--)
     875             :     {
     876           7 :       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
     877           7 :       for (j=i+2; j<=n; j++)
     878             :       {
     879           0 :         if (m & HIGHMASK) m %= p;
     880           0 :         m += ucoeff(A,i,j) * uel(u,j);
     881             :       }
     882           7 :       u[i] = Fl_neg(m % p, p);
     883             :     }
     884             :   else
     885          21 :     for (i--; i>0; i--)
     886             :     {
     887           7 :       ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
     888           7 :       for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
     889           7 :       u[i] = Fl_neg(m, p);
     890             :     }
     891          28 :   return u;
     892             : }
     893             : static GEN
     894          14 : Flm_inv_upper_1(GEN A, ulong p)
     895             : {
     896             :   long i, l;
     897          14 :   GEN B = cgetg_copy(A, &l);
     898          14 :   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
     899          14 :   return B;
     900             : }
     901             : /* destroy a, b */
     902             : static GEN
     903       34245 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
     904             : {
     905       34245 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     906       34245 :   ulong det = 1;
     907             :   GEN u;
     908             : 
     909       34245 :   li = nbrows(a);
     910       34245 :   bco = lg(b)-1;
     911      113093 :   for (i=1; i<=aco; i++)
     912             :   {
     913             :     ulong invpiv;
     914             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
     915      113093 :     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
     916      128402 :     for (k = i; k <= li; k++)
     917             :     {
     918      128395 :       ulong piv = ( ucoeff(a,k,i) %= p );
     919      128395 :       if (piv)
     920             :       {
     921      113086 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     922      113086 :         if (detp)
     923             :         {
     924           0 :           if (det & HIGHMASK) det %= p;
     925           0 :           det *= piv;
     926             :         }
     927      113086 :         break;
     928             :       }
     929             :     }
     930             :     /* found a pivot on line k */
     931      113093 :     if (k > li) return NULL;
     932      113086 :     if (k != i)
     933             :     { /* swap lines so that k = i */
     934       10046 :       s = -s;
     935       10046 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     936       10046 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     937             :     }
     938      113086 :     if (i == aco) break;
     939             : 
     940       78848 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
     941      267843 :     for (k=i+1; k<=li; k++)
     942             :     {
     943      188995 :       ulong m = ( ucoeff(a,k,i) %= p );
     944      188995 :       if (!m) continue;
     945             : 
     946       59154 :       m = Fl_mul(m, invpiv, p);
     947       59154 :       if (m == 1) {
     948       12635 :         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
     949       12635 :         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
     950             :       } else {
     951       46519 :         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
     952       46519 :         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
     953             :       }
     954             :     }
     955             :   }
     956       34238 :   if (detp)
     957             :   {
     958           0 :     det %=  p;
     959           0 :     if (s < 0 && det) det = p - det;
     960           0 :     *detp = det;
     961             :   }
     962       34238 :   u = cgetg(bco+1,t_MAT);
     963       34238 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
     964       34238 :   return u;
     965             : }
     966             : 
     967             : /* destroy a, b */
     968             : static GEN
     969      535609 : Flm_gauss_sp_i(GEN a, GEN b, ulong *detp, ulong p)
     970             : {
     971      535609 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     972      535609 :   ulong det = 1;
     973             :   GEN u;
     974             :   ulong pi;
     975      535609 :   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
     976      535609 :   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
     977      501364 :   pi = get_Fl_red(p);
     978      501364 :   li = nbrows(a);
     979      501364 :   bco = lg(b)-1;
     980     1369424 :   for (i=1; i<=aco; i++)
     981             :   {
     982             :     ulong invpiv;
     983             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
     984     1430465 :     for (k = i; k <= li; k++)
     985             :     {
     986     1430458 :       ulong piv = ucoeff(a,k,i);
     987     1430458 :       if (piv)
     988             :       {
     989     1369417 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     990     1369417 :         if (detp) det = Fl_mul_pre(det, piv, p, pi);
     991     1369417 :         break;
     992             :       }
     993             :     }
     994             :     /* found a pivot on line k */
     995     1369424 :     if (k > li) return NULL;
     996     1369417 :     if (k != i)
     997             :     { /* swap lines so that k = i */
     998       53543 :       s = -s;
     999       53543 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
    1000       53543 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
    1001             :     }
    1002     1369417 :     if (i == aco) break;
    1003             : 
    1004      868060 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    1005     2254339 :     for (k=i+1; k<=li; k++)
    1006             :     {
    1007     1386279 :       ulong m = ucoeff(a,k,i);
    1008     1386279 :       if (!m) continue;
    1009             : 
    1010     1169240 :       m = Fl_mul_pre(m, invpiv, p, pi);
    1011     1169240 :       if (m == 1) {
    1012       36163 :         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
    1013       36163 :         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
    1014             :       } else {
    1015     1133077 :         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
    1016     1133077 :         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
    1017             :       }
    1018             :     }
    1019             :   }
    1020      501357 :   if (detp)
    1021             :   {
    1022           0 :     if (s < 0 && det) det = p - det;
    1023           0 :     *detp = det;
    1024             :   }
    1025      501357 :   u = cgetg(bco+1,t_MAT);
    1026      501357 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
    1027      501357 :   return u;
    1028             : }
    1029             : 
    1030             : static GEN
    1031      212104 : Flm_gauss_from_CUP(GEN b, GEN R, GEN C, GEN U, GEN P, ulong p, ulong pi, ulong *detp)
    1032             : {
    1033      212104 :   GEN Y = Flm_rsolve_lower_unit_pre(rowpermute(C, R), rowpermute(b, R), p, pi);
    1034      212105 :   GEN X = rowpermute(Flm_rsolve_upper_pre(U, Y, p, pi), perm_inv(P));
    1035      212105 :   if (detp)
    1036             :   {
    1037      160452 :     ulong d = perm_sign(P) == 1? 1: p-1;
    1038      160449 :     long i, r = lg(R);
    1039     1837287 :     for (i = 1; i < r; i++)
    1040     1676835 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
    1041      160452 :     *detp = d;
    1042             :   }
    1043      212105 :   return X;
    1044             : }
    1045             : 
    1046             : static GEN
    1047       51667 : Flm_gauss_CUP(GEN a, GEN b, ulong *detp, ulong p) {
    1048             :   GEN R, C, U, P;
    1049       51667 :   long n = lg(a) - 1, r;
    1050       51667 :   ulong pi = get_Fl_red(p);
    1051       51667 :   if (nbrows(a) < n || (r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi)) < n)
    1052          14 :     return NULL;
    1053       51653 :   return Flm_gauss_from_CUP(b, R, C, U, P, p, pi, detp);
    1054             : }
    1055             : 
    1056             : GEN
    1057      551773 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p) {
    1058      551773 :   pari_sp av = avma;
    1059             :   GEN x;
    1060      551773 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1061       16164 :     x = Flm_gauss_CUP(a, b, detp, p);
    1062             :   else
    1063      535609 :     x = Flm_gauss_sp_i(a, b, detp, p);
    1064      551773 :   if (!x) return gc_NULL(av);
    1065      551759 :   return gerepileupto(av, x);
    1066             : }
    1067             : 
    1068             : GEN
    1069          63 : Flm_gauss(GEN a, GEN b, ulong p) {
    1070          63 :   pari_sp av = avma;
    1071             :   GEN x;
    1072          63 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1073          14 :     x = Flm_gauss_CUP(a, b, NULL, p);
    1074             :   else {
    1075          49 :     a = RgM_shallowcopy(a);
    1076          49 :     b = RgM_shallowcopy(b);
    1077          49 :     x = Flm_gauss_sp(a, b, NULL, p);
    1078             :   }
    1079          63 :   if (!x) return gc_NULL(av);
    1080          49 :   return gerepileupto(av, x);
    1081             : }
    1082             : 
    1083             : static GEN
    1084      547725 : Flm_inv_i(GEN a, ulong *detp, ulong p, long inplace) {
    1085      547725 :   pari_sp av = avma;
    1086      547725 :   long n = lg(a) - 1;
    1087             :   GEN b, x;
    1088      547725 :   if (n == 0) return cgetg(1, t_MAT);
    1089      547725 :   b = matid_Flm(nbrows(a));
    1090      547725 :   if (n >= Flm_CUP_LIMIT)
    1091       35489 :     x = Flm_gauss_CUP(a, b, detp, p);
    1092             :   else {
    1093      512236 :     if (!inplace)
    1094       10899 :       a = RgM_shallowcopy(a);
    1095      512236 :     x = Flm_gauss_sp(a, b, detp, p);
    1096             :   }
    1097      547725 :   if (!x) return gc_NULL(av);
    1098      547711 :   return gerepileupto(av, x);
    1099             : }
    1100             : 
    1101             : GEN
    1102      525402 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
    1103      525402 :   return Flm_inv_i(a, detp, p, 1);
    1104             : }
    1105             : 
    1106             : GEN
    1107       22323 : Flm_inv(GEN a, ulong p) {
    1108       22323 :   return Flm_inv_i(a, NULL, p, 0);
    1109             : }
    1110             : 
    1111             : GEN
    1112          21 : Flm_Flc_gauss(GEN a, GEN b, ulong p) {
    1113          21 :   pari_sp av = avma;
    1114          21 :   GEN z = Flm_gauss(a, mkmat(b), p);
    1115          21 :   if (!z) return gc_NULL(av);
    1116          14 :   if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
    1117          14 :   return gerepileuptoleaf(av, gel(z,1));
    1118             : }
    1119             : 
    1120             : GEN
    1121      160465 : Flm_adjoint(GEN A, ulong p)
    1122             : {
    1123      160465 :   pari_sp av = avma;
    1124             :   GEN R, C, U, P, C1, U1, v, c, d;
    1125      160465 :   long r, i, q, n = lg(A)-1, m;
    1126             :   ulong D;
    1127      160465 :   ulong pi = get_Fl_red(p);
    1128      160466 :   if (n == 0) return cgetg(1, t_MAT);
    1129      160466 :   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
    1130      160466 :   m = nbrows(A);
    1131      160466 :   if (r == n)
    1132             :   {
    1133      160452 :     GEN X = Flm_gauss_from_CUP(matid_Flm(m), R, C, U, P, p, pi, &D);
    1134      160450 :     return gerepileupto(av, Flm_Fl_mul_pre(X, D, p, pi));
    1135             :   }
    1136          14 :   if (r < n-1) return zero_Flm(n, m);
    1137          28 :   for (q = n, i = 1; i < n; i++)
    1138          14 :     if (R[i] != i) { q = i; break; }
    1139          14 :   C1 = matslice(C, 1, q-1, 1, q-1);
    1140          14 :   c = vecslice(Flm_row(C, q), 1, q-1);
    1141          14 :   c = Flm_lsolve_lower_unit_pre(C1, Flm_transpose(mkmat(c)), p, pi);
    1142          14 :   d = cgetg(m+1, t_VECSMALL);
    1143          14 :   for (i=1; i<q; i++)    uel(d,i) = ucoeff(c,1,i);
    1144          14 :   uel(d,q) = p-1;
    1145          14 :   for (i=q+1; i<=m; i++) uel(d,i) = 0;
    1146          14 :   U1 = vecslice(U, 1, n-1);
    1147          14 :   v = gel(Flm_rsolve_upper_pre(U1, mkmat(gel(U,n)), p, pi),1);
    1148          14 :   v = vecsmall_append(v, p-1);
    1149          14 :   D = perm_sign(P) != (odd(q+n)?-1:1) ? p-1 : 1;
    1150          28 :   for (i = 1; i <= n-1; i++)
    1151          14 :     D = Fl_mul_pre(D, ucoeff(U1, i, i), p, pi);
    1152          14 :   d = Flv_Fl_mul(d, D, p);
    1153          14 :   return rowpermute(Flc_Flv_mul(v, d, p),perm_inv(P));
    1154             : }
    1155             : 
    1156             : static GEN
    1157          21 : Flm_invimage_CUP(GEN A, GEN B, ulong p) {
    1158          21 :   pari_sp av = avma;
    1159             :   GEN R, Rc, C, U, P, B1, B2, C1, C2, X, Y, Z;
    1160             :   long r;
    1161          21 :   ulong pi = get_Fl_red(p);
    1162          21 :   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
    1163          21 :   Rc = indexcompl(R, nbrows(B));
    1164          21 :   C1 = rowpermute(C, R);
    1165          21 :   C2 = rowpermute(C, Rc);
    1166          21 :   B1 = rowpermute(B, R);
    1167          21 :   B2 = rowpermute(B, Rc);
    1168          21 :   Z = Flm_rsolve_lower_unit_pre(C1, B1, p, pi);
    1169          21 :   if (!gequal(Flm_mul_pre(C2, Z, p, pi), B2))
    1170          14 :     return NULL;
    1171          14 :   Y = vconcat(Flm_rsolve_upper_pre(vecslice(U, 1, r), Z, p, pi),
    1172          14 :               zero_Flm(lg(A) - 1 - r, lg(B) - 1));
    1173           7 :   X = rowpermute(Y, perm_inv(P));
    1174           7 :   return gerepileupto(av, X);
    1175             : }
    1176             : 
    1177             : GEN
    1178          42 : Flm_invimage_i(GEN A, GEN B, ulong p)
    1179             : {
    1180             :   GEN d, x, X, Y;
    1181          42 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    1182             : 
    1183          42 :   if (!nB) return cgetg(1, t_MAT);
    1184          42 :   if (nA + nB >= Flm_CUP_LIMIT && nbrows(B) >= Flm_CUP_LIMIT)
    1185          21 :     return Flm_invimage_CUP(A, B, p);
    1186             : 
    1187          21 :   x = Flm_ker_sp(shallowconcat(Flm_neg(A,p), B), p, 0);
    1188             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    1189             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    1190             :    * Y has at least nB columns and full rank */
    1191          21 :   nY = lg(x)-1;
    1192          21 :   if (nY < nB) return NULL;
    1193          21 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    1194          21 :   d = cgetg(nB+1, t_VECSMALL);
    1195          56 :   for (i = nB, j = nY; i >= 1; i--, j--)
    1196             :   {
    1197          49 :     for (; j>=1; j--)
    1198          42 :       if (coeff(Y,i,j)) { d[i] = j; break; }
    1199          42 :     if (!j) return NULL;
    1200             :   }
    1201             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    1202          14 :   Y = vecpermute(Y, d);
    1203          14 :   x = vecpermute(x, d);
    1204          14 :   X = rowslice(x, 1, nA);
    1205          14 :   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
    1206             : }
    1207             : GEN
    1208           0 : Flm_invimage(GEN A, GEN B, ulong p)
    1209             : {
    1210           0 :   pari_sp av = avma;
    1211           0 :   GEN X = Flm_invimage_i(A,B,p);
    1212           0 :   if (!X) return gc_NULL(av);
    1213           0 :   return gerepileupto(av, X);
    1214             : }
    1215             : 
    1216             : GEN
    1217       32042 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
    1218             : {
    1219       32042 :   pari_sp av = avma;
    1220       32042 :   long i, l = lg(A);
    1221             :   GEN M, x;
    1222             :   ulong t;
    1223             : 
    1224       32042 :   if (l==1) return NULL;
    1225       32042 :   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
    1226       32042 :   M = cgetg(l+1,t_MAT);
    1227       32042 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    1228       32042 :   gel(M,l) = y; M = Flm_ker(M,p);
    1229       32042 :   i = lg(M)-1; if (!i) return gc_NULL(av);
    1230             : 
    1231       32042 :   x = gel(M,i); t = x[l];
    1232       32042 :   if (!t) return gc_NULL(av);
    1233             : 
    1234       32035 :   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
    1235       32035 :   if (t!=1) x = Flv_Fl_mul(x, t, p);
    1236       32035 :   return gerepileuptoleaf(av, x);
    1237             : }

Generated by: LCOV version 1.13