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 - Flx.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.1 lcov report (development 28920-5f2f1e7126) Lines: 2402 2712 88.6 %
Date: 2023-12-09 07:54:41 Functions: 296 342 86.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2004  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             : /* Not so fast arithmetic with polynomials with small coefficients. */
      19             : 
      20             : static GEN
      21   898564430 : get_Flx_red(GEN T, GEN *B)
      22             : {
      23   898564430 :   if (typ(T)!=t_VEC) { *B=NULL; return T; }
      24      624650 :   *B = gel(T,1); return gel(T,2);
      25             : }
      26             : 
      27             : /***********************************************************************/
      28             : /**                              Flx                                  **/
      29             : /***********************************************************************/
      30             : /* Flx objects are defined as follows:
      31             :  * Let l an ulong. An Flx is a t_VECSMALL:
      32             :  * x[0] = codeword
      33             :  * x[1] = evalvarn(variable number)  (signe is not stored).
      34             :  * x[2] = a_0 x[3] = a_1, etc. with 0 <= a_i < l
      35             :  *
      36             :  * signe(x) is not valid. Use degpol(x)>0 instead. */
      37             : /***********************************************************************/
      38             : /**                      Conversion from Flx                          **/
      39             : /***********************************************************************/
      40             : 
      41             : GEN
      42    36153591 : Flx_to_ZX(GEN z)
      43             : {
      44    36153591 :   long i, l = lg(z);
      45    36153591 :   GEN x = cgetg(l,t_POL);
      46   238668424 :   for (i=2; i<l; i++) gel(x,i) = utoi(z[i]);
      47    36130933 :   x[1] = evalsigne(l-2!=0)| z[1]; return x;
      48             : }
      49             : 
      50             : GEN
      51       71332 : Flx_to_FlxX(GEN z, long sv)
      52             : {
      53       71332 :   long i, l = lg(z);
      54       71332 :   GEN x = cgetg(l,t_POL);
      55      277928 :   for (i=2; i<l; i++) gel(x,i) = Fl_to_Flx(z[i], sv);
      56       71333 :   x[1] = evalsigne(l-2!=0)| z[1]; return x;
      57             : }
      58             : 
      59             : /* same as Flx_to_ZX, in place */
      60             : GEN
      61    36344932 : Flx_to_ZX_inplace(GEN z)
      62             : {
      63    36344932 :   long i, l = lg(z);
      64   226917150 :   for (i=2; i<l; i++) gel(z,i) = utoi(z[i]);
      65    36304181 :   settyp(z, t_POL); z[1]=evalsigne(l-2!=0)|z[1]; return z;
      66             : }
      67             : 
      68             : /*Flx_to_Flv=zx_to_zv*/
      69             : GEN
      70    63073687 : Flx_to_Flv(GEN x, long N)
      71             : {
      72    63073687 :   GEN z = cgetg(N+1,t_VECSMALL);
      73    63061556 :   long i, l = lg(x)-1;
      74    63061556 :   x++;
      75   692392329 :   for (i=1; i<l ; i++) z[i]=x[i];
      76   323634608 :   for (   ; i<=N; i++) z[i]=0;
      77    63061556 :   return z;
      78             : }
      79             : 
      80             : /*Flv_to_Flx=zv_to_zx*/
      81             : GEN
      82    26786046 : Flv_to_Flx(GEN x, long sv)
      83             : {
      84    26786046 :   long i, l=lg(x)+1;
      85    26786046 :   GEN z = cgetg(l,t_VECSMALL); z[1]=sv;
      86    26775324 :   x--;
      87   291286389 :   for (i=2; i<l ; i++) z[i]=x[i];
      88    26775324 :   return Flx_renormalize(z,l);
      89             : }
      90             : 
      91             : /*Flm_to_FlxV=zm_to_zxV*/
      92             : GEN
      93        2268 : Flm_to_FlxV(GEN x, long sv)
      94        6153 : { pari_APPLY_type(t_VEC, Flv_to_Flx(gel(x,i), sv)) }
      95             : 
      96             : /*FlxC_to_ZXC=zxC_to_ZXC*/
      97             : GEN
      98      108837 : FlxC_to_ZXC(GEN x)
      99      561563 : { pari_APPLY_type(t_COL, Flx_to_ZX(gel(x,i))) }
     100             : 
     101             : /*FlxC_to_ZXC=zxV_to_ZXV*/
     102             : GEN
     103      593776 : FlxV_to_ZXV(GEN x)
     104     2406728 : { pari_APPLY_type(t_VEC, Flx_to_ZX(gel(x,i))) }
     105             : 
     106             : void
     107     2786157 : FlxV_to_ZXV_inplace(GEN v)
     108             : {
     109             :   long i;
     110     7443500 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
     111     2786066 : }
     112             : 
     113             : /*FlxM_to_ZXM=zxM_to_ZXM*/
     114             : GEN
     115        4570 : FlxM_to_ZXM(GEN x)
     116       15760 : { pari_APPLY_same(FlxC_to_ZXC(gel(x,i))) }
     117             : 
     118             : GEN
     119      394400 : FlxV_to_FlxX(GEN x, long v)
     120             : {
     121      394400 :   long i, l = lg(x)+1;
     122      394400 :   GEN z = cgetg(l,t_POL); z[1] = evalvarn(v);
     123      394400 :   x--;
     124     4961885 :   for (i=2; i<l ; i++) gel(z,i) = gel(x,i);
     125      394400 :   return FlxX_renormalize(z,l);
     126             : }
     127             : 
     128             : GEN
     129           0 : FlxM_to_FlxXV(GEN x, long v)
     130           0 : { pari_APPLY_type(t_COL, FlxV_to_FlxX(gel(x,i), v)) }
     131             : 
     132             : GEN
     133           0 : FlxM_Flx_add_shallow(GEN x, GEN y, ulong p)
     134             : {
     135           0 :   long l = lg(x), i, j;
     136           0 :   GEN z = cgetg(l,t_MAT);
     137             : 
     138           0 :   if (l==1) return z;
     139           0 :   if (l != lgcols(x)) pari_err_OP( "+", x, y);
     140           0 :   for (i=1; i<l; i++)
     141             :   {
     142           0 :     GEN zi = cgetg(l,t_COL), xi = gel(x,i);
     143           0 :     gel(z,i) = zi;
     144           0 :     for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
     145           0 :     gel(zi,i) = Flx_add(gel(zi,i), y, p);
     146             :   }
     147           0 :   return z;
     148             : }
     149             : 
     150             : /***********************************************************************/
     151             : /**                      Conversion to Flx                            **/
     152             : /***********************************************************************/
     153             : /* Take an integer and return a scalar polynomial mod p,  with evalvarn=vs */
     154             : GEN
     155    16798542 : Fl_to_Flx(ulong x, long sv) { return x? mkvecsmall2(sv, x): pol0_Flx(sv); }
     156             : 
     157             : /* a X^d */
     158             : GEN
     159      907729 : monomial_Flx(ulong a, long d, long vs)
     160             : {
     161             :   GEN P;
     162      907729 :   if (a==0) return pol0_Flx(vs);
     163      907729 :   P = const_vecsmall(d+2, 0);
     164      907742 :   P[1] = vs; P[d+2] = a; return P;
     165             : }
     166             : 
     167             : GEN
     168     2658272 : Z_to_Flx(GEN x, ulong p, long sv)
     169             : {
     170     2658272 :   long u = umodiu(x,p);
     171     2658424 :   return u? mkvecsmall2(sv, u): pol0_Flx(sv);
     172             : }
     173             : 
     174             : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
     175             : GEN
     176   158109992 : ZX_to_Flx(GEN x, ulong p)
     177             : {
     178   158109992 :   long i, lx = lg(x);
     179   158109992 :   GEN a = cgetg(lx, t_VECSMALL);
     180   157907221 :   a[1]=((ulong)x[1])&VARNBITS;
     181  1066354671 :   for (i=2; i<lx; i++) a[i] = umodiu(gel(x,i), p);
     182   158094522 :   return Flx_renormalize(a,lx);
     183             : }
     184             : 
     185             : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
     186             : GEN
     187     6031292 : zx_to_Flx(GEN x, ulong p)
     188             : {
     189     6031292 :   long i, lx = lg(x);
     190     6031292 :   GEN a = cgetg(lx, t_VECSMALL);
     191     6022118 :   a[1] = x[1];
     192    18498935 :   for (i=2; i<lx; i++) uel(a,i) = umodsu(x[i], p);
     193     6020562 :   return Flx_renormalize(a,lx);
     194             : }
     195             : 
     196             : ulong
     197    53652860 : Rg_to_Fl(GEN x, ulong p)
     198             : {
     199    53652860 :   switch(typ(x))
     200             :   {
     201    27962655 :     case t_INT: return umodiu(x, p);
     202      330726 :     case t_FRAC: {
     203      330726 :       ulong z = umodiu(gel(x,1), p);
     204      330727 :       if (!z) return 0;
     205      326548 :       return Fl_div(z, umodiu(gel(x,2), p), p);
     206             :     }
     207      205956 :     case t_PADIC: return padic_to_Fl(x, p);
     208    25153545 :     case t_INTMOD: {
     209    25153545 :       GEN q = gel(x,1), a = gel(x,2);
     210    25153545 :       if (absequaliu(q, p)) return itou(a);
     211           0 :       if (!dvdiu(q,p)) pari_err_MODULUS("Rg_to_Fl", q, utoipos(p));
     212           0 :       return umodiu(a, p);
     213             :     }
     214           0 :     default: pari_err_TYPE("Rg_to_Fl",x);
     215             :       return 0; /* LCOV_EXCL_LINE */
     216             :   }
     217             : }
     218             : 
     219             : ulong
     220     1664535 : Rg_to_F2(GEN x)
     221             : {
     222     1664535 :   switch(typ(x))
     223             :   {
     224      261497 :     case t_INT: return mpodd(x);
     225           0 :     case t_FRAC:
     226           0 :       if (!mpodd(gel(x,2))) (void)Fl_inv(0,2); /* error */
     227           0 :       return mpodd(gel(x,1));
     228           0 :     case t_PADIC:
     229           0 :       if (!absequaliu(gel(x,2),2)) pari_err_OP("",x, mkintmodu(1,2));
     230           0 :       if (valp(x) < 0) (void)Fl_inv(0,2);
     231           0 :       return valp(x) & 1;
     232     1403038 :     case t_INTMOD: {
     233     1403038 :       GEN q = gel(x,1), a = gel(x,2);
     234     1403038 :       if (mpodd(q)) pari_err_MODULUS("Rg_to_F2", q, gen_2);
     235     1403038 :       return mpodd(a);
     236             :     }
     237           0 :     default: pari_err_TYPE("Rg_to_F2",x);
     238             :       return 0; /* LCOV_EXCL_LINE */
     239             :   }
     240             : }
     241             : 
     242             : GEN
     243     2825977 : RgX_to_Flx(GEN x, ulong p)
     244             : {
     245     2825977 :   long i, lx = lg(x);
     246     2825977 :   GEN a = cgetg(lx, t_VECSMALL);
     247     2825977 :   a[1]=((ulong)x[1])&VARNBITS;
     248    22237186 :   for (i=2; i<lx; i++) a[i] = Rg_to_Fl(gel(x,i), p);
     249     2825977 :   return Flx_renormalize(a,lx);
     250             : }
     251             : 
     252             : GEN
     253           7 : RgXV_to_FlxV(GEN x, ulong p)
     254         175 : { pari_APPLY_type(t_VEC, RgX_to_Flx(gel(x,i), p)) }
     255             : 
     256             : /* If x is a POLMOD, assume modulus is a multiple of T. */
     257             : GEN
     258     3447389 : Rg_to_Flxq(GEN x, GEN T, ulong p)
     259             : {
     260     3447389 :   long ta, tx = typ(x), v = get_Flx_var(T);
     261             :   ulong pi;
     262             :   GEN a, b;
     263     3447395 :   if (is_const_t(tx))
     264             :   {
     265     3196901 :     if (tx == t_FFELT) return FF_to_Flxq(x);
     266     2465893 :     return Fl_to_Flx(Rg_to_Fl(x, p), v);
     267             :   }
     268      250498 :   switch(tx)
     269             :   {
     270        8576 :     case t_POLMOD:
     271        8576 :       b = gel(x,1);
     272        8576 :       a = gel(x,2); ta = typ(a);
     273        8576 :       if (is_const_t(ta)) return Fl_to_Flx(Rg_to_Fl(a, p), v);
     274        8422 :       b = RgX_to_Flx(b, p); if (b[1] != v) break;
     275        8422 :       a = RgX_to_Flx(a, p); if (Flx_equal(b,T)) return a;
     276           0 :       pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     277           0 :       if (lgpol(Flx_rem_pre(b,T,p,pi))==0) return Flx_rem_pre(a, T, p, pi);
     278           0 :       break;
     279      241922 :     case t_POL:
     280      241922 :       x = RgX_to_Flx(x,p);
     281      241922 :       if (x[1] != v) break;
     282      241922 :       return Flx_rem(x, T, p);
     283           0 :     case t_RFRAC:
     284           0 :       a = Rg_to_Flxq(gel(x,1), T,p);
     285           0 :       b = Rg_to_Flxq(gel(x,2), T,p);
     286           0 :       return Flxq_div(a,b, T,p);
     287             :   }
     288           0 :   pari_err_TYPE("Rg_to_Flxq",x);
     289             :   return NULL; /* LCOV_EXCL_LINE */
     290             : }
     291             : 
     292             : /***********************************************************************/
     293             : /**                   Basic operation on Flx                          **/
     294             : /***********************************************************************/
     295             : /* = zx_renormalize. Similar to normalizepol, in place */
     296             : GEN
     297  2011575084 : Flx_renormalize(GEN /*in place*/ x, long lx)
     298             : {
     299             :   long i;
     300  2268468815 :   for (i = lx-1; i>1; i--)
     301  2175686832 :     if (x[i]) break;
     302  2011575084 :   stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
     303  2009631027 :   setlg(x, i+1); return x;
     304             : }
     305             : 
     306             : GEN
     307     1885549 : Flx_red(GEN z, ulong p)
     308             : {
     309     1885549 :   long i, l = lg(z);
     310     1885549 :   GEN x = cgetg(l, t_VECSMALL);
     311     1885376 :   x[1] = z[1];
     312    34643416 :   for (i=2; i<l; i++) x[i] = uel(z,i)%p;
     313     1885376 :   return Flx_renormalize(x,l);
     314             : }
     315             : 
     316             : int
     317    26304907 : Flx_equal(GEN V, GEN W)
     318             : {
     319    26304907 :   long l = lg(V);
     320    26304907 :   if (lg(W) != l) return 0;
     321    27329354 :   while (--l > 1) /* do not compare variables, V[1] */
     322    26073551 :     if (V[l] != W[l]) return 0;
     323     1255803 :   return 1;
     324             : }
     325             : 
     326             : GEN
     327     2646277 : random_Flx(long d1, long vs, ulong p)
     328             : {
     329     2646277 :   long i, d = d1+2;
     330     2646277 :   GEN y = cgetg(d,t_VECSMALL); y[1] = vs;
     331    18323349 :   for (i=2; i<d; i++) y[i] = random_Fl(p);
     332     2646633 :   return Flx_renormalize(y,d);
     333             : }
     334             : 
     335             : static GEN
     336     7181202 : Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
     337             : {
     338             :   long i,lz;
     339             :   GEN z;
     340             : 
     341     7181202 :   if (ly>lx) swapspec(x,y, lx,ly);
     342     7181202 :   lz = lx+2; z = cgetg(lz, t_VECSMALL);
     343   107770597 :   for (i=0; i<ly; i++) z[i+2] = Fl_add(x[i], y[i], p);
     344    91187365 :   for (   ; i<lx; i++) z[i+2] = x[i];
     345     7181202 :   z[1] = 0; return Flx_renormalize(z, lz);
     346             : }
     347             : 
     348             : GEN
     349    60449592 : Flx_add(GEN x, GEN y, ulong p)
     350             : {
     351             :   long i,lz;
     352             :   GEN z;
     353    60449592 :   long lx=lg(x);
     354    60449592 :   long ly=lg(y);
     355    60449592 :   if (ly>lx) swapspec(x,y, lx,ly);
     356    60449592 :   lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
     357   559151104 :   for (i=2; i<ly; i++) z[i] = Fl_add(x[i], y[i], p);
     358   128239613 :   for (   ; i<lx; i++) z[i] = x[i];
     359    60404999 :   return Flx_renormalize(z, lz);
     360             : }
     361             : 
     362             : GEN
     363     9648174 : Flx_Fl_add(GEN y, ulong x, ulong p)
     364             : {
     365             :   GEN z;
     366             :   long lz, i;
     367     9648174 :   if (!lgpol(y))
     368      228298 :     return Fl_to_Flx(x,y[1]);
     369     9420913 :   lz=lg(y);
     370     9420913 :   z=cgetg(lz,t_VECSMALL);
     371     9419453 :   z[1]=y[1];
     372     9419453 :   z[2] = Fl_add(y[2],x,p);
     373    45985520 :   for(i=3;i<lz;i++)
     374    36565939 :     z[i] = y[i];
     375     9419581 :   if (lz==3) z = Flx_renormalize(z,lz);
     376     9419298 :   return z;
     377             : }
     378             : 
     379             : static GEN
     380      897417 : Flx_subspec(GEN x, GEN y, ulong p, long lx, long ly)
     381             : {
     382             :   long i,lz;
     383             :   GEN z;
     384             : 
     385      897417 :   if (ly <= lx)
     386             :   {
     387      897431 :     lz = lx+2; z = cgetg(lz, t_VECSMALL);
     388    54155631 :     for (i=0; i<ly; i++) z[i+2] = Fl_sub(x[i],y[i],p);
     389     1442930 :     for (   ; i<lx; i++) z[i+2] = x[i];
     390             :   }
     391             :   else
     392             :   {
     393           0 :     lz = ly+2; z = cgetg(lz, t_VECSMALL);
     394           0 :     for (i=0; i<lx; i++) z[i+2] = Fl_sub(x[i],y[i],p);
     395           0 :     for (   ; i<ly; i++) z[i+2] = Fl_neg(y[i],p);
     396             :   }
     397      896906 :   z[1] = 0; return Flx_renormalize(z, lz);
     398             : }
     399             : 
     400             : GEN
     401   137807930 : Flx_sub(GEN x, GEN y, ulong p)
     402             : {
     403   137807930 :   long i,lz,lx = lg(x), ly = lg(y);
     404             :   GEN z;
     405             : 
     406   137807930 :   if (ly <= lx)
     407             :   {
     408    86604193 :     lz = lx; z = cgetg(lz, t_VECSMALL);
     409   455177343 :     for (i=2; i<ly; i++) z[i] = Fl_sub(x[i],y[i],p);
     410   173033884 :     for (   ; i<lx; i++) z[i] = x[i];
     411             :   }
     412             :   else
     413             :   {
     414    51203737 :     lz = ly; z = cgetg(lz, t_VECSMALL);
     415   267388583 :     for (i=2; i<lx; i++) z[i] = Fl_sub(x[i],y[i],p);
     416   237423009 :     for (   ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
     417             :   }
     418   137745842 :   z[1]=x[1]; return Flx_renormalize(z, lz);
     419             : }
     420             : 
     421             : GEN
     422       43562 : Flx_Fl_sub(GEN y, ulong x, ulong p)
     423             : {
     424             :   GEN z;
     425       43562 :   long lz = lg(y), i;
     426       43562 :   if (lz==2)
     427         513 :     return Fl_to_Flx(Fl_neg(x, p),y[1]);
     428       43049 :   z = cgetg(lz, t_VECSMALL);
     429       43048 :   z[1] = y[1];
     430       43048 :   z[2] = Fl_sub(uel(y,2), x, p);
     431      281951 :   for(i=3; i<lz; i++)
     432      238903 :     z[i] = y[i];
     433       43048 :   if (lz==3) z = Flx_renormalize(z,lz);
     434       43048 :   return z;
     435             : }
     436             : 
     437             : static GEN
     438     3430553 : Flx_negspec(GEN x, ulong p, long l)
     439             : {
     440             :   long i;
     441     3430553 :   GEN z = cgetg(l+2, t_VECSMALL) + 2;
     442    21228304 :   for (i=0; i<l; i++) z[i] = Fl_neg(x[i], p);
     443     3430549 :   return z-2;
     444             : }
     445             : 
     446             : GEN
     447     3430623 : Flx_neg(GEN x, ulong p)
     448             : {
     449     3430623 :   GEN z = Flx_negspec(x+2, p, lgpol(x));
     450     3430898 :   z[1] = x[1];
     451     3430898 :   return z;
     452             : }
     453             : 
     454             : GEN
     455     1760375 : Flx_neg_inplace(GEN x, ulong p)
     456             : {
     457     1760375 :   long i, l = lg(x);
     458    52943181 :   for (i=2; i<l; i++)
     459    51182806 :     if (x[i]) x[i] = p - x[i];
     460     1760375 :   return x;
     461             : }
     462             : 
     463             : GEN
     464     2433686 : Flx_double(GEN y, ulong p)
     465             : {
     466             :   long i, l;
     467     2433686 :   GEN z = cgetg_copy(y, &l); z[1] = y[1];
     468    20624785 :   for(i=2; i<l; i++) z[i] = Fl_double(y[i], p);
     469     2433686 :   return Flx_renormalize(z, l);
     470             : }
     471             : GEN
     472     1040382 : Flx_triple(GEN y, ulong p)
     473             : {
     474             :   long i, l;
     475     1040382 :   GEN z = cgetg_copy(y, &l); z[1] = y[1];
     476     8294688 :   for(i=2; i<l; i++) z[i] = Fl_triple(y[i], p);
     477     1040382 :   return Flx_renormalize(z, l);
     478             : }
     479             : 
     480             : GEN
     481    18120076 : Flx_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
     482             : {
     483             :   GEN z;
     484             :   long i, l;
     485    18120076 :   if (!x) return pol0_Flx(y[1]);
     486    17431799 :   z = cgetg_copy(y, &l); z[1] = y[1];
     487    17430983 :   if (pi==0)
     488             :   {
     489    15477850 :     if (HIGHWORD(x | p))
     490           0 :       for(i=2; i<l; i++) z[i] = Fl_mul(uel(y,i), x, p);
     491             :     else
     492    91505343 :       for(i=2; i<l; i++) z[i] = (uel(y,i) * x) % p;
     493             :   } else
     494    14396663 :       for(i=2; i<l; i++) z[i] = Fl_mul_pre(uel(y,i), x, p, pi);
     495    17431778 :   return Flx_renormalize(z, l);
     496             : }
     497             : 
     498             : GEN
     499     6826316 : Flx_Fl_mul(GEN x, ulong y, ulong p)
     500     6826316 : { return Flx_Fl_mul_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
     501             : 
     502             : GEN
     503           0 : Flx_convol(GEN x, GEN y, ulong p)
     504             : {
     505           0 :   long lx = lg(x), ly = lg(y), i;
     506             :   GEN z;
     507           0 :   if (lx < ly) swapspec(x,y, lx,ly);
     508           0 :   z = cgetg(ly,t_VECSMALL); z[1] = x[1];
     509           0 :   for (i=2; i<ly; i++) uel(z,i) = Fl_mul(uel(x,i),uel(y,i), p);
     510           0 :   return Flx_renormalize(z, ly);
     511             : }
     512             : 
     513             : GEN
     514    11865590 : Flx_Fl_mul_to_monic(GEN y, ulong x, ulong p)
     515             : {
     516             :   GEN z;
     517             :   long i, l;
     518    11865590 :   z = cgetg_copy(y, &l); z[1] = y[1];
     519    11860975 :   if (HIGHWORD(x | p))
     520     5353209 :     for(i=2; i<l-1; i++) z[i] = Fl_mul(y[i], x, p);
     521             :   else
     522    26606237 :     for(i=2; i<l-1; i++) z[i] = (y[i] * x) % p;
     523    11860971 :   z[l-1] = 1; return z;
     524             : }
     525             : 
     526             : /* Return a*x^n if n>=0 and a\x^(-n) if n<0 */
     527             : GEN
     528    26185902 : Flx_shift(GEN a, long n)
     529             : {
     530    26185902 :   long i, l = lg(a);
     531             :   GEN  b;
     532    26185902 :   if (l==2 || !n) return Flx_copy(a);
     533    25842915 :   if (l+n<=2) return pol0_Flx(a[1]);
     534    25628326 :   b = cgetg(l+n, t_VECSMALL);
     535    25624589 :   b[1] = a[1];
     536    25624589 :   if (n < 0)
     537    71550817 :     for (i=2-n; i<l; i++) b[i+n] = a[i];
     538             :   else
     539             :   {
     540    49826118 :     for (i=0; i<n; i++) b[2+i] = 0;
     541   145009158 :     for (i=2; i<l; i++) b[i+n] = a[i];
     542             :   }
     543    25624589 :   return b;
     544             : }
     545             : 
     546             : GEN
     547    61050546 : Flx_normalize(GEN z, ulong p)
     548             : {
     549    61050546 :   long l = lg(z)-1;
     550    61050546 :   ulong p1 = z[l]; /* leading term */
     551    61050546 :   if (p1 == 1) return z;
     552    11844046 :   return Flx_Fl_mul_to_monic(z, Fl_inv(p1,p), p);
     553             : }
     554             : 
     555             : /* return (x * X^d) + y. Assume d > 0, shallow if x == 0*/
     556             : static GEN
     557     3689827 : Flx_addshift(GEN x, GEN y, ulong p, long d)
     558             : {
     559     3689827 :   GEN xd,yd,zd = (GEN)avma;
     560     3689827 :   long a,lz,ny = lgpol(y), nx = lgpol(x);
     561     3689827 :   long vs = x[1];
     562     3689827 :   if (nx == 0) return y;
     563     3687975 :   x += 2; y += 2; a = ny-d;
     564     3687975 :   if (a <= 0)
     565             :   {
     566       85237 :     lz = (a>nx)? ny+2: nx+d+2;
     567       85237 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
     568     1735135 :     while (xd > x) *--zd = *--xd;
     569       85237 :     x = zd + a;
     570      166914 :     while (zd > x) *--zd = 0;
     571             :   }
     572             :   else
     573             :   {
     574     3602738 :     xd = new_chunk(d); yd = y+d;
     575     3602738 :     x = Flx_addspec(x,yd,p, nx,a);
     576     3602738 :     lz = (a>nx)? ny+2: lg(x)+d;
     577   134239541 :     x += 2; while (xd > x) *--zd = *--xd;
     578             :   }
     579    60972427 :   while (yd > y) *--zd = *--yd;
     580     3687975 :   *--zd = vs;
     581     3687975 :   *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
     582             : }
     583             : 
     584             : /* shift polynomial + gerepile */
     585             : /* Do not set evalvarn*/
     586             : static GEN
     587   612907119 : Flx_shiftip(pari_sp av, GEN x, long v)
     588             : {
     589   612907119 :   long i, lx = lg(x), ly;
     590             :   GEN y;
     591   612907119 :   if (!v || lx==2) return gerepileuptoleaf(av, x);
     592   171383785 :   ly = lx + v; /* result length */
     593   171383785 :   (void)new_chunk(ly); /* check that result fits */
     594   171218233 :   x += lx; y = (GEN)av;
     595  1223470929 :   for (i = 2; i<lx; i++) *--y = *--x;
     596   695780704 :   for (i = 0; i< v; i++) *--y = 0;
     597   171218233 :   y -= 2; y[0] = evaltyp(t_VECSMALL) | evallg(ly);
     598   171614226 :   return gc_const((pari_sp)y, y);
     599             : }
     600             : 
     601             : static long
     602  2182170480 : get_Fl_threshold(ulong p, long mul, long mul2)
     603             : {
     604  2182170480 :   return SMALL_ULONG(p) ? mul: mul2;
     605             : }
     606             : 
     607             : #define BITS_IN_QUARTULONG (BITS_IN_HALFULONG >> 1)
     608             : #define QUARTMASK ((1UL<<BITS_IN_QUARTULONG)-1UL)
     609             : #define LLQUARTWORD(x) ((x) & QUARTMASK)
     610             : #define HLQUARTWORD(x) (((x) >> BITS_IN_QUARTULONG) & QUARTMASK)
     611             : #define LHQUARTWORD(x) (((x) >> (2*BITS_IN_QUARTULONG)) & QUARTMASK)
     612             : #define HHQUARTWORD(x) (((x) >> (3*BITS_IN_QUARTULONG)) & QUARTMASK)
     613             : INLINE long
     614     8043690 : maxbitcoeffpol(ulong p, long n)
     615             : {
     616     8043690 :   GEN z = muliu(sqru(p - 1), n);
     617     8038872 :   long b = expi(z) + 1;
     618             :   /* only do expensive bit-packing if it saves at least 1 limb */
     619     8041594 :   if (b <= BITS_IN_QUARTULONG)
     620             :   {
     621      934520 :     if (nbits2nlong(n*b) == (n + 3)>>2)
     622      107362 :       b = BITS_IN_QUARTULONG;
     623             :   }
     624     7107074 :   else if (b <= BITS_IN_HALFULONG)
     625             :   {
     626     1550909 :     if (nbits2nlong(n*b) == (n + 1)>>1)
     627        5590 :       b = BITS_IN_HALFULONG;
     628             :   }
     629             :   else
     630             :   {
     631     5556165 :     long l = lgefint(z) - 2;
     632     5556165 :     if (nbits2nlong(n*b) == n*l)
     633      306546 :       b = l*BITS_IN_LONG;
     634             :   }
     635     8041311 :   return b;
     636             : }
     637             : 
     638             : INLINE ulong
     639  3382329062 : Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
     640             : { /* Assume OK_ULONG*/
     641  3382329062 :   ulong p1 = 0;
     642             :   long i;
     643 16132704695 :   for (i=a; i<b; i++)
     644 12750375633 :     if (y[i])
     645             :     {
     646 10693217282 :       p1 += y[i] * x[-i];
     647 10693217282 :       if (p1 & HIGHBIT) p1 %= p;
     648             :     }
     649  3382329062 :   return p1 % p;
     650             : }
     651             : 
     652             : INLINE ulong
     653  1124129936 : Flx_mullimb(GEN x, GEN y, ulong p, ulong pi, long a, long b)
     654             : {
     655  1124129936 :   ulong p1 = 0;
     656             :   long i;
     657  3580606886 :   for (i=a; i<b; i++)
     658  2456486595 :     if (y[i])
     659  2414919603 :       p1 = Fl_addmul_pre(p1, y[i], x[-i], p, pi);
     660  1124120291 :   return p1;
     661             : }
     662             : 
     663             : /* assume nx >= ny > 0 */
     664             : static GEN
     665   328649908 : Flx_mulspec_basecase(GEN x, GEN y, ulong p, ulong pi, long nx, long ny)
     666             : {
     667             :   long i,lz,nz;
     668             :   GEN z;
     669             : 
     670   328649908 :   lz = nx+ny+1; nz = lz-2;
     671   328649908 :   z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
     672   328227675 :   if (!pi)
     673             :   {
     674  1139533964 :     for (i=0; i<ny; i++)z[i] = Flx_mullimb_ok(x+i,y,p,0,i+1);
     675   732075167 :     for (  ; i<nx; i++) z[i] = Flx_mullimb_ok(x+i,y,p,0,ny);
     676   888946863 :     for (  ; i<nz; i++) z[i] = Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
     677             :   }
     678             :   else
     679             :   {
     680   275785671 :     for (i=0; i<ny; i++)z[i] = Flx_mullimb(x+i,y,p,pi,0,i+1);
     681   194669867 :     for (  ; i<nx; i++) z[i] = Flx_mullimb(x+i,y,p,pi,0,ny);
     682   200347922 :     for (  ; i<nz; i++) z[i] = Flx_mullimb(x+i,y,p,pi,i-nx+1,ny);
     683             :   }
     684   328494280 :   z -= 2; return Flx_renormalize(z, lz);
     685             : }
     686             : 
     687             : static GEN
     688       12330 : int_to_Flx(GEN z, ulong p)
     689             : {
     690       12330 :   long i, l = lgefint(z);
     691       12330 :   GEN x = cgetg(l, t_VECSMALL);
     692     1061730 :   for (i=2; i<l; i++) x[i] = uel(z,i)%p;
     693       12322 :   return Flx_renormalize(x, l);
     694             : }
     695             : 
     696             : INLINE GEN
     697       10058 : Flx_mulspec_mulii(GEN a, GEN b, ulong p, long na, long nb)
     698             : {
     699       10058 :   GEN z=muliispec(a,b,na,nb);
     700       10065 :   return int_to_Flx(z,p);
     701             : }
     702             : 
     703             : static GEN
     704      505996 : Flx_to_int_halfspec(GEN a, long na)
     705             : {
     706             :   long j;
     707      505996 :   long n = (na+1)>>1UL;
     708      505996 :   GEN V = cgetipos(2+n);
     709             :   GEN w;
     710     1432388 :   for (w = int_LSW(V), j=0; j+1<na; j+=2, w=int_nextW(w))
     711      926392 :     *w = a[j]|(a[j+1]<<BITS_IN_HALFULONG);
     712      505996 :   if (j<na)
     713      338933 :     *w = a[j];
     714      505996 :   return V;
     715             : }
     716             : 
     717             : static GEN
     718      592412 : int_to_Flx_half(GEN z, ulong p)
     719             : {
     720             :   long i;
     721      592412 :   long lx = (lgefint(z)-2)*2+2;
     722      592412 :   GEN w, x = cgetg(lx, t_VECSMALL);
     723     2018042 :   for (w = int_LSW(z), i=2; i<lx; i+=2, w=int_nextW(w))
     724             :   {
     725     1425630 :     x[i]   = LOWWORD((ulong)*w)%p;
     726     1425630 :     x[i+1] = HIGHWORD((ulong)*w)%p;
     727             :   }
     728      592412 :   return Flx_renormalize(x, lx);
     729             : }
     730             : 
     731             : static GEN
     732        5454 : Flx_mulspec_halfmulii(GEN a, GEN b, ulong p, long na, long nb)
     733             : {
     734        5454 :   GEN A = Flx_to_int_halfspec(a,na);
     735        5454 :   GEN B = Flx_to_int_halfspec(b,nb);
     736        5454 :   GEN z = mulii(A,B);
     737        5454 :   return int_to_Flx_half(z,p);
     738             : }
     739             : 
     740             : static GEN
     741      204502 : Flx_to_int_quartspec(GEN a, long na)
     742             : {
     743             :   long j;
     744      204502 :   long n = (na+3)>>2UL;
     745      204502 :   GEN V = cgetipos(2+n);
     746             :   GEN w;
     747     4377862 :   for (w = int_LSW(V), j=0; j+3<na; j+=4, w=int_nextW(w))
     748     4173359 :     *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG))|(a[j+3]<<(3*BITS_IN_QUARTULONG));
     749      204503 :   switch (na-j)
     750             :   {
     751      116434 :   case 3:
     752      116434 :     *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG));
     753      116434 :     break;
     754       34450 :   case 2:
     755       34450 :     *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG);
     756       34450 :     break;
     757       27219 :   case 1:
     758       27219 :     *w = a[j];
     759       27219 :     break;
     760       26400 :   case 0:
     761       26400 :     break;
     762             :   }
     763      204503 :   return V;
     764             : }
     765             : 
     766             : static GEN
     767      107362 : int_to_Flx_quart(GEN z, ulong p)
     768             : {
     769             :   long i;
     770      107362 :   long lx = (lgefint(z)-2)*4+2;
     771      107362 :   GEN w, x = cgetg(lx, t_VECSMALL);
     772     4874009 :   for (w = int_LSW(z), i=2; i<lx; i+=4, w=int_nextW(w))
     773             :   {
     774     4766647 :     x[i]   = LLQUARTWORD((ulong)*w)%p;
     775     4766647 :     x[i+1] = HLQUARTWORD((ulong)*w)%p;
     776     4766647 :     x[i+2] = LHQUARTWORD((ulong)*w)%p;
     777     4766647 :     x[i+3] = HHQUARTWORD((ulong)*w)%p;
     778             :   }
     779      107362 :   return Flx_renormalize(x, lx);
     780             : }
     781             : 
     782             : static GEN
     783       97141 : Flx_mulspec_quartmulii(GEN a, GEN b, ulong p, long na, long nb)
     784             : {
     785       97141 :   GEN A = Flx_to_int_quartspec(a,na);
     786       97141 :   GEN B = Flx_to_int_quartspec(b,nb);
     787       97141 :   GEN z = mulii(A,B);
     788       97141 :   return int_to_Flx_quart(z,p);
     789             : }
     790             : 
     791             : /*Eval x in 2^(k*BIL) in linear time, k==2 or 3*/
     792             : static GEN
     793      580334 : Flx_eval2BILspec(GEN x, long k, long l)
     794             : {
     795      580334 :   long i, lz = k*l, ki;
     796      580334 :   GEN pz = cgetipos(2+lz);
     797    16327382 :   for (i=0; i < lz; i++)
     798    15747048 :     *int_W(pz,i) = 0UL;
     799     8453858 :   for (i=0, ki=0; i<l; i++, ki+=k)
     800     7873524 :     *int_W(pz,ki) = x[i];
     801      580334 :   return int_normalize(pz,0);
     802             : }
     803             : 
     804             : static GEN
     805      297153 : Z_mod2BIL_Flx_2(GEN x, long d, ulong p)
     806             : {
     807      297153 :   long i, offset, lm = lgefint(x)-2, l = d+3;
     808      297153 :   ulong pi = get_Fl_red(p);
     809      297153 :   GEN pol = cgetg(l, t_VECSMALL);
     810      297151 :   pol[1] = 0;
     811     7990113 :   for (i=0, offset=0; offset+1 < lm; i++, offset += 2)
     812     7692960 :     pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
     813      297153 :   if (offset < lm)
     814      224989 :     pol[i+2] = (*int_W(x,offset)) % p;
     815      297153 :   return Flx_renormalize(pol,l);
     816             : }
     817             : 
     818             : static GEN
     819           0 : Z_mod2BIL_Flx_3(GEN x, long d, ulong p)
     820             : {
     821           0 :   long i, offset, lm = lgefint(x)-2, l = d+3;
     822           0 :   ulong pi = get_Fl_red(p);
     823           0 :   GEN pol = cgetg(l, t_VECSMALL);
     824           0 :   pol[1] = 0;
     825           0 :   for (i=0, offset=0; offset+2 < lm; i++, offset += 3)
     826           0 :     pol[i+2] = remlll_pre(*int_W(x,offset+2), *int_W(x,offset+1),
     827           0 :                           *int_W(x,offset), p, pi);
     828           0 :   if (offset+1 < lm)
     829           0 :     pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
     830           0 :   else if (offset < lm)
     831           0 :     pol[i+2] = (*int_W(x,offset)) % p;
     832           0 :   return Flx_renormalize(pol,l);
     833             : }
     834             : 
     835             : static GEN
     836      294223 : Z_mod2BIL_Flx(GEN x, long bs, long d, ulong p)
     837             : {
     838      294223 :   return bs==2 ? Z_mod2BIL_Flx_2(x, d, p): Z_mod2BIL_Flx_3(x, d, p);
     839             : }
     840             : 
     841             : static GEN
     842      282722 : Flx_mulspec_mulii_inflate(GEN x, GEN y, long N, ulong p, long nx, long ny)
     843             : {
     844      282722 :   pari_sp av = avma;
     845      282722 :   GEN z = mulii(Flx_eval2BILspec(x,N,nx), Flx_eval2BILspec(y,N,ny));
     846      282722 :   return gerepileupto(av, Z_mod2BIL_Flx(z, N, nx+ny-2, p));
     847             : }
     848             : 
     849             : static GEN
     850    20165108 : kron_pack_Flx_spec_bits(GEN x, long b, long l) {
     851             :   GEN y;
     852             :   long i;
     853    20165108 :   if (l == 0)
     854     3426611 :     return gen_0;
     855    16738497 :   y = cgetg(l + 1, t_VECSMALL);
     856   810692214 :   for(i = 1; i <= l; i++)
     857   793963525 :     y[i] = x[l - i];
     858    16728689 :   return nv_fromdigits_2k(y, b);
     859             : }
     860             : 
     861             : /* assume b < BITS_IN_LONG */
     862             : static GEN
     863     5694929 : kron_unpack_Flx_bits_narrow(GEN z, long b, ulong p) {
     864     5694929 :   GEN v = binary_2k_nv(z, b), x;
     865     5694932 :   long i, l = lg(v) + 1;
     866     5694932 :   x = cgetg(l, t_VECSMALL);
     867   625500562 :   for (i = 2; i < l; i++)
     868   619805307 :     x[i] = v[l - i] % p;
     869     5695255 :   return Flx_renormalize(x, l);
     870             : }
     871             : 
     872             : static GEN
     873     5187322 : kron_unpack_Flx_bits_wide(GEN z, long b, ulong p, ulong pi) {
     874     5187322 :   GEN v = binary_2k(z, b), x, y;
     875     5188224 :   long i, l = lg(v) + 1, ly;
     876     5188224 :   x = cgetg(l, t_VECSMALL);
     877   227624147 :   for (i = 2; i < l; i++) {
     878   222434718 :     y = gel(v, l - i);
     879   222434718 :     ly = lgefint(y);
     880   222434718 :     switch (ly) {
     881     6278306 :     case 2: x[i] = 0; break;
     882    28852510 :     case 3: x[i] = *int_W_lg(y, 0, ly) % p; break;
     883   171445877 :     case 4: x[i] = remll_pre(*int_W_lg(y, 1, ly), *int_W_lg(y, 0, ly), p, pi); break;
     884    31715723 :     case 5: x[i] = remlll_pre(*int_W_lg(y, 2, ly), *int_W_lg(y, 1, ly),
     885    15858025 :                               *int_W_lg(y, 0, ly), p, pi); break;
     886           0 :     default: x[i] = umodiu(gel(v, l - i), p);
     887             :     }
     888             :   }
     889     5189429 :   return Flx_renormalize(x, l);
     890             : }
     891             : 
     892             : static GEN
     893     6932687 : Flx_mulspec_Kronecker(GEN A, GEN B, long b, ulong p, long lA, long lB)
     894             : {
     895             :   GEN C, D;
     896     6932687 :   pari_sp av = avma;
     897     6932687 :   A =  kron_pack_Flx_spec_bits(A, b, lA);
     898     6937710 :   B =  kron_pack_Flx_spec_bits(B, b, lB);
     899     6937859 :   C = gerepileuptoint(av, mulii(A, B));
     900     6935379 :   if (b < BITS_IN_LONG)
     901     2121802 :     D =  kron_unpack_Flx_bits_narrow(C, b, p);
     902             :   else
     903             :   {
     904     4813577 :     ulong pi = get_Fl_red(p);
     905     4811771 :     D = kron_unpack_Flx_bits_wide(C, b, p, pi);
     906             :   }
     907     6935158 :   return D;
     908             : }
     909             : 
     910             : static GEN
     911      688973 : Flx_sqrspec_Kronecker(GEN A, long b, ulong p, long lA)
     912             : {
     913             :   GEN C, D;
     914      688973 :   A =  kron_pack_Flx_spec_bits(A, b, lA);
     915      689053 :   C = sqri(A);
     916      689075 :   if (b < BITS_IN_LONG)
     917      480723 :     D =  kron_unpack_Flx_bits_narrow(C, b, p);
     918             :   else
     919             :   {
     920      208352 :     ulong pi = get_Fl_red(p);
     921      208345 :     D = kron_unpack_Flx_bits_wide(C, b, p, pi);
     922             :   }
     923      689025 :   return D;
     924             : }
     925             : 
     926             : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
     927             :  * b+2 were sent instead. na, nb = number of terms of a, b.
     928             :  * Only c, c0, c1, c2 are genuine GEN.
     929             :  */
     930             : static GEN
     931   365755536 : Flx_mulspec(GEN a, GEN b, ulong p, ulong pi, long na, long nb)
     932             : {
     933             :   GEN a0,c,c0;
     934   365755536 :   long n0, n0a, i, v = 0;
     935             :   pari_sp av;
     936             : 
     937   470222098 :   while (na && !a[0]) { a++; na--; v++; }
     938   551460889 :   while (nb && !b[0]) { b++; nb--; v++; }
     939   365755536 :   if (na < nb) swapspec(a,b, na,nb);
     940   365755536 :   if (!nb) return pol0_Flx(0);
     941             : 
     942   337573074 :   av = avma;
     943   337573074 :   if (nb >= get_Fl_threshold(p, Flx_MUL_MULII_LIMIT, Flx_MUL2_MULII_LIMIT))
     944             :   {
     945     7330969 :     long m = maxbitcoeffpol(p,nb);
     946     7327201 :     switch (m)
     947             :     {
     948       97141 :     case BITS_IN_QUARTULONG:
     949       97141 :       return Flx_shiftip(av,Flx_mulspec_quartmulii(a,b,p,na,nb), v);
     950        5454 :     case BITS_IN_HALFULONG:
     951        5454 :       return Flx_shiftip(av,Flx_mulspec_halfmulii(a,b,p,na,nb), v);
     952       10058 :     case BITS_IN_LONG:
     953       10058 :       return Flx_shiftip(av,Flx_mulspec_mulii(a,b,p,na,nb), v);
     954      282722 :     case 2*BITS_IN_LONG:
     955      282722 :       return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,2,p,na,nb), v);
     956           0 :     case 3*BITS_IN_LONG:
     957           0 :       return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,3,p,na,nb), v);
     958     6931826 :     default:
     959     6931826 :       return Flx_shiftip(av,Flx_mulspec_Kronecker(a,b,m,p,na,nb), v);
     960             :     }
     961             :   }
     962   330407143 :   if (nb < get_Fl_threshold(p, Flx_MUL_KARATSUBA_LIMIT, Flx_MUL2_KARATSUBA_LIMIT))
     963   328637486 :     return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,pi,na,nb), v);
     964     1814005 :   i=(na>>1); n0=na-i; na=i;
     965     1814005 :   a0=a+n0; n0a=n0;
     966     2591242 :   while (n0a && !a[n0a-1]) n0a--;
     967             : 
     968     1814005 :   if (nb > n0)
     969             :   {
     970             :     GEN b0,c1,c2;
     971             :     long n0b;
     972             : 
     973     1760375 :     nb -= n0; b0 = b+n0; n0b = n0;
     974     2846603 :     while (n0b && !b[n0b-1]) n0b--;
     975     1760375 :     c =  Flx_mulspec(a,b,p,pi,n0a,n0b);
     976     1760375 :     c0 = Flx_mulspec(a0,b0,p,pi,na,nb);
     977             : 
     978     1760375 :     c2 = Flx_addspec(a0,a,p,na,n0a);
     979     1760375 :     c1 = Flx_addspec(b0,b,p,nb,n0b);
     980             : 
     981     1760375 :     c1 = Flx_mul_pre(c1,c2,p,pi);
     982     1760375 :     c2 = Flx_add(c0,c,p);
     983             : 
     984     1760375 :     c2 = Flx_neg_inplace(c2,p);
     985     1760375 :     c2 = Flx_add(c1,c2,p);
     986     1760375 :     c0 = Flx_addshift(c0,c2 ,p, n0);
     987             :   }
     988             :   else
     989             :   {
     990       53630 :     c  = Flx_mulspec(a,b,p,pi,n0a,nb);
     991       53630 :     c0 = Flx_mulspec(a0,b,p,pi,na,nb);
     992             :   }
     993     1814005 :   c0 = Flx_addshift(c0,c,p,n0);
     994     1814005 :   return Flx_shiftip(av,c0, v);
     995             : }
     996             : 
     997             : GEN
     998   360084010 : Flx_mul_pre(GEN x, GEN y, ulong p, ulong pi)
     999             : {
    1000   360084010 :   GEN z = Flx_mulspec(x+2,y+2,p, pi, lgpol(x),lgpol(y));
    1001   360305829 :   z[1] = x[1]; return z;
    1002             : }
    1003             : GEN
    1004    27918486 : Flx_mul(GEN x, GEN y, ulong p)
    1005    27918486 : { return Flx_mul_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    1006             : 
    1007             : static GEN
    1008   276203208 : Flx_sqrspec_basecase(GEN x, ulong p, ulong pi, long nx)
    1009             : {
    1010             :   long i, lz, nz;
    1011             :   ulong p1;
    1012             :   GEN z;
    1013             : 
    1014   276203208 :   if (!nx) return pol0_Flx(0);
    1015   276203208 :   lz = (nx << 1) + 1, nz = lz-2;
    1016   276203208 :   z = cgetg(lz, t_VECSMALL) + 2;
    1017   275445157 :   if (!pi)
    1018             :   {
    1019   212449061 :     z[0] = x[0]*x[0]%p;
    1020   916275649 :     for (i=1; i<nx; i++)
    1021             :     {
    1022   703623536 :       p1 = Flx_mullimb_ok(x+i,x,p,0, (i+1)>>1);
    1023   703826588 :       p1 <<= 1;
    1024   703826588 :       if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
    1025   703826588 :       z[i] = p1 % p;
    1026             :     }
    1027   920617208 :     for (  ; i<nz; i++)
    1028             :     {
    1029   707521702 :       p1 = Flx_mullimb_ok(x+i,x,p,i-nx+1, (i+1)>>1);
    1030   707965095 :       p1 <<= 1;
    1031   707965095 :       if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
    1032   707965095 :       z[i] = p1 % p;
    1033             :     }
    1034             :   }
    1035             :   else
    1036             :   {
    1037    62996096 :     z[0] = Fl_sqr_pre(x[0], p, pi);
    1038   409287928 :     for (i=1; i<nx; i++)
    1039             :     {
    1040   346400990 :       p1 = Flx_mullimb(x+i,x,p,pi,0, (i+1)>>1);
    1041   346650900 :       p1 = Fl_add(p1, p1, p);
    1042   346181151 :       if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
    1043   346213292 :       z[i] = p1;
    1044             :     }
    1045   409396837 :     for (  ; i<nz; i++)
    1046             :     {
    1047   346405140 :       p1 = Flx_mullimb(x+i,x,p,pi,i-nx+1, (i+1)>>1);
    1048   347172882 :       p1 = Fl_add(p1, p1, p);
    1049   346782063 :       if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
    1050   346509899 :       z[i] = p1;
    1051             :     }
    1052             :   }
    1053   276087203 :   z -= 2; return Flx_renormalize(z, lz);
    1054             : }
    1055             : 
    1056             : static GEN
    1057        2265 : Flx_sqrspec_sqri(GEN a, ulong p, long na)
    1058             : {
    1059        2265 :   GEN z=sqrispec(a,na);
    1060        2265 :   return int_to_Flx(z,p);
    1061             : }
    1062             : 
    1063             : static GEN
    1064         136 : Flx_sqrspec_halfsqri(GEN a, ulong p, long na)
    1065             : {
    1066         136 :   GEN z = sqri(Flx_to_int_halfspec(a,na));
    1067         136 :   return int_to_Flx_half(z,p);
    1068             : }
    1069             : 
    1070             : static GEN
    1071       10221 : Flx_sqrspec_quartsqri(GEN a, ulong p, long na)
    1072             : {
    1073       10221 :   GEN z = sqri(Flx_to_int_quartspec(a,na));
    1074       10221 :   return int_to_Flx_quart(z,p);
    1075             : }
    1076             : 
    1077             : static GEN
    1078       11501 : Flx_sqrspec_sqri_inflate(GEN x, long N, ulong p, long nx)
    1079             : {
    1080       11501 :   pari_sp av = avma;
    1081       11501 :   GEN  z = sqri(Flx_eval2BILspec(x,N,nx));
    1082       11501 :   return gerepileupto(av, Z_mod2BIL_Flx(z, N, (nx-1)*2, p));
    1083             : }
    1084             : 
    1085             : static GEN
    1086   276685689 : Flx_sqrspec(GEN a, ulong p, ulong pi, long na)
    1087             : {
    1088             :   GEN a0, c, c0;
    1089   276685689 :   long n0, n0a, i, v = 0, m;
    1090             :   pari_sp av;
    1091             : 
    1092   395844136 :   while (na && !a[0]) { a++; na--; v += 2; }
    1093   276685689 :   if (!na) return pol0_Flx(0);
    1094             : 
    1095   276449334 :   av = avma;
    1096   276449334 :   if (na >= get_Fl_threshold(p, Flx_SQR_SQRI_LIMIT, Flx_SQR2_SQRI_LIMIT))
    1097             :   {
    1098      713112 :     m = maxbitcoeffpol(p,na);
    1099      713088 :     switch(m)
    1100             :     {
    1101       10221 :     case BITS_IN_QUARTULONG:
    1102       10221 :       return Flx_shiftip(av, Flx_sqrspec_quartsqri(a,p,na), v);
    1103         136 :     case BITS_IN_HALFULONG:
    1104         136 :       return Flx_shiftip(av, Flx_sqrspec_halfsqri(a,p,na), v);
    1105        2265 :     case BITS_IN_LONG:
    1106        2265 :       return Flx_shiftip(av, Flx_sqrspec_sqri(a,p,na), v);
    1107       11501 :     case 2*BITS_IN_LONG:
    1108       11501 :       return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,2,p,na), v);
    1109           0 :     case 3*BITS_IN_LONG:
    1110           0 :       return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,3,p,na), v);
    1111      688965 :     default:
    1112      688965 :       return Flx_shiftip(av, Flx_sqrspec_Kronecker(a,m,p,na), v);
    1113             :     }
    1114             :   }
    1115   276036384 :   if (na < get_Fl_threshold(p, Flx_SQR_KARATSUBA_LIMIT, Flx_SQR2_KARATSUBA_LIMIT))
    1116   275928300 :     return Flx_shiftip(av, Flx_sqrspec_basecase(a,p,pi,na), v);
    1117       57733 :   i=(na>>1); n0=na-i; na=i;
    1118       57733 :   a0=a+n0; n0a=n0;
    1119       72524 :   while (n0a && !a[n0a-1]) n0a--;
    1120             : 
    1121       57733 :   c = Flx_sqrspec(a,p,pi,n0a);
    1122       57733 :   c0= Flx_sqrspec(a0,p,pi,na);
    1123       57733 :   if (p == 2) n0 *= 2;
    1124             :   else
    1125             :   {
    1126       57714 :     GEN c1, t = Flx_addspec(a0,a,p,na,n0a);
    1127       57714 :     t = Flx_sqr_pre(t,p,pi);
    1128       57714 :     c1= Flx_add(c0,c, p);
    1129       57714 :     c1= Flx_sub(t, c1, p);
    1130       57714 :     c0 = Flx_addshift(c0,c1,p,n0);
    1131             :   }
    1132       57733 :   c0 = Flx_addshift(c0,c,p,n0);
    1133       57733 :   return Flx_shiftip(av,c0,v);
    1134             : }
    1135             : 
    1136             : GEN
    1137   276446154 : Flx_sqr_pre(GEN x, ulong p, ulong pi)
    1138             : {
    1139   276446154 :   GEN z = Flx_sqrspec(x+2,p, pi, lgpol(x));
    1140   277481864 :   z[1] = x[1]; return z;
    1141             : }
    1142             : GEN
    1143      356044 : Flx_sqr(GEN x, ulong p)
    1144      356044 : { return Flx_sqr_pre(x, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    1145             : 
    1146             : GEN
    1147        8284 : Flx_powu_pre(GEN x, ulong n, ulong p, ulong pi)
    1148             : {
    1149        8284 :   GEN y = pol1_Flx(x[1]), z;
    1150             :   ulong m;
    1151        8281 :   if (n == 0) return y;
    1152        8281 :   m = n; z = x;
    1153             :   for (;;)
    1154             :   {
    1155       31921 :     if (m&1UL) y = Flx_mul_pre(y,z, p, pi);
    1156       31919 :     m >>= 1; if (!m) return y;
    1157       23636 :     z = Flx_sqr_pre(z, p, pi);
    1158             :   }
    1159             : }
    1160             : GEN
    1161           0 : Flx_powu(GEN x, ulong n, ulong p)
    1162             : {
    1163           0 :   if (n == 0) return pol1_Flx(x[1]);
    1164           0 :   return Flx_powu_pre(x, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p));
    1165             : }
    1166             : 
    1167             : GEN
    1168       14237 : Flx_halve(GEN y, ulong p)
    1169             : {
    1170             :   GEN z;
    1171             :   long i, l;
    1172       14237 :   z = cgetg_copy(y, &l); z[1] = y[1];
    1173       59382 :   for(i=2; i<l; i++) uel(z,i) = Fl_halve(uel(y,i), p);
    1174       14237 :   return z;
    1175             : }
    1176             : 
    1177             : static GEN
    1178     7094938 : Flx_recipspec(GEN x, long l, long n)
    1179             : {
    1180             :   long i;
    1181     7094938 :   GEN z=cgetg(n+2,t_VECSMALL)+2;
    1182   116210387 :   for(i=0; i<l; i++)
    1183   109117114 :     z[n-i-1] = x[i];
    1184    15515021 :   for(   ; i<n; i++)
    1185     8421748 :     z[n-i-1] = 0;
    1186     7093273 :   return Flx_renormalize(z-2,n+2);
    1187             : }
    1188             : 
    1189             : GEN
    1190           0 : Flx_recip(GEN x)
    1191             : {
    1192           0 :   GEN z=Flx_recipspec(x+2,lgpol(x),lgpol(x));
    1193           0 :   z[1]=x[1];
    1194           0 :   return z;
    1195             : }
    1196             : 
    1197             : /* Return h^degpol(P) P(x / h) */
    1198             : GEN
    1199        1117 : Flx_rescale(GEN P, ulong h, ulong p)
    1200             : {
    1201        1117 :   long i, l = lg(P);
    1202        1117 :   GEN Q = cgetg(l,t_VECSMALL);
    1203        1117 :   ulong hi = h;
    1204        1117 :   Q[l-1] = P[l-1];
    1205       12538 :   for (i=l-2; i>=2; i--)
    1206             :   {
    1207       12538 :     Q[i] = Fl_mul(P[i], hi, p);
    1208       12538 :     if (i == 2) break;
    1209       11421 :     hi = Fl_mul(hi,h, p);
    1210             :   }
    1211        1117 :   Q[1] = P[1]; return Q;
    1212             : }
    1213             : 
    1214             : /* x/polrecip(P)+O(x^n); allow pi = 0 */
    1215             : static GEN
    1216      134093 : Flx_invBarrett_basecase(GEN T, ulong p, ulong pi)
    1217             : {
    1218      134093 :   long i, l=lg(T)-1, lr=l-1, k;
    1219      134093 :   GEN r=cgetg(lr,t_VECSMALL); r[1] = T[1];
    1220      134094 :   r[2] = 1;
    1221      134094 :   if (!pi)
    1222      764068 :     for (i=3;i<lr;i++)
    1223             :     {
    1224      757076 :       ulong u = uel(T, l-i+2);
    1225    45377003 :       for (k=3; k<i; k++)
    1226    44619927 :         { u += uel(T,l-i+k) * uel(r, k); if (u & HIGHBIT) u %= p; }
    1227      757076 :       r[i] = Fl_neg(u % p, p);
    1228             :     }
    1229             :   else
    1230     2107904 :     for (i=3;i<lr;i++)
    1231             :     {
    1232     1980804 :       ulong u = Fl_neg(uel(T,l-i+2), p);
    1233    59503840 :       for (k=3; k<i; k++)
    1234             :       {
    1235    57523038 :         ulong t = Fl_neg(uel(T,l-i+k), p);
    1236    57523033 :         u = Fl_addmul_pre(u, t, uel(r,k), p, pi);
    1237             :       }
    1238     1980802 :       r[i] = u;
    1239             :     }
    1240      134092 :   return Flx_renormalize(r,lr);
    1241             : }
    1242             : 
    1243             : /* Return new lgpol */
    1244             : static long
    1245     2130728 : Flx_lgrenormalizespec(GEN x, long lx)
    1246             : {
    1247             :   long i;
    1248     7413105 :   for (i = lx-1; i>=0; i--)
    1249     7412279 :     if (x[i]) break;
    1250     2130728 :   return i+1;
    1251             : }
    1252             : /* allow pi = 0 */
    1253             : static GEN
    1254       23029 : Flx_invBarrett_Newton(GEN T, ulong p, ulong pi)
    1255             : {
    1256       23029 :   long nold, lx, lz, lq, l = degpol(T), lQ;
    1257       23029 :   GEN q, y, z, x = zero_zv(l+1) + 2;
    1258       23029 :   ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
    1259             :   pari_sp av;
    1260             : 
    1261       23029 :   y = T+2;
    1262       23029 :   q = Flx_recipspec(y,l+1,l+1); lQ = lgpol(q); q+=2;
    1263       23028 :   av = avma;
    1264             :   /* We work on _spec_ Flx's, all the l[xzq12] below are lgpol's */
    1265             : 
    1266             :   /* initialize */
    1267       23028 :   x[0] = Fl_inv(q[0], p);
    1268       23028 :   if (lQ>1 && q[1])
    1269        5108 :   {
    1270        5108 :     ulong u = q[1];
    1271        5108 :     if (x[0] != 1) u = Fl_mul(u, Fl_sqr(x[0],p), p);
    1272        5108 :     x[1] = p - u; lx = 2;
    1273             :   }
    1274             :   else
    1275       17920 :     lx = 1;
    1276       23028 :   nold = 1;
    1277      158181 :   for (; mask > 1; set_avma(av))
    1278             :   { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
    1279      135163 :     long i, lnew, nnew = nold << 1;
    1280             : 
    1281      135163 :     if (mask & 1) nnew--;
    1282      135163 :     mask >>= 1;
    1283             : 
    1284      135163 :     lnew = nnew + 1;
    1285      135163 :     lq = Flx_lgrenormalizespec(q, minss(lQ, lnew));
    1286      135165 :     z = Flx_mulspec(x, q, p, pi, lx, lq); /* FIXME: high product */
    1287      135154 :     lz = lgpol(z); if (lz > lnew) lz = lnew;
    1288      135154 :     z += 2;
    1289             :     /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
    1290      290179 :     for (i = nold; i < lz; i++) if (z[i]) break;
    1291      135154 :     nold = nnew;
    1292      135154 :     if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
    1293             : 
    1294             :     /* z + i represents (x*q - 1) / t^i */
    1295      100740 :     lz = Flx_lgrenormalizespec (z+i, lz-i);
    1296      100742 :     z = Flx_mulspec(x, z+i, p, pi, lx, lz); /* FIXME: low product */
    1297      100740 :     lz = lgpol(z); z += 2;
    1298      100739 :     if (lz > lnew-i) lz = Flx_lgrenormalizespec(z, lnew-i);
    1299             : 
    1300      100739 :     lx = lz+ i;
    1301      100739 :     y  = x + i; /* x -= z * t^i, in place */
    1302      915334 :     for (i = 0; i < lz; i++) y[i] = Fl_neg(z[i], p);
    1303             :   }
    1304       23029 :   x -= 2; setlg(x, lx + 2); x[1] = T[1];
    1305       23029 :   return x;
    1306             : }
    1307             : 
    1308             : /* allow pi = 0 */
    1309             : static GEN
    1310      158513 : Flx_invBarrett_pre(GEN T, ulong p, ulong pi)
    1311             : {
    1312      158513 :   pari_sp ltop = avma;
    1313      158513 :   long l = lgpol(T);
    1314             :   GEN r;
    1315      158513 :   if (l < 3) return pol0_Flx(T[1]);
    1316      157123 :   if (l < get_Fl_threshold(p, Flx_INVBARRETT_LIMIT, Flx_INVBARRETT2_LIMIT))
    1317             :   {
    1318      134094 :     ulong c = T[l+1];
    1319      134094 :     if (c != 1)
    1320             :     {
    1321       98117 :       ulong ci = Fl_inv(c,p);
    1322       98118 :       T = Flx_Fl_mul_pre(T, ci, p, pi);
    1323       98116 :       r = Flx_invBarrett_basecase(T, p, pi);
    1324       98117 :       r = Flx_Fl_mul_pre(r, ci, p, pi);
    1325             :     }
    1326             :     else
    1327       35977 :       r = Flx_invBarrett_basecase(T, p, pi);
    1328             :   }
    1329             :   else
    1330       23029 :     r = Flx_invBarrett_Newton(T, p, pi);
    1331      157121 :   return gerepileuptoleaf(ltop, r);
    1332             : }
    1333             : GEN
    1334           0 : Flx_invBarrett(GEN T, ulong p)
    1335           0 : { return Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    1336             : 
    1337             : /* allow pi = 0 */
    1338             : GEN
    1339    97032577 : Flx_get_red_pre(GEN T, ulong p, ulong pi)
    1340             : {
    1341    97032577 :   if (typ(T)!=t_VECSMALL
    1342    97001389 :     || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
    1343             :                                        Flx_BARRETT2_LIMIT))
    1344    97003554 :     return T;
    1345        7611 :   retmkvec2(Flx_invBarrett_pre(T, p, pi),T);
    1346             : }
    1347             : GEN
    1348    13736479 : Flx_get_red(GEN T, ulong p)
    1349             : {
    1350    13736479 :   if (typ(T)!=t_VECSMALL
    1351    13736548 :     || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
    1352             :                                        Flx_BARRETT2_LIMIT))
    1353    13729776 :     return T;
    1354        5194 :   retmkvec2(Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)),T);
    1355             : }
    1356             : 
    1357             : /* separate from Flx_divrem for maximal speed. */
    1358             : static GEN
    1359   712002910 : Flx_rem_basecase(GEN x, GEN y, ulong p, ulong pi)
    1360             : {
    1361             :   pari_sp av;
    1362             :   GEN z, c;
    1363             :   long dx,dy,dy1,dz,i,j;
    1364             :   ulong p1,inv;
    1365   712002910 :   long vs=x[1];
    1366             : 
    1367   712002910 :   dy = degpol(y); if (!dy) return pol0_Flx(x[1]);
    1368   680774396 :   dx = degpol(x);
    1369   680818294 :   dz = dx-dy; if (dz < 0) return Flx_copy(x);
    1370   680818294 :   x += 2; y += 2;
    1371   680818294 :   inv = y[dy];
    1372   680818294 :   if (inv != 1UL) inv = Fl_inv(inv,p);
    1373   823084634 :   for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
    1374             : 
    1375   682464180 :   c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
    1376   680483640 :   z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
    1377             : 
    1378   679230197 :   if (!pi)
    1379             :   {
    1380   472581373 :     z[dz] = (inv*x[dx]) % p;
    1381  1792779974 :     for (i=dx-1; i>=dy; --i)
    1382             :     {
    1383  1320198601 :       p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
    1384 10528032771 :       for (j=i-dy1; j<=i && j<=dz; j++)
    1385             :       {
    1386  9207834170 :         p1 += z[j]*y[i-j];
    1387  9207834170 :         if (p1 & HIGHBIT) p1 %= p;
    1388             :       }
    1389  1320198601 :       p1 %= p;
    1390  1320198601 :       z[i-dy] = p1? ((p - p1)*inv) % p: 0;
    1391             :     }
    1392  3256809063 :     for (i=0; i<dy; i++)
    1393             :     {
    1394  2783174216 :       p1 = z[0]*y[i];
    1395 14523431176 :       for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
    1396             :       {
    1397 11740256960 :         p1 += z[j]*y[i-j];
    1398 11740256960 :         if (p1 & HIGHBIT) p1 %= p;
    1399             :       }
    1400  2783503635 :       c[i] = Fl_sub(x[i], p1%p, p);
    1401             :     }
    1402             :   }
    1403             :   else
    1404             :   {
    1405   206648824 :     z[dz] = Fl_mul_pre(inv, x[dx], p, pi);
    1406   710966888 :     for (i=dx-1; i>=dy; --i)
    1407             :     {
    1408   504142351 :       p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
    1409  2234561520 :       for (j=i-dy1; j<=i && j<=dz; j++)
    1410  1729836609 :         p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
    1411   504724911 :       z[i-dy] = p1? Fl_mul_pre(p - p1, inv, p, pi): 0;
    1412             :     }
    1413  1519802759 :     for (i=0; i<dy; i++)
    1414             :     {
    1415  1313681262 :       p1 = Fl_mul_pre(z[0],y[i],p,pi);
    1416  3900103217 :       for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
    1417  2580109324 :         p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
    1418  1301990311 :       c[i] = Fl_sub(x[i], p1, p);
    1419             :     }
    1420             :   }
    1421   844865239 :   i = dy-1; while (i>=0 && !c[i]) i--;
    1422   679756344 :   set_avma(av); return Flx_renormalize(c-2, i+3);
    1423             : }
    1424             : 
    1425             : /* as FpX_divrem but working only on ulong types.
    1426             :  * if relevant, *pr is the last object on stack */
    1427             : static GEN
    1428    62219835 : Flx_divrem_basecase(GEN x, GEN y, ulong p, ulong pi, GEN *pr)
    1429             : {
    1430             :   GEN z,q,c;
    1431             :   long dx,dy,dy1,dz,i,j;
    1432             :   ulong p1,inv;
    1433    62219835 :   long sv=x[1];
    1434             : 
    1435    62219835 :   dy = degpol(y);
    1436    62216914 :   if (dy<0) pari_err_INV("Flx_divrem",y);
    1437    62217384 :   if (pr == ONLY_REM) return Flx_rem_basecase(x, y, p, pi);
    1438    62216986 :   if (!dy)
    1439             :   {
    1440     7508750 :     if (pr && pr != ONLY_DIVIDES) *pr = pol0_Flx(sv);
    1441     7508576 :     if (y[2] == 1UL) return Flx_copy(x);
    1442     5316791 :     return Flx_Fl_mul_pre(x, Fl_inv(y[2], p), p, pi);
    1443             :   }
    1444    54708236 :   dx = degpol(x);
    1445    54714732 :   dz = dx-dy;
    1446    54714732 :   if (dz < 0)
    1447             :   {
    1448     1017086 :     q = pol0_Flx(sv);
    1449     1017079 :     if (pr && pr != ONLY_DIVIDES) *pr = Flx_copy(x);
    1450     1017079 :     return q;
    1451             :   }
    1452    53697646 :   x += 2;
    1453    53697646 :   y += 2;
    1454    53697646 :   z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
    1455    53685305 :   inv = uel(y, dy);
    1456    53685305 :   if (inv != 1UL) inv = Fl_inv(inv,p);
    1457    79051534 :   for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
    1458             : 
    1459    53691362 :   if (SMALL_ULONG(p))
    1460             :   {
    1461    51822482 :     z[dz] = (inv*x[dx]) % p;
    1462   131777328 :     for (i=dx-1; i>=dy; --i)
    1463             :     {
    1464    79954846 :       p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
    1465   258002094 :       for (j=i-dy1; j<=i && j<=dz; j++)
    1466             :       {
    1467   178047248 :         p1 += z[j]*y[i-j];
    1468   178047248 :         if (p1 & HIGHBIT) p1 %= p;
    1469             :       }
    1470    79954846 :       p1 %= p;
    1471    79954846 :       z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
    1472             :     }
    1473             :   }
    1474             :   else
    1475             :   {
    1476     1868880 :     z[dz] = Fl_mul(inv, x[dx], p);
    1477     9216166 :     for (i=dx-1; i>=dy; --i)
    1478             :     { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
    1479     7348211 :       p1 = p - uel(x,i);
    1480    26307546 :       for (j=i-dy1; j<=i && j<=dz; j++)
    1481    18959334 :         p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
    1482     7348212 :       z[i-dy] = p1? Fl_mul(p - p1, inv, p): 0;
    1483             :     }
    1484             :   }
    1485    53690437 :   q = Flx_renormalize(z-2, dz+3);
    1486    53697231 :   if (!pr) return q;
    1487             : 
    1488    27192259 :   c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
    1489    27195268 :   if (SMALL_ULONG(p))
    1490             :   {
    1491   233540105 :     for (i=0; i<dy; i++)
    1492             :     {
    1493   207977891 :       p1 = (ulong)z[0]*y[i];
    1494   487393064 :       for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
    1495             :       {
    1496   279415173 :         p1 += (ulong)z[j]*y[i-j];
    1497   279415173 :         if (p1 & HIGHBIT) p1 %= p;
    1498             :       }
    1499   207977558 :       c[i] = Fl_sub(x[i], p1%p, p);
    1500             :     }
    1501             :   }
    1502             :   else
    1503             :   {
    1504    16088277 :     for (i=0; i<dy; i++)
    1505             :     {
    1506    14454920 :       p1 = Fl_mul(z[0],y[i],p);
    1507    50348738 :       for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
    1508    35893815 :         p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
    1509    14454925 :       c[i] = Fl_sub(x[i], p1, p);
    1510             :     }
    1511             :   }
    1512    36588804 :   i=dy-1; while (i>=0 && !c[i]) i--;
    1513    27195571 :   c = Flx_renormalize(c-2, i+3);
    1514    27195406 :   if (pr == ONLY_DIVIDES)
    1515         417 :   { if (lg(c) != 2) return NULL; }
    1516             :   else
    1517    27194989 :     *pr = c;
    1518    27195266 :   return q;
    1519             : }
    1520             : 
    1521             : /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
    1522             :  * and mg is the Barrett inverse of T. */
    1523             : static GEN
    1524      904887 : Flx_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
    1525             : {
    1526             :   GEN q, r;
    1527      904887 :   long lt = degpol(T); /*We discard the leading term*/
    1528             :   long ld, lm, lT, lmg;
    1529      904839 :   ld = l-lt;
    1530      904839 :   lm = minss(ld, lgpol(mg));
    1531      905144 :   lT  = Flx_lgrenormalizespec(T+2,lt);
    1532      905232 :   lmg = Flx_lgrenormalizespec(mg+2,lm);
    1533      905081 :   q = Flx_recipspec(x+lt,ld,ld);               /* q = rec(x)      lz<=ld*/
    1534      904380 :   q = Flx_mulspec(q+2,mg+2,p,pi,lgpol(q),lmg); /* q = rec(x) * mg lz<=ld+lm*/
    1535      905242 :   q = Flx_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lz<=ld*/
    1536      904392 :   if (!pr) return q;
    1537      896704 :   r = Flx_mulspec(q+2,T+2,p,pi,lgpol(q),lT);   /* r = q*pol      lz<=ld+lt*/
    1538      897560 :   r = Flx_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - q*pol lz<=lt */
    1539      897145 :   if (pr == ONLY_REM) return r;
    1540      424965 :   *pr = r; return q;
    1541             : }
    1542             : 
    1543             : static GEN
    1544      607316 : Flx_divrem_Barrett(GEN x, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
    1545             : {
    1546      607316 :   GEN q = NULL, r = Flx_copy(x);
    1547      607364 :   long l = lgpol(x), lt = degpol(T), lm = 2*lt-1, v = T[1];
    1548             :   long i;
    1549      607357 :   if (l <= lt)
    1550             :   {
    1551           0 :     if (pr == ONLY_REM) return Flx_copy(x);
    1552           0 :     if (pr == ONLY_DIVIDES) return lgpol(x)? NULL: pol0_Flx(v);
    1553           0 :     if (pr) *pr = Flx_copy(x);
    1554           0 :     return pol0_Flx(v);
    1555             :   }
    1556      607357 :   if (lt <= 1)
    1557        1390 :     return Flx_divrem_basecase(x,T,p,pi,pr);
    1558      605967 :   if (pr != ONLY_REM && l>lm)
    1559       28889 :   { q = zero_zv(l-lt+1); q[1] = T[1]; }
    1560      906393 :   while (l>lm)
    1561             :   {
    1562      300558 :     GEN zr, zq = Flx_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,pi,&zr);
    1563      300475 :     long lz = lgpol(zr);
    1564      300426 :     if (pr != ONLY_REM)
    1565             :     {
    1566       57929 :       long lq = lgpol(zq);
    1567      872174 :       for(i=0; i<lq; i++) q[2+l-lm+i] = zq[2+i];
    1568             :     }
    1569     4369683 :     for(i=0; i<lz; i++)   r[2+l-lm+i] = zr[2+i];
    1570      300426 :     l = l-lm+lz;
    1571             :   }
    1572      605835 :   if (pr == ONLY_REM)
    1573             :   {
    1574      472232 :     if (l > lt)
    1575      472190 :       r = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi,ONLY_REM);
    1576             :     else
    1577          42 :       r = Flx_renormalize(r, l+2);
    1578      472216 :     r[1] = v; return r;
    1579             :   }
    1580      133603 :   if (l > lt)
    1581             :   {
    1582      132159 :     GEN zq = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi, pr ? &r: NULL);
    1583      132159 :     if (!q) q = zq;
    1584             :     else
    1585             :     {
    1586       27315 :       long lq = lgpol(zq);
    1587      158603 :       for(i=0; i<lq; i++) q[2+i] = zq[2+i];
    1588             :     }
    1589             :   }
    1590        1444 :   else if (pr)
    1591        1535 :     r = Flx_renormalize(r, l+2);
    1592      133603 :   q[1] = v; q = Flx_renormalize(q, lg(q));
    1593      133733 :   if (pr == ONLY_DIVIDES) return lgpol(r)? NULL: q;
    1594      133733 :   if (pr) { r[1] = v; *pr = r; }
    1595      133733 :   return q;
    1596             : }
    1597             : 
    1598             : /* allow pi = 0 (SMALL_ULONG) */
    1599             : GEN
    1600    79611669 : Flx_divrem_pre(GEN x, GEN T, ulong p, ulong pi, GEN *pr)
    1601             : {
    1602             :   GEN B, y;
    1603             :   long dy, dx, d;
    1604    79611669 :   if (pr==ONLY_REM) return Flx_rem_pre(x, T, p, pi);
    1605    62343082 :   y = get_Flx_red(T, &B);
    1606    62363942 :   dy = degpol(y); dx = degpol(x); d = dx-dy;
    1607    62353789 :   if (!B && d+3 < get_Fl_threshold(p, Flx_DIVREM_BARRETT_LIMIT,Flx_DIVREM2_BARRETT_LIMIT))
    1608    62217516 :     return Flx_divrem_basecase(x,y,p,pi,pr);
    1609             :   else
    1610             :   {
    1611      134725 :     pari_sp av = avma;
    1612      134725 :     GEN mg = B? B: Flx_invBarrett_pre(y, p, pi);
    1613      134725 :     GEN q1 = Flx_divrem_Barrett(x,mg,y,p,pi,pr);
    1614      134725 :     if (!q1) return gc_NULL(av);
    1615      134725 :     if (!pr || pr==ONLY_DIVIDES) return gerepileuptoleaf(av, q1);
    1616      126437 :     return gc_all(av, 2, &q1, pr);
    1617             :   }
    1618             : }
    1619             : GEN
    1620    29985937 : Flx_divrem(GEN x, GEN T, ulong p, GEN *pr)
    1621    29985937 : { return Flx_divrem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p), pr); }
    1622             : 
    1623             : GEN
    1624   834807486 : Flx_rem_pre(GEN x, GEN T, ulong p, ulong pi)
    1625             : {
    1626   834807486 :   GEN B, y = get_Flx_red(T, &B);
    1627   834784229 :   long d = degpol(x) - degpol(y);
    1628   834641106 :   if (d < 0) return Flx_copy(x);
    1629   712744540 :   if (!B && d+3 < get_Fl_threshold(p, Flx_REM_BARRETT_LIMIT,Flx_REM2_BARRETT_LIMIT))
    1630   712039834 :     return Flx_rem_basecase(x,y,p, pi);
    1631             :   else
    1632             :   {
    1633      472595 :     pari_sp av=avma;
    1634      472595 :     GEN mg = B ? B: Flx_invBarrett_pre(y, p, pi);
    1635      472594 :     GEN r  = Flx_divrem_Barrett(x, mg, y, p, pi, ONLY_REM);
    1636      472612 :     return gerepileuptoleaf(av, r);
    1637             :   }
    1638             : }
    1639             : GEN
    1640    41006460 : Flx_rem(GEN x, GEN T, ulong p)
    1641    41006460 : { return Flx_rem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    1642             : 
    1643             : /* reduce T mod (X^n - 1, p). Shallow function */
    1644             : GEN
    1645     5090608 : Flx_mod_Xnm1(GEN T, ulong n, ulong p)
    1646             : {
    1647     5090608 :   long i, j, L = lg(T), l = n+2;
    1648             :   GEN S;
    1649     5090608 :   if (L <= l || n & ~LGBITS) return T;
    1650        3361 :   S = cgetg(l, t_VECSMALL);
    1651        3361 :   S[1] = T[1];
    1652       13657 :   for (i = 2; i < l; i++) S[i] = T[i];
    1653        9198 :   for (j = 2; i < L; i++) {
    1654        5837 :     S[j] = Fl_add(S[j], T[i], p);
    1655        5837 :     if (++j == l) j = 2;
    1656             :   }
    1657        3361 :   return Flx_renormalize(S, l);
    1658             : }
    1659             : /* reduce T mod (X^n + 1, p). Shallow function */
    1660             : GEN
    1661       29713 : Flx_mod_Xn1(GEN T, ulong n, ulong p)
    1662             : {
    1663       29713 :   long i, j, L = lg(T), l = n+2;
    1664             :   GEN S;
    1665       29713 :   if (L <= l || n & ~LGBITS) return T;
    1666        2593 :   S = cgetg(l, t_VECSMALL);
    1667        2593 :   S[1] = T[1];
    1668       10991 :   for (i = 2; i < l; i++) S[i] = T[i];
    1669        6752 :   for (j = 2; i < L; i++) {
    1670        4159 :     S[j] = Fl_sub(S[j], T[i], p);
    1671        4159 :     if (++j == l) j = 2;
    1672             :   }
    1673        2593 :   return Flx_renormalize(S, l);
    1674             : }
    1675             : 
    1676             : struct _Flxq {
    1677             :   GEN aut, T;
    1678             :   ulong p, pi;
    1679             : };
    1680             : /* allow pi = 0 */
    1681             : static void
    1682    69602003 : set_Flxq_pre(struct _Flxq *D, GEN T, ulong p, ulong pi)
    1683             : {
    1684    69602003 :   D->p = p;
    1685    69602003 :   D->pi = pi;
    1686    69602003 :   D->T = Flx_get_red_pre(T, p, pi);
    1687    69608423 : }
    1688             : static void
    1689       73048 : set_Flxq(struct _Flxq *D, GEN T, ulong p)
    1690       73048 : { set_Flxq_pre(D, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    1691             : 
    1692             : static GEN
    1693           0 : _Flx_divrem(void * E, GEN x, GEN y, GEN *r)
    1694             : {
    1695           0 :   struct _Flxq *D = (struct _Flxq*) E;
    1696           0 :   return Flx_divrem_pre(x, y, D->p, D->pi, r);
    1697             : }
    1698             : static GEN
    1699      578248 : _Flx_add(void * E, GEN x, GEN y) {
    1700      578248 :   struct _Flxq *D = (struct _Flxq*) E;
    1701      578248 :   return Flx_add(x, y, D->p);
    1702             : }
    1703             : static GEN
    1704    10446397 : _Flx_mul(void *E, GEN x, GEN y) {
    1705    10446397 :   struct _Flxq *D = (struct _Flxq*) E;
    1706    10446397 :   return Flx_mul_pre(x, y, D->p, D->pi);
    1707             : }
    1708             : static GEN
    1709           0 : _Flx_sqr(void *E, GEN x) {
    1710           0 :   struct _Flxq *D = (struct _Flxq*) E;
    1711           0 :   return Flx_sqr_pre(x, D->p, D->pi);
    1712             : }
    1713             : 
    1714             : static struct bb_ring Flx_ring = { _Flx_add,_Flx_mul,_Flx_sqr };
    1715             : 
    1716             : GEN
    1717           0 : Flx_digits(GEN x, GEN T, ulong p)
    1718             : {
    1719             :   struct _Flxq D;
    1720           0 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
    1721           0 :   D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    1722           0 :   return gen_digits(x,T,n,(void *)&D, &Flx_ring, _Flx_divrem);
    1723             : }
    1724             : 
    1725             : GEN
    1726           0 : FlxV_Flx_fromdigits(GEN x, GEN T, ulong p)
    1727             : {
    1728             :   struct _Flxq D;
    1729           0 :   D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    1730           0 :   return gen_fromdigits(x,T,(void *)&D, &Flx_ring);
    1731             : }
    1732             : 
    1733             : long
    1734     3922426 : Flx_val(GEN x)
    1735             : {
    1736     3922426 :   long i, l=lg(x);
    1737     3922426 :   if (l==2)  return LONG_MAX;
    1738     3931286 :   for (i=2; i<l && x[i]==0; i++) /*empty*/;
    1739     3922426 :   return i-2;
    1740             : }
    1741             : long
    1742    26236379 : Flx_valrem(GEN x, GEN *Z)
    1743             : {
    1744    26236379 :   long v, i, l=lg(x);
    1745             :   GEN y;
    1746    26236379 :   if (l==2) { *Z = Flx_copy(x); return LONG_MAX; }
    1747    28402486 :   for (i=2; i<l && x[i]==0; i++) /*empty*/;
    1748    26236379 :   v = i-2;
    1749    26236379 :   if (v == 0) { *Z = x; return 0; }
    1750     1016643 :   l -= v;
    1751     1016643 :   y = cgetg(l, t_VECSMALL); y[1] = x[1];
    1752     2607082 :   for (i=2; i<l; i++) y[i] = x[i+v];
    1753     1016508 :   *Z = y; return v;
    1754             : }
    1755             : 
    1756             : GEN
    1757    18508960 : Flx_deriv(GEN z, ulong p)
    1758             : {
    1759    18508960 :   long i,l = lg(z)-1;
    1760             :   GEN x;
    1761    18508960 :   if (l < 2) l = 2;
    1762    18508960 :   x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
    1763    18500817 :   if (HIGHWORD(l | p))
    1764    53710142 :     for (i=2; i<l; i++) x[i] = Fl_mul((ulong)i-1, z[i], p);
    1765             :   else
    1766    75344125 :     for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
    1767    18501822 :   return Flx_renormalize(x,l);
    1768             : }
    1769             : 
    1770             : static GEN
    1771      422659 : Flx_integXn(GEN x, long n, ulong p)
    1772             : {
    1773      422659 :   long i, lx = lg(x);
    1774             :   GEN y;
    1775      422659 :   if (lx == 2) return Flx_copy(x);
    1776      412846 :   y = cgetg(lx, t_VECSMALL); y[1] = x[1];
    1777     2096458 :   for (i=2; i<lx; i++)
    1778             :   {
    1779     1683070 :     ulong xi = uel(x,i);
    1780     1683070 :     if (xi == 0)
    1781       13345 :       uel(y,i) = 0;
    1782             :     else
    1783             :     {
    1784     1669725 :       ulong j = n+i-1;
    1785     1669725 :       ulong d = ugcd(j, xi);
    1786     1669623 :       if (d==1)
    1787     1018170 :         uel(y,i) = Fl_div(xi, j, p);
    1788             :       else
    1789      651453 :         uel(y,i) = Fl_div(xi/d, j/d, p);
    1790             :     }
    1791             :   }
    1792      413388 :   return Flx_renormalize(y, lx);;
    1793             : }
    1794             : 
    1795             : GEN
    1796           0 : Flx_integ(GEN x, ulong p)
    1797             : {
    1798           0 :   long i, lx = lg(x);
    1799             :   GEN y;
    1800           0 :   if (lx == 2) return Flx_copy(x);
    1801           0 :   y = cgetg(lx+1, t_VECSMALL); y[1] = x[1];
    1802           0 :   uel(y,2) = 0;
    1803           0 :   for (i=3; i<=lx; i++)
    1804           0 :     uel(y,i) = uel(x,i-1) ? Fl_div(uel(x,i-1), (i-2)%p, p): 0UL;
    1805           0 :   return Flx_renormalize(y, lx+1);;
    1806             : }
    1807             : 
    1808             : /* assume p prime */
    1809             : GEN
    1810       13447 : Flx_diff1(GEN P, ulong p)
    1811             : {
    1812       13447 :   return Flx_sub(Flx_translate1(P, p), P, p);
    1813             : }
    1814             : 
    1815             : GEN
    1816      343993 : Flx_deflate(GEN x0, long d)
    1817             : {
    1818             :   GEN z, y, x;
    1819      343993 :   long i,id, dy, dx = degpol(x0);
    1820      343993 :   if (d == 1 || dx <= 0) return Flx_copy(x0);
    1821      280853 :   dy = dx/d;
    1822      280853 :   y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
    1823      280850 :   z = y + 2;
    1824      280850 :   x = x0+ 2;
    1825      929260 :   for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
    1826      280850 :   return y;
    1827             : }
    1828             : 
    1829             : GEN
    1830      163932 : Flx_inflate(GEN x0, long d)
    1831             : {
    1832      163932 :   long i, id, dy, dx = degpol(x0);
    1833      163934 :   GEN x = x0 + 2, z, y;
    1834      163934 :   if (dx <= 0) return Flx_copy(x0);
    1835      162795 :   dy = dx*d;
    1836      162795 :   y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
    1837      162759 :   z = y + 2;
    1838     9228345 :   for (i=0; i<=dy; i++) z[i] = 0;
    1839     4496944 :   for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
    1840      162759 :   return y;
    1841             : }
    1842             : 
    1843             : /* write p(X) = a_0(X^k) + X*a_1(X^k) + ... + X^(k-1)*a_{k-1}(X^k) */
    1844             : GEN
    1845      148911 : Flx_splitting(GEN p, long k)
    1846             : {
    1847      148911 :   long n = degpol(p), v = p[1], m, i, j, l;
    1848             :   GEN r;
    1849             : 
    1850      148909 :   m = n/k;
    1851      148909 :   r = cgetg(k+1,t_VEC);
    1852      684015 :   for(i=1; i<=k; i++)
    1853             :   {
    1854      535108 :     gel(r,i) = cgetg(m+3, t_VECSMALL);
    1855      535101 :     mael(r,i,1) = v;
    1856             :   }
    1857     4659463 :   for (j=1, i=0, l=2; i<=n; i++)
    1858             :   {
    1859     4510556 :     mael(r,j,l) = p[2+i];
    1860     4510556 :     if (j==k) { j=1; l++; } else j++;
    1861             :   }
    1862      684030 :   for(i=1; i<=k; i++)
    1863      535145 :     gel(r,i) = Flx_renormalize(gel(r,i),i<j?l+1:l);
    1864      148885 :   return r;
    1865             : }
    1866             : 
    1867             : /* ux + vy */
    1868             : static GEN
    1869      431285 : Flx_addmulmul(GEN u, GEN v, GEN x, GEN y, ulong p, ulong pi)
    1870      431285 : { return Flx_add(Flx_mul_pre(u,x, p,pi), Flx_mul_pre(v,y, p,pi), p); }
    1871             : 
    1872             : static GEN
    1873       25344 : FlxM_Flx_mul2(GEN M, GEN x, GEN y, ulong p, ulong pi)
    1874             : {
    1875       25344 :   GEN res = cgetg(3, t_COL);
    1876       25343 :   gel(res, 1) = Flx_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p, pi);
    1877       25344 :   gel(res, 2) = Flx_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p, pi);
    1878       25344 :   return res;
    1879             : }
    1880             : 
    1881             : #if 0
    1882             : static GEN
    1883             : FlxM_mul2_old(GEN M, GEN N, ulong p)
    1884             : {
    1885             :   GEN res = cgetg(3, t_MAT);
    1886             :   gel(res, 1) = FlxM_Flx_mul2(M,gcoeff(N,1,1),gcoeff(N,2,1),p);
    1887             :   gel(res, 2) = FlxM_Flx_mul2(M,gcoeff(N,1,2),gcoeff(N,2,2),p);
    1888             :   return res;
    1889             : }
    1890             : #endif
    1891             : /* A,B are 2x2 matrices, Flx entries. Return A x B using Strassen 7M formula */
    1892             : static GEN
    1893        6813 : FlxM_mul2(GEN A, GEN B, ulong p, ulong pi)
    1894             : {
    1895        6813 :   GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
    1896        6813 :   GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
    1897        6813 :   GEN M1 = Flx_mul_pre(Flx_add(A11,A22, p), Flx_add(B11,B22, p), p, pi);
    1898        6813 :   GEN M2 = Flx_mul_pre(Flx_add(A21,A22, p), B11, p, pi);
    1899        6813 :   GEN M3 = Flx_mul_pre(A11, Flx_sub(B12,B22, p), p, pi);
    1900        6813 :   GEN M4 = Flx_mul_pre(A22, Flx_sub(B21,B11, p), p, pi);
    1901        6813 :   GEN M5 = Flx_mul_pre(Flx_add(A11,A12, p), B22, p, pi);
    1902        6813 :   GEN M6 = Flx_mul_pre(Flx_sub(A21,A11, p), Flx_add(B11,B12, p), p, pi);
    1903        6813 :   GEN M7 = Flx_mul_pre(Flx_sub(A12,A22, p), Flx_add(B21,B22, p), p, pi);
    1904        6813 :   GEN T1 = Flx_add(M1,M4, p), T2 = Flx_sub(M7,M5, p);
    1905        6813 :   GEN T3 = Flx_sub(M1,M2, p), T4 = Flx_add(M3,M6, p);
    1906        6813 :   retmkmat22(Flx_add(T1,T2, p), Flx_add(M3,M5, p),
    1907             :              Flx_add(M2,M4, p), Flx_add(T3,T4, p));
    1908             : }
    1909             : 
    1910             : /* Return [0,1;1,-q]*M */
    1911             : static GEN
    1912        6641 : Flx_FlxM_qmul(GEN q, GEN M, ulong p, ulong pi)
    1913             : {
    1914        6641 :   GEN u = Flx_mul_pre(gcoeff(M,2,1), q, p, pi);
    1915        6641 :   GEN v = Flx_mul_pre(gcoeff(M,2,2), q, p, pi);
    1916        6641 :   retmkmat22(gcoeff(M,2,1), gcoeff(M,2,2),
    1917             :     Flx_sub(gcoeff(M,1,1), u, p), Flx_sub(gcoeff(M,1,2), v, p));
    1918             : }
    1919             : 
    1920             : static GEN
    1921         895 : matid2_FlxM(long v)
    1922         895 : { retmkmat22(pol1_Flx(v),pol0_Flx(v),pol0_Flx(v),pol1_Flx(v)); }
    1923             : 
    1924             : static GEN
    1925          13 : matJ2_FlxM(long v)
    1926          13 : { retmkmat22(pol0_Flx(v),pol1_Flx(v),pol1_Flx(v),pol0_Flx(v)); }
    1927             : 
    1928             : struct Flx_res
    1929             : {
    1930             :    ulong res, lc;
    1931             :    long deg0, deg1, off;
    1932             : };
    1933             : 
    1934             : INLINE void
    1935        9405 : Flx_halfres_update_pre(long da, long db, long dr, ulong p, ulong pi, struct Flx_res *res)
    1936             : {
    1937        9405 :   if (dr >= 0)
    1938             :   {
    1939        9405 :     if (res->lc != 1)
    1940             :     {
    1941        7596 :       if (pi)
    1942             :       {
    1943        3127 :         res->lc  = Fl_powu_pre(res->lc, da - dr, p, pi);
    1944        3127 :         res->res = Fl_mul_pre(res->res, res->lc, p, pi);
    1945             :       } else
    1946             :       {
    1947        4469 :         res->lc  = Fl_powu(res->lc, da - dr, p);
    1948        4469 :         res->res = Fl_mul(res->res, res->lc, p);
    1949             :       }
    1950             :     }
    1951        9405 :     if (both_odd(da + res->off, db + res->off))
    1952          63 :       res->res = Fl_neg(res->res, p);
    1953             :   } else
    1954             :   {
    1955           0 :     if (db == 0)
    1956             :     {
    1957           0 :       if (res->lc != 1)
    1958             :       {
    1959           0 :         if (pi)
    1960             :         {
    1961           0 :           res->lc  = Fl_powu_pre(res->lc, da, p, pi);
    1962           0 :           res->res = Fl_mul_pre(res->res, res->lc, p, pi);
    1963             :         } else
    1964             :         {
    1965           0 :           res->lc  = Fl_powu(res->lc, da, p);
    1966           0 :           res->res = Fl_mul(res->res, res->lc, p);
    1967             :         }
    1968             :       }
    1969             :     } else
    1970           0 :       res->res = 0;
    1971             :   }
    1972        9405 : }
    1973             : 
    1974             : static GEN
    1975     1115393 : Flx_halfres_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *pa, GEN *pb, struct Flx_res *res)
    1976             : {
    1977     1115393 :   pari_sp av = avma;
    1978             :   GEN u, u1, v, v1, M;
    1979     1115393 :   long vx = a[1], n = lgpol(a)>>1;
    1980     1115390 :   u1 = v = pol0_Flx(vx);
    1981     1115376 :   u = v1 = pol1_Flx(vx);
    1982     6992023 :   while (lgpol(b)>n)
    1983             :   {
    1984             :     GEN r, q;
    1985     5876659 :     q = Flx_divrem_pre(a,b,p,pi, &r);
    1986     5876666 :     if (res)
    1987             :     {
    1988        8362 :       long da = degpol(a), db=degpol(b), dr = degpol(r);
    1989        8362 :       res->lc = b[db+2];
    1990        8362 :       if (dr >= n)
    1991        7133 :         Flx_halfres_update_pre(da, db, dr, p, pi, res);
    1992             :       else
    1993             :       {
    1994        1229 :         res->deg0 = da;
    1995        1229 :         res->deg1 = db;
    1996             :       }
    1997             :     }
    1998     5876666 :     a = b; b = r; swap(u,u1); swap(v,v1);
    1999     5876666 :     u1 = Flx_sub(u1, Flx_mul(u, q, p), p);
    2000     5876592 :     v1 = Flx_sub(v1, Flx_mul(v, q, p), p);
    2001     5876642 :     if (gc_needed(av,2))
    2002             :     {
    2003           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_halfgcd (d = %ld)",degpol(b));
    2004           0 :       gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
    2005             :     }
    2006             :   }
    2007     1115201 :   M = mkmat22(u,v,u1,v1); *pa = a; *pb = b;
    2008     1115361 :   return gc_all(av,3, &M, pa, pb);
    2009             : }
    2010             : 
    2011             : static GEN Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res);
    2012             : 
    2013             : static GEN
    2014       19578 : Flx_halfres_split(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
    2015             : {
    2016       19578 :   pari_sp av = avma;
    2017             :   GEN R, S, T, V1, V2;
    2018             :   GEN x1, y1, r, q;
    2019       19578 :   long l = lgpol(x), n = l>>1, k;
    2020       19578 :   if (lgpol(y) <= n)
    2021         855 :     { *a = Flx_copy(x); *b = Flx_copy(y); return matid2_FlxM(x[1]); }
    2022       18722 :   if (res)
    2023             :   {
    2024        3262 :      res->lc = Flx_lead(y);
    2025        3262 :      res->deg0 -= n;
    2026        3262 :      res->deg1 -= n;
    2027        3262 :      res->off += n;
    2028             :   }
    2029       18722 :   R = Flx_halfres_i(Flx_shift(x,-n),Flx_shift(y,-n),p,pi,a,b,res);
    2030       18722 :   if (res)
    2031             :   {
    2032        3262 :     res->off -= n;
    2033        3262 :     res->deg0 += n;
    2034        3262 :     res->deg1 += n;
    2035             :   }
    2036       18722 :   V1 = FlxM_Flx_mul2(R, Flxn_red(x,n), Flxn_red(y,n), p, pi);
    2037       18723 :   x1 = Flx_add(Flx_shift(*a,n), gel(V1,1), p);
    2038       18723 :   y1 = Flx_add(Flx_shift(*b,n), gel(V1,2), p);
    2039       18723 :   if (lgpol(y1) <= n)
    2040       12102 :     { *a = x1; *b = y1; return gc_all(av, 3, &R, a, b); }
    2041        6621 :   k = 2*n-degpol(y1);
    2042        6621 :   q = Flx_divrem_pre(x1, y1, p, pi, &r);
    2043        6621 :   if (res)
    2044             :   {
    2045        1043 :     long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
    2046        1043 :     if (dy1 < degpol(y))
    2047         185 :       Flx_halfres_update_pre(res->deg0, res->deg1, dy1, p, pi, res);
    2048        1043 :     res->lc = uel(y1, dy1+2);
    2049        1043 :     res->deg0 = dx1;
    2050        1043 :     res->deg1 = dy1;
    2051        1043 :     if (dr >= n)
    2052             :     {
    2053        1043 :       Flx_halfres_update_pre(dx1, dy1, dr, p, pi, res);
    2054        1043 :       res->deg0 = dy1;
    2055        1043 :       res->deg1 = dr;
    2056             :     }
    2057        1043 :     res->deg0 -= k;
    2058        1043 :     res->deg1 -= k;
    2059        1043 :     res->off += k;
    2060             :   }
    2061        6621 :   S = Flx_halfres_i(Flx_shift(y1,-k), Flx_shift(r,-k), p, pi, a, b, res);
    2062        6621 :   if (res)
    2063             :   {
    2064        1043 :     res->deg0 += k;
    2065        1043 :     res->deg1 += k;
    2066        1043 :     res->off -= k;
    2067             :   }
    2068        6621 :   T = FlxM_mul2(S, Flx_FlxM_qmul(q, R, p,pi), p, pi);
    2069        6621 :   V2 = FlxM_Flx_mul2(S, Flxn_red(y1,k), Flxn_red(r,k), p, pi);
    2070        6621 :   *a = Flx_add(Flx_shift(*a,k), gel(V2,1), p);
    2071        6621 :   *b = Flx_add(Flx_shift(*b,k), gel(V2,2), p);
    2072        6621 :   return gc_all(av, 3, &T, a, b);
    2073             : }
    2074             : 
    2075             : static GEN
    2076     1134974 : Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
    2077             : {
    2078     1134974 :   if (lgpol(x) < get_Fl_threshold(p, Flx_HALFGCD_LIMIT, Flx_HALFGCD2_LIMIT))
    2079     1115393 :     return Flx_halfres_basecase(x, y, p, pi, a, b, res);
    2080       19578 :   return Flx_halfres_split(x, y, p, pi, a, b, res);
    2081             : }
    2082             : 
    2083             : static GEN
    2084     1108587 : Flx_halfgcd_all_i(GEN x, GEN y, ulong p, ulong pi, GEN *pa, GEN *pb)
    2085             : {
    2086             :   GEN a, b, R;
    2087     1108587 :   R = Flx_halfres_i(x, y, p, pi, &a, &b, NULL);
    2088     1108587 :   if (pa) *pa = a;
    2089     1108587 :   if (pb) *pb = b;
    2090     1108587 :   return R;
    2091             : }
    2092             : 
    2093             : /* Return M in GL_2(Fl[X]) such that:
    2094             : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
    2095             : */
    2096             : 
    2097             : GEN
    2098     1108592 : Flx_halfgcd_all_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b)
    2099             : {
    2100             :   pari_sp av;
    2101             :   GEN R, q ,r;
    2102     1108592 :   long lx = lgpol(x), ly = lgpol(y);
    2103     1108587 :   if (!lx)
    2104             :   {
    2105           0 :     if (a) *a = Flx_copy(y);
    2106           0 :     if (b) *b = Flx_copy(x);
    2107           0 :     return matJ2_FlxM(x[1]);
    2108             :   }
    2109     1108587 :   if (ly < lx) return Flx_halfgcd_all_i(x, y, p, pi, a, b);
    2110       10779 :   av = avma;
    2111       10779 :   q = Flx_divrem(y,x,p,&r);
    2112       10779 :   R = Flx_halfgcd_all_i(x, r, p, pi, a, b);
    2113       10779 :   gcoeff(R,1,1) = Flx_sub(gcoeff(R,1,1), Flx_mul_pre(q,gcoeff(R,1,2), p,pi), p);
    2114       10779 :   gcoeff(R,2,1) = Flx_sub(gcoeff(R,2,1), Flx_mul_pre(q,gcoeff(R,2,2), p,pi), p);
    2115       10779 :   return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
    2116             : }
    2117             : 
    2118             : GEN
    2119         126 : Flx_halfgcd_all(GEN x, GEN y, ulong p, GEN *a, GEN *b)
    2120         126 : { return Flx_halfgcd_all_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), a, b); }
    2121             : 
    2122             : GEN
    2123      841435 : Flx_halfgcd_pre(GEN x, GEN y, ulong p, ulong pi)
    2124      841435 : { return Flx_halfgcd_all_pre(x, y, p, pi, NULL, NULL); }
    2125             : 
    2126             : GEN
    2127           0 : Flx_halfgcd(GEN x, GEN y, ulong p)
    2128           0 : { return Flx_halfgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2129             : 
    2130             : /*Do not garbage collect*/
    2131             : static GEN
    2132    76752700 : Flx_gcd_basecase(GEN a, GEN b, ulong p, ulong pi)
    2133             : {
    2134    76752700 :   pari_sp av = avma;
    2135    76752700 :   ulong iter = 0;
    2136    76752700 :   if (lg(b) > lg(a)) swap(a, b);
    2137   266552726 :   while (lgpol(b))
    2138             :   {
    2139   189156096 :     GEN c = Flx_rem_pre(a,b,p,pi);
    2140   189800026 :     iter++; a = b; b = c;
    2141   189800026 :     if (gc_needed(av,2))
    2142             :     {
    2143           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (d = %ld)",degpol(c));
    2144           0 :       gerepileall(av,2, &a,&b);
    2145             :     }
    2146             :   }
    2147    76682279 :   return iter < 2 ? Flx_copy(a) : a;
    2148             : }
    2149             : 
    2150             : GEN
    2151    78363315 : Flx_gcd_pre(GEN x, GEN y, ulong p, ulong pi)
    2152             : {
    2153    78363315 :   pari_sp av = avma;
    2154             :   long lim;
    2155    78363315 :   if (!lgpol(x)) return Flx_copy(y);
    2156    76756233 :   lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
    2157    76767858 :   while (lgpol(y) >= lim)
    2158             :   {
    2159         150 :     if (lgpol(y)<=(lgpol(x)>>1))
    2160             :     {
    2161           0 :       GEN r = Flx_rem_pre(x, y, p, pi);
    2162           0 :       x = y; y = r;
    2163             :     }
    2164         150 :     (void) Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
    2165         150 :     if (gc_needed(av,2))
    2166             :     {
    2167           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (y = %ld)",degpol(y));
    2168           0 :       gerepileall(av,2,&x,&y);
    2169             :     }
    2170             :   }
    2171    76748948 :   return gerepileuptoleaf(av, Flx_gcd_basecase(x,y,p,pi));
    2172             : }
    2173             : GEN
    2174    27655217 : Flx_gcd(GEN x, GEN y, ulong p)
    2175    27655217 : { return Flx_gcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2176             : 
    2177             : int
    2178     6285168 : Flx_is_squarefree(GEN z, ulong p)
    2179             : {
    2180     6285168 :   pari_sp av = avma;
    2181     6285168 :   GEN d = Flx_gcd(z, Flx_deriv(z,p) , p);
    2182     6285161 :   return gc_bool(av, degpol(d) == 0);
    2183             : }
    2184             : 
    2185             : static long
    2186      139846 : Flx_is_smooth_squarefree(GEN f, long r, ulong p, ulong pi)
    2187             : {
    2188      139846 :   pari_sp av = avma;
    2189             :   long i;
    2190      139846 :   GEN sx = polx_Flx(f[1]), a = sx;
    2191      595403 :   for(i=1;;i++)
    2192             :   {
    2193      595403 :     if (degpol(f)<=r) return gc_long(av,1);
    2194      569105 :     a = Flxq_powu_pre(Flx_rem_pre(a,f,p,pi), p, f, p, pi);
    2195      573117 :     if (Flx_equal(a, sx)) return gc_long(av,1);
    2196      568777 :     if (i==r) return gc_long(av,0);
    2197      456947 :     f = Flx_div_pre(f, Flx_gcd_pre(Flx_sub(a,sx,p),f,p,pi),p,pi);
    2198             :   }
    2199             : }
    2200             : 
    2201             : static long
    2202        9187 : Flx_is_l_pow(GEN x, ulong p)
    2203             : {
    2204        9187 :   ulong i, lx = lgpol(x);
    2205       18496 :   for (i=1; i<lx; i++)
    2206       16519 :     if (x[i+2] && i%p) return 0;
    2207        1977 :   return 1;
    2208             : }
    2209             : 
    2210             : int
    2211      139819 : Flx_is_smooth_pre(GEN g, long r, ulong p, ulong pi)
    2212             : {
    2213             :   while (1)
    2214        9188 :   {
    2215      139819 :     GEN f = Flx_gcd_pre(g, Flx_deriv(g, p), p, pi);
    2216      139734 :     if (!Flx_is_smooth_squarefree(Flx_div_pre(g, f, p, pi), r, p, pi))
    2217      111801 :       return 0;
    2218       28149 :     if (degpol(f)==0) return 1;
    2219        9167 :     g = Flx_is_l_pow(f,p) ? Flx_deflate(f, p): f;
    2220             :   }
    2221             : }
    2222             : int
    2223       74256 : Flx_is_smooth(GEN g, long r, ulong p)
    2224       74256 : { return Flx_is_smooth_pre(g, r, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2225             : 
    2226             : static GEN
    2227     6619243 : Flx_extgcd_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptu, GEN *ptv)
    2228             : {
    2229     6619243 :   pari_sp av=avma;
    2230             :   GEN u,v,u1,v1;
    2231     6619243 :   long vx = a[1];
    2232     6619243 :   v = pol0_Flx(vx); v1 = pol1_Flx(vx);
    2233     6618962 :   if (ptu) { u = pol1_Flx(vx); u1 = pol0_Flx(vx); }
    2234    29277864 :   while (lgpol(b))
    2235             :   {
    2236    22656950 :     GEN r, q = Flx_divrem_pre(a,b,p,pi, &r);
    2237    22658665 :     a = b; b = r;
    2238    22658665 :     if (ptu)
    2239             :     {
    2240     2455915 :       swap(u,u1);
    2241     2455915 :       u1 = Flx_sub(u1, Flx_mul_pre(u, q, p, pi), p);
    2242             :     }
    2243    22658660 :     swap(v,v1);
    2244    22658660 :     v1 = Flx_sub(v1, Flx_mul_pre(v, q, p, pi), p);
    2245    22658919 :     if (gc_needed(av,2))
    2246             :     {
    2247           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_extgcd (d = %ld)",degpol(a));
    2248           0 :       gerepileall(av,ptu ? 6: 4, &a,&b,&v,&v1,&u,&u1);
    2249             :     }
    2250             :   }
    2251     6619089 :   if (ptu) *ptu = u;
    2252     6619089 :   *ptv = v;
    2253     6619089 :   return a;
    2254             : }
    2255             : 
    2256             : static GEN
    2257      151994 : Flx_extgcd_halfgcd(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
    2258             : {
    2259             :   GEN u, v;
    2260      151994 :   long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
    2261      151994 :   GEN V = cgetg(expu(lgpol(y))+2,t_VEC);
    2262      151994 :   long i, n = 0, vs = x[1];
    2263      414967 :   while (lgpol(y) >= lim)
    2264             :   {
    2265      262973 :     if (lgpol(y)<=(lgpol(x)>>1))
    2266             :     {
    2267          26 :       GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
    2268          26 :       x = y; y = r;
    2269          26 :       gel(V,++n) = mkmat22(pol0_Flx(vs),pol1_Flx(vs),pol1_Flx(vs),Flx_neg(q,p));
    2270             :     } else
    2271      262947 :       gel(V,++n) = Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
    2272             :   }
    2273      151994 :   y = Flx_extgcd_basecase(x,y,p,pi,&u,&v);
    2274      262972 :   for (i = n; i>1; i--)
    2275             :   {
    2276      110979 :     GEN R = gel(V,i);
    2277      110979 :     GEN u1 = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
    2278      110979 :     GEN v1 = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
    2279      110979 :     u = u1; v = v1;
    2280             :   }
    2281             :   {
    2282      151993 :     GEN R = gel(V,1);
    2283      151993 :     if (ptu)
    2284        6543 :       *ptu = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
    2285      151993 :     *ptv   = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
    2286             :   }
    2287      151994 :   return y;
    2288             : }
    2289             : 
    2290             : /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
    2291             :  * ux + vy = gcd (mod p) */
    2292             : GEN
    2293     6619220 : Flx_extgcd_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
    2294             : {
    2295     6619220 :   pari_sp av = avma;
    2296             :   GEN d;
    2297     6619220 :   long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
    2298     6619245 :   if (lgpol(y) >= lim)
    2299      151994 :     d = Flx_extgcd_halfgcd(x, y, p, pi, ptu, ptv);
    2300             :   else
    2301     6467243 :     d = Flx_extgcd_basecase(x, y, p, pi, ptu, ptv);
    2302     6619085 :   return gc_all(av, ptu?3:2, &d, ptv, ptu);
    2303             : }
    2304             : GEN
    2305      841317 : Flx_extgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
    2306      841317 : { return Flx_extgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptu, ptv); }
    2307             : 
    2308             : static GEN
    2309        1044 : Flx_halfres_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, ulong *r)
    2310             : {
    2311             :   struct Flx_res res;
    2312             :   GEN R;
    2313             :   long dB;
    2314             : 
    2315        1044 :   res.res  = *r;
    2316        1044 :   res.lc   = Flx_lead(y);
    2317        1044 :   res.deg0 = degpol(x);
    2318        1044 :   res.deg1 = degpol(y);
    2319        1044 :   res.off = 0;
    2320        1044 :   R = Flx_halfres_i(x, y, p, pi, a, b, &res);
    2321        1044 :   dB = degpol(*b);
    2322        1044 :   if (dB < degpol(y))
    2323        1044 :     Flx_halfres_update_pre(res.deg0, res.deg1, dB, p, pi, &res);
    2324        1044 :   *r = res.res;
    2325        1044 :   return R;
    2326             : }
    2327             : 
    2328             : static ulong
    2329     7076787 : Flx_resultant_basecase_pre(GEN a, GEN b, ulong p, ulong pi)
    2330             : {
    2331             :   pari_sp av;
    2332             :   long da,db,dc;
    2333     7076787 :   ulong lb, res = 1UL;
    2334             :   GEN c;
    2335             : 
    2336     7076787 :   da = degpol(a);
    2337     7076662 :   db = degpol(b);
    2338     7076853 :   if (db > da)
    2339             :   {
    2340           0 :     swapspec(a,b, da,db);
    2341           0 :     if (both_odd(da,db)) res = p-res;
    2342             :   }
    2343     7076853 :   else if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
    2344     7076853 :   av = avma;
    2345    54681976 :   while (db)
    2346             :   {
    2347    47622453 :     lb = b[db+2];
    2348    47622453 :     c = Flx_rem_pre(a,b, p,pi);
    2349    47424031 :     a = b; b = c; dc = degpol(c);
    2350    47417410 :     if (dc < 0) return gc_long(av,0);
    2351             : 
    2352    47404845 :     if (both_odd(da,db)) res = p - res;
    2353    47409334 :     if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, da - dc, p, pi), p);
    2354    47605615 :     if (gc_needed(av,2))
    2355             :     {
    2356           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant (da = %ld)",da);
    2357           0 :       gerepileall(av,2, &a,&b);
    2358             :     }
    2359    47605123 :     da = db; /* = degpol(a) */
    2360    47605123 :     db = dc; /* = degpol(b) */
    2361             :   }
    2362     7059523 :   return gc_ulong(av, Fl_mul(res, Fl_powu_pre(b[2], da, p, pi), p));
    2363             : }
    2364             : 
    2365             : ulong
    2366     7078662 : Flx_resultant_pre(GEN x, GEN y, ulong p, ulong pi)
    2367             : {
    2368     7078662 :   pari_sp av = avma;
    2369             :   long lim;
    2370     7078662 :   ulong res = 1;
    2371     7078662 :   long dx = degpol(x), dy = degpol(y);
    2372     7078334 :   if (dx < 0 || dy < 0) return 0;
    2373     7076892 :   if (dx < dy)
    2374             :   {
    2375     1064518 :     swap(x,y);
    2376     1064518 :     if (both_odd(dx, dy))
    2377        1906 :       res = Fl_neg(res, p);
    2378             :   }
    2379     7076895 :   lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
    2380     7077736 :   while (lgpol(y) >= lim)
    2381             :   {
    2382         852 :     if (lgpol(y)<=(lgpol(x)>>1))
    2383             :     {
    2384           0 :       GEN r = Flx_rem_pre(x, y, p, pi);
    2385           0 :       long dx = degpol(x), dy = degpol(y), dr = degpol(r);
    2386           0 :       ulong ly = y[dy+2];
    2387           0 :       if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
    2388           0 :       if (both_odd(dx, dy))
    2389           0 :         res = Fl_neg(res, p);
    2390           0 :       x = y; y = r;
    2391             :     }
    2392         852 :     (void) Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
    2393         852 :     if (gc_needed(av,2))
    2394             :     {
    2395           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_res (y = %ld)",degpol(y));
    2396           0 :       gerepileall(av,2,&x,&y);
    2397             :     }
    2398             :   }
    2399     7076784 :   return gc_ulong(av, Fl_mul(res, Flx_resultant_basecase_pre(x, y, p, pi), p));
    2400             : }
    2401             : 
    2402             : ulong
    2403     4810437 : Flx_resultant(GEN a, GEN b, ulong p)
    2404     4810437 : { return Flx_resultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2405             : 
    2406             : /* If resultant is 0, *ptU and *ptV are not set */
    2407             : static ulong
    2408          53 : Flx_extresultant_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptU, GEN *ptV)
    2409             : {
    2410          53 :   GEN z,q,u,v, x = a, y = b;
    2411          53 :   ulong lb, res = 1UL;
    2412          53 :   pari_sp av = avma;
    2413             :   long dx, dy, dz;
    2414          53 :   long vs = a[1];
    2415             : 
    2416          53 :   u = pol0_Flx(vs);
    2417          53 :   v = pol1_Flx(vs); /* v = 1 */
    2418          53 :   dx = degpol(x);
    2419          53 :   dy = degpol(y);
    2420         764 :   while (dy)
    2421             :   { /* b u = x (a), b v = y (a) */
    2422         711 :     lb = y[dy+2];
    2423         711 :     q = Flx_divrem_pre(x,y, p, pi, &z);
    2424         711 :     x = y; y = z; /* (x,y) = (y, x - q y) */
    2425         711 :     dz = degpol(z); if (dz < 0) return gc_ulong(av,0);
    2426         711 :     z = Flx_sub(u, Flx_mul_pre(q,v, p, pi), p);
    2427         711 :     u = v; v = z; /* (u,v) = (v, u - q v) */
    2428             : 
    2429         711 :     if (both_odd(dx,dy)) res = p - res;
    2430         711 :     if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, dx-dz, p, pi), p);
    2431         711 :     dx = dy; /* = degpol(x) */
    2432         711 :     dy = dz; /* = degpol(y) */
    2433             :   }
    2434          53 :   res = Fl_mul(res, Fl_powu_pre(y[2], dx, p, pi), p);
    2435          53 :   lb = Fl_mul(res, Fl_inv(y[2],p), p);
    2436          53 :   v = gerepileuptoleaf(av, Flx_Fl_mul_pre(v, lb, p, pi));
    2437          53 :   av = avma;
    2438          53 :   u = Flx_sub(Fl_to_Flx(res,vs), Flx_mul_pre(b,v,p,pi), p);
    2439          53 :   u = gerepileuptoleaf(av, Flx_div_pre(u,a,p,pi)); /* = (res - b v) / a */
    2440          53 :   *ptU = u;
    2441          53 :   *ptV = v; return res;
    2442             : }
    2443             : 
    2444             : ulong
    2445          53 : Flx_extresultant_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptU, GEN *ptV)
    2446             : {
    2447          53 :   pari_sp av=avma;
    2448             :   GEN u, v, R;
    2449          53 :   long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
    2450          53 :   ulong res = 1, res1;
    2451          53 :   long dx = degpol(x), dy = degpol(y);
    2452          53 :   if (dy > dx)
    2453             :   {
    2454          13 :     swap(x,y); lswap(dx,dy);
    2455          13 :     if (both_odd(dx,dy)) res = p-res;
    2456          13 :     R = matJ2_FlxM(x[1]);
    2457          40 :   } else R = matid2_FlxM(x[1]);
    2458          53 :   if (dy < 0) return 0;
    2459         245 :   while (lgpol(y) >= lim)
    2460             :   {
    2461             :     GEN M;
    2462         192 :     if (lgpol(y)<=(lgpol(x)>>1))
    2463             :     {
    2464          20 :       GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
    2465          20 :       long dx = degpol(x), dy = degpol(y), dr = degpol(r);
    2466          20 :       ulong ly = y[dy+2];
    2467          20 :       if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
    2468          20 :       if (both_odd(dx, dy))
    2469           0 :         res = Fl_neg(res, p);
    2470          20 :       x = y; y = r;
    2471          20 :       R = Flx_FlxM_qmul(q, R, p,pi);
    2472             :     }
    2473         192 :     M = Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
    2474         192 :     if (!res) return gc_ulong(av, 0);
    2475         192 :     R = FlxM_mul2(M, R, p, pi);
    2476         192 :     gerepileall(av,3,&x,&y,&R);
    2477             :   }
    2478          53 :   res1 = Flx_extresultant_basecase(x,y,p,pi,&u,&v);
    2479          53 :   if (!res1) return gc_ulong(av, 0);
    2480          53 :   *ptU = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi), res, p, pi);
    2481          53 :   *ptV = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi), res, p, pi);
    2482          53 :   gerepileall(av, 2, ptU, ptV);
    2483          53 :   return Fl_mul(res1,res,p);
    2484             : }
    2485             : 
    2486             : ulong
    2487          53 : Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
    2488          53 : { return Flx_extresultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptU, ptV); }
    2489             : 
    2490             : /* allow pi = 0 (SMALL_ULONG) */
    2491             : ulong
    2492    43677028 : Flx_eval_powers_pre(GEN x, GEN y, ulong p, ulong pi)
    2493             : {
    2494    43677028 :   ulong l0, l1, h0, h1, v1,  i = 1, lx = lg(x)-1;
    2495             : 
    2496    43677028 :   if (lx == 1) return 0;
    2497    40924795 :   x++;
    2498    40924795 :   if (pi)
    2499             :   {
    2500             :     LOCAL_OVERFLOW;
    2501             :     LOCAL_HIREMAINDER;
    2502    40860301 :     l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
    2503    97668025 :     while (++i < lx)
    2504             :     {
    2505    56807724 :       l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
    2506    56807724 :       l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
    2507             :     }
    2508       81117 :     return v1? remlll_pre(v1, h1, l1, p, pi)
    2509    40941418 :              : remll_pre(h1, l1, p, pi);
    2510             :   }
    2511             :   else
    2512             :   {
    2513       64494 :     l1 = x[i] * y[i];
    2514    30912057 :     while (++i < lx) { l1 += x[i] * y[i]; if (l1 & HIGHBIT) l1 %= p; }
    2515       64494 :     return l1 % p;
    2516             :   }
    2517             : }
    2518             : 
    2519             : /* allow pi = 0 (SMALL_ULONG) */
    2520             : ulong
    2521    44789845 : Flx_eval_pre(GEN x, ulong y, ulong p, ulong pi)
    2522             : {
    2523    44789845 :   long i, n = degpol(x);
    2524             :   ulong t;
    2525    44788301 :   if (n <= 0) return n? 0: x[2];
    2526    29607428 :   if (n > 15)
    2527             :   {
    2528      180389 :     pari_sp av = avma;
    2529      180389 :     GEN v = Fl_powers_pre(y, n, p, pi);
    2530      180380 :     return gc_ulong(av, Flx_eval_powers_pre(x, v, p, pi));
    2531             :   }
    2532    29427039 :   i = n+2; t = x[i];
    2533    29427039 :   if (pi)
    2534             :   {
    2535   116244231 :     for (i--; i>=2; i--) t = Fl_addmul_pre(uel(x, i), t, y, p, pi);
    2536    28315687 :     return t;
    2537             :   }
    2538     2667584 :   for (i--; i>=2; i--) t = (t * y + x[i]) % p;
    2539     1114866 :   return t %= p;
    2540             : }
    2541             : ulong
    2542    20366552 : Flx_eval(GEN x, ulong y, ulong p)
    2543    20366552 : { return Flx_eval_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2544             : 
    2545             : ulong
    2546        3255 : Flv_prod_pre(GEN x, ulong p, ulong pi)
    2547             : {
    2548        3255 :   pari_sp ltop = avma;
    2549             :   GEN v;
    2550        3255 :   long i,k,lx = lg(x);
    2551        3255 :   if (lx == 1) return 1UL;
    2552        3255 :   if (lx == 2) return uel(x,1);
    2553        3003 :   v = cgetg(1+(lx << 1), t_VECSMALL);
    2554        3003 :   k = 1;
    2555       28593 :   for (i=1; i<lx-1; i+=2)
    2556       25590 :     uel(v,k++) = Fl_mul_pre(uel(x,i), uel(x,i+1), p, pi);
    2557        3003 :   if (i < lx) uel(v,k++) = uel(x,i);
    2558       13529 :   while (k > 2)
    2559             :   {
    2560       10526 :     lx = k; k = 1;
    2561       36116 :     for (i=1; i<lx-1; i+=2)
    2562       25590 :       uel(v,k++) = Fl_mul_pre(uel(v,i), uel(v,i+1), p, pi);
    2563       10526 :     if (i < lx) uel(v,k++) = uel(v,i);
    2564             :   }
    2565        3003 :   return gc_ulong(ltop, uel(v,1));
    2566             : }
    2567             : 
    2568             : ulong
    2569           0 : Flv_prod(GEN v, ulong p)
    2570             : {
    2571           0 :   return Flv_prod_pre(v, p, get_Fl_red(p));
    2572             : }
    2573             : 
    2574             : GEN
    2575           0 : FlxV_prod(GEN V, ulong p)
    2576             : {
    2577             :   struct _Flxq D;
    2578           0 :   D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2579           0 :   return gen_product(V, (void *)&D, &_Flx_mul);
    2580             : }
    2581             : 
    2582             : /* compute prod (x - a[i]) */
    2583             : GEN
    2584      725830 : Flv_roots_to_pol(GEN a, ulong p, long vs)
    2585             : {
    2586             :   struct _Flxq D;
    2587      725830 :   long i,k,lx = lg(a);
    2588             :   GEN p1;
    2589      725830 :   if (lx == 1) return pol1_Flx(vs);
    2590      725830 :   p1 = cgetg(lx, t_VEC);
    2591    11824298 :   for (k=1,i=1; i<lx-1; i+=2)
    2592    11097874 :     gel(p1,k++) = mkvecsmall4(vs, Fl_mul(a[i], a[i+1], p),
    2593    11098928 :                               Fl_neg(Fl_add(a[i],a[i+1],p),p), 1);
    2594      725370 :   if (i < lx)
    2595       58112 :     gel(p1,k++) = mkvecsmall3(vs, Fl_neg(a[i],p), 1);
    2596      725364 :   D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2597      725356 :   setlg(p1, k); return gen_product(p1, (void *)&D, _Flx_mul);
    2598             : }
    2599             : 
    2600             : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for "large" p */
    2601             : INLINE void
    2602    18670197 : Flv_inv_pre_indir(GEN w, GEN v, ulong p, ulong pi)
    2603             : {
    2604    18670197 :   pari_sp av = avma;
    2605    18670197 :   long n = lg(w), i;
    2606             :   ulong u;
    2607             :   GEN c;
    2608             : 
    2609    18670197 :   if (n == 1) return;
    2610    18670197 :   c = cgetg(n, t_VECSMALL); c[1] = w[1];
    2611    77057565 :   for (i = 2; i < n; ++i) c[i] = Fl_mul_pre(w[i], c[i-1], p, pi);
    2612    18817901 :   i = n-1; u = Fl_inv(c[i], p);
    2613    77454265 :   for ( ; i > 1; --i)
    2614             :   {
    2615    58605037 :     ulong t = Fl_mul_pre(u, c[i-1], p, pi);
    2616    58589517 :     u = Fl_mul_pre(u, w[i], p, pi); v[i] = t;
    2617             :   }
    2618    18849228 :   v[1] = u; set_avma(av);
    2619             : }
    2620             : 
    2621             : void
    2622    18244877 : Flv_inv_pre_inplace(GEN v, ulong p, ulong pi) { Flv_inv_pre_indir(v,v, p, pi); }
    2623             : 
    2624             : GEN
    2625       10849 : Flv_inv_pre(GEN w, ulong p, ulong pi)
    2626       10849 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_pre_indir(w, v, p, pi); return v; }
    2627             : 
    2628             : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for SMALL_ULONG p */
    2629             : INLINE void
    2630       48739 : Flv_inv_indir(GEN w, GEN v, ulong p)
    2631             : {
    2632       48739 :   pari_sp av = avma;
    2633       48739 :   long n = lg(w), i;
    2634             :   ulong u;
    2635             :   GEN c;
    2636             : 
    2637       48739 :   if (n == 1) return;
    2638       48739 :   c = cgetg(n, t_VECSMALL); c[1] = w[1];
    2639     1709586 :   for (i = 2; i < n; ++i) c[i] = Fl_mul(w[i], c[i-1], p);
    2640       48743 :   i = n-1; u = Fl_inv(c[i], p);
    2641     1709668 :   for ( ; i > 1; --i)
    2642             :   {
    2643     1660921 :     ulong t = Fl_mul(u, c[i-1], p);
    2644     1660919 :     u = Fl_mul(u, w[i], p); v[i] = t;
    2645             :   }
    2646       48747 :   v[1] = u; set_avma(av);
    2647             : }
    2648             : static void
    2649      443269 : Flv_inv_i(GEN v, GEN w, ulong p)
    2650             : {
    2651      443269 :   if (SMALL_ULONG(p)) Flv_inv_indir(w, v, p);
    2652      394530 :   else Flv_inv_pre_indir(w, v, p, get_Fl_red(p));
    2653      443281 : }
    2654             : void
    2655       12017 : Flv_inv_inplace(GEN v, ulong p) { Flv_inv_i(v, v, p); }
    2656             : GEN
    2657      431259 : Flv_inv(GEN w, ulong p)
    2658      431259 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_i(v, w, p); return v; }
    2659             : 
    2660             : GEN
    2661    32895388 : Flx_div_by_X_x(GEN a, ulong x, ulong p, ulong *rem)
    2662             : {
    2663    32895388 :   long l = lg(a), i;
    2664             :   GEN a0, z0, z;
    2665    32895388 :   if (l <= 3)
    2666             :   {
    2667           0 :     if (rem) *rem = l == 2? 0: a[2];
    2668           0 :     return zero_Flx(a[1]);
    2669             :   }
    2670    32895388 :   z = cgetg(l-1,t_VECSMALL); z[1] = a[1];
    2671    32732466 :   a0 = a + l-1;
    2672    32732466 :   z0 = z + l-2; *z0 = *a0--;
    2673    32732466 :   if (SMALL_ULONG(p))
    2674             :   {
    2675    79370668 :     for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
    2676             :     {
    2677    58836050 :       ulong t = (*a0-- + x *  *z0--) % p;
    2678    58836050 :       *z0 = (long)t;
    2679             :     }
    2680    20534618 :     if (rem) *rem = (*a0 + x *  *z0) % p;
    2681             :   }
    2682             :   else
    2683             :   {
    2684    48192657 :     for (i=l-3; i>1; i--)
    2685             :     {
    2686    35938146 :       ulong t = Fl_add((ulong)*a0--, Fl_mul(x, *z0--, p), p);
    2687    35994809 :       *z0 = (long)t;
    2688             :     }
    2689    12254511 :     if (rem) *rem = Fl_add((ulong)*a0, Fl_mul(x, *z0, p), p);
    2690             :   }
    2691    32835062 :   return z;
    2692             : }
    2693             : 
    2694             : /* xa, ya = t_VECSMALL */
    2695             : static GEN
    2696      432456 : Flv_producttree(GEN xa, GEN s, ulong p, ulong pi, long vs)
    2697             : {
    2698      432456 :   long n = lg(xa)-1;
    2699      432456 :   long m = n==1 ? 1: expu(n-1)+1;
    2700      432456 :   long i, j, k, ls = lg(s);
    2701      432456 :   GEN T = cgetg(m+1, t_VEC);
    2702      432447 :   GEN t = cgetg(ls, t_VEC);
    2703     4572490 :   for (j=1, k=1; j<ls; k+=s[j++])
    2704     4139864 :     gel(t, j) = s[j] == 1 ?
    2705     4140038 :              mkvecsmall3(vs, Fl_neg(xa[k], p), 1):
    2706     1304185 :              mkvecsmall4(vs, Fl_mul(xa[k], xa[k+1], p),
    2707     1304180 :                  Fl_neg(Fl_add(xa[k],xa[k+1],p),p), 1);
    2708      432452 :   gel(T,1) = t;
    2709     1395464 :   for (i=2; i<=m; i++)
    2710             :   {
    2711      963066 :     GEN u = gel(T, i-1);
    2712      963066 :     long n = lg(u)-1;
    2713      963066 :     GEN t = cgetg(((n+1)>>1)+1, t_VEC);
    2714     4670047 :     for (j=1, k=1; k<n; j++, k+=2)
    2715     3707035 :       gel(t, j) = Flx_mul_pre(gel(u, k), gel(u, k+1), p, pi);
    2716      963012 :     gel(T, i) = t;
    2717             :   }
    2718      432398 :   return T;
    2719             : }
    2720             : 
    2721             : static GEN
    2722      472762 : Flx_Flv_multieval_tree(GEN P, GEN xa, GEN T, ulong p, ulong pi)
    2723             : {
    2724             :   long i,j,k;
    2725      472762 :   long m = lg(T)-1;
    2726      472762 :   GEN R = cgetg(lg(xa), t_VECSMALL);
    2727      472756 :   GEN Tp = cgetg(m+1, t_VEC), t;
    2728      472751 :   gel(Tp, m) = mkvec(P);
    2729     1621160 :   for (i=m-1; i>=1; i--)
    2730             :   {
    2731     1148403 :     GEN u = gel(T, i), v = gel(Tp, i+1);
    2732     1148403 :     long n = lg(u)-1;
    2733     1148403 :     t = cgetg(n+1, t_VEC);
    2734     5885549 :     for (j=1, k=1; k<n; j++, k+=2)
    2735             :     {
    2736     4737139 :       gel(t, k)   = Flx_rem_pre(gel(v, j), gel(u, k), p, pi);
    2737     4737048 :       gel(t, k+1) = Flx_rem_pre(gel(v, j), gel(u, k+1), p, pi);
    2738             :     }
    2739     1148410 :     gel(Tp, i) = t;
    2740             :   }
    2741             :   {
    2742      472757 :     GEN u = gel(T, i+1), v = gel(Tp, i+1);
    2743      472757 :     long n = lg(u)-1;
    2744     5685026 :     for (j=1, k=1; j<=n; j++)
    2745             :     {
    2746     5212212 :       long c, d = degpol(gel(u,j));
    2747    11977570 :       for (c=1; c<=d; c++, k++) R[k] = Flx_eval_pre(gel(v, j), xa[k], p, pi);
    2748             :     }
    2749      472814 :     return gc_const((pari_sp)R, R);
    2750             :   }
    2751             : }
    2752             : 
    2753             : static GEN
    2754     1184534 : FlvV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, ulong p, ulong pi, long vs)
    2755             : {
    2756     1184534 :   pari_sp av = avma;
    2757     1184534 :   long m = lg(T)-1;
    2758     1184534 :   long i, j, k, ls = lg(s);
    2759     1184534 :   GEN Tp = cgetg(m+1, t_VEC);
    2760     1183995 :   GEN t = cgetg(ls, t_VEC);
    2761    21563137 :   for (j=1, k=1; j<ls; k+=s[j++])
    2762    20379342 :     if (s[j]==2)
    2763             :     {
    2764     6709863 :       ulong a = Fl_mul(ya[k], R[k], p);
    2765     6711159 :       ulong b = Fl_mul(ya[k+1], R[k+1], p);
    2766     6720687 :       gel(t, j) = mkvecsmall3(vs, Fl_neg(Fl_add(Fl_mul(xa[k], b, p ),
    2767     6716757 :                   Fl_mul(xa[k+1], a, p), p), p), Fl_add(a, b, p));
    2768     6713009 :       gel(t, j) = Flx_renormalize(gel(t, j), 4);
    2769             :     }
    2770             :     else
    2771    13669479 :       gel(t, j) = Fl_to_Flx(Fl_mul(ya[k], R[k], p), vs);
    2772     1183795 :   gel(Tp, 1) = t;
    2773     5386891 :   for (i=2; i<=m; i++)
    2774             :   {
    2775     4202889 :     GEN u = gel(T, i-1);
    2776     4202889 :     GEN t = cgetg(lg(gel(T,i)), t_VEC);
    2777     4199744 :     GEN v = gel(Tp, i-1);
    2778     4199744 :     long n = lg(v)-1;
    2779    23351165 :     for (j=1, k=1; k<n; j++, k+=2)
    2780    19162008 :       gel(t, j) = Flx_add(Flx_mul_pre(gel(u, k), gel(v, k+1), p, pi),
    2781    19148069 :                           Flx_mul_pre(gel(u, k+1), gel(v, k), p, pi), p);
    2782     4203096 :     gel(Tp, i) = t;
    2783             :   }
    2784     1184002 :   return gerepileuptoleaf(av, gmael(Tp,m,1));
    2785             : }
    2786             : 
    2787             : GEN
    2788           0 : Flx_Flv_multieval(GEN P, GEN xa, ulong p)
    2789             : {
    2790           0 :   pari_sp av = avma;
    2791           0 :   GEN s = producttree_scheme(lg(xa)-1);
    2792           0 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2793           0 :   GEN T = Flv_producttree(xa, s, p, pi, P[1]);
    2794           0 :   return gerepileuptoleaf(av, Flx_Flv_multieval_tree(P, xa, T, p, pi));
    2795             : }
    2796             : 
    2797             : static GEN
    2798        2471 : FlxV_Flv_multieval_tree(GEN x, GEN xa, GEN T, ulong p, ulong pi)
    2799       45248 : { pari_APPLY_same(Flx_Flv_multieval_tree(gel(x,i), xa, T, p, pi)) }
    2800             : 
    2801             : GEN
    2802        2471 : FlxV_Flv_multieval(GEN P, GEN xa, ulong p)
    2803             : {
    2804        2471 :   pari_sp av = avma;
    2805        2471 :   GEN s = producttree_scheme(lg(xa)-1);
    2806        2471 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2807        2471 :   GEN T = Flv_producttree(xa, s, p, pi, P[1]);
    2808        2471 :   return gerepileupto(av, FlxV_Flv_multieval_tree(P, xa, T, p, pi));
    2809             : }
    2810             : 
    2811             : GEN
    2812      177136 : Flv_polint(GEN xa, GEN ya, ulong p, long vs)
    2813             : {
    2814      177136 :   pari_sp av = avma;
    2815      177136 :   GEN s = producttree_scheme(lg(xa)-1);
    2816      177138 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2817      177136 :   GEN T = Flv_producttree(xa, s, p, pi, vs);
    2818      177140 :   long m = lg(T)-1;
    2819      177140 :   GEN P = Flx_deriv(gmael(T, m, 1), p);
    2820      177142 :   GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
    2821      177142 :   return gerepileuptoleaf(av, FlvV_polint_tree(T, R, s, xa, ya, p, pi, vs));
    2822             : }
    2823             : 
    2824             : GEN
    2825       99861 : Flv_Flm_polint(GEN xa, GEN ya, ulong p, long vs)
    2826             : {
    2827       99861 :   pari_sp av = avma;
    2828       99861 :   GEN s = producttree_scheme(lg(xa)-1);
    2829       99862 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2830       99862 :   GEN T = Flv_producttree(xa, s, p, pi, vs);
    2831       99858 :   long i, m = lg(T)-1, l = lg(ya)-1;
    2832       99858 :   GEN P = Flx_deriv(gmael(T, m, 1), p);
    2833       99857 :   GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
    2834       99863 :   GEN M = cgetg(l+1, t_VEC);
    2835     1107092 :   for (i=1; i<=l; i++)
    2836     1007230 :     gel(M,i) = FlvV_polint_tree(T, R, s, xa, gel(ya,i), p, pi, vs);
    2837       99862 :   return gerepileupto(av, M);
    2838             : }
    2839             : 
    2840             : GEN
    2841      152987 : Flv_invVandermonde(GEN L, ulong den, ulong p)
    2842             : {
    2843      152987 :   pari_sp av = avma;
    2844      152987 :   long i, n = lg(L);
    2845             :   GEN M, R;
    2846      152987 :   GEN s = producttree_scheme(n-1);
    2847      152987 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2848      152987 :   GEN tree = Flv_producttree(L, s, p, pi, 0);
    2849      152987 :   long m = lg(tree)-1;
    2850      152987 :   GEN T = gmael(tree, m, 1);
    2851      152987 :   R = Flv_inv(Flx_Flv_multieval_tree(Flx_deriv(T, p), L, tree, p, pi), p);
    2852      152987 :   if (den!=1) R = Flv_Fl_mul(R, den, p);
    2853      152987 :   M = cgetg(n, t_MAT);
    2854      600481 :   for (i = 1; i < n; i++)
    2855             :   {
    2856      447494 :     GEN P = Flx_Fl_mul(Flx_div_by_X_x(T, uel(L,i), p, NULL), uel(R,i), p);
    2857      447494 :     gel(M,i) = Flx_to_Flv(P, n-1);
    2858             :   }
    2859      152987 :   return gerepilecopy(av, M);
    2860             : }
    2861             : 
    2862             : /***********************************************************************/
    2863             : /**                               Flxq                                **/
    2864             : /***********************************************************************/
    2865             : /* Flxq objects are Flx modulo another Flx called q. */
    2866             : 
    2867             : /* Product of y and x in Z/pZ[X]/(T), as t_VECSMALL. */
    2868             : GEN
    2869   186694403 : Flxq_mul_pre(GEN x,GEN y,GEN T,ulong p,ulong pi)
    2870   186694403 : { return Flx_rem_pre(Flx_mul_pre(x,y,p,pi),T,p,pi); }
    2871             : GEN
    2872    13149316 : Flxq_mul(GEN x,GEN y,GEN T,ulong p)
    2873    13149316 : { return Flxq_mul_pre(x,y,T,p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2874             : 
    2875             : GEN
    2876   275387934 : Flxq_sqr_pre(GEN x,GEN T,ulong p,ulong pi)
    2877   275387934 : { return Flx_rem_pre(Flx_sqr_pre(x, p,pi), T, p,pi); }
    2878             : /* Square of y in Z/pZ[X]/(T), as t_VECSMALL. */
    2879             : GEN
    2880     2725270 : Flxq_sqr(GEN x,GEN T,ulong p)
    2881     2725270 : { return Flxq_sqr_pre(x,T,p,SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2882             : 
    2883             : static GEN
    2884     1767764 : _Flxq_red(void *E, GEN x)
    2885     1767764 : { struct _Flxq *s = (struct _Flxq *)E;
    2886     1767764 :   return Flx_rem_pre(x, s->T, s->p, s->pi); }
    2887             : #if 0
    2888             : static GEN
    2889             : _Flx_sub(void *E, GEN x, GEN y)
    2890             : { struct _Flxq *s = (struct _Flxq *)E;
    2891             :   return Flx_sub(x,y,s->p); }
    2892             : #endif
    2893             : static GEN
    2894   268117428 : _Flxq_sqr(void *data, GEN x)
    2895             : {
    2896   268117428 :   struct _Flxq *D = (struct _Flxq*)data;
    2897   268117428 :   return Flxq_sqr_pre(x, D->T, D->p, D->pi);
    2898             : }
    2899             : static GEN
    2900   145325698 : _Flxq_mul(void *data, GEN x, GEN y)
    2901             : {
    2902   145325698 :   struct _Flxq *D = (struct _Flxq*)data;
    2903   145325698 :   return Flxq_mul_pre(x,y, D->T, D->p, D->pi);
    2904             : }
    2905             : static GEN
    2906    22499912 : _Flxq_one(void *data)
    2907             : {
    2908    22499912 :   struct _Flxq *D = (struct _Flxq*)data;
    2909    22499912 :   return pol1_Flx(get_Flx_var(D->T));
    2910             : }
    2911             : 
    2912             : static GEN
    2913    22158814 : _Flxq_powu_i(struct _Flxq *D, GEN x, ulong n)
    2914    22158814 : { return gen_powu_i(x, n, (void*)D, &_Flxq_sqr, &_Flxq_mul); }
    2915             : static GEN
    2916          68 : _Flxq_powu(struct _Flxq *D, GEN x, ulong n)
    2917          68 : { pari_sp av = avma; return gerepileuptoleaf(av, _Flxq_powu_i(D, x, n)); }
    2918             : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
    2919             : GEN
    2920    22887540 : Flxq_powu_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
    2921             : {
    2922             :   pari_sp av;
    2923             :   struct _Flxq D;
    2924    22887540 :   switch(n)
    2925             :   {
    2926           0 :     case 0: return pol1_Flx(get_Flx_var(T));
    2927      141663 :     case 1: return Flx_copy(x);
    2928      593653 :     case 2: return Flxq_sqr_pre(x, T, p, pi);
    2929             :   }
    2930    22152224 :   av = avma; set_Flxq_pre(&D, T, p, pi);
    2931    22158785 :   return gerepileuptoleaf(av, _Flxq_powu_i(&D, x, n));
    2932             : }
    2933             : GEN
    2934      486582 : Flxq_powu(GEN x, ulong n, GEN T, ulong p)
    2935      486582 : { return Flxq_powu_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2936             : 
    2937             : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
    2938             : GEN
    2939    24763748 : Flxq_pow_pre(GEN x, GEN n, GEN T, ulong p, ulong pi)
    2940             : {
    2941    24763748 :   pari_sp av = avma;
    2942             :   struct _Flxq D;
    2943             :   GEN y;
    2944    24763748 :   long s = signe(n);
    2945    24763748 :   if (!s) return pol1_Flx(get_Flx_var(T));
    2946    24528042 :   if (s < 0) x = Flxq_inv_pre(x,T,p,pi);
    2947    24528037 :   if (is_pm1(n)) return s < 0 ? x : Flx_copy(x);
    2948    23352306 :   set_Flxq_pre(&D, T, p, pi);
    2949    23352356 :   y = gen_pow_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
    2950    23352182 :   return gerepileuptoleaf(av, y);
    2951             : }
    2952             : GEN
    2953      871141 : Flxq_pow(GEN x, GEN n, GEN T, ulong p)
    2954      871141 : { return Flxq_pow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2955             : 
    2956             : GEN
    2957          28 : Flxq_pow_init_pre(GEN x, GEN n, long k, GEN T, ulong p, ulong pi)
    2958             : {
    2959          28 :   struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
    2960          28 :   return gen_pow_init(x, n, k, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
    2961             : }
    2962             : GEN
    2963           0 : Flxq_pow_init(GEN x, GEN n, long k, GEN T, ulong p)
    2964           0 : { return Flxq_pow_init_pre(x, n, k, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2965             : 
    2966             : GEN
    2967        4397 : Flxq_pow_table_pre(GEN R, GEN n, GEN T, ulong p, ulong pi)
    2968             : {
    2969        4397 :   struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
    2970        4397 :   return gen_pow_table(R, n, (void*)&D, &_Flxq_one, &_Flxq_mul);
    2971             : }
    2972             : GEN
    2973           0 : Flxq_pow_table(GEN R, GEN n, GEN T, ulong p)
    2974           0 : { return Flxq_pow_table_pre(R, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2975             : 
    2976             : /* Inverse of x in Z/lZ[X]/(T) or NULL if inverse doesn't exist
    2977             :  * not stack clean. */
    2978             : GEN
    2979     5777950 : Flxq_invsafe_pre(GEN x, GEN T, ulong p, ulong pi)
    2980             : {
    2981     5777950 :   GEN V, z = Flx_extgcd_pre(get_Flx_mod(T), x, p, pi, NULL, &V);
    2982             :   ulong iz;
    2983     5778060 :   if (degpol(z)) return NULL;
    2984     5777370 :   iz = Fl_inv(uel(z,2), p);
    2985     5777376 :   return Flx_Fl_mul_pre(V, iz, p, pi);
    2986             : }
    2987             : GEN
    2988      540462 : Flxq_invsafe(GEN x, GEN T, ulong p)
    2989      540462 : { return Flxq_invsafe_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2990             : 
    2991             : GEN
    2992     4793770 : Flxq_inv_pre(GEN x, GEN T, ulong p, ulong pi)
    2993             : {
    2994     4793770 :   pari_sp av=avma;
    2995     4793770 :   GEN U = Flxq_invsafe_pre(x, T, p, pi);
    2996     4793722 :   if (!U) pari_err_INV("Flxq_inv",Flx_to_ZX(x));
    2997     4793694 :   return gerepileuptoleaf(av, U);
    2998             : }
    2999             : GEN
    3000      345212 : Flxq_inv(GEN x, GEN T, ulong p)
    3001      345212 : { return Flxq_inv_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3002             : 
    3003             : GEN
    3004     2399860 : Flxq_div_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
    3005             : {
    3006     2399860 :   pari_sp av = avma;
    3007     2399860 :   return gerepileuptoleaf(av, Flxq_mul_pre(x,Flxq_inv_pre(y,T,p,pi),T,p,pi));
    3008             : }
    3009             : GEN
    3010      225649 : Flxq_div(GEN x, GEN y, GEN T, ulong p)
    3011      225649 : { return Flxq_div_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3012             : 
    3013             : GEN
    3014    22500717 : Flxq_powers_pre(GEN x, long l, GEN T, ulong p, ulong pi)
    3015             : {
    3016    22500717 :   int use_sqr = 2*degpol(x) >= get_Flx_degree(T);
    3017    22490541 :   struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
    3018    22488987 :   return gen_powers(x, l, use_sqr, (void*)&D, &_Flxq_sqr, &_Flxq_mul, &_Flxq_one);
    3019             : }
    3020             : GEN
    3021      232169 : Flxq_powers(GEN x, long l, GEN T, ulong p)
    3022      232169 : { return Flxq_powers_pre(x, l, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3023             : 
    3024             : GEN
    3025      167575 : Flxq_matrix_pow_pre(GEN y, long n, long m, GEN P, ulong l, ulong li)
    3026      167575 : { return FlxV_to_Flm(Flxq_powers_pre(y,m-1,P,l,li),n); }
    3027             : GEN
    3028         399 : Flxq_matrix_pow(GEN y, long n, long m, GEN P, ulong l)
    3029         399 : { return Flxq_matrix_pow_pre(y, n, m, P, l, SMALL_ULONG(l)? 0: get_Fl_red(l)); }
    3030             : 
    3031             : GEN
    3032    12625946 : Flx_Frobenius_pre(GEN T, ulong p, ulong pi)
    3033    12625946 : { return Flxq_powu_pre(polx_Flx(get_Flx_var(T)), p, T, p, pi); }
    3034             : GEN
    3035       86551 : Flx_Frobenius(GEN T, ulong p)
    3036       86551 : { return Flx_Frobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3037             : 
    3038             : GEN
    3039       85029 : Flx_matFrobenius_pre(GEN T, ulong p, ulong pi)
    3040             : {
    3041       85029 :   long n = get_Flx_degree(T);
    3042       85028 :   return Flxq_matrix_pow_pre(Flx_Frobenius_pre(T, p, pi), n, n, T, p, pi);
    3043             : }
    3044             : GEN
    3045           0 : Flx_matFrobenius(GEN T, ulong p)
    3046           0 : { return Flx_matFrobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3047             : 
    3048             : static GEN
    3049    12956501 : Flx_blocks_Flm(GEN P, long n, long m)
    3050             : {
    3051    12956501 :   GEN z = cgetg(m+1,t_MAT);
    3052    12953811 :   long i,j, k=2, l = lg(P);
    3053    38377191 :   for(i=1; i<=m; i++)
    3054             :   {
    3055    25427458 :     GEN zi = cgetg(n+1,t_VECSMALL);
    3056    25423380 :     gel(z,i) = zi;
    3057   113009140 :     for(j=1; j<=n; j++)
    3058    87585760 :       uel(zi, j) = k==l ? 0 : uel(P,k++);
    3059             :   }
    3060    12949733 :   return z;
    3061             : }
    3062             : 
    3063             : GEN
    3064      516562 : Flx_blocks(GEN P, long n, long m)
    3065             : {
    3066      516562 :   GEN z = cgetg(m+1,t_VEC);
    3067      516256 :   long i,j, k=2, l = lg(P);
    3068     1546858 :   for(i=1; i<=m; i++)
    3069             :   {
    3070     1030857 :     GEN zi = cgetg(n+2,t_VECSMALL);
    3071     1029895 :     zi[1] = P[1];
    3072     1029895 :     gel(z,i) = zi;
    3073     6462599 :     for(j=2; j<n+2; j++)
    3074     5432704 :       uel(zi, j) = k==l ? 0 : uel(P,k++);
    3075     1029895 :     zi = Flx_renormalize(zi, n+2);
    3076             :   }
    3077      516001 :   return z;
    3078             : }
    3079             : 
    3080             : static GEN
    3081    12956963 : FlxV_to_Flm_lg(GEN x, long m, long n)
    3082             : {
    3083             :   long i;
    3084    12956963 :   GEN y = cgetg(n+1, t_MAT);
    3085    58954057 :   for (i=1; i<=n; i++) gel(y,i) = Flx_to_Flv(gel(x,i), m);
    3086    12953229 :   return y;
    3087             : }
    3088             : 
    3089             : /* allow pi = 0 (SMALL_ULONG) */
    3090             : GEN
    3091    13154347 : Flx_FlxqV_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
    3092             : {
    3093    13154347 :   pari_sp btop, av = avma;
    3094    13154347 :   long sv = get_Flx_var(T), m = get_Flx_degree(T);
    3095    13154186 :   long i, l = lg(x)-1, lQ = lgpol(Q), n,  d;
    3096             :   GEN A, B, C, S, g;
    3097    13155175 :   if (lQ == 0) return pol0_Flx(sv);
    3098    12958162 :   if (lQ <= l)
    3099             :   {
    3100     5779844 :     n = l;
    3101     5779844 :     d = 1;
    3102             :   }
    3103             :   else
    3104             :   {
    3105     7178318 :     n = l-1;
    3106     7178318 :     d = (lQ+n-1)/n;
    3107             :   }
    3108    12958162 :   A = FlxV_to_Flm_lg(x, m, n);
    3109    12955836 :   B = Flx_blocks_Flm(Q, n, d);
    3110    12953260 :   C = gerepileupto(av, Flm_mul(A, B, p));
    3111    12959001 :   g = gel(x, l);
    3112    12959001 :   if (pi && SMALL_ULONG(p)) pi = 0;
    3113    12959001 :   T = Flx_get_red_pre(T, p, pi);
    3114    12957825 :   btop = avma;
    3115    12957825 :   S = Flv_to_Flx(gel(C, d), sv);
    3116    25436117 :   for (i = d-1; i>0; i--)
    3117             :   {
    3118    12480772 :     S = Flx_add(Flxq_mul_pre(S, g, T, p, pi), Flv_to_Flx(gel(C,i), sv), p);
    3119    12479962 :     if (gc_needed(btop,1))
    3120           0 :       S = gerepileuptoleaf(btop, S);
    3121             :   }
    3122    12955345 :   return gerepileuptoleaf(av, S);
    3123             : }
    3124             : GEN
    3125        5089 : Flx_FlxqV_eval(GEN Q, GEN x, GEN T, ulong p)
    3126        5089 : { return Flx_FlxqV_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3127             : 
    3128             : /* allow pi = 0 (SMALL_ULONG) */
    3129             : GEN
    3130     2622291 : Flx_Flxq_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
    3131             : {
    3132     2622291 :   pari_sp av = avma;
    3133             :   GEN z, V;
    3134     2622291 :   long d = degpol(Q), rtd;
    3135     2622293 :   if (d < 0) return pol0_Flx(get_Flx_var(T));
    3136     2622202 :   rtd = (long) sqrt((double)d);
    3137     2622202 :   T = Flx_get_red_pre(T, p, pi);
    3138     2622214 :   V = Flxq_powers_pre(x, rtd, T, p, pi);
    3139     2622260 :   z = Flx_FlxqV_eval_pre(Q, V, T, p, pi);
    3140     2622211 :   return gerepileupto(av, z);
    3141             : }
    3142             : GEN
    3143      786176 : Flx_Flxq_eval(GEN Q, GEN x, GEN T, ulong p)
    3144      786176 : { return Flx_Flxq_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3145             : 
    3146             : /* allow pi = 0 (SMALL_ULONG) */
    3147             : GEN
    3148           0 : FlxC_FlxqV_eval_pre(GEN x, GEN v, GEN T, ulong p, ulong pi)
    3149           0 : { pari_APPLY_type(t_COL, Flx_FlxqV_eval_pre(gel(x,i), v, T, p, pi)) }
    3150             : GEN
    3151           0 : FlxC_FlxqV_eval(GEN x, GEN v, GEN T, ulong p)
    3152           0 : { return FlxC_FlxqV_eval_pre(x, v, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3153             : 
    3154             : /* allow pi = 0 (SMALL_ULONG) */
    3155             : GEN
    3156           0 : FlxC_Flxq_eval_pre(GEN x, GEN F, GEN T, ulong p, ulong pi)
    3157             : {
    3158           0 :   long d = brent_kung_optpow(get_Flx_degree(T)-1,lg(x)-1,1);
    3159           0 :   GEN Fp = Flxq_powers_pre(F, d, T, p, pi);
    3160           0 :   return FlxC_FlxqV_eval_pre(x, Fp, T, p, pi);
    3161             : }
    3162             : GEN
    3163           0 : FlxC_Flxq_eval(GEN x, GEN F, GEN T, ulong p)
    3164           0 : { return FlxC_Flxq_eval_pre(x, F, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3165             : 
    3166             : #if 0
    3167             : static struct bb_algebra Flxq_algebra = { _Flxq_red, _Flx_add, _Flx_sub,
    3168             :               _Flxq_mul, _Flxq_sqr, _Flxq_one, _Flxq_zero};
    3169             : #endif
    3170             : 
    3171             : static GEN
    3172      386265 : Flxq_autpow_sqr(void *E, GEN x)
    3173             : {
    3174      386265 :   struct _Flxq *D = (struct _Flxq*)E;
    3175      386265 :   return Flx_Flxq_eval_pre(x, x, D->T, D->p, D->pi);
    3176             : }
    3177             : static GEN
    3178       20713 : Flxq_autpow_msqr(void *E, GEN x)
    3179             : {
    3180       20713 :   struct _Flxq *D = (struct _Flxq*)E;
    3181       20713 :   return Flx_FlxqV_eval_pre(Flxq_autpow_sqr(E, x), D->aut, D->T, D->p, D->pi);
    3182             : }
    3183             : 
    3184             : GEN
    3185      304445 : Flxq_autpow_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
    3186             : {
    3187      304445 :   pari_sp av = avma;
    3188             :   struct _Flxq D;
    3189             :   long d;
    3190      304445 :   if (n==0) return Flx_rem_pre(polx_Flx(x[1]), T, p, pi);
    3191      304438 :   if (n==1) return Flx_rem_pre(x, T, p, pi);
    3192      303955 :   set_Flxq_pre(&D, T, p, pi);
    3193      303955 :   d = brent_kung_optpow(get_Flx_degree(T), hammingl(n)-1, 1);
    3194      303955 :   D.aut = Flxq_powers_pre(x, d, T, p, D.pi);
    3195      303955 :   x = gen_powu_fold_i(x,n,(void*)&D,Flxq_autpow_sqr,Flxq_autpow_msqr);
    3196      303955 :   return gerepilecopy(av, x);
    3197             : }
    3198             : GEN
    3199           7 : Flxq_autpow(GEN x, ulong n, GEN T, ulong p)
    3200           7 : { return Flxq_autpow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3201             : 
    3202             : GEN
    3203        1733 : Flxq_autpowers(GEN x, ulong l, GEN T, ulong p)
    3204             : {
    3205        1733 :   long d, vT = get_Flx_var(T), dT = get_Flx_degree(T);
    3206             :   ulong i, pi;
    3207        1733 :   pari_sp av = avma;
    3208        1733 :   GEN xp, V = cgetg(l+2,t_VEC);
    3209        1733 :   gel(V,1) = polx_Flx(vT); if (l==0) return V;
    3210        1733 :   gel(V,2) = gcopy(x); if (l==1) return V;
    3211        1733 :   pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    3212        1733 :   T = Flx_get_red_pre(T, p, pi);
    3213        1733 :   d = brent_kung_optpow(dT-1, l-1, 1);
    3214        1733 :   xp = Flxq_powers_pre(x, d, T, p, pi);
    3215        7376 :   for(i = 3; i < l+2; i++)
    3216        5643 :     gel(V,i) = Flx_FlxqV_eval_pre(gel(V,i-1), xp, T, p, pi);
    3217        1733 :   return gerepilecopy(av, V);
    3218             : }
    3219             : 
    3220             : static GEN
    3221      588669 : Flxq_autsum_mul(void *E, GEN x, GEN y)
    3222             : {
    3223      588669 :   struct _Flxq *D = (struct _Flxq*)E;
    3224      588669 :   GEN T = D->T;
    3225      588669 :   ulong p = D->p, pi = D->pi;
    3226      588669 :   GEN phi1 = gel(x,1), a1 = gel(x,2);
    3227      588669 :   GEN phi2 = gel(y,1), a2 = gel(y,2);
    3228      588669 :   ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
    3229      588669 :   GEN V2 = Flxq_powers_pre(phi2, d, T, p, pi);
    3230      588669 :   GEN phi3 = Flx_FlxqV_eval_pre(phi1, V2, T, p, pi);
    3231      588669 :   GEN aphi = Flx_FlxqV_eval_pre(a1, V2, T, p, pi);
    3232      588669 :   GEN a3 = Flxq_mul_pre(aphi, a2, T, p, pi);
    3233      588669 :   return mkvec2(phi3, a3);
    3234             : }
    3235             : static GEN
    3236      379814 : Flxq_autsum_sqr(void *E, GEN x)
    3237      379814 : { return Flxq_autsum_mul(E, x, x); }
    3238             : 
    3239             : static GEN
    3240      315636 : Flxq_autsum_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
    3241             : {
    3242      315636 :   pari_sp av = avma;
    3243      315636 :   struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
    3244      315636 :   x = gen_powu_i(x,n,(void*)&D,Flxq_autsum_sqr,Flxq_autsum_mul);
    3245      315636 :   return gerepilecopy(av, x);
    3246             : }
    3247             : GEN
    3248           0 : Flxq_autsum(GEN x, ulong n, GEN T, ulong p)
    3249           0 : { return Flxq_autsum_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3250             : 
    3251             : static GEN
    3252      759919 : Flxq_auttrace_mul(void *E, GEN x, GEN y)
    3253             : {
    3254      759919 :   struct _Flxq *D = (struct _Flxq*)E;
    3255      759919 :   GEN T = D->T;
    3256      759919 :   ulong p = D->p, pi = D->pi;
    3257      759919 :   GEN phi1 = gel(x,1), a1 = gel(x,2);
    3258      759919 :   GEN phi2 = gel(y,1), a2 = gel(y,2);
    3259      759919 :   ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
    3260      759935 :   GEN V1 = Flxq_powers_pre(phi1, d, T, p, pi);
    3261      759867 :   GEN phi3 = Flx_FlxqV_eval_pre(phi2, V1, T, p, pi);
    3262      759908 :   GEN aphi = Flx_FlxqV_eval_pre(a2, V1, T, p, pi);
    3263      759916 :   GEN a3 = Flx_add(a1, aphi, p);
    3264      759917 :   return mkvec2(phi3, a3);
    3265             : }
    3266             : 
    3267             : static GEN
    3268      633111 : Flxq_auttrace_sqr(void *E, GEN x)
    3269      633111 : { return Flxq_auttrace_mul(E, x, x); }
    3270             : 
    3271             : GEN
    3272      929090 : Flxq_auttrace_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
    3273             : {
    3274      929090 :   pari_sp av = avma;
    3275             :   struct _Flxq D;
    3276      929090 :   set_Flxq_pre(&D, T, p, pi);
    3277      929088 :   x = gen_powu_i(x,n,(void*)&D,Flxq_auttrace_sqr,Flxq_auttrace_mul);
    3278      929070 :   return gerepilecopy(av, x);
    3279             : }
    3280             : GEN
    3281           0 : Flxq_auttrace(GEN x, ulong n, GEN T, ulong p)
    3282           0 : { return Flxq_auttrace_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3283             : 
    3284             : static long
    3285      849765 : bounded_order(ulong p, GEN b, long k)
    3286             : {
    3287      849765 :   GEN a = modii(utoipos(p), b);
    3288             :   long i;
    3289     2078418 :   for(i = 1; i < k; i++)
    3290             :   {
    3291     1612846 :     if (equali1(a)) return i;
    3292     1228656 :     a = modii(muliu(a,p),b);
    3293             :   }
    3294      465572 :   return 0;
    3295             : }
    3296             : 
    3297             : /*
    3298             :   n = (p^d-a)\b
    3299             :   b = bb*p^vb
    3300             :   p^k = 1 [bb]
    3301             :   d = m*k+r+vb
    3302             :   u = (p^k-1)/bb;
    3303             :   v = (p^(r+vb)-a)/b;
    3304             :   w = (p^(m*k)-1)/(p^k-1)
    3305             :   n = p^r*w*u+v
    3306             :   w*u = p^vb*(p^(m*k)-1)/b
    3307             :   n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b
    3308             : */
    3309             : 
    3310             : static GEN
    3311    24058081 : Flxq_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, ulong p, ulong pi)
    3312             : {
    3313    24058081 :   pari_sp av=avma;
    3314    24058081 :   long d = get_Flx_degree(T);
    3315    24058105 :   GEN an = absi_shallow(n), z, q;
    3316    24058075 :   if (abscmpiu(an,p)<0 || cmpis(an,d)<=0) return Flxq_pow_pre(x, n, T, p, pi);
    3317      850170 :   q = powuu(p, d);
    3318      850169 :   if (dvdii(q, n))
    3319             :   {
    3320         370 :     long vn = logint(an, utoipos(p));
    3321         370 :     GEN autvn = vn==1 ? aut: Flxq_autpow_pre(aut,vn,T,p,pi);
    3322         370 :     z = Flx_Flxq_eval_pre(x,autvn,T,p,pi);
    3323             :   } else
    3324             :   {
    3325      849796 :     GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
    3326             :     GEN bb, u, v, autk;
    3327      849794 :     long vb = Z_lvalrem(b,p,&bb);
    3328      849799 :     long m, r, k = is_pm1(bb)? 1: bounded_order(p,bb,d);
    3329      849795 :     if (!k || d-vb < k) return Flxq_pow_pre(x,n, T,p,pi);
    3330      384216 :     m = (d-vb)/k; r = (d-vb)%k;
    3331      384216 :     u = diviiexact(subiu(powuu(p,k),1),bb);
    3332      384216 :     v = diviiexact(subii(powuu(p,r+vb),a),b);
    3333      384216 :     autk = k==1 ? aut: Flxq_autpow_pre(aut,k,T,p,pi);
    3334      384216 :     if (r)
    3335             :     {
    3336       94534 :       GEN autr = r==1 ? aut: Flxq_autpow_pre(aut,r,T,p,pi);
    3337       94534 :       z = Flx_Flxq_eval_pre(x,autr,T,p,pi);
    3338      289682 :     } else z = x;
    3339      384216 :     if (m > 1) z = gel(Flxq_autsum_pre(mkvec2(autk, z), m, T, p, pi), 2);
    3340      384216 :     if (!is_pm1(u)) z = Flxq_pow_pre(z, u, T, p, pi);
    3341      384216 :     if (signe(v)) z = Flxq_mul_pre(z, Flxq_pow_pre(x, v, T, p, pi), T, p, pi);
    3342             :   }
    3343      384586 :   return gerepileupto(av,signe(n)>0 ? z : Flxq_inv_pre(z,T,p,pi));
    3344             : }
    3345             : 
    3346             : static GEN
    3347    24050604 : _Flxq_pow(void *data, GEN x, GEN n)
    3348             : {
    3349    24050604 :   struct _Flxq *D = (struct _Flxq*)data;
    3350    24050604 :   return Flxq_pow_Frobenius(x, n, D->aut, D->T, D->p, D->pi);
    3351             : }
    3352             : 
    3353             : static GEN
    3354      328965 : _Flxq_rand(void *data)
    3355             : {
    3356      328965 :   pari_sp av=avma;
    3357      328965 :   struct _Flxq *D = (struct _Flxq*)data;
    3358             :   GEN z;
    3359             :   do
    3360             :   {
    3361      330283 :     set_avma(av);
    3362      330284 :     z = random_Flx(get_Flx_degree(D->T),get_Flx_var(D->T),D->p);
    3363      330286 :   } while (lgpol(z)==0);
    3364      328968 :   return z;
    3365             : }
    3366             : 
    3367             : /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
    3368             : static GEN
    3369       35763 : Fl_Flxq_log(ulong a, GEN g, GEN o, GEN T, ulong p)
    3370             : {
    3371       35763 :   pari_sp av = avma;
    3372             :   GEN q,n_q,ord,ordp, op;
    3373             : 
    3374       35763 :   if (a == 1UL) return gen_0;
    3375             :   /* p > 2 */
    3376             : 
    3377       35763 :   ordp = utoi(p - 1);
    3378       35763 :   ord  = get_arith_Z(o);
    3379       35763 :   if (!ord) ord = T? subiu(powuu(p, get_FpX_degree(T)), 1): ordp;
    3380       35763 :   if (a == p - 1) /* -1 */
    3381        7718 :     return gerepileuptoint(av, shifti(ord,-1));
    3382       28045 :   ordp = gcdii(ordp, ord);
    3383       28045 :   op = typ(o)==t_MAT ? famat_Z_gcd(o, ordp) : ordp;
    3384             : 
    3385       28045 :   q = NULL;
    3386       28045 :   if (T)
    3387             :   { /* we want < g > = Fp^* */
    3388       28045 :     if (!equalii(ord,ordp)) {
    3389       11899 :       q = diviiexact(ord,ordp);
    3390       11899 :       g = Flxq_pow(g,q,T,p);
    3391             :     }
    3392             :   }
    3393       28045 :   n_q = Fp_log(utoi(a), utoipos(uel(g,2)), op, utoipos(p));
    3394       28045 :   if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
    3395       28045 :   if (q) n_q = mulii(q, n_q);
    3396       28045 :   return gerepileuptoint(av, n_q);
    3397             : }
    3398             : 
    3399             : static GEN
    3400      691194 : Flxq_easylog(void* E, GEN a, GEN g, GEN ord)
    3401             : {
    3402      691194 :   struct _Flxq *f = (struct _Flxq *)E;
    3403      691194 :   GEN T = f->T;
    3404      691194 :   ulong p = f->p;
    3405      691194 :   long d = get_Flx_degree(T);
    3406      691194 :   if (Flx_equal1(a)) return gen_0;
    3407      531092 :   if (Flx_equal(a,g)) return gen_1;
    3408      174146 :   if (!degpol(a))
    3409       35763 :     return Fl_Flxq_log(uel(a,2), g, ord, T, p);
    3410      138383 :   if (typ(ord)!=t_INT || d <= 4 || d == 6 || abscmpiu(ord,1UL<<27)<0)
    3411      138355 :     return NULL;
    3412          28 :   return Flxq_log_index(a, g, ord, T, p);
    3413             : }
    3414             : 
    3415             : static const struct bb_group Flxq_star={_Flxq_mul,_Flxq_pow,_Flxq_rand,hash_GEN,Flx_equal,Flx_equal1,Flxq_easylog};
    3416             : 
    3417             : const struct bb_group *
    3418      427699 : get_Flxq_star(void **E, GEN T, ulong p)
    3419             : {
    3420      427699 :   struct _Flxq *e = (struct _Flxq *) stack_malloc(sizeof(struct _Flxq));
    3421      427699 :   e->T = T; e->p  = p; e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    3422      427699 :   e->aut =  Flx_Frobenius_pre(T, p, e->pi);
    3423      427700 :   *E = (void*)e; return &Flxq_star;
    3424             : }
    3425             : 
    3426             : GEN
    3427       95205 : Flxq_order(GEN a, GEN ord, GEN T, ulong p)
    3428             : {
    3429             :   void *E;
    3430       95205 :   const struct bb_group *S = get_Flxq_star(&E,T,p);
    3431       95205 :   return gen_order(a,ord,E,S);
    3432             : }
    3433             : 
    3434             : GEN
    3435      164532 : Flxq_log(GEN a, GEN g, GEN ord, GEN T, ulong p)
    3436             : {
    3437             :   void *E;
    3438      164532 :   pari_sp av = avma;
    3439      164532 :   const struct bb_group *S = get_Flxq_star(&E,T,p);
    3440      164532 :   GEN v = get_arith_ZZM(ord), F = gmael(v,2,1);
    3441      164532 :   if (Flxq_log_use_index(veclast(F), T, p))
    3442       24311 :     v = mkvec2(gel(v, 1), ZM_famat_limit(gel(v, 2), int2n(27)));
    3443      164532 :   return gerepileuptoleaf(av, gen_PH_log(a, g, v, E, S));
    3444             : }
    3445             : 
    3446             : GEN
    3447      170916 : Flxq_sqrtn(GEN a, GEN n, GEN T, ulong p, GEN *zeta)
    3448             : {
    3449      170916 :   if (!lgpol(a))
    3450             :   {
    3451        2954 :     if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
    3452        2947 :     if (zeta)
    3453           0 :       *zeta=pol1_Flx(get_Flx_var(T));
    3454        2947 :     return pol0_Flx(get_Flx_var(T));
    3455             :   }
    3456             :   else
    3457             :   {
    3458             :     void *E;
    3459      167962 :     pari_sp av = avma;
    3460      167962 :     const struct bb_group *S = get_Flxq_star(&E,T,p);
    3461      167963 :     GEN o = subiu(powuu(p,get_Flx_degree(T)), 1);
    3462      167960 :     GEN s = gen_Shanks_sqrtn(a,n,o,zeta,E,S);
    3463      167962 :     if (!s) return gc_NULL(av);
    3464      167864 :     return gc_all(av, zeta?2:1, &s, zeta);
    3465             :   }
    3466             : }
    3467             : 
    3468             : GEN
    3469      218670 : Flxq_sqrt_pre(GEN z, GEN T, ulong p, ulong pi)
    3470             : {
    3471      218670 :   pari_sp av = avma;
    3472      218670 :   if (p==2)
    3473             :   {
    3474           0 :     GEN r = F2xq_sqrt(Flx_to_F2x(z), Flx_to_F2x(get_Flx_mod(T)));
    3475           0 :     return gerepileupto(av, F2x_to_Flx(r));
    3476             :   }
    3477      218670 :   if (get_Flx_degree(T)==2)
    3478             :   {
    3479       67856 :     GEN P = get_Flx_mod(T), s;
    3480       67856 :     ulong c = uel(P,2), b = uel(P,3), a = uel(P,4);
    3481       67856 :     ulong y = degpol(z)<1 ? 0: uel(z,3);
    3482       67856 :     if (a==1 && b==0)
    3483       14890 :     {
    3484       15670 :       ulong x = degpol(z)<1 ? Flx_constant(z): uel(z,2);
    3485       15670 :       GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), Fl_neg(c, p), p, pi);
    3486       15670 :       if (!r) return gc_NULL(av);
    3487       14890 :       s = mkvecsmall3(P[1], uel(r,1), uel(r,2));
    3488             :     }
    3489             :     else
    3490             :     {
    3491       52186 :       ulong b2 = Fl_halve(b, p), t = Fl_div(b2, a, p);
    3492       52186 :       ulong D = Fl_sub(Fl_sqr(b2, p), Fl_mul(a, c, p), p);
    3493       52186 :       ulong x = degpol(z)<1 ? Flx_constant(z): Fl_sub(uel(z,2), Fl_mul(uel(z,3), t, p), p);
    3494       52186 :       GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), D, p, pi);
    3495       52186 :       if (!r) return gc_NULL(av);
    3496       49792 :       s = mkvecsmall3(P[1], Fl_add(uel(r,1), Fl_mul(uel(r,2),t,p), p), uel(r,2));
    3497             :     }
    3498       64682 :     return gerepileuptoleaf(av, Flx_renormalize(s, 4));
    3499             :   }
    3500             :   else
    3501      150814 :     return Flxq_sqrtn(z, gen_2, T, p, NULL);
    3502             : }
    3503             : 
    3504             : GEN
    3505      218670 : Flxq_sqrt(GEN a, GEN T, ulong p)
    3506      218670 : { return Flxq_sqrt_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3507             : 
    3508             : /* assume T irreducible mod p */
    3509             : int
    3510      380306 : Flxq_issquare(GEN x, GEN T, ulong p)
    3511             : {
    3512      380306 :   if (lgpol(x) == 0 || p == 2) return 1;
    3513      377107 :   return krouu(Flxq_norm(x,T,p), p) == 1;
    3514             : }
    3515             : 
    3516             : /* assume T irreducible mod p */
    3517             : int
    3518           0 : Flxq_is2npower(GEN x, long n, GEN T, ulong p)
    3519             : {
    3520             :   pari_sp av;
    3521             :   GEN m;
    3522           0 :   if (n==1) return Flxq_issquare(x, T, p);
    3523           0 :   if (lgpol(x) == 0 || p == 2) return 1;
    3524           0 :   av = avma;
    3525           0 :   m = shifti(subiu(powuu(p, get_Flx_degree(T)), 1), -n);
    3526           0 :   return gc_bool(av, Flx_equal1(Flxq_pow(x, m, T, p)));
    3527             : }
    3528             : 
    3529             : GEN
    3530      113589 : Flxq_lroot_fast_pre(GEN a, GEN sqx, GEN T, long p, ulong pi)
    3531             : {
    3532      113589 :   pari_sp av=avma;
    3533      113589 :   GEN A = Flx_splitting(a,p);
    3534      113589 :   return gerepileuptoleaf(av, FlxqV_dotproduct_pre(A,sqx,T,p,pi));
    3535             : }
    3536             : GEN
    3537           0 : Flxq_lroot_fast(GEN a, GEN sqx, GEN T, long p)
    3538           0 : { return Flxq_lroot_fast_pre(a, sqx, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3539             : 
    3540             : GEN
    3541       25053 : Flxq_lroot_pre(GEN a, GEN T, long p, ulong pi)
    3542             : {
    3543       25053 :   pari_sp av=avma;
    3544       25053 :   long n = get_Flx_degree(T), d = degpol(a);
    3545             :   GEN sqx, V;
    3546       25053 :   if (n==1) return leafcopy(a);
    3547       25053 :   if (n==2) return Flxq_powu_pre(a, p, T, p, pi);
    3548       25053 :   sqx = Flxq_autpow_pre(Flx_Frobenius_pre(T, p, pi), n-1, T, p, pi);
    3549       25053 :   if (d==1 && a[2]==0 && a[3]==1) return gerepileuptoleaf(av, sqx);
    3550           0 :   if (d>=p)
    3551             :   {
    3552           0 :     V = Flxq_powers_pre(sqx,p-1,T,p,pi);
    3553           0 :     return gerepileuptoleaf(av, Flxq_lroot_fast_pre(a,V,T,p,pi));
    3554             :   } else
    3555           0 :     return gerepileuptoleaf(av, Flx_Flxq_eval_pre(a,sqx,T,p,pi));
    3556             : }
    3557             : GEN
    3558           0 : Flxq_lroot(GEN a, GEN T, long p)
    3559           0 : { return Flxq_lroot_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3560             : 
    3561             : ulong
    3562      422459 : Flxq_norm(GEN x, GEN TB, ulong p)
    3563             : {
    3564      422459 :   GEN T = get_Flx_mod(TB);
    3565      422459 :   ulong y = Flx_resultant(T, x, p), L = Flx_lead(T);
    3566      422459 :   if (L==1 || lgpol(x)==0) return y;
    3567           0 :   return Fl_div(y, Fl_powu(L, (ulong)degpol(x), p), p);
    3568             : }
    3569             : 
    3570             : ulong
    3571        3219 : Flxq_trace(GEN x, GEN TB, ulong p)
    3572             : {
    3573        3219 :   pari_sp av = avma;
    3574             :   ulong t;
    3575        3219 :   GEN T = get_Flx_mod(TB);
    3576        3219 :   long n = degpol(T)-1;
    3577        3219 :   GEN z = Flxq_mul(x, Flx_deriv(T, p), TB, p);
    3578        3219 :   t = degpol(z)<n ? 0 : Fl_div(z[2+n],T[3+n],p);
    3579        3219 :   return gc_ulong(av, t);
    3580             : }
    3581             : 
    3582             : /*x must be reduced*/
    3583             : GEN
    3584        3632 : Flxq_charpoly(GEN x, GEN TB, ulong p)
    3585             : {
    3586        3632 :   pari_sp ltop=avma;
    3587        3632 :   GEN T = get_Flx_mod(TB);
    3588        3632 :   long vs = evalvarn(fetch_var());
    3589        3632 :   GEN xm1 = deg1pol_shallow(pol1_Flx(x[1]),Flx_neg(x,p),vs);
    3590        3632 :   GEN r = Flx_FlxY_resultant(T, xm1, p);
    3591        3632 :   r[1] = x[1];
    3592        3632 :   (void)delete_var(); return gerepileupto(ltop, r);
    3593             : }
    3594             : 
    3595             : /* Computing minimal polynomial :                         */
    3596             : /* cf Shoup 'Efficient Computation of Minimal Polynomials */
    3597             : /*          in Algebraic Extensions of Finite Fields'     */
    3598             : 
    3599             : /* Let v a linear form, return the linear form z->v(tau*z)
    3600             :    that is, v*(M_tau) */
    3601             : 
    3602             : static GEN
    3603     1682837 : Flxq_transmul_init(GEN tau, GEN T, ulong p, ulong pi)
    3604             : {
    3605             :   GEN bht;
    3606     1682837 :   GEN h, Tp = get_Flx_red(T, &h);
    3607     1682833 :   long n = degpol(Tp), vT = Tp[1];
    3608     1682831 :   GEN ft = Flx_recipspec(Tp+2, n+1, n+1);
    3609     1682807 :   GEN bt = Flx_recipspec(tau+2, lgpol(tau), n);
    3610     1682810 :   ft[1] = vT; bt[1] = vT;
    3611     1682810 :   if (h)
    3612        2688 :     bht = Flxn_mul_pre(bt, h, n-1, p, pi);
    3613             :   else
    3614             :   {
    3615     1680122 :     GEN bh = Flx_div_pre(Flx_shift(tau, n-1), T, p, pi);
    3616     1680116 :     bht = Flx_recipspec(bh+2, lgpol(bh), n-1);
    3617     1680109 :     bht[1] = vT;
    3618             :   }
    3619     1682797 :   return mkvec3(bt, bht, ft);
    3620             : }
    3621             : 
    3622             : static GEN
    3623     4064455 : Flxq_transmul(GEN tau, GEN a, long n, ulong p, ulong pi)
    3624             : {
    3625     4064455 :   pari_sp ltop = avma;
    3626             :   GEN t1, t2, t3, vec;
    3627     4064455 :   GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
    3628     4064455 :   if (lgpol(a)==0) return pol0_Flx(a[1]);
    3629     4034256 :   t2  = Flx_shift(Flx_mul_pre(bt, a, p, pi),1-n);
    3630     4033987 :   if (lgpol(bht)==0) return gerepileuptoleaf(ltop, t2);
    3631     3047634 :   t1  = Flx_shift(Flx_mul_pre(ft, a, p, pi),-n);
    3632     3047640 :   t3  = Flxn_mul_pre(t1, bht, n-1, p, pi);
    3633     3047626 :   vec = Flx_sub(t2, Flx_shift(t3, 1), p);
    3634     3047651 :   return gerepileuptoleaf(ltop, vec);
    3635             : }
    3636             : 
    3637             : GEN
    3638      780335 : Flxq_minpoly_pre(GEN x, GEN T, ulong p, ulong pi)
    3639             : {
    3640      780335 :   pari_sp ltop = avma;
    3641      780335 :   long vT = get_Flx_var(T), n = get_Flx_degree(T);
    3642             :   GEN v_x;
    3643      780331 :   GEN g = pol1_Flx(vT), tau = pol1_Flx(vT);
    3644      780306 :   T = Flx_get_red_pre(T, p, pi);
    3645      780313 :   v_x = Flxq_powers_pre(x, usqrt(2*n), T, p, pi);
    3646     1621714 :   while (lgpol(tau) != 0)
    3647             :   {
    3648             :     long i, j, m, k1;
    3649             :     GEN M, v, tr, g_prime, c;
    3650      841409 :     if (degpol(g) == n) { tau = pol1_Flx(vT); g = pol1_Flx(vT); }
    3651      841410 :     v = random_Flx(n, vT, p);
    3652      841423 :     tr = Flxq_transmul_init(tau, T, p, pi);
    3653      841413 :     v = Flxq_transmul(tr, v, n, p, pi);
    3654      841429 :     m = 2*(n-degpol(g));
    3655      841427 :     k1 = usqrt(m);
    3656      841427 :     tr = Flxq_transmul_init(gel(v_x,k1+1), T, p, pi);
    3657      841407 :     c = cgetg(m+2,t_VECSMALL);
    3658      841357 :     c[1] = vT;
    3659     4064336 :     for (i=0; i<m; i+=k1)
    3660             :     {
    3661     3222913 :       long mj = minss(m-i, k1);
    3662    12599013 :       for (j=0; j<mj; j++)
    3663     9375638 :         uel(c,m+1-(i+j)) = Flx_dotproduct_pre(v, gel(v_x,j+1), p, pi);
    3664     3223375 :       v = Flxq_transmul(tr, v, n, p, pi);
    3665             :     }
    3666      841423 :     c = Flx_renormalize(c, m+2);
    3667             :     /* now c contains <v,x^i> , i = 0..m-1  */
    3668      841420 :     M = Flx_halfgcd_pre(monomial_Flx(1, m, vT), c, p, pi);
    3669      841430 :     g_prime = gmael(M, 2, 2);
    3670      841430 :     if (degpol(g_prime) < 1) continue;
    3671      829819 :     g = Flx_mul_pre(g, g_prime, p, pi);
    3672      829802 :     tau = Flxq_mul_pre(tau, Flx_FlxqV_eval_pre(g_prime, v_x, T,p,pi), T,p,pi);
    3673             :   }
    3674      780263 :   g = Flx_normalize(g,p);
    3675      780320 :   return gerepileuptoleaf(ltop,g);
    3676             : }
    3677             : GEN
    3678       44381 : Flxq_minpoly(GEN x, GEN T, ulong p)
    3679       44381 : { return Flxq_minpoly_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3680             : 
    3681             : GEN
    3682          20 : Flxq_conjvec(GEN x, GEN T, ulong p)
    3683             : {
    3684          20 :   long i, l = 1+get_Flx_degree(T);
    3685          20 :   GEN z = cgetg(l,t_COL);
    3686          20 :   struct _Flxq D; set_Flxq(&D, T, p);
    3687          20 :   gel(z,1) = Flx_copy(x);
    3688          88 :   for (i=2; i<l; i++) gel(z,i) = _Flxq_powu(&D, gel(z,i-1), p);
    3689          20 :   return z;
    3690             : }
    3691             : 
    3692             : GEN
    3693        7166 : gener_Flxq(GEN T, ulong p, GEN *po)
    3694             : {
    3695        7166 :   long i, j, vT = get_Flx_var(T), f = get_Flx_degree(T);
    3696             :   ulong p_1, pi;
    3697             :   GEN g, L, L2, o, q, F;
    3698             :   pari_sp av0, av;
    3699             : 
    3700        7166 :   if (f == 1) {
    3701             :     GEN fa;
    3702          28 :     o = utoipos(p-1);
    3703          28 :     fa = Z_factor(o);
    3704          28 :     L = gel(fa,1);
    3705          28 :     L = vecslice(L, 2, lg(L)-1); /* remove 2 for efficiency */
    3706          28 :     g = Fl_to_Flx(pgener_Fl_local(p, vec_to_vecsmall(L)), vT);
    3707          28 :     if (po) *po = mkvec2(o, fa);
    3708          28 :     return g;
    3709             :   }
    3710             : 
    3711        7138 :   av0 = avma; p_1 = p - 1;
    3712        7138 :   q = diviuexact(subiu(powuu(p,f), 1), p_1);
    3713             : 
    3714        7138 :   L = cgetg(1, t_VECSMALL);
    3715        7138 :   if (p > 3)
    3716             :   {
    3717        2343 :     ulong t = p_1 >> vals(p_1);
    3718        2343 :     GEN P = gel(factoru(t), 1);
    3719        2343 :     L = cgetg_copy(P, &i);
    3720        3745 :     while (--i) L[i] = p_1 / P[i];
    3721             :   }
    3722        7138 :   o = factor_pn_1(utoipos(p),f);
    3723        7138 :   L2 = leafcopy( gel(o, 1) );
    3724       19100 :   for (i = j = 1; i < lg(L2); i++)
    3725             :   {
    3726       11962 :     if (umodui(p_1, gel(L2,i)) == 0) continue;
    3727        6460 :     gel(L2,j++) = diviiexact(q, gel(L2,i));
    3728             :   }
    3729        7138 :   setlg(L2, j); pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    3730        7138 :   F = Flx_Frobenius_pre(T, p, pi);
    3731       17612 :   for (av = avma;; set_avma(av))
    3732       10474 :   {
    3733             :     GEN tt;
    3734       17612 :     g = random_Flx(f, vT, p);
    3735       17612 :     if (degpol(g) < 1) continue;
    3736       11945 :     if (p == 2) tt = g;
    3737             :     else
    3738             :     {
    3739        8746 :       ulong t = Flxq_norm(g, T, p);
    3740        8746 :       if (t == 1 || !is_gener_Fl(t, p, p_1, L)) continue;
    3741        4792 :       tt = Flxq_powu_pre(g, p_1>>1, T, p, pi);
    3742             :     }
    3743       14570 :     for (i = 1; i < j; i++)
    3744             :     {
    3745        7432 :       GEN a = Flxq_pow_Frobenius(tt, gel(L2,i), F, T, p, pi);
    3746        7432 :       if (!degpol(a) && uel(a,2) == p_1) break;
    3747             :     }
    3748        7991 :     if (i == j) break;
    3749             :   }
    3750        7138 :   if (!po)
    3751             :   {
    3752         187 :     set_avma((pari_sp)g);
    3753         187 :     g = gerepileuptoleaf(av0, g);
    3754             :   }
    3755             :   else {
    3756        6951 :     *po = mkvec2(subiu(powuu(p,f), 1), o);
    3757        6951 :     gerepileall(av0, 2, &g, po);
    3758             :   }
    3759        7138 :   return g;
    3760             : }
    3761             : 
    3762             : static GEN
    3763      553283 : _Flxq_neg(void *E, GEN x)
    3764      553283 : { struct _Flxq *s = (struct _Flxq *)E;
    3765      553283 :   return Flx_neg(x,s->p); }
    3766             : 
    3767             : static GEN
    3768     1572233 : _Flxq_rmul(void *E, GEN x, GEN y)
    3769     1572233 : { struct _Flxq *s = (struct _Flxq *)E;
    3770     1572233 :   return Flx_mul_pre(x,y,s->p,s->pi); }
    3771             : 
    3772             : static GEN
    3773       18987 : _Flxq_inv(void *E, GEN x)
    3774       18987 : { struct _Flxq *s = (struct _Flxq *)E;
    3775       18987 :   return Flxq_inv(x,s->T,s->p); }
    3776             : 
    3777             : static int
    3778      165019 : _Flxq_equal0(GEN x) { return lgpol(x)==0; }
    3779             : 
    3780             : static GEN
    3781       25164 : _Flxq_s(void *E, long x)
    3782       25164 : { struct _Flxq *s = (struct _Flxq *)E;
    3783       25164 :   ulong u = x<0 ? s->p+x: (ulong)x;
    3784       25164 :   return Fl_to_Flx(u, get_Flx_var(s->T));
    3785             : }
    3786             : 
    3787             : static const struct bb_field Flxq_field={_Flxq_red,_Flx_add,_Flxq_rmul,_Flxq_neg,
    3788             :                                          _Flxq_inv,_Flxq_equal0,_Flxq_s};
    3789             : 
    3790       73028 : const struct bb_field *get_Flxq_field(void **E, GEN T, ulong p)
    3791             : {
    3792       73028 :   GEN z = new_chunk(sizeof(struct _Flxq));
    3793       73028 :   set_Flxq((struct _Flxq *)z, T, p); *E = (void*)z; return &Flxq_field;
    3794             : }
    3795             : 
    3796             : /***********************************************************************/
    3797             : /**                               Flxn                                **/
    3798             : /***********************************************************************/
    3799             : 
    3800             : GEN
    3801       54331 : Flx_invLaplace(GEN x, ulong p)
    3802             : {
    3803       54331 :   long i, d = degpol(x);
    3804             :   ulong t;
    3805             :   GEN y;
    3806       54324 :   if (d <= 1) return Flx_copy(x);
    3807       54324 :   t = Fl_inv(factorial_Fl(d, p), p);
    3808       54387 :   y = cgetg(d+3, t_VECSMALL);
    3809       54300 :   y[1] = x[1];
    3810     1328321 :   for (i=d; i>=2; i--)
    3811             :   {
    3812     1273949 :     uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
    3813     1273989 :     t = Fl_mul(t, i, p);
    3814             :   }
    3815       54372 :   uel(y,3) = uel(x,3);
    3816       54372 :   uel(y,2) = uel(x,2);
    3817       54372 :   return y;
    3818             : }
    3819             : 
    3820             : GEN
    3821       27331 : Flx_Laplace(GEN x, ulong p)
    3822             : {
    3823       27331 :   long i, d = degpol(x);
    3824       27328 :   ulong t = 1;
    3825             :   GEN y;
    3826       27328 :   if (d <= 1) return Flx_copy(x);
    3827       27328 :   y = cgetg(d+3, t_VECSMALL);
    3828       27313 :   y[1] = x[1];
    3829       27313 :   uel(y,2) = uel(x,2);
    3830       27313 :   uel(y,3) = uel(x,3);
    3831      757546 :   for (i=2; i<=d; i++)
    3832             :   {
    3833      730201 :     t = Fl_mul(t, i%p, p);
    3834      730238 :     uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
    3835             :   }
    3836       27345 :   return y;
    3837             : }
    3838             : 
    3839             : GEN
    3840     6219385 : Flxn_red(GEN a, long n)
    3841             : {
    3842     6219385 :   long i, L, l = lg(a);
    3843             :   GEN  b;
    3844     6219385 :   if (l == 2 || !n) return zero_Flx(a[1]);
    3845     5829734 :   L = n+2; if (L > l) L = l;
    3846     5829734 :   b = cgetg(L, t_VECSMALL); b[1] = a[1];
    3847    58569510 :   for (i=2; i<L; i++) b[i] = a[i];
    3848     5826237 :   return Flx_renormalize(b,L);
    3849             : }
    3850             : 
    3851             : GEN
    3852     5050646 : Flxn_mul_pre(GEN a, GEN b, long n, ulong p, ulong pi)
    3853     5050646 : { return Flxn_red(Flx_mul_pre(a, b, p, pi), n); }
    3854             : GEN
    3855       75201 : Flxn_mul(GEN a, GEN b, long n, ulong p)
    3856       75201 : { return Flxn_mul_pre(a, b, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3857             : 
    3858             : GEN
    3859           0 : Flxn_sqr_pre(GEN a, long n, ulong p, ulong pi)
    3860           0 : { return Flxn_red(Flx_sqr_pre(a, p, pi), n); }
    3861             : GEN
    3862           0 : Flxn_sqr(GEN a, long n, ulong p)
    3863           0 : { return Flxn_sqr_pre(a, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3864             : 
    3865             : /* (f*g) \/ x^n */
    3866             : static GEN
    3867      938060 : Flx_mulhigh_i(GEN f, GEN g, long n, ulong p, ulong pi)
    3868      938060 : { return Flx_shift(Flx_mul_pre(f, g, p, pi),-n); }
    3869             : 
    3870             : static GEN
    3871      516509 : Flxn_mulhigh(GEN f, GEN g, long n2, long n, ulong p, ulong pi)
    3872             : {
    3873      516509 :   GEN F = Flx_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    3874      516136 :   return Flx_add(Flx_mulhigh_i(fl, g, n2, p, pi),
    3875             :                  Flxn_mul_pre(fh, g, n - n2, p, pi), p);
    3876             : }
    3877             : 
    3878             : /* g==NULL -> assume g==1 */
    3879             : GEN
    3880       55148 : Flxn_div_pre(GEN g, GEN f, long e, ulong p, ulong pi)
    3881             : {
    3882       55148 :   pari_sp av = avma, av2;
    3883             :   ulong mask;
    3884             :   GEN W;
    3885       55148 :   long n = 1;
    3886       55148 :   if (lg(f) <= 2) pari_err_INV("Flxn_inv",f);
    3887       55148 :   W = Fl_to_Flx(Fl_inv(uel(f,2),p), f[1]);
    3888       55145 :   mask = quadratic_prec_mask(e);
    3889       55169 :   av2 = avma;
    3890      258767 :   for (;mask>1;)
    3891             :   {
    3892             :     GEN u, fr;
    3893      203600 :     long n2 = n;
    3894      203600 :     n<<=1; if (mask & 1) n--;
    3895      203600 :     mask >>= 1;
    3896      203600 :     fr = Flxn_red(f, n);
    3897      203434 :     if (mask>1 || !g)
    3898             :     {
    3899      149405 :       u = Flxn_mul_pre(W, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
    3900      149634 :       W = Flx_sub(W, Flx_shift(u, n2), p);
    3901             :     } else
    3902             :     {
    3903       54029 :       GEN y = Flxn_mul_pre(g, W, n, p, pi), yt =  Flxn_red(y, n-n2);
    3904       54052 :       u = Flxn_mul_pre(yt, Flxn_mulhigh(fr,  W, n2, n, p, pi), n-n2, p, pi);
    3905       54061 :       W = Flx_sub(y, Flx_shift(u, n2), p);
    3906             :     }
    3907      203566 :     if (gc_needed(av2,2))
    3908             :     {
    3909           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Flxn_div, e = %ld", n);
    3910           0 :       W = gerepileupto(av2, W);
    3911             :     }
    3912             :   }
    3913       55167 :   return gerepileupto(av, W);
    3914             : }
    3915             : GEN
    3916       55124 : Flxn_div(GEN g, GEN f, long e, ulong p)
    3917       55124 : { return Flxn_div_pre(g, f, e, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    3918             : 
    3919             : GEN
    3920        1039 : Flxn_inv(GEN f, long e, ulong p)
    3921        1039 : { return Flxn_div(NULL, f, e, p); }
    3922             : 
    3923             : GEN
    3924      109372 : Flxn_expint(GEN h, long e, ulong p)
    3925             : {
    3926      109372 :   pari_sp av = avma, av2;
    3927      109372 :   long v = h[1], n=1;
    3928      109372 :   GEN f = pol1_Flx(v), g = pol1_Flx(v);
    3929      109351 :   ulong mask = quadratic_prec_mask(e), pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    3930      109360 :   av2 = avma;
    3931      422706 :   for (;mask>1;)
    3932             :   {
    3933             :     GEN u, w;
    3934      422667 :     long n2 = n;
    3935      422667 :     n<<=1; if (mask & 1) n--;
    3936      422667 :     mask >>= 1;
    3937      422667 :     u = Flxn_mul_pre(g, Flx_mulhigh_i(f, Flxn_red(h, n2-1), n2-1, p,pi), n-n2, p,pi);
    3938      422642 :     u = Flx_add(u, Flx_shift(Flxn_red(h, n-1), 1-n2), p);
    3939      422654 :     w = Flxn_mul_pre(f, Flx_integXn(u, n2-1, p), n-n2, p, pi);
    3940      422582 :     f = Flx_add(f, Flx_shift(w, n2), p);
    3941      422726 :     if (mask<=1) break;
    3942      313358 :     u = Flxn_mul_pre(g, Flxn_mulhigh(f, g, n2, n, p, pi), n-n2, p, pi);
    3943      313323 :     g = Flx_sub(g, Flx_shift(u, n2), p);
    3944      313346 :     if (gc_needed(av2,2))
    3945             :     {
    3946           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flxn_exp, e = %ld", n);
    3947           0 :       gerepileall(av2, 2, &f, &g);
    3948             :     }
    3949             :   }
    3950      109407 :   return gerepileupto(av, f);
    3951             : }
    3952             : 
    3953             : GEN
    3954           0 : Flxn_exp(GEN h, long e, ulong p)
    3955             : {
    3956           0 :   if (degpol(h)<1 || uel(h,2)!=0)
    3957           0 :     pari_err_DOMAIN("Flxn_exp","valuation", "<", gen_1, h);
    3958           0 :   return Flxn_expint(Flx_deriv(h, p), e, p);
    3959             : }
    3960             : 
    3961             : INLINE GEN
    3962      217176 : Flxn_recip(GEN x, long n)
    3963             : {
    3964      217176 :   GEN z=Flx_recipspec(x+2,lgpol(x),n);
    3965      217003 :   z[1]=x[1];
    3966      217003 :   return z;
    3967             : }
    3968             : 
    3969             : GEN
    3970       54056 : Flx_Newton(GEN P, long n, ulong p)
    3971             : {
    3972       54056 :   pari_sp av = avma;
    3973       54056 :   long d = degpol(P);
    3974       54053 :   GEN dP = Flxn_recip(Flx_deriv(P, p), d);
    3975       53971 :   GEN Q = Flxn_div(dP, Flxn_recip(P, d+1), n, p);
    3976       54027 :   return gerepileuptoleaf(av, Q);
    3977             : }
    3978             : 
    3979             : GEN
    3980      109380 : Flx_fromNewton(GEN P, ulong p)
    3981             : {
    3982      109380 :   pari_sp av = avma;
    3983      109380 :   ulong n = Flx_constant(P)+1;
    3984      109380 :   GEN z = Flx_neg(Flx_shift(P, -1), p);
    3985      109373 :   GEN Q = Flxn_recip(Flxn_expint(z, n, p), n);
    3986      109361 :   return gerepileuptoleaf(av, Q);
    3987             : }
    3988             : 
    3989             : static void
    3990       12514 : init_invlaplace(long d, ulong p, GEN *pt_P, GEN *pt_V)
    3991             : {
    3992             :   long i;
    3993             :   ulong e;
    3994       12514 :   GEN P = cgetg(d+1, t_VECSMALL);
    3995       12514 :   GEN V = cgetg(d+1, t_VECSMALL);
    3996     1396581 :   for (i=1, e=1; i<=d; i++, e++)
    3997             :   {
    3998     1384067 :     if (e==p)
    3999             :     {
    4000      459153 :       e = 0;
    4001      459153 :       V[i] = u_lvalrem(i, p, &uel(P,i));
    4002             :     } else
    4003             :     {
    4004      924914 :       V[i] = 0; uel(P,i) = i;
    4005             :     }
    4006             :   }
    4007       12514 :   *pt_P = P; *pt_V = V;
    4008       12514 : }
    4009             : 
    4010             : /* return p^val * FpX_invLaplace(1+x+...x^(n-1), q), with q a power of p and
    4011             :  * val large enough to compensate for the power of p in the factorials */
    4012             : 
    4013             : static GEN
    4014         497 : ZpX_invLaplace_init(long n, GEN q, ulong p, long v, long sv)
    4015             : {
    4016         497 :   pari_sp av = avma;
    4017         497 :   long i, d = n-1, w;
    4018             :   GEN y, W, E, t;
    4019         497 :   init_invlaplace(d, p, &E, &W);
    4020         497 :   t = Fp_inv(FpV_prod(Flv_to_ZV(E), q), q);
    4021         497 :   w = zv_sum(W);
    4022         497 :   if (v > w) t = Fp_mul(t, powuu(p, v-w), q);
    4023         497 :   y = cgetg(d+3,t_POL);
    4024         497 :   y[1] = evalsigne(1) | sv;
    4025       28882 :   for (i=d; i>=1; i--)
    4026             :   {
    4027       28385 :     gel(y,i+2) = t;
    4028       28385 :     t = Fp_mulu(t, uel(E,i), q);
    4029       28385 :     if (uel(W,i)) t = Fp_mul(t, powuu(p, uel(W,i)), q);
    4030             :   }
    4031         497 :   gel(y,2) = t;
    4032         497 :   return gerepilecopy(av, ZX_renormalize(y, d+3));
    4033             : }
    4034             : 
    4035             : GEN
    4036       27529 : Flx_composedsum(GEN P, GEN Q, ulong p)
    4037             : {
    4038       27529 :   pari_sp av = avma;
    4039       27529 :   long n = 1 + degpol(P)*degpol(Q);
    4040       27530 :   ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
    4041       27525 :                       Fl_powu(Flx_lead(Q), degpol(P), p), p);
    4042             :   GEN R;
    4043       27533 :   if (p >= (ulong)n)
    4044             :   {
    4045       27036 :     GEN Pl = Flx_invLaplace(Flx_Newton(P,n,p), p);
    4046       27028 :     GEN Ql = Flx_invLaplace(Flx_Newton(Q,n,p), p);
    4047       27037 :     GEN L  = Flx_Laplace(Flxn_mul(Pl, Ql, n, p), p);
    4048       27038 :     R = Flx_fromNewton(L, p);
    4049             :   } else
    4050             :   {
    4051         497 :     long v = factorial_lval(n-1, p);
    4052         497 :     long w = 1 + ulogint(n-1, p);
    4053         497 :     GEN pv = powuu(p, v);
    4054         497 :     GEN qf = powuu(p, w), q = mulii(pv, qf), q2 = mulii(q, pv);
    4055         497 :     GEN iL = ZpX_invLaplace_init(n, q, p, v, P[1]);
    4056         497 :     GEN Pl = FpX_convol(iL, FpX_Newton(Flx_to_ZX(P), n, qf), q);
    4057         497 :     GEN Ql = FpX_convol(iL, FpX_Newton(Flx_to_ZX(Q), n, qf), q);
    4058         497 :     GEN Ln = ZX_Z_divexact(FpXn_mul(Pl, Ql, n, q2), pv);
    4059         497 :     GEN L  = ZX_Z_divexact(FpX_Laplace(Ln, q), pv);
    4060         497 :     R = ZX_to_Flx(FpX_fromNewton(L, qf), p);
    4061             :   }
    4062       27521 :   return gerepileuptoleaf(av, Flx_Fl_mul(R, lead, p));
    4063             : }
    4064             : 
    4065             : static GEN
    4066        3826 : _Flx_composedsum(void *E, GEN a, GEN b)
    4067        3826 : { return Flx_composedsum(a, b, (ulong)E); }
    4068             : 
    4069             : GEN
    4070       28884 : FlxV_composedsum(GEN V, ulong p)
    4071       28884 : { return gen_product(V, (void *)p, &_Flx_composedsum); }
    4072             : 
    4073             : GEN
    4074           0 : Flx_composedprod(GEN P, GEN Q, ulong p)
    4075             : {
    4076           0 :   pari_sp av = avma;
    4077           0 :   long n = 1+ degpol(P)*degpol(Q);
    4078           0 :   ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
    4079           0 :                       Fl_powu(Flx_lead(Q), degpol(P), p), p);
    4080             :   GEN R;
    4081           0 :   if (p >= (ulong)n)
    4082             :   {
    4083           0 :     GEN L = Flx_convol(Flx_Newton(P,n,p), Flx_Newton(Q,n,p), p);
    4084           0 :     R = Flx_fromNewton(L, p);
    4085             :   } else
    4086             :   {
    4087           0 :     long w = 1 + ulogint(n, p);
    4088           0 :     GEN qf = powuu(p, w);
    4089           0 :     GEN Pl = FpX_convol(FpX_Newton(Flx_to_ZX(P), n, qf), FpX_Newton(Flx_to_ZX(Q), n, qf), qf);
    4090           0 :     R = ZX_to_Flx(FpX_fromNewton(Pl, qf), p);
    4091             :   }
    4092           0 :   return gerepileuptoleaf(av, Flx_Fl_mul(R, lead, p));
    4093             : 
    4094             : }
    4095             : 
    4096             : /* (x+1)^n mod p; assume 2 <= n < 2p prime */
    4097             : static GEN
    4098           0 : Fl_Xp1_powu(ulong n, ulong p, long v)
    4099             : {
    4100           0 :   ulong k, d = (n + 1) >> 1;
    4101           0 :   GEN C, V = identity_zv(d);
    4102             : 
    4103           0 :   Flv_inv_inplace(V, p); /* could restrict to odd integers in [3,d] */
    4104           0 :   C = cgetg(n+3, t_VECSMALL);
    4105           0 :   C[1] = v;
    4106           0 :   uel(C,2) = 1UL;
    4107           0 :   uel(C,3) = n%p;
    4108           0 :   uel(C,4) = Fl_mul(odd(n)? n: n-1, n >> 1, p);
    4109             :     /* binom(n,k) = binom(n,k-1) * (n-k+1) / k */
    4110           0 :   if (SMALL_ULONG(p))
    4111           0 :     for (k = 3; k <= d; k++)
    4112           0 :       uel(C,k+2) = Fl_mul(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p);
    4113             :   else
    4114             :   {
    4115           0 :     ulong pi  = get_Fl_red(p);
    4116           0 :     for (k = 3; k <= d; k++)
    4117           0 :       uel(C,k+2) = Fl_mul_pre(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p, pi);
    4118             :   }
    4119           0 :   for (   ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
    4120           0 :   return C; /* normalized */
    4121             : }
    4122             : 
    4123             : /* p arbitrary */
    4124             : GEN
    4125       28194 : Flx_translate1_basecase(GEN P, ulong p)
    4126             : {
    4127       28194 :   GEN R = Flx_copy(P);
    4128       28194 :   long i, k, n = degpol(P);
    4129      654732 :   for (i = 1; i <= n; i++)
    4130    14846488 :     for (k = n-i; k < n; k++) uel(R,k+2) = Fl_add(uel(R,k+2), uel(R,k+3), p);
    4131       28194 :   return R;
    4132             : }
    4133             : 
    4134             : static int
    4135       41359 : translate_basecase(long n, ulong p)
    4136             : {
    4137             : #ifdef LONG_IS_64BIT
    4138       36066 :   if (p <= 19) return n < 40;
    4139       29910 :   if (p < 1UL<<30) return n < 58;
    4140           0 :   if (p < 1UL<<59) return n < 100;
    4141           0 :   if (p < 1UL<<62) return n < 120;
    4142           0 :   if (p < 1UL<<63) return n < 240;
    4143           0 :   return n < 250;
    4144             : #else
    4145        5293 :   if (p <= 13) return n < 18;
    4146        4136 :   if (p <= 17) return n < 22;
    4147        4078 :   if (p <= 29) return n < 39;
    4148        3886 :   if (p <= 67) return n < 69;
    4149        3667 :   if (p < 1UL<< 15) return n < 80;
    4150        2047 :   if (p < 1UL<< 16) return n < 100;
    4151           0 :   if (p < 1UL<< 28) return n < 300;
    4152           0 :   return n < 650;
    4153             : #endif
    4154             : }
    4155             : /* assume p prime */
    4156             : GEN
    4157       16100 : Flx_translate1(GEN P, ulong p)
    4158             : {
    4159       16100 :   long d, n = degpol(P);
    4160             :   GEN R, Q, S;
    4161       16100 :   if (translate_basecase(n, p)) return Flx_translate1_basecase(P, p);
    4162             :   /* n > 0 */
    4163        1148 :   d = n >> 1;
    4164        1148 :   if ((ulong)n < p)
    4165             :   {
    4166           0 :     R = Flx_translate1(Flxn_red(P, d), p);
    4167           0 :     Q = Flx_translate1(Flx_shift(P, -d), p);
    4168           0 :     S = Fl_Xp1_powu(d, p, P[1]);
    4169           0 :     return Flx_add(Flx_mul(Q, S, p), R, p);
    4170             :   }
    4171             :   else
    4172             :   {
    4173             :     ulong q;
    4174        1148 :     if ((ulong)d > p) (void)ulogintall(d, p, &q); else q = p;
    4175        1148 :     R = Flx_translate1(Flxn_red(P, q), p);
    4176        1148 :     Q = Flx_translate1(Flx_shift(P, -q), p);
    4177        1148 :     S = Flx_add(Flx_shift(Q, q), Q, p);
    4178        1148 :     return Flx_add(S, R, p); /* P(x+1) = Q(x+1) (x^q+1) + R(x+1) */
    4179             :   }
    4180             : }
    4181             : 
    4182             : static GEN
    4183       12017 : zl_Xp1_powu(ulong n, ulong p, ulong q, long e, long vs)
    4184             : {
    4185       12017 :   ulong k, d = n >> 1, c, v = 0;
    4186       12017 :   GEN C, V, W, U = upowers(p, e-1);
    4187       12017 :   init_invlaplace(d, p, &V, &W);
    4188       12017 :   Flv_inv_inplace(V, q);
    4189       12017 :   C = cgetg(n+3, t_VECSMALL);
    4190       12017 :   C[1] = vs;
    4191       12017 :   uel(C,2) = 1UL;
    4192       12017 :   uel(C,3) = n%q;
    4193       12017 :   v = u_lvalrem(n, p, &c);
    4194     1355682 :   for (k = 2; k <= d; k++)
    4195             :   {
    4196             :     ulong w;
    4197     1343665 :     v += u_lvalrem(n-k+1, p, &w) - W[k];
    4198     1343665 :     c = Fl_mul(Fl_mul(w%q, c, q), uel(V,k), q);
    4199     1343665 :     uel(C,2+k) = v >= (ulong)e ? 0: v==0 ? c : Fl_mul(c, uel(U, v+1), q);
    4200             :   }
    4201     1374521 :   for (   ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
    4202       12017 :   return C; /* normalized */
    4203             : }
    4204             : 
    4205             : GEN
    4206       25259 : zlx_translate1(GEN P, ulong p, long e)
    4207             : {
    4208       25259 :   ulong d, q = upowuu(p,e), n = degpol(P);
    4209             :   GEN R, Q, S;
    4210       25259 :   if (translate_basecase(n, q)) return Flx_translate1_basecase(P, q);
    4211             :   /* n > 0 */
    4212       12017 :   d = n >> 1;
    4213       12017 :   R = zlx_translate1(Flxn_red(P, d), p, e);
    4214       12017 :   Q = zlx_translate1(Flx_shift(P, -d), p, e);
    4215       12017 :   S = zl_Xp1_powu(d, p, q, e, P[1]);
    4216       12017 :   return Flx_add(Flx_mul(Q, S, q), R, q);
    4217             : }
    4218             : 
    4219             : /***********************************************************************/
    4220             : /**                               Fl2                                 **/
    4221             : /***********************************************************************/
    4222             : /* Fl2 objects are Flv of length 2 [a,b] representing a+bsqrt(D) for
    4223             :  * a nonsquare D. */
    4224             : 
    4225             : INLINE GEN
    4226     7146946 : mkF2(ulong a, ulong b) { return mkvecsmall2(a,b); }
    4227             : 
    4228             : /* allow pi = 0 */
    4229             : GEN
    4230     1905026 : Fl2_mul_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
    4231             : {
    4232             :   ulong xaya, xbyb, Db2, mid, z1, z2;
    4233     1905026 :   ulong x1 = x[1], x2 = x[2], y1 = y[1], y2 = y[2];
    4234     1905026 :   if (pi)
    4235             :   {
    4236     1905355 :     xaya = Fl_mul_pre(x1,y1,p,pi);
    4237     1905673 :     if (x2==0 && y2==0) return mkF2(xaya,0);
    4238     1837340 :     if (x2==0) return mkF2(xaya,Fl_mul_pre(x1,y2,p,pi));
    4239     1812941 :     if (y2==0) return mkF2(xaya,Fl_mul_pre(x2,y1,p,pi));
    4240     1812639 :     xbyb = Fl_mul_pre(x2,y2,p,pi);
    4241     1812539 :     mid = Fl_mul_pre(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p,pi);
    4242     1812707 :     Db2 = Fl_mul_pre(D, xbyb, p,pi);
    4243             :   }
    4244           0 :   else if (p & HIGHMASK)
    4245             :   {
    4246           0 :     xaya = Fl_mul(x1,y1,p);
    4247           0 :     if (x2==0 && y2==0) return mkF2(xaya,0);
    4248           0 :     if (x2==0) return mkF2(xaya,Fl_mul(x1,y2,p));
    4249           0 :     if (y2==0) return mkF2(xaya,Fl_mul(x2,y1,p));
    4250           0 :     xbyb = Fl_mul(x2,y2,p);
    4251           0 :     mid = Fl_mul(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p);
    4252           0 :     Db2 = Fl_mul(D, xbyb, p);
    4253             :   }
    4254             :   else
    4255             :   {
    4256           0 :     xaya = (x1 * y1) % p;
    4257           0 :     if (x2==0 && y2==0) return mkF2(xaya,0);
    4258           0 :     if (x2==0) return mkF2(xaya, (x1 * y2) % p);
    4259           0 :     if (y2==0) return mkF2(xaya, (x2 * y1) % p);
    4260           0 :     xbyb = (x2 * y2) % p;
    4261           0 :     mid = (Fl_add(x1,x2,p) * Fl_add(y1,y2,p)) % p;
    4262           0 :     Db2 = (D * xbyb) % p;
    4263             :   }
    4264     1812606 :   z1 = Fl_add(xaya,Db2,p);
    4265     1812527 :   z2 = Fl_sub(mid,Fl_add(xaya,xbyb,p),p);
    4266     1812495 :   return mkF2(z1,z2);
    4267             : }
    4268             : 
    4269             : /* allow pi = 0 */
    4270             : GEN
    4271     4791931 : Fl2_sqr_pre(GEN x, ulong D, ulong p, ulong pi)
    4272             : {
    4273     4791931 :   ulong a = x[1], b = x[2];
    4274             :   ulong a2, Db2, ab;
    4275     4791931 :   if (pi)
    4276             :   {
    4277     4792552 :     a2 = Fl_sqr_pre(a,p,pi);
    4278     4795295 :     if (b==0) return mkF2(a2,0);
    4279     4586126 :     Db2= Fl_mul_pre(D, Fl_sqr_pre(b,p,pi), p,pi);
    4280     4586331 :     ab = Fl_mul_pre(a,b,p,pi);
    4281             :   }
    4282           0 :   else if (p & HIGHMASK)
    4283             :   {
    4284           0 :     a2 = Fl_sqr(a,p);
    4285           0 :     if (b==0) return mkF2(a2,0);
    4286           0 :     Db2= Fl_mul(D, Fl_sqr(b,p), p);
    4287           0 :     ab = Fl_mul(a,b,p);
    4288             :   }
    4289             :   else
    4290             :   {
    4291           0 :     a2 = (a * a) % p;
    4292           0 :     if (b==0) return mkF2(a2,0);
    4293           0 :     Db2= (D * ((b * b) % p)) % p;
    4294           0 :     ab = (a * b) % p;
    4295             :   }
    4296     4585548 :   return mkF2(Fl_add(a2,Db2,p), Fl_double(ab,p));
    4297             : }
    4298             : 
    4299             : /* allow pi = 0 */
    4300             : ulong
    4301      125563 : Fl2_norm_pre(GEN x, ulong D, ulong p, ulong pi)
    4302             : {
    4303      125563 :   ulong a = x[1], b = x[2], a2;
    4304      125563 :   if (pi)
    4305             :   {
    4306       71778 :     a2 = Fl_sqr_pre(a,p,pi);
    4307       71778 :     return b? Fl_sub(a2, Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p): a2;
    4308             :   }
    4309       53785 :   else if (p & HIGHMASK)
    4310             :   {
    4311           0 :     a2 = Fl_sqr(a,p);
    4312           0 :     return b? Fl_sub(a2, Fl_mul(D, Fl_sqr(b, p), p), p): a2;
    4313             :   }
    4314             :   else
    4315             :   {
    4316       53785 :     a2 = (a * a) % p;
    4317       53785 :     return b? Fl_sub(a2, (D * ((b * b) % p)) % p, p): a2;
    4318             :   }
    4319             : }
    4320             : 
    4321             : /* allow pi = 0 */
    4322             : GEN
    4323      192181 : Fl2_inv_pre(GEN x, ulong D, ulong p, ulong pi)
    4324             : {
    4325      192181 :   ulong a = x[1], b = x[2], n, ni;
    4326      192181 :   if (b == 0) return mkF2(Fl_inv(a,p), 0);
    4327      161633 :   b = Fl_neg(b, p);
    4328      161635 :   if (pi)
    4329             :   {
    4330      161635 :     n = Fl_sub(Fl_sqr_pre(a, p,pi),
    4331             :                Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p);
    4332      161634 :     ni = Fl_inv(n,p);
    4333      161630 :     return mkF2(Fl_mul_pre(a, ni, p,pi), Fl_mul_pre(b, ni, p,pi));
    4334             :   }
    4335           0 :   else if (p & HIGHMASK)
    4336             :   {
    4337           0 :     n = Fl_sub(Fl_sqr(a, p), Fl_mul(D, Fl_sqr(b, p), p), p);
    4338           0 :     ni = Fl_inv(n,p);
    4339           0 :     return mkF2(Fl_mul(a, ni, p), Fl_mul(b, ni, p));
    4340             :   }
    4341             :   else
    4342             :   {
    4343           0 :     n = Fl_sub((a * a) % p, (D * ((b * b) % p)) % p, p);
    4344           0 :     ni = Fl_inv(n,p);
    4345           0 :     return mkF2((a * ni) % p, (b * ni) % p);
    4346             :   }
    4347             : }
    4348             : 
    4349             : int
    4350      437435 : Fl2_equal1(GEN x) { return x[1]==1 && x[2]==0; }
    4351             : 
    4352             : struct _Fl2 {
    4353             :   ulong p, pi, D;
    4354             : };
    4355             : 
    4356             : static GEN
    4357     4791374 : _Fl2_sqr(void *data, GEN x)
    4358             : {
    4359     4791374 :   struct _Fl2 *D = (struct _Fl2*)data;
    4360     4791374 :   return Fl2_sqr_pre(x, D->D, D->p, D->pi);
    4361             : }
    4362             : static GEN
    4363     1876753 : _Fl2_mul(void *data, GEN x, GEN y)
    4364             : {
    4365     1876753 :   struct _Fl2 *D = (struct _Fl2*)data;
    4366     1876753 :   return Fl2_mul_pre(x,y, D->D, D->p, D->pi);
    4367             : }
    4368             : 
    4369             : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL; allow pi = 0 */
    4370             : GEN
    4371      653762 : Fl2_pow_pre(GEN x, GEN n, ulong D, ulong p, ulong pi)
    4372             : {
    4373      653762 :   pari_sp av = avma;
    4374             :   struct _Fl2 d;
    4375             :   GEN y;
    4376      653762 :   long s = signe(n);
    4377      653762 :   if (!s) return mkF2(1,0);
    4378      579572 :   if (s < 0)
    4379      192181 :     x = Fl2_inv_pre(x,D,p,pi);
    4380      579561 :   if (is_pm1(n)) return s < 0 ? x : zv_copy(x);
    4381      426694 :   d.p = p; d.pi = pi; d.D=D;
    4382      426694 :   y = gen_pow_i(x, n, (void*)&d, &_Fl2_sqr, &_Fl2_mul);
    4383      426697 :   return gerepileuptoleaf(av, y);
    4384             : }
    4385             : 
    4386             : static GEN
    4387      653747 : _Fl2_pow(void *data, GEN x, GEN n)
    4388             : {
    4389      653747 :   struct _Fl2 *D = (struct _Fl2*)data;
    4390      653747 :   return Fl2_pow_pre(x, n, D->D, D->p, D->pi);
    4391             : }
    4392             : 
    4393             : static GEN
    4394      110072 : _Fl2_rand(void *data)
    4395             : {
    4396      110072 :   struct _Fl2 *D = (struct _Fl2*)data;
    4397      110072 :   ulong a = random_Fl(D->p), b=random_Fl(D->p-1)+1;
    4398      110074 :   return mkF2(a,b);
    4399             : }
    4400             : 
    4401             : GEN
    4402       67856 : Fl2_sqrt_pre(GEN z, ulong D, ulong p, ulong pi)
    4403             : {
    4404       67856 :   ulong a = uel(z,1), b = uel(z,2), as2, u, v, s;
    4405       67856 :   ulong y = Fl_2gener_pre_i(D, p, pi);
    4406       67856 :   if (b == 0)
    4407       19504 :     return krouu(a, p)==1 ? mkF2(Fl_sqrt_pre_i(a, y, p, pi), 0)
    4408       19504 :                           : mkF2(0, Fl_sqrt_pre_i(Fl_div(a, D, p), y, p, pi));
    4409       54615 :   s = Fl_sqrt_pre_i(Fl2_norm_pre(z, D, p, pi), y, p, pi);
    4410       54615 :   if (s==~0UL) return NULL;
    4411       51441 :   as2 = Fl_halve(Fl_add(a, s, p), p);
    4412       51441 :   if (krouu(as2, p)==-1) as2 = Fl_sub(as2, s, p);
    4413       51441 :   u = Fl_sqrt_pre_i(as2, y, p, pi);
    4414       51441 :   v = Fl_div(b, Fl_double(u, p), p);
    4415       51441 :   return mkF2(u,v);
    4416             : }
    4417             : 
    4418             : static const struct bb_group Fl2_star={_Fl2_mul, _Fl2_pow, _Fl2_rand,
    4419             :        hash_GEN, zv_equal, Fl2_equal1, NULL};
    4420             : 
    4421             : /* allow pi = 0 */
    4422             : GEN
    4423       74189 : Fl2_sqrtn_pre(GEN a, GEN n, ulong D, ulong p, ulong pi, GEN *zeta)
    4424             : {
    4425             :   struct _Fl2 E;
    4426             :   GEN o;
    4427       74189 :   if (a[1]==0 && a[2]==0)
    4428             :   {
    4429           0 :     if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
    4430           0 :     if (zeta) *zeta=mkF2(1,0);
    4431           0 :     return zv_copy(a);
    4432             :   }
    4433       74189 :   E.p=p; E.pi = pi; E.D = D;
    4434       74189 :   o = subiu(powuu(p,2), 1);
    4435       74188 :   return gen_Shanks_sqrtn(a,n,o,zeta,(void*)&E,&Fl2_star);
    4436             : }
    4437             : 
    4438             : /* allow pi = 0 */
    4439             : GEN
    4440       10402 : Flx_Fl2_eval_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
    4441             : {
    4442             :   GEN p1;
    4443       10402 :   long i = lg(x)-1;
    4444       10402 :   if (i <= 2)
    4445        2065 :     return mkF2(i == 2? x[2]: 0, 0);
    4446        8337 :   p1 = mkF2(x[i], 0);
    4447       36456 :   for (i--; i>=2; i--)
    4448             :   {
    4449       28119 :     p1 = Fl2_mul_pre(p1, y, D, p, pi);
    4450       28119 :     uel(p1,1) = Fl_add(uel(p1,1), uel(x,i), p);
    4451             :   }
    4452        8337 :   return p1;
    4453             : }
    4454             : 
    4455             : /***********************************************************************/
    4456             : /**                               FlxV                                **/
    4457             : /***********************************************************************/
    4458             : /* FlxV are t_VEC with Flx coefficients. */
    4459             : 
    4460             : GEN
    4461       34482 : FlxV_Flc_mul(GEN V, GEN W, ulong p)
    4462             : {
    4463       34482 :   pari_sp ltop=avma;
    4464             :   long i;
    4465       34482 :   GEN z = Flx_Fl_mul(gel(V,1),W[1],p);
    4466      257068 :   for(i=2;i<lg(V);i++)
    4467      222586 :     z=Flx_add(z,Flx_Fl_mul(gel(V,i),W[i],p),p);
    4468       34482 :   return gerepileuptoleaf(ltop,z);
    4469             : }
    4470             : 
    4471             : GEN
    4472           0 : ZXV_to_FlxV(GEN x, ulong p)
    4473           0 : { pari_APPLY_type(t_VEC, ZX_to_Flx(gel(x,i), p)) }
    4474             : 
    4475             : GEN
    4476     3737392 : ZXT_to_FlxT(GEN x, ulong p)
    4477             : {
    4478     3737392 :   if (typ(x) == t_POL)
    4479     3677485 :     return ZX_to_Flx(x, p);
    4480             :   else
    4481      196523 :     pari_APPLY_type(t_VEC, ZXT_to_FlxT(gel(x,i), p))
    4482             : }
    4483             : 
    4484             : GEN
    4485      168781 : FlxV_to_Flm(GEN x, long n)
    4486      916895 : { pari_APPLY_type(t_MAT, Flx_to_Flv(gel(x,i), n)) }
    4487             : 
    4488             : GEN
    4489           0 : FlxV_red(GEN x, ulong p)
    4490           0 : { pari_APPLY_type(t_VEC, Flx_red(gel(x,i), p)) }
    4491             : 
    4492             : GEN
    4493      308194 : FlxT_red(GEN x, ulong p)
    4494             : {
    4495      308194 :   if (typ(x) == t_VECSMALL)
    4496      207371 :     return Flx_red(x, p);
    4497             :   else
    4498      338252 :     pari_APPLY_type(t_VEC, FlxT_red(gel(x,i), p))
    4499             : }
    4500             : 
    4501             : GEN
    4502      113589 : FlxqV_dotproduct_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
    4503             : {
    4504      113589 :   long i, lx = lg(x);
    4505             :   pari_sp av;
    4506             :   GEN c;
    4507      113589 :   if (lx == 1) return pol0_Flx(get_Flx_var(T));
    4508      113589 :   av = avma; c = Flx_mul_pre(gel(x,1),gel(y,1), p, pi);
    4509      464499 :   for (i=2; i<lx; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
    4510      113589 :   return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
    4511             : }
    4512             : GEN
    4513           0 : FlxqV_dotproduct(GEN x, GEN y, GEN T, ulong p)
    4514           0 : { return FlxqV_dotproduct_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    4515             : 
    4516             : GEN
    4517        2336 : FlxqX_dotproduct(GEN x, GEN y, GEN T, ulong p)
    4518             : {
    4519        2336 :   long i, l = minss(lg(x), lg(y));
    4520             :   ulong pi;
    4521             :   pari_sp av;
    4522             :   GEN c;
    4523        2336 :   if (l == 2) return pol0_Flx(get_Flx_var(T));
    4524        2265 :   av = avma; pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    4525        2265 :   c = Flx_mul_pre(gel(x,2),gel(y,2), p, pi);
    4526        6830 :   for (i=3; i<l; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
    4527        2265 :   return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
    4528             : }
    4529             : 
    4530             : /* allow pi = 0 */
    4531             : GEN
    4532      250783 : FlxC_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
    4533             : {
    4534      250783 :   long i, l = lg(z);
    4535      250783 :   GEN y = cgetg(l, t_VECSMALL);
    4536    12740263 :   for (i=1; i<l; i++) uel(y,i) = Flx_eval_powers_pre(gel(z,i), x, p, pi);
    4537      250795 :   return y;
    4538             : }
    4539             : 
    4540             : /***********************************************************************/
    4541             : /**                               FlxM                                **/
    4542             : /***********************************************************************/
    4543             : /* allow pi = 0 */
    4544             : GEN
    4545       19404 : FlxM_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
    4546             : {
    4547       19404 :   long i, l = lg(z);
    4548       19404 :   GEN y = cgetg(l, t_MAT);
    4549      270187 :   for (i=1; i<l; i++) gel(y,i) = FlxC_eval_powers_pre(gel(z,i), x, p, pi);
    4550       19404 :   return y;
    4551             : }
    4552             : 
    4553             : GEN
    4554           0 : zero_FlxC(long n, long sv)
    4555             : {
    4556           0 :   GEN x = cgetg(n + 1, t_COL), z = zero_Flx(sv);
    4557             :   long i;
    4558           0 :   for (i = 1; i <= n; i++) gel(x, i) = z;
    4559           0 :   return x;
    4560             : }
    4561             : 
    4562             : GEN
    4563           0 : FlxC_neg(GEN x, ulong p)
    4564           0 : { pari_APPLY_type(t_COL, Flx_neg(gel(x, i), p)) }
    4565             : 
    4566             : GEN
    4567           0 : FlxC_sub(GEN x, GEN y, ulong p)
    4568           0 : { pari_APPLY_type(t_COL, Flx_sub(gel(x, i), gel(y, i), p)) }
    4569             : 
    4570             : GEN
    4571           0 : zero_FlxM(long r, long c, long sv)
    4572             : {
    4573           0 :   GEN x = cgetg(c + 1, t_MAT), z = zero_FlxC(r, sv);
    4574             :   long j;
    4575           0 :   for (j = 1; j <= c; j++) gel(x, j) = z;
    4576           0 :   return x;
    4577             : }
    4578             : 
    4579             : GEN
    4580           0 : FlxM_neg(GEN x, ulong p)
    4581           0 : { pari_APPLY_same(FlxC_neg(gel(x, i), p)) }
    4582             : 
    4583             : GEN
    4584           0 : FlxM_sub(GEN x, GEN y, ulong p)
    4585           0 : { pari_APPLY_same(FlxC_sub(gel(x, i), gel(y,i), p)) }
    4586             : 
    4587             : GEN
    4588           0 : FlxqC_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
    4589           0 : { pari_APPLY_type(t_COL, Flxq_mul(gel(x, i), y, T, p)) }
    4590             : 
    4591             : GEN
    4592           0 : FlxqM_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
    4593           0 : { pari_APPLY_same(FlxqC_Flxq_mul(gel(x, i), y, T, p)) }
    4594             : 
    4595             : static GEN
    4596       58112 : FlxM_pack_ZM(GEN M, GEN (*pack)(GEN, long)) {
    4597             :   long i, j, l, lc;
    4598       58112 :   GEN N = cgetg_copy(M, &l), x;
    4599       58112 :   if (l == 1)
    4600           0 :     return N;
    4601       58112 :   lc = lgcols(M);
    4602      278301 :   for (j = 1; j < l; j++) {
    4603      220189 :     gel(N, j) = cgetg(lc, t_COL);
    4604     1131044 :     for (i = 1; i < lc; i++) {
    4605      910855 :       x = gcoeff(M, i, j);
    4606      910855 :       gcoeff(N, i, j) = pack(x + 2, lgpol(x));
    4607             :     }
    4608             :   }
    4609       58112 :   return N;
    4610             : }
    4611             : 
    4612             : static GEN
    4613      853644 : kron_pack_Flx_spec_half(GEN x, long l) {
    4614      853644 :   if (l == 0) return gen_0;
    4615      494952 :   return Flx_to_int_halfspec(x, l);
    4616             : }
    4617             : 
    4618             : static GEN
    4619       53822 : kron_pack_Flx_spec(GEN x, long l) {
    4620             :   long i;
    4621             :   GEN w, y;
    4622       53822 :   if (l == 0)
    4623       10075 :     return gen_0;
    4624       43747 :   y = cgetipos(l + 2);
    4625      159446 :   for (i = 0, w = int_LSW(y); i < l; i++, w = int_nextW(w))
    4626      115699 :     *w = x[i];
    4627       43747 :   return y;
    4628             : }
    4629             : 
    4630             : static GEN
    4631        3389 : kron_pack_Flx_spec_2(GEN x, long l) { return Flx_eval2BILspec(x, 2, l); }
    4632             : 
    4633             : static GEN
    4634           0 : kron_pack_Flx_spec_3(GEN x, long l) { return Flx_eval2BILspec(x, 3, l); }
    4635             : 
    4636             : static GEN
    4637       42982 : kron_unpack_Flx(GEN z, ulong p)
    4638             : {
    4639       42982 :   long i, l = lgefint(z);
    4640       42982 :   GEN x = cgetg(l, t_VECSMALL), w;
    4641      202068 :   for (w = int_LSW(z), i = 2; i < l; w = int_nextW(w), i++)
    4642      159086 :     x[i] = ((ulong) *w) % p;
    4643       42982 :   return Flx_renormalize(x, l);
    4644             : }
    4645             : 
    4646             : static GEN
    4647        2930 : kron_unpack_Flx_2(GEN x, ulong p) {
    4648        2930 :   long d = (lgefint(x)-1)/2 - 1;
    4649        2930 :   return Z_mod2BIL_Flx_2(x, d, p);
    4650             : }
    4651             : 
    4652             : static GEN
    4653           0 : kron_unpack_Flx_3(GEN x, ulong p) {
    4654           0 :   long d = lgefint(x)/3 - 1;
    4655           0 :   return Z_mod2BIL_Flx_3(x, d, p);
    4656             : }
    4657             : 
    4658             : static GEN
    4659      118719 : FlxM_pack_ZM_bits(GEN M, long b)
    4660             : {
    4661             :   long i, j, l, lc;
    4662      118719 :   GEN N = cgetg_copy(M, &l), x;
    4663      118719 :   if (l == 1)
    4664           0 :     return N;
    4665      118719 :   lc = lgcols(M);
    4666      491688 :   for (j = 1; j < l; j++) {
    4667      372969 :     gel(N, j) = cgetg(lc, t_COL);
    4668     5983078 :     for (i = 1; i < lc; i++) {
    4669     5610109 :       x = gcoeff(M, i, j);
    4670     5610109 :       gcoeff(N, i, j) = kron_pack_Flx_spec_bits(x + 2, b, lgpol(x));
    4671             :     }
    4672             :   }
    4673      118719 :   return N;
    4674             : }
    4675             : 
    4676             : static GEN
    4677       29056 : ZM_unpack_FlxqM(GEN M, GEN T, ulong p, ulong pi, GEN (*unpack)(GEN, ulong))
    4678             : {
    4679       29056 :   long i, j, l, lc, sv = get_Flx_var(T);
    4680       29056 :   GEN N = cgetg_copy(M, &l), x;
    4681       29056 :   if (l == 1)
    4682           0 :     return N;
    4683       29056 :   lc = lgcols(M);
    4684      163328 :   for (j = 1; j < l; j++) {
    4685      134272 :     gel(N, j) = cgetg(lc, t_COL);
    4686      767006 :     for (i = 1; i < lc; i++) {
    4687      632734 :       x = unpack(gcoeff(M, i, j), p);
    4688      632734 :       x[1] = sv;
    4689      632734 :       gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
    4690             :     }
    4691             :   }
    4692       29056 :   return N;
    4693             : }
    4694             : 
    4695             : static GEN
    4696       59400 : ZM_unpack_FlxqM_bits(GEN M, long b, GEN T, ulong p, ulong pi)
    4697             : {
    4698       59400 :   long i, j, l, lc, sv = get_Flx_var(T);
    4699       59400 :   GEN N = cgetg_copy(M, &l), x;
    4700       59400 :   if (l == 1)
    4701           0 :     return N;
    4702       59400 :   lc = lgcols(M);
    4703       59400 :   if (b < BITS_IN_LONG) {
    4704      200902 :     for (j = 1; j < l; j++) {
    4705      143161 :       gel(N, j) = cgetg(lc, t_COL);
    4706     3235607 :       for (i = 1; i < lc; i++) {
    4707     3092446 :         x = kron_unpack_Flx_bits_narrow(gcoeff(M, i, j), b, p);
    4708     3092446 :         x[1] = sv;
    4709     3092446 :         gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
    4710             :       }
    4711             :     }
    4712             :   } else {
    4713        1659 :     ulong pi = get_Fl_red(p);
    4714        9796 :     for (j = 1; j < l; j++) {
    4715        8137 :       gel(N, j) = cgetg(lc, t_COL);
    4716      175265 :       for (i = 1; i < lc; i++) {
    4717      167128 :         x = kron_unpack_Flx_bits_wide(gcoeff(M, i, j), b, p, pi);
    4718      167128 :         x[1] = sv;
    4719      167128 :         gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
    4720             :       }
    4721             :     }
    4722             :   }
    4723       59400 :   return N;
    4724             : }
    4725             : 
    4726             : GEN
    4727       88456 : FlxqM_mul_Kronecker(GEN A, GEN B, GEN T, ulong p)
    4728             : {
    4729       88456 :   pari_sp av = avma;
    4730       88456 :   long b, d = get_Flx_degree(T), n = lg(A) - 1;
    4731             :   GEN C, D, z;
    4732             :   GEN (*pack)(GEN, long), (*unpack)(GEN, ulong);
    4733       88456 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    4734       88456 :   int is_sqr = A==B;
    4735             : 
    4736       88456 :   z = muliu(muliu(sqru(p - 1), d), n);
    4737       88456 :   b = expi(z) + 1;
    4738             :   /* only do expensive bit-packing if it saves at least 1 limb */
    4739       88456 :   if (b <= BITS_IN_HALFULONG)
    4740       84053 :   { if (nbits2nlong(d*b) == (d + 1)/2) b = BITS_IN_HALFULONG; }
    4741             :   else
    4742             :   {
    4743        4403 :     long l = lgefint(z) - 2;
    4744        4403 :     if (nbits2nlong(d*b) == d*l) b = l*BITS_IN_LONG;
    4745             :   }
    4746       88456 :   set_avma(av);
    4747             : 
    4748       88456 :   switch (b) {
    4749       27995 :   case BITS_IN_HALFULONG:
    4750       27995 :     pack = kron_pack_Flx_spec_half;
    4751       27995 :     unpack = int_to_Flx_half;
    4752       27995 :     break;
    4753        1012 :   case BITS_IN_LONG:
    4754        1012 :     pack = kron_pack_Flx_spec;
    4755        1012 :     unpack = kron_unpack_Flx;
    4756        1012 :     break;
    4757          49 :   case 2*BITS_IN_LONG:
    4758          49 :     pack = kron_pack_Flx_spec_2;
    4759          49 :     unpack = kron_unpack_Flx_2;
    4760          49 :     break;
    4761           0 :   case 3*BITS_IN_LONG:
    4762           0 :     pack = kron_pack_Flx_spec_3;
    4763           0 :     unpack = kron_unpack_Flx_3;
    4764           0 :     break;
    4765       59400 :   default:
    4766       59400 :     A = FlxM_pack_ZM_bits(A, b);
    4767       59400 :     B = is_sqr? A: FlxM_pack_ZM_bits(B, b);
    4768       59400 :     C = ZM_mul(A, B);
    4769       59400 :     D = ZM_unpack_FlxqM_bits(C, b, T, p, pi);
    4770       59400 :     return gerepilecopy(av, D);
    4771             :   }
    4772       29056 :   A = FlxM_pack_ZM(A, pack);
    4773       29056 :   B = is_sqr? A: FlxM_pack_ZM(B, pack);
    4774       29056 :   C = ZM_mul(A, B);
    4775       29056 :   D = ZM_unpack_FlxqM(C, T, p, pi, unpack);
    4776       29056 :   return gerepilecopy(av, D);
    4777             : }

Generated by: LCOV version 1.14