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 - polarit2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 2170 2397 90.5 %
Date: 2022-07-03 07:33:15 Functions: 193 203 95.1 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  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             : /**                         (second part)                             **/
      19             : /**                                                                   **/
      20             : /***********************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_pol
      25             : 
      26             : /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
      27             :  * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
      28             :  * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
      29             :  * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
      30             :  * Not memory clean in the latter case */
      31             : GEN
      32      127982 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      33             : {
      34      127982 :   long dP=degpol(P), i, k, m;
      35             :   pari_sp av1, av2;
      36             :   GEN s,y,P_lead;
      37             : 
      38      127982 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      39      127982 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      40      127982 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      41      127982 :   y = cgetg(n+2,t_COL);
      42      127982 :   if (y0)
      43             :   {
      44       13237 :     if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
      45       13237 :     m = lg(y0)-1;
      46       63987 :     for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
      47             :   }
      48             :   else
      49             :   {
      50      114745 :     m = 1;
      51      114745 :     gel(y,1) = stoi(dP);
      52             :   }
      53      127983 :   P += 2; /* strip codewords */
      54             : 
      55      127983 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      56      127983 :   if (P_lead)
      57             :   {
      58           7 :     if (N) P_lead = Fq_inv(P_lead,T,N);
      59           7 :     else if (T) P_lead = QXQ_inv(P_lead,T);
      60             :   }
      61      380104 :   for (k=m; k<=n; k++)
      62             :   {
      63      252120 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      64      725537 :     for (i=1; i<k && i<=dP; i++)
      65      473433 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      66      252104 :     if (N)
      67             :     {
      68       18760 :       s = Fq_red(s, T, N);
      69       18760 :       if (P_lead) s = Fq_mul(s, P_lead, T, N);
      70             :     }
      71      233344 :     else if (T)
      72             :     {
      73           0 :       s = grem(s, T);
      74           0 :       if (P_lead) s = grem(gmul(s, P_lead), T);
      75             :     }
      76             :     else
      77      233344 :       if (P_lead) s = gdiv(s, P_lead);
      78      252104 :     av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
      79             :   }
      80      127984 :   return y;
      81             : }
      82             : 
      83             : GEN
      84      110567 : polsym(GEN x, long n)
      85             : {
      86      110567 :   return polsym_gen(x, NULL, n, NULL,NULL);
      87             : }
      88             : 
      89             : /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
      90             : GEN
      91    85008097 : centermodii(GEN x, GEN p, GEN po2)
      92             : {
      93    85008097 :   GEN y = remii(x, p);
      94    85181398 :   switch(signe(y))
      95             :   {
      96    10524046 :     case 0: break;
      97    51485772 :     case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
      98    51346177 :       break;
      99    23504771 :     case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
     100    23304761 :       break;
     101             :   }
     102    84841793 :   return y;
     103             : }
     104             : 
     105             : static long
     106           0 : s_centermod(long x, ulong pp, ulong pps2)
     107             : {
     108           0 :   long y = x % (long)pp;
     109           0 :   if (y < 0) y += pp;
     110           0 :   return Fl_center(y, pp,pps2);
     111             : }
     112             : 
     113             : /* for internal use */
     114             : GEN
     115    12648404 : centermod_i(GEN x, GEN p, GEN ps2)
     116             : {
     117             :   long i, lx;
     118             :   pari_sp av;
     119             :   GEN y;
     120             : 
     121    12648404 :   if (!ps2) ps2 = shifti(p,-1);
     122    12646122 :   switch(typ(x))
     123             :   {
     124     1123830 :     case t_INT: return centermodii(x,p,ps2);
     125             : 
     126     6143289 :     case t_POL: lx = lg(x);
     127     6143289 :       y = cgetg(lx,t_POL); y[1] = x[1];
     128    42161020 :       for (i=2; i<lx; i++)
     129             :       {
     130    36015603 :         av = avma;
     131    36015603 :         gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
     132             :       }
     133     6145417 :       return normalizepol_lg(y, lx);
     134             : 
     135     5377932 :     case t_COL: lx = lg(x);
     136     5377932 :       y = cgetg(lx,t_COL);
     137    25451437 :       for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
     138     5376867 :       return y;
     139             : 
     140        1071 :     case t_MAT: lx = lg(x);
     141        1071 :       y = cgetg(lx,t_MAT);
     142       13531 :       for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
     143        1071 :       return y;
     144             : 
     145           0 :     case t_VECSMALL: lx = lg(x);
     146             :     {
     147           0 :       ulong pp = itou(p), pps2 = itou(ps2);
     148           0 :       y = cgetg(lx,t_VECSMALL);
     149           0 :       for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
     150           0 :       return y;
     151             :     }
     152             :   }
     153           0 :   return x;
     154             : }
     155             : 
     156             : GEN
     157     9069201 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
     158             : 
     159             : static GEN
     160         343 : RgX_Frobenius_deflate(GEN S, ulong p)
     161             : {
     162         343 :   if (degpol(S)%p)
     163           0 :     return NULL;
     164             :   else
     165             :   {
     166         343 :     GEN F = RgX_deflate(S, p);
     167         343 :     long i, l = lg(F);
     168        1043 :     for (i=2; i<l; i++)
     169             :     {
     170         721 :       GEN Fi = gel(F,i), R;
     171         721 :       if (typ(Fi)==t_POL)
     172             :       {
     173         259 :         if (signe(RgX_deriv(Fi))==0)
     174         238 :           gel(F,i) = RgX_Frobenius_deflate(gel(F, i), p);
     175          21 :         else return NULL;
     176             :       }
     177         462 :       else if (ispower(Fi, utoi(p), &R))
     178         462 :         gel(F,i) = R;
     179           0 :       else return NULL;
     180             :     }
     181         322 :     return F;
     182             :   }
     183             : }
     184             : 
     185             : static GEN
     186         245 : RgXY_squff(GEN f)
     187             : {
     188         245 :   long i, q, n = degpol(f);
     189         245 :   ulong p = itos_or_0(characteristic(f));
     190         245 :   GEN u = const_vec(n+1, pol_1(varn(f)));
     191         245 :   for(q = 1;;q *= p)
     192          84 :   {
     193         329 :     GEN t, v, tv, r = RgX_gcd(f, RgX_deriv(f));
     194         329 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
     195         126 :     t = RgX_div(f, r);
     196         126 :     if (degpol(t) > 0)
     197             :     {
     198             :       long j;
     199          28 :       for(j = 1;;j++)
     200             :       {
     201         140 :         v = RgX_gcd(r, t);
     202         140 :         tv = RgX_div(t, v);
     203         140 :         if (degpol(tv) > 0) gel(u, j*q) = tv;
     204         140 :         if (degpol(v) <= 0) break;
     205         112 :         r = RgX_div(r, v);
     206         112 :         t = v;
     207             :       }
     208          28 :       if (degpol(r) == 0) break;
     209             :     }
     210         105 :     if (!p) break;
     211         105 :     f = RgX_Frobenius_deflate(r, p);
     212         105 :     if (!f) { gel(u, q) = r; break; }
     213             :   }
     214         931 :   for (i = n; i; i--)
     215         931 :     if (degpol(gel(u,i))) break;
     216         245 :   setlg(u,i+1); return u;
     217             : }
     218             : 
     219             : /* Lmod contains modular factors of *F (NULL codes an empty slot: used factor)
     220             :  * Lfac accumulates irreducible factors as they are found.
     221             :  * p is a product of modular factors in Lmod[1..i-1] (NULL for p = 1), not
     222             :  * a rational factor of *F
     223             :  * Find an irreducible factor of *F divisible by p (by including
     224             :  * exhaustively further factors from Lmod[i..]); return 0 on failure, else 1.
     225             :  * Update Lmod, Lfac and *F */
     226             : static int
     227        9121 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
     228             : {
     229             :   GEN q;
     230        9121 :   if (i == lg(Lmod)) return 0;
     231        4704 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
     232        4501 :   if (!gel(Lmod,i)) return 0;
     233        4487 :   p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
     234        4487 :   q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
     235        4487 :   if (degpol(q))
     236             :   {
     237        4130 :     GEN R, Q = RgX_divrem(*F, q, &R);
     238        4130 :     if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
     239             :   }
     240        4158 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
     241        3892 :   return 0;
     242             : }
     243             : 
     244             : static GEN factor_domain(GEN x, GEN flag);
     245             : 
     246             : static GEN
     247         434 : ok_bloc(GEN f, GEN BLOC, ulong c)
     248             : {
     249         434 :   GEN F = poleval(f, BLOC);
     250         434 :   return issquarefree(c ? gmul(F,mkintmodu(1,c)): F)? F: NULL;
     251             : }
     252             : static GEN
     253         112 : random_FpX_monic(long n, long v, GEN p)
     254             : {
     255         112 :   long i, d = n + 2;
     256         112 :   GEN y = cgetg(d + 1, t_POL); y[1] = evalsigne(1) | evalvarn(v);
     257         364 :   for (i = 2; i < d; i++) gel(y,i) = randomi(p);
     258         112 :   gel(y,i) = gen_1; return y;
     259             : }
     260             : static GEN
     261         273 : RgXY_factor_squarefree(GEN f, GEN dom)
     262             : {
     263         273 :   pari_sp av = avma;
     264         273 :   ulong i, c = itou_or_0(residual_characteristic(f));
     265         273 :   long vy = gvar2(f), val = RgX_valrem(f, &f), n = RgXY_degreex(f);
     266         273 :   GEN y, Lmod, F = NULL, BLOC = NULL, Lfac = coltrunc_init(degpol(f)+2);
     267         273 :   GEN gc = c? utoipos(c): NULL;
     268         273 :   if (val)
     269             :   {
     270          35 :     GEN x = pol_x(varn(f));
     271          35 :     if (dom)
     272             :     {
     273          14 :       GEN c = Rg_get_1(dom);
     274          14 :       if (typ(c) != t_INT) x = RgX_Rg_mul(x, c);
     275             :     }
     276          35 :     vectrunc_append(Lfac, x); if (!degpol(f)) return Lfac;
     277             :   }
     278         259 :   y = pol_x(vy);
     279             :   for(;;)
     280             :   {
     281         322 :     for (i = 0; !c || i < c; i++)
     282             :     {
     283         322 :       BLOC = gpowgs(gaddgs(y, i), n+1);
     284         322 :       if ((F = ok_bloc(f, BLOC, c))) break;
     285         147 :       if (c)
     286             :       {
     287         112 :         BLOC = random_FpX_monic(n, vy, gc);
     288         112 :         if ((F = ok_bloc(f, BLOC, c))) break;
     289             :       }
     290             :     }
     291         259 :     if (!c || i < c) break;
     292           0 :     n++;
     293             :   }
     294         259 :   if (DEBUGLEVEL >= 2)
     295           0 :     err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
     296         259 :   Lmod = gel(factor_domain(F,dom),1);
     297         259 :   if (DEBUGLEVEL >= 2)
     298           0 :     err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
     299         259 :   (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
     300         259 :   if (degpol(f)) vectrunc_append(Lfac, f);
     301         259 :   return gerepilecopy(av, Lfac);
     302             : }
     303             : 
     304             : static GEN
     305         245 : FE_matconcat(GEN F, GEN E, long l)
     306             : {
     307         245 :   setlg(E,l); E = shallowconcat1(E);
     308         245 :   setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
     309             : }
     310             : 
     311             : static int
     312         385 : gen_cmp_RgXY(void *data, GEN x, GEN y)
     313             : {
     314         385 :   long vx = varn(x), vy = varn(y);
     315         385 :   return (vx == vy)? gen_cmp_RgX(data, x, y): -varncmp(vx, vy);
     316             : }
     317             : static GEN
     318         245 : RgXY_factor(GEN f, GEN dom)
     319             : {
     320         245 :   pari_sp av = avma;
     321             :   GEN C, F, E, cf, V;
     322             :   long i, j, l;
     323         245 :   if (dom) { GEN c = Rg_get_1(dom); if (typ(c) != t_INT) f = RgX_Rg_mul(f,c); }
     324         245 :   cf = content(f);
     325         245 :   V = RgXY_squff(gdiv(f, cf)); l = lg(V);
     326         245 :   C = factor_domain(cf, dom);
     327         245 :   F = cgetg(l+1, t_VEC); gel(F,1) = gel(C,1);
     328         245 :   E = cgetg(l+1, t_VEC); gel(E,1) = gel(C,2);
     329         756 :   for (i=1, j=2; i < l; i++)
     330             :   {
     331         511 :     GEN v = gel(V,i);
     332         511 :     if (degpol(v))
     333             :     {
     334         273 :       gel(F,j) = v = RgXY_factor_squarefree(v, dom);
     335         273 :       gel(E,j) = const_col(lg(v)-1, utoipos(i));
     336         273 :       j++;
     337             :     }
     338             :   }
     339         245 :   f = FE_matconcat(F,E,j);
     340         245 :   (void)sort_factor(f,(void*)cmp_universal, &gen_cmp_RgXY);
     341         245 :   return gerepilecopy(av, f);
     342             : }
     343             : 
     344             : /***********************************************************************/
     345             : /**                                                                   **/
     346             : /**                          FACTORIZATION                            **/
     347             : /**                                                                   **/
     348             : /***********************************************************************/
     349             : static long RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var);
     350             : #define assign_or_fail(x,y) { GEN __x = x;\
     351             :   if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
     352             : }
     353             : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
     354             : 
     355             : static const long tsh = 6;
     356             : #define code(t1,t2) ((t1 << 6) | t2)
     357             : void
     358       89456 : RgX_type_decode(long x, long *t1, long *t2)
     359             : {
     360       89456 :   *t1 = x >> tsh;
     361       89456 :   *t2 = (x & ((1L<<tsh)-1));
     362       89456 : }
     363             : int
     364     9112494 : RgX_type_is_composite(long t) { return t >= tsh; }
     365             : 
     366             : static int
     367  2139029264 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     368             : {
     369             :   long j;
     370  2139029264 :   switch(typ(c))
     371             :   {
     372  1603450618 :     case t_INT:
     373  1603450618 :       break;
     374    25922044 :     case t_FRAC:
     375    25922044 :       t[1]=1; break;
     376             :       break;
     377   234776474 :     case t_REAL:
     378   234776474 :       update_prec(precision(c), pa);
     379   234775668 :       t[2]=1; break;
     380    23503244 :     case t_INTMOD:
     381    23503244 :       assign_or_fail(gel(c,1),p);
     382    23503244 :       t[3]=1; break;
     383     1870168 :     case t_FFELT:
     384     1870168 :       if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
     385     1870168 :       assign_or_fail(FF_p_i(c),p);
     386     1870168 :       t[5]=1; break;
     387   215557010 :     case t_COMPLEX:
     388   646668530 :       for (j=1; j<=2; j++)
     389             :       {
     390   431111641 :         GEN d = gel(c,j);
     391   431111641 :         switch(typ(d))
     392             :         {
     393     1643416 :           case t_INT: case t_FRAC:
     394     1643416 :             if (!*t2) *t2 = t_COMPLEX;
     395     1643416 :             t[1]=1; break;
     396   429468190 :           case t_REAL:
     397   429468190 :             update_prec(precision(d), pa);
     398   429468076 :             if (!*t2) *t2 = t_COMPLEX;
     399   429468076 :             t[2]=1; break;
     400          14 :           case t_INTMOD:
     401          14 :             assign_or_fail(gel(d,1),p);
     402          14 :             if (!signe(*p) || mod4(*p) != 3) return 0;
     403           7 :             if (!*t2) *t2 = t_COMPLEX;
     404           7 :             t[3]=1; break;
     405          21 :           case t_PADIC:
     406          21 :             update_prec(precp(d)+valp(d), pa);
     407          21 :             assign_or_fail(gel(d,2),p);
     408          21 :             if (!*t2) *t2 = t_COMPLEX;
     409          21 :             t[7]=1; break;
     410           0 :           default: return 0;
     411             :         }
     412             :       }
     413   215556889 :       if (!t[2]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
     414   215556889 :       break;
     415     3296524 :     case t_PADIC:
     416     3296524 :       update_prec(precp(c)+valp(c), pa);
     417     3296524 :       assign_or_fail(gel(c,2),p);
     418     3296524 :       t[7]=1; break;
     419        1491 :     case t_QUAD:
     420        1491 :       assign_or_fail(gel(c,1),pol);
     421        4473 :       for (j=2; j<=3; j++)
     422             :       {
     423        2982 :         GEN d = gel(c,j);
     424        2982 :         switch(typ(d))
     425             :         {
     426        2947 :           case t_INT: case t_FRAC:
     427        2947 :             t[8]=1; break;
     428          28 :           case t_INTMOD:
     429          28 :             assign_or_fail(gel(d,1),p);
     430          28 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     431          28 :             t[3]=1; break;
     432           7 :           case t_PADIC:
     433           7 :             update_prec(precp(d)+valp(d), pa);
     434           7 :             assign_or_fail(gel(d,2),p);
     435           7 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     436           7 :             t[7]=1; break;
     437           0 :           default: return 0;
     438             :         }
     439             :       }
     440        1491 :       break;
     441     3694434 :     case t_POLMOD:
     442     3694434 :       assign_or_fail(gel(c,1),pol);
     443     3694219 :       if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1))) return 0;
     444    11055280 :       for (j=1; j<=2; j++)
     445             :       {
     446             :         GEN pbis, polbis;
     447             :         long pabis;
     448     7376858 :         *t2 = t_POLMOD;
     449     7376858 :         switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
     450             :         {
     451     3944226 :           case t_INT: break;
     452      844426 :           case t_FRAC:   t[1]=1; break;
     453     2576875 :           case t_INTMOD: t[3]=1; break;
     454           7 :           case t_PADIC:  t[7]=1; update_prec(pabis,pa); break;
     455       11337 :           default: return 0;
     456             :         }
     457     7365534 :         if (pbis) assign_or_fail(pbis,p);
     458     7365534 :         if (polbis) assign_or_fail(polbis,pol);
     459             :       }
     460     3678422 :       break;
     461     5325148 :     case t_RFRAC: t[10] = 1;
     462     5325148 :       if (!settype(gel(c,1),t,p,pol,pa,ff,t2,var)) return 0;
     463     5325148 :       c = gel(c,2); /* fall through */
     464    26955556 :     case t_POL: t[10] = 1;
     465    26955556 :       if (!RgX_settype(c,t,p,pol,pa,ff,t2,var)) return 0;
     466    26968484 :       if (*var == NO_VARIABLE) { *var = varn(c); break; }
     467             :       /* if more than one free var, ensure varn() == *var fails. FIXME: should
     468             :        * keep the list of all variables, later t_POLMOD may cancel them */
     469    19264296 :       if (*var != varn(c)) *var = MAXVARN+1;
     470    19264296 :       break;
     471        1701 :     default: return 0;
     472             :   }
     473  2139023552 :   return 1;
     474             : }
     475             : /* t[0] unused. Other values, if set, indicate a coefficient of type
     476             :  * t[1] : t_FRAC
     477             :  * t[2] : t_REAL
     478             :  * t[3] : t_INTMOD
     479             :  * t[4] : Unused
     480             :  * t[5] : t_FFELT
     481             :  * t[6] : Unused
     482             :  * t[7] : t_PADIC
     483             :  * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
     484             :  * t[9]:  Unused
     485             :  * t[10]: t_POL (recursive factorisation) */
     486             : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
     487             :  * given by t) */
     488             : static long
     489   155400011 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
     490             : {
     491   155400011 :   if (t[10] && (!*pol || var!=varn(*pol))) return t_POL;
     492   147707717 :   if (t2) /* polmod/quad/complex of intmod/padic */
     493             :   {
     494    20875416 :     if (t[2] && (t[3]||t[7])) return 0;
     495    20875416 :     if (t[3]) return code(t2,t_INTMOD);
     496    20846163 :     if (t[7]) return code(t2,t_PADIC);
     497    20846114 :     if (t[2]) return t_COMPLEX;
     498      487747 :     if (t[1]) return code(t2,t_FRAC);
     499      223391 :     return code(t2,t_INT);
     500             :   }
     501   126832301 :   if (t[5]) /* ffelt */
     502             :   {
     503      219745 :     if (t[2]||t[8]||t[9]) return 0;
     504      219745 :     *pol=ff; return t_FFELT;
     505             :   }
     506   126612556 :   if (t[2]) /* inexact, real */
     507             :   {
     508    41580476 :     if (t[3]||t[7]||t[9]) return 0;
     509    41580508 :     return t_REAL;
     510             :   }
     511    85032080 :   if (t[10]) return t_POL;
     512    85032080 :   if (t[8]) return code(t_QUAD,t_INT);
     513    85031443 :   if (t[3]) return t_INTMOD;
     514    81339728 :   if (t[7]) return t_PADIC;
     515    80743370 :   if (t[1]) return t_FRAC;
     516    74884609 :   return t_INT;
     517             : }
     518             : 
     519             : static long
     520   229484839 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     521             : {
     522   229484839 :   long i, lx = lg(x);
     523   883317205 :   for (i=2; i<lx; i++)
     524   653854533 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     525   229462672 :   return 1;
     526             : }
     527             : 
     528             : static long
     529   243040370 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     530             : {
     531   243040370 :   long i, l = lg(x);
     532  1722290471 :   for (i = 1; i<l; i++)
     533  1479252631 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     534   243037840 :   return 1;
     535             : }
     536             : 
     537             : static long
     538    49477089 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     539             : {
     540    49477089 :   long i, l = lg(x);
     541   250928309 :   for (i = 1; i < l; i++)
     542   201453673 :     if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     543    49474636 :   return 1;
     544             : }
     545             : 
     546             : long
     547    16489358 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
     548             : {
     549    16489358 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0};
     550    16489358 :   long t2 = 0, var = NO_VARIABLE;
     551    16489358 :   GEN ff = NULL;
     552    16489358 :   *p = *pol = NULL; *pa = LONG_MAX;
     553    16489358 :   switch(typ(x))
     554             :   {
     555      602080 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
     556             :     case t_COMPLEX: case t_PADIC: case t_QUAD:
     557      602080 :       if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     558      602082 :       break;
     559    15304640 :     case t_POL: case t_SER:
     560    15304640 :       if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     561    15304650 :       break;
     562          14 :     case t_VEC: case t_COL:
     563          14 :       if(!RgC_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
     564          14 :       break;
     565         126 :     case t_MAT:
     566         126 :       if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
     567         126 :       break;
     568      582498 :     default: return 0;
     569             :   }
     570    15906872 :   return choosetype(t,t2,ff,pol,var);
     571             : }
     572             : 
     573             : long
     574     1256868 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     575             : {
     576     1256868 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     577     1256868 :   long t2 = 0, var = NO_VARIABLE;
     578     1256868 :   GEN ff = NULL;
     579     1256868 :   *p = *pol = NULL; *pa = LONG_MAX;
     580     1256868 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     581     1256820 :   return choosetype(t,t2,ff,pol,var);
     582             : }
     583             : 
     584             : long
     585         294 : RgX_Rg_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     586             : {
     587         294 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     588         294 :   long t2 = 0, var = NO_VARIABLE;
     589         294 :   GEN ff = NULL;
     590         294 :   *p = *pol = NULL; *pa = LONG_MAX;
     591         294 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     592         294 :   if (!settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     593         294 :   return choosetype(t,t2,ff,pol,var);
     594             : }
     595             : 
     596             : long
     597    91722912 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     598             : {
     599    91722912 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     600    91722912 :   long t2 = 0, var = NO_VARIABLE;
     601    91722912 :   GEN ff = NULL;
     602    91722912 :   *p = *pol = NULL; *pa = LONG_MAX;
     603   183434474 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     604    91723922 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     605    91711444 :   return choosetype(t,t2,ff,pol,var);
     606             : }
     607             : 
     608             : long
     609      839214 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
     610             : {
     611      839214 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     612      839214 :   long t2 = 0, var = NO_VARIABLE;
     613      839214 :   GEN ff = NULL;
     614      839214 :   *p = *pol = NULL; *pa = LONG_MAX;
     615     1675937 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     616     1673445 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
     617      839215 :       !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
     618      836720 :   return choosetype(t,t2,ff,pol,var);
     619             : }
     620             : 
     621             : long
     622      316430 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
     623             : {
     624      316430 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     625      316430 :   long t2 = 0, var = NO_VARIABLE;
     626      316430 :   GEN ff = NULL;
     627      316430 :   *p = *pol = NULL; *pa = LONG_MAX;
     628      316430 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     629      315459 :   return choosetype(t,t2,ff,pol,var);
     630             : }
     631             : 
     632             : long
     633          70 : RgV_type(GEN x, GEN *p, GEN *pol, long *pa)
     634             : {
     635          70 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     636          70 :   long t2 = 0, var = NO_VARIABLE;
     637          70 :   GEN ff = NULL;
     638          70 :   *p = *pol = NULL; *pa = LONG_MAX;
     639          70 :   if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     640          70 :   return choosetype(t,t2,ff,pol,var);
     641             : }
     642             : 
     643             : long
     644         203 : RgV_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     645             : {
     646         203 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     647         203 :   long t2 = 0, var = NO_VARIABLE;
     648         203 :   GEN ff = NULL;
     649         203 :   *p = *pol = NULL; *pa = LONG_MAX;
     650         406 :   if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     651         203 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     652         203 :   return choosetype(t,t2,ff,pol,var);
     653             : }
     654             : 
     655             : long
     656    41587557 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     657             : {
     658    41587557 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     659    41587557 :   long t2 = 0, var = NO_VARIABLE;
     660    41587557 :   GEN ff = NULL;
     661    41587557 :   *p = *pol = NULL; *pa = LONG_MAX;
     662    83174773 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     663    41588439 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     664    41586417 :   return choosetype(t,t2,ff,pol,var);
     665             : }
     666             : 
     667             : long
     668     3786675 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     669             : {
     670     3786675 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     671     3786675 :   long t2 = 0, var = NO_VARIABLE;
     672     3786675 :   GEN ff = NULL;
     673     3786675 :   *p = *pol = NULL; *pa = LONG_MAX;
     674     7572893 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     675     3786801 :       !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     676     3786113 :   return choosetype(t,t2,ff,pol,var);
     677             : }
     678             : 
     679             : GEN
     680       59326 : factor0(GEN x, GEN flag)
     681             : {
     682             :   ulong B;
     683       59326 :   long tx = typ(x);
     684       59326 :   if (!flag) return factor(x);
     685         224 :   if ((tx != t_INT && tx!=t_FRAC) || typ(flag) != t_INT)
     686         175 :     return factor_domain(x, flag);
     687          49 :   if (signe(flag) < 0) pari_err_FLAG("factor");
     688          49 :   switch(lgefint(flag))
     689             :   {
     690           7 :     case 2: B = 0; break;
     691          42 :     case 3: B = flag[2]; break;
     692           0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     693             :              return NULL; /*LCOV_EXCL_LINE*/
     694             :   }
     695          49 :   return boundfact(x, B);
     696             : }
     697             : 
     698             : GEN
     699      134185 : deg1_from_roots(GEN L, long v)
     700             : {
     701      134185 :   long i, l = lg(L);
     702      134185 :   GEN z = cgetg(l,t_COL);
     703      411278 :   for (i=1; i<l; i++)
     704      277094 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     705      134184 :   return z;
     706             : }
     707             : GEN
     708       63852 : roots_from_deg1(GEN x)
     709             : {
     710       63852 :   long i,l = lg(x);
     711       63852 :   GEN r = cgetg(l,t_VEC);
     712      392640 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     713       63853 :   return r;
     714             : }
     715             : 
     716             : static GEN
     717          42 : gauss_factor_p(GEN p)
     718             : {
     719          42 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     720          42 :   return mkcomplex(a, b);
     721             : }
     722             : 
     723             : static GEN
     724          49 : gauss_primpart(GEN x, GEN *c)
     725             : {
     726          49 :   GEN a = real_i(x), b = imag_i(x), n = gcdii(a, b);
     727          49 :   *c = n; if (n == gen_1) return x;
     728          49 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     729             : }
     730             : 
     731             : static GEN
     732          70 : gauss_primpart_try(GEN x, GEN c)
     733             : {
     734             :   GEN r, y;
     735          70 :   if (typ(x) == t_INT)
     736             :   {
     737          42 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     738             :   }
     739             :   else
     740             :   {
     741          28 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     742          28 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     743          14 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     744             :   }
     745          56 :   return y;
     746             : }
     747             : 
     748             : static int
     749          91 : gauss_cmp(GEN x, GEN y)
     750             : {
     751             :   int v;
     752          91 :   if (typ(x) != t_COMPLEX)
     753           0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     754          91 :   if (typ(y) != t_COMPLEX) return 1;
     755          63 :   v = cmpii(gel(x,2), gel(y,2));
     756          63 :   if (v) return v;
     757          28 :   return gcmp(gel(x,1), gel(y,1));
     758             : }
     759             : 
     760             : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
     761             : static GEN
     762         469 : gauss_normal(GEN x)
     763             : {
     764         469 :   if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
     765         392 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     766         392 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     767         392 :   return x;
     768             : }
     769             : 
     770             : static GEN
     771          49 : gauss_factor(GEN x)
     772             : {
     773          49 :   pari_sp av = avma;
     774          49 :   GEN a = real_i(x), b = imag_i(x), d = gen_1, n, y, fa, P, E, P2, E2;
     775          49 :   long t1 = typ(a);
     776          49 :   long t2 = typ(b), i, j, l, exp = 0;
     777          49 :   if (t1 == t_FRAC) d = gel(a,2);
     778          49 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     779          49 :   if (d == gen_1) y = x;
     780             :   else
     781             :   {
     782          21 :     y = gmul(x, d);
     783          21 :     a = real_i(y); t1 = typ(a);
     784          21 :     b = imag_i(y); t2 = typ(b);
     785             :   }
     786          49 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     787          49 :   y = gauss_primpart(y, &n);
     788          49 :   fa = factor(cxnorm(y));
     789          49 :   P = gel(fa,1);
     790          49 :   E = gel(fa,2); l = lg(P);
     791          49 :   P2 = cgetg(l, t_COL);
     792          49 :   E2 = cgetg(l, t_COL);
     793         105 :   for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
     794             :   { /* either p = 2 (ramified) or those factors split in Q(i) */
     795          56 :     GEN p = gel(P,i), w, w2, t, we, pe;
     796          56 :     long v, e = itos(gel(E,i));
     797          56 :     int is2 = absequaliu(p, 2);
     798          56 :     w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     799          56 :     w2 = gauss_normal( conj_i(w) );
     800             :     /* w * w2 * I^3 = p, w2 = conj(w) * I */
     801          56 :     pe = powiu(p, e);
     802          56 :     we = gpowgs(w, e);
     803          56 :     t = gauss_primpart_try( gmul(y, conj_i(we)), pe );
     804          56 :     if (t) y = t; /* y /= w^e */
     805             :     else {
     806             :       /* y /= conj(w)^e, should be y /= w2^e */
     807          14 :       y = gauss_primpart_try( gmul(y, we), pe );
     808          14 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     809             :     }
     810          56 :     gel(P,i) = w;
     811          56 :     v = Z_pvalrem(n, p, &n);
     812          56 :     if (v) {
     813           7 :       exp -= v; /* += 3*v mod 4 */
     814           7 :       if (is2) v <<= 1; /* 2 = w^2 I^3 */
     815             :       else {
     816           0 :         gel(P2,j) = w2;
     817           0 :         gel(E2,j) = utoipos(v); j++;
     818             :       }
     819           7 :       gel(E,i) = stoi(e + v);
     820             :     }
     821          56 :     v = Z_pvalrem(d, p, &d);
     822          56 :     if (v) {
     823           7 :       exp += v; /* -= 3*v mod 4 */
     824           7 :       if (is2) v <<= 1; /* 2 is ramified */
     825             :       else {
     826           7 :         gel(P2,j) = w2;
     827           7 :         gel(E2,j) = utoineg(v); j++;
     828             :       }
     829           7 :       gel(E,i) = stoi(e - v);
     830             :     }
     831          56 :     exp &= 3;
     832             :   }
     833          49 :   if (j > 1) {
     834           7 :     long k = 1;
     835           7 :     GEN P1 = cgetg(l, t_COL);
     836           7 :     GEN E1 = cgetg(l, t_COL);
     837             :     /* remove factors with exponent 0 */
     838          14 :     for (i = 1; i < l; i++)
     839           7 :       if (signe(gel(E,i)))
     840             :       {
     841           0 :         gel(P1,k) = gel(P,i);
     842           0 :         gel(E1,k) = gel(E,i);
     843           0 :         k++;
     844             :       }
     845           7 :     setlg(P1, k); setlg(E1, k);
     846           7 :     setlg(P2, j); setlg(E2, j);
     847           7 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     848             :   }
     849          49 :   if (!equali1(n) || !equali1(d))
     850             :   {
     851          28 :     GEN Fa = factor(Qdivii(n, d));
     852          28 :     P = gel(Fa,1); l = lg(P);
     853          28 :     E = gel(Fa,2);
     854          70 :     for (i = 1; i < l; i++)
     855             :     {
     856          42 :       GEN w, p = gel(P,i);
     857             :       long e;
     858             :       int is2;
     859          42 :       switch(mod4(p))
     860             :       {
     861          14 :         case 3: continue;
     862          14 :         case 2: is2 = 1; break;
     863          14 :         default:is2 = 0; break;
     864             :       }
     865          28 :       e = itos(gel(E,i));
     866          28 :       w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     867          28 :       gel(P,i) = w;
     868          28 :       if (is2)
     869          14 :         gel(E,i) = stoi(2*e);
     870             :       else
     871             :       {
     872          14 :         P = vec_append(P, gauss_normal( conj_i(w) ));
     873          14 :         E = vec_append(E, gel(E,i));
     874             :       }
     875          28 :       exp -= e; /* += 3*e mod 4 */
     876          28 :       exp &= 3;
     877             :     }
     878          28 :     gel(Fa,1) = P;
     879          28 :     gel(Fa,2) = E;
     880          28 :     fa = famat_mul_shallow(fa, Fa);
     881             :   }
     882          49 :   fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
     883             : 
     884          49 :   y = gmul(y, powIs(exp));
     885          49 :   if (!gequal1(y)) {
     886          35 :     gel(fa,1) = vec_prepend(gel(fa,1), y);
     887          35 :     gel(fa,2) = vec_prepend(gel(fa,2), gen_1);
     888             :   }
     889          49 :   return gerepilecopy(av, fa);
     890             : }
     891             : 
     892             : GEN
     893        9754 : Q_factor_limit(GEN x, ulong lim)
     894             : {
     895        9754 :   pari_sp av = avma;
     896             :   GEN a, b;
     897        9754 :   if (typ(x) == t_INT) return Z_factor_limit(x, lim);
     898        5292 :   a = Z_factor_limit(gel(x,1), lim);
     899        5292 :   b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
     900        5292 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     901             : }
     902             : GEN
     903       20812 : Q_factor(GEN x)
     904             : {
     905       20812 :   pari_sp av = avma;
     906             :   GEN a, b;
     907       20812 :   if (typ(x) == t_INT) return Z_factor(x);
     908          33 :   a = Z_factor(gel(x,1));
     909          35 :   b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
     910          35 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     911             : }
     912             : /* replace t_QUAD/t_COMPLEX coeffs by t_POLMOD in T */
     913             : static GEN
     914         126 : RgX_fix_quad(GEN x, GEN T)
     915             : {
     916         126 :   long i, l, v = varn(T);
     917         126 :   GEN y = cgetg_copy(x,&l);
     918         630 :   for (i = 2; i < l; i++)
     919             :   {
     920         504 :     GEN c = gel(x,i);
     921         504 :     switch(typ(c))
     922             :     {
     923          56 :       case t_QUAD: c++;/* fall through */
     924          98 :       case t_COMPLEX: c = deg1pol_shallow(gel(c,2),gel(c,1),v);
     925             :     }
     926         504 :     gel(y,i) = c;
     927             :   }
     928         126 :   y[1] = x[1]; return y;
     929             : }
     930             : 
     931             : static GEN
     932       13629 : RgX_factor(GEN x, GEN dom)
     933             : {
     934             :   pari_sp av;
     935             :   long pa, v, lx, r1, i;
     936             :   GEN  p, pol, y, p1, p2;
     937       13629 :   long tx = dom ? RgX_Rg_type(x,dom,&p,&pol,&pa): RgX_type(x,&p,&pol,&pa);
     938       13629 :   switch(tx)
     939             :   {
     940           7 :     case 0: pari_err_IMPL("factor for general polynomials");
     941         245 :     case t_POL: return RgXY_factor(x, dom);
     942       12600 :     case t_INT: return ZX_factor(x);
     943           7 :     case t_FRAC: return QX_factor(x);
     944         329 :     case t_INTMOD: return factmod(x,p);
     945          42 :     case t_PADIC: return factorpadic(x,p,pa);
     946          98 :     case t_FFELT: return FFX_factor(x,pol);
     947             : 
     948          21 :     case t_COMPLEX: y = cgetg(3,t_MAT);
     949          21 :       av = avma; p1 = deg1_from_roots(roots(x,pa), varn(x));
     950          21 :       gel(y,1) = p1 = gerepileupto(av, p1);
     951          21 :       gel(y,2) = const_col(lg(p1)-1, gen_1); return y;
     952             : 
     953          28 :     case t_REAL: y=cgetg(3,t_MAT); v=varn(x);
     954          28 :       av=avma; p1=cleanroots(x,pa);
     955          28 :       lx = lg(p1);
     956          70 :       for (r1 = 1; r1 < lx; r1++)
     957          49 :         if (typ(gel(p1,r1)) == t_COMPLEX) break;
     958          28 :       lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
     959          70 :       for (i = 1; i < r1; i++)
     960          42 :         gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
     961          35 :       for (   ; i < lx; i++)
     962             :       {
     963           7 :         GEN a = gel(p1,2*i-r1);
     964           7 :         p = cgetg(5, t_POL); gel(p2,i) = p;
     965           7 :         p[1] = x[1];
     966           7 :         gel(p,2) = gnorm(a);
     967           7 :         gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
     968           7 :         gel(p,4) = gen_1;
     969             :       }
     970          28 :       gel(y,1) = gerepileupto(av,p2);
     971          28 :       gel(y,2) = const_col(lx-1, gen_1); return y;
     972             : 
     973         252 :     default:
     974             :     {
     975         252 :       GEN w = NULL, T = pol;
     976             :       long t1, t2;
     977         252 :       av = avma;
     978         252 :       RgX_type_decode(tx, &t1, &t2);
     979         252 :       if (t1 == t_COMPLEX) w = gen_I();
     980         203 :       else if (t1 == t_QUAD) w = mkquad(pol,gen_0,gen_1);
     981         252 :       if (w)
     982             :       { /* substitute I or w by t_POLMOD */
     983         126 :         T = leafcopy(pol); setvarn(T, fetch_var());
     984         126 :         x = RgX_fix_quad(x, T);
     985             :       }
     986         252 :       switch (t2)
     987             :       {
     988         161 :         case t_INT: case t_FRAC: p1 = nffactor(T,x); break;
     989          56 :         case t_INTMOD:
     990          56 :           T = RgX_to_FpX(T,p);
     991          56 :           if (FpX_is_irred(T,p)) { p1 = factmod(x,mkvec2(p,T)); break; }
     992             :         /*fall through*/
     993             :         default:
     994          56 :           if (w) (void)delete_var();
     995          56 :           pari_err_IMPL("factor for general polynomial");
     996             :           return NULL; /* LCOV_EXCL_LINE */
     997             :       }
     998         196 :       if (t1 == t_POLMOD) return gerepileupto(av, p1);
     999             :       /* substitute back I or w */
    1000          98 :       gel(p1,1) = gsubst(liftpol_shallow(gel(p1,1)), varn(T), w);
    1001          98 :       (void)delete_var(); return gerepilecopy(av, p1);
    1002             :     }
    1003             :   }
    1004             : }
    1005             : 
    1006             : static GEN
    1007       62791 : factor_domain(GEN x, GEN dom)
    1008             : {
    1009       62791 :   long tx = typ(x);
    1010       62791 :   long tdom = dom ? typ(dom): 0;
    1011             :   pari_sp av;
    1012             : 
    1013       62791 :   if (gequal0(x))
    1014          63 :     switch(tx)
    1015             :     {
    1016          63 :       case t_INT:
    1017             :       case t_COMPLEX:
    1018             :       case t_POL:
    1019          63 :       case t_RFRAC: return prime_fact(x);
    1020           0 :       default: pari_err_TYPE("factor",x);
    1021             :     }
    1022       62728 :   av = avma;
    1023       62728 :   switch(tx)
    1024             :   {
    1025        2429 :     case t_POL: return RgX_factor(x, dom);
    1026          35 :     case t_RFRAC: {
    1027          35 :       GEN a = gel(x,1), b = gel(x,2);
    1028          35 :       GEN y = famat_inv_shallow(RgX_factor(b, dom));
    1029          35 :       if (typ(a)==t_POL) y = famat_mul_shallow(RgX_factor(a, dom), y);
    1030          35 :       return gerepilecopy(av, sort_factor_pol(y, cmp_universal));
    1031             :     }
    1032       60194 :     case t_INT:  if (tdom==0 || tdom==t_INT) return Z_factor(x);
    1033          28 :     case t_FRAC: if (tdom==0 || tdom==t_INT) return Q_factor(x);
    1034             :     case t_COMPLEX: /* fall through */
    1035          49 :       if (tdom==0 || tdom==t_COMPLEX)
    1036             :       {
    1037          49 :         GEN y = gauss_factor(x); if (y) return y;
    1038             :       }
    1039             :       /* fall through */
    1040             :   }
    1041           0 :   pari_err_TYPE("factor",x);
    1042             :   return NULL; /* LCOV_EXCL_LINE */
    1043             : }
    1044             : 
    1045             : GEN
    1046       62112 : factor(GEN x) { return factor_domain(x, NULL); }
    1047             : 
    1048             : /*******************************************************************/
    1049             : /*                                                                 */
    1050             : /*                     ROOTS --> MONIC POLYNOMIAL                  */
    1051             : /*                                                                 */
    1052             : /*******************************************************************/
    1053             : static GEN
    1054     1684321 : normalized_mul(void *E, GEN x, GEN y)
    1055             : {
    1056     1684321 :   long a = gel(x,1)[1], b = gel(y,1)[1];
    1057             :   (void) E;
    1058     1684308 :   return mkvec2(mkvecsmall(a + b),
    1059     1684321 :                 RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
    1060             : }
    1061             : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
    1062             : static GEN
    1063      931554 : normalized_to_RgX(GEN L)
    1064             : {
    1065      931554 :   long i, a = gel(L,1)[1];
    1066      931554 :   GEN A = gel(L,2);
    1067      931554 :   GEN z = cgetg(a + 3, t_POL);
    1068      931554 :   z[1] = evalsigne(1) | evalvarn(varn(A));
    1069     5439245 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
    1070      934012 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
    1071      931555 :   gel(z,i) = gen_1; return z;
    1072             : }
    1073             : 
    1074             : /* compute prod (x - a[i]) */
    1075             : GEN
    1076      671978 : roots_to_pol(GEN a, long v)
    1077             : {
    1078      671978 :   pari_sp av = avma;
    1079      671978 :   long i, k, lx = lg(a);
    1080             :   GEN L;
    1081      671978 :   if (lx == 1) return pol_1(v);
    1082      671930 :   L = cgetg(lx, t_VEC);
    1083     1426080 :   for (k=1,i=1; i<lx-1; i+=2)
    1084             :   {
    1085      754150 :     GEN s = gel(a,i), t = gel(a,i+1);
    1086      754150 :     GEN x0 = gmul(s,t);
    1087      754149 :     GEN x1 = gneg(gadd(s,t));
    1088      754150 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1089             :   }
    1090     1321654 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
    1091      649724 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
    1092      671930 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1093      671930 :   return gerepileupto(av, normalized_to_RgX(L));
    1094             : }
    1095             : 
    1096             : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
    1097             : GEN
    1098      259625 : roots_to_pol_r1(GEN a, long v, long r1)
    1099             : {
    1100      259625 :   pari_sp av = avma;
    1101      259625 :   long i, k, lx = lg(a);
    1102             :   GEN L;
    1103      259625 :   if (lx == 1) return pol_1(v);
    1104      259625 :   L = cgetg(lx, t_VEC);
    1105      775223 :   for (k=1,i=1; i<r1; i+=2)
    1106             :   {
    1107      515596 :     GEN s = gel(a,i), t = gel(a,i+1);
    1108      515596 :     GEN x0 = gmul(s,t);
    1109      515597 :     GEN x1 = gneg(gadd(s,t));
    1110      515598 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1111             :   }
    1112      331499 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
    1113       71872 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
    1114      884135 :   for (i=r1+1; i<lx; i++)
    1115             :   {
    1116      624507 :     GEN s = gel(a,i);
    1117      624507 :     GEN x0 = gnorm(s);
    1118      624496 :     GEN x1 = gneg(gtrace(s));
    1119      624500 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1120             :   }
    1121      259628 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1122      259625 :   return gerepileupto(av, normalized_to_RgX(L));
    1123             : }
    1124             : 
    1125             : /*******************************************************************/
    1126             : /*                                                                 */
    1127             : /*                          FACTORBACK                             */
    1128             : /*                                                                 */
    1129             : /*******************************************************************/
    1130             : static GEN
    1131    53921832 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
    1132             : static GEN
    1133    78371063 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
    1134             : static GEN
    1135    25240512 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
    1136             : static GEN
    1137      205261 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
    1138             : 
    1139             : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
    1140             : GEN
    1141    31982723 : gen_factorback(GEN L, GEN e, void *data, GEN (*_mul)(void*,GEN,GEN),
    1142             :                GEN (*_pow)(void*,GEN,GEN), GEN (*_one)(void*))
    1143             : {
    1144    31982723 :   pari_sp av = avma;
    1145             :   long k, l, lx;
    1146             :   GEN p,x;
    1147             : 
    1148    31982723 :   if (e) /* supplied vector of exponents */
    1149     1309823 :     p = L;
    1150             :   else
    1151             :   {
    1152    30672900 :     switch(typ(L)) {
    1153     6951031 :       case t_VEC:
    1154             :       case t_COL: /* product of the L[i] */
    1155     6951031 :         if (lg(L)==1) return _one? _one(data): gen_1;
    1156     6877482 :         return gerepileupto(av, gen_product(L, data, _mul));
    1157    23721881 :       case t_MAT: /* genuine factorization */
    1158    23721881 :         l = lg(L);
    1159    23721881 :         if (l == 3) break;
    1160             :         /*fall through*/
    1161             :       default:
    1162           6 :         pari_err_TYPE("factorback [not a factorization]", L);
    1163             :     }
    1164    23721874 :     p = gel(L,1);
    1165    23721874 :     e = gel(L,2);
    1166             :   }
    1167             :   /* p = elts, e = expo */
    1168    25031697 :   lx = lg(p);
    1169             :   /* check whether e is an integral vector of correct length */
    1170    25031697 :   switch(typ(e))
    1171             :   {
    1172       94584 :     case t_VECSMALL:
    1173       94584 :       if (lx != lg(e))
    1174           0 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1175       94584 :       if (lx == 1) return _one? _one(data): gen_1;
    1176       87241 :       x = cgetg(lx,t_VEC);
    1177     1149868 :       for (l=1,k=1; k<lx; k++)
    1178     1062627 :         if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
    1179       87241 :       break;
    1180    24937111 :     case t_VEC: case t_COL:
    1181    24937111 :       if (lx != lg(e) || !RgV_is_ZV(e))
    1182           7 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1183    24937104 :       if (lx == 1) return _one? _one(data): gen_1;
    1184    24836120 :       x = cgetg(lx,t_VEC);
    1185   104368292 :       for (l=1,k=1; k<lx; k++)
    1186    79532175 :         if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
    1187    24836117 :       break;
    1188           2 :     default:
    1189           2 :       pari_err_TYPE("factorback [not an exponent vector]", e);
    1190             :       return NULL;/*LCOV_EXCL_LINE*/
    1191             :   }
    1192    24923358 :   if (l==1) return gerepileupto(av, _one? _one(data): gen_1);
    1193    24831065 :   x[0] = evaltyp(t_VEC) | _evallg(l);
    1194    24831065 :   return gerepileupto(av, gen_product(x, data, _mul));
    1195             : }
    1196             : 
    1197             : GEN
    1198     7104765 : FpV_factorback(GEN L, GEN e, GEN p)
    1199     7104765 : { return gen_factorback(L, e, (void*)p, &Fpmul, &Fppow, NULL); }
    1200             : 
    1201             : ulong
    1202       89887 : Flv_factorback(GEN L, GEN e, ulong p)
    1203             : {
    1204       89887 :   long i, l = lg(e);
    1205       89887 :   ulong r = 1UL, ri = 1UL;
    1206      437549 :   for (i = 1; i < l; i++)
    1207             :   {
    1208      347662 :     long c = e[i];
    1209      347662 :     if (!c) continue;
    1210      124383 :     if (c < 0)
    1211           0 :       ri = Fl_mul(ri, Fl_powu(L[i],-c,p), p);
    1212             :     else
    1213      124383 :       r = Fl_mul(r, Fl_powu(L[i],c,p), p);
    1214             :   }
    1215       89887 :   if (ri != 1UL) r = Fl_div(r, ri, p);
    1216       89887 :   return r;
    1217             : }
    1218             : GEN
    1219        2467 : FlxqV_factorback(GEN L, GEN e, GEN Tp, ulong p)
    1220             : {
    1221        2467 :   pari_sp av = avma;
    1222        2467 :   GEN Hi = NULL, H = NULL;
    1223        2467 :   long i, l = lg(L), v = get_Flx_var(Tp);
    1224      170338 :   for (i = 1; i < l; i++)
    1225             :   {
    1226      167805 :     GEN x, ei = gel(e,i);
    1227      167805 :     long s = signe(ei);
    1228      167805 :     if (!s) continue;
    1229      159572 :     x = Flxq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
    1230      159564 :     if (s > 0)
    1231       82484 :       H = H? Flxq_mul(H, x, Tp, p): x;
    1232             :     else
    1233       77080 :       Hi = Hi? Flxq_mul(Hi, x, Tp, p): x;
    1234             :   }
    1235        2533 :   if (!Hi)
    1236             :   {
    1237           0 :     if (!H) { set_avma(av); return mkvecsmall2(v,1); }
    1238           0 :     return gerepileuptoleaf(av, H);
    1239             :   }
    1240        2533 :   Hi = Flxq_inv(Hi, Tp, p);
    1241        2467 :   return gerepileuptoleaf(av, H? Flxq_mul(H,Hi,Tp,p): Hi);
    1242             : }
    1243             : GEN
    1244          14 : FqV_factorback(GEN L, GEN e, GEN Tp, GEN p)
    1245             : {
    1246          14 :   pari_sp av = avma;
    1247          14 :   GEN Hi = NULL, H = NULL;
    1248          14 :   long i, l = lg(L), small = typ(e) == t_VECSMALL;
    1249        1744 :   for (i = 1; i < l; i++)
    1250             :   {
    1251             :     GEN x;
    1252             :     long s;
    1253        1730 :     if (small)
    1254             :     {
    1255           0 :       s = e[i]; if (!s) continue;
    1256           0 :       x = Fq_powu(gel(L,i), labs(s), Tp, p);
    1257             :     }
    1258             :     else
    1259             :     {
    1260        1730 :       GEN ei = gel(e,i);
    1261        1730 :       s = signe(ei); if (!s) continue;
    1262        1730 :       x = Fq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
    1263             :     }
    1264        1730 :     if (s > 0)
    1265         846 :       H = H? Fq_mul(H, x, Tp, p): x;
    1266             :     else
    1267         884 :       Hi = Hi? Fq_mul(Hi, x, Tp, p): x;
    1268             :   }
    1269          14 :   if (Hi)
    1270             :   {
    1271           7 :     Hi = Fq_inv(Hi, Tp, p);
    1272           7 :     H = H? Fq_mul(H,Hi,Tp,p): Hi;
    1273             :   }
    1274           7 :   else if (!H) return gc_const(av, gen_1);
    1275          14 :   return gerepileupto(av, H);
    1276             : }
    1277             : 
    1278             : GEN
    1279    24490718 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, NULL, &mul, &powi, NULL); }
    1280             : GEN
    1281     1163779 : factorback(GEN fa) { return factorback2(fa, NULL); }
    1282             : 
    1283             : GEN
    1284          70 : vecprod(GEN v)
    1285             : {
    1286          70 :   pari_sp av = avma;
    1287          70 :   if (!is_vec_t(typ(v)))
    1288           0 :     pari_err_TYPE("vecprod", v);
    1289          70 :   if (lg(v) == 1) return gen_1;
    1290          63 :   return gerepilecopy(av, gen_product(v, NULL, mul));
    1291             : }
    1292             : 
    1293             : static int
    1294       11165 : RgX_is_irred_i(GEN x)
    1295             : {
    1296             :   GEN y, p, pol;
    1297       11165 :   long l = lg(x), pa;
    1298             : 
    1299       11165 :   if (!signe(x) || l <= 3) return 0;
    1300       11165 :   switch(RgX_type(x,&p,&pol,&pa))
    1301             :   {
    1302          21 :     case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
    1303           0 :     case t_COMPLEX: return l == 4;
    1304           0 :     case t_REAL:
    1305           0 :       if (l == 4) return 1;
    1306           0 :       if (l > 5) return 0;
    1307           0 :       return gsigne(RgX_disc(x)) > 0;
    1308             :   }
    1309       11144 :   y = RgX_factor(x, NULL);
    1310       11144 :   return (lg(gcoeff(y,1,1))==l);
    1311             : }
    1312             : static int
    1313       11165 : RgX_is_irred(GEN x)
    1314       11165 : { pari_sp av = avma; return gc_bool(av, RgX_is_irred_i(x)); }
    1315             : long
    1316       11165 : polisirreducible(GEN x)
    1317             : {
    1318       11165 :   long tx = typ(x);
    1319       11165 :   if (tx == t_POL) return RgX_is_irred(x);
    1320           0 :   if (!is_scalar_t(tx)) pari_err_TYPE("polisirreducible",x);
    1321           0 :   return 0;
    1322             : }
    1323             : 
    1324             : /*******************************************************************/
    1325             : /*                                                                 */
    1326             : /*                         GENERIC GCD                             */
    1327             : /*                                                                 */
    1328             : /*******************************************************************/
    1329             : /* x is a COMPLEX or a QUAD */
    1330             : static GEN
    1331        2555 : triv_cont_gcd(GEN x, GEN y)
    1332             : {
    1333        2555 :   pari_sp av = avma;
    1334             :   GEN c;
    1335        2555 :   if (typ(x)==t_COMPLEX)
    1336             :   {
    1337        2233 :     GEN a = gel(x,1), b = gel(x,2);
    1338        2233 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
    1339          21 :     c = ggcd(a,b);
    1340             :   }
    1341             :   else
    1342         322 :     c = ggcd(gel(x,2),gel(x,3));
    1343         343 :   return gerepileupto(av, ggcd(c,y));
    1344             : }
    1345             : 
    1346             : /* y is a PADIC, x a rational number or an INTMOD */
    1347             : static GEN
    1348        1407 : padic_gcd(GEN x, GEN y)
    1349             : {
    1350        1407 :   GEN p = gel(y,2);
    1351        1407 :   long v = gvaluation(x,p), w = valp(y);
    1352        1407 :   if (w < v) v = w;
    1353        1407 :   return powis(p, v);
    1354             : }
    1355             : 
    1356             : /* x,y in Z[i], at least one of which is t_COMPLEX */
    1357             : static GEN
    1358         399 : gauss_gcd(GEN x, GEN y)
    1359             : {
    1360         399 :   pari_sp av = avma;
    1361             :   GEN dx, dy;
    1362         399 :   dx = denom_i(x); x = gmul(x, dx);
    1363         399 :   dy = denom_i(y); y = gmul(y, dy);
    1364         847 :   while (!gequal0(y))
    1365             :   {
    1366         448 :     GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
    1367         448 :     x = y; y = z;
    1368             :   }
    1369         399 :   x = gauss_normal(x);
    1370         399 :   if (typ(x) == t_COMPLEX)
    1371             :   {
    1372         196 :     if      (gequal0(gel(x,2))) x = gel(x,1);
    1373         196 :     else if (gequal0(gel(x,1))) x = gel(x,2);
    1374             :   }
    1375         399 :   return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
    1376             : }
    1377             : 
    1378             : static int
    1379        2590 : c_is_rational(GEN x)
    1380        2590 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
    1381             : static GEN
    1382        1267 : c_zero_gcd(GEN c)
    1383             : {
    1384        1267 :   GEN x = gel(c,1), y = gel(c,2);
    1385        1267 :   long tx = typ(x), ty = typ(y);
    1386        1267 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
    1387          42 :   if (tx == t_PADIC || tx == t_INTMOD
    1388          42 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
    1389          35 :   return gauss_gcd(c, gen_0);
    1390             : }
    1391             : 
    1392             : /* gcd(x, 0) */
    1393             : static GEN
    1394     8244858 : zero_gcd(GEN x)
    1395             : {
    1396             :   pari_sp av;
    1397     8244858 :   switch(typ(x))
    1398             :   {
    1399       42322 :     case t_INT: return absi(x);
    1400       13850 :     case t_FRAC: return absfrac(x);
    1401        1267 :     case t_COMPLEX: return c_zero_gcd(x);
    1402         616 :     case t_REAL: return gen_1;
    1403         728 :     case t_PADIC: return powis(gel(x,2), valp(x));
    1404         245 :     case t_SER: return pol_xnall(valp(x), varn(x));
    1405        3172 :     case t_POLMOD: {
    1406        3172 :       GEN d = gel(x,2);
    1407        3172 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
    1408         406 :       return isinexact(d)? zero_gcd(d): gcopy(d);
    1409             :     }
    1410     7920140 :     case t_POL:
    1411     7920140 :       if (!isinexact(x)) break;
    1412          14 :       av = avma;
    1413          14 :       return gerepileupto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
    1414             : 
    1415      217067 :     case t_RFRAC:
    1416      217067 :       if (!isinexact(x)) break;
    1417           0 :       av = avma;
    1418           0 :       return gerepileupto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
    1419             :   }
    1420     8182644 :   return gcopy(x);
    1421             : }
    1422             : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
    1423             : static GEN
    1424     9357573 : zero_gcd2(GEN y, GEN z)
    1425             : {
    1426             :   pari_sp av;
    1427     9357573 :   switch(typ(z))
    1428             :   {
    1429     8229817 :     case t_INT: return zero_gcd(y);
    1430     1126706 :     case t_INTMOD:
    1431     1126706 :       av = avma;
    1432     1126706 :       return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
    1433        1050 :     case t_FFELT:
    1434        1050 :       av = avma;
    1435        1050 :       return gerepileupto(av, gmul(y, FF_1(z)));
    1436           0 :     default:
    1437           0 :       pari_err_TYPE("zero_gcd", z);
    1438             :       return NULL;/*LCOV_EXCL_LINE*/
    1439             :   }
    1440             : }
    1441             : static GEN
    1442     2786273 : cont_gcd_pol_i(GEN x, GEN y) { return scalarpol(ggcd(content(x),y), varn(x));}
    1443             : /* tx = t_POL, y considered as constant */
    1444             : static GEN
    1445     2786266 : cont_gcd_pol(GEN x, GEN y)
    1446     2786266 : { pari_sp av = avma; return gerepileupto(av, cont_gcd_pol_i(x,y)); }
    1447             : /* tx = t_RFRAC, y considered as constant */
    1448             : static GEN
    1449      864932 : cont_gcd_rfrac(GEN x, GEN y)
    1450             : {
    1451      864932 :   pari_sp av = avma;
    1452      864932 :   GEN cx; x = primitive_part(x, &cx);
    1453             :   /* e.g. Mod(1,2) / (2*y+1) => primitive_part = Mod(1,2)*y^0 */
    1454      864932 :   if (typ(x) != t_RFRAC) x = cont_gcd_pol_i(x, y);
    1455      864925 :   else x = gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2));
    1456      864932 :   return gerepileupto(av, x);
    1457             : }
    1458             : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
    1459             : static GEN
    1460        2492 : cont_gcd_gen(GEN x, GEN y)
    1461             : {
    1462        2492 :   pari_sp av = avma;
    1463        2492 :   return gerepileupto(av, ggcd(content(x),y));
    1464             : }
    1465             : /* !is_const(tx), y considered as constant */
    1466             : static GEN
    1467     2799678 : cont_gcd(GEN x, long tx, GEN y)
    1468             : {
    1469     2799678 :   switch(tx)
    1470             :   {
    1471       10934 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1472     2786252 :     case t_POL: return cont_gcd_pol(x,y);
    1473        2492 :     default: return cont_gcd_gen(x,y);
    1474             :   }
    1475             : }
    1476             : static GEN
    1477    11202039 : gcdiq(GEN x, GEN y)
    1478             : {
    1479             :   GEN z;
    1480    11202039 :   if (!signe(x)) return Q_abs(y);
    1481     4113473 :   z = cgetg(3,t_FRAC);
    1482     4113495 :   gel(z,1) = gcdii(x,gel(y,1));
    1483     4113453 :   gel(z,2) = icopy(gel(y,2));
    1484     4113478 :   return z;
    1485             : }
    1486             : static GEN
    1487    25534678 : gcdqq(GEN x, GEN y)
    1488             : {
    1489    25534678 :   GEN z = cgetg(3,t_FRAC);
    1490    25534679 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1491    25534264 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1492    25534604 :   return z;
    1493             : }
    1494             : /* assume x,y t_INT or t_FRAC */
    1495             : GEN
    1496   239025833 : Q_gcd(GEN x, GEN y)
    1497             : {
    1498   239025833 :   long tx = typ(x), ty = typ(y);
    1499   239025833 :   if (tx == t_INT)
    1500   204890426 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1501             :   else
    1502    34135407 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1503             : }
    1504             : 
    1505             : GEN
    1506    32492512 : ggcd(GEN x, GEN y)
    1507             : {
    1508    32492512 :   long vx, vy, tx = typ(x), ty = typ(y);
    1509             :   pari_sp av, tetpil;
    1510             :   GEN p1,z;
    1511             : 
    1512    64985024 :   if (is_noncalc_t(tx) || is_matvec_t(tx) ||
    1513    64985024 :       is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
    1514    32492513 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1515             :   /* tx <= ty */
    1516    32492513 :   z = gisexactzero(x); if (z) return zero_gcd2(y,z);
    1517    29413611 :   z = gisexactzero(y); if (z) return zero_gcd2(x,z);
    1518    23134940 :   if (is_const_t(tx))
    1519             :   {
    1520    15171318 :     if (ty == tx) switch(tx)
    1521             :     {
    1522     8132995 :       case t_INT:
    1523     8132995 :         return gcdii(x,y);
    1524             : 
    1525     4084710 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1526     4084710 :         if (equalii(gel(x,1),gel(y,1)))
    1527     4084703 :           gel(z,1) = icopy(gel(x,1));
    1528             :         else
    1529           7 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1530     4084710 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1531             :         else
    1532             :         {
    1533     4084710 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1534     4084710 :           if (!equali1(p1))
    1535             :           {
    1536           7 :             p1 = gcdii(p1,gel(y,2));
    1537           7 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1538           7 :             else p1 = gerepileuptoint(av, p1);
    1539             :           }
    1540     4084710 :           gel(z,2) = p1;
    1541             :         }
    1542     4084710 :         return z;
    1543             : 
    1544      158716 :       case t_FRAC:
    1545      158716 :         return gcdqq(x,y);
    1546             : 
    1547        5397 :       case t_FFELT:
    1548        5397 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1549        5397 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1550             : 
    1551          21 :       case t_COMPLEX:
    1552          21 :         if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
    1553           7 :         return triv_cont_gcd(y,x);
    1554             : 
    1555          14 :       case t_PADIC:
    1556          14 :         if (!equalii(gel(x,2),gel(y,2))) return gen_1;
    1557           7 :         return powis(gel(y,2), minss(valp(x), valp(y)));
    1558             : 
    1559         133 :       case t_QUAD:
    1560         133 :         av=avma; p1=gdiv(x,y);
    1561         133 :         if (gequal0(gel(p1,3)))
    1562             :         {
    1563          14 :           p1=gel(p1,2);
    1564          14 :           if (typ(p1)==t_INT) { set_avma(av); return gcopy(y); }
    1565           7 :           tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
    1566             :         }
    1567         119 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {set_avma(av); return gcopy(y);}
    1568         112 :         p1 = ginv(p1); set_avma(av);
    1569         112 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1570         105 :         return triv_cont_gcd(y,x);
    1571             : 
    1572           0 :       default: return gen_1; /* t_REAL */
    1573             :     }
    1574     2789332 :     if (is_const_t(ty)) switch(tx)
    1575             :     {
    1576       16786 :       case t_INT:
    1577       16786 :         switch(ty)
    1578             :         {
    1579          70 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1580          70 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1581          70 :             p1 = gcdii(gel(y,1),gel(y,2));
    1582          70 :             if (!equali1(p1)) {
    1583          14 :               p1 = gcdii(x,p1);
    1584          14 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1585             :               else
    1586          14 :                 p1 = gerepileuptoint(av, p1);
    1587             :             }
    1588          70 :             gel(z,2) = p1; return z;
    1589             : 
    1590        8162 :           case t_REAL: return gen_1;
    1591             : 
    1592        4452 :           case t_FRAC:
    1593        4452 :             return gcdiq(x,y);
    1594             : 
    1595        2471 :           case t_COMPLEX:
    1596        2471 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1597        2128 :             return triv_cont_gcd(y,x);
    1598             : 
    1599         112 :           case t_FFELT:
    1600         112 :             if (!FF_equal0(y)) return FF_1(y);
    1601           0 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1602             : 
    1603        1393 :           case t_PADIC:
    1604        1393 :             return padic_gcd(x,y);
    1605             : 
    1606         126 :           case t_QUAD:
    1607         126 :             return triv_cont_gcd(y,x);
    1608           0 :           default:
    1609           0 :             pari_err_TYPE2("gcd",x,y);
    1610             :         }
    1611             : 
    1612          14 :       case t_REAL:
    1613          14 :         switch(ty)
    1614             :         {
    1615          14 :           case t_INTMOD:
    1616             :           case t_FFELT:
    1617          14 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1618           0 :           default: return gen_1;
    1619             :         }
    1620             : 
    1621          49 :       case t_INTMOD:
    1622          49 :         switch(ty)
    1623             :         {
    1624          14 :           case t_FRAC:
    1625          14 :             av = avma; p1=gcdii(gel(x,1),gel(y,2)); set_avma(av);
    1626          14 :             if (!equali1(p1)) pari_err_OP("gcd",x,y);
    1627           7 :             return ggcd(gel(y,1), x);
    1628             : 
    1629          14 :           case t_FFELT:
    1630             :           {
    1631          14 :             GEN p = gel(y,4);
    1632          14 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1633           7 :             if (!FF_equal0(y)) return FF_1(y);
    1634           0 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1635             :           }
    1636             : 
    1637          14 :           case t_COMPLEX: case t_QUAD:
    1638          14 :             return triv_cont_gcd(y,x);
    1639             : 
    1640           7 :           case t_PADIC:
    1641           7 :             return padic_gcd(x,y);
    1642             : 
    1643           0 :           default: pari_err_TYPE2("gcd",x,y);
    1644             :         }
    1645             : 
    1646         210 :       case t_FRAC:
    1647         210 :         switch(ty)
    1648             :         {
    1649          84 :           case t_COMPLEX:
    1650          84 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1651             :           case t_QUAD:
    1652         154 :             return triv_cont_gcd(y,x);
    1653          42 :           case t_FFELT:
    1654             :           {
    1655          42 :             GEN p = gel(y,4);
    1656          42 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1657          21 :             if (!FF_equal0(y)) return FF_1(y);
    1658           0 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1659             :           }
    1660             : 
    1661           7 :           case t_PADIC:
    1662           7 :             return padic_gcd(x,y);
    1663             : 
    1664           0 :           default: pari_err_TYPE2("gcd",x,y);
    1665             :         }
    1666          70 :       case t_FFELT:
    1667          70 :         switch(ty)
    1668             :         {
    1669          42 :           case t_PADIC:
    1670             :           {
    1671          42 :             GEN p = gel(y,2);
    1672          42 :             long v = valp(y);
    1673          42 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1674          14 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1675             :           }
    1676          28 :           default: pari_err_TYPE2("gcd",x,y);
    1677             :         }
    1678             : 
    1679          14 :       case t_COMPLEX:
    1680          14 :         switch(ty)
    1681             :         {
    1682          14 :           case t_PADIC:
    1683          14 :           case t_QUAD: return triv_cont_gcd(x,y);
    1684           0 :           default: pari_err_TYPE2("gcd",x,y);
    1685             :         }
    1686             : 
    1687           7 :       case t_PADIC:
    1688           7 :         switch(ty)
    1689             :         {
    1690           7 :           case t_QUAD: return triv_cont_gcd(y,x);
    1691           0 :           default: pari_err_TYPE2("gcd",x,y);
    1692             :         }
    1693             : 
    1694           0 :       default: return gen_1; /* tx = t_REAL */
    1695             :     }
    1696     2772182 :     return cont_gcd(y,ty, x);
    1697             :   }
    1698             : 
    1699     7963622 :   if (tx == t_POLMOD)
    1700             :   {
    1701        6119 :     if (ty == t_POLMOD)
    1702             :     {
    1703        6063 :       GEN T = gel(x,1);
    1704        6063 :       z = cgetg(3,t_POLMOD);
    1705        6063 :       T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
    1706        6063 :       gel(z,1) = T;
    1707        6063 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1708             :       else
    1709             :       {
    1710             :         GEN X, Y, d;
    1711        6063 :         av = avma; X = gel(x,2); Y = gel(y,2);
    1712        6063 :         d = ggcd(content(X), content(Y));
    1713        6063 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1714        6063 :         p1 = ggcd(T, X);
    1715        6063 :         gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
    1716             :       }
    1717        6063 :       return z;
    1718             :     }
    1719          56 :     vx = varn(gel(x,1));
    1720          56 :     switch(ty)
    1721             :     {
    1722          28 :       case t_POL:
    1723          28 :         vy = varn(y);
    1724          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1725          14 :         z = cgetg(3,t_POLMOD);
    1726          14 :         gel(z,1) = RgX_copy(gel(x,1));
    1727          14 :         av = avma; p1 = ggcd(gel(x,1),gel(x,2));
    1728          14 :         gel(z,2) = gerepileupto(av, ggcd(p1,y));
    1729          14 :         return z;
    1730             : 
    1731          28 :       case t_RFRAC:
    1732          28 :         vy = varn(gel(y,2));
    1733          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1734          28 :         av = avma;
    1735          28 :         p1 = ggcd(gel(x,1),gel(y,2));
    1736          28 :         if (degpol(p1)) pari_err_OP("gcd",x,y);
    1737          21 :         set_avma(av); return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1738             :     }
    1739             :   }
    1740             : 
    1741     7957503 :   vx = gvar(x);
    1742     7957503 :   vy = gvar(y);
    1743     7957503 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1744     7945484 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1745             : 
    1746             :   /* vx = vy: same main variable */
    1747     7930007 :   switch(tx)
    1748             :   {
    1749     7769073 :     case t_POL:
    1750             :       switch(ty)
    1751             :       {
    1752     6915054 :         case t_POL: return RgX_gcd(x,y);
    1753          21 :         case t_SER:
    1754          21 :           z = ggcd(content(x), content(y));
    1755          21 :           return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
    1756      853998 :         case t_RFRAC: return cont_gcd_rfrac(y, x);
    1757             :       }
    1758           0 :       break;
    1759             : 
    1760          14 :     case t_SER:
    1761          14 :       z = ggcd(content(x), content(y));
    1762             :       switch(ty)
    1763             :       {
    1764           7 :         case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
    1765           7 :         case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
    1766             :       }
    1767           0 :       break;
    1768             : 
    1769      160920 :     case t_RFRAC:
    1770             :     {
    1771      160920 :       GEN xd = gel(x,2), yd = gel(y,2);
    1772      160920 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1773      160920 :       z = cgetg(3,t_RFRAC); av = avma;
    1774      160920 :       gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1775      160920 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1776             :     }
    1777             :   }
    1778           0 :   pari_err_TYPE2("gcd",x,y);
    1779             :   return NULL; /* LCOV_EXCL_LINE */
    1780             : }
    1781             : GEN
    1782     5203504 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
    1783             : 
    1784             : static GEN
    1785        3568 : fix_lcm(GEN x)
    1786             : {
    1787             :   GEN t;
    1788        3568 :   switch(typ(x))
    1789             :   {
    1790        3470 :     case t_INT: if (signe(x)<0) x = negi(x);
    1791        3470 :       break;
    1792          98 :     case t_POL:
    1793          98 :       if (lg(x) <= 2) break;
    1794          98 :       t = leading_coeff(x);
    1795          98 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1796             :   }
    1797        3568 :   return x;
    1798             : }
    1799             : GEN
    1800        2877 : glcm0(GEN x, GEN y)
    1801             : {
    1802        2877 :   if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
    1803        2828 :   return glcm(x,y);
    1804             : }
    1805             : GEN
    1806        3470 : ZV_lcm(GEN x) { return fix_lcm(gassoc_proto(lcmii,x,NULL)); }
    1807             : GEN
    1808        3171 : glcm(GEN x, GEN y)
    1809             : {
    1810             :   pari_sp av;
    1811             :   GEN z;
    1812        3171 :   if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
    1813          63 :   av = avma; z = ggcd(x,y);
    1814          63 :   if (!gequal1(z))
    1815             :   {
    1816          63 :     if (gequal0(z)) { set_avma(av); return gmul(x,y); }
    1817          49 :     y = gdiv(y,z);
    1818             :   }
    1819          49 :   return gerepileupto(av, fix_lcm(gmul(x,y)));
    1820             : }
    1821             : 
    1822             : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
    1823             : static int
    1824           0 : pol_approx0(GEN r, GEN x, int exact)
    1825             : {
    1826             :   long i, l;
    1827           0 :   if (exact) return !signe(r);
    1828           0 :   l = minss(lg(x), lg(r));
    1829           0 :   for (i = 2; i < l; i++)
    1830           0 :     if (!cx_approx0(gel(r,i), gel(x,i))) return 0;
    1831           0 :   return 1;
    1832             : }
    1833             : 
    1834             : GEN
    1835           0 : RgX_gcd_simple(GEN x, GEN y)
    1836             : {
    1837           0 :   pari_sp av1, av = avma;
    1838           0 :   GEN r, yorig = y;
    1839           0 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1840             : 
    1841             :   for(;;)
    1842             :   {
    1843           0 :     av1 = avma; r = RgX_rem(x,y);
    1844           0 :     if (pol_approx0(r, x, exact))
    1845             :     {
    1846           0 :       set_avma(av1);
    1847           0 :       if (y == yorig) return RgX_copy(y);
    1848           0 :       y = normalizepol_approx(y, lg(y));
    1849           0 :       if (lg(y) == 3) { set_avma(av); return pol_1(varn(x)); }
    1850           0 :       return gerepileupto(av,y);
    1851             :     }
    1852           0 :     x = y; y = r;
    1853           0 :     if (gc_needed(av,1)) {
    1854           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1855           0 :       gerepileall(av,2, &x,&y);
    1856             :     }
    1857             :   }
    1858             : }
    1859             : GEN
    1860           0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
    1861             : {
    1862           0 :   pari_sp av = avma;
    1863             :   GEN q, r, d, d1, u, v, v1;
    1864           0 :   int exact = !(isinexactreal(a) || isinexactreal(b));
    1865             : 
    1866           0 :   d = a; d1 = b; v = gen_0; v1 = gen_1;
    1867             :   for(;;)
    1868             :   {
    1869           0 :     if (pol_approx0(d1, a, exact)) break;
    1870           0 :     q = poldivrem(d,d1, &r);
    1871           0 :     v = gsub(v, gmul(q,v1));
    1872           0 :     u=v; v=v1; v1=u;
    1873           0 :     u=r; d=d1; d1=u;
    1874             :   }
    1875           0 :   u = gsub(d, gmul(b,v));
    1876           0 :   u = RgX_div(u,a);
    1877             : 
    1878           0 :   gerepileall(av, 3, &u,&v,&d);
    1879           0 :   *pu = u;
    1880           0 :   *pv = v; return d;
    1881             : }
    1882             : 
    1883             : GEN
    1884          42 : ghalfgcd(GEN x, GEN y)
    1885             : {
    1886          42 :   long tx = typ(x), ty = typ(y);
    1887          42 :   if (tx==t_INT && ty==t_INT) return halfgcdii(x, y);
    1888          35 :   if (tx==t_POL && ty==t_POL && varn(x)==varn(y)) return RgX_halfgcd(x, y);
    1889           0 :   pari_err_OP("halfgcd", x, y);
    1890             :   return NULL; /* LCOV_EXCL_LINE */
    1891             : }
    1892             : 
    1893             : /*******************************************************************/
    1894             : /*                                                                 */
    1895             : /*                    CONTENT / PRIMITIVE PART                     */
    1896             : /*                                                                 */
    1897             : /*******************************************************************/
    1898             : 
    1899             : GEN
    1900    74074595 : content(GEN x)
    1901             : {
    1902    74074595 :   long lx, i, t, tx = typ(x);
    1903    74074595 :   pari_sp av = avma;
    1904             :   GEN c;
    1905             : 
    1906    74074595 :   if (is_scalar_t(tx)) return zero_gcd(x);
    1907    74062355 :   switch(tx)
    1908             :   {
    1909      864967 :     case t_RFRAC:
    1910             :     {
    1911      864967 :       GEN n = gel(x,1), d = gel(x,2);
    1912             :       /* -- varncmp(vn, vd) < 0 can't happen
    1913             :        * -- if n is POLMOD, its main variable (in the sense of gvar2)
    1914             :        *    has lower priority than denominator */
    1915      864967 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    1916      824622 :         n = isinexact(n)? zero_gcd(n): gcopy(n);
    1917             :       else
    1918       40345 :         n = content(n);
    1919      864967 :       return gerepileupto(av, gdiv(n, content(d)));
    1920             :     }
    1921             : 
    1922     2568326 :     case t_VEC: case t_COL:
    1923     2568326 :       lx = lg(x); if (lx==1) return gen_0;
    1924     2568319 :       break;
    1925             : 
    1926        1652 :     case t_MAT:
    1927             :     {
    1928             :       long hx, j;
    1929        1652 :       lx = lg(x);
    1930        1652 :       if (lx == 1) return gen_0;
    1931        1645 :       hx = lgcols(x);
    1932        1645 :       if (hx == 1) return gen_0;
    1933        1638 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    1934        1638 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    1935        1638 :       c = content(gel(x,1));
    1936        3276 :       for (j=2; j<lx; j++)
    1937        4914 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    1938        1638 :       if (typ(c) == t_INTMOD || isinexact(c)) return gc_const(av, gen_1);
    1939        1638 :       return gerepileupto(av,c);
    1940             :     }
    1941             : 
    1942    70627207 :     case t_POL: case t_SER:
    1943    70627207 :       lx = lg(x); if (lx == 2) return gen_0;
    1944    70604023 :       break;
    1945          21 :     case t_VECSMALL: return utoi(zv_content(x));
    1946         182 :     case t_QFB:
    1947         182 :       lx = 4; break;
    1948             : 
    1949           0 :     default: pari_err_TYPE("content",x);
    1950             :       return NULL; /* LCOV_EXCL_LINE */
    1951             :   }
    1952   215793778 :   for (i=lontyp[tx]; i<lx; i++)
    1953   155812306 :     if (typ(gel(x,i)) != t_INT) break;
    1954    73172525 :   lx--; c = gel(x,lx);
    1955    73172525 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    1956    73172528 :   if (i > lx)
    1957             :   { /* integer coeffs */
    1958    62875072 :     while (lx-- > lontyp[tx])
    1959             :     {
    1960    61381162 :       c = gcdii(c, gel(x,lx));
    1961    61381129 :       if (equali1(c)) return gc_const(av, gen_1);
    1962             :     }
    1963             :   }
    1964             :   else
    1965             :   {
    1966    13191053 :     if (isinexact(c)) c = zero_gcd(c);
    1967    36639912 :     while (lx-- > lontyp[tx])
    1968             :     {
    1969    23448859 :       GEN d = gel(x,lx);
    1970    23448859 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    1971    23448859 :       c = ggcd(c, d);
    1972             :     }
    1973    13191053 :     if (isinexact(c)) return gc_const(av, gen_1);
    1974             :   }
    1975    14684963 :   switch(typ(c))
    1976             :   {
    1977     1498825 :     case t_INT:
    1978     1498825 :       if (signe(c) < 0) c = negi(c);
    1979     1498825 :       break;
    1980           0 :     case t_VEC: case t_COL: case t_MAT:
    1981           0 :       pari_err_TYPE("content",x);
    1982             :   }
    1983             : 
    1984    14684977 :   return av==avma? gcopy(c): gerepileupto(av,c);
    1985             : }
    1986             : 
    1987             : GEN
    1988     1700922 : primitive_part(GEN x, GEN *ptc)
    1989             : {
    1990     1700922 :   pari_sp av = avma;
    1991     1700922 :   GEN c = content(x);
    1992     1700908 :   if (gequal1(c)) { set_avma(av); c = NULL; }
    1993       99807 :   else if (!gequal0(c)) x = gdiv(x,c);
    1994     1700913 :   if (ptc) *ptc = c;
    1995     1700913 :   return x;
    1996             : }
    1997             : GEN
    1998        2233 : primpart(GEN x) { return primitive_part(x, NULL); }
    1999             : 
    2000             : static GEN
    2001    71608248 : Q_content_v(GEN x, long imin, long l)
    2002             : {
    2003    71608248 :   pari_sp av = avma;
    2004    71608248 :   long i = l-1;
    2005    71608248 :   GEN d = Q_content_safe(gel(x,i));
    2006    71606833 :   if (!d) return NULL;
    2007   310577513 :   for (i--; i >= imin; i--)
    2008             :   {
    2009   239001496 :     GEN c = Q_content_safe(gel(x,i));
    2010   239012440 :     if (!c) return NULL;
    2011   239012426 :     d = Q_gcd(d, c);
    2012   238972050 :     if (gc_needed(av,1)) d = gerepileupto(av, d);
    2013             :   }
    2014    71576017 :   return gerepileupto(av, d);
    2015             : }
    2016             : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    2017             :  * of Q(x2,...,xn)[x1] */
    2018             : GEN
    2019   359467160 : Q_content_safe(GEN x)
    2020             : {
    2021             :   long l;
    2022   359467160 :   switch(typ(x))
    2023             :   {
    2024   253822740 :     case t_INT:  return absi(x);
    2025    33686428 :     case t_FRAC: return absfrac(x);
    2026    35839147 :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2027    35839147 :       l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
    2028    36100585 :     case t_POL:
    2029    36100585 :       l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
    2030       32683 :     case t_POLMOD: return Q_content_safe(gel(x,2));
    2031          70 :     case t_RFRAC:
    2032             :     {
    2033             :       GEN a, b;
    2034          70 :       a = Q_content_safe(gel(x,1)); if (!a) return NULL;
    2035          63 :       b = Q_content_safe(gel(x,2)); if (!b) return NULL;
    2036          63 :       return gdiv(a, b);
    2037             :     }
    2038             :   }
    2039         366 :   return NULL;
    2040             : }
    2041             : GEN
    2042      308318 : Q_content(GEN x)
    2043             : {
    2044      308318 :   GEN c = Q_content_safe(x);
    2045      308319 :   if (!c) pari_err_TYPE("Q_content",x);
    2046      308319 :   return c;
    2047             : }
    2048             : 
    2049             : GEN
    2050       13146 : ZX_content(GEN x)
    2051             : {
    2052       13146 :   long i, l = lg(x);
    2053             :   GEN d;
    2054             :   pari_sp av;
    2055             : 
    2056       13146 :   if (l == 2) return gen_0;
    2057       13146 :   d = gel(x,2);
    2058       13146 :   if (l == 3) return absi(d);
    2059        9170 :   av = avma;
    2060       18977 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    2061        9170 :   if (signe(d) < 0) d = negi(d);
    2062        9170 :   return gerepileuptoint(av, d);
    2063             : }
    2064             : 
    2065             : static GEN
    2066       43169 : Z_content_v(GEN x, long i, long l)
    2067             : {
    2068       43169 :   pari_sp av = avma;
    2069       43169 :   GEN d = Z_content(gel(x,i));
    2070       43169 :   if (!d) return NULL;
    2071      273098 :   for (i++; i<l; i++)
    2072             :   {
    2073      235060 :     GEN c = Z_content(gel(x,i));
    2074      235060 :     if (!c) return NULL;
    2075      234878 :     d = gcdii(d, c); if (equali1(d)) return NULL;
    2076      234752 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    2077             :   }
    2078       38038 :   return gerepileuptoint(av, d);
    2079             : }
    2080             : /* return NULL for 1 */
    2081             : GEN
    2082      280049 : Z_content(GEN x)
    2083             : {
    2084             :   long l;
    2085      280049 :   switch(typ(x))
    2086             :   {
    2087      217875 :     case t_INT:
    2088      217875 :       if (is_pm1(x)) return NULL;
    2089      216692 :       return absi(x);
    2090        9002 :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2091        9002 :       l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
    2092       53172 :     case t_POL:
    2093       53172 :       l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
    2094           0 :     case t_POLMOD: return Z_content(gel(x,2));
    2095             :   }
    2096           0 :   pari_err_TYPE("Z_content", x);
    2097             :   return NULL; /* LCOV_EXCL_LINE */
    2098             : }
    2099             : 
    2100             : static GEN
    2101    31737440 : Q_denom_v(GEN x, long i, long l)
    2102             : {
    2103    31737440 :   pari_sp av = avma;
    2104    31737440 :   GEN d = Q_denom_safe(gel(x,i));
    2105    31737283 :   if (!d) return NULL;
    2106   123542005 :   for (i++; i<l; i++)
    2107             :   {
    2108    91804885 :     GEN D = Q_denom_safe(gel(x,i));
    2109    91804651 :     if (!D) return NULL;
    2110    91804651 :     if (D != gen_1) d = lcmii(d, D);
    2111    91804563 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    2112             :   }
    2113    31737120 :   return gerepileuptoint(av, d);
    2114             : }
    2115             : /* NOT MEMORY CLEAN (because of t_FRAC).
    2116             :  * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    2117             :  * of Q(x2,...,xn)[x1] */
    2118             : GEN
    2119   161606191 : Q_denom_safe(GEN x)
    2120             : {
    2121             :   long l;
    2122   161606191 :   switch(typ(x))
    2123             :   {
    2124   109233985 :     case t_INT: return gen_1;
    2125          28 :     case t_PADIC: l = valp(x); return l < 0? powiu(gel(x,2), -l): gen_1;
    2126    20447020 :     case t_FRAC: return gel(x,2);
    2127         504 :     case t_QUAD: return Q_denom_v(x, 2, 4);
    2128    21426506 :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2129    21426506 :       l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
    2130    10401308 :     case t_POL: case t_SER:
    2131    10401308 :       l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
    2132       93282 :     case t_POLMOD: return Q_denom(gel(x,2));
    2133        8134 :     case t_RFRAC:
    2134             :     {
    2135             :       GEN a, b;
    2136        8134 :       a = Q_content(gel(x,1)); if (!a) return NULL;
    2137        8134 :       b = Q_content(gel(x,2)); if (!b) return NULL;
    2138        8134 :       return Q_denom(gdiv(a, b));
    2139             :     }
    2140             :   }
    2141          66 :   return NULL;
    2142             : }
    2143             : GEN
    2144     2160891 : Q_denom(GEN x)
    2145             : {
    2146     2160891 :   GEN d = Q_denom_safe(x);
    2147     2160882 :   if (!d) pari_err_TYPE("Q_denom",x);
    2148     2160882 :   return d;
    2149             : }
    2150             : 
    2151             : GEN
    2152    35906154 : Q_remove_denom(GEN x, GEN *ptd)
    2153             : {
    2154    35906154 :   GEN d = Q_denom_safe(x);
    2155    35905642 :   if (d) { if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d); }
    2156    35905592 :   if (ptd) *ptd = d;
    2157    35905592 :   return x;
    2158             : }
    2159             : 
    2160             : /* return y = x * d, assuming x rational, and d,y integral */
    2161             : GEN
    2162   106140852 : Q_muli_to_int(GEN x, GEN d)
    2163             : {
    2164             :   long i, l;
    2165             :   GEN y, xn, xd;
    2166             :   pari_sp av;
    2167             : 
    2168   106140852 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    2169   106145830 :   switch (typ(x))
    2170             :   {
    2171    34458058 :     case t_INT:
    2172    34458058 :       return mulii(x,d);
    2173             : 
    2174    50555480 :     case t_FRAC:
    2175    50555480 :       xn = gel(x,1);
    2176    50555480 :       xd = gel(x,2); av = avma;
    2177    50555480 :       y = mulii(xn, diviiexact(d, xd));
    2178    50546708 :       return gerepileuptoint(av, y);
    2179           7 :     case t_COMPLEX:
    2180           7 :       y = cgetg(3,t_COMPLEX);
    2181           7 :       gel(y,1) = Q_muli_to_int(gel(x,1),d);
    2182           7 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2183           7 :       return y;
    2184          14 :     case t_PADIC:
    2185          14 :       y = gcopy(x); if (!isint1(d)) setvalp(y, 0);
    2186          14 :       return y;
    2187         175 :     case t_QUAD:
    2188         175 :       y = cgetg(4,t_QUAD);
    2189         175 :       gel(y,1) = ZX_copy(gel(x,1));
    2190         175 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2191         175 :       gel(y,3) = Q_muli_to_int(gel(x,3),d); return y;
    2192             : 
    2193    12222528 :     case t_VEC: case t_COL: case t_MAT:
    2194    12222528 :       y = cgetg_copy(x, &l);
    2195    70559820 :       for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2196    12221941 :       return y;
    2197             : 
    2198     8859315 :     case t_POL: case t_SER:
    2199     8859315 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2200    41950708 :       for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2201     8858117 :       return y;
    2202             : 
    2203       50232 :     case t_POLMOD:
    2204       50232 :       retmkpolmod(Q_muli_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2205          21 :     case t_RFRAC:
    2206          21 :       return gmul(x, d);
    2207             :   }
    2208           0 :   pari_err_TYPE("Q_muli_to_int",x);
    2209             :   return NULL; /* LCOV_EXCL_LINE */
    2210             : }
    2211             : 
    2212             : static void
    2213     5866092 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
    2214             : {
    2215             :   long e, i;
    2216     5866092 :   switch(typ(c))
    2217             :   {
    2218     4214257 :     case t_REAL:
    2219     4214257 :       *exact = 0;
    2220     4214257 :       if (!signe(c)) return;
    2221     4187362 :       e = expo(c) + 1 - bit_prec(c);
    2222     4215344 :       for (i = lg(c)-1; i > 2; i--, e += BITS_IN_LONG)
    2223     3976083 :         if (c[i]) break;
    2224     4187363 :       e += vals(c[i]); break; /* e[2] != 0 */
    2225     1647315 :     case t_INT:
    2226     1647315 :       if (!signe(c)) return;
    2227      726583 :       e = expi(c);
    2228      726583 :       break;
    2229        4494 :     case t_FRAC:
    2230        4494 :       e = expi(gel(c,1)) - expi(gel(c,2));
    2231        4494 :       if (*exact) *D = lcmii(*D, gel(c,2));
    2232        4494 :       break;
    2233          26 :     default:
    2234          26 :       pari_err_TYPE("rescale_to_int",c);
    2235             :       return; /* LCOV_EXCL_LINE */
    2236             :   }
    2237     4918438 :   if (e < *emin) *emin = e;
    2238             : }
    2239             : GEN
    2240      561797 : RgM_rescale_to_int(GEN x)
    2241             : {
    2242      561797 :   long lx = lg(x), i,j, hx, emin;
    2243             :   GEN D;
    2244             :   int exact;
    2245             : 
    2246      561797 :   if (lx == 1) return cgetg(1,t_MAT);
    2247      561797 :   hx = lgcols(x);
    2248      561798 :   exact = 1;
    2249      561798 :   emin = HIGHEXPOBIT;
    2250      561798 :   D = gen_1;
    2251     2168066 :   for (j = 1; j < lx; j++)
    2252     7280648 :     for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
    2253      561769 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2254      561699 :   return grndtoi(gmul2n(x, -emin), NULL);
    2255             : }
    2256             : GEN
    2257       37387 : RgX_rescale_to_int(GEN x)
    2258             : {
    2259       37387 :   long lx = lg(x), i, emin;
    2260             :   GEN D;
    2261             :   int exact;
    2262       37387 :   if (lx == 2) return gcopy(x); /* rare */
    2263       37387 :   exact = 1;
    2264       37387 :   emin = HIGHEXPOBIT;
    2265       37387 :   D = gen_1;
    2266      229100 :   for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
    2267       37387 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2268       36281 :   return grndtoi(gmul2n(x, -emin), NULL);
    2269             : }
    2270             : 
    2271             : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    2272             : static GEN
    2273     9532643 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    2274             : {
    2275             :   long i, l;
    2276             :   GEN y, xn, xd;
    2277             :   pari_sp av;
    2278             : 
    2279     9532643 :   switch(typ(x))
    2280             :   {
    2281     2781146 :     case t_INT:
    2282     2781146 :       av = avma; y = diviiexact(x,d);
    2283     2781146 :       return gerepileuptoint(av, mulii(y,n));
    2284             : 
    2285     4586763 :     case t_FRAC:
    2286     4586763 :       xn = gel(x,1);
    2287     4586763 :       xd = gel(x,2); av = avma;
    2288     4586763 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    2289     4586763 :       return gerepileuptoint(av, y);
    2290             : 
    2291      453057 :     case t_VEC: case t_COL: case t_MAT:
    2292      453057 :       y = cgetg_copy(x, &l);
    2293     4095019 :       for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2294      453057 :       return y;
    2295             : 
    2296     1711677 :     case t_POL:
    2297     1711677 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2298     5808911 :       for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2299     1711677 :       return y;
    2300             : 
    2301           0 :     case t_RFRAC:
    2302           0 :       av = avma;
    2303           0 :       return gerepileupto(av, gmul(x,mkfrac(n,d)));
    2304             : 
    2305           0 :     case t_POLMOD:
    2306           0 :       retmkpolmod(Q_divmuli_to_int(gel(x,2), d,n), RgX_copy(gel(x,1)));
    2307             :   }
    2308           0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    2309             :   return NULL; /* LCOV_EXCL_LINE */
    2310             : }
    2311             : 
    2312             : /* return x / d. x: rational; d,result: integral. */
    2313             : static GEN
    2314    38174797 : Q_divi_to_int(GEN x, GEN d)
    2315             : {
    2316             :   long i, l;
    2317             :   GEN y;
    2318             : 
    2319    38174797 :   switch(typ(x))
    2320             :   {
    2321    30440551 :     case t_INT:
    2322    30440551 :       return diviiexact(x,d);
    2323             : 
    2324     4044844 :     case t_VEC: case t_COL: case t_MAT:
    2325     4044844 :       y = cgetg_copy(x, &l);
    2326    23393518 :       for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2327     4043675 :       return y;
    2328             : 
    2329     3684820 :     case t_POL:
    2330     3684820 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2331    16868963 :       for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2332     3684771 :       return y;
    2333             : 
    2334           7 :     case t_RFRAC:
    2335           7 :       return gdiv(x,d);
    2336             : 
    2337        5915 :     case t_POLMOD:
    2338        5915 :       retmkpolmod(Q_divi_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2339             :   }
    2340           0 :   pari_err_TYPE("Q_divi_to_int",x);
    2341             :   return NULL; /* LCOV_EXCL_LINE */
    2342             : }
    2343             : /* c t_FRAC */
    2344             : static GEN
    2345     8310854 : Q_divq_to_int(GEN x, GEN c)
    2346             : {
    2347     8310854 :   GEN n = gel(c,1), d = gel(c,2);
    2348     8310854 :   if (is_pm1(n)) {
    2349     6517404 :     GEN y = Q_muli_to_int(x,d);
    2350     6517382 :     if (signe(n) < 0) y = gneg(y);
    2351     6517382 :     return y;
    2352             :   }
    2353     1793447 :   return Q_divmuli_to_int(x, n,d);
    2354             : }
    2355             : 
    2356             : /* return y = x / c, assuming x,c rational, and y integral */
    2357             : GEN
    2358        4746 : Q_div_to_int(GEN x, GEN c)
    2359             : {
    2360        4746 :   switch(typ(c))
    2361             :   {
    2362        3283 :     case t_INT:  return Q_divi_to_int(x, c);
    2363        1463 :     case t_FRAC: return Q_divq_to_int(x, c);
    2364             :   }
    2365           0 :   pari_err_TYPE("Q_div_to_int",c);
    2366             :   return NULL; /* LCOV_EXCL_LINE */
    2367             : }
    2368             : /* return y = x * c, assuming x,c rational, and y integral */
    2369             : GEN
    2370           0 : Q_mul_to_int(GEN x, GEN c)
    2371             : {
    2372             :   GEN d, n;
    2373           0 :   switch(typ(c))
    2374             :   {
    2375           0 :     case t_INT: return Q_muli_to_int(x, c);
    2376           0 :     case t_FRAC:
    2377           0 :       n = gel(c,1);
    2378           0 :       d = gel(c,2);
    2379           0 :       return Q_divmuli_to_int(x, d,n);
    2380             :   }
    2381           0 :   pari_err_TYPE("Q_mul_to_int",c);
    2382             :   return NULL; /* LCOV_EXCL_LINE */
    2383             : }
    2384             : 
    2385             : GEN
    2386    48512047 : Q_primitive_part(GEN x, GEN *ptc)
    2387             : {
    2388    48512047 :   pari_sp av = avma;
    2389    48512047 :   GEN c = Q_content_safe(x);
    2390    48510386 :   if (c)
    2391             :   {
    2392    48510234 :     if (typ(c) == t_INT)
    2393             :     {
    2394    40200842 :       if (equali1(c)) { set_avma(av); c = NULL; }
    2395     6061636 :       else if (signe(c)) x = Q_divi_to_int(x, c);
    2396             :     }
    2397     8309392 :     else x = Q_divq_to_int(x, c);
    2398             :   }
    2399    48509237 :   if (ptc) *ptc = c;
    2400    48509237 :   return x;
    2401             : }
    2402             : GEN
    2403     5006522 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    2404             : GEN
    2405      110801 : vec_Q_primpart(GEN x)
    2406      626427 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
    2407             : GEN
    2408       16128 : row_Q_primpart(GEN M)
    2409       16128 : { return shallowtrans(vec_Q_primpart(shallowtrans(M))); }
    2410             : 
    2411             : /*******************************************************************/
    2412             : /*                                                                 */
    2413             : /*                           SUBRESULTANT                          */
    2414             : /*                                                                 */
    2415             : /*******************************************************************/
    2416             : /* for internal use */
    2417             : GEN
    2418    20027971 : gdivexact(GEN x, GEN y)
    2419             : {
    2420             :   long i,lx;
    2421             :   GEN z;
    2422    20027971 :   if (gequal1(y)) return x;
    2423    20060893 :   if (typ(y) == t_POLMOD) return gmul(x, ginv(y));
    2424    20060816 :   switch(typ(x))
    2425             :   {
    2426    17187252 :     case t_INT:
    2427    17187252 :       if (typ(y)==t_INT) return diviiexact(x,y);
    2428          75 :       if (!signe(x)) return gen_0;
    2429          40 :       break;
    2430        8316 :     case t_INTMOD:
    2431             :     case t_FFELT:
    2432        8316 :     case t_POLMOD: return gmul(x,ginv(y));
    2433     2867736 :     case t_POL:
    2434     2867736 :       switch(typ(y))
    2435             :       {
    2436         714 :         case t_INTMOD:
    2437             :         case t_FFELT:
    2438         714 :         case t_POLMOD: return gmul(x,ginv(y));
    2439       36386 :         case t_POL: { /* not stack-clean */
    2440             :           long v;
    2441       36386 :           if (varn(x)!=varn(y)) break;
    2442       35420 :           v = RgX_valrem(y,&y);
    2443       35420 :           if (v) x = RgX_shift_shallow(x,-v);
    2444       35420 :           if (!degpol(y)) { y = gel(y,2); break; }
    2445       31564 :           return RgX_div(x,y);
    2446             :         }
    2447          42 :         case t_RFRAC:
    2448          42 :           if (varn(gel(y,2)) != varn(x)) break;
    2449          42 :           return gdiv(x, y);
    2450             :       }
    2451     2835416 :       return RgX_Rg_divexact(x, y);
    2452        4844 :     case t_VEC: case t_COL: case t_MAT:
    2453        4844 :       lx = lg(x); z = new_chunk(lx);
    2454       53060 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    2455        4844 :       z[0] = x[0]; return z;
    2456             :   }
    2457          36 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    2458          36 :   return gdiv(x,y);
    2459             : }
    2460             : 
    2461             : static GEN
    2462      399603 : init_resultant(GEN x, GEN y)
    2463             : {
    2464      399603 :   long tx = typ(x), ty = typ(y), vx, vy;
    2465      399603 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    2466             :   {
    2467          14 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    2468          14 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    2469           0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    2470           0 :     return gen_1;
    2471             :   }
    2472      399588 :   if (tx!=t_POL) pari_err_TYPE("resultant",x);
    2473      399588 :   if (ty!=t_POL) pari_err_TYPE("resultant",y);
    2474      399588 :   if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
    2475      399511 :   vx = varn(x);
    2476      399511 :   vy = varn(y); if (vx == vy) return NULL;
    2477           7 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    2478             : }
    2479             : 
    2480             : /* x an RgX, y a scalar */
    2481             : static GEN
    2482           0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    2483             : {
    2484           0 :   *V = gpowgs(y,degpol(x)-1);
    2485           0 :   *U = gen_0; return gmul(y, *V);
    2486             : }
    2487             : 
    2488             : /* return 0 if the subresultant chain can be interrupted.
    2489             :  * Set u = NULL if the resultant is 0. */
    2490             : static int
    2491       11439 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    2492             : {
    2493       11439 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    2494             :   long du, dv, dr, degq;
    2495             : 
    2496       11439 :   if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
    2497       11439 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    2498       11292 :   du = degpol(*u);
    2499       11292 :   dv = degpol(*v);
    2500       11292 :   degq = du - dv;
    2501       11292 :   if (*um1 == gen_1)
    2502        6266 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    2503        5026 :   else if (*um1 == gen_0)
    2504        2275 :     u0 = gen_0;
    2505             :   else /* except in those 2 cases, um1 is an RgX */
    2506        2751 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2507             : 
    2508       11292 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2509        6266 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2510             :   else
    2511        5026 :     u0 = gsub(u0, gmul(q,*uze));
    2512             : 
    2513       11292 :   *um1 = *uze;
    2514       11292 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2515             : 
    2516       11292 :   *u = *v; c = *g; *g  = leading_coeff(*u);
    2517       11292 :   switch(degq)
    2518             :   {
    2519        1666 :     case 0: break;
    2520        7833 :     case 1:
    2521        7833 :       c = gmul(*h,c); *h = *g; break;
    2522        1793 :     default:
    2523        1793 :       c = gmul(gpowgs(*h,degq), c);
    2524        1793 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2525             :   }
    2526       11292 :   if (typ(c) == t_POLMOD)
    2527             :   {
    2528         918 :     c = ginv(c);
    2529         918 :     *v = RgX_Rg_mul(r,c);
    2530         918 :     *uze = RgX_Rg_mul(*uze,c);
    2531             :   }
    2532             :   else
    2533             :   {
    2534       10374 :     *v  = RgX_Rg_divexact(r,c);
    2535       10374 :     *uze= RgX_Rg_divexact(*uze,c);
    2536             :   }
    2537       11292 :   if (both_odd(du, dv)) *signh = -*signh;
    2538       11292 :   return (dr > 3);
    2539             : }
    2540             : 
    2541             : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2542             : static GEN
    2543        2395 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2544             : {
    2545             :   pari_sp av, av2;
    2546        2395 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2547             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2548             : 
    2549        2395 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2550        2395 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2551        2395 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2552        2395 :   if (tx != t_POL) {
    2553           0 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2554           0 :     return scalar_res(y,x,V,U);
    2555             :   }
    2556        2395 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2557        2395 :   if (varn(x) != varn(y))
    2558           0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2559           0 :                                         : scalar_res(y,x,V,U);
    2560        2395 :   if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
    2561        2395 :   if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
    2562        2395 :   dx = degpol(x);
    2563        2395 :   dy = degpol(y);
    2564        2395 :   signh = 1;
    2565        2395 :   if (dx < dy)
    2566             :   {
    2567         904 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2568         904 :     if (both_odd(dx, dy)) signh = -signh;
    2569             :   }
    2570        2395 :   if (dy == 0)
    2571             :   {
    2572           0 :     *V = gpowgs(gel(y,2),dx-1);
    2573           0 :     *U = gen_0; return gmul(*V,gel(y,2));
    2574             :   }
    2575        2395 :   av = avma;
    2576        2395 :   u = x = primitive_part(x, &cu);
    2577        2395 :   v = y = primitive_part(y, &cv);
    2578        2395 :   g = h = gen_1; av2 = avma;
    2579        2395 :   um1 = gen_1; uze = gen_0;
    2580             :   for(;;)
    2581             :   {
    2582        7064 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2583        4669 :     if (gc_needed(av2,1))
    2584             :     {
    2585           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2586           0 :       gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    2587             :     }
    2588             :   }
    2589             :   /* uze an RgX */
    2590        2395 :   if (!u) { *U = *V = gen_0; return gc_const(av, gen_0); }
    2591        2388 :   z = gel(v,2); du = degpol(u);
    2592        2388 :   if (du > 1)
    2593             :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2594         245 :     p1 = gpowgs(gdiv(z,h),du-1);
    2595         245 :     z = gmul(z,p1);
    2596         245 :     uze = RgX_Rg_mul(uze, p1);
    2597             :   }
    2598        2388 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2599             : 
    2600        2388 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2601        2388 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2602             :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2603        2388 :   p1 = gen_1;
    2604        2388 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2605        2388 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2606        2388 :   cu = cu? gdiv(p1,cu): p1;
    2607        2388 :   cv = cv? gdiv(p1,cv): p1;
    2608        2388 :   z = gmul(z,p1);
    2609        2388 :   *U = RgX_Rg_mul(uze,cu);
    2610        2388 :   *V = RgX_Rg_mul(vze,cv);
    2611        2388 :   return z;
    2612             : }
    2613             : GEN
    2614           0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2615             : {
    2616           0 :   pari_sp av = avma;
    2617           0 :   GEN z = subresext_i(x, y, U, V);
    2618           0 :   return gc_all(av, 3, &z, U, V);
    2619             : }
    2620             : 
    2621             : static GEN
    2622          28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2623             : {
    2624          28 :   GEN x=content(y);
    2625          28 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2626             : }
    2627             : 
    2628             : static int
    2629        6167 : must_negate(GEN x)
    2630             : {
    2631        6167 :   GEN t = leading_coeff(x);
    2632        6167 :   switch(typ(t))
    2633             :   {
    2634        4263 :     case t_INT: case t_REAL:
    2635        4263 :       return (signe(t) < 0);
    2636           0 :     case t_FRAC:
    2637           0 :       return (signe(gel(t,1)) < 0);
    2638             :   }
    2639        1904 :   return 0;
    2640             : }
    2641             : 
    2642             : static GEN
    2643          63 : gc_gcdext(pari_sp av, GEN r, GEN *u, GEN *v)
    2644             : {
    2645          63 :   if (!u && !v) return gerepileupto(av, r);
    2646          63 :   if (u  &&  v) return gc_all(av, 3, &r, u, v);
    2647           0 :   return gc_all(av, 2, &r, u ? u: v);
    2648             : }
    2649             : 
    2650             : static GEN
    2651          56 : RgX_extgcd_FpX(GEN x, GEN y, GEN p, GEN *u, GEN *v)
    2652             : {
    2653          56 :   pari_sp av = avma;
    2654          56 :   GEN r = FpX_extgcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, u, v);
    2655          56 :   if (u) *u = FpX_to_mod(*u, p);
    2656          56 :   if (v) *v = FpX_to_mod(*v, p);
    2657          56 :   return gc_gcdext(av, FpX_to_mod(r, p), u, v);
    2658             : }
    2659             : 
    2660             : static GEN
    2661           7 : RgX_extgcd_FpXQX(GEN x, GEN y, GEN pol, GEN p, GEN *U, GEN *V)
    2662             : {
    2663           7 :   pari_sp av = avma;
    2664           7 :   GEN r, T = RgX_to_FpX(pol, p);
    2665           7 :   r = FpXQX_extgcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p, U, V);
    2666           7 :   return gc_gcdext(av, FpXQX_to_mod(r, T, p), U, V);
    2667             : }
    2668             : 
    2669             : static GEN
    2670        4221 : RgX_extgcd_fast(GEN x, GEN y, GEN *U, GEN *V)
    2671             : {
    2672             :   GEN p, pol;
    2673             :   long pa;
    2674        4221 :   long t = RgX_type(x, &p,&pol,&pa);
    2675        4221 :   switch(t)
    2676             :   {
    2677          21 :     case t_FFELT:  return FFX_extgcd(x, y, pol, U, V);
    2678          56 :     case t_INTMOD: return RgX_extgcd_FpX(x, y, p, U, V);
    2679           7 :     case code(t_POLMOD, t_INTMOD):
    2680           7 :                    return RgX_extgcd_FpXQX(x, y, pol, p, U, V);
    2681        4137 :     default:       return NULL;
    2682             :   }
    2683             : }
    2684             : 
    2685             : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2686             : GEN
    2687        4256 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2688             : {
    2689             :   pari_sp av, av2, tetpil;
    2690             :   long signh; /* junk */
    2691        4256 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2692             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
    2693             : 
    2694        4256 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2695        4256 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2696        4256 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2697        4256 :   vx=varn(x);
    2698        4256 :   if (!signe(x))
    2699             :   {
    2700          14 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2701           7 :     *U = pol_0(vx); *V = pol_0(vx);
    2702           7 :     return pol_0(vx);
    2703             :   }
    2704        4242 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2705        4221 :   r = RgX_extgcd_fast(x, y, U, V);
    2706        4221 :   if (r) return r;
    2707        4137 :   dx = degpol(x); dy = degpol(y);
    2708        4137 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2709        4137 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2710             : 
    2711        3948 :   av = avma;
    2712        3948 :   u = x = primitive_part(x, &cu);
    2713        3948 :   v = y = primitive_part(y, &cv);
    2714        3948 :   g = h = gen_1; av2 = avma;
    2715        3948 :   um1 = gen_1; uze = gen_0;
    2716             :   for(;;)
    2717             :   {
    2718        4088 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2719         140 :     if (gc_needed(av2,1))
    2720             :     {
    2721           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2722           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2723             :     }
    2724             :   }
    2725        3948 :   if (uze != gen_0) {
    2726             :     GEN r;
    2727        3808 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2728        3808 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2729        3808 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2730        3808 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2731        3808 :     p1 = ginv(content(v));
    2732             :   }
    2733             :   else /* y | x */
    2734             :   {
    2735         140 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2736         140 :     uze = pol_0(vx);
    2737         140 :     p1 = gen_1;
    2738             :   }
    2739        3948 :   if (must_negate(v)) p1 = gneg(p1);
    2740        3948 :   tetpil = avma;
    2741        3948 :   z = RgX_Rg_mul(v,p1);
    2742        3948 :   *U = RgX_Rg_mul(uze,p1);
    2743        3948 :   *V = RgX_Rg_mul(vze,p1);
    2744        3948 :   gptr[0] = &z;
    2745        3948 :   gptr[1] = U;
    2746        3948 :   gptr[2] = V;
    2747        3948 :   gerepilemanysp(av,tetpil,gptr,3); return z;
    2748             : }
    2749             : 
    2750             : static GEN
    2751          14 : RgX_halfgcd_i(GEN a, GEN b)
    2752             : {
    2753          14 :   pari_sp av=avma;
    2754          14 :   long m = degpol(a), va = varn(a);
    2755             :   GEN u,u1,v,v1;
    2756          14 :   u1 = v = pol_0(va);
    2757          14 :   u = v1 = pol_1(va);
    2758          14 :   if (degpol(a)<degpol(b))
    2759             :   {
    2760           0 :     swap(a,b);
    2761           0 :     swap(u,v); swap(u1,v1);
    2762             :   }
    2763          42 :   while (2*degpol(b) >= m)
    2764             :   {
    2765          28 :     GEN r, q = RgX_pseudodivrem(a,b,&r);
    2766          28 :     GEN l = gpowgs(leading_coeff(b), degpol(a)-degpol(b)+1);
    2767          28 :     GEN g = ggcd(l, content(r));
    2768          28 :     q = RgX_Rg_div(q, g);
    2769          28 :     r = RgX_Rg_div(r, g);
    2770          28 :     l = gdiv(l, g);
    2771          28 :     a = b; b = r; swap(u,v); swap(u1,v1);
    2772          28 :     v  = RgX_sub(gmul(l,v),  RgX_mul(u, q));
    2773          28 :     v1 = RgX_sub(gmul(l,v1), RgX_mul(u1, q));
    2774          28 :     if (gc_needed(av,2))
    2775             :     {
    2776           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"halfgcd (d = %ld)",degpol(b));
    2777           0 :       gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
    2778             :     }
    2779             :   }
    2780          14 :   return gerepilecopy(av, mkvec2(mkmat22(u,u1,v,v1), mkcol2(a, b)));
    2781             : }
    2782             : 
    2783             : static GEN
    2784           7 : RgX_halfgcd_FFX(GEN x, GEN y, GEN fa)
    2785             : {
    2786           7 :   pari_sp av = avma;
    2787           7 :   GEN M = FFX_halfgcd(x, y, fa);
    2788           7 :   GEN a = FFX_add(FFX_mul(gcoeff(M,1,1), x, fa),
    2789           7 :                   FFX_mul(gcoeff(M,1,2), y, fa), fa);
    2790           7 :   GEN b = FFX_add(FFX_mul(gcoeff(M,2,1), x, fa),
    2791           7 :                   FFX_mul(gcoeff(M,2,2), y, fa), fa);
    2792           7 :   return gerepilecopy(av, mkvec2(M, mkcol2(a, b)));
    2793             : }
    2794             : 
    2795             : static GEN
    2796          14 : RgX_halfgcd_FpX(GEN x, GEN y, GEN p)
    2797             : {
    2798          14 :   pari_sp av = avma;
    2799             :   GEN M, V, a, b;
    2800          14 :   if (lgefint(p) == 3)
    2801             :   {
    2802           7 :     ulong pp = uel(p, 2);
    2803           7 :     GEN xp = RgX_to_Flx(x, pp), yp = RgX_to_Flx(y, pp);
    2804           7 :     M = Flx_halfgcd(xp, yp, pp);
    2805           7 :     a = Flx_add(Flx_mul(gcoeff(M,1,1), xp, pp),
    2806           7 :                 Flx_mul(gcoeff(M,1,2), yp, pp), pp);
    2807           7 :     b = Flx_add(Flx_mul(gcoeff(M,2,1), xp, pp),
    2808           7 :                 Flx_mul(gcoeff(M,2,2), yp, pp), pp);
    2809           7 :     M = FlxM_to_ZXM(M); a = Flx_to_ZX(a); b = Flx_to_ZX(b);
    2810             :   }
    2811             :   else
    2812             :   {
    2813           7 :     x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
    2814           7 :     M = FpX_halfgcd(x, y, p);
    2815           7 :     a = FpX_add(FpX_mul(gcoeff(M,1,1), x, p),
    2816           7 :                 FpX_mul(gcoeff(M,1,2), y, p), p);
    2817           7 :     b = FpX_add(FpX_mul(gcoeff(M,2,1), x, p),
    2818           7 :                 FpX_mul(gcoeff(M,2,2), y, p), p);
    2819             :   }
    2820          14 :   V = mkcol2(a, b);
    2821          14 :   return gerepilecopy(av, mkvec2(FpXM_to_mod(M, p), FpXC_to_mod(V, p)));
    2822             : }
    2823             : 
    2824             : static GEN
    2825           0 : RgX_halfgcd_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    2826             : {
    2827           0 :   pari_sp av = avma;
    2828           0 :   GEN M, V, a, b, T = RgX_to_FpX(pol, p);
    2829           0 :   if (signe(T)==0) pari_err_OP("halfgcd", x, y);
    2830           0 :   x = RgX_to_FpXQX(x, T, p); y = RgX_to_FpXQX(y, T, p);
    2831           0 :   M = FpXQX_halfgcd(x, y, T, p);
    2832           0 :   a = FpXX_add(FpXQX_mul(gcoeff(M,1,1), x, T, p), FpXQX_mul(gcoeff(M,1,2), y, T, p), p);
    2833           0 :   b = FpXX_add(FpXQX_mul(gcoeff(M,2,1), x, T, p), FpXQX_mul(gcoeff(M,2,2), y, T, p), p);
    2834           0 :   V = mkcol2(a, b);
    2835           0 :   return gerepilecopy(av, mkvec2(FqXM_to_mod(M, T, p), FqXC_to_mod(V, T, p)));
    2836             : }
    2837             : 
    2838             : static GEN
    2839          35 : RgX_halfgcd_fast(GEN x, GEN y)
    2840             : {
    2841             :   GEN p, pol;
    2842             :   long pa;
    2843          35 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    2844          35 :   switch(t)
    2845             :   {
    2846           7 :     case t_FFELT:  return RgX_halfgcd_FFX(x, y, pol);
    2847          14 :     case t_INTMOD: return RgX_halfgcd_FpX(x, y, p);
    2848           0 :     case code(t_POLMOD, t_INTMOD):
    2849           0 :                    return RgX_halfgcd_FpXQX(x, y, pol, p);
    2850          14 :     default:       return NULL;
    2851             :   }
    2852             : }
    2853             : 
    2854             : GEN
    2855          35 : RgX_halfgcd(GEN a, GEN b)
    2856             : {
    2857          35 :   GEN z = RgX_halfgcd_fast(a,b);
    2858          35 :   if (z) return z;
    2859          14 :   return RgX_halfgcd_i(a, b);
    2860             : }
    2861             : 
    2862             : int
    2863          70 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2864             : {
    2865          70 :   pari_sp av = avma, av2, tetpil;
    2866             :   long signh; /* junk */
    2867             :   long vx;
    2868             :   GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
    2869             : 
    2870          70 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2871          70 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2872          70 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2873          70 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2874          70 :   if (!signe(T)) {
    2875           0 :     if (degpol(x) <= amax) {
    2876           0 :       *P = RgX_copy(x);
    2877           0 :       *Q = pol_1(varn(x));
    2878           0 :       return 1;
    2879             :     }
    2880           0 :     return 0;
    2881             :   }
    2882          70 :   if (amax+bmax >= degpol(T))
    2883           0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2884             :                     mkvec3(stoi(amax), stoi(bmax), T));
    2885          70 :   vx = varn(T);
    2886          70 :   u = x = primitive_part(x, &cu);
    2887          70 :   v = T = primitive_part(T, &cv);
    2888          70 :   g = h = gen_1; av2 = avma;
    2889          70 :   um1 = gen_1; uze = gen_0;
    2890             :   for(;;)
    2891             :   {
    2892         287 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2893         287 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) return gc_bool(av,0);
    2894         287 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2895         217 :     if (gc_needed(av2,1))
    2896             :     {
    2897           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2898           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2899             :     }
    2900             :   }
    2901          70 :   if (uze == gen_0)
    2902             :   {
    2903           0 :     set_avma(av); *P = pol_0(vx); *Q = pol_1(vx);
    2904           0 :     return 1;
    2905             :   }
    2906          70 :   if (cu) uze = RgX_Rg_div(uze,cu);
    2907          70 :   p1 = ginv(content(v));
    2908          70 :   if (must_negate(v)) p1 = gneg(p1);
    2909          70 :   tetpil = avma;
    2910          70 :   *P = RgX_Rg_mul(v,p1);
    2911          70 :   *Q = RgX_Rg_mul(uze,p1);
    2912          70 :   gptr[0] = P;
    2913          70 :   gptr[1] = Q;
    2914          70 :   gerepilemanysp(av,tetpil,gptr,2); return 1;
    2915             : }
    2916             : 
    2917             : GEN
    2918           0 : RgX_chinese_coprime(GEN x, GEN y, GEN Tx, GEN Ty, GEN Tz)
    2919             : {
    2920           0 :   pari_sp av = avma;
    2921           0 :   GEN ax = RgX_mul(RgXQ_inv(Tx,Ty), Tx);
    2922           0 :   GEN p1 = RgX_mul(ax, RgX_sub(y,x));
    2923           0 :   p1 = RgX_add(x,p1);
    2924           0 :   if (!Tz) Tz = RgX_mul(Tx,Ty);
    2925           0 :   p1 = RgX_rem(p1, Tz);
    2926           0 :   return gerepileupto(av,p1);
    2927             : }
    2928             : 
    2929             : /*******************************************************************/
    2930             : /*                                                                 */
    2931             : /*                RESULTANT USING DUCOS VARIANT                    */
    2932             : /*                                                                 */
    2933             : /*******************************************************************/
    2934             : /* x^n / y^(n-1), assume n > 0 */
    2935             : static GEN
    2936      129121 : Lazard(GEN x, GEN y, long n)
    2937             : {
    2938             :   long a;
    2939             :   GEN c;
    2940             : 
    2941      129121 :   if (n == 1) return x;
    2942        8147 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    2943        8147 :   c=x; n-=a;
    2944       16918 :   while (a>1)
    2945             :   {
    2946        8770 :     a>>=1; c=gdivexact(gsqr(c),y);
    2947        8771 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    2948             :   }
    2949        8148 :   return c;
    2950             : }
    2951             : 
    2952             : /* F (x/y)^(n-1), assume n >= 1 */
    2953             : static GEN
    2954       99226 : Lazard2(GEN F, GEN x, GEN y, long n)
    2955             : {
    2956       99226 :   if (n == 1) return F;
    2957        1228 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    2958             : }
    2959             : 
    2960             : static GEN
    2961       99226 : RgX_neg_i(GEN x, long lx)
    2962             : {
    2963             :   long i;
    2964       99226 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    2965      220511 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    2966       99226 :   return y;
    2967             : }
    2968             : static GEN
    2969      259913 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    2970             : {
    2971             :   long i;
    2972             :   GEN z;
    2973      259913 :   if (isrationalzero(x)) return pol_0(varn(y));
    2974      259853 :   z = cgetg(ly,t_POL); z[1] = y[1];
    2975      578408 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2976      259846 :   return z;
    2977             : }
    2978             : static long
    2979      256593 : reductum_lg(GEN x, long lx)
    2980             : {
    2981      256593 :   long i = lx-2;
    2982      305535 :   while (i > 1 && gequal0(gel(x,i))) i--;
    2983      256593 :   return i+1;
    2984             : }
    2985             : 
    2986             : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
    2987             : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    2988             :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    2989             : static GEN
    2990       99226 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    2991             : {
    2992       99226 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
    2993             :   long p, q, j, lP, lQ;
    2994             :   pari_sp av;
    2995             : 
    2996       99226 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    2997       99225 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    2998             :   /* p > q. Very often p - 1 = q */
    2999       99226 :   av = avma;
    3000             :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    3001       99226 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    3002             : 
    3003       99226 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    3004      103561 :   for (j = q+1; j < p; j++)
    3005             :   {
    3006        4336 :     if (degpol(H) == q-1)
    3007             :     { /* h0 = coeff of degree q-1 = leading coeff */
    3008        3713 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    3009        3713 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    3010             :     }
    3011             :     else
    3012         623 :       H = RgX_shift_shallow(H, 1);
    3013        4336 :     if (j+2 < lP)
    3014             :     {
    3015        3249 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    3016        3249 :       A = A? RgX_add(A, TMP): TMP;
    3017             :     }
    3018        4336 :     if (gc_needed(av,1))
    3019             :     {
    3020         147 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    3021         147 :       gerepileall(av,A?2:1,&H,&A);
    3022             :     }
    3023             :   }
    3024       99225 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    3025       99225 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    3026       99223 :   A = A? RgX_add(A, TMP): TMP;
    3027       99222 :   A = RgX_Rg_divexact(A, p0);
    3028       99225 :   if (degpol(H) == q-1)
    3029             :   {
    3030       98833 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    3031       98832 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    3032             :   }
    3033             :   else
    3034         392 :     A = RgX_Rg_mul(addshift(H,A), q0);
    3035       99226 :   return RgX_Rg_divexact(A, s);
    3036             : }
    3037             : #undef addshift
    3038             : 
    3039             : /* Ducos's subresultant */
    3040             : GEN
    3041      128541 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    3042             : {
    3043             :   pari_sp av, av2;
    3044      128541 :   long dP, dQ, delta, sig = 1;
    3045             :   GEN cP, cQ, Z, s;
    3046             : 
    3047      128541 :   dP = degpol(P);
    3048      128541 :   dQ = degpol(Q); delta = dP - dQ;
    3049      128541 :   if (delta < 0)
    3050             :   {
    3051         595 :     if (both_odd(dP, dQ)) sig = -1;
    3052         595 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    3053             :   }
    3054      128541 :   if (sol) *sol = gen_0;
    3055      128541 :   av = avma;
    3056      128541 :   if (dQ <= 0)
    3057             :   {
    3058         651 :     if (dQ < 0) return Rg_get_0(P);
    3059         651 :     s = gpowgs(gel(Q,2), dP);
    3060         651 :     if (sig == -1) s = gerepileupto(av, gneg(s));
    3061         651 :     return s;
    3062             :   }
    3063             :   /* primitive_part is also possible here, but possibly very costly,
    3064             :    * and hardly ever worth it */
    3065      127890 :   P = Q_primitive_part(P, &cP);
    3066      127893 :   Q = Q_primitive_part(Q, &cQ);
    3067      127890 :   av2 = avma;
    3068      127890 :   s = gpowgs(leading_coeff(Q),delta);
    3069      127892 :   if (both_odd(dP, dQ)) sig = -sig;
    3070      127892 :   Z = Q;
    3071      127892 :   Q = RgX_pseudorem(P, Q);
    3072      127893 :   P = Z;
    3073      227119 :   while(degpol(Q) > 0)
    3074             :   {
    3075       99225 :     delta = degpol(P) - degpol(Q); /* > 0 */
    3076       99225 :     Z = Lazard2(Q, leading_coeff(Q), s, delta);
    3077       99226 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    3078       99226 :     Q = nextSousResultant(P, Q, Z, s);
    3079       99226 :     P = Z;
    3080       99226 :     if (gc_needed(av,1))
    3081             :     {
    3082          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    3083          13 :       gerepileall(av2,2,&P,&Q);
    3084             :     }
    3085       99226 :     s = leading_coeff(P);
    3086             :   }
    3087      127894 :   if (!signe(Q)) { set_avma(av); return Rg_get_0(Q); }
    3088      127894 :   s = Lazard(leading_coeff(Q), s, degpol(P));
    3089      127894 :   if (sig == -1) s = gneg(s);
    3090      127894 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    3091      127894 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    3092      127894 :   if (!sol) return gerepilecopy(av, s);
    3093        9499 :   *sol = P; return gc_all(av, 2, &s, sol);
    3094             : }
    3095             : 
    3096             : static GEN
    3097          14 : RgX_resultant_FpX(GEN x, GEN y, GEN p)
    3098             : {
    3099          14 :   pari_sp av = avma;
    3100             :   GEN r;
    3101          14 :   if (lgefint(p) == 3)
    3102             :   {
    3103           7 :     ulong pp = uel(p, 2);
    3104           7 :     r = utoi(Flx_resultant(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), pp));
    3105             :   }
    3106             :   else
    3107           7 :     r = FpX_resultant(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3108          14 :   return gerepileupto(av, Fp_to_mod(r, p));
    3109             : }
    3110             : 
    3111             : static GEN
    3112          21 : RgX_resultant_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3113             : {
    3114          21 :   pari_sp av = avma;
    3115          21 :   GEN r, T = RgX_to_FpX(pol, p);
    3116          21 :   r = FpXQX_resultant(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3117          21 :   return gerepileupto(av, FpX_to_mod(r, p));
    3118             : }
    3119             : 
    3120             : static GEN
    3121      399581 : resultant_fast(GEN x, GEN y)
    3122             : {
    3123             :   GEN p, pol;
    3124             :   long pa, t;
    3125      399581 :   p = init_resultant(x,y);
    3126      399581 :   if (p) return p;
    3127      399483 :   t = RgX_type2(x,y, &p,&pol,&pa);
    3128      399484 :   switch(t)
    3129             :   {
    3130        2660 :     case t_INT:    return ZX_resultant(x,y);
    3131          56 :     case t_FRAC:   return QX_resultant(x,y);
    3132          21 :     case t_FFELT:  return FFX_resultant(x,y,pol);
    3133          14 :     case t_INTMOD: return RgX_resultant_FpX(x, y, p);
    3134          21 :     case code(t_POLMOD, t_INTMOD):
    3135          21 :                    return RgX_resultant_FpXQX(x, y, pol, p);
    3136      396712 :     default:       return NULL;
    3137             :   }
    3138             : }
    3139             : 
    3140             : static GEN
    3141      278334 : RgX_resultant_sylvester(GEN x, GEN y)
    3142             : {
    3143      278334 :   pari_sp av = avma;
    3144      278334 :   return gerepileupto(av, det(RgX_sylvestermatrix(x,y)));
    3145             : }
    3146             : 
    3147             : /* Return resultant(P,Q).
    3148             :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    3149             :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    3150             :  * in the "generic" case. */
    3151             : GEN
    3152      399580 : resultant(GEN P, GEN Q)
    3153             : {
    3154      399580 :   GEN z = resultant_fast(P,Q);
    3155      399582 :   if (z) return z;
    3156      396712 :   if (isinexact(P) || isinexact(Q)) return RgX_resultant_sylvester(P,Q);
    3157      118399 :   return RgX_resultant_all(P, Q, NULL);
    3158             : }
    3159             : 
    3160             : /*******************************************************************/
    3161             : /*                                                                 */
    3162             : /*               RESULTANT USING SYLVESTER MATRIX                  */
    3163             : /*                                                                 */
    3164             : /*******************************************************************/
    3165             : static GEN
    3166      531622 : syl_RgC(GEN x, long j, long d, long D, long cp)
    3167             : {
    3168      531622 :   GEN c = cgetg(d+1,t_COL);
    3169             :   long i;
    3170     1380827 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    3171     2978318 :   for (   ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
    3172     1380827 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    3173      531622 :   return c;
    3174             : }
    3175             : static GEN
    3176      278341 : syl_RgM(GEN x, GEN y, long cp)
    3177             : {
    3178      278341 :   long j, d, dx = degpol(x), dy = degpol(y);
    3179             :   GEN M;
    3180      278341 :   if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
    3181      278341 :   if (dy < 0) return zeromat(dx,dx);
    3182      278341 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    3183      684537 :   for (j=1; j<=dy; j++) gel(M,j)    = syl_RgC(x,j,d,j+dx, cp);
    3184      403767 :   for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
    3185      278341 :   return M;
    3186             : }
    3187             : GEN
    3188      278334 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
    3189             : GEN
    3190           7 : sylvestermatrix(GEN x, GEN y)
    3191             : {
    3192           7 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    3193           7 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    3194           7 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    3195           7 :   return syl_RgM(x,y,1);
    3196             : }
    3197             : 
    3198             : GEN
    3199          21 : resultant2(GEN x, GEN y)
    3200             : {
    3201          21 :   GEN r = init_resultant(x,y);
    3202          21 :   return r? r: RgX_resultant_sylvester(x,y);
    3203             : }
    3204             : 
    3205             : /* let vx = main variable of x, v0 a variable of highest priority;
    3206             :  * return a t_POL in variable v0:
    3207             :  * if vx <= v, return subst(x, v, pol_x(v0))
    3208             :  * if vx >  v, return scalarpol(x, v0) */
    3209             : static GEN
    3210         329 : fix_pol(GEN x, long v, long v0)
    3211             : {
    3212         329 :   long vx, tx = typ(x);
    3213         329 :   if (tx != t_POL)
    3214          35 :     vx = gvar(x);
    3215             :   else
    3216             :   { /* shortcut: almost nothing to do */
    3217         294 :     vx = varn(x);
    3218         294 :     if (v == vx)
    3219             :     {
    3220         119 :       if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
    3221         119 :       return x;
    3222             :     }
    3223             :   }
    3224         210 :   if (varncmp(v, vx) > 0)
    3225             :   {
    3226         203 :     x = gsubst(x, v, pol_x(v0));
    3227         203 :     if (typ(x) != t_POL) vx = gvar(x);
    3228             :     else
    3229             :     {
    3230         196 :       vx = varn(x);
    3231         196 :       if (vx == v0) return x;
    3232             :     }
    3233             :   }
    3234          49 :   if (varncmp(vx, v0) <= 0) pari_err_TYPE("polresultant", x);
    3235          42 :   return scalarpol_shallow(x, v0);
    3236             : }
    3237             : 
    3238             : /* resultant of x and y with respect to variable v, or with respect to their
    3239             :  * main variable if v < 0. */
    3240             : GEN
    3241         476 : polresultant0(GEN x, GEN y, long v, long flag)
    3242             : {
    3243         476 :   pari_sp av = avma;
    3244             : 
    3245         476 :   if (v >= 0)
    3246             :   {
    3247         140 :     long v0 = fetch_var_higher();
    3248         140 :     x = fix_pol(x,v, v0);
    3249         140 :     y = fix_pol(y,v, v0);
    3250             :   }
    3251         476 :   switch(flag)
    3252             :   {
    3253         469 :     case 2:
    3254         469 :     case 0: x=resultant(x,y); break;
    3255           7 :     case 1: x=resultant2(x,y); break;
    3256           0 :     default: pari_err_FLAG("polresultant");
    3257             :   }
    3258         476 :   if (v >= 0) (void)delete_var();
    3259         476 :   return gerepileupto(av,x);
    3260             : }
    3261             : GEN
    3262        1491 : polresultantext0(GEN x, GEN y, long v)
    3263             : {
    3264             :   GEN R, U, V;
    3265        1491 :   pari_sp av = avma;
    3266             : 
    3267        1491 :   if (v >= 0)
    3268             :   {
    3269          14 :     long v0 = fetch_var_higher();
    3270          14 :     x = fix_pol(x,v, v0);
    3271          14 :     y = fix_pol(y,v, v0);
    3272             :   }
    3273        1491 :   R = subresext_i(x,y, &U,&V);
    3274        1491 :   if (v >= 0)
    3275             :   {
    3276          14 :     (void)delete_var();
    3277          14 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    3278          14 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    3279             :   }
    3280        1491 :   return gerepilecopy(av, mkvec3(U,V,R));
    3281             : }
    3282             : GEN
    3283        1463 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    3284             : 
    3285             : /*******************************************************************/
    3286             : /*                                                                 */
    3287             : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    3288             : /*                                                                 */
    3289             : /*******************************************************************/
    3290             : 
    3291             : /* (v - x)^d */
    3292             : static GEN
    3293         140 : caract_const(pari_sp av, GEN x, long v, long d)
    3294         140 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    3295             : 
    3296             : /* return caract(Mod(x,T)) in variable v */
    3297             : GEN
    3298      116851 : RgXQ_charpoly(GEN x, GEN T, long v)
    3299             : {
    3300      116851 :   pari_sp av = avma;
    3301      116851 :   long d = degpol(T), dx, vx, vp, v0;
    3302             :   GEN ch, L;
    3303             : 
    3304      116852 :   if (typ(x) != t_POL) return caract_const(av, x, v, d);
    3305      116838 :   vx = varn(x);
    3306      116838 :   vp = varn(T);
    3307      116838 :   if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
    3308      116838 :   if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
    3309      116838 :   dx = degpol(x);
    3310      116838 :   if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
    3311      116838 :   if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
    3312             : 
    3313      116768 :   v0 = fetch_var_higher();
    3314      116768 :   x = RgX_neg(x);
    3315      116767 :   gel(x,2) = gadd(gel(x,2), pol_x(v));
    3316      116766 :   setvarn(x, v0);
    3317      116766 :   T = leafcopy(T); setvarn(T, v0);
    3318      116766 :   ch = resultant(T, x);
    3319      116771 :   (void)delete_var();
    3320             :   /* test for silly input: x mod (deg 0 polynomial) */
    3321      116771 :   if (typ(ch) != t_POL)
    3322           7 :     pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
    3323      116764 :   L = leading_coeff(ch);
    3324      116764 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    3325      116764 :   return gerepileupto(av, ch);
    3326             : }
    3327             : 
    3328             : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    3329             :  * algebra nf[t]/(Q(t)) */
    3330             : GEN
    3331         224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    3332             : {
    3333         224 :   const char *f = "rnfcharpoly";
    3334         224 :   long dQ = degpol(Q);
    3335         224 :   pari_sp av = avma;
    3336             :   GEN T;
    3337             : 
    3338         224 :   if (v < 0) v = 0;
    3339         224 :   nf = checknf(nf); T = nf_get_pol(nf);
    3340         224 :   Q = RgX_nffix(f, T,Q,0);
    3341         224 :   switch(typ(x))
    3342             :   {
    3343          28 :     case t_INT:
    3344          28 :     case t_FRAC: return caract_const(av, x, v, dQ);
    3345          91 :     case t_POLMOD:
    3346          91 :       x = polmod_nffix2(f,T,Q, x,0);
    3347          56 :       break;
    3348          56 :     case t_POL:
    3349          56 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    3350          42 :       break;
    3351          49 :     default: pari_err_TYPE(f,x);
    3352             :   }
    3353          98 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    3354             :   /* x a t_POL in variable vQ */
    3355          56 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    3356          56 :   if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
    3357          56 :   return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    3358             : }
    3359             : 
    3360             : /*******************************************************************/
    3361             : /*                                                                 */
    3362             : /*                  GCD USING SUBRESULTANT                         */
    3363             : /*                                                                 */
    3364             : /*******************************************************************/
    3365             : static int inexact(GEN x, int *simple);
    3366             : static int
    3367       38367 : isinexactall(GEN x, int *simple)
    3368             : {
    3369       38367 :   long i, lx = lg(x);
    3370      139020 :   for (i=2; i<lx; i++)
    3371      100667 :     if (inexact(gel(x,i), simple)) return 1;
    3372       38353 :   return 0;
    3373             : }
    3374             : /* return 1 if coeff explosion is not possible */
    3375             : static int
    3376      100919 : inexact(GEN x, int *simple)
    3377             : {
    3378      100919 :   int junk = 0;
    3379      100919 :   switch(typ(x))
    3380             :   {
    3381       64729 :     case t_INT: case t_FRAC: return 0;
    3382             : 
    3383           7 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    3384             : 
    3385        4011 :     case t_INTMOD:
    3386             :     case t_FFELT:
    3387        4011 :       if (!*simple) *simple = 1;
    3388        4011 :       return 0;
    3389             : 
    3390          77 :     case t_COMPLEX:
    3391          77 :       return inexact(gel(x,1), simple)
    3392          77 :           || inexact(gel(x,2), simple);
    3393           0 :     case t_QUAD:
    3394           0 :       *simple = 0;
    3395           0 :       return inexact(gel(x,2), &junk)
    3396           0 :           || inexact(gel(x,3), &junk);
    3397             : 
    3398         819 :     case t_POLMOD:
    3399         819 :       return isinexactall(gel(x,1), simple);
    3400       31227 :     case t_POL:
    3401       31227 :       *simple = -1;
    3402       31227 :       return isinexactall(x, &junk);
    3403          49 :     case t_RFRAC:
    3404          49 :       *simple = -1;
    3405          49 :       return inexact(gel(x,1), &junk)
    3406          49 :           || inexact(gel(x,2), &junk);
    3407             :   }
    3408           0 :   *simple = -1; return 0;
    3409             : }
    3410             : 
    3411             : /* x monomial, y t_POL in the same variable */
    3412             : static GEN
    3413     1706468 : gcdmonome(GEN x, GEN y)
    3414             : {
    3415     1706468 :   pari_sp av = avma;
    3416     1706468 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    3417     1706468 :   long i, l = lg(y);
    3418     1706468 :   GEN t, v = cgetg(l, t_VEC);
    3419     1706468 :   gel(v,1) = gel(x,dx+2);
    3420     3426099 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    3421     1706468 :   t = content(v); /* gcd(lc(x), cont(y)) */
    3422     1706468 :   t = simplify_shallow(t);
    3423     1706468 :   if (dx < e) e = dx;
    3424     1706468 :   return gerepileupto(av, monomialcopy(t, e, varn(x)));
    3425             : }
    3426             : 
    3427             : static GEN
    3428      198912 : RgX_gcd_FpX(GEN x, GEN y, GEN p)
    3429             : {
    3430      198912 :   pari_sp av = avma;
    3431             :   GEN r;
    3432      198912 :   if (lgefint(p) == 3)
    3433             :   {
    3434      198898 :     ulong pp = uel(p, 2);
    3435      198898 :     r = Flx_to_ZX_inplace(Flx_gcd(RgX_to_Flx(x, pp),
    3436             :                                   RgX_to_Flx(y, pp), pp));
    3437             :   }
    3438             :   else
    3439          14 :     r = FpX_gcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3440      198912 :   return gerepileupto(av, FpX_to_mod(r, p));
    3441             : }
    3442             : 
    3443             : static GEN
    3444           7 : RgX_gcd_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3445             : {
    3446           7 :   pari_sp av = avma;
    3447           7 :   GEN r, T = RgX_to_FpX(pol, p);
    3448           7 :   if (signe(T)==0) pari_err_OP("gcd", x, y);
    3449           7 :   r = FpXQX_gcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3450           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3451             : }
    3452             : 
    3453             : static GEN
    3454        7392 : RgX_liftred(GEN x, GEN T)
    3455        7392 : { return RgXQX_red(liftpol_shallow(x), T); }
    3456             : 
    3457             : static GEN
    3458        2268 : RgX_gcd_ZXQX(GEN x, GEN y, GEN T)
    3459             : {
    3460        2268 :   pari_sp av = avma;
    3461        2268 :   GEN r = ZXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3462        2268 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    3463             : }
    3464             : 
    3465             : static GEN
    3466        1428 : RgX_gcd_QXQX(GEN x, GEN y, GEN T)
    3467             : {
    3468        1428 :   pari_sp av = avma;
    3469        1428 :   GEN r = QXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3470        1428 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    3471             : }
    3472             : 
    3473             : static GEN
    3474    11324206 : RgX_gcd_fast(GEN x, GEN y)
    3475             : {
    3476             :   GEN p, pol;
    3477             :   long pa;
    3478    11324206 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3479    11324206 :   switch(t)
    3480             :   {
    3481     9406842 :     case t_INT:    return ZX_gcd(x, y);
    3482        2499 :     case t_FRAC:   return QX_gcd(x, y);
    3483        2520 :     case t_FFELT:  return FFX_gcd(x, y, pol);
    3484      198912 :     case t_INTMOD: return RgX_gcd_FpX(x, y, p);
    3485           7 :     case code(t_POLMOD, t_INTMOD):
    3486           7 :                    return RgX_gcd_FpXQX(x, y, pol, p);
    3487        2275 :     case code(t_POLMOD, t_INT):
    3488        2275 :                    return ZX_is_monic(pol)? RgX_gcd_ZXQX(x,y,pol): NULL;
    3489        1442 :     case code(t_POLMOD, t_FRAC):
    3490        2884 :                    return RgX_is_ZX(pol) && ZX_is_monic(pol) ?
    3491        2884 :                                             RgX_gcd_QXQX(x,y,pol): NULL;
    3492     1709709 :     default:       return NULL;
    3493             :   }
    3494             : }
    3495             : 
    3496             : /* x, y are t_POL in the same variable */
    3497             : GEN
    3498    11324206 : RgX_gcd(GEN x, GEN y)
    3499             : {
    3500             :   long dx, dy;
    3501             :   pari_sp av, av1;
    3502             :   GEN d, g, h, p1, p2, u, v;
    3503    11324206 :   int simple = 0;
    3504    11324206 :   GEN z = RgX_gcd_fast(x, y);
    3505    11324206 :   if (z) return z;
    3506     1709730 :   if (isexactzero(y)) return RgX_copy(x);
    3507     1709632 :   if (isexactzero(x)) return RgX_copy(y);
    3508     1709632 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    3509        4270 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    3510        3164 :   if (isinexactall(x,&simple) || isinexactall(y,&simple))
    3511             :   {
    3512           7 :     av = avma; u = ggcd(content(x), content(y));
    3513           7 :     return gerepileupto(av, scalarpol(u, varn(x)));
    3514             :   }
    3515             : 
    3516        3157 :   av = avma;
    3517        3157 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    3518             :   else
    3519             :   {
    3520        3157 :     dx = lg(x); dy = lg(y);
    3521        3157 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    3522        3157 :     if (dy==3)
    3523             :     {
    3524           0 :       d = ggcd(gel(y,2), content(x));
    3525           0 :       return gerepileupto(av, scalarpol(d, varn(x)));
    3526             :     }
    3527        3157 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    3528        3157 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    3529        3157 :     d = ggcd(p1,p2);
    3530        3157 :     av1 = avma;
    3531        3157 :     g = h = gen_1;
    3532             :     for(;;)
    3533        1197 :     {
    3534        4354 :       GEN r = RgX_pseudorem(u,v);
    3535        4354 :       long degq, du, dv, dr = lg(r);
    3536             : 
    3537        4354 :       if (!signe(r)) break;
    3538        2205 :       if (dr <= 3)
    3539             :       {
    3540        1008 :         set_avma(av1);
    3541        1008 :         return gerepileupto(av, scalarpol(d, varn(x)));
    3542             :       }
    3543        1197 :       du = lg(u); dv = lg(v); degq = du-dv;
    3544        1197 :       u = v; p1 = g; g = leading_coeff(u);
    3545        1197 :       switch(degq)
    3546             :       {
    3547         189 :         case 0: break;
    3548         917 :         case 1:
    3549         917 :           p1 = gmul(h,p1); h = g; break;
    3550          91 :         default:
    3551          91 :           p1 = gmul(gpowgs(h,degq), p1);
    3552          91 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    3553             :       }
    3554        1197 :       v = RgX_Rg_div(r,p1);
    3555        1197 :       if (gc_needed(av1,1))
    3556             :       {
    3557           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd, dr = %ld", degpol(r));
    3558           0 :         gerepileall(av1,4, &u,&v,&g,&h);
    3559             :       }
    3560             :     }
    3561        2149 :     x = RgX_Rg_mul(primpart(v), d);
    3562             :   }
    3563        2149 :   if (must_negate(x)) x = RgX_neg(x);
    3564        2149 :   return gerepileupto(av,x);
    3565             : }
    3566             : 
    3567             : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
    3568             : static GEN
    3569         308 : RgX_disc_i(GEN P)
    3570             : {
    3571         308 :   long n = degpol(P), dd;
    3572             :   GEN N, D, L, y;
    3573         308 :   if (!signe(P) || !n) return Rg_get_0(P);
    3574         301 :   if (n == 1) return Rg_get_1(P);
    3575         301 :   if (n == 2) {
    3576         105 :     GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
    3577         105 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    3578             :   }
    3579         196 :   y = RgX_deriv(P);
    3580         196 :   N = characteristic(P);
    3581         196 :   if (signe(N)) y = gmul(y, mkintmod(gen_1,N));
    3582         196 :   if (!signe(y)) return Rg_get_0(y);
    3583         196 :   dd = n - 2 - degpol(y);
    3584         196 :   if (isinexact(P))
    3585          14 :     D = resultant2(P,y);
    3586             :   else
    3587             :   {
    3588         182 :     D = RgX_resultant_all(P, y, NULL);
    3589         182 :     if (D == gen_0) return Rg_get_0(y);
    3590             :   }
    3591         196 :   L = leading_coeff(P);
    3592         196 :   if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
    3593         196 :   if (n & 2) D = gneg(D);
    3594         196 :   return D;
    3595             : }
    3596             : 
    3597             : static GEN
    3598          42 : RgX_disc_FpX(GEN x, GEN p)
    3599             : {
    3600          42 :   pari_sp av = avma;
    3601          42 :   GEN r = FpX_disc(RgX_to_FpX(x, p), p);
    3602          42 :   return gerepileupto(av, Fp_to_mod(r, p));
    3603             : }
    3604             : 
    3605             : static GEN
    3606          28 : RgX_disc_FpXQX(GEN x, GEN pol, GEN p)
    3607             : {
    3608          28 :   pari_sp av = avma;
    3609          28 :   GEN r, T = RgX_to_FpX(pol, p);
    3610          28 :   r = FpXQX_disc(RgX_to_FpXQX(x, T, p), T, p);
    3611          28 :   return gerepileupto(av, FpX_to_mod(r, p));
    3612             : }
    3613             : 
    3614             : static GEN
    3615      122902 : RgX_disc_fast(GEN x)
    3616             : {
    3617             :   GEN p, pol;
    3618             :   long pa;
    3619      122902 :   long t = RgX_type(x, &p,&pol,&pa);
    3620      122902 :   switch(t)
    3621             :   {
    3622      122482 :     case t_INT:    return ZX_disc(x);
    3623           7 :     case t_FRAC:   return QX_disc(x);
    3624          35 :     case t_FFELT:  return FFX_disc(x, pol);
    3625          42 :     case t_INTMOD: return RgX_disc_FpX(x, p);
    3626          28 :     case code(t_POLMOD, t_INTMOD):
    3627          28 :                    return RgX_disc_FpXQX(x, pol, p);
    3628         308 :     default:       return NULL;
    3629             :   }
    3630             : }
    3631             : 
    3632             : GEN
    3633      122902 : RgX_disc(GEN x)
    3634             : {
    3635             :   pari_sp av;
    3636      122902 :   GEN z = RgX_disc_fast(x);
    3637      122902 :   if (z) return z;
    3638         308 :   av = avma;
    3639         308 :   return gerepileupto(av, RgX_disc_i(x));
    3640             : }
    3641             : 
    3642             : GEN
    3643        4672 : poldisc0(GEN x, long v)
    3644             : {
    3645        4672 :   long v0, tx = typ(x);
    3646             :   pari_sp av;
    3647             :   GEN D;
    3648        4672 :   if (tx == t_POL && (v < 0 || v == varn(x))) return RgX_disc(x);
    3649          35 :   switch(tx)
    3650             :   {
    3651           0 :     case t_QUAD:
    3652           0 :       return quad_disc(x);
    3653           0 :     case t_POLMOD:
    3654           0 :       if (v >= 0 && varn(gel(x,1)) != v) break;
    3655           0 :       return RgX_disc(gel(x,1));
    3656          14 :     case t_QFB:
    3657          14 :       return icopy(qfb_disc(x));
    3658           0 :     case t_VEC: case t_COL: case t_MAT:
    3659             :     {
    3660             :       long i;
    3661           0 :       GEN z = cgetg_copy(x, &i);
    3662           0 :       for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
    3663           0 :       return z;
    3664             :     }
    3665             :   }
    3666          21 :   if (v < 0) pari_err_TYPE("poldisc",x);
    3667          21 :   av = avma; v0 = fetch_var_higher();
    3668          21 :   x = fix_pol(x,v, v0);
    3669          14 :   D = RgX_disc(x); (void)delete_var();
    3670          14 :   return gerepileupto(av, D);
    3671             : }
    3672             : 
    3673             : GEN
    3674           7 : reduceddiscsmith(GEN x)
    3675             : {
    3676           7 :   long j, n = degpol(x);
    3677           7 :   pari_sp av = avma;
    3678             :   GEN xp, M;
    3679             : 
    3680           7 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    3681           7 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    3682           7 :   RgX_check_ZX(x,"poldiscreduced");
    3683           7 :   if (!gequal1(gel(x,n+2)))
    3684           0 :     pari_err_IMPL("nonmonic polynomial in poldiscreduced");
    3685           7 :   M = cgetg(n+1,t_MAT);
    3686           7 :   xp = ZX_deriv(x);
    3687          28 :   for (j=1; j<=n; j++)
    3688             :   {
    3689          21 :     gel(M,j) = RgX_to_RgC(xp, n);
    3690          21 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    3691             :   }
    3692           7 :   return gerepileupto(av, ZM_snf(M));
    3693             : }
    3694             : 
    3695             : /***********************************************************************/
    3696             : /**                                                                   **/
    3697             : /**                       STURM ALGORITHM                             **/
    3698             : /**              (number of real roots of x in [a,b])                 **/
    3699             : /**                                                                   **/
    3700             : /***********************************************************************/
    3701             : static GEN
    3702        1169 : R_to_Q_up(GEN x)
    3703             : {
    3704             :   long e;
    3705        1169 :   switch(typ(x))
    3706             :   {
    3707        1169 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3708           0 :     case t_REAL:
    3709           0 :       x = mantissa_real(x,&e);
    3710           0 :       return gmul2n(addiu(x,1), -e);
    3711           0 :     default: pari_err_TYPE("R_to_Q_up", x);
    3712             :       return NULL; /* LCOV_EXCL_LINE */
    3713             :   }
    3714             : }
    3715             : static GEN
    3716        1169 : R_to_Q_down(GEN x)
    3717             : {
    3718             :   long e;
    3719        1169 :   switch(typ(x))
    3720             :   {
    3721        1155 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3722          14 :     case t_REAL:
    3723          14 :       x = mantissa_real(x,&e);
    3724          14 :       return gmul2n(subiu(x,1), -e);
    3725           0 :     default: pari_err_TYPE("R_to_Q_down", x);
    3726             :       return NULL; /* LCOV_EXCL_LINE */
    3727             :   }
    3728             : }
    3729             : 
    3730             : static long
    3731        1169 : sturmpart_i(GEN x, GEN ab)
    3732             : {
    3733        1169 :   long tx = typ(x);
    3734        1169 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    3735        1169 :   if (tx != t_POL)
    3736             :   {
    3737           0 :     if (is_real_t(tx)) return 0;
    3738           0 :     pari_err_TYPE("sturm",x);
    3739             :   }
    3740        1169 :   if (lg(x) == 3) return 0;
    3741        1169 :   if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
    3742        1169 :   (void)ZX_gcd_all(x, ZX_deriv(x), &x);
    3743        1169 :   if (ab)
    3744             :   {
    3745             :     GEN A, B;
    3746        1169 :     if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
    3747        1169 :     A = R_to_Q_down(gel(ab,1));
    3748        1169 :     B = R_to_Q_up(gel(ab,2));
    3749        1169 :     ab = mkvec2(A, B);
    3750             :   }
    3751        1169 :   return ZX_sturmpart(x, ab);
    3752             : }
    3753             : /* Deprecated: RgX_sturmpart() should be preferred  */
    3754             : long
    3755        1169 : sturmpart(GEN x, GEN a, GEN b)
    3756             : {
    3757        1169 :   pari_sp av = avma;
    3758        1169 :   if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
    3759        1029 :   if (!a) a = mkmoo();
    3760        1029 :   if (!b) b = mkoo();
    3761        1029 :   return gc_long(av, sturmpart_i(x, mkvec2(a,b)));
    3762             : }
    3763             : long
    3764         140 : RgX_sturmpart(GEN x, GEN ab)
    3765         140 : { pari_sp av = avma; return gc_long(av, sturmpart_i(x, ab)); }
    3766             : 
    3767             : /***********************************************************************/
    3768             : /**                                                                   **/
    3769             : /**                        GENERIC EXTENDED GCD                       **/
    3770             : /**                                                                   **/
    3771             : /***********************************************************************/
    3772             : /* assume typ(x) = typ(y) = t_POL */
    3773             : static GEN
    3774         904 : RgXQ_inv_i(GEN x, GEN y)
    3775             : {
    3776         904 :   long vx=varn(x), vy=varn(y);
    3777             :   pari_sp av;
    3778             :   GEN u, v, d;
    3779             : 
    3780         904 :   while (vx != vy)
    3781             :   {
    3782           0 :     if (varncmp(vx,vy) > 0)
    3783             :     {
    3784           0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3785           0 :       return scalarpol(d, vy);
    3786             :     }
    3787           0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3788           0 :     x = gel(x,2); vx = gvar(x);
    3789             :   }
    3790         904 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3791         904 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3792         904 :   d = gdiv(u,d);
    3793         904 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3794         904 :   return gerepileupto(av, d);
    3795             : }
    3796             : 
    3797             : /*Assume x is a polynomial and y is not */
    3798             : static GEN
    3799         112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3800             : {
    3801         112 :   long vx = varn(x);
    3802         112 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    3803         112 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    3804          84 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    3805          56 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    3806             : }
    3807             : /* Assume x==0, y!=0 */
    3808             : static GEN
    3809          63 : zero_bezout(GEN y, GEN *U, GEN *V)
    3810             : {
    3811          63 :   *U=gen_0; *V = ginv(y); return gen_1;
    3812             : }
    3813             : 
    3814             : GEN
    3815         350 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    3816             : {
    3817         350 :   long tx=typ(x), ty=typ(y), vx;
    3818         350 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    3819         315 :   if (tx != t_POL)
    3820             :   {
    3821         140 :     if (ty == t_POL)
    3822          56 :       return scalar_bezout(y,x,v,u);
    3823             :     else
    3824             :     {
    3825          84 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    3826          84 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    3827          63 :       if (xis0) return zero_bezout(y,u,v);
    3828          42 :       else return zero_bezout(x,v,u);
    3829             :     }
    3830             :   }
    3831         175 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    3832         119 :   vx = varn(x);
    3833         119 :   if (vx != varn(y))
    3834           0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    3835           0 :                                    : scalar_bezout(y,x,v,u);
    3836         119 :   return RgX_extgcd(x,y,u,v);
    3837             : }
    3838             : 
    3839             : GEN
    3840         350 : gcdext0(GEN x, GEN y)
    3841             : {
    3842         350 :   GEN z=cgetg(4,t_VEC);
    3843         350 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    3844         350 :   return z;
    3845             : }
    3846             : 
    3847             : /*******************************************************************/
    3848             : /*                                                                 */
    3849             : /*                    GENERIC (modular) INVERSE                    */
    3850             : /*                                                                 */
    3851             : /*******************************************************************/
    3852             : 
    3853             : GEN
    3854       27251 : ginvmod(GEN x, GEN y)
    3855             : {
    3856       27251 :   long tx=typ(x);
    3857             : 
    3858       27251 :   switch(typ(y))
    3859             :   {
    3860       27251 :     case t_POL:
    3861       27251 :       if (tx==t_POL) return RgXQ_inv(x,y);
    3862       13727 :       if (is_scalar_t(tx)) return ginv(x);
    3863           0 :       break;
    3864           0 :     case t_INT:
    3865           0 :       if (tx==t_INT) return Fp_inv(x,y);
    3866           0 :       if (tx==t_POL) return gen_0;
    3867             :   }
    3868           0 :   pari_err_TYPE2("ginvmod",x,y);
    3869             :   return NULL; /* LCOV_EXCL_LINE */
    3870             : }
    3871             : 
    3872             : /***********************************************************************/
    3873             : /**                                                                   **/
    3874             : /**                         NEWTON POLYGON                            **/
    3875             : /**                                                                   **/
    3876             : /***********************************************************************/
    3877             : 
    3878             : /* assume leading coeff of x is nonzero */
    3879             : GEN
    3880          28 : newtonpoly(GEN x, GEN p)
    3881             : {
    3882          28 :   pari_sp av = avma;
    3883             :   long n, ind, a, b;
    3884             :   GEN y, vval;
    3885             : 
    3886          28 :   if (typ(x) != t_POL) pari_err_TYPE("newtonpoly",x);
    3887          28 :   n = degpol(x); if (n<=0) return cgetg(1,t_VEC);
    3888          28 :   vval = new_chunk(n+1);
    3889          28 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    3890         168 :   for (a = 0; a <= n; a++) vval[a] = gvaluation(gel(x,a),p);
    3891          42 :   for (a = 0, ind = 1; a < n; a++)
    3892             :   {
    3893          42 :     if (vval[a] != LONG_MAX) break;
    3894          14 :     gel(y,ind++) = mkoo();
    3895             :   }
    3896          84 :   for (b = a+1; b <= n; a = b, b = a+1)
    3897             :   {
    3898             :     long u1, u2, c;
    3899          70 :     while (vval[b] == LONG_MAX) b++;
    3900          56 :     u1 = vval[a] - vval[b];
    3901          56 :     u2 = b - a;
    3902         154 :     for (c = b+1; c <= n; c++)
    3903             :     {
    3904             :       long r1, r2;
    3905          98 :       if (vval[c] == LONG_MAX) continue;
    3906          70 :       r1 = vval[a] - vval[c];
    3907          70 :       r2 = c - a;
    3908          70 :       if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
    3909             :     }
    3910         154 :     while (ind <= b) gel(y,ind++) = sstoQ(u1,u2);
    3911             :   }
    3912          28 :   stackdummy((pari_sp)vval, av); return y;
    3913             : }
    3914             : 
    3915             : static GEN
    3916      274288 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
    3917             : {
    3918      274288 :   pari_sp av = avma;
    3919             :   GEN r;
    3920      274288 :   if (lgefint(p) == 3)
    3921             :   {
    3922      152381 :     ulong pp = uel(p, 2);
    3923      152381 :     r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
    3924             :                 RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
    3925             :   }
    3926             :   else
    3927      121907 :     r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
    3928      274288 :   return gerepileupto(av, FpX_to_mod(r, p));
    3929             : }
    3930             : 
    3931             : static GEN
    3932          14 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
    3933             : {
    3934          14 :   pari_sp av = avma;
    3935             :   GEN r;
    3936          14 :   if (lgefint(p) == 3)
    3937             :   {
    3938           7 :     ulong pp = uel(p, 2);
    3939           7 :     r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
    3940             :                                    RgX_to_Flx(y, pp), pp));
    3941             :   }
    3942             :   else
    3943           7 :     r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3944          14 :   return gerepileupto(av, FpX_to_mod(r, p));
    3945             : }
    3946             : 
    3947             : static GEN
    3948       12054 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
    3949             : {
    3950       12054 :   pari_sp av = avma;
    3951             :   GEN r;
    3952       12054 :   if (lgefint(p) == 3)
    3953             :   {
    3954        6088 :     ulong pp = uel(p, 2);
    3955        6088 :     r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
    3956             :                                    RgX_to_Flx(y, pp), pp));
    3957             :   }
    3958             :   else
    3959        5966 :     r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3960       12054 :   return gerepileupto(av, FpX_to_mod(r, p));
    3961             : }
    3962             : 
    3963             : static GEN
    3964         385 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
    3965             : {
    3966         385 :   pari_sp av = avma;
    3967             :   GEN r;
    3968         385 :   GEN T = RgX_to_FpX(pol, p);
    3969         385 :   if (signe(T)==0) pari_err_OP("*",x,y);
    3970         385 :   if (lgefint(p) == 3)
    3971             :   {
    3972         241 :     ulong pp = uel(p, 2);
    3973         241 :     GEN Tp = ZX_to_Flx(T, pp);
    3974         241 :     r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
    3975             :                                RgX_to_FlxqX(y, Tp, pp),
    3976             :                                RgX_to_FlxqX(S, Tp, pp), Tp, pp));
    3977             :   }
    3978             :   else
    3979         144 :     r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
    3980             :                    RgX_to_FpXQX(S, T, p), T, p);
    3981         385 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3982             : }
    3983             : 
    3984             : static GEN
    3985           0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    3986             : {
    3987           0 :   pari_sp av = avma;
    3988             :   GEN r;
    3989           0 :   GEN T = RgX_to_FpX(pol, p);
    3990           0 :   if (signe(T)==0) pari_err_OP("*",x,x);
    3991           0 :   if (lgefint(p) == 3)
    3992             :   {
    3993           0 :     ulong pp = uel(p, 2);
    3994           0 :     GEN Tp = ZX_to_Flx(T, pp);
    3995           0 :     r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
    3996             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3997             :   }
    3998             :   else
    3999           0 :     r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    4000           0 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    4001             : }
    4002             : 
    4003             : static GEN
    4004           7 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    4005             : {
    4006           7 :   pari_sp av = avma;
    4007             :   GEN r;
    4008           7 :   GEN T = RgX_to_FpX(pol, p);
    4009           7 :   if (signe(T)==0) pari_err_OP("^",x,gen_m1);
    4010           7 :   if (lgefint(p) == 3)
    4011             :   {
    4012           7 :     ulong pp = uel(p, 2);
    4013           7 :     GEN Tp = ZX_to_Flx(T, pp);
    4014           7 :     r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
    4015             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    4016             :   }
    4017             :   else
    4018           0 :     r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    4019           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    4020             : }
    4021             : 
    4022             : static GEN
    4023      839214 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
    4024             : {
    4025             :   GEN p, pol;
    4026             :   long pa;
    4027      839214 :   long t = RgX_type3(x,y,T, &p,&pol,&pa);
    4028      839212 :   switch(t)
    4029             :   {
    4030      252004 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
    4031      277755 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
    4032         105 :     case t_FFELT:  return FFXQ_mul(x, y, T, pol);
    4033      274288 :     case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
    4034         385 :     case code(t_POLMOD, t_INTMOD):
    4035         385 :                    return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
    4036       34675 :     default:       return NULL;
    4037             :   }
    4038             : }
    4039             : 
    4040             : GEN
    4041      839215 : RgXQ_mul(GEN x, GEN y, GEN T)
    4042             : {
    4043      839215 :   GEN z = RgXQ_mul_fast(x, y, T);
    4044      839216 :   if (!z) z = RgX_rem(RgX_mul(x, y), T);
    4045      839216 :   return z;
    4046             : }
    4047             : 
    4048             : static GEN
    4049      104391 : RgXQ_sqr_fast(GEN x, GEN T)
    4050             : {
    4051             :   GEN p, pol;
    4052             :   long pa;
    4053      104391 :   long t = RgX_type2(x, T, &p,&pol,&pa);
    4054      104391 :   switch(t)
    4055             :   {
    4056       78176 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
    4057       19308 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
    4058           7 :     case t_FFELT:  return FFXQ_sqr(x, T, pol);
    4059          14 :     case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
    4060           0 :     case code(t_POLMOD, t_INTMOD):
    4061           0 :                    return RgXQ_sqr_FpXQXQ(x, T, pol, p);
    4062        6886 :     default:       return NULL;
    4063             :   }
    4064             : }
    4065             : 
    4066             : GEN
    4067      104391 : RgXQ_sqr(GEN x, GEN T)
    4068             : {
    4069      104391 :   GEN z = RgXQ_sqr_fast(x, T);
    4070      104391 :   if (!z) z = RgX_rem(RgX_sqr(x), T);
    4071      104391 :   return z;
    4072             : }
    4073             : 
    4074             : static GEN
    4075       19709 : RgXQ_inv_fast(GEN x, GEN y)
    4076             : {
    4077             :   GEN p, pol;
    4078             :   long pa;
    4079       19709 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    4080       19709 :   switch(t)
    4081             :   {
    4082        5869 :     case t_INT:    return QXQ_inv(x,y);
    4083         868 :     case t_FRAC:   return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
    4084          14 :     case t_FFELT:  return FFXQ_inv(x, y, pol);
    4085       12054 :     case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
    4086           7 :     case code(t_POLMOD, t_INTMOD):
    4087           7 :                    return RgXQ_inv_FpXQXQ(x, y, pol, p);
    4088         897 :     default:       return NULL;
    4089             :   }
    4090             : }
    4091             : 
    4092             : GEN
    4093       19709 : RgXQ_inv(GEN x, GEN y)
    4094             : {
    4095       19709 :   GEN z = RgXQ_inv_fast(x, y);
    4096       19695 :   if (!z) z = RgXQ_inv_i(x, y);
    4097       19695 :   return z;
    4098             : }

Generated by: LCOV version 1.13