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.0 lcov report (development 23886-801e231f7) Lines: 700 903 77.5 %
Date: 2019-05-27 05:43:56 Functions: 105 127 82.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. 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    11072999 : FpC_red(GEN x, GEN p)
      25    11072999 : { 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      207830 : FpV_red(GEN x, GEN p)
      30      207830 : { pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
      31             : 
      32             : GEN
      33     1378073 : FpC_center(GEN x, GEN p, GEN pov2)
      34     1378073 : { 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      142275 : { 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     2215399 : FpM_red(GEN x, GEN p)
      56     2215399 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
      57             : 
      58             : GEN
      59      164472 : FpM_center(GEN x, GEN p, GEN pov2)
      60      164472 : { 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          14 :   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        9947 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
      71             : 
      72             : GEN
      73         126 : random_FpV(long d, GEN p)
      74             : {
      75             :   long i;
      76         126 :   GEN y = cgetg(d+1,t_VEC);
      77         126 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
      78         126 :   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         413 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
      87         413 :   return y;
      88             : }
      89             : 
      90             : GEN
      91        4177 : random_Flv(long n, ulong p)
      92             : {
      93        4177 :   GEN y = cgetg(n+1, t_VECSMALL);
      94             :   long i;
      95        4177 :   for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
      96        4177 :   return y;
      97             : }
      98             : 
      99             : /********************************************************************/
     100             : /**                                                                **/
     101             : /**                           ADD, SUB                             **/
     102             : /**                                                                **/
     103             : /********************************************************************/
     104             : GEN
     105      186102 : FpC_add(GEN x, GEN y, GEN p)
     106      186102 : { 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      589804 : Flv_add(GEN x, GEN y, ulong p)
     118             : {
     119      589804 :   if (p==2)
     120       77300 :     pari_APPLY_ulong(x[i]^y[i])
     121             :   else
     122      512504 :     pari_APPLY_ulong(Fl_add(x[i], y[i], p))
     123             : }
     124             : 
     125             : void
     126      175071 : Flv_add_inplace(GEN x, GEN y, ulong p)
     127             : {
     128      175071 :   long i, l = lg(x);
     129      175071 :   if (p==2)
     130      171935 :     for (i = 1; i < l; i++) x[i] ^= y[i];
     131             :   else
     132        3136 :     for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
     133      175071 : }
     134             : 
     135             : ulong
     136        3823 : Flv_sum(GEN x, ulong p)
     137             : {
     138        3823 :   long i, l = lg(x);
     139        3823 :   ulong s = 0;
     140        3823 :   if (p==2)
     141        3823 :     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        3823 :   return s;
     145             : }
     146             : 
     147             : GEN
     148      481649 : FpC_sub(GEN x, GEN y, GEN p)
     149      481649 : { 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           0 : FpM_sub(GEN x, GEN y, GEN p)
     157           0 : { pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
     158             : 
     159             : GEN
     160    77915292 : Flv_sub(GEN x, GEN y, ulong p)
     161    77915292 : { 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       14572 : Flm_Fl_add(GEN x, ulong y, ulong p)
     172             : {
     173       14572 :   long l = lg(x), i, j;
     174       14572 :   GEN z = cgetg(l,t_MAT);
     175             : 
     176       14572 :   if (l==1) return z;
     177       14362 :   if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
     178       55280 :   for (i=1; i<l; i++)
     179             :   {
     180       40918 :     GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
     181       40918 :     gel(z,i) = zi;
     182       40918 :     for (j=1; j<l; j++) zi[j] = xi[j];
     183       40918 :     zi[i] = Fl_add(zi[i], y, p);
     184             :   }
     185       14362 :   return z;
     186             : }
     187             : 
     188             : GEN
     189        8715 : Flm_add(GEN x, GEN y, ulong p)
     190        8715 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
     191             : 
     192             : GEN
     193     8640860 : Flm_sub(GEN x, GEN y, ulong p)
     194     8640860 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
     195             : 
     196             : /********************************************************************/
     197             : /**                                                                **/
     198             : /**                           MULTIPLICATION                       **/
     199             : /**                                                                **/
     200             : /********************************************************************/
     201             : GEN
     202      737708 : FpC_Fp_mul(GEN x, GEN y, GEN p)
     203      737708 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
     204             : 
     205             : GEN
     206      553734 : Flv_Fl_mul(GEN x, ulong y, ulong p)
     207      553734 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
     208             : 
     209             : GEN
     210         464 : Flv_Fl_div(GEN x, ulong y, ulong p)
     211             : {
     212         464 :   return Flv_Fl_mul(x, Fl_inv(y, p), p);
     213             : }
     214             : void
     215           0 : Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
     216             : {
     217           0 :   Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
     218           0 : }
     219             : GEN
     220      124049 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
     221      124049 :   long i, j, h, l = lg(X);
     222      124049 :   GEN A = cgetg(l, t_MAT);
     223      124049 :   if (l == 1) return A;
     224      124049 :   h = lgcols(X);
     225     1149344 :   for (j=1; j<l; j++)
     226             :   {
     227     1025295 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     228     1025295 :     for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
     229     1025295 :     gel(A,j) = a;
     230             :   }
     231      124049 :   return A;
     232             : }
     233             : 
     234             : /* x *= y */
     235             : void
     236     1801631 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
     237             : {
     238             :   long i;
     239     1801631 :   for (i=1;i<=l;i++) x[i] = Fl_mul(x[i], y, p);
     240     1801631 : }
     241             : void
     242           0 : Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
     243           0 : { Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
     244             : 
     245             : /* set y *= x */
     246             : void
     247           0 : Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
     248             : {
     249           0 :   long i, j, m, l = lg(y);
     250           0 :   if (l == 1) return;
     251           0 :   m = lgcols(y);
     252           0 :   if (HIGHWORD(x | p))
     253           0 :     for(j=1; j<l; j++)
     254           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
     255             :   else
     256           0 :     for(j=1; j<l; j++)
     257           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
     258             : }
     259             : 
     260             : /* return x * y */
     261             : static GEN
     262     4810038 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
     263             : {
     264     4810038 :   long i, j, m, l = lg(y);
     265     4810038 :   GEN z = cgetg(l, t_MAT);
     266     4809148 :   if (l == 1) return z;
     267     4809148 :   m = lgcols(y);
     268    49603607 :   for(j=1; j<l; j++) {
     269    44793295 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     270    44700324 :     for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
     271             :   }
     272     4810312 :   return z;
     273             : }
     274             : 
     275             : /* return x * y */
     276             : static GEN
     277     3071058 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
     278             : {
     279     3071058 :   long i, j, m, l = lg(y);
     280     3071058 :   GEN z = cgetg(l, t_MAT);
     281     3071058 :   if (l == 1) return z;
     282     3071058 :   m = lgcols(y);
     283    18567948 :   for(j=1; j<l; j++) {
     284    15496890 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     285    15496890 :     for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
     286             :   }
     287     3071058 :   return z;
     288             : }
     289             : 
     290             : /* return x * y */
     291             : GEN
     292     7847640 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
     293             : {
     294     7847640 :   if (HIGHWORD(p))
     295     4810161 :     return Flm_Fl_mul_pre_i(y, x, p, pi);
     296             :   else
     297     3037479 :     return Flm_Fl_mul_OK(y, x, p);
     298             : }
     299             : 
     300             : /* return x * y */
     301             : GEN
     302       33579 : Flm_Fl_mul(GEN y, ulong x, ulong p)
     303             : {
     304       33579 :   if (HIGHWORD(p))
     305           0 :     return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
     306             :   else
     307       33579 :     return Flm_Fl_mul_OK(y, x, p);
     308             : }
     309             : 
     310             : GEN
     311     1697074 : Flv_neg(GEN x, ulong p)
     312     1697074 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
     313             : 
     314             : void
     315        5764 : Flv_neg_inplace(GEN v, ulong p)
     316             : {
     317             :   long i;
     318      184358 :   for (i = 1; i < lg(v); ++i)
     319      178594 :     v[i] = Fl_neg(v[i], p);
     320        5764 : }
     321             : 
     322             : GEN
     323      136463 : Flm_neg(GEN x, ulong p)
     324      136463 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
     325             : 
     326             : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
     327             : static GEN
     328    12921770 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
     329             : {
     330    12921770 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1));
     331             :   long k;
     332   196540234 :   for (k = 2; k < lx; k++)
     333             :   {
     334   183618181 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     335   183610751 :     if (signe(t)) c = addii(c, t);
     336             :   }
     337    12922053 :   return c;
     338             : }
     339             : 
     340             : static long
     341    68774300 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
     342             : {
     343             :   long k;
     344    68774300 :   long c = coeff(x,i,1) * y[1];
     345  1042862765 :   for (k = 2; k < lx; k++)
     346   974088465 :     c += coeff(x,i,k) * y[k];
     347    68774300 :   return c;
     348             : }
     349             : 
     350             : GEN
     351     4606882 : zm_zc_mul(GEN x, GEN y)
     352             : {
     353     4606882 :   long lx = lg(x), l, i;
     354             :   GEN z;
     355     4606882 :   if (lx == 1) return cgetg(1, t_VECSMALL);
     356     4606882 :   l = lg(gel(x,1));
     357     4606882 :   z = cgetg(l,t_VECSMALL);
     358     4606882 :   for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
     359     4606882 :   return z;
     360             : }
     361             : 
     362             : GEN
     363        1771 : zm_mul(GEN x, GEN y)
     364             : {
     365        1771 :   long i,j,lx=lg(x), ly=lg(y);
     366             :   GEN z;
     367        1771 :   if (ly==1) return cgetg(1,t_MAT);
     368        1771 :   z = cgetg(ly,t_MAT);
     369        1771 :   if (lx==1)
     370             :   {
     371           0 :     for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
     372           0 :     return z;
     373             :   }
     374       24402 :   for (j=1; j<ly; j++)
     375       22631 :     gel(z,j) = zm_zc_mul(x, gel(y,j));
     376        1771 :   return z;
     377             : }
     378             : 
     379             : static ulong
     380   339067378 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
     381             : {
     382   339067378 :   ulong c = ucoeff(x,i,1) * uel(y,1);
     383             :   long k;
     384  4793879579 :   for (k = 2; k < lx; k++) {
     385  4454812201 :     c += ucoeff(x,i,k) * uel(y,k);
     386  4454812201 :     if (c & HIGHBIT) c %= p;
     387             :   }
     388   339067378 :   return c % p;
     389             : }
     390             : 
     391             : static ulong
     392   173903534 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
     393             : {
     394             :   ulong l0, l1, v1, h0, h1;
     395   173903534 :   long k = 1;
     396             :   LOCAL_OVERFLOW;
     397             :   LOCAL_HIREMAINDER;
     398   173903534 :   l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
     399  2146383519 :   while (++k < lx) {
     400  1798576451 :     l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
     401  1798576451 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     402             :   }
     403   173903534 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     404     5088530 :   else return remlll_pre(v1, h1, l1, p, pi);
     405             : }
     406             : 
     407             : static GEN
     408      144346 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
     409             : {
     410             :   long i,j;
     411      144346 :   GEN z = NULL;
     412             : 
     413     1506060 :   for (j=1; j<lx; j++)
     414             :   {
     415     1361714 :     if (!y[j]) continue;
     416      390420 :     if (!z) z = Flv_copy(gel(x,j));
     417      264890 :     else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
     418             :   }
     419      144346 :   if (!z) z = zero_zv(l-1);
     420      144346 :   return z;
     421             : }
     422             : 
     423             : static GEN
     424     1076821 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
     425             : {
     426     1076821 :   GEN z = cgetg(l,t_COL);
     427             :   long i;
     428    12980610 :   for (i = 1; i < l; i++)
     429             :   {
     430    11903815 :     pari_sp av = avma;
     431    11903815 :     GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
     432    11904108 :     gel(z,i) = gerepileuptoint(av, modii(c,p));
     433             :   }
     434     1076795 :   return z;
     435             : }
     436             : 
     437             : static void
     438    37238131 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
     439             : {
     440             :   long i;
     441    37238131 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
     442    37254911 : }
     443             : static GEN
     444    34137793 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
     445             : {
     446    34137793 :   GEN z = cgetg(l,t_VECSMALL);
     447    34137282 :   __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     448    34138028 :   return z;
     449             : }
     450             : 
     451             : static void
     452    32844376 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     453             : {
     454             :   long i;
     455    32844376 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
     456    33046056 : }
     457             : static GEN
     458    28520669 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     459             : {
     460    28520669 :   GEN z = cgetg(l,t_VECSMALL);
     461    28461298 :   __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     462    28533755 :   return z;
     463             : }
     464             : 
     465             : GEN
     466      581723 : FpM_mul(GEN x, GEN y, GEN p)
     467             : {
     468      581723 :   long lx=lg(x), ly=lg(y);
     469             :   GEN z;
     470      581723 :   pari_sp av = avma;
     471      581723 :   if (ly==1) return cgetg(1,t_MAT);
     472      581723 :   if (lx==1) return zeromat(0, ly-1);
     473      581723 :   if (lgefint(p) == 3)
     474             :   {
     475      554514 :     ulong pp = uel(p,2);
     476      554514 :     if (pp == 2)
     477             :     {
     478      100316 :       x = ZM_to_F2m(x);
     479      100316 :       y = ZM_to_F2m(y);
     480      100316 :       z = F2m_to_ZM(F2m_mul(x,y));
     481             :     }
     482             :     else
     483             :     {
     484      454198 :       x = ZM_to_Flm(x, pp);
     485      454198 :       y = ZM_to_Flm(y, pp);
     486      454198 :       z = Flm_to_ZM(Flm_mul(x,y, pp));
     487             :     }
     488             :   } else
     489       27209 :     z = FpM_red(ZM_mul(x, y), p);
     490      581723 :   return gerepileupto(av, z);
     491             : }
     492             : 
     493             : static GEN
     494    11445232 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     495             : {
     496             :   long j;
     497    11445232 :   GEN z = cgetg(ly, t_MAT);
     498    11442252 :   if (SMALL_ULONG(p))
     499    41814332 :     for (j=1; j<ly; j++)
     500    34005454 :       gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
     501             :   else
     502    32162521 :     for (j=1; j<ly; j++)
     503    28525630 :       gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
     504    11445769 :   return z;
     505             : }
     506             : 
     507             : /*
     508             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     509             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     510             : */
     511             : static void
     512        7424 : add_slices_ip(long m, long n,
     513             :            GEN A, long ma, long da, long na, long ea,
     514             :            GEN B, long mb, long db, long nb, long eb,
     515             :            GEN M, long dy, long dx, ulong p)
     516             : {
     517        7424 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     518             :   GEN C;
     519             : 
     520      277078 :   for (j = 1; j <= min_e; j++) {
     521      269654 :     C = gel(M, j + dx) + dy;
     522    11953946 :     for (i = 1; i <= min_d; i++)
     523    23368584 :       uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
     524    11684292 :                         ucoeff(B, mb + i, nb + j), p);
     525      275954 :     for (; i <= da; i++)
     526        6300 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     527      269654 :     for (; i <= db; i++)
     528           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     529      269654 :     for (; i <= m; i++)
     530           0 :       uel(C, i) = 0;
     531             :   }
     532        8016 :   for (; j <= ea; j++) {
     533         592 :     C = gel(M, j + dx) + dy;
     534       27488 :     for (i = 1; i <= da; i++)
     535       26896 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     536         592 :     for (; i <= m; i++)
     537           0 :       uel(C, i) = 0;
     538             :   }
     539        7424 :   for (; j <= eb; j++) {
     540           0 :     C = gel(M, j + dx) + dy;
     541           0 :     for (i = 1; i <= db; i++)
     542           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     543           0 :     for (; i <= m; i++)
     544           0 :       uel(C, i) = 0;
     545             :   }
     546        7424 :   for (; j <= n; j++) {
     547           0 :     C = gel(M, j + dx) + dy;
     548           0 :     for (i = 1; i <= m; i++)
     549           0 :       uel(C, i) = 0;
     550             :   }
     551        7424 : }
     552             : 
     553             : static GEN
     554        3712 : add_slices(long m, long n,
     555             :            GEN A, long ma, long da, long na, long ea,
     556             :            GEN B, long mb, long db, long nb, long eb, ulong p)
     557             : {
     558             :   GEN M;
     559             :   long j;
     560        3712 :   M = cgetg(n + 1, t_MAT);
     561      137166 :   for (j = 1; j <= n; j++)
     562      133454 :     gel(M, j) = cgetg(m + 1, t_VECSMALL);
     563        3712 :   add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
     564        3712 :   return M;
     565             : }
     566             : 
     567             : /*
     568             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     569             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     570             : */
     571             : static GEN
     572        6496 : subtract_slices(long m, long n,
     573             :                 GEN A, long ma, long da, long na, long ea,
     574             :                 GEN B, long mb, long db, long nb, long eb, ulong p)
     575             : {
     576        6496 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     577        6496 :   GEN M = cgetg(n + 1, t_MAT), C;
     578             : 
     579      235060 :   for (j = 1; j <= min_e; j++) {
     580      228564 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     581     9725482 :     for (i = 1; i <= min_d; i++)
     582    18993836 :       uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
     583     9496918 :                         ucoeff(B, mb + i, nb + j), p);
     584      247150 :     for (; i <= da; i++)
     585       18586 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     586      234432 :     for (; i <= db; i++)
     587        5868 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     588      228564 :     for (; i <= m; i++)
     589           0 :       uel(C, i) = 0;
     590             :   }
     591        6496 :   for (; j <= ea; j++) {
     592           0 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     593           0 :     for (i = 1; i <= da; i++)
     594           0 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     595           0 :     for (; i <= m; i++)
     596           0 :       uel(C, i) = 0;
     597             :   }
     598        6640 :   for (; j <= eb; j++) {
     599         144 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     600        5616 :     for (i = 1; i <= db; i++)
     601        5472 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     602         144 :     for (; i <= m; i++)
     603           0 :       uel(C, i) = 0;
     604             :   }
     605        6640 :   for (; j <= n; j++)
     606         144 :     gel(M, j) = zero_Flv(m);
     607        6496 :   return M;
     608             : }
     609             : 
     610             : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
     611             : 
     612             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     613             : static GEN
     614         928 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
     615             : {
     616             :   pari_sp av;
     617         928 :   long m1 = (m + 1)/2, m2 = m/2,
     618         928 :     n1 = (n + 1)/2, n2 = n/2,
     619         928 :     p1 = (p + 1)/2, p2 = p/2;
     620             :   GEN A11, A12, A22, B11, B21, B22,
     621             :     S1, S2, S3, S4, T1, T2, T3, T4,
     622             :     M1, M2, M3, M4, M5, M6, M7,
     623             :     V1, V2, V3, C;
     624             :   long j;
     625         928 :   C = cgetg(p + 1, t_MAT);
     626       69324 :   for (j = 1; j <= p; j++)
     627       68396 :     gel(C, j) = cgetg(m + 1, t_VECSMALL);
     628         928 :   av = avma;
     629         928 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
     630         928 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
     631         928 :   M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
     632         928 :   if (gc_needed(av, 1))
     633           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     634         928 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
     635         928 :   if (gc_needed(av, 1))
     636           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     637         928 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
     638         928 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
     639         928 :   M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
     640         928 :   if (gc_needed(av, 1))
     641           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     642         928 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
     643         928 :   if (gc_needed(av, 1))
     644           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     645         928 :   A11 = matslice(A, 1, m1, 1, n1);
     646         928 :   B11 = matslice(B, 1, n1, 1, p1);
     647         928 :   M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
     648         928 :   if (gc_needed(av, 1))
     649           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     650         928 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     651         928 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     652         928 :   M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
     653         928 :   if (gc_needed(av, 1))
     654           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     655         928 :   add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
     656         928 :   if (gc_needed(av, 1))
     657           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy M4 */
     658         928 :   M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
     659         928 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
     660         928 :   if (gc_needed(av, 1))
     661           0 :     gerepileall(av, 6, &M2, &T3, &M3, &M1, &M5, &S4);  /* destroy S3 */
     662         928 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
     663         928 :   if (gc_needed(av, 1))
     664           0 :     gerepileall(av, 6, &M2, &M3, &M1, &M5, &S4, &T4);  /* destroy T3 */
     665         928 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
     666         928 :   if (gc_needed(av, 1))
     667           0 :     gerepileall(av, 5, &M2, &M3, &S4, &T4, &V1);  /* destroy M1, M5 */
     668         928 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     669         928 :   M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
     670         928 :   if (gc_needed(av, 1))
     671           0 :     gerepileall(av, 5, &M2, &M3, &T4, &V1, &M6);  /* destroy S4, B22 */
     672         928 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     673         928 :   M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
     674         928 :   if (gc_needed(av, 1))
     675           0 :     gerepileall(av, 5, &M2, &M3, &V1, &M6, &M7);  /* destroy A22, T4 */
     676         928 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
     677         928 :   add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
     678         928 :   if (gc_needed(av, 1))
     679           0 :     gerepileall(av, 4, &M2, &M3, &V1, &M7);  /* destroy V3, M6 */
     680         928 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
     681         928 :   if (gc_needed(av, 1))
     682           0 :     gerepileall(av, 3, &M3, &M7, &V2);  /* destroy V1, M2 */
     683         928 :   add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
     684         928 :   add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
     685         928 :   set_avma(av); return C;
     686             : }
     687             : 
     688             : /* Strassen-Winograd used for dim >= ZM_sw_bound */
     689             : static GEN
     690    11445960 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     691             : {
     692    11445960 :   ulong e = expu(p);
     693             : #ifdef LONG_IS_64BIT
     694     9404350 :   long ZM_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
     695             : #else
     696     2041737 :   long ZM_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
     697             : #endif
     698    11446087 :   if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
     699    11445159 :     return Flm_mul_classical(x, y, l, lx, ly, p, pi);
     700             :   else
     701         928 :     return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
     702             : }
     703             : 
     704             : GEN
     705     4708453 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     706             : {
     707     4708453 :   long lx=lg(x), ly=lg(y);
     708     4708453 :   if (ly==1) return cgetg(1,t_MAT);
     709     4708453 :   if (lx==1) return zero_Flm(0, ly-1);
     710     4708453 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
     711             : }
     712             : 
     713             : GEN
     714     6730741 : Flm_mul(GEN x, GEN y, ulong p)
     715             : {
     716     6730741 :   long lx=lg(x), ly=lg(y);
     717     6730741 :   if (ly==1) return cgetg(1,t_MAT);
     718     6730580 :   if (lx==1) return zero_Flm(0, ly-1);
     719     6730580 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
     720             : }
     721             : 
     722             : struct _Flm
     723             : {
     724             :   ulong p;
     725             :   long n;
     726             : };
     727             : 
     728             : static GEN
     729        8421 : _Flm_mul(void *E , GEN x, GEN y)
     730        8421 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
     731             : static GEN
     732       43393 : _Flm_sqr(void *E, GEN x)
     733       43393 : { return Flm_mul(x,x,((struct _Flm*)E)->p); }
     734             : static GEN
     735         266 : _Flm_one(void *E)
     736         266 : { return matid_Flm(((struct _Flm*)E)->n); }
     737             : GEN
     738       24528 : Flm_powu(GEN x, ulong n, ulong p)
     739             : {
     740             :   struct _Flm d;
     741       24528 :   if (!n) return matid(lg(x)-1);
     742       24528 :   d.p = p;
     743       24528 :   return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
     744             : }
     745             : GEN
     746         266 : Flm_powers(GEN x, ulong n, ulong p)
     747             : {
     748             :   struct _Flm d;
     749         266 :   d.p = p;
     750         266 :   d.n = lg(x)-1;
     751         266 :   return gen_powers(x, n, 1, (void*)&d,
     752             :                     &_Flm_sqr, &_Flm_mul, &_Flm_one);
     753             : }
     754             : 
     755             : static GEN
     756           0 : _FpM_mul(void *p , GEN x, GEN y)
     757           0 : { return FpM_mul(x,y,(GEN)p); }
     758             : static GEN
     759           0 : _FpM_sqr(void *p, GEN x)
     760           0 : { return FpM_mul(x,x,(GEN)p); }
     761             : GEN
     762           0 : FpM_powu(GEN x, ulong n, GEN p)
     763             : {
     764           0 :   if (!n) return matid(lg(x)-1);
     765           0 :   if (lgefint(p) == 3)
     766             :   {
     767           0 :     pari_sp av = avma;
     768           0 :     ulong pp = uel(p,2);
     769             :     GEN z;
     770           0 :     if (pp == 2)
     771           0 :       z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
     772             :     else
     773           0 :       z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
     774           0 :     return gerepileupto(av, z);
     775             :   }
     776           0 :   return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
     777             : }
     778             : 
     779             : /*Multiple a column vector by a line vector to make a matrix*/
     780             : GEN
     781          14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
     782             : {
     783          14 :   long i,j, lx=lg(x), ly=lg(y);
     784             :   GEN z;
     785          14 :   if (ly==1) return cgetg(1,t_MAT);
     786          14 :   z = cgetg(ly, t_MAT);
     787          49 :   for (j=1; j < ly; j++)
     788             :   {
     789          35 :     GEN zj = cgetg(lx,t_VECSMALL);
     790         112 :     for (i=1; i<lx; i++)
     791          77 :       uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
     792          35 :     gel(z,j) = zj;
     793             :   }
     794          14 :   return z;
     795             : }
     796             : 
     797             : /*Multiple a column vector by a line vector to make a matrix*/
     798             : GEN
     799        1877 : FpC_FpV_mul(GEN x, GEN y, GEN p)
     800             : {
     801        1877 :   long i,j, lx=lg(x), ly=lg(y);
     802             :   GEN z;
     803        1877 :   if (ly==1) return cgetg(1,t_MAT);
     804        1877 :   z = cgetg(ly,t_MAT);
     805       30557 :   for (j=1; j < ly; j++)
     806             :   {
     807       28680 :     GEN zj = cgetg(lx,t_COL);
     808       28680 :     for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
     809       28680 :     gel(z, j) = zj;
     810             :   }
     811        1877 :   return z;
     812             : }
     813             : 
     814             : /* Multiply a line vector by a column and return a scalar (t_INT) */
     815             : GEN
     816     1942598 : FpV_dotproduct(GEN x, GEN y, GEN p)
     817             : {
     818     1942598 :   long i, lx = lg(x);
     819             :   pari_sp av;
     820             :   GEN c;
     821     1942598 :   if (lx == 1) return gen_0;
     822     1942017 :   av = avma; c = mulii(gel(x,1),gel(y,1));
     823     1941927 :   for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
     824     1941981 :   return gerepileuptoint(av, modii(c,p));
     825             : }
     826             : GEN
     827           0 : FpV_dotsquare(GEN x, GEN p)
     828             : {
     829           0 :   long i, lx = lg(x);
     830             :   pari_sp av;
     831             :   GEN c;
     832           0 :   if (lx == 1) return gen_0;
     833           0 :   av = avma; c = sqri(gel(x,1));
     834           0 :   for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
     835           0 :   return gerepileuptoint(av, modii(c,p));
     836             : }
     837             : 
     838             : INLINE ulong
     839     2984439 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
     840             : {
     841     2984439 :   ulong c = uel(x,0) * uel(y,0);
     842             :   long k;
     843    22768903 :   for (k = 1; k < lx; k++) {
     844    19784464 :     c += uel(x,k) * uel(y,k);
     845    19784464 :     if (c & HIGHBIT) c %= p;
     846             :   }
     847     2984439 :   return c % p;
     848             : }
     849             : 
     850             : INLINE ulong
     851        8314 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
     852             : {
     853             :   ulong l0, l1, v1, h0, h1;
     854        8314 :   long i = 0;
     855             :   LOCAL_OVERFLOW;
     856             :   LOCAL_HIREMAINDER;
     857        8314 :   l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
     858      183673 :   while (++i < lx) {
     859      167045 :     l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
     860      167045 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     861             :   }
     862        8314 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     863           3 :   else return remlll_pre(v1, h1, l1, p, pi);
     864             : }
     865             : 
     866             : ulong
     867      322252 : Flv_dotproduct(GEN x, GEN y, ulong p)
     868             : {
     869      322252 :   long lx = lg(x)-1;
     870      322252 :   if (lx == 0) return 0;
     871      322252 :   if (SMALL_ULONG(p))
     872      322252 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     873             :   else
     874           0 :     return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
     875             : }
     876             : 
     877             : ulong
     878       57184 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
     879             : {
     880       57184 :   long lx = lg(x)-1;
     881       57184 :   if (lx == 0) return 0;
     882       57184 :   if (SMALL_ULONG(p))
     883       51692 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     884             :   else
     885        5492 :     return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
     886             : }
     887             : 
     888             : ulong
     889     2637616 : Flx_dotproduct(GEN x, GEN y, ulong p)
     890             : {
     891     2637616 :   long lx = minss(lgpol(x), lgpol(y));
     892     2637685 :   if (lx == 0) return 0;
     893     2613057 :   if (SMALL_ULONG(p))
     894     2610235 :     return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
     895             :   else
     896        2822 :     return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
     897             : }
     898             : 
     899             : GEN
     900     1076822 : FpM_FpC_mul(GEN x, GEN y, GEN p)
     901             : {
     902     1076822 :   long lx = lg(x);
     903     1076822 :   return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
     904             : }
     905             : GEN
     906      276792 : Flm_Flc_mul(GEN x, GEN y, ulong p)
     907             : {
     908      276792 :   long l, lx = lg(x);
     909      276792 :   if (lx==1) return cgetg(1,t_VECSMALL);
     910      276792 :   l = lgcols(x);
     911      276792 :   if (p==2)
     912      144346 :     return Flm_Flc_mul_i_2(x, y, lx, l);
     913      132446 :   else if (SMALL_ULONG(p))
     914      132332 :     return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
     915             :   else
     916         114 :     return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
     917             : }
     918             : 
     919             : GEN
     920        4816 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     921             : {
     922        4816 :   long l, lx = lg(x);
     923             :   GEN z;
     924        4816 :   if (lx==1) return cgetg(1,t_VECSMALL);
     925        4816 :   l = lgcols(x);
     926        4816 :   z = cgetg(l, t_VECSMALL);
     927        4815 :   if (SMALL_ULONG(p))
     928        4127 :     __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     929             :   else
     930         688 :     __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     931        4816 :   return z;
     932             : }
     933             : 
     934             : GEN
     935     7604287 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
     936             : {
     937     7604287 :   long l, lx = lg(x);
     938             :   GEN z;
     939     7604287 :   if (lx==1) return pol0_Flx(sv);
     940     7604287 :   l = lgcols(x);
     941     7601360 :   z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
     942     7603645 :   if (SMALL_ULONG(p))
     943     3104331 :     __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
     944             :   else
     945     4499314 :     __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
     946     7625956 :   return Flx_renormalize(z, l + 1);
     947             : }
     948             : 
     949             : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
     950             : GEN
     951      354805 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
     952             : {
     953      354805 :   long i, l, lx = lg(x);
     954             :   GEN z;
     955      354805 :   if (lx==1) return pol_0(v);
     956      354805 :   l = lgcols(x);
     957      354805 :   z = new_chunk(l+1);
     958      759425 :   for (i=l-1; i; i--)
     959             :   {
     960      604534 :     pari_sp av = avma;
     961      604534 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
     962      604534 :     p1 = modii(p1, p);
     963      604534 :     if (signe(p1))
     964             :     {
     965      199914 :       if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
     966      199914 :       gel(z,i+1) = gerepileuptoint(av, p1);
     967      199914 :       break;
     968             :     }
     969      404620 :     set_avma(av);
     970             :   }
     971      354805 :   if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
     972      199914 :   z[0] = evaltyp(t_POL) | evallg(i+2);
     973      199914 :   z[1] = evalsigne(1) | evalvarn(v);
     974      613320 :   for (; i; i--)
     975             :   {
     976      413406 :     pari_sp av = avma;
     977      413406 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
     978      413406 :     gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
     979             :   }
     980      199914 :   return z;
     981             : }
     982             : 
     983             : /********************************************************************/
     984             : /**                                                                **/
     985             : /**                           TRANSPOSITION                        **/
     986             : /**                                                                **/
     987             : /********************************************************************/
     988             : 
     989             : /* == zm_transpose */
     990             : GEN
     991      306691 : Flm_transpose(GEN x)
     992             : {
     993      306691 :   long i, dx, lx = lg(x);
     994             :   GEN y;
     995      306691 :   if (lx == 1) return cgetg(1,t_MAT);
     996      306481 :   dx = lgcols(x); y = cgetg(dx,t_MAT);
     997      306490 :   for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
     998      306499 :   return y;
     999             : }
    1000             : 
    1001             : /********************************************************************/
    1002             : /**                                                                **/
    1003             : /**                           SCALAR MATRICES                      **/
    1004             : /**                                                                **/
    1005             : /********************************************************************/
    1006             : 
    1007             : GEN
    1008        3762 : gen_matid(long n, void *E, const struct bb_field *S)
    1009             : {
    1010        3762 :   GEN y = cgetg(n+1,t_MAT), _0, _1;
    1011             :   long i;
    1012        3762 :   if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
    1013        3762 :   _0 = S->s(E,0);
    1014        3762 :   _1 = S->s(E,1);
    1015       16760 :   for (i=1; i<=n; i++)
    1016             :   {
    1017       12998 :     GEN z = const_col(n, _0); gel(z,i) = _1;
    1018       12998 :     gel(y, i) = z;
    1019             :   }
    1020        3762 :   return y;
    1021             : }
    1022             : 
    1023             : GEN
    1024          35 : matid_F2xqM(long n, GEN T)
    1025             : {
    1026             :   void *E;
    1027          35 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1028          35 :   return gen_matid(n, E, S);
    1029             : }
    1030             : GEN
    1031          56 : matid_FlxqM(long n, GEN T, ulong p)
    1032             : {
    1033             :   void *E;
    1034          56 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1035          56 :   return gen_matid(n, E, S);
    1036             : }
    1037             : 
    1038             : GEN
    1039     1007678 : matid_Flm(long n)
    1040             : {
    1041     1007678 :   GEN y = cgetg(n+1,t_MAT);
    1042             :   long i;
    1043     1007677 :   if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
    1044     1007668 :   for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
    1045     1007685 :   return y;
    1046             : }
    1047             : 
    1048             : GEN
    1049          42 : scalar_Flm(long s, long n)
    1050             : {
    1051             :   long i;
    1052          42 :   GEN y = cgetg(n+1,t_MAT);
    1053          42 :   for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
    1054          42 :   return y;
    1055             : }
    1056             : 
    1057             : /********************************************************************/
    1058             : /**                                                                **/
    1059             : /**                           CONVERSIONS                          **/
    1060             : /**                                                                **/
    1061             : /********************************************************************/
    1062             : GEN
    1063    23045228 : ZV_to_Flv(GEN x, ulong p)
    1064    23045228 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
    1065             : 
    1066             : GEN
    1067     4254118 : ZM_to_Flm(GEN x, ulong p)
    1068     4254118 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
    1069             : 
    1070             : GEN
    1071        2317 : ZMV_to_FlmV(GEN x, ulong m)
    1072        2317 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
    1073             : 
    1074             : /*                          TO INTMOD                        */
    1075             : static GEN
    1076    12069312 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
    1077             : static GEN
    1078       19376 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
    1079             : 
    1080             : GEN
    1081      302127 : Fp_to_mod(GEN z, GEN p)
    1082             : {
    1083      302127 :   retmkintmod(modii(z, p), icopy(p));
    1084             : }
    1085             : 
    1086             : /* z in Z[X], return z * Mod(1,p), normalized*/
    1087             : GEN
    1088     1864035 : FpX_to_mod(GEN z, GEN p)
    1089             : {
    1090     1864035 :   long i,l = lg(z);
    1091             :   GEN x;
    1092     1864035 :   if (l == 2)
    1093             :   {
    1094      525637 :     x = cgetg(3,t_POL); x[1] = z[1];
    1095      525637 :     gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
    1096             :   }
    1097     1338398 :   x = cgetg(l,t_POL);
    1098     1338398 :   if (l >2) p = icopy(p);
    1099     1338398 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1100     1338398 :   x[1] = z[1]; return normalizepol_lg(x,l);
    1101             : }
    1102             : 
    1103             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1104             : GEN
    1105        7196 : FpV_to_mod(GEN z, GEN p)
    1106             : {
    1107        7196 :   long i,l = lg(z);
    1108        7196 :   GEN x = cgetg(l, t_VEC);
    1109        7196 :   if (l == 1) return x;
    1110        7196 :   p = icopy(p);
    1111        7196 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1112        7196 :   return x;
    1113             : }
    1114             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1115             : GEN
    1116          77 : FpC_to_mod(GEN z, GEN p)
    1117             : {
    1118          77 :   long i, l = lg(z);
    1119          77 :   GEN x = cgetg(l, t_COL);
    1120          77 :   if (l == 1) return x;
    1121          63 :   p = icopy(p);
    1122          63 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1123          63 :   return x;
    1124             : }
    1125             : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
    1126             : GEN
    1127       56392 : FpM_to_mod(GEN z, GEN p)
    1128             : {
    1129       56392 :   long i, j, m, l = lg(z);
    1130       56392 :   GEN  x = cgetg(l,t_MAT), y, zi;
    1131       56392 :   if (l == 1) return x;
    1132       56371 :   m = lgcols(z);
    1133       56371 :   p = icopy(p);
    1134      181986 :   for (i=1; i<l; i++)
    1135             :   {
    1136      125615 :     gel(x,i) = cgetg(m,t_COL);
    1137      125615 :     y = gel(x,i); zi= gel(z,i);
    1138      125615 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1139             :   }
    1140       56371 :   return x;
    1141             : }
    1142             : GEN
    1143          28 : Flc_to_mod(GEN z, ulong pp)
    1144             : {
    1145          28 :   long i, l = lg(z);
    1146          28 :   GEN p, x = cgetg(l, t_COL);
    1147          28 :   if (l == 1) return x;
    1148          28 :   p = utoipos(pp);
    1149          28 :   for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
    1150          28 :   return x;
    1151             : }
    1152             : GEN
    1153         189 : Flm_to_mod(GEN z, ulong pp)
    1154             : {
    1155         189 :   long i, j, m, l = lg(z);
    1156         189 :   GEN p, x = cgetg(l,t_MAT), y, zi;
    1157         189 :   if (l == 1) return x;
    1158         161 :   m = lgcols(z);
    1159         161 :   p = utoipos(pp);
    1160        1414 :   for (i=1; i<l; i++)
    1161             :   {
    1162        1253 :     gel(x,i) = cgetg(m,t_COL);
    1163        1253 :     y = gel(x,i); zi= gel(z,i);
    1164        1253 :     for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
    1165             :   }
    1166         161 :   return x;
    1167             : }
    1168             : 
    1169             : GEN
    1170         574 : FpVV_to_mod(GEN z, GEN p)
    1171             : {
    1172         574 :   long i, j, m, l = lg(z);
    1173         574 :   GEN  x = cgetg(l,t_VEC), y, zi;
    1174         574 :   if (l == 1) return x;
    1175         574 :   m = lgcols(z);
    1176         574 :   p = icopy(p);
    1177        1246 :   for (i=1; i<l; i++)
    1178             :   {
    1179         672 :     gel(x,i) = cgetg(m,t_VEC);
    1180         672 :     y = gel(x,i); zi= gel(z,i);
    1181         672 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1182             :   }
    1183         574 :   return x;
    1184             : }
    1185             : 
    1186             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1187             : GEN
    1188           7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
    1189             : {
    1190           7 :   long i,l = lg(z);
    1191           7 :   GEN x = cgetg(l, t_COL);
    1192           7 :   if (l == 1) return x;
    1193           7 :   T = FpX_to_mod(T, p);
    1194          21 :   for (i=1; i<l; i++)
    1195          14 :     gel(x,i) = mkpolmod(FpX_to_mod(gel(z,i), p), T);
    1196           7 :   return x;
    1197             : }
    1198             : 
    1199             : static GEN
    1200      584437 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
    1201             : {
    1202      584437 :   GEN z = typ(x)==t_INT ? Fp_to_mod(x, p): FpX_to_mod(x, p);
    1203      584437 :   return mkpolmod(z, T);
    1204             : }
    1205             : 
    1206             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1207             : static GEN
    1208      108115 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
    1209      108115 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
    1210             : 
    1211             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1212             : GEN
    1213       26159 : FqM_to_mod(GEN z, GEN T, GEN p)
    1214             : {
    1215             :   GEN x;
    1216       26159 :   long i,l = lg(z);
    1217       26159 :   if (!T) return FpM_to_mod(z, p);
    1218       26159 :   x = cgetg(l, t_MAT);
    1219       26159 :   if (l == 1) return x;
    1220       26131 :   T = FpX_to_mod(T, p);
    1221      134246 :   for (i=1; i<l; i++)
    1222      108115 :     gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
    1223       26131 :   return x;
    1224             : }
    1225             : 
    1226             : /********************************************************************/
    1227             : /*                                                                  */
    1228             : /*                     Blackbox linear algebra                      */
    1229             : /*                                                                  */
    1230             : /********************************************************************/
    1231             : 
    1232             : /* A sparse column (zCs) v is a t_COL with two components C and E which are
    1233             :  * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
    1234             :  * (e_j) is the canonical basis.
    1235             :  * A sparse matrix (zMs) is a t_VEC of zCs */
    1236             : 
    1237             : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
    1238             :  * integer representing an element of Fp. This is important since p can be
    1239             :  * large and p+E[i] would not fit in a C long.  */
    1240             : 
    1241             : /* RgCs and RgMs are similar, except that the type of the component is
    1242             :  * unspecified. Functions handling RgCs/RgMs must be independent of the type
    1243             :  * of E. */
    1244             : 
    1245             : /* Most functions take an argument nbrow which is the number of lines of the
    1246             :  * column/matrix, which cannot be derived from the data. */
    1247             : 
    1248             : GEN
    1249           0 : zCs_to_ZC(GEN R, long nbrow)
    1250             : {
    1251           0 :   GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
    1252           0 :   long j, l = lg(C);
    1253           0 :   for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
    1254           0 :   return c;
    1255             : }
    1256             : 
    1257             : GEN
    1258           0 : zMs_to_ZM(GEN x, long nbrow)
    1259           0 : { pari_APPLY_same(zCs_to_ZC(gel(x, i), nbrow)) }
    1260             : 
    1261             : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
    1262             :  * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
    1263             : GEN
    1264         126 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
    1265             : {
    1266         126 :   pari_sp ltop = avma;
    1267         126 :   long col = 0, n = lg(B)-1, m = 2*n+1;
    1268         126 :   if (ZV_equal0(B)) return zerocol(n);
    1269         252 :   while (++col <= n)
    1270             :   {
    1271         126 :     pari_sp btop = avma, av;
    1272             :     long i, lQ;
    1273         126 :     GEN V, Q, M, W = B;
    1274         126 :     GEN b = cgetg(m+2, t_POL);
    1275         126 :     b[1] = evalsigne(1)|evalvarn(0);
    1276         126 :     gel(b, 2) = gel(W, col);
    1277       51524 :     for (i = 3; i<m+2; i++)
    1278       51398 :       gel(b, i) = cgeti(lgefint(p));
    1279         126 :     av = avma;
    1280       51524 :     for (i = 3; i<m+2; i++)
    1281             :     {
    1282       51398 :       W = f(E, W);
    1283       51398 :       affii(gel(W, col),gel(b, i));
    1284       51398 :       if (gc_needed(av,1))
    1285             :       {
    1286        2000 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
    1287        2000 :         W = gerepileupto(av, W);
    1288             :       }
    1289             :     }
    1290         126 :     b = FpX_renormalize(b, m+2);
    1291         126 :     if (lgpol(b)==0) {set_avma(btop); continue; }
    1292         126 :     M = FpX_halfgcd(b, pol_xn(m, 0), p);
    1293         126 :     Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
    1294         126 :     W = B; lQ =lg(Q);
    1295         126 :     if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
    1296         126 :     V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
    1297         126 :     av = avma;
    1298       25279 :     for (i = lQ-3; i > 1; i--)
    1299             :     {
    1300       25153 :       W = f(E, W);
    1301       25153 :       V = ZC_lincomb(gen_1, gel(Q,i), V, W);
    1302       25153 :       if (gc_needed(av,1))
    1303             :       {
    1304        1583 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
    1305        1583 :         gerepileall(av, 2, &V, &W);
    1306             :       }
    1307             :     }
    1308         126 :     V = FpC_red(V, p);
    1309         126 :     W = FpC_sub(f(E,V), B, p);
    1310         252 :     if (ZV_equal0(W)) return gerepilecopy(ltop, V);
    1311           0 :     av = avma;
    1312           0 :     for (i = 1; i <= n; ++i)
    1313             :     {
    1314           0 :       V = W;
    1315           0 :       W = f(E, V);
    1316           0 :       if (ZV_equal0(W))
    1317           0 :         return gerepilecopy(ltop, shallowtrans(V));
    1318           0 :       gerepileall(av, 2, &V, &W);
    1319             :     }
    1320           0 :     set_avma(btop);
    1321             :   }
    1322           0 :   return NULL;
    1323             : }
    1324             : 
    1325             : GEN
    1326           0 : zMs_ZC_mul(GEN M, GEN B)
    1327             : {
    1328             :   long i, j;
    1329           0 :   long n = lg(B)-1;
    1330           0 :   GEN V = zerocol(n);
    1331           0 :   for (i = 1; i <= n; ++i)
    1332           0 :     if (signe(gel(B, i)))
    1333             :     {
    1334           0 :       GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1335           0 :       long l = lg(C);
    1336           0 :       for (j = 1; j < l; ++j)
    1337             :       {
    1338           0 :         long k = C[j];
    1339           0 :         switch(E[j])
    1340             :         {
    1341             :         case 1:
    1342           0 :           gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
    1343           0 :           break;
    1344             :         case -1:
    1345           0 :           gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
    1346           0 :           break;
    1347             :         default:
    1348           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]));
    1349           0 :           break;
    1350             :         }
    1351             :       }
    1352             :     }
    1353           0 :   return V;
    1354             : }
    1355             : 
    1356             : GEN
    1357           0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
    1358             : 
    1359             : GEN
    1360       76803 : ZV_zMs_mul(GEN B, GEN M)
    1361             : {
    1362             :   long i, j;
    1363       76803 :   long m = lg(M)-1;
    1364       76803 :   GEN V = cgetg(m+1,t_VEC);
    1365    39361216 :   for (i = 1; i <= m; ++i)
    1366             :   {
    1367    39284413 :     GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1368    39284413 :     long l = lg(C);
    1369    39284413 :     GEN z = mulis(gel(B, C[1]), E[1]);
    1370   330271731 :     for (j = 2; j < l; ++j)
    1371             :     {
    1372   290987318 :       long k = C[j];
    1373   290987318 :       switch(E[j])
    1374             :       {
    1375             :       case 1:
    1376   196103761 :         z = addii(z, gel(B,k));
    1377   196103761 :         break;
    1378             :       case -1:
    1379    79342341 :         z = subii(z, gel(B,k));
    1380    79342341 :         break;
    1381             :       default:
    1382    15541216 :         z = addii(z, mulis(gel(B,k), E[j]));
    1383    15541216 :         break;
    1384             :       }
    1385             :     }
    1386    39284413 :     gel(V,i) = z;
    1387             :   }
    1388       76803 :   return V;
    1389             : }
    1390             : 
    1391             : GEN
    1392         126 : FpV_FpMs_mul(GEN B, GEN M, GEN p) { return FpV_red(ZV_zMs_mul(B, M), p); }
    1393             : 
    1394             : GEN
    1395           0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
    1396             : {
    1397           0 :   pari_sp av = avma, av2;
    1398           0 :   GEN xi, xb, pi = gen_1, P;
    1399             :   long i;
    1400           0 :   if (!C) {
    1401           0 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    1402           0 :     if (!C) pari_err_INV("ZlM_gauss", a);
    1403             :   }
    1404           0 :   P = utoipos(p);
    1405           0 :   av2 = avma;
    1406           0 :   xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1407           0 :   xb = Flm_to_ZM(xi);
    1408           0 :   for (i = 2; i <= e; i++)
    1409             :   {
    1410           0 :     pi = muliu(pi, p); /* = p^(i-1) */
    1411           0 :     b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
    1412           0 :     if (gc_needed(av,2))
    1413             :     {
    1414           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
    1415           0 :       gerepileall(av2,3, &pi,&b,&xb);
    1416             :     }
    1417           0 :     xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1418           0 :     xb = ZM_add(xb, nm_Z_mul(xi, pi));
    1419             :   }
    1420           0 :   return gerepileupto(av, xb);
    1421             : }
    1422             : 
    1423             : struct wrapper_modp_s {
    1424             :   GEN (*f)(void*, GEN);
    1425             :   void *E;
    1426             :   GEN p;
    1427             : };
    1428             : 
    1429             : /* compute f(x) mod p */
    1430             : static GEN
    1431       76677 : wrap_relcomb_modp(void *E, GEN x)
    1432             : {
    1433       76677 :   struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
    1434       76677 :   return FpC_red(W->f(W->E, x), W->p);
    1435             : }
    1436             : static GEN
    1437           0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
    1438             : 
    1439             : static GEN
    1440       76677 : wrap_relker(void*E, GEN x) { return ZV_zMs_mul(x, (GEN)E); }
    1441             : 
    1442             : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
    1443             : GEN
    1444           0 : gen_ZpM_Dixon(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
    1445             : {
    1446             :   struct wrapper_modp_s W;
    1447           0 :   pari_sp av = avma;
    1448           0 :   GEN xb, xi, pi = gen_1;
    1449             :   long i;
    1450           0 :   W.E = E;
    1451           0 :   W.f = f;
    1452           0 :   W.p = p;
    1453           0 :   xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
    1454           0 :   if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
    1455           0 :   xb = xi;
    1456           0 :   for (i = 2; i <= e; i++)
    1457             :   {
    1458           0 :     pi = mulii(pi, p); /* = p^(i-1) */
    1459           0 :     B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
    1460           0 :     if (gc_needed(av,2))
    1461             :     {
    1462           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon. i=%ld",i);
    1463           0 :       gerepileall(av,3, &pi,&B,&xb);
    1464             :     }
    1465           0 :     xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
    1466           0 :     if (!xi) return NULL;
    1467           0 :     if (typ(xi) == t_VEC) return gerepileupto(av, xi);
    1468           0 :     xb = ZC_add(xb, ZC_Z_mul(xi, pi));
    1469             :   }
    1470           0 :   return gerepileupto(av, xb);
    1471             : }
    1472             : 
    1473             : static GEN
    1474       25699 : vecprow(GEN A, GEN prow)
    1475             : {
    1476       25699 :   return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
    1477             : }
    1478             : 
    1479             : /* Solve the equation MX = A. Return either a solution as a t_COL,
    1480             :  * or the index of a column which is linearly dependent from the others as a
    1481             :  * t_VECSMALL with a single component. */
    1482             : GEN
    1483           0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
    1484             : {
    1485           0 :   pari_sp av = avma;
    1486             :   GEN pcol, prow;
    1487           0 :   long nbi=lg(M)-1, lR;
    1488             :   long i, n;
    1489             :   GEN Mp, Ap, Rp;
    1490             :   pari_timer ti;
    1491           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1492           0 :   RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
    1493           0 :   if (!pcol) return gc_NULL(av);
    1494           0 :   if (DEBUGLEVEL)
    1495           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
    1496           0 :   n = lg(pcol)-1;
    1497           0 :   Mp = cgetg(n+1, t_MAT);
    1498           0 :   for(i=1; i<=n; i++)
    1499           0 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1500           0 :   Ap = zCs_to_ZC(vecprow(A, prow), n);
    1501           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1502           0 :   Rp = gen_ZpM_Dixon((void*)Mp,wrap_relcomb, Ap, p, e);
    1503           0 :   if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
    1504           0 :   if (!Rp) return gc_NULL(av);
    1505           0 :   lR = lg(Rp)-1;
    1506           0 :   if (typ(Rp) == t_COL)
    1507             :   {
    1508           0 :     GEN R = zerocol(nbi+1);
    1509           0 :     for (i=1; i<=lR; i++)
    1510           0 :        gel(R,pcol[i]) = gel(Rp,i);
    1511           0 :     return gerepilecopy(av, R);
    1512             :   }
    1513           0 :   for (i = 1; i <= lR; ++i)
    1514           0 :     if (signe(gel(Rp, i)))
    1515           0 :       return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
    1516             :   return NULL; /* LCOV_EXCL_LINE */
    1517             : }
    1518             : 
    1519             : GEN
    1520           0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
    1521             : {
    1522           0 :   return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1523             : }
    1524             : 
    1525             : GEN
    1526           0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
    1527             : {
    1528             :   GEN res;
    1529           0 :   pari_CATCH(e_INV)
    1530             :   {
    1531           0 :     GEN E = pari_err_last();
    1532           0 :     GEN x = err_get_compo(E,2);
    1533           0 :     if (typ(x) != t_INTMOD) pari_err(0,E);
    1534           0 :     if (DEBUGLEVEL)
    1535           0 :       pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
    1536           0 :     res = NULL;
    1537             :   } pari_TRY {
    1538           0 :     res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1539           0 :   } pari_ENDCATCH
    1540           0 :   return res;
    1541             : }
    1542             : 
    1543             : static GEN
    1544         126 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
    1545             : {
    1546         126 :   long i, j, oldf = 0, f = 0;
    1547         126 :   long lrow = lg(prow), lM = lg(M);
    1548         126 :   GEN W = const_vecsmall(lM-1,1);
    1549         126 :   GEN R = cgetg(lrow, t_VEC);
    1550      237272 :   for (i=1; i<lrow; i++)
    1551      237146 :     gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
    1552             :   do
    1553             :   {
    1554         414 :     oldf = f;
    1555      443945 :     for(i=1; i<lM; i++)
    1556      443531 :       if (W[i])
    1557             :       {
    1558      150721 :         GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
    1559      150721 :         long c=0, cj=0, lF = lg(F);
    1560     1317607 :         for(j=1; j<lF; j++)
    1561     1166886 :           if (!gel(R,F[j])) { c++; cj=j; }
    1562      150721 :         if (c>=2) continue;
    1563      129590 :         if (c==1)
    1564             :         {
    1565       41694 :           pari_sp av = avma;
    1566       41694 :           GEN s = gen_0;
    1567      362307 :           for(j=1; j<lF; j++)
    1568      320613 :             if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
    1569             :           /* s /= -E[cj] mod p */
    1570       41694 :           s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
    1571       41694 :           gel(R,F[cj]) = gerepileuptoint(av, s);
    1572       41694 :           f++;
    1573             :         }
    1574      129590 :         W[i]=0;
    1575             :       }
    1576         414 :   } while(oldf!=f);
    1577      237272 :   for (i=1; i<lrow; i++)
    1578      237146 :     if (!gel(R,i)) gel(R,i) = gen_0;
    1579         126 :   return R;
    1580             : }
    1581             : 
    1582             : /* Return a linear form Y such that YM is essentially 0 */
    1583             : GEN
    1584         126 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
    1585             : {
    1586         126 :   pari_sp av = avma, av2;
    1587             :   GEN pcol, prow;
    1588             :   long i, n;
    1589             :   GEN Mp, B, MB, R, Rp;
    1590             :   pari_timer ti;
    1591             :   struct wrapper_modp_s W;
    1592         126 :   if (DEBUGLEVEL) timer_start(&ti);
    1593         126 :   RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
    1594         126 :   if (!pcol) return gc_NULL(av);
    1595         126 :   if (DEBUGLEVEL)
    1596           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
    1597         126 :   n = lg(pcol)-1;
    1598         126 :   Mp = cgetg(n+1, t_MAT);
    1599       25825 :   for(i=1; i<=n; i++)
    1600       25699 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1601         126 :   W.E = (void*) Mp;
    1602         126 :   W.f = wrap_relker;
    1603         126 :   W.p = p;
    1604         126 :   av2 = avma;
    1605             :   for(;;)
    1606             :   {
    1607         126 :     set_avma(av2);
    1608         126 :     B = random_FpV(n, p);
    1609         126 :     MB = FpV_FpMs_mul(B, Mp, p);
    1610         126 :     if (DEBUGLEVEL) timer_start(&ti);
    1611         126 :     pari_CATCH(e_INV)
    1612             :     {
    1613           0 :       GEN E = pari_err_last();
    1614           0 :       GEN x = err_get_compo(E,2);
    1615           0 :       if (typ(x) != t_INTMOD) pari_err(0,E);
    1616           0 :       if (DEBUGLEVEL)
    1617           0 :         pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
    1618           0 :       Rp = NULL;
    1619             :     } pari_TRY {
    1620         126 :       Rp = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, MB, p);
    1621         126 :     } pari_ENDCATCH
    1622         126 :     if (!Rp) continue;
    1623         126 :     if (typ(Rp)==t_COL)
    1624             :     {
    1625         126 :       Rp = FpC_sub(Rp,B,p);
    1626         126 :       if (ZV_equal0(Rp)) continue;
    1627             :     }
    1628         126 :     R = FpMs_structelim_back(M, Rp, prow, p);
    1629         126 :     if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
    1630         126 :     return gerepilecopy(av, R);
    1631             :   }
    1632             : }
    1633             : 
    1634             : GEN
    1635          42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
    1636             : {
    1637          42 :   return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
    1638             : }

Generated by: LCOV version 1.13