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 - FpX.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30005-fc14bb602a) Lines: 1761 1926 91.4 %
Date: 2025-02-18 09:22:46 Functions: 189 202 93.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2007  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 over Fp */
      19             : 
      20             : static GEN
      21    85525043 : get_FpX_red(GEN T, GEN *B)
      22             : {
      23    85525043 :   if (typ(T)!=t_VEC) { *B=NULL; return T; }
      24      860481 :   *B = gel(T,1); return gel(T,2);
      25             : }
      26             : 
      27             : /***********************************************************************/
      28             : /**                                                                   **/
      29             : /**                              FpX                                  **/
      30             : /**                                                                   **/
      31             : /***********************************************************************/
      32             : 
      33             : /* FpX are polynomials over Z/pZ represented as t_POL with
      34             :  * t_INT coefficients.
      35             :  * 1) Coefficients should belong to {0,...,p-1}, though nonreduced
      36             :  * coefficients should work but be slower.
      37             :  *
      38             :  * 2) p is not assumed to be prime, but it is assumed that impossible divisions
      39             :  *    will not happen.
      40             :  * 3) Theses functions let some garbage on the stack, but are gerepileupto
      41             :  * compatible.
      42             :  */
      43             : 
      44             : static ulong
      45    44518432 : to_Flx(GEN *P, GEN *Q, GEN p)
      46             : {
      47    44518432 :   ulong pp = uel(p,2);
      48    44518432 :   *P = ZX_to_Flx(*P, pp);
      49    44525122 :   if(Q) *Q = ZX_to_Flx(*Q, pp);
      50    44524329 :   return pp;
      51             : }
      52             : 
      53             : static ulong
      54     2136853 : to_Flxq(GEN *P, GEN *T, GEN p)
      55             : {
      56     2136853 :   ulong pp = uel(p,2);
      57     2136853 :   if (P) *P = ZX_to_Flx(*P, pp);
      58     2136864 :   *T = ZXT_to_FlxT(*T, pp); return pp;
      59             : }
      60             : 
      61             : GEN
      62        1726 : Z_to_FpX(GEN a, GEN p, long v)
      63             : {
      64        1726 :   pari_sp av = avma;
      65        1726 :   GEN z = cgetg(3, t_POL);
      66        1726 :   GEN x = modii(a, p);
      67        1726 :   if (!signe(x)) { set_avma(av); return pol_0(v); }
      68        1726 :   z[1] = evalsigne(1) | evalvarn(v);
      69        1726 :   gel(z,2) = x; return z;
      70             : }
      71             : 
      72             : /* z in Z[X], return lift(z * Mod(1,p)), normalized*/
      73             : GEN
      74    95587423 : FpX_red(GEN z, GEN p)
      75             : {
      76    95587423 :   long i, l = lg(z);
      77    95587423 :   GEN x = cgetg(l, t_POL);
      78  1404671339 :   for (i=2; i<l; i++) gel(x,i) = modii(gel(z,i),p);
      79    95200958 :   x[1] = z[1]; return FpX_renormalize(x,l);
      80             : }
      81             : 
      82             : GEN
      83      404469 : FpXV_red(GEN x, GEN p)
      84     1901610 : { pari_APPLY_type(t_VEC, FpX_red(gel(x,i), p)) }
      85             : 
      86             : GEN
      87     1663671 : FpXT_red(GEN x, GEN p)
      88             : {
      89     1663671 :   if (typ(x) == t_POL)
      90     1575739 :     return FpX_red(x, p);
      91             :   else
      92      391520 :     pari_APPLY_type(t_VEC, FpXT_red(gel(x,i), p))
      93             : }
      94             : 
      95             : GEN
      96     1841389 : FpX_normalize(GEN z, GEN p)
      97             : {
      98     1841389 :   GEN p1 = leading_coeff(z);
      99     1841394 :   if (lg(z) == 2 || equali1(p1)) return z;
     100      126170 :   return FpX_Fp_mul_to_monic(z, Fp_inv(p1,p), p);
     101             : }
     102             : 
     103             : GEN
     104      312733 : FpX_center(GEN T, GEN p, GEN pov2)
     105             : {
     106      312733 :   long i, l = lg(T);
     107      312733 :   GEN P = cgetg(l,t_POL);
     108     1415125 :   for(i=2; i<l; i++) gel(P,i) = Fp_center(gel(T,i), p, pov2);
     109      312739 :   P[1] = T[1]; return P;
     110             : }
     111             : GEN
     112     1242273 : FpX_center_i(GEN T, GEN p, GEN pov2)
     113             : {
     114     1242273 :   long i, l = lg(T);
     115     1242273 :   GEN P = cgetg(l,t_POL);
     116     5646633 :   for(i=2; i<l; i++) gel(P,i) = Fp_center_i(gel(T,i), p, pov2);
     117     1242287 :   P[1] = T[1]; return P;
     118             : }
     119             : 
     120             : GEN
     121    16056331 : FpX_add(GEN x,GEN y,GEN p)
     122             : {
     123    16056331 :   long lx = lg(x), ly = lg(y), i;
     124             :   GEN z;
     125    16056331 :   if (lx < ly) swapspec(x,y, lx,ly);
     126    16056331 :   z = cgetg(lx,t_POL); z[1] = x[1];
     127   147206956 :   for (i=2; i<ly; i++) gel(z,i) = Fp_add(gel(x,i),gel(y,i), p);
     128    34908734 :   for (   ; i<lx; i++) gel(z,i) = modii(gel(x,i), p);
     129    16056289 :   z = ZX_renormalize(z, lx);
     130    16056334 :   if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); return pol_0(varn(x)); }
     131    15738761 :   return z;
     132             : }
     133             : 
     134             : static GEN
     135       21506 : Fp_red_FpX(GEN x, GEN p, long v)
     136             : {
     137             :   GEN z;
     138       21506 :   if (!signe(x)) return pol_0(v);
     139       14595 :   z = cgetg(3, t_POL);
     140       14595 :   gel(z,2) = Fp_red(x,p);
     141       14595 :   z[1] = evalvarn(v);
     142       14595 :   return FpX_renormalize(z, 3);
     143             : }
     144             : 
     145             : static GEN
     146         935 : Fp_neg_FpX(GEN x, GEN p, long v)
     147             : {
     148             :   GEN z;
     149         935 :   if (!signe(x)) return pol_0(v);
     150         794 :   z = cgetg(3, t_POL);
     151         794 :   gel(z,2) = Fp_neg(x,p);
     152         794 :   z[1] = evalvarn(v);
     153         794 :   return FpX_renormalize(z, 3);
     154             : }
     155             : 
     156             : GEN
     157      908011 : FpX_Fp_add(GEN y,GEN x,GEN p)
     158             : {
     159      908011 :   long i, lz = lg(y);
     160             :   GEN z;
     161      908011 :   if (lz == 2) return Fp_red_FpX(x,p,varn(y));
     162      886505 :   z = cgetg(lz,t_POL); z[1] = y[1];
     163      886505 :   gel(z,2) = Fp_add(gel(y,2),x, p);
     164      886505 :   if (lz == 3) z = FpX_renormalize(z,lz);
     165             :   else
     166     2285985 :     for(i=3;i<lz;i++) gel(z,i) = icopy(gel(y,i));
     167      886505 :   return z;
     168             : }
     169             : GEN
     170           0 : FpX_Fp_add_shallow(GEN y,GEN x,GEN p)
     171             : {
     172           0 :   long i, lz = lg(y);
     173             :   GEN z;
     174           0 :   if (lz == 2) return scalar_ZX_shallow(x,varn(y));
     175           0 :   z = cgetg(lz,t_POL); z[1] = y[1];
     176           0 :   gel(z,2) = Fp_add(gel(y,2),x, p);
     177           0 :   if (lz == 3) z = FpX_renormalize(z,lz);
     178             :   else
     179           0 :     for(i=3;i<lz;i++) gel(z,i) = gel(y,i);
     180           0 :   return z;
     181             : }
     182             : GEN
     183      588244 : FpX_Fp_sub(GEN y,GEN x,GEN p)
     184             : {
     185      588244 :   long i, lz = lg(y);
     186             :   GEN z;
     187      588244 :   if (lz == 2) return Fp_neg_FpX(x,p,varn(y));
     188      587309 :   z = cgetg(lz,t_POL); z[1] = y[1];
     189      587309 :   gel(z,2) = Fp_sub(gel(y,2),x, p);
     190      587310 :   if (lz == 3) z = FpX_renormalize(z,lz);
     191             :   else
     192     1349859 :     for(i=3;i<lz;i++) gel(z,i) = icopy(gel(y,i));
     193      587310 :   return z;
     194             : }
     195             : GEN
     196       11146 : FpX_Fp_sub_shallow(GEN y,GEN x,GEN p)
     197             : {
     198       11146 :   long i, lz = lg(y);
     199             :   GEN z;
     200       11146 :   if (lz == 2) return Fp_neg_FpX(x,p,varn(y));
     201       11146 :   z = cgetg(lz,t_POL); z[1] = y[1];
     202       11146 :   gel(z,2) = Fp_sub(gel(y,2),x, p);
     203       11146 :   if (lz == 3) z = FpX_renormalize(z,lz);
     204             :   else
     205       37357 :     for(i=3;i<lz;i++) gel(z,i) = gel(y,i);
     206       11146 :   return z;
     207             : }
     208             : 
     209             : GEN
     210      467994 : FpX_neg(GEN x,GEN p)
     211     5011729 : { pari_APPLY_ZX(Fp_neg(gel(x,i), p)); }
     212             : 
     213             : static GEN
     214    15414396 : FpX_subspec(GEN x,GEN y,GEN p, long nx, long ny)
     215             : {
     216             :   long i, lz;
     217             :   GEN z;
     218    15414396 :   if (nx >= ny)
     219             :   {
     220    10923364 :     lz = nx+2;
     221    10923364 :     z = cgetg(lz,t_POL); z[1] = 0; z += 2;
     222   118373465 :     for (i=0; i<ny; i++) gel(z,i) = Fp_sub(gel(x,i),gel(y,i), p);
     223    11672337 :     for (   ; i<nx; i++) gel(z,i) = modii(gel(x,i), p);
     224             :   }
     225             :   else
     226             :   {
     227     4491032 :     lz = ny+2;
     228     4491032 :     z = cgetg(lz,t_POL); z[1] = 0; z += 2;
     229    23437079 :     for (i=0; i<nx; i++) gel(z,i) = Fp_sub(gel(x,i),gel(y,i), p);
     230    14762274 :     for (   ; i<ny; i++) gel(z,i) = Fp_neg(gel(y,i), p);
     231             :   }
     232    15411500 :   z = FpX_renormalize(z-2, lz);
     233    15414409 :   if (!lgpol(z)) { set_avma((pari_sp)(z + lz)); return pol_0(0); }
     234    15156343 :   return z;
     235             : }
     236             : 
     237             : GEN
     238    14598141 : FpX_sub(GEN x,GEN y,GEN p)
     239             : {
     240    14598141 :   GEN z = FpX_subspec(x+2,y+2,p,lgpol(x),lgpol(y));
     241    14598148 :   setvarn(z, varn(x));
     242    14598148 :   return z;
     243             : }
     244             : 
     245             : GEN
     246       25690 : Fp_FpX_sub(GEN x, GEN y, GEN p)
     247             : {
     248       25690 :   long ly = lg(y), i;
     249             :   GEN z;
     250       25690 :   if (ly <= 3) {
     251         482 :     z = cgetg(3, t_POL);
     252         482 :     x = (ly == 3)? Fp_sub(x, gel(y,2), p): modii(x, p);
     253         482 :     if (!signe(x)) { set_avma((pari_sp)(z + 3)); return pol_0(varn(y)); }
     254         399 :     z[1] = evalsigne(1)|y[1]; gel(z,2) = x; return z;
     255             :   }
     256       25208 :   z = cgetg(ly,t_POL);
     257       25208 :   gel(z,2) = Fp_sub(x, gel(y,2), p);
     258       93612 :   for (i = 3; i < ly; i++) gel(z,i) = Fp_neg(gel(y,i), p);
     259       25208 :   z = ZX_renormalize(z, ly);
     260       25208 :   if (!lgpol(z)) { set_avma((pari_sp)(z + ly)); return pol_0(varn(x)); }
     261       25208 :   z[1] = y[1]; return z;
     262             : }
     263             : 
     264             : GEN
     265         994 : FpX_convol(GEN x, GEN y, GEN p)
     266             : {
     267         994 :   long lx = lg(x), ly = lg(y), i;
     268             :   GEN z;
     269         994 :   if (lx < ly) swapspec(x,y, lx,ly);
     270         994 :   z = cgetg(ly,t_POL); z[1] = x[1];
     271       58751 :   for (i=2; i<ly; i++) gel(z,i) = Fp_mul(gel(x,i),gel(y,i), p);
     272         994 :   z = ZX_renormalize(z, ly);
     273         994 :   if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); return pol_0(varn(x)); }
     274         994 :   return z;
     275             : }
     276             : 
     277             : GEN
     278    26401263 : FpX_mul(GEN x,GEN y,GEN p)
     279             : {
     280    26401263 :   if (lgefint(p) == 3)
     281             :   {
     282    13651202 :     ulong pp = to_Flx(&x, &y, p);
     283    13652171 :     return Flx_to_ZX(Flx_mul(x, y, pp));
     284             :   }
     285    12750061 :   return FpX_red(ZX_mul(x, y), p);
     286             : }
     287             : 
     288             : GEN
     289     9202398 : FpX_mulspec(GEN a, GEN b, GEN p, long na, long nb)
     290     9202398 : { return FpX_red(ZX_mulspec(a, b, na, nb), p); }
     291             : 
     292             : GEN
     293     6565454 : FpX_sqr(GEN x,GEN p)
     294             : {
     295     6565454 :   if (lgefint(p) == 3)
     296             :   {
     297      355008 :     ulong pp = to_Flx(&x, NULL, p);
     298      355005 :     return Flx_to_ZX(Flx_sqr(x, pp));
     299             :   }
     300     6210446 :   return FpX_red(ZX_sqr(x), p);
     301             : }
     302             : 
     303             : GEN
     304     1060651 : FpX_mulu(GEN x, ulong t,GEN p)
     305             : {
     306     1060651 :   t = umodui(t, p); if (!t) return zeropol(varn(x));
     307     6637510 :   pari_APPLY_ZX(Fp_mulu(gel(x,i), t, p));
     308             : }
     309             : 
     310             : GEN
     311        8610 : FpX_divu(GEN y, ulong x, GEN p)
     312        8610 : { return FpX_Fp_div(y, utoi(umodui(x, p)), p); }
     313             : 
     314             : GEN
     315     5647536 : FpX_Fp_mulspec(GEN y,GEN x,GEN p,long ly)
     316             : {
     317             :   GEN z;
     318             :   long i;
     319     5647536 :   if (!signe(x)) return pol_0(0);
     320     5617654 :   z = cgetg(ly+2,t_POL); z[1] = evalsigne(1);
     321    33283912 :   for(i=0; i<ly; i++) gel(z,i+2) = Fp_mul(gel(y,i), x, p);
     322     5615593 :   return ZX_renormalize(z, ly+2);
     323             : }
     324             : 
     325             : GEN
     326     5632979 : FpX_Fp_mul(GEN y,GEN x,GEN p)
     327             : {
     328     5632979 :   GEN z = FpX_Fp_mulspec(y+2,x,p,lgpol(y));
     329     5632672 :   setvarn(z, varn(y)); return z;
     330             : }
     331             : 
     332             : GEN
     333      613071 : FpX_Fp_div(GEN y, GEN x, GEN p)
     334             : {
     335      613071 :   return FpX_Fp_mul(y, Fp_inv(x, p), p);
     336             : }
     337             : 
     338             : GEN
     339      134390 : FpX_Fp_mul_to_monic(GEN y,GEN x,GEN p)
     340             : {
     341             :   GEN z;
     342             :   long i, l;
     343      134390 :   z = cgetg_copy(y, &l); z[1] = y[1];
     344      577947 :   for(i=2; i<l-1; i++) gel(z,i) = Fp_mul(gel(y,i), x, p);
     345      134390 :   gel(z,l-1) = gen_1; return z;
     346             : }
     347             : 
     348             : struct _FpXQ {
     349             :   GEN T, p, aut;
     350             : };
     351             : 
     352             : struct _FpX
     353             : {
     354             :   GEN p;
     355             :   long v;
     356             : };
     357             : 
     358             : static GEN
     359      371535 : _FpX_mul(void* E, GEN x, GEN y)
     360      371535 : { struct _FpX *D = (struct _FpX *)E; return FpX_mul(x, y, D->p); }
     361             : static GEN
     362       86390 : _FpX_sqr(void *E, GEN x)
     363       86390 : { struct _FpX *D = (struct _FpX *)E; return FpX_sqr(x, D->p); }
     364             : 
     365             : GEN
     366      318508 : FpX_powu(GEN x, ulong n, GEN p)
     367             : {
     368             :   struct _FpX D;
     369      318508 :   if (n==0) return pol_1(varn(x));
     370       61720 :   D.p = p;
     371       61720 :   return gen_powu(x, n, (void *)&D, _FpX_sqr, _FpX_mul);
     372             : }
     373             : 
     374             : GEN
     375      309968 : FpXV_prod(GEN V, GEN p)
     376             : {
     377             :   struct _FpX D;
     378      309968 :   D.p = p;
     379      309968 :   return gen_product(V, (void *)&D, &_FpX_mul);
     380             : }
     381             : 
     382             : static GEN
     383       35659 : _FpX_pow(void* E, GEN x, GEN y)
     384       35659 : { struct _FpX *D = (struct _FpX *)E; return FpX_powu(x, itou(y), D->p); }
     385             : static GEN
     386           0 : _FpX_one(void *E)
     387           0 : { struct _FpX *D = (struct _FpX *)E; return pol_1(D->v); }
     388             : 
     389             : GEN
     390       23300 : FpXV_factorback(GEN f, GEN e, GEN p, long v)
     391             : {
     392             :   struct _FpX D;
     393       23300 :   D.p = p; D.v = v;
     394       23300 :   return gen_factorback(f, e, (void *)&D, &_FpX_mul, &_FpX_pow, &_FpX_one);
     395             : }
     396             : 
     397             : GEN
     398       94911 : FpX_halve(GEN x, GEN p)
     399      281519 : { pari_APPLY_pol_normalized(Fp_halve(gel(x,i), p)); }
     400             : 
     401             : static GEN
     402    65775116 : FpX_divrem_basecase(GEN x, GEN y, GEN p, GEN *pr)
     403             : {
     404             :   long vx, dx, dy, dy1, dz, i, j, sx, lr;
     405             :   pari_sp av0, av;
     406             :   GEN z,p1,rem,lead;
     407             : 
     408    65775116 :   if (!signe(y)) pari_err_INV("FpX_divrem",y);
     409    65775116 :   vx = varn(x);
     410    65775116 :   dy = degpol(y);
     411    65774411 :   dx = degpol(x);
     412    65774897 :   if (dx < dy)
     413             :   {
     414      126102 :     if (pr)
     415             :     {
     416      125550 :       av0 = avma; x = FpX_red(x, p);
     417      125550 :       if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: pol_0(vx); }
     418      125550 :       if (pr == ONLY_REM) return x;
     419      125550 :       *pr = x;
     420             :     }
     421      126102 :     return pol_0(vx);
     422             :   }
     423    65648795 :   lead = leading_coeff(y);
     424    65651571 :   if (!dy) /* y is constant */
     425             :   {
     426      607609 :     if (pr && pr != ONLY_DIVIDES)
     427             :     {
     428      603874 :       if (pr == ONLY_REM) return pol_0(vx);
     429      585612 :       *pr = pol_0(vx);
     430             :     }
     431      589347 :     av0 = avma;
     432      589347 :     if (equali1(lead)) return FpX_red(x, p);
     433      584740 :     else return gerepileupto(av0, FpX_Fp_div(x, lead, p));
     434             :   }
     435    65043962 :   av0 = avma; dz = dx-dy;
     436    65043962 :   if (lgefint(p) == 3)
     437             :   { /* assume ab != 0 mod p */
     438    28293761 :     ulong pp = to_Flx(&x, &y, p);
     439    28298003 :     z = Flx_divrem(x, y, pp, pr);
     440    28285042 :     set_avma(av0); /* HACK: assume pr last on stack, then z */
     441    28283867 :     if (!z) return NULL;
     442    28283727 :     z = leafcopy(z);
     443    28284844 :     if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
     444             :     {
     445     5544118 :       *pr = leafcopy(*pr);
     446     5544176 :       *pr = Flx_to_ZX_inplace(*pr);
     447             :     }
     448    28284909 :     return Flx_to_ZX_inplace(z);
     449             :   }
     450    36750201 :   lead = equali1(lead)? NULL: gclone(Fp_inv(lead,p));
     451    36749940 :   set_avma(av0);
     452    36749952 :   z=cgetg(dz+3,t_POL); z[1] = x[1];
     453    36749944 :   x += 2; y += 2; z += 2;
     454    39905015 :   for (dy1=dy-1; dy1>=0 && !signe(gel(y, dy1)); dy1--);
     455             : 
     456    36749944 :   p1 = gel(x,dx); av = avma;
     457    36749944 :   gel(z,dz) = lead? gerepileuptoint(av, Fp_mul(p1,lead, p)): icopy(p1);
     458   111974999 :   for (i=dx-1; i>=dy; i--)
     459             :   {
     460    75224018 :     av=avma; p1=gel(x,i);
     461   959961222 :     for (j=i-dy1; j<=i && j<=dz; j++)
     462   884745516 :       p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
     463    75215706 :     if (lead) p1 = mulii(p1,lead);
     464    75215706 :     gel(z,i-dy) = gerepileuptoint(av,modii(p1, p));
     465             :   }
     466    36750981 :   if (!pr) { guncloneNULL(lead); return z-2; }
     467             : 
     468    36670627 :   rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
     469    40609100 :   for (sx=0; ; i--)
     470             :   {
     471    40609100 :     p1 = gel(x,i);
     472   224940973 :     for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
     473   184332649 :       p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
     474    40608240 :     p1 = modii(p1,p); if (signe(p1)) { sx = 1; break; }
     475     4075304 :     if (!i) break;
     476     3938469 :     set_avma(av);
     477             :   }
     478    36669142 :   if (pr == ONLY_DIVIDES)
     479             :   {
     480           0 :     guncloneNULL(lead);
     481           0 :     if (sx) return gc_NULL(av0);
     482           0 :     return gc_const((pari_sp)rem, z-2);
     483             :   }
     484    36669142 :   lr=i+3; rem -= lr;
     485    36669142 :   rem[0] = evaltyp(t_POL) | _evallg(lr);
     486    36669142 :   rem[1] = z[-1];
     487    36669142 :   p1 = gerepileuptoint((pari_sp)rem, p1);
     488    36670480 :   rem += 2; gel(rem,i) = p1;
     489   165745032 :   for (i--; i>=0; i--)
     490             :   {
     491   129074601 :     av=avma; p1 = gel(x,i);
     492  1088835790 :     for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
     493   959782006 :       p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
     494   129053535 :     gel(rem,i) = gerepileuptoint(av, modii(p1,p));
     495             :   }
     496    36670431 :   rem -= 2;
     497    36670431 :   guncloneNULL(lead);
     498    36670384 :   if (!sx) (void)FpX_renormalize(rem, lr);
     499    36670391 :   if (pr == ONLY_REM) return gerepileupto(av0,rem);
     500     2537131 :   *pr = rem; return z-2;
     501             : }
     502             : 
     503             : GEN
     504      165124 : FpX_div_by_X_x(GEN a, GEN x, GEN p, GEN *r)
     505             : {
     506      165124 :   long l = lg(a), i;
     507             :   GEN z;
     508      165124 :   if (l <= 3)
     509             :   {
     510           0 :     if (r) *r = l == 2? gen_0: icopy(gel(a,2));
     511           0 :     return pol_0(varn(a));
     512             :   }
     513      165124 :   l--; z = cgetg(l, t_POL); z[1] = a[1];
     514      165124 :   gel(z, l-1) = gel(a,l);
     515     2481992 :   for (i = l-2; i > 1; i--) /* z[i] = a[i+1] + x*z[i+1] */
     516     2316874 :     gel(z,i) = Fp_addmul(gel(a,i+1), x, gel(z,i+1), p);
     517      165118 :   if (r) *r = Fp_addmul(gel(a,2), x, gel(z,2), p);
     518      165118 :   return z;
     519             : }
     520             : 
     521             : static GEN
     522      134778 : _FpX_divrem(void * E, GEN x, GEN y, GEN *r)
     523             : {
     524      134778 :   struct _FpX *D = (struct _FpX*) E;
     525      134778 :   return FpX_divrem(x, y, D->p, r);
     526             : }
     527             : static GEN
     528       20062 : _FpX_add(void * E, GEN x, GEN y) {
     529       20062 :   struct _FpX *D = (struct _FpX*) E;
     530       20062 :   return FpX_add(x, y, D->p);
     531             : }
     532             : 
     533             : static struct bb_ring FpX_ring = { _FpX_add,_FpX_mul,_FpX_sqr };
     534             : 
     535             : GEN
     536       11403 : FpX_digits(GEN x, GEN T, GEN p)
     537             : {
     538             :   struct _FpX D;
     539       11403 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
     540       11403 :   D.p = p;
     541       11403 :   return gen_digits(x,T,n,(void *)&D, &FpX_ring, _FpX_divrem);
     542             : }
     543             : 
     544             : GEN
     545        4564 : FpXV_FpX_fromdigits(GEN x, GEN T, GEN p)
     546             : {
     547             :   struct _FpX D;
     548        4564 :   D.p = p;
     549        4564 :   return gen_fromdigits(x,T,(void *)&D, &FpX_ring);
     550             : }
     551             : 
     552             : long
     553      254647 : FpX_valrem(GEN x, GEN t, GEN p, GEN *py)
     554             : {
     555      254647 :   pari_sp av=avma;
     556             :   long k;
     557             :   GEN r, y;
     558             : 
     559      254647 :   for (k=0; ; k++)
     560             :   {
     561      650739 :     y = FpX_divrem(x, t, p, &r);
     562      650744 :     if (signe(r)) break;
     563      396092 :     x = y;
     564             :   }
     565      254652 :   *py = gerepilecopy(av,x);
     566      254655 :   return k;
     567             : }
     568             : 
     569             : static GEN
     570       87863 : FpX_addmulmul(GEN u, GEN v, GEN x, GEN y, GEN p)
     571             : {
     572       87863 :   return FpX_add(FpX_mul(u, x, p),FpX_mul(v, y, p), p);
     573             : }
     574             : 
     575             : static GEN
     576       36106 : FpXM_FpX_mul2(GEN M, GEN x, GEN y, GEN p)
     577             : {
     578       36106 :   GEN res = cgetg(3, t_COL);
     579       36106 :   gel(res, 1) = FpX_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p);
     580       36106 :   gel(res, 2) = FpX_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p);
     581       36106 :   return res;
     582             : }
     583             : 
     584             : static GEN
     585       17253 : FpXM_mul2(GEN A, GEN B, GEN p)
     586             : {
     587       17253 :   GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
     588       17253 :   GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
     589       17253 :   GEN M1 = FpX_mul(FpX_add(A11,A22, p), FpX_add(B11,B22, p), p);
     590       17253 :   GEN M2 = FpX_mul(FpX_add(A21,A22, p), B11, p);
     591       17253 :   GEN M3 = FpX_mul(A11, FpX_sub(B12,B22, p), p);
     592       17253 :   GEN M4 = FpX_mul(A22, FpX_sub(B21,B11, p), p);
     593       17253 :   GEN M5 = FpX_mul(FpX_add(A11,A12, p), B22, p);
     594       17253 :   GEN M6 = FpX_mul(FpX_sub(A21,A11, p), FpX_add(B11,B12, p), p);
     595       17253 :   GEN M7 = FpX_mul(FpX_sub(A12,A22, p), FpX_add(B21,B22, p), p);
     596       17253 :   GEN T1 = FpX_add(M1,M4, p), T2 = FpX_sub(M7,M5, p);
     597       17253 :   GEN T3 = FpX_sub(M1,M2, p), T4 = FpX_add(M3,M6, p);
     598       17253 :   retmkmat22(FpX_add(T1,T2, p), FpX_add(M3,M5, p),
     599             :              FpX_add(M2,M4, p), FpX_add(T3,T4, p));
     600             : }
     601             : 
     602             : /* Return [0,1;1,-q]*M */
     603             : static GEN
     604       17162 : FpX_FpXM_qmul(GEN q, GEN M, GEN p)
     605             : {
     606       17162 :   GEN u = FpX_mul(gcoeff(M,2,1), q, p);
     607       17162 :   GEN v = FpX_mul(gcoeff(M,2,2), q, p);
     608       17162 :   retmkmat22(gcoeff(M,2,1), gcoeff(M,2,2),
     609             :     FpX_sub(gcoeff(M,1,1), u, p), FpX_sub(gcoeff(M,1,2), v, p));
     610             : }
     611             : 
     612             : static GEN
     613          24 : matid2_FpXM(long v)
     614          24 : { retmkmat22(pol_1(v), pol_0(v), pol_0(v), pol_1(v)); }
     615             : 
     616             : static GEN
     617           8 : matJ2_FpXM(long v)
     618           8 : { retmkmat22(pol_0(v), pol_1(v), pol_1(v), pol_0(v)); }
     619             : 
     620             : INLINE GEN
     621      989051 : FpX_shift(GEN a, long n) { return RgX_shift_shallow(a, n); }
     622             : 
     623             : INLINE GEN
     624      205706 : FpXn_red(GEN a, long n) { return RgXn_red_shallow(a, n); }
     625             : 
     626             : /* Fast resultant formula from William Hart in Flint <http://flintlib.org/> */
     627             : 
     628             : struct FpX_res
     629             : {
     630             :    GEN res, lc;
     631             :    long deg0, deg1, off;
     632             : };
     633             : 
     634             : INLINE void
     635        3749 : FpX_halfres_update(long da, long db, long dr, GEN p, struct FpX_res *res)
     636             : {
     637        3749 :   if (dr >= 0)
     638             :   {
     639        3749 :     if (!equali1(res->lc))
     640             :     {
     641        3749 :       res->lc  = Fp_powu(res->lc, da - dr, p);
     642        3749 :       res->res = Fp_mul(res->res, res->lc, p);
     643             :     }
     644        3749 :     if (both_odd(da + res->off, db + res->off))
     645           0 :       res->res = Fp_neg(res->res, p);
     646             :   } else
     647             :   {
     648           0 :     if (db == 0)
     649             :     {
     650           0 :       if (!equali1(res->lc))
     651             :       {
     652           0 :           res->lc  = Fp_powu(res->lc, da, p);
     653           0 :           res->res = Fp_mul(res->res, res->lc, p);
     654             :       }
     655             :     } else
     656           0 :       res->res = gen_0;
     657             :   }
     658        3749 : }
     659             : 
     660             : static GEN
     661       33795 : FpX_halfres_basecase(GEN a, GEN b, GEN p, GEN *pa, GEN *pb, struct FpX_res *res)
     662             : {
     663       33795 :   pari_sp av=avma;
     664             :   GEN u,u1,v,v1, M;
     665       33795 :   long vx = varn(a), n = lgpol(a)>>1;
     666       33795 :   u1 = v = pol_0(vx);
     667       33795 :   u = v1 = pol_1(vx);
     668      473549 :   while (lgpol(b)>n)
     669             :   {
     670             :     GEN r, q;
     671      439754 :     q = FpX_divrem(a,b,p, &r);
     672      439754 :     if (res)
     673             :     {
     674        3625 :       long da = degpol(a), db=degpol(b), dr = degpol(r);
     675        3625 :       res->lc = leading_coeff(b);
     676        3625 :       if (dr >= n)
     677        3403 :         FpX_halfres_update(da,db,dr,p,res);
     678             :       else
     679             :       {
     680         222 :         res->deg0 = da;
     681         222 :         res->deg1 = db;
     682             :       }
     683             :     }
     684      439754 :     a = b; b = r; swap(u,u1); swap(v,v1);
     685      439754 :     u1 = FpX_sub(u1, FpX_mul(u, q, p), p);
     686      439754 :     v1 = FpX_sub(v1, FpX_mul(v, q, p), p);
     687      439754 :     if (gc_needed(av,2))
     688             :     {
     689           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpX_halfgcd (d = %ld)",degpol(b));
     690           0 :       if (res)
     691           0 :         gerepileall(av, 8, &a,&b,&u1,&v1,&u,&v,&res->res,&res->lc);
     692             :       else
     693           0 :         gerepileall(av, 6, &a,&b,&u1,&v1,&u,&v);
     694             :     }
     695             :   }
     696       33795 :   M = mkmat22(u,v,u1,v1); *pa = a; *pb = b;
     697         222 :   return res ? gc_all(av, 5, &M, pa, pb, &res->res, &res->lc)
     698       34017 :              : gc_all(av, 3, &M, pa, pb);
     699             : }
     700             : 
     701             : static GEN FpX_halfres_i(GEN x, GEN y, GEN p, GEN *a, GEN *b, struct FpX_res *res);
     702             : 
     703             : static GEN
     704       18960 : FpX_halfres_split(GEN x, GEN y, GEN p, GEN *a, GEN *b, struct FpX_res *res)
     705             : {
     706       18960 :   pari_sp av = avma;
     707             :   GEN R, S, T, V1, V2;
     708             :   GEN x1, y1, r, q;
     709       18960 :   long l = lgpol(x), n = l>>1, k;
     710       18960 :   if (lgpol(y) <= n)
     711           8 :     { *a = RgX_copy(x); *b = RgX_copy(y); return matid2_FpXM(varn(x)); }
     712       18952 :   if (res)
     713             :   {
     714         166 :      res->lc = leading_coeff(y);
     715         166 :      res->deg0 -= n;
     716         166 :      res->deg1 -= n;
     717         166 :      res->off += n;
     718             :   }
     719       18952 :   R = FpX_halfres_i(FpX_shift(x,-n), FpX_shift(y,-n), p, a, b, res);
     720       18952 :   if (res)
     721             :   {
     722         166 :     res->off -= n;
     723         166 :     res->deg0 += n;
     724         166 :     res->deg1 += n;
     725             :   }
     726       18952 :   V1 = FpXM_FpX_mul2(R, FpXn_red(x,n), FpXn_red(y,n), p);
     727       18952 :   x1 = FpX_add(FpX_shift(*a,n), gel(V1,1), p);
     728       18952 :   y1 = FpX_add(FpX_shift(*b,n), gel(V1,2), p);
     729       18952 :   if (lgpol(y1) <= n)
     730             :   {
     731        1798 :     *a = x1; *b = y1;
     732          42 :     return res ? gc_all(av, 5, &R, a, b, &res->res, &res->lc)
     733        1840 :                : gc_all(av, 3, &R, a, b);
     734             :   }
     735       17154 :   k = 2*n-degpol(y1);
     736       17154 :   q = FpX_divrem(x1, y1, p, &r);
     737       17154 :   if (res)
     738             :   {
     739         124 :     long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
     740         124 :     if (dy1 < degpol(y))
     741         116 :       FpX_halfres_update(res->deg0, res->deg1, dy1, p,res);
     742         124 :     res->lc = gel(y1, dy1+2);
     743         124 :     res->deg0 = dx1;
     744         124 :     res->deg1 = dy1;
     745         124 :     if (dr >= n)
     746             :     {
     747         124 :       FpX_halfres_update(dx1, dy1, dr, p,res);
     748         124 :       res->deg0 = dy1;
     749         124 :       res->deg1 = dr;
     750             :     }
     751         124 :     res->deg0 -= k;
     752         124 :     res->deg1 -= k;
     753         124 :     res->off += k;
     754             :   }
     755       17154 :   S = FpX_halfres_i(FpX_shift(y1,-k), FpX_shift(r,-k), p, a, b, res);
     756       17154 :   if (res)
     757             :   {
     758         124 :     res->deg0 += k;
     759         124 :     res->deg1 += k;
     760         124 :     res->off -= k;
     761             :   }
     762       17154 :   T = FpXM_mul2(S, FpX_FpXM_qmul(q, R, p), p);
     763       17154 :   V2 = FpXM_FpX_mul2(S, FpXn_red(y1,k), FpXn_red(r,k), p);
     764       17154 :   *a = FpX_add(FpX_shift(*a,k), gel(V2,1), p);
     765       17154 :   *b = FpX_add(FpX_shift(*b,k), gel(V2,2), p);
     766         124 :   return res ? gc_all(av, 5, &T, a, b, &res->res, &res->lc)
     767       17278 :              : gc_all(av, 3, &T, a, b);
     768             : }
     769             : 
     770             : static GEN
     771       52755 : FpX_halfres_i(GEN x, GEN y, GEN p, GEN *a, GEN *b, struct FpX_res *res)
     772             : {
     773       52755 :   if (lgpol(x) < FpX_HALFGCD_LIMIT)
     774       33795 :     return FpX_halfres_basecase(x, y, p, a, b, res);
     775       18960 :   return FpX_halfres_split(x, y, p, a, b, res);
     776             : }
     777             : 
     778             : static GEN
     779       16543 : FpX_halfgcd_all_i(GEN x, GEN y, GEN p, GEN *pa, GEN *pb)
     780             : {
     781             :   GEN a, b;
     782       16543 :   GEN R = FpX_halfres_i(x, y, p, &a, &b, NULL);
     783       16543 :   if (pa) *pa = a;
     784       16543 :   if (pb) *pb = b;
     785       16543 :   return R;
     786             : }
     787             : 
     788             : /* Return M in GL_2(Fp[X]) such that:
     789             : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
     790             : */
     791             : 
     792             : GEN
     793       16683 : FpX_halfgcd_all(GEN x, GEN y, GEN p, GEN *a, GEN *b)
     794             : {
     795       16683 :   pari_sp av = avma;
     796             :   GEN R, q, r;
     797       16683 :   if (lgefint(p)==3)
     798             :   {
     799         140 :     ulong pp = to_Flx(&x, &y, p);
     800         140 :     R = Flx_halfgcd_all(x, y, pp, a, b);
     801         140 :     R = FlxM_to_ZXM(R);
     802         140 :     if (a) *a = Flx_to_ZX(*a);
     803         140 :     if (b) *b = Flx_to_ZX(*b);
     804         140 :     return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
     805             :   }
     806       16543 :   if (!signe(x))
     807             :   {
     808           0 :     if (a) *a = RgX_copy(y);
     809           0 :     if (b) *b = RgX_copy(x);
     810           0 :     return matJ2_FpXM(varn(x));
     811             :   }
     812       16543 :   if (degpol(y)<degpol(x)) return FpX_halfgcd_all_i(x, y, p, a, b);
     813         389 :   q = FpX_divrem(y,x,p,&r);
     814         389 :   R = FpX_halfgcd_all_i(x, r, p, a, b);
     815         389 :   gcoeff(R,1,1) = FpX_sub(gcoeff(R,1,1), FpX_mul(q, gcoeff(R,1,2), p), p);
     816         389 :   gcoeff(R,2,1) = FpX_sub(gcoeff(R,2,1), FpX_mul(q, gcoeff(R,2,2), p), p);
     817         389 :   return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
     818             : }
     819             : 
     820             : GEN
     821         970 : FpX_halfgcd(GEN x, GEN y, GEN p)
     822         970 : { return FpX_halfgcd_all(x, y, p, NULL, NULL); }
     823             : 
     824             : static GEN
     825       52811 : FpX_gcd_basecase(GEN a, GEN b, GEN p)
     826             : {
     827       52811 :   pari_sp av = avma, av0=avma;
     828      451787 :   while (signe(b))
     829             :   {
     830             :     GEN c;
     831      399264 :     if (gc_needed(av0,2))
     832             :     {
     833           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpX_gcd (d = %ld)",degpol(b));
     834           0 :       gerepileall(av0,2, &a,&b);
     835             :     }
     836      399264 :     av = avma; c = FpX_rem(a,b,p); a=b; b=c;
     837             :   }
     838       52523 :   return gc_const(av, a);
     839             : }
     840             : 
     841             : GEN
     842     1011206 : FpX_gcd(GEN x, GEN y, GEN p)
     843             : {
     844     1011206 :   pari_sp av = avma;
     845     1011206 :   if (lgefint(p)==3)
     846             :   {
     847             :     ulong pp;
     848      957965 :     (void)new_chunk((lg(x) + lg(y)) << 2); /* scratch space */
     849      957965 :     pp = to_Flx(&x, &y, p);
     850      957965 :     x = Flx_gcd(x, y, pp);
     851      957965 :     set_avma(av); return Flx_to_ZX(x);
     852             :   }
     853       53241 :   x = FpX_red(x, p);
     854       53241 :   y = FpX_red(y, p);
     855       53241 :   if (!signe(x)) return gerepileupto(av, y);
     856       53858 :   while (lgpol(y) >= FpX_GCD_LIMIT)
     857             :   {
     858        1047 :     if (lgpol(y)<=(lgpol(x)>>1))
     859             :     {
     860           0 :       GEN r = FpX_rem(x, y, p);
     861           0 :       x = y; y = r;
     862             :     }
     863        1047 :     (void) FpX_halfgcd_all(x, y, p, &x, &y);
     864        1047 :     if (gc_needed(av,2))
     865             :     {
     866           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpX_gcd (y = %ld)",degpol(y));
     867           0 :       gerepileall(av,2,&x,&y);
     868             :     }
     869             :   }
     870       52811 :   return gerepileupto(av, FpX_gcd_basecase(x,y,p));
     871             : }
     872             : 
     873             : /* Return NULL if gcd can be computed else return a factor of p */
     874             : GEN
     875         813 : FpX_gcd_check(GEN x, GEN y, GEN p)
     876             : {
     877         813 :   pari_sp av = avma;
     878             :   GEN a,b,c;
     879             : 
     880         813 :   a = FpX_red(x, p);
     881         813 :   b = FpX_red(y, p);
     882        9034 :   while (signe(b))
     883             :   {
     884             :     GEN g;
     885        8284 :     if (!invmod(leading_coeff(b), p, &g)) return gerepileuptoint(av,g);
     886        8221 :     b = FpX_Fp_mul_to_monic(b, g, p);
     887        8221 :     c = FpX_rem(a, b, p); a = b; b = c;
     888        8221 :     if (gc_needed(av,1))
     889             :     {
     890           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpX_gcd_check (d = %ld)",degpol(b));
     891           0 :       gerepileall(av,2,&a,&b);
     892             :     }
     893             :   }
     894         750 :   return gc_NULL(av);
     895             : }
     896             : 
     897             : static GEN
     898      585609 : FpX_extgcd_basecase(GEN a, GEN b, GEN p, GEN *ptu, GEN *ptv)
     899             : {
     900      585609 :   pari_sp av=avma;
     901      585609 :   GEN v,v1, A = a, B = b;
     902      585609 :   long vx = varn(a);
     903      585609 :   if (!lgpol(b))
     904             :   {
     905           0 :     if (ptu) *ptu = pol_1(vx);
     906           0 :     *ptv = pol_0(vx);
     907           0 :     return RgX_copy(a);
     908             :   }
     909      585609 :   v = pol_0(vx); v1 = pol_1(vx);
     910             :   while (1)
     911     1555599 :   {
     912     2141208 :     GEN r, q = FpX_divrem(a,b,p, &r);
     913     2141208 :     a = b; b = r;
     914     2141208 :     swap(v,v1);
     915     2141208 :     if (!lgpol(b)) break;
     916     1555599 :     v1 = FpX_sub(v1, FpX_mul(v, q, p), p);
     917     1555599 :     if (gc_needed(av,2))
     918             :     {
     919           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpX_extgcd (d = %ld)",degpol(a));
     920           0 :       gerepileall(av,4,&a,&b,&v,&v1);
     921             :     }
     922             :   }
     923      585609 :   if (ptu) *ptu = FpX_div(FpX_sub(a,FpX_mul(B,v,p),p),A,p);
     924      585609 :   *ptv = v;
     925      585609 :   return a;
     926             : }
     927             : 
     928             : static GEN
     929       13757 : FpX_extgcd_halfgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
     930             : {
     931             :   GEN u, v;
     932       13757 :   GEN V = cgetg(expu(lgpol(y))+2,t_VEC);
     933       13757 :   long i, n = 0, vs = varn(x);
     934       28417 :   while (lgpol(y) >= FpX_EXTGCD_LIMIT)
     935             :   {
     936       14660 :     if (lgpol(y)<=(lgpol(x)>>1))
     937             :     {
     938           8 :       GEN r, q = FpX_divrem(x, y, p, &r);
     939           8 :       x = y; y = r;
     940           8 :       gel(V,++n) = mkmat22(pol_0(vs),pol_1(vs),pol_1(vs),FpX_neg(q,p));
     941             :     } else
     942       14652 :       gel(V,++n) = FpX_halfgcd_all(x, y, p, &x, &y);
     943             :   }
     944       13757 :   y = FpX_extgcd_basecase(x, y, p, &u, &v);
     945       14660 :   for (i = n; i>1; i--)
     946             :   {
     947         903 :     GEN R = gel(V,i);
     948         903 :     GEN u1 = FpX_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p);
     949         903 :     GEN v1 = FpX_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p);
     950         903 :     u = u1; v = v1;
     951             :   }
     952             :   {
     953       13757 :     GEN R = gel(V,1);
     954       13757 :     if (ptu)
     955          40 :       *ptu = FpX_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p);
     956       13757 :     *ptv   = FpX_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p);
     957             :   }
     958       13757 :   return y;
     959             : }
     960             : 
     961             : /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
     962             :  * ux + vy = gcd (mod p) */
     963             : GEN
     964     1440738 : FpX_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
     965             : {
     966     1440738 :   pari_sp av = avma;
     967             :   GEN d;
     968     1440738 :   if (lgefint(p)==3)
     969             :   {
     970      855129 :     ulong pp = to_Flx(&x, &y, p);
     971      855131 :     d = Flx_extgcd(x,y, pp, ptu,ptv);
     972      855154 :     d = Flx_to_ZX(d);
     973      855121 :     if (ptu) *ptu = Flx_to_ZX(*ptu);
     974      855122 :     *ptv = Flx_to_ZX(*ptv);
     975             :   }
     976             :   else
     977             :   {
     978      585609 :     x = FpX_red(x, p);
     979      585609 :     y = FpX_red(y, p);
     980      585609 :     if (lgpol(y) >= FpX_EXTGCD_LIMIT)
     981       13757 :       d = FpX_extgcd_halfgcd(x, y, p, ptu, ptv);
     982             :     else
     983      571852 :       d = FpX_extgcd_basecase(x, y, p, ptu, ptv);
     984             :   }
     985     1440722 :   return gc_all(av, ptu?3:2, &d, ptv, ptu);
     986             : }
     987             : 
     988             : static GEN
     989         106 : FpX_halfres(GEN x, GEN y, GEN p, GEN *a, GEN *b, GEN *r)
     990             : {
     991             :   struct FpX_res res;
     992             :   GEN V;
     993             :   long dB;
     994             : 
     995         106 :   res.res  = *r;
     996         106 :   res.lc   = leading_coeff(y);
     997         106 :   res.deg0 = degpol(x);
     998         106 :   res.deg1 = degpol(y);
     999         106 :   res.off = 0;
    1000         106 :   V = FpX_halfres_i(x, y, p, a, b, &res);
    1001         106 :   dB = degpol(*b);
    1002         106 :   if (dB < degpol(y))
    1003         106 :     FpX_halfres_update(res.deg0, res.deg1, dB, p, &res);
    1004         106 :   *r = res.res;
    1005         106 :   return V;
    1006             : }
    1007             : 
    1008             : static GEN
    1009        4002 : FpX_resultant_basecase(GEN a, GEN b, GEN p)
    1010             : {
    1011        4002 :   pari_sp av = avma;
    1012             :   long da,db,dc;
    1013        4002 :   GEN c, lb, res = gen_1;
    1014             : 
    1015        4002 :   if (!signe(a) || !signe(b)) return pol_0(varn(a));
    1016             : 
    1017        4002 :   da = degpol(a);
    1018        4002 :   db = degpol(b);
    1019        4002 :   if (db > da)
    1020             :   {
    1021           0 :     swapspec(a,b, da,db);
    1022           0 :     if (both_odd(da,db)) res = subii(p, res);
    1023             :   }
    1024        4002 :   if (!da) return gc_const(av, gen_1); /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
    1025       11324 :   while (db)
    1026             :   {
    1027        7322 :     lb = gel(b,db+2);
    1028        7322 :     c = FpX_rem(a,b, p);
    1029        7322 :     a = b; b = c; dc = degpol(c);
    1030        7322 :     if (dc < 0) return gc_const(av, gen_0);
    1031             : 
    1032        7322 :     if (both_odd(da,db)) res = subii(p, res);
    1033        7322 :     if (!equali1(lb)) res = Fp_mul(res, Fp_powu(lb, da - dc, p), p);
    1034        7322 :     if (gc_needed(av,2))
    1035             :     {
    1036           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpX_resultant (da = %ld)",da);
    1037           0 :       gerepileall(av,3, &a,&b,&res);
    1038             :     }
    1039        7322 :     da = db; /* = degpol(a) */
    1040        7322 :     db = dc; /* = degpol(b) */
    1041             :   }
    1042        4002 :   return gerepileuptoint(av, Fp_mul(res, Fp_powu(gel(b,2), da, p), p));
    1043             : }
    1044             : 
    1045             : GEN
    1046      416662 : FpX_resultant(GEN x, GEN y, GEN p)
    1047             : {
    1048      416662 :   pari_sp av = avma;
    1049             :   long dx, dy;
    1050      416662 :   GEN res = gen_1;
    1051      416662 :   if (!signe(x) || !signe(y)) return gen_0;
    1052      416662 :   if (lgefint(p) == 3)
    1053             :   {
    1054      412660 :     pari_sp av = avma;
    1055      412660 :     ulong pp = to_Flx(&x, &y, p);
    1056      412660 :     ulong res = Flx_resultant(x, y, pp);
    1057      412660 :     return gc_utoi(av, res);
    1058             :   }
    1059        4002 :   dx = degpol(x); dy = degpol(y);
    1060        4002 :   if (dx < dy)
    1061             :   {
    1062           0 :     swap(x,y);
    1063           0 :     if (both_odd(dx, dy))
    1064           0 :       res = Fp_neg(res, p);
    1065             :   }
    1066        4009 :   while (lgpol(y) >= FpX_GCD_LIMIT)
    1067             :   {
    1068           7 :     if (lgpol(y)<=(lgpol(x)>>1))
    1069             :     {
    1070           0 :       GEN r = FpX_rem(x, y, p);
    1071           0 :       long dx = degpol(x), dy = degpol(y), dr = degpol(r);
    1072           0 :       GEN ly = gel(y,dy+2);
    1073           0 :       if (!equali1(ly)) res = Fp_mul(res, Fp_powu(ly, dx - dr, p), p);
    1074           0 :       if (both_odd(dx, dy))
    1075           0 :         res = Fp_neg(res, p);
    1076           0 :       x = y; y = r;
    1077             :     }
    1078           7 :     (void) FpX_halfres(x, y, p, &x, &y, &res);
    1079           7 :     if (gc_needed(av,2))
    1080             :     {
    1081           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpX_res (y = %ld)",degpol(y));
    1082           0 :       gerepileall(av,3,&x,&y,&res);
    1083             :     }
    1084             :   }
    1085        4002 :   return gerepileuptoint(av, Fp_mul(res, FpX_resultant_basecase(x, y, p), p));
    1086             : }
    1087             : 
    1088             : /* If resultant is 0, *ptU and *ptV are not set */
    1089             : static GEN
    1090          24 : FpX_extresultant_basecase(GEN a, GEN b, GEN p, GEN *ptU, GEN *ptV)
    1091             : {
    1092          24 :   pari_sp av = avma;
    1093          24 :   GEN z,q,u,v, x = a, y = b;
    1094          24 :   GEN lb, res = gen_1;
    1095             :   long dx, dy, dz;
    1096          24 :   long vs = varn(a);
    1097             : 
    1098          24 :   u = pol_0(vs);
    1099          24 :   v = pol_1(vs); /* v = 1 */
    1100          24 :   dx = degpol(x);
    1101          24 :   dy = degpol(y);
    1102         281 :   while (dy)
    1103             :   { /* b u = x (a), b v = y (a) */
    1104         257 :     lb = gel(y,dy+2);
    1105         257 :     q = FpX_divrem(x,y, p, &z);
    1106         257 :     x = y; y = z; /* (x,y) = (y, x - q y) */
    1107         257 :     dz = degpol(z); if (dz < 0) return gc_const(av,gen_0);
    1108         257 :     z = FpX_sub(u, FpX_mul(q,v, p), p);
    1109         257 :     u = v; v = z; /* (u,v) = (v, u - q v) */
    1110             : 
    1111         257 :     if (both_odd(dx,dy)) res = Fp_neg(res, p);
    1112         257 :     if (!equali1(lb)) res = Fp_mul(res, Fp_powu(lb, dx-dz, p), p);
    1113         257 :     dx = dy; /* = degpol(x) */
    1114         257 :     dy = dz; /* = degpol(y) */
    1115             :   }
    1116          24 :   res = Fp_mul(res, Fp_powu(gel(y,2), dx, p), p);
    1117          24 :   lb = Fp_mul(res, Fp_inv(gel(y,2),p), p);
    1118          24 :   v = FpX_Fp_mul(v, lb, p);
    1119          24 :   u = Fp_FpX_sub(res, FpX_mul(b,v,p), p);
    1120          24 :   u = FpX_div(u,a,p); /* = (res - b v) / a */
    1121          24 :   *ptU = u;
    1122          24 :   *ptV = v;
    1123          24 :   return res;
    1124             : }
    1125             : 
    1126             : GEN
    1127          77 : FpX_extresultant(GEN x, GEN y, GEN p, GEN *ptU, GEN *ptV)
    1128             : {
    1129          77 :   pari_sp av=avma;
    1130             :   GEN u, v, R;
    1131          77 :   GEN res = gen_1, res1;
    1132          77 :   long dx = degpol(x), dy = degpol(y);
    1133          77 :   if (lgefint(p) == 3)
    1134             :   {
    1135          53 :     pari_sp av = avma;
    1136          53 :     ulong pp = to_Flx(&x, &y, p);
    1137          53 :     ulong resp = Flx_extresultant(x, y, pp, &u, &v);
    1138          53 :     if (!resp) return gc_const(av, gen_0);
    1139          53 :     res = utoi(resp);
    1140          53 :     *ptU = Flx_to_ZX(u); *ptV = Flx_to_ZX(v);
    1141          53 :     return gc_all(av, 3, &res, ptU, ptV);
    1142             :   }
    1143          24 :   if (dy > dx)
    1144             :   {
    1145           8 :     swap(x,y); lswap(dx,dy);
    1146           8 :     if (both_odd(dx,dy)) res = Fp_neg(res,p);
    1147           8 :     R = matJ2_FpXM(x[1]);
    1148          16 :   } else R = matid2_FpXM(x[1]);
    1149          24 :   if (dy < 0) return gen_0;
    1150         123 :   while (lgpol(y) >= FpX_EXTGCD_LIMIT)
    1151             :   {
    1152             :     GEN M;
    1153          99 :     if (lgpol(y)<=(lgpol(x)>>1))
    1154             :     {
    1155           8 :       GEN r, q = FpX_divrem(x, y, p, &r);
    1156           8 :       long dx = degpol(x), dy = degpol(y), dr = degpol(r);
    1157           8 :       GEN ly = gel(y,dy+2);
    1158           8 :       if (!equali1(ly)) res = Fp_mul(res, Fp_powu(ly, dx - dr, p), p);
    1159           8 :       if (both_odd(dx, dy))
    1160           0 :         res = Fp_neg(res, p);
    1161           8 :       x = y; y = r;
    1162           8 :       R = FpX_FpXM_qmul(q, R, p);
    1163             :     }
    1164          99 :     M = FpX_halfres(x, y, p, &x, &y, &res);
    1165          99 :     if (!signe(res)) return gc_const(av, gen_0);
    1166          99 :     R = FpXM_mul2(M, R, p);
    1167          99 :     gerepileall(av,4,&x,&y,&R,&res);
    1168             :   }
    1169          24 :   res1 = FpX_extresultant_basecase(x,y,p,&u,&v);
    1170          24 :   if (!signe(res1)) return gc_const(av, gen_0);
    1171          24 :   *ptU = FpX_Fp_mul(FpX_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p), res, p);
    1172          24 :   *ptV = FpX_Fp_mul(FpX_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p), res, p);
    1173          24 :   res = Fp_mul(res1,res,p);
    1174          24 :   return gc_all(av, 3, &res, ptU, ptV);
    1175             : }
    1176             : 
    1177             : GEN
    1178      178891 : FpX_rescale(GEN P, GEN h, GEN p)
    1179             : {
    1180      178891 :   long i, l = lg(P);
    1181      178891 :   GEN Q = cgetg(l,t_POL), hi = h;
    1182      178888 :   gel(Q,l-1) = gel(P,l-1);
    1183      365712 :   for (i=l-2; i>=2; i--)
    1184             :   {
    1185      365710 :     gel(Q,i) = Fp_mul(gel(P,i), hi, p);
    1186      365713 :     if (i == 2) break;
    1187      186821 :     hi = Fp_mul(hi,h, p);
    1188             :   }
    1189      178894 :   Q[1] = P[1]; return Q;
    1190             : }
    1191             : 
    1192             : GEN
    1193     1633255 : FpX_deriv(GEN x, GEN p) { return FpX_red(ZX_deriv(x), p); }
    1194             : 
    1195             : /* Compute intformal(x^n*S)/x^(n+1) */
    1196             : static GEN
    1197       56163 : FpX_integXn(GEN x, long n, GEN p)
    1198             : {
    1199       56163 :   long i, lx = lg(x);
    1200             :   GEN y;
    1201       56163 :   if (lx == 2) return ZX_copy(x);
    1202       54898 :   y = cgetg(lx, t_POL); y[1] = x[1];
    1203      199691 :   for (i=2; i<lx; i++)
    1204             :   {
    1205      144793 :     GEN xi = gel(x,i);
    1206      144793 :     if (!signe(xi))
    1207           0 :       gel(y,i) = gen_0;
    1208             :     else
    1209             :     {
    1210      144793 :       ulong j = n+i-1, d = ugcdiu(xi, j);
    1211      144793 :       if (d==1)
    1212       93207 :         gel(y,i) = Fp_divu(xi, j, p);
    1213             :       else
    1214       51586 :         gel(y,i) = Fp_divu(diviuexact(xi, d), j/d, p);
    1215             :     }
    1216             :   }
    1217       54898 :   return ZX_renormalize(y, lx);;
    1218             : }
    1219             : 
    1220             : GEN
    1221           0 : FpX_integ(GEN x, GEN p)
    1222             : {
    1223           0 :   long i, lx = lg(x);
    1224             :   GEN y;
    1225           0 :   if (lx == 2) return ZX_copy(x);
    1226           0 :   y = cgetg(lx+1, t_POL); y[1] = x[1];
    1227           0 :   gel(y,2) = gen_0;
    1228           0 :   for (i=3; i<=lx; i++)
    1229           0 :     gel(y,i) = signe(gel(x,i-1))? Fp_divu(gel(x,i-1), i-2, p): gen_0;
    1230           0 :   return ZX_renormalize(y, lx+1);;
    1231             : }
    1232             : 
    1233             : INLINE GEN
    1234      535526 : FpXn_recip(GEN P, long n)
    1235      535526 : { return RgXn_recip_shallow(P, n); }
    1236             : 
    1237             : GEN
    1238      524577 : FpX_Newton(GEN P, long n, GEN p)
    1239             : {
    1240      524577 :   pari_sp av = avma;
    1241      524577 :   GEN dP = FpX_deriv(P, p);
    1242      524569 :   GEN Q = FpXn_recip(FpX_div(FpX_shift(dP,n), P, p), n);
    1243      524585 :   return gerepilecopy(av, Q);
    1244             : }
    1245             : 
    1246             : GEN
    1247       11446 : FpX_fromNewton(GEN P, GEN p)
    1248             : {
    1249       11446 :   pari_sp av = avma;
    1250       11446 :   if (lgefint(p)==3)
    1251             :   {
    1252         497 :     ulong pp = p[2];
    1253         497 :     GEN Q = Flx_fromNewton(ZX_to_Flx(P, pp), pp);
    1254         497 :     return gerepileupto(av, Flx_to_ZX(Q));
    1255             :   } else
    1256             :   {
    1257       10949 :     long n = itos(modii(constant_coeff(P), p))+1;
    1258       10949 :     GEN z = FpX_neg(FpX_shift(P,-1),p);
    1259       10949 :     GEN Q = FpXn_recip(FpXn_expint(z, n, p), n);
    1260       10949 :     return gerepilecopy(av, Q);
    1261             :   }
    1262             : }
    1263             : 
    1264             : GEN
    1265         382 : FpX_invLaplace(GEN x, GEN p)
    1266             : {
    1267         382 :   pari_sp av = avma;
    1268         382 :   long i, d = degpol(x);
    1269             :   GEN t, y;
    1270         382 :   if (d <= 1) return gcopy(x);
    1271         382 :   t = Fp_inv(factorial_Fp(d, p), p);
    1272         382 :   y = cgetg(d+3, t_POL);
    1273         382 :   y[1] = x[1];
    1274       13424 :   for (i=d; i>=2; i--)
    1275             :   {
    1276       13042 :     gel(y,i+2) = Fp_mul(gel(x,i+2), t, p);
    1277       13042 :     t = Fp_mulu(t, i, p);
    1278             :   }
    1279         382 :   gel(y,3) = gel(x,3);
    1280         382 :   gel(y,2) = gel(x,2);
    1281         382 :   return gerepilecopy(av, y);
    1282             : }
    1283             : 
    1284             : GEN
    1285         688 : FpX_Laplace(GEN x, GEN p)
    1286             : {
    1287         688 :   pari_sp av = avma;
    1288         688 :   long i, d = degpol(x);
    1289         688 :   GEN t = gen_1;
    1290             :   GEN y;
    1291         688 :   if (d <= 1) return gcopy(x);
    1292         688 :   y = cgetg(d+3, t_POL);
    1293         688 :   y[1] = x[1];
    1294         688 :   gel(y,2) = gel(x,2);
    1295         688 :   gel(y,3) = gel(x,3);
    1296       35097 :   for (i=2; i<=d; i++)
    1297             :   {
    1298       34409 :     t = Fp_mulu(t, i, p);
    1299       34409 :     gel(y,i+2) = Fp_mul(gel(x,i+2), t, p);
    1300             :   }
    1301         688 :   return gerepilecopy(av, y);
    1302             : }
    1303             : 
    1304             : int
    1305       40968 : FpX_is_squarefree(GEN f, GEN p)
    1306             : {
    1307       40968 :   pari_sp av = avma;
    1308       40968 :   GEN z = FpX_gcd(f,FpX_deriv(f,p),p);
    1309       40969 :   set_avma(av);
    1310       40969 :   return degpol(z)==0;
    1311             : }
    1312             : 
    1313             : GEN
    1314      256175 : random_FpX(long d1, long v, GEN p)
    1315             : {
    1316      256175 :   long i, d = d1+2;
    1317      256175 :   GEN y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
    1318      934749 :   for (i=2; i<d; i++) gel(y,i) = randomi(p);
    1319      256175 :   return FpX_renormalize(y,d);
    1320             : }
    1321             : 
    1322             : GEN
    1323       61132 : FpX_dotproduct(GEN x, GEN y, GEN p)
    1324             : {
    1325       61132 :   long i, l = minss(lg(x), lg(y));
    1326             :   pari_sp av;
    1327             :   GEN c;
    1328       61132 :   if (l == 2) return gen_0;
    1329       61055 :   av = avma; c = mulii(gel(x,2),gel(y,2));
    1330     5010203 :   for (i=3; i<l; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
    1331       61055 :   return gerepileuptoint(av, modii(c,p));
    1332             : }
    1333             : 
    1334             : /* Evaluation in Fp
    1335             :  * x a ZX and y an Fp, return x(y) mod p
    1336             :  *
    1337             :  * If p is very large (several longs) and x has small coefficients(<<p),
    1338             :  * then Brent & Kung algorithm is faster. */
    1339             : GEN
    1340      964318 : FpX_eval(GEN x,GEN y,GEN p)
    1341             : {
    1342             :   pari_sp av;
    1343             :   GEN p1,r,res;
    1344      964318 :   long j, i=lg(x)-1;
    1345      964318 :   if (i<=2 || !signe(y))
    1346      181867 :     return (i==1)? gen_0: modii(gel(x,2),p);
    1347      782451 :   res=cgeti(lgefint(p));
    1348      782454 :   av=avma; p1=gel(x,i);
    1349             :   /* specific attention to sparse polynomials (see poleval)*/
    1350             :   /*You've guessed it! It's a copy-paste(tm)*/
    1351     3404151 :   for (i--; i>=2; i=j-1)
    1352             :   {
    1353     3700489 :     for (j=i; !signe(gel(x,j)); j--)
    1354     1078780 :       if (j==2)
    1355             :       {
    1356      162424 :         if (i!=j) y = Fp_powu(y,i-j+1,p);
    1357      162424 :         p1=mulii(p1,y);
    1358      162411 :         goto fppoleval;/*sorry break(2) no implemented*/
    1359             :       }
    1360     2621709 :     r = (i==j)? y: Fp_powu(y,i-j+1,p);
    1361     2621705 :     p1 = Fp_addmul(gel(x,j), p1, r, p);
    1362     2621697 :     if ((i & 7) == 0) { affii(p1, res); p1 = res; set_avma(av); }
    1363             :   }
    1364      620018 :  fppoleval:
    1365      782429 :   modiiz(p1,p,res); return gc_const(av, res);
    1366             : }
    1367             : 
    1368             : /* Tz=Tx*Ty where Tx and Ty coprime
    1369             :  * return lift(chinese(Mod(x*Mod(1,p),Tx*Mod(1,p)),Mod(y*Mod(1,p),Ty*Mod(1,p))))
    1370             :  * if Tz is NULL it is computed
    1371             :  * As we do not return it, and the caller will frequently need it,
    1372             :  * it must compute it and pass it.
    1373             :  */
    1374             : GEN
    1375           0 : FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,GEN p)
    1376             : {
    1377           0 :   pari_sp av = avma;
    1378             :   GEN ax,p1;
    1379           0 :   ax = FpX_mul(FpXQ_inv(Tx,Ty,p), Tx,p);
    1380           0 :   p1 = FpX_mul(ax, FpX_sub(y,x,p),p);
    1381           0 :   p1 = FpX_add(x,p1,p);
    1382           0 :   if (!Tz) Tz=FpX_mul(Tx,Ty,p);
    1383           0 :   p1 = FpX_rem(p1,Tz,p);
    1384           0 :   return gerepileupto(av,p1);
    1385             : }
    1386             : 
    1387             : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
    1388             : GEN
    1389          42 : FpX_disc(GEN P, GEN p)
    1390             : {
    1391          42 :   pari_sp av = avma;
    1392          42 :   GEN L, dP = FpX_deriv(P,p), D = FpX_resultant(P, dP, p);
    1393             :   long dd;
    1394          42 :   if (!signe(D)) return gen_0;
    1395          35 :   dd = degpol(P) - 2 - degpol(dP); /* >= -1; > -1 iff p | deg(P) */
    1396          35 :   L = leading_coeff(P);
    1397          35 :   if (dd && !equali1(L))
    1398           7 :     D = (dd == -1)? Fp_div(D,L,p): Fp_mul(D, Fp_powu(L, dd, p), p);
    1399          35 :   if (degpol(P) & 2) D = Fp_neg(D ,p);
    1400          35 :   return gerepileuptoint(av, D);
    1401             : }
    1402             : 
    1403             : GEN
    1404       93400 : FpV_roots_to_pol(GEN V, GEN p, long v)
    1405             : {
    1406       93400 :   pari_sp ltop=avma;
    1407             :   long i;
    1408       93400 :   GEN g=cgetg(lg(V),t_VEC);
    1409      404119 :   for(i=1;i<lg(V);i++)
    1410      310719 :     gel(g,i) = deg1pol_shallow(gen_1,modii(negi(gel(V,i)),p),v);
    1411       93400 :   return gerepileupto(ltop,FpXV_prod(g,p));
    1412             : }
    1413             : 
    1414             : /* invert all elements of x mod p using Montgomery's multi-inverse trick.
    1415             :  * Not stack-clean. */
    1416             : GEN
    1417       34108 : FpV_inv(GEN x, GEN p)
    1418             : {
    1419       34108 :   long i, lx = lg(x);
    1420       34108 :   GEN u, y = cgetg(lx, t_VEC);
    1421             : 
    1422       34108 :   gel(y,1) = gel(x,1);
    1423      471634 :   for (i=2; i<lx; i++) gel(y,i) = Fp_mul(gel(y,i-1), gel(x,i), p);
    1424             : 
    1425       34108 :   u = Fp_inv(gel(y,--i), p);
    1426      471633 :   for ( ; i > 1; i--)
    1427             :   {
    1428      437526 :     gel(y,i) = Fp_mul(u, gel(y,i-1), p);
    1429      437526 :     u = Fp_mul(u, gel(x,i), p); /* u = 1 / (x[1] ... x[i-1]) */
    1430             :   }
    1431       34107 :   gel(y,1) = u; return y;
    1432             : }
    1433             : GEN
    1434           0 : FqV_inv(GEN x, GEN T, GEN p)
    1435             : {
    1436           0 :   long i, lx = lg(x);
    1437           0 :   GEN u, y = cgetg(lx, t_VEC);
    1438             : 
    1439           0 :   gel(y,1) = gel(x,1);
    1440           0 :   for (i=2; i<lx; i++) gel(y,i) = Fq_mul(gel(y,i-1), gel(x,i), T,p);
    1441             : 
    1442           0 :   u = Fq_inv(gel(y,--i), T,p);
    1443           0 :   for ( ; i > 1; i--)
    1444             :   {
    1445           0 :     gel(y,i) = Fq_mul(u, gel(y,i-1), T,p);
    1446           0 :     u = Fq_mul(u, gel(x,i), T,p); /* u = 1 / (x[1] ... x[i-1]) */
    1447             :   }
    1448           0 :   gel(y,1) = u; return y;
    1449             : }
    1450             : 
    1451             : /***********************************************************************/
    1452             : /**                                                                   **/
    1453             : /**                      Barrett reduction                            **/
    1454             : /**                                                                   **/
    1455             : /***********************************************************************/
    1456             : 
    1457             : static GEN
    1458        3975 : FpX_invBarrett_basecase(GEN T, GEN p)
    1459             : {
    1460        3975 :   long i, l=lg(T)-1, lr = l-1, k;
    1461        3975 :   GEN r=cgetg(lr, t_POL); r[1]=T[1];
    1462        3975 :   gel(r,2) = gen_1;
    1463      226419 :   for (i=3; i<lr; i++)
    1464             :   {
    1465      222444 :     pari_sp av = avma;
    1466      222444 :     GEN u = gel(T,l-i+2);
    1467     7066499 :     for (k=3; k<i; k++)
    1468     6844055 :       u = addii(u, mulii(gel(T,l-i+k), gel(r,k)));
    1469      222444 :     gel(r,i) = gerepileupto(av, modii(negi(u), p));
    1470             :   }
    1471        3975 :   return FpX_renormalize(r,lr);
    1472             : }
    1473             : 
    1474             : /* Return new lgpol */
    1475             : static long
    1476     1662841 : ZX_lgrenormalizespec(GEN x, long lx)
    1477             : {
    1478             :   long i;
    1479     2028414 :   for (i = lx-1; i>=0; i--)
    1480     2028415 :     if (signe(gel(x,i))) break;
    1481     1662841 :   return i+1;
    1482             : }
    1483             : 
    1484             : INLINE GEN
    1485     1638348 : FpX_recipspec(GEN x, long l, long n)
    1486             : {
    1487     1638348 :   return RgX_recipspec_shallow(x, l, n);
    1488             : }
    1489             : 
    1490             : static GEN
    1491        1486 : FpX_invBarrett_Newton(GEN T, GEN p)
    1492             : {
    1493        1486 :   pari_sp av = avma;
    1494        1486 :   long nold, lx, lz, lq, l = degpol(T), i, lQ;
    1495        1486 :   GEN q, y, z, x = cgetg(l+2, t_POL) + 2;
    1496        1486 :   ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
    1497      595486 :   for (i=0;i<l;i++) gel(x,i) = gen_0;
    1498        1486 :   q = FpX_recipspec(T+2,l+1,l+1); lQ = lgpol(q); q+=2;
    1499             :   /* We work on _spec_ FpX's, all the l[xzq] below are lgpol's */
    1500             : 
    1501             :   /* initialize */
    1502        1486 :   gel(x,0) = Fp_inv(gel(q,0), p);
    1503        1486 :   if (lQ>1) gel(q,1) = Fp_red(gel(q,1), p);
    1504        1486 :   if (lQ>1 && signe(gel(q,1)))
    1505        1099 :   {
    1506        1099 :     GEN u = gel(q, 1);
    1507        1099 :     if (!equali1(gel(x,0))) u = Fp_mul(u, Fp_sqr(gel(x,0), p), p);
    1508        1099 :     gel(x,1) = Fp_neg(u, p); lx = 2;
    1509             :   }
    1510             :   else
    1511         387 :     lx = 1;
    1512        1486 :   nold = 1;
    1513       13401 :   for (; mask > 1; )
    1514             :   { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
    1515       11932 :     long i, lnew, nnew = nold << 1;
    1516             : 
    1517       11932 :     if (mask & 1) nnew--;
    1518       11932 :     mask >>= 1;
    1519             : 
    1520       11932 :     lnew = nnew + 1;
    1521       11932 :     lq = ZX_lgrenormalizespec(q, minss(lQ,lnew));
    1522       11931 :     z = FpX_mulspec(x, q, p, lx, lq); /* FIXME: high product */
    1523       11933 :     lz = lgpol(z); if (lz > lnew) lz = lnew;
    1524       11933 :     z += 2;
    1525             :     /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
    1526       83869 :     for (i = nold; i < lz; i++) if (signe(gel(z,i))) break;
    1527       11933 :     nold = nnew;
    1528       11933 :     if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
    1529             : 
    1530             :     /* z + i represents (x*q - 1) / t^i */
    1531        9390 :     lz = ZX_lgrenormalizespec (z+i, lz-i);
    1532        9390 :     z = FpX_mulspec(x, z+i, p, lx, lz); /* FIXME: low product */
    1533        9390 :     lz = lgpol(z); z += 2;
    1534        9390 :     if (lz > lnew-i) lz = ZX_lgrenormalizespec(z, lnew-i);
    1535             : 
    1536        9390 :     lx = lz+ i;
    1537        9390 :     y  = x + i; /* x -= z * t^i, in place */
    1538      428240 :     for (i = 0; i < lz; i++) gel(y,i) = Fp_neg(gel(z,i), p);
    1539             :   }
    1540        1469 :   x -= 2; setlg(x, lx + 2); x[1] = T[1];
    1541        1486 :   return gerepilecopy(av, x);
    1542             : }
    1543             : 
    1544             : /* 1/polrecip(T)+O(x^(deg(T)-1)) */
    1545             : GEN
    1546        5514 : FpX_invBarrett(GEN T, GEN p)
    1547             : {
    1548        5514 :   pari_sp ltop = avma;
    1549        5514 :   long l = lg(T);
    1550             :   GEN r;
    1551        5514 :   if (l<5) return pol_0(varn(T));
    1552        5461 :   if (l<=FpX_INVBARRETT_LIMIT)
    1553             :   {
    1554        3975 :     GEN c = gel(T,l-1), ci=gen_1;
    1555        3975 :     if (!equali1(c))
    1556             :     {
    1557          14 :       ci = Fp_inv(c, p);
    1558          14 :       T = FpX_Fp_mul(T, ci, p);
    1559          14 :       r = FpX_invBarrett_basecase(T, p);
    1560          14 :       r = FpX_Fp_mul(r, ci, p);
    1561             :     } else
    1562        3961 :       r = FpX_invBarrett_basecase(T, p);
    1563             :   }
    1564             :   else
    1565        1486 :     r = FpX_invBarrett_Newton(T, p);
    1566        5461 :   return gerepileupto(ltop, r);
    1567             : }
    1568             : 
    1569             : GEN
    1570     1067268 : FpX_get_red(GEN T, GEN p)
    1571             : {
    1572     1067268 :   if (typ(T)==t_POL && lg(T)>FpX_BARRETT_LIMIT)
    1573        4677 :     retmkvec2(FpX_invBarrett(T,p),T);
    1574     1062591 :   return T;
    1575             : }
    1576             : 
    1577             : /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
    1578             :  * and mg is the Barrett inverse of T. */
    1579             : static GEN
    1580      816265 : FpX_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, GEN p, GEN *pr)
    1581             : {
    1582             :   GEN q, r;
    1583      816265 :   long lt = degpol(T); /*We discard the leading term*/
    1584             :   long ld, lm, lT, lmg;
    1585      816265 :   ld = l-lt;
    1586      816265 :   lm = minss(ld, lgpol(mg));
    1587      816265 :   lT  = ZX_lgrenormalizespec(T+2,lt);
    1588      816265 :   lmg = ZX_lgrenormalizespec(mg+2,lm);
    1589      816265 :   q = FpX_recipspec(x+lt,ld,ld);              /* q = rec(x)     lq<=ld*/
    1590      816266 :   q = FpX_mulspec(q+2,mg+2,p,lgpol(q),lmg);    /* q = rec(x) * mg lq<=ld+lm*/
    1591      816265 :   q = FpX_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lq<=ld*/
    1592      816265 :   if (!pr) return q;
    1593      816265 :   r = FpX_mulspec(q+2,T+2,p,lgpol(q),lT);      /* r = q*pol        lr<=ld+lt*/
    1594      816263 :   r = FpX_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - r   lr<=lt */
    1595      816266 :   if (pr == ONLY_REM) return r;
    1596        1330 :   *pr = r; return q;
    1597             : }
    1598             : 
    1599             : static GEN
    1600      815520 : FpX_divrem_Barrett(GEN x, GEN mg, GEN T, GEN p, GEN *pr)
    1601             : {
    1602      815520 :   GEN q = NULL, r = FpX_red(x, p);
    1603      815519 :   long l = lgpol(r), lt = degpol(T), lm = 2*lt-1, v = varn(T);
    1604             :   long i;
    1605      815519 :   if (l <= lt)
    1606             :   {
    1607           0 :     if (pr == ONLY_REM) return r;
    1608           0 :     if (pr == ONLY_DIVIDES) return signe(r)? NULL: pol_0(v);
    1609           0 :     if (pr) *pr = r;
    1610           0 :     return pol_0(v);
    1611             :   }
    1612      815519 :   if (lt <= 1)
    1613          53 :     return FpX_divrem_basecase(r,T,p,pr);
    1614      815466 :   if (pr != ONLY_REM && l>lm)
    1615             :   {
    1616         497 :     q = cgetg(l-lt+2, t_POL); q[1] = T[1];
    1617      905007 :     for (i=0;i<l-lt;i++) gel(q+2,i) = gen_0;
    1618             :   }
    1619      816266 :   while (l>lm)
    1620             :   {
    1621         800 :     GEN zr, zq = FpX_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,&zr);
    1622         800 :     long lz = lgpol(zr);
    1623         800 :     if (pr != ONLY_REM)
    1624             :     {
    1625         626 :       long lq = lgpol(zq);
    1626      464768 :       for(i=0; i<lq; i++) gel(q+2+l-lm,i) = gel(zq,2+i);
    1627             :     }
    1628      475648 :     for(i=0; i<lz; i++) gel(r+2+l-lm,i) = gel(zr,2+i);
    1629         800 :     l = l-lm+lz;
    1630             :   }
    1631      815466 :   if (pr == ONLY_REM)
    1632             :   {
    1633      814935 :     if (l > lt)
    1634      814935 :       r = FpX_divrem_Barrettspec(r+2, l, mg, T, p, ONLY_REM);
    1635             :     else
    1636           0 :       r = FpX_renormalize(r, l+2);
    1637      814936 :     setvarn(r, v); return r;
    1638             :   }
    1639         531 :   if (l > lt)
    1640             :   {
    1641         530 :     GEN zq = FpX_divrem_Barrettspec(r+2,l,mg,T,p, pr? &r: NULL);
    1642         530 :     if (!q) q = zq;
    1643             :     else
    1644             :     {
    1645         496 :       long lq = lgpol(zq);
    1646      440483 :       for(i=0; i<lq; i++) gel(q+2,i) = gel(zq,2+i);
    1647             :     }
    1648             :   }
    1649           1 :   else if (pr)
    1650           1 :     r = FpX_renormalize(r, l+2);
    1651         531 :   setvarn(q, v); q = FpX_renormalize(q, lg(q));
    1652         531 :   if (pr == ONLY_DIVIDES) return signe(r)? NULL: q;
    1653         531 :   if (pr) { setvarn(r, v); *pr = r; }
    1654         531 :   return q;
    1655             : }
    1656             : 
    1657             : GEN
    1658    14311118 : FpX_divrem(GEN x, GEN T, GEN p, GEN *pr)
    1659             : {
    1660             :   GEN B, y;
    1661             :   long dy, dx, d;
    1662    14311118 :   if (pr==ONLY_REM) return FpX_rem(x, T, p);
    1663    14311118 :   y = get_FpX_red(T, &B);
    1664    14311108 :   dy = degpol(y); dx = degpol(x); d = dx-dy;
    1665    14311059 :   if (!B && d+3 < FpX_DIVREM_BARRETT_LIMIT)
    1666    14309177 :     return FpX_divrem_basecase(x,y,p,pr);
    1667        1882 :   else if (lgefint(p)==3)
    1668             :   {
    1669        1318 :     pari_sp av = avma;
    1670        1318 :     ulong pp = to_Flxq(&x, &T, p);
    1671        1318 :     GEN z = Flx_divrem(x, T, pp, pr);
    1672        1318 :     if (!z) return gc_NULL(av);
    1673        1318 :     if (!pr || pr == ONLY_DIVIDES)
    1674          59 :       return Flx_to_ZX_inplace(gerepileuptoleaf(av, z));
    1675        1259 :     z = Flx_to_ZX(z);
    1676        1259 :     *pr = Flx_to_ZX(*pr);
    1677        1259 :     return gc_all(av, 2, &z, pr);
    1678             :   } else
    1679             :   {
    1680         564 :     pari_sp av = avma;
    1681         564 :     GEN mg = B? B: FpX_invBarrett(y, p);
    1682         564 :     GEN z = FpX_divrem_Barrett(x,mg,y,p,pr);
    1683         564 :     if (!z) return gc_NULL(av);
    1684         564 :     if (!pr || pr==ONLY_DIVIDES) return gerepilecopy(av, z);
    1685         564 :     return gc_all(av, 2, &z, pr);
    1686             :   }
    1687             : }
    1688             : 
    1689             : GEN
    1690    71201191 : FpX_rem(GEN x, GEN T, GEN p)
    1691             : {
    1692    71201191 :   GEN B, y = get_FpX_red(T, &B);
    1693    71224449 :   long dy = degpol(y), dx = degpol(x), d = dx-dy;
    1694    71244598 :   if (d < 0) return FpX_red(x,p);
    1695    52318759 :   if (!B && d+3 < FpX_REM_BARRETT_LIMIT)
    1696    51464106 :     return FpX_divrem_basecase(x,y,p,ONLY_REM);
    1697      854653 :   else if (lgefint(p)==3)
    1698             :   {
    1699       39697 :     pari_sp av = avma;
    1700       39697 :     ulong pp = to_Flxq(&x, &T, p);
    1701       39699 :     return Flx_to_ZX_inplace(gerepileuptoleaf(av, Flx_rem(x, T, pp)));
    1702             :   } else
    1703             :   {
    1704      814956 :     pari_sp av = avma;
    1705      814956 :     GEN mg = B? B: FpX_invBarrett(y, p);
    1706      814956 :     return gerepileupto(av, FpX_divrem_Barrett(x, mg, y, p, ONLY_REM));
    1707             :   }
    1708             : }
    1709             : 
    1710             : static GEN
    1711       32269 : FpXV_producttree_dbl(GEN t, long n, GEN p)
    1712             : {
    1713       32269 :   long i, j, k, m = n==1 ? 1: expu(n-1)+1;
    1714       32269 :   GEN T = cgetg(m+1, t_VEC);
    1715       32269 :   gel(T,1) = t;
    1716       63649 :   for (i=2; i<=m; i++)
    1717             :   {
    1718       31380 :     GEN u = gel(T, i-1);
    1719       31380 :     long n = lg(u)-1;
    1720       31380 :     GEN t = cgetg(((n+1)>>1)+1, t_VEC);
    1721      102394 :     for (j=1, k=1; k<n; j++, k+=2)
    1722       71014 :       gel(t, j) = FpX_mul(gel(u, k), gel(u, k+1), p);
    1723       31380 :     gel(T, i) = t;
    1724             :   }
    1725       32269 :   return T;
    1726             : }
    1727             : 
    1728             : static GEN
    1729       31681 : FpV_producttree(GEN xa, GEN s, GEN p, long vs)
    1730             : {
    1731       31681 :   long n = lg(xa)-1;
    1732       31681 :   long j, k, ls = lg(s);
    1733       31681 :   GEN t = cgetg(ls, t_VEC);
    1734      132108 :   for (j=1, k=1; j<ls; k+=s[j++])
    1735      100427 :     gel(t, j) = s[j] == 1 ?
    1736      100427 :              deg1pol_shallow(gen_1, Fp_neg(gel(xa,k), p), vs):
    1737       61829 :              deg2pol_shallow(gen_1,
    1738       61829 :                Fp_neg(Fp_add(gel(xa,k), gel(xa,k+1), p), p),
    1739       61829 :                Fp_mul(gel(xa,k), gel(xa,k+1), p), vs);
    1740       31681 :   return FpXV_producttree_dbl(t, n, p);
    1741             : }
    1742             : 
    1743             : static GEN
    1744       32269 : FpX_FpXV_multirem_dbl_tree(GEN P, GEN T, GEN p)
    1745             : {
    1746             :   long i,j,k;
    1747       32269 :   long m = lg(T)-1;
    1748             :   GEN t;
    1749       32269 :   GEN Tp = cgetg(m+1, t_VEC);
    1750       32269 :   gel(Tp, m) = mkvec(P);
    1751       63649 :   for (i=m-1; i>=1; i--)
    1752             :   {
    1753       31380 :     GEN u = gel(T, i);
    1754       31380 :     GEN v = gel(Tp, i+1);
    1755       31380 :     long n = lg(u)-1;
    1756       31380 :     t = cgetg(n+1, t_VEC);
    1757      102394 :     for (j=1, k=1; k<n; j++, k+=2)
    1758             :     {
    1759       71014 :       gel(t, k)   = FpX_rem(gel(v, j), gel(u, k), p);
    1760       71014 :       gel(t, k+1) = FpX_rem(gel(v, j), gel(u, k+1), p);
    1761             :     }
    1762       31380 :     gel(Tp, i) = t;
    1763             :   }
    1764       32269 :   return Tp;
    1765             : }
    1766             : 
    1767             : static GEN
    1768       31681 : FpX_FpV_multieval_tree(GEN P, GEN xa, GEN T, GEN p)
    1769             : {
    1770       31681 :   pari_sp av = avma;
    1771             :   long j,k;
    1772       31681 :   GEN Tp = FpX_FpXV_multirem_dbl_tree(P, T, p);
    1773       31681 :   GEN R = cgetg(lg(xa), t_VEC);
    1774       31681 :   GEN u = gel(T, 1);
    1775       31681 :   GEN v = gel(Tp, 1);
    1776       31681 :   long n = lg(u)-1;
    1777      132108 :   for (j=1, k=1; j<=n; j++)
    1778             :   {
    1779      100427 :     long c, d = degpol(gel(u,j));
    1780      262680 :     for (c=1; c<=d; c++, k++)
    1781      162253 :       gel(R,k) = FpX_eval(gel(v, j), gel(xa,k), p);
    1782             :   }
    1783       31681 :   return gerepileupto(av, R);
    1784             : }
    1785             : 
    1786             : static GEN
    1787          15 : FpVV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, GEN p, long vs)
    1788             : {
    1789          15 :   pari_sp av = avma;
    1790          15 :   long m = lg(T)-1;
    1791          15 :   long i, j, k, ls = lg(s);
    1792          15 :   GEN Tp = cgetg(m+1, t_VEC);
    1793          15 :   GEN t = cgetg(ls, t_VEC);
    1794         241 :   for (j=1, k=1; j<ls; k+=s[j++])
    1795         226 :     if (s[j]==2)
    1796             :     {
    1797          58 :       GEN a = Fp_mul(gel(ya,k), gel(R,k), p);
    1798          58 :       GEN b = Fp_mul(gel(ya,k+1), gel(R,k+1), p);
    1799          58 :       gel(t, j) = deg1pol_shallow(Fp_add(a, b, p),
    1800          58 :               Fp_neg(Fp_add(Fp_mul(gel(xa,k), b, p ),
    1801          58 :               Fp_mul(gel(xa,k+1), a, p), p), p), vs);
    1802             :     }
    1803             :     else
    1804         168 :       gel(t, j) = scalarpol(Fp_mul(gel(ya,k), gel(R,k), p), vs);
    1805          15 :   gel(Tp, 1) = t;
    1806          72 :   for (i=2; i<=m; i++)
    1807             :   {
    1808          57 :     GEN u = gel(T, i-1);
    1809          57 :     GEN t = cgetg(lg(gel(T,i)), t_VEC);
    1810          57 :     GEN v = gel(Tp, i-1);
    1811          57 :     long n = lg(v)-1;
    1812         268 :     for (j=1, k=1; k<n; j++, k+=2)
    1813         211 :       gel(t, j) = FpX_add(ZX_mul(gel(u, k), gel(v, k+1)),
    1814         211 :                           ZX_mul(gel(u, k+1), gel(v, k)), p);
    1815          57 :     gel(Tp, i) = t;
    1816             :   }
    1817          15 :   return gerepilecopy(av, gmael(Tp,m,1));
    1818             : }
    1819             : 
    1820             : GEN
    1821           0 : FpX_FpV_multieval(GEN P, GEN xa, GEN p)
    1822             : {
    1823           0 :   pari_sp av = avma;
    1824           0 :   GEN s = producttree_scheme(lg(xa)-1);
    1825           0 :   GEN T = FpV_producttree(xa, s, p, varn(P));
    1826           0 :   return gerepileupto(av, FpX_FpV_multieval_tree(P, xa, T, p));
    1827             : }
    1828             : 
    1829             : GEN
    1830          22 : FpV_polint(GEN xa, GEN ya, GEN p, long vs)
    1831             : {
    1832          22 :   pari_sp av = avma;
    1833             :   GEN s, T, P, R;
    1834             :   long m;
    1835          22 :   if (lgefint(p) == 3)
    1836             :   {
    1837           7 :     ulong pp = p[2];
    1838           7 :     P = Flv_polint(ZV_to_Flv(xa, pp), ZV_to_Flv(ya, pp), pp, evalvarn(vs));
    1839           7 :     return gerepileupto(av, Flx_to_ZX(P));
    1840             :   }
    1841          15 :   s = producttree_scheme(lg(xa)-1);
    1842          15 :   T = FpV_producttree(xa, s, p, vs);
    1843          15 :   m = lg(T)-1;
    1844          15 :   P = FpX_deriv(gmael(T, m, 1), p);
    1845          15 :   R = FpV_inv(FpX_FpV_multieval_tree(P, xa, T, p), p);
    1846          15 :   return gerepileupto(av, FpVV_polint_tree(T, R, s, xa, ya, p, vs));
    1847             : }
    1848             : 
    1849             : GEN
    1850           0 : FpV_FpM_polint(GEN xa, GEN ya, GEN p, long vs)
    1851             : {
    1852           0 :   pari_sp av = avma;
    1853           0 :   GEN s = producttree_scheme(lg(xa)-1);
    1854           0 :   GEN T = FpV_producttree(xa, s, p, vs);
    1855           0 :   long i, m = lg(T)-1, l = lg(ya)-1;
    1856           0 :   GEN P = FpX_deriv(gmael(T, m, 1), p);
    1857           0 :   GEN R = FpV_inv(FpX_FpV_multieval_tree(P, xa, T, p), p);
    1858           0 :   GEN M = cgetg(l+1, t_VEC);
    1859           0 :   for (i=1; i<=l; i++)
    1860           0 :     gel(M,i) = FpVV_polint_tree(T, R, s, xa, gel(ya,i), p, vs);
    1861           0 :   return gerepileupto(av, M);
    1862             : }
    1863             : 
    1864             : GEN
    1865       31666 : FpV_invVandermonde(GEN L, GEN den, GEN p)
    1866             : {
    1867       31666 :   pari_sp av = avma;
    1868       31666 :   long i, n = lg(L);
    1869             :   GEN M, R;
    1870       31666 :   GEN s = producttree_scheme(n-1);
    1871       31666 :   GEN tree = FpV_producttree(L, s, p, 0);
    1872       31666 :   long m = lg(tree)-1;
    1873       31666 :   GEN T = gmael(tree, m, 1);
    1874       31666 :   R = FpV_inv(FpX_FpV_multieval_tree(FpX_deriv(T, p), L, tree, p), p);
    1875       31665 :   if (den) R = FpC_Fp_mul(R, den, p);
    1876       31666 :   M = cgetg(n, t_MAT);
    1877      193636 :   for (i = 1; i < n; i++)
    1878             :   {
    1879      161970 :     GEN P = FpX_Fp_mul(FpX_div_by_X_x(T, gel(L,i), p, NULL), gel(R,i), p);
    1880      161969 :     gel(M,i) = RgX_to_RgC(P, n-1);
    1881             :   }
    1882       31666 :   return gerepilecopy(av, M);
    1883             : }
    1884             : 
    1885             : static GEN
    1886         588 : FpXV_producttree(GEN xa, GEN s, GEN p)
    1887             : {
    1888         588 :   long n = lg(xa)-1;
    1889         588 :   long j, k, ls = lg(s);
    1890         588 :   GEN t = cgetg(ls, t_VEC);
    1891        3444 :   for (j=1, k=1; j<ls; k+=s[j++])
    1892        2856 :     gel(t, j) = s[j] == 1 ?
    1893        2856 :              gel(xa,k): FpX_mul(gel(xa,k),gel(xa,k+1),p);
    1894         588 :   return FpXV_producttree_dbl(t, n, p);
    1895             : }
    1896             : 
    1897             : static GEN
    1898         588 : FpX_FpXV_multirem_tree(GEN P, GEN xa, GEN T, GEN s, GEN p)
    1899             : {
    1900         588 :   pari_sp av = avma;
    1901         588 :   long j, k, ls = lg(s);
    1902         588 :   GEN Tp = FpX_FpXV_multirem_dbl_tree(P, T, p);
    1903         588 :   GEN R = cgetg(lg(xa), t_VEC);
    1904         588 :   GEN v = gel(Tp, 1);
    1905        3444 :   for (j=1, k=1; j<ls; k+=s[j++])
    1906             :   {
    1907        2856 :     gel(R,k) = FpX_rem(gel(v, j), gel(xa,k), p);
    1908        2856 :     if (s[j] == 2)
    1909        1050 :       gel(R,k+1) = FpX_rem(gel(v, j), gel(xa,k+1), p);
    1910             :   }
    1911         588 :   return gerepileupto(av, R);
    1912             : }
    1913             : 
    1914             : GEN
    1915           0 : FpX_FpXV_multirem(GEN P, GEN xa, GEN p)
    1916             : {
    1917           0 :   pari_sp av = avma;
    1918           0 :   GEN s = producttree_scheme(lg(xa)-1);
    1919           0 :   GEN T = FpXV_producttree(xa, s, p);
    1920           0 :   return gerepileupto(av, FpX_FpXV_multirem_tree(P, xa, T, s, p));
    1921             : }
    1922             : 
    1923             : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
    1924             : static GEN
    1925         588 : FpXV_chinese_tree(GEN A, GEN P, GEN T, GEN R, GEN s, GEN p)
    1926             : {
    1927         588 :   long m = lg(T)-1, ls = lg(s);
    1928             :   long i,j,k;
    1929         588 :   GEN Tp = cgetg(m+1, t_VEC);
    1930         588 :   GEN M = gel(T, 1);
    1931         588 :   GEN t = cgetg(lg(M), t_VEC);
    1932        3444 :   for (j=1, k=1; j<ls; k+=s[j++])
    1933        2856 :     if (s[j] == 2)
    1934             :     {
    1935        1050 :       pari_sp av = avma;
    1936        1050 :       GEN a = FpX_mul(gel(A,k), gel(R,k), p), b = FpX_mul(gel(A,k+1), gel(R,k+1), p);
    1937        1050 :       GEN tj = FpX_rem(FpX_add(FpX_mul(gel(P,k), b, p),
    1938        1050 :             FpX_mul(gel(P,k+1), a, p), p), gel(M,j), p);
    1939        1050 :       gel(t, j) = gerepileupto(av, tj);
    1940             :     }
    1941             :     else
    1942        1806 :       gel(t, j) = FpX_rem(FpX_mul(gel(A,k), gel(R,k), p), gel(M, j), p);
    1943         588 :   gel(Tp, 1) = t;
    1944        1890 :   for (i=2; i<=m; i++)
    1945             :   {
    1946        1302 :     GEN u = gel(T, i-1), M = gel(T, i);
    1947        1302 :     GEN t = cgetg(lg(M), t_VEC);
    1948        1302 :     GEN v = gel(Tp, i-1);
    1949        1302 :     long n = lg(v)-1;
    1950        3570 :     for (j=1, k=1; k<n; j++, k+=2)
    1951             :     {
    1952        2268 :       pari_sp av = avma;
    1953        2268 :       gel(t, j) = gerepileupto(av, FpX_rem(FpX_add(FpX_mul(gel(u, k), gel(v, k+1), p),
    1954        2268 :               FpX_mul(gel(u, k+1), gel(v, k), p), p), gel(M, j), p));
    1955             :     }
    1956        1302 :     if (k==n) gel(t, j) = gel(v, k);
    1957        1302 :     gel(Tp, i) = t;
    1958             :   }
    1959         588 :   return gmael(Tp,m,1);
    1960             : }
    1961             : 
    1962             : static GEN
    1963         588 : FpXV_sqr(GEN x, GEN p)
    1964        4494 : { pari_APPLY_type(t_VEC, FpX_sqr(gel(x,i), p)) }
    1965             : 
    1966             : static GEN
    1967        7602 : FpXT_sqr(GEN x, GEN p)
    1968             : {
    1969        7602 :   if (typ(x) == t_POL)
    1970        5124 :     return FpX_sqr(x, p);
    1971        9492 :   pari_APPLY_type(t_VEC, FpXT_sqr(gel(x,i), p))
    1972             : }
    1973             : 
    1974             : static GEN
    1975         588 : FpXV_invdivexact(GEN x, GEN y, GEN p)
    1976        4494 : { pari_APPLY_type(t_VEC, FpXQ_inv(FpX_div(gel(x,i), gel(y,i),p), gel(y,i),p)) }
    1977             : 
    1978             : static GEN
    1979         588 : FpXV_chinesetree(GEN P, GEN T, GEN s, GEN p)
    1980             : {
    1981         588 :   GEN T2 = FpXT_sqr(T, p), P2 = FpXV_sqr(P, p);
    1982         588 :   GEN mod = gmael(T,lg(T)-1,1);
    1983         588 :   return FpXV_invdivexact(FpX_FpXV_multirem_tree(mod, P2, T2, s, p), P, p);
    1984             : }
    1985             : 
    1986             : static GEN
    1987         588 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
    1988             : {
    1989         588 :   if (!pt_mod)
    1990         588 :     return gerepileupto(av, a);
    1991             :   else
    1992             :   {
    1993           0 :     GEN mod = gmael(T, lg(T)-1, 1);
    1994           0 :     gerepileall(av, 2, &a, &mod);
    1995           0 :     *pt_mod = mod;
    1996           0 :     return a;
    1997             :   }
    1998             : }
    1999             : 
    2000             : GEN
    2001         588 : FpXV_chinese(GEN A, GEN P, GEN p, GEN *pt_mod)
    2002             : {
    2003         588 :   pari_sp av = avma;
    2004         588 :   GEN s = producttree_scheme(lg(P)-1);
    2005         588 :   GEN T = FpXV_producttree(P, s, p);
    2006         588 :   GEN R = FpXV_chinesetree(P, T, s, p);
    2007         588 :   GEN a = FpXV_chinese_tree(A, P, T, R, s, p);
    2008         588 :   return gc_chinese(av, T, a, pt_mod);
    2009             : }
    2010             : 
    2011             : /***********************************************************************/
    2012             : /**                                                                   **/
    2013             : /**                              FpXQ                                 **/
    2014             : /**                                                                   **/
    2015             : /***********************************************************************/
    2016             : 
    2017             : /* FpXQ are elements of Fp[X]/(T), represented by FpX*/
    2018             : 
    2019             : GEN
    2020    17805973 : FpXQ_red(GEN x, GEN T, GEN p)
    2021             : {
    2022    17805973 :   GEN z = FpX_red(x,p);
    2023    17775626 :   return FpX_rem(z, T,p);
    2024             : }
    2025             : 
    2026             : GEN
    2027    11074575 : FpXQ_mul(GEN x,GEN y,GEN T,GEN p)
    2028             : {
    2029    11074575 :   GEN z = FpX_mul(x,y,p);
    2030    11074934 :   return FpX_rem(z, T, p);
    2031             : }
    2032             : 
    2033             : GEN
    2034     6373038 : FpXQ_sqr(GEN x, GEN T, GEN p)
    2035             : {
    2036     6373038 :   GEN z = FpX_sqr(x,p);
    2037     6371828 :   return FpX_rem(z, T, p);
    2038             : }
    2039             : 
    2040             : /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
    2041             :  * return lift(1 / (x mod (p,pol))) */
    2042             : GEN
    2043     1098641 : FpXQ_invsafe(GEN x, GEN y, GEN p)
    2044             : {
    2045     1098641 :   GEN V, z = FpX_extgcd(get_FpX_mod(y), x, p, NULL, &V);
    2046     1098656 :   if (degpol(z)) return NULL;
    2047     1098656 :   z = Fp_invsafe(gel(z,2), p);
    2048     1098606 :   if (!z) return NULL;
    2049     1098606 :   return FpX_Fp_mul(V, z, p);
    2050             : }
    2051             : 
    2052             : GEN
    2053     1098641 : FpXQ_inv(GEN x,GEN T,GEN p)
    2054             : {
    2055     1098641 :   pari_sp av = avma;
    2056     1098641 :   GEN U = FpXQ_invsafe(x, T, p);
    2057     1098593 :   if (!U) pari_err_INV("FpXQ_inv",x);
    2058     1098593 :   return gerepileupto(av, U);
    2059             : }
    2060             : 
    2061             : GEN
    2062      531956 : FpXQ_div(GEN x,GEN y,GEN T,GEN p)
    2063             : {
    2064      531956 :   pari_sp av = avma;
    2065      531956 :   return gerepileupto(av, FpXQ_mul(x,FpXQ_inv(y,T,p),T,p));
    2066             : }
    2067             : 
    2068             : static GEN
    2069     3041295 : _FpXQ_add(void *data, GEN x, GEN y)
    2070             : {
    2071             :   (void) data;
    2072     3041295 :   return ZX_add(x, y);
    2073             : }
    2074             : static GEN
    2075       52941 : _FpXQ_sub(void *data, GEN x, GEN y)
    2076             : {
    2077             :   (void) data;
    2078       52941 :   return ZX_sub(x, y);
    2079             : }
    2080             : static GEN
    2081     3458820 : _FpXQ_cmul(void *data, GEN P, long a, GEN x)
    2082             : {
    2083             :   (void) data;
    2084     3458820 :   return ZX_Z_mul(x, gel(P,a+2));
    2085             : }
    2086             : static GEN
    2087     5457753 : _FpXQ_sqr(void *data, GEN x)
    2088             : {
    2089     5457753 :   struct _FpXQ *D = (struct _FpXQ*)data;
    2090     5457753 :   return FpXQ_sqr(x, D->T, D->p);
    2091             : }
    2092             : static GEN
    2093     1945193 : _FpXQ_mul(void *data, GEN x, GEN y)
    2094             : {
    2095     1945193 :   struct _FpXQ *D = (struct _FpXQ*)data;
    2096     1945193 :   return FpXQ_mul(x,y, D->T, D->p);
    2097             : }
    2098             : static GEN
    2099        4095 : _FpXQ_zero(void *data)
    2100             : {
    2101        4095 :   struct _FpXQ *D = (struct _FpXQ*)data;
    2102        4095 :   return pol_0(get_FpX_var(D->T));
    2103             : }
    2104             : static GEN
    2105      933495 : _FpXQ_one(void *data)
    2106             : {
    2107      933495 :   struct _FpXQ *D = (struct _FpXQ*)data;
    2108      933495 :   return pol_1(get_FpX_var(D->T));
    2109             : }
    2110             : static GEN
    2111      939104 : _FpXQ_red(void *data, GEN x)
    2112             : {
    2113      939104 :   struct _FpXQ *D = (struct _FpXQ*)data;
    2114      939104 :   return FpX_red(x,D->p);
    2115             : }
    2116             : 
    2117             : static struct bb_algebra FpXQ_algebra = { _FpXQ_red, _FpXQ_add, _FpXQ_sub,
    2118             :        _FpXQ_mul, _FpXQ_sqr, _FpXQ_one, _FpXQ_zero };
    2119             : 
    2120             : const struct bb_algebra *
    2121       10199 : get_FpXQ_algebra(void **E, GEN T, GEN p)
    2122             : {
    2123       10199 :   GEN z = new_chunk(sizeof(struct _FpXQ));
    2124       10199 :   struct _FpXQ *e = (struct _FpXQ *) z;
    2125       10199 :   e->T = FpX_get_red(T, p);
    2126       10199 :   e->p  = p; *E = (void*)e;
    2127       10199 :   return &FpXQ_algebra;
    2128             : }
    2129             : 
    2130             : static GEN
    2131           0 : _FpX_red(void *E, GEN x)
    2132           0 : { struct _FpX *D = (struct _FpX*)E; return FpX_red(x,D->p); }
    2133             : 
    2134             : static GEN
    2135           0 : _FpX_zero(void *E)
    2136           0 : { struct _FpX *D = (struct _FpX *)E; return pol_0(D->v); }
    2137             : 
    2138             : 
    2139             : static struct bb_algebra FpX_algebra = { _FpX_red, _FpXQ_add, _FpXQ_sub,
    2140             :        _FpX_mul, _FpX_sqr, _FpX_one, _FpX_zero };
    2141             : 
    2142             : const struct bb_algebra *
    2143           0 : get_FpX_algebra(void **E, GEN p, long v)
    2144             : {
    2145           0 :   GEN z = new_chunk(sizeof(struct _FpX));
    2146           0 :   struct _FpX *e = (struct _FpX *) z;
    2147           0 :   e->p  = p; e->v = v; *E = (void*)e;
    2148           0 :   return &FpX_algebra;
    2149             : }
    2150             : 
    2151             : /* x,pol in Z[X], p in Z, n in Z, compute lift(x^n mod (p, pol)) */
    2152             : GEN
    2153     1018502 : FpXQ_pow(GEN x, GEN n, GEN T, GEN p)
    2154             : {
    2155             :   struct _FpXQ D;
    2156             :   pari_sp av;
    2157     1018502 :   long s = signe(n);
    2158             :   GEN y;
    2159     1018502 :   if (!s) return pol_1(varn(x));
    2160     1017898 :   if (is_pm1(n)) /* +/- 1 */
    2161       36032 :     return (s < 0)? FpXQ_inv(x,T,p): FpXQ_red(x,T,p);
    2162      981864 :   av = avma;
    2163      981864 :   if (!is_bigint(p))
    2164             :   {
    2165      659598 :     ulong pp = to_Flxq(&x, &T, p);
    2166      659601 :     y = Flxq_pow(x, n, T, pp);
    2167      659584 :     return Flx_to_ZX_inplace(gerepileuptoleaf(av, y));
    2168             :   }
    2169      322267 :   if (s < 0) x = FpXQ_inv(x,T,p);
    2170      322267 :   D.p = p; D.T = FpX_get_red(T,p);
    2171      322267 :   y = gen_pow_i(x, n, (void*)&D, &_FpXQ_sqr, &_FpXQ_mul);
    2172      322267 :   return gerepilecopy(av, y);
    2173             : }
    2174             : 
    2175             : GEN /*Assume n is very small*/
    2176      608794 : FpXQ_powu(GEN x, ulong n, GEN T, GEN p)
    2177             : {
    2178             :   struct _FpXQ D;
    2179             :   pari_sp av;
    2180             :   GEN y;
    2181      608794 :   if (!n) return pol_1(varn(x));
    2182      608794 :   if (n==1) return FpXQ_red(x,T,p);
    2183      206442 :   av = avma;
    2184      206442 :   if (!is_bigint(p))
    2185             :   {
    2186      197927 :     ulong pp = to_Flxq(&x, &T, p);
    2187      197926 :     y = Flxq_powu(x, n, T, pp);
    2188      197923 :     return Flx_to_ZX_inplace(gerepileuptoleaf(av, y));
    2189             :   }
    2190        8524 :   D.T = FpX_get_red(T, p); D.p = p;
    2191        8524 :   y = gen_powu_i(x, n, (void*)&D, &_FpXQ_sqr, &_FpXQ_mul);
    2192        8524 :   return gerepilecopy(av, y);
    2193             : }
    2194             : 
    2195             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2196             : GEN
    2197      383518 : FpXQ_powers(GEN x, long l, GEN T, GEN p)
    2198             : {
    2199             :   struct _FpXQ D;
    2200             :   int use_sqr;
    2201      383518 :   if (l>2 && lgefint(p) == 3) {
    2202      209547 :     pari_sp av = avma;
    2203      209547 :     ulong pp = to_Flxq(&x, &T, p);
    2204      209549 :     GEN z = FlxV_to_ZXV(Flxq_powers(x, l, T, pp));
    2205      209550 :     return gerepileupto(av, z);
    2206             :   }
    2207      173971 :   use_sqr = 2*degpol(x)>=get_FpX_degree(T);
    2208      173972 :   D.T = FpX_get_red(T,p); D.p = p;
    2209      173972 :   return gen_powers(x, l, use_sqr, (void*)&D, &_FpXQ_sqr, &_FpXQ_mul,&_FpXQ_one);
    2210             : }
    2211             : 
    2212             : GEN
    2213       66290 : FpXQ_matrix_pow(GEN y, long n, long m, GEN P, GEN l)
    2214             : {
    2215       66290 :   return RgXV_to_RgM(FpXQ_powers(y,m-1,P,l),n);
    2216             : }
    2217             : 
    2218             : GEN
    2219      445445 : FpX_Frobenius(GEN T, GEN p)
    2220             : {
    2221      445445 :   return FpXQ_pow(pol_x(get_FpX_var(T)), p, T, p);
    2222             : }
    2223             : 
    2224             : GEN
    2225       31494 : FpX_matFrobenius(GEN T, GEN p)
    2226             : {
    2227       31494 :   long n = get_FpX_degree(T);
    2228       31494 :   return FpXQ_matrix_pow(FpX_Frobenius(T, p), n, n, T, p);
    2229             : }
    2230             : 
    2231             : GEN
    2232      411618 : FpX_FpXQV_eval(GEN Q, GEN x, GEN T, GEN p)
    2233             : {
    2234             :   struct _FpXQ D;
    2235      411618 :   D.T = FpX_get_red(T,p); D.p = p;
    2236      411632 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)&D,&FpXQ_algebra,_FpXQ_cmul);
    2237             : }
    2238             : 
    2239             : GEN
    2240      796051 : FpX_FpXQ_eval(GEN Q, GEN x, GEN T, GEN p)
    2241             : {
    2242             :   struct _FpXQ D;
    2243             :   int use_sqr;
    2244      796051 :   if (lgefint(p) == 3)
    2245             :   {
    2246      789403 :     pari_sp av = avma;
    2247      789403 :     ulong pp = to_Flxq(&x, &T, p);
    2248      789426 :     GEN z = Flx_Flxq_eval(ZX_to_Flx(Q, pp), x, T, pp);
    2249      789423 :     return Flx_to_ZX_inplace(gerepileuptoleaf(av, z));
    2250             :   }
    2251        6648 :   use_sqr = 2*degpol(x) >= get_FpX_degree(T);
    2252        6648 :   D.T = FpX_get_red(T,p); D.p = p;
    2253        6648 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&D,&FpXQ_algebra,_FpXQ_cmul);
    2254             : }
    2255             : 
    2256             : GEN
    2257        1470 : FpXC_FpXQV_eval(GEN x, GEN v, GEN T, GEN p)
    2258        8316 : { pari_APPLY_type(t_COL, FpX_FpXQV_eval(gel(x,i), v, T, p)) }
    2259             : 
    2260             : GEN
    2261         315 : FpXM_FpXQV_eval(GEN x, GEN v, GEN T, GEN p)
    2262        1197 : { pari_APPLY_same(FpXC_FpXQV_eval(gel(x,i), v, T, p)) }
    2263             : 
    2264             : GEN
    2265         588 : FpXC_FpXQ_eval(GEN x, GEN F, GEN T, GEN p)
    2266             : {
    2267         588 :   long d = brent_kung_optpow(RgXV_maxdegree(x), lg(x)-1, 1);
    2268         588 :   GEN Fp = FpXQ_powers(F, d, T, p);
    2269         588 :   return FpXC_FpXQV_eval(x, Fp, T, p);
    2270             : }
    2271             : 
    2272             : GEN
    2273        1771 : FpXQ_autpowers(GEN aut, long f, GEN T, GEN p)
    2274             : {
    2275        1771 :   pari_sp av = avma;
    2276        1771 :   long n = get_FpX_degree(T);
    2277        1771 :   long i, nautpow = brent_kung_optpow(n-1,f-2,1);
    2278        1771 :   long v = get_FpX_var(T);
    2279             :   GEN autpow, V;
    2280        1771 :   T = FpX_get_red(T, p);
    2281        1771 :   autpow = FpXQ_powers(aut, nautpow,T,p);
    2282        1771 :   V = cgetg(f + 2, t_VEC);
    2283        1771 :   gel(V,1) = pol_x(v); if (f==0) return gerepileupto(av, V);
    2284        1771 :   gel(V,2) = gcopy(aut);
    2285        6272 :   for (i = 3; i <= f+1; i++)
    2286        4501 :     gel(V,i) = FpX_FpXQV_eval(gel(V,i-1),autpow,T,p);
    2287        1771 :   return gerepileupto(av, V);
    2288             : }
    2289             : 
    2290             : static GEN
    2291        4558 : FpXQ_autpow_sqr(void *E, GEN x)
    2292             : {
    2293        4558 :   struct _FpXQ *D = (struct _FpXQ*)E;
    2294        4558 :   return FpX_FpXQ_eval(x, x, D->T, D->p);
    2295             : }
    2296             : 
    2297             : static GEN
    2298          42 : FpXQ_autpow_msqr(void *E, GEN x)
    2299             : {
    2300          42 :   struct _FpXQ *D = (struct _FpXQ*)E;
    2301          42 :   return FpX_FpXQV_eval(FpXQ_autpow_sqr(E, x), D->aut, D->T, D->p);
    2302             : }
    2303             : 
    2304             : GEN
    2305        4464 : FpXQ_autpow(GEN x, ulong n, GEN T, GEN p)
    2306             : {
    2307        4464 :   pari_sp av = avma;
    2308             :   struct _FpXQ D;
    2309             :   long d;
    2310        4464 :   if (n==0) return FpX_rem(pol_x(varn(x)), T, p);
    2311        4464 :   if (n==1) return FpX_rem(x, T, p);
    2312        4243 :   D.T = FpX_get_red(T, p); D.p = p;
    2313        4243 :   d = brent_kung_optpow(degpol(T), hammingl(n)-1, 1);
    2314        4243 :   D.aut = FpXQ_powers(x, d, T, p);
    2315        4243 :   x = gen_powu_fold(x,n,(void*)&D,FpXQ_autpow_sqr,FpXQ_autpow_msqr);
    2316        4243 :   return gerepilecopy(av, x);
    2317             : }
    2318             : 
    2319             : static GEN
    2320         360 : FpXQ_auttrace_mul(void *E, GEN x, GEN y)
    2321             : {
    2322         360 :   struct _FpXQ *D = (struct _FpXQ*)E;
    2323         360 :   GEN T = D->T, p = D->p;
    2324         360 :   GEN phi1 = gel(x,1), a1 = gel(x,2);
    2325         360 :   GEN phi2 = gel(y,1), a2 = gel(y,2);
    2326         360 :   ulong d = brent_kung_optpow(maxss(degpol(phi2),degpol(a2)),2,1);
    2327         360 :   GEN V1 = FpXQ_powers(phi1, d, T, p);
    2328         360 :   GEN phi3 = FpX_FpXQV_eval(phi2, V1, T, p);
    2329         360 :   GEN aphi = FpX_FpXQV_eval(a2, V1, T, p);
    2330         360 :   GEN a3 = FpX_add(a1, aphi, p);
    2331         360 :   return mkvec2(phi3, a3);
    2332             : }
    2333             : 
    2334             : static GEN
    2335         317 : FpXQ_auttrace_sqr(void *E, GEN x)
    2336         317 : { return FpXQ_auttrace_mul(E, x, x); }
    2337             : 
    2338             : GEN
    2339         427 : FpXQ_auttrace(GEN x, ulong n, GEN T, GEN p)
    2340             : {
    2341         427 :   pari_sp av = avma;
    2342             :   struct _FpXQ D;
    2343         427 :   D.T = FpX_get_red(T, p); D.p = p;
    2344         427 :   x = gen_powu_i(x,n,(void*)&D,FpXQ_auttrace_sqr,FpXQ_auttrace_mul);
    2345         427 :   return gerepilecopy(av, x);
    2346             : }
    2347             : 
    2348             : static GEN
    2349        3078 : FpXQ_autsum_mul(void *E, GEN x, GEN y)
    2350             : {
    2351        3078 :   struct _FpXQ *D = (struct _FpXQ*)E;
    2352        3078 :   GEN T = D->T, p = D->p;
    2353        3078 :   GEN phi1 = gel(x,1), a1 = gel(x,2);
    2354        3078 :   GEN phi2 = gel(y,1), a2 = gel(y,2);
    2355        3078 :   ulong d = brent_kung_optpow(maxss(degpol(phi2),degpol(a2)),2,1);
    2356        3078 :   GEN V1 = FpXQ_powers(phi1, d, T, p);
    2357        3078 :   GEN phi3 = FpX_FpXQV_eval(phi2, V1, T, p);
    2358        3078 :   GEN aphi = FpX_FpXQV_eval(a2, V1, T, p);
    2359        3078 :   GEN a3 = FpXQ_mul(a1, aphi, T, p);
    2360        3078 :   return mkvec2(phi3, a3);
    2361             : }
    2362             : static GEN
    2363        2918 : FpXQ_autsum_sqr(void *E, GEN x)
    2364        2918 : { return FpXQ_autsum_mul(E, x, x); }
    2365             : 
    2366             : GEN
    2367        2792 : FpXQ_autsum(GEN x, ulong n, GEN T, GEN p)
    2368             : {
    2369        2792 :   pari_sp av = avma;
    2370             :   struct _FpXQ D;
    2371        2792 :   D.T = FpX_get_red(T, p); D.p = p;
    2372        2792 :   x = gen_powu_i(x,n,(void*)&D,FpXQ_autsum_sqr,FpXQ_autsum_mul);
    2373        2792 :   return gerepilecopy(av, x);
    2374             : }
    2375             : 
    2376             : static GEN
    2377         315 : FpXQM_autsum_mul(void *E, GEN x, GEN y)
    2378             : {
    2379         315 :   struct _FpXQ *D = (struct _FpXQ*)E;
    2380         315 :   GEN T = D->T, p = D->p;
    2381         315 :   GEN phi1 = gel(x,1), a1 = gel(x,2);
    2382         315 :   GEN phi2 = gel(y,1), a2 = gel(y,2);
    2383         315 :   long g = lg(a2)-1, dT = get_FpX_degree(T);
    2384         315 :   ulong d = brent_kung_optpow(dT-1, g*g+1, 1);
    2385         315 :   GEN V1 = FpXQ_powers(phi1, d, T, p);
    2386         315 :   GEN phi3 = FpX_FpXQV_eval(phi2, V1, T, p);
    2387         315 :   GEN aphi = FpXM_FpXQV_eval(a2, V1, T, p);
    2388         315 :   GEN a3 = FqM_mul(a1, aphi, T, p);
    2389         315 :   return mkvec2(phi3, a3);
    2390             : }
    2391             : static GEN
    2392         217 : FpXQM_autsum_sqr(void *E, GEN x)
    2393         217 : { return FpXQM_autsum_mul(E, x, x); }
    2394             : 
    2395             : GEN
    2396         147 : FpXQM_autsum(GEN x, ulong n, GEN T, GEN p)
    2397             : {
    2398         147 :   pari_sp av = avma;
    2399             :   struct _FpXQ D;
    2400         147 :   D.T = FpX_get_red(T, p); D.p = p;
    2401         147 :   x = gen_powu_i(x, n, (void*)&D, FpXQM_autsum_sqr, FpXQM_autsum_mul);
    2402         147 :   return gerepilecopy(av, x);
    2403             : }
    2404             : 
    2405             : static long
    2406        4414 : bounded_order(GEN p, GEN b, long k)
    2407             : {
    2408             :   long i;
    2409        4414 :   GEN a=modii(p,b);
    2410        9538 :   for(i=1;i<k;i++)
    2411             :   {
    2412        7916 :     if (equali1(a))
    2413        2792 :       return i;
    2414        5124 :     a = Fp_mul(a,p,b);
    2415             :   }
    2416        1622 :   return 0;
    2417             : }
    2418             : 
    2419             : /* n = (p^d-a)\b
    2420             :  * b = bb*p^vb
    2421             :  * p^k = 1 [bb]
    2422             :  * d = m*k+r+vb
    2423             :  * u = (p^k-1)/bb;
    2424             :  * v = (p^(r+vb)-a)/b;
    2425             :  * w = (p^(m*k)-1)/(p^k-1)
    2426             :  * n = p^r*w*u+v
    2427             :  * w*u = p^vb*(p^(m*k)-1)/b
    2428             :  * n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b */
    2429             : static GEN
    2430      296988 : FpXQ_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, GEN p)
    2431             : {
    2432      296988 :   pari_sp av=avma;
    2433      296988 :   long d = get_FpX_degree(T);
    2434      296988 :   GEN an = absi_shallow(n), z, q;
    2435      296988 :   if (cmpii(an,p)<0 || cmpis(an,d)<=0) return FpXQ_pow(x, n, T, p);
    2436        4421 :   q = powiu(p, d);
    2437        4421 :   if (dvdii(q, n))
    2438             :   {
    2439           0 :     long vn = logint(an,p);
    2440           0 :     GEN autvn = vn==1 ? aut: FpXQ_autpow(aut,vn,T,p);
    2441           0 :     z = FpX_FpXQ_eval(x,autvn,T,p);
    2442             :   } else
    2443             :   {
    2444        4421 :     GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
    2445             :     GEN bb, u, v, autk;
    2446        4421 :     long vb = Z_pvalrem(b,p,&bb);
    2447        4421 :     long m, r, k = is_pm1(bb) ? 1 : bounded_order(p,bb,d);
    2448        4421 :     if (!k || d-vb<k) return FpXQ_pow(x,n, T, p);
    2449        2799 :     m = (d-vb)/k; r = (d-vb)%k;
    2450        2799 :     u = diviiexact(subiu(powiu(p,k),1),bb);
    2451        2799 :     v = diviiexact(subii(powiu(p,r+vb),a),b);
    2452        2799 :     autk = k==1 ? aut: FpXQ_autpow(aut,k,T,p);
    2453        2799 :     if (r)
    2454             :     {
    2455          14 :       GEN autr = r==1 ? aut: FpXQ_autpow(aut,r,T,p);
    2456          14 :       z = FpX_FpXQ_eval(x,autr,T,p);
    2457        2785 :     } else z = x;
    2458        2799 :     if (m > 1) z = gel(FpXQ_autsum(mkvec2(autk, z), m, T, p), 2);
    2459        2799 :     if (!is_pm1(u)) z = FpXQ_pow(z, u, T, p);
    2460        2799 :     if (signe(v)) z = FpXQ_mul(z, FpXQ_pow(x, v, T, p), T, p);
    2461             :   }
    2462        2799 :   return gerepileupto(av,signe(n)>0 ? z : FpXQ_inv(z,T,p));
    2463             : }
    2464             : 
    2465             : /* assume T irreducible mod p */
    2466             : int
    2467      401380 : FpXQ_issquare(GEN x, GEN T, GEN p)
    2468             : {
    2469             :   pari_sp av;
    2470      401380 :   if (lg(x) == 2 || absequalui(2, p)) return 1;
    2471      401366 :   if (lg(x) == 3) return Fq_issquare(gel(x,2), T, p);
    2472      363835 :   av = avma; /* Ng = g^((q-1)/(p-1)) */
    2473      363835 :   return gc_bool(av, kronecker(FpXQ_norm(x,T,p), p) != -1);
    2474             : }
    2475             : int
    2476     1336138 : Fp_issquare(GEN x, GEN p)
    2477     1336138 : { return absequalui(2, p) || kronecker(x, p) != -1; }
    2478             : /* assume T irreducible mod p */
    2479             : int
    2480     1630934 : Fq_issquare(GEN x, GEN T, GEN p)
    2481             : {
    2482     1630934 :   if (typ(x) != t_INT) return FpXQ_issquare(x, T, p);
    2483     1233642 :   return (T && ! odd(get_FpX_degree(T))) || Fp_issquare(x, p);
    2484             : }
    2485             : 
    2486             : long
    2487          70 : Fq_ispower(GEN x, GEN K, GEN T, GEN p)
    2488             : {
    2489          70 :   pari_sp av = avma;
    2490             :   long d;
    2491             :   GEN Q;
    2492          70 :   if (equaliu(K,2)) return Fq_issquare(x, T, p);
    2493           0 :   if (!T) return Fp_ispower(x, K, p);
    2494           0 :   d = get_FpX_degree(T);
    2495           0 :   if (typ(x) == t_INT && !umodui(d, K)) return 1;
    2496           0 :   Q = subiu(powiu(p,d), 1);
    2497           0 :   Q = diviiexact(Q, gcdii(Q, K));
    2498           0 :   d = gequal1(Fq_pow(x, Q, T,p));
    2499           0 :   return gc_long(av, d);
    2500             : }
    2501             : 
    2502             : /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
    2503             : GEN
    2504      544143 : Fp_FpXQ_log(GEN a, GEN g, GEN o, GEN T, GEN p)
    2505             : {
    2506      544143 :   pari_sp av = avma;
    2507             :   GEN q,n_q,ord,ordp, op;
    2508             : 
    2509      544143 :   if (equali1(a)) return gen_0;
    2510             :   /* p > 2 */
    2511             : 
    2512        6974 :   ordp = subiu(p, 1); /* even */
    2513        6974 :   ord  = get_arith_Z(o);
    2514        6946 :   if (!ord) ord = T? subiu(powiu(p, get_FpX_degree(T)), 1): ordp;
    2515        6946 :   if (equalii(a, ordp)) /* -1 */
    2516        5021 :     return gerepileuptoint(av, shifti(ord,-1));
    2517        1925 :   ordp = gcdii(ordp,ord);
    2518        1925 :   op = typ(o)==t_MAT ? famat_Z_gcd(o,ordp) : ordp;
    2519             : 
    2520        1925 :   q = NULL;
    2521        1925 :   if (T)
    2522             :   { /* we want < g > = Fp^* */
    2523        1925 :     if (!equalii(ord,ordp)) {
    2524        1903 :       q = diviiexact(ord,ordp);
    2525        1903 :       g = FpXQ_pow(g,q,T,p);
    2526             :     }
    2527        1925 :     g = constant_coeff(g);
    2528             :   }
    2529        1925 :   n_q = Fp_log(a,g,op,p);
    2530        1925 :   if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
    2531        1925 :   if (q) n_q = mulii(q, n_q);
    2532        1925 :   return gerepileuptoint(av, n_q);
    2533             : }
    2534             : 
    2535             : static GEN
    2536      281774 : _FpXQ_pow(void *data, GEN x, GEN n)
    2537             : {
    2538      281774 :   struct _FpXQ *D = (struct _FpXQ*)data;
    2539      281774 :   return FpXQ_pow_Frobenius(x,n, D->aut, D->T, D->p);
    2540             : }
    2541             : 
    2542             : static GEN
    2543         738 : _FpXQ_rand(void *data)
    2544             : {
    2545         738 :   pari_sp av=avma;
    2546         738 :   struct _FpXQ *D = (struct _FpXQ*)data;
    2547             :   GEN z;
    2548             :   do
    2549             :   {
    2550         738 :     set_avma(av);
    2551         738 :     z=random_FpX(get_FpX_degree(D->T),get_FpX_var(D->T),D->p);
    2552         738 :   } while (!signe(z));
    2553         738 :   return z;
    2554             : }
    2555             : 
    2556             : static GEN
    2557         420 : _FpXQ_easylog(void *E, GEN a, GEN g, GEN ord)
    2558             : {
    2559         420 :   struct _FpXQ *s=(struct _FpXQ*) E;
    2560         420 :   if (degpol(a)) return NULL;
    2561         193 :   return Fp_FpXQ_log(constant_coeff(a),g,ord,s->T,s->p);
    2562             : }
    2563             : 
    2564             : static const struct bb_group FpXQ_star={_FpXQ_mul,_FpXQ_pow,_FpXQ_rand,hash_GEN,ZX_equal,ZX_equal1,_FpXQ_easylog};
    2565             : 
    2566             : const struct bb_group *
    2567        2809 : get_FpXQ_star(void **E, GEN T, GEN p)
    2568             : {
    2569        2809 :   struct _FpXQ *e = (struct _FpXQ *) stack_malloc(sizeof(struct _FpXQ));
    2570        2809 :   e->T = T; e->p  = p; e->aut =  FpX_Frobenius(T, p);
    2571        2809 :   *E = (void*)e; return &FpXQ_star;
    2572             : }
    2573             : 
    2574             : GEN
    2575        1816 : FpXQ_order(GEN a, GEN ord, GEN T, GEN p)
    2576             : {
    2577        1816 :   if (lgefint(p)==3)
    2578             :   {
    2579           0 :     pari_sp av=avma;
    2580           0 :     ulong pp = to_Flxq(&a, &T, p);
    2581           0 :     GEN z = Flxq_order(a, ord, T, pp);
    2582           0 :     return gerepileuptoint(av,z);
    2583             :   }
    2584             :   else
    2585             :   {
    2586             :     void *E;
    2587        1816 :     const struct bb_group *S = get_FpXQ_star(&E,T,p);
    2588        1816 :     return gen_order(a,ord,E,S);
    2589             :   }
    2590             : }
    2591             : 
    2592             : GEN
    2593      708205 : FpXQ_log(GEN a, GEN g, GEN ord, GEN T, GEN p)
    2594             : {
    2595      708205 :   pari_sp av=avma;
    2596      708205 :   if (lgefint(p)==3)
    2597             :   {
    2598      708070 :     if (uel(p,2) == 2)
    2599             :     {
    2600      543690 :       GEN z = F2xq_log(ZX_to_F2x(a), ZX_to_F2x(g), ord,
    2601             :                                      ZX_to_F2x(get_FpX_mod(T)));
    2602      543690 :       return gerepileuptoleaf(av, z);
    2603             :     }
    2604             :     else
    2605             :     {
    2606      164380 :       ulong pp = to_Flxq(&a, &T, p);
    2607      164380 :       GEN z = Flxq_log(a, ZX_to_Flx(g, pp), ord, T, pp);
    2608      164380 :       return gerepileuptoleaf(av, z);
    2609             :     }
    2610             :   }
    2611             :   else
    2612             :   {
    2613             :     void *E;
    2614         135 :     const struct bb_group *S = get_FpXQ_star(&E,T,p);
    2615         135 :     GEN z = gen_PH_log(a,g,ord,E,S);
    2616         107 :     return gerepileuptoleaf(av, z);
    2617             :   }
    2618             : }
    2619             : 
    2620             : GEN
    2621     2193909 : Fq_log(GEN a, GEN g, GEN ord, GEN T, GEN p)
    2622             : {
    2623     2193909 :   if (!T) return Fp_log(a,g,ord,p);
    2624     1252103 :   if (typ(g) == t_INT)
    2625             :   {
    2626           0 :     if (typ(a) == t_POL)
    2627             :     {
    2628           0 :       if (degpol(a)) return cgetg(1,t_VEC);
    2629           0 :       a = gel(a,2);
    2630             :     }
    2631           0 :     return Fp_log(a,g,ord,p);
    2632             :   }
    2633     1252103 :   return typ(a) == t_INT? Fp_FpXQ_log(a,g,ord,T,p): FpXQ_log(a,g,ord,T,p);
    2634             : }
    2635             : 
    2636             : static GEN
    2637        2264 : FpXQ_sumautsum_sqr(void *E, GEN xzd)
    2638             : {
    2639        2264 :   struct _FpXQ *D = (struct _FpXQ*)E;
    2640        2264 :   pari_sp av = avma;
    2641             :   GEN xi, zeta, delta, xi2, zeta2, delta2, temp, xipow;
    2642        2264 :   GEN T = D->T, p = D-> p;
    2643             :   ulong d;
    2644        2264 :   xi = gel(xzd, 1); zeta = gel(xzd, 2); delta = gel(xzd, 3);
    2645             : 
    2646        2264 :   d = brent_kung_optpow(get_FpX_degree(T)-1,3,1);
    2647        2264 :   xipow = FpXQ_powers(xi, d, T, p);
    2648             : 
    2649        2264 :   xi2 = FpX_FpXQV_eval(xi, xipow, T, p);
    2650        2264 :   zeta2 = FpXQ_mul(zeta, FpX_FpXQV_eval(zeta,  xipow, T, p), T, p);
    2651        2264 :   temp  = FpXQ_mul(zeta, FpX_FpXQV_eval(delta, xipow, T, p), T, p);
    2652        2264 :   delta2 = FpX_add(delta, temp, p);
    2653        2264 :   return gerepilecopy(av, mkvec3(xi2, zeta2, delta2));
    2654             : }
    2655             : 
    2656             : static GEN
    2657        1123 : FpXQ_sumautsum_msqr(void *E, GEN xzd)
    2658             : {
    2659        1123 :   struct _FpXQ *D = (struct _FpXQ*)E;
    2660        1123 :   pari_sp av = avma;
    2661             :   GEN xii, zetai, deltai, xzd2;
    2662        1123 :   GEN T = D->T, p = D-> p, xi0pow = gel(D->aut, 1), zeta0 = gel(D->aut, 2);
    2663        1123 :   xzd2 = FpXQ_sumautsum_sqr(E, xzd);
    2664        1123 :   xii = FpX_FpXQV_eval(gel(xzd2, 1), xi0pow, T, p);
    2665        1123 :   zetai = FpXQ_mul(zeta0, FpX_FpXQV_eval(gel(xzd2, 2), xi0pow, T, p), T, p);
    2666        1123 :   deltai = FpX_add(gel(xzd2, 3), zetai, p);
    2667             : 
    2668        1123 :   return gerepilecopy(av, mkvec3(xii, zetai, deltai));
    2669             : }
    2670             : 
    2671             : /*returns a + a^(1+s) + a^(1+s+2s) + ... + a^(1+s+...+is)
    2672             :   where ax = [a,s] with s an automorphism */
    2673             : static GEN
    2674        1266 : FpXQ_sumautsum(GEN ax, long i, GEN T, GEN p) {
    2675        1266 :   pari_sp av = avma;
    2676             :   GEN a, xi, zeta, vec, res;
    2677             :   struct _FpXQ D;
    2678             :   ulong d;
    2679        1266 :   D.T = FpX_get_red(T, p); D.p = p;
    2680        1266 :   a = gel(ax, 1); xi = gel(ax,2);
    2681        1266 :   d = brent_kung_optpow(get_FpX_degree(T)-1,2*(hammingl(i)-1),1);
    2682        1266 :   zeta = FpX_FpXQ_eval(a, xi, T, p);
    2683        1266 :   D.aut = mkvec2(FpXQ_powers(xi, d, T, p), zeta);
    2684             : 
    2685        1266 :   vec = gen_powu_fold(mkvec3(xi, zeta, zeta), i, (void *)&D, FpXQ_sumautsum_sqr, FpXQ_sumautsum_msqr);
    2686        1266 :   res = FpXQ_mul(a, FpX_add(pol_1(get_FpX_var(T)), gel(vec, 3), p), T, p);
    2687             : 
    2688        1266 :   return gerepilecopy(av, res);
    2689             : }
    2690             : 
    2691             : /*algorithm from
    2692             : Doliskani, J., & Schost, E. (2014).
    2693             : Taking roots over high extensions of finite fields*/
    2694             : static GEN
    2695         543 : FpXQ_sqrtl_spec(GEN z, GEN n, GEN T, GEN p, GEN *zetan)
    2696             : {
    2697         543 :   pari_sp av = avma;
    2698             :   GEN psn, c, b, new_z, beta, x, y, w, ax, g, zeta;
    2699         543 :   long s, l, v = get_FpX_var(T), d = get_FpX_degree(T);
    2700         543 :   if(!isprime(n)) pari_err_PRIME("FpXQ_sqrtn", n);
    2701         543 :   s = itos(Fp_order(p, subiu(n,1), n));
    2702         543 :   if(s >= d || d % s != 0)
    2703           0 :     pari_err(e_MISC, "expected p's order mod n to divide the degree of T");
    2704         543 :   l = d/s;
    2705         543 :   if (!signe(z)) return pol_0(varn(z));
    2706         543 :   T = FpX_get_red(T, p);
    2707         543 :   ax = mkvec2(NULL, FpXQ_autpow(FpX_Frobenius(T,p), s, T, p));
    2708         543 :   psn = diviiexact(subii(powiu(p, s), gen_1), n);
    2709             :   do {
    2710         543 :     do c = random_FpX(d, v, p); while (!signe(c));
    2711         543 :     new_z = FpXQ_mul(z, FpXQ_pow(c, n, T, p), T, p);
    2712         543 :     gel(ax,1) = FpXQ_pow(new_z, psn, T, p);
    2713             : 
    2714             :     /*If l == 2, b has to be 1 + a^((p^s-1)/n)*/
    2715         543 :     if(l == 2) y = gel(ax, 1);
    2716         543 :     else y = FpXQ_sumautsum(ax, l-2, T, p);
    2717         543 :     b = FpX_Fp_add(y, gen_1, p);
    2718         543 :   } while (!signe(b));
    2719             : 
    2720         543 :   x = FpXQ_mul(new_z, FpXQ_pow(b, n, T, p), T, p);
    2721         543 :   if(s == 1) {
    2722         221 :     if (degpol(x) > 0) return gc_NULL(av);
    2723         186 :     beta = Fp_sqrtn(constant_coeff(x), n, p, &zeta);
    2724         186 :     if (!beta) return gc_NULL(av);
    2725         186 :     if(zetan) *zetan = scalarpol(zeta, varn(z));
    2726         186 :     w = FpX_Fp_mul(FpXQ_inv(FpXQ_mul(b, c, T, p), T, p), beta, p);
    2727         186 :     gerepileall(av, zetan? 2: 1, &w, zetan);
    2728         186 :     return w;
    2729             :   }
    2730         322 :   g = FpXQ_minpoly(x, T, p);
    2731         322 :   if (degpol(g) != s) return gc_NULL(av);
    2732             : 
    2733         322 :   beta = FpXQ_sqrtn(pol_x(varn(z)), n, g, p, &zeta);
    2734         322 :   if (!beta) return gc_NULL(av);
    2735             : 
    2736         322 :   if(zetan) *zetan = FpX_FpXQ_eval(zeta, x, T, p);
    2737         322 :   beta = FpX_FpXQ_eval(beta, x, T, p);
    2738         322 :   w = FpXQ_mul(FpXQ_inv(FpXQ_mul(b, c, T, p), T, p), beta, T, p);
    2739         322 :   gerepileall(av, zetan? 2: 1, &w, zetan);
    2740         322 :   return w;
    2741             : }
    2742             : 
    2743             : static GEN
    2744         858 : FpXQ_sqrtn_spec(GEN a, GEN n, GEN T, GEN p, GEN q, GEN *zetan)
    2745             : {
    2746         858 :   pari_sp ltop = avma;
    2747             :   GEN z, m, u1, u2;
    2748             :   int is_1;
    2749         858 :   if (is_pm1(n))
    2750             :   {
    2751         581 :     if (zetan) *zetan = pol_1(varn(a));
    2752         581 :     return signe(n) < 0? FpXQ_inv(a, T, p): gcopy(a);
    2753             :   }
    2754         277 :   is_1 = gequal1(a);
    2755         277 :   if (is_1 && !zetan) return gcopy(a);
    2756         270 :   z = pol_1(varn(a));
    2757         270 :   m = bezout(n,q,&u1,&u2);
    2758         270 :   if (!is_pm1(m))
    2759             :   {
    2760         270 :     GEN F = Z_factor(m);
    2761         270 :     long i, j, j2 = 0;
    2762             :     GEN y, l;
    2763         270 :     pari_sp av1 = avma;
    2764         505 :     for (i = nbrows(F); i; i--)
    2765             :     {
    2766         270 :       l = gcoeff(F,i,1);
    2767         270 :       j = itos(gcoeff(F,i,2));
    2768         270 :       if(zetan) {
    2769          71 :         a = FpXQ_sqrtl_spec(a,l,T,p,&y);
    2770         106 :         if (!a) return gc_NULL(ltop);
    2771          71 :         j--;
    2772          71 :         j2 = j;
    2773             :       }
    2774         270 :       if (!is_1 && j > 0) {
    2775             :         do
    2776             :         {
    2777         395 :           a = FpXQ_sqrtl_spec(a,l,T,p,NULL);
    2778         395 :           if (!a) return gc_NULL(ltop);
    2779         360 :         } while (--j);
    2780             :       }
    2781             :       /*This is below finding a's root,
    2782             :       so we don't spend time doing this, if a is not n-th root*/
    2783         235 :       if(zetan) {
    2784         148 :         for(; j2>0; j2--) y = FpXQ_sqrtl_spec(y, l, T, p, NULL);
    2785          71 :         z = FpXQ_mul(z, y, T, p);
    2786             :       }
    2787         235 :       if (gc_needed(ltop,1))
    2788             :       { /* n can have lots of prime factors*/
    2789           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"FpXQ_sqrtn_spec");
    2790           0 :         gerepileall(av1, zetan? 2: 1, &a, &z);
    2791             :       }
    2792             :     }
    2793             :   }
    2794             : 
    2795         235 :   if (!equalii(m, n))
    2796          42 :     a = FpXQ_pow(a,modii(u1,q), T, p);
    2797         235 :   if (zetan)
    2798             :   {
    2799          71 :     *zetan = z;
    2800          71 :     gerepileall(ltop,2,&a,zetan);
    2801             :   }
    2802             :   else /* is_1 is 0: a was modified above -> gerepileupto valid */
    2803         164 :     a = gerepileupto(ltop, a);
    2804         235 :   return a;
    2805             : }
    2806             : GEN
    2807        1222 : FpXQ_sqrtn(GEN a, GEN n, GEN T, GEN p, GEN *zeta)
    2808             : {
    2809        1222 :   pari_sp av = avma;
    2810             :   GEN z;
    2811        1222 :   if (!signe(a))
    2812             :   {
    2813         154 :     long v=varn(a);
    2814         154 :     if (signe(n) < 0) pari_err_INV("FpXQ_sqrtn",a);
    2815         147 :     if (zeta) *zeta=pol_1(v);
    2816         147 :     return pol_0(v);
    2817             :   }
    2818        1068 :   if (lgefint(p)==3)
    2819             :   {
    2820         210 :     if (uel(p,2) == 2)
    2821             :     {
    2822          14 :       z = F2xq_sqrtn(ZX_to_F2x(a), n, ZX_to_F2x(get_FpX_mod(T)), zeta);
    2823          14 :       if (!z) return NULL;
    2824          14 :       z = F2x_to_ZX(z);
    2825          14 :       if (!zeta) return gerepileuptoleaf(av, z);
    2826           7 :       *zeta=F2x_to_ZX(*zeta);
    2827             :     } else
    2828             :     {
    2829         196 :       ulong pp = to_Flxq(&a, &T, p);
    2830         196 :       z = Flxq_sqrtn(a, n, T, pp, zeta);
    2831         196 :       if (!z) return NULL;
    2832         189 :       if (!zeta) return Flx_to_ZX_inplace(gerepileuptoleaf(av, z));
    2833          63 :       z = Flx_to_ZX(z);
    2834          63 :       *zeta=Flx_to_ZX(*zeta);
    2835             :     }
    2836             :   }
    2837             :   else
    2838             :   {
    2839             :     void *E;
    2840         858 :     const struct bb_group *S = get_FpXQ_star(&E,T,p);
    2841         858 :     GEN o = subiu(powiu(p,get_FpX_degree(T)),1);
    2842             :     GEN m, u1, u2, l, zeta2, F, n2;
    2843         858 :     long i, s, d = degpol(T);
    2844             : 
    2845         858 :     m = bezout(n,o,&u1,&u2);
    2846         858 :     F = Z_factor(m);
    2847        1632 :     for (i = nbrows(F); i; i--)
    2848             :     {
    2849         774 :       l = gcoeff(F,i,1);
    2850         774 :       s = itos(Fp_order(p, subiu(l, 1), l));
    2851             :       /*FpXQ_sqrtn_spec only works if d > s and s | d
    2852             :       for those factors of m we use FpXQ_sqrtn_spec
    2853             :       for the other factor we stay with gen_Shanks_sqrtn*/
    2854         774 :       if(d <= s || d % s != 0) {
    2855         497 :         gcoeff(F,i,2) = gen_0;
    2856             :       }
    2857         277 :       else gcoeff(F,i,2) = stoi(Z_pval(n,l));
    2858             :     }
    2859         858 :     F = factorback(F);
    2860         858 :     z = FpXQ_sqrtn_spec(a,F,T, p, o,zeta);
    2861        1260 :     if(!z) return gc_NULL(av);
    2862         823 :     n2 = diviiexact(n, F);
    2863         823 :     if(!gequal1(n2)) {
    2864         602 :       if(zeta) zeta2 = gcopy(*zeta);
    2865         602 :       z = gen_Shanks_sqrtn(z, n2, o, zeta, E, S);
    2866         602 :       if (!z) return gc_NULL(av);
    2867         602 :       if(zeta) *zeta = FpXQ_mul(*zeta, zeta2, T, p);
    2868             :     }
    2869         823 :     if (!zeta) return gerepileupto(av, z);
    2870             :   }
    2871         491 :   return gc_all(av, 2, &z,zeta);
    2872             : }
    2873             : 
    2874             : static GEN
    2875       19390 : Fp2_norm(GEN x, GEN D, GEN p)
    2876             : {
    2877       19390 :   GEN a = gel(x,1), b = gel(x,2);
    2878       19390 :   if (signe(b)==0) return Fp_sqr(a,p);
    2879       19390 :   return Fp_sub(sqri(a), mulii(D, Fp_sqr(b, p)), p);
    2880             : }
    2881             : 
    2882             : static GEN
    2883       19821 : Fp2_sqrt(GEN z, GEN D, GEN p)
    2884             : {
    2885       19821 :   GEN a = gel(z,1), b = gel(z,2), as2, u, v, s;
    2886       19821 :   GEN y = Fp_2gener_i(D, p);
    2887       19821 :   if (signe(b)==0)
    2888         431 :     return kronecker(a, p)==1 ? mkvec2(Fp_sqrt_i(a, y, p), gen_0)
    2889         431 :                               : mkvec2(gen_0,Fp_sqrt_i(Fp_div(a, D, p), y, p));
    2890       19390 :   s = Fp_sqrt_i(Fp2_norm(z, D, p), y, p);
    2891       19390 :   if(!s) return NULL;
    2892       18980 :   as2 = Fp_halve(Fp_add(a, s, p), p);
    2893       18980 :   if (kronecker(as2, p)==-1) as2 = Fp_sub(as2,s,p);
    2894       18980 :   u = Fp_sqrt_i(as2, y, p);
    2895       18980 :   v = Fp_div(b, Fp_double(u, p), p);
    2896       18980 :   return mkvec2(u,v);
    2897             : }
    2898             : 
    2899             : 
    2900             : GEN
    2901       80860 : FpXQ_sqrt(GEN z, GEN T, GEN p)
    2902             : {
    2903       80860 :   pari_sp av = avma;
    2904             :   long d;
    2905       80860 :   if (lgefint(p)==3)
    2906             :   {
    2907       60308 :     if (uel(p,2) == 2)
    2908             :     {
    2909        5362 :       GEN r = F2xq_sqrt(ZX_to_F2x(z), ZX_to_F2x(get_FpX_mod(T)));
    2910        5362 :       return gerepileupto(av, F2x_to_ZX(r));
    2911             :     } else
    2912             :     {
    2913       54946 :       ulong pp = to_Flxq(&z, &T, p);
    2914       54946 :       z = Flxq_sqrt(z, T, pp);
    2915       54946 :       if (!z) return NULL;
    2916       52178 :       return gerepileupto(av, Flx_to_ZX(z));
    2917             :     }
    2918             :   }
    2919       20552 :   d = get_FpX_degree(T);
    2920       20552 :   if (d==2)
    2921             :   {
    2922       19821 :     GEN P = get_FpX_mod(T);
    2923       19821 :     GEN c = gel(P,2), b = gel(P,3), a = gel(P,4), b2 = Fp_halve(b, p);
    2924       19821 :     GEN t = Fp_div(b2, a, p);
    2925       19821 :     GEN D = Fp_sub(Fp_sqr(b2, p), Fp_mul(a, c, p), p);
    2926       20252 :     GEN x = degpol(z)<1 ? constant_coeff(z)
    2927       19821 :                         : Fp_sub(gel(z,2), Fp_mul(gel(z,3), t, p), p);
    2928       19821 :     GEN y = degpol(z)<1 ? gen_0: gel(z,3);
    2929       19821 :     GEN r = Fp2_sqrt(mkvec2(x, y), D, p), s;
    2930       19821 :     if (!r) return gc_NULL(av);
    2931       19411 :     s = deg1pol_shallow(gel(r,2),Fp_add(gel(r,1), Fp_mul(gel(r,2),t,p), p), varn(P));
    2932       19411 :     return gerepilecopy(av, s);
    2933             :   }
    2934         731 :   if (lgpol(z)<=1 && odd(d))
    2935             :   {
    2936           8 :     pari_sp av = avma;
    2937           8 :     GEN s = Fp_sqrt(constant_coeff(z), p);
    2938           8 :     if (!s) return gc_NULL(av);
    2939           8 :     return gerepilecopy(av, scalarpol_shallow(s, get_FpX_var(T)));
    2940             :   } else
    2941             :   {
    2942             :     GEN p2, c, b, new_z, beta, x, y, w, ax;
    2943         723 :     long v = get_FpX_var(T);
    2944         723 :     if (!signe(z)) return pol_0(varn(z));
    2945         723 :     T = FpX_get_red(T, p);
    2946         723 :     ax = mkvec2(NULL, FpX_Frobenius(T,p));
    2947         723 :     p2 = shifti(p, -1); /* (p - 1) / 2 */
    2948             :     do {
    2949         723 :       do c = random_FpX(d, v, p); while (!signe(c));
    2950         723 :       new_z = FpXQ_mul(z, FpXQ_sqr(c, T, p), T, p);
    2951         723 :       gel(ax,1) = FpXQ_pow(new_z, p2, T, p);
    2952         723 :       y = FpXQ_sumautsum(ax, d-2, T, p); /* d > 2 */
    2953         723 :       b = FpX_Fp_add(y, gen_1, p);
    2954         723 :     } while (!signe(b));
    2955             : 
    2956         723 :     x = FpXQ_mul(new_z, FpXQ_sqr(b, T, p), T, p);
    2957         723 :     if (degpol(x) > 0) return gc_NULL(av);
    2958         709 :     beta = Fp_sqrt(constant_coeff(x), p);
    2959         709 :     if (!beta) return gc_NULL(av);
    2960         709 :     w = FpX_Fp_mul(FpXQ_inv(FpXQ_mul(b, c, T, p), T, p), beta, p);
    2961         709 :     return gerepilecopy(av, w);
    2962             :   }
    2963             : }
    2964             : 
    2965             : GEN
    2966      363843 : FpXQ_norm(GEN x, GEN TB, GEN p)
    2967             : {
    2968      363843 :   pari_sp av = avma;
    2969      363843 :   GEN T = get_FpX_mod(TB);
    2970      363843 :   GEN y = FpX_resultant(T, x, p);
    2971      363843 :   GEN L = leading_coeff(T);
    2972      363843 :   if (gequal1(L) || signe(x)==0) return y;
    2973           0 :   return gerepileupto(av, Fp_div(y, Fp_pows(L, degpol(x), p), p));
    2974             : }
    2975             : 
    2976             : GEN
    2977       21063 : FpXQ_trace(GEN x, GEN TB, GEN p)
    2978             : {
    2979       21063 :   pari_sp av = avma;
    2980       21063 :   GEN T = get_FpX_mod(TB);
    2981       21063 :   GEN dT = FpX_deriv(T,p);
    2982       21063 :   long n = degpol(dT);
    2983       21063 :   GEN z = FpXQ_mul(x, dT, TB, p);
    2984       21063 :   if (degpol(z)<n) return gc_const(av, gen_0);
    2985       19880 :   return gerepileuptoint(av, Fp_div(gel(z,2+n), gel(T,3+n),p));
    2986             : }
    2987             : 
    2988             : GEN
    2989          15 : FpXQ_charpoly(GEN x, GEN T, GEN p)
    2990             : {
    2991          15 :   pari_sp ltop=avma;
    2992          15 :   long vT, v = fetch_var();
    2993             :   GEN R;
    2994          15 :   T = leafcopy(get_FpX_mod(T));
    2995          15 :   vT = varn(T); setvarn(T, v);
    2996          15 :   x = leafcopy(x); setvarn(x, v);
    2997          15 :   R = FpX_FpXY_resultant(T, deg1pol_shallow(gen_1,FpX_neg(x,p),vT),p);
    2998          15 :   (void)delete_var(); return gerepileupto(ltop,R);
    2999             : }
    3000             : 
    3001             : /* Computing minimal polynomial :                         */
    3002             : /* cf Shoup 'Efficient Computation of Minimal Polynomials */
    3003             : /*          in Algebraic Extensions of Finite Fields'     */
    3004             : 
    3005             : /* Let v a linear form, return the linear form z->v(tau*z)
    3006             :    that is, v*(M_tau) */
    3007             : 
    3008             : static GEN
    3009        1646 : FpXQ_transmul_init(GEN tau, GEN T, GEN p)
    3010             : {
    3011             :   GEN bht;
    3012        1646 :   GEN h, Tp = get_FpX_red(T, &h);
    3013        1646 :   long n = degpol(Tp), vT = varn(Tp);
    3014        1646 :   GEN ft = FpX_recipspec(Tp+2, n+1, n+1);
    3015        1646 :   GEN bt = FpX_recipspec(tau+2, lgpol(tau), n);
    3016        1646 :   setvarn(ft, vT); setvarn(bt, vT);
    3017        1646 :   if (h)
    3018         602 :     bht = FpXn_mul(bt, h, n-1, p);
    3019             :   else
    3020             :   {
    3021        1044 :     GEN bh = FpX_div(FpX_shift(tau, n-1), T, p);
    3022        1044 :     bht = FpX_recipspec(bh+2, lgpol(bh), n-1);
    3023        1044 :     setvarn(bht, vT);
    3024             :   }
    3025        1646 :   return mkvec3(bt, bht, ft);
    3026             : }
    3027             : 
    3028             : static GEN
    3029        7141 : FpXQ_transmul(GEN tau, GEN a, long n, GEN p)
    3030             : {
    3031        7141 :   pari_sp ltop = avma;
    3032             :   GEN t1, t2, t3, vec;
    3033        7141 :   GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
    3034        7141 :   if (signe(a)==0) return pol_0(varn(a));
    3035        7106 :   t2 = FpX_shift(FpX_mul(bt, a, p),1-n);
    3036        7106 :   if (signe(bht)==0) return gerepilecopy(ltop, t2);
    3037        6262 :   t1 = FpX_shift(FpX_mul(ft, a, p),-n);
    3038        6262 :   t3 = FpXn_mul(t1, bht, n-1, p);
    3039        6262 :   vec = FpX_sub(t2, FpX_shift(t3, 1), p);
    3040        6262 :   return gerepileupto(ltop, vec);
    3041             : }
    3042             : 
    3043             : GEN
    3044       13701 : FpXQ_minpoly(GEN x, GEN T, GEN p)
    3045             : {
    3046       13701 :   pari_sp ltop = avma;
    3047             :   long vT, n;
    3048             :   GEN v_x, g, tau;
    3049       13701 :   if (lgefint(p)==3)
    3050             :   {
    3051       12878 :     ulong pp = to_Flxq(&x, &T, p);
    3052       12878 :     GEN g = Flxq_minpoly(x, T, pp);
    3053       12878 :     return gerepileupto(ltop, Flx_to_ZX(g));
    3054             :   }
    3055         823 :   vT = get_FpX_var(T);
    3056         823 :   n = get_FpX_degree(T);
    3057         823 :   g = pol_1(vT);
    3058         823 :   tau = pol_1(vT);
    3059         823 :   T = FpX_get_red(T, p);
    3060         823 :   x = FpXQ_red(x, T, p);
    3061         823 :   v_x = FpXQ_powers(x, usqrt(2*n), T, p);
    3062        1646 :   while(signe(tau) != 0)
    3063             :   {
    3064             :     long i, j, m, k1;
    3065             :     GEN M, v, tr;
    3066             :     GEN g_prime, c;
    3067         823 :     if (degpol(g) == n) { tau = pol_1(vT); g = pol_1(vT); }
    3068         823 :     v = random_FpX(n, vT, p);
    3069         823 :     tr = FpXQ_transmul_init(tau, T, p);
    3070         823 :     v = FpXQ_transmul(tr, v, n, p);
    3071         823 :     m = 2*(n-degpol(g));
    3072         823 :     k1 = usqrt(m);
    3073         823 :     tr = FpXQ_transmul_init(gel(v_x,k1+1), T, p);
    3074         823 :     c = cgetg(m+2,t_POL);
    3075         823 :     c[1] = evalsigne(1)|evalvarn(vT);
    3076        7141 :     for (i=0; i<m; i+=k1)
    3077             :     {
    3078        6318 :       long mj = minss(m-i, k1);
    3079       67450 :       for (j=0; j<mj; j++)
    3080       61132 :         gel(c,m+1-(i+j)) = FpX_dotproduct(v, gel(v_x,j+1), p);
    3081        6318 :       v = FpXQ_transmul(tr, v, n, p);
    3082             :     }
    3083         823 :     c = FpX_renormalize(c, m+2);
    3084             :     /* now c contains <v,x^i> , i = 0..m-1  */
    3085         823 :     M = FpX_halfgcd(pol_xn(m, vT), c, p);
    3086         823 :     g_prime = gmael(M, 2, 2);
    3087         823 :     if (degpol(g_prime) < 1) continue;
    3088         823 :     g = FpX_mul(g, g_prime, p);
    3089         823 :     tau = FpXQ_mul(tau, FpX_FpXQV_eval(g_prime, v_x, T, p), T, p);
    3090             :   }
    3091         823 :   g = FpX_normalize(g,p);
    3092         823 :   return gerepilecopy(ltop,g);
    3093             : }
    3094             : 
    3095             : GEN
    3096           8 : FpXQ_conjvec(GEN x, GEN T, GEN p)
    3097             : {
    3098           8 :   pari_sp av=avma;
    3099             :   long i;
    3100           8 :   long n = get_FpX_degree(T), v = varn(x);
    3101           8 :   GEN M = FpX_matFrobenius(T, p);
    3102           8 :   GEN z = cgetg(n+1,t_COL);
    3103           8 :   gel(z,1) = RgX_to_RgC(x,n);
    3104          17 :   for (i=2; i<=n; i++) gel(z,i) = FpM_FpC_mul(M,gel(z,i-1),p);
    3105           8 :   gel(z,1) = x;
    3106          17 :   for (i=2; i<=n; i++) gel(z,i) = RgV_to_RgX(gel(z,i),v);
    3107           8 :   return gerepilecopy(av,z);
    3108             : }
    3109             : 
    3110             : /* p prime, p_1 = p-1, q = p^deg T, Lp = cofactors of some prime divisors
    3111             :  * l_p of p-1, Lq = cofactors of some prime divisors l_q of q-1, return a
    3112             :  * g in Fq such that
    3113             :  * - Ng generates all l_p-Sylows of Fp^*
    3114             :  * - g generates all l_q-Sylows of Fq^* */
    3115             : static GEN
    3116       83398 : gener_FpXQ_i(GEN T, GEN p, GEN p_1, GEN Lp, GEN Lq)
    3117             : {
    3118             :   pari_sp av;
    3119       83398 :   long vT = varn(T), f = degpol(T), l = lg(Lq);
    3120       83398 :   GEN F = FpX_Frobenius(T, p);
    3121       83400 :   int p_is_2 = is_pm1(p_1);
    3122      167381 :   for (av = avma;; set_avma(av))
    3123       83981 :   {
    3124      167381 :     GEN t, g = random_FpX(f, vT, p);
    3125             :     long i;
    3126      167381 :     if (degpol(g) < 1) continue;
    3127      108167 :     if (p_is_2)
    3128       55688 :       t = g;
    3129             :     else
    3130             :     {
    3131       52479 :       t = FpX_resultant(T, g, p); /* Ng = g^((q-1)/(p-1)), assuming T monic */
    3132       52479 :       if (kronecker(t, p) == 1) continue;
    3133       31306 :       if (lg(Lp) > 1 && !is_gener_Fp(t, p, p_1, Lp)) continue;
    3134       30101 :       t = FpXQ_pow(g, shifti(p_1,-1), T, p);
    3135             :     }
    3136       98614 :     for (i = 1; i < l; i++)
    3137             :     {
    3138       15214 :       GEN a = FpXQ_pow_Frobenius(t, gel(Lq,i), F, T, p);
    3139       15214 :       if (!degpol(a) && equalii(gel(a,2), p_1)) break;
    3140             :     }
    3141       85789 :     if (i == l) return g;
    3142             :   }
    3143             : }
    3144             : 
    3145             : GEN
    3146        7030 : gener_FpXQ(GEN T, GEN p, GEN *po)
    3147             : {
    3148        7030 :   long i, j, f = get_FpX_degree(T);
    3149             :   GEN g, Lp, Lq, p_1, q_1, N, o;
    3150        7030 :   pari_sp av = avma;
    3151             : 
    3152        7030 :   p_1 = subiu(p,1);
    3153        7030 :   if (f == 1) {
    3154             :     GEN Lp, fa;
    3155           7 :     o = p_1;
    3156           7 :     fa = Z_factor(o);
    3157           7 :     Lp = gel(fa,1);
    3158           7 :     Lp = vecslice(Lp, 2, lg(Lp)-1); /* remove 2 for efficiency */
    3159             : 
    3160           7 :     g = cgetg(3, t_POL);
    3161           7 :     g[1] = evalsigne(1) | evalvarn(get_FpX_var(T));
    3162           7 :     gel(g,2) = pgener_Fp_local(p, Lp);
    3163           7 :     if (po) *po = mkvec2(o, fa);
    3164           7 :     return g;
    3165             :   }
    3166        7023 :   if (lgefint(p) == 3)
    3167             :   {
    3168        6986 :     ulong pp = to_Flxq(NULL, &T, p);
    3169        6986 :     g = gener_Flxq(T, pp, po);
    3170        6986 :     if (!po) return Flx_to_ZX_inplace(gerepileuptoleaf(av, g));
    3171        6986 :     g = Flx_to_ZX(g); return gc_all(av, 2, &g, po);
    3172             :   }
    3173             :   /* p now odd */
    3174          37 :   q_1 = subiu(powiu(p,f), 1);
    3175          37 :   N = diviiexact(q_1, p_1);
    3176          37 :   Lp = odd_prime_divisors(p_1);
    3177         168 :   for (i=lg(Lp)-1; i; i--) gel(Lp,i) = diviiexact(p_1, gel(Lp,i));
    3178          37 :   o = factor_pn_1(p,f);
    3179          37 :   Lq = leafcopy( gel(o, 1) );
    3180         353 :   for (i = j = 1; i < lg(Lq); i++)
    3181             :   {
    3182         316 :     if (dvdii(p_1, gel(Lq,i))) continue;
    3183         148 :     gel(Lq,j++) = diviiexact(N, gel(Lq,i));
    3184             :   }
    3185          37 :   setlg(Lq, j);
    3186          37 :   g = gener_FpXQ_i(get_FpX_mod(T), p, p_1, Lp, Lq);
    3187          37 :   if (!po) g = gerepilecopy(av, g);
    3188             :   else {
    3189          21 :     *po = mkvec2(q_1, o);
    3190          21 :     gerepileall(av, 2, &g, po);
    3191             :   }
    3192          37 :   return g;
    3193             : }
    3194             : 
    3195             : GEN
    3196       83362 : gener_FpXQ_local(GEN T, GEN p, GEN L)
    3197             : {
    3198             :   GEN Lp, Lq, p_1, q_1, N, Q;
    3199             :   long f, i, ip, iq, l;
    3200             : 
    3201       83362 :   T = get_FpX_mod(T);
    3202       83362 :   f = degpol(T);
    3203       83362 :   if (f == 1) return pgener_Fp_local(p, L);
    3204       83362 :   l = lg(L);
    3205       83362 :   p_1 = subiu(p,1);
    3206       83360 :   q_1 = subiu(powiu(p,f), 1);
    3207       83361 :   N = diviiexact(q_1, p_1);
    3208             : 
    3209       83360 :   Q = is_pm1(p_1)? gen_1: shifti(p_1,-1);
    3210       83361 :   Lp = cgetg(l, t_VEC); ip = 1;
    3211       83360 :   Lq = cgetg(l, t_VEC); iq = 1;
    3212       98873 :   for (i=1; i < l; i++)
    3213             :   {
    3214       15512 :     GEN a, b, ell = gel(L,i);
    3215       15512 :     if (absequaliu(ell,2)) continue;
    3216       15232 :     a = dvmdii(Q, ell, &b);
    3217       15232 :     if (b == gen_0)
    3218        2555 :       gel(Lp,ip++) = a;
    3219             :     else
    3220       12677 :       gel(Lq,iq++) = diviiexact(N,ell);
    3221             :   }
    3222       83361 :   setlg(Lp, ip);
    3223       83361 :   setlg(Lq, iq);
    3224       83361 :   return gener_FpXQ_i(T, p, p_1, Lp, Lq);
    3225             : }
    3226             : 
    3227             : /***********************************************************************/
    3228             : /**                                                                   **/
    3229             : /**                              FpXn                                 **/
    3230             : /**                                                                   **/
    3231             : /***********************************************************************/
    3232             : 
    3233             : GEN
    3234     2567038 : FpXn_mul(GEN a, GEN b, long n, GEN p)
    3235             : {
    3236     2567038 :   return FpX_red(ZXn_mul(a, b, n), p);
    3237             : }
    3238             : 
    3239             : GEN
    3240           0 : FpXn_sqr(GEN a, long n, GEN p)
    3241             : {
    3242           0 :   return FpX_red(ZXn_sqr(a, n), p);
    3243             : }
    3244             : 
    3245             : /* (f*g) \/ x^n */
    3246             : static GEN
    3247      116133 : FpX_mulhigh_i(GEN f, GEN g, long n, GEN p)
    3248             : {
    3249      116133 :   return FpX_shift(FpX_mul(f,g, p),-n);
    3250             : }
    3251             : 
    3252             : static GEN
    3253       59970 : FpXn_mulhigh(GEN f, GEN g, long n2, long n, GEN p)
    3254             : {
    3255       59970 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    3256       59970 :   return FpX_add(FpX_mulhigh_i(fl, g, n2, p), FpXn_mul(fh, g, n - n2, p), p);
    3257             : }
    3258             : 
    3259             : GEN
    3260        6412 : FpXn_div(GEN g, GEN f, long e, GEN p)
    3261             : {
    3262        6412 :   pari_sp av = avma, av2;
    3263             :   ulong mask;
    3264             :   GEN W, a;
    3265        6412 :   long v = varn(f), n = 1;
    3266             : 
    3267        6412 :   if (!signe(f)) pari_err_INV("FpXn_inv",f);
    3268        6412 :   a = Fp_inv(gel(f,2), p);
    3269        6412 :   if (e == 1 && !g) return scalarpol(a, v);
    3270        6412 :   else if (e == 2 && !g)
    3271             :   {
    3272             :     GEN b;
    3273           0 :     if (degpol(f) <= 0) return scalarpol(a, v);
    3274           0 :     b = Fp_neg(gel(f,3),p);
    3275           0 :     if (signe(b)==0) return scalarpol(a, v);
    3276           0 :     if (!is_pm1(a)) b = Fp_mul(b, Fp_sqr(a, p), p);
    3277           0 :     W = deg1pol_shallow(b, a, v);
    3278           0 :     return gerepilecopy(av, W);
    3279             :   }
    3280        6412 :   W = scalarpol_shallow(Fp_inv(gel(f,2), p),v);
    3281        6412 :   mask = quadratic_prec_mask(e);
    3282        6412 :   av2 = avma;
    3283       27580 :   for (;mask>1;)
    3284             :   {
    3285             :     GEN u, fr;
    3286       21168 :     long n2 = n;
    3287       21168 :     n<<=1; if (mask & 1) n--;
    3288       21168 :     mask >>= 1;
    3289       21168 :     fr = FpXn_red(f, n);
    3290       21168 :     if (mask>1 || !g)
    3291             :     {
    3292       21168 :       u = FpXn_mul(W, FpXn_mulhigh(fr, W, n2, n, p), n-n2, p);
    3293       21168 :       W = FpX_sub(W, FpX_shift(u, n2), p);
    3294             :     }
    3295             :     else
    3296             :     {
    3297           0 :       GEN y = FpXn_mul(g, W, n, p), yt =  FpXn_red(y, n-n2);
    3298           0 :       u = FpXn_mul(yt, FpXn_mulhigh(fr,  W, n2, n, p), n-n2, p);
    3299           0 :       W = FpX_sub(y, FpX_shift(u, n2), p);
    3300             :     }
    3301       21168 :     if (gc_needed(av2,2))
    3302             :     {
    3303           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpXn_inv, e = %ld", n);
    3304           0 :       W = gerepileupto(av2, W);
    3305             :     }
    3306             :   }
    3307        6412 :   return gerepileupto(av, W);
    3308             : }
    3309             : 
    3310             : GEN
    3311        6412 : FpXn_inv(GEN f, long e, GEN p)
    3312        6412 : { return FpXn_div(NULL, f, e, p); }
    3313             : 
    3314             : GEN
    3315       17361 : FpXn_expint(GEN h, long e, GEN p)
    3316             : {
    3317       17361 :   pari_sp av = avma, av2;
    3318       17361 :   long v = varn(h), n=1;
    3319       17361 :   GEN f = pol_1(v), g = pol_1(v);
    3320       17361 :   ulong mask = quadratic_prec_mask(e);
    3321       17361 :   av2 = avma;
    3322       56163 :   for (;mask>1;)
    3323             :   {
    3324             :     GEN u, w;
    3325       56163 :     long n2 = n;
    3326       56163 :     n<<=1; if (mask & 1) n--;
    3327       56163 :     mask >>= 1;
    3328       56163 :     u = FpXn_mul(g, FpX_mulhigh_i(f, FpXn_red(h, n2-1), n2-1, p), n-n2, p);
    3329       56163 :     u = FpX_add(u, FpX_shift(FpXn_red(h, n-1), 1-n2), p);
    3330       56163 :     w = FpXn_mul(f, FpX_integXn(u, n2-1, p), n-n2, p);
    3331       56163 :     f = FpX_add(f, FpX_shift(w, n2), p);
    3332       56163 :     if (mask<=1) break;
    3333       38802 :     u = FpXn_mul(g, FpXn_mulhigh(f, g, n2, n, p), n-n2, p);
    3334       38802 :     g = FpX_sub(g, FpX_shift(u, n2), p);
    3335       38802 :     if (gc_needed(av2,2))
    3336             :     {
    3337           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpXn_exp, e = %ld", n);
    3338           0 :       gerepileall(av2, 2, &f, &g);
    3339             :     }
    3340             :   }
    3341       17361 :   return gerepileupto(av, f);
    3342             : }
    3343             : 
    3344             : GEN
    3345           0 : FpXn_exp(GEN h, long e, GEN p)
    3346             : {
    3347           0 :   if (signe(h)==0 || degpol(h)<1 || !gequal0(gel(h,2)))
    3348           0 :     pari_err_DOMAIN("FpXn_exp","valuation", "<", gen_1, h);
    3349           0 :   return FpXn_expint(FpX_deriv(h, p), e, p);
    3350             : }

Generated by: LCOV version 1.16