Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - polarit3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.13.0 lcov report (development 25825-edc109b529) Lines: 1722 1930 89.2 %
Date: 2020-09-21 06:08:33 Functions: 180 195 92.3 %
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        1239 : char_update_prime(struct charact *S, GEN p)
      34             : {
      35        1239 :   if (!S->isprime) { S->isprime = 1; S->q = p; }
      36        1239 :   if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
      37        1232 : }
      38             : static void
      39        5194 : char_update_int(struct charact *S, GEN n)
      40             : {
      41        5194 :   if (S->isprime)
      42             :   {
      43           7 :     if (dvdii(n, S->q)) return;
      44           7 :     pari_err_MODULUS("characteristic", S->q, n);
      45             :   }
      46        5187 :   S->q = gcdii(S->q, n);
      47             : }
      48             : static void
      49      596463 : charact(struct charact *S, GEN x)
      50             : {
      51      596463 :   const long tx = typ(x);
      52             :   long i, l;
      53      596463 :   switch(tx)
      54             :   {
      55        4221 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      56        1148 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      57       15722 :     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       15722 :       l = lg(x);
      61      601468 :       for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
      62       15708 :       break;
      63           7 :     case t_LIST:
      64           7 :       x = list_data(x);
      65           7 :       if (x) charact(S, x);
      66           7 :       break;
      67             :   }
      68      596435 : }
      69             : static void
      70       33315 : charact_res(struct charact *S, GEN x)
      71             : {
      72       33315 :   const long tx = typ(x);
      73             :   long i, l;
      74       33315 :   switch(tx)
      75             :   {
      76         973 :     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       10380 :     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       10380 :       l = lg(x);
      83       41184 :       for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
      84       10380 :       break;
      85           0 :     case t_LIST:
      86           0 :       x = list_data(x);
      87           0 :       if (x) charact_res(S, x);
      88           0 :       break;
      89             :   }
      90       33315 : }
      91             : GEN
      92       10703 : characteristic(GEN x)
      93             : {
      94             :   struct charact S;
      95       10703 :   S.q = gen_0; S.isprime = 0;
      96       10703 :   charact(&S, x); return S.q;
      97             : }
      98             : GEN
      99        2511 : residual_characteristic(GEN x)
     100             : {
     101             :   struct charact S;
     102        2511 :   S.q = gen_0; S.isprime = 0;
     103        2511 :   charact_res(&S, x); return S.q;
     104             : }
     105             : 
     106             : int
     107    60821536 : Rg_is_Fp(GEN x, GEN *pp)
     108             : {
     109             :   GEN mod;
     110    60821536 :   switch(typ(x))
     111             :   {
     112     4256875 :   case t_INTMOD:
     113     4256875 :     mod = gel(x,1);
     114     4256875 :     if (!*pp) *pp = mod;
     115     4021955 :     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     4256875 :     return 1;
     121    52328704 :   case t_INT:
     122    52328704 :     return 1;
     123     4235957 :   default: return 0;
     124             :   }
     125             : }
     126             : 
     127             : int
     128    20550159 : RgX_is_FpX(GEN x, GEN *pp)
     129             : {
     130    20550159 :   long i, lx = lg(x);
     131    77109943 :   for (i=2; i<lx; i++)
     132    60795741 :     if (!Rg_is_Fp(gel(x, i), pp))
     133     4235957 :       return 0;
     134    16314202 :   return 1;
     135             : }
     136             : 
     137             : int
     138           0 : RgV_is_FpV(GEN x, GEN *pp)
     139             : {
     140           0 :   long i, lx = lg(x);
     141           0 :   for (i=1; i<lx; i++)
     142           0 :     if (!Rg_is_Fp(gel(x,i), pp)) return 0;
     143           0 :   return 1;
     144             : }
     145             : 
     146             : int
     147           0 : RgM_is_FpM(GEN x, GEN *pp)
     148             : {
     149           0 :   long i, lx = lg(x);
     150           0 :   for (i=1; i<lx; i++)
     151           0 :     if (!RgV_is_FpV(gel(x, i), pp)) return 0;
     152           0 :   return 1;
     153             : }
     154             : 
     155             : int
     156       58345 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
     157             : {
     158             :   GEN pol, mod, p;
     159       58345 :   switch(typ(x))
     160             :   {
     161       25788 :   case t_INTMOD:
     162       25788 :     return Rg_is_Fp(x, pp);
     163        6538 :   case t_INT:
     164        6538 :     return 1;
     165          42 :   case t_POL:
     166          42 :     return RgX_is_FpX(x, pp);
     167       21350 :   case t_FFELT:
     168       21350 :     mod = x; p = FF_p_i(x);
     169       21350 :     if (!*pp) *pp = p;
     170       21350 :     if (!*pT) *pT = mod;
     171       19824 :     else if (typ(*pT)!=t_FFELT || !FF_samefield(*pT,mod))
     172             :     {
     173          42 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     174          42 :       return 0;
     175             :     }
     176       21308 :     return 1;
     177        4543 :   case t_POLMOD:
     178        4543 :     mod = gel(x,1); pol = gel(x, 2);
     179        4543 :     if (!RgX_is_FpX(mod, pp)) return 0;
     180        4543 :     if (typ(pol)==t_POL)
     181             :     {
     182        4536 :       if (!RgX_is_FpX(pol, pp)) return 0;
     183             :     }
     184           7 :     else if (!Rg_is_Fp(pol, pp)) return 0;
     185        4543 :     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        4543 :     return 1;
     192             : 
     193          84 :   default: return 0;
     194             :   }
     195             : }
     196             : 
     197             : int
     198        3045 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
     199             : {
     200        3045 :   long i, lx = lg(x);
     201       60837 :   for (i = 2; i < lx; i++)
     202       57834 :     if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
     203        3003 :   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    32995458 : Rg_to_Fp(GEN x, GEN p)
     216             : {
     217    32995458 :   if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
     218     3073226 :   switch(typ(x))
     219             :   {
     220      179955 :     case t_INT: return modii(x, p);
     221         121 :     case t_FRAC: {
     222         121 :       pari_sp av = avma;
     223         121 :       GEN z = modii(gel(x,1), p);
     224         121 :       if (z == gen_0) return gen_0;
     225         121 :       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     2893164 :     case t_INTMOD: {
     229     2893164 :       GEN q = gel(x,1), a = gel(x,2);
     230     2893164 :       if (equalii(q, p)) return icopy(a);
     231          12 :       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     1292383 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
     241             : {
     242     1292383 :   long ta, tx = typ(x), v = get_FpX_var(T);
     243             :   GEN a, b;
     244     1292383 :   if (is_const_t(tx))
     245             :   {
     246       57913 :     if (tx == t_FFELT)
     247             :     {
     248       17085 :       GEN z = FF_to_FpXQ(x);
     249       17085 :       setvarn(z, v);
     250       17085 :       return z;
     251             :     }
     252       40828 :     return scalar_ZX(degpol(T)? Rg_to_Fp(x, p): gen_0, v);
     253             :   }
     254     1234470 :   switch(tx)
     255             :   {
     256     1229619 :     case t_POLMOD:
     257     1229619 :       b = gel(x,1);
     258     1229619 :       a = gel(x,2); ta = typ(a);
     259     1229619 :       if (is_const_t(ta))
     260        4025 :         return scalar_ZX(degpol(T)? Rg_to_Fp(a, p): gen_0, v);
     261     1225594 :       b = RgX_to_FpX(b, p); if (varn(b) != v) break;
     262     1225594 :       a = RgX_to_FpX(a, p);
     263     1225594 :       if (ZX_equal(b,get_FpX_mod(T)) || signe(FpX_rem(b,T,p))==0)
     264     1225594 :         return FpX_rem(a, T, p);
     265           0 :       break;
     266        4851 :     case t_POL:
     267        4851 :       if (varn(x) != v) break;
     268        4851 :       return FpX_rem(RgX_to_FpX(x,p), T, p);
     269           0 :     case t_RFRAC:
     270           0 :       a = Rg_to_FpXQ(gel(x,1), T,p);
     271           0 :       b = Rg_to_FpXQ(gel(x,2), T,p);
     272           0 :       return FpXQ_div(a,b, T,p);
     273             :   }
     274           0 :   pari_err_TYPE("Rg_to_FpXQ",x);
     275             :   return NULL; /* LCOV_EXCL_LINE */
     276             : }
     277             : GEN
     278     3432517 : RgX_to_FpX(GEN x, GEN p)
     279             : {
     280             :   long i, l;
     281     3432517 :   GEN z = cgetg_copy(x, &l); z[1] = x[1];
     282    15101997 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     283     3432517 :   return FpX_renormalize(z, l);
     284             : }
     285             : 
     286             : GEN
     287         126 : RgV_to_FpV(GEN x, GEN p)
     288         385 : { pari_APPLY_type(t_VEC, Rg_to_Fp(gel(x,i), p)) }
     289             : 
     290             : GEN
     291      925732 : RgC_to_FpC(GEN x, GEN p)
     292    11881302 : { pari_APPLY_type(t_COL, Rg_to_Fp(gel(x,i), p)) }
     293             : 
     294             : GEN
     295      129453 : RgM_to_FpM(GEN x, GEN p)
     296     1055143 : { pari_APPLY_same(RgC_to_FpC(gel(x,i), p)) }
     297             : 
     298             : GEN
     299      282099 : RgV_to_Flv(GEN x, ulong p)
     300     1715382 : { pari_APPLY_ulong(Rg_to_Fl(gel(x,i), p)) }
     301             : 
     302             : GEN
     303      114250 : RgM_to_Flm(GEN x, ulong p)
     304      395831 : { pari_APPLY_same(RgV_to_Flv(gel(x,i), p)) }
     305             : 
     306             : GEN
     307        5000 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
     308             : {
     309        5000 :   long i, l = lg(x);
     310        5000 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     311       42855 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
     312        5000 :   return FpXQX_renormalize(z, l);
     313             : }
     314             : GEN
     315        3108 : RgX_to_FqX(GEN x, GEN T, GEN p)
     316             : {
     317        3108 :   long i, l = lg(x);
     318        3108 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     319        3108 :   if (T)
     320       10990 :     for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
     321             :   else
     322       58436 :     for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     323        3108 :   return FpXQX_renormalize(z, l);
     324             : }
     325             : 
     326             : GEN
     327      219184 : RgC_to_FqC(GEN x, GEN T, GEN p)
     328             : {
     329      219184 :   long i, l = lg(x);
     330      219184 :   GEN z = cgetg(l, t_COL);
     331      219184 :   if (T)
     332     1431710 :     for (i = 1; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
     333             :   else
     334           0 :     for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     335      219184 :   return z;
     336             : }
     337             : 
     338             : GEN
     339       52332 : RgM_to_FqM(GEN x, GEN T, GEN p)
     340      271502 : { pari_APPLY_same(RgC_to_FqC(gel(x, i), T, p)) }
     341             : 
     342             : /* lg(V) > 1 */
     343             : GEN
     344      849765 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
     345             : {
     346      849765 :   pari_sp av = avma;
     347      849765 :   long i, l = lg(V);
     348      849765 :   GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
     349     4181499 :   for(i=2; i<l; i++)
     350             :   {
     351     3331734 :     z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
     352     3331734 :     if ((i & 7) == 0) z = gerepileupto(av, z);
     353             :   }
     354      849765 :   return gerepileupto(av, FpX_red(z,p));
     355             : }
     356             : 
     357             : GEN
     358       52976 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
     359             : {
     360       52976 :   long i, lz = lg(y);
     361             :   GEN z;
     362       52976 :   if (!T) return FpX_Fp_add(y, x, p);
     363        5012 :   if (lz == 2) return scalarpol(x, varn(y));
     364        4599 :   z = cgetg(lz,t_POL); z[1] = y[1];
     365        4599 :   gel(z,2) = Fq_add(gel(y,2),x, T, p);
     366        4599 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     367             :   else
     368       18725 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     369        4599 :   return z;
     370             : }
     371             : 
     372             : GEN
     373        1048 : FqX_Fq_sub(GEN y, GEN x, GEN T, GEN p)
     374             : {
     375        1048 :   long i, lz = lg(y);
     376             :   GEN z;
     377        1048 :   if (!T) return FpX_Fp_sub(y, x, p);
     378        1048 :   if (lz == 2) return scalarpol(x, varn(y));
     379        1048 :   z = cgetg(lz,t_POL); z[1] = y[1];
     380        1048 :   gel(z,2) = Fq_sub(gel(y,2), x, T, p);
     381        1048 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     382             :   else
     383       10288 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     384        1048 :   return z;
     385             : }
     386             : 
     387             : GEN
     388      148818 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
     389             : {
     390             :   long i, lP;
     391      148818 :   GEN res = cgetg_copy(P, &lP); res[1] = P[1];
     392      918380 :   for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
     393      148818 :   gel(res,lP-1) = gen_1; return res;
     394             : }
     395             : 
     396             : GEN
     397        4877 : FpXQX_normalize(GEN z, GEN T, GEN p)
     398             : {
     399             :   GEN lc;
     400        4877 :   if (lg(z) == 2) return z;
     401        4863 :   lc = leading_coeff(z);
     402        4863 :   if (typ(lc) == t_POL)
     403             :   {
     404        1947 :     if (lg(lc) > 3) /* non-constant */
     405        1677 :       return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
     406             :     /* constant */
     407         270 :     lc = gel(lc,2);
     408         270 :     z = shallowcopy(z);
     409         270 :     gel(z, lg(z)-1) = lc;
     410             :   }
     411             :   /* lc a t_INT */
     412        3186 :   if (equali1(lc)) return z;
     413          78 :   return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
     414             : }
     415             : 
     416             : GEN
     417      127820 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
     418             : {
     419             :   pari_sp av;
     420             :   GEN p1, r;
     421      127820 :   long j, i=lg(x)-1;
     422      127820 :   if (i<=2)
     423       26327 :     return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
     424      101493 :   av=avma; p1=gel(x,i);
     425             :   /* specific attention to sparse polynomials (see poleval)*/
     426             :   /*You've guessed it! It's a copy-paste(tm)*/
     427      299894 :   for (i--; i>=2; i=j-1)
     428             :   {
     429      198905 :     for (j=i; !signe(gel(x,j)); j--)
     430         504 :       if (j==2)
     431             :       {
     432         315 :         if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
     433         315 :         return gerepileupto(av, Fq_mul(p1,y, T, p));
     434             :       }
     435      198401 :     r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
     436      198401 :     p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
     437             :   }
     438      101178 :   return gerepileupto(av, p1);
     439             : }
     440             : 
     441             : GEN
     442       31598 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
     443             : {
     444       31598 :   long i, lb = lg(Q);
     445             :   GEN z;
     446       31598 :   if (!T) return FpXY_evalx(Q, x, p);
     447       20993 :   z = cgetg(lb, t_POL); z[1] = Q[1];
     448      117068 :   for (i=2; i<lb; i++)
     449             :   {
     450       96075 :     GEN q = gel(Q,i);
     451       96075 :     gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
     452             :   }
     453       20993 :   return FpXQX_renormalize(z, lb);
     454             : }
     455             : 
     456             : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
     457             : GEN
     458       14392 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
     459             : {
     460       14392 :   pari_sp av = avma;
     461       14392 :   if (!T) return FpXY_eval(Q, y, x, p);
     462         588 :   return gerepileupto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
     463             : }
     464             : 
     465             : /* a X^d */
     466             : GEN
     467     7636941 : monomial(GEN a, long d, long v)
     468             : {
     469             :   long i, n;
     470             :   GEN P;
     471     7636941 :   if (d < 0) {
     472           0 :     if (isrationalzero(a)) return pol_0(v);
     473           0 :     retmkrfrac(a, pol_xn(-d, v));
     474             :   }
     475     7636941 :   if (gequal0(a))
     476             :   {
     477       11760 :     if (isexactzero(a)) return scalarpol_shallow(a,v);
     478           0 :     n = d+2; P = cgetg(n+1, t_POL);
     479           0 :     P[1] = evalsigne(0) | evalvarn(v);
     480             :   }
     481             :   else
     482             :   {
     483     7625181 :     n = d+2; P = cgetg(n+1, t_POL);
     484     7625181 :     P[1] = evalsigne(1) | evalvarn(v);
     485             :   }
     486    15592359 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     487     7625181 :   gel(P,i) = a; return P;
     488             : }
     489             : GEN
     490     1860804 : monomialcopy(GEN a, long d, long v)
     491             : {
     492             :   long i, n;
     493             :   GEN P;
     494     1860804 :   if (d < 0) {
     495          14 :     if (isrationalzero(a)) return pol_0(v);
     496          14 :     retmkrfrac(gcopy(a), pol_xn(-d, v));
     497             :   }
     498     1860790 :   if (gequal0(a))
     499             :   {
     500           7 :     if (isexactzero(a)) return scalarpol(a,v);
     501           0 :     n = d+2; P = cgetg(n+1, t_POL);
     502           0 :     P[1] = evalsigne(0) | evalvarn(v);
     503             :   }
     504             :   else
     505             :   {
     506     1860783 :     n = d+2; P = cgetg(n+1, t_POL);
     507     1860783 :     P[1] = evalsigne(1) | evalvarn(v);
     508             :   }
     509     3505075 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     510     1860783 :   gel(P,i) = gcopy(a); return P;
     511             : }
     512             : GEN
     513       25937 : pol_x_powers(long N, long v)
     514             : {
     515       25937 :   GEN L = cgetg(N+1,t_VEC);
     516             :   long i;
     517       88868 :   for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
     518       25937 :   return L;
     519             : }
     520             : 
     521             : GEN
     522           0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
     523             : {
     524           0 :   return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
     525             : }
     526             : 
     527             : GEN
     528           0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
     529             : {
     530           0 :   return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
     531             : }
     532             : 
     533             : /*******************************************************************/
     534             : /*                                                                 */
     535             : /*                             Fq                                  */
     536             : /*                                                                 */
     537             : /*******************************************************************/
     538             : 
     539             : GEN
     540     4441698 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     541             : {
     542             :   (void)T;
     543     4441698 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     544             :   {
     545      541302 :     case 0: return Fp_add(x,y,p);
     546      205324 :     case 1: return FpX_Fp_add(x,y,p);
     547       57842 :     case 2: return FpX_Fp_add(y,x,p);
     548     3637230 :     case 3: return FpX_add(x,y,p);
     549             :   }
     550             :   return NULL;/*LCOV_EXCL_LINE*/
     551             : }
     552             : 
     553             : GEN
     554     4761831 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
     555             : {
     556             :   (void)T;
     557     4761831 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     558             :   {
     559      194946 :     case 0: return Fp_sub(x,y,p);
     560        8601 :     case 1: return FpX_Fp_sub(x,y,p);
     561       10149 :     case 2: return Fp_FpX_sub(x,y,p);
     562     4548135 :     case 3: return FpX_sub(x,y,p);
     563             :   }
     564             :   return NULL;/*LCOV_EXCL_LINE*/
     565             : }
     566             : 
     567             : GEN
     568      562981 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
     569             : {
     570             :   (void)T;
     571      562981 :   return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
     572             : }
     573             : 
     574             : GEN
     575       12829 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
     576             : {
     577             :   (void)T;
     578       12829 :   return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
     579             : }
     580             : 
     581             : /* If T==NULL do not reduce*/
     582             : GEN
     583     5726009 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
     584             : {
     585     5726009 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     586             :   {
     587      664363 :     case 0: return Fp_mul(x,y,p);
     588       67085 :     case 1: return FpX_Fp_mul(x,y,p);
     589      123597 :     case 2: return FpX_Fp_mul(y,x,p);
     590     4870964 :     case 3: if (T) return FpXQ_mul(x,y,T,p);
     591     2841184 :             else return FpX_mul(x,y,p);
     592             :   }
     593             :   return NULL;/*LCOV_EXCL_LINE*/
     594             : }
     595             : 
     596             : /* If T==NULL do not reduce*/
     597             : GEN
     598      372644 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
     599             : {
     600             :   (void) T;
     601      372644 :   return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
     602             : }
     603             : 
     604             : /* y t_INT */
     605             : GEN
     606       33137 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
     607             : {
     608             :   (void)T;
     609        5956 :   return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
     610       39093 :                           : Fp_mul(x,y,p);
     611             : }
     612             : /* If T==NULL do not reduce*/
     613             : GEN
     614      255569 : Fq_sqr(GEN x, GEN T, GEN p)
     615             : {
     616      255569 :   if (typ(x) == t_POL)
     617             :   {
     618       12047 :     if (T) return FpXQ_sqr(x,T,p);
     619           0 :     else return FpX_sqr(x,p);
     620             :   }
     621             :   else
     622      243522 :     return Fp_sqr(x,p);
     623             : }
     624             : 
     625             : GEN
     626           0 : Fq_neg_inv(GEN x, GEN T, GEN p)
     627             : {
     628           0 :   if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
     629           0 :   return FpXQ_inv(FpX_neg(x,p),T,p);
     630             : }
     631             : 
     632             : GEN
     633           0 : Fq_invsafe(GEN x, GEN pol, GEN p)
     634             : {
     635           0 :   if (typ(x) == t_INT) return Fp_invsafe(x,p);
     636           0 :   return FpXQ_invsafe(x,pol,p);
     637             : }
     638             : 
     639             : GEN
     640       36051 : Fq_inv(GEN x, GEN pol, GEN p)
     641             : {
     642       36051 :   if (typ(x) == t_INT) return Fp_inv(x,p);
     643       29803 :   return FpXQ_inv(x,pol,p);
     644             : }
     645             : 
     646             : GEN
     647      303009 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
     648             : {
     649      303009 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     650             :   {
     651      287245 :     case 0: return Fp_div(x,y,p);
     652       10241 :     case 1: return FpX_Fp_div(x,y,p);
     653         280 :     case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
     654        5243 :     case 3: return FpXQ_div(x,y,pol,p);
     655             :   }
     656             :   return NULL;/*LCOV_EXCL_LINE*/
     657             : }
     658             : 
     659             : GEN
     660      144516 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
     661             : {
     662      144516 :   if (typ(x) == t_INT) return Fp_pow(x,n,p);
     663       34847 :   return FpXQ_pow(x,n,pol,p);
     664             : }
     665             : 
     666             : GEN
     667       14693 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
     668             : {
     669       14693 :   if (typ(x) == t_INT) return Fp_powu(x,n,p);
     670         749 :   return FpXQ_powu(x,n,pol,p);
     671             : }
     672             : 
     673             : GEN
     674      709980 : Fq_sqrt(GEN x, GEN T, GEN p)
     675             : {
     676      709980 :   if (typ(x) == t_INT)
     677             :   {
     678      699160 :     if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
     679         336 :     x = scalarpol_shallow(x, get_FpX_var(T));
     680             :   }
     681       11156 :   return FpXQ_sqrt(x,T,p);
     682             : }
     683             : GEN
     684       60193 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
     685             : {
     686       60193 :   if (typ(x) == t_INT)
     687             :   {
     688             :     long d;
     689       59927 :     if (!T) return Fp_sqrtn(x,n,p,zeta);
     690         147 :     d = get_FpX_degree(T);
     691         147 :     if (ugcdiu(n,d) == 1)
     692             :     {
     693          14 :       if (!zeta) return Fp_sqrtn(x,n,p,NULL);
     694             :       /* gcd(n,p-1)=gcd(n,q-1): same number of solutions in Fp and F_q */
     695           7 :       if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
     696           7 :         return Fp_sqrtn(x,n,p,zeta);
     697             :     }
     698         133 :     x = scalarpol(x, get_FpX_var(T)); /* left on stack */
     699             :   }
     700         399 :   return FpXQ_sqrtn(x,n,T,p,zeta);
     701             : }
     702             : 
     703             : struct _Fq_field
     704             : {
     705             :   GEN T, p;
     706             : };
     707             : 
     708             : static GEN
     709      302463 : _Fq_red(void *E, GEN x)
     710      302463 : { struct _Fq_field *s = (struct _Fq_field *)E;
     711      302463 :   return Fq_red(x, s->T, s->p);
     712             : }
     713             : 
     714             : static GEN
     715      357525 : _Fq_add(void *E, GEN x, GEN y)
     716             : {
     717             :   (void) E;
     718      357525 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     719             :   {
     720          14 :     case 0: return addii(x,y);
     721           0 :     case 1: return ZX_Z_add(x,y);
     722       15918 :     case 2: return ZX_Z_add(y,x);
     723      341593 :     default: return ZX_add(x,y);
     724             :   }
     725             : }
     726             : 
     727             : static GEN
     728      315028 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
     729             : 
     730             : static GEN
     731      237559 : _Fq_mul(void *E, GEN x, GEN y)
     732             : {
     733             :   (void) E;
     734      237559 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     735             :   {
     736         133 :     case 0: return mulii(x,y);
     737        1085 :     case 1: return ZX_Z_mul(x,y);
     738        1043 :     case 2: return ZX_Z_mul(y,x);
     739      235298 :     default: return ZX_mul(x,y);
     740             :   }
     741             : }
     742             : 
     743             : static GEN
     744        9331 : _Fq_inv(void *E, GEN x)
     745        9331 : { struct _Fq_field *s = (struct _Fq_field *)E;
     746        9331 :   return Fq_inv(x,s->T,s->p);
     747             : }
     748             : 
     749             : static int
     750       38388 : _Fq_equal0(GEN x) { return signe(x)==0; }
     751             : 
     752             : static GEN
     753       13965 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
     754             : 
     755             : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
     756             :                                        _Fq_inv,_Fq_equal0,_Fq_s};
     757             : 
     758        4165 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
     759             : {
     760        4165 :   if (!T)
     761           0 :     return get_Fp_field(E, p);
     762             :   else
     763             :   {
     764        4165 :     GEN z = new_chunk(sizeof(struct _Fq_field));
     765        4165 :     struct _Fq_field *e = (struct _Fq_field *) z;
     766        4165 :     e->T = T; e->p  = p; *E = (void*)e;
     767        4165 :     return &Fq_field;
     768             :   }
     769             : }
     770             : 
     771             : /*******************************************************************/
     772             : /*                                                                 */
     773             : /*                             Fq[X]                               */
     774             : /*                                                                 */
     775             : /*******************************************************************/
     776             : /* P(X + c) */
     777             : GEN
     778         266 : FpX_translate(GEN P, GEN c, GEN p)
     779             : {
     780         266 :   pari_sp av = avma;
     781             :   GEN Q, *R;
     782             :   long i, k, n;
     783             : 
     784         266 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     785         266 :   Q = leafcopy(P);
     786         266 :   R = (GEN*)(Q+2); n = degpol(P);
     787        3738 :   for (i=1; i<=n; i++)
     788             :   {
     789      118153 :     for (k=n-i; k<n; k++)
     790      114681 :       R[k] = Fp_add(R[k], Fp_mul(c, R[k+1], p), p);
     791             : 
     792        3472 :     if (gc_needed(av,2))
     793             :     {
     794           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpX_translate, i = %ld/%ld", i,n);
     795           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     796             :     }
     797             :   }
     798         266 :   return gerepilecopy(av, FpX_renormalize(Q, lg(Q)));
     799             : }
     800             : /* P(X + c), c an Fq */
     801             : GEN
     802       33880 : FqX_translate(GEN P, GEN c, GEN T, GEN p)
     803             : {
     804       33880 :   pari_sp av = avma;
     805             :   GEN Q, *R;
     806             :   long i, k, n;
     807             : 
     808             :   /* signe works for t_(INT|POL) */
     809       33880 :   if (!signe(P) || !signe(c)) return RgX_copy(P);
     810       33880 :   Q = leafcopy(P);
     811       33880 :   R = (GEN*)(Q+2); n = degpol(P);
     812      150059 :   for (i=1; i<=n; i++)
     813             :   {
     814      433559 :     for (k=n-i; k<n; k++)
     815      317380 :       R[k] = Fq_add(R[k], Fq_mul(c, R[k+1], T, p), T, p);
     816             : 
     817      116179 :     if (gc_needed(av,2))
     818             :     {
     819           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FqX_translate, i = %ld/%ld", i,n);
     820           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     821             :     }
     822             :   }
     823       33880 :   return gerepilecopy(av, FpXQX_renormalize(Q, lg(Q)));
     824             : }
     825             : 
     826             : GEN
     827        5579 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
     828             : {
     829        5579 :   pari_sp ltop = avma;
     830             :   long k;
     831             :   GEN W;
     832        5579 :   if (lgefint(p) == 3)
     833             :   {
     834         900 :     ulong pp = p[2];
     835         900 :     GEN Tl = ZX_to_Flx(T, pp);
     836         900 :     GEN Vl = FqV_to_FlxV(V, T, p);
     837         900 :     Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
     838         900 :     return gerepileupto(ltop, FlxX_to_ZXX(Tl));
     839             :   }
     840        4679 :   W = cgetg(lg(V),t_VEC);
     841       60157 :   for(k=1; k < lg(V); k++)
     842       55478 :     gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
     843        4679 :   return gerepileupto(ltop, FpXQXV_prod(W, T, p));
     844             : }
     845             : 
     846             : GEN
     847      124393 : FqV_red(GEN x, GEN T, GEN p)
     848     1102300 : { pari_APPLY_same(Fq_red(gel(x,i), T, p)) }
     849             : 
     850             : GEN
     851           0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
     852             : {
     853           0 :   if (!T) return FpC_add(x, y, p);
     854           0 :   pari_APPLY_type(t_COL, Fq_add(gel(x,i), gel(y,i), T, p))
     855             : }
     856             : 
     857             : GEN
     858           0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
     859             : {
     860           0 :   if (!T) return FpC_sub(x, y, p);
     861           0 :   pari_APPLY_type(t_COL, Fq_sub(gel(x,i), gel(y,i), T, p))
     862             : }
     863             : 
     864             : GEN
     865           0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
     866             : {
     867           0 :   if (!T) return FpC_Fp_mul(x, y, p);
     868           0 :   pari_APPLY_type(t_COL, Fq_mul(gel(x,i),y,T,p))
     869             : }
     870             : 
     871             : GEN
     872          77 : FqC_FqV_mul(GEN x, GEN y, GEN T, GEN p)
     873             : {
     874          77 :   long i,j, lx=lg(x), ly=lg(y);
     875             :   GEN z;
     876          77 :   if (ly==1) return cgetg(1,t_MAT);
     877          77 :   z = cgetg(ly,t_MAT);
     878         623 :   for (j=1; j < ly; j++)
     879             :   {
     880         546 :     GEN zj = cgetg(lx,t_COL);
     881        3360 :     for (i=1; i<lx; i++) gel(zj,i) = Fq_mul(gel(x,i),gel(y,j), T, p);
     882         546 :     gel(z, j) = zj;
     883             :   }
     884          77 :   return z;
     885             : }
     886             : 
     887             : GEN
     888         900 : FqV_to_FlxV(GEN x, GEN T, GEN pp)
     889             : {
     890         900 :   long vT = evalvarn(get_FpX_var(T));
     891         900 :   ulong p = pp[2];
     892        3914 :   pari_APPLY_type(t_VEC, typ(gel(x,i))==t_INT?  Z_to_Flx(gel(x,i), p, vT)
     893             :                                              : ZX_to_Flx(gel(x,i), p))
     894             : }
     895             : 
     896             : GEN
     897       58489 : FqC_to_FlxC(GEN x, GEN T, GEN pp)
     898             : {
     899       58489 :   long vT = evalvarn(get_FpX_var(T));
     900       58489 :   ulong p = pp[2];
     901     1645462 :   pari_APPLY_type(t_COL, typ(gel(x,i))==t_INT?  Z_to_Flx(gel(x,i), p, vT)
     902             :                                              : ZX_to_Flx(gel(x,i), p))
     903             : }
     904             : 
     905             : GEN
     906       10121 : FqM_to_FlxM(GEN x, GEN T, GEN p)
     907       68610 : { pari_APPLY_same(FqC_to_FlxC(gel(x,i), T, p)) }
     908             : 
     909             : GEN
     910        7009 : FpXC_center(GEN x, GEN p, GEN pov2)
     911       75000 : { pari_APPLY_type(t_COL, FpX_center(gel(x,i), p, pov2)) }
     912             : 
     913             : GEN
     914        1705 : FpXM_center(GEN x, GEN p, GEN pov2)
     915        8714 : { pari_APPLY_same(FpXC_center(gel(x,i), p, pov2)) }
     916             : 
     917             : /*******************************************************************/
     918             : /*                                                                 */
     919             : /*                          GENERIC CRT                            */
     920             : /*                                                                 */
     921             : /*******************************************************************/
     922             : static GEN
     923     3793005 : primelist(forprime_t *S, long n, GEN dB)
     924             : {
     925     3793005 :   GEN P = cgetg(n+1, t_VECSMALL);
     926     3793004 :   long i = 1;
     927             :   ulong p;
     928     9116826 :   while (i <= n && (p = u_forprime_next(S)))
     929     5323822 :     if (!dB || umodiu(dB, p)) P[i++] = p;
     930     3793006 :   return P;
     931             : }
     932             : 
     933             : void
     934     3598569 : gen_inccrt_i(const char *str, GEN worker, GEN dB, long n, long mmin,
     935             :              forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
     936             :              GEN center(GEN, GEN, GEN))
     937             : {
     938     3598569 :   long m = mmin? minss(mmin, n): usqrt(n);
     939             :   GEN  H, P, mod;
     940             :   pari_timer ti;
     941     3598570 :   if (DEBUGLEVEL > 4)
     942             :   {
     943           0 :     timer_start(&ti);
     944           0 :     err_printf("%s: nb primes: %ld\n",str, n);
     945             :   }
     946     3598570 :   if (m == 1)
     947             :   {
     948     3531311 :     GEN P = primelist(S, n, dB);
     949     3531311 :     GEN done = closure_callgen1(worker, P);
     950     3531311 :     H = gel(done,1);
     951     3531311 :     mod = gel(done,2);
     952     3531311 :     if (!*pH && center) H = center(H, mod, shifti(mod,-1));
     953     3531311 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
     954             :   }
     955             :   else
     956             :   {
     957       67259 :     long i, s = (n+m-1)/m, r = m - (m*s-n), di = 0;
     958             :     struct pari_mt pt;
     959       67259 :     long pending = 0;
     960       67259 :     H = cgetg(m+1, t_VEC); P = cgetg(m+1, t_VEC);
     961       67257 :     mt_queue_start_lim(&pt, worker, m);
     962      351791 :     for (i=1; i<=m || pending; i++)
     963             :     {
     964             :       GEN done;
     965      284531 :       GEN pr = i <= m ? mkvec(primelist(S, i<=r ? s: s-1, dB)): NULL;
     966      284532 :       mt_queue_submit(&pt, i, pr);
     967      284533 :       done = mt_queue_get(&pt, NULL, &pending);
     968      284533 :       if (done)
     969             :       {
     970      261697 :         di++;
     971      261697 :         gel(H, di) = gel(done,1);
     972      261697 :         gel(P, di) = gel(done,2);
     973      261697 :         if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
     974             :       }
     975             :     }
     976       67260 :     mt_queue_end(&pt);
     977       67260 :     if (DEBUGLEVEL>5) err_printf("\n");
     978       67260 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
     979       67260 :     H = crt(H, P, &mod);
     980       67260 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: chinese", str);
     981             :   }
     982     3598571 :   if (*pH) H = crt(mkvec2(*pH, H), mkvec2(*pmod, mod), &mod);
     983     3598571 :   *pH = H; *pmod = mod;
     984     3598571 : }
     985             : void
     986      272670 : gen_inccrt(const char *str, GEN worker, GEN dB, long n, long mmin,
     987             :            forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
     988             :            GEN center(GEN, GEN, GEN))
     989             : {
     990      272670 :   pari_sp av = avma;
     991      272670 :   gen_inccrt_i(str, worker, dB, n, mmin, S, pH, pmod, crt, center);
     992      272671 :   gerepileall(av, 2, pH, pmod);
     993      272671 : }
     994             : 
     995             : GEN
     996      135186 : gen_crt(const char *str, GEN worker, forprime_t *S, GEN dB, ulong bound, long mmin, GEN *pmod,
     997             :         GEN crt(GEN, GEN, GEN*), GEN center(GEN, GEN, GEN))
     998             : {
     999      135186 :   GEN mod = gen_1, H = NULL;
    1000             :   ulong e;
    1001             : 
    1002      135186 :   bound++;
    1003      270373 :   while (bound > (e = expi(mod)))
    1004             :   {
    1005      135186 :     long n = (bound - e) / expu(S->p) + 1;
    1006      135186 :     gen_inccrt(str, worker, dB, n, mmin, S, &H, &mod, crt, center);
    1007             :   }
    1008      135187 :   if (pmod) *pmod = mod;
    1009      135187 :   return H;
    1010             : }
    1011             : 
    1012             : /*******************************************************************/
    1013             : /*                                                                 */
    1014             : /*                          MODULAR GCD                            */
    1015             : /*                                                                 */
    1016             : /*******************************************************************/
    1017             : /* return z = a mod q, b mod p (p,q) = 1; qinv = 1/q mod p; a in ]-q,q] */
    1018             : static GEN
    1019      452445 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq, GEN pq2)
    1020             : {
    1021      452445 :   ulong d, amod = umodiu(a, p);
    1022      452445 :   pari_sp av = avma;
    1023             :   GEN ax;
    1024             : 
    1025      452445 :   if (b == amod) return NULL;
    1026      217145 :   d = Fl_mul(Fl_sub(b, amod, p), qinv, p); /* != 0 */
    1027      217145 :   if (d >= 1 + (p>>1))
    1028      107473 :     ax = subii(a, mului(p-d, q));
    1029             :   else
    1030             :   {
    1031      109672 :     ax = addii(a, mului(d, q)); /* in ]0, pq[ assuming a in ]-q,q[ */
    1032      109672 :     if (cmpii(ax,pq2) > 0) ax = subii(ax,pq);
    1033             :   }
    1034      217145 :   return gerepileuptoint(av, ax);
    1035             : }
    1036             : GEN
    1037         364 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
    1038             : GEN
    1039        2289 : ZX_init_CRT(GEN Hp, ulong p, long v)
    1040             : {
    1041        2289 :   long i, l = lg(Hp), lim = (long)(p>>1);
    1042        2289 :   GEN H = cgetg(l, t_POL);
    1043        2289 :   H[1] = evalsigne(1) | evalvarn(v);
    1044       29311 :   for (i=2; i<l; i++)
    1045       27022 :     gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
    1046        2289 :   return ZX_renormalize(H,l);
    1047             : }
    1048             : 
    1049             : GEN
    1050        2597 : ZM_init_CRT(GEN Hp, ulong p)
    1051             : {
    1052        2597 :   long i,j, m, l = lg(Hp), lim = (long)(p>>1);
    1053        2597 :   GEN c, cp, H = cgetg(l, t_MAT);
    1054        2597 :   if (l==1) return H;
    1055        2597 :   m = lgcols(Hp);
    1056        8988 :   for (j=1; j<l; j++)
    1057             :   {
    1058        6391 :     cp = gel(Hp,j);
    1059        6391 :     c = cgetg(m, t_COL);
    1060        6391 :     gel(H,j) = c;
    1061       76860 :     for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
    1062             :   }
    1063        2597 :   return H;
    1064             : }
    1065             : 
    1066             : int
    1067        7511 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
    1068             : {
    1069        7511 :   GEN h, q = *ptq, qp = muliu(q,p);
    1070        7511 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1071        7511 :   int stable = 1;
    1072        7511 :   h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp,shifti(qp,-1));
    1073        7511 :   if (h) { *H = h; stable = 0; }
    1074        7511 :   *ptq = qp; return stable;
    1075             : }
    1076             : 
    1077             : static int
    1078        8922 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
    1079             : {
    1080        8922 :   GEN H = *ptH, h, qp2 = shifti(qp,-1);
    1081        8922 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1082        8922 :   long i, l = lg(H), lp = lg(Hp);
    1083        8922 :   int stable = 1;
    1084             : 
    1085        8922 :   if (l < lp)
    1086             :   { /* degree increases */
    1087           0 :     GEN x = cgetg(lp, t_POL);
    1088           0 :     for (i=1; i<l; i++)  x[i] = H[i];
    1089           0 :     for (   ; i<lp; i++) gel(x,i) = gen_0;
    1090           0 :     *ptH = H = x;
    1091           0 :     stable = 0;
    1092        8922 :   } else if (l > lp)
    1093             :   { /* degree decreases */
    1094           0 :     GEN x = cgetg(l, t_VECSMALL);
    1095           0 :     for (i=1; i<lp; i++)  x[i] = Hp[i];
    1096           0 :     for (   ; i<l; i++) x[i] = 0;
    1097           0 :     Hp = x; lp = l;
    1098             :   }
    1099      133197 :   for (i=2; i<lp; i++)
    1100             :   {
    1101      124275 :     h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp,qp2);
    1102      124275 :     if (h) { gel(H,i) = h; stable = 0; }
    1103             :   }
    1104        8922 :   (void)ZX_renormalize(H,lp);
    1105        8922 :   return stable;
    1106             : }
    1107             : 
    1108             : int
    1109           0 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
    1110             : {
    1111           0 :   GEN q = *ptq, qp = muliu(q,p);
    1112           0 :   int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
    1113           0 :   *ptq = qp; return stable;
    1114             : }
    1115             : 
    1116             : int
    1117        4484 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1118             : {
    1119        4484 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1120        4484 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1121        4484 :   long i,j, l = lg(H), m = lgcols(H);
    1122        4484 :   int stable = 1;
    1123       16385 :   for (j=1; j<l; j++)
    1124      136502 :     for (i=1; i<m; i++)
    1125             :     {
    1126      124601 :       h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp,qp2);
    1127      124601 :       if (h) { gcoeff(H,i,j) = h; stable = 0; }
    1128             :     }
    1129        4484 :   *ptq = qp; return stable;
    1130             : }
    1131             : 
    1132             : GEN
    1133         952 : ZXM_init_CRT(GEN Hp, long deg, ulong p)
    1134             : {
    1135             :   long i, j, k;
    1136             :   GEN H;
    1137         952 :   long m, l = lg(Hp), lim = (long)(p>>1), n;
    1138         952 :   H = cgetg(l, t_MAT);
    1139         952 :   if (l==1) return H;
    1140         952 :   m = lgcols(Hp);
    1141         952 :   n = deg + 3;
    1142        3850 :   for (j=1; j<l; j++)
    1143             :   {
    1144        2898 :     GEN cp = gel(Hp,j);
    1145        2898 :     GEN c = cgetg(m, t_COL);
    1146        2898 :     gel(H,j) = c;
    1147       42413 :     for (i=1; i<m; i++)
    1148             :     {
    1149       39515 :       GEN dp = gel(cp, i);
    1150       39515 :       long l = lg(dp);
    1151       39515 :       GEN d = cgetg(n, t_POL);
    1152       39515 :       gel(c, i) = d;
    1153       39515 :       d[1] = dp[1] | evalsigne(1);
    1154       79905 :       for (k=2; k<l; k++)
    1155       40390 :         gel(d,k) = stoi(Fl_center(dp[k], p, lim));
    1156       78547 :       for (   ; k<n; k++)
    1157       39032 :         gel(d,k) = gen_0;
    1158             :     }
    1159             :   }
    1160         952 :   return H;
    1161             : }
    1162             : 
    1163             : int
    1164         930 : ZXM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1165             : {
    1166         930 :   GEN v, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1167         930 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1168         930 :   long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
    1169         930 :   int stable = 1;
    1170        5727 :   for (j=1; j<l; j++)
    1171      102816 :     for (i=1; i<m; i++)
    1172             :     {
    1173       98019 :       GEN h = gmael(H,j,i), hp = gmael(Hp,j,i);
    1174       98019 :       long lh = lg(hp);
    1175      207083 :       for (k=2; k<lh; k++)
    1176             :       {
    1177      109064 :         v = Fl_chinese_coprime(gel(h,k),uel(hp,k),q,p,qinv,qp,qp2);
    1178      109064 :         if (v) { gel(h,k) = v; stable = 0; }
    1179             :       }
    1180      185013 :       for (; k<n; k++)
    1181             :       {
    1182       86994 :         v = Fl_chinese_coprime(gel(h,k),0,q,p,qinv,qp,qp2);
    1183       86994 :         if (v) { gel(h,k) = v; stable = 0; }
    1184             :       }
    1185             :     }
    1186         930 :   *ptq = qp; return stable;
    1187             : }
    1188             : 
    1189             : /* record the degrees of Euclidean remainders (make them as large as
    1190             :  * possible : smaller values correspond to a degenerate sequence) */
    1191             : static void
    1192        3255 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
    1193             : {
    1194             :   long da,db,dc, ind;
    1195        3255 :   pari_sp av = avma;
    1196             : 
    1197        3255 :   if (lgpol(a)==0 || lgpol(b)==0) return;
    1198        3213 :   da = degpol(a);
    1199        3213 :   db = degpol(b);
    1200        3213 :   if (db > da)
    1201           0 :   { swapspec(a,b, da,db); }
    1202        3213 :   else if (!da) return;
    1203        3213 :   ind = 0;
    1204       17108 :   while (db)
    1205             :   {
    1206       13895 :     GEN c = Flx_rem(a,b, p);
    1207       13895 :     a = b; b = c; dc = degpol(c);
    1208       13895 :     if (dc < 0) break;
    1209             : 
    1210       13895 :     ind++;
    1211       13895 :     if (dc > dglist[ind]) dglist[ind] = dc;
    1212       13895 :     if (gc_needed(av,2))
    1213             :     {
    1214           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1215           0 :       gerepileall(av, 2, &a,&b);
    1216             :     }
    1217       13895 :     db = dc; /* = degpol(b) */
    1218             :   }
    1219        3213 :   if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
    1220        3213 :   set_avma(av);
    1221             : }
    1222             : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
    1223             :  * "generic" degree sequence as given by dglist, set *Ci and return
    1224             :  * resultant(a,b). Modular version of Collins's subresultant */
    1225             : static ulong
    1226       61965 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
    1227             : {
    1228             :   long da,db,dc, ind;
    1229       61965 :   ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
    1230       61965 :   int s = 1;
    1231       61965 :   pari_sp av = avma;
    1232             : 
    1233       61965 :   *C0 = 1; *C1 = 0;
    1234       61965 :   if (lgpol(a)==0 || lgpol(b)==0) return 0;
    1235       61881 :   da = degpol(a);
    1236       61881 :   db = degpol(b);
    1237       61881 :   if (db > da)
    1238             :   {
    1239           0 :     swapspec(a,b, da,db);
    1240           0 :     if (both_odd(da,db)) s = -s;
    1241             :   }
    1242       61881 :   else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
    1243       61881 :   ind = 0;
    1244      584019 :   while (db)
    1245             :   { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
    1246             :      * da = deg a, db = deg b */
    1247      522764 :     GEN c = Flx_rem(a,b, p);
    1248      522764 :     long delta = da - db;
    1249             : 
    1250      522764 :     if (both_odd(da,db)) s = -s;
    1251      522764 :     lb = Fl_mul(b[db+2], cb, p);
    1252      522764 :     a = b; b = c; dc = degpol(c);
    1253      522764 :     ind++;
    1254      522764 :     if (dc != dglist[ind]) return gc_ulong(av,0); /* degenerates */
    1255      522138 :     if (g == h)
    1256             :     { /* frequent */
    1257      516479 :       ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
    1258      516479 :       ca = cb;
    1259      516479 :       cb = cc;
    1260             :     }
    1261             :     else
    1262             :     {
    1263        5659 :       ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
    1264        5659 :       ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
    1265        5659 :       ca = cb;
    1266        5659 :       cb = Fl_div(cc, ghdelta, p);
    1267             :     }
    1268      522138 :     da = db; /* = degpol(a) */
    1269      522138 :     db = dc; /* = degpol(b) */
    1270             : 
    1271      522138 :     g = lb;
    1272      522138 :     if (delta == 1)
    1273      504849 :       h = g; /* frequent */
    1274             :     else
    1275       17289 :       h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
    1276             : 
    1277      522138 :     if (gc_needed(av,2))
    1278             :     {
    1279           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1280           0 :       gerepileall(av, 2, &a,&b);
    1281             :     }
    1282             :   }
    1283       61255 :   if (da > 1) return 0; /* Failure */
    1284             :   /* last non-constant polynomial has degree 1 */
    1285       61255 :   *C0 = Fl_mul(ca, a[2], p);
    1286       61255 :   *C1 = Fl_mul(ca, a[3], p);
    1287       61255 :   res = Fl_mul(cb, b[2], p);
    1288       61255 :   if (s == -1) res = p - res;
    1289       61255 :   return gc_ulong(av,res);
    1290             : }
    1291             : 
    1292             : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
    1293             :  * Return 0 in case of degree drop. */
    1294             : static GEN
    1295       65220 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
    1296             : {
    1297             :   GEN z;
    1298       65220 :   long i, lb = lg(Q);
    1299       65220 :   ulong leadz = Flx_eval(leading_coeff(Q), x, p);
    1300       65220 :   long vs=mael(Q,2,1);
    1301       65220 :   if (!leadz) return zero_Flx(vs);
    1302             : 
    1303       65094 :   z = cgetg(lb, t_VECSMALL); z[1] = vs;
    1304      629369 :   for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
    1305       65094 :   z[i] = leadz; return z;
    1306             : }
    1307             : 
    1308             : GEN
    1309        1232 : FpXY_FpXQ_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
    1310             : {
    1311        1232 :   pari_sp av = avma;
    1312        1232 :   long i, lb = lg(Q);
    1313             :   GEN z;
    1314        1232 :   if (lb == 2) return pol_0(vx);
    1315        1232 :   z = gel(Q, lb-1);
    1316        1232 :   if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
    1317             : 
    1318        1232 :   if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
    1319       26572 :   for (i=lb-2; i>=2; i--)
    1320             :   {
    1321       25340 :     GEN c = gel(Q,i);
    1322       25340 :     z = FqX_Fq_mul(z, y, T, p);
    1323       25340 :     z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
    1324             :   }
    1325        1232 :   return gerepileupto(av, z);
    1326             : }
    1327             : 
    1328             : static GEN
    1329       21763 : ZX_norml1(GEN x)
    1330             : {
    1331       21763 :   long i, l = lg(x);
    1332             :   GEN s;
    1333             : 
    1334       21763 :   if (l == 2) return gen_0;
    1335       12908 :   s = gel(x, l-1); /* != 0 */
    1336       47784 :   for (i = l-2; i > 1; i--) {
    1337       34876 :     GEN xi = gel(x,i);
    1338       34876 :     if (!signe(x)) continue;
    1339       34876 :     s = addii_sign(s,1, xi,1);
    1340             :   }
    1341       12908 :   return s;
    1342             : }
    1343             : 
    1344             : static GEN
    1345        2635 : L2_bound(GEN nf, GEN den, GEN *pt_roots)
    1346             : {
    1347        2635 :   GEN M, L, prep, T = nf_get_pol(nf), tozk = nf_get_invzk(nf);
    1348        2635 :   long prec = ZM_max_lg(tozk) + ZX_max_lg(T) + nbits2prec(degpol(T));
    1349        2635 :   (void)initgaloisborne(nf, den? den: gen_1, prec, &L, &prep, NULL);
    1350        2635 :   M = vandermondeinverse(L, RgX_gtofp(T,prec), den, prep);
    1351        2635 :   *pt_roots = L;
    1352        2635 :   return RgM_fpnorml2(RgM_mul(tozk,M), DEFAULTPREC);
    1353             : }
    1354             : 
    1355             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1356             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1357             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1358             :  * Return e such that Res(A, B) < 2^e */
    1359             : static GEN
    1360        5043 : RgX_RgXY_ResBound(GEN A, GEN B, long prec)
    1361             : {
    1362        5043 :   pari_sp av = avma, av2;
    1363        5043 :   GEN a = gen_0, b = gen_0, bnd;
    1364        5043 :   long i , lA = lg(A), lB = lg(B);
    1365       25303 :   for (i=2; i<lA; i++)
    1366             :   {
    1367       20260 :     a = gadd(a, gabs(gnorm(gel(A,i)), prec));
    1368       20260 :     if (gc_needed(av,1))
    1369             :     {
    1370           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
    1371           0 :       a = gerepileupto(av, a);
    1372             :     }
    1373             :   }
    1374        5043 :   av2 = avma;
    1375       20911 :   for (i=2; i<lB; i++)
    1376             :   {
    1377       15868 :     GEN t = gel(B,i);
    1378       15868 :     if (typ(t) == t_POL) t = gnorml1(t, prec);
    1379       15868 :     b = gadd(b, gabs(gsqr(t), prec));
    1380       15868 :     if (gc_needed(av,1))
    1381             :     {
    1382           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
    1383           0 :       b = gerepileupto(av2, b);
    1384             :     }
    1385             :   }
    1386        5043 :   bnd = gsqrt(gmul(gpowgs(a, degpol(B)), gpowgs(b, degpol(A))), prec);
    1387        5043 :   return gerepileupto(av, bnd);
    1388             : }
    1389             : 
    1390             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1391             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1392             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1393             :  * Return e such that Res(A, B) < 2^e */
    1394             : ulong
    1395       48728 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
    1396             : {
    1397       48728 :   pari_sp av = avma;
    1398       48728 :   GEN a = gen_0, b = gen_0;
    1399       48728 :   long i , lA = lg(A), lB = lg(B);
    1400             :   double loga, logb;
    1401      405359 :   for (i=2; i<lA; i++)
    1402             :   {
    1403      356631 :     a = addii(a, sqri(gel(A,i)));
    1404      356631 :     if (gc_needed(av,1))
    1405             :     {
    1406           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1407           0 :       a = gerepileupto(av, a);
    1408             :     }
    1409             :   }
    1410       48728 :   loga = dbllog2(a); set_avma(av);
    1411      363964 :   for (i=2; i<lB; i++)
    1412             :   {
    1413      315236 :     GEN t = gel(B,i);
    1414      315236 :     if (typ(t) == t_POL) t = ZX_norml1(t);
    1415      315236 :     b = addii(b, sqri(t));
    1416      315236 :     if (gc_needed(av,1))
    1417             :     {
    1418           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1419           0 :       b = gerepileupto(av, b);
    1420             :     }
    1421             :   }
    1422       48728 :   logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
    1423       48728 :   i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
    1424       48728 :   return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
    1425             : }
    1426             : /* special case B = A' */
    1427             : static ulong
    1428       60611 : ZX_discbound(GEN A)
    1429             : {
    1430       60611 :   pari_sp av = avma;
    1431       60611 :   GEN a = gen_0, b = gen_0;
    1432       60611 :   long i , lA = lg(A), dA = degpol(A);
    1433             :   double loga, logb;
    1434      804271 :   for (i = 2; i < lA; i++)
    1435             :   {
    1436      743666 :     GEN c = sqri(gel(A,i));
    1437      743671 :     a = addii(a, c);
    1438      743677 :     if (i > 2) b = addii(b, mulii(c, sqru(i-2)));
    1439      743659 :     if (gc_needed(av,1))
    1440             :     {
    1441           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_discbound i = %ld",i);
    1442           0 :       gerepileall(av, 2, &a, &b);
    1443             :     }
    1444             :   }
    1445       60605 :   loga = dbllog2(a);
    1446       60611 :   logb = dbllog2(b); set_avma(av);
    1447       60611 :   i = (long)(((dA-1) * loga + dA * logb) / 2);
    1448       60611 :   return (i <= 0)? 1: 1 + (ulong)i;
    1449             : }
    1450             : 
    1451             : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
    1452             : static ulong
    1453      304703 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong la)
    1454             : {
    1455      304703 :   GEN ev = FlxY_evalx(b, n, p);
    1456      304720 :   long drop = lg(b) - lg(ev);
    1457      304720 :   ulong r = Flx_resultant(a, ev, p);
    1458      304699 :   if (drop && la != 1) r = Fl_mul(r, Fl_powu(la, drop,p),p);
    1459      304699 :   return r;
    1460             : }
    1461             : static GEN
    1462           4 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
    1463             : {
    1464           4 :   GEN ev = FpXY_evaly(b, n, p, vX);
    1465           4 :   long drop = db-degpol(ev);
    1466           4 :   GEN r = FpX_resultant(a, ev, p);
    1467           4 :   if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
    1468           4 :   return r;
    1469             : }
    1470             : 
    1471             : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1472             : /* Return a Fly */
    1473             : static GEN
    1474       17778 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, long dres, long sx)
    1475             : {
    1476             :   long i;
    1477       17778 :   ulong n, la = Flx_lead(a);
    1478       17778 :   GEN  x = cgetg(dres+2, t_VECSMALL);
    1479       17778 :   GEN  y = cgetg(dres+2, t_VECSMALL);
    1480             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1481             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1482      162737 :   for (i=0,n = 1; i < dres; n++)
    1483             :   {
    1484      144959 :     x[++i] = n;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1485      144962 :     x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1486             :   }
    1487       17778 :   if (i == dres)
    1488             :   {
    1489       14864 :     x[++i] = 0;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1490             :   }
    1491       17778 :   return Flv_polint(x,y, p, sx);
    1492             : }
    1493             : 
    1494             : static GEN
    1495        8446 : FlxX_pseudorem(GEN x, GEN y, ulong p)
    1496             : {
    1497        8446 :   long vx = varn(x), dx, dy, dz, i, lx, dp;
    1498        8446 :   pari_sp av = avma, av2;
    1499             : 
    1500        8446 :   if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
    1501        8446 :   (void)new_chunk(2);
    1502        8445 :   dx=degpol(x); x = RgX_recip_shallow(x)+2;
    1503        8449 :   dy=degpol(y); y = RgX_recip_shallow(y)+2; dz=dx-dy; dp = dz+1;
    1504        8448 :   av2 = avma;
    1505             :   for (;;)
    1506             :   {
    1507       69592 :     gel(x,0) = Flx_neg(gel(x,0), p); dp--;
    1508      260120 :     for (i=1; i<=dy; i++)
    1509      188829 :       gel(x,i) = Flx_add( Flx_mul(gel(y,0), gel(x,i), p),
    1510      190432 :                               Flx_mul(gel(x,0), gel(y,i), p), p );
    1511     1218904 :     for (   ; i<=dx; i++)
    1512     1143786 :       gel(x,i) = Flx_mul(gel(y,0), gel(x,i), p);
    1513       73913 :     do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
    1514       75118 :     if (dx < dy) break;
    1515       66683 :     if (gc_needed(av2,1))
    1516             :     {
    1517           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
    1518           0 :       gerepilecoeffs(av2,x,dx+1);
    1519             :     }
    1520             :   }
    1521        8435 :   if (dx < 0) return zero_Flx(0);
    1522        8435 :   lx = dx+3; x -= 2;
    1523        8435 :   x[0]=evaltyp(t_POL) | evallg(lx);
    1524        8434 :   x[1]=evalsigne(1) | evalvarn(vx);
    1525        8434 :   x = RgX_recip_shallow(x);
    1526        8451 :   if (dp)
    1527             :   { /* multiply by y[0]^dp   [beware dummy vars from FpX_FpXY_resultant] */
    1528        2235 :     GEN t = Flx_powu(gel(y,0), dp, p);
    1529        8932 :     for (i=2; i<lx; i++)
    1530        6698 :       gel(x,i) = Flx_mul(gel(x,i), t, p);
    1531             :   }
    1532        8450 :   return gerepilecopy(av, x);
    1533             : }
    1534             : 
    1535             : /* return a Flx */
    1536             : GEN
    1537        2827 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
    1538             : {
    1539        2827 :   pari_sp av = avma, av2;
    1540             :   long degq,dx,dy,du,dv,dr,signh;
    1541             :   GEN z,g,h,r,p1;
    1542             : 
    1543        2827 :   dx=degpol(u); dy=degpol(v); signh=1;
    1544        2827 :   if (dx < dy)
    1545             :   {
    1546           7 :     swap(u,v); lswap(dx,dy);
    1547           7 :     if (both_odd(dx, dy)) signh = -signh;
    1548             :   }
    1549        2827 :   if (dy < 0) return zero_Flx(sx);
    1550        2827 :   if (dy==0) return gerepileupto(av, Flx_powu(gel(v,2),dx,p));
    1551             : 
    1552        2827 :   g = h = pol1_Flx(sx); av2 = avma;
    1553             :   for(;;)
    1554             :   {
    1555        8446 :     r = FlxX_pseudorem(u,v,p); dr = lg(r);
    1556        8452 :     if (dr == 2) { set_avma(av); return zero_Flx(sx); }
    1557        8452 :     du = degpol(u); dv = degpol(v); degq = du-dv;
    1558        8451 :     u = v; p1 = g; g = leading_coeff(u);
    1559        8451 :     switch(degq)
    1560             :     {
    1561           0 :       case 0: break;
    1562        6202 :       case 1:
    1563        6202 :         p1 = Flx_mul(h,p1, p); h = g; break;
    1564        2249 :       default:
    1565        2249 :         p1 = Flx_mul(Flx_powu(h,degq,p), p1, p);
    1566        2247 :         h = Flx_div(Flx_powu(g,degq,p), Flx_powu(h,degq-1,p), p);
    1567             :     }
    1568        8443 :     if (both_odd(du,dv)) signh = -signh;
    1569        8439 :     v = FlxY_Flx_div(r, p1, p);
    1570        8440 :     if (dr==3) break;
    1571        5615 :     if (gc_needed(av2,1))
    1572             :     {
    1573           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_resultant, dr = %ld",dr);
    1574           0 :       gerepileall(av2,4, &u, &v, &g, &h);
    1575             :     }
    1576             :   }
    1577        2825 :   z = gel(v,2);
    1578        2825 :   if (dv > 1) z = Flx_div(Flx_powu(z,dv,p), Flx_powu(h,dv-1,p), p);
    1579        2825 :   if (signh < 0) z = Flx_neg(z,p);
    1580        2825 :   return gerepileupto(av, z);
    1581             : }
    1582             : 
    1583             : /* Warning:
    1584             :  * This function switches between valid and invalid variable ordering*/
    1585             : 
    1586             : static GEN
    1587        6079 : FlxY_to_FlyX(GEN b, long sv)
    1588             : {
    1589        6079 :   long i, n=-1;
    1590        6079 :   long sw = b[1]&VARNBITS;
    1591       21067 :   for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
    1592        6078 :   return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
    1593             : }
    1594             : 
    1595             : /* Return a Fly*/
    1596             : GEN
    1597        6079 : Flx_FlxY_resultant(GEN a, GEN b, ulong pp)
    1598             : {
    1599        6079 :   pari_sp ltop=avma;
    1600        6079 :   long dres = degpol(a)*degpol(b);
    1601        6079 :   long sx=a[1], sy=b[1]&VARNBITS;
    1602             :   GEN z;
    1603        6079 :   b = FlxY_to_FlyX(b,sx);
    1604        6079 :   if ((ulong)dres >= pp)
    1605        2825 :     z = FlxX_resultant(Fly_to_FlxY(a, sy), b, pp, sx);
    1606             :   else
    1607        3254 :     z = Flx_FlxY_resultant_polint(a, b, pp, (ulong)dres, sy);
    1608        6080 :   return gerepileupto(ltop,z);
    1609             : }
    1610             : 
    1611             : /* return a t_POL (in variable v >= 0) whose coeffs are the coeffs of b,
    1612             :  * in variable v. This is an incorrect PARI object if initially varn(b) << v.
    1613             :  * We could return a vector of coeffs, but it is convenient to have degpol()
    1614             :  * and friends available. Even in that case, it will behave nicely with all
    1615             :  * functions treating a polynomial as a vector of coeffs (eg poleval).
    1616             :  * FOR INTERNAL USE! */
    1617             : GEN
    1618       11887 : swap_vars(GEN b0, long v)
    1619             : {
    1620       11887 :   long i, n = RgX_degree(b0, v);
    1621             :   GEN b, x;
    1622       11887 :   if (n < 0) return pol_0(v);
    1623       11887 :   b = cgetg(n+3, t_POL); x = b + 2;
    1624       11887 :   b[1] = evalsigne(1) | evalvarn(v);
    1625       78229 :   for (i=0; i<=n; i++) gel(x,i) = polcoef_i(b0, i, v);
    1626       11887 :   return b;
    1627             : }
    1628             : 
    1629             : /* assume varn(b) << varn(a) */
    1630             : /* return a FpY*/
    1631             : GEN
    1632           1 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
    1633             : {
    1634           1 :   long i,n,dres, db, vY = varn(b), vX = varn(a);
    1635             :   GEN la,x,y;
    1636             : 
    1637           1 :   if (lgefint(p) == 3)
    1638             :   {
    1639           0 :     ulong pp = uel(p,2);
    1640           0 :     b = ZXX_to_FlxX(b, pp, vX);
    1641           0 :     a = ZX_to_Flx(a, pp);
    1642           0 :     x = Flx_FlxY_resultant(a, b, pp);
    1643           0 :     return Flx_to_ZX(x);
    1644             :   }
    1645           1 :   db = RgXY_degreex(b);
    1646           1 :   dres = degpol(a)*degpol(b);
    1647           1 :   la = leading_coeff(a);
    1648           1 :   x = cgetg(dres+2, t_VEC);
    1649           1 :   y = cgetg(dres+2, t_VEC);
    1650             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1651             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1652           3 :   for (i=0,n = 1; i < dres; n++)
    1653             :   {
    1654           2 :     gel(x,++i) = utoipos(n);
    1655           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1656           2 :     gel(x,++i) = subiu(p,n);
    1657           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1658             :   }
    1659           1 :   if (i == dres)
    1660             :   {
    1661           0 :     gel(x,++i) = gen_0;
    1662           0 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
    1663             :   }
    1664           1 :   return FpV_polint(x,y, p, vY);
    1665             : }
    1666             : 
    1667             : static GEN
    1668          30 : FpX_composedsum(GEN P, GEN Q, GEN p)
    1669             : {
    1670          30 :   long n = 1+ degpol(P)*degpol(Q);
    1671          30 :   GEN Pl = FpX_invLaplace(FpX_Newton(P,n,p), p);
    1672          30 :   GEN Ql = FpX_invLaplace(FpX_Newton(Q,n,p), p);
    1673          30 :   GEN L = FpX_Laplace(FpXn_mul(Pl, Ql, n, p), p);
    1674          30 :   GEN lead = Fp_mul(Fp_powu(leading_coeff(P),degpol(Q), p),
    1675          30 :                     Fp_powu(leading_coeff(Q),degpol(P), p), p);
    1676          30 :   GEN R = FpX_fromNewton(L, p);
    1677          30 :   return FpX_Fp_mul(R, lead, p);
    1678             : }
    1679             : 
    1680             : #if 0
    1681             : GEN
    1682             : FpX_composedprod(GEN P, GEN Q, GEN p)
    1683             : {
    1684             :   long n = 1+ degpol(P)*degpol(Q);
    1685             :   GEN L=FpX_convol(FpX_Newton(P,n,p), FpX_Newton(Q,n,p), p);
    1686             :   return FpX_fromNewton(L, p);
    1687             : }
    1688             : #endif
    1689             : 
    1690             : GEN
    1691          30 : FpX_direct_compositum(GEN a, GEN b, GEN p)
    1692             : {
    1693          30 :   if (lgefint(p)==3)
    1694             :   {
    1695           0 :     pari_sp av = avma;
    1696           0 :     ulong pp = p[2];
    1697           0 :     GEN z = Flx_direct_compositum(ZX_to_Flx(a, pp), ZX_to_Flx(b, pp), pp);
    1698           0 :     return gerepileupto(av, Flx_to_ZX(z));
    1699             :   }
    1700          30 :   return FpX_composedsum(a, b, p);
    1701             : }
    1702             : 
    1703             : static GEN
    1704          30 : _FpX_direct_compositum(void *E, GEN a, GEN b)
    1705          30 : { return FpX_direct_compositum(a,b, (GEN)E); }
    1706             : 
    1707             : GEN
    1708         497 : FpXV_direct_compositum(GEN V, GEN p)
    1709             : {
    1710         497 :   if (lgefint(p)==3)
    1711             :   {
    1712           0 :     ulong pp = p[2];
    1713           0 :     return Flx_to_ZX(FlxV_direct_compositum(ZXV_to_FlxV(V, pp), pp));
    1714             :   }
    1715         497 :   return gen_product(V, (void *)p, &_FpX_direct_compositum);
    1716             : }
    1717             : 
    1718             : /* 0, 1, -1, 2, -2, ... */
    1719             : #define next_lambda(a) (a>0 ? -a : 1-a)
    1720             : GEN
    1721           0 : FpX_compositum(GEN a, GEN b, GEN p)
    1722             : {
    1723           0 :   long k, v = fetch_var_higher();
    1724           0 :   for (k = 1;; k = next_lambda(k))
    1725           0 :   {
    1726           0 :     GEN x = deg1pol_shallow(gen_1, gmulsg(k, pol_x(v)), 0); /* x + k y */
    1727           0 :     GEN C = FpX_FpXY_resultant(a, poleval(b,x),p);
    1728           0 :     if (FpX_is_squarefree(C, p)) { (void)delete_var(); return C; }
    1729             :   }
    1730             : }
    1731             : 
    1732             : /* Assume A in Z[Y], B in Q[Y][X], B squarefree in (Q[Y]/(A))[X] and
    1733             :  * Res_Y(A, B) in Z[X]. Find a small lambda (start from *lambda, use
    1734             :  * next_lambda successively) such that C(X) = Res_Y(A(Y), B(X + lambda Y))
    1735             :  * is squarefree, reset *lambda to the chosen value and return C. Set LERS to
    1736             :  * the Last non-constant polynomial in the Euclidean Remainder Sequence */
    1737             : static GEN
    1738        3871 : ZX_ZXY_resultant_LERS(GEN A, GEN B0, long *plambda, GEN *LERS)
    1739             : {
    1740             :   ulong bound, dp;
    1741        3871 :   pari_sp av = avma, av2 = 0;
    1742        3871 :   long lambda = *plambda, degA = degpol(A), dres = degA*degpol(B0);
    1743             :   long stable, checksqfree, i,n, cnt, degB;
    1744        3871 :   long v, vX = varn(B0), vY = varn(A); /* vY < vX */
    1745             :   GEN x, y, dglist, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
    1746             :   forprime_t S;
    1747             : 
    1748        3871 :   if (degA == 1)
    1749             :   {
    1750         630 :     GEN a1 = gel(A,3), a0 = gel(A,2);
    1751         630 :     B = lambda? RgX_translate(B0, monomial(stoi(lambda), 1, vY)): B0;
    1752         630 :     H = gsubst(B, vY, gdiv(gneg(a0),a1));
    1753         630 :    if (!equali1(a1)) H = RgX_Rg_mul(H, powiu(a1, poldegree(B,vY)));
    1754         630 :     *LERS = mkvec2(scalarpol_shallow(a0,vX), scalarpol_shallow(a1,vX));
    1755         630 :     gerepileall(av, 2, &H, LERS);
    1756         630 :     return H;
    1757             :   }
    1758             : 
    1759        3241 :   dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
    1760        3241 :   C0 = cgetg(dres+2, t_VECSMALL);
    1761        3241 :   C1 = cgetg(dres+2, t_VECSMALL);
    1762        3241 :   dglist = cgetg(dres+1, t_VECSMALL);
    1763        3241 :   x = cgetg(dres+2, t_VECSMALL);
    1764        3241 :   y = cgetg(dres+2, t_VECSMALL);
    1765        3241 :   B0 = leafcopy(B0);
    1766        3241 :   A = leafcopy(A);
    1767        3241 :   B = B0;
    1768        3241 :   v = fetch_var_higher(); setvarn(A,v);
    1769             :   /* make sure p large enough */
    1770        3829 : INIT:
    1771             :   /* always except the first time */
    1772        3829 :   if (av2) { set_avma(av2); lambda = next_lambda(lambda); }
    1773        3829 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    1774        3829 :   B = swap_vars(B, vY); setvarn(B,v);
    1775             :   /* B0(lambda v + x, v) */
    1776        3829 :   if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    1777        3829 :   av2 = avma;
    1778             : 
    1779        3829 :   if (degA <= 3)
    1780             :   { /* sub-resultant faster for small degrees */
    1781        2975 :     H = RgX_resultant_all(A,B,&q);
    1782        2975 :     if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
    1783        2520 :     H0 = gel(q,2);
    1784        2520 :     if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
    1785        2520 :     H1 = gel(q,3);
    1786        2520 :     if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
    1787        2520 :     if (!ZX_is_squarefree(H)) goto INIT;
    1788        2478 :     goto END;
    1789             :   }
    1790             : 
    1791         854 :   H = H0 = H1 = NULL;
    1792         854 :   degB = degpol(B);
    1793         854 :   bound = ZX_ZXY_ResBound(A, B, NULL);
    1794         854 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    1795         854 :   dp = 1;
    1796         854 :   init_modular_big(&S);
    1797         854 :   for(cnt = 0, checksqfree = 1;;)
    1798        2974 :   {
    1799        3828 :     ulong p = u_forprime_next(&S);
    1800             :     GEN Hi;
    1801        3828 :     a = ZX_to_Flx(A, p);
    1802        3828 :     b = ZXX_to_FlxX(B, p, varn(A));
    1803        3828 :     if (degpol(a) < degA || degpol(b) < degB) continue; /* p | lc(A)lc(B) */
    1804        3828 :     if (checksqfree)
    1805             :     { /* find degree list for generic Euclidean Remainder Sequence */
    1806         854 :       long goal = minss(degpol(a), degpol(b)); /* longest possible */
    1807        5355 :       for (n=1; n <= goal; n++) dglist[n] = 0;
    1808         854 :       setlg(dglist, 1);
    1809        3423 :       for (n=0; n <= dres; n++)
    1810             :       {
    1811        3255 :         ev = FlxY_evalx_drop(b, n, p);
    1812        3255 :         Flx_resultant_set_dglist(a, ev, dglist, p);
    1813        3255 :         if (lg(dglist)-1 == goal) break;
    1814             :       }
    1815             :       /* last pol in ERS has degree > 1 ? */
    1816         854 :       goal = lg(dglist)-1;
    1817         854 :       if (degpol(B) == 1) { if (!goal) goto INIT; }
    1818             :       else
    1819             :       {
    1820         805 :         if (goal <= 1) goto INIT;
    1821         756 :         if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
    1822             :       }
    1823         805 :       if (DEBUGLEVEL>4)
    1824           0 :         err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
    1825             :     }
    1826             : 
    1827       65744 :     for (i=0,n = 0; i <= dres; n++)
    1828             :     {
    1829       61965 :       ev = FlxY_evalx_drop(b, n, p);
    1830       61965 :       x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
    1831       61965 :       if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
    1832             :     }
    1833        3779 :     Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
    1834        3779 :     Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
    1835        3779 :     if (!H && degpol(Hp) != dres) continue;
    1836        3779 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1837        3779 :     if (checksqfree) {
    1838         805 :       if (!Flx_is_squarefree(Hp, p)) goto INIT;
    1839         763 :       if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    1840         763 :       checksqfree = 0;
    1841             :     }
    1842             : 
    1843        3737 :     if (!H)
    1844             :     { /* initialize */
    1845         763 :       q = utoipos(p); stable = 0;
    1846         763 :       H = ZX_init_CRT(Hp, p,vX);
    1847         763 :       H0= ZX_init_CRT(H0p, p,vX);
    1848         763 :       H1= ZX_init_CRT(H1p, p,vX);
    1849             :     }
    1850             :     else
    1851             :     {
    1852        2974 :       GEN qp = muliu(q,p);
    1853        2974 :       stable  = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
    1854        2974 :               & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
    1855        2974 :               & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
    1856        2974 :       q = qp;
    1857             :     }
    1858             :     /* could make it probabilistic for H ? [e.g if stable twice, etc]
    1859             :      * Probabilistic anyway for H0, H1 */
    1860        3737 :     if (DEBUGLEVEL>5 && (stable ||  ++cnt==100))
    1861           0 :     { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
    1862        3737 :     if (stable && (ulong)expi(q) >= bound) break; /* DONE */
    1863        2974 :     if (gc_needed(av,2))
    1864             :     {
    1865           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
    1866           0 :       gerepileall(av2, 4, &H, &q, &H0, &H1);
    1867             :     }
    1868             :   }
    1869        3241 : END:
    1870        3241 :   if (DEBUGLEVEL>5) err_printf(" done\n");
    1871        3241 :   setvarn(H, vX); (void)delete_var();
    1872        3241 :   *LERS = mkvec2(H0,H1);
    1873        3241 :   gerepileall(av, 2, &H, LERS);
    1874        3241 :   *plambda = lambda; return H;
    1875             : }
    1876             : 
    1877             : GEN
    1878        5012 : ZX_ZXY_resultant_all(GEN A, GEN B, long *plambda, GEN *LERS)
    1879             : {
    1880        5012 :   if (LERS)
    1881             :   {
    1882        3871 :     if (!plambda)
    1883           0 :       pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
    1884        3871 :     return ZX_ZXY_resultant_LERS(A, B, plambda, LERS);
    1885             :   }
    1886        1141 :   return ZX_ZXY_rnfequation(A, B, plambda);
    1887             : }
    1888             : 
    1889             : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
    1890             :  * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
    1891             :  * squarefree */
    1892             : GEN
    1893        1961 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
    1894             : {
    1895        1961 :   pari_sp av = avma;
    1896             :   GEN R, a;
    1897             :   long dA;
    1898             :   int delvar;
    1899             : 
    1900        1961 :   if (v < 0) v = 0;
    1901        1961 :   switch (typ(A))
    1902             :   {
    1903        1961 :     case t_POL: dA = degpol(A); if (dA > 0) break;
    1904           0 :       A = constant_coeff(A);
    1905           0 :     default:
    1906           0 :       if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
    1907           0 :       return gerepileupto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
    1908             :   }
    1909        1961 :   delvar = 0;
    1910        1961 :   if (varn(T) == 0)
    1911             :   {
    1912        1870 :     long v0 = fetch_var(); delvar = 1;
    1913        1870 :     T = leafcopy(T); setvarn(T,v0);
    1914        1870 :     A = leafcopy(A); setvarn(A,v0);
    1915             :   }
    1916        1961 :   R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
    1917        1961 :   if (delvar) (void)delete_var();
    1918        1961 :   setvarn(R, v); a = leading_coeff(T);
    1919        1961 :   if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
    1920        1961 :   return gerepileupto(av, R);
    1921             : }
    1922             : 
    1923             : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
    1924             : GEN
    1925       14927 : ZXQ_charpoly(GEN A, GEN T, long v)
    1926             : {
    1927       14927 :   return (degpol(T) < 16) ? RgXQ_charpoly(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
    1928             : }
    1929             : 
    1930             : GEN
    1931         812 : QXQ_charpoly(GEN A, GEN T, long v)
    1932             : {
    1933         812 :   pari_sp av = avma;
    1934         812 :   GEN den, B = Q_remove_denom(A, &den);
    1935         812 :   GEN P = ZXQ_charpoly(B, T, v);
    1936         812 :   return gerepilecopy(av, den ? RgX_rescale(P, ginv(den)): P);
    1937             : }
    1938             : 
    1939             : static ulong
    1940     1610381 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
    1941             : {
    1942     1610381 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    1943             :   ulong H, dp;
    1944     1610209 :   if (dropa && dropb) return 0; /* p | lc(A), p | lc(B) */
    1945     1610209 :   H = Flx_resultant(a, b, p);
    1946     1609941 :   if (dropa)
    1947             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    1948           0 :     ulong c = b[degB+2]; /* lc(B) */
    1949           0 :     if (odd(degB)) c = p - c;
    1950           0 :     c = Fl_powu(c, dropa, p);
    1951           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1952             :   }
    1953     1609941 :   else if (dropb)
    1954             :   { /* multiply by lc(A)^(deg B - deg b) */
    1955           0 :     ulong c = a[degA+2]; /* lc(A) */
    1956           0 :     c = Fl_powu(c, dropb, p);
    1957           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1958             :   }
    1959     1609938 :   dp = dB ? umodiu(dB, p): 1;
    1960     1609938 :   if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1961     1609937 :   return H;
    1962             : }
    1963             : 
    1964             : /* If B=NULL, assume B=A' */
    1965             : static GEN
    1966      273110 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
    1967             : {
    1968      273110 :   pari_sp av = avma, av2;
    1969      273110 :   long degA, degB, i, n = lg(P)-1;
    1970             :   GEN H, T;
    1971             : 
    1972      273110 :   degA = degpol(A);
    1973      273110 :   degB = B? degpol(B): degA - 1;
    1974      273110 :   if (n == 1)
    1975             :   {
    1976       43302 :     ulong Hp, p = uel(P,1);
    1977       43302 :     GEN a = ZX_to_Flx(A, p), b = B? ZX_to_Flx(B, p): Flx_deriv(a, p);
    1978       43302 :     Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1979       43302 :     set_avma(av); *mod = utoipos(p); return utoi(Hp);
    1980             :   }
    1981      229808 :   T = ZV_producttree(P);
    1982      229810 :   A = ZX_nv_mod_tree(A, P, T);
    1983      229781 :   if (B) B = ZX_nv_mod_tree(B, P, T);
    1984      229781 :   H = cgetg(n+1, t_VECSMALL); av2 = avma;
    1985     1796541 :   for(i=1; i <= n; i++, set_avma(av2))
    1986             :   {
    1987     1566754 :     ulong p = P[i];
    1988     1566754 :     GEN a = gel(A,i), b = B? gel(B,i): Flx_deriv(a, p);
    1989     1567111 :     H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1990             :   }
    1991      229787 :   H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
    1992      229801 :   *mod = gmael(T, lg(T)-1, 1);
    1993      229801 :   gerepileall(av, 2, &H, mod);
    1994      229806 :   return H;
    1995             : }
    1996             : 
    1997             : GEN
    1998      273110 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
    1999             : {
    2000      273110 :   GEN V = cgetg(3, t_VEC);
    2001      273110 :   if (typ(B) == t_INT) B = NULL;
    2002      273110 :   if (!signe(dB)) dB = NULL;
    2003      273110 :   gel(V,1) = ZX_resultant_slice(A, B, dB, P, &gel(V,2));
    2004      273108 :   return V;
    2005             : }
    2006             : 
    2007             : /* Compute Res(A, B/dB) in Z, assuming A,B in Z[X], dB in Z or NULL (= 1)
    2008             :  * If B=NULL, take B = A' and assume deg A > 1 and 'bound' is set */
    2009             : GEN
    2010      107888 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
    2011             : {
    2012      107888 :   pari_sp av = avma;
    2013             :   forprime_t S;
    2014             :   GEN  H, worker;
    2015      107888 :   if (B)
    2016             :   {
    2017       43623 :     long a = degpol(A), b = degpol(B);
    2018       43623 :     if (a < 0 || b < 0) return gen_0;
    2019       43592 :     if (!a) return powiu(gel(A,2), b);
    2020       43592 :     if (!b) return powiu(gel(B,2), a);
    2021       42441 :     if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
    2022             :   }
    2023      106706 :   worker = snm_closure(is_entry("_ZX_resultant_worker"),
    2024             :                        mkvec3(A, B? B: gen_0, dB? dB: gen_0));
    2025      106706 :   init_modular_big(&S);
    2026      106705 :   H = gen_crt("ZX_resultant_all", worker, &S, dB, bound, 0, NULL,
    2027             :               ZV_chinese_center, Fp_center);
    2028      106706 :   return gerepileuptoint(av, H);
    2029             : }
    2030             : 
    2031             : /* A0 and B0 in Q[X] */
    2032             : GEN
    2033          56 : QX_resultant(GEN A0, GEN B0)
    2034             : {
    2035             :   GEN s, a, b, A, B;
    2036          56 :   pari_sp av = avma;
    2037             : 
    2038          56 :   A = Q_primitive_part(A0, &a);
    2039          56 :   B = Q_primitive_part(B0, &b);
    2040          56 :   s = ZX_resultant(A, B);
    2041          56 :   if (!signe(s)) { set_avma(av); return gen_0; }
    2042          56 :   if (a) s = gmul(s, gpowgs(a,degpol(B)));
    2043          56 :   if (b) s = gmul(s, gpowgs(b,degpol(A)));
    2044          56 :   return gerepileupto(av, s);
    2045             : }
    2046             : 
    2047             : GEN
    2048        2982 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
    2049             : 
    2050             : GEN
    2051           0 : QXQ_intnorm(GEN A, GEN B)
    2052             : {
    2053             :   GEN c, n, R, lB;
    2054           0 :   long dA = degpol(A), dB = degpol(B);
    2055           0 :   pari_sp av = avma;
    2056           0 :   if (dA < 0) return gen_0;
    2057           0 :   A = Q_primitive_part(A, &c);
    2058           0 :   if (!c || typ(c) == t_INT) {
    2059           0 :     n = c;
    2060           0 :     R = ZX_resultant(B, A);
    2061             :   } else {
    2062           0 :     n = gel(c,1);
    2063           0 :     R = ZX_resultant_all(B, A, gel(c,2), 0);
    2064             :   }
    2065           0 :   if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
    2066           0 :   lB = leading_coeff(B);
    2067           0 :   if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
    2068           0 :   return gerepileuptoint(av, R);
    2069             : }
    2070             : 
    2071             : GEN
    2072           0 : QXQ_norm(GEN A, GEN B)
    2073             : {
    2074             :   GEN c, R, lB;
    2075           0 :   long dA = degpol(A), dB = degpol(B);
    2076           0 :   pari_sp av = avma;
    2077           0 :   if (dA < 0) return gen_0;
    2078           0 :   A = Q_primitive_part(A, &c);
    2079           0 :   R = ZX_resultant(B, A);
    2080           0 :   if (c) R = gmul(R, gpowgs(c, dB));
    2081           0 :   lB = leading_coeff(B);
    2082           0 :   if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
    2083           0 :   return gerepileupto(av, R);
    2084             : }
    2085             : 
    2086             : /* assume x has integral coefficients */
    2087             : GEN
    2088       66039 : ZX_disc_all(GEN x, ulong bound)
    2089             : {
    2090       66039 :   pari_sp av = avma;
    2091       66039 :   long s, d = degpol(x);
    2092             :   GEN l, R;
    2093             : 
    2094       66039 :   if (d <= 1) return d == 1? gen_1: gen_0;
    2095       64265 :   s = (d & 2) ? -1: 1;
    2096       64265 :   l = leading_coeff(x);
    2097       64265 :   if (!bound) bound = ZX_discbound(x);
    2098       64265 :   R = ZX_resultant_all(x, NULL, NULL, bound);
    2099       64265 :   if (is_pm1(l))
    2100       61402 :   { if (signe(l) < 0) s = -s; }
    2101             :   else
    2102        2863 :     R = diviiexact(R,l);
    2103       64265 :   if (s == -1) togglesign_safe(&R);
    2104       64265 :   return gerepileuptoint(av,R);
    2105             : }
    2106             : 
    2107             : GEN
    2108       62350 : ZX_disc(GEN x) { return ZX_disc_all(x,0); }
    2109             : 
    2110             : static GEN
    2111        5998 : ZXQX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, GEN T, ulong p)
    2112             : {
    2113        5998 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    2114             :   GEN H, dp;
    2115        5998 :   if (dropa && dropb) return pol0_Flx(T[1]); /* p | lc(A), p | lc(B) */
    2116        5998 :   H = FlxqX_saferesultant(a, b, T, p);
    2117        5999 :   if (!H) return NULL;
    2118        5999 :   if (dropa)
    2119             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2120           0 :     GEN c = gel(b,degB+2); /* lc(B) */
    2121           0 :     if (odd(degB)) c = Flx_neg(c, p);
    2122           0 :     c = Flxq_powu(c, dropa, T, p);
    2123           0 :     if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
    2124             :   }
    2125        5999 :   else if (dropb)
    2126             :   { /* multiply by lc(A)^(deg B - deg b) */
    2127           0 :     GEN c = gel(a,degA+2); /* lc(A) */
    2128           0 :     c = Flxq_powu(c, dropb, T, p);
    2129           0 :     if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
    2130             :   }
    2131        5999 :   dp = dB ? ZX_to_Flx(dB, p): pol1_Flx(T[1]);
    2132        5999 :   if (!Flx_equal1(dp))
    2133             :   {
    2134           0 :     GEN idp = Flxq_invsafe(dp, T, p);
    2135           0 :     if (!idp) return NULL;
    2136           0 :     H = Flxq_mul(H, Flxq_powu(idp, degA, T, p), T, p);
    2137             :   }
    2138        5999 :   return H;
    2139             : }
    2140             : 
    2141             : /* If B=NULL, assume B=A' */
    2142             : static GEN
    2143        2894 : ZXQX_resultant_slice(GEN A, GEN B, GEN U, GEN dB, GEN P, GEN *mod)
    2144             : {
    2145        2894 :   pari_sp av = avma;
    2146        2894 :   long degA, degB, i, n = lg(P)-1;
    2147             :   GEN H, T;
    2148        2894 :   long v = varn(U), redo = 0;
    2149             : 
    2150        2894 :   degA = degpol(A);
    2151        2894 :   degB = B? degpol(B): degA - 1;
    2152        2894 :   if (n == 1)
    2153             :   {
    2154        1856 :     ulong p = uel(P,1);
    2155        1856 :     GEN a = ZXX_to_FlxX(A, p, v), b = B? ZXX_to_FlxX(B, p, v): FlxX_deriv(a, p);
    2156        1856 :     GEN u = ZX_to_Flx(U, p);
    2157        1856 :     GEN Hp = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
    2158        1856 :     if (!Hp) { set_avma(av); *mod = gen_1; return pol_0(v); }
    2159        1856 :     set_avma(av); *mod = utoipos(p); return Flx_to_ZX(Hp);
    2160             :   }
    2161        1038 :   T = ZV_producttree(P);
    2162        1038 :   A = ZXX_nv_mod_tree(A, P, T, v);
    2163        1038 :   if (B) B = ZXX_nv_mod_tree(B, P, T, v);
    2164        1038 :   U = ZX_nv_mod_tree(U, P, T);
    2165        1038 :   H = cgetg(n+1, t_VEC);
    2166        5181 :   for(i=1; i <= n; i++)
    2167             :   {
    2168        4143 :     ulong p = P[i];
    2169        4143 :     GEN a = gel(A,i), b = B? gel(B,i): FlxX_deriv(a, p), u = gel(U, i);
    2170        4142 :     GEN h = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
    2171        4143 :     if (!h)
    2172             :     {
    2173           0 :       gel(H,i) = pol_0(v);
    2174           0 :       P[i] = 1; redo = 1;
    2175             :     }
    2176             :     else
    2177        4143 :       gel(H,i) = h;
    2178             :   }
    2179        1038 :   if (redo) T = ZV_producttree(P);
    2180        1038 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2181        1038 :   *mod = gmael(T, lg(T)-1, 1);
    2182        1038 :   gerepileall(av, 2, &H, mod);
    2183        1038 :   return H;
    2184             : }
    2185             : 
    2186             : GEN
    2187        2894 : ZXQX_resultant_worker(GEN P, GEN A, GEN B, GEN T, GEN dB)
    2188             : {
    2189        2894 :   GEN V = cgetg(3, t_VEC);
    2190        2894 :   if (isintzero(B)) B = NULL;
    2191        2894 :   if (!signe(dB)) dB = NULL;
    2192        2894 :   gel(V,1) = ZXQX_resultant_slice(A, B, T, dB, P, &gel(V,2));
    2193        2894 :   return V;
    2194             : }
    2195             : 
    2196             : static ulong
    2197        2411 : ZXQX_resultant_bound(GEN nf, GEN A, GEN B)
    2198             : {
    2199        2411 :   pari_sp av = avma;
    2200        2411 :   GEN r, M = L2_bound(nf, NULL, &r);
    2201        2411 :   long v = nf_get_varn(nf), i, l = lg(r);
    2202        2411 :   GEN a = cgetg(l, t_COL);
    2203        6887 :   for (i = 1; i < l; i++)
    2204        4476 :     gel(a, i) =  RgX_RgXY_ResBound(gsubst(A, v, gel(r,i)),
    2205        4476 :                                    gsubst(B, v, gel(r,i)), DEFAULTPREC);
    2206        2411 :   return gc_ulong(av, (ulong) dbllog2(gmul(M,RgC_fpnorml2(a, DEFAULTPREC))));
    2207             : }
    2208             : 
    2209             : /* Compute Res(A, B/dB) in Z, assuming A,B in Z[X], dB in Z or NULL (= 1)
    2210             :  * If B=NULL, take B = A' and assume deg A > 1 and 'bound' is set */
    2211             : static GEN
    2212        2408 : ZXQX_resultant_all(GEN A, GEN B, GEN T, GEN dB, ulong bound)
    2213             : {
    2214        2408 :   pari_sp av = avma;
    2215             :   forprime_t S;
    2216             :   GEN  H, worker;
    2217        2408 :   if (B)
    2218             :   {
    2219          70 :     long a = degpol(A), b = degpol(B);
    2220          70 :     if (a < 0 || b < 0) return gen_0;
    2221          70 :     if (!a) return gpowgs(gel(A,2), b);
    2222          70 :     if (!b) return gpowgs(gel(B,2), a);
    2223          49 :     if (!bound) bound = ZXQX_resultant_bound(nfinit(T, DEFAULTPREC), A, B);
    2224             :   } else
    2225        2338 :     if (!bound) bound = ZXQX_resultant_bound(nfinit(T, DEFAULTPREC), A, RgX_deriv(A));
    2226        2387 :   worker = snm_closure(is_entry("_ZXQX_resultant_worker"),
    2227             :                        mkvec4(A, B? B: gen_0, T, dB? dB: gen_0));
    2228        2387 :   init_modular_big(&S);
    2229        2387 :   H = gen_crt("ZXQX_resultant_all", worker, &S, dB, bound, 0, NULL,
    2230             :               nxV_chinese_center, FpX_center);
    2231        2387 :   if (DEBUGLEVEL)
    2232           0 :     err_printf("ZXQX_resultant_all: a priori bound: %lu, a posteriori: %lu\n",
    2233             :                bound, expi(gsupnorm(H, DEFAULTPREC)));
    2234        2387 :   return gerepileupto(av, H);
    2235             : }
    2236             : 
    2237             : GEN
    2238          70 : nfX_resultant(GEN nf, GEN x, GEN y)
    2239             : {
    2240          70 :   pari_sp av = avma;
    2241          70 :   GEN cx, cy, D, T = nf_get_pol(nf);
    2242             :   ulong bound;
    2243          70 :   long d = degpol(x), v = varn(T);
    2244          70 :   if (d <= 1) return d == 1? pol_1(v): pol_0(v);
    2245          70 :   x = Q_primitive_part(x, &cx);
    2246          70 :   y = Q_primitive_part(y, &cy);
    2247          70 :   bound = ZXQX_resultant_bound(nf, x, y);
    2248          70 :   D = ZXQX_resultant_all(x, y, T, NULL, bound);
    2249          70 :   if (cx) D = gmul(D, gpowgs(cx, degpol(y)));
    2250          70 :   if (cy) D = gmul(D, gpowgs(cy, degpol(x)));
    2251          70 :   return gerepileupto(av, D);
    2252             : }
    2253             : 
    2254             : static GEN
    2255         147 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
    2256             : 
    2257             : static GEN
    2258        2338 : ZXQX_disc_all(GEN x, GEN T, ulong bound)
    2259             : {
    2260        2338 :   pari_sp av = avma;
    2261        2338 :   long s, d = degpol(x), v = varn(T);
    2262             :   GEN l, R;
    2263             : 
    2264        2338 :   if (d <= 1) return d == 1? pol_1(v): pol_0(v);
    2265        2338 :   s = (d & 2) ? -1: 1;
    2266        2338 :   l = leading_coeff(x);
    2267        2338 :   R = ZXQX_resultant_all(x, NULL, T, NULL, bound);
    2268        2338 :   if (!gequal1(l)) R = QXQ_div(R, to_ZX(l,v), T);
    2269        2338 :   if (s == -1) R = RgX_neg(R);
    2270        2338 :   return gerepileupto(av, R);
    2271             : }
    2272             : 
    2273             : GEN
    2274           7 : QX_disc(GEN x)
    2275             : {
    2276           7 :   pari_sp av = avma;
    2277           7 :   GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
    2278           7 :   if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
    2279           7 :   return gerepileupto(av, d);
    2280             : }
    2281             : 
    2282             : GEN
    2283        2471 : nfX_disc(GEN nf, GEN x)
    2284             : {
    2285        2471 :   pari_sp av = avma;
    2286        2471 :   GEN c, D, T = nf_get_pol(nf);
    2287             :   ulong bound;
    2288        2471 :   long d = degpol(x), v = varn(T);
    2289        2471 :   if (d <= 1) return d == 1? pol_1(v): pol_0(v);
    2290        2338 :   x = Q_primitive_part(x, &c);
    2291        2338 :   bound = ZXQX_resultant_bound(nf, x, RgX_deriv(x));
    2292        2338 :   D = ZXQX_disc_all(x, T, bound);
    2293        2338 :   if (c) D = gmul(D, gpowgs(c, 2*d - 2));
    2294        2338 :   return gerepileupto(av, D);
    2295             : }
    2296             : 
    2297             : GEN
    2298       63487 : QXQ_mul(GEN x, GEN y, GEN T)
    2299             : {
    2300       63487 :   GEN dx, nx = Q_primitive_part(x, &dx);
    2301       63487 :   GEN dy, ny = Q_primitive_part(y, &dy);
    2302       63487 :   GEN z = ZXQ_mul(nx, ny, T);
    2303       63487 :   if (dx || dy)
    2304             :   {
    2305       63417 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
    2306       63417 :     if (!gequal1(d)) z = ZX_Q_mul(z, d);
    2307             :   }
    2308       63487 :   return z;
    2309             : }
    2310             : 
    2311             : GEN
    2312       11565 : QXQ_sqr(GEN x, GEN T)
    2313             : {
    2314       11565 :   GEN dx, nx = Q_primitive_part(x, &dx);
    2315       11565 :   GEN z = ZXQ_sqr(nx, T);
    2316       11565 :   if (dx)
    2317       11565 :     z = ZX_Q_mul(z, gsqr(dx));
    2318       11565 :   return z;
    2319             : }
    2320             : 
    2321             : static GEN
    2322       47616 : QXQ_inv_slice(GEN A, GEN B, GEN P, GEN *mod)
    2323             : {
    2324       47616 :   pari_sp av = avma;
    2325       47616 :   long i, n = lg(P)-1, v = varn(A), redo = 0;
    2326             :   GEN H, T;
    2327       47616 :   if (n == 1)
    2328             :   {
    2329       35093 :     ulong p = uel(P,1);
    2330       35093 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
    2331       35093 :     GEN U = Flxq_invsafe(a, b, p);
    2332       35093 :     if (!U)
    2333             :     {
    2334          24 :       set_avma(av);
    2335          24 :       *mod = gen_1; return pol_0(v);
    2336             :     }
    2337       35069 :     H = gerepilecopy(av, Flx_to_ZX(U));
    2338       35069 :     *mod = utoi(p);
    2339       35069 :     return H;
    2340             :   }
    2341       12523 :   T = ZV_producttree(P);
    2342       12524 :   A = ZX_nv_mod_tree(A, P, T);
    2343       12521 :   B = ZX_nv_mod_tree(B, P, T);
    2344       12522 :   H = cgetg(n+1, t_VEC);
    2345       67667 :   for(i=1; i <= n; i++)
    2346             :   {
    2347       55145 :     ulong p = P[i];
    2348       55145 :     GEN a = gel(A,i), b = gel(B,i);
    2349       55145 :     GEN U = Flxq_invsafe(a, b, p);
    2350       55129 :     if (!U)
    2351             :     {
    2352         603 :       gel(H,i) = pol_0(v);
    2353         600 :       P[i] = 1; redo = 1;
    2354             :     }
    2355             :     else
    2356       54526 :       gel(H,i) = U;
    2357             :   }
    2358       12522 :   if (redo) T = ZV_producttree(P);
    2359       12522 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2360       12524 :   *mod = gmael(T, lg(T)-1, 1);
    2361       12524 :   gerepileall(av, 2, &H, mod);
    2362       12524 :   return H;
    2363             : }
    2364             : 
    2365             : GEN
    2366       47616 : QXQ_inv_worker(GEN P, GEN A, GEN B)
    2367             : {
    2368       47616 :   GEN V = cgetg(3, t_VEC);
    2369       47616 :   gel(V,1) = QXQ_inv_slice(A, B, P, &gel(V,2));
    2370       47617 :   return V;
    2371             : }
    2372             : 
    2373             : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
    2374             : GEN
    2375       27652 : QXQ_inv(GEN A, GEN B)
    2376             : {
    2377             :   GEN D, Ap, Bp;
    2378             :   ulong pp;
    2379       27652 :   pari_sp av2, av = avma;
    2380             :   forprime_t S;
    2381       27652 :   GEN worker, U, H = NULL, mod = gen_1;
    2382             :   pari_timer ti;
    2383             :   long k, dA, dB;
    2384       27652 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2385             :   /* A a QX, B a ZX */
    2386       27652 :   A = Q_primitive_part(A, &D);
    2387       27652 :   dA = degpol(A); dB= degpol(B);
    2388             :   /* A, B in Z[X] */
    2389       27652 :   init_modular_small(&S);
    2390             :   do {
    2391       27652 :     pp = u_forprime_next(&S);
    2392       27652 :     Ap = ZX_to_Flx(A, pp);
    2393       27652 :     Bp = ZX_to_Flx(B, pp);
    2394       27652 :   } while (degpol(Ap) != dA || degpol(Bp) != dB);
    2395       27652 :   if (degpol(Flx_gcd(Ap, Bp, pp)) != 0 && degpol(ZX_gcd(A,B))!=0)
    2396          14 :     pari_err_INV("QXQ_inv",mkpolmod(A,B));
    2397       27638 :   worker = snm_closure(is_entry("_QXQ_inv_worker"), mkvec2(A, B));
    2398       27638 :   av2 = avma;
    2399       27638 :   for (k = 1; ;k *= 2)
    2400       14082 :   {
    2401             :     GEN res, b, N, den;
    2402       41720 :     gen_inccrt_i("QXQ_inv", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
    2403             :                  nxV_chinese_center, FpX_center);
    2404       41720 :     gerepileall(av2, 2, &H, &mod);
    2405       41720 :     b = sqrti(shifti(mod,-1));
    2406       41720 :     if (DEBUGLEVEL>5) timer_start(&ti);
    2407       41720 :     U = FpX_ratlift(H, mod, b, b, NULL);
    2408       41720 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: ratlift");
    2409       43196 :     if (!U) continue;
    2410       29114 :     N = Q_remove_denom(U, &den); if (!den) den = gen_1;
    2411       29114 :     res = Flx_rem(Flx_Fl_sub(Flx_mul(Ap, ZX_to_Flx(N,pp), pp),
    2412             :                   umodiu(den, pp), pp), Bp, pp);
    2413       29114 :     if (degpol(res) >= 0) continue;
    2414       27638 :     res = ZX_Z_sub(ZX_mul(A, N), den);
    2415       27638 :     res = ZX_is_monic(B) ? ZX_rem(res, B): RgX_pseudorem(res, B);
    2416       27638 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: final check");
    2417       27638 :     if (degpol(res)<0)
    2418             :     {
    2419       27638 :       if (D) U = RgX_Rg_div(U, D);
    2420       27638 :       return gerepilecopy(av, U);
    2421             :     }
    2422             :   }
    2423             : }
    2424             : 
    2425             : static GEN
    2426        7479 : QXQ_div_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
    2427             : {
    2428        7479 :   pari_sp av = avma;
    2429        7479 :   long i, n = lg(P)-1, v = varn(A), redo = 0;
    2430             :   GEN H, T;
    2431        7479 :   if (n == 1)
    2432             :   {
    2433        4710 :     ulong p = uel(P,1);
    2434        4710 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p), c = ZX_to_Flx(C, p);
    2435        4710 :     GEN bi = Flxq_invsafe(b, c, p), U;
    2436        4710 :     if (!bi)
    2437             :     {
    2438           0 :       set_avma(av);
    2439           0 :       *mod = gen_1; return pol_0(v);
    2440             :     }
    2441        4710 :     U = Flxq_mul(a, bi, c, p);
    2442        4710 :     H = gerepilecopy(av, Flx_to_ZX(U));
    2443        4710 :     *mod = utoi(p);
    2444        4710 :     return H;
    2445             :   }
    2446        2769 :   T = ZV_producttree(P);
    2447        2769 :   A = ZX_nv_mod_tree(A, P, T);
    2448        2769 :   B = ZX_nv_mod_tree(B, P, T);
    2449        2769 :   C = ZX_nv_mod_tree(C, P, T);
    2450        2769 :   H = cgetg(n+1, t_VEC);
    2451       10306 :   for(i=1; i <= n; i++)
    2452             :   {
    2453        7537 :     ulong p = P[i];
    2454        7537 :     GEN a = gel(A,i), b = gel(B,i), c = gel(C, i);
    2455        7537 :     GEN bi = Flxq_invsafe(b, c, p);
    2456        7538 :     if (!bi)
    2457             :     {
    2458           0 :       gel(H,i) = pol_0(v);
    2459           0 :       P[i] = 1; redo = 1;
    2460             :     }
    2461             :     else
    2462        7538 :       gel(H,i) = Flxq_mul(a, bi, c, p);
    2463             :   }
    2464        2769 :   if (redo) T = ZV_producttree(P);
    2465        2769 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2466        2769 :   *mod = gmael(T, lg(T)-1, 1);
    2467        2769 :   gerepileall(av, 2, &H, mod);
    2468        2769 :   return H;
    2469             : }
    2470             : 
    2471             : GEN
    2472        7479 : QXQ_div_worker(GEN P, GEN A, GEN B, GEN C)
    2473             : {
    2474        7479 :   GEN V = cgetg(3, t_VEC);
    2475        7479 :   gel(V,1) = QXQ_div_slice(A, B, C, P, &gel(V,2));
    2476        7479 :   return V;
    2477             : }
    2478             : 
    2479             : /* lift(Mod(A/B, C)). C a ZX, A, B a scalar or a QX */
    2480             : GEN
    2481        3640 : QXQ_div(GEN A, GEN B, GEN C)
    2482             : {
    2483             :   GEN DA, DB, Ap, Bp, Cp;
    2484             :   ulong pp;
    2485        3640 :   pari_sp av2, av = avma;
    2486             :   forprime_t S;
    2487        3640 :   GEN worker, U, H = NULL, mod = gen_1;
    2488             :   pari_timer ti;
    2489             :   long k, dA, dB, dC;
    2490        3640 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2491             :   /* A a QX, B a ZX */
    2492        3640 :   A = Q_primitive_part(A, &DA);
    2493        3640 :   B = Q_primitive_part(B, &DB);
    2494        3640 :   dA = degpol(A); dB = degpol(B); dC = degpol(C);
    2495             :   /* A, B in Z[X] */
    2496        3640 :   init_modular_small(&S);
    2497             :   do {
    2498        3640 :     pp = u_forprime_next(&S);
    2499        3640 :     Ap = ZX_to_Flx(A, pp);
    2500        3640 :     Bp = ZX_to_Flx(B, pp);
    2501        3640 :     Cp = ZX_to_Flx(C, pp);
    2502        3640 :   } while (degpol(Ap) != dA || degpol(Bp) != dB || degpol(Cp) != dC);
    2503        3640 :   if (degpol(Flx_gcd(Bp, Cp, pp)) != 0 && degpol(ZX_gcd(B,C))!=0)
    2504           0 :     pari_err_INV("QXQ_div",mkpolmod(B,C));
    2505        3640 :   worker = snm_closure(is_entry("_QXQ_div_worker"), mkvec3(A, B, C));
    2506        3640 :   av2 = avma;
    2507        3640 :   for (k = 1; ;k *= 2)
    2508        2724 :   {
    2509             :     GEN res, b, N, den;
    2510        6364 :     gen_inccrt_i("QXQ_div", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
    2511             :                  nxV_chinese_center, FpX_center);
    2512        6364 :     gerepileall(av2, 2, &H, &mod);
    2513        6364 :     b = sqrti(shifti(mod,-1));
    2514        6364 :     if (DEBUGLEVEL>5) timer_start(&ti);
    2515        6364 :     U = FpX_ratlift(H, mod, b, b, NULL);
    2516        6364 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: ratlift");
    2517        6506 :     if (!U) continue;
    2518        3782 :     N = Q_remove_denom(U, &den); if (!den) den = gen_1;
    2519        3782 :     res = Flx_rem(Flx_sub(Flx_mul(Bp, ZX_to_Flx(N,pp), pp),
    2520             :                           Flx_Fl_mul(Ap, umodiu(den, pp), pp), pp), Cp, pp);
    2521        3782 :     if (degpol(res) >= 0) continue;
    2522        3640 :     res = ZX_sub(ZX_mul(B, N), ZX_Z_mul(A,den));
    2523        3640 :     res = ZX_is_monic(C) ? ZX_rem(res, C): RgX_pseudorem(res, C);
    2524        3640 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: final check");
    2525        3640 :     if (degpol(res)<0)
    2526             :     {
    2527        3640 :       if (DA && DB) U = RgX_Rg_mul(U, gdiv(DA,DB));
    2528        2765 :       else if (DA) U = RgX_Rg_mul(U, DA);
    2529        2233 :       else if (DB) U = RgX_Rg_div(U, DB);
    2530        3640 :       return gerepilecopy(av, U);
    2531             :     }
    2532             :   }
    2533             : }
    2534             : 
    2535             : /************************************************************************
    2536             :  *                                                                      *
    2537             :  *                           ZXQ_minpoly                                *
    2538             :  *                                                                      *
    2539             :  ************************************************************************/
    2540             : 
    2541             : static GEN
    2542        3524 : ZXQ_minpoly_slice(GEN A, GEN B, long d, GEN P, GEN *mod)
    2543             : {
    2544        3524 :   pari_sp av = avma;
    2545        3524 :   long i, n = lg(P)-1, v = evalvarn(varn(B));
    2546             :   GEN H, T;
    2547        3524 :   if (n == 1)
    2548             :   {
    2549         483 :     ulong p = uel(P,1);
    2550         483 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
    2551         483 :     GEN Hp = Flxq_minpoly(a, b, p);
    2552         483 :     if (degpol(Hp) != d) { p = 1; Hp = pol0_Flx(v); }
    2553         483 :     H = Flx_to_ZX(Hp);
    2554         483 :     *mod = utoi(p);
    2555         483 :     gerepileall(av, 2, &H, mod);
    2556         483 :     return H;
    2557             :   }
    2558        3041 :   T = ZV_producttree(P);
    2559        3041 :   A = ZX_nv_mod_tree(A, P, T);
    2560        3041 :   B = ZX_nv_mod_tree(B, P, T);
    2561        3041 :   H = cgetg(n+1, t_VEC);
    2562       18615 :   for(i=1; i <= n; i++)
    2563             :   {
    2564       15574 :     ulong p = P[i];
    2565       15574 :     GEN a = gel(A,i), b = gel(B,i);
    2566       15574 :     GEN m = Flxq_minpoly(a, b, p);
    2567       15574 :     if (degpol(m) != d) { P[i] = 1; m = pol0_Flx(v); }
    2568       15574 :     gel(H, i) = m;
    2569             :   }
    2570        3041 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2571        3041 :   *mod = gmael(T, lg(T)-1, 1);
    2572        3041 :   gerepileall(av, 2, &H, mod);
    2573        3041 :   return H;
    2574             : }
    2575             : 
    2576             : GEN
    2577        3524 : ZXQ_minpoly_worker(GEN P, GEN A, GEN B, long d)
    2578             : {
    2579        3524 :   GEN V = cgetg(3, t_VEC);
    2580        3524 :   gel(V,1) = ZXQ_minpoly_slice(A, B, d, P, &gel(V,2));
    2581        3524 :   return V;
    2582             : }
    2583             : 
    2584             : GEN
    2585        1596 : ZXQ_minpoly(GEN A, GEN B, long d, ulong bound)
    2586             : {
    2587        1596 :   pari_sp av = avma;
    2588             :   GEN worker, H, dB;
    2589             :   forprime_t S;
    2590        1596 :   B = Q_remove_denom(B, &dB);
    2591        1596 :   worker = strtoclosure("_ZXQ_minpoly_worker", 3, A, B, stoi(d));
    2592        1596 :   init_modular_big(&S);
    2593        1596 :   H = gen_crt("ZXQ_minpoly", worker, &S, dB, bound, 0, NULL,
    2594             :                nxV_chinese_center, FpX_center_i);
    2595        1596 :   return gerepilecopy(av, H);
    2596             : }
    2597             : 
    2598             : /************************************************************************
    2599             :  *                                                                      *
    2600             :  *                   ZX_ZXY_resultant                                   *
    2601             :  *                                                                      *
    2602             :  ************************************************************************/
    2603             : 
    2604             : static GEN
    2605       14524 : ZX_ZXY_resultant_prime(GEN a, GEN b, ulong dp, ulong p,
    2606             :                        long degA, long degB, long dres, long sX)
    2607             : {
    2608       14524 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    2609       14524 :   GEN Hp = Flx_FlxY_resultant_polint(a, b, p, dres, sX);
    2610       14524 :   if (dropa && dropb)
    2611           0 :     Hp = zero_Flx(sX);
    2612             :   else {
    2613       14524 :     if (dropa)
    2614             :     { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2615           0 :       GEN c = gel(b,degB+2); /* lc(B) */
    2616           0 :       if (odd(degB)) c = Flx_neg(c, p);
    2617           0 :       if (!Flx_equal1(c)) {
    2618           0 :         c = Flx_powu(c, dropa, p);
    2619           0 :         if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
    2620             :       }
    2621             :     }
    2622       14524 :     else if (dropb)
    2623             :     { /* multiply by lc(A)^(deg B - deg b) */
    2624           0 :       ulong c = uel(a, degA+2); /* lc(A) */
    2625           0 :       c = Fl_powu(c, dropb, p);
    2626           0 :       if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
    2627             :     }
    2628             :   }
    2629       14524 :   if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2630       14524 :   return Hp;
    2631             : }
    2632             : 
    2633             : static GEN
    2634        5781 : ZX_ZXY_resultant_slice(GEN A, GEN B, GEN dB, long degA, long degB, long dres,
    2635             :                        GEN P, GEN *mod, long sX, long vY)
    2636             : {
    2637        5781 :   pari_sp av = avma;
    2638        5781 :   long i, n = lg(P)-1;
    2639             :   GEN H, T, D;
    2640        5781 :   if (n == 1)
    2641             :   {
    2642        2334 :     ulong p = uel(P,1);
    2643        2334 :     ulong dp = dB ? umodiu(dB, p): 1;
    2644        2334 :     GEN a = ZX_to_Flx(A, p), b = ZXX_to_FlxX(B, p, vY);
    2645        2334 :     GEN Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2646        2334 :     H = Flx_to_ZX(Hp);
    2647        2334 :     *mod = utoi(p);
    2648        2334 :     gerepileall(av, 2, &H, mod);
    2649        2334 :     return H;
    2650             :   }
    2651        3447 :   T = ZV_producttree(P);
    2652        3447 :   A = ZX_nv_mod_tree(A, P, T);
    2653        3447 :   B = ZXX_nv_mod_tree(B, P, T, vY);
    2654        3447 :   D = dB ? Z_ZV_mod_tree(dB, P, T): NULL;
    2655        3447 :   H = cgetg(n+1, t_VEC);
    2656       12480 :   for(i=1; i <= n; i++)
    2657             :   {
    2658        9033 :     ulong p = P[i];
    2659        9033 :     GEN a = gel(A,i), b = gel(B,i);
    2660        9033 :     ulong dp = D ? uel(D, i): 1;
    2661        9033 :     gel(H,i) = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2662             :   }
    2663        3447 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2664        3447 :   *mod = gmael(T, lg(T)-1, 1);
    2665        3447 :   gerepileall(av, 2, &H, mod);
    2666        3447 :   return H;
    2667             : }
    2668             : 
    2669             : GEN
    2670        5781 : ZX_ZXY_resultant_worker(GEN P, GEN A, GEN B, GEN dB, GEN v)
    2671             : {
    2672        5781 :   GEN V = cgetg(3, t_VEC);
    2673        5781 :   if (isintzero(dB)) dB = NULL;
    2674        5781 :   gel(V,1) = ZX_ZXY_resultant_slice(A, B, dB, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
    2675        5781 :   return V;
    2676             : }
    2677             : 
    2678             : GEN
    2679        4866 : ZX_ZXY_resultant(GEN A, GEN B)
    2680             : {
    2681        4866 :   pari_sp av = avma;
    2682             :   forprime_t S;
    2683             :   ulong bound;
    2684        4866 :   long v = fetch_var_higher();
    2685        4866 :   long degA = degpol(A), degB, dres = degA * degpol(B);
    2686        4866 :   long vX = varn(B), vY = varn(A); /* assume vY has lower priority */
    2687        4866 :   long sX = evalvarn(vX);
    2688             :   GEN worker, H, dB;
    2689        4866 :   B = Q_remove_denom(B, &dB);
    2690        4866 :   if (!dB) B = leafcopy(B);
    2691        4866 :   A = leafcopy(A); setvarn(A,v);
    2692        4866 :   B = swap_vars(B, vY); setvarn(B,v); degB = degpol(B);
    2693        4866 :   bound = ZX_ZXY_ResBound(A, B, dB);
    2694        4866 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    2695        9732 :   worker = snm_closure(is_entry("_ZX_ZXY_resultant_worker"),
    2696        4866 :                        mkvec4(A, B, dB? dB: gen_0,
    2697             :                               mkvecsmall5(degA, degB,dres, vY, sX)));
    2698        4866 :   init_modular_big(&S);
    2699        4866 :   H = gen_crt("ZX_ZXY_resultant_all", worker, &S, dB, bound, 0, NULL,
    2700             :                nxV_chinese_center, FpX_center_i);
    2701        4866 :   setvarn(H, vX); (void)delete_var();
    2702        4866 :   return gerepilecopy(av, H);
    2703             : }
    2704             : 
    2705             : static long
    2706        2583 : ZX_ZXY_rnfequation_lambda(GEN A, GEN B0, long lambda)
    2707             : {
    2708        2583 :   pari_sp av = avma;
    2709        2583 :   long degA = degpol(A), degB, dres = degA*degpol(B0);
    2710        2583 :   long v = fetch_var_higher();
    2711        2583 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    2712        2583 :   long sX = evalvarn(vX);
    2713             :   GEN dB, B, a, b, Hp;
    2714             :   forprime_t S;
    2715             : 
    2716        2583 :   B0 = Q_remove_denom(B0, &dB);
    2717        2583 :   if (!dB) B0 = leafcopy(B0);
    2718        2583 :   A = leafcopy(A);
    2719        2583 :   B = B0;
    2720        2583 :   setvarn(A,v);
    2721        3157 : INIT:
    2722        3157 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    2723        3157 :   B = swap_vars(B, vY); setvarn(B,v);
    2724             :   /* B0(lambda v + x, v) */
    2725        3157 :   if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    2726             : 
    2727        3157 :   degB = degpol(B);
    2728        3157 :   init_modular_big(&S);
    2729             :   while (1)
    2730           0 :   {
    2731        3157 :     ulong p = u_forprime_next(&S);
    2732        3157 :     ulong dp = dB ? umodiu(dB, p): 1;
    2733        3157 :     if (!dp) continue;
    2734        3157 :     a = ZX_to_Flx(A, p);
    2735        3157 :     b = ZXX_to_FlxX(B, p, v);
    2736        3157 :     Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2737        3157 :     if (degpol(Hp) != dres) continue;
    2738        3157 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2739        3157 :     if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); goto INIT; }
    2740        2583 :     if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    2741        2583 :     (void)delete_var(); return gc_long(av,lambda);
    2742             :   }
    2743             : }
    2744             : 
    2745             : GEN
    2746        3256 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
    2747             : {
    2748        3256 :   if (lambda)
    2749             :   {
    2750        2583 :     *lambda = ZX_ZXY_rnfequation_lambda(A, B, *lambda);
    2751        2583 :     B = RgX_translate(B, monomial(stoi(*lambda), 1, varn(A)));
    2752             :   }
    2753        3256 :   return ZX_ZXY_resultant(A,B);
    2754             : }
    2755             : 
    2756             : static GEN
    2757         668 : ZX_direct_compositum_slice(GEN A, GEN B, GEN P, GEN *mod)
    2758             : {
    2759         668 :   pari_sp av = avma;
    2760         668 :   long i, n = lg(P)-1;
    2761             :   GEN H, T;
    2762         668 :   if (n == 1)
    2763             :   {
    2764         499 :     ulong p = uel(P,1);
    2765         499 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
    2766         499 :     GEN Hp = Flx_direct_compositum(a, b, p);
    2767         499 :     H = gerepileupto(av, Flx_to_ZX(Hp));
    2768         499 :     *mod = utoi(p);
    2769         499 :     return H;
    2770             :   }
    2771         169 :   T = ZV_producttree(P);
    2772         169 :   A = ZX_nv_mod_tree(A, P, T);
    2773         169 :   B = ZX_nv_mod_tree(B, P, T);
    2774         169 :   H = cgetg(n+1, t_VEC);
    2775         843 :   for(i=1; i <= n; i++)
    2776             :   {
    2777         674 :     ulong p = P[i];
    2778         674 :     GEN a = gel(A,i), b = gel(B,i);
    2779         674 :     gel(H,i) = Flx_direct_compositum(a, b, p);
    2780             :   }
    2781         169 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2782         169 :   *mod = gmael(T, lg(T)-1, 1);
    2783         169 :   gerepileall(av, 2, &H, mod);
    2784         169 :   return H;
    2785             : }
    2786             : 
    2787             : GEN
    2788         668 : ZX_direct_compositum_worker(GEN P, GEN A, GEN B)
    2789             : {
    2790         668 :   GEN V = cgetg(3, t_VEC);
    2791         668 :   gel(V,1) = ZX_direct_compositum_slice(A, B, P, &gel(V,2));
    2792         668 :   return V;
    2793             : }
    2794             : 
    2795             : /* Assume A,B irreducible (in particular squarefree) and define linearly
    2796             :  * disjoint extensions: no factorisation needed */
    2797             : static GEN
    2798         567 : ZX_direct_compositum(GEN A, GEN B, GEN lead)
    2799             : {
    2800         567 :   pari_sp av = avma;
    2801             :   forprime_t S;
    2802             :   ulong bound;
    2803             :   GEN H, worker, mod;
    2804         567 :   bound = ZX_ZXY_ResBound(A, poleval(B,deg1pol(gen_1,pol_x(1),0)), NULL);
    2805         567 :   worker = snm_closure(is_entry("_ZX_direct_compositum_worker"), mkvec2(A,B));
    2806         567 :   init_modular_big(&S);
    2807         567 :   H = gen_crt("ZX_direct_compositum", worker, &S, lead, bound, 0, &mod,
    2808             :               nxV_chinese_center, FpX_center);
    2809         567 :   return gerepileupto(av, H);
    2810             : }
    2811             : 
    2812             : static long
    2813         182 : ZX_compositum_lambda(GEN A, GEN B, GEN lead, long lambda)
    2814             : {
    2815         182 :   pari_sp av = avma;
    2816             :   forprime_t S;
    2817             :   ulong p;
    2818         182 :   init_modular_big(&S);
    2819         182 :   p = u_forprime_next(&S);
    2820             :   while (1)
    2821          56 :   {
    2822             :     GEN Hp, a;
    2823         238 :     if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    2824         238 :     if (lead && dvdiu(lead,p)) { p = u_forprime_next(&S); continue; }
    2825         231 :     a = ZX_to_Flx(ZX_rescale(A, stoi(-lambda)), p);
    2826         231 :     Hp = Flx_direct_compositum(a, ZX_to_Flx(B, p), p);
    2827         231 :     if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); continue; }
    2828         182 :     if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    2829         182 :     return gc_long(av, lambda);
    2830             :   }
    2831             : }
    2832             : 
    2833             : GEN
    2834         567 : ZX_compositum(GEN A, GEN B, long *lambda)
    2835             : {
    2836         567 :   GEN lead  = mulii(leading_coeff(A),leading_coeff(B));
    2837         567 :   if (lambda)
    2838             :   {
    2839         182 :     *lambda = ZX_compositum_lambda(A, B, lead, *lambda);
    2840         182 :     A = ZX_rescale(A, stoi(-*lambda));
    2841             :   }
    2842         567 :   return ZX_direct_compositum(A, B, lead);
    2843             : }
    2844             : 
    2845             : GEN
    2846         385 : ZX_compositum_disjoint(GEN A, GEN B)
    2847         385 : { return ZX_compositum(A, B, NULL); }
    2848             : 
    2849             : static GEN
    2850         243 : ZXQX_direct_compositum_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
    2851             : {
    2852         243 :   pari_sp av = avma;
    2853         243 :   long i, n = lg(P)-1, dC = degpol(C), v = varn(C);
    2854             :   GEN H, T;
    2855         243 :   if (n == 1)
    2856             :   {
    2857         170 :     ulong p = uel(P,1);
    2858         170 :     GEN a = ZXX_to_FlxX(A, p, v), b = ZXX_to_FlxX(B, p, v);
    2859         170 :     GEN c = ZX_to_Flx(C, p);
    2860         170 :     GEN Hp = FlxX_to_Flm(FlxqX_direct_compositum(a, b, c, p), dC);
    2861         170 :     H = gerepileupto(av, Flm_to_ZM(Hp));
    2862         170 :     *mod = utoi(p);
    2863         170 :     return H;
    2864             :   }
    2865          73 :   T = ZV_producttree(P);
    2866          73 :   A = ZXX_nv_mod_tree(A, P, T, v);
    2867          73 :   B = ZXX_nv_mod_tree(B, P, T, v);
    2868          73 :   C = ZX_nv_mod_tree(C, P, T);
    2869          73 :   H = cgetg(n+1, t_VEC);
    2870         279 :   for(i=1; i <= n; i++)
    2871             :   {
    2872         206 :     ulong p = P[i];
    2873         206 :     GEN a = gel(A,i), b = gel(B,i), c = gel(C,i);
    2874         206 :     gel(H,i) = FlxX_to_Flm(FlxqX_direct_compositum(a, b, c, p), dC);
    2875             :   }
    2876          73 :   H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
    2877          73 :   *mod = gmael(T, lg(T)-1, 1);
    2878          73 :   gerepileall(av, 2, &H, mod);
    2879          73 :   return H;
    2880             : }
    2881             : 
    2882             : GEN
    2883         243 : ZXQX_direct_compositum_worker(GEN P, GEN A, GEN B, GEN C)
    2884             : {
    2885         243 :   GEN V = cgetg(3, t_VEC);
    2886         243 :   gel(V,1) = ZXQX_direct_compositum_slice(A, B, C, P, &gel(V,2));
    2887         243 :   return V;
    2888             : }
    2889             : 
    2890             : static GEN
    2891         224 : ZXQX_direct_compositum(GEN A, GEN B, GEN T, ulong bound)
    2892             : {
    2893         224 :   pari_sp av = avma;
    2894             :   forprime_t S;
    2895             :   GEN H, worker, mod;
    2896         224 :   GEN lead = mulii(Q_content(leading_coeff(A)), Q_content(leading_coeff(B)));
    2897         224 :   worker = snm_closure(is_entry("_ZXQX_direct_compositum_worker")
    2898             :                       , mkvec3(A,B,T));
    2899         224 :   init_modular_big(&S);
    2900         224 :   H = gen_crt("ZXQX_direct_compositum", worker, &S, lead, bound, 0, &mod,
    2901             :               nmV_chinese_center, FpM_center);
    2902         224 :   if (DEBUGLEVEL > 4)
    2903           0 :     err_printf("nfcompositum: a priori bound: %lu, a posteriori: %lu\n",
    2904             :                bound, expi(gsupnorm(H, DEFAULTPREC)));
    2905         224 :   return gerepilecopy(av, RgM_to_RgXX(H, varn(A), varn(T)));
    2906             : }
    2907             : 
    2908             : static long
    2909         224 : ZXQX_direct_compositum_bound(GEN nf, GEN A, GEN B)
    2910             : {
    2911         224 :   pari_sp av = avma;
    2912         224 :   GEN r, M = L2_bound(nf, NULL, &r);
    2913         224 :   long v = nf_get_varn(nf), i, l = lg(r);
    2914         224 :   GEN a = cgetg(l, t_COL);
    2915         791 :   for (i = 1; i < l; i++)
    2916        1134 :     gel(a, i) =  RgX_RgXY_ResBound(gsubst(A, v, gel(r,i)),
    2917         567 :                  poleval(gsubst(B, v, gel(r,i)),
    2918             :                          deg1pol(gen_1, pol_x(1), 0)), DEFAULTPREC);
    2919         224 :   return gc_long(av, (long) dbllog2(gmul(M,RgC_fpnorml2(a, DEFAULTPREC))));
    2920             : }
    2921             : 
    2922             : GEN
    2923         224 : nf_direct_compositum(GEN nf, GEN A, GEN B)
    2924             : {
    2925         224 :   ulong bnd = ZXQX_direct_compositum_bound(nf, A, B);
    2926         224 :   return ZXQX_direct_compositum(A, B, nf_get_pol(nf), bnd);
    2927             : }
    2928             : 
    2929             : /************************************************************************
    2930             :  *                                                                      *
    2931             :  *                   IRREDUCIBLE POLYNOMIAL / Fp                        *
    2932             :  *                                                                      *
    2933             :  ************************************************************************/
    2934             : 
    2935             : /* irreducible (unitary) polynomial of degree n over Fp */
    2936             : GEN
    2937           0 : ffinit_rand(GEN p,long n)
    2938             : {
    2939           0 :   for(;;) {
    2940           0 :     pari_sp av = avma;
    2941           0 :     GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
    2942           0 :     if (FpX_is_irred(pol, p)) return pol;
    2943           0 :     set_avma(av);
    2944             :   }
    2945             : }
    2946             : 
    2947             : /* return an extension of degree 2^l of F_2, assume l > 0
    2948             :  * Not stack clean. */
    2949             : static GEN
    2950         627 : ffinit_Artin_Schreier_2(long l)
    2951             : {
    2952             :   GEN Q, T, S;
    2953             :   long i, v;
    2954             : 
    2955         627 :   if (l == 1) return mkvecsmall4(0,1,1,1); /*x^2 + x + 1*/
    2956         592 :   v = fetch_var_higher();
    2957         592 :   S = mkvecsmall5(0, 0, 0, 1, 1); /* y(y^2 + y) */
    2958         591 :   Q = mkpoln(3, pol1_Flx(0), pol1_Flx(0), S); /* x^2 + x + y(y^2+y) */
    2959         591 :   setvarn(Q, v);
    2960             : 
    2961             :   /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
    2962         591 :   T = mkvecsmalln(6,evalvarn(v),1UL,1UL,0UL,0UL,1UL);
    2963             :   /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
    2964             :    * ==> x^2 + x + a(y) b irred. over K for any root b of Q
    2965             :    * ==> x^2 + x + (b^2+b)b */
    2966        3395 :   for (i=2; i<l; i++) T = Flx_FlxY_resultant(T, Q, 2); /* minpoly(b) / F2*/
    2967         592 :   (void)delete_var(); T[1] = 0; return T;
    2968             : }
    2969             : 
    2970             : /* return an extension of degree p^l of F_p, assume l > 0
    2971             :  * Not stack clean. */
    2972             : GEN
    2973         823 : ffinit_Artin_Schreier(ulong p, long l)
    2974             : {
    2975             :   long i, v;
    2976             :   GEN Q, R, S, T, xp;
    2977         823 :   if (p==2) return ffinit_Artin_Schreier_2(l);
    2978         196 :   xp = polxn_Flx(p,0); /* x^p */
    2979         196 :   T = Flx_sub(xp, mkvecsmall3(0,1,1),p); /* x^p - x - 1 */
    2980         196 :   if (l == 1) return T;
    2981             : 
    2982           7 :   v = evalvarn(fetch_var_higher());
    2983           7 :   xp[1] = v;
    2984           7 :   R = Flx_sub(polxn_Flx(2*p-1,0), polxn_Flx(p,0),p);
    2985           7 :   S = Flx_sub(xp, polx_Flx(0), p);
    2986           7 :   Q = FlxX_Flx_sub(Flx_to_FlxX(S, v), R, p); /* x^p - x - (y^(2p-1)-y^p) */
    2987          14 :   for (i = 2; i <= l; ++i) T = Flx_FlxY_resultant(T, Q, p);
    2988           7 :   (void)delete_var(); T[1] = 0; return T;
    2989             : }
    2990             : 
    2991             : static long
    2992       20853 : flinit_check(ulong p, long n, long l)
    2993             : {
    2994             :   ulong q;
    2995       20853 :   if (!uisprime(n)) return 0;
    2996       12848 :   q = p % n; if (!q) return 0;
    2997       10769 :   return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    2998             : }
    2999             : 
    3000             : static GEN
    3001        5024 : flinit(ulong p, long l)
    3002             : {
    3003        5024 :   ulong n = 1+l;
    3004       13569 :   while (!flinit_check(p,n,l)) n += l;
    3005        5024 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    3006        5024 :   return ZX_to_Flx(polsubcyclo(n,l,0), p);
    3007             : }
    3008             : 
    3009             : static GEN
    3010        5219 : ffinit_fact_Flx(ulong p, long n)
    3011             : {
    3012        5219 :   GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
    3013        5219 :   long i, l = lg(Fm);
    3014        5219 :   P = cgetg(l, t_VEC);
    3015       11066 :   for (i = 1; i < l; ++i)
    3016        5847 :     gel(P,i) = p==uel(Fp,i) ?
    3017         823 :                  ffinit_Artin_Schreier(uel(Fp,i), Fe[i])
    3018        5847 :                : flinit(p, uel(Fm,i));
    3019        5219 :   return FlxV_direct_compositum(P, p);
    3020             : }
    3021             : 
    3022             : static GEN
    3023        7284 : init_Flxq_i(ulong p, long n, long sv)
    3024             : {
    3025             :   GEN P;
    3026        7284 :   if (n == 1) return polx_Flx(sv);
    3027        7284 :   if (flinit_check(p, n+1, n))
    3028             :   {
    3029        2065 :     P = const_vecsmall(n+2,1);
    3030        2065 :     P[1] = sv; return P;
    3031             :   }
    3032        5219 :   P = ffinit_fact_Flx(p,n);
    3033        5219 :   P[1] = sv; return P;
    3034             : }
    3035             : 
    3036             : GEN
    3037           0 : init_Flxq(ulong p, long n, long v)
    3038             : {
    3039           0 :   pari_sp av = avma;
    3040           0 :   return gerepileupto(av, init_Flxq_i(p, n, v));
    3041             : }
    3042             : 
    3043             : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
    3044             : static long
    3045        2188 : fpinit_check(GEN p, long n, long l)
    3046             : {
    3047             :   ulong q;
    3048        2188 :   if (!uisprime(n)) return 0;
    3049        1537 :   q = umodiu(p,n); if (!q) return 0;
    3050        1537 :   return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    3051             : }
    3052             : 
    3053             : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
    3054             :  * Return an irreducible polynomial of degree l over F_p.
    3055             :  * Variant of Adleman and Lenstra "Finding irreducible polynomials over
    3056             :  * finite fields", ACM, 1986 (5) 350--355.
    3057             :  * Not stack clean */
    3058             : static GEN
    3059         527 : fpinit(GEN p, long l)
    3060             : {
    3061         527 :   ulong n = 1+l;
    3062        1600 :   while (!fpinit_check(p,n,l)) n += l;
    3063         527 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    3064         527 :   return FpX_red(polsubcyclo(n,l,0),p);
    3065             : }
    3066             : 
    3067             : static GEN
    3068         497 : ffinit_fact(GEN p, long n)
    3069             : {
    3070         497 :   GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
    3071         497 :   long i, l = lg(Fm);
    3072         497 :   P = cgetg(l, t_VEC);
    3073        1024 :   for (i = 1; i < l; ++i)
    3074        1054 :     gel(P,i) = absequaliu(p, Fp[i]) ?
    3075           0 :                  Flx_to_ZX(ffinit_Artin_Schreier(Fp[i], Fe[i]))
    3076         527 :                : fpinit(p, Fm[i]);
    3077         497 :   return FpXV_direct_compositum(P, p);
    3078             : }
    3079             : 
    3080             : static GEN
    3081        8124 : init_Fq_i(GEN p, long n, long v)
    3082             : {
    3083             :   GEN P;
    3084        8124 :   if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
    3085        8124 :   if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
    3086        8124 :   if (signe(p) <= 0) pari_err_PRIME("ffinit",p);
    3087        8124 :   if (v < 0) v = 0;
    3088        8124 :   if (n == 1) return pol_x(v);
    3089        7872 :   if (lgefint(p) == 3)
    3090        7284 :     return Flx_to_ZX(init_Flxq_i(p[2], n, evalvarn(v)));
    3091         588 :   if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
    3092         497 :   P = ffinit_fact(p,n);
    3093         497 :   setvarn(P, v); return P;
    3094             : }
    3095             : GEN
    3096        7669 : init_Fq(GEN p, long n, long v)
    3097             : {
    3098        7669 :   pari_sp av = avma;
    3099        7669 :   return gerepileupto(av, init_Fq_i(p, n, v));
    3100             : }
    3101             : GEN
    3102         455 : ffinit(GEN p, long n, long v)
    3103             : {
    3104         455 :   pari_sp av = avma;
    3105         455 :   return gerepileupto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
    3106             : }
    3107             : 
    3108             : GEN
    3109        3178 : ffnbirred(GEN p, long n)
    3110             : {
    3111        3178 :   pari_sp av = avma;
    3112        3178 :   GEN s = powiu(p,n), F = factoru(n), D = divisorsu_moebius(gel(F, 1));
    3113        3178 :   long j, l = lg(D);
    3114        6797 :   for (j = 2; j < l; j++) /* skip d = 1 */
    3115             :   {
    3116        3619 :     long md = D[j]; /* mu(d) * d, d squarefree */
    3117        3619 :     GEN pd = powiu(p, n / labs(md)); /* p^{n/d} */
    3118        3619 :     s = md > 0? addii(s, pd): subii(s,pd);
    3119             :   }
    3120        3178 :   return gerepileuptoint(av, diviuexact(s, n));
    3121             : }
    3122             : 
    3123             : GEN
    3124         420 : ffsumnbirred(GEN p, long n)
    3125             : {
    3126         420 :   pari_sp av = avma, av2;
    3127         420 :   GEN q, t = p, v = vecfactoru(1, n);
    3128             :   long i;
    3129         420 :   q = cgetg(n+1,t_VEC); gel(q,1) = p;
    3130        1568 :   for (i=2; i<=n; i++) gel(q,i) = mulii(gel(q,i-1), p);
    3131         420 :   av2 = avma;
    3132        1568 :   for (i=2; i<=n; i++)
    3133             :   {
    3134        1148 :     GEN s = gel(q,i), F = gel(v,i), D = divisorsu_moebius(gel(F,1));
    3135        1148 :     long j, l = lg(D);
    3136        2534 :     for (j = 2; j < l; j++) /* skip 1 */
    3137             :     {
    3138        1386 :       long md = D[j];
    3139        1386 :       GEN pd = gel(q, i / labs(md)); /* p^{i/d} */
    3140        1386 :       s = md > 0? addii(s, pd): subii(s, pd);
    3141             :     }
    3142        1148 :     t = gerepileuptoint(av2, addii(t, diviuexact(s, i)));
    3143             :   }
    3144         420 :   return gerepileuptoint(av, t);
    3145             : }
    3146             : 
    3147             : GEN
    3148         140 : ffnbirred0(GEN p, long n, long flag)
    3149             : {
    3150         140 :   if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
    3151         140 :   if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
    3152         140 :   switch(flag)
    3153             :   {
    3154          70 :     case 0: return ffnbirred(p, n);
    3155          70 :     case 1: return ffsumnbirred(p, n);
    3156             :   }
    3157           0 :   pari_err_FLAG("ffnbirred");
    3158             :   return NULL; /* LCOV_EXCL_LINE */
    3159             : }
    3160             : 
    3161             : static void
    3162        2254 : checkmap(GEN m, const char *s)
    3163             : {
    3164        2254 :   if (typ(m)!=t_VEC || lg(m)!=3 || typ(gel(m,1))!=t_FFELT)
    3165           0 :     pari_err_TYPE(s,m);
    3166        2254 : }
    3167             : 
    3168             : GEN
    3169         182 : ffembed(GEN a, GEN b)
    3170             : {
    3171         182 :   pari_sp av = avma;
    3172         182 :   GEN p, Ta, Tb, g, r = NULL;
    3173         182 :   if (typ(a)!=t_FFELT) pari_err_TYPE("ffembed",a);
    3174         182 :   if (typ(b)!=t_FFELT) pari_err_TYPE("ffembed",b);
    3175         182 :   p = FF_p_i(a); g = FF_gen(a);
    3176         182 :   if (!equalii(p, FF_p_i(b))) pari_err_MODULUS("ffembed",a,b);
    3177         182 :   Ta = FF_mod(a);
    3178         182 :   Tb = FF_mod(b);
    3179         182 :   if (degpol(Tb)%degpol(Ta)!=0)
    3180           7 :     pari_err_DOMAIN("ffembed",GENtostr_raw(a),"is not a subfield of",b,a);
    3181         175 :   r = gel(FFX_roots(Ta, b), 1);
    3182         175 :   return gerepilecopy(av, mkvec2(g,r));
    3183             : }
    3184             : 
    3185             : GEN
    3186          91 : ffextend(GEN a, GEN P, long v)
    3187             : {
    3188          91 :   pari_sp av = avma;
    3189             :   long n;
    3190             :   GEN p, T, R, g, m;
    3191          91 :   if (typ(a)!=t_FFELT) pari_err_TYPE("ffextend",a);
    3192          91 :   T = a; p = FF_p_i(a);
    3193          91 :   if (typ(P)!=t_POL || !RgX_is_FpXQX(P,&T,&p)) pari_err_TYPE("ffextend", P);
    3194          49 :   if (!FF_samefield(a, T)) pari_err_MODULUS("ffextend",a,T);
    3195          49 :   if (v < 0) v = varn(P);
    3196          49 :   n = FF_f(T) * degpol(P); R = ffinit(p, n, v); g = ffgen(R, v);
    3197          49 :   m = ffembed(a, g);
    3198          49 :   R = FFX_roots(ffmap(m, P),g);
    3199          49 :   return gerepilecopy(av, mkvec2(gel(R,1), m));
    3200             : }
    3201             : 
    3202             : GEN
    3203          42 : fffrobenius(GEN a, long n)
    3204             : {
    3205          42 :   if (typ(a)!=t_FFELT) pari_err_TYPE("fffrobenius",a);
    3206          42 :   retmkvec2(FF_gen(a), FF_Frobenius(a, n));
    3207             : }
    3208             : 
    3209             : GEN
    3210         133 : ffinvmap(GEN m)
    3211             : {
    3212         133 :   pari_sp av = avma;
    3213             :   long i, l;
    3214         133 :   GEN T, F, a, g, r, f = NULL;
    3215         133 :   checkmap(m, "ffinvmap");
    3216         133 :   a = gel(m,1); r = gel(m,2);
    3217         133 :   if (typ(r) != t_FFELT)
    3218           7 :    pari_err_TYPE("ffinvmap", m);
    3219         126 :   g = FF_gen(a);
    3220         126 :   T = FF_mod(r);
    3221         126 :   F = gel(FFX_factor(T, a), 1);
    3222         126 :   l = lg(F);
    3223         490 :   for(i=1; i<l; i++)
    3224             :   {
    3225         490 :     GEN s = FFX_rem(FF_to_FpXQ_i(r), gel(F, i), a);
    3226         490 :     if (degpol(s)==0 && gequal(constant_coeff(s),g)) { f = gel(F, i); break; }
    3227             :   }
    3228         126 :   if (f==NULL) pari_err_TYPE("ffinvmap", m);
    3229         126 :   if (degpol(f)==1) f = FF_neg_i(gel(f,2));
    3230         126 :   return gerepilecopy(av, mkvec2(FF_gen(r),f));
    3231             : }
    3232             : 
    3233             : static GEN
    3234        1260 : ffpartmapimage(const char *s, GEN r)
    3235             : {
    3236        1260 :    GEN a = NULL, p = NULL;
    3237        1260 :    if (typ(r)==t_POL && degpol(r) >= 1
    3238        1260 :       && RgX_is_FpXQX(r,&a,&p) && a && typ(a)==t_FFELT) return a;
    3239           0 :    pari_err_TYPE(s, r);
    3240             :    return NULL; /* LCOV_EXCL_LINE */
    3241             : }
    3242             : 
    3243             : static GEN
    3244        2702 : ffeltmap_i(GEN m, GEN x)
    3245             : {
    3246        2702 :    GEN r = gel(m,2);
    3247        2702 :    if (!FF_samefield(x, gel(m,1)))
    3248          84 :      pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
    3249        2618 :    if (typ(r)==t_FFELT)
    3250        1652 :      return FF_map(r, x);
    3251             :    else
    3252         966 :      return FFX_preimage(x, r, ffpartmapimage("ffmap", r));
    3253             : }
    3254             : 
    3255             : static GEN
    3256        4452 : ffmap_i(GEN m, GEN x)
    3257             : {
    3258             :   GEN y;
    3259        4452 :   long i, lx, tx = typ(x);
    3260        4452 :   switch(tx)
    3261             :   {
    3262        2534 :     case t_FFELT:
    3263        2534 :       return ffeltmap_i(m, x);
    3264        1267 :     case t_POL: case t_RFRAC: case t_SER:
    3265             :     case t_VEC: case t_COL: case t_MAT:
    3266        1267 :       y = cgetg_copy(x, &lx);
    3267        1988 :       for (i=1; i<lontyp[tx]; i++) y[i] = x[1];
    3268        4564 :       for (i=lontyp[tx]; i<lx; i++)
    3269             :       {
    3270        3339 :         GEN yi = ffmap_i(m, gel(x,i));
    3271        3297 :         if (!yi) return NULL;
    3272        3297 :         gel(y,i) = yi;
    3273             :       }
    3274        1225 :       return y;
    3275             :   }
    3276         651 :   return gcopy(x);
    3277             : }
    3278             : 
    3279             : GEN
    3280        1029 : ffmap(GEN m, GEN x)
    3281             : {
    3282        1029 :   pari_sp ltop = avma;
    3283             :   GEN y;
    3284        1029 :   checkmap(m, "ffmap");
    3285        1029 :   y = ffmap_i(m, x);
    3286        1029 :   if (y) return y;
    3287          42 :   set_avma(ltop); return cgetg(1,t_VEC);
    3288             : }
    3289             : 
    3290             : static GEN
    3291         252 : ffeltmaprel_i(GEN m, GEN x)
    3292             : {
    3293         252 :    GEN g = gel(m,1), r = gel(m,2);
    3294         252 :    if (!FF_samefield(x, g))
    3295           0 :      pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
    3296         252 :    if (typ(r)==t_FFELT)
    3297          84 :      retmkpolmod(FF_map(r, x), pol_x(FF_var(g)));
    3298             :    else
    3299         168 :      retmkpolmod(FFX_preimagerel(x, r, ffpartmapimage("ffmap", r)), gcopy(r));
    3300             : }
    3301             : 
    3302             : static GEN
    3303         252 : ffmaprel_i(GEN m, GEN x)
    3304             : {
    3305             :   GEN y;
    3306         252 :   long i, lx, tx = typ(x);
    3307         252 :   switch(tx)
    3308             :   {
    3309         252 :     case t_FFELT:
    3310         252 :       return ffeltmaprel_i(m, x);
    3311           0 :     case t_POL: case t_RFRAC: case t_SER:
    3312             :     case t_VEC: case t_COL: case t_MAT:
    3313           0 :       y = cgetg_copy(x, &lx);
    3314           0 :       for (i=1; i<lontyp[tx]; i++) y[i] = x[1];
    3315           0 :       for (i=lontyp[tx]; i<lx; i++)
    3316           0 :         gel(y,i) = ffmaprel_i(m, gel(x,i));
    3317           0 :       return y;
    3318             :   }
    3319           0 :   return gcopy(x);
    3320             : }
    3321             : 
    3322             : GEN
    3323         252 : ffmaprel(GEN m, GEN x)
    3324             : {
    3325         252 :   checkmap(m, "ffmaprel");
    3326         252 :   return ffmaprel_i(m, x);
    3327             : }
    3328             : 
    3329             : static void
    3330          84 : err_compo(GEN m, GEN n)
    3331          84 : { pari_err_DOMAIN("ffcompomap","m","domain does not contain codomain of",n,m); }
    3332             : 
    3333             : GEN
    3334         420 : ffcompomap(GEN m, GEN n)
    3335             : {
    3336         420 :   pari_sp av = avma;
    3337         420 :   GEN g = gel(n,1), r, m2, n2;
    3338         420 :   checkmap(m, "ffcompomap");
    3339         420 :   checkmap(n, "ffcompomap");
    3340         420 :   m2 = gel(m,2); n2 = gel(n,2);
    3341         420 :   switch((typ(m2)==t_POL)|((typ(n2)==t_POL)<<1))
    3342             :   {
    3343          84 :     case 0:
    3344          84 :       if (!FF_samefield(gel(m,1),n2)) err_compo(m,n);
    3345          42 :       r = FF_map(gel(m,2), n2);
    3346          42 :       break;
    3347          84 :     case 2:
    3348          84 :       r = ffmap_i(m, n2);
    3349          42 :       if (lg(r) == 1) err_compo(m,n);
    3350          42 :       break;
    3351         168 :     case 1:
    3352         168 :       r = ffeltmap_i(m, n2);
    3353         126 :       if (!r)
    3354             :       {
    3355             :         GEN a, A, R, M;
    3356             :         long dm, dn;
    3357          42 :         a = ffpartmapimage("ffcompomap",m2);
    3358          42 :         A = FF_to_FpXQ_i(FF_neg(n2));
    3359          42 :         setvarn(A, 1);
    3360          42 :         R = deg1pol(gen_1, A, 0);
    3361          42 :         setvarn(R, 0);
    3362          42 :         M = gcopy(m2);
    3363          42 :         setvarn(M, 1);
    3364          42 :         r = polresultant0(R, M, 1, 0);
    3365          42 :         dm = FF_f(gel(m,1)); dn = FF_f(gel(n,1));
    3366          42 :         if (dm % dn || !FFX_ispower(r, dm/dn, a, &r)) err_compo(m,n);
    3367          42 :         setvarn(r, varn(FF_mod(g)));
    3368             :       }
    3369         126 :       break;
    3370          84 :     case 3:
    3371             :     {
    3372             :       GEN M, R, T, p, a;
    3373          84 :       a = ffpartmapimage("ffcompomap",n2);
    3374          84 :       if (!FF_samefield(a, gel(m,1))) err_compo(m,n);
    3375          42 :       p = FF_p_i(gel(n,1));
    3376          42 :       T = FF_mod(gel(n,1));
    3377          42 :       setvarn(T, 1);
    3378          42 :       R = RgX_to_FpXQX(n2,T,p);
    3379          42 :       setvarn(R, 0);
    3380          42 :       M = gcopy(m2);
    3381          42 :       setvarn(M, 1);
    3382          42 :       r = polresultant0(R, M, 1, 0);
    3383          42 :       setvarn(r, varn(n2));
    3384             :     }
    3385             :   }
    3386         252 :   return gerepilecopy(av, mkvec2(g,r));
    3387             : }

Generated by: LCOV version 1.13