Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - polarit3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21196-f12677d) Lines: 1274 1490 85.5 %
Date: 2017-10-22 06:23:24 Functions: 133 147 90.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000-2005  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /***********************************************************************/
      15             : /**                                                                   **/
      16             : /**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
      17             : /**                         (third part)                              **/
      18             : /**                                                                   **/
      19             : /***********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /************************************************************************
      24             :  **                                                                    **
      25             :  **                      Ring membership                               **
      26             :  **                                                                    **
      27             :  ************************************************************************/
      28             : struct charact {
      29             :   GEN q;
      30             :   int isprime;
      31             : };
      32             : static void
      33         637 : char_update_prime(struct charact *S, GEN p)
      34             : {
      35         637 :   if (!S->isprime) { S->isprime = 1; S->q = p; }
      36         637 :   if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
      37         630 : }
      38             : static void
      39        1365 : char_update_int(struct charact *S, GEN n)
      40             : {
      41        1365 :   if (S->isprime)
      42             :   {
      43        1365 :     if (dvdii(n, S->q)) return;
      44           7 :     pari_err_MODULUS("characteristic", S->q, n);
      45             :   }
      46        1358 :   S->q = gcdii(S->q, n);
      47             : }
      48             : static void
      49      578340 : charact(struct charact *S, GEN x)
      50             : {
      51      578340 :   const long tx = typ(x);
      52             :   long i, l;
      53      578340 :   switch(tx)
      54             :   {
      55         777 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      56         546 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      57             :     case t_COMPLEX: case t_QUAD:
      58             :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      59             :     case t_VEC: case t_COL: case t_MAT:
      60       10941 :       l = lg(x);
      61       10941 :       for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
      62       10927 :       break;
      63             :     case t_LIST:
      64           7 :       x = list_data(x);
      65           7 :       if (x) charact(S, x);
      66           7 :       break;
      67             :   }
      68      578312 : }
      69             : static void
      70       32340 : charact_res(struct charact *S, GEN x)
      71             : {
      72       32340 :   const long tx = typ(x);
      73             :   long i, l;
      74       32340 :   switch(tx)
      75             :   {
      76         588 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      77           0 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      78          91 :     case t_PADIC: char_update_prime(S, gel(x,2)); break;
      79             :     case t_COMPLEX: case t_QUAD:
      80             :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      81             :     case t_VEC: case t_COL: case t_MAT:
      82        9919 :       l = lg(x);
      83        9919 :       for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
      84        9919 :       break;
      85             :     case t_LIST:
      86           0 :       x = list_data(x);
      87           0 :       if (x) charact_res(S, x);
      88           0 :       break;
      89             :   }
      90       32340 : }
      91             : GEN
      92        8449 : characteristic(GEN x)
      93             : {
      94             :   struct charact S;
      95        8449 :   S.q = gen_0; S.isprime = 0;
      96        8449 :   charact(&S, x); return S.q;
      97             : }
      98             : GEN
      99        2415 : residual_characteristic(GEN x)
     100             : {
     101             :   struct charact S;
     102        2415 :   S.q = gen_0; S.isprime = 0;
     103        2415 :   charact_res(&S, x); return S.q;
     104             : }
     105             : 
     106             : int
     107    69590439 : Rg_is_Fp(GEN x, GEN *pp)
     108             : {
     109             :   GEN mod;
     110    69590439 :   switch(typ(x))
     111             :   {
     112             :   case t_INTMOD:
     113     6066823 :     mod = gel(x,1);
     114     6066823 :     if (!*pp) *pp = mod;
     115     5903429 :     else if (mod != *pp && !equalii(mod, *pp))
     116             :     {
     117           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_Fp");
     118           0 :       return 0;
     119             :     }
     120     6066823 :     return 1;
     121             :   case t_INT:
     122    59192096 :     return 1;
     123     4331520 :   default: return 0;
     124             :   }
     125             : }
     126             : 
     127             : int
     128    19349087 : RgX_is_FpX(GEN x, GEN *pp)
     129             : {
     130    19349087 :   long i, lx = lg(x);
     131    74994776 :   for (i=2; i<lx; i++)
     132    59131301 :     if (!Rg_is_Fp(gel(x, i), pp))
     133     3485612 :       return 0;
     134    15863475 :   return 1;
     135             : }
     136             : 
     137             : int
     138     2252615 : RgV_is_FpV(GEN x, GEN *pp)
     139             : {
     140     2252615 :   long i, lx = lg(x);
     141    11840379 :   for (i=1; i<lx; i++)
     142    10433672 :     if (!Rg_is_Fp(gel(x,i), pp)) return 0;
     143     1406707 :   return 1;
     144             : }
     145             : 
     146             : int
     147     1048339 : RgM_is_FpM(GEN x, GEN *pp)
     148             : {
     149     1048339 :   long i, lx = lg(x);
     150     2452645 :   for (i=1; i<lx; i++)
     151     2250193 :     if (!RgV_is_FpV(gel(x, i), pp)) return 0;
     152      202452 :   return 1;
     153             : }
     154             : 
     155             : int
     156       54782 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
     157             : {
     158             :   GEN pol, mod, p;
     159       54782 :   switch(typ(x))
     160             :   {
     161             :   case t_INTMOD:
     162       25459 :     return Rg_is_Fp(x, pp);
     163             :   case t_INT:
     164        7000 :     return 1;
     165             :   case t_POL:
     166          21 :     return RgX_is_FpX(x, pp);
     167             :   case t_FFELT:
     168       17353 :     mod = x; p = FF_p_i(x);
     169       17353 :     if (!*pp) *pp = p;
     170       17353 :     if (!*pT) *pT = mod;
     171       16520 :     else if (typ(*pT)!=t_FFELT || !FF_samefield(*pT,mod))
     172             :     {
     173          21 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     174          21 :       return 0;
     175             :     }
     176       17332 :     return 1;
     177             :   case t_POLMOD:
     178        4865 :     mod = gel(x,1); pol = gel(x, 2);
     179        4865 :     if (!RgX_is_FpX(mod, pp)) return 0;
     180        4865 :     if (typ(pol)==t_POL)
     181             :     {
     182        4858 :       if (!RgX_is_FpX(pol, pp)) return 0;
     183             :     }
     184           7 :     else if (!Rg_is_Fp(pol, pp)) return 0;
     185        4865 :     if (!*pT) *pT = mod;
     186        4431 :     else if (mod != *pT && !gequal(mod, *pT))
     187             :     {
     188           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     189           0 :       return 0;
     190             :     }
     191        4865 :     return 1;
     192             : 
     193          84 :   default: return 0;
     194             :   }
     195             : }
     196             : 
     197             : int
     198        2555 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
     199             : {
     200        2555 :   long i, lx = lg(x);
     201       56847 :   for (i = 2; i < lx; i++)
     202       54313 :     if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
     203        2534 :   return 1;
     204             : }
     205             : 
     206             : /************************************************************************
     207             :  **                                                                    **
     208             :  **                      Ring conversion                               **
     209             :  **                                                                    **
     210             :  ************************************************************************/
     211             : 
     212             : /* p > 0 a t_INT, return lift(x * Mod(1,p)).
     213             :  * If x is an INTMOD, assume modulus is a multiple of p. */
     214             : GEN
     215    21608028 : Rg_to_Fp(GEN x, GEN p)
     216             : {
     217    21608028 :   if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
     218      496427 :   switch(typ(x))
     219             :   {
     220       38541 :     case t_INT: return modii(x, p);
     221             :     case t_FRAC: {
     222          60 :       pari_sp av = avma;
     223          60 :       GEN z = modii(gel(x,1), p);
     224          60 :       if (z == gen_0) return gen_0;
     225          60 :       return gerepileuptoint(av, remii(mulii(z, Fp_inv(gel(x,2), p)), p));
     226             :     }
     227           0 :     case t_PADIC: return padic_to_Fp(x, p);
     228             :     case t_INTMOD: {
     229      457826 :       GEN q = gel(x,1), a = gel(x,2);
     230      457826 :       if (equalii(q, p)) return icopy(a);
     231          14 :       if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
     232           0 :       return remii(a, p);
     233             :     }
     234           0 :     default: pari_err_TYPE("Rg_to_Fp",x);
     235             :       return NULL; /* LCOV_EXCL_LINE */
     236             :   }
     237             : }
     238             : /* If x is a POLMOD, assume modulus is a multiple of T. */
     239             : GEN
     240       50443 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
     241             : {
     242       50443 :   long ta, tx = typ(x), v = get_FpX_var(T);
     243             :   GEN a, b;
     244       50443 :   if (is_const_t(tx))
     245             :   {
     246       38193 :     if (tx == t_FFELT)
     247             :     {
     248       18725 :       GEN z = FF_to_FpXQ(x);
     249       18725 :       setvarn(z, v);
     250       18725 :       return z;
     251             :     }
     252       19468 :     return scalar_ZX(Rg_to_Fp(x, p), v);
     253             :   }
     254       12250 :   switch(tx)
     255             :   {
     256             :     case t_POLMOD:
     257       11235 :       b = gel(x,1);
     258       11235 :       a = gel(x,2); ta = typ(a);
     259       11235 :       if (is_const_t(ta)) return scalar_ZX(Rg_to_Fp(a, p), v);
     260       11144 :       b = RgX_to_FpX(b, p); if (varn(b) != v) break;
     261       11144 :       a = RgX_to_FpX(a, p); if (ZX_equal(b,get_FpX_mod(T))) return a;
     262           0 :       if (signe(FpX_rem(b,T,p))==0) return FpX_rem(a, T, p);
     263           0 :       break;
     264             :     case t_POL:
     265        1015 :       if (varn(x) != v) break;
     266        1015 :       return FpX_rem(RgX_to_FpX(x,p), T, p);
     267             :     case t_RFRAC:
     268           0 :       a = Rg_to_FpXQ(gel(x,1), T,p);
     269           0 :       b = Rg_to_FpXQ(gel(x,2), T,p);
     270           0 :       return FpXQ_div(a,b, T,p);
     271             :   }
     272           0 :   pari_err_TYPE("Rg_to_FpXQ",x);
     273             :   return NULL; /* LCOV_EXCL_LINE */
     274             : }
     275             : GEN
     276      393228 : RgX_to_FpX(GEN x, GEN p)
     277             : {
     278             :   long i, l;
     279      393228 :   GEN z = cgetg_copy(x, &l); z[1] = x[1];
     280      393228 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     281      393228 :   return FpX_renormalize(z, l);
     282             : }
     283             : 
     284             : GEN
     285        1099 : RgV_to_FpV(GEN x, GEN p)
     286             : {
     287        1099 :   long i, l = lg(x);
     288        1099 :   GEN z = cgetg(l, t_VEC);
     289        1099 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     290        1099 :   return z;
     291             : }
     292             : 
     293             : GEN
     294      465233 : RgC_to_FpC(GEN x, GEN p)
     295             : {
     296      465233 :   long i, l = lg(x);
     297      465233 :   GEN z = cgetg(l, t_COL);
     298      465233 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     299      465233 :   return z;
     300             : }
     301             : 
     302             : GEN
     303       71004 : RgM_to_FpM(GEN x, GEN p)
     304             : {
     305       71004 :   long i, l = lg(x);
     306       71004 :   GEN z = cgetg(l, t_MAT);
     307       71004 :   for (i = 1; i < l; i++) gel(z,i) = RgC_to_FpC(gel(x,i), p);
     308       71004 :   return z;
     309             : }
     310             : GEN
     311       17116 : RgV_to_Flv(GEN x, ulong p)
     312             : {
     313       17116 :   long l = lg(x), i;
     314       17116 :   GEN a = cgetg(l, t_VECSMALL);
     315       17116 :   for (i=1; i<l; i++) a[i] = Rg_to_Fl(gel(x,i), p);
     316       17116 :   return a;
     317             : }
     318             : GEN
     319        2034 : RgM_to_Flm(GEN x, ulong p)
     320             : {
     321             :   long l, i;
     322        2034 :   GEN a = cgetg_copy(x, &l);
     323        2034 :   for (i=1; i<l; i++) gel(a,i) = RgV_to_Flv(gel(x,i), p);
     324        2034 :   return a;
     325             : }
     326             : 
     327             : GEN
     328         441 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
     329             : {
     330         441 :   long i, l = lg(x);
     331         441 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     332         441 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
     333         441 :   return FpXQX_renormalize(z, l);
     334             : }
     335             : GEN
     336         777 : RgX_to_FqX(GEN x, GEN T, GEN p)
     337             : {
     338         777 :   long i, l = lg(x);
     339         777 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     340         777 :   if (T)
     341       10878 :     for (i = 2; i < l; i++)
     342       10269 :       gel(z,i) = simplify_shallow(Rg_to_FpXQ(gel(x,i), T, p));
     343             :   else
     344       17976 :     for (i = 2; i < l; i++)
     345       17808 :       gel(z,i) = Rg_to_Fp(gel(x,i), p);
     346         777 :   return FpXQX_renormalize(z, l);
     347             : }
     348             : 
     349             : /* lg(V) > 1 */
     350             : GEN
     351      849765 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
     352             : {
     353      849765 :   pari_sp av = avma;
     354      849765 :   long i, l = lg(V);
     355      849765 :   GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
     356     4181499 :   for(i=2; i<l; i++)
     357             :   {
     358     3331734 :     z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
     359     3331734 :     if ((i & 7) == 0) z = gerepileupto(av, z);
     360             :   }
     361      849765 :   return gerepileupto(av, FpX_red(z,p));
     362             : }
     363             : 
     364             : GEN
     365        1260 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
     366             : {
     367        1260 :   long i, lz = lg(y);
     368             :   GEN z;
     369        1260 :   if (!T) return FpX_Fp_add(y, x, p);
     370        1260 :   if (lz == 2) return scalarpol(x, varn(y));
     371        1260 :   z = cgetg(lz,t_POL); z[1] = y[1];
     372        1260 :   gel(z,2) = Fq_add(gel(y,2),x, T, p);
     373        1260 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     374             :   else
     375         350 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     376        1260 :   return z;
     377             : }
     378             : 
     379             : GEN
     380        1048 : FqX_Fq_sub(GEN y, GEN x, GEN T, GEN p)
     381             : {
     382        1048 :   long i, lz = lg(y);
     383             :   GEN z;
     384        1048 :   if (!T) return FpX_Fp_sub(y, x, p);
     385        1048 :   if (lz == 2) return scalarpol(x, varn(y));
     386        1048 :   z = cgetg(lz,t_POL); z[1] = y[1];
     387        1048 :   gel(z,2) = Fq_sub(gel(y,2), x, T, p);
     388        1048 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     389             :   else
     390         926 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     391        1048 :   return z;
     392             : }
     393             : 
     394             : GEN
     395       63567 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
     396             : {
     397             :   long i, lP;
     398       63567 :   GEN res = cgetg_copy(P, &lP); res[1] = P[1];
     399       63567 :   for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
     400       63567 :   gel(res,lP-1) = gen_1; return res;
     401             : }
     402             : 
     403             : GEN
     404        3083 : FpXQX_normalize(GEN z, GEN T, GEN p)
     405             : {
     406             :   GEN lc;
     407        3083 :   if (lg(z) == 2) return z;
     408        3069 :   lc = leading_coeff(z);
     409        3069 :   if (typ(lc) == t_POL)
     410             :   {
     411        1550 :     if (lg(lc) > 3) /* non-constant */
     412        1442 :       return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
     413             :     /* constant */
     414         108 :     lc = gel(lc,2);
     415         108 :     z = shallowcopy(z);
     416         108 :     gel(z, lg(z)-1) = lc;
     417             :   }
     418             :   /* lc a t_INT */
     419        1627 :   if (equali1(lc)) return z;
     420          35 :   return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
     421             : }
     422             : 
     423             : GEN
     424      123403 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
     425             : {
     426             :   pari_sp av;
     427             :   GEN p1, r;
     428      123403 :   long j, i=lg(x)-1;
     429      123403 :   if (i<=2)
     430       24122 :     return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
     431       99281 :   av=avma; p1=gel(x,i);
     432             :   /* specific attention to sparse polynomials (see poleval)*/
     433             :   /*You've guessed it! It's a copy-paste(tm)*/
     434      291235 :   for (i--; i>=2; i=j-1)
     435             :   {
     436      192283 :     for (j=i; !signe(gel(x,j)); j--)
     437         329 :       if (j==2)
     438             :       {
     439         182 :         if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
     440         182 :         return gerepileupto(av, Fq_mul(p1,y, T, p));
     441             :       }
     442      191954 :     r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
     443      191954 :     p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
     444             :   }
     445       99099 :   return gerepileupto(av, p1);
     446             : }
     447             : 
     448             : GEN
     449       30380 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
     450             : {
     451       30380 :   long i, lb = lg(Q);
     452             :   GEN z;
     453       30380 :   if (!T) return FpXY_evalx(Q, x, p);
     454       20720 :   z = cgetg(lb, t_POL); z[1] = Q[1];
     455      115969 :   for (i=2; i<lb; i++)
     456             :   {
     457       95249 :     GEN q = gel(Q,i);
     458       95249 :     gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
     459             :   }
     460       20720 :   return FpXQX_renormalize(z, lb);
     461             : }
     462             : 
     463             : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
     464             : GEN
     465       12733 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
     466             : {
     467       12733 :   pari_sp av = avma;
     468       12733 :   if (!T) return FpXY_eval(Q, y, x, p);
     469         336 :   return gerepileupto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
     470             : }
     471             : 
     472             : /* a X^d */
     473             : GEN
     474      575897 : monomial(GEN a, long d, long v)
     475             : {
     476             :   long i, n;
     477             :   GEN P;
     478      575897 :   if (d < 0) {
     479           0 :     if (isrationalzero(a)) return pol_0(v);
     480           0 :     retmkrfrac(a, pol_xn(-d, v));
     481             :   }
     482      575897 :   if (gequal0(a))
     483             :   {
     484        8617 :     if (isexactzero(a)) return scalarpol_shallow(a,v);
     485           0 :     n = d+2; P = cgetg(n+1, t_POL);
     486           0 :     P[1] = evalsigne(0) | evalvarn(v);
     487             :   }
     488             :   else
     489             :   {
     490      567280 :     n = d+2; P = cgetg(n+1, t_POL);
     491      567280 :     P[1] = evalsigne(1) | evalvarn(v);
     492             :   }
     493      567280 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     494      567280 :   gel(P,i) = a; return P;
     495             : }
     496             : GEN
     497     7598948 : monomialcopy(GEN a, long d, long v)
     498             : {
     499             :   long i, n;
     500             :   GEN P;
     501     7598948 :   if (d < 0) {
     502           7 :     if (isrationalzero(a)) return pol_0(v);
     503           7 :     retmkrfrac(gcopy(a), pol_xn(-d, v));
     504             :   }
     505     7598941 :   if (gequal0(a))
     506             :   {
     507           7 :     if (isexactzero(a)) return scalarpol(a,v);
     508           0 :     n = d+2; P = cgetg(n+1, t_POL);
     509           0 :     P[1] = evalsigne(0) | evalvarn(v);
     510             :   }
     511             :   else
     512             :   {
     513     7598934 :     n = d+2; P = cgetg(n+1, t_POL);
     514     7598934 :     P[1] = evalsigne(1) | evalvarn(v);
     515             :   }
     516     7598934 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     517     7598934 :   gel(P,i) = gcopy(a); return P;
     518             : }
     519             : GEN
     520       19908 : pol_x_powers(long N, long v)
     521             : {
     522       19908 :   GEN L = cgetg(N+1,t_VEC);
     523             :   long i;
     524       19908 :   for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
     525       19908 :   return L;
     526             : }
     527             : 
     528             : GEN
     529           0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
     530             : {
     531           0 :   return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
     532             : }
     533             : 
     534             : GEN
     535           0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
     536             : {
     537           0 :   return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
     538             : }
     539             : 
     540             : /*******************************************************************/
     541             : /*                                                                 */
     542             : /*                             Fq                                  */
     543             : /*                                                                 */
     544             : /*******************************************************************/
     545             : 
     546             : GEN
     547     6951205 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     548             : {
     549             :   (void)T;
     550     6951205 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     551             :   {
     552     2382232 :     case 0: return Fp_add(x,y,p);
     553      204169 :     case 1: return FpX_Fp_add(x,y,p);
     554      338059 :     case 2: return FpX_Fp_add(y,x,p);
     555     4026745 :     case 3: return FpX_add(x,y,p);
     556             :   }
     557           0 :   return NULL;
     558             : }
     559             : 
     560             : GEN
     561     4850376 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
     562             : {
     563             :   (void)T;
     564     4850376 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     565             :   {
     566      155494 :     case 0: return Fp_sub(x,y,p);
     567        2370 :     case 1: return FpX_Fp_sub(x,y,p);
     568        9836 :     case 2: return Fp_FpX_sub(x,y,p);
     569     4682676 :     case 3: return FpX_sub(x,y,p);
     570             :   }
     571           0 :   return NULL;
     572             : }
     573             : 
     574             : GEN
     575      467999 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
     576             : {
     577             :   (void)T;
     578      467999 :   return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
     579             : }
     580             : 
     581             : GEN
     582       11499 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
     583             : {
     584             :   (void)T;
     585       11499 :   return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
     586             : }
     587             : 
     588             : /* If T==NULL do not reduce*/
     589             : GEN
     590    42526547 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
     591             : {
     592    42526547 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     593             :   {
     594     2380768 :     case 0: return Fp_mul(x,y,p);
     595      187587 :     case 1: return FpX_Fp_mul(x,y,p);
     596      143777 :     case 2: return FpX_Fp_mul(y,x,p);
     597    39814415 :     case 3: if (T) return FpXQ_mul(x,y,T,p);
     598     2904964 :             else return FpX_mul(x,y,p);
     599             :   }
     600           0 :   return NULL;
     601             : }
     602             : 
     603             : /* If T==NULL do not reduce*/
     604             : GEN
     605      731527 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
     606             : {
     607             :   (void) T;
     608      731527 :   return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
     609             : }
     610             : 
     611             : /* y t_INT */
     612             : GEN
     613       51729 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
     614             : {
     615             :   (void)T;
     616      103458 :   return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
     617       51729 :                           : Fp_mul(x,y,p);
     618             : }
     619             : /* If T==NULL do not reduce*/
     620             : GEN
     621      261491 : Fq_sqr(GEN x, GEN T, GEN p)
     622             : {
     623      261491 :   if (typ(x) == t_POL)
     624             :   {
     625       11241 :     if (T) return FpXQ_sqr(x,T,p);
     626           0 :     else return FpX_sqr(x,p);
     627             :   }
     628             :   else
     629      250250 :     return Fp_sqr(x,p);
     630             : }
     631             : 
     632             : GEN
     633           0 : Fq_neg_inv(GEN x, GEN T, GEN p)
     634             : {
     635           0 :   if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
     636           0 :   return FpXQ_inv(FpX_neg(x,p),T,p);
     637             : }
     638             : 
     639             : GEN
     640           0 : Fq_invsafe(GEN x, GEN pol, GEN p)
     641             : {
     642           0 :   if (typ(x) == t_INT) return Fp_invsafe(x,p);
     643           0 :   return FpXQ_invsafe(x,pol,p);
     644             : }
     645             : 
     646             : GEN
     647       29434 : Fq_inv(GEN x, GEN pol, GEN p)
     648             : {
     649       29434 :   if (typ(x) == t_INT) return Fp_inv(x,p);
     650       24674 :   return FpXQ_inv(x,pol,p);
     651             : }
     652             : 
     653             : GEN
     654      479248 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
     655             : {
     656      479248 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     657             :   {
     658      451976 :     case 0: return Fp_div(x,y,p);
     659       22806 :     case 1: return FpX_Fp_mul(x,Fp_inv(y,p),p);
     660         196 :     case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
     661        4270 :     case 3: return FpXQ_div(x,y,pol,p);
     662             :   }
     663           0 :   return NULL;
     664             : }
     665             : 
     666             : GEN
     667       19803 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
     668             : {
     669       19803 :   if (typ(x) == t_INT) return Fp_pow(x,n,p);
     670        8743 :   return FpXQ_pow(x,n,pol,p);
     671             : }
     672             : 
     673             : GEN
     674       12985 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
     675             : {
     676       12985 :   if (typ(x) == t_INT) return Fp_powu(x,n,p);
     677         553 :   return FpXQ_powu(x,n,pol,p);
     678             : }
     679             : 
     680             : GEN
     681      709098 : Fq_sqrt(GEN x, GEN T, GEN p)
     682             : {
     683      709098 :   if (typ(x) == t_INT)
     684             :   {
     685      698544 :     if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
     686         287 :     x = scalarpol_shallow(x, get_FpX_var(T));
     687             :   }
     688       10841 :   return FpXQ_sqrt(x,T,p);
     689             : }
     690             : GEN
     691       60411 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
     692             : {
     693       60411 :   if (typ(x) == t_INT)
     694             :   {
     695             :     long d;
     696       60215 :     if (!T) return Fp_sqrtn(x,n,p,zeta);
     697         477 :     d = get_FpX_degree(T);
     698         477 :     if (ugcd(umodiu(n,d),d) == 1)
     699             :     {
     700         414 :       if (!zeta)
     701           7 :         return Fp_sqrtn(x,n,p,NULL);
     702             :       else
     703             :       {
     704             :         /* gcd(n,p-1)=gcd(n,p^d-1) <=> same number of solutions if Fp and F_{p^d} */
     705         407 :         if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
     706         386 :           return Fp_sqrtn(x,n,p,zeta);
     707             :       }
     708             :     }
     709          84 :     x = scalarpol(x, get_FpX_var(T)); /* left on stack */
     710             :   }
     711         280 :   return FpXQ_sqrtn(x,n,T,p,zeta);
     712             : }
     713             : 
     714             : struct _Fq_field
     715             : {
     716             :   GEN T, p;
     717             : };
     718             : 
     719             : static GEN
     720        2268 : _Fq_red(void *E, GEN x)
     721        2268 : { struct _Fq_field *s = (struct _Fq_field *)E;
     722        2268 :   return Fq_red(x, s->T, s->p);
     723             : }
     724             : 
     725             : static GEN
     726        1344 : _Fq_add(void *E, GEN x, GEN y)
     727             : {
     728             :   (void) E;
     729        1344 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     730             :   {
     731          14 :     case 0: return addii(x,y);
     732           0 :     case 1: return ZX_Z_add(x,y);
     733         210 :     case 2: return ZX_Z_add(y,x);
     734        1120 :     default: return ZX_add(x,y);
     735             :   }
     736             : }
     737             : 
     738             : static GEN
     739         665 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
     740             : 
     741             : static GEN
     742        2219 : _Fq_mul(void *E, GEN x, GEN y)
     743             : {
     744             :   (void) E;
     745        2219 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     746             :   {
     747          63 :     case 0: return mulii(x,y);
     748         315 :     case 1: return ZX_Z_mul(x,y);
     749          56 :     case 2: return ZX_Z_mul(y,x);
     750        1785 :     default: return ZX_mul(x,y);
     751             :   }
     752             : }
     753             : 
     754             : static GEN
     755         322 : _Fq_inv(void *E, GEN x)
     756         322 : { struct _Fq_field *s = (struct _Fq_field *)E;
     757         322 :   return Fq_inv(x,s->T,s->p);
     758             : }
     759             : 
     760             : static int
     761         714 : _Fq_equal0(GEN x) { return signe(x)==0; }
     762             : 
     763             : static GEN
     764         448 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
     765             : 
     766             : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
     767             :                                        _Fq_inv,_Fq_equal0,_Fq_s};
     768             : 
     769         161 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
     770             : {
     771         161 :   GEN z = new_chunk(sizeof(struct _Fq_field));
     772         161 :   struct _Fq_field *e = (struct _Fq_field *) z;
     773         161 :   e->T = T; e->p  = p; *E = (void*)e;
     774         161 :   return &Fq_field;
     775             : }
     776             : 
     777             : /*******************************************************************/
     778             : /*                                                                 */
     779             : /*                             Fq[X]                               */
     780             : /*                                                                 */
     781             : /*******************************************************************/
     782             : /* P(X + c) */
     783             : GEN
     784           0 : FpX_translate(GEN P, GEN c, GEN p)
     785             : {
     786           0 :   pari_sp av = avma;
     787             :   GEN Q, *R;
     788             :   long i, k, n;
     789             : 
     790           0 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     791           0 :   Q = leafcopy(P);
     792           0 :   R = (GEN*)(Q+2); n = degpol(P);
     793           0 :   for (i=1; i<=n; i++)
     794             :   {
     795           0 :     for (k=n-i; k<n; k++)
     796           0 :       R[k] = Fp_add(R[k], Fp_mul(c, R[k+1], p), p);
     797             : 
     798           0 :     if (gc_needed(av,2))
     799             :     {
     800           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpX_translate, i = %ld/%ld", i,n);
     801           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     802             :     }
     803             :   }
     804           0 :   return gerepilecopy(av, FpX_renormalize(Q, lg(Q)));
     805             : }
     806             : /* P(X + c), c an Fq */
     807             : GEN
     808       43589 : FqX_translate(GEN P, GEN c, GEN T, GEN p)
     809             : {
     810       43589 :   pari_sp av = avma;
     811             :   GEN Q, *R;
     812             :   long i, k, n;
     813             : 
     814             :   /* signe works for t_(INT|POL) */
     815       43589 :   if (!signe(P) || !signe(c)) return RgX_copy(P);
     816       43589 :   Q = leafcopy(P);
     817       43589 :   R = (GEN*)(Q+2); n = degpol(P);
     818      194922 :   for (i=1; i<=n; i++)
     819             :   {
     820      553028 :     for (k=n-i; k<n; k++)
     821      401695 :       R[k] = Fq_add(R[k], Fq_mul(c, R[k+1], T, p), T, p);
     822             : 
     823      151333 :     if (gc_needed(av,2))
     824             :     {
     825           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FqX_translate, i = %ld/%ld", i,n);
     826           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     827             :     }
     828             :   }
     829       43589 :   return gerepilecopy(av, FpXQX_renormalize(Q, lg(Q)));
     830             : }
     831             : 
     832             : GEN
     833         658 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
     834             : {
     835         658 :   pari_sp ltop = avma;
     836             :   long k;
     837             :   GEN W;
     838         658 :   if (lgefint(p) == 3)
     839             :   {
     840         591 :     ulong pp = p[2];
     841         591 :     GEN Tl = ZX_to_Flx(T, pp);
     842         591 :     GEN Vl = FqV_to_FlxV(V, T, p);
     843         591 :     Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
     844         591 :     return gerepileupto(ltop, FlxX_to_ZXX(Tl));
     845             :   }
     846          67 :   W = cgetg(lg(V),t_VEC);
     847         381 :   for(k=1; k < lg(V); k++)
     848         314 :     gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
     849          67 :   return gerepileupto(ltop, FpXQXV_prod(W, T, p));
     850             : }
     851             : 
     852             : GEN
     853      123277 : FqV_red(GEN z, GEN T, GEN p)
     854             : {
     855      123277 :   long i, l = lg(z);
     856      123277 :   GEN res = cgetg(l, typ(z));
     857      123277 :   for(i=1;i<l;i++) gel(res,i) = Fq_red(gel(z,i),T,p);
     858      123277 :   return res;
     859             : }
     860             : 
     861             : GEN
     862           0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
     863             : {
     864           0 :   long i, lx = lg(x);
     865             :   GEN z;
     866           0 :   if (!T) return FpC_add(x, y, p);
     867           0 :   z = cgetg(lx, t_COL);
     868           0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_add(gel(x, i), gel(y, i), T, p);
     869           0 :   return z;
     870             : }
     871             : 
     872             : GEN
     873           0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
     874             : {
     875           0 :   long i, lx = lg(x);
     876             :   GEN z;
     877           0 :   if (!T) return FpC_sub(x, y, p);
     878           0 :   z = cgetg(lx, t_COL);
     879           0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_sub(gel(x, i), gel(y, i), T, p);
     880           0 :   return z;
     881             : }
     882             : 
     883             : GEN
     884           0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
     885             : {
     886           0 :   long i, l = lg(x);
     887             :   GEN z;
     888           0 :   if (!T) return FpC_Fp_mul(x, y, p);
     889           0 :   z = cgetg(l, t_COL);
     890           0 :   for (i=1;i<l;i++) gel(z,i) = Fq_mul(gel(x,i),y,T,p);
     891           0 :   return z;
     892             : }
     893             : 
     894             : GEN
     895         591 : FqV_to_FlxV(GEN v, GEN T, GEN pp)
     896             : {
     897         591 :   long j, N = lg(v);
     898         591 :   long vT = evalvarn(get_FpX_var(T));
     899         591 :   ulong p = pp[2];
     900         591 :   GEN y = cgetg(N, t_VEC);
     901        2874 :   for (j=1; j<N; j++)
     902        4566 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     903        2283 :                                     : ZX_to_Flx(gel(v,j), p));
     904         591 :   return y;
     905             : }
     906             : 
     907             : GEN
     908       41673 : FqC_to_FlxC(GEN v, GEN T, GEN pp)
     909             : {
     910       41673 :   long j, N = lg(v);
     911       41673 :   long vT = evalvarn(get_FpX_var(T));
     912       41673 :   ulong p = pp[2];
     913       41673 :   GEN y = cgetg(N, t_COL);
     914     1097112 :   for (j=1; j<N; j++)
     915     2346458 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     916     1291019 :                                     : ZX_to_Flx(gel(v,j), p));
     917       41673 :   return y;
     918             : }
     919             : 
     920             : GEN
     921        8010 : FqM_to_FlxM(GEN x, GEN T, GEN pp)
     922             : {
     923        8010 :   long j, n = lg(x);
     924        8010 :   GEN y = cgetg(n,t_MAT);
     925        8010 :   if (n == 1) return y;
     926       49683 :   for (j=1; j<n; j++)
     927       41673 :     gel(y,j) = FqC_to_FlxC(gel(x,j), T, pp);
     928        8010 :   return y;
     929             : }
     930             : 
     931             : /*******************************************************************/
     932             : /*                                                                 */
     933             : /*                          GENERIC CRT                            */
     934             : /*                                                                 */
     935             : /*******************************************************************/
     936             : 
     937             : static long
     938      216107 : get_nbprimes(ulong bound, ulong *pt_start)
     939             : {
     940             : #ifdef LONG_IS_64BIT
     941      185106 :   ulong pstart = 4611686018427388039UL;
     942             : #else
     943       31001 :   ulong pstart = 1073741827UL;
     944             : #endif
     945      216107 :   if (pt_start) *pt_start = pstart;
     946      216107 :   return (bound/expu(pstart))+1;
     947             : }
     948             : 
     949             : static GEN
     950      817862 : primelist_disc(ulong *p, long n, GEN dB)
     951             : {
     952      817862 :   GEN P = cgetg(n+1, t_VECSMALL);
     953             :   long i;
     954     2430262 :   for (i=1; i <= n; i++, *p = unextprime(*p+1))
     955             :   {
     956     1612400 :     if (dB && umodiu(dB, *p)==0) { i--; continue; }
     957     1612400 :     P[i] = *p;
     958             :   }
     959      817862 :   return P;
     960             : }
     961             : 
     962             : GEN
     963      216107 : gen_crt(const char *str, GEN worker, GEN dB, ulong bound, long mmin, GEN *pt_mod,
     964             :         GEN crt(GEN, GEN, GEN*), GEN center(GEN, GEN, GEN))
     965             : {
     966      216107 :   pari_sp av = avma;
     967             :   ulong p;
     968             :   long n, m;
     969             :   GEN  H, P, mod;
     970             :   pari_timer ti;
     971      216107 :   bound++; /* +1 to account for sign */
     972      216107 :   n = get_nbprimes(bound, &p);
     973      216107 :   m = minss(mmin, n);
     974      216107 :   if (DEBUGLEVEL > 4) timer_start(&ti);
     975      216107 :   if (m == 1)
     976             :   {
     977      132697 :     GEN P = primelist_disc(&p, n, dB);
     978      132697 :     GEN done = closure_callgen1(worker, P);
     979      132697 :     H = gel(done,1);
     980      132697 :     mod = gel(done,2);
     981      132697 :     if (center) H = center(H, mod, shifti(mod,-1));
     982      132697 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
     983             :   }
     984             :   else
     985             :   {
     986       83410 :     long i, s = (n+m-1)/m, r = m - (m*s-n), di = 0;
     987             :     struct pari_mt pt;
     988             :     long pending;
     989       83410 :     if (DEBUGLEVEL > 4)
     990           0 :       err_printf("%s: bound 2^%ld, nb primes: %ld\n",str, bound, n);
     991       83410 :     H = cgetg(m+1, t_VEC); P = cgetg(m+1, t_VEC);
     992       83410 :     mt_queue_start_lim(&pt, worker, m);
     993      832992 :     for (i=1; i<=m || pending; i++)
     994             :     {
     995             :       GEN done;
     996      749582 :       GEN pr = i <= m ? mkvec(primelist_disc(&p, i<=r ? s: s-1, dB)): NULL;
     997      749582 :       mt_queue_submit(&pt, i, pr);
     998      749582 :       done = mt_queue_get(&pt, NULL, &pending);
     999      749582 :       if (done)
    1000             :       {
    1001      685165 :         di++;
    1002      685165 :         gel(H, di) = gel(done,1);
    1003      685165 :         gel(P, di) = gel(done,2);
    1004      685165 :         if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
    1005             :       }
    1006             :     }
    1007       83410 :     mt_queue_end(&pt);
    1008       83410 :     if (DEBUGLEVEL>5) err_printf("\n");
    1009       83410 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
    1010       83410 :     H = crt(H, P, &mod);
    1011       83410 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: chinese", str);
    1012             :   }
    1013      432214 :   while ((ulong)expi(mod) < bound)
    1014             :   {
    1015           0 :     long s = get_nbprimes(bound-expi(mod), NULL);
    1016             :     GEN P, done, Hi, modi;
    1017           0 :     if (DEBUGLEVEL > 4) err_printf("%s: need %ld extra primes\n",str, s);
    1018           0 :     P = primelist_disc(&p, s, dB);
    1019           0 :     done = closure_callgen1(worker, P);
    1020           0 :     Hi = gel(done,1); modi = gel(done,2);
    1021           0 :     H = crt(mkvec2(H, Hi), mkvec2(mod, modi), &mod);
    1022           0 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: extra primes", str);
    1023             :   }
    1024      216107 :   if (pt_mod) *pt_mod = mod;
    1025      216107 :   gerepileall(av, pt_mod? 2: 1, &H, pt_mod);
    1026      216107 :   return H;
    1027             : }
    1028             : 
    1029             : /*******************************************************************/
    1030             : /*                                                                 */
    1031             : /*                          MODULAR GCD                            */
    1032             : /*                                                                 */
    1033             : /*******************************************************************/
    1034             : /* return z = a mod q, b mod p (p,q) = 1; qinv = 1/q mod p; a in ]-q,q] */
    1035             : static GEN
    1036     3126942 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq, GEN pq2)
    1037             : {
    1038     3126942 :   ulong d, amod = umodiu(a, p);
    1039     3126942 :   pari_sp av = avma;
    1040             :   GEN ax;
    1041             : 
    1042     3126942 :   if (b == amod) return NULL;
    1043     2169305 :   d = Fl_mul(Fl_sub(b, amod, p), qinv, p); /* != 0 */
    1044     2169305 :   if (d >= 1 + (p>>1))
    1045     1081530 :     ax = subii(a, mului(p-d, q));
    1046             :   else
    1047             :   {
    1048     1087775 :     ax = addii(a, mului(d, q)); /* in ]0, pq[ assuming a in ]-q,q[ */
    1049     1087775 :     if (cmpii(ax,pq2) > 0) ax = subii(ax,pq);
    1050             :   }
    1051     2169305 :   return gerepileuptoint(av, ax);
    1052             : }
    1053             : GEN
    1054         315 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
    1055             : GEN
    1056     3162840 : ZX_init_CRT(GEN Hp, ulong p, long v)
    1057             : {
    1058     3162840 :   long i, l = lg(Hp), lim = (long)(p>>1);
    1059     3162840 :   GEN H = cgetg(l, t_POL);
    1060     3162840 :   H[1] = evalsigne(1) | evalvarn(v);
    1061    11205571 :   for (i=2; i<l; i++)
    1062     8042731 :     gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
    1063     3162840 :   return H;
    1064             : }
    1065             : 
    1066             : GEN
    1067       96103 : ZM_init_CRT(GEN Hp, ulong p)
    1068             : {
    1069       96103 :   long i,j, m, l = lg(Hp), lim = (long)(p>>1);
    1070       96103 :   GEN c, cp, H = cgetg(l, t_MAT);
    1071       96103 :   if (l==1) return H;
    1072       52696 :   m = lgcols(Hp);
    1073      169000 :   for (j=1; j<l; j++)
    1074             :   {
    1075      116304 :     cp = gel(Hp,j);
    1076      116304 :     c = cgetg(m, t_COL);
    1077      116304 :     gel(H,j) = c;
    1078      116304 :     for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
    1079             :   }
    1080       52696 :   return H;
    1081             : }
    1082             : 
    1083             : int
    1084        6951 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
    1085             : {
    1086        6951 :   GEN h, q = *ptq, qp = muliu(q,p);
    1087        6951 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1088        6951 :   int stable = 1;
    1089        6951 :   h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp,shifti(qp,-1));
    1090        6951 :   if (h) { *H = h; stable = 0; }
    1091        6951 :   *ptq = qp; return stable;
    1092             : }
    1093             : 
    1094             : static int
    1095      171945 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
    1096             : {
    1097      171945 :   GEN H = *ptH, h, qp2 = shifti(qp,-1);
    1098      171945 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1099      171945 :   long i, l = lg(H), lp = lg(Hp);
    1100      171945 :   int stable = 1;
    1101             : 
    1102      171945 :   if (l < lp)
    1103             :   { /* degree increases */
    1104           0 :     GEN x = cgetg(lp, t_POL);
    1105           0 :     for (i=1; i<l; i++)  x[i] = H[i];
    1106           0 :     for (   ; i<lp; i++) gel(x,i) = gen_0;
    1107           0 :     *ptH = H = x;
    1108           0 :     stable = 0;
    1109      171945 :   } else if (l > lp)
    1110             :   { /* degree decreases */
    1111           0 :     GEN x = cgetg(l, t_VECSMALL);
    1112           0 :     for (i=1; i<lp; i++)  x[i] = Hp[i];
    1113           0 :     for (   ; i<l; i++) x[i] = 0;
    1114           0 :     Hp = x; lp = l;
    1115             :   }
    1116     1409881 :   for (i=2; i<lp; i++)
    1117             :   {
    1118     1237936 :     h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp,qp2);
    1119     1237936 :     if (h) { gel(H,i) = h; stable = 0; }
    1120             :   }
    1121      171945 :   return stable;
    1122             : }
    1123             : 
    1124             : int
    1125        1168 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
    1126             : {
    1127        1168 :   GEN q = *ptq, qp = muliu(q,p);
    1128        1168 :   int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
    1129        1168 :   *ptq = qp; return stable;
    1130             : }
    1131             : 
    1132             : int
    1133       17803 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1134             : {
    1135       17803 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1136       17803 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1137       17803 :   long i,j, l = lg(H), m = lgcols(H);
    1138       17803 :   int stable = 1;
    1139       71739 :   for (j=1; j<l; j++)
    1140     1790703 :     for (i=1; i<m; i++)
    1141             :     {
    1142     1736767 :       h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp,qp2);
    1143     1736767 :       if (h) { gcoeff(H,i,j) = h; stable = 0; }
    1144             :     }
    1145       17803 :   *ptq = qp; return stable;
    1146             : }
    1147             : 
    1148             : GEN
    1149        1449 : ZVM_init_CRT(GEN Hp, ulong p)
    1150             : {
    1151             :   long i, j, k;
    1152             :   GEN c, cp, d, dp, H;
    1153        1449 :   long m, l = lg(Hp), lim = (long)(p>>1), n;
    1154        1449 :   H = cgetg(l, t_MAT);
    1155        1449 :   if (l==1) return H;
    1156        1449 :   m = lgcols(Hp);
    1157        1449 :   n = lg(gmael(Hp,1,1));
    1158        6097 :   for (j=1; j<l; j++)
    1159             :   {
    1160        4648 :     cp = gel(Hp,j);
    1161        4648 :     c = cgetg(m, t_COL);
    1162        4648 :     gel(H,j) = c;
    1163       67963 :     for (i=1; i<m; i++)
    1164             :     {
    1165       63315 :       dp = gel(cp, i);
    1166       63315 :       d = cgetg(n, t_VEC);
    1167       63315 :       gel(c, i) = d;
    1168      190827 :       for (k=1; k<n; k++)
    1169      127512 :         gel(d,k) = stoi(Fl_center(dp[k], p, lim));
    1170             :     }
    1171             :   }
    1172        1449 :   return H;
    1173             : }
    1174             : 
    1175             : int
    1176         590 : ZVM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1177             : {
    1178         590 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1179         590 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1180         590 :   long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
    1181         590 :   int stable = 1;
    1182        3873 :   for (j=1; j<l; j++)
    1183       75927 :     for (i=1; i<m; i++)
    1184      217932 :       for (k=1; k<n; k++)
    1185             :       {
    1186      145288 :         h = Fl_chinese_coprime(gmael3(H,j,i,k),mael3(Hp,j,i,k),q,p,qinv,qp,qp2);
    1187      145288 :         if (h) { gmael3(H,j,i,k) = h; stable = 0; }
    1188             :       }
    1189         590 :   *ptq = qp; return stable;
    1190             : }
    1191             : 
    1192             : /* record the degrees of Euclidean remainders (make them as large as
    1193             :  * possible : smaller values correspond to a degenerate sequence) */
    1194             : static void
    1195        1561 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
    1196             : {
    1197             :   long da,db,dc, ind;
    1198        1561 :   pari_sp av = avma;
    1199             : 
    1200        1561 :   if (lgpol(a)==0 || lgpol(b)==0) return;
    1201        1561 :   da = degpol(a);
    1202        1561 :   db = degpol(b);
    1203        1561 :   if (db > da)
    1204           0 :   { swapspec(a,b, da,db); }
    1205        1561 :   else if (!da) return;
    1206        1561 :   ind = 0;
    1207        9814 :   while (db)
    1208             :   {
    1209        6692 :     GEN c = Flx_rem(a,b, p);
    1210        6692 :     a = b; b = c; dc = degpol(c);
    1211        6692 :     if (dc < 0) break;
    1212             : 
    1213        6692 :     ind++;
    1214        6692 :     if (dc > dglist[ind]) dglist[ind] = dc;
    1215        6692 :     if (gc_needed(av,2))
    1216             :     {
    1217           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1218           0 :       gerepileall(av, 2, &a,&b);
    1219             :     }
    1220        6692 :     db = dc; /* = degpol(b) */
    1221             :   }
    1222        1561 :   if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
    1223        1561 :   avma = av; return;
    1224             : }
    1225             : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
    1226             :  * "generic" degree sequence as given by dglist, set *Ci and return
    1227             :  * resultant(a,b). Modular version of Collins's subresultant */
    1228             : static ulong
    1229        6996 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
    1230             : {
    1231             :   long da,db,dc, ind;
    1232        6996 :   ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
    1233        6996 :   int s = 1;
    1234        6996 :   pari_sp av = avma;
    1235             : 
    1236        6996 :   *C0 = 1; *C1 = 0;
    1237        6996 :   if (lgpol(a)==0 || lgpol(b)==0) return 0;
    1238        6996 :   da = degpol(a);
    1239        6996 :   db = degpol(b);
    1240        6996 :   if (db > da)
    1241             :   {
    1242           0 :     swapspec(a,b, da,db);
    1243           0 :     if (both_odd(da,db)) s = -s;
    1244             :   }
    1245        6996 :   else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
    1246        6996 :   ind = 0;
    1247       40727 :   while (db)
    1248             :   { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
    1249             :      * da = deg a, db = deg b */
    1250       27099 :     GEN c = Flx_rem(a,b, p);
    1251       27099 :     long delta = da - db;
    1252             : 
    1253       27099 :     if (both_odd(da,db)) s = -s;
    1254       27099 :     lb = Fl_mul(b[db+2], cb, p);
    1255       27099 :     a = b; b = c; dc = degpol(c);
    1256       27099 :     ind++;
    1257       27099 :     if (dc != dglist[ind]) { avma = av; return 0; } /* degenerates */
    1258       26735 :     if (g == h)
    1259             :     { /* frequent */
    1260       24621 :       ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
    1261       24621 :       ca = cb;
    1262       24621 :       cb = cc;
    1263             :     }
    1264             :     else
    1265             :     {
    1266        2114 :       ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
    1267        2114 :       ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
    1268        2114 :       ca = cb;
    1269        2114 :       cb = Fl_div(cc, ghdelta, p);
    1270             :     }
    1271       26735 :     da = db; /* = degpol(a) */
    1272       26735 :     db = dc; /* = degpol(b) */
    1273             : 
    1274       26735 :     g = lb;
    1275       26735 :     if (delta == 1)
    1276       17482 :       h = g; /* frequent */
    1277             :     else
    1278        9253 :       h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
    1279             : 
    1280       26735 :     if (gc_needed(av,2))
    1281             :     {
    1282           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1283           0 :       gerepileall(av, 2, &a,&b);
    1284             :     }
    1285             :   }
    1286        6632 :   if (da > 1) return 0; /* Failure */
    1287             :   /* last non-constant polynomial has degree 1 */
    1288        6632 :   *C0 = Fl_mul(ca, a[2], p);
    1289        6632 :   *C1 = Fl_mul(ca, a[3], p);
    1290        6632 :   res = Fl_mul(cb, b[2], p);
    1291        6632 :   if (s == -1) res = p - res;
    1292        6632 :   avma = av; return res;
    1293             : }
    1294             : 
    1295             : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
    1296             :  * Return 0 in case of degree drop. */
    1297             : static GEN
    1298        8557 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
    1299             : {
    1300             :   GEN z;
    1301        8557 :   long i, lb = lg(Q);
    1302        8557 :   ulong leadz = Flx_eval(leading_coeff(Q), x, p);
    1303        8557 :   long vs=mael(Q,2,1);
    1304        8557 :   if (!leadz) return zero_Flx(vs);
    1305             : 
    1306        8557 :   z = cgetg(lb, t_VECSMALL); z[1] = vs;
    1307        8557 :   for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
    1308        8557 :   z[i] = leadz; return z;
    1309             : }
    1310             : 
    1311             : GEN
    1312       17836 : FpXY_Fq_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
    1313             : {
    1314       17836 :   pari_sp av = avma;
    1315       17836 :   long i, lb = lg(Q);
    1316             :   GEN z;
    1317       17836 :   if (!T) return FpXY_evaly(Q, y, p, vx);
    1318        1148 :   if (lb == 2) return pol_0(vx);
    1319        1148 :   z = gel(Q, lb-1);
    1320        1148 :   if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
    1321             : 
    1322        1148 :   if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
    1323       28084 :   for (i=lb-2; i>=2; i--)
    1324             :   {
    1325       26936 :     GEN c = gel(Q,i);
    1326       26936 :     z = FqX_Fq_mul(z, y, T, p);
    1327       26936 :     z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
    1328             :   }
    1329        1148 :   return gerepileupto(av, z);
    1330             : }
    1331             : 
    1332             : static GEN
    1333       15246 : ZX_norml1(GEN x)
    1334             : {
    1335       15246 :   long i, l = lg(x);
    1336             :   GEN s;
    1337             : 
    1338       15246 :   if (l == 2) return gen_0;
    1339        8624 :   s = gel(x, l-1); /* != 0 */
    1340       31556 :   for (i = l-2; i > 1; i--) {
    1341       22932 :     GEN xi = gel(x,i);
    1342       22932 :     if (!signe(x)) continue;
    1343       22932 :     s = addii_sign(s,1, xi,1);
    1344             :   }
    1345        8624 :   return s;
    1346             : }
    1347             : 
    1348             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1349             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1350             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1351             :  * Return e such that Res(A, B) < 2^e */
    1352             : ulong
    1353       75297 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
    1354             : {
    1355       75297 :   pari_sp av = avma, av2;
    1356       75297 :   GEN a = gen_0, b = gen_0;
    1357       75297 :   long i , lA = lg(A), lB = lg(B);
    1358             :   double loga, logb;
    1359      856911 :   for (i=2; i<lA; i++)
    1360             :   {
    1361      781614 :     a = addii(a, sqri(gel(A,i)));
    1362      781614 :     if (gc_needed(av,1))
    1363             :     {
    1364           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1365           0 :       a = gerepileupto(av, a);
    1366             :     }
    1367             :   }
    1368       75297 :   a = gerepileuptoint(av, a);
    1369       75297 :   av2 = avma;
    1370      787756 :   for (i=2; i<lB; i++)
    1371             :   {
    1372      712459 :     GEN t = gel(B,i);
    1373      712459 :     if (typ(t) == t_POL) t = ZX_norml1(t);
    1374      712459 :     b = addii(b, sqri(t));
    1375      712459 :     if (gc_needed(av2,1))
    1376             :     {
    1377           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1378           0 :       b = gerepileupto(av2, b);
    1379             :     }
    1380             :   }
    1381       75297 :   loga = dbllog2(a);
    1382       75297 :   logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
    1383       75297 :   i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
    1384       75297 :   avma = av; return (i <= 0)? 1: 1 + (ulong)i;
    1385             : }
    1386             : 
    1387             : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
    1388             : static ulong
    1389      258606 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong la)
    1390             : {
    1391      258606 :   GEN ev = FlxY_evalx(b, n, p);
    1392      258553 :   long drop = lg(b) - lg(ev);
    1393      258553 :   ulong r = Flx_resultant(a, ev, p);
    1394      258589 :   if (drop && la != 1) r = Fl_mul(r, Fl_powu(la, drop,p),p);
    1395      258589 :   return r;
    1396             : }
    1397             : static GEN
    1398         137 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
    1399             : {
    1400         137 :   GEN ev = FpXY_evaly(b, n, p, vX);
    1401         137 :   long drop = db-degpol(ev);
    1402         137 :   GEN r = FpX_resultant(a, ev, p);
    1403         137 :   if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
    1404         137 :   return r;
    1405             : }
    1406             : 
    1407             : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1408             : /* Return a Fly */
    1409             : static GEN
    1410       12663 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, long dres, long sx)
    1411             : {
    1412             :   long i;
    1413       12663 :   ulong n, la = Flx_lead(a);
    1414       12662 :   GEN  x = cgetg(dres+2, t_VECSMALL);
    1415       12661 :   GEN  y = cgetg(dres+2, t_VECSMALL);
    1416             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1417             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1418      137079 :   for (i=0,n = 1; i < dres; n++)
    1419             :   {
    1420      124416 :     x[++i] = n;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1421      124404 :     x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1422             :   }
    1423       12663 :   if (i == dres)
    1424             :   {
    1425        9867 :     x[++i] = 0;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1426             :   }
    1427       12664 :   return Flv_polint(x,y, p, sx);
    1428             : }
    1429             : 
    1430             : static GEN
    1431        9435 : FlxX_pseudorem(GEN x, GEN y, ulong p)
    1432             : {
    1433        9435 :   long vx = varn(x), dx, dy, dz, i, lx, dp;
    1434        9435 :   pari_sp av = avma, av2;
    1435             : 
    1436        9435 :   if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
    1437        9435 :   (void)new_chunk(2);
    1438        9437 :   dx=degpol(x); x = RgX_recip_shallow(x)+2;
    1439        9443 :   dy=degpol(y); y = RgX_recip_shallow(y)+2; dz=dx-dy; dp = dz+1;
    1440        9442 :   av2 = avma;
    1441             :   for (;;)
    1442             :   {
    1443       72411 :     gel(x,0) = Flx_neg(gel(x,0), p); dp--;
    1444      278213 :     for (i=1; i<=dy; i++)
    1445      409400 :       gel(x,i) = Flx_add( Flx_mul(gel(y,0), gel(x,i), p),
    1446      204700 :                               Flx_mul(gel(x,0), gel(y,i), p), p );
    1447     1273650 :     for (   ; i<=dx; i++)
    1448     1201360 :       gel(x,i) = Flx_mul(gel(y,0), gel(x,i), p);
    1449       76615 :     do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
    1450       72290 :     if (dx < dy) break;
    1451       62850 :     if (gc_needed(av2,1))
    1452             :     {
    1453           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
    1454           0 :       gerepilecoeffs(av2,x,dx+1);
    1455             :     }
    1456       62969 :   }
    1457        9440 :   if (dx < 0) return zero_Flx(0);
    1458        9440 :   lx = dx+3; x -= 2;
    1459        9440 :   x[0]=evaltyp(t_POL) | evallg(lx);
    1460        9438 :   x[1]=evalsigne(1) | evalvarn(vx);
    1461        9438 :   x = RgX_recip_shallow(x);
    1462        9439 :   if (dp)
    1463             :   { /* multiply by y[0]^dp   [beware dummy vars from FpX_FpXY_resultant] */
    1464        2188 :     GEN t = Flx_powu(gel(y,0), dp, p);
    1465        8747 :     for (i=2; i<lx; i++)
    1466        6561 :       gel(x,i) = Flx_mul(gel(x,i), t, p);
    1467             :   }
    1468        9437 :   return gerepilecopy(av, x);
    1469             : }
    1470             : 
    1471             : /* return a Flx */
    1472             : GEN
    1473        3091 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
    1474             : {
    1475        3091 :   pari_sp av = avma, av2;
    1476             :   long degq,dx,dy,du,dv,dr,signh;
    1477             :   GEN z,g,h,r,p1;
    1478             : 
    1479        3091 :   dx=degpol(u); dy=degpol(v); signh=1;
    1480        3094 :   if (dx < dy)
    1481             :   {
    1482          42 :     swap(u,v); lswap(dx,dy);
    1483          42 :     if (both_odd(dx, dy)) signh = -signh;
    1484             :   }
    1485        3094 :   if (dy < 0) return zero_Flx(sx);
    1486        3094 :   if (dy==0) return gerepileupto(av, Flx_powu(gel(v,2),dx,p));
    1487             : 
    1488        3094 :   g = h = pol1_Flx(sx); av2 = avma;
    1489             :   for(;;)
    1490             :   {
    1491        9439 :     r = FlxX_pseudorem(u,v,p); dr = lg(r);
    1492        9437 :     if (dr == 2) { avma = av; return zero_Flx(sx); }
    1493        9437 :     du = degpol(u); dv = degpol(v); degq = du-dv;
    1494        9440 :     u = v; p1 = g; g = leading_coeff(u);
    1495        9443 :     switch(degq)
    1496             :     {
    1497           0 :       case 0: break;
    1498             :       case 1:
    1499        6975 :         p1 = Flx_mul(h,p1, p); h = g; break;
    1500             :       default:
    1501        2468 :         p1 = Flx_mul(Flx_powu(h,degq,p), p1, p);
    1502        2468 :         h = Flx_div(Flx_powu(g,degq,p), Flx_powu(h,degq-1,p), p);
    1503             :     }
    1504        9437 :     if (both_odd(du,dv)) signh = -signh;
    1505        9438 :     v = FlxY_Flx_div(r, p1, p);
    1506        9440 :     if (dr==3) break;
    1507        6346 :     if (gc_needed(av2,1))
    1508             :     {
    1509           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_resultant, dr = %ld",dr);
    1510           0 :       gerepileall(av2,4, &u, &v, &g, &h);
    1511             :     }
    1512        6346 :   }
    1513        3094 :   z = gel(v,2);
    1514        3094 :   if (dv > 1) z = Flx_div(Flx_powu(z,dv,p), Flx_powu(h,dv-1,p), p);
    1515        3094 :   if (signh < 0) z = Flx_neg(z,p);
    1516        3094 :   return gerepileupto(av, z);
    1517             : }
    1518             : 
    1519             : /* Warning:
    1520             :  * This function switches between valid and invalid variable ordering*/
    1521             : 
    1522             : static GEN
    1523        3199 : FlxY_to_FlyX(GEN b, long sv)
    1524             : {
    1525        3199 :   long i, n=-1;
    1526        3199 :   long sw = b[1]&VARNBITS;
    1527        3199 :   for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
    1528        3198 :   return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
    1529             : }
    1530             : 
    1531             : /* Return a Fly*/
    1532             : GEN
    1533        3195 : Flx_FlxY_resultant(GEN a, GEN b, ulong pp)
    1534             : {
    1535        3195 :   pari_sp ltop=avma;
    1536        3195 :   long dres = degpol(a)*degpol(b);
    1537        3197 :   long sx=a[1], sy=b[1]&VARNBITS;
    1538             :   GEN z;
    1539        3197 :   b = FlxY_to_FlyX(b,sx);
    1540        3195 :   if ((ulong)dres >= pp)
    1541        3091 :     z = FlxX_resultant(Fly_to_FlxY(a, sy), b, pp, sx);
    1542             :   else
    1543         104 :     z = Flx_FlxY_resultant_polint(a, b, pp, (ulong)dres, sy);
    1544        3199 :   return gerepileupto(ltop,z);
    1545             : }
    1546             : 
    1547             : /* return a t_POL (in variable v >= 0) whose coeffs are the coeffs of b,
    1548             :  * in variable v. This is an incorrect PARI object if initially varn(b) << v.
    1549             :  * We could return a vector of coeffs, but it is convenient to have degpol()
    1550             :  * and friends available. Even in that case, it will behave nicely with all
    1551             :  * functions treating a polynomial as a vector of coeffs (eg poleval).
    1552             :  * FOR INTERNAL USE! */
    1553             : GEN
    1554        9478 : swap_vars(GEN b0, long v)
    1555             : {
    1556        9478 :   long i, n = RgX_degree(b0, v);
    1557             :   GEN b, x;
    1558        9478 :   if (n < 0) return pol_0(v);
    1559        9478 :   b = cgetg(n+3, t_POL); x = b + 2;
    1560        9478 :   b[1] = evalsigne(1) | evalvarn(v);
    1561        9478 :   for (i=0; i<=n; i++) gel(x,i) = polcoeff_i(b0, i, v);
    1562        9478 :   return b;
    1563             : }
    1564             : 
    1565             : /* assume varn(b) << varn(a) */
    1566             : /* return a FpY*/
    1567             : GEN
    1568        3174 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
    1569             : {
    1570        3174 :   long i,n,dres, db, vY = varn(b), vX = varn(a);
    1571             :   GEN la,x,y;
    1572             : 
    1573        3174 :   if (lgefint(p) == 3)
    1574             :   {
    1575        3166 :     ulong pp = uel(p,2);
    1576        3166 :     b = ZXX_to_FlxX(b, pp, vX);
    1577        3170 :     a = ZX_to_Flx(a, pp);
    1578        3168 :     x = Flx_FlxY_resultant(a, b, pp);
    1579        3172 :     return Flx_to_ZX(x);
    1580             :   }
    1581           8 :   db = RgXY_degreex(b);
    1582           8 :   dres = degpol(a)*degpol(b);
    1583           8 :   la = leading_coeff(a);
    1584           8 :   x = cgetg(dres+2, t_VEC);
    1585           8 :   y = cgetg(dres+2, t_VEC);
    1586             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1587             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1588          73 :   for (i=0,n = 1; i < dres; n++)
    1589             :   {
    1590          65 :     gel(x,++i) = utoipos(n);
    1591          65 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1592          65 :     gel(x,++i) = subiu(p,n);
    1593          65 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1594             :   }
    1595           8 :   if (i == dres)
    1596             :   {
    1597           7 :     gel(x,++i) = gen_0;
    1598           7 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
    1599             :   }
    1600           8 :   return FpV_polint(x,y, p, vY);
    1601             : }
    1602             : 
    1603             : GEN
    1604         490 : FpX_direct_compositum(GEN a, GEN b, GEN p)
    1605             : {
    1606         490 :   GEN x = deg1pol_shallow(gen_1, pol_x(varn(a)), fetch_var_higher()); /* x+y */
    1607         490 :   x = FpX_FpXY_resultant(a, poleval(b,x),p);
    1608         490 :   (void)delete_var(); return x;
    1609             : }
    1610             : 
    1611             : /* 0, 1, -1, 2, -2, ... */
    1612             : #define next_lambda(a) (a>0 ? -a : 1-a)
    1613             : GEN
    1614           0 : FpX_compositum(GEN a, GEN b, GEN p)
    1615             : {
    1616           0 :   long k, v = fetch_var_higher();
    1617           0 :   for (k = 1;; k = next_lambda(k))
    1618             :   {
    1619           0 :     GEN x = deg1pol_shallow(gen_1, gmulsg(k, pol_x(v)), 0); /* x + k y */
    1620           0 :     GEN C = FpX_FpXY_resultant(a, poleval(b,x),p);
    1621           0 :     if (FpX_is_squarefree(C, p)) { (void)delete_var(); return C; }
    1622           0 :   }
    1623             : }
    1624             : 
    1625             : /* Assume A in Z[Y], B in Q[Y][X], and Res_Y(A, B) in Z[X].
    1626             :  * If lambda = NULL, return Res_Y(A,B).
    1627             :  * Otherwise, find a small lambda (start from *lambda, use the sequence above)
    1628             :  * such that R(X) = Res_Y(A(Y), B(X + lambda Y)) is squarefree, reset *lambda
    1629             :  * to the chosen value and return R
    1630             :  *
    1631             :  * If LERS is non-NULL, set it to the Last non-constant polynomial in the
    1632             :  * Euclidean Remainder Sequence */
    1633             : static GEN
    1634        1582 : ZX_ZXY_resultant_LERS(GEN A, GEN B0, long *plambda, GEN *LERS)
    1635             : {
    1636        1582 :   int checksqfree = plambda? 1: 0, stable;
    1637        1582 :   long lambda = plambda? *plambda: 0, cnt = 0;
    1638             :   ulong bound, dp;
    1639        1582 :   pari_sp av = avma, av2 = 0;
    1640        1582 :   long i,n, degA = degpol(A), degB, dres = degA*degpol(B0);
    1641        1582 :   long v = fetch_var_higher();
    1642        1582 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    1643        1582 :   long sX = evalvarn(vX);
    1644             :   GEN x, y, dglist, dB, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
    1645             :   forprime_t S;
    1646             : 
    1647        1582 :   dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
    1648        1582 :   if (LERS)
    1649             :   {
    1650        1582 :     if (!checksqfree)
    1651           0 :       pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
    1652        1582 :     C0 = cgetg(dres+2, t_VECSMALL);
    1653        1582 :     C1 = cgetg(dres+2, t_VECSMALL);
    1654        1582 :     dglist = cgetg(dres+1, t_VECSMALL);
    1655             :   }
    1656        1582 :   x = cgetg(dres+2, t_VECSMALL);
    1657        1582 :   y = cgetg(dres+2, t_VECSMALL);
    1658        1582 :   B0 = Q_remove_denom(B0, &dB);
    1659        1582 :   if (!dB) B0 = leafcopy(B0);
    1660        1582 :   A = leafcopy(A);
    1661        1582 :   B = B0;
    1662        1582 :   setvarn(A,v);
    1663             :   /* make sure p large enough */
    1664             : INIT:
    1665             :   /* always except the first time */
    1666        2317 :   if (av2) { avma = av2; lambda = next_lambda(lambda); }
    1667        2317 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    1668        2317 :   B = swap_vars(B, vY); setvarn(B,v);
    1669             :   /* B0(lambda v + x, v) */
    1670        2317 :   if (DEBUGLEVEL>4 && checksqfree) err_printf("Trying lambda = %ld\n", lambda);
    1671        2317 :   av2 = avma;
    1672             : 
    1673        2317 :   if (degA <= 3)
    1674             :   { /* sub-resultant faster for small degrees */
    1675        2121 :     if (LERS)
    1676             :     { /* implies checksqfree */
    1677        2121 :       H = RgX_resultant_all(A,B,&q);
    1678        2121 :       if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
    1679        1435 :       H0 = gel(q,2);
    1680        1435 :       if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
    1681        1435 :       H1 = gel(q,3);
    1682        1435 :       if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
    1683             :     }
    1684             :     else
    1685           0 :       H = resultant(A,B);
    1686        1435 :     if (checksqfree && !ZX_is_squarefree(H)) goto INIT;
    1687        1393 :     if (dB) H = ZX_Z_divexact(H, powiu(dB, degA));
    1688        1393 :     goto END;
    1689             :   }
    1690             : 
    1691         196 :   H = H0 = H1 = NULL;
    1692         196 :   degB = degpol(B);
    1693         196 :   bound = ZX_ZXY_ResBound(A, B, dB);
    1694         196 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    1695         196 :   dp = 1;
    1696         196 :   init_modular_big(&S);
    1697             :   while (1)
    1698             :   {
    1699         397 :     ulong p = u_forprime_next(&S);
    1700         397 :     if (dB) { dp = umodiu(dB, p); if (!dp) continue; }
    1701             : 
    1702         397 :     a = ZX_to_Flx(A, p);
    1703         397 :     b = ZXX_to_FlxX(B, p, varn(A));
    1704         397 :     if (LERS)
    1705             :     {
    1706             :       GEN Hi;
    1707         397 :       if (degpol(a) < degA || degpol(b) < degB) continue; /* p | lc(A)lc(B) */
    1708         397 :       if (checksqfree)
    1709             :       { /* find degree list for generic Euclidean Remainder Sequence */
    1710         196 :         long goal = minss(degpol(a), degpol(b)); /* longest possible */
    1711         196 :         for (n=1; n <= goal; n++) dglist[n] = 0;
    1712         196 :         setlg(dglist, 1);
    1713        1645 :         for (n=0; n <= dres; n++)
    1714             :         {
    1715        1561 :           ev = FlxY_evalx_drop(b, n, p);
    1716        1561 :           Flx_resultant_set_dglist(a, ev, dglist, p);
    1717        1561 :           if (lg(dglist)-1 == goal) break;
    1718             :         }
    1719             :         /* last pol in ERS has degree > 1 ? */
    1720         196 :         goal = lg(dglist)-1;
    1721         196 :         if (degpol(B) == 1) { if (!goal) goto INIT; }
    1722             :         else
    1723             :         {
    1724         189 :           if (goal <= 1) goto INIT;
    1725         182 :           if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
    1726             :         }
    1727         189 :         if (DEBUGLEVEL>4)
    1728           0 :           err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
    1729             :       }
    1730             : 
    1731        7386 :       for (i=0,n = 0; i <= dres; n++)
    1732             :       {
    1733        6996 :         ev = FlxY_evalx_drop(b, n, p);
    1734        6996 :         x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
    1735        6996 :         if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
    1736             :       }
    1737         390 :       Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
    1738         390 :       Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
    1739             :     }
    1740             :     else
    1741             :     {
    1742           0 :       long dropa = degA - degpol(a), dropb = degB - degpol(b);
    1743           0 :       Hp = Flx_FlxY_resultant_polint(a, b, p, (ulong)dres, sX);
    1744           0 :       if (dropa && dropb)
    1745           0 :         Hp = zero_Flx(sX);
    1746             :       else {
    1747           0 :         if (dropa)
    1748             :         { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    1749           0 :           GEN c = gel(b,degB+2); /* lc(B) */
    1750           0 :           if (odd(degB)) c = Flx_neg(c, p);
    1751           0 :           if (!Flx_equal1(c)) {
    1752           0 :             c = Flx_powu(c, dropa, p);
    1753           0 :             if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
    1754             :           }
    1755             :         }
    1756           0 :         else if (dropb)
    1757             :         { /* multiply by lc(A)^(deg B - deg b) */
    1758           0 :           ulong c = a[degA+2]; /* lc(A) */
    1759           0 :           c = Fl_powu(c, dropb, p);
    1760           0 :           if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
    1761             :         }
    1762             :       }
    1763             :     }
    1764         390 :     if (!H && degpol(Hp) != dres) continue;
    1765         390 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1766         390 :     if (checksqfree) {
    1767         189 :       if (!Flx_is_squarefree(Hp, p)) goto INIT;
    1768         189 :       if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    1769         189 :       checksqfree = 0;
    1770             :     }
    1771             : 
    1772         390 :     if (!H)
    1773             :     { /* initialize */
    1774         189 :       q = utoipos(p); stable = 0;
    1775         189 :       H = ZX_init_CRT(Hp, p,vX);
    1776         189 :       if (LERS) {
    1777         189 :         H0= ZX_init_CRT(H0p, p,vX);
    1778         189 :         H1= ZX_init_CRT(H1p, p,vX);
    1779             :       }
    1780             :     }
    1781             :     else
    1782             :     {
    1783         201 :       if (LERS) {
    1784         201 :         GEN qp = muliu(q,p);
    1785         402 :         stable  = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
    1786         201 :                 & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
    1787         201 :                 & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
    1788         201 :         q = qp;
    1789             :       }
    1790             :       else
    1791           0 :         stable = ZX_incremental_CRT(&H, Hp, &q, p);
    1792             :     }
    1793             :     /* could make it probabilistic for H ? [e.g if stable twice, etc]
    1794             :      * Probabilistic anyway for H0, H1 */
    1795         390 :     if (DEBUGLEVEL>5 && (stable ||  ++cnt==100))
    1796           0 :     { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
    1797         390 :     if (stable && (ulong)expi(q) >= bound) break; /* DONE */
    1798         201 :     if (gc_needed(av,2))
    1799             :     {
    1800           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
    1801           0 :       gerepileall(av2, LERS? 4: 2, &H, &q, &H0, &H1);
    1802             :     }
    1803         201 :   }
    1804             : END:
    1805        1582 :   if (DEBUGLEVEL>5) err_printf(" done\n");
    1806        1582 :   setvarn(H, vX); (void)delete_var();
    1807        1582 :   if (plambda) *plambda = lambda;
    1808        1582 :   if (LERS)
    1809             :   {
    1810        1582 :     *LERS = mkvec2(H0,H1);
    1811        1582 :     gerepileall(av, 2, &H, LERS);
    1812        1582 :     return H;
    1813             :   }
    1814           0 :   return gerepilecopy(av, H);
    1815             : }
    1816             : 
    1817             : GEN
    1818        2590 : ZX_ZXY_resultant_all(GEN A, GEN B, long *plambda, GEN *LERS)
    1819             : {
    1820        2590 :   if (LERS)
    1821        1582 :     return ZX_ZXY_resultant_LERS(A, B, plambda, LERS);
    1822        1008 :   return ZX_ZXY_rnfequation(A, B, plambda);
    1823             : }
    1824             : 
    1825             : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
    1826             :  * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
    1827             :  * squarefree */
    1828             : GEN
    1829        1862 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
    1830             : {
    1831        1862 :   pari_sp av = avma;
    1832             :   GEN R, a;
    1833             :   long dA;
    1834             :   int delvar;
    1835             : 
    1836        1862 :   if (v < 0) v = 0;
    1837        1862 :   switch (typ(A))
    1838             :   {
    1839        1862 :     case t_POL: dA = degpol(A); if (dA > 0) break;
    1840           0 :       A = constant_coeff(A);
    1841             :     default:
    1842           0 :       if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
    1843           0 :       return gerepileupto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
    1844             :   }
    1845        1862 :   delvar = 0;
    1846        1862 :   if (varn(T) == 0)
    1847             :   {
    1848        1806 :     long v0 = fetch_var(); delvar = 1;
    1849        1806 :     T = leafcopy(T); setvarn(T,v0);
    1850        1806 :     A = leafcopy(A); setvarn(A,v0);
    1851             :   }
    1852        1862 :   R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
    1853        1862 :   if (delvar) (void)delete_var();
    1854        1862 :   setvarn(R, v); a = leading_coeff(T);
    1855        1862 :   if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
    1856        1862 :   return gerepileupto(av, R);
    1857             : }
    1858             : 
    1859             : 
    1860             : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
    1861             : GEN
    1862       11479 : ZXQ_charpoly(GEN A, GEN T, long v)
    1863             : {
    1864       11479 :   return (degpol(T) < 16) ? RgXQ_charpoly(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
    1865             : }
    1866             : 
    1867             : GEN
    1868         819 : QXQ_charpoly(GEN A, GEN T, long v)
    1869             : {
    1870         819 :   pari_sp av = avma;
    1871         819 :   GEN den, B = Q_remove_denom(A, &den);
    1872         819 :   GEN P = ZXQ_charpoly(B, T, v);
    1873         819 :   return gerepilecopy(av, den ? RgX_rescale(P, ginv(den)): P);
    1874             : }
    1875             : 
    1876             : static GEN
    1877      155555 : trivial_case(GEN A, GEN B)
    1878             : {
    1879             :   long d;
    1880      155555 :   if (typ(A) == t_INT) return powiu(A, degpol(B));
    1881      148078 :   d = degpol(A);
    1882      148078 :   if (d == 0) return trivial_case(gel(A,2),B);
    1883      145011 :   if (d < 0) return gen_0;
    1884      144996 :   return NULL;
    1885             : }
    1886             : 
    1887             : static ulong
    1888     1287642 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
    1889             : {
    1890     1287642 :   pari_sp av = avma;
    1891             :   ulong H;
    1892             :   long dropa, dropb;
    1893     1287642 :   ulong dp = dB ? umodiu(dB, p): 1;
    1894     1287686 :   if (!b) b = Flx_deriv(a, p);
    1895     1287622 :   dropa = degA - degpol(a);
    1896     1287640 :   dropb = degB - degpol(b);
    1897     1287621 :   if (dropa && dropb) /* p | lc(A), p | lc(B) */
    1898           0 :   { avma = av; return 0; }
    1899     1287621 :   H = Flx_resultant(a, b, p);
    1900     1287467 :   if (dropa)
    1901             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    1902           0 :     ulong c = b[degB+2]; /* lc(B) */
    1903           0 :     if (odd(degB)) c = p - c;
    1904           0 :     c = Fl_powu(c, dropa, p);
    1905           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1906             :   }
    1907     1287467 :   else if (dropb)
    1908             :   { /* multiply by lc(A)^(deg B - deg b) */
    1909           0 :     ulong c = a[degA+2]; /* lc(A) */
    1910           0 :     c = Fl_powu(c, dropb, p);
    1911           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1912             :   }
    1913     1287455 :   if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1914     1287458 :   avma = av; return H;
    1915             : }
    1916             : 
    1917             : /* If B=NULL, assume B=A' */
    1918             : static GEN
    1919      537036 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
    1920             : {
    1921      537036 :   pari_sp av = avma;
    1922      537036 :   long degA, degB, i, n = lg(P)-1;
    1923             :   GEN H, T;
    1924             : 
    1925      537036 :   degA = degpol(A);
    1926      537030 :   degB = B ? degpol(B): degA - 1;
    1927      537066 :   if (n == 1)
    1928             :   {
    1929      158263 :     ulong Hp, p = uel(P,1);
    1930             :     GEN a, b;
    1931      158263 :     a = ZX_to_Flx(A, p), b = B ? ZX_to_Flx(B, p): NULL;
    1932      158258 :     Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1933      158266 :     avma = av;
    1934      158266 :     *mod = utoi(p); return utoi(Hp);
    1935             :   }
    1936      378803 :   T = ZV_producttree(P);
    1937      378801 :   A = ZX_nv_mod_tree(A, P, T);
    1938      378790 :   if (B) B = ZX_nv_mod_tree(B, P, T);
    1939      378791 :   H = cgetg(n+1, t_VECSMALL);
    1940     1508113 :   for(i=1; i <= n; i++)
    1941             :   {
    1942     1129373 :     ulong p = P[i];
    1943     1129373 :     GEN a = gel(A,i), b = B? gel(B,i): NULL;
    1944     1129373 :     H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1945             :   }
    1946      378740 :   H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
    1947      378779 :   *mod = gmael(T, lg(T)-1, 1);
    1948      378779 :   gerepileall(av, 2, &H, mod);
    1949      378790 :   return H;
    1950             : }
    1951             : 
    1952             : GEN
    1953      537054 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
    1954             : {
    1955      537054 :   GEN V = cgetg(3, t_VEC);
    1956      537089 :   if (isintzero(B)) B = NULL;
    1957      537063 :   if (isintzero(dB)) dB = NULL;
    1958      537051 :   gel(V,1) = ZX_resultant_slice(A,B,dB,P,&gel(V,2));
    1959      536960 :   return V;
    1960             : }
    1961             : 
    1962             : /* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */
    1963             : /* if B=NULL, take B = A' */
    1964             : GEN
    1965       79415 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
    1966             : {
    1967       79415 :   pari_sp av = avma;
    1968             :   long m;
    1969             :   GEN  H, worker;
    1970       79415 :   int is_disc = !B;
    1971       79415 :   if (is_disc) B = ZX_deriv(A);
    1972       79415 :   if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H;
    1973       71923 :   if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
    1974       71923 :   if (is_disc)
    1975       47700 :     B = NULL;
    1976       71923 :   worker = strtoclosure("_ZX_resultant_worker", 3, A, B?B:gen_0, dB?dB:gen_0);
    1977       71923 :   m = degpol(A)+(B ? degpol(B): 0);
    1978       71923 :   H = gen_crt("ZX_resultant_all", worker, dB, bound, m, NULL,
    1979             :                ZV_chinese_center, Fp_center);
    1980       71923 :   return gerepileuptoint(av, H);
    1981             : }
    1982             : 
    1983             : /* A0 and B0 in Q[X] */
    1984             : GEN
    1985       10530 : QX_resultant(GEN A0, GEN B0)
    1986             : {
    1987             :   GEN s, a, b, A, B;
    1988       10530 :   pari_sp av = avma;
    1989             : 
    1990       10530 :   A = Q_primitive_part(A0, &a);
    1991       10530 :   B = Q_primitive_part(B0, &b);
    1992       10530 :   s = ZX_resultant(A, B);
    1993       10530 :   if (!signe(s)) { avma = av; return gen_0; }
    1994       10530 :   if (a) s = gmul(s, gpowgs(a,degpol(B)));
    1995       10530 :   if (b) s = gmul(s, gpowgs(b,degpol(A)));
    1996       10530 :   return gerepileupto(av, s);
    1997             : }
    1998             : 
    1999             : GEN
    2000       31070 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
    2001             : 
    2002             : GEN
    2003           0 : QXQ_intnorm(GEN A, GEN B)
    2004             : {
    2005             :   GEN c, n, R, lB;
    2006           0 :   long dA = degpol(A), dB = degpol(B);
    2007           0 :   pari_sp av = avma;
    2008           0 :   if (dA < 0) return gen_0;
    2009           0 :   A = Q_primitive_part(A, &c);
    2010           0 :   if (!c || typ(c) == t_INT) {
    2011           0 :     n = c;
    2012           0 :     R = ZX_resultant(B, A);
    2013             :   } else {
    2014           0 :     n = gel(c,1);
    2015           0 :     R = ZX_resultant_all(B, A, gel(c,2), 0);
    2016             :   }
    2017           0 :   if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
    2018           0 :   lB = leading_coeff(B);
    2019           0 :   if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
    2020           0 :   return gerepileuptoint(av, R);
    2021             : }
    2022             : 
    2023             : GEN
    2024           0 : QXQ_norm(GEN A, GEN B)
    2025             : {
    2026             :   GEN c, R, lB;
    2027           0 :   long dA = degpol(A), dB = degpol(B);
    2028           0 :   pari_sp av = avma;
    2029           0 :   if (dA < 0) return gen_0;
    2030           0 :   A = Q_primitive_part(A, &c);
    2031           0 :   R = ZX_resultant(B, A);
    2032           0 :   if (c) R = gmul(R, gpowgs(c, dB));
    2033           0 :   lB = leading_coeff(B);
    2034           0 :   if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
    2035           0 :   return gerepileupto(av, R);
    2036             : }
    2037             : 
    2038             : /* assume x has integral coefficients */
    2039             : GEN
    2040       49037 : ZX_disc_all(GEN x, ulong bound)
    2041             : {
    2042       49037 :   pari_sp av = avma;
    2043             :   GEN l, R;
    2044       49037 :   long s, d = degpol(x);
    2045       49037 :   if (d <= 1) return d ? gen_1: gen_0;
    2046       47700 :   s = (d & 2) ? -1: 1;
    2047       47700 :   l = leading_coeff(x);
    2048       47700 :   R = ZX_resultant_all(x, NULL, NULL, bound);
    2049       47700 :   if (is_pm1(l))
    2050       44893 :   { if (signe(l) < 0) s = -s; }
    2051             :   else
    2052        2807 :     R = diviiexact(R,l);
    2053       47700 :   if (s == -1) togglesign_safe(&R);
    2054       47700 :   return gerepileuptoint(av,R);
    2055             : }
    2056       48008 : GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); }
    2057             : 
    2058             : GEN
    2059           0 : QX_disc(GEN x)
    2060             : {
    2061           0 :   pari_sp av = avma;
    2062           0 :   GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
    2063           0 :   if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
    2064           0 :   return gerepileupto(av, d);
    2065             : }
    2066             : 
    2067             : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
    2068             : GEN
    2069       24049 : QXQ_inv(GEN A, GEN B)
    2070             : {
    2071             :   GEN D, cU, q, U, V;
    2072             :   ulong p;
    2073       24049 :   pari_sp av2, av = avma;
    2074             :   forprime_t S;
    2075             :   pari_timer ti;
    2076       24049 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2077             :   /* A a QX, B a ZX */
    2078       24049 :   A = Q_primitive_part(A, &D);
    2079             :   /* A, B in Z[X] */
    2080       24049 :   init_modular_small(&S);
    2081       24049 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2082       24049 :   av2 = avma; U = NULL;
    2083      133185 :   while ((p = u_forprime_next(&S)))
    2084             :   {
    2085             :     GEN a, b, qp, Up, Vp;
    2086             :     int stable;
    2087             : 
    2088      109136 :     a = ZX_to_Flx(A, p);
    2089      109136 :     b = ZX_to_Flx(B, p);
    2090             :     /* if p | Res(A/G, B/G), discard */
    2091      133178 :     if (!Flx_extresultant(b,a,p, &Vp,&Up)) continue;
    2092             : 
    2093      109129 :     if (!U)
    2094             :     { /* First time */
    2095       24042 :       U = ZX_init_CRT(Up,p,varn(A));
    2096       24042 :       V = ZX_init_CRT(Vp,p,varn(A));
    2097       24042 :       q = utoipos(p); continue;
    2098             :     }
    2099       85087 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: mod %ld (bound 2^%ld)", p,expi(q));
    2100       85087 :     qp = muliu(q,p);
    2101      170174 :     stable = ZX_incremental_CRT_raw(&U, Up, q,qp, p)
    2102       85087 :            & ZX_incremental_CRT_raw(&V, Vp, q,qp, p);
    2103       85087 :     if (stable)
    2104             :     { /* all stable: check divisibility */
    2105       24042 :       GEN res = ZX_add(ZX_mul(A,U), ZX_mul(B,V));
    2106       24042 :       if (degpol(res) == 0) {
    2107       24042 :         res = gel(res,2);
    2108       24042 :         D = D? gmul(D, res): res;
    2109       48084 :         break;
    2110             :       } /* DONE */
    2111           0 :       if (DEBUGLEVEL) err_printf("QXQ_inv: char 0 check failed");
    2112             :     }
    2113       61045 :     q = qp;
    2114       61045 :     if (gc_needed(av,1))
    2115             :     {
    2116           8 :       if (DEBUGMEM>1) pari_warn(warnmem,"QXQ_inv");
    2117           8 :       gerepileall(av2, 3, &q,&U,&V);
    2118             :     }
    2119             :   }
    2120       24042 :   if (!p) pari_err_OVERFLOW("QXQ_inv [ran out of primes]");
    2121       24042 :   cU = ZX_content(U);
    2122       24042 :   if (!is_pm1(cU)) { U = Q_div_to_int(U, cU); D = gdiv(D, cU); }
    2123       24042 :   return gerepileupto(av, RgX_Rg_div(U, D));
    2124             : }
    2125             : 
    2126             : /************************************************************************
    2127             :  *                                                                      *
    2128             :  *                   ZX_ZXY_resultant                                   *
    2129             :  *                                                                      *
    2130             :  ************************************************************************/
    2131             : 
    2132             : static GEN
    2133       12559 : ZX_ZXY_resultant_prime(GEN a, GEN b, ulong dp, ulong p,
    2134             :                        long degA, long degB, long dres, long sX)
    2135             : {
    2136       12559 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    2137       12560 :   GEN Hp = Flx_FlxY_resultant_polint(a, b, p, dres, sX);
    2138       12560 :   if (dropa && dropb)
    2139           0 :     Hp = zero_Flx(sX);
    2140             :   else {
    2141       12560 :     if (dropa)
    2142             :     { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2143           0 :       GEN c = gel(b,degB+2); /* lc(B) */
    2144           0 :       if (odd(degB)) c = Flx_neg(c, p);
    2145           0 :       if (!Flx_equal1(c)) {
    2146           0 :         c = Flx_powu(c, dropa, p);
    2147           0 :         if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
    2148             :       }
    2149             :     }
    2150       12560 :     else if (dropb)
    2151             :     { /* multiply by lc(A)^(deg B - deg b) */
    2152           0 :       ulong c = uel(a, degA+2); /* lc(A) */
    2153           0 :       c = Fl_powu(c, dropb, p);
    2154           0 :       if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
    2155             :     }
    2156             :   }
    2157       12560 :   if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2158       12560 :   return Hp;
    2159             : }
    2160             : 
    2161             : static GEN
    2162        8989 : ZX_ZXY_resultant_slice(GEN A, GEN B, GEN dB, long degA, long degB, long dres,
    2163             :                        GEN P, GEN *mod, long sX, long vY)
    2164             : {
    2165        8989 :   pari_sp av = avma;
    2166        8989 :   long i, n = lg(P)-1;
    2167             :   GEN H, T, D;
    2168        8989 :   if (n == 1)
    2169             :   {
    2170        8566 :     ulong p = uel(P,1);
    2171        8566 :     ulong dp = dB ? umodiu(dB, p): 1;
    2172        8566 :     GEN a = ZX_to_Flx(A, p), b = ZXX_to_FlxX(B, p, vY);
    2173        8567 :     GEN Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2174        8568 :     H = Flx_to_ZX(Hp);
    2175        8568 :     *mod = utoi(p);
    2176        8567 :     gerepileall(av, 2, &H, mod);
    2177        8564 :     return H;
    2178             :   }
    2179         423 :   T = ZV_producttree(P);
    2180         423 :   A = ZX_nv_mod_tree(A, P, T);
    2181         423 :   B = ZXX_nv_mod_tree(B, P, T, vY);
    2182         423 :   D = dB ? Z_ZV_mod_tree(dB, P, T): NULL;
    2183         423 :   H = cgetg(n+1, t_VEC);
    2184        1489 :   for(i=1; i <= n; i++)
    2185             :   {
    2186        1066 :     ulong p = P[i];
    2187        1066 :     GEN a = gel(A,i), b = gel(B,i);
    2188        1066 :     ulong dp = D ? uel(D, i): 1;
    2189        1066 :     gel(H,i) = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2190             :   }
    2191         423 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2192         423 :   *mod = gmael(T, lg(T)-1, 1);
    2193         423 :   gerepileall(av, 2, &H, mod);
    2194         423 :   return H;
    2195             : }
    2196             : 
    2197             : GEN
    2198        8989 : ZX_ZXY_resultant_worker(GEN P, GEN A, GEN B, GEN dB, GEN v)
    2199             : {
    2200        8989 :   GEN V = cgetg(3, t_VEC);
    2201        8991 :   if (isintzero(dB)) dB = NULL;
    2202        8991 :   gel(V,1) = ZX_ZXY_resultant_slice(A, B, dB, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
    2203        8987 :   return V;
    2204             : }
    2205             : 
    2206             : GEN
    2207        4200 : ZX_ZXY_resultant(GEN A, GEN B)
    2208             : {
    2209        4200 :   pari_sp av = avma;
    2210             :   ulong bound;
    2211        4200 :   long v = fetch_var_higher();
    2212        4200 :   long degA = degpol(A), degB, dres = degA * degpol(B);
    2213        4200 :   long vX = varn(B), vY = varn(A); /* assume vY has lower priority */
    2214        4200 :   long sX = evalvarn(vX);
    2215             :   GEN worker, H, dB;
    2216        4200 :   B = Q_remove_denom(B, &dB);
    2217        4200 :   if (!dB) B = leafcopy(B);
    2218        4200 :   A = leafcopy(A); setvarn(A,v);
    2219        4200 :   B = swap_vars(B, vY); setvarn(B,v); degB = degpol(B);
    2220        4200 :   bound = ZX_ZXY_ResBound(A, B, dB);
    2221        4200 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    2222        4200 :   worker = strtoclosure("_ZX_ZXY_resultant_worker", 4, A, B, dB?dB:gen_0,
    2223             :                         mkvecsmall5(degA, degB,dres, vY, sX));
    2224        4200 :   H = gen_crt("ZX_ZXY_resultant_all", worker, dB, bound, degpol(A)+degpol(B), NULL,
    2225             :                nxV_chinese_center, FpX_center);
    2226        4200 :   setvarn(H, vX); (void)delete_var();
    2227        4200 :   return gerepilecopy(av, H);
    2228             : }
    2229             : 
    2230             : static long
    2231        2443 : ZX_ZXY_rnfequation_lambda(GEN A, GEN B0, long lambda)
    2232             : {
    2233        2443 :   pari_sp av = avma;
    2234        2443 :   long degA = degpol(A), degB, dres = degA*degpol(B0);
    2235        2443 :   long v = fetch_var_higher();
    2236        2443 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    2237        2443 :   long sX = evalvarn(vX);
    2238             :   GEN dB, B, a, b, Hp;
    2239             :   forprime_t S;
    2240             : 
    2241        2443 :   B0 = Q_remove_denom(B0, &dB);
    2242        2443 :   if (!dB) B0 = leafcopy(B0);
    2243        2443 :   A = leafcopy(A);
    2244        2443 :   B = B0;
    2245        2443 :   setvarn(A,v);
    2246             : INIT:
    2247        2926 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    2248        2926 :   B = swap_vars(B, vY); setvarn(B,v);
    2249             :   /* B0(lambda v + x, v) */
    2250        2926 :   if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    2251             : 
    2252        2926 :   degB = degpol(B);
    2253        2926 :   init_modular_big(&S);
    2254             :   while (1)
    2255             :   {
    2256        2926 :     ulong p = u_forprime_next(&S);
    2257        2926 :     ulong dp = dB ? umodiu(dB, p): 1;
    2258        2926 :     if (!dp) continue;
    2259        2926 :     a = ZX_to_Flx(A, p);
    2260        2926 :     b = ZXX_to_FlxX(B, p, v);
    2261        2926 :     Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2262        2926 :     if (degpol(Hp) != dres) continue;
    2263        2926 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2264        2926 :     if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); goto INIT; }
    2265        2443 :     if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    2266        4886 :     avma = av; (void)delete_var(); return lambda;
    2267           0 :   }
    2268             : }
    2269             : 
    2270             : GEN
    2271        2996 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
    2272             : {
    2273        2996 :   if (lambda)
    2274             :   {
    2275        2443 :     *lambda = ZX_ZXY_rnfequation_lambda(A, B, *lambda);
    2276        2443 :     B = RgX_translate(B, monomial(stoi(*lambda), 1, varn(A)));
    2277             :   }
    2278        2996 :   return ZX_ZXY_resultant(A,B);
    2279             : }
    2280             : 
    2281             : /************************************************************************
    2282             :  *                                                                      *
    2283             :  *                   IRREDUCIBLE POLYNOMIAL / Fp                        *
    2284             :  *                                                                      *
    2285             :  ************************************************************************/
    2286             : 
    2287             : /* irreducible (unitary) polynomial of degree n over Fp */
    2288             : GEN
    2289           0 : ffinit_rand(GEN p,long n)
    2290             : {
    2291             :   for(;;) {
    2292           0 :     pari_sp av = avma;
    2293           0 :     GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
    2294           0 :     if (FpX_is_irred(pol, p)) return pol;
    2295           0 :     avma = av;
    2296           0 :   }
    2297             : }
    2298             : 
    2299             : /* return an extension of degree 2^l of F_2, assume l > 0
    2300             :  * Not stack clean. */
    2301             : static GEN
    2302         596 : f2init(long l)
    2303             : {
    2304             :   GEN Q, T, S;
    2305             :   long i, v;
    2306             : 
    2307         596 :   if (l == 1) return polcyclo(3, 0);
    2308         561 :   v = fetch_var_higher();
    2309         562 :   S = mkpoln(4, gen_1,gen_1,gen_0,gen_0); /* y(y^2 + y) */
    2310         562 :   Q = mkpoln(3, gen_1,gen_1, S); /* x^2 + x + y(y^2+y) */
    2311         562 :   setvarn(Q, v);
    2312             : 
    2313             :   /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
    2314         562 :   T = mkpoln(5, gen_1,gen_0,gen_0,gen_1,gen_1);
    2315         561 :   setvarn(T, v);
    2316             :   /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
    2317             :    * ==> x^2 + x + a(y) b irred. over K for any root b of Q
    2318             :    * ==> x^2 + x + (b^2+b)b */
    2319         561 :   for (i=2; i<l; i++) T = FpX_FpXY_resultant(T, Q, gen_2); /* minpoly(b) / F2*/
    2320         563 :   (void)delete_var(); setvarn(T,0); return T;
    2321             : }
    2322             : 
    2323             : /* return an extension of degree p^l of F_p, assume l > 0
    2324             :  * Not stack clean. */
    2325             : GEN
    2326           0 : ffinit_Artin_Shreier(GEN ip, long l)
    2327             : {
    2328           0 :   long i, v, p = itos(ip);
    2329           0 :   GEN T, Q, xp = pol_xn(p,0); /* x^p */
    2330           0 :   T = ZX_sub(xp, deg1pol_shallow(gen_1,gen_1,0)); /* x^p - x - 1 */
    2331           0 :   if (l == 1) return T;
    2332             : 
    2333           0 :   v = fetch_var_higher();
    2334           0 :   setvarn(xp, v);
    2335           0 :   Q = ZX_sub(pol_xn(2*p-1,0), pol_xn(p,0));
    2336           0 :   Q = gsub(xp, deg1pol_shallow(gen_1, Q, v)); /* x^p - x - (y^(2p-1)-y^p) */
    2337           0 :   for (i = 2; i <= l; ++i) T = FpX_FpXY_resultant(T, Q, ip);
    2338           0 :   (void)delete_var(); setvarn(T,0); return T;
    2339             : }
    2340             : 
    2341             : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
    2342             : static long
    2343       12740 : fpinit_check(GEN p, long n, long l)
    2344             : {
    2345             :   ulong q;
    2346       12740 :   if (!uisprime(n)) return 0;
    2347        6139 :   q = umodiu(p,n); if (!q) return 0;
    2348        5558 :   return cgcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    2349             : }
    2350             : 
    2351             : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
    2352             :  * Return an irreducible polynomial of degree l over F_p.
    2353             :  * Variant of Adleman and Lenstra "Finding irreducible polynomials over
    2354             :  * finite fields", ACM, 1986 (5) 350--355.
    2355             :  * Not stack clean */
    2356             : static GEN
    2357        3206 : fpinit(GEN p, long l)
    2358             : {
    2359        3206 :   ulong n = 1+l;
    2360        3206 :   while (!fpinit_check(p,n,l)) n += l;
    2361        3206 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    2362        3206 :   return FpX_red(polsubcyclo(n,l,0),p);
    2363             : }
    2364             : 
    2365             : static GEN
    2366        3185 : ffinit_fact(GEN p, long n)
    2367             : {
    2368        3185 :   GEN P, F = gel(factoru_pow(n),3);
    2369             :   long i;
    2370        3187 :   if (!odd(n) && absequaliu(p, 2))
    2371         596 :     P = f2init(vals(n)); /* if n is even, F[1] = 2^vals(n)*/
    2372             :   else
    2373        2590 :     P = fpinit(p, F[1]);
    2374        3678 :   for (i = 2; i < lg(F); ++i)
    2375         490 :     P = FpX_direct_compositum(fpinit(p, F[i]), P, p);
    2376        3188 :   return P;
    2377             : }
    2378             : 
    2379             : static GEN
    2380         126 : ffinit_nofact(GEN p, long n)
    2381             : {
    2382         126 :   GEN P, Q = NULL;
    2383         126 :   if (lgefint(p)==3)
    2384             :   {
    2385           0 :     ulong pp = p[2], q;
    2386           0 :     long v = u_lvalrem(n,pp,&q);
    2387           0 :     if (v>0)
    2388             :     {
    2389           0 :       Q = (pp == 2)? f2init(v): fpinit(p,n/q);
    2390           0 :       n = q;
    2391             :     }
    2392             :   }
    2393             :   /* n coprime to p */
    2394         126 :   if (n==1) P = Q;
    2395             :   else
    2396             :   {
    2397         126 :     P = fpinit(p, n);
    2398         126 :     if (Q) P = FpX_direct_compositum(P, Q, p);
    2399             :   }
    2400         126 :   return P;
    2401             : }
    2402             : 
    2403             : static GEN
    2404        4116 : init_Fq_i(GEN p, long n, long v)
    2405             : {
    2406             :   GEN P;
    2407        4116 :   if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
    2408        4117 :   if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
    2409        4117 :   if (signe(p) <= 0) pari_err_PRIME("ffinit",p);
    2410        4117 :   if (v < 0) v = 0;
    2411        4117 :   if (n == 1) return pol_x(v);
    2412        3928 :   if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
    2413        3312 :   if (lgefint(p)-2 <= expu(n))
    2414        3185 :     P = ffinit_fact(p,n);
    2415             :   else
    2416         126 :     P = ffinit_nofact(p,n);
    2417        3314 :   setvarn(P, v); return P;
    2418             : }
    2419             : GEN
    2420        3961 : init_Fq(GEN p, long n, long v)
    2421             : {
    2422        3961 :   pari_sp av = avma;
    2423        3961 :   return gerepileupto(av, init_Fq_i(p, n, v));
    2424             : }
    2425             : GEN
    2426         154 : ffinit(GEN p, long n, long v)
    2427             : {
    2428         154 :   pari_sp av = avma;
    2429         154 :   return gerepileupto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
    2430             : }
    2431             : 
    2432             : GEN
    2433        3178 : ffnbirred(GEN p, long n)
    2434             : {
    2435        3178 :   pari_sp av = avma;
    2436             :   long j;
    2437        3178 :   GEN s = gen_0, dk, pd;
    2438        3178 :   dk = divisorsu(n);
    2439       10535 :   for (j = 1; j < lg(dk); ++j)
    2440             :   {
    2441        7357 :     long d = dk[j];
    2442        7357 :     long m = moebiusu(d);
    2443        7357 :     if (!m) continue;
    2444        6797 :     pd = powiu(p, n/d);
    2445        6797 :     s = m>0 ? addii(s, pd): subii(s,pd);
    2446             :   }
    2447        3178 :   return gerepileuptoint(av, divis(s, n));
    2448             : }
    2449             : 
    2450             : GEN
    2451         434 : ffsumnbirred(GEN p, long n)
    2452             : {
    2453         434 :   pari_sp av = avma;
    2454             :   long i,j;
    2455         434 :   GEN v,q, t = gen_0;
    2456         434 :   v = cgetg(n+1,t_VECSMALL); v[1] = 1;
    2457         434 :   q = cgetg(n+1,t_VEC); gel(q,1) = p;
    2458        1547 :   for(i=2; i<=n; i++)
    2459             :   {
    2460        1113 :     v[i] = moebiusu(i);
    2461        1113 :     gel(q,i) = mulii(gel(q,i-1), p);
    2462             :   }
    2463        1981 :   for(i=1; i<=n; i++)
    2464             :   {
    2465        1547 :     GEN s = gen_0;
    2466        1547 :     GEN dk = divisorsu(i);
    2467        4725 :     for (j = 1; j < lg(dk); ++j)
    2468             :     {
    2469        3178 :       long d = dk[j], m = v[d];
    2470        3178 :       if (!m) continue;
    2471        2884 :       s = m>0 ? addii(s, gel(q, i/d)): subii(s, gel(q, i/d));
    2472             :     }
    2473        1547 :     t = addii(t, divis(s, i));
    2474             :   }
    2475         434 :   return gerepileuptoint(av, t);
    2476             : }
    2477             : 
    2478             : GEN
    2479         140 : ffnbirred0(GEN p, long n, long flag)
    2480             : {
    2481         140 :   if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
    2482         140 :   if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
    2483         140 :   switch(flag)
    2484             :   {
    2485             :     case 0:
    2486          70 :       return ffnbirred(p, n);
    2487             :     case 1:
    2488          70 :       return ffsumnbirred(p, n);
    2489             :     default:
    2490           0 :       pari_err_FLAG("ffnbirred");
    2491             :   }
    2492             :   return NULL; /* LCOV_EXCL_LINE */
    2493             : }
    2494             : 
    2495             : static void
    2496         994 : checkmap(GEN m, const char *s)
    2497             : {
    2498         994 :   if (typ(m)!=t_VEC || lg(m)!=3 || typ(gel(m,1))!=t_FFELT)
    2499           0 :     pari_err_TYPE(s,m);
    2500         994 : }
    2501             : 
    2502             : GEN
    2503          84 : ffembed(GEN a, GEN b)
    2504             : {
    2505          84 :   pari_sp av = avma;
    2506          84 :   GEN p, Ta, Tb, g, r = NULL;
    2507          84 :   if (typ(a)!=t_FFELT)
    2508           0 :     pari_err_TYPE("ffembed",a);
    2509          84 :   if (typ(b)!=t_FFELT)
    2510           0 :     pari_err_TYPE("ffembed",b);
    2511          84 :   p = FF_p_i(a); g = FF_gen(a);
    2512          84 :   if (!equalii(p, FF_p_i(b)))
    2513           0 :     pari_err_MODULUS("ffembed",a,b);
    2514          84 :   Ta = FF_mod(a);
    2515          84 :   Tb = FF_mod(b);
    2516          84 :   if (degpol(Tb)%degpol(Ta)!=0)
    2517           0 :     pari_err_DOMAIN("ffembed","a","in a subfield of",b,a);
    2518          84 :   r = gel(FFX_roots(Ta, b), 1);
    2519          84 :   return gerepilecopy(av, mkvec2(g,r));
    2520             : }
    2521             : 
    2522             : GEN
    2523          42 : ffextend(GEN a, GEN P, long v)
    2524             : {
    2525          42 :   pari_sp av = avma;
    2526             :   long n;
    2527             :   GEN p, T, R, g, m;
    2528          42 :   if (typ(a)!=t_FFELT)
    2529           0 :     pari_err_TYPE("ffextend",a);
    2530          42 :   T = a; p = FF_p_i(a);
    2531          42 :   if (typ(P)!=t_POL || !RgX_is_FpXQX(P,&T,&p))
    2532          21 :     pari_err_TYPE("ffextend", P);
    2533          21 :   if (!FF_samefield(a, T))
    2534           0 :     pari_err_MODULUS("ffextend",a,T);
    2535          21 :   if (v < 0) v = varn(P);
    2536          21 :   n = FF_f(T) * degpol(P); R = ffinit(p, n, v); g = ffgen(R, v);
    2537          21 :   m = ffembed(a, g);
    2538          21 :   R = FFX_roots(ffmap(m, P),g);
    2539          21 :   return gerepilecopy(av, mkvec2(gel(R,1), m));
    2540             : }
    2541             : 
    2542             : GEN
    2543          21 : fffrobenius(GEN a, long n)
    2544             : {
    2545             :   GEN g;
    2546          21 :   if (typ(a)!=t_FFELT) pari_err_TYPE("fffrobenius",a);
    2547          21 :   retmkvec2(g=FF_gen(a), FF_Frobenius(g, n));
    2548             : }
    2549             : 
    2550             : GEN
    2551          63 : ffinvmap(GEN m)
    2552             : {
    2553          63 :   pari_sp av = avma;
    2554             :   long i, l;
    2555          63 :   GEN T, F, a, g, r, f = NULL;
    2556          63 :   checkmap(m, "ffinvmap");
    2557          63 :   a = gel(m,1); r = gel(m,2);
    2558          63 :   g = FF_gen(a);
    2559          63 :   T = FF_mod(r);
    2560          63 :   F = gel(FFX_factor(T, a), 1);
    2561          63 :   l = lg(F);
    2562         406 :   for(i=1; i<l; i++)
    2563             :   {
    2564         406 :     GEN s = FFX_rem(FF_to_FpXQ_i(r), gel(F, i), a);
    2565         406 :     if (degpol(s)==0 && gequal(constant_term(s),g))
    2566          63 :       { f = gel(F, i); break; }
    2567             :   }
    2568          63 :   if (f==NULL) pari_err_TYPE("ffinvmap", m);
    2569          63 :   if (degpol(f)==1) f = FF_neg_i(gel(f,2));
    2570          63 :   return gerepilecopy(av, mkvec2(FF_gen(r),f));
    2571             : }
    2572             : 
    2573             : static GEN
    2574         567 : ffpartmapimage(const char *s, GEN r)
    2575             : {
    2576         567 :    GEN a = NULL, p = NULL;
    2577         567 :    if (typ(r)==t_POL && degpol(r) >= 1
    2578         567 :       && RgX_is_FpXQX(r,&a,&p) && a && typ(a)==t_FFELT)
    2579         567 :      return a;
    2580           0 :    pari_err_TYPE(s, r);
    2581             :    return NULL; /* LCOV_EXCL_LINE */
    2582             : }
    2583             : 
    2584             : static GEN
    2585        1393 : ffeltmap_i(GEN m, GEN x)
    2586             : {
    2587        1393 :    GEN r = gel(m,2);
    2588        1393 :    if (!FF_samefield(x, gel(m,1)))
    2589          42 :      pari_err_DOMAIN("ffmap","m","domain do not contain", x, r);
    2590        1351 :    if (typ(r)==t_FFELT)
    2591         847 :      return FF_map(r, x);
    2592             :    else
    2593         504 :      return FFX_preimage(x, r, ffpartmapimage("ffmap", r));
    2594             : }
    2595             : 
    2596             : static GEN
    2597        2212 : ffmap_i(GEN m, GEN x)
    2598             : {
    2599             :   GEN y;
    2600        2212 :   long i, lx, tx = typ(x);
    2601        2212 :   switch(tx)
    2602             :   {
    2603             :     case t_FFELT:
    2604        1309 :       return ffeltmap_i(m, x);
    2605             :     case t_POL: case t_RFRAC: case t_SER:
    2606             :     case t_VEC: case t_COL: case t_MAT:
    2607         630 :       y = cgetg_copy(x, &lx);
    2608         630 :       for (i=1; i<lontyp[tx]; i++) y[i] = x[1];
    2609        2268 :       for (i=lontyp[tx]; i<lx; i++)
    2610             :       {
    2611        1659 :         GEN yi = ffmap_i(m, gel(x,i));
    2612        1638 :         if (!yi) return NULL;
    2613        1638 :         gel(y,i) = yi;
    2614             :       }
    2615         609 :       return y;
    2616             :   }
    2617         273 :   return gcopy(x);
    2618             : }
    2619             : 
    2620             : GEN
    2621         511 : ffmap(GEN m, GEN x)
    2622             : {
    2623         511 :   pari_sp ltop = avma;
    2624             :   GEN y;
    2625         511 :   checkmap(m, "ffmap");
    2626         511 :   y = ffmap_i(m, x);
    2627         511 :   if (y) return y;
    2628          21 :   avma = ltop; return cgetg(1,t_VEC);
    2629             : }
    2630             : 
    2631             : static void
    2632          42 : err_compo(GEN m, GEN n)
    2633             : {
    2634          42 :   pari_err_DOMAIN("ffcompomap","m","domain do not contain codomain of",n,m);
    2635           0 : }
    2636             : 
    2637             : GEN
    2638         210 : ffcompomap(GEN m, GEN n)
    2639             : {
    2640         210 :   pari_sp av = avma;
    2641         210 :   GEN g = gel(n,1), r, m2, n2;
    2642         210 :   checkmap(m, "ffcompomap");
    2643         210 :   checkmap(n, "ffcompomap");
    2644         210 :   m2 = gel(m,2); n2 = gel(n,2);
    2645         210 :   switch((typ(m2)==t_POL)|((typ(n2)==t_POL)<<1))
    2646             :   {
    2647             :     case 0:
    2648          42 :       if (!FF_samefield(gel(m,1),n2)) err_compo(m,n);
    2649          21 :       r = FF_map(gel(m,2), n2);
    2650          21 :       break;
    2651             :     case 2:
    2652          42 :       r = ffmap_i(m, n2);
    2653          21 :       if (lg(r) == 1) err_compo(m,n);
    2654          21 :       break;
    2655             :     case 1:
    2656          84 :       r = ffeltmap_i(m, n2);
    2657          63 :       if (!r)
    2658             :       {
    2659             :         GEN a, A, R, M;
    2660             :         long dm, dn;
    2661          21 :         a = ffpartmapimage("ffcompomap",m2);
    2662          21 :         A = FF_to_FpXQ_i(FF_neg(n2));
    2663          21 :         setvarn(A, 1);
    2664          21 :         R = deg1pol(gen_1, A, 0);
    2665          21 :         setvarn(R, 0);
    2666          21 :         M = gcopy(m2);
    2667          21 :         setvarn(M, 1);
    2668          21 :         r = polresultant0(R, M, 1, 0);
    2669          21 :         dm = FF_f(gel(m,1)); dn = FF_f(gel(n,1));
    2670          21 :         if (dm % dn || !FFX_ispower(r, dm/dn, a, &r)) err_compo(m,n);
    2671          21 :         setvarn(r, varn(FF_mod(g)));
    2672             :       }
    2673          63 :       break;
    2674             :     case 3:
    2675             :     {
    2676             :       GEN M, R, T, p, a;
    2677          42 :       a = ffpartmapimage("ffcompomap",n2);
    2678          42 :       if (!FF_samefield(a, gel(m,1))) err_compo(m,n);
    2679          21 :       p = FF_p_i(gel(n,1));
    2680          21 :       T = FF_mod(gel(n,1));
    2681          21 :       setvarn(T, 1);
    2682          21 :       R = RgX_to_FpXQX(n2,T,p);
    2683          21 :       setvarn(R, 0);
    2684          21 :       M = gcopy(m2);
    2685          21 :       setvarn(M, 1);
    2686          21 :       r = polresultant0(R, M, 1, 0);
    2687          21 :       setvarn(r, varn(n2));
    2688             :     }
    2689             :   }
    2690         126 :   return gerepilecopy(av, mkvec2(g,r));
    2691             : }

Generated by: LCOV version 1.11