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.18.0 lcov report (development 29785-cb4530e674) Lines: 737 758 97.2 %
Date: 2024-12-02 09:09:34 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      388511 : Flv_to_ZV(GEN x)
      20     4096924 : { pari_APPLY_type(t_VEC, utoi(x[i])) }
      21             : 
      22             : GEN
      23    28810332 : Flc_to_ZC(GEN x)
      24   240461824 : { pari_APPLY_type(t_COL, utoi(x[i])) }
      25             : 
      26             : GEN
      27    11411100 : Flm_to_ZM(GEN x)
      28    39273114 : { 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     1797621 : Flm_solve_upper_1(GEN U, GEN B, ulong p, ulong pi)
      49     1797621 : { return Flm_Fl_mul_pre(B, Fl_inv(ucoeff(U, 1, 1), p), p, pi); }
      50             : 
      51             : static GEN
      52     2790762 : Flm_rsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      53             : {
      54     2790762 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      55     2790762 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      56     2790932 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      57     2790911 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
      58     2790888 :   GEN B1 = rowslice(B, 1, 1);
      59     2790715 :   GEN B2 = rowslice(B, 2, 2);
      60     2790705 :   GEN X2 = Flm_Fl_mul_pre(B2, dinv, p, pi);
      61     2790861 :   GEN X1 = Flm_Fl_mul_pre(Flm_sub(B1, Flm_Fl_mul_pre(X2, b, p, pi), p),
      62             :                       ainv, p, pi);
      63     2790884 :   return vconcat(X1, X2);
      64             : }
      65             : 
      66             : /* solve U*X = B,  U upper triangular and invertible */
      67             : static GEN
      68     6687350 : Flm_rsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
      69             : {
      70     6687350 :   long n = lg(U) - 1, n1;
      71             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
      72     6687350 :   pari_sp av = avma;
      73             : 
      74     6687350 :   if (n == 0) return B;
      75     6687343 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
      76     5591040 :   if (n == 2) return Flm_rsolve_upper_2(U, B, p, pi);
      77     2800282 :   n1 = (n + 1)/2;
      78     2800282 :   U2 = vecslice(U, n1 + 1, n);
      79     2801175 :   U22 = rowslice(U2, n1 + 1, n);
      80     2801078 :   B2 = rowslice(B, n1 + 1, n);
      81     2801023 :   X2 = Flm_rsolve_upper_pre(U22, B2, p, pi);
      82     2801369 :   U12 = rowslice(U2, 1, n1);
      83     2801121 :   B1 = rowslice(B, 1, n1);
      84     2801094 :   B1 = Flm_sub(B1, Flm_mul_pre(U12, X2, p, pi), p);
      85     2801104 :   if (gc_needed(av, 1)) gerepileall(av, 2, &B1, &X2);
      86     2801104 :   U11 = matslice(U, 1,n1, 1,n1);
      87     2801061 :   X1 = Flm_rsolve_upper_pre(U11, B1, p, pi);
      88     2801415 :   X = vconcat(X1, X2);
      89     2801425 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
      90     2801425 :   return X;
      91             : }
      92             : 
      93             : static GEN
      94     2537679 : Flm_lsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
      95             : {
      96     2537679 :   ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
      97     2537679 :   ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
      98     2537751 :   ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
      99     2537737 :   ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
     100     2537724 :   GEN B1 = vecslice(B, 1, 1);
     101     2537688 :   GEN B2 = vecslice(B, 2, 2);
     102     2537640 :   GEN X1 = Flm_Fl_mul_pre(B1, ainv, p, pi);
     103     2537693 :   GEN X2 = Flm_Fl_mul_pre(Flm_sub(B2, Flm_Fl_mul_pre(X1, b, p, pi), p),
     104             :                       dinv, p, pi);
     105     2537709 :   return shallowconcat(X1, X2);
     106             : }
     107             : 
     108             : /* solve X*U = B,  U upper triangular and invertible */
     109             : static GEN
     110     5466759 : Flm_lsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
     111             : {
     112     5466759 :   long n = lg(U) - 1, n1;
     113             :   GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
     114     5466759 :   pari_sp av = avma;
     115             : 
     116     5466759 :   if (n == 0) return B;
     117     5466759 :   if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
     118     4765393 :   if (n == 2) return Flm_lsolve_upper_2(U, B, p, pi);
     119     2227711 :   n1 = (n + 1)/2;
     120     2227711 :   U2 = vecslice(U, n1 + 1, n);
     121     2227991 :   U11 = matslice(U, 1,n1, 1,n1);
     122     2227961 :   U12 = rowslice(U2, 1, n1);
     123     2227957 :   U22 = rowslice(U2, n1 + 1, n);
     124     2227971 :   B1 = vecslice(B, 1, n1);
     125     2227932 :   B2 = vecslice(B, n1 + 1, n);
     126     2227934 :   X1 = Flm_lsolve_upper_pre(U11, B1, p, pi);
     127     2227996 :   B2 = Flm_sub(B2, Flm_mul_pre(X1, U12, p, pi), p);
     128     2227955 :   if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
     129     2227955 :   X2 = Flm_lsolve_upper_pre(U22, B2, p, pi);
     130     2228024 :   X = shallowconcat(X1, X2);
     131     2228015 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     132     2228016 :   return X;
     133             : }
     134             : 
     135             : static GEN
     136     4795759 : Flm_rsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     137             : {
     138     4795759 :   GEN X1 = rowslice(A, 1, 1);
     139     4795841 :   GEN X2 = Flm_sub(rowslice(A, 2, 2),
     140     4795684 :                    Flm_Fl_mul_pre(X1, ucoeff(L, 2, 1), p, pi), p);
     141     4795607 :   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    10981711 : Flm_rsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     148             : {
     149    10981711 :   long m = lg(L) - 1, m1, n;
     150             :   GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
     151    10981711 :   pari_sp av = avma;
     152             : 
     153    10981711 :   if (m == 0) return zero_Flm(0, lg(A) - 1);
     154    10981711 :   if (m == 1) return rowslice(A, 1, 1);
     155     9386263 :   if (m == 2) return Flm_rsolve_lower_unit_2(L, A, p, pi);
     156     4590511 :   m1 = (m + 1)/2;
     157     4590511 :   n = nbrows(L);
     158     4591319 :   L1 = vecslice(L, 1, m1);
     159     4591359 :   L11 = rowslice(L1, 1, m1);
     160     4591266 :   L21 = rowslice(L1, m1 + 1, n);
     161     4591223 :   A1 = rowslice(A, 1, m1);
     162     4591272 :   X1 = Flm_rsolve_lower_unit_pre(L11, A1, p, pi);
     163     4591518 :   A2 = rowslice(A, m1 + 1, n);
     164     4591158 :   A2 = Flm_sub(A2, Flm_mul_pre(L21, X1, p, pi), p);
     165     4591161 :   if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
     166     4591161 :   L22 = matslice(L, m1+1,n, m1+1,m);
     167     4591180 :   X2 = Flm_rsolve_lower_unit_pre(L22, A2, p, pi);
     168     4591347 :   X = vconcat(X1, X2);
     169     4591554 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     170     4591551 :   return X;
     171             : }
     172             : 
     173             : static GEN
     174     1938664 : Flm_lsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
     175             : {
     176     1938664 :   GEN X2 = vecslice(A, 2, 2);
     177     1938662 :   GEN X1 = Flm_sub(vecslice(A, 1, 1),
     178     1938662 :                    Flm_Fl_mul_pre(X2, ucoeff(L, 2, 1), p, pi), p);
     179     1938663 :   return shallowconcat(X1, X2);
     180             : }
     181             : 
     182             : /* solve L*X = A,  L square lower triangular with ones on the diagonal */
     183             : static GEN
     184     4503551 : Flm_lsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
     185             : {
     186     4503551 :   long m = lg(L) - 1, m1;
     187             :   GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
     188     4503551 :   pari_sp av = avma;
     189             : 
     190     4503551 :   if (m <= 1) return A;
     191     4015026 :   if (m == 2) return Flm_lsolve_lower_unit_2(L, A, p, pi);
     192     2076361 :   m1 = (m + 1)/2;
     193     2076361 :   L2 = vecslice(L, m1 + 1, m);
     194     2076380 :   L22 = rowslice(L2, m1 + 1, m);
     195     2076376 :   A2 = vecslice(A, m1 + 1, m);
     196     2076373 :   X2 = Flm_lsolve_lower_unit_pre(L22, A2, p, pi);
     197     2076384 :   if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
     198     2076384 :   L1 = vecslice(L, 1, m1);
     199     2076384 :   L21 = rowslice(L1, m1 + 1, m);
     200     2076384 :   A1 = vecslice(A, 1, m1);
     201     2076381 :   A1 = Flm_sub(A1, Flm_mul_pre(X2, L21, p, pi), p);
     202     2076381 :   L11 = rowslice(L1, 1, m1);
     203     2076384 :   if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
     204     2076384 :   X1 = Flm_lsolve_lower_unit_pre(L11, A1, p, pi);
     205     2076383 :   X = shallowconcat(X1, X2);
     206     2076383 :   if (gc_needed(av, 1)) X = gerepilecopy(av, X);
     207     2076383 :   return X;
     208             : }
     209             : 
     210             : /* destroy A */
     211             : static long
     212     3584969 : Flm_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     213             : {
     214     3584969 :   long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
     215             :   ulong u, v;
     216             : 
     217     3584969 :   if (P) *P = identity_perm(n);
     218     3584963 :   *R = cgetg(m + 1, t_VECSMALL);
     219    17885328 :   for (j = 1, pr = 0; j <= n; j++)
     220             :   {
     221    20517194 :     for (pr++, pc = 0; pr <= m; pr++)
     222             :     {
     223   138795868 :       for (k = j; k <= n; k++) { v = ucoeff(A, pr, k); if (!pc && v) pc = k; }
     224    19126003 :       if (pc) break;
     225             :     }
     226    15691742 :     if (!pc) break;
     227    14300452 :     (*R)[j] = pr;
     228    14300452 :     if (pc != j)
     229             :     {
     230     2293318 :       swap(gel(A, j), gel(A, pc));
     231     2293318 :       if (P) lswap((*P)[j], (*P)[pc]);
     232             :     }
     233    14300452 :     u = Fl_inv(ucoeff(A, pr, j), p);
     234    97137671 :     for (i = pr + 1; i <= m; i++)
     235             :     {
     236    82837302 :       v = Fl_mul_pre(ucoeff(A, i, j), u, p, pi);
     237    82836731 :       ucoeff(A, i, j) = v;
     238    82836731 :       v = Fl_neg(v, p);
     239   383277525 :       for (k = j + 1; k <= n; k++)
     240   300445670 :         ucoeff(A, i, k) = Fl_addmul_pre(ucoeff(A, i, k),
     241   300440477 :                                         ucoeff(A, pr, k), v, p, pi);
     242             :     }
     243             :   }
     244     3584876 :   setlg(*R, j);
     245     3584978 :   *C = vecslice(A, 1, j - 1);
     246     3584958 :   if (U) *U = rowpermute(A, *R);
     247     3584948 :   return j - 1;
     248             : }
     249             : 
     250             : static const long Flm_CUP_LIMIT = 8;
     251             : 
     252             : static long
     253     3473893 : Flm_CUP_pre(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
     254             : {
     255     3473893 :   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     3473889 :   pari_sp av = avma;
     259             : 
     260     3473889 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     261             :     /* destroy A; not called at the outermost recursion level */
     262     2435604 :     return Flm_CUP_basecase(A, R, C, U, P, p, pi);
     263     1038285 :   m1 = (minss(m, n) + 1)/2;
     264     1038285 :   A1 = rowslice(A, 1, m1);
     265     1038281 :   A2 = rowslice(A, m1 + 1, m);
     266     1038274 :   r1 = Flm_CUP_pre(A1, &R1, &C1, &U1, &P1, p, pi);
     267     1038271 :   if (r1 == 0)
     268             :   {
     269       27197 :     r2 = Flm_CUP_pre(A2, &R2, &C2, &U2, &P2, p, pi);
     270       27197 :     *R = cgetg(r2 + 1, t_VECSMALL);
     271       36915 :     for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
     272       27197 :     *C = vconcat(zero_Flm(m1, r2), C2);
     273       27197 :     *U = U2;
     274       27197 :     *P = P2;
     275       27197 :     r = r2;
     276             :   }
     277             :   else
     278             :   {
     279     1011074 :     U11 = vecslice(U1, 1, r1);
     280     1011074 :     U12 = vecslice(U1, r1 + 1, n);
     281     1011078 :     T21 = vecslicepermute(A2, P1, 1, r1);
     282     1011078 :     T22 = vecslicepermute(A2, P1, r1 + 1, n);
     283     1011076 :     C21 = Flm_lsolve_upper_pre(U11, T21, p, pi);
     284     1011082 :     if (gc_needed(av, 1))
     285           0 :       gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
     286     1011082 :     B2 = Flm_sub(T22, Flm_mul_pre(C21, U12, p, pi), p);
     287     1011076 :     r2 = Flm_CUP_pre(B2, &R2, &C2, &U2, &P2, p, pi);
     288     1011082 :     r = r1 + r2;
     289     1011082 :     *R = cgetg(r + 1, t_VECSMALL);
     290     6787959 :     for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
     291     6319306 :     for (;      i <= r; i++)  (*R)[i] = R2[i - r1] + m1;
     292     1011078 :     *C = shallowconcat(vconcat(C1, C21),
     293             :                        vconcat(zero_Flm(m1, r2), C2));
     294     1011086 :     *U = shallowconcat(vconcat(U11, zero_Flm(r2, r1)),
     295             :                        vconcat(vecpermute(U12, P2), U2));
     296     1011092 :     *P = cgetg(n + 1, t_VECSMALL);
     297     6788039 :     for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
     298    12643799 :     for (; i <= n; i++)       (*P)[i] = P1[P2[i - r1] + r1];
     299             :   }
     300     1038287 :   if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
     301     1038287 :   return r;
     302             : }
     303             : 
     304             : static long
     305     1149327 : Flm_echelon_gauss(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     306     1149327 : { 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     1065040 : indexcompl(GEN v, long n)
     311             : {
     312     1065040 :   long i, j, k, m = lg(v) - 1;
     313     1065040 :   GEN w = cgetg(n - m + 1, t_VECSMALL);
     314    23071727 :   for (i = j = k = 1; i <= n; i++)
     315    22006687 :     if (j <= m && v[j] == i) j++; else w[k++] = i;
     316     1065040 :   return w;
     317             : }
     318             : 
     319             : /* column echelon form */
     320             : static long
     321     1943292 : Flm_echelon_pre(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
     322             : {
     323     1943292 :   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     1943292 :   pari_sp av = avma;
     327             : 
     328     1943292 :   if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
     329     1149327 :     return Flm_echelon_gauss(Flm_copy(A), R, C, p, pi);
     330             : 
     331      793965 :   n1 = (n + 1)/2;
     332      793965 :   A1 = vecslice(A, 1, n1);
     333      793964 :   A2 = vecslice(A, n1 + 1, n);
     334      793964 :   r1 = Flm_echelon_pre(A1, &R1, &C1, p, pi);
     335      793965 :   if (!r1) return Flm_echelon_pre(A2, R, C, p, pi);
     336      719086 :   if (r1 == m) { *R = R1; *C = C1; return r1; }
     337             : 
     338      713961 :   R1c = indexcompl(R1, m);
     339      713961 :   C11 = rowpermute(C1, R1);
     340      713960 :   C21 = rowpermute(C1, R1c);
     341      713960 :   A12 = rowpermute(A2, R1);
     342      713960 :   A22 = rowpermute(A2, R1c);
     343      713959 :   M12 = Flm_rsolve_lower_unit_pre(C11, A12, p, pi);
     344      713960 :   B2 = Flm_sub(A22, Flm_mul_pre(C21, M12, p, pi), p);
     345      713961 :   r2 = Flm_echelon_pre(B2, &R2, &C2, p, pi);
     346      713960 :   if (!r2) { *R = R1; *C = C1; r = r1; }
     347             :   else
     348             :   {
     349      668507 :     R2 = perm_mul(R1c, R2);
     350      668507 :     C2 = rowpermute(vconcat(zero_Flm(r1, r2), C2),
     351             :                     perm_inv(vecsmall_concat(R1, R1c)));
     352      668506 :     r = r1 + r2;
     353      668506 :     *R = cgetg(r + 1, t_VECSMALL);
     354      668506 :     *C = cgetg(r + 1, t_MAT);
     355     8380000 :     for (j = j1 = j2 = 1; j <= r; j++)
     356     7711494 :       if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
     357             :       {
     358     4083177 :         gel(*C, j) = gel(C1, j1);
     359     4083177 :         (*R)[j] = R1[j1++];
     360             :       }
     361             :       else
     362             :       {
     363     3628317 :         gel(*C, j) = gel(C2, j2);
     364     3628317 :         (*R)[j] = R2[j2++];
     365             :       }
     366             :   }
     367      713959 :   if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
     368      713959 :   return r;
     369             : }
     370             : 
     371             : static void /* assume m < p */
     372    11845343 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
     373             : {
     374    11845343 :   uel(b,k) = Fl_addmul_pre(uel(b, k), m, uel(b, i), p, pi);
     375    11845443 : }
     376             : static void /* same m = 1 */
     377      667145 : _Fl_add(GEN b, long k, long i, ulong p)
     378             : {
     379      667145 :   uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
     380      667145 : }
     381             : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     382     4427017 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
     383             : {
     384     4427017 :   uel(b,k) += m * uel(b,i);
     385     4427017 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     386     4427017 : }
     387             : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
     388     2196332 : _Fl_add_OK(GEN b, long k, long i, ulong p)
     389             : {
     390     2196332 :   uel(b,k) += uel(b,i);
     391     2196332 :   if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
     392     2196332 : }
     393             : 
     394             : /* assume 0 <= a[i,j] < p */
     395             : static GEN
     396      440333 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
     397             : {
     398      440333 :   GEN u = cgetg(li+1,t_VECSMALL);
     399      440333 :   ulong m = uel(b,li) % p;
     400             :   long i,j;
     401             : 
     402      440333 :   uel(u,li) = (m * ucoeff(a,li,li)) % p;
     403     1700368 :   for (i = li-1; i > 0; i--)
     404             :   {
     405     1260035 :     m = p - uel(b,i)%p;
     406     4040035 :     for (j = i+1; j <= li; j++) {
     407     2780000 :       if (m & HIGHBIT) m %= p;
     408     2780000 :       m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
     409             :     }
     410     1260035 :     m %= p;
     411     1260035 :     if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
     412     1260035 :     uel(u,i) = m;
     413             :   }
     414      440333 :   return u;
     415             : }
     416             : static GEN
     417     2164000 : Fl_get_col(GEN a, GEN b, long li, ulong p)
     418             : {
     419     2164000 :   GEN u = cgetg(li+1,t_VECSMALL);
     420     2163996 :   ulong m = uel(b,li) % p;
     421             :   long i,j;
     422             : 
     423     2163996 :   uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
     424     5726429 :   for (i=li-1; i>0; i--)
     425             :   {
     426     3562433 :     m = b[i]%p;
     427     9117954 :     for (j = i+1; j <= li; j++)
     428     5555526 :       m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
     429     3562428 :     if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
     430     3562427 :     uel(u,i) = m;
     431             :   }
     432     2163996 :   return u;
     433             : }
     434             : 
     435             : static GEN
     436      667199 : 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      667199 :   n = lg(x)-1;
     443      667199 :   m=nbrows(x); r=0;
     444             : 
     445      667199 :   c = zero_zv(m);
     446      667198 :   d = cgetg(n+1, t_VECSMALL);
     447      667196 :   a = 0; /* for gcc -Wall */
     448     3277108 :   for (k=1; k<=n; k++)
     449             :   {
     450     9459337 :     for (j=1; j<=m; j++)
     451     8271505 :       if (!c[j])
     452             :       {
     453     4441670 :         a = ucoeff(x,j,k) % p;
     454     4441670 :         if (a) break;
     455             :       }
     456     2663795 :     if (j > m)
     457             :     {
     458     1187842 :       if (deplin==1) {
     459       53889 :         c = cgetg(n+1, t_VECSMALL);
     460      176160 :         for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
     461      105209 :         c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
     462       53889 :         return c;
     463             :       }
     464     1133953 :       r++; d[k] = 0;
     465             :     }
     466             :     else
     467             :     {
     468     1475953 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     469     1475959 :       c[j] = k; d[k] = j;
     470     1475959 :       ucoeff(x,j,k) = p-1;
     471     1475959 :       if (piv != 1)
     472     3181310 :         for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
     473     7590756 :       for (t=1; t<=m; t++)
     474             :       {
     475     6114797 :         if (t == j) continue;
     476             : 
     477     4638840 :         piv = ( ucoeff(x,t,k) %= p );
     478     4638840 :         if (!piv) continue;
     479     2169294 :         if (piv == 1)
     480     2548653 :           for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
     481             :         else
     482     4888275 :           for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
     483             :       }
     484             :     }
     485             :   }
     486      613313 :   if (deplin==1) return NULL;
     487             : 
     488      613306 :   y = cgetg(r+1, t_MAT);
     489     1747258 :   for (j=k=1; j<=r; j++,k++)
     490             :   {
     491     1133952 :     GEN C = cgetg(n+1, t_VECSMALL);
     492             : 
     493     2202345 :     gel(y,j) = C; while (d[k]) k++;
     494     9357855 :     for (i=1; i<k; i++)
     495     8223903 :       if (d[i])
     496     2337213 :         uel(C,i) = ucoeff(x,d[i],k) % p;
     497             :       else
     498     5886690 :         uel(C,i) = 0UL;
     499     7745571 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     500             :   }
     501      613306 :   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      613306 :   return y;
     508             : }
     509             : 
     510             : /* in place, destroy x */
     511             : static GEN
     512      814518 : 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      814518 :   n = lg(x)-1;
     518      814518 :   if (!n) return cgetg(1,t_MAT);
     519      814511 :   if (SMALL_ULONG(p)) return Flm_ker_gauss_OK(x, p, deplin);
     520      147313 :   pi = get_Fl_red(p);
     521             : 
     522      147312 :   m=nbrows(x); r=0;
     523             : 
     524      147312 :   c = zero_zv(m);
     525      147312 :   d = cgetg(n+1, t_VECSMALL);
     526      147312 :   a = 0; /* for gcc -Wall */
     527      474560 :   for (k=1; k<=n; k++)
     528             :   {
     529      854525 :     for (j=1; j<=m; j++)
     530      751176 :       if (!c[j])
     531             :       {
     532      567251 :         a = ucoeff(x,j,k);
     533      567251 :         if (a) break;
     534             :       }
     535      327255 :     if (j > m)
     536             :     {
     537      103349 :       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      103342 :       r++; d[k] = 0;
     544             :     }
     545             :     else
     546             :     {
     547      223906 :       ulong piv = p - Fl_inv(a, p); /* -1/a */
     548      223906 :       c[j] = k; d[k] = j;
     549      223906 :       ucoeff(x,j,k) = p-1;
     550      223906 :       if (piv != 1)
     551      405220 :         for (i=k+1; i<=n; i++)
     552      187237 :           ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
     553      971929 :       for (t=1; t<=m; t++)
     554             :       {
     555      748023 :         if (t == j) continue;
     556             : 
     557      524117 :         piv = ucoeff(x,t,k);
     558      524117 :         if (!piv) continue;
     559      308616 :         if (piv == 1)
     560       46124 :           for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
     561             :         else
     562      710347 :           for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
     563             :       }
     564             :     }
     565             :   }
     566      147305 :   if (deplin==1) return NULL;
     567             : 
     568      147298 :   y = cgetg(r+1, t_MAT);
     569      250640 :   for (j=k=1; j<=r; j++,k++)
     570             :   {
     571      103342 :     GEN C = cgetg(n+1, t_VECSMALL);
     572             : 
     573      166193 :     gel(y,j) = C; while (d[k]) k++;
     574      218087 :     for (i=1; i<k; i++)
     575      114745 :       if (d[i])
     576       79019 :         uel(C,i) = ucoeff(x,d[i],k);
     577             :       else
     578       35726 :         uel(C,i) = 0UL;
     579      166977 :     uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
     580             :   }
     581      147298 :   if (deplin == 2) {
     582      142444 :     GEN pc = cgetg(n - r + 1, t_VECSMALL);  /* indices of pivot columns */
     583      442500 :     for (i = j = 1; j <= n; j++)
     584      300056 :       if (d[j]) pc[i++] = j;
     585      142444 :     return mkvec2(y, pc);
     586             :   }
     587        4854 :   return y;
     588             : }
     589             : 
     590             : GEN
     591      271604 : Flm_intersect_i(GEN x, GEN y, ulong p)
     592             : {
     593      271604 :   long j, lx = lg(x);
     594             :   GEN z;
     595             : 
     596      271604 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
     597      271590 :   z = Flm_ker(shallowconcat(x,y), p);
     598     1470136 :   for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
     599      271588 :   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      319996 : Flm_ker_echelon(GEN x, ulong p, long pivots) {
     610      319996 :   pari_sp av = avma;
     611             :   GEN R, Rc, C, C1, C2, S, K;
     612      319996 :   long n = lg(x) - 1, r;
     613      319996 :   ulong pi = get_Fl_red(p);
     614      319996 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     615      319996 :   Rc = indexcompl(R, n);
     616      319997 :   C1 = rowpermute(C, R);
     617      319996 :   C2 = rowpermute(C, Rc);
     618      319996 :   S = Flm_lsolve_lower_unit_pre(C1, C2, p, pi);
     619      319997 :   K = vecpermute(shallowconcat(Flm_neg(S, p), matid_Flm(n - r)),
     620             :                  perm_inv(vecsmall_concat(R, Rc)));
     621      319997 :   K = Flm_transpose(K);
     622      319997 :   if (pivots)
     623       32219 :     return gerepilecopy(av, mkvec2(K, R));
     624      287778 :   return gerepileupto(av, K);
     625             : }
     626             : 
     627             : static GEN
     628       30803 : Flm_deplin_echelon(GEN x, ulong p) {
     629       30803 :   pari_sp av = avma;
     630             :   GEN R, Rc, C, C1, C2, s, v;
     631       30803 :   long i, n = lg(x) - 1, r;
     632       30803 :   ulong pi = get_Fl_red(p);
     633       30803 :   r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
     634       30803 :   if (r == n) return gc_NULL(av);
     635       30796 :   Rc = indexcompl(R, n);
     636       30796 :   i = Rc[1];
     637       30796 :   C1 = rowpermute(C, R);
     638       30796 :   C2 = rowslice(C, i, i);
     639       30796 :   s = Flm_row(Flm_lsolve_lower_unit_pre(C1, C2, p, pi), 1);
     640       30796 :   v = vecsmallpermute(vecsmall_concat(Flv_neg(s, p), vecsmall_ei(n - r, 1)),
     641             :                  perm_inv(vecsmall_concat(R, Rc)));
     642       30796 :   return gerepileuptoleaf(av, v);
     643             : }
     644             : 
     645             : static GEN
     646     1165319 : Flm_ker_i(GEN x, ulong p, long deplin, long inplace) {
     647     1165319 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     648      350799 :     switch(deplin) {
     649      287777 :     case 0: return Flm_ker_echelon(x, p, 0);
     650       30803 :     case 1: return Flm_deplin_echelon(x, p);
     651       32219 :     case 2: return Flm_ker_echelon(x, p, 1);
     652             :     }
     653      814520 :   return Flm_ker_gauss(inplace? x: Flm_copy(x), p, deplin);
     654             : }
     655             : 
     656             : GEN
     657      615278 : Flm_ker_sp(GEN x, ulong p, long deplin) { return Flm_ker_i(x, p, deplin, 1); }
     658             : GEN
     659      550041 : 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        2310 : Flm_det_gauss_OK(GEN a, long nbco, ulong p)
     666             : {
     667        2310 :   long i,j,k, s = 1;
     668        2310 :   ulong q, x = 1;
     669             : 
     670        9534 :   for (i=1; i<nbco; i++)
     671             :   {
     672        8841 :     for(k=i; k<=nbco; k++)
     673             :     {
     674        8659 :       ulong c = ucoeff(a,k,i) % p;
     675        8659 :       ucoeff(a,k,i) = c;
     676        8659 :       if (c) break;
     677             :     }
     678       23191 :     for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
     679        7406 :     if (k > nbco) return ucoeff(a,i,i);
     680        7224 :     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        7224 :     q = ucoeff(a,i,i);
     686             : 
     687        7224 :     if (x & HIGHMASK) x %= p;
     688        7224 :     x *= q;
     689        7224 :     q = Fl_inv(q,p);
     690       24080 :     for (k=i+1; k<=nbco; k++)
     691             :     {
     692       16856 :       ulong m = ucoeff(a,i,k) % p;
     693       16856 :       if (!m) continue;
     694             : 
     695        9618 :       m = p - ((m*q)%p);
     696       37576 :       for (j=i+1; j<=nbco; j++)
     697             :       {
     698       27958 :         ulong c = ucoeff(a,j,k);
     699       27958 :         if (c & HIGHMASK) c %= p;
     700       27958 :         ucoeff(a,j,k) = c  + m*ucoeff(a,j,i);
     701             :       }
     702             :     }
     703             :   }
     704        2128 :   if (x & HIGHMASK) x %= p;
     705        2128 :   q = ucoeff(a,nbco,nbco);
     706        2128 :   if (q & HIGHMASK) q %= p;
     707        2128 :   x = (x*q) % p;
     708        2128 :   if (s < 0 && x) x = p - x;
     709        2128 :   return x;
     710             : }
     711             : 
     712             : /* in place, destroy a */
     713             : static ulong
     714       53345 : Flm_det_gauss(GEN a, ulong p)
     715             : {
     716       53345 :   long i,j,k, s = 1, nbco = lg(a)-1;
     717       53345 :   ulong pi, q, x = 1;
     718             : 
     719       53345 :   if (SMALL_ULONG(p)) return Flm_det_gauss_OK(a, nbco, p);
     720       51035 :   pi = get_Fl_red(p);
     721      325543 :   for (i=1; i<nbco; i++)
     722             :   {
     723      274911 :     for(k=i; k<=nbco; k++)
     724      274911 :       if (ucoeff(a,k,i)) break;
     725      274508 :     if (k > nbco) return ucoeff(a,i,i);
     726      274508 :     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      274508 :     q = ucoeff(a,i,i);
     732             : 
     733      274508 :     x = Fl_mul_pre(x, q, p, pi);
     734      274508 :     q = Fl_inv(q,p);
     735     1176224 :     for (k=i+1; k<=nbco; k++)
     736             :     {
     737      901716 :       ulong m = ucoeff(a,i,k);
     738      901716 :       if (!m) continue;
     739             : 
     740      842414 :       m = Fl_mul_pre(m, q, p, pi);
     741     4329292 :       for (j=i+1; j<=nbco; j++)
     742     3486878 :         ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
     743             :     }
     744             :   }
     745       51035 :   if (s < 0) x = Fl_neg(x, p);
     746       51035 :   return Fl_mul(x, ucoeff(a,nbco,nbco), p);
     747             : }
     748             : 
     749             : static ulong
     750       37377 : Flm_det_CUP(GEN a, ulong p) {
     751             :   GEN R, C, U, P;
     752       37377 :   long i, n = lg(a) - 1, r;
     753             :   ulong d;
     754       37377 :   ulong pi = get_Fl_red(p);
     755       37377 :   r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi);
     756       37378 :   if (r < n)
     757          42 :     d = 0;
     758             :   else {
     759       37336 :     d = perm_sign(P) == 1? 1: p-1;
     760      429910 :     for (i = 1; i <= n; i++)
     761      392574 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
     762             :   }
     763       37378 :   return d;
     764             : }
     765             : 
     766             : static ulong
     767       90722 : Flm_det_i(GEN x, ulong p, long inplace) {
     768       90722 :   pari_sp av = avma;
     769             :   ulong d;
     770       90722 :   if (lg(x) - 1 >= Flm_CUP_LIMIT)
     771       37377 :     d = Flm_det_CUP(x, p);
     772             :   else
     773       53345 :     d = Flm_det_gauss(inplace? x: Flm_copy(x), p);
     774       90723 :   return gc_ulong(av, d);
     775             : }
     776             : 
     777             : ulong
     778       90722 : 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     4397340 : 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     4397340 :   n=lg(x)-1; if (!n) { *rr=0; return NULL; }
     790             : 
     791     4396031 :   m=nbrows(x); r=0;
     792     4396029 :   d=cgetg(n+1,t_VECSMALL);
     793     4396028 :   c = zero_zv(m);
     794    18907579 :   for (k=1; k<=n; k++)
     795             :   {
     796    36108527 :     for (j=1; j<=m; j++)
     797    34308059 :       if (!c[j])
     798             :       {
     799    18418920 :         ucoeff(x,j,k) %= p;
     800    18418920 :         if (ucoeff(x,j,k)) break;
     801             :       }
     802    14511522 :     if (j>m) { r++; d[k]=0; }
     803             :     else
     804             :     {
     805    12711020 :       ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
     806    12711085 :       c[j]=k; d[k]=j;
     807    29213769 :       for (i=k+1; i<=n; i++)
     808    16502686 :         ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
     809    60432498 :       for (t=1; t<=m; t++)
     810    47721533 :         if (!c[t]) /* no pivot on that line yet */
     811             :         {
     812    20430948 :           piv = ucoeff(x,t,k);
     813    20430948 :           if (piv)
     814             :           {
     815    12029107 :             ucoeff(x,t,k) = 0;
     816    37274104 :             for (i=k+1; i<=n; i++)
     817    25245108 :               ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
     818    25245115 :                                      Fl_mul(piv,ucoeff(x,j,i),p),p);
     819             :           }
     820             :         }
     821    41924658 :       for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
     822             :     }
     823             :   }
     824     4396057 :   *rr = r; return gc_const((pari_sp)d, d);
     825             : }
     826             : 
     827             : static GEN
     828      274039 : Flm_pivots_CUP(GEN x, ulong p, long *rr)
     829             : {
     830      274039 :   long i, n = lg(x) - 1, r;
     831      274039 :   GEN R, C, U, P, d = zero_zv(n);
     832      274039 :   ulong pi = get_Fl_red(p);
     833      274039 :   r = Flm_CUP_pre(x, &R, &C, &U, &P, p, pi);
     834     3018406 :   for(i = 1; i <= r; i++)
     835     2744367 :     d[P[i]] = R[i];
     836      274039 :   *rr = n - r; return gc_const((pari_sp)d, d);
     837             : }
     838             : 
     839             : GEN
     840     4671379 : Flm_pivots(GEN x, ulong p, long *rr, long inplace)
     841             : {
     842     4671379 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
     843      274039 :     return Flm_pivots_CUP(x, p, rr);
     844     4397340 :   return Flm_gauss_pivot(inplace? x: Flm_copy(x), p, rr);
     845             : }
     846             : 
     847             : long
     848       39844 : Flm_rank(GEN x, ulong p)
     849             : {
     850       39844 :   pari_sp av = avma;
     851             :   long r;
     852       39844 :   if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT) {
     853             :     GEN R, C;
     854        9688 :     ulong pi = get_Fl_red(p);
     855        9688 :     return gc_long(av, Flm_echelon_pre(x, &R, &C, p, pi));
     856             :   }
     857       30156 :   (void) Flm_pivots(x, p, &r, 0);
     858       30156 :   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        1274 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
     865             : {
     866        1274 :   long n = lg(A)-1, i = index, j;
     867        1274 :   GEN u = const_vecsmall(n, 0);
     868        1274 :   u[i] = 1;
     869        1274 :   if (SMALL_ULONG(p))
     870        3269 :     for (i--; i>0; i--)
     871             :     {
     872        2009 :       ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
     873        7035 :       for (j=i+2; j<=n; j++)
     874             :       {
     875        5026 :         if (m & HIGHMASK) m %= p;
     876        5026 :         m += ucoeff(A,i,j) * uel(u,j);
     877             :       }
     878        2009 :       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        1274 :   return u;
     888             : }
     889             : static GEN
     890         490 : Flm_inv_upper_1(GEN A, ulong p)
     891             : {
     892             :   long i, l;
     893         490 :   GEN B = cgetg_copy(A, &l);
     894        1764 :   for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
     895         490 :   return B;
     896             : }
     897             : /* destroy a, b */
     898             : static GEN
     899      118581 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
     900             : {
     901      118581 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     902      118581 :   ulong det = 1;
     903             :   GEN u;
     904             : 
     905      118581 :   li = nbrows(a);
     906      118581 :   bco = lg(b)-1;
     907      416974 :   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     1028305 :     for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
     912      473705 :     for (k = i; k <= li; k++)
     913             :     {
     914      473698 :       ulong piv = ( ucoeff(a,k,i) %= p );
     915      473698 :       if (piv)
     916             :       {
     917      416967 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     918      416967 :         if (detp)
     919             :         {
     920           0 :           if (det & HIGHMASK) det %= p;
     921           0 :           det *= piv;
     922             :         }
     923      416967 :         break;
     924             :       }
     925             :     }
     926             :     /* found a pivot on line k */
     927      416974 :     if (k > li) return NULL;
     928      416967 :     if (k != i)
     929             :     { /* swap lines so that k = i */
     930       38059 :       s = -s;
     931      157332 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     932      206065 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     933             :     }
     934      416967 :     if (i == aco) break;
     935             : 
     936      298393 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
     937      947153 :     for (k=i+1; k<=li; k++)
     938             :     {
     939      648760 :       ulong m = ( ucoeff(a,k,i) %= p );
     940      648760 :       if (!m) continue;
     941             : 
     942      183353 :       m = Fl_mul(m, invpiv, p);
     943      183353 :       if (m == 1) {
     944      142397 :         for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
     945      225202 :         for (j=1;   j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
     946             :       } else {
     947      509175 :         for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
     948      845658 :         for (j=1;   j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
     949             :       }
     950             :     }
     951             :   }
     952      118574 :   if (detp)
     953             :   {
     954           0 :     det %=  p;
     955           0 :     if (s < 0 && det) det = p - det;
     956           0 :     *detp = det;
     957             :   }
     958      118574 :   u = cgetg(bco+1,t_MAT);
     959      558907 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
     960      118574 :   return u;
     961             : }
     962             : 
     963             : /* destroy a, b */
     964             : static GEN
     965     2157980 : Flm_gauss_sp_i(GEN a, GEN b, ulong *detp, ulong p)
     966             : {
     967     2157980 :   long i, j, k, li, bco, aco = lg(a)-1, s = 1;
     968     2157980 :   ulong det = 1;
     969             :   GEN u;
     970             :   ulong pi;
     971     2157980 :   if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
     972     2157980 :   if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
     973     2039399 :   pi = get_Fl_red(p);
     974     2039400 :   li = nbrows(a);
     975     2039399 :   bco = lg(b)-1;
     976     5308061 :   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     5547214 :     for (k = i; k <= li; k++)
     981             :     {
     982     5547214 :       ulong piv = ucoeff(a,k,i);
     983     5547214 :       if (piv)
     984             :       {
     985     5308035 :         ucoeff(a,k,i) = Fl_inv(piv, p);
     986     5308038 :         if (detp) det = Fl_mul_pre(det, piv, p, pi);
     987     5308036 :         break;
     988             :       }
     989             :     }
     990             :     /* found a pivot on line k */
     991     5308036 :     if (k > li) return NULL;
     992     5308036 :     if (k != i)
     993             :     { /* swap lines so that k = i */
     994      182115 :       s = -s;
     995      678000 :       for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
     996      369965 :       for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
     997             :     }
     998     5308036 :     if (i == aco) break;
     999             : 
    1000     3268631 :     invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
    1001     8247634 :     for (k=i+1; k<=li; k++)
    1002             :     {
    1003     4978972 :       ulong m = ucoeff(a,k,i);
    1004     4978972 :       if (!m) continue;
    1005             : 
    1006     4141229 :       m = Fl_mul_pre(m, invpiv, p, pi);
    1007     4141223 :       if (m == 1) {
    1008      637894 :         for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
    1009      443503 :         for (j=1;   j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
    1010             :       } else {
    1011    11369912 :         for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
    1012     7895955 :         for (j=1;   j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
    1013             :       }
    1014             :     }
    1015             :   }
    1016     2039431 :   if (detp)
    1017             :   {
    1018           0 :     if (s < 0 && det) det = p - det;
    1019           0 :     *detp = det;
    1020             :   }
    1021     2039431 :   u = cgetg(bco+1,t_MAT);
    1022     4203400 :   for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
    1023     2039400 :   return u;
    1024             : }
    1025             : 
    1026             : static GEN
    1027     1085637 : Flm_gauss_from_CUP(GEN b, GEN R, GEN C, GEN U, GEN P, ulong p, ulong pi, ulong *detp)
    1028             : {
    1029     1085637 :   GEN Y = Flm_rsolve_lower_unit_pre(rowpermute(C, R), rowpermute(b, R), p, pi);
    1030     1085644 :   GEN X = rowpermute(Flm_rsolve_upper_pre(U, Y, p, pi), perm_inv(P));
    1031     1085642 :   if (detp)
    1032             :   {
    1033      918379 :     ulong d = perm_sign(P) == 1? 1: p-1;
    1034      918377 :     long i, r = lg(R);
    1035     5671375 :     for (i = 1; i < r; i++)
    1036     4752994 :       d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
    1037      918381 :     *detp = d;
    1038             :   }
    1039     1085644 :   return X;
    1040             : }
    1041             : 
    1042             : static GEN
    1043      167278 : Flm_gauss_CUP(GEN a, GEN b, ulong *detp, ulong p) {
    1044             :   GEN R, C, U, P;
    1045      167278 :   long n = lg(a) - 1, r;
    1046      167278 :   ulong pi = get_Fl_red(p);
    1047      167278 :   if (nbrows(a) < n || (r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi)) < n)
    1048          14 :     return NULL;
    1049      167264 :   return Flm_gauss_from_CUP(b, R, C, U, P, p, pi, detp);
    1050             : }
    1051             : 
    1052             : GEN
    1053     2205166 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p) {
    1054     2205166 :   pari_sp av = avma;
    1055             :   GEN x;
    1056     2205166 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1057       47186 :     x = Flm_gauss_CUP(a, b, detp, p);
    1058             :   else
    1059     2157980 :     x = Flm_gauss_sp_i(a, b, detp, p);
    1060     2205165 :   if (!x) return gc_NULL(av);
    1061     2205158 :   return gerepileupto(av, x);
    1062             : }
    1063             : 
    1064             : GEN
    1065     2143751 : Flm_gauss(GEN a, GEN b, ulong p) {
    1066     2143751 :   pari_sp av = avma;
    1067             :   GEN x;
    1068     2143751 :   if (lg(a) - 1 >= Flm_CUP_LIMIT)
    1069      104360 :     x = Flm_gauss_CUP(a, b, NULL, p);
    1070             :   else {
    1071     2039391 :     a = RgM_shallowcopy(a);
    1072     2039390 :     b = RgM_shallowcopy(b);
    1073     2039390 :     x = Flm_gauss_sp(a, b, NULL, p);
    1074             :   }
    1075     2143749 :   if (!x) return gc_NULL(av);
    1076     2143735 :   return gerepileupto(av, x);
    1077             : }
    1078             : 
    1079             : static GEN
    1080       63628 : Flm_inv_i(GEN a, ulong *detp, ulong p, long inplace) {
    1081       63628 :   pari_sp av = avma;
    1082       63628 :   long n = lg(a) - 1;
    1083             :   GEN b, x;
    1084       63628 :   if (n == 0) return cgetg(1, t_MAT);
    1085       63628 :   b = matid_Flm(nbrows(a));
    1086       63628 :   if (n >= Flm_CUP_LIMIT)
    1087       15732 :     x = Flm_gauss_CUP(a, b, detp, p);
    1088             :   else {
    1089       47896 :     if (!inplace)
    1090       46300 :       a = RgM_shallowcopy(a);
    1091       47896 :     x = Flm_gauss_sp(a, b, detp, p);
    1092             :   }
    1093       63628 :   if (!x) return gc_NULL(av);
    1094       63621 :   return gerepileupto(av, x);
    1095             : }
    1096             : 
    1097             : GEN
    1098        1827 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
    1099        1827 :   return Flm_inv_i(a, detp, p, 1);
    1100             : }
    1101             : 
    1102             : GEN
    1103       61801 : Flm_inv(GEN a, ulong p) {
    1104       61801 :   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      918391 : Flm_adjoint(GEN A, ulong p)
    1118             : {
    1119      918391 :   pari_sp av = avma;
    1120             :   GEN R, C, U, P, C1, U1, v, c, d;
    1121      918391 :   long r, i, q, n = lg(A)-1, m;
    1122             :   ulong D;
    1123      918391 :   ulong pi = get_Fl_red(p);
    1124      918389 :   if (n == 0) return cgetg(1, t_MAT);
    1125      918389 :   r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
    1126      918387 :   m = nbrows(A);
    1127      918387 :   if (r == n)
    1128             :   {
    1129      918373 :     GEN X = Flm_gauss_from_CUP(matid_Flm(m), R, C, U, P, p, pi, &D);
    1130      918380 :     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         931 : Flm_invimage_i(GEN A, GEN B, ulong p)
    1175             : {
    1176             :   GEN d, x, X, Y;
    1177         931 :   long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
    1178             : 
    1179         931 :   if (!nB) return cgetg(1, t_MAT);
    1180         784 :   if (maxss(nA, nB) >= Flm_CUP_LIMIT && nbrows(B) >= Flm_CUP_LIMIT)
    1181         287 :     return Flm_invimage_CUP(A, B, p);
    1182             : 
    1183         497 :   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         497 :   nY = lg(x)-1;
    1188         497 :   if (nY < nB) return NULL;
    1189         497 :   Y = rowslice(x, nA+1, nA+nB); /* nB rows */
    1190         497 :   d = cgetg(nB+1, t_VECSMALL);
    1191        1778 :   for (i = nB, j = nY; i >= 1; i--, j--)
    1192             :   {
    1193        1295 :     for (; j>=1; j--)
    1194        1288 :       if (coeff(Y,i,j)) { d[i] = j; break; }
    1195        1288 :     if (!j) return NULL;
    1196             :   }
    1197             :   /* reduce to the case Y square, upper triangular with 1s on diagonal */
    1198         490 :   Y = vecpermute(Y, d);
    1199         490 :   x = vecpermute(x, d);
    1200         490 :   X = rowslice(x, 1, nA);
    1201         490 :   return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
    1202             : }
    1203             : GEN
    1204         889 : Flm_invimage(GEN A, GEN B, ulong p)
    1205             : {
    1206         889 :   pari_sp av = avma;
    1207         889 :   GEN X = Flm_invimage_i(A,B,p);
    1208         889 :   if (!X) return gc_NULL(av);
    1209         889 :   return gerepileupto(av, X);
    1210             : }
    1211             : 
    1212             : GEN
    1213      131611 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
    1214             : {
    1215      131611 :   pari_sp av = avma;
    1216      131611 :   long i, l = lg(A);
    1217             :   GEN M, x;
    1218             :   ulong t;
    1219             : 
    1220      131611 :   if (l==1) return NULL;
    1221      131471 :   if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
    1222      131471 :   M = cgetg(l+1,t_MAT);
    1223     1080326 :   for (i=1; i<l; i++) gel(M,i) = gel(A,i);
    1224      131471 :   gel(M,l) = y; M = Flm_ker(M,p);
    1225      131472 :   i = lg(M)-1; if (!i) return gc_NULL(av);
    1226             : 
    1227      128239 :   x = gel(M,i); t = x[l];
    1228      128239 :   if (!t) return gc_NULL(av);
    1229             : 
    1230      128232 :   setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
    1231      128231 :   if (t!=1) x = Flv_Fl_mul(x, t, p);
    1232      128232 :   return gerepileuptoleaf(av, x);
    1233             : }

Generated by: LCOV version 1.16