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 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.10.0 lcov report (development 21501-1931cb9) Lines: 722 936 77.1 %
Date: 2017-12-16 06:20:36 Functions: 105 131 80.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /********************************************************************/
      18             : /**                                                                **/
      19             : /**                           REDUCTION                            **/
      20             : /**                                                                **/
      21             : /********************************************************************/
      22             : /* z in Z^n, return lift(Col(z) * Mod(1,p)) */
      23             : GEN
      24     9240466 : FpC_red(GEN x, GEN p)
      25     9240466 : { pari_APPLY_type(t_COL, modii(gel(x,i), p)) }
      26             : 
      27             : /* z in Z^n, return lift(Vec(z) * Mod(1,p)) */
      28             : GEN
      29      207823 : FpV_red(GEN x, GEN p)
      30      207823 : { pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
      31             : 
      32             : GEN
      33     1028290 : FpC_center(GEN x, GEN p, GEN pov2)
      34     1028290 : { 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      111692 : Flv_center(GEN x, ulong p, ulong ps2)
      51      111692 : { 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     1773976 : FpM_red(GEN x, GEN p)
      56     1773976 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
      57             : 
      58             : GEN
      59      134926 : FpM_center(GEN x, GEN p, GEN pov2)
      60      134926 : { 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        7896 : Flm_center(GEN x, ulong p, ulong ps2)
      70        7896 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
      71             : 
      72             : GEN
      73         119 : random_FpV(long d, GEN p)
      74             : {
      75             :   long i;
      76         119 :   GEN y = cgetg(d+1,t_VEC);
      77         119 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
      78         119 :   return y;
      79             : }
      80             : 
      81             : GEN
      82         182 : random_FpC(long d, GEN p)
      83             : {
      84             :   long i;
      85         182 :   GEN y = cgetg(d+1,t_COL);
      86         182 :   for (i=1; i<=d; i++) gel(y,i) = randomi(p);
      87         182 :   return y;
      88             : }
      89             : 
      90             : GEN
      91        3416 : random_Flv(long n, ulong p)
      92             : {
      93        3416 :   GEN y = cgetg(n+1, t_VECSMALL);
      94             :   long i;
      95        3416 :   for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
      96        3416 :   return y;
      97             : }
      98             : 
      99             : /********************************************************************/
     100             : /**                                                                **/
     101             : /**                           ADD, SUB                             **/
     102             : /**                                                                **/
     103             : /********************************************************************/
     104             : GEN
     105      137144 : FpC_add(GEN x, GEN y, GEN p)
     106      137144 : { 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      629007 : Flv_add(GEN x, GEN y, ulong p)
     118             : {
     119      629007 :   if (p==2)
     120       87885 :     pari_APPLY_ulong(x[i]^y[i])
     121             :   else
     122      541122 :     pari_APPLY_ulong(Fl_add(x[i], y[i], p))
     123             : }
     124             : 
     125             : void
     126      168344 : Flv_add_inplace(GEN x, GEN y, ulong p)
     127             : {
     128      168344 :   long i, l = lg(x);
     129      168344 :   if (p==2)
     130      165978 :     for (i = 1; i < l; i++) x[i] ^= y[i];
     131             :   else
     132        2366 :     for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
     133      168344 : }
     134             : 
     135             : ulong
     136        3822 : Flv_sum(GEN x, ulong p)
     137             : {
     138        3822 :   long i, l = lg(x);
     139        3822 :   ulong s = 0;
     140        3822 :   if (p==2)
     141        3822 :     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        3822 :   return s;
     145             : }
     146             : 
     147             : GEN
     148      199507 : FpC_sub(GEN x, GEN y, GEN p)
     149      199507 : { 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    84083200 : Flv_sub(GEN x, GEN y, ulong p)
     161    84083200 : { 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        8092 : Flm_Fl_add(GEN x, ulong y, ulong p)
     172             : {
     173        8092 :   long l = lg(x), i, j;
     174        8092 :   GEN z = cgetg(l,t_MAT);
     175             : 
     176        8092 :   if (l==1) return z;
     177        7987 :   if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
     178       32606 :   for (i=1; i<l; i++)
     179             :   {
     180       24619 :     GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
     181       24619 :     gel(z,i) = zi;
     182       24619 :     for (j=1; j<l; j++) zi[j] = xi[j];
     183       24619 :     zi[i] = Fl_add(zi[i], y, p);
     184             :   }
     185        7987 :   return z;
     186             : }
     187             : 
     188             : GEN
     189        6993 : Flm_add(GEN x, GEN y, ulong p)
     190        6993 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
     191             : 
     192             : GEN
     193     8711620 : Flm_sub(GEN x, GEN y, ulong p)
     194     8711620 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
     195             : 
     196             : /********************************************************************/
     197             : /**                                                                **/
     198             : /**                           MULTIPLICATION                       **/
     199             : /**                                                                **/
     200             : /********************************************************************/
     201             : GEN
     202      622335 : FpC_Fp_mul(GEN x, GEN y, GEN p)
     203      622335 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
     204             : 
     205             : GEN
     206      433668 : Flv_Fl_mul(GEN x, ulong y, ulong p)
     207      433668 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
     208             : 
     209             : GEN
     210         422 : Flv_Fl_div(GEN x, ulong y, ulong p)
     211             : {
     212         422 :   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        1610 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
     221        1610 :   long i, j, h, l = lg(X);
     222        1610 :   GEN A = cgetg(l, t_MAT);
     223        1610 :   if (l == 1) return A;
     224        1610 :   h = lgcols(X);
     225       12579 :   for (j=1; j<l; j++)
     226             :   {
     227       10969 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     228       10969 :     for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
     229       10969 :     gel(A,j) = a;
     230             :   }
     231        1610 :   return A;
     232             : }
     233             : 
     234             : /* x *= y */
     235             : void
     236     1539717 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
     237             : {
     238             :   long i;
     239     1539717 :   for (i=1;i<=l;i++) x[i] = Fl_mul(x[i], y, p);
     240     1539717 : }
     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             : GEN
     262     7861776 : Flm_Fl_mul(GEN y, ulong x, ulong p)
     263             : {
     264     7861776 :   long i, j, m, l = lg(y);
     265     7861776 :   GEN z = cgetg(l, t_MAT);
     266     7860940 :   if (l == 1) return z;
     267     7860940 :   m = lgcols(y);
     268     7861809 :   if (HIGHWORD(x | p))
     269    60234377 :     for(j=1; j<l; j++) {
     270    55216321 :       GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     271    55038717 :       for(i=1; i<m; i++) c[i] = Fl_mul(ucoeff(y,i,j), x, p);
     272             :     }
     273             :   else
     274    15719285 :     for(j=1; j<l; j++) {
     275    12875315 :       GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
     276    12875315 :       for(i=1; i<m; i++) c[i] = (ucoeff(y,i,j) * x) % p;
     277             :     }
     278     7862026 :   return z;
     279             : }
     280             : 
     281             : GEN
     282     1812190 : Flv_neg(GEN x, ulong p)
     283     1812190 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
     284             : 
     285             : void
     286        5289 : Flv_neg_inplace(GEN v, ulong p)
     287             : {
     288             :   long i;
     289      167660 :   for (i = 1; i < lg(v); ++i)
     290      162371 :     v[i] = Fl_neg(v[i], p);
     291        5289 : }
     292             : 
     293             : GEN
     294      168517 : Flm_neg(GEN x, ulong p)
     295      168517 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
     296             : 
     297             : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
     298             : static GEN
     299    11490393 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
     300             : {
     301    11490393 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1));
     302             :   long k;
     303   174017300 :   for (k = 2; k < lx; k++)
     304             :   {
     305   162526907 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     306   162526907 :     if (signe(t)) c = addii(c, t);
     307             :   }
     308    11490393 :   return c;
     309             : }
     310             : 
     311             : static long
     312    21635523 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
     313             : {
     314             :   long k;
     315    21635523 :   long c = coeff(x,i,1) * y[1];
     316   289227596 :   for (k = 2; k < lx; k++)
     317   267592073 :     c += coeff(x,i,k) * y[k];
     318    21635523 :   return c;
     319             : }
     320             : 
     321             : GEN
     322     1653225 : zm_zc_mul(GEN x, GEN y)
     323             : {
     324     1653225 :   long lx = lg(x), l, i;
     325             :   GEN z;
     326     1653225 :   if (lx == 1) return cgetg(1, t_VECSMALL);
     327     1653225 :   l = lg(gel(x,1));
     328     1653225 :   z = cgetg(l,t_VECSMALL);
     329     1653225 :   for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
     330     1653225 :   return z;
     331             : }
     332             : 
     333             : GEN
     334         252 : zm_mul(GEN x, GEN y)
     335             : {
     336         252 :   long i,j,lx=lg(x), ly=lg(y);
     337             :   GEN z;
     338         252 :   if (ly==1) return cgetg(1,t_MAT);
     339         252 :   z = cgetg(ly,t_MAT);
     340         252 :   if (lx==1)
     341             :   {
     342           0 :     for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
     343           0 :     return z;
     344             :   }
     345        2128 :   for (j=1; j<ly; j++)
     346        1876 :     gel(z,j) = zm_zc_mul(x, gel(y,j));
     347         252 :   return z;
     348             : }
     349             : 
     350             : static ulong
     351   194883587 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
     352             : {
     353   194883587 :   ulong c = ucoeff(x,i,1) * uel(y,1);
     354             :   long k;
     355  3201147172 :   for (k = 2; k < lx; k++) {
     356  3006263585 :     c += ucoeff(x,i,k) * uel(y,k);
     357  3006263585 :     if (c & HIGHBIT) c %= p;
     358             :   }
     359   194883587 :   return c % p;
     360             : }
     361             : 
     362             : static ulong
     363   188902030 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
     364             : {
     365             :   ulong l0, l1, v1, h0, h1;
     366   188902030 :   long k = 1;
     367             :   LOCAL_OVERFLOW;
     368             :   LOCAL_HIREMAINDER;
     369   188902030 :   l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
     370  2093327253 :   while (++k < lx) {
     371  1715523193 :     l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
     372  1715523193 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     373             :   }
     374   188902030 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     375     4083577 :   else return remlll_pre(v1, h1, l1, p, pi);
     376             : }
     377             : 
     378             : static GEN
     379       82593 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
     380             : {
     381             :   long i,j;
     382       82593 :   GEN z = NULL;
     383             : 
     384     1064119 :   for (j=1; j<lx; j++)
     385             :   {
     386      981526 :     if (!y[j]) continue;
     387      306377 :     if (!z) z = Flv_copy(gel(x,j));
     388      231449 :     else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
     389             :   }
     390       82593 :   if (!z) z = zero_zv(l-1);
     391       82593 :   return z;
     392             : }
     393             : 
     394             : static GEN
     395      941465 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
     396             : {
     397      941465 :   GEN z = cgetg(l,t_COL);
     398             :   long i;
     399    11501439 :   for (i = 1; i < l; i++)
     400             :   {
     401    10559974 :     pari_sp av = avma;
     402    10559974 :     GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
     403    10559974 :     gel(z,i) = gerepileuptoint(av, modii(c,p));
     404             :   }
     405      941465 :   return z;
     406             : }
     407             : 
     408             : static void
     409    23460465 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
     410             : {
     411             :   long i;
     412    23460465 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
     413    23451361 : }
     414             : static GEN
     415    20314877 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
     416             : {
     417    20314877 :   GEN z = cgetg(l,t_VECSMALL);
     418    20314877 :   __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     419    20314877 :   return z;
     420             : }
     421             : 
     422             : static void
     423    37723202 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     424             : {
     425             :   long i;
     426    37723202 :   for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
     427    37498265 : }
     428             : static GEN
     429    33312567 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
     430             : {
     431    33312567 :   GEN z = cgetg(l,t_VECSMALL);
     432    33303849 :   __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     433    33294340 :   return z;
     434             : }
     435             : 
     436             : INLINE GEN
     437      292091 : F2m_F2c_mul_i(GEN x, GEN y, long lx, long l)
     438             : {
     439             :   long j;
     440      292091 :   GEN z = NULL;
     441             : 
     442     3445271 :   for (j=1; j<lx; j++)
     443             :   {
     444     3153180 :     if (!F2v_coeff(y,j)) continue;
     445      706630 :     if (!z) z = vecsmall_copy(gel(x,j));
     446      435035 :     else F2v_add_inplace(z,gel(x,j));
     447             :   }
     448      292091 :   if (!z) z = zero_F2v(l);
     449      292091 :   return z;
     450             : }
     451             : 
     452             : GEN
     453      376632 : FpM_mul(GEN x, GEN y, GEN p)
     454             : {
     455      376632 :   long j, l, lx=lg(x), ly=lg(y);
     456             :   GEN z;
     457      376632 :   if (ly==1) return cgetg(1,t_MAT);
     458      376632 :   if (lx==1) return zeromat(0, ly-1);
     459      376632 :   if (lgefint(p) == 3)
     460             :   {
     461      372966 :     pari_sp av = avma;
     462      372966 :     ulong pp = uel(p,2);
     463      372966 :     if (pp == 2)
     464             :     {
     465       61058 :       x = ZM_to_F2m(x);
     466       61058 :       y = ZM_to_F2m(y);
     467       61058 :       z = F2m_to_ZM(F2m_mul(x,y));
     468             :     }
     469             :     else
     470             :     {
     471      311908 :       x = ZM_to_Flm(x, pp);
     472      311908 :       y = ZM_to_Flm(y, pp);
     473      311908 :       z = Flm_to_ZM(Flm_mul(x,y, pp));
     474             :     }
     475      372966 :     return gerepileupto(av, z);
     476             :   }
     477        3666 :   l = lgcols(x); z = cgetg(ly,t_MAT);
     478        3666 :   for (j=1; j<ly; j++) gel(z,j) = FpM_FpC_mul_i(x, gel(y,j), lx, l, p);
     479        3666 :   return z;
     480             : }
     481             : 
     482             : static GEN
     483     6363966 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     484             : {
     485             :   long j;
     486     6363966 :   GEN z = cgetg(ly, t_MAT);
     487     6363661 :   if (SMALL_ULONG(p))
     488    23041784 :     for (j=1; j<ly; j++)
     489    20188920 :       gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
     490             :   else
     491    36800734 :     for (j=1; j<ly; j++)
     492    33289693 :       gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
     493     6363905 :   return z;
     494             : }
     495             : 
     496             : /*
     497             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     498             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     499             : */
     500             : static void
     501        7480 : add_slices_ip(long m, long n,
     502             :            GEN A, long ma, long da, long na, long ea,
     503             :            GEN B, long mb, long db, long nb, long eb,
     504             :            GEN M, long dy, long dx, ulong p)
     505             : {
     506        7480 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     507             :   GEN C;
     508             : 
     509      284155 :   for (j = 1; j <= min_e; j++) {
     510      276675 :     C = gel(M, j + dx) + dy;
     511    12840349 :     for (i = 1; i <= min_d; i++)
     512    25127348 :       uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
     513    12563674 :                         ucoeff(B, mb + i, nb + j), p);
     514      283850 :     for (; i <= da; i++)
     515        7175 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     516      276675 :     for (; i <= db; i++)
     517           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     518      276675 :     for (; i <= m; i++)
     519           0 :       uel(C, i) = 0;
     520             :   }
     521        8086 :   for (; j <= ea; j++) {
     522         606 :     C = gel(M, j + dx) + dy;
     523       29252 :     for (i = 1; i <= da; i++)
     524       28646 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     525         606 :     for (; i <= m; i++)
     526           0 :       uel(C, i) = 0;
     527             :   }
     528        7480 :   for (; j <= eb; j++) {
     529           0 :     C = gel(M, j + dx) + dy;
     530           0 :     for (i = 1; i <= db; i++)
     531           0 :       uel(C, i) = ucoeff(B, mb + i, nb + j);
     532           0 :     for (; i <= m; i++)
     533           0 :       uel(C, i) = 0;
     534             :   }
     535        7480 :   for (; j <= n; j++) {
     536           0 :     C = gel(M, j + dx) + dy;
     537           0 :     for (i = 1; i <= m; i++)
     538           0 :       uel(C, i) = 0;
     539             :   }
     540        7480 : }
     541             : 
     542             : static GEN
     543        3740 : add_slices(long m, long n,
     544             :            GEN A, long ma, long da, long na, long ea,
     545             :            GEN B, long mb, long db, long nb, long eb, ulong p)
     546             : {
     547             :   GEN M;
     548             :   long j;
     549        3740 :   M = cgetg(n + 1, t_MAT);
     550      140715 :   for (j = 1; j <= n; j++)
     551      136975 :     gel(M, j) = cgetg(m + 1, t_VECSMALL);
     552        3740 :   add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
     553        3740 :   return M;
     554             : }
     555             : 
     556             : /*
     557             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     558             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     559             : */
     560             : static GEN
     561        6545 : subtract_slices(long m, long n,
     562             :                 GEN A, long ma, long da, long na, long ea,
     563             :                 GEN B, long mb, long db, long nb, long eb, ulong p)
     564             : {
     565        6545 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     566        6545 :   GEN M = cgetg(n + 1, t_MAT), C;
     567             : 
     568      241255 :   for (j = 1; j <= min_e; j++) {
     569      234710 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     570    10503385 :     for (i = 1; i <= min_d; i++)
     571    20537350 :       uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
     572    10268675 :                         ucoeff(B, mb + i, nb + j), p);
     573      254171 :     for (; i <= da; i++)
     574       19461 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     575      241460 :     for (; i <= db; i++)
     576        6750 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     577      234710 :     for (; i <= m; i++)
     578           0 :       uel(C, i) = 0;
     579             :   }
     580        6545 :   for (; j <= ea; j++) {
     581           0 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     582           0 :     for (i = 1; i <= da; i++)
     583           0 :       uel(C, i) = ucoeff(A, ma + i, na + j);
     584           0 :     for (; i <= m; i++)
     585           0 :       uel(C, i) = 0;
     586             :   }
     587        6696 :   for (; j <= eb; j++) {
     588         151 :     gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
     589        6505 :     for (i = 1; i <= db; i++)
     590        6354 :       uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
     591         151 :     for (; i <= m; i++)
     592           0 :       uel(C, i) = 0;
     593             :   }
     594        6696 :   for (; j <= n; j++)
     595         151 :     gel(M, j) = zero_Flv(m);
     596        6545 :   return M;
     597             : }
     598             : 
     599             : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
     600             : 
     601             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     602             : static GEN
     603         935 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
     604             : {
     605             :   pari_sp av;
     606         935 :   long m1 = (m + 1)/2, m2 = m/2,
     607         935 :     n1 = (n + 1)/2, n2 = n/2,
     608         935 :     p1 = (p + 1)/2, p2 = p/2;
     609             :   GEN A11, A12, A22, B11, B21, B22,
     610             :     S1, S2, S3, S4, T1, T2, T3, T4,
     611             :     M1, M2, M3, M4, M5, M6, M7,
     612             :     V1, V2, V3, C;
     613             :   long j;
     614         935 :   C = cgetg(p + 1, t_MAT);
     615       71088 :   for (j = 1; j <= p; j++)
     616       70153 :     gel(C, j) = cgetg(m + 1, t_VECSMALL);
     617         935 :   av = avma;
     618         935 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
     619         935 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
     620         935 :   M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
     621         935 :   if (gc_needed(av, 1))
     622           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     623         935 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
     624         935 :   if (gc_needed(av, 1))
     625           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     626         935 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
     627         935 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
     628         935 :   M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
     629         935 :   if (gc_needed(av, 1))
     630           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     631         935 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
     632         935 :   if (gc_needed(av, 1))
     633           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     634         935 :   A11 = matslice(A, 1, m1, 1, n1);
     635         935 :   B11 = matslice(B, 1, n1, 1, p1);
     636         935 :   M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
     637         935 :   if (gc_needed(av, 1))
     638           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     639         935 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     640         935 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     641         935 :   M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
     642         935 :   if (gc_needed(av, 1))
     643           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     644         935 :   add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
     645         935 :   if (gc_needed(av, 1))
     646           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy M4 */
     647         935 :   M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
     648         935 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
     649         935 :   if (gc_needed(av, 1))
     650           0 :     gerepileall(av, 6, &M2, &T3, &M3, &M1, &M5, &S4);  /* destroy S3 */
     651         935 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
     652         935 :   if (gc_needed(av, 1))
     653           0 :     gerepileall(av, 6, &M2, &M3, &M1, &M5, &S4, &T4);  /* destroy T3 */
     654         935 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
     655         935 :   if (gc_needed(av, 1))
     656           0 :     gerepileall(av, 5, &M2, &M3, &S4, &T4, &V1);  /* destroy M1, M5 */
     657         935 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     658         935 :   M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
     659         935 :   if (gc_needed(av, 1))
     660           0 :     gerepileall(av, 5, &M2, &M3, &T4, &V1, &M6);  /* destroy S4, B22 */
     661         935 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     662         935 :   M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
     663         935 :   if (gc_needed(av, 1))
     664           0 :     gerepileall(av, 5, &M2, &M3, &V1, &M6, &M7);  /* destroy A22, T4 */
     665         935 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
     666         935 :   add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
     667         935 :   if (gc_needed(av, 1))
     668           6 :     gerepileall(av, 4, &M2, &M3, &V1, &M7);  /* destroy V3, M6 */
     669         935 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
     670         935 :   if (gc_needed(av, 1))
     671           0 :     gerepileall(av, 3, &M3, &M7, &V2);  /* destroy V1, M2 */
     672         935 :   add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
     673         935 :   add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
     674         935 :   avma = av; return C;
     675             : }
     676             : 
     677             : /* Strassen-Winograd used for dim >= ZM_sw_bound */
     678             : static GEN
     679     6364928 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
     680             : {
     681     6364928 :   ulong e = expu(p);
     682             : #ifdef LONG_IS_64BIT
     683     5011749 :   long ZM_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
     684             : #else
     685     1353210 :   long ZM_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
     686             : #endif
     687     6364959 :   if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
     688     6364024 :     return Flm_mul_classical(x, y, l, lx, ly, p, pi);
     689             :   else
     690         935 :     return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
     691             : }
     692             : 
     693             : GEN
     694     6358694 : Flm_mul(GEN x, GEN y, ulong p)
     695             : {
     696     6358694 :   long lx=lg(x), ly=lg(y);
     697     6358694 :   if (ly==1) return cgetg(1,t_MAT);
     698     6358589 :   if (lx==1) return zero_Flm(0, ly-1);
     699     6358589 :   return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
     700             : }
     701             : 
     702             : GEN
     703       61065 : F2m_mul(GEN x, GEN y)
     704             : {
     705       61065 :   long i,j,l,lx=lg(x), ly=lg(y);
     706             :   GEN z;
     707       61065 :   if (ly==1) return cgetg(1,t_MAT);
     708       61065 :   z = cgetg(ly,t_MAT);
     709       61065 :   if (lx==1)
     710             :   {
     711           0 :     for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
     712           0 :     return z;
     713             :   }
     714       61065 :   l = coeff(x,1,1);
     715       61065 :   for (j=1; j<ly; j++) gel(z,j) = F2m_F2c_mul_i(x, gel(y,j), lx, l);
     716       61065 :   return z;
     717             : }
     718             : 
     719             : struct _Flm
     720             : {
     721             :   ulong p;
     722             :   long n;
     723             : };
     724             : 
     725             : static GEN
     726        8099 : _Flm_mul(void *E , GEN x, GEN y)
     727        8099 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
     728             : static GEN
     729       27279 : _Flm_sqr(void *E, GEN x)
     730       27279 : { return Flm_mul(x,x,((struct _Flm*)E)->p); }
     731             : static GEN
     732         126 : _Flm_one(void *E)
     733         126 : { return matid_Flm(((struct _Flm*)E)->n); }
     734             : GEN
     735       14658 : Flm_powu(GEN x, ulong n, ulong p)
     736             : {
     737       14658 :   pari_sp av = avma;
     738             :   struct _Flm d;
     739       14658 :   if (!n) return matid(lg(x)-1);
     740       14658 :   d.p = p;
     741       14658 :   return gerepileupto(av, gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul));
     742             : }
     743             : GEN
     744         126 : Flm_powers(GEN x, ulong n, ulong p)
     745             : {
     746         126 :   pari_sp av = avma;
     747             :   struct _Flm d;
     748         126 :   d.p = p;
     749         126 :   d.n = lg(x)-1;
     750         126 :   return gerepileupto(av, gen_powers(x, n, 1, (void*)&d,
     751             :                           &_Flm_sqr, &_Flm_mul, &_Flm_one));
     752             : }
     753             : 
     754             : static GEN
     755           0 : _F2m_mul(void *data, GEN x, GEN y)
     756           0 : { (void) data; return F2m_mul(x,y); }
     757             : static GEN
     758           0 : _F2m_sqr(void *data, GEN x)
     759           0 : { (void) data; return F2m_mul(x,x); }
     760             : GEN
     761           0 : F2m_powu(GEN x, ulong n)
     762             : {
     763           0 :   pari_sp av = avma;
     764           0 :   if (!n) return matid(lg(x)-1);
     765           0 :   return gerepileupto(av, gen_powu(x, n,NULL, &_F2m_sqr, &_F2m_mul));
     766             : }
     767             : static GEN
     768           0 : _FpM_mul(void *p , GEN x, GEN y)
     769           0 : { return FpM_mul(x,y,(GEN)p); }
     770             : static GEN
     771           0 : _FpM_sqr(void *p, GEN x)
     772           0 : { return FpM_mul(x,x,(GEN)p); }
     773             : GEN
     774           0 : FpM_powu(GEN x, ulong n, GEN p)
     775             : {
     776           0 :   pari_sp av = avma;
     777           0 :   if (!n) return matid(lg(x)-1);
     778           0 :   if (lgefint(p) == 3)
     779             :   {
     780           0 :     pari_sp av = avma;
     781           0 :     ulong pp = uel(p,2);
     782             :     GEN z;
     783           0 :     if (pp == 2)
     784           0 :       z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
     785             :     else
     786           0 :       z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
     787           0 :     return gerepileupto(av, z);
     788             :   }
     789           0 :   return gerepileupto(av, gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul));
     790             : }
     791             : 
     792             : /*Multiple a column vector by a line vector to make a matrix*/
     793             : GEN
     794           7 : Flc_Flv_mul(GEN x, GEN y, ulong p)
     795             : {
     796           7 :   long i,j, lx=lg(x), ly=lg(y);
     797             :   GEN z;
     798           7 :   if (ly==1) return cgetg(1,t_MAT);
     799           7 :   z = cgetg(ly, t_MAT);
     800          28 :   for (j=1; j < ly; j++)
     801             :   {
     802          21 :     GEN zj = cgetg(lx,t_VECSMALL);
     803          84 :     for (i=1; i<lx; i++)
     804          63 :       uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
     805          21 :     gel(z,j) = zj;
     806             :   }
     807           7 :   return z;
     808             : }
     809             : 
     810             : /*Multiple a column vector by a line vector to make a matrix*/
     811             : GEN
     812        1540 : FpC_FpV_mul(GEN x, GEN y, GEN p)
     813             : {
     814        1540 :   long i,j, lx=lg(x), ly=lg(y);
     815             :   GEN z;
     816        1540 :   if (ly==1) return cgetg(1,t_MAT);
     817        1540 :   z = cgetg(ly,t_MAT);
     818       28252 :   for (j=1; j < ly; j++)
     819             :   {
     820       26712 :     GEN zj = cgetg(lx,t_COL);
     821       26712 :     for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
     822       26712 :     gel(z, j) = zj;
     823             :   }
     824        1540 :   return z;
     825             : }
     826             : 
     827             : /* Multiply a line vector by a column and return a scalar (t_INT) */
     828             : GEN
     829     1315363 : FpV_dotproduct(GEN x, GEN y, GEN p)
     830             : {
     831     1315363 :   long i, lx = lg(x);
     832             :   pari_sp av;
     833             :   GEN c;
     834     1315363 :   if (lx == 1) return gen_0;
     835     1313004 :   av = avma; c = mulii(gel(x,1),gel(y,1));
     836     1313004 :   for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
     837     1313004 :   return gerepileuptoint(av, modii(c,p));
     838             : }
     839             : GEN
     840           0 : FpV_dotsquare(GEN x, GEN p)
     841             : {
     842           0 :   long i, lx = lg(x);
     843             :   pari_sp av;
     844             :   GEN c;
     845           0 :   if (lx == 1) return gen_0;
     846           0 :   av = avma; c = sqri(gel(x,1));
     847           0 :   for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
     848           0 :   return gerepileuptoint(av, modii(c,p));
     849             : }
     850             : 
     851             : INLINE ulong
     852      692868 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
     853             : {
     854      692868 :   ulong c = uel(x,0) * uel(y,0);
     855             :   long k;
     856     6912991 :   for (k = 1; k < lx; k++) {
     857     6220123 :     c += uel(x,k) * uel(y,k);
     858     6220123 :     if (c & HIGHBIT) c %= p;
     859             :   }
     860      692868 :   return c % p;
     861             : }
     862             : 
     863             : INLINE ulong
     864        5542 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
     865             : {
     866             :   ulong l0, l1, v1, h0, h1;
     867        5542 :   long i = 0;
     868             :   LOCAL_OVERFLOW;
     869             :   LOCAL_HIREMAINDER;
     870        5542 :   l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
     871      166552 :   while (++i < lx) {
     872      155468 :     l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
     873      155468 :     l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
     874             :   }
     875        5542 :   if (v1 == 0) return remll_pre(h1, l1, p, pi);
     876           0 :   else return remlll_pre(v1, h1, l1, p, pi);
     877             : }
     878             : 
     879             : ulong
     880      254772 : Flv_dotproduct(GEN x, GEN y, ulong p)
     881             : {
     882      254772 :   long lx = lg(x)-1;
     883      254772 :   if (lx == 0) return 0;
     884      254772 :   if (SMALL_ULONG(p))
     885      254772 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     886             :   else
     887           0 :     return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
     888             : }
     889             : 
     890             : ulong
     891       53460 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
     892             : {
     893       53460 :   long lx = lg(x)-1;
     894       53460 :   if (lx == 0) return 0;
     895       53460 :   if (SMALL_ULONG(p))
     896       48002 :     return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
     897             :   else
     898        5458 :     return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
     899             : }
     900             : 
     901             : ulong
     902      400694 : Flx_dotproduct(GEN x, GEN y, ulong p)
     903             : {
     904      400694 :   long lx = minss(lgpol(x), lgpol(y));
     905      400694 :   if (lx == 0) return 0;
     906      390178 :   if (SMALL_ULONG(p))
     907      390094 :     return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
     908             :   else
     909          84 :     return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
     910             : }
     911             : 
     912             : ulong
     913           0 : F2v_dotproduct(GEN x, GEN y)
     914             : {
     915           0 :   long i, lx = lg(x);
     916             :   ulong c;
     917           0 :   if (lx == 1) return 0;
     918           0 :   c = uel(x,1) & uel(y,1);
     919           0 :   for (i=2; i<lx; i++) c ^= uel(x,i) & uel(y,i);
     920             : #ifdef LONG_IS_64BIT
     921           0 :   c ^= c >> 32;
     922             : #endif
     923           0 :   c ^= c >> 16;
     924           0 :   c ^= c >>  8;
     925           0 :   c ^= c >>  4;
     926           0 :   c ^= c >>  2;
     927           0 :   c ^= c >>  1;
     928           0 :   return c & 1;
     929             : }
     930             : 
     931             : GEN
     932      895831 : FpM_FpC_mul(GEN x, GEN y, GEN p)
     933             : {
     934      895831 :   long lx = lg(x);
     935      895831 :   return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
     936             : }
     937             : GEN
     938      208664 : Flm_Flc_mul(GEN x, GEN y, ulong p)
     939             : {
     940      208664 :   long l, lx = lg(x);
     941      208664 :   if (lx==1) return cgetg(1,t_VECSMALL);
     942      208664 :   l = lgcols(x);
     943      208664 :   if (p==2)
     944       82593 :     return Flm_Flc_mul_i_2(x, y, lx, l);
     945      126071 :   else if (SMALL_ULONG(p))
     946      125957 :     return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
     947             :   else
     948         114 :     return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
     949             : }
     950             : 
     951             : GEN
     952        4816 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     953             : {
     954        4816 :   long l, lx = lg(x);
     955             :   GEN z;
     956        4816 :   if (lx==1) return cgetg(1,t_VECSMALL);
     957        4816 :   l = lgcols(x);
     958        4815 :   z = cgetg(l, t_VECSMALL);
     959        4815 :   if (SMALL_ULONG(p))
     960        4127 :     __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
     961             :   else
     962         688 :     __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
     963        4816 :   return z;
     964             : }
     965             : 
     966             : GEN
     967     7304676 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
     968             : {
     969     7304676 :   long l, lx = lg(x);
     970             :   GEN z;
     971     7304676 :   if (lx==1) return pol0_Flx(sv);
     972     7304676 :   l = lgcols(x);
     973     7310229 :   z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
     974     7313949 :   if (SMALL_ULONG(p))
     975     3129389 :     __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
     976             :   else
     977     4184560 :     __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
     978     7318031 :   return Flx_renormalize(z, l + 1);
     979             : }
     980             : 
     981             : GEN
     982           0 : F2m_F2c_mul(GEN x, GEN y)
     983             : {
     984           0 :   long l, lx = lg(x);
     985           0 :   if (lx==1) return cgetg(1,t_VECSMALL);
     986           0 :   l = coeff(x,1,1);
     987           0 :   return F2m_F2c_mul_i(x, y, lx, l);
     988             : }
     989             : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
     990             : GEN
     991      325048 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
     992             : {
     993      325048 :   long i, l, lx = lg(x);
     994             :   GEN z;
     995      325048 :   if (lx==1) return pol_0(v);
     996      325048 :   l = lgcols(x);
     997      325048 :   z = new_chunk(l+1);
     998      709640 :   for (i=l-1; i; i--)
     999             :   {
    1000      561161 :     pari_sp av = avma;
    1001      561161 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
    1002      561161 :     p1 = modii(p1, p);
    1003      561161 :     if (signe(p1))
    1004             :     {
    1005      176569 :       if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
    1006      176569 :       gel(z,i+1) = gerepileuptoint(av, p1);
    1007      176569 :       break;
    1008             :     }
    1009      384592 :     avma = av;
    1010             :   }
    1011      325048 :   if (!i) { avma = (pari_sp)(z + l+1); return pol_0(v); }
    1012      176569 :   z[0] = evaltyp(t_POL) | evallg(i+2);
    1013      176569 :   z[1] = evalsigne(1) | evalvarn(v);
    1014      545827 :   for (; i; i--)
    1015             :   {
    1016      369258 :     pari_sp av = avma;
    1017      369258 :     GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
    1018      369258 :     gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
    1019             :   }
    1020      176569 :   return z;
    1021             : }
    1022             : 
    1023             : /********************************************************************/
    1024             : /**                                                                **/
    1025             : /**                           TRANSPOSITION                        **/
    1026             : /**                                                                **/
    1027             : /********************************************************************/
    1028             : 
    1029             : /* == zm_transpose */
    1030             : GEN
    1031      365718 : Flm_transpose(GEN x)
    1032             : {
    1033      365718 :   long i, dx, lx = lg(x);
    1034             :   GEN y;
    1035      365718 :   if (lx == 1) return cgetg(1,t_MAT);
    1036      365613 :   dx = lgcols(x); y = cgetg(dx,t_MAT);
    1037      365617 :   for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
    1038      365615 :   return y;
    1039             : }
    1040             : 
    1041             : /********************************************************************/
    1042             : /**                                                                **/
    1043             : /**                           SCALAR MATRICES                      **/
    1044             : /**                                                                **/
    1045             : /********************************************************************/
    1046             : 
    1047             : GEN
    1048        1897 : gen_matid(long n, void *E, const struct bb_field *S)
    1049             : {
    1050        1897 :   GEN y = cgetg(n+1,t_MAT), _0, _1;
    1051             :   long i;
    1052        1897 :   if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
    1053        1897 :   _0 = S->s(E,0);
    1054        1897 :   _1 = S->s(E,1);
    1055        8505 :   for (i=1; i<=n; i++)
    1056             :   {
    1057        6608 :     GEN z = const_col(n, _0); gel(z,i) = _1;
    1058        6608 :     gel(y, i) = z;
    1059             :   }
    1060        1897 :   return y;
    1061             : }
    1062             : 
    1063             : GEN
    1064       12866 : matid_F2m(long n)
    1065             : {
    1066       12866 :   GEN y = cgetg(n+1,t_MAT);
    1067             :   long i;
    1068       12866 :   if (n < 0) pari_err_DOMAIN("matid_F2m", "dimension","<",gen_0,stoi(n));
    1069       12866 :   for (i=1; i<=n; i++) { gel(y,i) = zero_F2v(n); F2v_set(gel(y,i),i); }
    1070       12866 :   return y;
    1071             : }
    1072             : 
    1073             : GEN
    1074          35 : matid_F2xqM(long n, GEN T)
    1075             : {
    1076             :   void *E;
    1077          35 :   const struct bb_field *S = get_F2xq_field(&E, T);
    1078          35 :   return gen_matid(n, E, S);
    1079             : }
    1080             : GEN
    1081        1862 : matid_FlxqM(long n, GEN T, ulong p)
    1082             : {
    1083             :   void *E;
    1084        1862 :   const struct bb_field *S = get_Flxq_field(&E, T, p);
    1085        1862 :   return gen_matid(n, E, S);
    1086             : }
    1087             : 
    1088             : GEN
    1089      879800 : matid_Flm(long n)
    1090             : {
    1091      879800 :   GEN y = cgetg(n+1,t_MAT);
    1092             :   long i;
    1093      879794 :   if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
    1094      879803 :   for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
    1095      879804 :   return y;
    1096             : }
    1097             : 
    1098             : GEN
    1099          42 : scalar_Flm(long s, long n)
    1100             : {
    1101             :   long i;
    1102          42 :   GEN y = cgetg(n+1,t_MAT);
    1103          42 :   for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
    1104          42 :   return y;
    1105             : }
    1106             : 
    1107             : /********************************************************************/
    1108             : /**                                                                **/
    1109             : /**                           CONVERSIONS                          **/
    1110             : /**                                                                **/
    1111             : /********************************************************************/
    1112             : GEN
    1113    17074645 : ZV_to_Flv(GEN x, ulong p)
    1114    17074645 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
    1115             : 
    1116             : GEN
    1117     3512388 : ZM_to_Flm(GEN x, ulong p)
    1118     3512388 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
    1119             : 
    1120             : GEN
    1121        1302 : ZMV_to_FlmV(GEN x, ulong m)
    1122        1302 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
    1123             : 
    1124             : /*                          TO INTMOD                        */
    1125             : static GEN
    1126     9299867 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
    1127             : static GEN
    1128       19348 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
    1129             : 
    1130             : GEN
    1131      289172 : Fp_to_mod(GEN z, GEN p)
    1132             : {
    1133      289172 :   retmkintmod(modii(z, p), icopy(p));
    1134             : }
    1135             : 
    1136             : /* z in Z[X], return z * Mod(1,p), normalized*/
    1137             : GEN
    1138     1121629 : FpX_to_mod(GEN z, GEN p)
    1139             : {
    1140     1121629 :   long i,l = lg(z);
    1141     1121629 :   GEN x = cgetg(l,t_POL);
    1142     1121629 :   if (l >2) p = icopy(p);
    1143     1121629 :   for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1144     1121629 :   x[1] = z[1]; return normalizepol_lg(x,l);
    1145             : }
    1146             : 
    1147             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1148             : GEN
    1149        4648 : FpV_to_mod(GEN z, GEN p)
    1150             : {
    1151        4648 :   long i,l = lg(z);
    1152        4648 :   GEN x = cgetg(l, t_VEC);
    1153        4648 :   if (l == 1) return x;
    1154        4648 :   p = icopy(p);
    1155        4648 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1156        4648 :   return x;
    1157             : }
    1158             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1159             : GEN
    1160          63 : FpC_to_mod(GEN z, GEN p)
    1161             : {
    1162          63 :   long i, l = lg(z);
    1163          63 :   GEN x = cgetg(l, t_COL);
    1164          63 :   if (l == 1) return x;
    1165          63 :   p = icopy(p);
    1166          63 :   for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
    1167          63 :   return x;
    1168             : }
    1169             : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
    1170             : GEN
    1171       56385 : FpM_to_mod(GEN z, GEN p)
    1172             : {
    1173       56385 :   long i, j, m, l = lg(z);
    1174       56385 :   GEN  x = cgetg(l,t_MAT), y, zi;
    1175       56385 :   if (l == 1) return x;
    1176       56364 :   m = lgcols(z);
    1177       56364 :   p = icopy(p);
    1178      181965 :   for (i=1; i<l; i++)
    1179             :   {
    1180      125601 :     gel(x,i) = cgetg(m,t_COL);
    1181      125601 :     y = gel(x,i); zi= gel(z,i);
    1182      125601 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1183             :   }
    1184       56364 :   return x;
    1185             : }
    1186             : GEN
    1187          28 : Flc_to_mod(GEN z, ulong pp)
    1188             : {
    1189          28 :   long i, l = lg(z);
    1190          28 :   GEN p, x = cgetg(l, t_COL);
    1191          28 :   if (l == 1) return x;
    1192          28 :   p = utoipos(pp);
    1193          28 :   for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
    1194          28 :   return x;
    1195             : }
    1196             : GEN
    1197         182 : Flm_to_mod(GEN z, ulong pp)
    1198             : {
    1199         182 :   long i, j, m, l = lg(z);
    1200         182 :   GEN p, x = cgetg(l,t_MAT), y, zi;
    1201         182 :   if (l == 1) return x;
    1202         154 :   m = lgcols(z);
    1203         154 :   p = utoipos(pp);
    1204        1393 :   for (i=1; i<l; i++)
    1205             :   {
    1206        1239 :     gel(x,i) = cgetg(m,t_COL);
    1207        1239 :     y = gel(x,i); zi= gel(z,i);
    1208        1239 :     for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
    1209             :   }
    1210         154 :   return x;
    1211             : }
    1212             : 
    1213             : GEN
    1214        1827 : FpVV_to_mod(GEN z, GEN p)
    1215             : {
    1216        1827 :   long i, j, m, l = lg(z);
    1217        1827 :   GEN  x = cgetg(l,t_VEC), y, zi;
    1218        1827 :   if (l == 1) return x;
    1219        1750 :   m = lgcols(z);
    1220        1750 :   p = icopy(p);
    1221        3661 :   for (i=1; i<l; i++)
    1222             :   {
    1223        1911 :     gel(x,i) = cgetg(m,t_VEC);
    1224        1911 :     y = gel(x,i); zi= gel(z,i);
    1225        1911 :     for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
    1226             :   }
    1227        1750 :   return x;
    1228             : }
    1229             : 
    1230             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1231             : GEN
    1232           7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
    1233             : {
    1234           7 :   long i,l = lg(z);
    1235           7 :   GEN x = cgetg(l, t_COL);
    1236           7 :   if (l == 1) return x;
    1237           7 :   T = FpX_to_mod(T, p);
    1238          21 :   for (i=1; i<l; i++)
    1239          14 :     gel(x,i) = mkpolmod(FpX_to_mod(gel(z,i), p), T);
    1240           7 :   return x;
    1241             : }
    1242             : 
    1243             : static GEN
    1244       61208 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
    1245             : {
    1246       61208 :   GEN z = typ(x)==t_INT ? Fp_to_mod(x, p): FpX_to_mod(x, p);
    1247       61208 :   return mkpolmod(z, T);
    1248             : }
    1249             : 
    1250             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1251             : static GEN
    1252        2660 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
    1253        2660 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
    1254             : 
    1255             : /* z in Z^n, return z * Mod(1,p), normalized*/
    1256             : GEN
    1257         196 : FqM_to_mod(GEN z, GEN T, GEN p)
    1258             : {
    1259             :   GEN x;
    1260         196 :   long i,l = lg(z);
    1261         196 :   if (!T) return FpM_to_mod(z, p);
    1262         196 :   x = cgetg(l, t_MAT);
    1263         196 :   if (l == 1) return x;
    1264         168 :   T = FpX_to_mod(T, p);
    1265        2828 :   for (i=1; i<l; i++)
    1266        2660 :     gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
    1267         168 :   return x;
    1268             : }
    1269             : 
    1270             : /********************************************************************/
    1271             : /*                                                                  */
    1272             : /*                     Blackbox linear algebra                      */
    1273             : /*                                                                  */
    1274             : /********************************************************************/
    1275             : 
    1276             : /* A sparse column (zCs) v is a t_COL with two components C and E which are
    1277             :  * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
    1278             :  * (e_j) is the canonical basis.
    1279             :  * A sparse matrix (zMs) is a t_VEC of zCs */
    1280             : 
    1281             : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
    1282             :  * integer representing an element of Fp. This is important since p can be
    1283             :  * large and p+E[i] would not fit in a C long.  */
    1284             : 
    1285             : /* RgCs and RgMs are similar, except that the type of the component is
    1286             :  * unspecified. Functions handling RgCs/RgMs must be independent of the type
    1287             :  * of E. */
    1288             : 
    1289             : /* Most functions take an argument nbrow which is the number of lines of the
    1290             :  * column/matrix, which cannot be derived from the data. */
    1291             : 
    1292             : GEN
    1293           0 : zCs_to_ZC(GEN R, long nbrow)
    1294             : {
    1295           0 :   GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
    1296           0 :   long j, l = lg(C);
    1297           0 :   for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
    1298           0 :   return c;
    1299             : }
    1300             : 
    1301             : GEN
    1302           0 : zMs_to_ZM(GEN x, long nbrow)
    1303           0 : { pari_APPLY_same(zCs_to_ZC(gel(x, i), nbrow)) }
    1304             : 
    1305             : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
    1306             :  * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
    1307             : GEN
    1308         119 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
    1309             : {
    1310         119 :   pari_sp ltop = avma;
    1311         119 :   long col = 0, n = lg(B)-1, m = 2*n+1;
    1312         119 :   if (ZV_equal0(B)) return zerocol(n);
    1313         238 :   while (++col <= n)
    1314             :   {
    1315         119 :     pari_sp btop = avma, av;
    1316             :     long i, lQ;
    1317         119 :     GEN V, Q, M, W = B;
    1318         119 :     GEN b = cgetg(m+2, t_POL);
    1319         119 :     b[1] = evalsigne(1)|evalvarn(0);
    1320         119 :     gel(b, 2) = gel(W, col);
    1321       50567 :     for (i = 3; i<m+2; i++)
    1322       50448 :       gel(b, i) = cgeti(lgefint(p));
    1323         119 :     av = avma;
    1324       50567 :     for (i = 3; i<m+2; i++)
    1325             :     {
    1326       50448 :       W = f(E, W);
    1327       50448 :       affii(gel(W, col),gel(b, i));
    1328       50448 :       if (gc_needed(av,1))
    1329             :       {
    1330        2000 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
    1331        2000 :         W = gerepileupto(av, W);
    1332             :       }
    1333             :     }
    1334         119 :     b = FpX_renormalize(b, m+2);
    1335         119 :     if (lgpol(b)==0) {avma = btop; continue; }
    1336         119 :     M = FpX_halfgcd(b, pol_xn(m, 0), p);
    1337         119 :     Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
    1338         119 :     W = B; lQ =lg(Q);
    1339         119 :     if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
    1340         119 :     V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
    1341         119 :     av = avma;
    1342       24811 :     for (i = lQ-3; i > 1; i--)
    1343             :     {
    1344       24692 :       W = f(E, W);
    1345       24692 :       V = ZC_lincomb(gen_1, gel(Q,i), V, W);
    1346       24692 :       if (gc_needed(av,1))
    1347             :       {
    1348        1583 :         if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
    1349        1583 :         gerepileall(av, 2, &V, &W);
    1350             :       }
    1351             :     }
    1352         119 :     V = FpC_red(V, p);
    1353         119 :     W = FpC_sub(f(E,V), B, p);
    1354         238 :     if (ZV_equal0(W)) return gerepilecopy(ltop, V);
    1355           0 :     av = avma;
    1356           0 :     for (i = 1; i <= n; ++i)
    1357             :     {
    1358           0 :       V = W;
    1359           0 :       W = f(E, V);
    1360           0 :       if (ZV_equal0(W))
    1361           0 :         return gerepilecopy(ltop, shallowtrans(V));
    1362           0 :       gerepileall(av, 2, &V, &W);
    1363             :     }
    1364           0 :     avma = btop;
    1365             :   }
    1366           0 :   return NULL;
    1367             : }
    1368             : 
    1369             : GEN
    1370           0 : zMs_ZC_mul(GEN M, GEN B)
    1371             : {
    1372             :   long i, j;
    1373           0 :   long n = lg(B)-1;
    1374           0 :   GEN V = zerocol(n);
    1375           0 :   for (i = 1; i <= n; ++i)
    1376           0 :     if (signe(gel(B, i)))
    1377             :     {
    1378           0 :       GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1379           0 :       long l = lg(C);
    1380           0 :       for (j = 1; j < l; ++j)
    1381             :       {
    1382           0 :         long k = C[j];
    1383           0 :         switch(E[j])
    1384             :         {
    1385             :         case 1:
    1386           0 :           gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
    1387           0 :           break;
    1388             :         case -1:
    1389           0 :           gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
    1390           0 :           break;
    1391             :         default:
    1392           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]));
    1393           0 :           break;
    1394             :         }
    1395             :       }
    1396             :     }
    1397           0 :   return V;
    1398             : }
    1399             : 
    1400             : GEN
    1401           0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
    1402             : 
    1403             : GEN
    1404       75378 : ZV_zMs_mul(GEN B, GEN M)
    1405             : {
    1406             :   long i, j;
    1407       75378 :   long m = lg(M)-1;
    1408       75378 :   GEN V = cgetg(m+1,t_VEC);
    1409    39262492 :   for (i = 1; i <= m; ++i)
    1410             :   {
    1411    39187114 :     GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
    1412    39187114 :     long l = lg(C);
    1413    39187114 :     GEN z = mulis(gel(B, C[1]), E[1]);
    1414   329595309 :     for (j = 2; j < l; ++j)
    1415             :     {
    1416   290408195 :       long k = C[j];
    1417   290408195 :       switch(E[j])
    1418             :       {
    1419             :       case 1:
    1420   196103761 :         z = addii(z, gel(B,k));
    1421   196103761 :         break;
    1422             :       case -1:
    1423    79106238 :         z = subii(z, gel(B,k));
    1424    79106238 :         break;
    1425             :       default:
    1426    15198196 :         z = addii(z, mulis(gel(B,k), E[j]));
    1427    15198196 :         break;
    1428             :       }
    1429             :     }
    1430    39187114 :     gel(V,i) = z;
    1431             :   }
    1432       75378 :   return V;
    1433             : }
    1434             : 
    1435             : GEN
    1436         119 : FpV_FpMs_mul(GEN B, GEN M, GEN p) { return FpV_red(ZV_zMs_mul(B, M), p); }
    1437             : 
    1438             : GEN
    1439      512611 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
    1440             : {
    1441      512611 :   pari_sp av = avma, av2;
    1442      512611 :   GEN xi, xb, pi = gen_1, P;
    1443             :   long i;
    1444      512611 :   if (!C) {
    1445           0 :     C = Flm_inv(ZM_to_Flm(a, p), p);
    1446           0 :     if (!C) pari_err_INV("ZlM_gauss", a);
    1447             :   }
    1448      512611 :   P = utoipos(p);
    1449      512611 :   av2 = avma;
    1450      512611 :   xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1451      512611 :   xb = Flm_to_ZM(xi);
    1452     1147449 :   for (i = 2; i <= e; i++)
    1453             :   {
    1454      634838 :     pi = muliu(pi, p); /* = p^(i-1) */
    1455      634838 :     b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
    1456      634838 :     if (gc_needed(av,2))
    1457             :     {
    1458          24 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
    1459          24 :       gerepileall(av2,3, &pi,&b,&xb);
    1460             :     }
    1461      634838 :     xi = Flm_mul(C, ZM_to_Flm(b, p), p);
    1462      634838 :     xb = ZM_add(xb, nm_Z_mul(xi, pi));
    1463             :   }
    1464      512611 :   return gerepileupto(av, xb);
    1465             : }
    1466             : 
    1467             : struct wrapper_modp_s {
    1468             :   GEN (*f)(void*, GEN);
    1469             :   void *E;
    1470             :   GEN p;
    1471             : };
    1472             : 
    1473             : /* compute f(x) mod p */
    1474             : static GEN
    1475       75259 : wrap_relcomb_modp(void *E, GEN x)
    1476             : {
    1477       75259 :   struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
    1478       75259 :   return FpC_red(W->f(W->E, x), W->p);
    1479             : }
    1480             : static GEN
    1481           0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
    1482             : 
    1483             : static GEN
    1484       75259 : wrap_relker(void*E, GEN x) { return ZV_zMs_mul(x, (GEN)E); }
    1485             : 
    1486             : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
    1487             : GEN
    1488           0 : gen_ZpM_Dixon(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
    1489             : {
    1490             :   struct wrapper_modp_s W;
    1491           0 :   pari_sp av = avma;
    1492           0 :   GEN xb, xi, pi = gen_1;
    1493             :   long i;
    1494           0 :   W.E = E;
    1495           0 :   W.f = f;
    1496           0 :   W.p = p;
    1497           0 :   xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
    1498           0 :   if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
    1499           0 :   xb = xi;
    1500           0 :   for (i = 2; i <= e; i++)
    1501             :   {
    1502           0 :     pi = mulii(pi, p); /* = p^(i-1) */
    1503           0 :     B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
    1504           0 :     if (gc_needed(av,2))
    1505             :     {
    1506           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon. i=%ld",i);
    1507           0 :       gerepileall(av,3, &pi,&B,&xb);
    1508             :     }
    1509           0 :     xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
    1510           0 :     if (!xi) return NULL;
    1511           0 :     if (typ(xi) == t_VEC) return gerepileupto(av, xi);
    1512           0 :     xb = ZC_add(xb, ZC_Z_mul(xi, pi));
    1513             :   }
    1514           0 :   return gerepileupto(av, xb);
    1515             : }
    1516             : 
    1517             : static GEN
    1518       25224 : vecprow(GEN A, GEN prow)
    1519             : {
    1520       25224 :   return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
    1521             : }
    1522             : 
    1523             : /* Solve the equation MX = A. Return either a solution as a t_COL,
    1524             :  * or the index of a column which is linearly dependent from the others as a
    1525             :  * t_VECSMALL with a single component. */
    1526             : GEN
    1527           0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
    1528             : {
    1529           0 :   pari_sp av = avma;
    1530             :   GEN pcol, prow;
    1531           0 :   long nbi=lg(M)-1, lR;
    1532             :   long i, n;
    1533             :   GEN Mp, Ap, Rp;
    1534             :   pari_timer ti;
    1535           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1536           0 :   RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
    1537           0 :   if (!pcol) { avma = av; return NULL; }
    1538           0 :   if (DEBUGLEVEL)
    1539           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
    1540           0 :   n = lg(pcol)-1;
    1541           0 :   Mp = cgetg(n+1, t_MAT);
    1542           0 :   for(i=1; i<=n; i++)
    1543           0 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1544           0 :   Ap = zCs_to_ZC(vecprow(A, prow), n);
    1545           0 :   if (DEBUGLEVEL) timer_start(&ti);
    1546           0 :   Rp = gen_ZpM_Dixon((void*)Mp,wrap_relcomb, Ap, p, e);
    1547           0 :   if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
    1548           0 :   if (!Rp) { avma = av; return NULL; }
    1549           0 :   lR = lg(Rp)-1;
    1550           0 :   if (typ(Rp) == t_COL)
    1551             :   {
    1552           0 :     GEN R = zerocol(nbi+1);
    1553           0 :     for (i=1; i<=lR; i++)
    1554           0 :        gel(R,pcol[i]) = gel(Rp,i);
    1555           0 :     return gerepilecopy(av, R);
    1556             :   }
    1557           0 :   for (i = 1; i <= lR; ++i)
    1558           0 :     if (signe(gel(Rp, i)))
    1559           0 :       return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
    1560             :   return NULL; /* LCOV_EXCL_LINE */
    1561             : }
    1562             : 
    1563             : GEN
    1564           0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
    1565             : {
    1566           0 :   return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1567             : }
    1568             : 
    1569             : GEN
    1570           0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
    1571             : {
    1572             :   GEN res;
    1573           0 :   pari_CATCH(e_INV)
    1574             :   {
    1575           0 :     GEN E = pari_err_last();
    1576           0 :     GEN x = err_get_compo(E,2);
    1577           0 :     if (typ(x) != t_INTMOD) pari_err(0,E);
    1578           0 :     if (DEBUGLEVEL)
    1579           0 :       pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
    1580           0 :     res = NULL;
    1581             :   } pari_TRY {
    1582           0 :     res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
    1583           0 :   } pari_ENDCATCH
    1584           0 :   return res;
    1585             : }
    1586             : 
    1587             : static GEN
    1588         119 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
    1589             : {
    1590         119 :   long i, j, oldf = 0, f = 0;
    1591         119 :   long lrow = lg(prow), lM = lg(M);
    1592         119 :   GEN W = const_vecsmall(lM-1,1);
    1593         119 :   GEN R = cgetg(lrow, t_VEC);
    1594      221956 :   for (i=1; i<lrow; i++)
    1595      221837 :     gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
    1596             :   do
    1597             :   {
    1598         386 :     oldf = f;
    1599      437057 :     for(i=1; i<lM; i++)
    1600      436671 :       if (W[i])
    1601             :       {
    1602      148516 :         GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
    1603      148516 :         long c=0, cj=0, lF = lg(F);
    1604     1300903 :         for(j=1; j<lF; j++)
    1605     1152387 :           if (!gel(R,F[j]))
    1606       86452 :           { c++; cj=j; }
    1607      148516 :         if (c>=2) continue;
    1608      127900 :         if (c==1)
    1609             :         {
    1610       40912 :           pari_sp av = avma;
    1611       40912 :           GEN s = gen_0;
    1612      356320 :           for(j=1; j<lF; j++)
    1613      315408 :             if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
    1614       40912 :           gel(R,F[cj]) = gerepileupto(av, Fp_div(Fp_neg(s, p), stoi(E[cj]), p));
    1615       40912 :           f++;
    1616             :         }
    1617      127900 :         W[i]=0;
    1618             :       }
    1619         386 :   } while(oldf!=f);
    1620      221956 :   for (i=1; i<lrow; i++)
    1621      221837 :     if (!gel(R,i)) gel(R,i) = gen_0;
    1622         119 :   return R;
    1623             : }
    1624             : 
    1625             : /* Return a linear form Y such that YM is essentially 0 */
    1626             : GEN
    1627         119 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
    1628             : {
    1629         119 :   pari_sp av = avma, av2;
    1630             :   GEN pcol, prow;
    1631             :   long i, n;
    1632             :   GEN Mp, B, MB, R, Rp;
    1633             :   pari_timer ti;
    1634             :   struct wrapper_modp_s W;
    1635         119 :   if (DEBUGLEVEL) timer_start(&ti);
    1636         119 :   RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
    1637         119 :   if (!pcol) { avma = av; return NULL; }
    1638         119 :   if (DEBUGLEVEL)
    1639           0 :     timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
    1640         119 :   n = lg(pcol)-1;
    1641         119 :   Mp = cgetg(n+1, t_MAT);
    1642       25343 :   for(i=1; i<=n; i++)
    1643       25224 :     gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
    1644         119 :   W.E = (void*) Mp;
    1645         119 :   W.f = wrap_relker;
    1646         119 :   W.p = p;
    1647         119 :   av2 = avma;
    1648             :   for(;;)
    1649             :   {
    1650         119 :     avma = av2;
    1651         119 :     B = random_FpV(n, p);
    1652         119 :     MB = FpV_FpMs_mul(B, Mp, p);
    1653         119 :     if (DEBUGLEVEL) timer_start(&ti);
    1654         119 :     pari_CATCH(e_INV)
    1655             :     {
    1656           0 :       GEN E = pari_err_last();
    1657           0 :       GEN x = err_get_compo(E,2);
    1658           0 :       if (typ(x) != t_INTMOD) pari_err(0,E);
    1659           0 :       if (DEBUGLEVEL)
    1660           0 :         pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
    1661           0 :       Rp = NULL;
    1662             :     } pari_TRY {
    1663         119 :       Rp = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, MB, p);
    1664         119 :     } pari_ENDCATCH
    1665         119 :     if (!Rp) continue;
    1666         119 :     if (typ(Rp)==t_COL)
    1667             :     {
    1668         119 :       Rp = FpC_sub(Rp,B,p);
    1669         119 :       if (ZV_equal0(Rp)) continue;
    1670             :     }
    1671         119 :     R = FpMs_structelim_back(M, Rp, prow, p);
    1672         119 :     if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
    1673         119 :     return gerepilecopy(av, R);
    1674           0 :   }
    1675             : }
    1676             : 
    1677             : GEN
    1678          35 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
    1679             : {
    1680          35 :   return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
    1681             : }

Generated by: LCOV version 1.11