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 - FpV.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30256-57c6c2b57f) Lines: 743 970 76.6 %
Date: 2025-05-09 09:19:30 Functions: 110 138 79.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  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             : #define DEBUGLEVEL DEBUGLEVEL_mat
      19             : 
      20             : /********************************************************************/
      21             : /**                                                                **/
      22             : /**                           REDUCTION                            **/
      23             : /**                                                                **/
      24             : /********************************************************************/
      25             : /* z in Z^n, return lift(Col(z) * Mod(1,p)) */
      26             : GEN
      27    29828819 : FpC_red(GEN x, GEN p)
      28   261375312 : { pari_APPLY_type(t_COL, modii(gel(x,i), p)) }
      29             : 
      30             : /* z in Z^n, return lift(Vec(z) * Mod(1,p)) */
      31             : GEN
      32         420 : FpV_red(GEN x, GEN p)
      33        2695 : { pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
      34             : 
      35             : GEN
      36     6863308 : FpC_center(GEN x, GEN p, GEN pov2)
      37    43259782 : { pari_APPLY_type(t_COL, Fp_center(gel(x,i), p, pov2)) }
      38             : 
      39             : GEN
      40      164710 : Flv_center(GEN x, ulong p, ulong ps2)
      41     2621724 : { pari_APPLY_ulong(Fl_center(uel(x,i),p,ps2)) }
      42             : 
      43             : /* z in Mat m,n(Z), return lift(z * Mod(1,p)) */
      44             : GEN
      45     5193365 : FpM_red(GEN x, GEN p)
      46    24397312 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
      47             : 
      48             : GEN
      49     2728228 : FpM_center(GEN x, GEN p, GEN pov2)
      50     8271567 : { pari_APPLY_same(FpC_center(gel(x,i), p,pov2)) }
      51             : 
      52             : /* HACK: assume p > 3, which ensures u = gen_[0-2] is never written to */
      53             : static void
      54         168 : Fp_center_inplace(GEN u, GEN p, GEN ps2)
      55             : {
      56         168 :   if (abscmpii(u, ps2) > 0)
      57             :   {
      58          28 :     pari_sp av = avma;
      59          28 :     affii(subii(u, p), u);
      60          28 :     set_avma(av);
      61             :   }
      62         168 : }
      63             : 
      64             : /* p > 3; assume entries in [0,p[ and ps2 = p>>1. */
      65             : static void
      66          42 : _FpC_center_inplace(GEN z, GEN p, GEN ps2)
      67             : {
      68          42 :   long i, l = lg(z);
      69         210 :   for (i = 1; i < l; i++) Fp_center_inplace(gel(z,i), p, ps2);
      70          42 : }
      71             : static void
      72           7 : _F3C_center_inplace(GEN z)
      73             : {
      74           7 :   long i, l = lg(z);
      75          35 :   for (i = 1; i < l; i++) /* z[i] = 0, 1 : no-op */
      76          28 :     if (equaliu(gel(z,i), 2)) gel(z,i) = gen_m1;
      77           7 : }
      78             : void
      79           0 : FpC_center_inplace(GEN z, GEN p, GEN ps2)
      80             : {
      81           0 :   switch (itou_or_0(p))
      82             :   {
      83           0 :     case 2: break;
      84           0 :     case 3:_F3C_center_inplace(z); break;
      85           0 :     default: _FpC_center_inplace(z, p, ps2);
      86             :   }
      87           0 : }
      88             : 
      89             : void
      90          42 : FpM_center_inplace(GEN z, GEN p, GEN pov2)
      91             : {
      92          42 :   long i, l = lg(z);
      93          42 :   switch (itou_or_0(p))
      94             :   {
      95          21 :     case 2: break;
      96           7 :     case 3:
      97          14 :       for (i = 1; i < l; i++) _F3C_center_inplace(gel(z,i));
      98           7 :       break;
      99          14 :     default:
     100          56 :       for (i = 1; i < l; i++) _FpC_center_inplace(gel(z,i), p, pov2);
     101             :   }
     102          42 : }
     103             : GEN
     104       11347 : Flm_center(GEN x, ulong p, ulong ps2)
     105      175875 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
     106             : 
     107             : GEN
     108         147 : random_FpV(long d, GEN p)
     109             : {
     110             :   long i;
     111         147 :   GEN y = cgetg(d+1,t_VEC);
     112       23191 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
     113         147 :   return y;
     114             : }
     115             : 
     116             : GEN
     117         924 : random_FpC(long d, GEN p)
     118             : {
     119             :   long i;
     120         924 :   GEN y = cgetg(d+1,t_COL);
     121        6188 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
     122         924 :   return y;
     123             : }
     124             : 
     125             : GEN
     126       41406 : random_Flv(long n, ulong p)
     127             : {
     128       41406 :   GEN y = cgetg(n+1, t_VECSMALL);
     129             :   long i;
     130      262114 :   for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
     131       41406 :   return y;
     132             : }
     133             : 
     134             : /********************************************************************/
     135             : /**                                                                **/
     136             : /**                           ADD, SUB                             **/
     137             : /**                                                                **/
     138             : /********************************************************************/
     139             : GEN
     140      283601 : FpC_add(GEN x, GEN y, GEN p)
     141     3551393 : { pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
     142             : 
     143             : GEN
     144           0 : FpV_add(GEN x, GEN y, GEN p)
     145           0 : { pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
     146             : 
     147             : GEN
     148           0 : FpM_add(GEN x, GEN y, GEN p)
     149           0 : { pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
     150             : 
     151             : GEN
     152     1812467 : Flv_add(GEN x, GEN y, ulong p)
     153             : {
     154     1812467 :   if (p==2)
     155     7096824 :     pari_APPLY_ulong(x[i]^y[i])
     156             :   else
     157    41998884 :     pari_APPLY_ulong(Fl_add(x[i], y[i], p))
     158             : }
     159             : 
     160             : void
     161      446634 : Flv_add_inplace(GEN x, GEN y, ulong p)
     162             : {
     163      446634 :   long i, l = lg(x);
     164      446634 :   if (p==2)
     165     1388909 :     for (i = 1; i < l; i++) x[i] ^= y[i];
     166             :   else
     167      164878 :     for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
     168      446634 : }
     169             : 
     170             : ulong
     171        3871 : Flv_sum(GEN x, ulong p)
     172             : {
     173        3871 :   long i, l = lg(x);
     174        3871 :   ulong s = 0;
     175        3871 :   if (p==2)
     176       17689 :     for (i = 1; i < l; i++) s ^= x[i];
     177             :   else
     178           0 :     for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
     179        3871 :   return s;
     180             : }
     181             : 
     182             : GEN
     183      965634 : FpC_sub(GEN x, GEN y, GEN p)
     184    18959611 : { pari_APPLY_type(t_COL, Fp_sub(gel(x,i), gel(y,i), p)) }
     185             : 
     186             : GEN
     187           0 : FpV_sub(GEN x, GEN y, GEN p)
     188           0 : { pari_APPLY_type(t_VEC, Fp_sub(gel(x,i), gel(y,i), p)) }
     189             : 
     190             : GEN
     191           0 : FpM_sub(GEN x, GEN y, GEN p)
     192           0 : { pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
     193             : 
     194             : GEN
     195   230929718 : Flv_sub(GEN x, GEN y, ulong p)
     196  1033965989 : { pari_APPLY_ulong(Fl_sub(uel(x,i), uel(y,i), p)) }
     197             : 
     198             : void
     199           0 : Flv_sub_inplace(GEN x, GEN y, ulong p)
     200             : {
     201           0 :   long i, l = lg(x);
     202           0 :   for (i = 1; i < l; i++) x[i] = Fl_sub(x[i], y[i], p);
     203           0 : }
     204             : 
     205             : GEN
     206       51191 : Flm_Fl_add(GEN x, ulong y, ulong p)
     207             : {
     208       51191 :   long l = lg(x), i, j;
     209       51191 :   GEN z = cgetg(l,t_MAT);
     210             : 
     211       51191 :   if (l==1) return z;
     212       51191 :   if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
     213      241292 :   for (i=1; i<l; i++)
     214             :   {
     215      190097 :     GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
     216      190099 :     gel(z,i) = zi;
     217     1331184 :     for (j=1; j<l; j++) zi[j] = xi[j];
     218      190099 :     zi[i] = Fl_add(zi[i], y, p);
     219             :   }
     220       51195 :   return z;
     221             : }
     222             : GEN
     223        3724 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
     224             : 
     225             : GEN
     226       29322 : Flm_add(GEN x, GEN y, ulong p)
     227      695951 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
     228             : 
     229             : GEN
     230    25765870 : Flm_sub(GEN x, GEN y, ulong p)
     231   256673757 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
     232             : 
     233             : /********************************************************************/
     234             : /**                                                                **/
     235             : /**                           MULTIPLICATION                       **/
     236             : /**                                                                **/
     237             : /********************************************************************/
     238             : GEN
     239      960453 : FpC_Fp_mul(GEN x, GEN y, GEN p)
     240    11704199 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
     241             : 
     242             : GEN
     243     1753294 : Flv_Fl_mul(GEN x, ulong y, ulong p)
     244    41051905 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
     245             : 
     246             : GEN
     247        4692 : Flv_Fl_div(GEN x, ulong y, ulong p)
     248             : {
     249        4692 :   return Flv_Fl_mul(x, Fl_inv(y, p), p);
     250             : }
     251             : void
     252           0 : Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
     253             : {
     254           0 :   Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
     255           0 : }
     256             : GEN
     257      765461 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
     258      765461 :   long i, j, h, l = lg(X);
     259      765461 :   GEN A = cgetg(l, t_MAT);
     260      765461 :   if (l == 1) return A;
     261      765461 :   h = lgcols(X);
     262     4075248 :   for (j=1; j<l; j++)
     263             :   {
     264     3309921 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     265    24311447 :     for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
     266     3309787 :     gel(A,j) = a;
     267             :   }
     268      765327 :   return A;
     269             : }
     270             : 
     271             : /* x *= y */
     272             : void
     273     7776052 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
     274             : {
     275             :   long i;
     276     7776052 :   if (HIGHWORD(y | p))
     277     8598169 :     for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
     278             :   else
     279    32083067 :     for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
     280     7776049 : }
     281             : void
     282           0 : Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
     283           0 : { Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
     284             : 
     285             : /* set y *= x */
     286             : void
     287           0 : Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
     288             : {
     289           0 :   long i, j, m, l = lg(y);
     290           0 :   if (l == 1) return;
     291           0 :   m = lgcols(y);
     292           0 :   if (HIGHWORD(x | p))
     293           0 :     for(j=1; j<l; j++)
     294           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
     295             :   else
     296           0 :     for(j=1; j<l; j++)
     297           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
     298             : }
     299             : 
     300             : /* return x * y */
     301             : static GEN
     302    17150082 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
     303             : {
     304    17150082 :   long i, j, m, l = lg(y);
     305    17150082 :   GEN z = cgetg(l, t_MAT);
     306    17147967 :   if (l == 1) return z;
     307    17147967 :   m = lgcols(y);
     308   153568615 :   for(j=1; j<l; j++) {
     309   136329334 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     310   366409040 :     for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
     311             :   }
     312    17239281 :   return z;
     313             : }
     314             : 
     315             : /* return x * y */
     316             : static GEN
     317     8639804 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
     318             : {
     319     8639804 :   long i, j, m, l = lg(y);
     320     8639804 :   GEN z = cgetg(l, t_MAT);
     321     8639794 :   if (l == 1) return z;
     322     8639794 :   m = lgcols(y);
     323    53929572 :   for(j=1; j<l; j++) {
     324    45289813 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     325   163431533 :     for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
     326             :   }
     327     8639759 :   return z;
     328             : }
     329             : 
     330             : /* return x * y */
     331             : GEN
     332    25707295 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
     333             : {
     334    25707295 :   if (HIGHWORD(p))
     335    17150144 :     return Flm_Fl_mul_pre_i(y, x, p, pi);
     336             :   else
     337     8557151 :     return Flm_Fl_mul_OK(y, x, p);
     338             : }
     339             : 
     340             : /* return x * y */
     341             : GEN
     342       82649 : Flm_Fl_mul(GEN y, ulong x, ulong p)
     343             : {
     344       82649 :   if (HIGHWORD(p))
     345          74 :     return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
     346             :   else
     347       82575 :     return Flm_Fl_mul_OK(y, x, p);
     348             : }
     349             : 
     350             : GEN
     351     4102228 : Flv_neg(GEN x, ulong p)
     352    31299537 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
     353             : 
     354             : void
     355        6509 : Flv_neg_inplace(GEN v, ulong p)
     356             : {
     357             :   long i;
     358      184269 :   for (i = 1; i < lg(v); ++i)
     359      177760 :     v[i] = Fl_neg(v[i], p);
     360        6509 : }
     361             : 
     362             : GEN
     363      324925 : Flm_neg(GEN x, ulong p)
     364     4396407 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
     365             : 
     366             : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
     367             : static GEN
     368    18918660 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
     369             : {
     370    18918660 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1));
     371             :   long k;
     372   230080292 :   for (k = 2; k < lx; k++)
     373             :   {
     374   211162265 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     375   211159963 :     if (signe(t)) c = addii(c, t);
     376             :   }
     377    18918027 :   return c;
     378             : }
     379             : 
     380             : static long
     381    97688616 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
     382             : {
     383             :   long k;
     384    97688616 :   long c = coeff(x,i,1) * y[1];
     385  1498281691 :   for (k = 2; k < lx; k++)
     386  1400593075 :     c += coeff(x,i,k) * y[k];
     387    97688616 :   return c;
     388             : }
     389             : 
     390             : GEN
     391     6450565 : zm_zc_mul(GEN x, GEN y)
     392             : {
     393     6450565 :   long lx = lg(x), l, i;
     394             :   GEN z;
     395     6450565 :   if (lx == 1) return cgetg(1, t_VECSMALL);
     396     6450565 :   l = lg(gel(x,1));
     397     6450565 :   z = cgetg(l,t_VECSMALL);
     398   104139181 :   for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
     399     6450565 :   return z;
     400             : }
     401             : 
     402             : GEN
     403        2114 : zm_mul(GEN x, GEN y)
     404             : {
     405        2114 :   long i,j,lx=lg(x), ly=lg(y);
     406             :   GEN z;
     407        2114 :   if (ly==1) return cgetg(1,t_MAT);
     408        2114 :   z = cgetg(ly,t_MAT);
     409        2114 :   if (lx==1)
     410             :   {
     411           0 :     for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
     412           0 :     return z;
     413             :   }
     414       30849 :   for (j=1; j<ly; j++)
     415       28735 :     gel(z,j) = zm_zc_mul(x, gel(y,j));
     416        2114 :   return z;
     417             : }
     418             : 
     419             : static ulong
     420   768193079 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
     421             : {
     422   768193079 :   ulong c = ucoeff(x,i,1) * uel(y,1);
     423             :   long k;
     424 12011511564 :   for (k = 2; k < lx; k++) {
     425 11243318485 :     c += ucoeff(x,i,k) * uel(y,k);
     426 11243318485 :     if (c & HIGHBIT) c %= p;
     427             :   }
     428   768193079 :   return c % p;
     429             : }
     430             : 
     431             : static ulong
     432   532947140 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
     433             : {
     434             :   ulong l0, l1, v1, h0, h1;
     435   532947140 :   long k = 1;
     436             :   LOCAL_OVERFLOW;
     437             :   LOCAL_HIREMAINDER;
     438   532947140 :   l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
     439  8865873870 :   while (++k < lx) {
     440  8332926730 :     l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
     441  8332926730 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     442             :   }
     443   532947140 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     444    54755728 :   else return remlll_pre(v1, h1, l1, p, pi);
     445             : }
     446             : 
     447             : static GEN
     448      329854 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
     449             : {
     450             :   long i,j;
     451      329854 :   GEN z = NULL;
     452             : 
     453     3725914 :   for (j=1; j<lx; j++)
     454             :   {
     455     3396060 :     if (!y[j]) continue;
     456     1106411 :     if (!z) z = Flv_copy(gel(x,j));
     457    29798064 :     else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
     458             :   }
     459      329854 :   if (!z) z = zero_zv(l-1);
     460      329854 :   return z;
     461             : }
     462             : 
     463             : static GEN
     464     1754860 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
     465             : {
     466     1754860 :   GEN z = cgetg(l,t_COL);
     467             :   long i;
     468    16782122 :   for (i = 1; i < l; i++)
     469             :   {
     470    15027301 :     pari_sp av = avma;
     471    15027301 :     GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
     472    15027128 :     gel(z,i) = gc_INT(av, modii(c,p));
     473             :   }
     474     1754821 :   return z;
     475             : }
     476             : 
     477             : static void
     478    86301249 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
     479             : {
     480             :   long i;
     481   854526833 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
     482    86315951 : }
     483             : static GEN
     484    82948204 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
     485             : {
     486    82948204 :   GEN z = cgetg(l,t_VECSMALL);
     487    82947318 :   __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     488    82948972 :   return z;
     489             : }
     490             : 
     491             : static void
     492    94842230 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     493             : {
     494             :   long i;
     495   629532243 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
     496    95101485 : }
     497             : static GEN
     498    90740920 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     499             : {
     500    90740920 :   GEN z = cgetg(l,t_VECSMALL);
     501    90560767 :   __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     502    90789939 :   return z;
     503             : }
     504             : 
     505             : GEN
     506     1634241 : FpM_mul(GEN x, GEN y, GEN p)
     507             : {
     508     1634241 :   long lx=lg(x), ly=lg(y);
     509             :   GEN z;
     510     1634241 :   pari_sp av = avma;
     511     1634241 :   if (ly==1) return cgetg(1,t_MAT);
     512     1634241 :   if (lx==1) return zeromat(0, ly-1);
     513     1634241 :   if (lgefint(p) == 3)
     514             :   {
     515     1588023 :     ulong pp = uel(p,2);
     516     1588023 :     if (pp == 2)
     517             :     {
     518      711588 :       x = ZM_to_F2m(x);
     519      711634 :       y = ZM_to_F2m(y);
     520      711626 :       z = F2m_to_ZM(F2m_mul(x,y));
     521             :     }
     522             :     else
     523             :     {
     524      876435 :       x = ZM_to_Flm(x, pp);
     525      876444 :       y = ZM_to_Flm(y, pp);
     526      876443 :       z = Flm_to_ZM(Flm_mul(x,y, pp));
     527             :     }
     528             :   } else
     529       46218 :     z = FpM_red(ZM_mul(x, y), p);
     530     1634290 :   return gc_upto(av, z);
     531             : }
     532             : 
     533             : static GEN
     534    28032561 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     535             : {
     536             :   long j;
     537    28032561 :   GEN z = cgetg(ly, t_MAT);
     538    28031287 :   if (SMALL_ULONG(p))
     539   102468350 :     for (j=1; j<ly; j++)
     540    82438758 :       gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
     541             :   else
     542    98534377 :     for (j=1; j<ly; j++)
     543    90529579 :       gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
     544    28034390 :   return z;
     545             : }
     546             : 
     547             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     548             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     549             : static void
     550       66303 : add_slices_ip(long m, long n,
     551             :            GEN A, long ma, long da, long na, long ea,
     552             :            GEN B, long mb, long db, long nb, long eb,
     553             :            GEN M, long dy, long dx, ulong p)
     554             : {
     555       66303 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     556             :   GEN C;
     557             : 
     558     2782382 :   for (j = 1; j <= min_e; j++) {
     559     2716079 :     C = gel(M, j + dx) + dy;
     560   123905915 :     for (i = 1; i <= min_d; i++)
     561   121189836 :       uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
     562   121189975 :                         ucoeff(B, mb + i, nb + j), p);
     563     2780500 :     for (; i <= da; i++)
     564       64560 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     565     2715940 :     for (; i <= db; i++)
     566           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     567     2715940 :     for (; i <= m; i++)
     568           0 :       uel(C, i) = 0;
     569             :   }
     570       69916 :   for (; j <= ea; j++) {
     571        3613 :     C = gel(M, j + dx) + dy;
     572      195159 :     for (i = 1; i <= da; i++)
     573      191546 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     574        3613 :     for (; i <= m; i++)
     575           0 :       uel(C, i) = 0;
     576             :   }
     577       66303 :   for (; j <= eb; j++) {
     578           0 :     C = gel(M, j + dx) + dy;
     579           0 :     for (i = 1; i <= db; i++)
     580           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     581           0 :     for (; i <= m; i++)
     582           0 :       uel(C, i) = 0;
     583             :   }
     584       66303 :   for (; j <= n; j++) {
     585           0 :     C = gel(M, j + dx) + dy;
     586           0 :     for (i = 1; i <= m; i++)
     587           0 :       uel(C, i) = 0;
     588             :   }
     589       66303 : }
     590             : 
     591             : static GEN
     592       33152 : add_slices(long m, long n,
     593             :            GEN A, long ma, long da, long na, long ea,
     594             :            GEN B, long mb, long db, long nb, long eb, ulong p)
     595             : {
     596             :   GEN M;
     597             :   long j;
     598       33152 :   M = cgetg(n + 1, t_MAT);
     599     1415674 :   for (j = 1; j <= n; j++)
     600     1382522 :     gel(M, j) = cgetg(m + 1, t_VECSMALL);
     601       33152 :   add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
     602       33152 :   return M;
     603             : }
     604             : 
     605             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     606             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     607             : static GEN
     608       58016 : subtract_slices(long m, long n,
     609             :                 GEN A, long ma, long da, long na, long ea,
     610             :                 GEN B, long mb, long db, long nb, long eb, ulong p)
     611             : {
     612       58016 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     613       58016 :   GEN M = cgetg(n + 1, t_MAT), C;
     614             : 
     615     2527514 :   for (j = 1; j <= min_e; j++) {
     616     2469498 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     617   116080294 :     for (i = 1; i <= min_d; i++)
     618   113610825 :       uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
     619   113610553 :                         ucoeff(B, mb + i, nb + j), p);
     620     2570189 :     for (; i <= da; i++)
     621      100448 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     622     2544784 :     for (; i <= db; i++)
     623       75043 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     624     2469741 :     for (; i <= m; i++)
     625           0 :       uel(C, i) = 0;
     626             :   }
     627       58016 :   for (; j <= ea; j++) {
     628           0 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     629           0 :     for (i = 1; i <= da; i++)
     630           0 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     631           0 :     for (; i <= m; i++)
     632           0 :       uel(C, i) = 0;
     633             :   }
     634       59435 :   for (; j <= eb; j++) {
     635        1419 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     636       82481 :     for (i = 1; i <= db; i++)
     637       81062 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     638        1419 :     for (; i <= m; i++)
     639           0 :       uel(C, i) = 0;
     640             :   }
     641       59435 :   for (; j <= n; j++)
     642        1419 :     gel(M, j) = zero_Flv(m);
     643       58016 :   return M;
     644             : }
     645             : 
     646             : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
     647             : 
     648             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     649             : static GEN
     650        8288 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
     651             : {
     652             :   pari_sp av;
     653        8288 :   long m1 = (m + 1)/2, m2 = m/2,
     654        8288 :     n1 = (n + 1)/2, n2 = n/2,
     655        8288 :     p1 = (p + 1)/2, p2 = p/2;
     656             :   GEN A11, A12, A22, B11, B21, B22,
     657             :     S1, S2, S3, S4, T1, T2, T3, T4,
     658             :     M1, M2, M3, M4, M5, M6, M7,
     659             :     V1, V2, V3, C;
     660             :   long j;
     661        8288 :   C = cgetg(p + 1, t_MAT);
     662      677039 :   for (j = 1; j <= p; j++)
     663      668751 :     gel(C, j) = cgetg(m + 1, t_VECSMALL);
     664        8288 :   av = avma;
     665        8288 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
     666        8288 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
     667        8288 :   M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
     668        8288 :   if (gc_needed(av, 1))
     669           0 :     (void)gc_all(av, 2, &T2, &M2);  /* destroy S1 */
     670        8288 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
     671        8288 :   if (gc_needed(av, 1))
     672           0 :     (void)gc_all(av, 2, &M2, &T3);  /* destroy T2 */
     673        8288 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
     674        8288 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
     675        8288 :   M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
     676        8288 :   if (gc_needed(av, 1))
     677           0 :     (void)gc_all(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     678        8288 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
     679        8288 :   if (gc_needed(av, 1))
     680           0 :     (void)gc_all(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     681        8288 :   A11 = matslice(A, 1, m1, 1, n1);
     682        8288 :   B11 = matslice(B, 1, n1, 1, p1);
     683        8288 :   M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
     684        8288 :   if (gc_needed(av, 1))
     685           0 :     (void)gc_all(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     686        8288 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     687        8288 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     688        8288 :   M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
     689        8288 :   if (gc_needed(av, 1))
     690           0 :     (void)gc_all(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     691        8288 :   add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
     692        8288 :   if (gc_needed(av, 1))
     693           0 :     (void)gc_all(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy M4 */
     694        8288 :   M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
     695        8288 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
     696        8288 :   if (gc_needed(av, 1))
     697           0 :     (void)gc_all(av, 6, &M2, &T3, &M3, &M1, &M5, &S4);  /* destroy S3 */
     698        8288 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
     699        8288 :   if (gc_needed(av, 1))
     700           0 :     (void)gc_all(av, 6, &M2, &M3, &M1, &M5, &S4, &T4);  /* destroy T3 */
     701        8288 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
     702        8288 :   if (gc_needed(av, 1))
     703           0 :     (void)gc_all(av, 5, &M2, &M3, &S4, &T4, &V1);  /* destroy M1, M5 */
     704        8288 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     705        8288 :   M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
     706        8288 :   if (gc_needed(av, 1))
     707           0 :     (void)gc_all(av, 5, &M2, &M3, &T4, &V1, &M6);  /* destroy S4, B22 */
     708        8288 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     709        8288 :   M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
     710        8288 :   if (gc_needed(av, 1))
     711           0 :     (void)gc_all(av, 5, &M2, &M3, &V1, &M6, &M7);  /* destroy A22, T4 */
     712        8288 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
     713        8288 :   add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
     714        8288 :   if (gc_needed(av, 1))
     715           0 :     (void)gc_all(av, 4, &M2, &M3, &V1, &M7);  /* destroy V3, M6 */
     716        8288 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
     717        8288 :   if (gc_needed(av, 1))
     718           0 :     (void)gc_all(av, 3, &M3, &M7, &V2);  /* destroy V1, M2 */
     719        8288 :   add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
     720        8288 :   add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
     721        8288 :   set_avma(av); return C;
     722             : }
     723             : 
     724             : /* Strassen-Winograd used for dim >= ZM_sw_bound */
     725             : static GEN
     726    28039675 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     727             : {
     728    28039675 :   ulong e = expu(p);
     729             : #ifdef LONG_IS_64BIT /* Beware to update ZM_mul_i if this changes */
     730    23548147 :   long Flm_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
     731             : #else
     732     4491927 :   long Flm_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
     733             : #endif
     734    28040074 :   if (l <= Flm_sw_bound || lx <= Flm_sw_bound || ly <= Flm_sw_bound)
     735    28031786 :     return Flm_mul_classical(x, y, l, lx, ly, p, pi);
     736             :   else
     737        8288 :     return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
     738             : }
     739             : 
     740             : GEN
     741    13794912 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     742             : {
     743    13794912 :   long lx=lg(x), ly=lg(y);
     744    13794912 :   if (ly==1) return cgetg(1,t_MAT);
     745    13794912 :   if (lx==1) return zero_Flm(0, ly-1);
     746    13794912 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
     747             : }
     748             : 
     749             : GEN
     750    14092419 : Flm_mul(GEN x, GEN y, ulong p)
     751             : {
     752    14092419 :   long lx=lg(x), ly=lg(y);
     753    14092419 :   if (ly==1) return cgetg(1,t_MAT);
     754    14092132 :   if (lx==1) return zero_Flm(0, ly-1);
     755    14092132 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
     756             : }
     757             : 
     758             : GEN
     759       92727 : Flm_sqr(GEN x, ulong p)
     760             : {
     761       92727 :   long lx = lg(x);
     762       92727 :   if (lx==1) return cgetg(1,t_MAT);
     763       92727 :   return Flm_mul_i(x, x, lx, lx, lx, p, get_Fl_red(p));
     764             : }
     765             : 
     766             : struct _Flm
     767             : {
     768             :   ulong p;
     769             :   long n;
     770             : };
     771             : 
     772             : static GEN
     773       21403 : _Flm_mul(void *E , GEN x, GEN y)
     774       21403 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
     775             : static GEN
     776       92727 : _Flm_sqr(void *E, GEN x)
     777       92727 : { return Flm_sqr(x,((struct _Flm*)E)->p); }
     778             : static GEN
     779         847 : _Flm_one(void *E)
     780         847 : { return matid_Flm(((struct _Flm*)E)->n); }
     781             : GEN
     782       51360 : Flm_powu(GEN x, ulong n, ulong p)
     783             : {
     784             :   struct _Flm d;
     785       51360 :   if (!n) return matid(lg(x)-1);
     786       51360 :   d.p = p;
     787       51360 :   return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
     788             : }
     789             : GEN
     790         847 : Flm_powers(GEN x, ulong n, ulong p)
     791             : {
     792             :   struct _Flm d;
     793         847 :   d.p = p;
     794         847 :   d.n = lg(x)-1;
     795         847 :   return gen_powers(x, n, 1, (void*)&d,
     796             :                     &_Flm_sqr, &_Flm_mul, &_Flm_one);
     797             : }
     798             : 
     799             : static GEN
     800           0 : _FpM_mul(void *p , GEN x, GEN y)
     801           0 : { return FpM_mul(x,y,(GEN)p); }
     802             : static GEN
     803           0 : _FpM_sqr(void *p, GEN x)
     804           0 : { return FpM_mul(x,x,(GEN)p); }
     805             : GEN
     806           0 : FpM_powu(GEN x, ulong n, GEN p)
     807             : {
     808           0 :   if (!n) return matid(lg(x)-1);
     809           0 :   if (lgefint(p) == 3)
     810             :   {
     811           0 :     pari_sp av = avma;
     812           0 :     ulong pp = uel(p,2);
     813             :     GEN z;
     814           0 :     if (pp == 2)
     815           0 :       z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
     816             :     else
     817           0 :       z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
     818           0 :     return gc_upto(av, z);
     819             :   }
     820           0 :   return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
     821             : }
     822             : 
     823             : /*Multiple a column vector by a line vector to make a matrix*/
     824             : GEN
     825          14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
     826             : {
     827          14 :   long i,j, lx=lg(x), ly=lg(y);
     828             :   GEN z;
     829          14 :   if (ly==1) return cgetg(1,t_MAT);
     830          14 :   z = cgetg(ly, t_MAT);
     831          49 :   for (j=1; j < ly; j++)
     832             :   {
     833          35 :     GEN zj = cgetg(lx,t_VECSMALL);
     834         112 :     for (i=1; i<lx; i++)
     835          77 :       uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
     836          35 :     gel(z,j) = zj;
     837             :   }
     838          14 :   return z;
     839             : }
     840             : 
     841             : /*Multiple a column vector by a line vector to make a matrix*/
     842             : GEN
     843        5651 : FpC_FpV_mul(GEN x, GEN y, GEN p)
     844             : {
     845        5651 :   long i,j, lx=lg(x), ly=lg(y);
     846             :   GEN z;
     847        5651 :   if (ly==1) return cgetg(1,t_MAT);
     848        5651 :   z = cgetg(ly,t_MAT);
     849       63794 :   for (j=1; j < ly; j++)
     850             :   {
     851       58143 :     GEN zj = cgetg(lx,t_COL);
     852     1379166 :     for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
     853       58143 :     gel(z, j) = zj;
     854             :   }
     855        5651 :   return z;
     856             : }
     857             : 
     858             : /* Multiply a line vector by a column and return a scalar (t_INT) */
     859             : GEN
     860     8500533 : FpV_dotproduct(GEN x, GEN y, GEN p)
     861             : {
     862     8500533 :   long i, lx = lg(x);
     863             :   pari_sp av;
     864             :   GEN c;
     865     8500533 :   if (lx == 1) return gen_0;
     866     8498027 :   av = avma; c = mulii(gel(x,1),gel(y,1));
     867    36634457 :   for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
     868     8498024 :   return gc_INT(av, modii(c,p));
     869             : }
     870             : GEN
     871           0 : FpV_dotsquare(GEN x, GEN p)
     872             : {
     873           0 :   long i, lx = lg(x);
     874             :   pari_sp av;
     875             :   GEN c;
     876           0 :   if (lx == 1) return gen_0;
     877           0 :   av = avma; c = sqri(gel(x,1));
     878           0 :   for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
     879           0 :   return gc_INT(av, modii(c,p));
     880             : }
     881             : 
     882             : INLINE ulong
     883     9399763 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
     884             : {
     885     9399763 :   ulong c = uel(x,0) * uel(y,0);
     886             :   long k;
     887    72139997 :   for (k = 1; k < lx; k++) {
     888    62740234 :     c += uel(x,k) * uel(y,k);
     889    62740234 :     if (c & HIGHBIT) c %= p;
     890             :   }
     891     9399763 :   return c % p;
     892             : }
     893             : 
     894             : INLINE ulong
     895      741410 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
     896             : {
     897             :   ulong l0, l1, v1, h0, h1;
     898      741410 :   long i = 0;
     899             :   LOCAL_OVERFLOW;
     900             :   LOCAL_HIREMAINDER;
     901      741410 :   l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
     902    15418120 :   while (++i < lx) {
     903    14676710 :     l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
     904    14676710 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     905             :   }
     906      741410 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     907      464137 :   else return remlll_pre(v1, h1, l1, p, pi);
     908             : }
     909             : 
     910             : ulong
     911      730736 : Flv_dotproduct(GEN x, GEN y, ulong p)
     912             : {
     913      730736 :   long lx = lg(x)-1;
     914      730736 :   if (lx == 0) return 0;
     915      730736 :   if (SMALL_ULONG(p))
     916      730736 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     917             :   else
     918           0 :     return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
     919             : }
     920             : 
     921             : ulong
     922       61822 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
     923             : {
     924       61822 :   long lx = lg(x)-1;
     925       61822 :   if (lx == 0) return 0;
     926       61822 :   if (SMALL_ULONG(p))
     927       56060 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     928             :   else
     929        5762 :     return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
     930             : }
     931             : 
     932             : ulong
     933     9416234 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
     934             : {
     935     9416234 :   long lx = minss(lgpol(x), lgpol(y));
     936     9416292 :   if (lx == 0) return 0;
     937     9348627 :   if (pi)
     938      735648 :     return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
     939             :   else
     940     8612979 :     return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
     941             : }
     942             : ulong
     943           0 : Flx_dotproduct(GEN x, GEN y, ulong p)
     944           0 : { return Flx_dotproduct_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
     945             : 
     946             : GEN
     947     1754863 : FpM_FpC_mul(GEN x, GEN y, GEN p)
     948             : {
     949     1754863 :   long lx = lg(x);
     950     1754863 :   return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
     951             : }
     952             : GEN
     953     1077185 : Flm_Flc_mul(GEN x, GEN y, ulong p)
     954             : {
     955     1077185 :   long l, lx = lg(x);
     956     1077185 :   if (lx==1) return cgetg(1,t_VECSMALL);
     957     1077185 :   l = lgcols(x);
     958     1077185 :   if (p==2)
     959      329854 :     return Flm_Flc_mul_i_2(x, y, lx, l);
     960      747331 :   else if (SMALL_ULONG(p))
     961      509993 :     return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
     962             :   else
     963      237338 :     return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
     964             : }
     965             : 
     966             : /* allow pi = 0 */
     967             : GEN
     968        6790 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     969             : {
     970        6790 :   long l, lx = lg(x);
     971             :   GEN z;
     972        6790 :   if (lx==1) return cgetg(1,t_VECSMALL);
     973        6790 :   l = lgcols(x);
     974        6790 :   z = cgetg(l, t_VECSMALL);
     975        6790 :   if (SMALL_ULONG(p))
     976        4406 :     __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     977             :   else
     978        2384 :     __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     979        6791 :   return z;
     980             : }
     981             : 
     982             : /* allow pi = 0 */
     983             : GEN
     984     7616400 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
     985             : {
     986     7616400 :   long l, lx = lg(x);
     987             :   GEN z;
     988     7616400 :   if (lx==1) return pol0_Flx(sv);
     989     7616400 :   l = lgcols(x);
     990     7614585 :   z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
     991     7614025 :   if (SMALL_ULONG(p))
     992     3348759 :     __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
     993             :   else
     994     4265266 :     __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
     995     7624744 :   return Flx_renormalize(z, l + 1);
     996             : }
     997             : 
     998             : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
     999             : GEN
    1000     1234396 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
    1001             : {
    1002     1234396 :   long i, l, lx = lg(x);
    1003             :   GEN z;
    1004     1234396 :   if (lx==1) return pol_0(v);
    1005     1234396 :   l = lgcols(x);
    1006     1234396 :   z = new_chunk(l+1);
    1007     1688675 :   for (i=l-1; i; i--)
    1008             :   {
    1009     1631689 :     pari_sp av = avma;
    1010     1631689 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
    1011     1631686 :     p1 = modii(p1, p);
    1012     1631688 :     if (signe(p1))
    1013             :     {
    1014     1177409 :       if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
    1015     1177409 :       gel(z,i+1) = gc_INT(av, p1);
    1016     1177410 :       break;
    1017             :     }
    1018      454279 :     set_avma(av);
    1019             :   }
    1020     1234396 :   if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
    1021     1177410 :   z[0] = evaltyp(t_POL) | _evallg(i+2);
    1022     1177410 :   z[1] = evalsigne(1) | evalvarn(v);
    1023     3437103 :   for (; i; i--)
    1024             :   {
    1025     2259694 :     pari_sp av = avma;
    1026     2259694 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
    1027     2259690 :     gel(z,i+1) = gc_INT(av, modii(p1,p));
    1028             :   }
    1029     1177409 :   return z;
    1030             : }
    1031             : 
    1032             : /********************************************************************/
    1033             : /**                                                                **/
    1034             : /**                           TRANSPOSITION                        **/
    1035             : /**                                                                **/
    1036             : /********************************************************************/
    1037             : 
    1038             : /* == zm_transpose */
    1039             : GEN
    1040      735413 : Flm_transpose(GEN x)
    1041             : {
    1042      735413 :   long i, dx, lx = lg(x);
    1043             :   GEN y;
    1044      735413 :   if (lx == 1) return cgetg(1,t_MAT);
    1045      735280 :   dx = lgcols(x); y = cgetg(dx,t_MAT);
    1046     9445693 :   for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
    1047      735278 :   return y;
    1048             : }
    1049             : 
    1050             : /********************************************************************/
    1051             : /**                                                                **/
    1052             : /**                           SCALAR MATRICES                      **/
    1053             : /**                                                                **/
    1054             : /********************************************************************/
    1055             : 
    1056             : GEN
    1057        2464 : gen_matid(long n, void *E, const struct bb_field *S)
    1058             : {
    1059        2464 :   GEN y = cgetg(n+1,t_MAT), _0, _1;
    1060             :   long i;
    1061        2464 :   if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
    1062        2464 :   _0 = S->s(E,0);
    1063        2464 :   _1 = S->s(E,1);
    1064       11515 :   for (i=1; i<=n; i++)
    1065             :   {
    1066        9051 :     GEN z = const_col(n, _0); gel(z,i) = _1;
    1067        9051 :     gel(y, i) = z;
    1068             :   }
    1069        2464 :   return y;
    1070             : }
    1071             : 
    1072             : GEN
    1073          35 : matid_F2xqM(long n, GEN T)
    1074             : {
    1075             :   void *E;
    1076          35 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1077          35 :   return gen_matid(n, E, S);
    1078             : }
    1079             : GEN
    1080          56 : matid_FlxqM(long n, GEN T, ulong p)
    1081             : {
    1082             :   void *E;
    1083          56 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1084          56 :   return gen_matid(n, E, S);
    1085             : }
    1086             : 
    1087             : GEN
    1088     1435404 : matid_Flm(long n)
    1089             : {
    1090     1435404 :   GEN y = cgetg(n+1,t_MAT);
    1091             :   long i;
    1092     1435397 :   if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
    1093     9107761 :   for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
    1094     1435371 :   return y;
    1095             : }
    1096             : 
    1097             : GEN
    1098          42 : scalar_Flm(long s, long n)
    1099             : {
    1100             :   long i;
    1101          42 :   GEN y = cgetg(n+1,t_MAT);
    1102         560 :   for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
    1103          42 :   return y;
    1104             : }
    1105             : 
    1106             : /********************************************************************/
    1107             : /**                                                                **/
    1108             : /**                           CONVERSIONS                          **/
    1109             : /**                                                                **/
    1110             : /********************************************************************/
    1111             : GEN
    1112    56246745 : ZV_to_Flv(GEN x, ulong p)
    1113   709055582 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
    1114             : 
    1115             : GEN
    1116    12901364 : ZM_to_Flm(GEN x, ulong p)
    1117    68043873 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
    1118             : 
    1119             : GEN
    1120       14722 : ZVV_to_FlvV(GEN x, ulong p)
    1121      112431 : { pari_APPLY_type(t_VEC, ZV_to_Flv(gel(x,i), p)) }
    1122             : 
    1123             : GEN
    1124        4940 : ZMV_to_FlmV(GEN x, ulong m)
    1125       41848 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
    1126             : 
    1127             : /*                          TO INTMOD                        */
    1128             : static GEN
    1129    18978188 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
    1130             : static GEN
    1131      578850 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
    1132             : 
    1133             : GEN
    1134        3997 : Fp_to_mod(GEN z, GEN p)
    1135             : {
    1136        3997 :   retmkintmod(modii(z, p), icopy(p));
    1137             : }
    1138             : 
    1139             : GEN
    1140      997241 : FpX_to_mod_raw(GEN z, GEN p)
    1141             : {
    1142      997241 :   long i, l = lg(z);
    1143             :   GEN x;
    1144      997241 :   if (l == 2)
    1145             :   {
    1146      522424 :     x = cgetg(3,t_POL); x[1] = z[1];
    1147      522424 :     gel(x,2) = mkintmod(gen_0,p); return x;
    1148             :   }
    1149      474817 :   x = cgetg(l,t_POL); x[1] = z[1];
    1150     3918607 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1151      474817 :   return normalizepol_lg(x,l);
    1152             : }
    1153             : 
    1154             : /* z in Z[X], return z * Mod(1,p), normalized*/
    1155             : GEN
    1156     1260819 : FpX_to_mod(GEN z, GEN p)
    1157             : {
    1158     1260819 :   long i, l = lg(z);
    1159             :   GEN x;
    1160     1260819 :   if (l == 2)
    1161             :   {
    1162        1197 :     x = cgetg(3,t_POL); x[1] = z[1];
    1163        1197 :     gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
    1164             :   }
    1165     1259622 :   x = cgetg(l,t_POL); x[1] = z[1]; p = icopy(p);
    1166    16690366 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1167     1259622 :   return normalizepol_lg(x,l);
    1168             : }
    1169             : 
    1170             : GEN
    1171       84021 : FpXC_to_mod(GEN z, GEN p)
    1172             : {
    1173       84021 :   long i,l = lg(z);
    1174       84021 :   GEN x = cgetg(l, t_COL);
    1175       84021 :   p = icopy(p);
    1176      474159 :   for (i=1; i<l; i++)
    1177      390138 :     gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
    1178       84021 :   return x;
    1179             : }
    1180             : 
    1181             : static GEN
    1182           0 : FpXC_to_mod_raw(GEN x, GEN p)
    1183           0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
    1184             : 
    1185             : GEN
    1186           0 : FpXM_to_mod(GEN z, GEN p)
    1187             : {
    1188           0 :   long i,l = lg(z);
    1189           0 :   GEN x = cgetg(l, t_MAT);
    1190           0 :   p = icopy(p);
    1191           0 :   for (i=1; i<l; i++)
    1192           0 :     gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
    1193           0 :   return x;
    1194             : }
    1195             : 
    1196             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1197             : GEN
    1198       36481 : FpV_to_mod(GEN z, GEN p)
    1199             : {
    1200       36481 :   long i,l = lg(z);
    1201       36481 :   GEN x = cgetg(l, t_VEC);
    1202       36481 :   if (l == 1) return x;
    1203       36481 :   p = icopy(p);
    1204       73529 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1205       36481 :   return x;
    1206             : }
    1207             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1208             : GEN
    1209         105 : FpC_to_mod(GEN z, GEN p)
    1210             : {
    1211         105 :   long i, l = lg(z);
    1212         105 :   GEN x = cgetg(l, t_COL);
    1213         105 :   if (l == 1) return x;
    1214          91 :   p = icopy(p);
    1215         805 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1216          91 :   return x;
    1217             : }
    1218             : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
    1219             : GEN
    1220         226 : FpM_to_mod(GEN z, GEN p)
    1221             : {
    1222         226 :   long i, j, m, l = lg(z);
    1223         226 :   GEN  x = cgetg(l,t_MAT), y, zi;
    1224         226 :   if (l == 1) return x;
    1225         205 :   m = lgcols(z);
    1226         205 :   p = icopy(p);
    1227        2456 :   for (i=1; i<l; i++)
    1228             :   {
    1229        2251 :     gel(x,i) = cgetg(m,t_COL);
    1230        2251 :     y = gel(x,i); zi= gel(z,i);
    1231       66799 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1232             :   }
    1233         205 :   return x;
    1234             : }
    1235             : GEN
    1236          28 : Flc_to_mod(GEN z, ulong pp)
    1237             : {
    1238          28 :   long i, l = lg(z);
    1239          28 :   GEN p, x = cgetg(l, t_COL);
    1240          28 :   if (l == 1) return x;
    1241          28 :   p = utoipos(pp);
    1242         112 :   for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
    1243          28 :   return x;
    1244             : }
    1245             : GEN
    1246       59694 : Flm_to_mod(GEN z, ulong pp)
    1247             : {
    1248       59694 :   long i, j, m, l = lg(z);
    1249       59694 :   GEN p, x = cgetg(l,t_MAT), y, zi;
    1250       59694 :   if (l == 1) return x;
    1251       59666 :   m = lgcols(z);
    1252       59666 :   p = utoipos(pp);
    1253      201328 :   for (i=1; i<l; i++)
    1254             :   {
    1255      141662 :     gel(x,i) = cgetg(m,t_COL);
    1256      141662 :     y = gel(x,i); zi= gel(z,i);
    1257      720428 :     for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
    1258             :   }
    1259       59666 :   return x;
    1260             : }
    1261             : 
    1262             : GEN
    1263         574 : FpVV_to_mod(GEN z, GEN p)
    1264             : {
    1265         574 :   long i, j, m, l = lg(z);
    1266         574 :   GEN  x = cgetg(l,t_VEC), y, zi;
    1267         574 :   if (l == 1) return x;
    1268         574 :   m = lgcols(z);
    1269         574 :   p = icopy(p);
    1270        1246 :   for (i=1; i<l; i++)
    1271             :   {
    1272         672 :     gel(x,i) = cgetg(m,t_VEC);
    1273         672 :     y = gel(x,i); zi= gel(z,i);
    1274        2016 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1275             :   }
    1276         574 :   return x;
    1277             : }
    1278             : 
    1279             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1280             : GEN
    1281           7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
    1282             : {
    1283           7 :   long i,l = lg(z);
    1284           7 :   GEN x = cgetg(l, t_COL);
    1285           7 :   if (l == 1) return x;
    1286           7 :   p = icopy(p);
    1287           7 :   T = FpX_to_mod_raw(T, p);
    1288          21 :   for (i=1; i<l; i++)
    1289          14 :     gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
    1290           7 :   return x;
    1291             : }
    1292             : 
    1293             : static GEN
    1294      582533 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
    1295             : {
    1296      582533 :   GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
    1297      582533 :   return mkpolmod(z, T);
    1298             : }
    1299             : 
    1300             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1301             : GEN
    1302          28 : FqC_to_mod(GEN z, GEN T, GEN p)
    1303             : {
    1304             :   GEN x;
    1305          28 :   long i,l = lg(z);
    1306          28 :   if (!T) return FpC_to_mod(z, p);
    1307          28 :   x = cgetg(l, t_COL);
    1308          28 :   if (l == 1) return x;
    1309          28 :   p = icopy(p);
    1310          28 :   T = FpX_to_mod_raw(T, p);
    1311         504 :   for (i=1; i<l; i++)
    1312         476 :     gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
    1313          28 :   return x;
    1314             : }
    1315             : 
    1316             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1317             : static GEN
    1318      107975 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
    1319      690032 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
    1320             : 
    1321             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1322             : GEN
    1323       26145 : FqM_to_mod(GEN z, GEN T, GEN p)
    1324             : {
    1325             :   GEN x;
    1326       26145 :   long i,l = lg(z);
    1327       26145 :   if (!T) return FpM_to_mod(z, p);
    1328       26145 :   x = cgetg(l, t_MAT);
    1329       26145 :   if (l == 1) return x;
    1330       26117 :   p = icopy(p);
    1331       26117 :   T = FpX_to_mod_raw(T, p);
    1332      134092 :   for (i=1; i<l; i++)
    1333      107975 :     gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
    1334       26117 :   return x;
    1335             : }
    1336             : 
    1337             : /********************************************************************/
    1338             : /*                                                                  */
    1339             : /*                     Blackbox linear algebra                      */
    1340             : /*                                                                  */
    1341             : /********************************************************************/
    1342             : 
    1343             : /* A sparse column (zCs) v is a t_COL with two components C and E which are
    1344             :  * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
    1345             :  * (e_j) is the canonical basis.
    1346             :  * A sparse matrix (zMs) is a t_VEC of zCs */
    1347             : 
    1348             : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
    1349             :  * integer representing an element of Fp. This is important since p can be
    1350             :  * large and p+E[i] would not fit in a C long.  */
    1351             : 
    1352             : /* RgCs and RgMs are similar, except that the type of the component is
    1353             :  * unspecified. Functions handling RgCs/RgMs must be independent of the type
    1354             :  * of E. */
    1355             : 
    1356             : /* Most functions take an argument nbrow which is the number of lines of the
    1357             :  * column/matrix, which cannot be derived from the data. */
    1358             : 
    1359             : GEN
    1360           0 : zCs_to_ZC(GEN R, long nbrow)
    1361             : {
    1362           0 :   GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
    1363           0 :   long j, l = lg(C);
    1364           0 :   for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
    1365           0 :   return c;
    1366             : }
    1367             : 
    1368             : GEN
    1369           0 : zMs_to_ZM(GEN x, long nbrow)
    1370           0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
    1371             : 
    1372             : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
    1373             :  * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
    1374             : GEN
    1375         147 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
    1376             : {
    1377         147 :   pari_sp ltop = avma;
    1378         147 :   long col = 0, n = lg(B)-1, m = 2*n+1;
    1379         147 :   if (ZV_equal0(B)) return zerocol(n);
    1380         147 :   while (++col <= n)
    1381             :   {
    1382         147 :     pari_sp btop = avma, av;
    1383             :     long i, lQ;
    1384         147 :     GEN V, Q, M, W = B;
    1385         147 :     GEN b = cgetg(m+2, t_POL);
    1386         147 :     b[1] = evalsigne(1)|evalvarn(0);
    1387         147 :     gel(b, 2) = gel(W, col);
    1388       46235 :     for (i = 3; i<m+2; i++)
    1389       46088 :       gel(b, i) = cgeti(lgefint(p));
    1390         147 :     av = avma;
    1391       46235 :     for (i = 3; i<m+2; i++)
    1392             :     {
    1393       46088 :       W = f(E, W);
    1394       46088 :       affii(gel(W, col),gel(b, i));
    1395       46088 :       if (gc_needed(av,1))
    1396             :       {
    1397          72 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
    1398          72 :         W = gc_upto(av, W);
    1399             :       }
    1400             :     }
    1401         147 :     b = FpX_renormalize(b, m+2);
    1402         147 :     if (lgpol(b)==0) {set_avma(btop); continue; }
    1403         147 :     M = FpX_halfgcd(b, pol_xn(m, 0), p);
    1404         147 :     Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
    1405         147 :     W = B; lQ =lg(Q);
    1406         147 :     if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
    1407         147 :     V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
    1408         147 :     av = avma;
    1409       22750 :     for (i = lQ-3; i > 1; i--)
    1410             :     {
    1411       22603 :       W = f(E, W);
    1412       22603 :       V = ZC_lincomb(gen_1, gel(Q,i), V, W);
    1413       22603 :       if (gc_needed(av,1))
    1414             :       {
    1415         110 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
    1416         110 :         (void)gc_all(av, 2, &V, &W);
    1417             :       }
    1418             :     }
    1419         147 :     V = FpC_red(V, p);
    1420         147 :     W = FpC_sub(f(E,V), B, p);
    1421         147 :     if (ZV_equal0(W)) return gc_GEN(ltop, V);
    1422           0 :     av = avma;
    1423           0 :     for (i = 1; i <= n; ++i)
    1424             :     {
    1425           0 :       V = W;
    1426           0 :       W = f(E, V);
    1427           0 :       if (ZV_equal0(W))
    1428           0 :         return gc_GEN(ltop, shallowtrans(V));
    1429           0 :       (void)gc_all(av, 2, &V, &W);
    1430             :     }
    1431           0 :     set_avma(btop);
    1432             :   }
    1433           0 :   return NULL;
    1434             : }
    1435             : 
    1436             : GEN
    1437           0 : zMs_ZC_mul(GEN M, GEN B)
    1438             : {
    1439             :   long i, j;
    1440           0 :   long n = lg(B)-1;
    1441           0 :   GEN V = zerocol(n);
    1442           0 :   for (i = 1; i <= n; ++i)
    1443           0 :     if (signe(gel(B, i)))
    1444             :     {
    1445           0 :       GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1446           0 :       long l = lg(C);
    1447           0 :       for (j = 1; j < l; ++j)
    1448             :       {
    1449           0 :         long k = C[j];
    1450           0 :         switch(E[j])
    1451             :         {
    1452           0 :         case 1:
    1453           0 :           gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
    1454           0 :           break;
    1455           0 :         case -1:
    1456           0 :           gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
    1457           0 :           break;
    1458           0 :         default:
    1459           0 :           gel(V, k) = gel(V,k)==gen_0 ? mulis(gel(B, i), E[j]) : addii(gel(V, k), mulis(gel(B, i), E[j]));
    1460           0 :           break;
    1461             :         }
    1462             :       }
    1463             :     }
    1464           0 :   return V;
    1465             : }
    1466             : 
    1467             : GEN
    1468           0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
    1469             : 
    1470             : GEN
    1471       68985 : FpV_FpMs_mul(GEN B, GEN M, GEN p)
    1472             : {
    1473       68985 :   long i, j, lM = lg(M);
    1474       68985 :   GEN V = cgetg(lM,t_VEC);
    1475    28335639 :   for (i = 1; i < lM; ++i)
    1476             :   {
    1477    28266654 :     GEN z, R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1478    28266654 :     pari_sp av = avma;
    1479    28266654 :     long lC = lg(C);
    1480    28266654 :     if (lC == 1) { gel(V,i) = gen_0; continue; }
    1481    28266654 :     z = mulis(gel(B, C[1]), E[1]);
    1482   228756107 :     for (j = 2; j < lC; ++j)
    1483             :     {
    1484   200489453 :       GEN b = gel(B, C[j]);
    1485   200489453 :       switch(E[j])
    1486             :       {
    1487   131566519 :         case  1: z = addii(z, b); break;
    1488    57536025 :         case -1: z = subii(z, b); break;
    1489    11386909 :         default: z = addii(z, mulis(b, E[j])); break;
    1490             :       }
    1491             :     }
    1492    28266654 :     gel(V,i) = gc_INT(av, p? Fp_red(z, p): z);
    1493             :   }
    1494       68985 :   return V;
    1495             : }
    1496             : GEN
    1497           0 : ZV_zMs_mul(GEN B, GEN M) { return FpV_FpMs_mul(B, M, NULL); }
    1498             : 
    1499             : GEN
    1500           0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
    1501             : {
    1502           0 :   pari_sp av = avma, av2;
    1503           0 :   GEN xi, xb, pi = gen_1, P;
    1504             :   long i;
    1505           0 :   if (!C) {
    1506           0 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    1507           0 :     if (!C) pari_err_INV("ZlM_gauss", a);
    1508             :   }
    1509           0 :   P = utoipos(p);
    1510           0 :   av2 = avma;
    1511           0 :   xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1512           0 :   xb = Flm_to_ZM(xi);
    1513           0 :   for (i = 2; i <= e; i++)
    1514             :   {
    1515           0 :     pi = muliu(pi, p); /* = p^(i-1) */
    1516           0 :     b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
    1517           0 :     if (gc_needed(av,2))
    1518             :     {
    1519           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
    1520           0 :       (void)gc_all(av2,3, &pi,&b,&xb);
    1521             :     }
    1522           0 :     xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1523           0 :     xb = ZM_add(xb, nm_Z_mul(xi, pi));
    1524             :   }
    1525           0 :   return gc_upto(av, xb);
    1526             : }
    1527             : 
    1528             : struct wrapper_modp_s {
    1529             :   GEN (*f)(void*, GEN);
    1530             :   void *E;
    1531             :   GEN p;
    1532             : };
    1533             : 
    1534             : /* compute f(x) mod p */
    1535             : static GEN
    1536           0 : wrap_relcomb_modp(void *E, GEN x)
    1537             : {
    1538           0 :   struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
    1539           0 :   return FpC_red(W->f(W->E, x), W->p);
    1540             : }
    1541             : static GEN
    1542           0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
    1543             : 
    1544             : static GEN
    1545       68838 : wrap_relker(void*E, GEN x)
    1546             : {
    1547       68838 :   struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
    1548       68838 :   return FpV_FpMs_mul(x, (GEN) W->E, W->p);
    1549             : }
    1550             : 
    1551             : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
    1552             : GEN
    1553           0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
    1554             : {
    1555             :   struct wrapper_modp_s W;
    1556           0 :   pari_sp av = avma;
    1557           0 :   GEN xb, xi, pi = gen_1;
    1558             :   long i;
    1559           0 :   W.E = E;
    1560           0 :   W.f = f;
    1561           0 :   W.p = p;
    1562           0 :   xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
    1563           0 :   if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
    1564           0 :   xb = xi;
    1565           0 :   for (i = 2; i <= e; i++)
    1566             :   {
    1567           0 :     pi = mulii(pi, p); /* = p^(i-1) */
    1568           0 :     B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
    1569           0 :     if (gc_needed(av,2))
    1570             :     {
    1571           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
    1572           0 :       (void)gc_all(av,3, &pi,&B,&xb);
    1573             :     }
    1574           0 :     xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
    1575           0 :     if (!xi) return NULL;
    1576           0 :     if (typ(xi) == t_VEC) return gc_upto(av, xi);
    1577           0 :     xb = ZC_add(xb, ZC_Z_mul(xi, pi));
    1578             :   }
    1579           0 :   return gc_upto(av, xb);
    1580             : }
    1581             : 
    1582             : static GEN
    1583       23044 : vecprow(GEN A, GEN prow)
    1584             : {
    1585       23044 :   return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
    1586             : }
    1587             : 
    1588             : /* Solve the equation MX = A. Return either a solution as a t_COL,
    1589             :  * or the index of a column which is linearly dependent from the others as a
    1590             :  * t_VECSMALL with a single component. */
    1591             : GEN
    1592           0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
    1593             : {
    1594           0 :   pari_sp av = avma;
    1595             :   GEN pcol, prow;
    1596           0 :   long nbi=lg(M)-1, lR;
    1597             :   long i, n;
    1598             :   GEN Mp, Ap, Rp;
    1599             :   pari_timer ti;
    1600           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1601           0 :   RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
    1602           0 :   if (!pcol) return gc_NULL(av);
    1603           0 :   if (DEBUGLEVEL)
    1604           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
    1605           0 :   n = lg(pcol)-1;
    1606           0 :   Mp = cgetg(n+1, t_MAT);
    1607           0 :   for(i=1; i<=n; i++)
    1608           0 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1609           0 :   Ap = zCs_to_ZC(vecprow(A, prow), n);
    1610           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1611           0 :   Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
    1612           0 :   if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
    1613           0 :   if (!Rp) return gc_NULL(av);
    1614           0 :   lR = lg(Rp)-1;
    1615           0 :   if (typ(Rp) == t_COL)
    1616             :   {
    1617           0 :     GEN R = zerocol(nbi+1);
    1618           0 :     for (i=1; i<=lR; i++)
    1619           0 :        gel(R,pcol[i]) = gel(Rp,i);
    1620           0 :     return gc_GEN(av, R);
    1621             :   }
    1622           0 :   for (i = 1; i <= lR; ++i)
    1623           0 :     if (signe(gel(Rp, i)))
    1624           0 :       return gc_uptoleaf(av, mkvecsmall(pcol[i]));
    1625             :   return NULL; /* LCOV_EXCL_LINE */
    1626             : }
    1627             : 
    1628             : GEN
    1629           0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
    1630             : {
    1631           0 :   return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1632             : }
    1633             : 
    1634             : GEN
    1635           0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
    1636             : {
    1637             :   GEN res;
    1638           0 :   pari_CATCH(e_INV)
    1639             :   {
    1640           0 :     GEN E = pari_err_last();
    1641           0 :     GEN x = err_get_compo(E,2);
    1642           0 :     if (typ(x) != t_INTMOD) pari_err(0,E);
    1643           0 :     if (DEBUGLEVEL)
    1644           0 :       pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
    1645           0 :     res = NULL;
    1646             :   } pari_TRY {
    1647           0 :     res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1648           0 :   } pari_ENDCATCH
    1649           0 :   return res;
    1650             : }
    1651             : 
    1652             : static GEN
    1653         147 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
    1654             : {
    1655         147 :   long i, j, oldf = 0, f = 0;
    1656         147 :   long lrow = lg(prow), lM = lg(M);
    1657         147 :   GEN W = const_vecsmall(lM-1,1);
    1658         147 :   GEN R = cgetg(lrow, t_VEC);
    1659      228354 :   for (i=1; i<lrow; i++)
    1660      228207 :     gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
    1661             :   do
    1662             :   {
    1663         442 :     oldf = f;
    1664      374995 :     for(i=1; i<lM; i++)
    1665      374553 :       if (W[i])
    1666             :       {
    1667      131751 :         GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
    1668      131751 :         long c=0, cj=0, lF = lg(F);
    1669     1093815 :         for(j=1; j<lF; j++)
    1670      962064 :           if (!gel(R,F[j])) { c++; cj=j; }
    1671      131751 :         if (c>=2) continue;
    1672      112310 :         if (c==1)
    1673             :         {
    1674       32891 :           pari_sp av = avma;
    1675       32891 :           GEN s = gen_0;
    1676      272904 :           for(j=1; j<lF; j++)
    1677      240013 :             if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
    1678             :           /* s /= -E[cj] mod p */
    1679       32891 :           s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
    1680       32891 :           gel(R,F[cj]) = gc_INT(av, s);
    1681       32891 :           f++;
    1682             :         }
    1683      112310 :         W[i]=0;
    1684             :       }
    1685         442 :   } while(oldf!=f);
    1686      228354 :   for (i=1; i<lrow; i++)
    1687      228207 :     if (!gel(R,i)) gel(R,i) = gen_0;
    1688         147 :   return R;
    1689             : }
    1690             : 
    1691             : /* Return a linear form Y such that YM is essentially 0 */
    1692             : GEN
    1693         147 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
    1694             : {
    1695         147 :   pari_sp av = avma, av2;
    1696             :   GEN Mp, pcol, prow;
    1697             :   long i, n;
    1698             :   pari_timer ti;
    1699             :   struct wrapper_modp_s W;
    1700         147 :   if (DEBUGLEVEL) timer_start(&ti);
    1701         147 :   RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
    1702         147 :   if (!pcol) return gc_NULL(av);
    1703         147 :   if (DEBUGLEVEL)
    1704           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
    1705         147 :   n = lg(pcol)-1;
    1706         147 :   Mp = cgetg(n+1, t_MAT);
    1707       23191 :   for (i=1; i<=n; i++) gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1708         147 :   W.E = (void*)Mp;
    1709         147 :   W.p = p;
    1710         147 :   av2 = avma;
    1711           0 :   for(;; set_avma(av2))
    1712           0 :   {
    1713         147 :     GEN R, Rp, B = random_FpV(n, p), MB = FpV_FpMs_mul(B, Mp, p);
    1714         147 :     if (DEBUGLEVEL) timer_start(&ti);
    1715         147 :     pari_CATCH(e_INV)
    1716             :     {
    1717           0 :       GEN E = pari_err_last(), x = err_get_compo(E,2);
    1718           0 :       if (typ(x) != t_INTMOD) pari_err(0,E);
    1719           0 :       if (DEBUGLEVEL)
    1720           0 :         pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
    1721           0 :       Rp = NULL;
    1722             :     } pari_TRY {
    1723         147 :       Rp = gen_FpM_Wiedemann((void*)&W, wrap_relker, MB, p);
    1724         147 :     } pari_ENDCATCH
    1725         147 :     if (!Rp) continue;
    1726         147 :     if (typ(Rp)==t_COL)
    1727             :     {
    1728         147 :       Rp = FpC_sub(Rp,B,p);
    1729         147 :       if (ZV_equal0(Rp)) continue;
    1730             :     }
    1731         147 :     R = FpMs_structelim_back(M, Rp, prow, p);
    1732         147 :     if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
    1733         147 :     return gc_GEN(av, R);
    1734             :   }
    1735             : }
    1736             : 
    1737             : GEN
    1738          42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
    1739             : {
    1740          42 :   return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
    1741             : }

Generated by: LCOV version 1.16