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

Generated by: LCOV version 1.13