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.12.0 lcov report (development 23036-b751c0af5) Lines: 1305 1492 87.5 %
Date: 2018-09-26 05:46:06 Functions: 141 157 89.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000-2005  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /***********************************************************************/
      15             : /**                                                                   **/
      16             : /**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
      17             : /**                         (third part)                              **/
      18             : /**                                                                   **/
      19             : /***********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /************************************************************************
      24             :  **                                                                    **
      25             :  **                      Ring membership                               **
      26             :  **                                                                    **
      27             :  ************************************************************************/
      28             : struct charact {
      29             :   GEN q;
      30             :   int isprime;
      31             : };
      32             : static void
      33         637 : char_update_prime(struct charact *S, GEN p)
      34             : {
      35         637 :   if (!S->isprime) { S->isprime = 1; S->q = p; }
      36         637 :   if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
      37         630 : }
      38             : static void
      39        2219 : char_update_int(struct charact *S, GEN n)
      40             : {
      41        2219 :   if (S->isprime)
      42             :   {
      43           7 :     if (dvdii(n, S->q)) return;
      44           7 :     pari_err_MODULUS("characteristic", S->q, n);
      45             :   }
      46        2212 :   S->q = gcdii(S->q, n);
      47             : }
      48             : static void
      49      582967 : charact(struct charact *S, GEN x)
      50             : {
      51      582967 :   const long tx = typ(x);
      52             :   long i, l;
      53      582967 :   switch(tx)
      54             :   {
      55        1246 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      56         546 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      57             :     case t_COMPLEX: case t_QUAD:
      58             :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      59             :     case t_VEC: case t_COL: case t_MAT:
      60       12145 :       l = lg(x);
      61       12145 :       for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
      62       12131 :       break;
      63             :     case t_LIST:
      64           7 :       x = list_data(x);
      65           7 :       if (x) charact(S, x);
      66           7 :       break;
      67             :   }
      68      582939 : }
      69             : static void
      70       33215 : charact_res(struct charact *S, GEN x)
      71             : {
      72       33215 :   const long tx = typ(x);
      73             :   long i, l;
      74       33215 :   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             :     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       10332 :       l = lg(x);
      83       10332 :       for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
      84       10332 :       break;
      85             :     case t_LIST:
      86           0 :       x = list_data(x);
      87           0 :       if (x) charact_res(S, x);
      88           0 :       break;
      89             :   }
      90       33215 : }
      91             : GEN
      92        9744 : characteristic(GEN x)
      93             : {
      94             :   struct charact S;
      95        9744 :   S.q = gen_0; S.isprime = 0;
      96        9744 :   charact(&S, x); return S.q;
      97             : }
      98             : GEN
      99        2485 : residual_characteristic(GEN x)
     100             : {
     101             :   struct charact S;
     102        2485 :   S.q = gen_0; S.isprime = 0;
     103        2485 :   charact_res(&S, x); return S.q;
     104             : }
     105             : 
     106             : int
     107    57513268 : Rg_is_Fp(GEN x, GEN *pp)
     108             : {
     109             :   GEN mod;
     110    57513268 :   switch(typ(x))
     111             :   {
     112             :   case t_INTMOD:
     113     4217619 :     mod = gel(x,1);
     114     4217619 :     if (!*pp) *pp = mod;
     115     3982629 :     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     4217619 :     return 1;
     121             :   case t_INT:
     122    49672281 :     return 1;
     123     3623368 :   default: return 0;
     124             :   }
     125             : }
     126             : 
     127             : int
     128    19547297 : RgX_is_FpX(GEN x, GEN *pp)
     129             : {
     130    19547297 :   long i, lx = lg(x);
     131    73411682 :   for (i=2; i<lx; i++)
     132    57487753 :     if (!Rg_is_Fp(gel(x, i), pp))
     133     3623368 :       return 0;
     134    15923929 :   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       57288 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
     157             : {
     158             :   GEN pol, mod, p;
     159       57288 :   switch(typ(x))
     160             :   {
     161             :   case t_INTMOD:
     162       25508 :     return Rg_is_Fp(x, pp);
     163             :   case t_INT:
     164        6517 :     return 1;
     165             :   case t_POL:
     166          21 :     return RgX_is_FpX(x, pp);
     167             :   case t_FFELT:
     168       20615 :     mod = x; p = FF_p_i(x);
     169       20615 :     if (!*pp) *pp = p;
     170       20615 :     if (!*pT) *pT = mod;
     171       19257 :     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       20573 :     return 1;
     177             :   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        2849 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
     199             : {
     200        2849 :   long i, lx = lg(x);
     201       59619 :   for (i = 2; i < lx; i++)
     202       56812 :     if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
     203        2807 :   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    31658062 : Rg_to_Fp(GEN x, GEN p)
     216             : {
     217    31658062 :   if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
     218     2558565 :   switch(typ(x))
     219             :   {
     220      179948 :     case t_INT: return modii(x, p);
     221             :     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             :     case t_INTMOD: {
     229     2378496 :       GEN q = gel(x,1), a = gel(x,2);
     230     2378496 :       if (equalii(q, p)) return icopy(a);
     231          14 :       if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
     232           0 :       return remii(a, p);
     233             :     }
     234           0 :     default: pari_err_TYPE("Rg_to_Fp",x);
     235             :       return NULL; /* LCOV_EXCL_LINE */
     236             :   }
     237             : }
     238             : /* If x is a POLMOD, assume modulus is a multiple of T. */
     239             : GEN
     240     1256544 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
     241             : {
     242     1256544 :   long ta, tx = typ(x), v = get_FpX_var(T);
     243             :   GEN a, b;
     244     1256544 :   if (is_const_t(tx))
     245             :   {
     246       55043 :     if (tx == t_FFELT)
     247             :     {
     248       17085 :       GEN z = FF_to_FpXQ(x);
     249       17085 :       setvarn(z, v);
     250       17085 :       return z;
     251             :     }
     252       37958 :     return scalar_ZX(degpol(T)? Rg_to_Fp(x, p): gen_0, v);
     253             :   }
     254     1201501 :   switch(tx)
     255             :   {
     256             :     case t_POLMOD:
     257     1196888 :       b = gel(x,1);
     258     1196888 :       a = gel(x,2); ta = typ(a);
     259     1196888 :       if (is_const_t(ta))
     260        2744 :         return scalar_ZX(degpol(T)? Rg_to_Fp(a, p): gen_0, v);
     261     1194144 :       b = RgX_to_FpX(b, p); if (varn(b) != v) break;
     262     1194144 :       a = RgX_to_FpX(a, p);
     263     1194144 :       if (ZX_equal(b,get_FpX_mod(T)) || signe(FpX_rem(b,T,p))==0)
     264     1194144 :         return FpX_rem(a, T, p);
     265           0 :       break;
     266             :     case t_POL:
     267        4613 :       if (varn(x) != v) break;
     268        4613 :       return FpX_rem(RgX_to_FpX(x,p), T, p);
     269             :     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     3332712 : RgX_to_FpX(GEN x, GEN p)
     279             : {
     280             :   long i, l;
     281     3332712 :   GEN z = cgetg_copy(x, &l); z[1] = x[1];
     282     3332712 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     283     3332712 :   return FpX_renormalize(z, l);
     284             : }
     285             : 
     286             : GEN
     287        1022 : RgV_to_FpV(GEN x, GEN p)
     288        1022 : { pari_APPLY_type(t_VEC, Rg_to_Fp(gel(x,i), p)) }
     289             : 
     290             : GEN
     291      922209 : RgC_to_FpC(GEN x, GEN p)
     292      922209 : { pari_APPLY_type(t_COL, Rg_to_Fp(gel(x,i), p)) }
     293             : 
     294             : GEN
     295      129343 : RgM_to_FpM(GEN x, GEN p)
     296      129343 : { pari_APPLY_same(RgC_to_FpC(gel(x,i), p)) }
     297             : 
     298             : GEN
     299      281602 : RgV_to_Flv(GEN x, ulong p)
     300      281602 : { pari_APPLY_ulong(Rg_to_Fl(gel(x,i), p)) }
     301             : 
     302             : GEN
     303      114236 : RgM_to_Flm(GEN x, ulong p)
     304      114236 : { pari_APPLY_same(RgV_to_Flv(gel(x,i), p)) }
     305             : 
     306             : GEN
     307         476 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
     308             : {
     309         476 :   long i, l = lg(x);
     310         476 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     311         476 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
     312         476 :   return FpXQX_renormalize(z, l);
     313             : }
     314             : GEN
     315        1267 : RgX_to_FqX(GEN x, GEN T, GEN p)
     316             : {
     317        1267 :   long i, l = lg(x);
     318        1267 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     319        1267 :   if (T)
     320         602 :     for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
     321             :   else
     322         665 :     for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     323        1267 :   return FpXQX_renormalize(z, l);
     324             : }
     325             : 
     326             : GEN
     327      218862 : RgC_to_FqC(GEN x, GEN T, GEN p)
     328             : {
     329      218862 :   long i, l = lg(x);
     330      218862 :   GEN z = cgetg(l, t_COL);
     331      218862 :   if (T)
     332      218862 :     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      218862 :   return z;
     336             : }
     337             : 
     338             : GEN
     339       52318 : RgM_to_FqM(GEN x, GEN T, GEN p)
     340       52318 : { 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        1596 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
     359             : {
     360        1596 :   long i, lz = lg(y);
     361             :   GEN z;
     362        1596 :   if (!T) return FpX_Fp_add(y, x, p);
     363        1596 :   if (lz == 2) return scalarpol(x, varn(y));
     364        1596 :   z = cgetg(lz,t_POL); z[1] = y[1];
     365        1596 :   gel(z,2) = Fq_add(gel(y,2),x, T, p);
     366        1596 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     367             :   else
     368         287 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     369        1596 :   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         926 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     384        1048 :   return z;
     385             : }
     386             : 
     387             : GEN
     388      144652 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
     389             : {
     390             :   long i, lP;
     391      144652 :   GEN res = cgetg_copy(P, &lP); res[1] = P[1];
     392      144652 :   for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
     393      144652 :   gel(res,lP-1) = gen_1; return res;
     394             : }
     395             : 
     396             : GEN
     397        3840 : FpXQX_normalize(GEN z, GEN T, GEN p)
     398             : {
     399             :   GEN lc;
     400        3840 :   if (lg(z) == 2) return z;
     401        3826 :   lc = leading_coeff(z);
     402        3826 :   if (typ(lc) == t_POL)
     403             :   {
     404        1841 :     if (lg(lc) > 3) /* non-constant */
     405        1592 :       return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
     406             :     /* constant */
     407         249 :     lc = gel(lc,2);
     408         249 :     z = shallowcopy(z);
     409         249 :     gel(z, lg(z)-1) = lc;
     410             :   }
     411             :   /* lc a t_INT */
     412        2234 :   if (equali1(lc)) return z;
     413          50 :   return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
     414             : }
     415             : 
     416             : GEN
     417      127379 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
     418             : {
     419             :   pari_sp av;
     420             :   GEN p1, r;
     421      127379 :   long j, i=lg(x)-1;
     422      127379 :   if (i<=2)
     423       26348 :     return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
     424      101031 :   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      297171 :   for (i--; i>=2; i=j-1)
     428             :   {
     429      196588 :     for (j=i; !signe(gel(x,j)); j--)
     430         448 :       if (j==2)
     431             :       {
     432         301 :         if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
     433         301 :         return gerepileupto(av, Fq_mul(p1,y, T, p));
     434             :       }
     435      196140 :     r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
     436      196140 :     p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
     437             :   }
     438      100730 :   return gerepileupto(av, p1);
     439             : }
     440             : 
     441             : GEN
     442       31591 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
     443             : {
     444       31591 :   long i, lb = lg(Q);
     445             :   GEN z;
     446       31591 :   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       14497 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
     459             : {
     460       14497 :   pari_sp av = avma;
     461       14497 :   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      348474 : monomial(GEN a, long d, long v)
     468             : {
     469             :   long i, n;
     470             :   GEN P;
     471      348474 :   if (d < 0) {
     472           0 :     if (isrationalzero(a)) return pol_0(v);
     473           0 :     retmkrfrac(a, pol_xn(-d, v));
     474             :   }
     475      348474 :   if (gequal0(a))
     476             :   {
     477        8379 :     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      340095 :     n = d+2; P = cgetg(n+1, t_POL);
     484      340095 :     P[1] = evalsigne(1) | evalvarn(v);
     485             :   }
     486      340095 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     487      340095 :   gel(P,i) = a; return P;
     488             : }
     489             : GEN
     490     7949173 : monomialcopy(GEN a, long d, long v)
     491             : {
     492             :   long i, n;
     493             :   GEN P;
     494     7949173 :   if (d < 0) {
     495          14 :     if (isrationalzero(a)) return pol_0(v);
     496          14 :     retmkrfrac(gcopy(a), pol_xn(-d, v));
     497             :   }
     498     7949159 :   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     7949152 :     n = d+2; P = cgetg(n+1, t_POL);
     507     7949152 :     P[1] = evalsigne(1) | evalvarn(v);
     508             :   }
     509     7949152 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     510     7949152 :   gel(P,i) = gcopy(a); return P;
     511             : }
     512             : GEN
     513       20223 : pol_x_powers(long N, long v)
     514             : {
     515       20223 :   GEN L = cgetg(N+1,t_VEC);
     516             :   long i;
     517       20223 :   for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
     518       20223 :   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     6962290 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     541             : {
     542             :   (void)T;
     543     6962290 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     544             :   {
     545     2485111 :     case 0: return Fp_add(x,y,p);
     546      203791 :     case 1: return FpX_Fp_add(x,y,p);
     547      327041 :     case 2: return FpX_Fp_add(y,x,p);
     548     3946347 :     case 3: return FpX_add(x,y,p);
     549             :   }
     550             :   return NULL;/*LCOV_EXCL_LINE*/
     551             : }
     552             : 
     553             : GEN
     554     4694280 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
     555             : {
     556             :   (void)T;
     557     4694280 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     558             :   {
     559      167380 :     case 0: return Fp_sub(x,y,p);
     560        2191 :     case 1: return FpX_Fp_sub(x,y,p);
     561       10066 :     case 2: return Fp_FpX_sub(x,y,p);
     562     4514643 :     case 3: return FpX_sub(x,y,p);
     563             :   }
     564             :   return NULL;/*LCOV_EXCL_LINE*/
     565             : }
     566             : 
     567             : GEN
     568      471433 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
     569             : {
     570             :   (void)T;
     571      471433 :   return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
     572             : }
     573             : 
     574             : GEN
     575       12899 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
     576             : {
     577             :   (void)T;
     578       12899 :   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    42744209 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
     584             : {
     585    42744209 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     586             :   {
     587     2539901 :     case 0: return Fp_mul(x,y,p);
     588       71525 :     case 1: return FpX_Fp_mul(x,y,p);
     589      131034 :     case 2: return FpX_Fp_mul(y,x,p);
     590    40001749 :     case 3: if (T) return FpXQ_mul(x,y,T,p);
     591     2756744 :             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      774479 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
     599             : {
     600             :   (void) T;
     601      774479 :   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       56860 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
     607             : {
     608             :   (void)T;
     609       56860 :   return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
     610       56860 :                           : Fp_mul(x,y,p);
     611             : }
     612             : /* If T==NULL do not reduce*/
     613             : GEN
     614      270836 : Fq_sqr(GEN x, GEN T, GEN p)
     615             : {
     616      270836 :   if (typ(x) == t_POL)
     617             :   {
     618       11843 :     if (T) return FpXQ_sqr(x,T,p);
     619           0 :     else return FpX_sqr(x,p);
     620             :   }
     621             :   else
     622      258993 :     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       31281 : Fq_inv(GEN x, GEN pol, GEN p)
     641             : {
     642       31281 :   if (typ(x) == t_INT) return Fp_inv(x,p);
     643       25278 :   return FpXQ_inv(x,pol,p);
     644             : }
     645             : 
     646             : GEN
     647      516467 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
     648             : {
     649      516467 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     650             :   {
     651      486976 :     case 0: return Fp_div(x,y,p);
     652       23975 :     case 1: return FpX_Fp_mul(x,Fp_inv(y,p),p);
     653         280 :     case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
     654        5236 :     case 3: return FpXQ_div(x,y,pol,p);
     655             :   }
     656             :   return NULL;/*LCOV_EXCL_LINE*/
     657             : }
     658             : 
     659             : GEN
     660       20013 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
     661             : {
     662       20013 :   if (typ(x) == t_INT) return Fp_pow(x,n,p);
     663        8729 :   return FpXQ_pow(x,n,pol,p);
     664             : }
     665             : 
     666             : GEN
     667       14770 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
     668             : {
     669       14770 :   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      709301 : Fq_sqrt(GEN x, GEN T, GEN p)
     675             : {
     676      709301 :   if (typ(x) == t_INT)
     677             :   {
     678      698670 :     if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
     679         301 :     x = scalarpol_shallow(x, get_FpX_var(T));
     680             :   }
     681       10932 :   return FpXQ_sqrt(x,T,p);
     682             : }
     683             : GEN
     684       60538 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
     685             : {
     686       60538 :   if (typ(x) == t_INT)
     687             :   {
     688             :     long d;
     689       60300 :     if (!T) return Fp_sqrtn(x,n,p,zeta);
     690         562 :     d = get_FpX_degree(T);
     691         562 :     if (ugcdiu(n,d) == 1)
     692             :     {
     693         401 :       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         394 :       if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
     696         373 :         return Fp_sqrtn(x,n,p,zeta);
     697             :     }
     698         182 :     x = scalarpol(x, get_FpX_var(T)); /* left on stack */
     699             :   }
     700         420 :   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      632794 : _Fq_red(void *E, GEN x)
     710      632794 : { struct _Fq_field *s = (struct _Fq_field *)E;
     711      632794 :   return Fq_red(x, s->T, s->p);
     712             : }
     713             : 
     714             : static GEN
     715     1225798 : _Fq_add(void *E, GEN x, GEN y)
     716             : {
     717             :   (void) E;
     718     1225798 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     719             :   {
     720        3094 :     case 0: return addii(x,y);
     721           0 :     case 1: return ZX_Z_add(x,y);
     722       25620 :     case 2: return ZX_Z_add(y,x);
     723     1197084 :     default: return ZX_add(x,y);
     724             :   }
     725             : }
     726             : 
     727             : static GEN
     728      207669 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
     729             : 
     730             : static GEN
     731     1309546 : _Fq_mul(void *E, GEN x, GEN y)
     732             : {
     733             :   (void) E;
     734     1309546 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     735             :   {
     736        4137 :     case 0: return mulii(x,y);
     737       36897 :     case 1: return ZX_Z_mul(x,y);
     738         119 :     case 2: return ZX_Z_mul(y,x);
     739     1268393 :     default: return ZX_mul(x,y);
     740             :   }
     741             : }
     742             : 
     743             : static GEN
     744        6055 : _Fq_inv(void *E, GEN x)
     745        6055 : { struct _Fq_field *s = (struct _Fq_field *)E;
     746        6055 :   return Fq_inv(x,s->T,s->p);
     747             : }
     748             : 
     749             : static int
     750      114387 : _Fq_equal0(GEN x) { return signe(x)==0; }
     751             : 
     752             : static GEN
     753       31822 : _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        3290 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
     759             : {
     760        3290 :   GEN z = new_chunk(sizeof(struct _Fq_field));
     761        3290 :   struct _Fq_field *e = (struct _Fq_field *) z;
     762        3290 :   e->T = T; e->p  = p; *E = (void*)e;
     763        3290 :   return &Fq_field;
     764             : }
     765             : 
     766             : /*******************************************************************/
     767             : /*                                                                 */
     768             : /*                             Fq[X]                               */
     769             : /*                                                                 */
     770             : /*******************************************************************/
     771             : /* P(X + c) */
     772             : GEN
     773           0 : FpX_translate(GEN P, GEN c, GEN p)
     774             : {
     775           0 :   pari_sp av = avma;
     776             :   GEN Q, *R;
     777             :   long i, k, n;
     778             : 
     779           0 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     780           0 :   Q = leafcopy(P);
     781           0 :   R = (GEN*)(Q+2); n = degpol(P);
     782           0 :   for (i=1; i<=n; i++)
     783             :   {
     784           0 :     for (k=n-i; k<n; k++)
     785           0 :       R[k] = Fp_add(R[k], Fp_mul(c, R[k+1], p), p);
     786             : 
     787           0 :     if (gc_needed(av,2))
     788             :     {
     789           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpX_translate, i = %ld/%ld", i,n);
     790           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     791             :     }
     792             :   }
     793           0 :   return gerepilecopy(av, FpX_renormalize(Q, lg(Q)));
     794             : }
     795             : /* P(X + c), c an Fq */
     796             : GEN
     797       34167 : FqX_translate(GEN P, GEN c, GEN T, GEN p)
     798             : {
     799       34167 :   pari_sp av = avma;
     800             :   GEN Q, *R;
     801             :   long i, k, n;
     802             : 
     803             :   /* signe works for t_(INT|POL) */
     804       34167 :   if (!signe(P) || !signe(c)) return RgX_copy(P);
     805       34167 :   Q = leafcopy(P);
     806       34167 :   R = (GEN*)(Q+2); n = degpol(P);
     807      151781 :   for (i=1; i<=n; i++)
     808             :   {
     809      439299 :     for (k=n-i; k<n; k++)
     810      321685 :       R[k] = Fq_add(R[k], Fq_mul(c, R[k+1], T, p), T, p);
     811             : 
     812      117614 :     if (gc_needed(av,2))
     813             :     {
     814           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FqX_translate, i = %ld/%ld", i,n);
     815           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     816             :     }
     817             :   }
     818       34167 :   return gerepilecopy(av, FpXQX_renormalize(Q, lg(Q)));
     819             : }
     820             : 
     821             : GEN
     822         665 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
     823             : {
     824         665 :   pari_sp ltop = avma;
     825             :   long k;
     826             :   GEN W;
     827         665 :   if (lgefint(p) == 3)
     828             :   {
     829         591 :     ulong pp = p[2];
     830         591 :     GEN Tl = ZX_to_Flx(T, pp);
     831         591 :     GEN Vl = FqV_to_FlxV(V, T, p);
     832         591 :     Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
     833         591 :     return gerepileupto(ltop, FlxX_to_ZXX(Tl));
     834             :   }
     835          74 :   W = cgetg(lg(V),t_VEC);
     836         402 :   for(k=1; k < lg(V); k++)
     837         328 :     gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
     838          74 :   return gerepileupto(ltop, FpXQXV_prod(W, T, p));
     839             : }
     840             : 
     841             : GEN
     842      124593 : FqV_red(GEN x, GEN T, GEN p)
     843      124593 : { pari_APPLY_same(Fq_red(gel(x,i), T, p)) }
     844             : 
     845             : GEN
     846           0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
     847             : {
     848           0 :   if (!T) return FpC_add(x, y, p);
     849           0 :   pari_APPLY_type(t_COL, Fq_add(gel(x,i), gel(y,i), T, p))
     850             : }
     851             : 
     852             : GEN
     853           0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
     854             : {
     855           0 :   if (!T) return FpC_sub(x, y, p);
     856           0 :   pari_APPLY_type(t_COL, Fq_sub(gel(x,i), gel(y,i), T, p))
     857             : }
     858             : 
     859             : GEN
     860           0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
     861             : {
     862           0 :   if (!T) return FpC_Fp_mul(x, y, p);
     863           0 :   pari_APPLY_type(t_COL, Fq_mul(gel(x,i),y,T,p))
     864             : }
     865             : 
     866             : GEN
     867         591 : FqV_to_FlxV(GEN x, GEN T, GEN pp)
     868             : {
     869         591 :   long vT = evalvarn(get_FpX_var(T));
     870         591 :   ulong p = pp[2];
     871         591 :   pari_APPLY_type(t_VEC, typ(gel(x,i))==t_INT?  Z_to_Flx(gel(x,i), p, vT)
     872             :                                              : ZX_to_Flx(gel(x,i), p))
     873             : }
     874             : 
     875             : GEN
     876       59165 : FqC_to_FlxC(GEN x, GEN T, GEN pp)
     877             : {
     878       59165 :   long vT = evalvarn(get_FpX_var(T));
     879       59166 :   ulong p = pp[2];
     880       59166 :   pari_APPLY_type(t_COL, typ(gel(x,i))==t_INT?  Z_to_Flx(gel(x,i), p, vT)
     881             :                                              : ZX_to_Flx(gel(x,i), p))
     882             : }
     883             : 
     884             : GEN
     885        9463 : FqM_to_FlxM(GEN x, GEN T, GEN p)
     886        9463 : { pari_APPLY_same(FqC_to_FlxC(gel(x,i), T, p)) }
     887             : 
     888             : GEN
     889        2402 : FpXC_center(GEN x, GEN p, GEN pov2)
     890        2402 : { pari_APPLY_type(t_COL, FpX_center(gel(x,i), p, pov2)) }
     891             : 
     892             : GEN
     893        1081 : FpXM_center(GEN x, GEN p, GEN pov2)
     894        1081 : { pari_APPLY_same(FpXC_center(gel(x,i), p, pov2)) }
     895             : 
     896             : /*******************************************************************/
     897             : /*                                                                 */
     898             : /*                          GENERIC CRT                            */
     899             : /*                                                                 */
     900             : /*******************************************************************/
     901             : 
     902             : static long
     903      309396 : get_nbprimes(ulong bound, ulong *pt_start)
     904             : {
     905             : #ifdef LONG_IS_64BIT
     906      264936 :   ulong pstart = 4611686018427388039UL;
     907             : #else
     908       44460 :   ulong pstart = 1073741827UL;
     909             : #endif
     910      309396 :   if (pt_start) *pt_start = pstart;
     911      309396 :   return (bound/expu(pstart))+1;
     912             : }
     913             : 
     914             : static GEN
     915      749903 : primelist_disc(ulong *p, long n, GEN dB)
     916             : {
     917      749903 :   ulong u = 0;
     918      749903 :   GEN P = cgetg(n+1, t_VECSMALL);
     919             :   long i;
     920      749903 :   if (dB && typ(dB)==t_VECSMALL) { u = uel(dB,1); dB = NULL; }
     921     2295973 :   for (i=1; i <= n; i++, *p = unextprime(*p+1))
     922             :   {
     923     1546070 :     if (dB && umodiu(dB, *p)==0) { i--; continue; }
     924     1546070 :     if (u && *p%u!=1) { i--; continue; }
     925     1541312 :     P[i] = *p;
     926             :   }
     927      749903 :   return P;
     928             : }
     929             : 
     930             : void
     931      225003 : gen_inccrt(const char *str, GEN worker, GEN dB, long n, long mmin,
     932             :            ulong *p, GEN *pt_H, GEN *pt_mod, GEN crt(GEN, GEN, GEN*),
     933             :            GEN center(GEN, GEN, GEN))
     934             : {
     935      225003 :   pari_sp av = avma;
     936             :   long m;
     937             :   GEN  H, P, mod;
     938             :   pari_timer ti;
     939      225003 :   if (!*p) (void) get_nbprimes(1, p);
     940      225003 :   m = minss(mmin, n);
     941      225003 :   if (DEBUGLEVEL > 4)
     942             :   {
     943           0 :       timer_start(&ti);
     944           0 :       err_printf("%s: nb primes: %ld\n",str, n);
     945             :   }
     946      225003 :   if (m == 1)
     947             :   {
     948      164827 :     GEN P = primelist_disc(p, n, dB);
     949      164827 :     GEN done = closure_callgen1(worker, P);
     950      164827 :     H = gel(done,1);
     951      164827 :     mod = gel(done,2);
     952      164827 :     if (!*pt_H && center) H = center(H, mod, shifti(mod,-1));
     953      164827 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
     954             :   }
     955             :   else
     956             :   {
     957       60176 :     long i, s = (n+m-1)/m, r = m - (m*s-n), di = 0;
     958             :     struct pari_mt pt;
     959       60176 :     long pending = 0;
     960       60176 :     H = cgetg(m+1, t_VEC); P = cgetg(m+1, t_VEC);
     961       60176 :     mt_queue_start_lim(&pt, worker, m);
     962      716095 :     for (i=1; i<=m || pending; i++)
     963             :     {
     964             :       GEN done;
     965      655919 :       GEN pr = i <= m ? mkvec(primelist_disc(p, i<=r ? s: s-1, dB)): NULL;
     966      655919 :       mt_queue_submit(&pt, i, pr);
     967      655919 :       done = mt_queue_get(&pt, NULL, &pending);
     968      655919 :       if (done)
     969             :       {
     970      585076 :         di++;
     971      585076 :         gel(H, di) = gel(done,1);
     972      585076 :         gel(P, di) = gel(done,2);
     973      585076 :         if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
     974             :       }
     975             :     }
     976       60176 :     mt_queue_end(&pt);
     977       60176 :     if (DEBUGLEVEL>5) err_printf("\n");
     978       60176 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
     979       60176 :     H = crt(H, P, &mod);
     980       60176 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: chinese", str);
     981             :   }
     982      225003 :   if (*pt_H)
     983       14140 :     H = crt(mkvec2(*pt_H, H), mkvec2(*pt_mod, mod), &mod);
     984      225003 :   *pt_H = H;
     985      225003 :   *pt_mod = mod;
     986      225003 :   gerepileall(av, 2, pt_H, pt_mod);
     987      225003 : }
     988             : 
     989             : GEN
     990       98533 : gen_crt(const char *str, GEN worker, GEN dB, ulong bound, long mmin, GEN *pt_mod,
     991             :         GEN crt(GEN, GEN, GEN*), GEN center(GEN, GEN, GEN))
     992             : {
     993       98533 :   ulong p = 0;
     994       98533 :   GEN mod = gen_1, H = NULL;
     995       98533 :   bound++;
     996      295599 :   while ((ulong)expi(mod) < bound)
     997             :   {
     998       98533 :     long n = get_nbprimes(bound-expi(mod), NULL);
     999       98533 :     gen_inccrt(str, worker, dB, n, mmin, &p, &H, &mod, crt, center);
    1000             :   }
    1001       98533 :   if (pt_mod) *pt_mod = mod;
    1002       98533 :   return H;
    1003             : }
    1004             : 
    1005             : /*******************************************************************/
    1006             : /*                                                                 */
    1007             : /*                          MODULAR GCD                            */
    1008             : /*                                                                 */
    1009             : /*******************************************************************/
    1010             : /* return z = a mod q, b mod p (p,q) = 1; qinv = 1/q mod p; a in ]-q,q] */
    1011             : static GEN
    1012     2111224 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq, GEN pq2)
    1013             : {
    1014     2111224 :   ulong d, amod = umodiu(a, p);
    1015     2111224 :   pari_sp av = avma;
    1016             :   GEN ax;
    1017             : 
    1018     2111224 :   if (b == amod) return NULL;
    1019     1308503 :   d = Fl_mul(Fl_sub(b, amod, p), qinv, p); /* != 0 */
    1020     1308503 :   if (d >= 1 + (p>>1))
    1021      650577 :     ax = subii(a, mului(p-d, q));
    1022             :   else
    1023             :   {
    1024      657926 :     ax = addii(a, mului(d, q)); /* in ]0, pq[ assuming a in ]-q,q[ */
    1025      657926 :     if (cmpii(ax,pq2) > 0) ax = subii(ax,pq);
    1026             :   }
    1027     1308503 :   return gerepileuptoint(av, ax);
    1028             : }
    1029             : GEN
    1030         364 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
    1031             : GEN
    1032     3179533 : ZX_init_CRT(GEN Hp, ulong p, long v)
    1033             : {
    1034     3179533 :   long i, l = lg(Hp), lim = (long)(p>>1);
    1035     3179533 :   GEN H = cgetg(l, t_POL);
    1036     3179533 :   H[1] = evalsigne(1) | evalvarn(v);
    1037    11298053 :   for (i=2; i<l; i++)
    1038     8118520 :     gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
    1039     3179533 :   return ZX_renormalize(H,l);
    1040             : }
    1041             : 
    1042             : GEN
    1043       94507 : ZM_init_CRT(GEN Hp, ulong p)
    1044             : {
    1045       94507 :   long i,j, m, l = lg(Hp), lim = (long)(p>>1);
    1046       94507 :   GEN c, cp, H = cgetg(l, t_MAT);
    1047       94507 :   if (l==1) return H;
    1048       50330 :   m = lgcols(Hp);
    1049      128197 :   for (j=1; j<l; j++)
    1050             :   {
    1051       77867 :     cp = gel(Hp,j);
    1052       77867 :     c = cgetg(m, t_COL);
    1053       77867 :     gel(H,j) = c;
    1054       77867 :     for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
    1055             :   }
    1056       50330 :   return H;
    1057             : }
    1058             : 
    1059             : int
    1060        7511 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
    1061             : {
    1062        7511 :   GEN h, q = *ptq, qp = muliu(q,p);
    1063        7511 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1064        7511 :   int stable = 1;
    1065        7511 :   h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp,shifti(qp,-1));
    1066        7511 :   if (h) { *H = h; stable = 0; }
    1067        7511 :   *ptq = qp; return stable;
    1068             : }
    1069             : 
    1070             : static int
    1071      222335 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
    1072             : {
    1073      222335 :   GEN H = *ptH, h, qp2 = shifti(qp,-1);
    1074      222335 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1075      222335 :   long i, l = lg(H), lp = lg(Hp);
    1076      222335 :   int stable = 1;
    1077             : 
    1078      222335 :   if (l < lp)
    1079             :   { /* degree increases */
    1080           0 :     GEN x = cgetg(lp, t_POL);
    1081           0 :     for (i=1; i<l; i++)  x[i] = H[i];
    1082           0 :     for (   ; i<lp; i++) gel(x,i) = gen_0;
    1083           0 :     *ptH = H = x;
    1084           0 :     stable = 0;
    1085      222335 :   } else if (l > lp)
    1086             :   { /* degree decreases */
    1087           0 :     GEN x = cgetg(l, t_VECSMALL);
    1088           0 :     for (i=1; i<lp; i++)  x[i] = Hp[i];
    1089           0 :     for (   ; i<l; i++) x[i] = 0;
    1090           0 :     Hp = x; lp = l;
    1091             :   }
    1092     1807270 :   for (i=2; i<lp; i++)
    1093             :   {
    1094     1584935 :     h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp,qp2);
    1095     1584935 :     if (h) { gel(H,i) = h; stable = 0; }
    1096             :   }
    1097      222335 :   (void)ZX_renormalize(H,lp);
    1098      222335 :   return stable;
    1099             : }
    1100             : 
    1101             : int
    1102       13660 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
    1103             : {
    1104       13660 :   GEN q = *ptq, qp = muliu(q,p);
    1105       13660 :   int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
    1106       13660 :   *ptq = qp; return stable;
    1107             : }
    1108             : 
    1109             : int
    1110       17975 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1111             : {
    1112       17975 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1113       17975 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1114       17975 :   long i,j, l = lg(H), m = lgcols(H);
    1115       17975 :   int stable = 1;
    1116       46613 :   for (j=1; j<l; j++)
    1117      453696 :     for (i=1; i<m; i++)
    1118             :     {
    1119      425058 :       h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp,qp2);
    1120      425058 :       if (h) { gcoeff(H,i,j) = h; stable = 0; }
    1121             :     }
    1122       17975 :   *ptq = qp; return stable;
    1123             : }
    1124             : 
    1125             : GEN
    1126         686 : ZXM_init_CRT(GEN Hp, long deg, ulong p)
    1127             : {
    1128             :   long i, j, k;
    1129             :   GEN H;
    1130         686 :   long m, l = lg(Hp), lim = (long)(p>>1), n;
    1131         686 :   H = cgetg(l, t_MAT);
    1132         686 :   if (l==1) return H;
    1133         686 :   m = lgcols(Hp);
    1134         686 :   n = deg + 3;
    1135        2548 :   for (j=1; j<l; j++)
    1136             :   {
    1137        1862 :     GEN cp = gel(Hp,j);
    1138        1862 :     GEN c = cgetg(m, t_COL);
    1139        1862 :     gel(H,j) = c;
    1140       25690 :     for (i=1; i<m; i++)
    1141             :     {
    1142       23828 :       GEN dp = gel(cp, i);
    1143       23828 :       long l = lg(dp);
    1144       23828 :       GEN d = cgetg(n, t_POL);
    1145       23828 :       gel(c, i) = d;
    1146       23828 :       d[1] = dp[1];
    1147       47075 :       for (k=2; k<l; k++)
    1148       23247 :         gel(d,k) = stoi(Fl_center(dp[k], p, lim));
    1149       48643 :       for (   ; k<n; k++)
    1150       24815 :         gel(d,k) = gen_0;
    1151             :     }
    1152             :   }
    1153         686 :   return H;
    1154             : }
    1155             : 
    1156             : int
    1157         404 : ZXM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1158             : {
    1159         404 :   GEN v, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1160         404 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1161         404 :   long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
    1162         404 :   int stable = 1;
    1163        2522 :   for (j=1; j<l; j++)
    1164       48968 :     for (i=1; i<m; i++)
    1165             :     {
    1166       46850 :       GEN h = gmael(H,j,i), hp = gmael(Hp,j,i);
    1167       46850 :       long lh = lg(hp);
    1168       95573 :       for (k=2; k<lh; k++)
    1169             :       {
    1170       48723 :         v = Fl_chinese_coprime(gel(h,k),uel(hp,k),q,p,qinv,qp,qp2);
    1171       48723 :         if (v) { gel(h,k) = v; stable = 0; }
    1172             :       }
    1173       91847 :       for (; k<n; k++)
    1174             :       {
    1175       44997 :         v = Fl_chinese_coprime(gel(h,k),0,q,p,qinv,qp,qp2);
    1176       44997 :         if (v) { gel(h,k) = v; stable = 0; }
    1177             :       }
    1178             :     }
    1179         404 :   *ptq = qp; return stable;
    1180             : }
    1181             : 
    1182             : /* record the degrees of Euclidean remainders (make them as large as
    1183             :  * possible : smaller values correspond to a degenerate sequence) */
    1184             : static void
    1185        1722 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
    1186             : {
    1187             :   long da,db,dc, ind;
    1188        1722 :   pari_sp av = avma;
    1189             : 
    1190        1722 :   if (lgpol(a)==0 || lgpol(b)==0) return;
    1191        1701 :   da = degpol(a);
    1192        1701 :   db = degpol(b);
    1193        1701 :   if (db > da)
    1194           0 :   { swapspec(a,b, da,db); }
    1195        1701 :   else if (!da) return;
    1196        1701 :   ind = 0;
    1197       10521 :   while (db)
    1198             :   {
    1199        7119 :     GEN c = Flx_rem(a,b, p);
    1200        7119 :     a = b; b = c; dc = degpol(c);
    1201        7119 :     if (dc < 0) break;
    1202             : 
    1203        7119 :     ind++;
    1204        7119 :     if (dc > dglist[ind]) dglist[ind] = dc;
    1205        7119 :     if (gc_needed(av,2))
    1206             :     {
    1207           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1208           0 :       gerepileall(av, 2, &a,&b);
    1209             :     }
    1210        7119 :     db = dc; /* = degpol(b) */
    1211             :   }
    1212        1701 :   if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
    1213        1701 :   set_avma(av);
    1214             : }
    1215             : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
    1216             :  * "generic" degree sequence as given by dglist, set *Ci and return
    1217             :  * resultant(a,b). Modular version of Collins's subresultant */
    1218             : static ulong
    1219        8707 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
    1220             : {
    1221             :   long da,db,dc, ind;
    1222        8707 :   ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
    1223        8707 :   int s = 1;
    1224        8707 :   pari_sp av = avma;
    1225             : 
    1226        8707 :   *C0 = 1; *C1 = 0;
    1227        8707 :   if (lgpol(a)==0 || lgpol(b)==0) return 0;
    1228        8665 :   da = degpol(a);
    1229        8665 :   db = degpol(b);
    1230        8665 :   if (db > da)
    1231             :   {
    1232           0 :     swapspec(a,b, da,db);
    1233           0 :     if (both_odd(da,db)) s = -s;
    1234             :   }
    1235        8665 :   else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
    1236        8665 :   ind = 0;
    1237       51490 :   while (db)
    1238             :   { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
    1239             :      * da = deg a, db = deg b */
    1240       34552 :     GEN c = Flx_rem(a,b, p);
    1241       34552 :     long delta = da - db;
    1242             : 
    1243       34552 :     if (both_odd(da,db)) s = -s;
    1244       34552 :     lb = Fl_mul(b[db+2], cb, p);
    1245       34552 :     a = b; b = c; dc = degpol(c);
    1246       34552 :     ind++;
    1247       34552 :     if (dc != dglist[ind]) return gc_ulong(av,0); /* degenerates */
    1248       34160 :     if (g == h)
    1249             :     { /* frequent */
    1250       31878 :       ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
    1251       31878 :       ca = cb;
    1252       31878 :       cb = cc;
    1253             :     }
    1254             :     else
    1255             :     {
    1256        2282 :       ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
    1257        2282 :       ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
    1258        2282 :       ca = cb;
    1259        2282 :       cb = Fl_div(cc, ghdelta, p);
    1260             :     }
    1261       34160 :     da = db; /* = degpol(a) */
    1262       34160 :     db = dc; /* = degpol(b) */
    1263             : 
    1264       34160 :     g = lb;
    1265       34160 :     if (delta == 1)
    1266       24487 :       h = g; /* frequent */
    1267             :     else
    1268        9673 :       h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
    1269             : 
    1270       34160 :     if (gc_needed(av,2))
    1271             :     {
    1272           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1273           0 :       gerepileall(av, 2, &a,&b);
    1274             :     }
    1275             :   }
    1276        8273 :   if (da > 1) return 0; /* Failure */
    1277             :   /* last non-constant polynomial has degree 1 */
    1278        8273 :   *C0 = Fl_mul(ca, a[2], p);
    1279        8273 :   *C1 = Fl_mul(ca, a[3], p);
    1280        8273 :   res = Fl_mul(cb, b[2], p);
    1281        8273 :   if (s == -1) res = p - res;
    1282        8273 :   return gc_ulong(av,res);
    1283             : }
    1284             : 
    1285             : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
    1286             :  * Return 0 in case of degree drop. */
    1287             : static GEN
    1288       10429 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
    1289             : {
    1290             :   GEN z;
    1291       10429 :   long i, lb = lg(Q);
    1292       10429 :   ulong leadz = Flx_eval(leading_coeff(Q), x, p);
    1293       10429 :   long vs=mael(Q,2,1);
    1294       10429 :   if (!leadz) return zero_Flx(vs);
    1295             : 
    1296       10366 :   z = cgetg(lb, t_VECSMALL); z[1] = vs;
    1297       10366 :   for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
    1298       10366 :   z[i] = leadz; return z;
    1299             : }
    1300             : 
    1301             : GEN
    1302       20062 : FpXY_Fq_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
    1303             : {
    1304       20062 :   pari_sp av = avma;
    1305       20062 :   long i, lb = lg(Q);
    1306             :   GEN z;
    1307       20062 :   if (!T) return FpXY_evaly(Q, y, p, vx);
    1308        1232 :   if (lb == 2) return pol_0(vx);
    1309        1232 :   z = gel(Q, lb-1);
    1310        1232 :   if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
    1311             : 
    1312        1232 :   if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
    1313       26572 :   for (i=lb-2; i>=2; i--)
    1314             :   {
    1315       25340 :     GEN c = gel(Q,i);
    1316       25340 :     z = FqX_Fq_mul(z, y, T, p);
    1317       25340 :     z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
    1318             :   }
    1319        1232 :   return gerepileupto(av, z);
    1320             : }
    1321             : 
    1322             : static GEN
    1323       15225 : ZX_norml1(GEN x)
    1324             : {
    1325       15225 :   long i, l = lg(x);
    1326             :   GEN s;
    1327             : 
    1328       15225 :   if (l == 2) return gen_0;
    1329        8603 :   s = gel(x, l-1); /* != 0 */
    1330       31626 :   for (i = l-2; i > 1; i--) {
    1331       23023 :     GEN xi = gel(x,i);
    1332       23023 :     if (!signe(x)) continue;
    1333       23023 :     s = addii_sign(s,1, xi,1);
    1334             :   }
    1335        8603 :   return s;
    1336             : }
    1337             : 
    1338             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1339             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1340             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1341             :  * Return e such that Res(A, B) < 2^e */
    1342             : ulong
    1343       79332 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
    1344             : {
    1345       79332 :   pari_sp av = avma, av2;
    1346       79332 :   GEN a = gen_0, b = gen_0;
    1347       79332 :   long i , lA = lg(A), lB = lg(B);
    1348             :   double loga, logb;
    1349      880419 :   for (i=2; i<lA; i++)
    1350             :   {
    1351      801087 :     a = addii(a, sqri(gel(A,i)));
    1352      801087 :     if (gc_needed(av,1))
    1353             :     {
    1354           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1355           0 :       a = gerepileupto(av, a);
    1356             :     }
    1357             :   }
    1358       79332 :   a = gerepileuptoint(av, a);
    1359       79332 :   av2 = avma;
    1360      806886 :   for (i=2; i<lB; i++)
    1361             :   {
    1362      727554 :     GEN t = gel(B,i);
    1363      727554 :     if (typ(t) == t_POL) t = ZX_norml1(t);
    1364      727554 :     b = addii(b, sqri(t));
    1365      727554 :     if (gc_needed(av2,1))
    1366             :     {
    1367           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1368           0 :       b = gerepileupto(av2, b);
    1369             :     }
    1370             :   }
    1371       79332 :   loga = dbllog2(a);
    1372       79332 :   logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
    1373       79332 :   i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
    1374       79332 :   return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
    1375             : }
    1376             : 
    1377             : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
    1378             : static ulong
    1379      247272 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong la)
    1380             : {
    1381      247272 :   GEN ev = FlxY_evalx(b, n, p);
    1382      247371 :   long drop = lg(b) - lg(ev);
    1383      247371 :   ulong r = Flx_resultant(a, ev, p);
    1384      247277 :   if (drop && la != 1) r = Fl_mul(r, Fl_powu(la, drop,p),p);
    1385      247279 :   return r;
    1386             : }
    1387             : static GEN
    1388           4 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
    1389             : {
    1390           4 :   GEN ev = FpXY_evaly(b, n, p, vX);
    1391           4 :   long drop = db-degpol(ev);
    1392           4 :   GEN r = FpX_resultant(a, ev, p);
    1393           4 :   if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
    1394           4 :   return r;
    1395             : }
    1396             : 
    1397             : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1398             : /* Return a Fly */
    1399             : static GEN
    1400       12021 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, long dres, long sx)
    1401             : {
    1402             :   long i;
    1403       12021 :   ulong n, la = Flx_lead(a);
    1404       12021 :   GEN  x = cgetg(dres+2, t_VECSMALL);
    1405       12021 :   GEN  y = cgetg(dres+2, t_VECSMALL);
    1406             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1407             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1408      131114 :   for (i=0,n = 1; i < dres; n++)
    1409             :   {
    1410      119095 :     x[++i] = n;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1411      119086 :     x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1412             :   }
    1413       12019 :   if (i == dres)
    1414             :   {
    1415        9360 :     x[++i] = 0;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1416             :   }
    1417       12020 :   return Flv_polint(x,y, p, sx);
    1418             : }
    1419             : 
    1420             : static GEN
    1421        6024 : FlxX_pseudorem(GEN x, GEN y, ulong p)
    1422             : {
    1423        6024 :   long vx = varn(x), dx, dy, dz, i, lx, dp;
    1424        6024 :   pari_sp av = avma, av2;
    1425             : 
    1426        6024 :   if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
    1427        6024 :   (void)new_chunk(2);
    1428        6025 :   dx=degpol(x); x = RgX_recip_shallow(x)+2;
    1429        6032 :   dy=degpol(y); y = RgX_recip_shallow(y)+2; dz=dx-dy; dp = dz+1;
    1430        6031 :   av2 = avma;
    1431             :   for (;;)
    1432             :   {
    1433       80367 :     gel(x,0) = Flx_neg(gel(x,0), p); dp--;
    1434      168715 :     for (i=1; i<=dy; i++)
    1435      246786 :       gel(x,i) = Flx_add( Flx_mul(gel(y,0), gel(x,i), p),
    1436      123393 :                               Flx_mul(gel(x,0), gel(y,i), p), p );
    1437      664712 :     for (   ; i<=dx; i++)
    1438      621516 :       gel(x,i) = Flx_mul(gel(y,0), gel(x,i), p);
    1439       45639 :     do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
    1440       43196 :     if (dx < dy) break;
    1441       37171 :     if (gc_needed(av2,1))
    1442             :     {
    1443           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
    1444           0 :       gerepilecoeffs(av2,x,dx+1);
    1445             :     }
    1446             :   }
    1447        6025 :   if (dx < 0) return zero_Flx(0);
    1448        6025 :   lx = dx+3; x -= 2;
    1449        6025 :   x[0]=evaltyp(t_POL) | evallg(lx);
    1450        6026 :   x[1]=evalsigne(1) | evalvarn(vx);
    1451        6026 :   x = RgX_recip_shallow(x);
    1452        6033 :   if (dp)
    1453             :   { /* multiply by y[0]^dp   [beware dummy vars from FpX_FpXY_resultant] */
    1454        1262 :     GEN t = Flx_powu(gel(y,0), dp, p);
    1455        5048 :     for (i=2; i<lx; i++)
    1456        3785 :       gel(x,i) = Flx_mul(gel(x,i), t, p);
    1457             :   }
    1458        6034 :   return gerepilecopy(av, x);
    1459             : }
    1460             : 
    1461             : /* return a Flx */
    1462             : GEN
    1463        1955 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
    1464             : {
    1465        1955 :   pari_sp av = avma, av2;
    1466             :   long degq,dx,dy,du,dv,dr,signh;
    1467             :   GEN z,g,h,r,p1;
    1468             : 
    1469        1955 :   dx=degpol(u); dy=degpol(v); signh=1;
    1470        1956 :   if (dx < dy)
    1471             :   {
    1472           0 :     swap(u,v); lswap(dx,dy);
    1473           0 :     if (both_odd(dx, dy)) signh = -signh;
    1474             :   }
    1475        1956 :   if (dy < 0) return zero_Flx(sx);
    1476        1956 :   if (dy==0) return gerepileupto(av, Flx_powu(gel(v,2),dx,p));
    1477             : 
    1478        1956 :   g = h = pol1_Flx(sx); av2 = avma;
    1479             :   for(;;)
    1480             :   {
    1481       10089 :     r = FlxX_pseudorem(u,v,p); dr = lg(r);
    1482        6036 :     if (dr == 2) { set_avma(av); return zero_Flx(sx); }
    1483        6036 :     du = degpol(u); dv = degpol(v); degq = du-dv;
    1484        6033 :     u = v; p1 = g; g = leading_coeff(u);
    1485        6032 :     switch(degq)
    1486             :     {
    1487           0 :       case 0: break;
    1488             :       case 1:
    1489        4471 :         p1 = Flx_mul(h,p1, p); h = g; break;
    1490             :       default:
    1491        1561 :         p1 = Flx_mul(Flx_powu(h,degq,p), p1, p);
    1492        1561 :         h = Flx_div(Flx_powu(g,degq,p), Flx_powu(h,degq-1,p), p);
    1493             :     }
    1494        6019 :     if (both_odd(du,dv)) signh = -signh;
    1495        6015 :     v = FlxY_Flx_div(r, p1, p);
    1496        6017 :     if (dr==3) break;
    1497        4066 :     if (gc_needed(av2,1))
    1498             :     {
    1499           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_resultant, dr = %ld",dr);
    1500           0 :       gerepileall(av2,4, &u, &v, &g, &h);
    1501             :     }
    1502             :   }
    1503        1951 :   z = gel(v,2);
    1504        1951 :   if (dv > 1) z = Flx_div(Flx_powu(z,dv,p), Flx_powu(h,dv-1,p), p);
    1505        1951 :   if (signh < 0) z = Flx_neg(z,p);
    1506        1951 :   return gerepileupto(av, z);
    1507             : }
    1508             : 
    1509             : /* Warning:
    1510             :  * This function switches between valid and invalid variable ordering*/
    1511             : 
    1512             : static GEN
    1513        1969 : FlxY_to_FlyX(GEN b, long sv)
    1514             : {
    1515        1969 :   long i, n=-1;
    1516        1969 :   long sw = b[1]&VARNBITS;
    1517        1969 :   for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
    1518        1970 :   return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
    1519             : }
    1520             : 
    1521             : /* Return a Fly*/
    1522             : GEN
    1523        1968 : Flx_FlxY_resultant(GEN a, GEN b, ulong pp)
    1524             : {
    1525        1968 :   pari_sp ltop=avma;
    1526        1968 :   long dres = degpol(a)*degpol(b);
    1527        1968 :   long sx=a[1], sy=b[1]&VARNBITS;
    1528             :   GEN z;
    1529        1968 :   b = FlxY_to_FlyX(b,sx);
    1530        1967 :   if ((ulong)dres >= pp)
    1531        1954 :     z = FlxX_resultant(Fly_to_FlxY(a, sy), b, pp, sx);
    1532             :   else
    1533          13 :     z = Flx_FlxY_resultant_polint(a, b, pp, (ulong)dres, sy);
    1534        1971 :   return gerepileupto(ltop,z);
    1535             : }
    1536             : 
    1537             : /* return a t_POL (in variable v >= 0) whose coeffs are the coeffs of b,
    1538             :  * in variable v. This is an incorrect PARI object if initially varn(b) << v.
    1539             :  * We could return a vector of coeffs, but it is convenient to have degpol()
    1540             :  * and friends available. Even in that case, it will behave nicely with all
    1541             :  * functions treating a polynomial as a vector of coeffs (eg poleval).
    1542             :  * FOR INTERNAL USE! */
    1543             : GEN
    1544        8575 : swap_vars(GEN b0, long v)
    1545             : {
    1546        8575 :   long i, n = RgX_degree(b0, v);
    1547             :   GEN b, x;
    1548        8575 :   if (n < 0) return pol_0(v);
    1549        8575 :   b = cgetg(n+3, t_POL); x = b + 2;
    1550        8575 :   b[1] = evalsigne(1) | evalvarn(v);
    1551        8575 :   for (i=0; i<=n; i++) gel(x,i) = polcoef_i(b0, i, v);
    1552        8575 :   return b;
    1553             : }
    1554             : 
    1555             : /* assume varn(b) << varn(a) */
    1556             : /* return a FpY*/
    1557             : GEN
    1558        1942 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
    1559             : {
    1560        1942 :   long i,n,dres, db, vY = varn(b), vX = varn(a);
    1561             :   GEN la,x,y;
    1562             : 
    1563        1942 :   if (lgefint(p) == 3)
    1564             :   {
    1565        1941 :     ulong pp = uel(p,2);
    1566        1941 :     b = ZXX_to_FlxX(b, pp, vX);
    1567        1939 :     a = ZX_to_Flx(a, pp);
    1568        1941 :     x = Flx_FlxY_resultant(a, b, pp);
    1569        1945 :     return Flx_to_ZX(x);
    1570             :   }
    1571           1 :   db = RgXY_degreex(b);
    1572           1 :   dres = degpol(a)*degpol(b);
    1573           1 :   la = leading_coeff(a);
    1574           1 :   x = cgetg(dres+2, t_VEC);
    1575           1 :   y = cgetg(dres+2, t_VEC);
    1576             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1577             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1578           3 :   for (i=0,n = 1; i < dres; n++)
    1579             :   {
    1580           2 :     gel(x,++i) = utoipos(n);
    1581           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1582           2 :     gel(x,++i) = subiu(p,n);
    1583           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1584             :   }
    1585           1 :   if (i == dres)
    1586             :   {
    1587           0 :     gel(x,++i) = gen_0;
    1588           0 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
    1589             :   }
    1590           1 :   return FpV_polint(x,y, p, vY);
    1591             : }
    1592             : 
    1593             : static GEN
    1594         182 : FpX_diamondsum(GEN P, GEN Q, GEN p)
    1595             : {
    1596         182 :   long n = 1+ degpol(P)*degpol(Q);
    1597         182 :   GEN Pl = FpX_invLaplace(FpX_Newton(P,n,p), p);
    1598         182 :   GEN Ql = FpX_invLaplace(FpX_Newton(Q,n,p), p);
    1599         182 :   GEN L = FpX_Laplace(FpXn_mul(Pl, Ql, n, p), p);
    1600         182 :   return FpX_fromNewton(L, p);
    1601             : }
    1602             : 
    1603             : #if 0
    1604             : GEN
    1605             : FpX_diamondprod(GEN P, GEN Q, GEN p)
    1606             : {
    1607             :   long n = 1+ degpol(P)*degpol(Q);
    1608             :   GEN L=FpX_convol(FpX_Newton(P,n,p), FpX_Newton(Q,n,p), p);
    1609             :   return FpX_fromNewton(L, p);
    1610             : }
    1611             : #endif
    1612             : 
    1613             : GEN
    1614         588 : FpX_direct_compositum(GEN a, GEN b, GEN p)
    1615             : {
    1616         588 :   long da = degpol(a), db = degpol(b);
    1617         588 :   if (cmpis(p, da*db) > 0)
    1618         182 :     return FpX_diamondsum(a, b, p);
    1619             :   else
    1620             :   {
    1621         406 :     long v = varn(a), w = fetch_var_higher();
    1622         406 :     GEN mx = deg1pol_shallow(gen_m1, gen_0, v);
    1623         406 :     GEN r, ymx = deg1pol_shallow(gen_1, mx, w); /* Y-X */
    1624         406 :     if (degpol(a) < degpol(b)) swap(a,b);
    1625         406 :     r = FpX_FpXY_resultant(a, poleval(b,ymx),p);
    1626         406 :     setvarn(r, v); (void)delete_var(); return r;
    1627             :   }
    1628             : }
    1629             : 
    1630             : static GEN
    1631         588 : _FpX_direct_compositum(void *E, GEN a, GEN b)
    1632         588 : { return FpX_direct_compositum(a,b, (GEN)E); }
    1633             : 
    1634             : GEN
    1635        5261 : FpXV_direct_compositum(GEN V, GEN p)
    1636             : {
    1637        5261 :   return gen_product(V, (void *)p, &_FpX_direct_compositum);
    1638             : }
    1639             : 
    1640             : /* 0, 1, -1, 2, -2, ... */
    1641             : #define next_lambda(a) (a>0 ? -a : 1-a)
    1642             : GEN
    1643           0 : FpX_compositum(GEN a, GEN b, GEN p)
    1644             : {
    1645           0 :   long k, v = fetch_var_higher();
    1646           0 :   for (k = 1;; k = next_lambda(k))
    1647           0 :   {
    1648           0 :     GEN x = deg1pol_shallow(gen_1, gmulsg(k, pol_x(v)), 0); /* x + k y */
    1649           0 :     GEN C = FpX_FpXY_resultant(a, poleval(b,x),p);
    1650           0 :     if (FpX_is_squarefree(C, p)) { (void)delete_var(); return C; }
    1651             :   }
    1652             : }
    1653             : 
    1654             : /* Assume A in Z[Y], B in Q[Y][X], B squarefree in (Q[Y]/(A))[X] and
    1655             :  * Res_Y(A, B) in Z[X]. Find a small lambda (start from *lambda, use
    1656             :  * next_lambda successively) such that C(X) = Res_Y(A(Y), B(X + lambda Y))
    1657             :  * is squarefree, reset *lambda to the chosen value and return C. Set LERS to
    1658             :  * the Last non-constant polynomial in the Euclidean Remainder Sequence */
    1659             : static GEN
    1660        1995 : ZX_ZXY_resultant_LERS(GEN A, GEN B0, long *plambda, GEN *LERS)
    1661             : {
    1662             :   ulong bound, dp;
    1663        1995 :   pari_sp av = avma, av2 = 0;
    1664        1995 :   long lambda = *plambda, degA = degpol(A), dres = degA*degpol(B0);
    1665             :   long stable, checksqfree, i,n, cnt, degB;
    1666        1995 :   long v, vX = varn(B0), vY = varn(A); /* vY < vX */
    1667             :   GEN x, y, dglist, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
    1668             :   forprime_t S;
    1669             : 
    1670        1995 :   if (degA == 1)
    1671             :   {
    1672         504 :     GEN a1 = gel(A,3), a0 = gel(A,2);
    1673         504 :     B = lambda? RgX_translate(B0, monomial(stoi(lambda), 1, vY)): B0;
    1674         504 :     H = gsubst(B, vY, gdiv(gneg(a0),a1));
    1675         504 :    if (!equali1(a1)) H = RgX_Rg_mul(H, powiu(a1, poldegree(B,vY)));
    1676         504 :     *LERS = mkvec2(scalarpol_shallow(a0,vX), scalarpol_shallow(a1,vX));
    1677         504 :     gerepileall(av, 2, &H, LERS);
    1678         504 :     return H;
    1679             :   }
    1680             : 
    1681        1491 :   dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
    1682        1491 :   C0 = cgetg(dres+2, t_VECSMALL);
    1683        1491 :   C1 = cgetg(dres+2, t_VECSMALL);
    1684        1491 :   dglist = cgetg(dres+1, t_VECSMALL);
    1685        1491 :   x = cgetg(dres+2, t_VECSMALL);
    1686        1491 :   y = cgetg(dres+2, t_VECSMALL);
    1687        1491 :   B0 = leafcopy(B0);
    1688        1491 :   A = leafcopy(A);
    1689        1491 :   B = B0;
    1690        1491 :   v = fetch_var_higher(); setvarn(A,v);
    1691             :   /* make sure p large enough */
    1692             : INIT:
    1693             :   /* always except the first time */
    1694        1911 :   if (av2) { set_avma(av2); lambda = next_lambda(lambda); }
    1695        1911 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    1696        1911 :   B = swap_vars(B, vY); setvarn(B,v);
    1697             :   /* B0(lambda v + x, v) */
    1698        1911 :   if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    1699        1911 :   av2 = avma;
    1700             : 
    1701        1911 :   if (degA <= 3)
    1702             :   { /* sub-resultant faster for small degrees */
    1703        1659 :     H = RgX_resultant_all(A,B,&q);
    1704        1659 :     if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
    1705        1295 :     H0 = gel(q,2);
    1706        1295 :     if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
    1707        1295 :     H1 = gel(q,3);
    1708        1295 :     if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
    1709        1295 :     if (!ZX_is_squarefree(H)) goto INIT;
    1710        1253 :     goto END;
    1711             :   }
    1712             : 
    1713         252 :   H = H0 = H1 = NULL;
    1714         252 :   degB = degpol(B);
    1715         252 :   bound = ZX_ZXY_ResBound(A, B, NULL);
    1716         252 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    1717         252 :   dp = 1;
    1718         252 :   init_modular_big(&S);
    1719         252 :   for(cnt = 0, checksqfree = 1;;)
    1720         273 :   {
    1721         525 :     ulong p = u_forprime_next(&S);
    1722             :     GEN Hi;
    1723         525 :     a = ZX_to_Flx(A, p);
    1724         525 :     b = ZXX_to_FlxX(B, p, varn(A));
    1725         525 :     if (degpol(a) < degA || degpol(b) < degB) continue; /* p | lc(A)lc(B) */
    1726         525 :     if (checksqfree)
    1727             :     { /* find degree list for generic Euclidean Remainder Sequence */
    1728         252 :       long goal = minss(degpol(a), degpol(b)); /* longest possible */
    1729         252 :       for (n=1; n <= goal; n++) dglist[n] = 0;
    1730         252 :       setlg(dglist, 1);
    1731        1813 :       for (n=0; n <= dres; n++)
    1732             :       {
    1733        1722 :         ev = FlxY_evalx_drop(b, n, p);
    1734        1722 :         Flx_resultant_set_dglist(a, ev, dglist, p);
    1735        1722 :         if (lg(dglist)-1 == goal) break;
    1736             :       }
    1737             :       /* last pol in ERS has degree > 1 ? */
    1738         252 :       goal = lg(dglist)-1;
    1739         252 :       if (degpol(B) == 1) { if (!goal) goto INIT; }
    1740             :       else
    1741             :       {
    1742         245 :         if (goal <= 1) goto INIT;
    1743         231 :         if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
    1744             :       }
    1745         238 :       if (DEBUGLEVEL>4)
    1746           0 :         err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
    1747             :     }
    1748             : 
    1749        9218 :     for (i=0,n = 0; i <= dres; n++)
    1750             :     {
    1751        8707 :       ev = FlxY_evalx_drop(b, n, p);
    1752        8707 :       x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
    1753        8707 :       if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
    1754             :     }
    1755         511 :     Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
    1756         511 :     Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
    1757         511 :     if (!H && degpol(Hp) != dres) continue;
    1758         511 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1759         511 :     if (checksqfree) {
    1760         238 :       if (!Flx_is_squarefree(Hp, p)) goto INIT;
    1761         238 :       if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    1762         238 :       checksqfree = 0;
    1763             :     }
    1764             : 
    1765         511 :     if (!H)
    1766             :     { /* initialize */
    1767         238 :       q = utoipos(p); stable = 0;
    1768         238 :       H = ZX_init_CRT(Hp, p,vX);
    1769         238 :       H0= ZX_init_CRT(H0p, p,vX);
    1770         238 :       H1= ZX_init_CRT(H1p, p,vX);
    1771             :     }
    1772             :     else
    1773             :     {
    1774         273 :       GEN qp = muliu(q,p);
    1775         546 :       stable  = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
    1776         273 :               & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
    1777         273 :               & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
    1778         273 :       q = qp;
    1779             :     }
    1780             :     /* could make it probabilistic for H ? [e.g if stable twice, etc]
    1781             :      * Probabilistic anyway for H0, H1 */
    1782         511 :     if (DEBUGLEVEL>5 && (stable ||  ++cnt==100))
    1783           0 :     { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
    1784         511 :     if (stable && (ulong)expi(q) >= bound) break; /* DONE */
    1785         273 :     if (gc_needed(av,2))
    1786             :     {
    1787           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
    1788           0 :       gerepileall(av2, 4, &H, &q, &H0, &H1);
    1789             :     }
    1790             :   }
    1791             : END:
    1792        1491 :   if (DEBUGLEVEL>5) err_printf(" done\n");
    1793        1491 :   setvarn(H, vX); (void)delete_var();
    1794        1491 :   *LERS = mkvec2(H0,H1);
    1795        1491 :   gerepileall(av, 2, &H, LERS);
    1796        1491 :   *plambda = lambda; return H;
    1797             : }
    1798             : 
    1799             : GEN
    1800        2331 : ZX_ZXY_resultant_all(GEN A, GEN B, long *plambda, GEN *LERS)
    1801             : {
    1802        2331 :   if (LERS)
    1803             :   {
    1804        1995 :     if (!plambda)
    1805           0 :       pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
    1806        1995 :     return ZX_ZXY_resultant_LERS(A, B, plambda, LERS);
    1807             :   }
    1808         336 :   return ZX_ZXY_rnfequation(A, B, plambda);
    1809             : }
    1810             : 
    1811             : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
    1812             :  * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
    1813             :  * squarefree */
    1814             : GEN
    1815        1848 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
    1816             : {
    1817        1848 :   pari_sp av = avma;
    1818             :   GEN R, a;
    1819             :   long dA;
    1820             :   int delvar;
    1821             : 
    1822        1848 :   if (v < 0) v = 0;
    1823        1848 :   switch (typ(A))
    1824             :   {
    1825        1848 :     case t_POL: dA = degpol(A); if (dA > 0) break;
    1826           0 :       A = constant_coeff(A);
    1827             :     default:
    1828           0 :       if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
    1829           0 :       return gerepileupto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
    1830             :   }
    1831        1848 :   delvar = 0;
    1832        1848 :   if (varn(T) == 0)
    1833             :   {
    1834        1764 :     long v0 = fetch_var(); delvar = 1;
    1835        1764 :     T = leafcopy(T); setvarn(T,v0);
    1836        1764 :     A = leafcopy(A); setvarn(A,v0);
    1837             :   }
    1838        1848 :   R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
    1839        1848 :   if (delvar) (void)delete_var();
    1840        1848 :   setvarn(R, v); a = leading_coeff(T);
    1841        1848 :   if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
    1842        1848 :   return gerepileupto(av, R);
    1843             : }
    1844             : 
    1845             : 
    1846             : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
    1847             : GEN
    1848       12247 : ZXQ_charpoly(GEN A, GEN T, long v)
    1849             : {
    1850       12247 :   return (degpol(T) < 16) ? RgXQ_charpoly(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
    1851             : }
    1852             : 
    1853             : GEN
    1854         819 : QXQ_charpoly(GEN A, GEN T, long v)
    1855             : {
    1856         819 :   pari_sp av = avma;
    1857         819 :   GEN den, B = Q_remove_denom(A, &den);
    1858         819 :   GEN P = ZXQ_charpoly(B, T, v);
    1859         819 :   return gerepilecopy(av, den ? RgX_rescale(P, ginv(den)): P);
    1860             : }
    1861             : 
    1862             : static GEN
    1863      164290 : trivial_case(GEN A, GEN B)
    1864             : {
    1865             :   long d;
    1866      164290 :   if (typ(A) == t_INT) return powiu(A, degpol(B));
    1867      156540 :   d = degpol(A);
    1868      156540 :   if (d == 0) return trivial_case(gel(A,2),B);
    1869      153508 :   if (d < 0) return gen_0;
    1870      153486 :   return NULL;
    1871             : }
    1872             : 
    1873             : static ulong
    1874     1324067 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
    1875             : {
    1876     1324067 :   pari_sp av = avma;
    1877             :   ulong H;
    1878             :   long dropa, dropb;
    1879     1324067 :   ulong dp = dB ? umodiu(dB, p): 1;
    1880     1324067 :   if (!b) b = Flx_deriv(a, p);
    1881     1324167 :   dropa = degA - degpol(a);
    1882     1324152 :   dropb = degB - degpol(b);
    1883     1324109 :   if (dropa && dropb) return gc_ulong(av,0); /* p | lc(A), p | lc(B) */
    1884     1324109 :   H = Flx_resultant(a, b, p);
    1885     1323866 :   if (dropa)
    1886             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    1887           0 :     ulong c = b[degB+2]; /* lc(B) */
    1888           0 :     if (odd(degB)) c = p - c;
    1889           0 :     c = Fl_powu(c, dropa, p);
    1890           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1891             :   }
    1892     1323866 :   else if (dropb)
    1893             :   { /* multiply by lc(A)^(deg B - deg b) */
    1894           0 :     ulong c = a[degA+2]; /* lc(A) */
    1895           0 :     c = Fl_powu(c, dropb, p);
    1896           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1897             :   }
    1898     1323864 :   if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1899     1323865 :   return gc_ulong(av,H);
    1900             : }
    1901             : 
    1902             : /* If B=NULL, assume B=A' */
    1903             : static GEN
    1904      552748 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
    1905             : {
    1906      552748 :   pari_sp av = avma;
    1907      552748 :   long degA, degB, i, n = lg(P)-1;
    1908             :   GEN H, T;
    1909             : 
    1910      552748 :   degA = degpol(A);
    1911      552747 :   degB = B ? degpol(B): degA - 1;
    1912      552747 :   if (n == 1)
    1913             :   {
    1914      164008 :     ulong Hp, p = uel(P,1);
    1915             :     GEN a, b;
    1916      164008 :     a = ZX_to_Flx(A, p), b = B ? ZX_to_Flx(B, p): NULL;
    1917      164011 :     Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1918      164003 :     set_avma(av);
    1919      164002 :     *mod = utoi(p); return utoi(Hp);
    1920             :   }
    1921      388739 :   T = ZV_producttree(P);
    1922      388735 :   A = ZX_nv_mod_tree(A, P, T);
    1923      388739 :   if (B) B = ZX_nv_mod_tree(B, P, T);
    1924      388739 :   H = cgetg(n+1, t_VECSMALL);
    1925     1548689 :   for(i=1; i <= n; i++)
    1926             :   {
    1927     1160079 :     ulong p = P[i];
    1928     1160079 :     GEN a = gel(A,i), b = B? gel(B,i): NULL;
    1929     1160079 :     H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1930             :   }
    1931      388610 :   H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
    1932      388692 :   *mod = gmael(T, lg(T)-1, 1);
    1933      388692 :   gerepileall(av, 2, &H, mod);
    1934      388716 :   return H;
    1935             : }
    1936             : 
    1937             : GEN
    1938      552748 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
    1939             : {
    1940      552748 :   GEN V = cgetg(3, t_VEC);
    1941      552753 :   if (isintzero(B)) B = NULL;
    1942      552748 :   if (isintzero(dB)) dB = NULL;
    1943      552747 :   gel(V,1) = ZX_resultant_slice(A,B,dB,P,&gel(V,2));
    1944      552677 :   return V;
    1945             : }
    1946             : 
    1947             : /* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */
    1948             : /* if B=NULL, take B = A' */
    1949             : GEN
    1950       83933 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
    1951             : {
    1952       83933 :   pari_sp av = avma;
    1953             :   long m;
    1954             :   GEN  H, worker;
    1955       83933 :   int is_disc = !B;
    1956       83933 :   if (is_disc) B = ZX_deriv(A);
    1957       83933 :   if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H;
    1958       76161 :   if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
    1959       76161 :   if (is_disc)
    1960       47826 :     B = NULL;
    1961       76161 :   worker = strtoclosure("_ZX_resultant_worker", 3, A, B?B:gen_0, dB?dB:gen_0);
    1962       76161 :   m = degpol(A)+(B ? degpol(B): 0);
    1963       76161 :   H = gen_crt("ZX_resultant_all", worker, dB, bound, m, NULL,
    1964             :                ZV_chinese_center, Fp_center);
    1965       76161 :   return gerepileuptoint(av, H);
    1966             : }
    1967             : 
    1968             : /* A0 and B0 in Q[X] */
    1969             : GEN
    1970       13035 : QX_resultant(GEN A0, GEN B0)
    1971             : {
    1972             :   GEN s, a, b, A, B;
    1973       13035 :   pari_sp av = avma;
    1974             : 
    1975       13035 :   A = Q_primitive_part(A0, &a);
    1976       13035 :   B = Q_primitive_part(B0, &b);
    1977       13035 :   s = ZX_resultant(A, B);
    1978       13035 :   if (!signe(s)) { set_avma(av); return gen_0; }
    1979       13035 :   if (a) s = gmul(s, gpowgs(a,degpol(B)));
    1980       13035 :   if (b) s = gmul(s, gpowgs(b,degpol(A)));
    1981       13035 :   return gerepileupto(av, s);
    1982             : }
    1983             : 
    1984             : GEN
    1985       35441 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
    1986             : 
    1987             : GEN
    1988           0 : QXQ_intnorm(GEN A, GEN B)
    1989             : {
    1990             :   GEN c, n, R, lB;
    1991           0 :   long dA = degpol(A), dB = degpol(B);
    1992           0 :   pari_sp av = avma;
    1993           0 :   if (dA < 0) return gen_0;
    1994           0 :   A = Q_primitive_part(A, &c);
    1995           0 :   if (!c || typ(c) == t_INT) {
    1996           0 :     n = c;
    1997           0 :     R = ZX_resultant(B, A);
    1998             :   } else {
    1999           0 :     n = gel(c,1);
    2000           0 :     R = ZX_resultant_all(B, A, gel(c,2), 0);
    2001             :   }
    2002           0 :   if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
    2003           0 :   lB = leading_coeff(B);
    2004           0 :   if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
    2005           0 :   return gerepileuptoint(av, R);
    2006             : }
    2007             : 
    2008             : GEN
    2009           0 : QXQ_norm(GEN A, GEN B)
    2010             : {
    2011             :   GEN c, R, lB;
    2012           0 :   long dA = degpol(A), dB = degpol(B);
    2013           0 :   pari_sp av = avma;
    2014           0 :   if (dA < 0) return gen_0;
    2015           0 :   A = Q_primitive_part(A, &c);
    2016           0 :   R = ZX_resultant(B, A);
    2017           0 :   if (c) R = gmul(R, gpowgs(c, dB));
    2018           0 :   lB = leading_coeff(B);
    2019           0 :   if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
    2020           0 :   return gerepileupto(av, R);
    2021             : }
    2022             : 
    2023             : /* assume x has integral coefficients */
    2024             : GEN
    2025       49233 : ZX_disc_all(GEN x, ulong bound)
    2026             : {
    2027       49233 :   pari_sp av = avma;
    2028             :   GEN l, R;
    2029       49233 :   long s, d = degpol(x);
    2030       49233 :   if (d <= 1) return d ? gen_1: gen_0;
    2031       47826 :   s = (d & 2) ? -1: 1;
    2032       47826 :   l = leading_coeff(x);
    2033       47826 :   R = ZX_resultant_all(x, NULL, NULL, bound);
    2034       47826 :   if (is_pm1(l))
    2035       44977 :   { if (signe(l) < 0) s = -s; }
    2036             :   else
    2037        2849 :     R = diviiexact(R,l);
    2038       47826 :   if (s == -1) togglesign_safe(&R);
    2039       47826 :   return gerepileuptoint(av,R);
    2040             : }
    2041       48176 : GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); }
    2042             : 
    2043             : GEN
    2044           0 : QX_disc(GEN x)
    2045             : {
    2046           0 :   pari_sp av = avma;
    2047           0 :   GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
    2048           0 :   if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
    2049           0 :   return gerepileupto(av, d);
    2050             : }
    2051             : 
    2052             : GEN
    2053       43321 : QXQ_mul(GEN x, GEN y, GEN T)
    2054             : {
    2055       43321 :   GEN dx, nx = Q_primitive_part(x, &dx);
    2056       43321 :   GEN dy, ny = Q_primitive_part(y, &dy);
    2057       43321 :   GEN z = ZXQ_mul(nx, ny, T);
    2058       43321 :   if (dx || dy)
    2059             :   {
    2060       43321 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
    2061       43321 :     if (!gequal1(d)) z = ZX_Q_mul(z, d);
    2062             :   }
    2063       43321 :   return z;
    2064             : }
    2065             : 
    2066             : GEN
    2067       11431 : QXQ_sqr(GEN x, GEN T)
    2068             : {
    2069       11431 :   GEN dx, nx = Q_primitive_part(x, &dx);
    2070       11431 :   GEN z = ZXQ_sqr(nx, T);
    2071       11431 :   if (dx)
    2072       11431 :     z = ZX_Q_mul(z, gsqr(dx));
    2073       11431 :   return z;
    2074             : }
    2075             : 
    2076             : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
    2077             : GEN
    2078       29004 : QXQ_inv(GEN A, GEN B)
    2079             : {
    2080             :   GEN D, cU, q, U, V;
    2081             :   ulong p;
    2082       29004 :   pari_sp av2, av = avma;
    2083             :   forprime_t S;
    2084             :   pari_timer ti;
    2085       29004 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2086             :   /* A a QX, B a ZX */
    2087       29004 :   A = Q_primitive_part(A, &D);
    2088             :   /* A, B in Z[X] */
    2089       29004 :   init_modular_small(&S);
    2090       29004 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2091       29004 :   av2 = avma; U = NULL;
    2092      161240 :   while ((p = u_forprime_next(&S)))
    2093             :   {
    2094             :     GEN a, b, qp, Up, Vp;
    2095             :     int stable;
    2096             : 
    2097      132236 :     a = ZX_to_Flx(A, p);
    2098      132236 :     b = ZX_to_Flx(B, p);
    2099             :     /* if p | Res(A/G, B/G), discard */
    2100      161226 :     if (!Flx_extresultant(b,a,p, &Vp,&Up)) continue;
    2101             : 
    2102      132222 :     if (!U)
    2103             :     { /* First time */
    2104       28990 :       U = ZX_init_CRT(Up,p,varn(A));
    2105       28990 :       V = ZX_init_CRT(Vp,p,varn(A));
    2106       28990 :       q = utoipos(p); continue;
    2107             :     }
    2108      103232 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: mod %ld (bound 2^%ld)", p,expi(q));
    2109      103232 :     qp = muliu(q,p);
    2110      206464 :     stable = ZX_incremental_CRT_raw(&U, Up, q,qp, p)
    2111      103232 :            & ZX_incremental_CRT_raw(&V, Vp, q,qp, p);
    2112      103232 :     if (stable)
    2113             :     { /* all stable: check divisibility */
    2114       28990 :       GEN res = ZX_add(ZX_mul(A,U), ZX_mul(B,V));
    2115       28990 :       if (degpol(res) == 0) {
    2116       28990 :         res = gel(res,2);
    2117       28990 :         D = D? gmul(D, res): res;
    2118       57980 :         break;
    2119             :       } /* DONE */
    2120           0 :       if (DEBUGLEVEL) err_printf("QXQ_inv: char 0 check failed");
    2121             :     }
    2122       74242 :     q = qp;
    2123       74242 :     if (gc_needed(av,1))
    2124             :     {
    2125          33 :       if (DEBUGMEM>1) pari_warn(warnmem,"QXQ_inv");
    2126          33 :       gerepileall(av2, 3, &q,&U,&V);
    2127             :     }
    2128             :   }
    2129       28990 :   if (!p) pari_err_OVERFLOW("QXQ_inv [ran out of primes]");
    2130       28990 :   cU = ZX_content(U);
    2131       28990 :   if (!is_pm1(cU)) { U = Q_div_to_int(U, cU); D = gdiv(D, cU); }
    2132       28990 :   return gerepileupto(av, RgX_Rg_div(U, D));
    2133             : }
    2134             : 
    2135             : /* lift(C / Mod(A,B)). B monic ZX, A and C scalar or QX. Use when result is
    2136             :  * small */
    2137             : GEN
    2138         273 : QXQ_div_ratlift(GEN C, GEN A, GEN B)
    2139             : {
    2140             :   GEN dA, dC, q, U;
    2141             :   ulong p, ct, delay;
    2142         273 :   pari_sp av2, av = avma;
    2143             :   forprime_t S;
    2144             :   pari_timer ti;
    2145         273 :   if (is_scalar_t(typ(A)))
    2146             :   {
    2147           0 :     A = gdiv(C,A);
    2148           0 :     if (typ(A) != t_POL) A = scalarpol(A, varn(B));
    2149           0 :     return A;
    2150             :   }
    2151             :   /* A a QX, B a ZX */
    2152         273 :   A = Q_remove_denom(A, &dA);
    2153         273 :   C = Q_remove_denom(C, &dC);
    2154         273 :   if (typ(C) != t_POL) C = scalarpol_shallow(C, varn(B));
    2155         273 :   if (dA) C = ZX_Z_mul(C,dA);
    2156             :   /* A, B, C in Z[X] */
    2157         273 :   init_modular_small(&S);
    2158         273 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2159         273 :   av2 = avma; U = NULL; ct = 0; delay = 1;
    2160        1938 :   while ((p = u_forprime_next(&S)))
    2161             :   {
    2162             :     GEN a, b, Up, Ur;
    2163        1665 :     a = ZX_to_Flx(A, p);
    2164        1665 :     b = ZX_to_Flx(B, p);
    2165             :     /* if p | Res(A/G, B/G), discard */
    2166        1665 :     Up = Flxq_invsafe(a,b,p); if (!Up) continue;
    2167        1665 :     Up = Flxq_mul(Up, ZX_to_Flx(C,p), b, p);
    2168             : 
    2169        1665 :     if (!U)
    2170             :     { /* First time */
    2171         273 :       U = ZX_init_CRT(Up,p,varn(A));
    2172         273 :       q = utoipos(p);
    2173             :     }
    2174             :     else
    2175             :     {
    2176        1392 :       GEN qp = muliu(q,p);
    2177        1392 :       (void)ZX_incremental_CRT_raw(&U, Up, q,qp, p);
    2178        1392 :       q = qp;
    2179             :     }
    2180        1665 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: mod %ld (bound 2^%ld)", p,expi(q));
    2181        1665 :     b = sqrti(shifti(q,-1));
    2182        1665 :     Ur = FpX_ratlift(U,q,b,b,NULL);
    2183        1665 :     if (Ur && ++ct == delay)
    2184             :     { /* check divisibility */
    2185         287 :       GEN d, V = Q_remove_denom(Ur,&d), W = d? ZX_Z_mul(C,d): C;
    2186         287 :       if (!signe(ZX_rem(ZX_sub(ZX_mul(A,V), W), B))) { U = Ur; break; }
    2187          14 :       delay <<= 1;
    2188          14 :       if (DEBUGLEVEL) err_printf("QXQ_div: check failed, delay = %ld",delay);
    2189             :     }
    2190        1392 :     if (gc_needed(av,1))
    2191             :     {
    2192           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"QXQ_div");
    2193           0 :       gerepileall(av2, 2, &q,&U);
    2194             :     }
    2195             :   }
    2196         273 :   if (!p) pari_err_OVERFLOW("QXQ_div [ran out of primes]");
    2197         273 :   if (!dC) return gerepilecopy(av, U);
    2198           0 :   return gerepileupto(av, RgX_Rg_div(U, dC));
    2199             : }
    2200             : 
    2201             : /************************************************************************
    2202             :  *                                                                      *
    2203             :  *                   ZX_ZXY_resultant                                   *
    2204             :  *                                                                      *
    2205             :  ************************************************************************/
    2206             : 
    2207             : static GEN
    2208       12008 : ZX_ZXY_resultant_prime(GEN a, GEN b, ulong dp, ulong p,
    2209             :                        long degA, long degB, long dres, long sX)
    2210             : {
    2211       12008 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    2212       12008 :   GEN Hp = Flx_FlxY_resultant_polint(a, b, p, dres, sX);
    2213       12006 :   if (dropa && dropb)
    2214           0 :     Hp = zero_Flx(sX);
    2215             :   else {
    2216       12006 :     if (dropa)
    2217             :     { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2218           0 :       GEN c = gel(b,degB+2); /* lc(B) */
    2219           0 :       if (odd(degB)) c = Flx_neg(c, p);
    2220           0 :       if (!Flx_equal1(c)) {
    2221           0 :         c = Flx_powu(c, dropa, p);
    2222           0 :         if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
    2223             :       }
    2224             :     }
    2225       12006 :     else if (dropb)
    2226             :     { /* multiply by lc(A)^(deg B - deg b) */
    2227           0 :       ulong c = uel(a, degA+2); /* lc(A) */
    2228           0 :       c = Fl_powu(c, dropb, p);
    2229           0 :       if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
    2230             :     }
    2231             :   }
    2232       12005 :   if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2233       12005 :   return Hp;
    2234             : }
    2235             : 
    2236             : static GEN
    2237        8716 : ZX_ZXY_resultant_slice(GEN A, GEN B, GEN dB, long degA, long degB, long dres,
    2238             :                        GEN P, GEN *mod, long sX, long vY)
    2239             : {
    2240        8716 :   pari_sp av = avma;
    2241        8716 :   long i, n = lg(P)-1;
    2242             :   GEN H, T, D;
    2243        8716 :   if (n == 1)
    2244             :   {
    2245        8311 :     ulong p = uel(P,1);
    2246        8311 :     ulong dp = dB ? umodiu(dB, p): 1;
    2247        8311 :     GEN a = ZX_to_Flx(A, p), b = ZXX_to_FlxX(B, p, vY);
    2248        8311 :     GEN Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2249        8309 :     H = Flx_to_ZX(Hp);
    2250        8309 :     *mod = utoi(p);
    2251        8308 :     gerepileall(av, 2, &H, mod);
    2252        8309 :     return H;
    2253             :   }
    2254         405 :   T = ZV_producttree(P);
    2255         405 :   A = ZX_nv_mod_tree(A, P, T);
    2256         405 :   B = ZXX_nv_mod_tree(B, P, T, vY);
    2257         405 :   D = dB ? Z_ZV_mod_tree(dB, P, T): NULL;
    2258         405 :   H = cgetg(n+1, t_VEC);
    2259        1435 :   for(i=1; i <= n; i++)
    2260             :   {
    2261        1030 :     ulong p = P[i];
    2262        1030 :     GEN a = gel(A,i), b = gel(B,i);
    2263        1030 :     ulong dp = D ? uel(D, i): 1;
    2264        1030 :     gel(H,i) = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2265             :   }
    2266         405 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2267         405 :   *mod = gmael(T, lg(T)-1, 1);
    2268         405 :   gerepileall(av, 2, &H, mod);
    2269         405 :   return H;
    2270             : }
    2271             : 
    2272             : GEN
    2273        8716 : ZX_ZXY_resultant_worker(GEN P, GEN A, GEN B, GEN dB, GEN v)
    2274             : {
    2275        8716 :   GEN V = cgetg(3, t_VEC);
    2276        8716 :   if (isintzero(dB)) dB = NULL;
    2277        8716 :   gel(V,1) = ZX_ZXY_resultant_slice(A, B, dB, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
    2278        8714 :   return V;
    2279             : }
    2280             : 
    2281             : GEN
    2282        3962 : ZX_ZXY_resultant(GEN A, GEN B)
    2283             : {
    2284        3962 :   pari_sp av = avma;
    2285             :   ulong bound;
    2286        3962 :   long v = fetch_var_higher();
    2287        3962 :   long degA = degpol(A), degB, dres = degA * degpol(B);
    2288        3962 :   long vX = varn(B), vY = varn(A); /* assume vY has lower priority */
    2289        3962 :   long sX = evalvarn(vX);
    2290             :   GEN worker, H, dB;
    2291        3962 :   B = Q_remove_denom(B, &dB);
    2292        3962 :   if (!dB) B = leafcopy(B);
    2293        3962 :   A = leafcopy(A); setvarn(A,v);
    2294        3962 :   B = swap_vars(B, vY); setvarn(B,v); degB = degpol(B);
    2295        3962 :   bound = ZX_ZXY_ResBound(A, B, dB);
    2296        3962 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    2297        3962 :   worker = strtoclosure("_ZX_ZXY_resultant_worker", 4, A, B, dB?dB:gen_0,
    2298             :                         mkvecsmall5(degA, degB,dres, vY, sX));
    2299        3962 :   H = gen_crt("ZX_ZXY_resultant_all", worker, dB, bound, degpol(A)+degpol(B), NULL,
    2300             :                nxV_chinese_center, FpX_center_i);
    2301        3962 :   setvarn(H, vX); (void)delete_var();
    2302        3962 :   return gerepilecopy(av, H);
    2303             : }
    2304             : 
    2305             : static long
    2306        2198 : ZX_ZXY_rnfequation_lambda(GEN A, GEN B0, long lambda)
    2307             : {
    2308        2198 :   pari_sp av = avma;
    2309        2198 :   long degA = degpol(A), degB, dres = degA*degpol(B0);
    2310        2198 :   long v = fetch_var_higher();
    2311        2198 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    2312        2198 :   long sX = evalvarn(vX);
    2313             :   GEN dB, B, a, b, Hp;
    2314             :   forprime_t S;
    2315             : 
    2316        2198 :   B0 = Q_remove_denom(B0, &dB);
    2317        2198 :   if (!dB) B0 = leafcopy(B0);
    2318        2198 :   A = leafcopy(A);
    2319        2198 :   B = B0;
    2320        2198 :   setvarn(A,v);
    2321             : INIT:
    2322        2667 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    2323        2667 :   B = swap_vars(B, vY); setvarn(B,v);
    2324             :   /* B0(lambda v + x, v) */
    2325        2667 :   if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    2326             : 
    2327        2667 :   degB = degpol(B);
    2328        2667 :   init_modular_big(&S);
    2329             :   while (1)
    2330           0 :   {
    2331        2667 :     ulong p = u_forprime_next(&S);
    2332        2667 :     ulong dp = dB ? umodiu(dB, p): 1;
    2333        2667 :     if (!dp) continue;
    2334        2667 :     a = ZX_to_Flx(A, p);
    2335        2667 :     b = ZXX_to_FlxX(B, p, v);
    2336        2667 :     Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2337        2667 :     if (degpol(Hp) != dres) continue;
    2338        2667 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2339        2667 :     if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); goto INIT; }
    2340        2198 :     if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    2341        4396 :     (void)delete_var(); return gc_long(av,lambda);
    2342             :   }
    2343             : }
    2344             : 
    2345             : GEN
    2346        2758 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
    2347             : {
    2348        2758 :   if (lambda)
    2349             :   {
    2350        2198 :     *lambda = ZX_ZXY_rnfequation_lambda(A, B, *lambda);
    2351        2198 :     B = RgX_translate(B, monomial(stoi(*lambda), 1, varn(A)));
    2352             :   }
    2353        2758 :   return ZX_ZXY_resultant(A,B);
    2354             : }
    2355             : 
    2356             : /************************************************************************
    2357             :  *                                                                      *
    2358             :  *                   IRREDUCIBLE POLYNOMIAL / Fp                        *
    2359             :  *                                                                      *
    2360             :  ************************************************************************/
    2361             : 
    2362             : /* irreducible (unitary) polynomial of degree n over Fp */
    2363             : GEN
    2364           0 : ffinit_rand(GEN p,long n)
    2365             : {
    2366           0 :   for(;;) {
    2367           0 :     pari_sp av = avma;
    2368           0 :     GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
    2369           0 :     if (FpX_is_irred(pol, p)) return pol;
    2370           0 :     set_avma(av);
    2371             :   }
    2372             : }
    2373             : 
    2374             : /* return an extension of degree 2^l of F_2, assume l > 0
    2375             :  * Not stack clean. */
    2376             : static GEN
    2377         365 : f2init(long l)
    2378             : {
    2379             :   GEN Q, T, S;
    2380             :   long i, v;
    2381             : 
    2382         365 :   if (l == 1) return polcyclo(3, 0);
    2383         330 :   v = fetch_var_higher();
    2384         329 :   S = mkpoln(4, gen_1,gen_1,gen_0,gen_0); /* y(y^2 + y) */
    2385         330 :   Q = mkpoln(3, gen_1,gen_1, S); /* x^2 + x + y(y^2+y) */
    2386         331 :   setvarn(Q, v);
    2387             : 
    2388             :   /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
    2389         331 :   T = mkpoln(5, gen_1,gen_0,gen_0,gen_1,gen_1);
    2390         331 :   setvarn(T, v);
    2391             :   /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
    2392             :    * ==> x^2 + x + a(y) b irred. over K for any root b of Q
    2393             :    * ==> x^2 + x + (b^2+b)b */
    2394         331 :   for (i=2; i<l; i++) T = FpX_FpXY_resultant(T, Q, gen_2); /* minpoly(b) / F2*/
    2395         333 :   (void)delete_var(); setvarn(T,0); return T;
    2396             : }
    2397             : 
    2398             : /* return an extension of degree p^l of F_p, assume l > 0
    2399             :  * Not stack clean. */
    2400             : GEN
    2401           0 : ffinit_Artin_Shreier(GEN ip, long l)
    2402             : {
    2403           0 :   long i, v, p = itos(ip);
    2404           0 :   GEN T, Q, xp = pol_xn(p,0); /* x^p */
    2405           0 :   T = ZX_sub(xp, deg1pol_shallow(gen_1,gen_1,0)); /* x^p - x - 1 */
    2406           0 :   if (l == 1) return T;
    2407             : 
    2408           0 :   v = fetch_var_higher();
    2409           0 :   setvarn(xp, v);
    2410           0 :   Q = ZX_sub(pol_xn(2*p-1,0), pol_xn(p,0));
    2411           0 :   Q = gsub(xp, deg1pol_shallow(gen_1, Q, v)); /* x^p - x - (y^(2p-1)-y^p) */
    2412           0 :   for (i = 2; i <= l; ++i) T = FpX_FpXY_resultant(T, Q, ip);
    2413           0 :   (void)delete_var(); setvarn(T,0); return T;
    2414             : }
    2415             : 
    2416             : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
    2417             : static long
    2418       22318 : fpinit_check(GEN p, long n, long l)
    2419             : {
    2420             :   ulong q;
    2421       22318 :   if (!uisprime(n)) return 0;
    2422       14112 :   q = umodiu(p,n); if (!q) return 0;
    2423       12075 :   return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    2424             : }
    2425             : 
    2426             : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
    2427             :  * Return an irreducible polynomial of degree l over F_p.
    2428             :  * Variant of Adleman and Lenstra "Finding irreducible polynomials over
    2429             :  * finite fields", ACM, 1986 (5) 350--355.
    2430             :  * Not stack clean */
    2431             : static GEN
    2432        5481 : fpinit(GEN p, long l)
    2433             : {
    2434        5481 :   ulong n = 1+l;
    2435        5481 :   while (!fpinit_check(p,n,l)) n += l;
    2436        5481 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    2437        5481 :   return FpX_red(polsubcyclo(n,l,0),p);
    2438             : }
    2439             : 
    2440             : static GEN
    2441        5259 : ffinit_fact(GEN p, long n)
    2442             : {
    2443        5259 :   GEN P, F = gel(factoru_pow(n),3);
    2444        5257 :   long i, l = lg(F);
    2445        5257 :   P= cgetg(l, t_VEC);
    2446        5259 :   if (!odd(n) && absequaliu(p, 2))
    2447         365 :     gel(P,1) = f2init(vals(n)); /* if n is even, F[1] = 2^vals(n)*/
    2448             :   else
    2449        4893 :     gel(P,1) = fpinit(p, F[1]);
    2450        5849 :   for (i = 2; i < l; ++i)
    2451         588 :     gel(P,i) = fpinit(p, F[i]);
    2452        5261 :   return FpXV_direct_compositum(P, p);
    2453             : }
    2454             : 
    2455             : static GEN
    2456        7632 : init_Fq_i(GEN p, long n, long v)
    2457             : {
    2458             :   GEN P;
    2459        7632 :   if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
    2460        7632 :   if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
    2461        7632 :   if (signe(p) <= 0) pari_err_PRIME("ffinit",p);
    2462        7632 :   if (v < 0) v = 0;
    2463        7632 :   if (n == 1) return pol_x(v);
    2464        7380 :   if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
    2465        5258 :   P = ffinit_fact(p,n);
    2466        5261 :   setvarn(P, v); return P;
    2467             : }
    2468             : GEN
    2469        7282 : init_Fq(GEN p, long n, long v)
    2470             : {
    2471        7282 :   pari_sp av = avma;
    2472        7282 :   return gerepileupto(av, init_Fq_i(p, n, v));
    2473             : }
    2474             : GEN
    2475         350 : ffinit(GEN p, long n, long v)
    2476             : {
    2477         350 :   pari_sp av = avma;
    2478         350 :   return gerepileupto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
    2479             : }
    2480             : 
    2481             : GEN
    2482        3178 : ffnbirred(GEN p, long n)
    2483             : {
    2484        3178 :   pari_sp av = avma;
    2485        3178 :   GEN s = powiu(p,n), F = factoru(n), D = divisorsu_moebius(gel(F, 1));
    2486        3178 :   long j, l = lg(D);
    2487        6797 :   for (j = 2; j < l; j++) /* skip d = 1 */
    2488             :   {
    2489        3619 :     long md = D[j]; /* mu(d) * d, d squarefree */
    2490        3619 :     GEN pd = powiu(p, n / labs(md)); /* p^{n/d} */
    2491        3619 :     s = md > 0? addii(s, pd): subii(s,pd);
    2492             :   }
    2493        3178 :   return gerepileuptoint(av, diviuexact(s, n));
    2494             : }
    2495             : 
    2496             : GEN
    2497         441 : ffsumnbirred(GEN p, long n)
    2498             : {
    2499         441 :   pari_sp av = avma, av2;
    2500         441 :   GEN q, t = p, v = vecfactoru(1, n);
    2501             :   long i;
    2502         441 :   q = cgetg(n+1,t_VEC); gel(q,1) = p;
    2503         441 :   for (i=2; i<=n; i++) gel(q,i) = mulii(gel(q,i-1), p);
    2504         441 :   av2 = avma;
    2505        1589 :   for (i=2; i<=n; i++)
    2506             :   {
    2507        1148 :     GEN s = gel(q,i), F = gel(v,i), D = divisorsu_moebius(gel(F,1));
    2508        1148 :     long j, l = lg(D);
    2509        2534 :     for (j = 2; j < l; j++) /* skip 1 */
    2510             :     {
    2511        1386 :       long md = D[j];
    2512        1386 :       GEN pd = gel(q, i / labs(md)); /* p^{i/d} */
    2513        1386 :       s = md > 0? addii(s, pd): subii(s, pd);
    2514             :     }
    2515        1148 :     t = gerepileuptoint(av2, addii(t, diviuexact(s, i)));
    2516             :   }
    2517         441 :   return gerepileuptoint(av, t);
    2518             : }
    2519             : 
    2520             : GEN
    2521         140 : ffnbirred0(GEN p, long n, long flag)
    2522             : {
    2523         140 :   if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
    2524         140 :   if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
    2525         140 :   switch(flag)
    2526             :   {
    2527          70 :     case 0: return ffnbirred(p, n);
    2528          70 :     case 1: return ffsumnbirred(p, n);
    2529             :   }
    2530           0 :   pari_err_FLAG("ffnbirred");
    2531             :   return NULL; /* LCOV_EXCL_LINE */
    2532             : }
    2533             : 
    2534             : static void
    2535        1988 : checkmap(GEN m, const char *s)
    2536             : {
    2537        1988 :   if (typ(m)!=t_VEC || lg(m)!=3 || typ(gel(m,1))!=t_FFELT)
    2538           0 :     pari_err_TYPE(s,m);
    2539        1988 : }
    2540             : 
    2541             : GEN
    2542         175 : ffembed(GEN a, GEN b)
    2543             : {
    2544         175 :   pari_sp av = avma;
    2545         175 :   GEN p, Ta, Tb, g, r = NULL;
    2546         175 :   if (typ(a)!=t_FFELT) pari_err_TYPE("ffembed",a);
    2547         175 :   if (typ(b)!=t_FFELT) pari_err_TYPE("ffembed",b);
    2548         175 :   p = FF_p_i(a); g = FF_gen(a);
    2549         175 :   if (!equalii(p, FF_p_i(b))) pari_err_MODULUS("ffembed",a,b);
    2550         175 :   Ta = FF_mod(a);
    2551         175 :   Tb = FF_mod(b);
    2552         175 :   if (degpol(Tb)%degpol(Ta)!=0)
    2553           7 :     pari_err_DOMAIN("ffembed",GENtostr_raw(a),"is not a subfield of",b,a);
    2554         168 :   r = gel(FFX_roots(Ta, b), 1);
    2555         168 :   return gerepilecopy(av, mkvec2(g,r));
    2556             : }
    2557             : 
    2558             : GEN
    2559          84 : ffextend(GEN a, GEN P, long v)
    2560             : {
    2561          84 :   pari_sp av = avma;
    2562             :   long n;
    2563             :   GEN p, T, R, g, m;
    2564          84 :   if (typ(a)!=t_FFELT) pari_err_TYPE("ffextend",a);
    2565          84 :   T = a; p = FF_p_i(a);
    2566          84 :   if (typ(P)!=t_POL || !RgX_is_FpXQX(P,&T,&p)) pari_err_TYPE("ffextend", P);
    2567          42 :   if (!FF_samefield(a, T)) pari_err_MODULUS("ffextend",a,T);
    2568          42 :   if (v < 0) v = varn(P);
    2569          42 :   n = FF_f(T) * degpol(P); R = ffinit(p, n, v); g = ffgen(R, v);
    2570          42 :   m = ffembed(a, g);
    2571          42 :   R = FFX_roots(ffmap(m, P),g);
    2572          42 :   return gerepilecopy(av, mkvec2(gel(R,1), m));
    2573             : }
    2574             : 
    2575             : GEN
    2576          42 : fffrobenius(GEN a, long n)
    2577             : {
    2578             :   GEN g;
    2579          42 :   if (typ(a)!=t_FFELT) pari_err_TYPE("fffrobenius",a);
    2580          42 :   retmkvec2(g=FF_gen(a), FF_Frobenius(g, n));
    2581             : }
    2582             : 
    2583             : GEN
    2584         126 : ffinvmap(GEN m)
    2585             : {
    2586         126 :   pari_sp av = avma;
    2587             :   long i, l;
    2588         126 :   GEN T, F, a, g, r, f = NULL;
    2589         126 :   checkmap(m, "ffinvmap");
    2590         126 :   a = gel(m,1); r = gel(m,2);
    2591         126 :   g = FF_gen(a);
    2592         126 :   T = FF_mod(r);
    2593         126 :   F = gel(FFX_factor(T, a), 1);
    2594         126 :   l = lg(F);
    2595         532 :   for(i=1; i<l; i++)
    2596             :   {
    2597         532 :     GEN s = FFX_rem(FF_to_FpXQ_i(r), gel(F, i), a);
    2598         532 :     if (degpol(s)==0 && gequal(constant_term(s),g)) { f = gel(F, i); break; }
    2599             :   }
    2600         126 :   if (f==NULL) pari_err_TYPE("ffinvmap", m);
    2601         126 :   if (degpol(f)==1) f = FF_neg_i(gel(f,2));
    2602         126 :   return gerepilecopy(av, mkvec2(FF_gen(r),f));
    2603             : }
    2604             : 
    2605             : static GEN
    2606        1092 : ffpartmapimage(const char *s, GEN r)
    2607             : {
    2608        1092 :    GEN a = NULL, p = NULL;
    2609        1092 :    if (typ(r)==t_POL && degpol(r) >= 1
    2610        1092 :       && RgX_is_FpXQX(r,&a,&p) && a && typ(a)==t_FFELT) return a;
    2611           0 :    pari_err_TYPE(s, r);
    2612             :    return NULL; /* LCOV_EXCL_LINE */
    2613             : }
    2614             : 
    2615             : static GEN
    2616        2695 : ffeltmap_i(GEN m, GEN x)
    2617             : {
    2618        2695 :    GEN r = gel(m,2);
    2619        2695 :    if (!FF_samefield(x, gel(m,1)))
    2620          84 :      pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
    2621        2611 :    if (typ(r)==t_FFELT)
    2622        1645 :      return FF_map(r, x);
    2623             :    else
    2624         966 :      return FFX_preimage(x, r, ffpartmapimage("ffmap", r));
    2625             : }
    2626             : 
    2627             : static GEN
    2628        4424 : ffmap_i(GEN m, GEN x)
    2629             : {
    2630             :   GEN y;
    2631        4424 :   long i, lx, tx = typ(x);
    2632        4424 :   switch(tx)
    2633             :   {
    2634             :     case t_FFELT:
    2635        2527 :       return ffeltmap_i(m, x);
    2636             :     case t_POL: case t_RFRAC: case t_SER:
    2637             :     case t_VEC: case t_COL: case t_MAT:
    2638        1260 :       y = cgetg_copy(x, &lx);
    2639        1260 :       for (i=1; i<lontyp[tx]; i++) y[i] = x[1];
    2640        4536 :       for (i=lontyp[tx]; i<lx; i++)
    2641             :       {
    2642        3318 :         GEN yi = ffmap_i(m, gel(x,i));
    2643        3276 :         if (!yi) return NULL;
    2644        3276 :         gel(y,i) = yi;
    2645             :       }
    2646        1218 :       return y;
    2647             :   }
    2648         637 :   return gcopy(x);
    2649             : }
    2650             : 
    2651             : GEN
    2652        1022 : ffmap(GEN m, GEN x)
    2653             : {
    2654        1022 :   pari_sp ltop = avma;
    2655             :   GEN y;
    2656        1022 :   checkmap(m, "ffmap");
    2657        1022 :   y = ffmap_i(m, x);
    2658        1022 :   if (y) return y;
    2659          42 :   set_avma(ltop); return cgetg(1,t_VEC);
    2660             : }
    2661             : 
    2662             : static void
    2663          84 : err_compo(GEN m, GEN n)
    2664          84 : { pari_err_DOMAIN("ffcompomap","m","domain does not contain codomain of",n,m); }
    2665             : 
    2666             : GEN
    2667         420 : ffcompomap(GEN m, GEN n)
    2668             : {
    2669         420 :   pari_sp av = avma;
    2670         420 :   GEN g = gel(n,1), r, m2, n2;
    2671         420 :   checkmap(m, "ffcompomap");
    2672         420 :   checkmap(n, "ffcompomap");
    2673         420 :   m2 = gel(m,2); n2 = gel(n,2);
    2674         420 :   switch((typ(m2)==t_POL)|((typ(n2)==t_POL)<<1))
    2675             :   {
    2676             :     case 0:
    2677          84 :       if (!FF_samefield(gel(m,1),n2)) err_compo(m,n);
    2678          42 :       r = FF_map(gel(m,2), n2);
    2679          42 :       break;
    2680             :     case 2:
    2681          84 :       r = ffmap_i(m, n2);
    2682          42 :       if (lg(r) == 1) err_compo(m,n);
    2683          42 :       break;
    2684             :     case 1:
    2685         168 :       r = ffeltmap_i(m, n2);
    2686         126 :       if (!r)
    2687             :       {
    2688             :         GEN a, A, R, M;
    2689             :         long dm, dn;
    2690          42 :         a = ffpartmapimage("ffcompomap",m2);
    2691          42 :         A = FF_to_FpXQ_i(FF_neg(n2));
    2692          42 :         setvarn(A, 1);
    2693          42 :         R = deg1pol(gen_1, A, 0);
    2694          42 :         setvarn(R, 0);
    2695          42 :         M = gcopy(m2);
    2696          42 :         setvarn(M, 1);
    2697          42 :         r = polresultant0(R, M, 1, 0);
    2698          42 :         dm = FF_f(gel(m,1)); dn = FF_f(gel(n,1));
    2699          42 :         if (dm % dn || !FFX_ispower(r, dm/dn, a, &r)) err_compo(m,n);
    2700          42 :         setvarn(r, varn(FF_mod(g)));
    2701             :       }
    2702         126 :       break;
    2703             :     case 3:
    2704             :     {
    2705             :       GEN M, R, T, p, a;
    2706          84 :       a = ffpartmapimage("ffcompomap",n2);
    2707          84 :       if (!FF_samefield(a, gel(m,1))) err_compo(m,n);
    2708          42 :       p = FF_p_i(gel(n,1));
    2709          42 :       T = FF_mod(gel(n,1));
    2710          42 :       setvarn(T, 1);
    2711          42 :       R = RgX_to_FpXQX(n2,T,p);
    2712          42 :       setvarn(R, 0);
    2713          42 :       M = gcopy(m2);
    2714          42 :       setvarn(M, 1);
    2715          42 :       r = polresultant0(R, M, 1, 0);
    2716          42 :       setvarn(r, varn(n2));
    2717             :     }
    2718             :   }
    2719         252 :   return gerepilecopy(av, mkvec2(g,r));
    2720             : }

Generated by: LCOV version 1.13