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.14.0 lcov report (development 27775-aca467eab2) Lines: 737 758 97.2 %
Date: 2022-07-03 07:33:15 Functions: 58 62 93.5 %
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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : GEN
      19      301880 : Flv_to_ZV(GEN x)
      20     3603616 : { pari_APPLY_type(t_VEC, utoi(x[i])) }
      21             : 
      22             : GEN
      23    31203401 : Flc_to_ZC(GEN x)
      24   229003487 : { pari_APPLY_type(t_COL, utoi(x[i])) }
      25             : 
      26             : GEN
      27    11096179 : Flm_to_ZM(GEN x)
      28    41385227 : { pari_APPLY_type(t_MAT, Flc_to_ZC(gel(x,i))) }
      29             : 
      30             : GEN
      31          14 : Flc_to_ZC_inplace(GEN z)
      32             : {
      33          14 :   long i, l = lg(z);
      34         252 :   for (i=1; i<l; i++) gel(z,i) = utoi(z[i]);
      35          14 :   settyp(z, t_COL);
      36          14 :   return z;
      37             : }
      38             : 
      39             : GEN
      40           0 : Flm_to_ZM_inplace(GEN z)
      41             : {
      42           0 :   long i, l = lg(z);
      43           0 :   for (i=1; i<l; i++) Flc_to_ZC_inplace(gel(z, i));
      44           0 :   return z;
      45             : }
      46             : 
      47             : static GEN
      48     1463912 : Flm_solve_upper_1(GEN U, GEN B, ulong p, ulong pi)
      49     1463912 : { return Flm_Fl_mul_pre(B, Fl_inv(ucoeff(U, 1, 1), p), p, pi); }
      50             : 
      51             : static GEN
      52     2023456 : Flm_rsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      53             : {
      54     2023456 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      55     2023456 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      56     2023583 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      57     2023570 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
      58     2023558 :   GEN B1 = rowslice(B, 1, 1);
      59     2023405 :   GEN B2 = rowslice(B, 2, 2);
      60     2023401 :   GEN X2 = Flm_Fl_mul_pre(B2, dinv, p, pi);
      61     2023511 :   GEN X1 = Flm_Fl_mul_pre(Flm_sub(B1, Flm_Fl_mul_pre(X2, b, p, pi), p),
      62             :                       ainv, p, pi);
      63     2023525 :   return vconcat(X1, X2);
      64             : }
      65             : 
      66             : /* solve U*X = B,  U upper triangular and invertible */
      67             : static GEN
      68     5039779 : Flm_rsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
      69             : {
      70     5039779 :   long n = lg(U) - 1, n1;
      71             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
      72     5039779 :   pari_sp av = avma;
      73             : 
      74     5039779 :   if (n == 0) return B;
      75     5039772 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
      76     4160154 :   if (n == 2) return Flm_rsolve_upper_2(U, B, p, pi);
      77     2136730 :   n1 = (n + 1)/2;
      78     2136730 :   U2 = vecslice(U, n1 + 1, n);
      79     2137494 :   U22 = rowslice(U2, n1 + 1, n);
      80     2137397 :   B2 = rowslice(B, n1 + 1, n);
      81     2137390 :   X2 = Flm_rsolve_upper_pre(U22, B2, p, pi);
      82     2137693 :   U12 = rowslice(U2, 1, n1);
      83     2137466 :   B1 = rowslice(B, 1, n1);
      84     2137461 :   B1 = Flm_sub(B1, Flm_mul_pre(U12, X2, p, pi), p);
      85     2137427 :   if (gc_needed(av, 1)) gerepileall(av, 2, &B1, &X2);
      86     2137427 :   U11 = matslice(U, 1,n1, 1,n1);
      87     2137372 :   X1 = Flm_rsolve_upper_pre(U11, B1, p, pi);
      88     2137741 :   X = vconcat(X1, X2);
      89     2137770 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
      90     2137774 :   return X;
      91             : }
      92             : 
      93             : static GEN
      94     1913852 : Flm_lsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      95             : {
      96     1913852 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      97     1913852 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      98     1913881 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      99     1913872 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
     100     1913872 :   GEN B1 = vecslice(B, 1, 1);
     101     1913838 :   GEN B2 = vecslice(B, 2, 2);
     102     1913793 :   GEN X1 = Flm_Fl_mul_pre(B1, ainv, p, pi);
     103     1913855 :   GEN X2 = Flm_Fl_mul_pre(Flm_sub(B2, Flm_Fl_mul_pre(X1, b, p, pi), p),
     104             :                       dinv, p, pi);
     105     1913852 :   return shallowconcat(X1, X2);
     106             : }
     107             : 
     108             : /* solve X*U = B,  U upper triangular and invertible */
     109             : static GEN
     110     4243121 : Flm_lsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
     111             : {
     112     4243121 :   long n = lg(U) - 1, n1;
     113             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
     114     4243121 :   pari_sp av = avma;
     115             : 
     116     4243121 :   if (n == 0) return B;
     117     4243121 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
     118     3658812 :   if (n == 2) return Flm_lsolve_upper_2(U, B, p, pi);
     119     1744966 :   n1 = (n + 1)/2;
     120     1744966 :   U2 = vecslice(U, n1 + 1, n);
     121     1745255 :   U11 = matslice(U, 1,n1, 1,n1);
     122     1745205 :   U12 = rowslice(U2, 1, n1);
     123     1745202 :   U22 = rowslice(U2, n1 + 1, n);
     124     1745192 :   B1 = vecslice(B, 1, n1);
     125     1745185 :   B2 = vecslice(B, n1 + 1, n);
     126     1745175 :   X1 = Flm_lsolve_upper_pre(U11, B1, p, pi);
     127     1745255 :   B2 = Flm_sub(B2, Flm_mul_pre(X1, U12, p, pi), p);
     128     1745229 :   if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
     129     1745229 :   X2 = Flm_lsolve_upper_pre(U22, B2, p, pi);
     130     1745261 :   X = shallowconcat(X1, X2);
     131     1745270 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     132     1745268 :   return X;
     133             : }
     134             : 
     135             : static GEN
     136     4477934 : Flm_rsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     137             : {
     138     4477934 :   GEN X1 = rowslice(A, 1, 1);
     139     4477915 :   GEN X2 = Flm_sub(rowslice(A, 2, 2),
     140     4477854 :                    Flm_Fl_mul_pre(X1, ucoeff(L, 2, 1), p, pi), p);
     141     4477745 :   return vconcat(X1, X2);
     142             : }
     143             : 
     144             : /* solve L*X = A,  L lower triangular with ones on the diagonal
     145             : * (at least as many rows as columns) */
     146             : static GEN
     147    10153652 : Flm_rsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     148             : {
     149    10153652 :   long m = lg(L) - 1, m1, n;
     150             :   GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
     151    10153652 :   pari_sp av = avma;
     152             : 
     153    10153652 :   if (m == 0) return zero_Flm(0, lg(A) - 1);
     154    10153652 :   if (m == 1) return rowslice(A, 1, 1);
     155     8718339 :   if (m == 2) return Flm_rsolve_lower_unit_2(L, A, p, pi);
     156     4240454 :   m1 = (m + 1)/2;
     157     4240454 :   n = nbrows(L);
     158     4241323 :   L1 = vecslice(L, 1, m1);
     159     4241261 :   L11 = rowslice(L1, 1, m1);
     160     4241162 :   L21 = rowslice(L1, m1 + 1, n);
     161     4241137 :   A1 = rowslice(A, 1, m1);
     162     4241206 :   X1 = Flm_rsolve_lower_unit_pre(L11, A1, p, pi);
     163     4241502 :   A2 = rowslice(A, m1 + 1, n);
     164     4241114 :   A2 = Flm_sub(A2, Flm_mul_pre(L21, X1, p, pi), p);
     165     4241132 :   if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
     166     4241132 :   L22 = matslice(L, m1+1,n, m1+1,m);
     167     4241110 :   X2 = Flm_rsolve_lower_unit_pre(L22, A2, p, pi);
     168     4241311 :   X = vconcat(X1, X2);
     169     4241537 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     170     4241537 :   return X;
     171             : }
     172             : 
     173             : static GEN
     174     2560500 : Flm_lsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     175             : {
     176     2560500 :   GEN X2 = vecslice(A, 2, 2);
     177     2560501 :   GEN X1 = Flm_sub(vecslice(A, 1, 1),
     178     2560500 :                    Flm_Fl_mul_pre(X2, ucoeff(L, 2, 1), p, pi), p);
     179     2560500 :   return shallowconcat(X1, X2);
     180             : }
     181             : 
     182             : /* solve L*X = A,  L square lower triangular with ones on the diagonal */
     183             : static GEN
     184     5945587 : Flm_lsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     185             : {
     186     5945587 :   long m = lg(L) - 1, m1;
     187             :   GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
     188     5945587 :   pari_sp av = avma;
     189             : 
     190     5945587 :   if (m <= 1) return A;
     191     5299293 :   if (m == 2) return Flm_lsolve_lower_unit_2(L, A, p, pi);
     192     2738793 :   m1 = (m + 1)/2;
     193     2738793 :   L2 = vecslice(L, m1 + 1, m);
     194     2738799 :   L22 = rowslice(L2, m1 + 1, m);
     195     2738802 :   A2 = vecslice(A, m1 + 1, m);
     196     2738802 :   X2 = Flm_lsolve_lower_unit_pre(L22, A2, p, pi);
     197     2738801 :   if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
     198     2738801 :   L1 = vecslice(L, 1, m1);
     199     2738802 :   L21 = rowslice(L1, m1 + 1, m);
     200     2738799 :   A1 = vecslice(A, 1, m1);
     201     2738798 :   A1 = Flm_sub(A1, Flm_mul_pre(X2, L21, p, pi), p);
     202     2738800 :   L11 = rowslice(L1, 1, m1);
     203     2738802 :   if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
     204     2738802 :   X1 = Flm_lsolve_lower_unit_pre(L11, A1, p, pi);
     205     2738804 :   X = shallowconcat(X1, X2);
     206     2738804 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     207     2738804 :   return X;
     208             : }
     209             : 
     210             : /* destroy A */
     211             : static long
     212     3180029 : Flm_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     213             : {
     214     3180029 :   long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
     215             :   ulong u, v;
     216             : 
     217     3180033 :   if (P) *P = identity_perm(n);
     218     3179999 :   *R = cgetg(m + 1, t_VECSMALL);
     219    16010126 :   for (j = 1, pr = 0; j <= n; j++)
     220             :   {
     221    18705406 :     for (pr++, pc = 0; pr <= m; pr++)
     222             :     {
     223   125633711 :       for (k = j; k <= n; k++) { v = ucoeff(A, pr, k); if (!pc && v) pc = k; }
     224    17573519 :       if (pc) break;
     225             :     }
     226    13962233 :     if (!pc) break;
     227    12830130 :     (*R)[j] = pr;
     228    12830130 :     if (pc != j)
     229             :     {
     230     1893316 :       swap(gel(A, j), gel(A, pc));
     231     1893316 :       if (P) lswap((*P)[j], (*P)[pc]);
     232             :     }
     233    12830130 :     u = Fl_inv(ucoeff(A, pr, j), p);
     234   103161448 :     for (i = pr + 1; i <= m; i++)
     235             :     {
     236    90331312 :       v = Fl_mul_pre(ucoeff(A, i, j), u, p, pi);
     237    90330995 :       ucoeff(A, i, j) = v;
     238    90330995 :       v = Fl_neg(v, p);
     239   391700941 :       for (k = j + 1; k <= n; k++)
     240   301377312 :         ucoeff(A, i, k) = Fl_addmul_pre(ucoeff(A, i, k),
     241   301369675 :                                         ucoeff(A, pr, k), v, p, pi);
     242             :     }
     243             :   }
     244     3179996 :   setlg(*R, j);
     245     3180035 :   *C = vecslice(A, 1, j - 1);
     246     3179998 :   if (U) *U = rowpermute(A, *R);
     247     3180012 :   return j - 1;
     248             : }
     249             : 
     250             : static const long Flm_CUP_LIMIT = 8;
     251             : 
     252             : static long
     253     2511854 : Flm_CUP_pre(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     254             : {
     255     2511854 :   long m = nbrows(A), m1, n = lg(A) - 1, i, r1, r2, r;
     256             :   GEN R1, C1, U1, P1, R2, C2, U2, P2;
     257             :   GEN A1, A2, B2, C21, U11, U12, T21, T22;
     258     2511856 :   pari_sp av = avma;
     259             : 
     260     2511856 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     261             :     /* destroy A; not called at the outermost recursion level */
     262     1740360 :     return Flm_CUP_basecase(A, R, C, U, P, p, pi);
     263      771496 :   m1 = (minss(m, n) + 1)/2;
     264      771495 :   A1 = rowslice(A, 1, m1);
     265      771489 :   A2 = rowslice(A, m1 + 1, m);
     266      771482 :   r1 = Flm_CUP_pre(A1, &R1, &C1, &U1, &P1, p, pi);
     267      771494 :   if (r1 == 0)
     268             :   {
     269       18566 :     r2 = Flm_CUP_pre(A2, &R2, &C2, &U2, &P2, p, pi);
     270       18566 :     *R = cgetg(r2 + 1, t_VECSMALL);
     271       26784 :     for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
     272       18566 :     *C = vconcat(zero_Flm(m1, r2), C2);
     273       18566 :     *U = U2;
     274       18566 :     *P = P2;
     275       18566 :     r = r2;
     276             :   }
     277             :   else
     278             :   {
     279      752928 :     U11 = vecslice(U1, 1, r1);
     280      752921 :     U12 = vecslice(U1, r1 + 1, n);
     281      752915 :     T21 = vecslicepermute(A2, P1, 1, r1);
     282      752918 :     T22 = vecslicepermute(A2, P1, r1 + 1, n);
     283      752912 :     C21 = Flm_lsolve_upper_pre(U11, T21, p, pi);
     284      752923 :     if (gc_needed(av, 1))
     285           0 :       gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
     286      752923 :     B2 = Flm_sub(T22, Flm_mul_pre(C21, U12, p, pi), p);
     287      752920 :     r2 = Flm_CUP_pre(B2, &R2, &C2, &U2, &P2, p, pi);
     288      752923 :     r = r1 + r2;
     289      752923 :     *R = cgetg(r + 1, t_VECSMALL);
     290     5165011 :     for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
     291     4659854 :     for (;      i <= r; i++)  (*R)[i] = R2[i - r1] + m1;
     292      752913 :     *C = shallowconcat(vconcat(C1, C21),
     293             :                        vconcat(zero_Flm(m1, r2), C2));
     294      752937 :     *U = shallowconcat(vconcat(U11, zero_Flm(r2, r1)),
     295             :                        vconcat(vecpermute(U12, P2), U2));
     296      752934 :     *P = cgetg(n + 1, t_VECSMALL);
     297     5165099 :     for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
     298    11107186 :     for (; i <= n; i++)       (*P)[i] = P1[P2[i - r1] + r1];
     299             :   }
     300      771493 :   if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
     301      771500 :   return r;
     302             : }
     303             : 
     304             : static long
     305     1439635 : Flm_echelon_gauss(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     306     1439635 : { return Flm_CUP_basecase(A, R, C, NULL, NULL, p, pi); }
     307             : 
     308             : /* complement of a strictly increasing subsequence of (1, 2, ..., n) */
     309             : static GEN
     310     1374645 : indexcompl(GEN v, long n)
     311             : {
     312     1374645 :   long i, j, k, m = lg(v) - 1;
     313     1374645 :   GEN w = cgetg(n - m + 1, t_VECSMALL);
     314    28057688 :   for (i = j = k = 1; i <= n; i++)
     315    26683043 :     if (j <= m && v[j] == i) j++; else w[k++] = i;
     316     1374645 :   return w;
     317             : }
     318             : 
     319             : /* column echelon form */
     320             : static long
     321     2406444 : Flm_echelon_pre(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     322             : {
     323     2406444 :   long j, j1, j2, m = nbrows(A), n = lg(A) - 1, n1, r, r1, r2;
     324             :   GEN A1, A2, R1, R1c, C1, R2, C2;
     325             :   GEN A12, A22, B2, C11, C21, M12;
     326     2406444 :   pari_sp av = avma;
     327             : 
     328     2406444 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     329     1439634 :     return Flm_echelon_gauss(Flm_copy(A), R, C, p, pi);
     330             : 
     331      966810 :   n1 = (n + 1)/2;
     332      966810 :   A1 = vecslice(A, 1, n1);
     333      966810 :   A2 = vecslice(A, n1 + 1, n);
     334      966810 :   r1 = Flm_echelon_pre(A1, &R1, &C1, p, pi);
     335      966810 :   if (!r1) return Flm_echelon_pre(A2, R, C, p, pi);
     336      911025 :   if (r1 == m) { *R = R1; *C = C1; return r1; }
     337             : 
     338      906380 :   R1c = indexcompl(R1, m);
     339      906380 :   C11 = rowpermute(C1, R1);
     340      906380 :   C21 = rowpermute(C1, R1c);
     341      906380 :   A12 = rowpermute(A2, R1);
     342      906380 :   A22 = rowpermute(A2, R1c);
     343      906379 :   M12 = Flm_rsolve_lower_unit_pre(C11, A12, p, pi);
     344      906380 :   B2 = Flm_sub(A22, Flm_mul_pre(C21, M12, p, pi), p);
     345      906379 :   r2 = Flm_echelon_pre(B2, &R2, &C2, p, pi);
     346      906380 :   if (!r2) { *R = R1; *C = C1; r = r1; }
     347             :   else
     348             :   {
     349      865672 :     R2 = perm_mul(R1c, R2);
     350      865672 :     C2 = rowpermute(vconcat(zero_Flm(r1, r2), C2),
     351             :                     perm_inv(vecsmall_concat(R1, R1c)));
     352      865672 :     r = r1 + r2;
     353      865672 :     *R = cgetg(r + 1, t_VECSMALL);
     354      865672 :     *C = cgetg(r + 1, t_MAT);
     355    10599773 :     for (j = j1 = j2 = 1; j <= r; j++)
     356     9734101 :       if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
     357             :       {
     358     5090115 :         gel(*C, j) = gel(C1, j1);
     359     5090115 :         (*R)[j] = R1[j1++];
     360             :       }
     361             :       else
     362             :       {
     363     4643986 :         gel(*C, j) = gel(C2, j2);
     364     4643986 :         (*R)[j] = R2[j2++];
     365             :       }
     366             :   }
     367      906380 :   if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
     368      906380 :   return r;
     369             : }
     370             : 
     371             : static void /* assume m < p */
     372    13004170 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
     373             : {
     374    13004170 :   uel(b,k) = Fl_addmul_pre(uel(b, k), m, uel(b, i), p, pi);
     375    13004190 : }
     376             : static void /* same m = 1 */
     377      548132 : _Fl_add(GEN b, long k, long i, ulong p)
     378             : {
     379      548132 :   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
     380      548132 : }
     381             : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     382     4249357 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
     383             : {
     384     4249357 :   uel(b,k) += m * uel(b,i);
     385     4249357 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     386     4249357 : }
     387             : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     388     2307056 : _Fl_add_OK(GEN b, long k, long i, ulong p)
     389             : {
     390     2307056 :   uel(b,k) += uel(b,i);
     391     2307056 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     392     2307056 : }
     393             : 
     394             : /* assume 0 <= a[i,j] < p */
     395             : static GEN
     396      372232 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
     397             : {
     398      372232 :   GEN u = cgetg(li+1,t_VECSMALL);
     399      372230 :   ulong m = uel(b,li) % p;
     400             :   long i,j;
     401             : 
     402      372230 :   uel(u,li) = (m * ucoeff(a,li,li)) % p;
     403     1458809 :   for (i = li-1; i > 0; i--)
     404             :   {
     405     1086579 :     m = p - uel(b,i)%p;
     406     3460771 :     for (j = i+1; j <= li; j++) {
     407     2374192 :       if (m & HIGHBIT) m %= p;
     408     2374192 :       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
     409             :     }
     410     1086579 :     m %= p;
     411     1086579 :     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
     412     1086579 :     uel(u,i) = m;
     413             :   }
     414      372230 :   return u;
     415             : }
     416             : static GEN
     417     2613154 : Fl_get_col(GEN a, GEN b, long li, ulong p)
     418             : {
     419     2613154 :   GEN u = cgetg(li+1,t_VECSMALL);
     420     2613144 :   ulong m = uel(b,li) % p;
     421             :   long i,j;
     422             : 
     423     2613144 :   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
     424     8373248 :   for (i=li-1; i>0; i--)
     425             :   {
     426     5760102 :     m = b[i]%p;
     427    16307175 :     for (j = i+1; j <= li; j++)
     428    10547089 :       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
     429     5760086 :     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
     430     5760099 :     uel(u,i) = m;
     431             :   }
     432     2613146 :   return u;
     433             : }
     434             : 
     435             : static GEN
     436      622232 : Flm_ker_gauss_OK(GEN x, ulong p, long deplin)
     437             : {
     438             :   GEN y, c, d;
     439             :   long i, j, k, r, t, m, n;
     440             :   ulong a;
     441             : 
     442      622232 :   n = lg(x)-1;
     443      622232 :   m=nbrows(x); r=0;
     444             : 
     445      622231 :   c = zero_zv(m);
     446      622231 :   d = cgetg(n+1, t_VECSMALL);
     447      622232 :   a = 0; /* for gcc -Wall */
     448     3129725 :   for (k=1; k<=n; k++)
     449             :   {
     450     9248531 :     for (j=1; j<=m; j++)
     451     8074647 :       if (!c[j])
     452             :       {
     453     4211707 :         a = ucoeff(x,j,k) % p;
     454     4211707 :         if (a) break;
     455             :       }
     456     2564040 :     if (j > m)
     457             :     {
     458     1173890 :       if (deplin==1) {
     459       56546 :         c = cgetg(n+1, t_VECSMALL);
     460      191987 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
     461      114663 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     462       56546 :         return c;
     463             :       }
     464     1117344 :       r++; d[k] = 0;
     465             :     }
     466             :     else
     467             :     {
     468     1390150 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     469     1390152 :       c[j] = k; d[k] = j;
     470     1390152 :       ucoeff(x,j,k) = p-1;
     471     1390152 :       if (piv != 1)
     472     3070977 :         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
     473     7164568 :       for (t=1; t<=m; t++)
     474             :       {
     475     5774419 :         if (t == j) continue;
     476             : 
     477     4384266 :         piv = ( ucoeff(x,t,k) %= p );
     478     4384266 :         if (!piv) continue;
     479     1984359 :         if (piv == 1)
     480     2671137 :           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
     481             :         else
     482     4702730 :           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
     483             :       }
     484             :     }
     485             :   }
     486      565685 :   if (deplin==1) return NULL;
     487             : 
     488      565678 :   y = cgetg(r+1, t_MAT);
     489     1683024 :   for (j=k=1; j<=r; j++,k++)
     490             :   {
     491     1117345 :     GEN C = cgetg(n+1, t_VECSMALL);
     492             : 
     493     2204375 :     gel(y,j) = C; while (d[k]) k++;
     494     9006096 :     for (i=1; i<k; i++)
     495     7888751 :       if (d[i])
     496     2367358 :         uel(C,i) = ucoeff(x,d[i],k) % p;
     497             :       else
     498     5521393 :         uel(C,i) = 0UL;
     499     7331589 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     500             :   }
     501      565679 :   if (deplin == 2) {
     502           0 :     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
     503           0 :     for (i = j = 1; j <= n; j++)
     504           0 :       if (d[j]) pc[i++] = j;
     505           0 :     return mkvec2(y, pc);
     506             :   }
     507      565679 :   return y;
     508             : }
     509             : 
     510             : /* in place, destroy x */
     511             : static GEN
     512      767225 : Flm_ker_gauss(GEN x, ulong p, long deplin)
     513             : {
     514             :   GEN y, c, d;
     515             :   long i, j, k, r, t, m, n;
     516             :   ulong a, pi;
     517      767225 :   n = lg(x)-1;
     518      767225 :   if (!n) return cgetg(1,t_MAT);
     519      767218 :   if (SMALL_ULONG(p)) return Flm_ker_gauss_OK(x, p, deplin);
     520      144986 :   pi = get_Fl_red(p);
     521             : 
     522      144986 :   m=nbrows(x); r=0;
     523             : 
     524      144985 :   c = zero_zv(m);
     525      144985 :   d = cgetg(n+1, t_VECSMALL);
     526      144969 :   a = 0; /* for gcc -Wall */
     527      467192 :   for (k=1; k<=n; k++)
     528             :   {
     529      831571 :     for (j=1; j<=m; j++)
     530      728703 :       if (!c[j])
     531             :       {
     532      547214 :         a = ucoeff(x,j,k);
     533      547214 :         if (a) break;
     534             :       }
     535      322219 :     if (j > m)
     536             :     {
     537      102875 :       if (deplin==1) {
     538           7 :         c = cgetg(n+1, t_VECSMALL);
     539          21 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k);
     540           7 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     541           7 :         return c;
     542             :       }
     543      102868 :       r++; d[k] = 0;
     544             :     }
     545             :     else
     546             :     {
     547      219344 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     548      219352 :       c[j] = k; d[k] = j;
     549      219352 :       ucoeff(x,j,k) = p-1;
     550      219352 :       if (piv != 1)
     551      401401 :         for (i=k+1; i<=n; i++)
     552      185552 :           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
     553      885354 :       for (t=1; t<=m; t++)
     554             :       {
     555      665999 :         if (t == j) continue;
     556             : 
     557      446645 :         piv = ucoeff(x,t,k);
     558      446645 :         if (!piv) continue;
     559      272470 :         if (piv == 1)
     560       30499 :           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
     561             :         else
     562      663510 :           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
     563             :       }
     564             :     }
     565             :   }
     566      144973 :   if (deplin==1) return NULL;
     567             : 
     568      144966 :   y = cgetg(r+1, t_MAT);
     569      247825 :   for (j=k=1; j<=r; j++,k++)
     570             :   {
     571      102865 :     GEN C = cgetg(n+1, t_VECSMALL);
     572             : 
     573      165411 :     gel(y,j) = C; while (d[k]) k++;
     574      216490 :     for (i=1; i<k; i++)
     575      113629 :       if (d[i])
     576       78404 :         uel(C,i) = ucoeff(x,d[i],k);
     577             :       else
     578       35225 :         uel(C,i) = 0UL;
     579      165966 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     580             :   }
     581      144960 :   if (deplin == 2) {
     582      140237 :     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
     583      435922 :     for (i = j = 1; j <= n; j++)
     584      295687 :       if (d[j]) pc[i++] = j;
     585      140235 :     return mkvec2(y, pc);
     586             :   }
     587        4723 :   return y;
     588             : }
     589             : 
     590             : GEN
     591      362176 : Flm_intersect_i(GEN x, GEN y, ulong p)
     592             : {
     593      362176 :   long j, lx = lg(x);
     594             :   GEN z;
     595             : 
     596      362176 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     597      362162 :   z = Flm_ker(shallowconcat(x,y), p);
     598     1920345 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     599      362163 :   return Flm_mul(x,z,p);
     600             : }
     601             : GEN
     602           0 : Flm_intersect(GEN x, GEN y, ulong p)
     603             : {
     604           0 :   pari_sp av = avma;
     605           0 :   return gerepileupto(av, Flm_image(Flm_intersect_i(x, y, p), p));
     606             : }
     607             : 
     608             : static GEN
     609      420907 : Flm_ker_echelon(GEN x, ulong p, long pivots) {
     610      420907 :   pari_sp av = avma;
     611             :   GEN R, Rc, C, C1, C2, S, K;
     612      420907 :   long n = lg(x) - 1, r;
     613      420907 :   ulong pi = get_Fl_red(p);
     614      420907 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     615      420907 :   Rc = indexcompl(R, n);
     616      420907 :   C1 = rowpermute(C, R);
     617      420907 :   C2 = rowpermute(C, Rc);
     618      420907 :   S = Flm_lsolve_lower_unit_pre(C1, C2, p, pi);
     619      420907 :   K = vecpermute(shallowconcat(Flm_neg(S, p), matid_Flm(n - r)),
     620             :                  perm_inv(vecsmall_concat(R, Rc)));
     621      420907 :   K = Flm_transpose(K);
     622      420907 :   if (pivots)
     623       32201 :     return gerepilecopy(av, mkvec2(K, R));
     624      388706 :   return gerepileupto(av, K);
     625             : }
     626             : 
     627             : static GEN
     628       47078 : Flm_deplin_echelon(GEN x, ulong p) {
     629       47078 :   pari_sp av = avma;
     630             :   GEN R, Rc, C, C1, C2, s, v;
     631       47078 :   long i, n = lg(x) - 1, r;
     632       47078 :   ulong pi = get_Fl_red(p);
     633       47078 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     634       47078 :   if (r == n) return gc_NULL(av);
     635       47071 :   Rc = indexcompl(R, n);
     636       47071 :   i = Rc[1];
     637       47071 :   C1 = rowpermute(C, R);
     638       47071 :   C2 = rowslice(C, i, i);
     639       47071 :   s = Flm_row(Flm_lsolve_lower_unit_pre(C1, C2, p, pi), 1);
     640       47071 :   v = vecsmallpermute(vecsmall_concat(Flv_neg(s, p), vecsmall_ei(n - r, 1)),
     641             :                  perm_inv(vecsmall_concat(R, Rc)));
     642       47071 :   return gerepileuptoleaf(av, v);
     643             : }
     644             : 
     645             : static GEN
     646     1235209 : Flm_ker_i(GEN x, ulong p, long deplin, long inplace) {
     647     1235209 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     648      467985 :     switch(deplin) {
     649      388706 :     case 0: return Flm_ker_echelon(x, p, 0);
     650       47078 :     case 1: return Flm_deplin_echelon(x, p);
     651       32201 :     case 2: return Flm_ker_echelon(x, p, 1);
     652             :     }
     653      767224 :   return Flm_ker_gauss(inplace? x: Flm_copy(x), p, deplin);
     654             : }
     655             : 
     656             : GEN
     657      585885 : Flm_ker_sp(GEN x, ulong p, long deplin) { return Flm_ker_i(x, p, deplin, 1); }
     658             : GEN
     659      649324 : Flm_ker(GEN x, ulong p) { return Flm_ker_i(x, p, 0, 0); }
     660             : GEN
     661           0 : Flm_deplin(GEN x, ulong p) { return Flm_ker_i(x, p, 1, 0); }
     662             : 
     663             : /* in place, destroy a, SMALL_ULONG(p) is TRUE */
     664             : static ulong
     665        1743 : Flm_det_gauss_OK(GEN a, long nbco, ulong p)
     666             : {
     667        1743 :   long i,j,k, s = 1;
     668        1743 :   ulong q, x = 1;
     669             : 
     670        7266 :   for (i=1; i<nbco; i++)
     671             :   {
     672        7140 :     for(k=i; k<=nbco; k++)
     673             :     {
     674        6958 :       ulong c = ucoeff(a,k,i) % p;
     675        6958 :       ucoeff(a,k,i) = c;
     676        6958 :       if (c) break;
     677             :     }
     678       17528 :     for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
     679        5705 :     if (k > nbco) return ucoeff(a,i,i);
     680        5523 :     if (k != i)
     681             :     { /* exchange the lines s.t. k = i */
     682        3122 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     683         784 :       s = -s;
     684             :     }
     685        5523 :     q = ucoeff(a,i,i);
     686             : 
     687        5523 :     if (x & HIGHMASK) x %= p;
     688        5523 :     x *= q;
     689        5523 :     q = Fl_inv(q,p);
     690       18417 :     for (k=i+1; k<=nbco; k++)
     691             :     {
     692       12894 :       ulong m = ucoeff(a,i,k) % p;
     693       12894 :       if (!m) continue;
     694             : 
     695        8708 :       m = p - ((m*q)%p);
     696       34489 :       for (j=i+1; j<=nbco; j++)
     697             :       {
     698       25781 :         ulong c = ucoeff(a,j,k);
     699       25781 :         if (c & HIGHMASK) c %= p;
     700       25781 :         ucoeff(a,j,k) = c  + m*ucoeff(a,j,i);
     701             :       }
     702             :     }
     703             :   }
     704        1561 :   if (x & HIGHMASK) x %= p;
     705        1561 :   q = ucoeff(a,nbco,nbco);
     706        1561 :   if (q & HIGHMASK) q %= p;
     707        1561 :   x = (x*q) % p;
     708        1561 :   if (s < 0 && x) x = p - x;
     709        1561 :   return x;
     710             : }
     711             : 
     712             : /* in place, destroy a */
     713             : static ulong
     714       51835 : Flm_det_gauss(GEN a, ulong p)
     715             : {
     716       51835 :   long i,j,k, s = 1, nbco = lg(a)-1;
     717       51835 :   ulong pi, q, x = 1;
     718             : 
     719       51835 :   if (SMALL_ULONG(p)) return Flm_det_gauss_OK(a, nbco, p);
     720       50092 :   pi = get_Fl_red(p);
     721      321727 :   for (i=1; i<nbco; i++)
     722             :   {
     723      272038 :     for(k=i; k<=nbco; k++)
     724      272039 :       if (ucoeff(a,k,i)) break;
     725      271635 :     if (k > nbco) return ucoeff(a,i,i);
     726      271635 :     if (k != i)
     727             :     { /* exchange the lines s.t. k = i */
     728        1090 :       for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
     729         212 :       s = -s;
     730             :     }
     731      271635 :     q = ucoeff(a,i,i);
     732             : 
     733      271635 :     x = Fl_mul_pre(x, q, p, pi);
     734      271633 :     q = Fl_inv(q,p);
     735     1167449 :     for (k=i+1; k<=nbco; k++)
     736             :     {
     737      895814 :       ulong m = ucoeff(a,i,k);
     738      895814 :       if (!m) continue;
     739             : 
     740      841037 :       m = Fl_mul_pre(m, q, p, pi);
     741     4326566 :       for (j=i+1; j<=nbco; j++)
     742     3485514 :         ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
     743             :     }
     744             :   }
     745       50092 :   if (s < 0) x = Fl_neg(x, p);
     746       50092 :   return Fl_mul(x, ucoeff(a,nbco,nbco), p);
     747             : }
     748             : 
     749             : static ulong
     750       36042 : Flm_det_CUP(GEN a, ulong p) {
     751             :   GEN R, C, U, P;
     752       36042 :   long i, n = lg(a) - 1, r;
     753             :   ulong d;
     754       36042 :   ulong pi = get_Fl_red(p);
     755       36043 :   r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi);
     756       36045 :   if (r < n)
     757          42 :     d = 0;
     758             :   else {
     759       36003 :     d = perm_sign(P) == 1? 1: p-1;
     760      381039 :     for (i = 1; i <= n; i++)
     761      345053 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
     762             :   }
     763       36028 :   return d;
     764             : }
     765             : 
     766             : static ulong
     767       87877 : Flm_det_i(GEN x, ulong p, long inplace) {
     768       87877 :   pari_sp av = avma;
     769             :   ulong d;
     770       87877 :   if (lg(x) - 1 >= Flm_CUP_LIMIT)
     771       36042 :     d = Flm_det_CUP(x, p);
     772             :   else
     773       51835 :     d = Flm_det_gauss(inplace? x: Flm_copy(x), p);
     774       87881 :   return gc_ulong(av, d);
     775             : }
     776             : 
     777             : ulong
     778       87877 : Flm_det_sp(GEN x, ulong p) { return Flm_det_i(x, p, 1); }
     779             : ulong
     780           0 : Flm_det(GEN x, ulong p) { return Flm_det_i(x, p, 0); }
     781             : 
     782             : /* Destroy x */
     783             : static GEN
     784      619177 : Flm_gauss_pivot(GEN x, ulong p, long *rr)
     785             : {
     786             :   GEN c,d;
     787             :   long i,j,k,r,t,n,m;
     788             : 
     789      619177 :   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
     790             : 
     791      617868 :   m=nbrows(x); r=0;
     792      617869 :   d=cgetg(n+1,t_VECSMALL);
     793      617867 :   c = zero_zv(m);
     794     3740744 :   for (k=1; k<=n; k++)
     795             :   {
     796    11513014 :     for (j=1; j<=m; j++)
     797    10099033 :       if (!c[j])
     798             :       {
     799     6462447 :         ucoeff(x,j,k) %= p;
     800     6462447 :         if (ucoeff(x,j,k)) break;
     801             :       }
     802     3122874 :     if (j>m) { r++; d[k]=0; }
     803             :     else
     804             :     {
     805     1708891 :       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
     806     1708887 :       c[j]=k; d[k]=j;
     807     5556430 :       for (i=k+1; i<=n; i++)
     808     3847549 :         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
     809    13121797 :       for (t=1; t<=m; t++)
     810    11412932 :         if (!c[t]) /* no pivot on that line yet */
     811             :         {
     812     7464021 :           piv = ucoeff(x,t,k);
     813     7464021 :           if (piv)
     814             :           {
     815     4458828 :             ucoeff(x,t,k) = 0;
     816    13925463 :             for (i=k+1; i<=n; i++)
     817     9466661 :               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
     818     9466651 :                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
     819             :           }
     820             :         }
     821     7265310 :       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
     822             :     }
     823             :   }
     824      617870 :   *rr = r; return gc_const((pari_sp)d, d);
     825             : }
     826             : 
     827             : static GEN
     828      167284 : Flm_pivots_CUP(GEN x, ulong p, long *rr)
     829             : {
     830      167284 :   long i, n = lg(x) - 1, r;
     831      167284 :   GEN R, C, U, P, d = zero_zv(n);
     832      167284 :   ulong pi = get_Fl_red(p);
     833      167284 :   r = Flm_CUP_pre(x, &R, &C, &U, &P, p, pi);
     834     1840189 :   for(i = 1; i <= r; i++)
     835     1672905 :     d[P[i]] = R[i];
     836      167284 :   *rr = n - r; return gc_const((pari_sp)d, d);
     837             : }
     838             : 
     839             : GEN
     840      786459 : Flm_pivots(GEN x, ulong p, long *rr, long inplace)
     841             : {
     842      786459 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     843      167284 :     return Flm_pivots_CUP(x, p, rr);
     844      619175 :   return Flm_gauss_pivot(inplace? x: Flm_copy(x), p, rr);
     845             : }
     846             : 
     847             : long
     848       39179 : Flm_rank(GEN x, ulong p)
     849             : {
     850       39179 :   pari_sp av = avma;
     851             :   long r;
     852       39179 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT) {
     853             :     GEN R, C;
     854        9485 :     ulong pi = get_Fl_red(p);
     855        9485 :     return gc_long(av, Flm_echelon_pre(x, &R, &C, p, pi));
     856             :   }
     857       29694 :   (void) Flm_pivots(x, p, &r, 0);
     858       29694 :   return gc_long(av, lg(x)-1 - r);
     859             : }
     860             : 
     861             : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
     862             :  * reduced mod p */
     863             : static GEN
     864        1225 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
     865             : {
     866        1225 :   long n = lg(A)-1, i = index, j;
     867        1225 :   GEN u = const_vecsmall(n, 0);
     868        1225 :   u[i] = 1;
     869        1225 :   if (SMALL_ULONG(p))
     870        3178 :     for (i--; i>0; i--)
     871             :     {
     872        1967 :       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
     873        6937 :       for (j=i+2; j<=n; j++)
     874             :       {
     875        4970 :         if (m & HIGHMASK) m %= p;
     876        4970 :         m += ucoeff(A,i,j) * uel(u,j);
     877             :       }
     878        1967 :       u[i] = Fl_neg(m % p, p);
     879             :     }
     880             :   else
     881          21 :     for (i--; i>0; i--)
     882             :     {
     883           7 :       ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
     884           7 :       for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
     885           7 :       u[i] = Fl_neg(m, p);
     886             :     }
     887        1225 :   return u;
     888             : }
     889             : static GEN
     890         462 : Flm_inv_upper_1(GEN A, ulong p)
     891             : {
     892             :   long i, l;
     893         462 :   GEN B = cgetg_copy(A, &l);
     894        1687 :   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
     895         462 :   return B;
     896             : }
     897             : /* destroy a, b */
     898             : static GEN
     899       98938 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
     900             : {
     901       98938 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     902       98938 :   ulong det = 1;
     903             :   GEN u;
     904             : 
     905       98938 :   li = nbrows(a);
     906       98938 :   bco = lg(b)-1;
     907      358780 :   for (i=1; i<=aco; i++)
     908             :   {
     909             :     ulong invpiv;
     910             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
     911      891157 :     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
     912      411018 :     for (k = i; k <= li; k++)
     913             :     {
     914      411012 :       ulong piv = ( ucoeff(a,k,i) %= p );
     915      411012 :       if (piv)
     916             :       {
     917      358780 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     918      358777 :         if (detp)
     919             :         {
     920           0 :           if (det & HIGHMASK) det %= p;
     921           0 :           det *= piv;
     922             :         }
     923      358777 :         break;
     924             :       }
     925             :     }
     926             :     /* found a pivot on line k */
     927      358783 :     if (k > li) return NULL;
     928      358776 :     if (k != i)
     929             :     { /* swap lines so that k = i */
     930       35531 :       s = -s;
     931      145560 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     932      189074 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     933             :     }
     934      358776 :     if (i == aco) break;
     935             : 
     936      259845 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
     937      814112 :     for (k=i+1; k<=li; k++)
     938             :     {
     939      554270 :       ulong m = ( ucoeff(a,k,i) %= p );
     940      554270 :       if (!m) continue;
     941             : 
     942      158652 :       m = Fl_mul(m, invpiv, p);
     943      158652 :       if (m == 1) {
     944      135678 :         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
     945      211721 :         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
     946             :       } else {
     947      424772 :         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
     948      712031 :         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
     949             :       }
     950             :     }
     951             :   }
     952       98925 :   if (detp)
     953             :   {
     954           0 :     det %=  p;
     955           0 :     if (s < 0 && det) det = p - det;
     956           0 :     *detp = det;
     957             :   }
     958       98925 :   u = cgetg(bco+1,t_MAT);
     959      471162 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
     960       98931 :   return u;
     961             : }
     962             : 
     963             : /* destroy a, b */
     964             : static GEN
     965      998824 : Flm_gauss_sp_i(GEN a, GEN b, ulong *detp, ulong p)
     966             : {
     967      998824 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     968      998824 :   ulong det = 1;
     969             :   GEN u;
     970             :   ulong pi;
     971      998824 :   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
     972      998824 :   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
     973      899886 :   pi = get_Fl_red(p);
     974      899886 :   li = nbrows(a);
     975      899886 :   bco = lg(b)-1;
     976     2613167 :   for (i=1; i<=aco; i++)
     977             :   {
     978             :     ulong invpiv;
     979             :     /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
     980     2718621 :     for (k = i; k <= li; k++)
     981             :     {
     982     2718616 :       ulong piv = ucoeff(a,k,i);
     983     2718616 :       if (piv)
     984             :       {
     985     2613164 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     986     2613167 :         if (detp) det = Fl_mul_pre(det, piv, p, pi);
     987     2613166 :         break;
     988             :       }
     989             :     }
     990             :     /* found a pivot on line k */
     991     2613171 :     if (k > li) return NULL;
     992     2613164 :     if (k != i)
     993             :     { /* swap lines so that k = i */
     994       93877 :       s = -s;
     995      364926 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     996      415893 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     997             :     }
     998     2613164 :     if (i == aco) break;
     999             : 
    1000     1713269 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    1001     4593411 :     for (k=i+1; k<=li; k++)
    1002             :     {
    1003     2880130 :       ulong m = ucoeff(a,k,i);
    1004     2880130 :       if (!m) continue;
    1005             : 
    1006     2318617 :       m = Fl_mul_pre(m, invpiv, p, pi);
    1007     2318617 :       if (m == 1) {
    1008      293176 :         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
    1009      417562 :         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
    1010             :       } else {
    1011     6771246 :         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
    1012    10286039 :         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
    1013             :       }
    1014             :     }
    1015             :   }
    1016      899893 :   if (detp)
    1017             :   {
    1018           0 :     if (s < 0 && det) det = p - det;
    1019           0 :     *detp = det;
    1020             :   }
    1021      899893 :   u = cgetg(bco+1,t_MAT);
    1022     3513038 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
    1023      899885 :   return u;
    1024             : }
    1025             : 
    1026             : static GEN
    1027      765257 : Flm_gauss_from_CUP(GEN b, GEN R, GEN C, GEN U, GEN P, ulong p, ulong pi, ulong *detp)
    1028             : {
    1029      765257 :   GEN Y = Flm_rsolve_lower_unit_pre(rowpermute(C, R), rowpermute(b, R), p, pi);
    1030      765270 :   GEN X = rowpermute(Flm_rsolve_upper_pre(U, Y, p, pi), perm_inv(P));
    1031      765273 :   if (detp)
    1032             :   {
    1033      648073 :     ulong d = perm_sign(P) == 1? 1: p-1;
    1034      648077 :     long i, r = lg(R);
    1035     4245263 :     for (i = 1; i < r; i++)
    1036     3597193 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
    1037      648070 :     *detp = d;
    1038             :   }
    1039      765270 :   return X;
    1040             : }
    1041             : 
    1042             : static GEN
    1043      117204 : Flm_gauss_CUP(GEN a, GEN b, ulong *detp, ulong p) {
    1044             :   GEN R, C, U, P;
    1045      117204 :   long n = lg(a) - 1, r;
    1046      117204 :   ulong pi = get_Fl_red(p);
    1047      117204 :   if (nbrows(a) < n || (r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi)) < n)
    1048          14 :     return NULL;
    1049      117190 :   return Flm_gauss_from_CUP(b, R, C, U, P, p, pi, detp);
    1050             : }
    1051             : 
    1052             : GEN
    1053     1052634 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p) {
    1054     1052634 :   pari_sp av = avma;
    1055             :   GEN x;
    1056     1052634 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1057       53809 :     x = Flm_gauss_CUP(a, b, detp, p);
    1058             :   else
    1059      998825 :     x = Flm_gauss_sp_i(a, b, detp, p);
    1060     1052650 :   if (!x) return gc_NULL(av);
    1061     1052636 :   return gerepileupto(av, x);
    1062             : }
    1063             : 
    1064             : GEN
    1065          63 : Flm_gauss(GEN a, GEN b, ulong p) {
    1066          63 :   pari_sp av = avma;
    1067             :   GEN x;
    1068          63 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1069          14 :     x = Flm_gauss_CUP(a, b, NULL, p);
    1070             :   else {
    1071          49 :     a = RgM_shallowcopy(a);
    1072          49 :     b = RgM_shallowcopy(b);
    1073          49 :     x = Flm_gauss_sp(a, b, NULL, p);
    1074             :   }
    1075          63 :   if (!x) return gc_NULL(av);
    1076          49 :   return gerepileupto(av, x);
    1077             : }
    1078             : 
    1079             : static GEN
    1080     1006679 : Flm_inv_i(GEN a, ulong *detp, ulong p, long inplace) {
    1081     1006679 :   pari_sp av = avma;
    1082     1006679 :   long n = lg(a) - 1;
    1083             :   GEN b, x;
    1084     1006679 :   if (n == 0) return cgetg(1, t_MAT);
    1085     1006679 :   b = matid_Flm(nbrows(a));
    1086     1006678 :   if (n >= Flm_CUP_LIMIT)
    1087       63381 :     x = Flm_gauss_CUP(a, b, detp, p);
    1088             :   else {
    1089      943297 :     if (!inplace)
    1090       43442 :       a = RgM_shallowcopy(a);
    1091      943298 :     x = Flm_gauss_sp(a, b, detp, p);
    1092             :   }
    1093     1006689 :   if (!x) return gc_NULL(av);
    1094     1006675 :   return gerepileupto(av, x);
    1095             : }
    1096             : 
    1097             : GEN
    1098      948467 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
    1099      948467 :   return Flm_inv_i(a, detp, p, 1);
    1100             : }
    1101             : 
    1102             : GEN
    1103       58210 : Flm_inv(GEN a, ulong p) {
    1104       58210 :   return Flm_inv_i(a, NULL, p, 0);
    1105             : }
    1106             : 
    1107             : GEN
    1108          21 : Flm_Flc_gauss(GEN a, GEN b, ulong p) {
    1109          21 :   pari_sp av = avma;
    1110          21 :   GEN z = Flm_gauss(a, mkmat(b), p);
    1111          21 :   if (!z) return gc_NULL(av);
    1112          14 :   if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
    1113          14 :   return gerepileuptoleaf(av, gel(z,1));
    1114             : }
    1115             : 
    1116             : GEN
    1117      648088 : Flm_adjoint(GEN A, ulong p)
    1118             : {
    1119      648088 :   pari_sp av = avma;
    1120             :   GEN R, C, U, P, C1, U1, v, c, d;
    1121      648088 :   long r, i, q, n = lg(A)-1, m;
    1122             :   ulong D;
    1123      648088 :   ulong pi = get_Fl_red(p);
    1124      648088 :   if (n == 0) return cgetg(1, t_MAT);
    1125      648088 :   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
    1126      648086 :   m = nbrows(A);
    1127      648085 :   if (r == n)
    1128             :   {
    1129      648071 :     GEN X = Flm_gauss_from_CUP(matid_Flm(m), R, C, U, P, p, pi, &D);
    1130      648078 :     return gerepileupto(av, Flm_Fl_mul_pre(X, D, p, pi));
    1131             :   }
    1132          14 :   if (r < n-1) return zero_Flm(n, m);
    1133          28 :   for (q = n, i = 1; i < n; i++)
    1134          14 :     if (R[i] != i) { q = i; break; }
    1135          14 :   C1 = matslice(C, 1, q-1, 1, q-1);
    1136          14 :   c = vecslice(Flm_row(C, q), 1, q-1);
    1137          14 :   c = Flm_lsolve_lower_unit_pre(C1, Flm_transpose(mkmat(c)), p, pi);
    1138          14 :   d = cgetg(m+1, t_VECSMALL);
    1139          28 :   for (i=1; i<q; i++)    uel(d,i) = ucoeff(c,1,i);
    1140          14 :   uel(d,q) = p-1;
    1141          21 :   for (i=q+1; i<=m; i++) uel(d,i) = 0;
    1142          14 :   U1 = vecslice(U, 1, n-1);
    1143          14 :   v = gel(Flm_rsolve_upper_pre(U1, mkmat(gel(U,n)), p, pi),1);
    1144          14 :   v = vecsmall_append(v, p-1);
    1145          14 :   D = perm_sign(P) != (odd(q+n)?-1:1) ? p-1 : 1;
    1146          28 :   for (i = 1; i <= n-1; i++)
    1147          14 :     D = Fl_mul_pre(D, ucoeff(U1, i, i), p, pi);
    1148          14 :   d = Flv_Fl_mul(d, D, p);
    1149          14 :   return rowpermute(Flc_Flv_mul(v, d, p),perm_inv(P));
    1150             : }
    1151             : 
    1152             : static GEN
    1153         287 : Flm_invimage_CUP(GEN A, GEN B, ulong p) {
    1154         287 :   pari_sp av = avma;
    1155             :   GEN R, Rc, C, U, P, B1, B2, C1, C2, X, Y, Z;
    1156             :   long r;
    1157         287 :   ulong pi = get_Fl_red(p);
    1158         287 :   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
    1159         287 :   Rc = indexcompl(R, nbrows(B));
    1160         287 :   C1 = rowpermute(C, R);
    1161         287 :   C2 = rowpermute(C, Rc);
    1162         287 :   B1 = rowpermute(B, R);
    1163         287 :   B2 = rowpermute(B, Rc);
    1164         287 :   Z = Flm_rsolve_lower_unit_pre(C1, B1, p, pi);
    1165         287 :   if (!gequal(Flm_mul_pre(C2, Z, p, pi), B2))
    1166          14 :     return NULL;
    1167         273 :   Y = vconcat(Flm_rsolve_upper_pre(vecslice(U, 1, r), Z, p, pi),
    1168         273 :               zero_Flm(lg(A) - 1 - r, lg(B) - 1));
    1169         273 :   X = rowpermute(Y, perm_inv(P));
    1170         273 :   return gerepileupto(av, X);
    1171             : }
    1172             : 
    1173             : GEN
    1174         896 : Flm_invimage_i(GEN A, GEN B, ulong p)
    1175             : {
    1176             :   GEN d, x, X, Y;
    1177         896 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    1178             : 
    1179         896 :   if (!nB) return cgetg(1, t_MAT);
    1180         756 :   if (maxss(nA, nB) >= Flm_CUP_LIMIT && nbrows(B) >= Flm_CUP_LIMIT)
    1181         287 :     return Flm_invimage_CUP(A, B, p);
    1182             : 
    1183         469 :   x = Flm_ker(shallowconcat(Flm_neg(A,p), B), p);
    1184             :   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    1185             :    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
    1186             :    * Y has at least nB columns and full rank */
    1187         469 :   nY = lg(x)-1;
    1188         469 :   if (nY < nB) return NULL;
    1189         469 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    1190         469 :   d = cgetg(nB+1, t_VECSMALL);
    1191        1701 :   for (i = nB, j = nY; i >= 1; i--, j--)
    1192             :   {
    1193        1246 :     for (; j>=1; j--)
    1194        1239 :       if (coeff(Y,i,j)) { d[i] = j; break; }
    1195        1239 :     if (!j) return NULL;
    1196             :   }
    1197             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    1198         462 :   Y = vecpermute(Y, d);
    1199         462 :   x = vecpermute(x, d);
    1200         462 :   X = rowslice(x, 1, nA);
    1201         462 :   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
    1202             : }
    1203             : GEN
    1204         854 : Flm_invimage(GEN A, GEN B, ulong p)
    1205             : {
    1206         854 :   pari_sp av = avma;
    1207         854 :   GEN X = Flm_invimage_i(A,B,p);
    1208         854 :   if (!X) return gc_NULL(av);
    1209         854 :   return gerepileupto(av, X);
    1210             : }
    1211             : 
    1212             : GEN
    1213      147806 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
    1214             : {
    1215      147806 :   pari_sp av = avma;
    1216      147806 :   long i, l = lg(A);
    1217             :   GEN M, x;
    1218             :   ulong t;
    1219             : 
    1220      147806 :   if (l==1) return NULL;
    1221      147680 :   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
    1222      147680 :   M = cgetg(l+1,t_MAT);
    1223     1313951 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    1224      147680 :   gel(M,l) = y; M = Flm_ker(M,p);
    1225      147680 :   i = lg(M)-1; if (!i) return gc_NULL(av);
    1226             : 
    1227      144502 :   x = gel(M,i); t = x[l];
    1228      144502 :   if (!t) return gc_NULL(av);
    1229             : 
    1230      144495 :   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
    1231      144495 :   if (t!=1) x = Flv_Fl_mul(x, t, p);
    1232      144495 :   return gerepileuptoleaf(av, x);
    1233             : }

Generated by: LCOV version 1.13