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.16.2 lcov report (development 29115-f22e516b23) Lines: 741 964 76.9 %
Date: 2024-03-18 08:03:28 Functions: 109 135 80.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_mat
      19             : 
      20             : /********************************************************************/
      21             : /**                                                                **/
      22             : /**                           REDUCTION                            **/
      23             : /**                                                                **/
      24             : /********************************************************************/
      25             : /* z in Z^n, return lift(Col(z) * Mod(1,p)) */
      26             : GEN
      27    27294784 : FpC_red(GEN x, GEN p)
      28   225933296 : { pari_APPLY_type(t_COL, modii(gel(x,i), p)) }
      29             : 
      30             : /* z in Z^n, return lift(Vec(z) * Mod(1,p)) */
      31             : GEN
      32         539 : FpV_red(GEN x, GEN p)
      33       27730 : { pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
      34             : 
      35             : GEN
      36     6811971 : FpC_center(GEN x, GEN p, GEN pov2)
      37    43079862 : { pari_APPLY_type(t_COL, Fp_center(gel(x,i), p, pov2)) }
      38             : 
      39             : GEN
      40      164710 : Flv_center(GEN x, ulong p, ulong ps2)
      41     2621724 : { pari_APPLY_ulong(Fl_center(uel(x,i),p,ps2)) }
      42             : 
      43             : /* z in Mat m,n(Z), return lift(z * Mod(1,p)) */
      44             : GEN
      45     5046265 : FpM_red(GEN x, GEN p)
      46    21749492 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
      47             : 
      48             : GEN
      49     2740673 : FpM_center(GEN x, GEN p, GEN pov2)
      50     8237223 : { pari_APPLY_same(FpC_center(gel(x,i), p,pov2)) }
      51             : 
      52             : /* p != 3; assume entries in [0,p[ and ps2 = p>>1. */
      53             : static void
      54          77 : _FpC_center_inplace(GEN z, GEN p, GEN ps2)
      55             : {
      56          77 :   long i, l = lg(z);
      57         357 :   for (i = 1; i < l; i++)
      58             :   { /* HACK: assume p != 3, which ensures u = gen_[0-2] is never written to */
      59         280 :     GEN u = gel(z,i);
      60         280 :     if (abscmpii(u, ps2) > 0) subiiz(u, p, u);
      61             :   }
      62          77 : }
      63             : static void
      64           7 : _F3C_center_inplace(GEN z)
      65             : {
      66           7 :   long i, l = lg(z);
      67          35 :   for (i = 1; i < l; i++) /* z[i] = 0, 1 : no-op */
      68          28 :     if (equaliu(gel(z,i), 2)) gel(z,i) = gen_m1;
      69           7 : }
      70             : void
      71           0 : FpC_center_inplace(GEN z, GEN p, GEN ps2)
      72             : {
      73           0 :   if (equaliu(p, 3))
      74           0 :     _F3C_center_inplace(z);
      75             :   else
      76           0 :     _FpC_center_inplace(z, p, ps2);
      77           0 : }
      78             : 
      79             : void
      80          42 : FpM_center_inplace(GEN z, GEN p, GEN pov2)
      81             : {
      82          42 :   long i, l = lg(z);
      83          42 :   if (equaliu(p, 3))
      84          14 :     for (i = 1; i < l; i++) _F3C_center_inplace(gel(z,i));
      85             :   else
      86         112 :     for (i = 1; i < l; i++) _FpC_center_inplace(gel(z,i), p, pov2);
      87          42 : }
      88             : GEN
      89       11347 : Flm_center(GEN x, ulong p, ulong ps2)
      90      175875 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
      91             : 
      92             : GEN
      93         119 : random_FpV(long d, GEN p)
      94             : {
      95             :   long i;
      96         119 :   GEN y = cgetg(d+1,t_VEC);
      97       25035 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
      98         119 :   return y;
      99             : }
     100             : 
     101             : GEN
     102         413 : random_FpC(long d, GEN p)
     103             : {
     104             :   long i;
     105         413 :   GEN y = cgetg(d+1,t_COL);
     106        3409 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
     107         413 :   return y;
     108             : }
     109             : 
     110             : GEN
     111       41262 : random_Flv(long n, ulong p)
     112             : {
     113       41262 :   GEN y = cgetg(n+1, t_VECSMALL);
     114             :   long i;
     115      261123 :   for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
     116       41262 :   return y;
     117             : }
     118             : 
     119             : /********************************************************************/
     120             : /**                                                                **/
     121             : /**                           ADD, SUB                             **/
     122             : /**                                                                **/
     123             : /********************************************************************/
     124             : GEN
     125      242687 : FpC_add(GEN x, GEN y, GEN p)
     126     3154231 : { pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
     127             : 
     128             : GEN
     129           0 : FpV_add(GEN x, GEN y, GEN p)
     130           0 : { pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
     131             : 
     132             : GEN
     133           0 : FpM_add(GEN x, GEN y, GEN p)
     134           0 : { pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
     135             : 
     136             : GEN
     137      773173 : Flv_add(GEN x, GEN y, ulong p)
     138             : {
     139      773173 :   if (p==2)
     140     1349984 :     pari_APPLY_ulong(x[i]^y[i])
     141             :   else
     142     8233261 :     pari_APPLY_ulong(Fl_add(x[i], y[i], p))
     143             : }
     144             : 
     145             : void
     146      393758 : Flv_add_inplace(GEN x, GEN y, ulong p)
     147             : {
     148      393758 :   long i, l = lg(x);
     149      393758 :   if (p==2)
     150     1235569 :     for (i = 1; i < l; i++) x[i] ^= y[i];
     151             :   else
     152      164878 :     for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
     153      393758 : }
     154             : 
     155             : ulong
     156        3934 : Flv_sum(GEN x, ulong p)
     157             : {
     158        3934 :   long i, l = lg(x);
     159        3934 :   ulong s = 0;
     160        3934 :   if (p==2)
     161       17738 :     for (i = 1; i < l; i++) s ^= x[i];
     162             :   else
     163           0 :     for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
     164        3934 :   return s;
     165             : }
     166             : 
     167             : GEN
     168      493826 : FpC_sub(GEN x, GEN y, GEN p)
     169     6679131 : { pari_APPLY_type(t_COL, Fp_sub(gel(x,i), gel(y,i), p)) }
     170             : 
     171             : GEN
     172           0 : FpV_sub(GEN x, GEN y, GEN p)
     173           0 : { pari_APPLY_type(t_VEC, Fp_sub(gel(x,i), gel(y,i), p)) }
     174             : 
     175             : GEN
     176           0 : FpM_sub(GEN x, GEN y, GEN p)
     177           0 : { pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
     178             : 
     179             : GEN
     180   211151659 : Flv_sub(GEN x, GEN y, ulong p)
     181   958637399 : { pari_APPLY_ulong(Fl_sub(uel(x,i), uel(y,i), p)) }
     182             : 
     183             : void
     184           0 : Flv_sub_inplace(GEN x, GEN y, ulong p)
     185             : {
     186           0 :   long i, l = lg(x);
     187           0 :   for (i = 1; i < l; i++) x[i] = Fl_sub(x[i], y[i], p);
     188           0 : }
     189             : 
     190             : GEN
     191       49654 : Flm_Fl_add(GEN x, ulong y, ulong p)
     192             : {
     193       49654 :   long l = lg(x), i, j;
     194       49654 :   GEN z = cgetg(l,t_MAT);
     195             : 
     196       49654 :   if (l==1) return z;
     197       49654 :   if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
     198      236093 :   for (i=1; i<l; i++)
     199             :   {
     200      186436 :     GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
     201      186439 :     gel(z,i) = zi;
     202     1306820 :     for (j=1; j<l; j++) zi[j] = xi[j];
     203      186439 :     zi[i] = Fl_add(zi[i], y, p);
     204             :   }
     205       49657 :   return z;
     206             : }
     207             : GEN
     208        3689 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
     209             : 
     210             : GEN
     211        9224 : Flm_add(GEN x, GEN y, ulong p)
     212      158026 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
     213             : 
     214             : GEN
     215    25032282 : Flm_sub(GEN x, GEN y, ulong p)
     216   236169860 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
     217             : 
     218             : /********************************************************************/
     219             : /**                                                                **/
     220             : /**                           MULTIPLICATION                       **/
     221             : /**                                                                **/
     222             : /********************************************************************/
     223             : GEN
     224     1117353 : FpC_Fp_mul(GEN x, GEN y, GEN p)
     225    12180233 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
     226             : 
     227             : GEN
     228      855692 : Flv_Fl_mul(GEN x, ulong y, ulong p)
     229    11181405 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
     230             : 
     231             : GEN
     232         907 : Flv_Fl_div(GEN x, ulong y, ulong p)
     233             : {
     234         907 :   return Flv_Fl_mul(x, Fl_inv(y, p), p);
     235             : }
     236             : void
     237           0 : Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
     238             : {
     239           0 :   Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
     240           0 : }
     241             : GEN
     242      755076 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
     243      755076 :   long i, j, h, l = lg(X);
     244      755076 :   GEN A = cgetg(l, t_MAT);
     245      755076 :   if (l == 1) return A;
     246      755076 :   h = lgcols(X);
     247     3960889 :   for (j=1; j<l; j++)
     248             :   {
     249     3205926 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     250    23160799 :     for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
     251     3205811 :     gel(A,j) = a;
     252             :   }
     253      754963 :   return A;
     254             : }
     255             : 
     256             : /* x *= y */
     257             : void
     258     7730458 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
     259             : {
     260             :   long i;
     261     7730458 :   if (HIGHWORD(y | p))
     262     8564682 :     for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
     263             :   else
     264    31975374 :     for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
     265     7730457 : }
     266             : void
     267           0 : Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
     268           0 : { Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
     269             : 
     270             : /* set y *= x */
     271             : void
     272           0 : Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
     273             : {
     274           0 :   long i, j, m, l = lg(y);
     275           0 :   if (l == 1) return;
     276           0 :   m = lgcols(y);
     277           0 :   if (HIGHWORD(x | p))
     278           0 :     for(j=1; j<l; j++)
     279           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
     280             :   else
     281           0 :     for(j=1; j<l; j++)
     282           0 :       for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
     283             : }
     284             : 
     285             : /* return x * y */
     286             : static GEN
     287    16913374 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
     288             : {
     289    16913374 :   long i, j, m, l = lg(y);
     290    16913374 :   GEN z = cgetg(l, t_MAT);
     291    16911081 :   if (l == 1) return z;
     292    16911081 :   m = lgcols(y);
     293   150689207 :   for(j=1; j<l; j++) {
     294   133641773 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     295   358592269 :     for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
     296             :   }
     297    17047434 :   return z;
     298             : }
     299             : 
     300             : /* return x * y */
     301             : static GEN
     302     8214402 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
     303             : {
     304     8214402 :   long i, j, m, l = lg(y);
     305     8214402 :   GEN z = cgetg(l, t_MAT);
     306     8214397 :   if (l == 1) return z;
     307     8214397 :   m = lgcols(y);
     308    44408035 :   for(j=1; j<l; j++) {
     309    36193659 :     GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     310   116320780 :     for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
     311             :   }
     312     8214376 :   return z;
     313             : }
     314             : 
     315             : /* return x * y */
     316             : GEN
     317    25091844 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
     318             : {
     319    25091844 :   if (HIGHWORD(p))
     320    16913506 :     return Flm_Fl_mul_pre_i(y, x, p, pi);
     321             :   else
     322     8178338 :     return Flm_Fl_mul_OK(y, x, p);
     323             : }
     324             : 
     325             : /* return x * y */
     326             : GEN
     327       36107 : Flm_Fl_mul(GEN y, ulong x, ulong p)
     328             : {
     329       36107 :   if (HIGHWORD(p))
     330          74 :     return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
     331             :   else
     332       36033 :     return Flm_Fl_mul_OK(y, x, p);
     333             : }
     334             : 
     335             : GEN
     336     3983393 : Flv_neg(GEN x, ulong p)
     337    30320974 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
     338             : 
     339             : void
     340        6187 : Flv_neg_inplace(GEN v, ulong p)
     341             : {
     342             :   long i;
     343      187303 :   for (i = 1; i < lg(v); ++i)
     344      181116 :     v[i] = Fl_neg(v[i], p);
     345        6187 : }
     346             : 
     347             : GEN
     348      314137 : Flm_neg(GEN x, ulong p)
     349     4266949 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
     350             : 
     351             : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
     352             : static GEN
     353    18706540 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
     354             : {
     355    18706540 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1));
     356             :   long k;
     357   220803094 :   for (k = 2; k < lx; k++)
     358             :   {
     359   202098075 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     360   202095014 :     if (signe(t)) c = addii(c, t);
     361             :   }
     362    18705019 :   return c;
     363             : }
     364             : 
     365             : static long
     366    97688656 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
     367             : {
     368             :   long k;
     369    97688656 :   long c = coeff(x,i,1) * y[1];
     370  1498282053 :   for (k = 2; k < lx; k++)
     371  1400593397 :     c += coeff(x,i,k) * y[k];
     372    97688656 :   return c;
     373             : }
     374             : 
     375             : GEN
     376     6450502 : zm_zc_mul(GEN x, GEN y)
     377             : {
     378     6450502 :   long lx = lg(x), l, i;
     379             :   GEN z;
     380     6450502 :   if (lx == 1) return cgetg(1, t_VECSMALL);
     381     6450502 :   l = lg(gel(x,1));
     382     6450502 :   z = cgetg(l,t_VECSMALL);
     383   104139158 :   for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
     384     6450502 :   return z;
     385             : }
     386             : 
     387             : GEN
     388        2114 : zm_mul(GEN x, GEN y)
     389             : {
     390        2114 :   long i,j,lx=lg(x), ly=lg(y);
     391             :   GEN z;
     392        2114 :   if (ly==1) return cgetg(1,t_MAT);
     393        2114 :   z = cgetg(ly,t_MAT);
     394        2114 :   if (lx==1)
     395             :   {
     396           0 :     for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
     397           0 :     return z;
     398             :   }
     399       30849 :   for (j=1; j<ly; j++)
     400       28735 :     gel(z,j) = zm_zc_mul(x, gel(y,j));
     401        2114 :   return z;
     402             : }
     403             : 
     404             : static ulong
     405   634139048 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
     406             : {
     407   634139048 :   ulong c = ucoeff(x,i,1) * uel(y,1);
     408             :   long k;
     409  7669222536 :   for (k = 2; k < lx; k++) {
     410  7035083488 :     c += ucoeff(x,i,k) * uel(y,k);
     411  7035083488 :     if (c & HIGHBIT) c %= p;
     412             :   }
     413   634139048 :   return c % p;
     414             : }
     415             : 
     416             : static ulong
     417   537899014 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
     418             : {
     419             :   ulong l0, l1, v1, h0, h1;
     420   537899014 :   long k = 1;
     421             :   LOCAL_OVERFLOW;
     422             :   LOCAL_HIREMAINDER;
     423   537899014 :   l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
     424  9910548468 :   while (++k < lx) {
     425  9372649454 :     l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
     426  9372649454 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     427             :   }
     428   537899014 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     429    61821132 :   else return remlll_pre(v1, h1, l1, p, pi);
     430             : }
     431             : 
     432             : static GEN
     433      179842 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
     434             : {
     435             :   long i,j;
     436      179842 :   GEN z = NULL;
     437             : 
     438     1610964 :   for (j=1; j<lx; j++)
     439             :   {
     440     1431122 :     if (!y[j]) continue;
     441      443882 :     if (!z) z = Flv_copy(gel(x,j));
     442    12370688 :     else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
     443             :   }
     444      179842 :   if (!z) z = zero_zv(l-1);
     445      179842 :   return z;
     446             : }
     447             : 
     448             : static GEN
     449     1659221 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
     450             : {
     451     1659221 :   GEN z = cgetg(l,t_COL);
     452             :   long i;
     453    16028594 :   for (i = 1; i < l; i++)
     454             :   {
     455    14369412 :     pari_sp av = avma;
     456    14369412 :     GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
     457    14369206 :     gel(z,i) = gerepileuptoint(av, modii(c,p));
     458             :   }
     459     1659182 :   return z;
     460             : }
     461             : 
     462             : static void
     463    73415178 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
     464             : {
     465             :   long i;
     466   707554204 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
     467    73428248 : }
     468             : static GEN
     469    70105975 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
     470             : {
     471    70105975 :   GEN z = cgetg(l,t_VECSMALL);
     472    70104893 :   __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     473    70106330 :   return z;
     474             : }
     475             : 
     476             : static void
     477    93121282 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     478             : {
     479             :   long i;
     480   632763953 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
     481    93396605 : }
     482             : static GEN
     483    89006123 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     484             : {
     485    89006123 :   GEN z = cgetg(l,t_VECSMALL);
     486    88843561 :   __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     487    89059694 :   return z;
     488             : }
     489             : 
     490             : GEN
     491     1353927 : FpM_mul(GEN x, GEN y, GEN p)
     492             : {
     493     1353927 :   long lx=lg(x), ly=lg(y);
     494             :   GEN z;
     495     1353927 :   pari_sp av = avma;
     496     1353927 :   if (ly==1) return cgetg(1,t_MAT);
     497     1353927 :   if (lx==1) return zeromat(0, ly-1);
     498     1353927 :   if (lgefint(p) == 3)
     499             :   {
     500     1308251 :     ulong pp = uel(p,2);
     501     1308251 :     if (pp == 2)
     502             :     {
     503      674672 :       x = ZM_to_F2m(x);
     504      674736 :       y = ZM_to_F2m(y);
     505      674742 :       z = F2m_to_ZM(F2m_mul(x,y));
     506             :     }
     507             :     else
     508             :     {
     509      633579 :       x = ZM_to_Flm(x, pp);
     510      633588 :       y = ZM_to_Flm(y, pp);
     511      633589 :       z = Flm_to_ZM(Flm_mul(x,y, pp));
     512             :     }
     513             :   } else
     514       45676 :     z = FpM_red(ZM_mul(x, y), p);
     515     1353997 :   return gerepileupto(av, z);
     516             : }
     517             : 
     518             : static GEN
     519    27467459 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     520             : {
     521             :   long j;
     522    27467459 :   GEN z = cgetg(ly, t_MAT);
     523    27466492 :   if (SMALL_ULONG(p))
     524    89423961 :     for (j=1; j<ly; j++)
     525    69823511 :       gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
     526             :   else
     527    96670800 :     for (j=1; j<ly; j++)
     528    88803630 :       gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
     529    27467620 :   return z;
     530             : }
     531             : 
     532             : /*
     533             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     534             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     535             : */
     536             : static void
     537       66664 : add_slices_ip(long m, long n,
     538             :            GEN A, long ma, long da, long na, long ea,
     539             :            GEN B, long mb, long db, long nb, long eb,
     540             :            GEN M, long dy, long dx, ulong p)
     541             : {
     542       66664 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     543             :   GEN C;
     544             : 
     545     2807189 :   for (j = 1; j <= min_e; j++) {
     546     2740525 :     C = gel(M, j + dx) + dy;
     547   125600386 :     for (i = 1; i <= min_d; i++)
     548   122859861 :       uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
     549   122859876 :                         ucoeff(B, mb + i, nb + j), p);
     550     2808130 :     for (; i <= da; i++)
     551       67620 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     552     2740510 :     for (; i <= db; i++)
     553           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     554     2740510 :     for (; i <= m; i++)
     555           0 :       uel(C, i) = 0;
     556             :   }
     557       70367 :   for (; j <= ea; j++) {
     558        3703 :     C = gel(M, j + dx) + dy;
     559      201369 :     for (i = 1; i <= da; i++)
     560      197666 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     561        3703 :     for (; i <= m; i++)
     562           0 :       uel(C, i) = 0;
     563             :   }
     564       66664 :   for (; j <= eb; j++) {
     565           0 :     C = gel(M, j + dx) + dy;
     566           0 :     for (i = 1; i <= db; i++)
     567           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     568           0 :     for (; i <= m; i++)
     569           0 :       uel(C, i) = 0;
     570             :   }
     571       66664 :   for (; j <= n; j++) {
     572           0 :     C = gel(M, j + dx) + dy;
     573           0 :     for (i = 1; i <= m; i++)
     574           0 :       uel(C, i) = 0;
     575             :   }
     576       66664 : }
     577             : 
     578             : static GEN
     579       33332 : add_slices(long m, long n,
     580             :            GEN A, long ma, long da, long na, long ea,
     581             :            GEN B, long mb, long db, long nb, long eb, ulong p)
     582             : {
     583             :   GEN M;
     584             :   long j;
     585       33332 :   M = cgetg(n + 1, t_MAT);
     586     1428179 :   for (j = 1; j <= n; j++)
     587     1394850 :     gel(M, j) = cgetg(m + 1, t_VECSMALL);
     588       33329 :   add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
     589       33332 :   return M;
     590             : }
     591             : 
     592             : /*
     593             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     594             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     595             : */
     596             : static GEN
     597       58331 : subtract_slices(long m, long n,
     598             :                 GEN A, long ma, long da, long na, long ea,
     599             :                 GEN B, long mb, long db, long nb, long eb, ulong p)
     600             : {
     601       58331 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     602       58331 :   GEN M = cgetg(n + 1, t_MAT), C;
     603             : 
     604     2549527 :   for (j = 1; j <= min_e; j++) {
     605     2491197 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     606   117619925 :     for (i = 1; i <= min_d; i++)
     607   115128747 :       uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
     608   115128567 :                         ucoeff(B, mb + i, nb + j), p);
     609     2594866 :     for (; i <= da; i++)
     610      103508 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     611     2569507 :     for (; i <= db; i++)
     612       78149 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     613     2491358 :     for (; i <= m; i++)
     614           0 :       uel(C, i) = 0;
     615             :   }
     616       58330 :   for (; j <= ea; j++) {
     617           0 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     618           0 :     for (i = 1; i <= da; i++)
     619           0 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     620           0 :     for (; i <= m; i++)
     621           0 :       uel(C, i) = 0;
     622             :   }
     623       59794 :   for (; j <= eb; j++) {
     624        1464 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     625       85631 :     for (i = 1; i <= db; i++)
     626       84167 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     627        1464 :     for (; i <= m; i++)
     628           0 :       uel(C, i) = 0;
     629             :   }
     630       59794 :   for (; j <= n; j++)
     631        1464 :     gel(M, j) = zero_Flv(m);
     632       58330 :   return M;
     633             : }
     634             : 
     635             : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
     636             : 
     637             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     638             : static GEN
     639        8333 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
     640             : {
     641             :   pari_sp av;
     642        8333 :   long m1 = (m + 1)/2, m2 = m/2,
     643        8333 :     n1 = (n + 1)/2, n2 = n/2,
     644        8333 :     p1 = (p + 1)/2, p2 = p/2;
     645             :   GEN A11, A12, A22, B11, B21, B22,
     646             :     S1, S2, S3, S4, T1, T2, T3, T4,
     647             :     M1, M2, M3, M4, M5, M6, M7,
     648             :     V1, V2, V3, C;
     649             :   long j;
     650        8333 :   C = cgetg(p + 1, t_MAT);
     651      683213 :   for (j = 1; j <= p; j++)
     652      674880 :     gel(C, j) = cgetg(m + 1, t_VECSMALL);
     653        8333 :   av = avma;
     654        8333 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
     655        8333 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
     656        8333 :   M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
     657        8333 :   if (gc_needed(av, 1))
     658           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     659        8333 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
     660        8333 :   if (gc_needed(av, 1))
     661           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     662        8333 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
     663        8333 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
     664        8333 :   M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
     665        8333 :   if (gc_needed(av, 1))
     666           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     667        8333 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
     668        8333 :   if (gc_needed(av, 1))
     669           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     670        8333 :   A11 = matslice(A, 1, m1, 1, n1);
     671        8333 :   B11 = matslice(B, 1, n1, 1, p1);
     672        8333 :   M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
     673        8333 :   if (gc_needed(av, 1))
     674           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     675        8333 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     676        8333 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     677        8333 :   M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
     678        8333 :   if (gc_needed(av, 1))
     679           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     680        8333 :   add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
     681        8333 :   if (gc_needed(av, 1))
     682           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy M4 */
     683        8333 :   M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
     684        8333 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
     685        8333 :   if (gc_needed(av, 1))
     686           0 :     gerepileall(av, 6, &M2, &T3, &M3, &M1, &M5, &S4);  /* destroy S3 */
     687        8333 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
     688        8333 :   if (gc_needed(av, 1))
     689           0 :     gerepileall(av, 6, &M2, &M3, &M1, &M5, &S4, &T4);  /* destroy T3 */
     690        8333 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
     691        8333 :   if (gc_needed(av, 1))
     692           0 :     gerepileall(av, 5, &M2, &M3, &S4, &T4, &V1);  /* destroy M1, M5 */
     693        8333 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     694        8333 :   M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
     695        8333 :   if (gc_needed(av, 1))
     696           0 :     gerepileall(av, 5, &M2, &M3, &T4, &V1, &M6);  /* destroy S4, B22 */
     697        8333 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     698        8333 :   M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
     699        8333 :   if (gc_needed(av, 1))
     700           0 :     gerepileall(av, 5, &M2, &M3, &V1, &M6, &M7);  /* destroy A22, T4 */
     701        8333 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
     702        8333 :   add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
     703        8333 :   if (gc_needed(av, 1))
     704           0 :     gerepileall(av, 4, &M2, &M3, &V1, &M7);  /* destroy V3, M6 */
     705        8333 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
     706        8333 :   if (gc_needed(av, 1))
     707           0 :     gerepileall(av, 3, &M3, &M7, &V2);  /* destroy V1, M2 */
     708        8333 :   add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
     709        8333 :   add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
     710        8333 :   set_avma(av); return C;
     711             : }
     712             : 
     713             : /* Strassen-Winograd used for dim >= ZM_sw_bound */
     714             : static GEN
     715    27474264 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     716             : {
     717    27474264 :   ulong e = expu(p);
     718             : #ifdef LONG_IS_64BIT
     719    23073913 :   long ZM_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
     720             : #else
     721     4400954 :   long ZM_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
     722             : #endif
     723    27474867 :   if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
     724    27466534 :     return Flm_mul_classical(x, y, l, lx, ly, p, pi);
     725             :   else
     726        8333 :     return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
     727             : }
     728             : 
     729             : GEN
     730    13379405 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     731             : {
     732    13379405 :   long lx=lg(x), ly=lg(y);
     733    13379405 :   if (ly==1) return cgetg(1,t_MAT);
     734    13379405 :   if (lx==1) return zero_Flm(0, ly-1);
     735    13379405 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
     736             : }
     737             : 
     738             : GEN
     739    14035005 : Flm_mul(GEN x, GEN y, ulong p)
     740             : {
     741    14035005 :   long lx=lg(x), ly=lg(y);
     742    14035005 :   if (ly==1) return cgetg(1,t_MAT);
     743    14034718 :   if (lx==1) return zero_Flm(0, ly-1);
     744    14034718 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
     745             : }
     746             : 
     747             : struct _Flm
     748             : {
     749             :   ulong p;
     750             :   long n;
     751             : };
     752             : 
     753             : static GEN
     754        8470 : _Flm_mul(void *E , GEN x, GEN y)
     755        8470 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
     756             : static GEN
     757       43946 : _Flm_sqr(void *E, GEN x)
     758       43946 : { return Flm_mul(x,x,((struct _Flm*)E)->p); }
     759             : static GEN
     760         861 : _Flm_one(void *E)
     761         861 : { return matid_Flm(((struct _Flm*)E)->n); }
     762             : GEN
     763       24902 : Flm_powu(GEN x, ulong n, ulong p)
     764             : {
     765             :   struct _Flm d;
     766       24902 :   if (!n) return matid(lg(x)-1);
     767       24902 :   d.p = p;
     768       24902 :   return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
     769             : }
     770             : GEN
     771         861 : Flm_powers(GEN x, ulong n, ulong p)
     772             : {
     773             :   struct _Flm d;
     774         861 :   d.p = p;
     775         861 :   d.n = lg(x)-1;
     776         861 :   return gen_powers(x, n, 1, (void*)&d,
     777             :                     &_Flm_sqr, &_Flm_mul, &_Flm_one);
     778             : }
     779             : 
     780             : static GEN
     781           0 : _FpM_mul(void *p , GEN x, GEN y)
     782           0 : { return FpM_mul(x,y,(GEN)p); }
     783             : static GEN
     784           0 : _FpM_sqr(void *p, GEN x)
     785           0 : { return FpM_mul(x,x,(GEN)p); }
     786             : GEN
     787           0 : FpM_powu(GEN x, ulong n, GEN p)
     788             : {
     789           0 :   if (!n) return matid(lg(x)-1);
     790           0 :   if (lgefint(p) == 3)
     791             :   {
     792           0 :     pari_sp av = avma;
     793           0 :     ulong pp = uel(p,2);
     794             :     GEN z;
     795           0 :     if (pp == 2)
     796           0 :       z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
     797             :     else
     798           0 :       z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
     799           0 :     return gerepileupto(av, z);
     800             :   }
     801           0 :   return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
     802             : }
     803             : 
     804             : /*Multiple a column vector by a line vector to make a matrix*/
     805             : GEN
     806          14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
     807             : {
     808          14 :   long i,j, lx=lg(x), ly=lg(y);
     809             :   GEN z;
     810          14 :   if (ly==1) return cgetg(1,t_MAT);
     811          14 :   z = cgetg(ly, t_MAT);
     812          49 :   for (j=1; j < ly; j++)
     813             :   {
     814          35 :     GEN zj = cgetg(lx,t_VECSMALL);
     815         112 :     for (i=1; i<lx; i++)
     816          77 :       uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
     817          35 :     gel(z,j) = zj;
     818             :   }
     819          14 :   return z;
     820             : }
     821             : 
     822             : /*Multiple a column vector by a line vector to make a matrix*/
     823             : GEN
     824        5670 : FpC_FpV_mul(GEN x, GEN y, GEN p)
     825             : {
     826        5670 :   long i,j, lx=lg(x), ly=lg(y);
     827             :   GEN z;
     828        5670 :   if (ly==1) return cgetg(1,t_MAT);
     829        5670 :   z = cgetg(ly,t_MAT);
     830       63861 :   for (j=1; j < ly; j++)
     831             :   {
     832       58191 :     GEN zj = cgetg(lx,t_COL);
     833     1379126 :     for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
     834       58191 :     gel(z, j) = zj;
     835             :   }
     836        5670 :   return z;
     837             : }
     838             : 
     839             : /* Multiply a line vector by a column and return a scalar (t_INT) */
     840             : GEN
     841     8617102 : FpV_dotproduct(GEN x, GEN y, GEN p)
     842             : {
     843     8617102 :   long i, lx = lg(x);
     844             :   pari_sp av;
     845             :   GEN c;
     846     8617102 :   if (lx == 1) return gen_0;
     847     8616584 :   av = avma; c = mulii(gel(x,1),gel(y,1));
     848    27924766 :   for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
     849     8616550 :   return gerepileuptoint(av, modii(c,p));
     850             : }
     851             : GEN
     852           0 : FpV_dotsquare(GEN x, GEN p)
     853             : {
     854           0 :   long i, lx = lg(x);
     855             :   pari_sp av;
     856             :   GEN c;
     857           0 :   if (lx == 1) return gen_0;
     858           0 :   av = avma; c = sqri(gel(x,1));
     859           0 :   for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
     860           0 :   return gerepileuptoint(av, modii(c,p));
     861             : }
     862             : 
     863             : INLINE ulong
     864     8988408 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
     865             : {
     866     8988408 :   ulong c = uel(x,0) * uel(y,0);
     867             :   long k;
     868    64212770 :   for (k = 1; k < lx; k++) {
     869    55224362 :     c += uel(x,k) * uel(y,k);
     870    55224362 :     if (c & HIGHBIT) c %= p;
     871             :   }
     872     8988408 :   return c % p;
     873             : }
     874             : 
     875             : INLINE ulong
     876      739816 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
     877             : {
     878             :   ulong l0, l1, v1, h0, h1;
     879      739816 :   long i = 0;
     880             :   LOCAL_OVERFLOW;
     881             :   LOCAL_HIREMAINDER;
     882      739816 :   l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
     883    15409934 :   while (++i < lx) {
     884    14670118 :     l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
     885    14670118 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     886             :   }
     887      739816 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     888      464436 :   else return remlll_pre(v1, h1, l1, p, pi);
     889             : }
     890             : 
     891             : ulong
     892      356660 : Flv_dotproduct(GEN x, GEN y, ulong p)
     893             : {
     894      356660 :   long lx = lg(x)-1;
     895      356660 :   if (lx == 0) return 0;
     896      356660 :   if (SMALL_ULONG(p))
     897      356660 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     898             :   else
     899           0 :     return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
     900             : }
     901             : 
     902             : ulong
     903       58817 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
     904             : {
     905       58817 :   long lx = lg(x)-1;
     906       58817 :   if (lx == 0) return 0;
     907       58817 :   if (SMALL_ULONG(p))
     908       53259 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     909             :   else
     910        5558 :     return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
     911             : }
     912             : 
     913             : ulong
     914     9379910 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
     915             : {
     916     9379910 :   long lx = minss(lgpol(x), lgpol(y));
     917     9379860 :   if (lx == 0) return 0;
     918     9312752 :   if (pi)
     919      734258 :     return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
     920             :   else
     921     8578494 :     return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
     922             : }
     923             : ulong
     924           0 : Flx_dotproduct(GEN x, GEN y, ulong p)
     925           0 : { return Flx_dotproduct_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
     926             : 
     927             : GEN
     928     1659223 : FpM_FpC_mul(GEN x, GEN y, GEN p)
     929             : {
     930     1659223 :   long lx = lg(x);
     931     1659223 :   return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
     932             : }
     933             : GEN
     934      699335 : Flm_Flc_mul(GEN x, GEN y, ulong p)
     935             : {
     936      699335 :   long l, lx = lg(x);
     937      699335 :   if (lx==1) return cgetg(1,t_VECSMALL);
     938      699335 :   l = lgcols(x);
     939      699335 :   if (p==2)
     940      179842 :     return Flm_Flc_mul_i_2(x, y, lx, l);
     941      519493 :   else if (SMALL_ULONG(p))
     942      282580 :     return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
     943             :   else
     944      236913 :     return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
     945             : }
     946             : 
     947             : /* allow pi = 0 */
     948             : GEN
     949        6693 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     950             : {
     951        6693 :   long l, lx = lg(x);
     952             :   GEN z;
     953        6693 :   if (lx==1) return cgetg(1,t_VECSMALL);
     954        6693 :   l = lgcols(x);
     955        6693 :   z = cgetg(l, t_VECSMALL);
     956        6693 :   if (SMALL_ULONG(p))
     957        4407 :     __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     958             :   else
     959        2286 :     __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     960        6693 :   return z;
     961             : }
     962             : 
     963             : /* allow pi = 0 */
     964             : GEN
     965     7535058 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
     966             : {
     967     7535058 :   long l, lx = lg(x);
     968             :   GEN z;
     969     7535058 :   if (lx==1) return pol0_Flx(sv);
     970     7535058 :   l = lgcols(x);
     971     7533167 :   z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
     972     7532895 :   if (SMALL_ULONG(p))
     973     3304941 :     __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
     974             :   else
     975     4227954 :     __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
     976     7549341 :   return Flx_renormalize(z, l + 1);
     977             : }
     978             : 
     979             : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
     980             : GEN
     981     1418108 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
     982             : {
     983     1418108 :   long i, l, lx = lg(x);
     984             :   GEN z;
     985     1418108 :   if (lx==1) return pol_0(v);
     986     1418108 :   l = lgcols(x);
     987     1418108 :   z = new_chunk(l+1);
     988     2158029 :   for (i=l-1; i; i--)
     989             :   {
     990     1974759 :     pari_sp av = avma;
     991     1974759 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
     992     1974755 :     p1 = modii(p1, p);
     993     1974758 :     if (signe(p1))
     994             :     {
     995     1234837 :       if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
     996     1234837 :       gel(z,i+1) = gerepileuptoint(av, p1);
     997     1234838 :       break;
     998             :     }
     999      739921 :     set_avma(av);
    1000             :   }
    1001     1418108 :   if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
    1002     1234838 :   z[0] = evaltyp(t_POL) | _evallg(i+2);
    1003     1234838 :   z[1] = evalsigne(1) | evalvarn(v);
    1004     3597220 :   for (; i; i--)
    1005             :   {
    1006     2362382 :     pari_sp av = avma;
    1007     2362382 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
    1008     2362379 :     gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
    1009             :   }
    1010     1234838 :   return z;
    1011             : }
    1012             : 
    1013             : /********************************************************************/
    1014             : /**                                                                **/
    1015             : /**                           TRANSPOSITION                        **/
    1016             : /**                                                                **/
    1017             : /********************************************************************/
    1018             : 
    1019             : /* == zm_transpose */
    1020             : GEN
    1021      713524 : Flm_transpose(GEN x)
    1022             : {
    1023      713524 :   long i, dx, lx = lg(x);
    1024             :   GEN y;
    1025      713524 :   if (lx == 1) return cgetg(1,t_MAT);
    1026      713391 :   dx = lgcols(x); y = cgetg(dx,t_MAT);
    1027     8471691 :   for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
    1028      713392 :   return y;
    1029             : }
    1030             : 
    1031             : /********************************************************************/
    1032             : /**                                                                **/
    1033             : /**                           SCALAR MATRICES                      **/
    1034             : /**                                                                **/
    1035             : /********************************************************************/
    1036             : 
    1037             : GEN
    1038        4022 : gen_matid(long n, void *E, const struct bb_field *S)
    1039             : {
    1040        4022 :   GEN y = cgetg(n+1,t_MAT), _0, _1;
    1041             :   long i;
    1042        4022 :   if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
    1043        4022 :   _0 = S->s(E,0);
    1044        4022 :   _1 = S->s(E,1);
    1045       17903 :   for (i=1; i<=n; i++)
    1046             :   {
    1047       13881 :     GEN z = const_col(n, _0); gel(z,i) = _1;
    1048       13881 :     gel(y, i) = z;
    1049             :   }
    1050        4022 :   return y;
    1051             : }
    1052             : 
    1053             : GEN
    1054          35 : matid_F2xqM(long n, GEN T)
    1055             : {
    1056             :   void *E;
    1057          35 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1058          35 :   return gen_matid(n, E, S);
    1059             : }
    1060             : GEN
    1061          56 : matid_FlxqM(long n, GEN T, ulong p)
    1062             : {
    1063             :   void *E;
    1064          56 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1065          56 :   return gen_matid(n, E, S);
    1066             : }
    1067             : 
    1068             : GEN
    1069     1406522 : matid_Flm(long n)
    1070             : {
    1071     1406522 :   GEN y = cgetg(n+1,t_MAT);
    1072             :   long i;
    1073     1406522 :   if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
    1074     8818601 :   for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
    1075     1406500 :   return y;
    1076             : }
    1077             : 
    1078             : GEN
    1079          42 : scalar_Flm(long s, long n)
    1080             : {
    1081             :   long i;
    1082          42 :   GEN y = cgetg(n+1,t_MAT);
    1083         560 :   for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
    1084          42 :   return y;
    1085             : }
    1086             : 
    1087             : /********************************************************************/
    1088             : /**                                                                **/
    1089             : /**                           CONVERSIONS                          **/
    1090             : /**                                                                **/
    1091             : /********************************************************************/
    1092             : GEN
    1093    49457378 : ZV_to_Flv(GEN x, ulong p)
    1094   513968766 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
    1095             : 
    1096             : GEN
    1097    12312563 : ZM_to_Flm(GEN x, ulong p)
    1098    60768317 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
    1099             : 
    1100             : GEN
    1101        2418 : ZMV_to_FlmV(GEN x, ulong m)
    1102       18920 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
    1103             : 
    1104             : /*                          TO INTMOD                        */
    1105             : static GEN
    1106    20250441 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
    1107             : static GEN
    1108      383844 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
    1109             : 
    1110             : GEN
    1111        3439 : Fp_to_mod(GEN z, GEN p)
    1112             : {
    1113        3439 :   retmkintmod(modii(z, p), icopy(p));
    1114             : }
    1115             : 
    1116             : GEN
    1117      997241 : FpX_to_mod_raw(GEN z, GEN p)
    1118             : {
    1119      997241 :   long i,l = lg(z);
    1120             :   GEN x;
    1121      997241 :   if (l == 2)
    1122             :   {
    1123      522424 :     x = cgetg(3,t_POL); x[1] = z[1];
    1124      522424 :     gel(x,2) = mkintmod(gen_0,p); return x;
    1125             :   }
    1126      474817 :   x = cgetg(l,t_POL);
    1127     3918481 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1128      474817 :   x[1] = z[1]; return normalizepol_lg(x,l);
    1129             : }
    1130             : 
    1131             : /* z in Z[X], return z * Mod(1,p), normalized*/
    1132             : GEN
    1133     1668749 : FpX_to_mod(GEN z, GEN p)
    1134             : {
    1135     1668749 :   long i,l = lg(z);
    1136             :   GEN x;
    1137     1668749 :   if (l == 2)
    1138             :   {
    1139        3157 :     x = cgetg(3,t_POL); x[1] = z[1];
    1140        3157 :     gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
    1141             :   }
    1142     1665592 :   x = cgetg(l,t_POL);
    1143     1665592 :   if (l >2) p = icopy(p);
    1144    18369436 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1145     1665592 :   x[1] = z[1]; return normalizepol_lg(x,l);
    1146             : }
    1147             : 
    1148             : GEN
    1149       84021 : FpXC_to_mod(GEN z, GEN p)
    1150             : {
    1151       84021 :   long i,l = lg(z);
    1152       84021 :   GEN x = cgetg(l, t_COL);
    1153       84021 :   p = icopy(p);
    1154      474166 :   for (i=1; i<l; i++)
    1155      390145 :     gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
    1156       84021 :   return x;
    1157             : }
    1158             : 
    1159             : static GEN
    1160           0 : FpXC_to_mod_raw(GEN x, GEN p)
    1161           0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
    1162             : 
    1163             : GEN
    1164           0 : FpXM_to_mod(GEN z, GEN p)
    1165             : {
    1166           0 :   long i,l = lg(z);
    1167           0 :   GEN x = cgetg(l, t_MAT);
    1168           0 :   p = icopy(p);
    1169           0 :   for (i=1; i<l; i++)
    1170           0 :     gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
    1171           0 :   return x;
    1172             : }
    1173             : 
    1174             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1175             : GEN
    1176       35928 : FpV_to_mod(GEN z, GEN p)
    1177             : {
    1178       35928 :   long i,l = lg(z);
    1179       35928 :   GEN x = cgetg(l, t_VEC);
    1180       35928 :   if (l == 1) return x;
    1181       35928 :   p = icopy(p);
    1182       72255 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1183       35928 :   return x;
    1184             : }
    1185             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1186             : GEN
    1187         105 : FpC_to_mod(GEN z, GEN p)
    1188             : {
    1189         105 :   long i, l = lg(z);
    1190         105 :   GEN x = cgetg(l, t_COL);
    1191         105 :   if (l == 1) return x;
    1192          91 :   p = icopy(p);
    1193         805 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1194          91 :   return x;
    1195             : }
    1196             : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
    1197             : GEN
    1198         226 : FpM_to_mod(GEN z, GEN p)
    1199             : {
    1200         226 :   long i, j, m, l = lg(z);
    1201         226 :   GEN  x = cgetg(l,t_MAT), y, zi;
    1202         226 :   if (l == 1) return x;
    1203         205 :   m = lgcols(z);
    1204         205 :   p = icopy(p);
    1205        2456 :   for (i=1; i<l; i++)
    1206             :   {
    1207        2251 :     gel(x,i) = cgetg(m,t_COL);
    1208        2251 :     y = gel(x,i); zi= gel(z,i);
    1209       66799 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1210             :   }
    1211         205 :   return x;
    1212             : }
    1213             : GEN
    1214          28 : Flc_to_mod(GEN z, ulong pp)
    1215             : {
    1216          28 :   long i, l = lg(z);
    1217          28 :   GEN p, x = cgetg(l, t_COL);
    1218          28 :   if (l == 1) return x;
    1219          28 :   p = utoipos(pp);
    1220         112 :   for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
    1221          28 :   return x;
    1222             : }
    1223             : GEN
    1224       56292 : Flm_to_mod(GEN z, ulong pp)
    1225             : {
    1226       56292 :   long i, j, m, l = lg(z);
    1227       56292 :   GEN p, x = cgetg(l,t_MAT), y, zi;
    1228       56292 :   if (l == 1) return x;
    1229       56264 :   m = lgcols(z);
    1230       56264 :   p = utoipos(pp);
    1231      179348 :   for (i=1; i<l; i++)
    1232             :   {
    1233      123084 :     gel(x,i) = cgetg(m,t_COL);
    1234      123084 :     y = gel(x,i); zi= gel(z,i);
    1235      506844 :     for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
    1236             :   }
    1237       56264 :   return x;
    1238             : }
    1239             : 
    1240             : GEN
    1241         574 : FpVV_to_mod(GEN z, GEN p)
    1242             : {
    1243         574 :   long i, j, m, l = lg(z);
    1244         574 :   GEN  x = cgetg(l,t_VEC), y, zi;
    1245         574 :   if (l == 1) return x;
    1246         574 :   m = lgcols(z);
    1247         574 :   p = icopy(p);
    1248        1246 :   for (i=1; i<l; i++)
    1249             :   {
    1250         672 :     gel(x,i) = cgetg(m,t_VEC);
    1251         672 :     y = gel(x,i); zi= gel(z,i);
    1252        2016 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1253             :   }
    1254         574 :   return x;
    1255             : }
    1256             : 
    1257             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1258             : GEN
    1259           7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
    1260             : {
    1261           7 :   long i,l = lg(z);
    1262           7 :   GEN x = cgetg(l, t_COL);
    1263           7 :   if (l == 1) return x;
    1264           7 :   p = icopy(p);
    1265           7 :   T = FpX_to_mod_raw(T, p);
    1266          21 :   for (i=1; i<l; i++)
    1267          14 :     gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
    1268           7 :   return x;
    1269             : }
    1270             : 
    1271             : static GEN
    1272      582533 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
    1273             : {
    1274      582533 :   GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
    1275      582533 :   return mkpolmod(z, T);
    1276             : }
    1277             : 
    1278             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1279             : GEN
    1280          28 : FqC_to_mod(GEN z, GEN T, GEN p)
    1281             : {
    1282             :   GEN x;
    1283          28 :   long i,l = lg(z);
    1284          28 :   if (!T) return FpC_to_mod(z, p);
    1285          28 :   x = cgetg(l, t_COL);
    1286          28 :   if (l == 1) return x;
    1287          28 :   p = icopy(p);
    1288          28 :   T = FpX_to_mod_raw(T, p);
    1289         504 :   for (i=1; i<l; i++)
    1290         476 :     gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
    1291          28 :   return x;
    1292             : }
    1293             : 
    1294             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1295             : static GEN
    1296      107975 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
    1297      690032 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
    1298             : 
    1299             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1300             : GEN
    1301       26145 : FqM_to_mod(GEN z, GEN T, GEN p)
    1302             : {
    1303             :   GEN x;
    1304       26145 :   long i,l = lg(z);
    1305       26145 :   if (!T) return FpM_to_mod(z, p);
    1306       26145 :   x = cgetg(l, t_MAT);
    1307       26145 :   if (l == 1) return x;
    1308       26117 :   p = icopy(p);
    1309       26117 :   T = FpX_to_mod_raw(T, p);
    1310      134092 :   for (i=1; i<l; i++)
    1311      107975 :     gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
    1312       26117 :   return x;
    1313             : }
    1314             : 
    1315             : /********************************************************************/
    1316             : /*                                                                  */
    1317             : /*                     Blackbox linear algebra                      */
    1318             : /*                                                                  */
    1319             : /********************************************************************/
    1320             : 
    1321             : /* A sparse column (zCs) v is a t_COL with two components C and E which are
    1322             :  * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
    1323             :  * (e_j) is the canonical basis.
    1324             :  * A sparse matrix (zMs) is a t_VEC of zCs */
    1325             : 
    1326             : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
    1327             :  * integer representing an element of Fp. This is important since p can be
    1328             :  * large and p+E[i] would not fit in a C long.  */
    1329             : 
    1330             : /* RgCs and RgMs are similar, except that the type of the component is
    1331             :  * unspecified. Functions handling RgCs/RgMs must be independent of the type
    1332             :  * of E. */
    1333             : 
    1334             : /* Most functions take an argument nbrow which is the number of lines of the
    1335             :  * column/matrix, which cannot be derived from the data. */
    1336             : 
    1337             : GEN
    1338           0 : zCs_to_ZC(GEN R, long nbrow)
    1339             : {
    1340           0 :   GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
    1341           0 :   long j, l = lg(C);
    1342           0 :   for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
    1343           0 :   return c;
    1344             : }
    1345             : 
    1346             : GEN
    1347           0 : zMs_to_ZM(GEN x, long nbrow)
    1348           0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
    1349             : 
    1350             : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
    1351             :  * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
    1352             : GEN
    1353         119 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
    1354             : {
    1355         119 :   pari_sp ltop = avma;
    1356         119 :   long col = 0, n = lg(B)-1, m = 2*n+1;
    1357         119 :   if (ZV_equal0(B)) return zerocol(n);
    1358         119 :   while (++col <= n)
    1359             :   {
    1360         119 :     pari_sp btop = avma, av;
    1361             :     long i, lQ;
    1362         119 :     GEN V, Q, M, W = B;
    1363         119 :     GEN b = cgetg(m+2, t_POL);
    1364         119 :     b[1] = evalsigne(1)|evalvarn(0);
    1365         119 :     gel(b, 2) = gel(W, col);
    1366       49951 :     for (i = 3; i<m+2; i++)
    1367       49832 :       gel(b, i) = cgeti(lgefint(p));
    1368         119 :     av = avma;
    1369       49951 :     for (i = 3; i<m+2; i++)
    1370             :     {
    1371       49832 :       W = f(E, W);
    1372       49832 :       affii(gel(W, col),gel(b, i));
    1373       49832 :       if (gc_needed(av,1))
    1374             :       {
    1375        1084 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
    1376        1084 :         W = gerepileupto(av, W);
    1377             :       }
    1378             :     }
    1379         119 :     b = FpX_renormalize(b, m+2);
    1380         119 :     if (lgpol(b)==0) {set_avma(btop); continue; }
    1381         119 :     M = FpX_halfgcd(b, pol_xn(m, 0), p);
    1382         119 :     Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
    1383         119 :     W = B; lQ =lg(Q);
    1384         119 :     if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
    1385         119 :     V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
    1386         119 :     av = avma;
    1387       24587 :     for (i = lQ-3; i > 1; i--)
    1388             :     {
    1389       24468 :       W = f(E, W);
    1390       24468 :       V = ZC_lincomb(gen_1, gel(Q,i), V, W);
    1391       24468 :       if (gc_needed(av,1))
    1392             :       {
    1393         864 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
    1394         864 :         gerepileall(av, 2, &V, &W);
    1395             :       }
    1396             :     }
    1397         119 :     V = FpC_red(V, p);
    1398         119 :     W = FpC_sub(f(E,V), B, p);
    1399         119 :     if (ZV_equal0(W)) return gerepilecopy(ltop, V);
    1400           0 :     av = avma;
    1401           0 :     for (i = 1; i <= n; ++i)
    1402             :     {
    1403           0 :       V = W;
    1404           0 :       W = f(E, V);
    1405           0 :       if (ZV_equal0(W))
    1406           0 :         return gerepilecopy(ltop, shallowtrans(V));
    1407           0 :       gerepileall(av, 2, &V, &W);
    1408             :     }
    1409           0 :     set_avma(btop);
    1410             :   }
    1411           0 :   return NULL;
    1412             : }
    1413             : 
    1414             : GEN
    1415           0 : zMs_ZC_mul(GEN M, GEN B)
    1416             : {
    1417             :   long i, j;
    1418           0 :   long n = lg(B)-1;
    1419           0 :   GEN V = zerocol(n);
    1420           0 :   for (i = 1; i <= n; ++i)
    1421           0 :     if (signe(gel(B, i)))
    1422             :     {
    1423           0 :       GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1424           0 :       long l = lg(C);
    1425           0 :       for (j = 1; j < l; ++j)
    1426             :       {
    1427           0 :         long k = C[j];
    1428           0 :         switch(E[j])
    1429             :         {
    1430           0 :         case 1:
    1431           0 :           gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
    1432           0 :           break;
    1433           0 :         case -1:
    1434           0 :           gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
    1435           0 :           break;
    1436           0 :         default:
    1437           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]));
    1438           0 :           break;
    1439             :         }
    1440             :       }
    1441             :     }
    1442           0 :   return V;
    1443             : }
    1444             : 
    1445             : GEN
    1446           0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
    1447             : 
    1448             : GEN
    1449       74538 : ZV_zMs_mul(GEN B, GEN M)
    1450             : {
    1451             :   long i, j;
    1452       74538 :   long m = lg(M)-1;
    1453       74538 :   GEN V = cgetg(m+1,t_VEC);
    1454    38968562 :   for (i = 1; i <= m; ++i)
    1455             :   {
    1456    38894024 :     GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1457    38894024 :     long l = lg(C);
    1458             :     GEN z;
    1459    38894024 :     if (l == 1)
    1460             :     {
    1461           0 :       gel(V,i) = gen_0;
    1462           0 :       continue;
    1463             :     }
    1464    38894024 :     z = mulis(gel(B, C[1]), E[1]);
    1465   327053616 :     for (j = 2; j < l; ++j)
    1466             :     {
    1467   288159592 :       long k = C[j];
    1468   288159592 :       switch(E[j])
    1469             :       {
    1470   194218066 :       case 1:
    1471   194218066 :         z = addii(z, gel(B,k));
    1472   194218066 :         break;
    1473    78567695 :       case -1:
    1474    78567695 :         z = subii(z, gel(B,k));
    1475    78567695 :         break;
    1476    15373831 :       default:
    1477    15373831 :         z = addii(z, mulis(gel(B,k), E[j]));
    1478    15373831 :         break;
    1479             :       }
    1480             :     }
    1481    38894024 :     gel(V,i) = z;
    1482             :   }
    1483       74538 :   return V;
    1484             : }
    1485             : 
    1486             : GEN
    1487         119 : FpV_FpMs_mul(GEN B, GEN M, GEN p) { return FpV_red(ZV_zMs_mul(B, M), p); }
    1488             : 
    1489             : GEN
    1490           0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
    1491             : {
    1492           0 :   pari_sp av = avma, av2;
    1493           0 :   GEN xi, xb, pi = gen_1, P;
    1494             :   long i;
    1495           0 :   if (!C) {
    1496           0 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    1497           0 :     if (!C) pari_err_INV("ZlM_gauss", a);
    1498             :   }
    1499           0 :   P = utoipos(p);
    1500           0 :   av2 = avma;
    1501           0 :   xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1502           0 :   xb = Flm_to_ZM(xi);
    1503           0 :   for (i = 2; i <= e; i++)
    1504             :   {
    1505           0 :     pi = muliu(pi, p); /* = p^(i-1) */
    1506           0 :     b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
    1507           0 :     if (gc_needed(av,2))
    1508             :     {
    1509           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
    1510           0 :       gerepileall(av2,3, &pi,&b,&xb);
    1511             :     }
    1512           0 :     xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1513           0 :     xb = ZM_add(xb, nm_Z_mul(xi, pi));
    1514             :   }
    1515           0 :   return gerepileupto(av, xb);
    1516             : }
    1517             : 
    1518             : struct wrapper_modp_s {
    1519             :   GEN (*f)(void*, GEN);
    1520             :   void *E;
    1521             :   GEN p;
    1522             : };
    1523             : 
    1524             : /* compute f(x) mod p */
    1525             : static GEN
    1526       74419 : wrap_relcomb_modp(void *E, GEN x)
    1527             : {
    1528       74419 :   struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
    1529       74419 :   return FpC_red(W->f(W->E, x), W->p);
    1530             : }
    1531             : static GEN
    1532           0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
    1533             : 
    1534             : static GEN
    1535       74419 : wrap_relker(void*E, GEN x) { return ZV_zMs_mul(x, (GEN)E); }
    1536             : 
    1537             : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
    1538             : GEN
    1539           0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
    1540             : {
    1541             :   struct wrapper_modp_s W;
    1542           0 :   pari_sp av = avma;
    1543           0 :   GEN xb, xi, pi = gen_1;
    1544             :   long i;
    1545           0 :   W.E = E;
    1546           0 :   W.f = f;
    1547           0 :   W.p = p;
    1548           0 :   xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
    1549           0 :   if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
    1550           0 :   xb = xi;
    1551           0 :   for (i = 2; i <= e; i++)
    1552             :   {
    1553           0 :     pi = mulii(pi, p); /* = p^(i-1) */
    1554           0 :     B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
    1555           0 :     if (gc_needed(av,2))
    1556             :     {
    1557           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
    1558           0 :       gerepileall(av,3, &pi,&B,&xb);
    1559             :     }
    1560           0 :     xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
    1561           0 :     if (!xi) return NULL;
    1562           0 :     if (typ(xi) == t_VEC) return gerepileupto(av, xi);
    1563           0 :     xb = ZC_add(xb, ZC_Z_mul(xi, pi));
    1564             :   }
    1565           0 :   return gerepileupto(av, xb);
    1566             : }
    1567             : 
    1568             : static GEN
    1569       24916 : vecprow(GEN A, GEN prow)
    1570             : {
    1571       24916 :   return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
    1572             : }
    1573             : 
    1574             : /* Solve the equation MX = A. Return either a solution as a t_COL,
    1575             :  * or the index of a column which is linearly dependent from the others as a
    1576             :  * t_VECSMALL with a single component. */
    1577             : GEN
    1578           0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
    1579             : {
    1580           0 :   pari_sp av = avma;
    1581             :   GEN pcol, prow;
    1582           0 :   long nbi=lg(M)-1, lR;
    1583             :   long i, n;
    1584             :   GEN Mp, Ap, Rp;
    1585             :   pari_timer ti;
    1586           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1587           0 :   RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
    1588           0 :   if (!pcol) return gc_NULL(av);
    1589           0 :   if (DEBUGLEVEL)
    1590           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
    1591           0 :   n = lg(pcol)-1;
    1592           0 :   Mp = cgetg(n+1, t_MAT);
    1593           0 :   for(i=1; i<=n; i++)
    1594           0 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1595           0 :   Ap = zCs_to_ZC(vecprow(A, prow), n);
    1596           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1597           0 :   Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
    1598           0 :   if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
    1599           0 :   if (!Rp) return gc_NULL(av);
    1600           0 :   lR = lg(Rp)-1;
    1601           0 :   if (typ(Rp) == t_COL)
    1602             :   {
    1603           0 :     GEN R = zerocol(nbi+1);
    1604           0 :     for (i=1; i<=lR; i++)
    1605           0 :        gel(R,pcol[i]) = gel(Rp,i);
    1606           0 :     return gerepilecopy(av, R);
    1607             :   }
    1608           0 :   for (i = 1; i <= lR; ++i)
    1609           0 :     if (signe(gel(Rp, i)))
    1610           0 :       return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
    1611             :   return NULL; /* LCOV_EXCL_LINE */
    1612             : }
    1613             : 
    1614             : GEN
    1615           0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
    1616             : {
    1617           0 :   return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1618             : }
    1619             : 
    1620             : GEN
    1621           0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
    1622             : {
    1623             :   GEN res;
    1624           0 :   pari_CATCH(e_INV)
    1625             :   {
    1626           0 :     GEN E = pari_err_last();
    1627           0 :     GEN x = err_get_compo(E,2);
    1628           0 :     if (typ(x) != t_INTMOD) pari_err(0,E);
    1629           0 :     if (DEBUGLEVEL)
    1630           0 :       pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
    1631           0 :     res = NULL;
    1632             :   } pari_TRY {
    1633           0 :     res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1634           0 :   } pari_ENDCATCH
    1635           0 :   return res;
    1636             : }
    1637             : 
    1638             : static GEN
    1639         119 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
    1640             : {
    1641         119 :   long i, j, oldf = 0, f = 0;
    1642         119 :   long lrow = lg(prow), lM = lg(M);
    1643         119 :   GEN W = const_vecsmall(lM-1,1);
    1644         119 :   GEN R = cgetg(lrow, t_VEC);
    1645      235298 :   for (i=1; i<lrow; i++)
    1646      235179 :     gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
    1647             :   do
    1648             :   {
    1649         393 :     oldf = f;
    1650      436630 :     for(i=1; i<lM; i++)
    1651      436237 :       if (W[i])
    1652             :       {
    1653      147864 :         GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
    1654      147864 :         long c=0, cj=0, lF = lg(F);
    1655     1294834 :         for(j=1; j<lF; j++)
    1656     1146970 :           if (!gel(R,F[j])) { c++; cj=j; }
    1657      147864 :         if (c>=2) continue;
    1658      127151 :         if (c==1)
    1659             :         {
    1660       41421 :           pari_sp av = avma;
    1661       41421 :           GEN s = gen_0;
    1662      359165 :           for(j=1; j<lF; j++)
    1663      317744 :             if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
    1664             :           /* s /= -E[cj] mod p */
    1665       41421 :           s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
    1666       41421 :           gel(R,F[cj]) = gerepileuptoint(av, s);
    1667       41421 :           f++;
    1668             :         }
    1669      127151 :         W[i]=0;
    1670             :       }
    1671         393 :   } while(oldf!=f);
    1672      235298 :   for (i=1; i<lrow; i++)
    1673      235179 :     if (!gel(R,i)) gel(R,i) = gen_0;
    1674         119 :   return R;
    1675             : }
    1676             : 
    1677             : /* Return a linear form Y such that YM is essentially 0 */
    1678             : GEN
    1679         119 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
    1680             : {
    1681         119 :   pari_sp av = avma, av2;
    1682             :   GEN pcol, prow;
    1683             :   long i, n;
    1684             :   GEN Mp, B, MB, R, Rp;
    1685             :   pari_timer ti;
    1686             :   struct wrapper_modp_s W;
    1687         119 :   if (DEBUGLEVEL) timer_start(&ti);
    1688         119 :   RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
    1689         119 :   if (!pcol) return gc_NULL(av);
    1690         119 :   if (DEBUGLEVEL)
    1691           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
    1692         119 :   n = lg(pcol)-1;
    1693         119 :   Mp = cgetg(n+1, t_MAT);
    1694       25035 :   for(i=1; i<=n; i++)
    1695       24916 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1696         119 :   W.E = (void*) Mp;
    1697         119 :   W.f = wrap_relker;
    1698         119 :   W.p = p;
    1699         119 :   av2 = avma;
    1700             :   for(;;)
    1701             :   {
    1702         119 :     set_avma(av2);
    1703         119 :     B = random_FpV(n, p);
    1704         119 :     MB = FpV_FpMs_mul(B, Mp, p);
    1705         119 :     if (DEBUGLEVEL) timer_start(&ti);
    1706         119 :     pari_CATCH(e_INV)
    1707             :     {
    1708           0 :       GEN E = pari_err_last();
    1709           0 :       GEN x = err_get_compo(E,2);
    1710           0 :       if (typ(x) != t_INTMOD) pari_err(0,E);
    1711           0 :       if (DEBUGLEVEL)
    1712           0 :         pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
    1713           0 :       Rp = NULL;
    1714             :     } pari_TRY {
    1715         119 :       Rp = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, MB, p);
    1716         119 :     } pari_ENDCATCH
    1717         119 :     if (!Rp) continue;
    1718         119 :     if (typ(Rp)==t_COL)
    1719             :     {
    1720         119 :       Rp = FpC_sub(Rp,B,p);
    1721         119 :       if (ZV_equal0(Rp)) continue;
    1722             :     }
    1723         119 :     R = FpMs_structelim_back(M, Rp, prow, p);
    1724         119 :     if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
    1725         119 :     return gerepilecopy(av, R);
    1726             :   }
    1727             : }
    1728             : 
    1729             : GEN
    1730          42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
    1731             : {
    1732          42 :   return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
    1733             : }

Generated by: LCOV version 1.14