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.12.1 lcov report (development 25819-e703fe1174) Lines: 746 951 78.4 %
Date: 2020-09-18 06:10:04 Functions: 112 133 84.2 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /********************************************************************/
      18             : /**                                                                **/
      19             : /**                           REDUCTION                            **/
      20             : /**                                                                **/
      21             : /********************************************************************/
      22             : /* z in Z^n, return lift(Col(z) * Mod(1,p)) */
      23             : GEN
      24    17844834 : FpC_red(GEN x, GEN p)
      25   182696651 : { pari_APPLY_type(t_COL, modii(gel(x,i), p)) }
      26             : 
      27             : /* z in Z^n, return lift(Vec(z) * Mod(1,p)) */
      28             : GEN
      29      207823 : FpV_red(GEN x, GEN p)
      30     2428554 : { pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
      31             : 
      32             : GEN
      33     1491517 : FpC_center(GEN x, GEN p, GEN pov2)
      34    18297011 : { pari_APPLY_type(t_COL, Fp_center(gel(x,i), p, pov2)) }
      35             : 
      36             : /* assume 0 <= u < p and ps2 = p>>1 */
      37             : INLINE void
      38          28 : Fp_center_inplace(GEN u, GEN p, GEN ps2)
      39          28 : { if (abscmpii(u,ps2) > 0) subiiz(u,p,u); }
      40             : 
      41             : void
      42          14 : FpC_center_inplace(GEN z, GEN p, GEN ps2)
      43             : {
      44          14 :   long i,l = lg(z);
      45          42 :   for (i=1; i<l; i++)
      46          28 :     Fp_center_inplace(gel(z,i), p, ps2);
      47          14 : }
      48             : 
      49             : GEN
      50      142275 : Flv_center(GEN x, ulong p, ulong ps2)
      51     2243563 : { pari_APPLY_ulong(Fl_center(uel(x,i),p,ps2)) }
      52             : 
      53             : /* z in Mat m,n(Z), return lift(z * Mod(1,p)) */
      54             : GEN
      55     4446169 : FpM_red(GEN x, GEN p)
      56    18838903 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
      57             : 
      58             : GEN
      59      170362 : FpM_center(GEN x, GEN p, GEN pov2)
      60     1346233 : { pari_APPLY_same(FpC_center(gel(x,i), p,pov2)) }
      61             : 
      62             : void
      63          14 : FpM_center_inplace(GEN z, GEN p, GEN pov2)
      64             : {
      65          14 :   long i, l = lg(z);
      66          28 :   for (i=1; i<l; i++) FpC_center_inplace(gel(z,i), p, pov2);
      67          14 : }
      68             : GEN
      69        9947 : Flm_center(GEN x, ulong p, ulong ps2)
      70      152187 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
      71             : 
      72             : GEN
      73         119 : random_FpV(long d, GEN p)
      74             : {
      75             :   long i;
      76         119 :   GEN y = cgetg(d+1,t_VEC);
      77       25048 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
      78         119 :   return y;
      79             : }
      80             : 
      81             : GEN
      82         413 : random_FpC(long d, GEN p)
      83             : {
      84             :   long i;
      85         413 :   GEN y = cgetg(d+1,t_COL);
      86        3409 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
      87         413 :   return y;
      88             : }
      89             : 
      90             : GEN
      91        6022 : random_Flv(long n, ulong p)
      92             : {
      93        6022 :   GEN y = cgetg(n+1, t_VECSMALL);
      94             :   long i;
      95       49646 :   for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
      96        6022 :   return y;
      97             : }
      98             : 
      99             : /********************************************************************/
     100             : /**                                                                **/
     101             : /**                           ADD, SUB                             **/
     102             : /**                                                                **/
     103             : /********************************************************************/
     104             : GEN
     105      185815 : FpC_add(GEN x, GEN y, GEN p)
     106     2672957 : { pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
     107             : 
     108             : GEN
     109           0 : FpV_add(GEN x, GEN y, GEN p)
     110           0 : { pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
     111             : 
     112             : GEN
     113           0 : FpM_add(GEN x, GEN y, GEN p)
     114           0 : { pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
     115             : 
     116             : GEN
     117      600939 : Flv_add(GEN x, GEN y, ulong p)
     118             : {
     119      600939 :   if (p==2)
     120     1116646 :     pari_APPLY_ulong(x[i]^y[i])
     121             :   else
     122     6878748 :     pari_APPLY_ulong(Fl_add(x[i], y[i], p))
     123             : }
     124             : 
     125             : void
     126      212200 : Flv_add_inplace(GEN x, GEN y, ulong p)
     127             : {
     128      212200 :   long i, l = lg(x);
     129      212200 :   if (p==2)
     130      666916 :     for (i = 1; i < l; i++) x[i] ^= y[i];
     131             :   else
     132      167174 :     for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
     133      212200 : }
     134             : 
     135             : ulong
     136        3907 : Flv_sum(GEN x, ulong p)
     137             : {
     138        3907 :   long i, l = lg(x);
     139        3907 :   ulong s = 0;
     140        3907 :   if (p==2)
     141       17629 :     for (i = 1; i < l; i++) s ^= x[i];
     142             :   else
     143           0 :     for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
     144        3907 :   return s;
     145             : }
     146             : 
     147             : GEN
     148     1834244 : FpC_sub(GEN x, GEN y, GEN p)
     149    14402717 : { pari_APPLY_type(t_COL, Fp_sub(gel(x,i), gel(y,i), p)) }
     150             : 
     151             : GEN
     152           0 : FpV_sub(GEN x, GEN y, GEN p)
     153           0 : { pari_APPLY_type(t_VEC, Fp_sub(gel(x,i), gel(y,i), p)) }
     154             : 
     155             : GEN
     156      436961 : FpM_sub(GEN x, GEN y, GEN p)
     157     1790263 : { pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
     158             : 
     159             : GEN
     160   156244468 : Flv_sub(GEN x, GEN y, ulong p)
     161   680240804 : { pari_APPLY_ulong(Fl_sub(uel(x,i), uel(y,i), p)) }
     162             : 
     163             : void
     164           0 : Flv_sub_inplace(GEN x, GEN y, ulong p)
     165             : {
     166           0 :   long i, l = lg(x);
     167           0 :   for (i = 1; i < l; i++) x[i] = Fl_sub(x[i], y[i], p);
     168           0 : }
     169             : 
     170             : GEN
     171       15428 : Flm_Fl_add(GEN x, ulong y, ulong p)
     172             : {
     173       15428 :   long l = lg(x), i, j;
     174       15428 :   GEN z = cgetg(l,t_MAT);
     175             : 
     176       15428 :   if (l==1) return z;
     177       15428 :   if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
     178       62125 :   for (i=1; i<l; i++)
     179             :   {
     180       46697 :     GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
     181       46697 :     gel(z,i) = zi;
     182      548422 :     for (j=1; j<l; j++) zi[j] = xi[j];
     183       46697 :     zi[i] = Fl_add(zi[i], y, p);
     184             :   }
     185       15428 :   return z;
     186             : }
     187             : GEN
     188         868 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
     189             : 
     190             : GEN
     191        8918 : Flm_add(GEN x, GEN y, ulong p)
     192      150745 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
     193             : 
     194             : GEN
     195    12418022 : Flm_sub(GEN x, GEN y, ulong p)
     196   168651898 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
     197             : 
     198             : /********************************************************************/
     199             : /**                                                                **/
     200             : /**                           MULTIPLICATION                       **/
     201             : /**                                                                **/
     202             : /********************************************************************/
     203             : GEN
     204      763297 : FpC_Fp_mul(GEN x, GEN y, GEN p)
     205    10269040 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
     206             : 
     207             : GEN
     208      571850 : Flv_Fl_mul(GEN x, ulong y, ulong p)
     209     9084902 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
     210             : 
     211             : GEN
     212         854 : Flv_Fl_div(GEN x, ulong y, ulong p)
     213             : {
     214         854 :   return Flv_Fl_mul(x, Fl_inv(y, p), p);
     215             : }
     216             : void
     217           0 : Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
     218             : {
     219           0 :   Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
     220           0 : }
     221             : GEN
     222      125039 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
     223      125039 :   long i, j, h, l = lg(X);
     224      125039 :   GEN A = cgetg(l, t_MAT);
     225      125039 :   if (l == 1) return A;
     226      125039 :   h = lgcols(X);
     227     1165156 :   for (j=1; j<l; j++)
     228             :   {
     229     1040117 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     230    12316551 :     for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
     231     1040117 :     gel(A,j) = a;
     232             :   }
     233      125039 :   return A;
     234             : }
     235             : 
     236             : /* x *= y */
     237             : void
     238     2276472 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
     239             : {
     240             :   long i;
     241     2276472 :   if (HIGHWORD(y | p))
     242     3567580 :     for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
     243             :   else
     244    12658661 :     for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
     245     2276472 : }
     246             : void
     247           0 : Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
     248           0 : { Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
     249             : 
     250             : /* set y *= x */
     251             : void
     252           0 : Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
     253             : {
     254           0 :   long i, j, m, l = lg(y);
     255           0 :   if (l == 1) return;
     256           0 :   m = lgcols(y);
     257           0 :   if (HIGHWORD(x | p))
     258           0 :     for(j=1; j<l; j++)
     259           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
     260             :   else
     261           0 :     for(j=1; j<l; j++)
     262           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
     263             : }
     264             : 
     265             : /* return x * y */
     266             : static GEN
     267     7551565 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
     268             : {
     269     7551565 :   long i, j, m, l = lg(y);
     270     7551565 :   GEN z = cgetg(l, t_MAT);
     271     7550583 :   if (l == 1) return z;
     272     7550583 :   m = lgcols(y);
     273   109477432 :   for(j=1; j<l; j++) {
     274   101880376 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     275   274750780 :     for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
     276             :   }
     277     7597056 :   return z;
     278             : }
     279             : 
     280             : /* return x * y */
     281             : static GEN
     282     3833983 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
     283             : {
     284     3833983 :   long i, j, m, l = lg(y);
     285     3833983 :   GEN z = cgetg(l, t_MAT);
     286     3833983 :   if (l == 1) return z;
     287     3833983 :   m = lgcols(y);
     288    26885041 :   for(j=1; j<l; j++) {
     289    23051058 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     290    67695387 :     for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
     291             :   }
     292     3833983 :   return z;
     293             : }
     294             : 
     295             : /* return x * y */
     296             : GEN
     297    11351821 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
     298             : {
     299    11351821 :   if (HIGHWORD(p))
     300     7551644 :     return Flm_Fl_mul_pre_i(y, x, p, pi);
     301             :   else
     302     3800177 :     return Flm_Fl_mul_OK(y, x, p);
     303             : }
     304             : 
     305             : /* return x * y */
     306             : GEN
     307       33789 : Flm_Fl_mul(GEN y, ulong x, ulong p)
     308             : {
     309       33789 :   if (HIGHWORD(p))
     310           0 :     return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
     311             :   else
     312       33789 :     return Flm_Fl_mul_OK(y, x, p);
     313             : }
     314             : 
     315             : GEN
     316     2055944 : Flv_neg(GEN x, ulong p)
     317    13828796 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
     318             : 
     319             : void
     320        5826 : Flv_neg_inplace(GEN v, ulong p)
     321             : {
     322             :   long i;
     323      185300 :   for (i = 1; i < lg(v); ++i)
     324      179474 :     v[i] = Fl_neg(v[i], p);
     325        5826 : }
     326             : 
     327             : GEN
     328      150705 : Flm_neg(GEN x, ulong p)
     329     2183410 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
     330             : 
     331             : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
     332             : static GEN
     333    12823973 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
     334             : {
     335    12823973 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1));
     336             :   long k;
     337   197122284 :   for (k = 2; k < lx; k++)
     338             :   {
     339   184299652 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     340   184296793 :     if (signe(t)) c = addii(c, t);
     341             :   }
     342    12822632 :   return c;
     343             : }
     344             : 
     345             : static long
     346    68774349 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
     347             : {
     348             :   long k;
     349    68774349 :   long c = coeff(x,i,1) * y[1];
     350  1042862814 :   for (k = 2; k < lx; k++)
     351   974088465 :     c += coeff(x,i,k) * y[k];
     352    68774349 :   return c;
     353             : }
     354             : 
     355             : GEN
     356     4606917 : zm_zc_mul(GEN x, GEN y)
     357             : {
     358     4606917 :   long lx = lg(x), l, i;
     359             :   GEN z;
     360     4606917 :   if (lx == 1) return cgetg(1, t_VECSMALL);
     361     4606917 :   l = lg(gel(x,1));
     362     4606917 :   z = cgetg(l,t_VECSMALL);
     363    73381266 :   for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
     364     4606917 :   return z;
     365             : }
     366             : 
     367             : GEN
     368        1771 : zm_mul(GEN x, GEN y)
     369             : {
     370        1771 :   long i,j,lx=lg(x), ly=lg(y);
     371             :   GEN z;
     372        1771 :   if (ly==1) return cgetg(1,t_MAT);
     373        1771 :   z = cgetg(ly,t_MAT);
     374        1771 :   if (lx==1)
     375             :   {
     376           0 :     for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
     377           0 :     return z;
     378             :   }
     379       24402 :   for (j=1; j<ly; j++)
     380       22631 :     gel(z,j) = zm_zc_mul(x, gel(y,j));
     381        1771 :   return z;
     382             : }
     383             : 
     384             : static ulong
     385   392287608 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
     386             : {
     387   392287608 :   ulong c = ucoeff(x,i,1) * uel(y,1);
     388             :   long k;
     389  5649600764 :   for (k = 2; k < lx; k++) {
     390  5257313156 :     c += ucoeff(x,i,k) * uel(y,k);
     391  5257313156 :     if (c & HIGHBIT) c %= p;
     392             :   }
     393   392287608 :   return c % p;
     394             : }
     395             : 
     396             : static ulong
     397   488005508 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
     398             : {
     399             :   ulong l0, l1, v1, h0, h1;
     400   488005508 :   long k = 1;
     401             :   LOCAL_OVERFLOW;
     402             :   LOCAL_HIREMAINDER;
     403   488005508 :   l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
     404  9038437317 :   while (++k < lx) {
     405  8550431809 :     l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
     406  8550431809 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     407             :   }
     408   488005508 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     409    54142772 :   else return remlll_pre(v1, h1, l1, p, pi);
     410             : }
     411             : 
     412             : static GEN
     413      142077 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
     414             : {
     415             :   long i,j;
     416      142077 :   GEN z = NULL;
     417             : 
     418     1490243 :   for (j=1; j<lx; j++)
     419             :   {
     420     1348166 :     if (!y[j]) continue;
     421      390551 :     if (!z) z = Flv_copy(gel(x,j));
     422    12159359 :     else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
     423             :   }
     424      142077 :   if (!z) z = zero_zv(l-1);
     425      142077 :   return z;
     426             : }
     427             : 
     428             : static GEN
     429     1089691 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
     430             : {
     431     1089691 :   GEN z = cgetg(l,t_COL);
     432             :   long i;
     433    13109086 :   for (i = 1; i < l; i++)
     434             :   {
     435    12019399 :     pari_sp av = avma;
     436    12019399 :     GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
     437    12019435 :     gel(z,i) = gerepileuptoint(av, modii(c,p));
     438             :   }
     439     1089687 :   return z;
     440             : }
     441             : 
     442             : static void
     443    43039387 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
     444             : {
     445             :   long i;
     446   435326614 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
     447    43045027 : }
     448             : static GEN
     449    39935525 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
     450             : {
     451    39935525 :   GEN z = cgetg(l,t_VECSMALL);
     452    39935129 :   __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     453    39935636 :   return z;
     454             : }
     455             : 
     456             : static void
     457    77811021 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     458             : {
     459             :   long i;
     460   568035729 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
     461    77888111 : }
     462             : static GEN
     463    73342440 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     464             : {
     465    73342440 :   GEN z = cgetg(l,t_VECSMALL);
     466    73223075 :   __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     467    73374066 :   return z;
     468             : }
     469             : 
     470             : GEN
     471     1905625 : FpM_mul(GEN x, GEN y, GEN p)
     472             : {
     473     1905625 :   long lx=lg(x), ly=lg(y);
     474             :   GEN z;
     475     1905625 :   pari_sp av = avma;
     476     1905625 :   if (ly==1) return cgetg(1,t_MAT);
     477     1905625 :   if (lx==1) return zeromat(0, ly-1);
     478     1905625 :   if (lgefint(p) == 3)
     479             :   {
     480      989379 :     ulong pp = uel(p,2);
     481      989379 :     if (pp == 2)
     482             :     {
     483      114150 :       x = ZM_to_F2m(x);
     484      114150 :       y = ZM_to_F2m(y);
     485      114150 :       z = F2m_to_ZM(F2m_mul(x,y));
     486             :     }
     487             :     else
     488             :     {
     489      875229 :       x = ZM_to_Flm(x, pp);
     490      875229 :       y = ZM_to_Flm(y, pp);
     491      875229 :       z = Flm_to_ZM(Flm_mul(x,y, pp));
     492             :     }
     493             :   } else
     494      916246 :     z = FpM_red(ZM_mul(x, y), p);
     495     1905625 :   return gerepileupto(av, z);
     496             : }
     497             : 
     498             : static GEN
     499    13912066 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     500             : {
     501             :   long j;
     502    13912066 :   GEN z = cgetg(ly, t_MAT);
     503    13911604 :   if (SMALL_ULONG(p))
     504    49029357 :     for (j=1; j<ly; j++)
     505    39792194 :       gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
     506             :   else
     507    78030695 :     for (j=1; j<ly; j++)
     508    73359627 :       gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
     509    13908231 :   return z;
     510             : }
     511             : 
     512             : /*
     513             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     514             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     515             : */
     516             : static void
     517       67232 : add_slices_ip(long m, long n,
     518             :            GEN A, long ma, long da, long na, long ea,
     519             :            GEN B, long mb, long db, long nb, long eb,
     520             :            GEN M, long dy, long dx, ulong p)
     521             : {
     522       67232 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     523             :   GEN C;
     524             : 
     525     2768771 :   for (j = 1; j <= min_e; j++) {
     526     2701539 :     C = gel(M, j + dx) + dy;
     527   115342324 :     for (i = 1; i <= min_d; i++)
     528   112640785 :       uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
     529   112641008 :                         ucoeff(B, mb + i, nb + j), p);
     530     2727414 :     for (; i <= da; i++)
     531       26098 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     532     2701316 :     for (; i <= db; i++)
     533           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     534     2701316 :     for (; i <= m; i++)
     535           0 :       uel(C, i) = 0;
     536             :   }
     537       68489 :   for (; j <= ea; j++) {
     538        1257 :     C = gel(M, j + dx) + dy;
     539       68091 :     for (i = 1; i <= da; i++)
     540       66834 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     541        1257 :     for (; i <= m; i++)
     542           0 :       uel(C, i) = 0;
     543             :   }
     544       67232 :   for (; j <= eb; j++) {
     545           0 :     C = gel(M, j + dx) + dy;
     546           0 :     for (i = 1; i <= db; i++)
     547           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     548           0 :     for (; i <= m; i++)
     549           0 :       uel(C, i) = 0;
     550             :   }
     551       67232 :   for (; j <= n; j++) {
     552           0 :     C = gel(M, j + dx) + dy;
     553           0 :     for (i = 1; i <= m; i++)
     554           0 :       uel(C, i) = 0;
     555             :   }
     556       67232 : }
     557             : 
     558             : static GEN
     559       33616 : add_slices(long m, long n,
     560             :            GEN A, long ma, long da, long na, long ea,
     561             :            GEN B, long mb, long db, long nb, long eb, ulong p)
     562             : {
     563             :   GEN M;
     564             :   long j;
     565       33616 :   M = cgetg(n + 1, t_MAT);
     566     1393320 :   for (j = 1; j <= n; j++)
     567     1359705 :     gel(M, j) = cgetg(m + 1, t_VECSMALL);
     568       33615 :   add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
     569       33616 :   return M;
     570             : }
     571             : 
     572             : /*
     573             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     574             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     575             : */
     576             : static GEN
     577       58825 : subtract_slices(long m, long n,
     578             :                 GEN A, long ma, long da, long na, long ea,
     579             :                 GEN B, long mb, long db, long nb, long eb, ulong p)
     580             : {
     581       58825 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     582       58825 :   GEN M = cgetg(n + 1, t_MAT), C;
     583             : 
     584     2455604 :   for (j = 1; j <= min_e; j++) {
     585     2396776 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     586   103701023 :     for (i = 1; i <= min_d; i++)
     587   101304238 :       uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
     588   101303978 :                         ucoeff(B, mb + i, nb + j), p);
     589     2434551 :     for (; i <= da; i++)
     590       37506 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     591     2426382 :     for (; i <= db; i++)
     592       29337 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     593     2397045 :     for (; i <= m; i++)
     594           0 :       uel(C, i) = 0;
     595             :   }
     596       58828 :   for (; j <= ea; j++) {
     597           0 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     598           0 :     for (i = 1; i <= da; i++)
     599           0 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     600           0 :     for (; i <= m; i++)
     601           0 :       uel(C, i) = 0;
     602             :   }
     603       59311 :   for (; j <= eb; j++) {
     604         483 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     605       29585 :     for (i = 1; i <= db; i++)
     606       29102 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     607         483 :     for (; i <= m; i++)
     608           0 :       uel(C, i) = 0;
     609             :   }
     610       59311 :   for (; j <= n; j++)
     611         483 :     gel(M, j) = zero_Flv(m);
     612       58828 :   return M;
     613             : }
     614             : 
     615             : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
     616             : 
     617             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     618             : static GEN
     619        8404 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
     620             : {
     621             :   pari_sp av;
     622        8404 :   long m1 = (m + 1)/2, m2 = m/2,
     623        8404 :     n1 = (n + 1)/2, n2 = n/2,
     624        8404 :     p1 = (p + 1)/2, p2 = p/2;
     625             :   GEN A11, A12, A22, B11, B21, B22,
     626             :     S1, S2, S3, S4, T1, T2, T3, T4,
     627             :     M1, M2, M3, M4, M5, M6, M7,
     628             :     V1, V2, V3, C;
     629             :   long j;
     630        8404 :   C = cgetg(p + 1, t_MAT);
     631      680190 :   for (j = 1; j <= p; j++)
     632      671788 :     gel(C, j) = cgetg(m + 1, t_VECSMALL);
     633        8402 :   av = avma;
     634        8402 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
     635        8404 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
     636        8404 :   M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
     637        8404 :   if (gc_needed(av, 1))
     638           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     639        8404 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
     640        8404 :   if (gc_needed(av, 1))
     641           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     642        8404 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
     643        8404 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
     644        8404 :   M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
     645        8404 :   if (gc_needed(av, 1))
     646           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     647        8404 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
     648        8404 :   if (gc_needed(av, 1))
     649           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     650        8404 :   A11 = matslice(A, 1, m1, 1, n1);
     651        8404 :   B11 = matslice(B, 1, n1, 1, p1);
     652        8404 :   M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
     653        8404 :   if (gc_needed(av, 1))
     654           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     655        8404 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     656        8404 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     657        8404 :   M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
     658        8404 :   if (gc_needed(av, 1))
     659           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     660        8404 :   add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
     661        8404 :   if (gc_needed(av, 1))
     662           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy M4 */
     663        8404 :   M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
     664        8404 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
     665        8404 :   if (gc_needed(av, 1))
     666           0 :     gerepileall(av, 6, &M2, &T3, &M3, &M1, &M5, &S4);  /* destroy S3 */
     667        8404 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
     668        8404 :   if (gc_needed(av, 1))
     669           0 :     gerepileall(av, 6, &M2, &M3, &M1, &M5, &S4, &T4);  /* destroy T3 */
     670        8404 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
     671        8404 :   if (gc_needed(av, 1))
     672           1 :     gerepileall(av, 5, &M2, &M3, &S4, &T4, &V1);  /* destroy M1, M5 */
     673        8404 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     674        8404 :   M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
     675        8404 :   if (gc_needed(av, 1))
     676           0 :     gerepileall(av, 5, &M2, &M3, &T4, &V1, &M6);  /* destroy S4, B22 */
     677        8404 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     678        8404 :   M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
     679        8404 :   if (gc_needed(av, 1))
     680           0 :     gerepileall(av, 5, &M2, &M3, &V1, &M6, &M7);  /* destroy A22, T4 */
     681        8404 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
     682        8404 :   add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
     683        8404 :   if (gc_needed(av, 1))
     684           0 :     gerepileall(av, 4, &M2, &M3, &V1, &M7);  /* destroy V3, M6 */
     685        8404 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
     686        8404 :   if (gc_needed(av, 1))
     687           0 :     gerepileall(av, 3, &M3, &M7, &V2);  /* destroy V1, M2 */
     688        8404 :   add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
     689        8404 :   add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
     690        8404 :   set_avma(av); return C;
     691             : }
     692             : 
     693             : /* Strassen-Winograd used for dim >= ZM_sw_bound */
     694             : static GEN
     695    13920213 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     696             : {
     697    13920213 :   ulong e = expu(p);
     698             : #ifdef LONG_IS_64BIT
     699    11632560 :   long ZM_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
     700             : #else
     701     2287751 :   long ZM_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
     702             : #endif
     703    13920311 :   if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
     704    13911907 :     return Flm_mul_classical(x, y, l, lx, ly, p, pi);
     705             :   else
     706        8404 :     return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
     707             : }
     708             : 
     709             : GEN
     710     6831043 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     711             : {
     712     6831043 :   long lx=lg(x), ly=lg(y);
     713     6831043 :   if (ly==1) return cgetg(1,t_MAT);
     714     6831043 :   if (lx==1) return zero_Flm(0, ly-1);
     715     6831043 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
     716             : }
     717             : 
     718             : GEN
     719     7029903 : Flm_mul(GEN x, GEN y, ulong p)
     720             : {
     721     7029903 :   long lx=lg(x), ly=lg(y);
     722     7029903 :   if (ly==1) return cgetg(1,t_MAT);
     723     7029665 :   if (lx==1) return zero_Flm(0, ly-1);
     724     7029665 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
     725             : }
     726             : 
     727             : struct _Flm
     728             : {
     729             :   ulong p;
     730             :   long n;
     731             : };
     732             : 
     733             : static GEN
     734        8505 : _Flm_mul(void *E , GEN x, GEN y)
     735        8505 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
     736             : static GEN
     737       43274 : _Flm_sqr(void *E, GEN x)
     738       43274 : { return Flm_mul(x,x,((struct _Flm*)E)->p); }
     739             : static GEN
     740         315 : _Flm_one(void *E)
     741         315 : { return matid_Flm(((struct _Flm*)E)->n); }
     742             : GEN
     743       24472 : Flm_powu(GEN x, ulong n, ulong p)
     744             : {
     745             :   struct _Flm d;
     746       24472 :   if (!n) return matid(lg(x)-1);
     747       24472 :   d.p = p;
     748       24472 :   return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
     749             : }
     750             : GEN
     751         315 : Flm_powers(GEN x, ulong n, ulong p)
     752             : {
     753             :   struct _Flm d;
     754         315 :   d.p = p;
     755         315 :   d.n = lg(x)-1;
     756         315 :   return gen_powers(x, n, 1, (void*)&d,
     757             :                     &_Flm_sqr, &_Flm_mul, &_Flm_one);
     758             : }
     759             : 
     760             : static GEN
     761           0 : _FpM_mul(void *p , GEN x, GEN y)
     762           0 : { return FpM_mul(x,y,(GEN)p); }
     763             : static GEN
     764           0 : _FpM_sqr(void *p, GEN x)
     765           0 : { return FpM_mul(x,x,(GEN)p); }
     766             : GEN
     767           0 : FpM_powu(GEN x, ulong n, GEN p)
     768             : {
     769           0 :   if (!n) return matid(lg(x)-1);
     770           0 :   if (lgefint(p) == 3)
     771             :   {
     772           0 :     pari_sp av = avma;
     773           0 :     ulong pp = uel(p,2);
     774             :     GEN z;
     775           0 :     if (pp == 2)
     776           0 :       z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
     777             :     else
     778           0 :       z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
     779           0 :     return gerepileupto(av, z);
     780             :   }
     781           0 :   return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
     782             : }
     783             : 
     784             : /*Multiple a column vector by a line vector to make a matrix*/
     785             : GEN
     786          14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
     787             : {
     788          14 :   long i,j, lx=lg(x), ly=lg(y);
     789             :   GEN z;
     790          14 :   if (ly==1) return cgetg(1,t_MAT);
     791          14 :   z = cgetg(ly, t_MAT);
     792          49 :   for (j=1; j < ly; j++)
     793             :   {
     794          35 :     GEN zj = cgetg(lx,t_VECSMALL);
     795         112 :     for (i=1; i<lx; i++)
     796          77 :       uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
     797          35 :     gel(z,j) = zj;
     798             :   }
     799          14 :   return z;
     800             : }
     801             : 
     802             : /*Multiple a column vector by a line vector to make a matrix*/
     803             : GEN
     804        2520 : FpC_FpV_mul(GEN x, GEN y, GEN p)
     805             : {
     806        2520 :   long i,j, lx=lg(x), ly=lg(y);
     807             :   GEN z;
     808        2520 :   if (ly==1) return cgetg(1,t_MAT);
     809        2520 :   z = cgetg(ly,t_MAT);
     810       35896 :   for (j=1; j < ly; j++)
     811             :   {
     812       33376 :     GEN zj = cgetg(lx,t_COL);
     813     1073058 :     for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
     814       33376 :     gel(z, j) = zj;
     815             :   }
     816        2520 :   return z;
     817             : }
     818             : 
     819             : /* Multiply a line vector by a column and return a scalar (t_INT) */
     820             : GEN
     821     2069062 : FpV_dotproduct(GEN x, GEN y, GEN p)
     822             : {
     823     2069062 :   long i, lx = lg(x);
     824             :   pari_sp av;
     825             :   GEN c;
     826     2069062 :   if (lx == 1) return gen_0;
     827     2068117 :   av = avma; c = mulii(gel(x,1),gel(y,1));
     828     9464713 :   for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
     829     2068025 :   return gerepileuptoint(av, modii(c,p));
     830             : }
     831             : GEN
     832           0 : FpV_dotsquare(GEN x, GEN p)
     833             : {
     834           0 :   long i, lx = lg(x);
     835             :   pari_sp av;
     836             :   GEN c;
     837           0 :   if (lx == 1) return gen_0;
     838           0 :   av = avma; c = sqri(gel(x,1));
     839           0 :   for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
     840           0 :   return gerepileuptoint(av, modii(c,p));
     841             : }
     842             : 
     843             : INLINE ulong
     844     4441137 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
     845             : {
     846     4441137 :   ulong c = uel(x,0) * uel(y,0);
     847             :   long k;
     848    47000956 :   for (k = 1; k < lx; k++) {
     849    42559819 :     c += uel(x,k) * uel(y,k);
     850    42559819 :     if (c & HIGHBIT) c %= p;
     851             :   }
     852     4441137 :   return c % p;
     853             : }
     854             : 
     855             : INLINE ulong
     856      798694 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
     857             : {
     858             :   ulong l0, l1, v1, h0, h1;
     859      798694 :   long i = 0;
     860             :   LOCAL_OVERFLOW;
     861             :   LOCAL_HIREMAINDER;
     862      798694 :   l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
     863    16826596 :   while (++i < lx) {
     864    16027902 :     l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
     865    16027902 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     866             :   }
     867      798694 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     868      506711 :   else return remlll_pre(v1, h1, l1, p, pi);
     869             : }
     870             : 
     871             : ulong
     872      321279 : Flv_dotproduct(GEN x, GEN y, ulong p)
     873             : {
     874      321279 :   long lx = lg(x)-1;
     875      321279 :   if (lx == 0) return 0;
     876      321279 :   if (SMALL_ULONG(p))
     877      321279 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     878             :   else
     879           0 :     return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
     880             : }
     881             : 
     882             : ulong
     883       57184 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
     884             : {
     885       57184 :   long lx = lg(x)-1;
     886       57184 :   if (lx == 0) return 0;
     887       57184 :   if (SMALL_ULONG(p))
     888       51692 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     889             :   else
     890        5492 :     return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
     891             : }
     892             : 
     893             : ulong
     894     4889035 : Flx_dotproduct(GEN x, GEN y, ulong p)
     895             : {
     896     4889035 :   long lx = minss(lgpol(x), lgpol(y));
     897     4889077 :   if (lx == 0) return 0;
     898     4861338 :   if (SMALL_ULONG(p))
     899     4068136 :     return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
     900             :   else
     901      793202 :     return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
     902             : }
     903             : 
     904             : GEN
     905     1089691 : FpM_FpC_mul(GEN x, GEN y, GEN p)
     906             : {
     907     1089691 :   long lx = lg(x);
     908     1089691 :   return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
     909             : }
     910             : GEN
     911      288015 : Flm_Flc_mul(GEN x, GEN y, ulong p)
     912             : {
     913      288015 :   long l, lx = lg(x);
     914      288015 :   if (lx==1) return cgetg(1,t_VECSMALL);
     915      288015 :   l = lgcols(x);
     916      288015 :   if (p==2)
     917      142077 :     return Flm_Flc_mul_i_2(x, y, lx, l);
     918      145938 :   else if (SMALL_ULONG(p))
     919      143371 :     return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
     920             :   else
     921        2567 :     return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
     922             : }
     923             : 
     924             : GEN
     925        5138 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     926             : {
     927        5138 :   long l, lx = lg(x);
     928             :   GEN z;
     929        5138 :   if (lx==1) return cgetg(1,t_VECSMALL);
     930        5138 :   l = lgcols(x);
     931        5138 :   z = cgetg(l, t_VECSMALL);
     932        5138 :   if (SMALL_ULONG(p))
     933        4407 :     __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     934             :   else
     935         731 :     __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     936        5138 :   return z;
     937             : }
     938             : 
     939             : GEN
     940     7657183 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
     941             : {
     942     7657183 :   long l, lx = lg(x);
     943             :   GEN z;
     944     7657183 :   if (lx==1) return pol0_Flx(sv);
     945     7657183 :   l = lgcols(x);
     946     7654675 :   z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
     947     7656951 :   if (SMALL_ULONG(p))
     948     3099596 :     __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
     949             :   else
     950     4557355 :     __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
     951     7671702 :   return Flx_renormalize(z, l + 1);
     952             : }
     953             : 
     954             : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
     955             : GEN
     956      304443 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
     957             : {
     958      304443 :   long i, l, lx = lg(x);
     959             :   GEN z;
     960      304443 :   if (lx==1) return pol_0(v);
     961      304443 :   l = lgcols(x);
     962      304443 :   z = new_chunk(l+1);
     963      681119 :   for (i=l-1; i; i--)
     964             :   {
     965      526272 :     pari_sp av = avma;
     966      526272 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
     967      526272 :     p1 = modii(p1, p);
     968      526272 :     if (signe(p1))
     969             :     {
     970      149596 :       if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
     971      149596 :       gel(z,i+1) = gerepileuptoint(av, p1);
     972      149596 :       break;
     973             :     }
     974      376676 :     set_avma(av);
     975             :   }
     976      304443 :   if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
     977      149596 :   z[0] = evaltyp(t_POL) | evallg(i+2);
     978      149596 :   z[1] = evalsigne(1) | evalvarn(v);
     979      427873 :   for (; i; i--)
     980             :   {
     981      278277 :     pari_sp av = avma;
     982      278277 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
     983      278277 :     gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
     984             :   }
     985      149596 :   return z;
     986             : }
     987             : 
     988             : /********************************************************************/
     989             : /**                                                                **/
     990             : /**                           TRANSPOSITION                        **/
     991             : /**                                                                **/
     992             : /********************************************************************/
     993             : 
     994             : /* == zm_transpose */
     995             : GEN
     996      338014 : Flm_transpose(GEN x)
     997             : {
     998      338014 :   long i, dx, lx = lg(x);
     999             :   GEN y;
    1000      338014 :   if (lx == 1) return cgetg(1,t_MAT);
    1001      338014 :   dx = lgcols(x); y = cgetg(dx,t_MAT);
    1002     5353222 :   for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
    1003      338017 :   return y;
    1004             : }
    1005             : 
    1006             : /********************************************************************/
    1007             : /**                                                                **/
    1008             : /**                           SCALAR MATRICES                      **/
    1009             : /**                                                                **/
    1010             : /********************************************************************/
    1011             : 
    1012             : GEN
    1013        3762 : gen_matid(long n, void *E, const struct bb_field *S)
    1014             : {
    1015        3762 :   GEN y = cgetg(n+1,t_MAT), _0, _1;
    1016             :   long i;
    1017        3762 :   if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
    1018        3762 :   _0 = S->s(E,0);
    1019        3762 :   _1 = S->s(E,1);
    1020       16760 :   for (i=1; i<=n; i++)
    1021             :   {
    1022       12998 :     GEN z = const_col(n, _0); gel(z,i) = _1;
    1023       12998 :     gel(y, i) = z;
    1024             :   }
    1025        3762 :   return y;
    1026             : }
    1027             : 
    1028             : GEN
    1029          35 : matid_F2xqM(long n, GEN T)
    1030             : {
    1031             :   void *E;
    1032          35 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1033          35 :   return gen_matid(n, E, S);
    1034             : }
    1035             : GEN
    1036          56 : matid_FlxqM(long n, GEN T, ulong p)
    1037             : {
    1038             :   void *E;
    1039          56 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1040          56 :   return gen_matid(n, E, S);
    1041             : }
    1042             : 
    1043             : GEN
    1044      812164 : matid_Flm(long n)
    1045             : {
    1046      812164 :   GEN y = cgetg(n+1,t_MAT);
    1047             :   long i;
    1048      812163 :   if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
    1049     5383401 :   for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
    1050      812163 :   return y;
    1051             : }
    1052             : 
    1053             : GEN
    1054          42 : scalar_Flm(long s, long n)
    1055             : {
    1056             :   long i;
    1057          42 :   GEN y = cgetg(n+1,t_MAT);
    1058         560 :   for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
    1059          42 :   return y;
    1060             : }
    1061             : 
    1062             : /********************************************************************/
    1063             : /**                                                                **/
    1064             : /**                           CONVERSIONS                          **/
    1065             : /**                                                                **/
    1066             : /********************************************************************/
    1067             : GEN
    1068    23314070 : ZV_to_Flv(GEN x, ulong p)
    1069   419547940 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
    1070             : 
    1071             : GEN
    1072     3447675 : ZM_to_Flm(GEN x, ulong p)
    1073    26522880 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
    1074             : 
    1075             : GEN
    1076        2324 : ZMV_to_FlmV(GEN x, ulong m)
    1077       18438 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
    1078             : 
    1079             : /*                          TO INTMOD                        */
    1080             : static GEN
    1081    12058567 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
    1082             : static GEN
    1083       19376 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
    1084             : 
    1085             : GEN
    1086      305238 : Fp_to_mod(GEN z, GEN p)
    1087             : {
    1088      305238 :   retmkintmod(modii(z, p), icopy(p));
    1089             : }
    1090             : 
    1091             : GEN
    1092      617463 : FpX_to_mod_raw(GEN z, GEN p)
    1093             : {
    1094      617463 :   long i,l = lg(z);
    1095             :   GEN x;
    1096      617463 :   if (l == 2)
    1097             :   {
    1098      524566 :     x = cgetg(3,t_POL); x[1] = z[1];
    1099      524566 :     gel(x,2) = mkintmod(gen_0,p); return x;
    1100             :   }
    1101       92897 :   x = cgetg(l,t_POL);
    1102      434868 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1103       92897 :   x[1] = z[1]; return normalizepol_lg(x,l);
    1104             : }
    1105             : 
    1106             : /* z in Z[X], return z * Mod(1,p), normalized*/
    1107             : GEN
    1108     1242925 : FpX_to_mod(GEN z, GEN p)
    1109             : {
    1110     1242925 :   long i,l = lg(z);
    1111             :   GEN x;
    1112     1242925 :   if (l == 2)
    1113             :   {
    1114        1309 :     x = cgetg(3,t_POL); x[1] = z[1];
    1115        1309 :     gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
    1116             :   }
    1117     1241616 :   x = cgetg(l,t_POL);
    1118     1241616 :   if (l >2) p = icopy(p);
    1119    12239851 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1120     1241616 :   x[1] = z[1]; return normalizepol_lg(x,l);
    1121             : }
    1122             : 
    1123             : GEN
    1124         357 : FpXC_to_mod(GEN z, GEN p)
    1125             : {
    1126         357 :   long i,l = lg(z);
    1127         357 :   GEN x = cgetg(l, t_COL);
    1128         357 :   p = icopy(p);
    1129        8533 :   for (i=1; i<l; i++)
    1130        8176 :     gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
    1131         357 :   return x;
    1132             : }
    1133             : 
    1134             : static GEN
    1135          28 : FpXC_to_mod_raw(GEN x, GEN p)
    1136          84 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
    1137             : 
    1138             : GEN
    1139          14 : FpXM_to_mod(GEN z, GEN p)
    1140             : {
    1141          14 :   long i,l = lg(z);
    1142          14 :   GEN x = cgetg(l, t_MAT);
    1143          14 :   p = icopy(p);
    1144          42 :   for (i=1; i<l; i++)
    1145          28 :     gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
    1146          14 :   return x;
    1147             : }
    1148             : 
    1149             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1150             : GEN
    1151        6629 : FpV_to_mod(GEN z, GEN p)
    1152             : {
    1153        6629 :   long i,l = lg(z);
    1154        6629 :   GEN x = cgetg(l, t_VEC);
    1155        6629 :   if (l == 1) return x;
    1156        6629 :   p = icopy(p);
    1157       13587 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1158        6629 :   return x;
    1159             : }
    1160             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1161             : GEN
    1162         105 : FpC_to_mod(GEN z, GEN p)
    1163             : {
    1164         105 :   long i, l = lg(z);
    1165         105 :   GEN x = cgetg(l, t_COL);
    1166         105 :   if (l == 1) return x;
    1167          91 :   p = icopy(p);
    1168         805 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1169          91 :   return x;
    1170             : }
    1171             : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
    1172             : GEN
    1173       56392 : FpM_to_mod(GEN z, GEN p)
    1174             : {
    1175       56392 :   long i, j, m, l = lg(z);
    1176       56392 :   GEN  x = cgetg(l,t_MAT), y, zi;
    1177       56392 :   if (l == 1) return x;
    1178       56371 :   m = lgcols(z);
    1179       56371 :   p = icopy(p);
    1180      181986 :   for (i=1; i<l; i++)
    1181             :   {
    1182      125615 :     gel(x,i) = cgetg(m,t_COL);
    1183      125615 :     y = gel(x,i); zi= gel(z,i);
    1184      834960 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1185             :   }
    1186       56371 :   return x;
    1187             : }
    1188             : GEN
    1189          28 : Flc_to_mod(GEN z, ulong pp)
    1190             : {
    1191          28 :   long i, l = lg(z);
    1192          28 :   GEN p, x = cgetg(l, t_COL);
    1193          28 :   if (l == 1) return x;
    1194          28 :   p = utoipos(pp);
    1195         112 :   for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
    1196          28 :   return x;
    1197             : }
    1198             : GEN
    1199         189 : Flm_to_mod(GEN z, ulong pp)
    1200             : {
    1201         189 :   long i, j, m, l = lg(z);
    1202         189 :   GEN p, x = cgetg(l,t_MAT), y, zi;
    1203         189 :   if (l == 1) return x;
    1204         161 :   m = lgcols(z);
    1205         161 :   p = utoipos(pp);
    1206        1414 :   for (i=1; i<l; i++)
    1207             :   {
    1208        1253 :     gel(x,i) = cgetg(m,t_COL);
    1209        1253 :     y = gel(x,i); zi= gel(z,i);
    1210       20545 :     for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
    1211             :   }
    1212         161 :   return x;
    1213             : }
    1214             : 
    1215             : GEN
    1216         574 : FpVV_to_mod(GEN z, GEN p)
    1217             : {
    1218         574 :   long i, j, m, l = lg(z);
    1219         574 :   GEN  x = cgetg(l,t_VEC), y, zi;
    1220         574 :   if (l == 1) return x;
    1221         574 :   m = lgcols(z);
    1222         574 :   p = icopy(p);
    1223        1246 :   for (i=1; i<l; i++)
    1224             :   {
    1225         672 :     gel(x,i) = cgetg(m,t_VEC);
    1226         672 :     y = gel(x,i); zi= gel(z,i);
    1227        2016 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1228             :   }
    1229         574 :   return x;
    1230             : }
    1231             : 
    1232             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1233             : GEN
    1234           7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
    1235             : {
    1236           7 :   long i,l = lg(z);
    1237           7 :   GEN x = cgetg(l, t_COL);
    1238           7 :   if (l == 1) return x;
    1239           7 :   p = icopy(p);
    1240           7 :   T = FpX_to_mod_raw(T, p);
    1241          21 :   for (i=1; i<l; i++)
    1242          14 :     gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
    1243           7 :   return x;
    1244             : }
    1245             : 
    1246             : static GEN
    1247      584675 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
    1248             : {
    1249      584675 :   GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
    1250      584675 :   return mkpolmod(z, T);
    1251             : }
    1252             : 
    1253             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1254             : GEN
    1255          14 : FqC_to_mod(GEN z, GEN T, GEN p)
    1256             : {
    1257             :   GEN x;
    1258          14 :   long i,l = lg(z);
    1259          14 :   if (!T) return FpC_to_mod(z, p);
    1260          14 :   x = cgetg(l, t_COL);
    1261          14 :   if (l == 1) return x;
    1262          14 :   p = icopy(p);
    1263          14 :   T = FpX_to_mod_raw(T, p);
    1264         252 :   for (i=1; i<l; i++)
    1265         238 :     gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
    1266          14 :   return x;
    1267             : }
    1268             : 
    1269             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1270             : static GEN
    1271      108115 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
    1272      692552 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
    1273             : 
    1274             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1275             : GEN
    1276       26159 : FqM_to_mod(GEN z, GEN T, GEN p)
    1277             : {
    1278             :   GEN x;
    1279       26159 :   long i,l = lg(z);
    1280       26159 :   if (!T) return FpM_to_mod(z, p);
    1281       26159 :   x = cgetg(l, t_MAT);
    1282       26159 :   if (l == 1) return x;
    1283       26131 :   p = icopy(p);
    1284       26131 :   T = FpX_to_mod_raw(T, p);
    1285      134246 :   for (i=1; i<l; i++)
    1286      108115 :     gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
    1287       26131 :   return x;
    1288             : }
    1289             : 
    1290             : /********************************************************************/
    1291             : /*                                                                  */
    1292             : /*                     Blackbox linear algebra                      */
    1293             : /*                                                                  */
    1294             : /********************************************************************/
    1295             : 
    1296             : /* A sparse column (zCs) v is a t_COL with two components C and E which are
    1297             :  * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
    1298             :  * (e_j) is the canonical basis.
    1299             :  * A sparse matrix (zMs) is a t_VEC of zCs */
    1300             : 
    1301             : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
    1302             :  * integer representing an element of Fp. This is important since p can be
    1303             :  * large and p+E[i] would not fit in a C long.  */
    1304             : 
    1305             : /* RgCs and RgMs are similar, except that the type of the component is
    1306             :  * unspecified. Functions handling RgCs/RgMs must be independent of the type
    1307             :  * of E. */
    1308             : 
    1309             : /* Most functions take an argument nbrow which is the number of lines of the
    1310             :  * column/matrix, which cannot be derived from the data. */
    1311             : 
    1312             : GEN
    1313           0 : zCs_to_ZC(GEN R, long nbrow)
    1314             : {
    1315           0 :   GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
    1316           0 :   long j, l = lg(C);
    1317           0 :   for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
    1318           0 :   return c;
    1319             : }
    1320             : 
    1321             : GEN
    1322           0 : zMs_to_ZM(GEN x, long nbrow)
    1323           0 : { pari_APPLY_same(zCs_to_ZC(gel(x, i), nbrow)) }
    1324             : 
    1325             : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
    1326             :  * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
    1327             : GEN
    1328         119 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
    1329             : {
    1330         119 :   pari_sp ltop = avma;
    1331         119 :   long col = 0, n = lg(B)-1, m = 2*n+1;
    1332         119 :   if (ZV_equal0(B)) return zerocol(n);
    1333         119 :   while (++col <= n)
    1334             :   {
    1335         119 :     pari_sp btop = avma, av;
    1336             :     long i, lQ;
    1337         119 :     GEN V, Q, M, W = B;
    1338         119 :     GEN b = cgetg(m+2, t_POL);
    1339         119 :     b[1] = evalsigne(1)|evalvarn(0);
    1340         119 :     gel(b, 2) = gel(W, col);
    1341       49977 :     for (i = 3; i<m+2; i++)
    1342       49858 :       gel(b, i) = cgeti(lgefint(p));
    1343         119 :     av = avma;
    1344       49977 :     for (i = 3; i<m+2; i++)
    1345             :     {
    1346       49858 :       W = f(E, W);
    1347       49858 :       affii(gel(W, col),gel(b, i));
    1348       49858 :       if (gc_needed(av,1))
    1349             :       {
    1350        1990 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
    1351        1990 :         W = gerepileupto(av, W);
    1352             :       }
    1353             :     }
    1354         119 :     b = FpX_renormalize(b, m+2);
    1355         119 :     if (lgpol(b)==0) {set_avma(btop); continue; }
    1356         119 :     M = FpX_halfgcd(b, pol_xn(m, 0), p);
    1357         119 :     Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
    1358         119 :     W = B; lQ =lg(Q);
    1359         119 :     if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
    1360         119 :     V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
    1361         119 :     av = avma;
    1362       24600 :     for (i = lQ-3; i > 1; i--)
    1363             :     {
    1364       24481 :       W = f(E, W);
    1365       24481 :       V = ZC_lincomb(gen_1, gel(Q,i), V, W);
    1366       24481 :       if (gc_needed(av,1))
    1367             :       {
    1368        1848 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
    1369        1848 :         gerepileall(av, 2, &V, &W);
    1370             :       }
    1371             :     }
    1372         119 :     V = FpC_red(V, p);
    1373         119 :     W = FpC_sub(f(E,V), B, p);
    1374         119 :     if (ZV_equal0(W)) return gerepilecopy(ltop, V);
    1375           0 :     av = avma;
    1376           0 :     for (i = 1; i <= n; ++i)
    1377             :     {
    1378           0 :       V = W;
    1379           0 :       W = f(E, V);
    1380           0 :       if (ZV_equal0(W))
    1381           0 :         return gerepilecopy(ltop, shallowtrans(V));
    1382           0 :       gerepileall(av, 2, &V, &W);
    1383             :     }
    1384           0 :     set_avma(btop);
    1385             :   }
    1386           0 :   return NULL;
    1387             : }
    1388             : 
    1389             : GEN
    1390           0 : zMs_ZC_mul(GEN M, GEN B)
    1391             : {
    1392             :   long i, j;
    1393           0 :   long n = lg(B)-1;
    1394           0 :   GEN V = zerocol(n);
    1395           0 :   for (i = 1; i <= n; ++i)
    1396           0 :     if (signe(gel(B, i)))
    1397             :     {
    1398           0 :       GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1399           0 :       long l = lg(C);
    1400           0 :       for (j = 1; j < l; ++j)
    1401             :       {
    1402           0 :         long k = C[j];
    1403           0 :         switch(E[j])
    1404             :         {
    1405           0 :         case 1:
    1406           0 :           gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
    1407           0 :           break;
    1408           0 :         case -1:
    1409           0 :           gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
    1410           0 :           break;
    1411           0 :         default:
    1412           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]));
    1413           0 :           break;
    1414             :         }
    1415             :       }
    1416             :     }
    1417           0 :   return V;
    1418             : }
    1419             : 
    1420             : GEN
    1421           0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
    1422             : 
    1423             : GEN
    1424       74577 : ZV_zMs_mul(GEN B, GEN M)
    1425             : {
    1426             :   long i, j;
    1427       74577 :   long m = lg(M)-1;
    1428       74577 :   GEN V = cgetg(m+1,t_VEC);
    1429    38974016 :   for (i = 1; i <= m; ++i)
    1430             :   {
    1431    38899439 :     GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1432    38899439 :     long l = lg(C);
    1433             :     GEN z;
    1434    38899439 :     if (l == 1)
    1435             :     {
    1436           0 :       gel(V,i) = gen_0;
    1437           0 :       continue;
    1438             :     }
    1439    38899439 :     z = mulis(gel(B, C[1]), E[1]);
    1440   327098202 :     for (j = 2; j < l; ++j)
    1441             :     {
    1442   288198763 :       long k = C[j];
    1443   288198763 :       switch(E[j])
    1444             :       {
    1445   194218066 :       case 1:
    1446   194218066 :         z = addii(z, gel(B,k));
    1447   194218066 :         break;
    1448    78587783 :       case -1:
    1449    78587783 :         z = subii(z, gel(B,k));
    1450    78587783 :         break;
    1451    15392914 :       default:
    1452    15392914 :         z = addii(z, mulis(gel(B,k), E[j]));
    1453    15392914 :         break;
    1454             :       }
    1455             :     }
    1456    38899439 :     gel(V,i) = z;
    1457             :   }
    1458       74577 :   return V;
    1459             : }
    1460             : 
    1461             : GEN
    1462         119 : FpV_FpMs_mul(GEN B, GEN M, GEN p) { return FpV_red(ZV_zMs_mul(B, M), p); }
    1463             : 
    1464             : GEN
    1465           0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
    1466             : {
    1467           0 :   pari_sp av = avma, av2;
    1468           0 :   GEN xi, xb, pi = gen_1, P;
    1469             :   long i;
    1470           0 :   if (!C) {
    1471           0 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    1472           0 :     if (!C) pari_err_INV("ZlM_gauss", a);
    1473             :   }
    1474           0 :   P = utoipos(p);
    1475           0 :   av2 = avma;
    1476           0 :   xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1477           0 :   xb = Flm_to_ZM(xi);
    1478           0 :   for (i = 2; i <= e; i++)
    1479             :   {
    1480           0 :     pi = muliu(pi, p); /* = p^(i-1) */
    1481           0 :     b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
    1482           0 :     if (gc_needed(av,2))
    1483             :     {
    1484           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
    1485           0 :       gerepileall(av2,3, &pi,&b,&xb);
    1486             :     }
    1487           0 :     xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1488           0 :     xb = ZM_add(xb, nm_Z_mul(xi, pi));
    1489             :   }
    1490           0 :   return gerepileupto(av, xb);
    1491             : }
    1492             : 
    1493             : struct wrapper_modp_s {
    1494             :   GEN (*f)(void*, GEN);
    1495             :   void *E;
    1496             :   GEN p;
    1497             : };
    1498             : 
    1499             : /* compute f(x) mod p */
    1500             : static GEN
    1501       74458 : wrap_relcomb_modp(void *E, GEN x)
    1502             : {
    1503       74458 :   struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
    1504       74458 :   return FpC_red(W->f(W->E, x), W->p);
    1505             : }
    1506             : static GEN
    1507           0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
    1508             : 
    1509             : static GEN
    1510       74458 : wrap_relker(void*E, GEN x) { return ZV_zMs_mul(x, (GEN)E); }
    1511             : 
    1512             : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
    1513             : GEN
    1514           0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
    1515             : {
    1516             :   struct wrapper_modp_s W;
    1517           0 :   pari_sp av = avma;
    1518           0 :   GEN xb, xi, pi = gen_1;
    1519             :   long i;
    1520           0 :   W.E = E;
    1521           0 :   W.f = f;
    1522           0 :   W.p = p;
    1523           0 :   xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
    1524           0 :   if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
    1525           0 :   xb = xi;
    1526           0 :   for (i = 2; i <= e; i++)
    1527             :   {
    1528           0 :     pi = mulii(pi, p); /* = p^(i-1) */
    1529           0 :     B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
    1530           0 :     if (gc_needed(av,2))
    1531             :     {
    1532           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
    1533           0 :       gerepileall(av,3, &pi,&B,&xb);
    1534             :     }
    1535           0 :     xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
    1536           0 :     if (!xi) return NULL;
    1537           0 :     if (typ(xi) == t_VEC) return gerepileupto(av, xi);
    1538           0 :     xb = ZC_add(xb, ZC_Z_mul(xi, pi));
    1539             :   }
    1540           0 :   return gerepileupto(av, xb);
    1541             : }
    1542             : 
    1543             : static GEN
    1544       24929 : vecprow(GEN A, GEN prow)
    1545             : {
    1546       24929 :   return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
    1547             : }
    1548             : 
    1549             : /* Solve the equation MX = A. Return either a solution as a t_COL,
    1550             :  * or the index of a column which is linearly dependent from the others as a
    1551             :  * t_VECSMALL with a single component. */
    1552             : GEN
    1553           0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
    1554             : {
    1555           0 :   pari_sp av = avma;
    1556             :   GEN pcol, prow;
    1557           0 :   long nbi=lg(M)-1, lR;
    1558             :   long i, n;
    1559             :   GEN Mp, Ap, Rp;
    1560             :   pari_timer ti;
    1561           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1562           0 :   RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
    1563           0 :   if (!pcol) return gc_NULL(av);
    1564           0 :   if (DEBUGLEVEL)
    1565           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
    1566           0 :   n = lg(pcol)-1;
    1567           0 :   Mp = cgetg(n+1, t_MAT);
    1568           0 :   for(i=1; i<=n; i++)
    1569           0 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1570           0 :   Ap = zCs_to_ZC(vecprow(A, prow), n);
    1571           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1572           0 :   Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
    1573           0 :   if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
    1574           0 :   if (!Rp) return gc_NULL(av);
    1575           0 :   lR = lg(Rp)-1;
    1576           0 :   if (typ(Rp) == t_COL)
    1577             :   {
    1578           0 :     GEN R = zerocol(nbi+1);
    1579           0 :     for (i=1; i<=lR; i++)
    1580           0 :        gel(R,pcol[i]) = gel(Rp,i);
    1581           0 :     return gerepilecopy(av, R);
    1582             :   }
    1583           0 :   for (i = 1; i <= lR; ++i)
    1584           0 :     if (signe(gel(Rp, i)))
    1585           0 :       return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
    1586             :   return NULL; /* LCOV_EXCL_LINE */
    1587             : }
    1588             : 
    1589             : GEN
    1590           0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
    1591             : {
    1592           0 :   return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1593             : }
    1594             : 
    1595             : GEN
    1596           0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
    1597             : {
    1598             :   GEN res;
    1599           0 :   pari_CATCH(e_INV)
    1600             :   {
    1601           0 :     GEN E = pari_err_last();
    1602           0 :     GEN x = err_get_compo(E,2);
    1603           0 :     if (typ(x) != t_INTMOD) pari_err(0,E);
    1604           0 :     if (DEBUGLEVEL)
    1605           0 :       pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
    1606           0 :     res = NULL;
    1607             :   } pari_TRY {
    1608           0 :     res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1609           0 :   } pari_ENDCATCH
    1610           0 :   return res;
    1611             : }
    1612             : 
    1613             : static GEN
    1614         119 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
    1615             : {
    1616         119 :   long i, j, oldf = 0, f = 0;
    1617         119 :   long lrow = lg(prow), lM = lg(M);
    1618         119 :   GEN W = const_vecsmall(lM-1,1);
    1619         119 :   GEN R = cgetg(lrow, t_VEC);
    1620      235298 :   for (i=1; i<lrow; i++)
    1621      235179 :     gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
    1622             :   do
    1623             :   {
    1624         393 :     oldf = f;
    1625      436630 :     for(i=1; i<lM; i++)
    1626      436237 :       if (W[i])
    1627             :       {
    1628      147866 :         GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
    1629      147866 :         long c=0, cj=0, lF = lg(F);
    1630     1294916 :         for(j=1; j<lF; j++)
    1631     1147050 :           if (!gel(R,F[j])) { c++; cj=j; }
    1632      147866 :         if (c>=2) continue;
    1633      127154 :         if (c==1)
    1634             :         {
    1635       41415 :           pari_sp av = avma;
    1636       41415 :           GEN s = gen_0;
    1637      359135 :           for(j=1; j<lF; j++)
    1638      317720 :             if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
    1639             :           /* s /= -E[cj] mod p */
    1640       41415 :           s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
    1641       41415 :           gel(R,F[cj]) = gerepileuptoint(av, s);
    1642       41415 :           f++;
    1643             :         }
    1644      127154 :         W[i]=0;
    1645             :       }
    1646         393 :   } while(oldf!=f);
    1647      235298 :   for (i=1; i<lrow; i++)
    1648      235179 :     if (!gel(R,i)) gel(R,i) = gen_0;
    1649         119 :   return R;
    1650             : }
    1651             : 
    1652             : /* Return a linear form Y such that YM is essentially 0 */
    1653             : GEN
    1654         119 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
    1655             : {
    1656         119 :   pari_sp av = avma, av2;
    1657             :   GEN pcol, prow;
    1658             :   long i, n;
    1659             :   GEN Mp, B, MB, R, Rp;
    1660             :   pari_timer ti;
    1661             :   struct wrapper_modp_s W;
    1662         119 :   if (DEBUGLEVEL) timer_start(&ti);
    1663         119 :   RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
    1664         119 :   if (!pcol) return gc_NULL(av);
    1665         119 :   if (DEBUGLEVEL)
    1666           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
    1667         119 :   n = lg(pcol)-1;
    1668         119 :   Mp = cgetg(n+1, t_MAT);
    1669       25048 :   for(i=1; i<=n; i++)
    1670       24929 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1671         119 :   W.E = (void*) Mp;
    1672         119 :   W.f = wrap_relker;
    1673         119 :   W.p = p;
    1674         119 :   av2 = avma;
    1675             :   for(;;)
    1676             :   {
    1677         119 :     set_avma(av2);
    1678         119 :     B = random_FpV(n, p);
    1679         119 :     MB = FpV_FpMs_mul(B, Mp, p);
    1680         119 :     if (DEBUGLEVEL) timer_start(&ti);
    1681         119 :     pari_CATCH(e_INV)
    1682             :     {
    1683           0 :       GEN E = pari_err_last();
    1684           0 :       GEN x = err_get_compo(E,2);
    1685           0 :       if (typ(x) != t_INTMOD) pari_err(0,E);
    1686           0 :       if (DEBUGLEVEL)
    1687           0 :         pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
    1688           0 :       Rp = NULL;
    1689             :     } pari_TRY {
    1690         119 :       Rp = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, MB, p);
    1691         119 :     } pari_ENDCATCH
    1692         119 :     if (!Rp) continue;
    1693         119 :     if (typ(Rp)==t_COL)
    1694             :     {
    1695         119 :       Rp = FpC_sub(Rp,B,p);
    1696         119 :       if (ZV_equal0(Rp)) continue;
    1697             :     }
    1698         119 :     R = FpMs_structelim_back(M, Rp, prow, p);
    1699         119 :     if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
    1700         119 :     return gerepilecopy(av, R);
    1701             :   }
    1702             : }
    1703             : 
    1704             : GEN
    1705          42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
    1706             : {
    1707          42 :   return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
    1708             : }

Generated by: LCOV version 1.13