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.18.1 lcov report (development 30166-834405808f) Lines: 2263 2504 90.4 %
Date: 2025-04-21 09:18:53 Functions: 213 223 95.5 %
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      131022 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      33             : {
      34      131022 :   long dP=degpol(P), i, k, m;
      35             :   pari_sp av1, av2;
      36             :   GEN s,y,P_lead;
      37             : 
      38      131022 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      39      131022 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      40      131022 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      41      131022 :   y = cgetg(n+2,t_COL);
      42      131022 :   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      117785 :     m = 1;
      51      117785 :     gel(y,1) = stoi(dP);
      52             :   }
      53      131022 :   P += 2; /* strip codewords */
      54             : 
      55      131022 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      56      131022 :   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      389179 :   for (k=m; k<=n; k++)
      62             :   {
      63      258155 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      64      746712 :     for (i=1; i<k && i<=dP; i++)
      65      488575 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      66      258137 :     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      239377 :     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      239377 :       if (P_lead) s = gdiv(s, P_lead);
      78      258137 :     av2 = avma; gel(y,k+1) = gc_GEN_unsafe(av1,av2, gneg(s));
      79             :   }
      80      131024 :   return y;
      81             : }
      82             : 
      83             : GEN
      84      113606 : polsym(GEN x, long n)
      85             : {
      86      113606 :   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    87844500 : centermodii(GEN x, GEN p, GEN po2)
      92             : {
      93    87844500 :   GEN y = remii(x, p);
      94    87859141 :   switch(signe(y))
      95             :   {
      96    10615805 :     case 0: break;
      97    54436503 :     case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
      98    54261551 :       break;
      99    23201613 :     case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
     100    22981033 :       break;
     101             :   }
     102    87463609 :   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    14472812 : centermod_i(GEN x, GEN p, GEN ps2)
     116             : {
     117             :   long i, lx;
     118             :   pari_sp av;
     119             :   GEN y;
     120             : 
     121    14472812 :   if (!ps2) ps2 = shifti(p,-1);
     122    14472720 :   switch(typ(x))
     123             :   {
     124     1416901 :     case t_INT: return centermodii(x,p,ps2);
     125             : 
     126     6097462 :     case t_POL: lx = lg(x);
     127     6097462 :       y = cgetg(lx,t_POL); y[1] = x[1];
     128    42026857 :       for (i=2; i<lx; i++)
     129             :       {
     130    35928204 :         av = avma;
     131    35928204 :         gel(y,i) = gc_INT(av, centermodii(gel(x,i),p,ps2));
     132             :       }
     133     6098653 :       return normalizepol_lg(y, lx);
     134             : 
     135    29632953 :     case t_COL: pari_APPLY_same(centermodii(gel(x,i),p,ps2));
     136       13531 :     case t_MAT: pari_APPLY_same(centermod_i(gel(x,i),p,ps2));
     137             : 
     138           0 :     case t_VECSMALL: lx = lg(x);
     139             :     {
     140           0 :       ulong pp = itou(p), pps2 = itou(ps2);
     141           0 :       pari_APPLY_long(s_centermod(x[i], pp, pps2));
     142             :     }
     143             :   }
     144           0 :   return x;
     145             : }
     146             : 
     147             : GEN
     148    10959456 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
     149             : 
     150             : static GEN
     151         343 : RgX_Frobenius_deflate(GEN S, ulong p)
     152             : {
     153         343 :   if (degpol(S)%p)
     154           0 :     return NULL;
     155             :   else
     156             :   {
     157         343 :     GEN F = RgX_deflate(S, p);
     158         343 :     long i, l = lg(F);
     159        1043 :     for (i=2; i<l; i++)
     160             :     {
     161         721 :       GEN Fi = gel(F,i), R;
     162         721 :       if (typ(Fi)==t_POL)
     163             :       {
     164         259 :         if (signe(RgX_deriv(Fi))==0)
     165         238 :           gel(F,i) = RgX_Frobenius_deflate(gel(F, i), p);
     166          21 :         else return NULL;
     167             :       }
     168         462 :       else if (ispower(Fi, utoi(p), &R))
     169         462 :         gel(F,i) = R;
     170           0 :       else return NULL;
     171             :     }
     172         322 :     return F;
     173             :   }
     174             : }
     175             : 
     176             : static GEN
     177         259 : RgXY_squff(GEN f)
     178             : {
     179         259 :   long i, q, n = degpol(f);
     180         259 :   ulong p = itos_or_0(characteristic(f));
     181         259 :   GEN u = const_vec(n+1, pol_1(varn(f)));
     182         259 :   for(q = 1;;q *= p)
     183          84 :   {
     184         343 :     GEN t, v, tv, r = RgX_gcd(f, RgX_deriv(f));
     185         343 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
     186         126 :     t = RgX_div(f, r);
     187         126 :     if (degpol(t) > 0)
     188             :     {
     189             :       long j;
     190          28 :       for(j = 1;;j++)
     191             :       {
     192         140 :         v = RgX_gcd(r, t);
     193         140 :         tv = RgX_div(t, v);
     194         140 :         if (degpol(tv) > 0) gel(u, j*q) = tv;
     195         140 :         if (degpol(v) <= 0) break;
     196         112 :         r = RgX_div(r, v);
     197         112 :         t = v;
     198             :       }
     199          28 :       if (degpol(r) == 0) break;
     200             :     }
     201         105 :     if (!p) break;
     202         105 :     f = RgX_Frobenius_deflate(r, p);
     203         105 :     if (!f) { gel(u, q) = r; break; }
     204             :   }
     205         945 :   for (i = n; i; i--)
     206         931 :     if (degpol(gel(u,i))) break;
     207         259 :   setlg(u,i+1); return u;
     208             : }
     209             : 
     210             : /* Lmod contains modular factors of *F (NULL codes an empty slot: used factor)
     211             :  * Lfac accumulates irreducible factors as they are found.
     212             :  * p is a product of modular factors in Lmod[1..i-1] (NULL for p = 1), not
     213             :  * a rational factor of *F
     214             :  * Find an irreducible factor of *F divisible by p (by including
     215             :  * exhaustively further factors from Lmod[i..]); return 0 on failure, else 1.
     216             :  * Update Lmod, Lfac and *F */
     217             : static int
     218        6475 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
     219             : {
     220             :   pari_sp av;
     221             :   GEN q;
     222        6475 :   if (i == lg(Lmod)) return 0;
     223        3388 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
     224        3206 :   if (!gel(Lmod,i)) return 0;
     225        3157 :   p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
     226        3157 :   av = avma;
     227        3157 :   q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
     228        3157 :   if (degpol(q))
     229             :   {
     230        2779 :     GEN R, Q = RgX_divrem(*F, q, &R);
     231        2779 :     if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
     232             :   }
     233        2828 :   set_avma(av);
     234        2828 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
     235        2569 :   return 0;
     236             : }
     237             : 
     238             : static GEN factor_domain(GEN x, GEN flag);
     239             : 
     240             : static GEN
     241         434 : ok_bloc(GEN f, GEN BLOC, ulong c)
     242             : {
     243         434 :   GEN F = poleval(f, BLOC);
     244         434 :   return issquarefree(c ? gmul(F,mkintmodu(1,c)): F)? F: NULL;
     245             : }
     246             : static GEN
     247         112 : random_FpX_monic(long n, long v, GEN p)
     248             : {
     249         112 :   long i, d = n + 2;
     250         112 :   GEN y = cgetg(d + 1, t_POL); y[1] = evalsigne(1) | evalvarn(v);
     251         364 :   for (i = 2; i < d; i++) gel(y,i) = randomi(p);
     252         112 :   gel(y,i) = gen_1; return y;
     253             : }
     254             : static GEN
     255         273 : RgXY_factor_squarefree(GEN f, GEN dom)
     256             : {
     257         273 :   pari_sp av = avma;
     258         273 :   ulong i, c = itou_or_0(residual_characteristic(f));
     259         273 :   long vy = gvar2(f), val = RgX_valrem(f, &f), n = RgXY_degreex(f);
     260         273 :   GEN y, Lmod, F = NULL, BLOC = NULL, Lfac = coltrunc_init(degpol(f)+2);
     261         273 :   GEN gc = c? utoipos(c): NULL;
     262         273 :   if (val)
     263             :   {
     264          35 :     GEN x = pol_x(varn(f));
     265          35 :     if (dom)
     266             :     {
     267          14 :       GEN one = Rg_get_1(dom);
     268          14 :       if (typ(one) != t_INT) x = RgX_Rg_mul(x, one);
     269             :     }
     270          35 :     vectrunc_append(Lfac, x); if (!degpol(f)) return Lfac;
     271             :   }
     272         259 :   y = pol_x(vy);
     273             :   for(;;)
     274             :   {
     275         322 :     for (i = 0; !c || i < c; i++)
     276             :     {
     277         322 :       BLOC = gpowgs(gaddgs(y, i), n+1);
     278         322 :       if ((F = ok_bloc(f, BLOC, c))) break;
     279         147 :       if (c)
     280             :       {
     281         112 :         BLOC = random_FpX_monic(n, vy, gc);
     282         112 :         if ((F = ok_bloc(f, BLOC, c))) break;
     283             :       }
     284             :     }
     285         259 :     if (!c || i < c) break;
     286           0 :     n++;
     287             :   }
     288         259 :   if (DEBUGLEVEL >= 2)
     289           0 :     err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
     290         259 :   Lmod = gel(factor_domain(F,dom),1);
     291         259 :   if (DEBUGLEVEL >= 2)
     292           0 :     err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
     293         259 :   (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
     294         259 :   if (degpol(f)) vectrunc_append(Lfac, f);
     295         259 :   return gc_GEN(av, Lfac);
     296             : }
     297             : 
     298             : static GEN
     299         259 : FE_matconcat(GEN F, GEN E, long l)
     300             : {
     301         259 :   setlg(E,l); E = shallowconcat1(E);
     302         259 :   setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
     303             : }
     304             : 
     305             : static int
     306         385 : gen_cmp_RgXY(void *data, GEN x, GEN y)
     307             : {
     308         385 :   long vx = varn(x), vy = varn(y);
     309         385 :   return (vx == vy)? gen_cmp_RgX(data, x, y): -varncmp(vx, vy);
     310             : }
     311             : static GEN
     312         259 : RgXY_factor(GEN f, GEN dom)
     313             : {
     314         259 :   pari_sp av = avma;
     315             :   GEN C, F, E, cf, V;
     316             :   long i, j, l;
     317         259 :   if (dom) { GEN c = Rg_get_1(dom); if (typ(c) != t_INT) f = RgX_Rg_mul(f,c); }
     318         259 :   cf = content(f);
     319         259 :   V = RgXY_squff(gdiv(f, cf)); l = lg(V);
     320         259 :   C = factor_domain(cf, dom);
     321         259 :   F = cgetg(l+1, t_VEC); gel(F,1) = gel(C,1);
     322         259 :   E = cgetg(l+1, t_VEC); gel(E,1) = gel(C,2);
     323         770 :   for (i=1, j=2; i < l; i++)
     324             :   {
     325         511 :     GEN v = gel(V,i);
     326         511 :     if (degpol(v))
     327             :     {
     328         273 :       gel(F,j) = v = RgXY_factor_squarefree(v, dom);
     329         273 :       gel(E,j) = const_col(lg(v)-1, utoipos(i));
     330         273 :       j++;
     331             :     }
     332             :   }
     333         259 :   f = FE_matconcat(F,E,j);
     334         259 :   (void)sort_factor(f,(void*)cmp_universal, &gen_cmp_RgXY);
     335         259 :   return gc_GEN(av, f);
     336             : }
     337             : 
     338             : /***********************************************************************/
     339             : /**                                                                   **/
     340             : /**                          FACTORIZATION                            **/
     341             : /**                                                                   **/
     342             : /***********************************************************************/
     343             : static long RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var);
     344             : #define assign_or_fail(x,y) { GEN __x = x;\
     345             :   if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
     346             : }
     347             : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
     348             : 
     349             : static const long RgX_type_shift = 6;
     350             : void
     351    11719146 : RgX_type_decode(long x, long *t1, long *t2)
     352             : {
     353    11719146 :   *t1 = x >> RgX_type_shift;
     354    11719146 :   *t2 = (x & ((1L<<RgX_type_shift)-1));
     355    11719146 : }
     356             : int
     357   163889078 : RgX_type_is_composite(long t) { return t >= RgX_type_shift; }
     358             : 
     359             : static int
     360  3813640389 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     361             : {
     362             :   long j;
     363  3813640389 :   switch(typ(c))
     364             :   {
     365  2872673243 :     case t_INT:
     366  2872673243 :       break;
     367    33490629 :     case t_FRAC:
     368    33490629 :       t[1]=1; break;
     369             :       break;
     370   304715210 :     case t_REAL:
     371   304715210 :       update_prec(precision(c), pa);
     372   304714858 :       t[2]=1; break;
     373    30195032 :     case t_INTMOD:
     374    30195032 :       assign_or_fail(gel(c,1),p);
     375    30195032 :       t[3]=1; break;
     376     1843036 :     case t_FFELT:
     377     1843036 :       if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
     378     1843036 :       assign_or_fail(FF_p_i(c),p);
     379     1843036 :       t[5]=1; break;
     380   327545446 :     case t_COMPLEX:
     381   982634348 :       for (j=1; j<=2; j++)
     382             :       {
     383   655089327 :         GEN d = gel(c,j);
     384   655089327 :         switch(typ(d))
     385             :         {
     386     3171875 :           case t_INT: case t_FRAC:
     387     3171875 :             if (!*t2) *t2 = t_COMPLEX;
     388     3171875 :             t[1]=1; break;
     389   651917417 :           case t_REAL:
     390   651917417 :             update_prec(precision(d), pa);
     391   651916999 :             if (!*t2) *t2 = t_COMPLEX;
     392   651916999 :             t[2]=1; break;
     393          14 :           case t_INTMOD:
     394          14 :             assign_or_fail(gel(d,1),p);
     395          14 :             if (!signe(*p) || mod4(*p) != 3) return 0;
     396           7 :             if (!*t2) *t2 = t_COMPLEX;
     397           7 :             t[3]=1; break;
     398          21 :           case t_PADIC:
     399          21 :             update_prec(precp(d)+valp(d), pa);
     400          21 :             assign_or_fail(padic_p(d), p);
     401          21 :             if (!*t2) *t2 = t_COMPLEX;
     402          21 :             t[7]=1; break;
     403           0 :           default: return 0;
     404             :         }
     405             :       }
     406   327545021 :       if (!t[2]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
     407   327545014 :       break;
     408     2333411 :     case t_PADIC:
     409     2333411 :       update_prec(precp(c)+valp(c), pa);
     410     2333411 :       assign_or_fail(padic_p(c), p);
     411     2333411 :       t[7]=1; break;
     412        2023 :     case t_QUAD:
     413        2023 :       assign_or_fail(gel(c,1),pol);
     414        6069 :       for (j=2; j<=3; j++)
     415             :       {
     416        4046 :         GEN d = gel(c,j);
     417        4046 :         switch(typ(d))
     418             :         {
     419        4011 :           case t_INT: case t_FRAC:
     420        4011 :             t[8]=1; break;
     421          28 :           case t_INTMOD:
     422          28 :             assign_or_fail(gel(d,1),p);
     423          28 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     424          28 :             t[3]=1; break;
     425           7 :           case t_PADIC:
     426           7 :             update_prec(precp(d)+valp(d), pa);
     427           7 :             assign_or_fail(padic_p(d), p);
     428           7 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     429           7 :             t[7]=1; break;
     430           0 :           default: return 0;
     431             :         }
     432             :       }
     433        2023 :       break;
     434     4029137 :     case t_POLMOD:
     435     4029137 :       assign_or_fail(gel(c,1),pol);
     436     4028935 :       if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1))) return 0;
     437    12080176 :       for (j=1; j<=2; j++)
     438             :       {
     439             :         GEN pbis, polbis;
     440             :         long pabis;
     441     8055175 :         *t2 = t_POLMOD;
     442     8055175 :         switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
     443             :         {
     444     4511341 :           case t_INT: break;
     445      958223 :           case t_FRAC:   t[1]=1; break;
     446     2581677 :           case t_INTMOD: t[3]=1; break;
     447           7 :           case t_PADIC:  t[7]=1; update_prec(pabis,pa); break;
     448        3927 :           default: return 0;
     449             :         }
     450     8051248 :         if (pbis) assign_or_fail(pbis,p);
     451     8051248 :         if (polbis) assign_or_fail(polbis,pol);
     452             :       }
     453     4025001 :       break;
     454     6775146 :     case t_RFRAC: t[11] = 1;
     455     6775146 :       if (!settype(gel(c,1),t,p,pol,pa,ff,t2,var)) return 0;
     456     6775146 :       c = gel(c,2); /* fall through */
     457   236811206 :     case t_POL: t[10] = 1;
     458   236811206 :       if (!RgX_settype(c,t,p,pol,pa,ff,t2,var)) return 0;
     459   236835836 :       if (*var == NO_VARIABLE) { *var = varn(c); break; }
     460             :       /* if more than one free var, ensure varn() == *var fails. FIXME: should
     461             :        * keep the list of all variables, later t_POLMOD may cancel them */
     462   198683115 :       if (*var != varn(c)) *var = MAXVARN+1;
     463   198683115 :       break;
     464        2016 :     default: return 0;
     465             :   }
     466  3813658083 :   return 1;
     467             : }
     468             : /* t[0] unused. Other values, if set, indicate a coefficient of type
     469             :  * t[1] : t_FRAC
     470             :  * t[2] : t_REAL
     471             :  * t[3] : t_INTMOD
     472             :  * t[4] : Unused
     473             :  * t[5] : t_FFELT
     474             :  * t[6] : Unused
     475             :  * t[7] : t_PADIC
     476             :  * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
     477             :  * t[9]:  Unused
     478             :  * t[10]: t_POL (multivariate polynomials)
     479             :  * t[11]: t_RFRAC (recursive factorisation)  */
     480             : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
     481             :  * given by t) */
     482             : 
     483             : static long
     484    38135230 : choosesubtype(long *t, long t2)
     485             : {
     486    38135230 :   if (t2 || t[11]) return 0;
     487    35678133 :   if (t[2] && (t[3]||t[7]||t[5])) return 0;
     488    35678133 :   if (t[8]) return t_QUAD;
     489    35678105 :   if (t[7]) return t_PADIC;
     490    35678091 :   if (t[5]) return t_FFELT;
     491    35675361 :   if (t[3]) return t_INTMOD;
     492    35672309 :   if (t[2]) return t_REAL;
     493    35672225 :   if (t[1]) return t_FRAC;
     494    35429276 :   return t_INT;
     495             : }
     496             : 
     497             : static long
     498   372276332 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
     499             : {
     500   372276332 :   if (t[10] && (!*pol || var!=varn(*pol)))
     501             :   {
     502    38134478 :     long ts = choosesubtype(t, t2);
     503    38135229 :     if (ts==t_FFELT) *pol=ff;
     504    38135229 :     return ts == 0 ? t_POL: RgX_type_code(t_POL,ts);
     505             :   }
     506   334141854 :   if (t2) /* polmod/quad/complex of intmod/padic */
     507             :   {
     508    23249563 :     if (t[2] && (t[3]||t[7])) return 0;
     509    23249563 :     if (t[3]) return RgX_type_code(t2,t_INTMOD);
     510    23219708 :     if (t[7]) return RgX_type_code(t2,t_PADIC);
     511    23219659 :     if (t[2]) return t_COMPLEX;
     512      669633 :     if (t[1]) return RgX_type_code(t2,t_FRAC);
     513      238101 :     return RgX_type_code(t2,t_INT);
     514             :   }
     515   310892291 :   if (t[5]) /* ffelt */
     516             :   {
     517      225093 :     if (t[2]||t[8]||t[9]) return 0;
     518      225093 :     *pol=ff; return t_FFELT;
     519             :   }
     520   310667198 :   if (t[2]) /* inexact, real */
     521             :   {
     522    47678175 :     if (t[3]||t[7]||t[9]) return 0;
     523    47678182 :     return t_REAL;
     524             :   }
     525   262989023 :   if (t[10])
     526             :   {
     527           0 :     long ts = choosesubtype(t, t2);
     528           0 :     if (ts==t_FFELT) *pol=ff;
     529           0 :     return ts == 0 ? t_POL: RgX_type_code(t_POL,ts);
     530             :   }
     531   262989023 :   if (t[8]) return RgX_type_code(t_QUAD,t_INT);
     532   262988162 :   if (t[3]) return t_INTMOD;
     533   258838086 :   if (t[7]) return t_PADIC;
     534   258460094 :   if (t[1]) return t_FRAC;
     535   250479586 :   return t_INT;
     536             : }
     537             : 
     538             : static long
     539   664459550 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     540             : {
     541   664459550 :   long i, lx = lg(x);
     542  2334014842 :   for (i=2; i<lx; i++)
     543  1669557831 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     544   664457011 :   return 1;
     545             : }
     546             : 
     547             : static long
     548   269299742 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     549             : {
     550   269299742 :   long i, l = lg(x);
     551  2344993323 :   for (i = 1; i<l; i++)
     552  2075696113 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     553   269297210 :   return 1;
     554             : }
     555             : 
     556             : static long
     557    48294398 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     558             : {
     559    48294398 :   long i, l = lg(x);
     560   284452646 :   for (i = 1; i < l; i++)
     561   236160410 :     if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     562    48292236 :   return 1;
     563             : }
     564             : 
     565             : long
     566   171944271 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
     567             : {
     568   171944271 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     569   171944271 :   long t2 = 0, var = NO_VARIABLE;
     570   171944271 :   GEN ff = NULL;
     571   171944271 :   *p = *pol = NULL; *pa = LONG_MAX;
     572   171944271 :   switch(typ(x))
     573             :   {
     574    55556045 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
     575             :     case t_COMPLEX: case t_PADIC: case t_QUAD:
     576    55556045 :       if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     577    55556046 :       break;
     578   115846792 :     case t_POL: case t_SER:
     579   115846792 :       if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     580   115846513 :       break;
     581          21 :     case t_VEC: case t_COL:
     582          21 :       if(!RgC_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
     583          21 :       break;
     584         126 :     case t_MAT:
     585         126 :       if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
     586         126 :       break;
     587      541287 :     default: return 0;
     588             :   }
     589   171402706 :   return choosetype(t,t2,ff,pol,var);
     590             : }
     591             : 
     592             : long
     593     2546402 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     594             : {
     595     2546402 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     596     2546402 :   long t2 = 0, var = NO_VARIABLE;
     597     2546402 :   GEN ff = NULL;
     598     2546402 :   *p = *pol = NULL; *pa = LONG_MAX;
     599     2546402 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     600     2546349 :   return choosetype(t,t2,ff,pol,var);
     601             : }
     602             : 
     603             : long
     604     6063663 : RgX_Rg_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     605             : {
     606     6063663 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     607     6063663 :   long t2 = 0, var = NO_VARIABLE;
     608     6063663 :   GEN ff = NULL;
     609     6063663 :   *p = *pol = NULL; *pa = LONG_MAX;
     610     6063663 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     611     6063668 :   if (!settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     612     6063658 :   return choosetype(t,t2,ff,pol,var);
     613             : }
     614             : 
     615             : long
     616   149293380 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     617             : {
     618   149293380 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     619   149293380 :   long t2 = 0, var = NO_VARIABLE;
     620   149293380 :   GEN ff = NULL;
     621   149293380 :   *p = *pol = NULL; *pa = LONG_MAX;
     622   298584731 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     623   149294385 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     624   149291233 :   return choosetype(t,t2,ff,pol,var);
     625             : }
     626             : 
     627             : long
     628     1527549 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
     629             : {
     630     1527549 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     631     1527549 :   long t2 = 0, var = NO_VARIABLE;
     632     1527549 :   GEN ff = NULL;
     633     1527549 :   *p = *pol = NULL; *pa = LONG_MAX;
     634     3055099 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     635     3055100 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
     636     1527550 :       !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
     637     1527550 :   return choosetype(t,t2,ff,pol,var);
     638             : }
     639             : 
     640             : long
     641      687397 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
     642             : {
     643      687397 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     644      687397 :   long t2 = 0, var = NO_VARIABLE;
     645      687397 :   GEN ff = NULL;
     646      687397 :   *p = *pol = NULL; *pa = LONG_MAX;
     647      687397 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     648      686369 :   return choosetype(t,t2,ff,pol,var);
     649             : }
     650             : 
     651             : long
     652      774594 : RgV_type(GEN x, GEN *p, GEN *pol, long *pa)
     653             : {
     654      774594 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     655      774594 :   long t2 = 0, var = NO_VARIABLE;
     656      774594 :   GEN ff = NULL;
     657      774594 :   *p = *pol = NULL; *pa = LONG_MAX;
     658      774594 :   if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     659      774594 :   return choosetype(t,t2,ff,pol,var);
     660             : }
     661             : 
     662             : long
     663         203 : RgV_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     664             : {
     665         203 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     666         203 :   long t2 = 0, var = NO_VARIABLE;
     667         203 :   GEN ff = NULL;
     668         203 :   *p = *pol = NULL; *pa = LONG_MAX;
     669         406 :   if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     670         203 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     671         203 :   return choosetype(t,t2,ff,pol,var);
     672             : }
     673             : 
     674             : long
     675    32365305 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     676             : {
     677    32365305 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     678    32365305 :   long t2 = 0, var = NO_VARIABLE;
     679    32365305 :   GEN ff = NULL;
     680    32365305 :   *p = *pol = NULL; *pa = LONG_MAX;
     681    64730324 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     682    32366223 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     683    32364218 :   return choosetype(t,t2,ff,pol,var);
     684             : }
     685             : 
     686             : long
     687     7621013 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     688             : {
     689     7621013 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
     690     7621013 :   long t2 = 0, var = NO_VARIABLE;
     691     7621013 :   GEN ff = NULL;
     692     7621013 :   *p = *pol = NULL; *pa = LONG_MAX;
     693    15241606 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     694     7621187 :       !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     695     7620487 :   return choosetype(t,t2,ff,pol,var);
     696             : }
     697             : 
     698             : GEN
     699       59375 : factor0(GEN x, GEN flag)
     700             : {
     701             :   ulong B;
     702       59375 :   long tx = typ(x);
     703       59375 :   if (!flag) return factor(x);
     704         259 :   if ((tx != t_INT && tx!=t_FRAC) || typ(flag) != t_INT)
     705         175 :     return factor_domain(x, flag);
     706          84 :   if (signe(flag) < 0) pari_err_FLAG("factor");
     707          84 :   switch(lgefint(flag))
     708             :   {
     709          14 :     case 2: B = 0; break;
     710          70 :     case 3: B = flag[2]; break;
     711           0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     712             :              return NULL; /*LCOV_EXCL_LINE*/
     713             :   }
     714          84 :   return boundfact(x, B);
     715             : }
     716             : 
     717             : GEN
     718      156896 : deg1_from_roots(GEN L, long v)
     719             : {
     720      156896 :   long i, l = lg(L);
     721      156896 :   GEN z = cgetg(l,t_COL);
     722      462513 :   for (i=1; i<l; i++)
     723      305617 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     724      156896 :   return z;
     725             : }
     726             : GEN
     727       63896 : roots_from_deg1(GEN x)
     728             : {
     729       63896 :   long i,l = lg(x);
     730       63896 :   GEN r = cgetg(l,t_VEC);
     731      392970 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     732       63894 :   return r;
     733             : }
     734             : 
     735             : static GEN
     736          42 : Qi_factor_p(GEN p)
     737             : {
     738          42 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     739          42 :   return mkcomplex(a, b);
     740             : }
     741             : 
     742             : static GEN
     743          49 : Qi_primpart(GEN x, GEN *c)
     744             : {
     745          49 :   GEN a = real_i(x), b = imag_i(x), n = gcdii(a, b);
     746          49 :   *c = n; if (n == gen_1) return x;
     747          49 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     748             : }
     749             : 
     750             : static GEN
     751          70 : Qi_primpart_try(GEN x, GEN c)
     752             : {
     753             :   GEN r, y;
     754          70 :   if (typ(x) == t_INT)
     755             :   {
     756          42 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     757             :   }
     758             :   else
     759             :   {
     760          28 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     761          28 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     762          14 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     763             :   }
     764          56 :   return y;
     765             : }
     766             : 
     767             : static int
     768          91 : Qi_cmp(GEN x, GEN y)
     769             : {
     770             :   int v;
     771          91 :   if (typ(x) != t_COMPLEX)
     772           0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     773          91 :   if (typ(y) != t_COMPLEX) return 1;
     774          63 :   v = cmpii(gel(x,2), gel(y,2));
     775          63 :   if (v) return v;
     776          28 :   return gcmp(gel(x,1), gel(y,1));
     777             : }
     778             : 
     779             : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
     780             : static GEN
     781         469 : Qi_normal(GEN x)
     782             : {
     783         469 :   if (typ(x) != t_COMPLEX) return absi_shallow(x);
     784         469 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     785         469 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     786         469 :   return x;
     787             : }
     788             : 
     789             : static GEN
     790          49 : Qi_factor(GEN x)
     791             : {
     792          49 :   pari_sp av = avma;
     793          49 :   GEN a = real_i(x), b = imag_i(x), d = gen_1, n, y, fa, P, E, P2, E2;
     794          49 :   long t1 = typ(a);
     795          49 :   long t2 = typ(b), i, j, l, exp = 0;
     796          49 :   if (t1 == t_FRAC) d = gel(a,2);
     797          49 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     798          49 :   if (d == gen_1) y = x;
     799             :   else
     800             :   {
     801          21 :     y = gmul(x, d);
     802          21 :     a = real_i(y); t1 = typ(a);
     803          21 :     b = imag_i(y); t2 = typ(b);
     804             :   }
     805          49 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     806          49 :   y = Qi_primpart(y, &n);
     807          49 :   fa = factor(cxnorm(y));
     808          49 :   P = gel(fa,1);
     809          49 :   E = gel(fa,2); l = lg(P);
     810          49 :   P2 = cgetg(l, t_COL);
     811          49 :   E2 = cgetg(l, t_COL);
     812         105 :   for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
     813             :   { /* either p = 2 (ramified) or those factors split in Q(i) */
     814          56 :     GEN p = gel(P,i), w, w2, t, we, pe;
     815          56 :     long v, e = itos(gel(E,i));
     816          56 :     int is2 = absequaliu(p, 2);
     817          56 :     w = is2? mkcomplex(gen_1,gen_1): Qi_factor_p(p);
     818          56 :     w2 = Qi_normal( conj_i(w) );
     819             :     /* w * w2 * I^3 = p, w2 = conj(w) * I */
     820          56 :     pe = powiu(p, e);
     821          56 :     we = gpowgs(w, e);
     822          56 :     t = Qi_primpart_try( gmul(y, conj_i(we)), pe );
     823          56 :     if (t) y = t; /* y /= w^e */
     824             :     else {
     825             :       /* y /= conj(w)^e, should be y /= w2^e */
     826          14 :       y = Qi_primpart_try( gmul(y, we), pe );
     827          14 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     828             :     }
     829          56 :     gel(P,i) = w;
     830          56 :     v = Z_pvalrem(n, p, &n);
     831          56 :     if (v) {
     832           7 :       exp -= v; /* += 3*v mod 4 */
     833           7 :       if (is2) v <<= 1; /* 2 = w^2 I^3 */
     834             :       else {
     835           0 :         gel(P2,j) = w2;
     836           0 :         gel(E2,j) = utoipos(v); j++;
     837             :       }
     838           7 :       gel(E,i) = stoi(e + v);
     839             :     }
     840          56 :     v = Z_pvalrem(d, p, &d);
     841          56 :     if (v) {
     842           7 :       exp += v; /* -= 3*v mod 4 */
     843           7 :       if (is2) v <<= 1; /* 2 is ramified */
     844             :       else {
     845           7 :         gel(P2,j) = w2;
     846           7 :         gel(E2,j) = utoineg(v); j++;
     847             :       }
     848           7 :       gel(E,i) = stoi(e - v);
     849             :     }
     850          56 :     exp &= 3;
     851             :   }
     852          49 :   if (j > 1) {
     853           7 :     long k = 1;
     854           7 :     GEN P1 = cgetg(l, t_COL);
     855           7 :     GEN E1 = cgetg(l, t_COL);
     856             :     /* remove factors with exponent 0 */
     857          14 :     for (i = 1; i < l; i++)
     858           7 :       if (signe(gel(E,i)))
     859             :       {
     860           0 :         gel(P1,k) = gel(P,i);
     861           0 :         gel(E1,k) = gel(E,i);
     862           0 :         k++;
     863             :       }
     864           7 :     setlg(P1, k); setlg(E1, k);
     865           7 :     setlg(P2, j); setlg(E2, j);
     866           7 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     867             :   }
     868          49 :   if (!equali1(n) || !equali1(d))
     869             :   {
     870          28 :     GEN Fa = factor(Qdivii(n, d));
     871          28 :     P = gel(Fa,1); l = lg(P);
     872          28 :     E = gel(Fa,2);
     873          70 :     for (i = 1; i < l; i++)
     874             :     {
     875          42 :       GEN w, p = gel(P,i);
     876             :       long e;
     877             :       int is2;
     878          42 :       switch(mod4(p))
     879             :       {
     880          14 :         case 3: continue;
     881          14 :         case 2: is2 = 1; break;
     882          14 :         default:is2 = 0; break;
     883             :       }
     884          28 :       e = itos(gel(E,i));
     885          28 :       w = is2? mkcomplex(gen_1,gen_1): Qi_factor_p(p);
     886          28 :       gel(P,i) = w;
     887          28 :       if (is2)
     888          14 :         gel(E,i) = stoi(2*e);
     889             :       else
     890             :       {
     891          14 :         P = vec_append(P, Qi_normal( conj_i(w) ));
     892          14 :         E = vec_append(E, gel(E,i));
     893             :       }
     894          28 :       exp -= e; /* += 3*e mod 4 */
     895          28 :       exp &= 3;
     896             :     }
     897          28 :     gel(Fa,1) = P;
     898          28 :     gel(Fa,2) = E;
     899          28 :     fa = famat_mul_shallow(fa, Fa);
     900             :   }
     901          49 :   fa = sort_factor(fa, (void*)&Qi_cmp, &cmp_nodata);
     902             : 
     903          49 :   y = gmul(y, powIs(exp));
     904          49 :   if (!gequal1(y)) {
     905          35 :     gel(fa,1) = vec_prepend(gel(fa,1), y);
     906          35 :     gel(fa,2) = vec_prepend(gel(fa,2), gen_1);
     907             :   }
     908          49 :   return gc_GEN(av, fa);
     909             : }
     910             : 
     911             : GEN
     912        9698 : Q_factor_limit(GEN x, ulong lim)
     913             : {
     914        9698 :   pari_sp av = avma;
     915             :   GEN a, b;
     916        9698 :   if (typ(x) == t_INT) return Z_factor_limit(x, lim);
     917        4859 :   a = Z_factor_limit(gel(x,1), lim);
     918        4859 :   b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
     919        4859 :   return gc_GEN(av, ZM_merge_factor(a,b));
     920             : }
     921             : GEN
     922       21572 : Q_factor(GEN x)
     923             : {
     924       21572 :   pari_sp av = avma;
     925             :   GEN a, b;
     926       21572 :   if (typ(x) == t_INT) return Z_factor(x);
     927          34 :   a = Z_factor(gel(x,1));
     928          35 :   b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
     929          35 :   return gc_GEN(av, ZM_merge_factor(a,b));
     930             : }
     931             : 
     932             : /* replace quadratic number over Fp or Q by t_POL in v */
     933             : static GEN
     934         420 : quadratic_to_RgX(GEN z, long v)
     935             : {
     936             :   GEN a, b;
     937         420 :   switch(typ(z))
     938             :   {
     939         343 :     case t_INT: case t_FRAC: case t_INTMOD: return z;
     940          35 :     case t_COMPLEX: a = gel(z,2); b = gel(z,1); break;
     941          42 :     case t_QUAD: a = gel(z,3); b = gel(z,2); break;
     942           0 :     default: pari_err_IMPL("factor for general polynomials"); /* paranoia */
     943             :              return NULL; /* LCOV_EXCL_LINE */
     944             :   }
     945          77 :   return deg1pol_shallow(a, b, v);
     946             : }
     947             : /* replace t_QUAD/t_COMPLEX [of rationals] coeffs by t_POL in v */
     948             : static GEN
     949          98 : RgX_fix_quadratic(GEN x, long v)
     950         518 : { pari_APPLY_pol_normalized(quadratic_to_RgX(gel(x,i), v)); }
     951             : static GEN
     952         252 : RgXQ_factor_i(GEN x, GEN T, GEN p, long t1, long t2, long *pv)
     953             : {
     954         252 :   *pv = -1;
     955         252 :   if (t2 == t_PADIC) return NULL;
     956         217 :   if (t2 == t_INTMOD)
     957             :   {
     958          56 :     T = RgX_to_FpX(T,p);
     959          56 :     if (!FpX_is_irred(T,p)) return NULL;
     960             :   }
     961         196 :   if (t1 != t_POLMOD)
     962             :   { /* replace w in x by t_POL */
     963          98 :     if (t2 != t_INTMOD) T = leafcopy(T);
     964          98 :     *pv = fetch_var(); setvarn(T, *pv);
     965          98 :     x = RgX_fix_quadratic(x, *pv);
     966             :   }
     967         196 :   if (t2 == t_INTMOD) return factmod(x, mkvec2(p,T));
     968         161 :   return nffactor(T, x);
     969             : }
     970             : static GEN
     971         252 : RgXQ_factor(GEN x, GEN T, GEN p, long tx)
     972             : {
     973         252 :   pari_sp av = avma;
     974             :   long t1, t2, v;
     975             :   GEN w, y;
     976         252 :   RgX_type_decode(tx, &t1, &t2);
     977         252 :   y = RgXQ_factor_i(x, T, p, t1, t2, &v);
     978         252 :   if (!y) pari_err_IMPL("factor for general polynomials");
     979         196 :   if (v < 0) return gc_upto(av, y);
     980             :   /* substitute back w */
     981          98 :   w = (t1 == t_COMPLEX)? gen_I(): mkquad(T,gen_0,gen_1);
     982          98 :   gel(y,1) = gsubst(liftpol_shallow(gel(y,1)), v, w);
     983          98 :   (void)delete_var(); return gc_GEN(av, y);
     984             : }
     985             : 
     986             : static GEN
     987          28 : RX_factor(GEN x, long prec)
     988             : {
     989          28 :   GEN y = cgetg(3,t_MAT), R, P;
     990          28 :   pari_sp av = avma;
     991          28 :   long v = varn(x), i, l, r1;
     992             : 
     993          28 :   R = cleanroots(x, prec); l = lg(R);
     994          70 :   for (r1 = 1; r1 < l; r1++)
     995          49 :     if (typ(gel(R,r1)) == t_COMPLEX) break;
     996          28 :   l = (r1+l)>>1; P = cgetg(l,t_COL);
     997          70 :   for (i = 1; i < r1; i++)
     998          42 :     gel(P,i) = deg1pol_shallow(gen_1, negr(gel(R,i)), v);
     999          35 :   for (   ; i < l; i++)
    1000             :   {
    1001           7 :     GEN a = gel(R,2*i-r1), t;
    1002           7 :     t = gmul2n(gel(a,1), 1); togglesign(t);
    1003           7 :     gel(P,i) = deg2pol_shallow(gen_1, t, gnorm(a), v);
    1004             :   }
    1005          28 :   gel(y,1) = gc_upto(av, P);
    1006          28 :   gel(y,2) = const_col(l-1, gen_1); return y;
    1007             : }
    1008             : static GEN
    1009          21 : CX_factor(GEN x, long prec)
    1010             : {
    1011          21 :   GEN y = cgetg(3,t_MAT), R;
    1012          21 :   pari_sp av = avma;
    1013          21 :   long v = varn(x);
    1014             : 
    1015          21 :   R = roots(x, prec);
    1016          21 :   gel(y,1) = gc_upto(av, deg1_from_roots(R, v));
    1017          21 :   gel(y,2) = const_col(degpol(x), gen_1); return y;
    1018             : }
    1019             : 
    1020             : static GEN
    1021       13832 : RgX_factor(GEN x, GEN dom)
    1022             : {
    1023             :   GEN  p, T;
    1024       13832 :   long pa, tx = dom ? RgX_Rg_type(x,dom,&p,&T,&pa): RgX_type(x,&p,&T,&pa);
    1025       13832 :   if (tx>>RgX_type_shift==t_POL) tx = t_POL;
    1026       13832 :   switch(tx)
    1027             :   {
    1028           7 :     case 0: pari_err_IMPL("factor for general polynomials");
    1029         259 :     case t_POL: return RgXY_factor(x, dom);
    1030       12789 :     case t_INT: return ZX_factor(x);
    1031           7 :     case t_FRAC: return QX_factor(x);
    1032         329 :     case t_INTMOD: return factmod(x, p);
    1033          42 :     case t_PADIC: return factorpadic(x, p, pa);
    1034          98 :     case t_FFELT: return FFX_factor(x, T);
    1035          21 :     case t_COMPLEX: return CX_factor(x, pa);
    1036          28 :     case t_REAL: return RX_factor(x, pa);
    1037             :   }
    1038         252 :   return RgXQ_factor(x, T, p, tx);
    1039             : }
    1040             : 
    1041             : static GEN
    1042       63036 : factor_domain(GEN x, GEN dom)
    1043             : {
    1044       63036 :   long tx = typ(x), tdom = dom ? typ(dom): 0;
    1045             :   pari_sp av;
    1046             : 
    1047       63036 :   if (gequal0(x))
    1048          63 :     switch(tx)
    1049             :     {
    1050          63 :       case t_INT:
    1051             :       case t_COMPLEX:
    1052             :       case t_POL:
    1053          63 :       case t_RFRAC: return prime_fact(x);
    1054           0 :       default: pari_err_TYPE("factor",x);
    1055             :     }
    1056       62973 :   av = avma;
    1057       62973 :   switch(tx)
    1058             :   {
    1059        2632 :     case t_POL: return RgX_factor(x, dom);
    1060          35 :     case t_RFRAC: {
    1061          35 :       GEN a = gel(x,1), b = gel(x,2);
    1062          35 :       GEN y = famat_inv_shallow(RgX_factor(b, dom));
    1063          35 :       if (typ(a)==t_POL) y = famat_mul_shallow(RgX_factor(a, dom), y);
    1064          35 :       return gc_GEN(av, sort_factor_pol(y, cmp_universal));
    1065             :     }
    1066       60236 :     case t_INT:  if (tdom==0 || tdom==t_INT) return Z_factor(x);
    1067          28 :     case t_FRAC: if (tdom==0 || tdom==t_INT) return Q_factor(x);
    1068             :     case t_COMPLEX: /* fall through */
    1069          49 :       if (tdom==0 || tdom==t_COMPLEX)
    1070          49 :       { GEN y = Qi_factor(x); if (y) return y; }
    1071             :       /* fall through */
    1072             :   }
    1073           0 :   pari_err_TYPE("factor",x);
    1074             :   return NULL; /* LCOV_EXCL_LINE */
    1075             : }
    1076             : 
    1077             : GEN
    1078       62343 : factor(GEN x) { return factor_domain(x, NULL); }
    1079             : 
    1080             : /*******************************************************************/
    1081             : /*                                                                 */
    1082             : /*                     ROOTS --> MONIC POLYNOMIAL                  */
    1083             : /*                                                                 */
    1084             : /*******************************************************************/
    1085             : static GEN
    1086     1719017 : normalized_mul(void *E, GEN x, GEN y)
    1087             : {
    1088     1719017 :   long a = gel(x,1)[1], b = gel(y,1)[1];
    1089             :   (void) E;
    1090     1719012 :   return mkvec2(mkvecsmall(a + b),
    1091     1719017 :                 RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
    1092             : }
    1093             : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
    1094             : static GEN
    1095     1038941 : normalized_to_RgX(GEN L)
    1096             : {
    1097     1038941 :   long i, a = gel(L,1)[1];
    1098     1038941 :   GEN A = gel(L,2);
    1099     1038941 :   GEN z = cgetg(a + 3, t_POL);
    1100     1038941 :   z[1] = evalsigne(1) | evalvarn(varn(A));
    1101     5785701 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
    1102     1044092 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
    1103     1038940 :   gel(z,i) = gen_1; return z;
    1104             : }
    1105             : 
    1106             : static GEN
    1107          14 : roots_to_pol_FpV(GEN x, long v, GEN p)
    1108             : {
    1109          14 :   pari_sp av = avma;
    1110             :   GEN r;
    1111          14 :   if (lgefint(p) == 3)
    1112             :   {
    1113          14 :     ulong pp = uel(p, 2);
    1114          14 :     r = Flx_to_ZX_inplace(Flv_roots_to_pol(RgV_to_Flv(x, pp), pp, v<<VARNSHIFT));
    1115             :   }
    1116             :   else
    1117           0 :     r = FpV_roots_to_pol(RgV_to_FpV(x, p), p, v);
    1118          14 :   return gc_upto(av, FpX_to_mod(r, p));
    1119             : }
    1120             : 
    1121             : static GEN
    1122           7 : roots_to_pol_FqV(GEN x, long v, GEN pol, GEN p)
    1123             : {
    1124           7 :   pari_sp av = avma;
    1125           7 :   GEN r, T = RgX_to_FpX(pol, p);
    1126           7 :   if (signe(T)==0) pari_err_OP("/", x, pol);
    1127           7 :   r = FqV_roots_to_pol(RgC_to_FqC(x, T, p), T, p, v);
    1128           7 :   return gc_upto(av, FpXQX_to_mod(r, T, p));
    1129             : }
    1130             : 
    1131             : static GEN
    1132      774510 : roots_to_pol_fast(GEN x, long v)
    1133             : {
    1134             :   GEN p, pol;
    1135             :   long pa;
    1136      774510 :   long t = RgV_type(x, &p,&pol,&pa);
    1137      774510 :   switch(t)
    1138             :   {
    1139          14 :     case t_INTMOD: return roots_to_pol_FpV(x, v, p);
    1140          14 :     case t_FFELT:  return FFV_roots_to_pol(x, pol, v);
    1141           7 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    1142           7 :                    return roots_to_pol_FqV(x, v, pol, p);
    1143      774475 :     default:       return NULL;
    1144             :   }
    1145             : }
    1146             : 
    1147             : /* compute prod (x - a[i]) */
    1148             : GEN
    1149      774584 : roots_to_pol(GEN a, long v)
    1150             : {
    1151      774584 :   pari_sp av = avma;
    1152      774584 :   long i, k, lx = lg(a);
    1153             :   GEN L;
    1154      774584 :   if (lx == 1) return pol_1(v);
    1155      774510 :   L = roots_to_pol_fast(a, v);
    1156      774510 :   if (L) return L;
    1157      774475 :   L = cgetg(lx, t_VEC);
    1158     1661241 :   for (k=1,i=1; i<lx-1; i+=2)
    1159             :   {
    1160      886766 :     GEN s = gel(a,i), t = gel(a,i+1);
    1161      886766 :     GEN x0 = gmul(s,t);
    1162      886766 :     GEN x1 = gneg(gadd(s,t));
    1163      886766 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1164             :   }
    1165     1466431 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
    1166      691956 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
    1167      774475 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1168      774475 :   return gc_upto(av, normalized_to_RgX(L));
    1169             : }
    1170             : 
    1171             : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
    1172             : GEN
    1173      264464 : roots_to_pol_r1(GEN a, long v, long r1)
    1174             : {
    1175      264464 :   pari_sp av = avma;
    1176      264464 :   long i, k, lx = lg(a);
    1177             :   GEN L;
    1178      264464 :   if (lx == 1) return pol_1(v);
    1179      264464 :   L = cgetg(lx, t_VEC);
    1180      710937 :   for (k=1,i=1; i<r1; i+=2)
    1181             :   {
    1182      446474 :     GEN s = gel(a,i), t = gel(a,i+1);
    1183      446474 :     GEN x0 = gmul(s,t);
    1184      446474 :     GEN x1 = gneg(gadd(s,t));
    1185      446471 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1186             :   }
    1187      336516 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
    1188       72053 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
    1189      925155 :   for (i=r1+1; i<lx; i++)
    1190             :   {
    1191      660691 :     GEN s = gel(a,i);
    1192      660691 :     GEN x0 = gnorm(s);
    1193      660681 :     GEN x1 = gneg(gtrace(s));
    1194      660687 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1195             :   }
    1196      264464 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1197      264467 :   return gc_upto(av, normalized_to_RgX(L));
    1198             : }
    1199             : 
    1200             : GEN
    1201          56 : polfromroots(GEN a, long v)
    1202             : {
    1203          56 :   if (!is_vec_t(typ(a)))
    1204           0 :     pari_err_TYPE("polfromroots",a);
    1205          56 :   if (v < 0) v = 0;
    1206          56 :   if (varncmp(gvar(a), v) <= 0) pari_err_PRIORITY("polfromroots",a,"<=",v);
    1207          49 :   return roots_to_pol(a, v);
    1208             : }
    1209             : 
    1210             : /*******************************************************************/
    1211             : /*                                                                 */
    1212             : /*                          FACTORBACK                             */
    1213             : /*                                                                 */
    1214             : /*******************************************************************/
    1215             : static GEN
    1216    54631713 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
    1217             : static GEN
    1218    79122425 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
    1219             : static GEN
    1220    28781952 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
    1221             : static GEN
    1222      244383 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
    1223             : 
    1224             : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
    1225             : GEN
    1226    33804027 : gen_factorback(GEN L, GEN e, void *data, GEN (*_mul)(void*,GEN,GEN),
    1227             :                GEN (*_pow)(void*,GEN,GEN), GEN (*_one)(void*))
    1228             : {
    1229    33804027 :   pari_sp av = avma;
    1230             :   long k, l, lx;
    1231             :   GEN p,x;
    1232             : 
    1233    33804027 :   if (e) /* supplied vector of exponents */
    1234     1416592 :     p = L;
    1235             :   else
    1236             :   {
    1237    32387435 :     switch(typ(L)) {
    1238     8412124 :       case t_VEC:
    1239             :       case t_COL: /* product of the L[i] */
    1240     8412124 :         if (lg(L)==1) return _one? _one(data): gen_1;
    1241     8335610 :         return gc_upto(av, gen_product(L, data, _mul));
    1242    23975319 :       case t_MAT: /* genuine factorization */
    1243    23975319 :         l = lg(L);
    1244    23975319 :         if (l == 3) break;
    1245             :         /*fall through*/
    1246             :       default:
    1247           6 :         pari_err_TYPE("factorback [not a factorization]", L);
    1248             :     }
    1249    23975312 :     p = gel(L,1);
    1250    23975312 :     e = gel(L,2);
    1251             :   }
    1252    25391904 :   if (!is_vec_t(typ(p))) pari_err_TYPE("factorback [not a vector]", p);
    1253             :   /* p = elts, e = expo */
    1254    25391889 :   lx = lg(p);
    1255             :   /* check whether e is an integral vector of correct length */
    1256    25391889 :   switch(typ(e))
    1257             :   {
    1258      123198 :     case t_VECSMALL:
    1259      123198 :       if (lx != lg(e))
    1260           0 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1261      123198 :       if (lx == 1) return _one? _one(data): gen_1;
    1262      122939 :       x = cgetg(lx,t_VEC);
    1263     1249087 :       for (l=1,k=1; k<lx; k++)
    1264     1126148 :         if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
    1265      122939 :       break;
    1266    25268684 :     case t_VEC: case t_COL:
    1267    25268684 :       if (lx != lg(e) || !RgV_is_ZV(e))
    1268          14 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1269    25268670 :       if (lx == 1) return _one? _one(data): gen_1;
    1270    25158740 :       x = cgetg(lx,t_VEC);
    1271   105202573 :       for (l=1,k=1; k<lx; k++)
    1272    80043834 :         if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
    1273    25158739 :       break;
    1274           7 :     default:
    1275           7 :       pari_err_TYPE("factorback [not an exponent vector]", e);
    1276             :       return NULL;/*LCOV_EXCL_LINE*/
    1277             :   }
    1278    25281678 :   if (l==1) return gc_upto(av, _one? _one(data): gen_1);
    1279    25219427 :   x[0] = evaltyp(t_VEC) | _evallg(l);
    1280    25219427 :   return gc_upto(av, gen_product(x, data, _mul));
    1281             : }
    1282             : 
    1283             : GEN
    1284     8519895 : FpV_factorback(GEN L, GEN e, GEN p)
    1285     8519895 : { return gen_factorback(L, e, (void*)p, &Fpmul, &Fppow, NULL); }
    1286             : 
    1287             : ulong
    1288      108623 : Flv_factorback(GEN L, GEN e, ulong p)
    1289             : {
    1290      108623 :   long i, l = lg(e);
    1291      108623 :   ulong r = 1UL, ri = 1UL;
    1292      529663 :   for (i = 1; i < l; i++)
    1293             :   {
    1294      421040 :     long c = e[i];
    1295      421040 :     if (!c) continue;
    1296      178632 :     if (c < 0)
    1297           0 :       ri = Fl_mul(ri, Fl_powu(L[i],-c,p), p);
    1298             :     else
    1299      178632 :       r = Fl_mul(r, Fl_powu(L[i],c,p), p);
    1300             :   }
    1301      108623 :   if (ri != 1UL) r = Fl_div(r, ri, p);
    1302      108623 :   return r;
    1303             : }
    1304             : GEN
    1305        2499 : FlxqV_factorback(GEN L, GEN e, GEN Tp, ulong p)
    1306             : {
    1307        2499 :   pari_sp av = avma;
    1308        2499 :   GEN Hi = NULL, H = NULL;
    1309        2499 :   long i, l = lg(L), v = get_Flx_var(Tp);
    1310      168187 :   for (i = 1; i < l; i++)
    1311             :   {
    1312      165627 :     GEN x, ei = gel(e,i);
    1313      165627 :     long s = signe(ei);
    1314      165627 :     if (!s) continue;
    1315      157606 :     x = Flxq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
    1316      157608 :     if (s > 0)
    1317       79417 :       H = H? Flxq_mul(H, x, Tp, p): x;
    1318             :     else
    1319       78191 :       Hi = Hi? Flxq_mul(Hi, x, Tp, p): x;
    1320             :   }
    1321        2560 :   if (!Hi)
    1322             :   {
    1323           0 :     if (!H) { set_avma(av); return mkvecsmall2(v,1); }
    1324           0 :     return gc_uptoleaf(av, H);
    1325             :   }
    1326        2560 :   Hi = Flxq_inv(Hi, Tp, p);
    1327        2499 :   return gc_uptoleaf(av, H? Flxq_mul(H,Hi,Tp,p): Hi);
    1328             : }
    1329             : GEN
    1330          14 : FqV_factorback(GEN L, GEN e, GEN Tp, GEN p)
    1331             : {
    1332          14 :   pari_sp av = avma;
    1333          14 :   GEN Hi = NULL, H = NULL;
    1334          14 :   long i, l = lg(L), small = typ(e) == t_VECSMALL;
    1335        1554 :   for (i = 1; i < l; i++)
    1336             :   {
    1337             :     GEN x;
    1338             :     long s;
    1339        1540 :     if (small)
    1340             :     {
    1341           0 :       s = e[i]; if (!s) continue;
    1342           0 :       x = Fq_powu(gel(L,i), labs(s), Tp, p);
    1343             :     }
    1344             :     else
    1345             :     {
    1346        1540 :       GEN ei = gel(e,i);
    1347        1540 :       s = signe(ei); if (!s) continue;
    1348        1540 :       x = Fq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
    1349             :     }
    1350        1540 :     if (s > 0)
    1351         819 :       H = H? Fq_mul(H, x, Tp, p): x;
    1352             :     else
    1353         721 :       Hi = Hi? Fq_mul(Hi, x, Tp, p): x;
    1354             :   }
    1355          14 :   if (Hi)
    1356             :   {
    1357           7 :     Hi = Fq_inv(Hi, Tp, p);
    1358           7 :     H = H? Fq_mul(H,Hi,Tp,p): Hi;
    1359             :   }
    1360           7 :   else if (!H) return gc_const(av, gen_1);
    1361          14 :   return gc_upto(av, H);
    1362             : }
    1363             : 
    1364             : GEN
    1365    24951662 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, NULL, &mul, &powi, NULL); }
    1366             : GEN
    1367     1417189 : factorback(GEN fa) { return factorback2(fa, NULL); }
    1368             : 
    1369             : GEN
    1370       10003 : vecprod(GEN v)
    1371             : {
    1372       10003 :   pari_sp av = avma;
    1373       10003 :   if (!is_vec_t(typ(v)))
    1374           0 :     pari_err_TYPE("vecprod", v);
    1375       10003 :   if (lg(v) == 1) return gen_1;
    1376        9226 :   return gc_GEN(av, gen_product(v, NULL, mul));
    1377             : }
    1378             : 
    1379             : static int
    1380       11165 : RgX_is_irred_i(GEN x)
    1381             : {
    1382             :   GEN y, p, pol;
    1383       11165 :   long l = lg(x), pa;
    1384             : 
    1385       11165 :   if (!signe(x) || l <= 3) return 0;
    1386       11165 :   switch(RgX_type(x,&p,&pol,&pa))
    1387             :   {
    1388          21 :     case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
    1389           0 :     case t_COMPLEX: return l == 4;
    1390           0 :     case t_REAL:
    1391           0 :       if (l == 4) return 1;
    1392           0 :       if (l > 5) return 0;
    1393           0 :       return gsigne(RgX_disc(x)) > 0;
    1394             :   }
    1395       11144 :   y = RgX_factor(x, NULL);
    1396       11144 :   return (lg(gcoeff(y,1,1))==l);
    1397             : }
    1398             : static int
    1399       11165 : RgX_is_irred(GEN x)
    1400       11165 : { pari_sp av = avma; return gc_bool(av, RgX_is_irred_i(x)); }
    1401             : long
    1402       11165 : polisirreducible(GEN x)
    1403             : {
    1404       11165 :   long tx = typ(x);
    1405       11165 :   if (tx == t_POL) return RgX_is_irred(x);
    1406           0 :   if (!is_scalar_t(tx)) pari_err_TYPE("polisirreducible",x);
    1407           0 :   return 0;
    1408             : }
    1409             : 
    1410             : /*******************************************************************/
    1411             : /*                                                                 */
    1412             : /*                         GENERIC GCD                             */
    1413             : /*                                                                 */
    1414             : /*******************************************************************/
    1415             : /* x is a COMPLEX or a QUAD */
    1416             : static GEN
    1417        3234 : triv_cont_gcd(GEN x, GEN y)
    1418             : {
    1419        3234 :   pari_sp av = avma;
    1420             :   GEN c;
    1421        3234 :   if (typ(x)==t_COMPLEX)
    1422             :   {
    1423        2912 :     GEN a = gel(x,1), b = gel(x,2);
    1424        2912 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
    1425          21 :     c = ggcd(a,b);
    1426             :   }
    1427             :   else
    1428         322 :     c = ggcd(gel(x,2),gel(x,3));
    1429         343 :   return gc_upto(av, ggcd(c,y));
    1430             : }
    1431             : 
    1432             : /* y is a PADIC, x a rational number or an INTMOD */
    1433             : static GEN
    1434        1407 : padic_gcd(GEN x, GEN y)
    1435             : {
    1436        1407 :   GEN p = padic_p(y);
    1437        1407 :   long v = gvaluation(x,p), w = valp(y);
    1438        1407 :   if (w < v) v = w;
    1439        1407 :   return powis(p, v);
    1440             : }
    1441             : 
    1442             : static void
    1443         896 : Zi_mul3(GEN xr, GEN xi, GEN yr, GEN yi, GEN *zr, GEN *zi)
    1444             : {
    1445         896 :   GEN p3 = addii(xr,xi);
    1446         896 :   GEN p4 = addii(yr,yi);
    1447         896 :   GEN p1 = mulii(xr,yr);
    1448         896 :   GEN p2 = mulii(xi,yi);
    1449         896 :   p3 = mulii(p3,p4);
    1450         896 :   p4 = addii(p2,p1);
    1451         896 :   *zr = subii(p1,p2); *zi = subii(p3,p4);
    1452         896 : }
    1453             : 
    1454             : static GEN
    1455         448 : Zi_rem(GEN x, GEN y)
    1456             : {
    1457         448 :   GEN xr = real_i(x), xi = imag_i(x);
    1458         448 :   GEN yr = real_i(y), yi = imag_i(y);
    1459         448 :   GEN n = addii(sqri(yr), sqri(yi));
    1460             :   GEN ur, ui, zr, zi;
    1461         448 :   Zi_mul3(xr, xi, yr, negi(yi), &ur, &ui);
    1462         448 :   Zi_mul3(yr, yi, diviiround(ur, n), diviiround(ui, n), &zr, &zi);
    1463         448 :   return mkcomplex(subii(xr,zr), subii(xi,zi));
    1464             : }
    1465             : 
    1466             : static GEN
    1467         399 : Qi_gcd(GEN x, GEN y)
    1468             : {
    1469         399 :   pari_sp av = avma, btop;
    1470             :   GEN dx, dy;
    1471         399 :   x = Q_remove_denom(x, &dx);
    1472         399 :   y = Q_remove_denom(y, &dy);
    1473         399 :   btop = avma;
    1474         847 :   while (!gequal0(y))
    1475             :   {
    1476         448 :     GEN z = Zi_rem(x,y);
    1477         448 :     x = y; y = z;
    1478         448 :     if (gc_needed(btop,1)) {
    1479           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Qi_gcd");
    1480           0 :       (void)gc_all(btop,2, &x,&y);
    1481             :     }
    1482             :   }
    1483         399 :   x = Qi_normal(x);
    1484         399 :   if (typ(x) == t_COMPLEX)
    1485             :   {
    1486         280 :     if      (gequal0(gel(x,2))) x = gel(x,1);
    1487         203 :     else if (gequal0(gel(x,1))) x = gel(x,2);
    1488             :   }
    1489         399 :   if (!dx && !dy) return gc_GEN(av, x);
    1490          35 :   return gc_upto(av, gdiv(x, dx? (dy? lcmii(dx, dy): dx): dy));
    1491             : }
    1492             : 
    1493             : static int
    1494        3269 : c_is_rational(GEN x)
    1495        3269 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
    1496             : static GEN
    1497        1386 : c_zero_gcd(GEN c)
    1498             : {
    1499        1386 :   GEN x = gel(c,1), y = gel(c,2);
    1500        1386 :   long tx = typ(x), ty = typ(y);
    1501        1386 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
    1502          42 :   if (tx == t_PADIC || tx == t_INTMOD
    1503          42 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
    1504          35 :   return Qi_gcd(c, gen_0);
    1505             : }
    1506             : 
    1507             : /* gcd(x, 0) */
    1508             : static GEN
    1509     8258205 : zero_gcd(GEN x)
    1510             : {
    1511             :   pari_sp av;
    1512     8258205 :   switch(typ(x))
    1513             :   {
    1514       47217 :     case t_INT: return absi(x);
    1515       46646 :     case t_FRAC: return absfrac(x);
    1516        1386 :     case t_COMPLEX: return c_zero_gcd(x);
    1517         707 :     case t_REAL: return gen_1;
    1518         728 :     case t_PADIC: return powis(padic_p(x), valp(x));
    1519         252 :     case t_SER: return pol_xnall(valser(x), varn(x));
    1520        3363 :     case t_POLMOD: {
    1521        3363 :       GEN d = gel(x,2);
    1522        3363 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
    1523         476 :       return isinexact(d)? zero_gcd(d): gcopy(d);
    1524             :     }
    1525     7913273 :     case t_POL:
    1526     7913273 :       if (!isinexact(x)) break;
    1527          14 :       av = avma;
    1528          14 :       return gc_upto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
    1529             : 
    1530      217452 :     case t_RFRAC:
    1531      217452 :       if (!isinexact(x)) break;
    1532           0 :       av = avma;
    1533           0 :       return gc_upto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
    1534             :   }
    1535     8157892 :   return gcopy(x);
    1536             : }
    1537             : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
    1538             : static GEN
    1539     8879691 : zero_gcd2(GEN y, GEN z)
    1540             : {
    1541             :   pari_sp av;
    1542     8879691 :   switch(typ(z))
    1543             :   {
    1544     8238708 :     case t_INT: return zero_gcd(y);
    1545      637476 :     case t_INTMOD:
    1546      637476 :       av = avma;
    1547      637476 :       return gc_upto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
    1548        3507 :     case t_FFELT:
    1549        3507 :       av = avma;
    1550        3507 :       return gc_upto(av, gmul(y, FF_1(z)));
    1551           0 :     default:
    1552           0 :       pari_err_TYPE("zero_gcd", z);
    1553             :       return NULL;/*LCOV_EXCL_LINE*/
    1554             :   }
    1555             : }
    1556             : static GEN
    1557     1912226 : cont_gcd_pol_i(GEN x, GEN y) { return scalarpol(ggcd(content(x),y), varn(x));}
    1558             : /* tx = t_POL, y considered as constant */
    1559             : static GEN
    1560     1912226 : cont_gcd_pol(GEN x, GEN y)
    1561     1912226 : { pari_sp av = avma; return gc_upto(av, cont_gcd_pol_i(x,y)); }
    1562             : /* tx = t_RFRAC, y considered as constant */
    1563             : static GEN
    1564       10115 : cont_gcd_rfrac(GEN x, GEN y)
    1565             : {
    1566       10115 :   pari_sp av = avma;
    1567       10115 :   GEN cx; x = primitive_part(x, &cx);
    1568             :   /* e.g. Mod(1,2) / (2*y+1) => primitive_part = Mod(1,2)*y^0 */
    1569       10115 :   if (typ(x) != t_RFRAC) x = cont_gcd_pol_i(x, y);
    1570       10115 :   else x = gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2));
    1571       10115 :   return gc_upto(av, x);
    1572             : }
    1573             : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
    1574             : static GEN
    1575        2635 : cont_gcd_gen(GEN x, GEN y)
    1576             : {
    1577        2635 :   pari_sp av = avma;
    1578        2635 :   return gc_upto(av, ggcd(content(x),y));
    1579             : }
    1580             : /* !is_const(tx), y considered as constant */
    1581             : static GEN
    1582     1924955 : cont_gcd(GEN x, long tx, GEN y)
    1583             : {
    1584     1924955 :   switch(tx)
    1585             :   {
    1586       10115 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1587     1912205 :     case t_POL: return cont_gcd_pol(x,y);
    1588        2635 :     default: return cont_gcd_gen(x,y);
    1589             :   }
    1590             : }
    1591             : static GEN
    1592    12569788 : gcdiq(GEN x, GEN y)
    1593             : {
    1594             :   GEN z;
    1595    12569788 :   if (!signe(x)) return Q_abs(y);
    1596     4662586 :   z = cgetg(3,t_FRAC);
    1597     4662612 :   gel(z,1) = gcdii(x,gel(y,1));
    1598     4662585 :   gel(z,2) = icopy(gel(y,2));
    1599     4662581 :   return z;
    1600             : }
    1601             : static GEN
    1602    27545103 : gcdqq(GEN x, GEN y)
    1603             : {
    1604    27545103 :   GEN z = cgetg(3,t_FRAC);
    1605    27545090 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1606    27544910 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1607    27545003 :   return z;
    1608             : }
    1609             : /* assume x,y t_INT or t_FRAC */
    1610             : GEN
    1611  1007365579 : Q_gcd(GEN x, GEN y)
    1612             : {
    1613  1007365579 :   long tx = typ(x), ty = typ(y);
    1614  1007365579 :   if (tx == t_INT)
    1615   970252708 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1616             :   else
    1617    37112871 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1618             : }
    1619             : static GEN
    1620        6005 : gcd3(GEN x, GEN y, GEN z) { return ggcd(ggcd(x, y), z); }
    1621             : 
    1622             : GEN
    1623    27700544 : ggcd(GEN x, GEN y)
    1624             : {
    1625    27700544 :   long vx, vy, tx = typ(x), ty = typ(y);
    1626             :   pari_sp av, tetpil;
    1627             :   GEN p1,z;
    1628             : 
    1629    55401088 :   if (is_noncalc_t(tx) || is_matvec_t(tx) ||
    1630    55401088 :       is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
    1631    27700544 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1632             :   /* tx <= ty */
    1633    27700544 :   z = gisexactzero(x); if (z) return zero_gcd2(y,z);
    1634    24612124 :   z = gisexactzero(y); if (z) return zero_gcd2(x,z);
    1635    18820853 :   if (is_const_t(tx))
    1636             :   {
    1637    11701654 :     if (ty == tx) switch(tx)
    1638             :     {
    1639     7272469 :       case t_INT:
    1640     7272469 :         return gcdii(x,y);
    1641             : 
    1642     2178029 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1643     2178029 :         if (equalii(gel(x,1),gel(y,1)))
    1644     2178022 :           gel(z,1) = icopy(gel(x,1));
    1645             :         else
    1646           7 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1647     2178029 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1648             :         else
    1649             :         {
    1650     2178029 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1651     2178029 :           if (!equali1(p1))
    1652             :           {
    1653           7 :             p1 = gcdii(p1,gel(y,2));
    1654           7 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1655           7 :             else p1 = gc_INT(av, p1);
    1656             :           }
    1657     2178029 :           gel(z,2) = p1;
    1658             :         }
    1659     2178029 :         return z;
    1660             : 
    1661      261957 :       case t_FRAC:
    1662      261957 :         return gcdqq(x,y);
    1663             : 
    1664        5845 :       case t_FFELT:
    1665        5845 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1666        5845 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1667             : 
    1668          21 :       case t_COMPLEX:
    1669          21 :         if (c_is_rational(x) && c_is_rational(y)) return Qi_gcd(x,y);
    1670           7 :         return triv_cont_gcd(y,x);
    1671             : 
    1672          14 :       case t_PADIC:
    1673          14 :         if (!equalii(padic_p(x), padic_p(y))) return gen_1;
    1674           7 :         return powis(padic_p(x), minss(valp(x), valp(y)));
    1675             : 
    1676         133 :       case t_QUAD:
    1677         133 :         av=avma; p1=gdiv(x,y);
    1678         133 :         if (gequal0(gel(p1,3)))
    1679             :         {
    1680          14 :           p1=gel(p1,2);
    1681          14 :           if (typ(p1)==t_INT) { set_avma(av); return gcopy(y); }
    1682           7 :           tetpil=avma; return gc_GEN_unsafe(av,tetpil, gdiv(y,gel(p1,2)));
    1683             :         }
    1684         119 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {set_avma(av); return gcopy(y);}
    1685         112 :         p1 = ginv(p1); set_avma(av);
    1686         112 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1687         105 :         return triv_cont_gcd(y,x);
    1688             : 
    1689           0 :       default: return gen_1; /* t_REAL */
    1690             :     }
    1691     1983186 :     if (is_const_t(ty)) switch(tx)
    1692             :     {
    1693       75269 :       case t_INT:
    1694       75269 :         switch(ty)
    1695             :         {
    1696          63 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1697          63 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1698          63 :             p1 = gcdii(gel(y,1),gel(y,2));
    1699          63 :             if (!equali1(p1)) {
    1700          14 :               p1 = gcdii(x,p1);
    1701          14 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1702             :               else
    1703          14 :                 p1 = gc_INT(av, p1);
    1704             :             }
    1705          63 :             gel(z,2) = p1; return z;
    1706             : 
    1707        8722 :           case t_REAL: return gen_1;
    1708             : 
    1709       61752 :           case t_FRAC:
    1710       61752 :             return gcdiq(x,y);
    1711             : 
    1712        3143 :           case t_COMPLEX:
    1713        3143 :             if (c_is_rational(y)) return Qi_gcd(x,y);
    1714        2800 :             return triv_cont_gcd(y,x);
    1715             : 
    1716          70 :           case t_FFELT:
    1717          70 :             if (!FF_equal0(y)) return FF_1(y);
    1718           0 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1719             : 
    1720        1393 :           case t_PADIC:
    1721        1393 :             return padic_gcd(x,y);
    1722             : 
    1723         126 :           case t_QUAD:
    1724         126 :             return triv_cont_gcd(y,x);
    1725           0 :           default:
    1726           0 :             pari_err_TYPE2("gcd",x,y);
    1727             :         }
    1728             : 
    1729          14 :       case t_REAL:
    1730          14 :         switch(ty)
    1731             :         {
    1732          14 :           case t_INTMOD:
    1733             :           case t_FFELT:
    1734          14 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1735           0 :           default: return gen_1;
    1736             :         }
    1737             : 
    1738          49 :       case t_INTMOD:
    1739          49 :         switch(ty)
    1740             :         {
    1741          14 :           case t_FRAC:
    1742          14 :             av = avma;
    1743          14 :             if (!equali1(gcdii(gel(x,1),gel(y,2)))) pari_err_OP("gcd",x,y);
    1744           7 :             set_avma(av); return ggcd(gel(y,1), x);
    1745             : 
    1746          14 :           case t_FFELT:
    1747             :           {
    1748          14 :             GEN p = gel(y,4);
    1749          14 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1750           7 :             if (!FF_equal0(y)) return FF_1(y);
    1751           0 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1752             :           }
    1753             : 
    1754          14 :           case t_COMPLEX: case t_QUAD:
    1755          14 :             return triv_cont_gcd(y,x);
    1756             : 
    1757           7 :           case t_PADIC:
    1758           7 :             return padic_gcd(x,y);
    1759             : 
    1760           0 :           default: pari_err_TYPE2("gcd",x,y);
    1761             :         }
    1762             : 
    1763         217 :       case t_FRAC:
    1764         217 :         switch(ty)
    1765             :         {
    1766          91 :           case t_COMPLEX:
    1767          91 :             if (c_is_rational(y)) return Qi_gcd(x,y);
    1768             :           case t_QUAD:
    1769         161 :             return triv_cont_gcd(y,x);
    1770          42 :           case t_FFELT:
    1771             :           {
    1772          42 :             GEN p = gel(y,4);
    1773          42 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1774          21 :             if (!FF_equal0(y)) return FF_1(y);
    1775           0 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1776             :           }
    1777             : 
    1778           7 :           case t_PADIC:
    1779           7 :             return padic_gcd(x,y);
    1780             : 
    1781           0 :           default: pari_err_TYPE2("gcd",x,y);
    1782             :         }
    1783          70 :       case t_FFELT:
    1784          70 :         switch(ty)
    1785             :         {
    1786          42 :           case t_PADIC:
    1787             :           {
    1788          42 :             GEN p = padic_p(y);
    1789          42 :             long v = valp(y);
    1790          42 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1791          14 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1792             :           }
    1793          28 :           default: pari_err_TYPE2("gcd",x,y);
    1794             :         }
    1795             : 
    1796          14 :       case t_COMPLEX:
    1797          14 :         switch(ty)
    1798             :         {
    1799          14 :           case t_PADIC:
    1800          14 :           case t_QUAD: return triv_cont_gcd(x,y);
    1801           0 :           default: pari_err_TYPE2("gcd",x,y);
    1802             :         }
    1803             : 
    1804           7 :       case t_PADIC:
    1805           7 :         switch(ty)
    1806             :         {
    1807           7 :           case t_QUAD: return triv_cont_gcd(y,x);
    1808           0 :           default: pari_err_TYPE2("gcd",x,y);
    1809             :         }
    1810             : 
    1811           0 :       default: return gen_1; /* tx = t_REAL */
    1812             :     }
    1813     1907546 :     return cont_gcd(y,ty, x);
    1814             :   }
    1815             : 
    1816     7119199 :   if (tx == t_POLMOD)
    1817             :   {
    1818        6054 :     if (ty == t_POLMOD)
    1819             :     {
    1820        5977 :       GEN T = gel(x,1), Ty = gel(y,1);
    1821        5977 :       vx = varn(T); vy = varn(Ty);
    1822        5977 :       z = cgetg(3,t_POLMOD);
    1823        5977 :       if (vx == vy)
    1824        5963 :         T = RgX_equal(T,Ty)? RgX_copy(T): RgX_gcd(T, Ty);
    1825             :       else
    1826          14 :         T = RgX_copy(varncmp(vx,vy) < 0? T: Ty);
    1827        5977 :       gel(z,1) = T;
    1828        5977 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1829             :       else
    1830             :       {
    1831        5977 :         GEN X = gel(x,2), Y = gel(y,2), d;
    1832        5977 :         av = avma; d = ggcd(content(X), content(Y));
    1833        5977 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1834        5977 :         gel(z,2) = gc_upto(av, gmul(d, gcd3(T, X, Y)));
    1835             :       }
    1836        5977 :       return z;
    1837             :     }
    1838          77 :     vx = varn(gel(x,1));
    1839          77 :     switch(ty)
    1840             :     {
    1841          49 :       case t_POL:
    1842          49 :         vy = varn(y);
    1843          49 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1844          28 :         z = cgetg(3,t_POLMOD);
    1845          28 :         gel(z,1) = RgX_copy(gel(x,1)); av = avma;
    1846          28 :         gel(z,2) = gc_upto(av, gcd3(gel(x,1), gel(x,2), y));
    1847          28 :         return z;
    1848             : 
    1849          28 :       case t_RFRAC:
    1850          28 :         vy = varn(gel(y,2));
    1851          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1852          28 :         av = avma;
    1853          28 :         if (degpol(ggcd(gel(x,1),gel(y,2)))) pari_err_OP("gcd",x,y);
    1854          21 :         set_avma(av); return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1855             :     }
    1856             :   }
    1857             : 
    1858     7113145 :   vx = gvar(x);
    1859     7113145 :   vy = gvar(y);
    1860     7113145 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1861     7106264 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1862             : 
    1863             :   /* vx = vy: same main variable */
    1864     7095736 :   switch(tx)
    1865             :   {
    1866     6934732 :     case t_POL:
    1867             :       switch(ty)
    1868             :       {
    1869             :         GEN cz, cx, cy;
    1870     6080832 :         case t_POL: return RgX_gcd(x,y);
    1871          28 :         case t_SER:
    1872          28 :           z = ggcd(content(x), content(y));
    1873          28 :           return monomialcopy(z, minss(valser(y),gval(x,vx)), vx);
    1874      853872 :         case t_RFRAC:
    1875      853872 :           av = avma;
    1876      853872 :           x = primitive_part(x, &cx);
    1877      853872 :           y = primitive_part(y, &cy);
    1878      853872 :           z = gred_rfrac_simple(ggcd(gel(y,1), x), gel(y,2));
    1879      853872 :           if (cx) cz = cy? ggcd(cx, cy): cx; else cz = cy? cy: NULL;
    1880      853872 :           if (cz) z = gmul(z, cz);
    1881      853872 :           return gc_upto(av, z);
    1882             :       }
    1883           0 :       break;
    1884             : 
    1885          14 :     case t_SER:
    1886          14 :       z = ggcd(content(x), content(y));
    1887             :       switch(ty)
    1888             :       {
    1889           7 :         case t_SER:   return monomialcopy(z, minss(valser(x),valser(y)), vx);
    1890           7 :         case t_RFRAC: return monomialcopy(z, minss(valser(x),gval(y,vx)), vx);
    1891             :       }
    1892           0 :       break;
    1893             : 
    1894      160990 :     case t_RFRAC:
    1895             :     {
    1896      160990 :       GEN xd = gel(x,2), yd = gel(y,2);
    1897      160990 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1898      160990 :       z = cgetg(3,t_RFRAC); av = avma;
    1899      160990 :       gel(z,2) = gc_upto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1900      160990 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1901             :     }
    1902             :   }
    1903           0 :   pari_err_TYPE2("gcd",x,y);
    1904             :   return NULL; /* LCOV_EXCL_LINE */
    1905             : }
    1906             : GEN
    1907     5203623 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
    1908             : 
    1909             : static GEN
    1910        3680 : fix_lcm(GEN x)
    1911             : {
    1912             :   GEN t;
    1913        3680 :   switch(typ(x))
    1914             :   {
    1915        3575 :     case t_INT:
    1916        3575 :       x = absi_shallow(x); break;
    1917          98 :     case t_POL:
    1918          98 :       if (lg(x) <= 2) break;
    1919          98 :       t = leading_coeff(x);
    1920          98 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1921             :   }
    1922        3680 :   return x;
    1923             : }
    1924             : GEN
    1925        2898 : glcm0(GEN x, GEN y)
    1926             : {
    1927        2898 :   if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
    1928        2849 :   return glcm(x,y);
    1929             : }
    1930             : GEN
    1931        3575 : ZV_lcm(GEN x) { return fix_lcm(gassoc_proto(lcmii,x,NULL)); }
    1932             : GEN
    1933        3283 : glcm(GEN x, GEN y)
    1934             : {
    1935             :   pari_sp av;
    1936             :   GEN z;
    1937        3283 :   if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
    1938          70 :   av = avma; z = ggcd(x,y);
    1939          70 :   if (!gequal1(z))
    1940             :   {
    1941          63 :     if (gequal0(z)) { set_avma(av); return gmul(x,y); }
    1942          49 :     y = gdiv(y,z);
    1943             :   }
    1944          56 :   return gc_upto(av, fix_lcm(gmul(x,y)));
    1945             : }
    1946             : 
    1947             : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
    1948             : static int
    1949           0 : pol_approx0(GEN r, GEN x, int exact)
    1950             : {
    1951             :   long i, l;
    1952           0 :   if (exact) return !signe(r);
    1953           0 :   l = minss(lg(x), lg(r));
    1954           0 :   for (i = 2; i < l; i++)
    1955           0 :     if (!cx_approx0(gel(r,i), gel(x,i))) return 0;
    1956           0 :   return 1;
    1957             : }
    1958             : 
    1959             : GEN
    1960           0 : RgX_gcd_simple(GEN x, GEN y)
    1961             : {
    1962           0 :   pari_sp av1, av = avma;
    1963           0 :   GEN r, yorig = y;
    1964           0 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1965             : 
    1966             :   for(;;)
    1967             :   {
    1968           0 :     av1 = avma; r = RgX_rem(x,y);
    1969           0 :     if (pol_approx0(r, x, exact))
    1970             :     {
    1971           0 :       set_avma(av1);
    1972           0 :       if (y == yorig) return RgX_copy(y);
    1973           0 :       y = normalizepol_approx(y, lg(y));
    1974           0 :       if (lg(y) == 3) { set_avma(av); return pol_1(varn(x)); }
    1975           0 :       return gc_upto(av,y);
    1976             :     }
    1977           0 :     x = y; y = r;
    1978           0 :     if (gc_needed(av,1)) {
    1979           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1980           0 :       (void)gc_all(av,2, &x,&y);
    1981             :     }
    1982             :   }
    1983             : }
    1984             : GEN
    1985           0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
    1986             : {
    1987           0 :   pari_sp av = avma;
    1988             :   GEN q, r, d, d1, u, v, v1;
    1989           0 :   int exact = !(isinexactreal(a) || isinexactreal(b));
    1990             : 
    1991           0 :   d = a; d1 = b; v = gen_0; v1 = gen_1;
    1992             :   for(;;)
    1993             :   {
    1994           0 :     if (pol_approx0(d1, a, exact)) break;
    1995           0 :     q = poldivrem(d,d1, &r);
    1996           0 :     v = gsub(v, gmul(q,v1));
    1997           0 :     u=v; v=v1; v1=u;
    1998           0 :     u=r; d=d1; d1=u;
    1999             :   }
    2000           0 :   u = gsub(d, gmul(b,v));
    2001           0 :   u = RgX_div(u,a);
    2002             : 
    2003           0 :   (void)gc_all(av, 3, &u,&v,&d);
    2004           0 :   *pu = u;
    2005           0 :   *pv = v; return d;
    2006             : }
    2007             : 
    2008             : GEN
    2009          91 : ghalfgcd(GEN x, GEN y)
    2010             : {
    2011          91 :   long tx = typ(x), ty = typ(y);
    2012          91 :   if (tx==t_INT && ty==t_INT) return halfgcdii(x, y);
    2013          63 :   if (tx==t_POL && ty==t_POL && varn(x)==varn(y))
    2014             :   {
    2015          63 :     pari_sp av = avma;
    2016          63 :     GEN a, b, M = RgX_halfgcd_all(x, y, &a, &b);
    2017          63 :     return gc_GEN(av, mkvec2(M, mkcol2(a,b)));
    2018             :   }
    2019           0 :   pari_err_OP("halfgcd", x, y);
    2020             :   return NULL; /* LCOV_EXCL_LINE */
    2021             : }
    2022             : 
    2023             : /*******************************************************************/
    2024             : /*                                                                 */
    2025             : /*                    CONTENT / PRIMITIVE PART                     */
    2026             : /*                                                                 */
    2027             : /*******************************************************************/
    2028             : 
    2029             : GEN
    2030    72467270 : content(GEN x)
    2031             : {
    2032    72467270 :   long lx, i, t, tx = typ(x);
    2033    72467270 :   pari_sp av = avma;
    2034             :   GEN c;
    2035             : 
    2036    72467270 :   if (is_scalar_t(tx)) return zero_gcd(x);
    2037    72450786 :   switch(tx)
    2038             :   {
    2039      864022 :     case t_RFRAC:
    2040             :     {
    2041      864022 :       GEN n = gel(x,1), d = gel(x,2);
    2042             :       /* -- varncmp(vn, vd) < 0 can't happen
    2043             :        * -- if n is POLMOD, its main variable (in the sense of gvar2)
    2044             :        *    has lower priority than denominator */
    2045      864022 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    2046      823598 :         n = isinexact(n)? zero_gcd(n): gcopy(n);
    2047             :       else
    2048       40424 :         n = content(n);
    2049      864022 :       return gc_upto(av, gdiv(n, content(d)));
    2050             :     }
    2051             : 
    2052     1117813 :     case t_VEC: case t_COL:
    2053     1117813 :       lx = lg(x); if (lx==1) return gen_0;
    2054     1117806 :       break;
    2055             : 
    2056          21 :     case t_MAT:
    2057             :     {
    2058             :       long hx, j;
    2059          21 :       lx = lg(x);
    2060          21 :       if (lx == 1) return gen_0;
    2061          14 :       hx = lgcols(x);
    2062          14 :       if (hx == 1) return gen_0;
    2063           7 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    2064           7 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    2065           7 :       c = content(gel(x,1));
    2066          14 :       for (j=2; j<lx; j++)
    2067          21 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    2068           7 :       if (typ(c) == t_INTMOD || isinexact(c)) return gc_const(av, gen_1);
    2069           7 :       return gc_upto(av,c);
    2070             :     }
    2071             : 
    2072    70468720 :     case t_POL: case t_SER:
    2073    70468720 :       lx = lg(x); if (lx == 2) return gen_0;
    2074    70445529 :       break;
    2075          21 :     case t_VECSMALL: return utoi(zv_content(x));
    2076         189 :     case t_QFB:
    2077         189 :       lx = 4; break;
    2078             : 
    2079           0 :     default: pari_err_TYPE("content",x);
    2080             :       return NULL; /* LCOV_EXCL_LINE */
    2081             :   }
    2082   213555818 :   for (i=lontyp[tx]; i<lx; i++)
    2083   153354943 :     if (typ(gel(x,i)) != t_INT) break;
    2084    71563524 :   lx--; c = gel(x,lx);
    2085    71563524 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    2086    71563526 :   if (i > lx)
    2087             :   { /* integer coeffs */
    2088    63070231 :     while (lx-- > lontyp[tx])
    2089             :     {
    2090    60793603 :       c = gcdii(c, gel(x,lx));
    2091    60793573 :       if (equali1(c)) return gc_const(av, gen_1);
    2092             :     }
    2093             :   }
    2094             :   else
    2095             :   {
    2096    11362649 :     if (isinexact(c)) c = zero_gcd(c);
    2097    30896912 :     while (lx-- > lontyp[tx])
    2098             :     {
    2099    19534263 :       GEN d = gel(x,lx);
    2100    19534263 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    2101    19534263 :       c = ggcd(c, d);
    2102             :     }
    2103    11362649 :     if (isinexact(c)) return gc_const(av, gen_1);
    2104             :   }
    2105    13639277 :   switch(typ(c))
    2106             :   {
    2107     2281663 :     case t_INT:
    2108     2281663 :       c = absi_shallow(c); break;
    2109           0 :     case t_VEC: case t_COL: case t_MAT:
    2110           0 :       pari_err_TYPE("content",x);
    2111             :   }
    2112             : 
    2113    13639284 :   return av==avma? gcopy(c): gc_upto(av,c);
    2114             : }
    2115             : 
    2116             : GEN
    2117     2800440 : primitive_part(GEN x, GEN *ptc)
    2118             : {
    2119     2800440 :   pari_sp av = avma;
    2120     2800440 :   GEN c = content(x);
    2121     2800402 :   if (gequal1(c)) { set_avma(av); c = NULL; }
    2122      203210 :   else if (!gequal0(c)) x = gdiv(x,c);
    2123     2800409 :   if (ptc) *ptc = c;
    2124     2800409 :   return x;
    2125             : }
    2126             : GEN
    2127         161 : primpart(GEN x) { return primitive_part(x, NULL); }
    2128             : 
    2129             : static GEN
    2130   179082002 : Q_content_v(GEN x, long imin, long l)
    2131             : {
    2132   179082002 :   pari_sp av = avma;
    2133   179082002 :   long i = l-1;
    2134   179082002 :   GEN d = Q_content_safe(gel(x,i));
    2135   179085849 :   if (!d) return NULL;
    2136  1186226079 :   for (i--; i >= imin; i--)
    2137             :   {
    2138  1007250112 :     GEN c = Q_content_safe(gel(x,i));
    2139  1007340909 :     if (!c) return NULL;
    2140  1007340867 :     d = Q_gcd(d, c);
    2141  1007143443 :     if (gc_needed(av,1)) d = gc_upto(av, d);
    2142             :   }
    2143   178975967 :   return gc_upto(av, d);
    2144             : }
    2145             : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    2146             :  * of Q(x2,...,xn)[x1] */
    2147             : GEN
    2148  1276547661 : Q_content_safe(GEN x)
    2149             : {
    2150             :   long l;
    2151  1276547661 :   switch(typ(x))
    2152             :   {
    2153  1058769978 :     case t_INT:  return absi(x);
    2154    37594427 :     case t_FRAC: return absfrac(x);
    2155   126325274 :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2156   126325274 :       l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
    2157    53901206 :     case t_POL:
    2158    53901206 :       l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
    2159       32687 :     case t_POLMOD: return Q_content_safe(gel(x,2));
    2160          21 :     case t_RFRAC:
    2161             :     {
    2162             :       GEN a, b;
    2163          21 :       a = Q_content_safe(gel(x,1)); if (!a) return NULL;
    2164          21 :       b = Q_content_safe(gel(x,2)); if (!b) return NULL;
    2165          21 :       return gdiv(a, b);
    2166             :     }
    2167             :   }
    2168         330 :   return NULL;
    2169             : }
    2170             : GEN
    2171     1844173 : Q_content(GEN x)
    2172             : {
    2173     1844173 :   GEN c = Q_content_safe(x);
    2174     1844171 :   if (!c) pari_err_TYPE("Q_content",x);
    2175     1844171 :   return c;
    2176             : }
    2177             : 
    2178             : GEN
    2179       13146 : ZX_content(GEN x)
    2180             : {
    2181       13146 :   long i, l = lg(x);
    2182             :   GEN d;
    2183             :   pari_sp av;
    2184             : 
    2185       13146 :   if (l == 2) return gen_0;
    2186       13146 :   d = gel(x,2);
    2187       13146 :   if (l == 3) return absi(d);
    2188        9170 :   av = avma;
    2189       18963 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    2190        9170 :   if (signe(d) < 0) d = negi(d);
    2191        9170 :   return gc_INT(av, d);
    2192             : }
    2193             : 
    2194             : static GEN
    2195     2282114 : Z_content_v(GEN x, long i, long l)
    2196             : {
    2197     2282114 :   pari_sp av = avma;
    2198     2282114 :   GEN d = Z_content(gel(x,i));
    2199     2282112 :   if (!d) return NULL;
    2200     5878976 :   for (i++; i<l; i++)
    2201             :   {
    2202     5359095 :     GEN c = Z_content(gel(x,i));
    2203     5359164 :     if (!c) return NULL;
    2204     4743574 :     d = gcdii(d, c); if (equali1(d)) return NULL;
    2205     3907693 :     if ((i & 255) == 0) d = gc_INT(av, d);
    2206             :   }
    2207      519881 :   return gc_INT(av, d);
    2208             : }
    2209             : /* return NULL for 1 */
    2210             : GEN
    2211     9881544 : Z_content(GEN x)
    2212             : {
    2213             :   long l;
    2214     9881544 :   switch(typ(x))
    2215             :   {
    2216     7580438 :     case t_INT:
    2217     7580438 :       if (is_pm1(x)) return NULL;
    2218     6657848 :       return absi(x);
    2219     2247765 :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2220     2247765 :       l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
    2221       53431 :     case t_POL:
    2222       53431 :       l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
    2223           0 :     case t_POLMOD: return Z_content(gel(x,2));
    2224             :   }
    2225           0 :   pari_err_TYPE("Z_content", x);
    2226             :   return NULL; /* LCOV_EXCL_LINE */
    2227             : }
    2228             : 
    2229             : static GEN
    2230    53793038 : Q_denom_v(GEN x, long i, long l)
    2231             : {
    2232    53793038 :   pari_sp av = avma;
    2233    53793038 :   GEN d = Q_denom_safe(gel(x,i));
    2234    53792712 :   if (!d) return NULL;
    2235   187079770 :   for (i++; i<l; i++)
    2236             :   {
    2237   133287127 :     GEN D = Q_denom_safe(gel(x,i));
    2238   133286935 :     if (!D) return NULL;
    2239   133286935 :     if (D != gen_1) d = lcmii(d, D);
    2240   133286801 :     if ((i & 255) == 0) d = gc_INT(av, d);
    2241             :   }
    2242    53792643 :   return gc_INT(av, d);
    2243             : }
    2244             : /* NOT MEMORY CLEAN (because of t_FRAC).
    2245             :  * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    2246             :  * of Q(x2,...,xn)[x1] */
    2247             : GEN
    2248   245873689 : Q_denom_safe(GEN x)
    2249             : {
    2250             :   long l;
    2251   245873689 :   switch(typ(x))
    2252             :   {
    2253   158093963 :     case t_INT: return gen_1;
    2254          28 :     case t_PADIC: l = valp(x); return l < 0? powiu(padic_p(x), -l): gen_1;
    2255    33723933 :     case t_FRAC: return gel(x,2);
    2256         504 :     case t_QUAD: return Q_denom_v(x, 2, 4);
    2257    41441937 :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2258    41441937 :       l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
    2259    12517725 :     case t_POL: case t_SER:
    2260    12517725 :       l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
    2261       93401 :     case t_POLMOD: return Q_denom(gel(x,2));
    2262        8134 :     case t_RFRAC:
    2263             :     {
    2264             :       GEN a, b;
    2265        8134 :       a = Q_content(gel(x,1)); if (!a) return NULL;
    2266        8134 :       b = Q_content(gel(x,2)); if (!b) return NULL;
    2267        8134 :       return Q_denom(gdiv(a, b));
    2268             :     }
    2269             :   }
    2270          66 :   return NULL;
    2271             : }
    2272             : GEN
    2273     3187988 : Q_denom(GEN x)
    2274             : {
    2275     3187988 :   GEN d = Q_denom_safe(x);
    2276     3187985 :   if (!d) pari_err_TYPE("Q_denom",x);
    2277     3187985 :   return d;
    2278             : }
    2279             : 
    2280             : GEN
    2281    55608970 : Q_remove_denom(GEN x, GEN *ptd)
    2282             : {
    2283    55608970 :   GEN d = Q_denom_safe(x);
    2284    55609013 :   if (d) { if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d); }
    2285    55608775 :   if (ptd) *ptd = d;
    2286    55608775 :   return x;
    2287             : }
    2288             : 
    2289             : /* return y = x * d, assuming x rational, and d,y integral */
    2290             : GEN
    2291   141444982 : Q_muli_to_int(GEN x, GEN d)
    2292             : {
    2293             :   GEN y, xn, xd;
    2294             :   pari_sp av;
    2295             : 
    2296   141444982 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    2297   141449517 :   switch (typ(x))
    2298             :   {
    2299    45216915 :     case t_INT:
    2300    45216915 :       return mulii(x,d);
    2301             : 
    2302    65476839 :     case t_FRAC:
    2303    65476839 :       xn = gel(x,1);
    2304    65476839 :       xd = gel(x,2); av = avma;
    2305    65476839 :       y = mulii(xn, diviiexact(d, xd));
    2306    65472181 :       return gc_INT(av, y);
    2307          42 :     case t_COMPLEX:
    2308          42 :       y = cgetg(3,t_COMPLEX);
    2309          42 :       gel(y,1) = Q_muli_to_int(gel(x,1),d);
    2310          42 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2311          42 :       return y;
    2312          14 :     case t_PADIC:
    2313          14 :       y = gcopy(x); if (!isint1(d)) setvalp(y, 0);
    2314          14 :       return y;
    2315         175 :     case t_QUAD:
    2316         175 :       y = cgetg(4,t_QUAD);
    2317         175 :       gel(y,1) = ZX_copy(gel(x,1));
    2318         175 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2319         175 :       gel(y,3) = Q_muli_to_int(gel(x,3),d); return y;
    2320             : 
    2321    20016380 :     case t_VEC:
    2322             :     case t_COL:
    2323   101723354 :     case t_MAT: pari_APPLY_same(Q_muli_to_int(gel(x,i), d));
    2324    48064284 :     case t_POL: pari_APPLY_pol_normalized(Q_muli_to_int(gel(x,i), d));
    2325          21 :     case t_SER: pari_APPLY_ser_normalized(Q_muli_to_int(gel(x,i), d));
    2326             : 
    2327       50539 :     case t_POLMOD:
    2328       50539 :       retmkpolmod(Q_muli_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2329          21 :     case t_RFRAC:
    2330          21 :       return gmul(x, d);
    2331             :   }
    2332           0 :   pari_err_TYPE("Q_muli_to_int",x);
    2333             :   return NULL; /* LCOV_EXCL_LINE */
    2334             : }
    2335             : 
    2336             : static void
    2337    29583940 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
    2338             : {
    2339             :   long e, i;
    2340    29583940 :   switch(typ(c))
    2341             :   {
    2342    20064753 :     case t_REAL:
    2343    20064753 :       *exact = 0;
    2344    20064753 :       if (!signe(c)) return;
    2345    19541445 :       e = expo(c) + 1 - bit_prec(c);
    2346    22057566 :       for (i = lg(c)-1; i > 2; i--, e += BITS_IN_LONG)
    2347    16674374 :         if (c[i]) break;
    2348    19541443 :       e += vals(c[i]); break; /* e[2] != 0 */
    2349     9514608 :     case t_INT:
    2350     9514608 :       if (!signe(c)) return;
    2351     1342276 :       e = expi(c);
    2352     1342287 :       break;
    2353        4545 :     case t_FRAC:
    2354        4545 :       e = expi(gel(c,1)) - expi(gel(c,2));
    2355        4545 :       if (*exact) *D = lcmii(*D, gel(c,2));
    2356        4545 :       break;
    2357          48 :     default:
    2358          48 :       pari_err_TYPE("rescale_to_int",c);
    2359             :       return; /* LCOV_EXCL_LINE */
    2360             :   }
    2361    20888274 :   if (e < *emin) *emin = e;
    2362             : }
    2363             : GEN
    2364     4607937 : RgM_rescale_to_int(GEN x)
    2365             : {
    2366     4607937 :   long lx = lg(x), i,j, hx, emin;
    2367             :   GEN D;
    2368             :   int exact;
    2369             : 
    2370     4607937 :   if (lx == 1) return cgetg(1,t_MAT);
    2371     4607937 :   hx = lgcols(x);
    2372     4607936 :   exact = 1;
    2373     4607936 :   emin = HIGHEXPOBIT;
    2374     4607936 :   D = gen_1;
    2375    15313109 :   for (j = 1; j < lx; j++)
    2376    40097003 :     for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
    2377     4607880 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2378     4607781 :   return grndtoi(gmul2n(x, -emin), NULL);
    2379             : }
    2380             : GEN
    2381       37478 : RgX_rescale_to_int(GEN x)
    2382             : {
    2383       37478 :   long lx = lg(x), i, emin;
    2384             :   GEN D;
    2385             :   int exact;
    2386       37478 :   if (lx == 2) return gcopy(x); /* rare */
    2387       37478 :   exact = 1;
    2388       37478 :   emin = HIGHEXPOBIT;
    2389       37478 :   D = gen_1;
    2390      229590 :   for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
    2391       37478 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2392       36351 :   return grndtoi(gmul2n(x, -emin), NULL);
    2393             : }
    2394             : 
    2395             : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    2396             : static GEN
    2397    11576772 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    2398             : {
    2399             :   GEN y, xn, xd;
    2400             :   pari_sp av;
    2401             : 
    2402    11576772 :   switch(typ(x))
    2403             :   {
    2404     2904071 :     case t_INT:
    2405     2904071 :       av = avma; y = diviiexact(x,d);
    2406     2904071 :       return gc_INT(av, mulii(y,n));
    2407             : 
    2408     5822988 :     case t_FRAC:
    2409     5822988 :       xn = gel(x,1);
    2410     5822988 :       xd = gel(x,2); av = avma;
    2411     5822988 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    2412     5822988 :       return gc_INT(av, y);
    2413             : 
    2414      451791 :     case t_VEC:
    2415             :     case t_COL:
    2416     4052134 :     case t_MAT: pari_APPLY_same(Q_divmuli_to_int(gel(x,i), d,n));
    2417     7892066 :     case t_POL: pari_APPLY_pol_normalized(Q_divmuli_to_int(gel(x,i), d,n));
    2418             : 
    2419           0 :     case t_RFRAC:
    2420           0 :       av = avma;
    2421           0 :       return gc_upto(av, gmul(x,mkfrac(n,d)));
    2422             : 
    2423           0 :     case t_POLMOD:
    2424           0 :       retmkpolmod(Q_divmuli_to_int(gel(x,2), d,n), RgX_copy(gel(x,1)));
    2425             :   }
    2426           0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    2427             :   return NULL; /* LCOV_EXCL_LINE */
    2428             : }
    2429             : 
    2430             : /* return x / d. x: rational; d,result: integral. */
    2431             : static GEN
    2432   172175421 : Q_divi_to_int(GEN x, GEN d)
    2433             : {
    2434   172175421 :   switch(typ(x))
    2435             :   {
    2436   145159289 :     case t_INT:
    2437   145159289 :       return diviiexact(x,d);
    2438             : 
    2439    20308432 :     case t_VEC:
    2440             :     case t_COL:
    2441   159494579 :     case t_MAT: pari_APPLY_same(Q_divi_to_int(gel(x,i), d));
    2442    25503922 :     case t_POL: pari_APPLY_pol_normalized(Q_divi_to_int(gel(x,i), d));
    2443             : 
    2444           0 :     case t_RFRAC:
    2445           0 :       return gdiv(x,d);
    2446             : 
    2447        5887 :     case t_POLMOD:
    2448        5887 :       retmkpolmod(Q_divi_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2449             :   }
    2450           0 :   pari_err_TYPE("Q_divi_to_int",x);
    2451             :   return NULL; /* LCOV_EXCL_LINE */
    2452             : }
    2453             : /* c t_FRAC */
    2454             : static GEN
    2455    10316881 : Q_divq_to_int(GEN x, GEN c)
    2456             : {
    2457    10316881 :   GEN n = gel(c,1), d = gel(c,2);
    2458    10316881 :   if (is_pm1(n)) {
    2459     7834595 :     GEN y = Q_muli_to_int(x,d);
    2460     7834519 :     if (signe(n) < 0) y = gneg(y);
    2461     7834519 :     return y;
    2462             :   }
    2463     2482285 :   return Q_divmuli_to_int(x, n,d);
    2464             : }
    2465             : 
    2466             : /* return y = x / c, assuming x,c rational, and y integral */
    2467             : GEN
    2468      486321 : Q_div_to_int(GEN x, GEN c)
    2469             : {
    2470      486321 :   switch(typ(c))
    2471             :   {
    2472      485386 :     case t_INT:  return Q_divi_to_int(x, c);
    2473         935 :     case t_FRAC: return Q_divq_to_int(x, c);
    2474             :   }
    2475           0 :   pari_err_TYPE("Q_div_to_int",c);
    2476             :   return NULL; /* LCOV_EXCL_LINE */
    2477             : }
    2478             : /* return y = x * c, assuming x,c rational, and y integral */
    2479             : GEN
    2480           0 : Q_mul_to_int(GEN x, GEN c)
    2481             : {
    2482             :   GEN d, n;
    2483           0 :   switch(typ(c))
    2484             :   {
    2485           0 :     case t_INT: return Q_muli_to_int(x, c);
    2486           0 :     case t_FRAC:
    2487           0 :       n = gel(c,1);
    2488           0 :       d = gel(c,2);
    2489           0 :       return Q_divmuli_to_int(x, d,n);
    2490             :   }
    2491           0 :   pari_err_TYPE("Q_mul_to_int",c);
    2492             :   return NULL; /* LCOV_EXCL_LINE */
    2493             : }
    2494             : 
    2495             : GEN
    2496    88425505 : Q_primitive_part(GEN x, GEN *ptc)
    2497             : {
    2498    88425505 :   pari_sp av = avma;
    2499    88425505 :   GEN c = Q_content_safe(x);
    2500    88423284 :   if (c)
    2501             :   {
    2502    88423389 :     if (typ(c) == t_INT)
    2503             :     {
    2504    78107442 :       if (equali1(c)) { set_avma(av); c = NULL; }
    2505    14143839 :       else if (signe(c)) x = Q_divi_to_int(x, c);
    2506             :     }
    2507    10315947 :     else x = Q_divq_to_int(x, c);
    2508             :   }
    2509    88420613 :   if (ptc) *ptc = c;
    2510    88420613 :   return x;
    2511             : }
    2512             : GEN
    2513     9492934 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    2514             : GEN
    2515      116491 : vec_Q_primpart(GEN x)
    2516      661784 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
    2517             : GEN
    2518       17738 : row_Q_primpart(GEN M)
    2519       17738 : { return shallowtrans(vec_Q_primpart(shallowtrans(M))); }
    2520             : 
    2521             : /*******************************************************************/
    2522             : /*                                                                 */
    2523             : /*                           SUBRESULTANT                          */
    2524             : /*                                                                 */
    2525             : /*******************************************************************/
    2526             : /* for internal use */
    2527             : GEN
    2528    27757821 : gdivexact(GEN x, GEN y)
    2529             : {
    2530             :   long i,lx;
    2531             :   GEN z;
    2532    27757821 :   if (gequal1(y)) return x;
    2533    27746167 :   if (typ(y) == t_POLMOD) return gmul(x, ginv(y));
    2534    27746069 :   switch(typ(x))
    2535             :   {
    2536    22929854 :     case t_INT:
    2537    22929854 :       if (typ(y)==t_INT) return diviiexact(x,y);
    2538          30 :       if (!signe(x)) return gen_0;
    2539           0 :       break;
    2540        8421 :     case t_INTMOD:
    2541             :     case t_FFELT:
    2542        8421 :     case t_POLMOD: return gmul(x,ginv(y));
    2543     4820819 :     case t_POL:
    2544     4820819 :       switch(typ(y))
    2545             :       {
    2546         714 :         case t_INTMOD:
    2547             :         case t_FFELT:
    2548         714 :         case t_POLMOD: return gmul(x,ginv(y));
    2549      222944 :         case t_POL: { /* not stack-clean */
    2550             :           long v;
    2551      222944 :           if (varn(x)!=varn(y)) break;
    2552      221978 :           v = RgX_valrem(y,&y);
    2553      221978 :           if (v) x = RgX_shift_shallow(x,-v);
    2554      221978 :           if (!degpol(y)) { y = gel(y,2); break; }
    2555      217740 :           return RgX_div(x,y);
    2556             :         }
    2557           0 :         case t_RFRAC:
    2558           0 :           if (varn(gel(y,2)) != varn(x)) break;
    2559           0 :           return gdiv(x, y);
    2560             :       }
    2561     4602365 :       return RgX_Rg_divexact(x, y);
    2562        4946 :     case t_VEC: case t_COL: case t_MAT:
    2563        4946 :       lx = lg(x); z = new_chunk(lx);
    2564       54182 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    2565        4946 :       z[0] = x[0]; return z;
    2566             :   }
    2567           0 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    2568           0 :   return gdiv(x,y);
    2569             : }
    2570             : 
    2571             : static GEN
    2572     1151598 : init_resultant(GEN x, GEN y)
    2573             : {
    2574     1151598 :   long tx = typ(x), ty = typ(y), vx, vy;
    2575     1151598 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    2576             :   {
    2577          14 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    2578          14 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    2579           0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    2580           0 :     return gen_1;
    2581             :   }
    2582     1151583 :   if (tx!=t_POL) pari_err_TYPE("resultant",x);
    2583     1151583 :   if (ty!=t_POL) pari_err_TYPE("resultant",y);
    2584     1151583 :   if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
    2585     1151506 :   vx = varn(x);
    2586     1151506 :   vy = varn(y); if (vx == vy) return NULL;
    2587           7 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    2588             : }
    2589             : 
    2590             : /* x an RgX, y a scalar */
    2591             : static GEN
    2592           7 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    2593             : {
    2594           7 :   *V = gpowgs(y,degpol(x)-1);
    2595           7 :   *U = gen_0; return gmul(y, *V);
    2596             : }
    2597             : 
    2598             : /* return 0 if the subresultant chain can be interrupted.
    2599             :  * Set u = NULL if the resultant is 0. */
    2600             : static int
    2601       11790 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    2602             : {
    2603       11790 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    2604             :   long du, dv, dr, degq;
    2605             : 
    2606       11790 :   if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
    2607       11790 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    2608       11545 :   du = degpol(*u);
    2609       11545 :   dv = degpol(*v);
    2610       11545 :   degq = du - dv;
    2611       11545 :   if (*um1 == gen_1)
    2612        6413 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    2613        5132 :   else if (*um1 == gen_0)
    2614        2318 :     u0 = gen_0;
    2615             :   else /* except in those 2 cases, um1 is an RgX */
    2616        2814 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2617             : 
    2618       11545 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2619        6413 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2620             :   else
    2621        5132 :     u0 = gsub(u0, gmul(q,*uze));
    2622             : 
    2623       11545 :   *um1 = *uze;
    2624       11545 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2625             : 
    2626       11545 :   *u = *v; c = *g; *g  = leading_coeff(*u);
    2627       11545 :   switch(degq)
    2628             :   {
    2629        1666 :     case 0: break;
    2630        8059 :     case 1:
    2631        8059 :       c = gmul(*h,c); *h = *g; break;
    2632        1820 :     default:
    2633        1820 :       c = gmul(gpowgs(*h,degq), c);
    2634        1820 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2635             :   }
    2636       11545 :   if (typ(c) == t_POLMOD)
    2637             :   {
    2638         904 :     c = ginv(c);
    2639         904 :     *v = RgX_Rg_mul(r,c);
    2640         904 :     *uze = RgX_Rg_mul(*uze,c);
    2641             :   }
    2642             :   else
    2643             :   {
    2644       10641 :     *v  = RgX_Rg_divexact(r,c);
    2645       10641 :     *uze= RgX_Rg_divexact(*uze,c);
    2646             :   }
    2647       11545 :   if (both_odd(du, dv)) *signh = -*signh;
    2648       11545 :   return (dr > 3);
    2649             : }
    2650             : 
    2651             : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2652             : static GEN
    2653        2346 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2654             : {
    2655             :   pari_sp av, av2;
    2656        2346 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2657             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2658             : 
    2659        2346 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2660        2346 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2661        2346 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2662        2346 :   if (tx != t_POL) {
    2663           7 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2664           7 :     return scalar_res(y,x,V,U);
    2665             :   }
    2666        2339 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2667        2339 :   if (varn(x) != varn(y))
    2668           0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2669           0 :                                         : scalar_res(y,x,V,U);
    2670        2339 :   if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
    2671        2339 :   if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
    2672        2339 :   dx = degpol(x);
    2673        2339 :   dy = degpol(y);
    2674        2339 :   signh = 1;
    2675        2339 :   if (dx < dy)
    2676             :   {
    2677         848 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2678         848 :     if (both_odd(dx, dy)) signh = -signh;
    2679             :   }
    2680        2339 :   if (dy == 0)
    2681             :   {
    2682           0 :     *V = gpowgs(gel(y,2),dx-1);
    2683           0 :     *U = gen_0; return gmul(*V,gel(y,2));
    2684             :   }
    2685        2339 :   av = avma;
    2686        2339 :   u = x = primitive_part(x, &cu);
    2687        2339 :   v = y = primitive_part(y, &cv);
    2688        2339 :   g = h = gen_1; av2 = avma;
    2689        2339 :   um1 = gen_1; uze = gen_0;
    2690             :   for(;;)
    2691             :   {
    2692        6995 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2693        4656 :     if (gc_needed(av2,1))
    2694             :     {
    2695           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2696           0 :       (void)gc_all(av2,6, &u,&v,&g,&h,&uze,&um1);
    2697             :     }
    2698             :   }
    2699             :   /* uze an RgX */
    2700        2339 :   if (!u) { *U = *V = gen_0; return gc_const(av, gen_0); }
    2701        2332 :   z = gel(v,2); du = degpol(u);
    2702        2332 :   if (du > 1)
    2703             :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2704         252 :     p1 = gpowgs(gdiv(z,h),du-1);
    2705         252 :     z = gmul(z,p1);
    2706         252 :     uze = RgX_Rg_mul(uze, p1);
    2707             :   }
    2708        2332 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2709             : 
    2710        2332 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2711        2332 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2712             :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2713        2332 :   p1 = gen_1;
    2714        2332 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2715        2332 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2716        2332 :   cu = cu? gdiv(p1,cu): p1;
    2717        2332 :   cv = cv? gdiv(p1,cv): p1;
    2718        2332 :   z = gmul(z,p1);
    2719        2332 :   *U = RgX_Rg_mul(uze,cu);
    2720        2332 :   *V = RgX_Rg_mul(vze,cv);
    2721        2332 :   return z;
    2722             : }
    2723             : GEN
    2724           0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2725             : {
    2726           0 :   pari_sp av = avma;
    2727           0 :   GEN z = subresext_i(x, y, U, V);
    2728           0 :   return gc_all(av, 3, &z, U, V);
    2729             : }
    2730             : 
    2731             : static GEN
    2732         434 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2733             : {
    2734         434 :   GEN x=content(y);
    2735         434 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2736             : }
    2737             : 
    2738             : static int
    2739        4368 : must_negate(GEN x)
    2740             : {
    2741        4368 :   GEN t = leading_coeff(x);
    2742        4368 :   switch(typ(t))
    2743             :   {
    2744        4263 :     case t_INT: case t_REAL:
    2745        4263 :       return (signe(t) < 0);
    2746           0 :     case t_FRAC:
    2747           0 :       return (signe(gel(t,1)) < 0);
    2748             :   }
    2749         105 :   return 0;
    2750             : }
    2751             : 
    2752             : static GEN
    2753         217 : gc_gcdext(pari_sp av, GEN r, GEN *u, GEN *v)
    2754             : {
    2755         217 :   if (!u && !v) return gc_upto(av, r);
    2756         217 :   if (u  &&  v) return gc_all(av, 3, &r, u, v);
    2757           0 :   return gc_all(av, 2, &r, u ? u: v);
    2758             : }
    2759             : 
    2760             : static GEN
    2761         133 : RgX_extgcd_FpX(GEN x, GEN y, GEN p, GEN *u, GEN *v)
    2762             : {
    2763         133 :   pari_sp av = avma;
    2764         133 :   GEN r = FpX_extgcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, u, v);
    2765         133 :   if (u) *u = FpX_to_mod(*u, p);
    2766         133 :   if (v) *v = FpX_to_mod(*v, p);
    2767         133 :   return gc_gcdext(av, FpX_to_mod(r, p), u, v);
    2768             : }
    2769             : 
    2770             : static GEN
    2771           7 : RgX_extgcd_FpXQX(GEN x, GEN y, GEN pol, GEN p, GEN *U, GEN *V)
    2772             : {
    2773           7 :   pari_sp av = avma;
    2774           7 :   GEN r, T = RgX_to_FpX(pol, p);
    2775           7 :   r = FpXQX_extgcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p, U, V);
    2776           7 :   return gc_gcdext(av, FpXQX_to_mod(r, T, p), U, V);
    2777             : }
    2778             : 
    2779             : static GEN
    2780        4529 : RgX_extgcd_fast(GEN x, GEN y, GEN *U, GEN *V)
    2781             : {
    2782             :   GEN p, pol;
    2783             :   long pa;
    2784        4529 :   long t = RgX_type(x, &p,&pol,&pa);
    2785        4529 :   switch(t)
    2786             :   {
    2787          21 :     case t_FFELT:  return FFX_extgcd(x, y, pol, U, V);
    2788         133 :     case t_INTMOD: return RgX_extgcd_FpX(x, y, p, U, V);
    2789           7 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    2790           7 :                    return RgX_extgcd_FpXQX(x, y, pol, p, U, V);
    2791        4368 :     default:       return NULL;
    2792             :   }
    2793             : }
    2794             : 
    2795             : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2796             : GEN
    2797        4970 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2798             : {
    2799             :   pari_sp av, av2, tetpil;
    2800             :   long signh; /* junk */
    2801        4970 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2802             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2803             : 
    2804        4970 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2805        4970 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2806        4970 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2807        4970 :   vx=varn(x);
    2808        4970 :   if (!signe(x))
    2809             :   {
    2810          14 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2811           7 :     *U = pol_0(vx); *V = pol_0(vx);
    2812           7 :     return pol_0(vx);
    2813             :   }
    2814        4956 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2815        4529 :   r = RgX_extgcd_fast(x, y, U, V);
    2816        4529 :   if (r) return r;
    2817        4368 :   dx = degpol(x); dy = degpol(y);
    2818        4368 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2819        4368 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2820             : 
    2821        4179 :   av = avma;
    2822        4179 :   u = x = primitive_part(x, &cu);
    2823        4179 :   v = y = primitive_part(y, &cv);
    2824        4179 :   g = h = gen_1; av2 = avma;
    2825        4179 :   um1 = gen_1; uze = gen_0;
    2826             :   for(;;)
    2827             :   {
    2828        4389 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2829         210 :     if (gc_needed(av2,1))
    2830             :     {
    2831           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2832           0 :       (void)gc_all(av2,6,&u,&v,&g,&h,&uze,&um1);
    2833             :     }
    2834             :   }
    2835        4179 :   if (uze != gen_0) {
    2836             :     GEN r;
    2837        3969 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2838        3969 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2839        3969 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2840        3969 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2841        3969 :     p1 = ginv(content(v));
    2842             :   }
    2843             :   else /* y | x */
    2844             :   {
    2845         210 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2846         210 :     uze = pol_0(vx);
    2847         210 :     p1 = gen_1;
    2848             :   }
    2849        4179 :   if (must_negate(v)) p1 = gneg(p1);
    2850        4179 :   tetpil = avma;
    2851        4179 :   z = RgX_Rg_mul(v,p1);
    2852        4179 :   *U = RgX_Rg_mul(uze,p1);
    2853        4179 :   *V = RgX_Rg_mul(vze,p1);
    2854        4179 :   return gc_all_unsafe(av,tetpil, 3, &z, U, V);
    2855             : }
    2856             : 
    2857             : static GEN
    2858          14 : RgX_halfgcd_all_i(GEN a, GEN b, GEN *pa, GEN *pb)
    2859             : {
    2860          14 :   pari_sp av=avma;
    2861          14 :   long m = degpol(a), va = varn(a);
    2862             :   GEN R, u,u1,v,v1;
    2863          14 :   u1 = v = pol_0(va);
    2864          14 :   u = v1 = pol_1(va);
    2865          14 :   if (degpol(a)<degpol(b))
    2866             :   {
    2867           0 :     swap(a,b);
    2868           0 :     swap(u,v); swap(u1,v1);
    2869             :   }
    2870          42 :   while (2*degpol(b) >= m)
    2871             :   {
    2872          28 :     GEN r, q = RgX_pseudodivrem(a,b,&r);
    2873          28 :     GEN l = gpowgs(leading_coeff(b), degpol(a)-degpol(b)+1);
    2874          28 :     GEN g = ggcd(l, content(r));
    2875          28 :     q = RgX_Rg_div(q, g);
    2876          28 :     r = RgX_Rg_div(r, g);
    2877          28 :     l = gdiv(l, g);
    2878          28 :     a = b; b = r; swap(u,v); swap(u1,v1);
    2879          28 :     v  = RgX_sub(gmul(l,v),  RgX_mul(u, q));
    2880          28 :     v1 = RgX_sub(gmul(l,v1), RgX_mul(u1, q));
    2881          28 :     if (gc_needed(av,2))
    2882             :     {
    2883           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"halfgcd (d = %ld)",degpol(b));
    2884           0 :       (void)gc_all(av,6, &a,&b,&u1,&v1,&u,&v);
    2885             :     }
    2886             :   }
    2887          14 :   if (pa) *pa = a;
    2888          14 :   if (pb) *pb = b;
    2889          14 :   R = mkmat22(u,u1,v,v1);
    2890          14 :   return !pa && pb ? gc_all(av, 2, &R, pb): gc_all(av, 1+!!pa+!!pb, &R, pa, pb);
    2891             : }
    2892             : 
    2893             : static GEN
    2894          28 : RgX_halfgcd_all_FpX(GEN x, GEN y, GEN p, GEN *a, GEN *b)
    2895             : {
    2896          28 :   pari_sp av = avma;
    2897             :   GEN M;
    2898          28 :   if (lgefint(p) == 3)
    2899             :   {
    2900          14 :     ulong pp = uel(p, 2);
    2901          14 :     GEN xp = RgX_to_Flx(x, pp), yp = RgX_to_Flx(y, pp);
    2902          14 :     M = Flx_halfgcd_all(xp, yp, pp, a, b);
    2903          14 :     M = FlxM_to_ZXM(M); *a = Flx_to_ZX(*a); *b = Flx_to_ZX(*b);
    2904             :   }
    2905             :   else
    2906             :   {
    2907          14 :     x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
    2908          14 :     M = FpX_halfgcd_all(x, y, p, a, b);
    2909             :   }
    2910          28 :   return !a && b ? gc_all(av, 2, &M, b): gc_all(av, 1+!!a+!!b, &M, a, b);
    2911             : }
    2912             : 
    2913             : static GEN
    2914           0 : RgX_halfgcd_all_FpXQX(GEN x, GEN y, GEN pol, GEN p, GEN *a, GEN *b)
    2915             : {
    2916           0 :   pari_sp av = avma;
    2917           0 :   GEN M, T = RgX_to_FpX(pol, p);
    2918           0 :   if (signe(T)==0) pari_err_OP("halfgcd", x, y);
    2919           0 :   x = RgX_to_FpXQX(x, T, p); y = RgX_to_FpXQX(y, T, p);
    2920           0 :   M = FpXQX_halfgcd_all(x, y, T, p, a, b);
    2921           0 :   if (a) *a = FqX_to_mod(*a, T, p);
    2922           0 :   if (b) *b = FqX_to_mod(*b, T, p);
    2923           0 :   M = FqXM_to_mod(M, T, p);
    2924           0 :   return !a && b ? gc_all(av, 2, &M, b): gc_all(av, 1+!!a+!!b, &M, a, b);
    2925             : }
    2926             : 
    2927             : static GEN
    2928          63 : RgX_halfgcd_all_fast(GEN x, GEN y, GEN *a, GEN *b)
    2929             : {
    2930             :   GEN p, pol;
    2931             :   long pa;
    2932          63 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    2933          63 :   switch(t)
    2934             :   {
    2935          21 :     case t_FFELT:  return FFX_halfgcd_all(x, y, pol, a, b);
    2936          28 :     case t_INTMOD: return RgX_halfgcd_all_FpX(x, y, p, a, b);
    2937           0 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    2938           0 :                    return RgX_halfgcd_all_FpXQX(x, y, pol, p, a, b);
    2939          14 :     default:       return NULL;
    2940             :   }
    2941             : }
    2942             : 
    2943             : GEN
    2944          63 : RgX_halfgcd_all(GEN x, GEN y, GEN *a, GEN *b)
    2945             : {
    2946          63 :   GEN z = RgX_halfgcd_all_fast(x, y, a, b);
    2947          63 :   if (z) return z;
    2948          14 :   return RgX_halfgcd_all_i(x, y, a, b);
    2949             : }
    2950             : 
    2951             : GEN
    2952           0 : RgX_halfgcd(GEN x, GEN y)
    2953           0 : { return  RgX_halfgcd_all(x, y, NULL, NULL); }
    2954             : 
    2955             : int
    2956         112 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2957             : {
    2958         112 :   pari_sp av = avma, av2, tetpil;
    2959             :   long signh; /* junk */
    2960             :   long vx;
    2961             :   GEN g, h, p1, cu, cv, u, v, um1, uze;
    2962             : 
    2963         112 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2964         112 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2965         112 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2966         112 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2967         112 :   if (!signe(T)) {
    2968           0 :     if (degpol(x) <= amax) {
    2969           0 :       *P = RgX_copy(x);
    2970           0 :       *Q = pol_1(varn(x));
    2971           0 :       return 1;
    2972             :     }
    2973           0 :     return 0;
    2974             :   }
    2975         112 :   if (amax+bmax >= degpol(T))
    2976           0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2977             :                     mkvec3(stoi(amax), stoi(bmax), T));
    2978         112 :   vx = varn(T);
    2979         112 :   u = x = primitive_part(x, &cu);
    2980         112 :   v = T = primitive_part(T, &cv);
    2981         112 :   g = h = gen_1; av2 = avma;
    2982         112 :   um1 = gen_1; uze = gen_0;
    2983             :   for(;;)
    2984             :   {
    2985         406 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2986         406 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) return gc_bool(av,0);
    2987         406 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2988         294 :     if (gc_needed(av2,1))
    2989             :     {
    2990           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2991           0 :       (void)gc_all(av2,6,&u,&v,&g,&h,&uze,&um1);
    2992             :     }
    2993             :   }
    2994         112 :   if (uze == gen_0)
    2995             :   {
    2996           0 :     set_avma(av); *P = pol_0(vx); *Q = pol_1(vx);
    2997           0 :     return 1;
    2998             :   }
    2999         112 :   if (cu) uze = RgX_Rg_div(uze,cu);
    3000         112 :   p1 = ginv(content(v));
    3001         112 :   if (must_negate(v)) p1 = gneg(p1);
    3002         112 :   tetpil = avma;
    3003         112 :   *P = RgX_Rg_mul(v,p1);
    3004         112 :   *Q = RgX_Rg_mul(uze,p1);
    3005         112 :   (void)gc_all_unsafe(av,tetpil,2,P,Q); return 1;
    3006             : }
    3007             : 
    3008             : GEN
    3009           0 : RgX_chinese_coprime(GEN x, GEN y, GEN Tx, GEN Ty, GEN Tz)
    3010             : {
    3011           0 :   pari_sp av = avma;
    3012           0 :   GEN ax = RgX_mul(RgXQ_inv(Tx,Ty), Tx);
    3013           0 :   GEN p1 = RgX_mul(ax, RgX_sub(y,x));
    3014           0 :   p1 = RgX_add(x,p1);
    3015           0 :   if (!Tz) Tz = RgX_mul(Tx,Ty);
    3016           0 :   p1 = RgX_rem(p1, Tz);
    3017           0 :   return gc_upto(av,p1);
    3018             : }
    3019             : 
    3020             : /*******************************************************************/
    3021             : /*                                                                 */
    3022             : /*                RESULTANT USING DUCOS VARIANT                    */
    3023             : /*                                                                 */
    3024             : /*******************************************************************/
    3025             : /* x^n / y^(n-1), assume n > 0 */
    3026             : static GEN
    3027      629317 : Lazard(GEN x, GEN y, long n)
    3028             : {
    3029             :   long a;
    3030             :   GEN c;
    3031             : 
    3032      629317 :   if (n == 1) return x;
    3033       10020 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    3034       10020 :   c=x; n-=a;
    3035       20663 :   while (a>1)
    3036             :   {
    3037       10643 :     a>>=1; c=gdivexact(gsqr(c),y);
    3038       10643 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    3039             :   }
    3040       10020 :   return c;
    3041             : }
    3042             : 
    3043             : /* F (x/y)^(n-1), assume n >= 1 */
    3044             : static GEN
    3045      834543 : Lazard2(GEN F, GEN x, GEN y, long n)
    3046             : {
    3047      834543 :   if (n == 1) return F;
    3048        4822 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    3049             : }
    3050             : 
    3051             : static GEN
    3052      834542 : RgX_neg_i(GEN x, long lx)
    3053             : {
    3054             :   long i;
    3055      834542 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    3056     2184609 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    3057      834543 :   return y;
    3058             : }
    3059             : static GEN
    3060     2451568 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    3061             : {
    3062             :   long i;
    3063             :   GEN z;
    3064     2451568 :   if (isrationalzero(x)) return pol_0(varn(y));
    3065     2451293 :   z = cgetg(ly,t_POL); z[1] = y[1];
    3066     6434633 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    3067     2451267 :   return z;
    3068             : }
    3069             : static long
    3070     2446461 : reductum_lg(GEN x, long lx)
    3071             : {
    3072     2446461 :   long i = lx-2;
    3073     2518491 :   while (i > 1 && gequal0(gel(x,i))) i--;
    3074     2446464 :   return i+1;
    3075             : }
    3076             : 
    3077             : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
    3078             : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    3079             :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    3080             : static GEN
    3081      834544 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    3082             : {
    3083      834544 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
    3084             :   long p, q, j, lP, lQ;
    3085             :   pari_sp av;
    3086             : 
    3087      834544 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    3088      834543 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    3089             :   /* p > q. Very often p - 1 = q */
    3090      834542 :   av = avma;
    3091             :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    3092      834542 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    3093             : 
    3094      834543 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    3095      843094 :   for (j = q+1; j < p; j++)
    3096             :   {
    3097        8558 :     if (degpol(H) == q-1)
    3098             :     { /* h0 = coeff of degree q-1 = leading coeff */
    3099        7354 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    3100        7354 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    3101             :     }
    3102             :     else
    3103        1204 :       H = RgX_shift_shallow(H, 1);
    3104        8558 :     if (j+2 < lP)
    3105             :     {
    3106        6515 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    3107        6515 :       A = A? RgX_add(A, TMP): TMP;
    3108             :     }
    3109        8558 :     if (gc_needed(av,1))
    3110             :     {
    3111         147 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    3112         147 :       (void)gc_all(av,A?2:1,&H,&A);
    3113             :     }
    3114             :   }
    3115      834536 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    3116      834535 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    3117      834542 :   A = A? RgX_add(A, TMP): TMP;
    3118      834543 :   A = RgX_Rg_divexact(A, p0);
    3119      834539 :   if (degpol(H) == q-1)
    3120             :   {
    3121      832293 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    3122      832294 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    3123             :   }
    3124             :   else
    3125        2245 :     A = RgX_Rg_mul(addshift(H,A), q0);
    3126      834542 :   return RgX_Rg_divexact(A, s);
    3127             : }
    3128             : #undef addshift
    3129             : 
    3130             : static GEN
    3131     1248991 : RgX_pseudodenom(GEN x)
    3132             : {
    3133     1248991 :   GEN m = NULL;
    3134     1248991 :   long l = lg(x), i;
    3135     6195352 :   for (i = 2; i < l; i++)
    3136             :   {
    3137     4946361 :     GEN xi = gel(x, i);
    3138     4946361 :     if (typ(xi) == t_RFRAC)
    3139             :     {
    3140          42 :       GEN d = denom_i(xi);
    3141          42 :       if (!m || signe(RgX_pseudorem(m, d)))
    3142          42 :         m = m ? gmul(m, d): d;
    3143             :     }
    3144             :   }
    3145     1248991 :   return m;
    3146             : }
    3147             : 
    3148             : /* Ducos's subresultant */
    3149             : GEN
    3150     1066704 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    3151             : {
    3152             :   pari_sp av, av2;
    3153     1066704 :   long dP, dQ, delta, sig = 1;
    3154             :   GEN DP, DQ, cP, cQ, Z, s;
    3155             : 
    3156     1066704 :   dP = degpol(P);
    3157     1066703 :   dQ = degpol(Q); delta = dP - dQ;
    3158     1066703 :   if (delta < 0)
    3159             :   {
    3160         910 :     if (both_odd(dP, dQ)) sig = -1;
    3161         910 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    3162             :   }
    3163     1066703 :   if (sol) *sol = gen_0;
    3164     1066703 :   av = avma;
    3165     1066703 :   if (dQ <= 0)
    3166             :   {
    3167         910 :     if (dQ < 0) return Rg_get_0(P);
    3168         910 :     s = gpowgs(gel(Q,2), dP);
    3169         910 :     if (sig == -1) s = gc_upto(av, gneg(s));
    3170         910 :     return s;
    3171             :   }
    3172     1065793 :   if (dQ == 1)
    3173             :   {
    3174      441297 :     if (sol) *sol = Q;
    3175      441297 :     s = RgX_homogenous_evalpow(P, gel(Q,2), gpowers(gneg(gel(Q,3)), dP));
    3176      441296 :     if (sig==-1) s = gneg(s);
    3177      441296 :     return gc_all(av, sol ? 2: 1, &s, sol);
    3178             :   }
    3179             :   /* primitive_part is also possible here, but possibly very costly,
    3180             :    * and hardly ever worth it */
    3181             : 
    3182      624496 :   DP = RgX_pseudodenom(P); if (DP) P = gmul(P,DP);
    3183      624496 :   DQ = RgX_pseudodenom(Q); if (DQ) Q = gmul(Q,DQ);
    3184      624496 :   P = Q_primitive_part(P, &cP);
    3185      624494 :   Q = Q_primitive_part(Q, &cQ);
    3186      624495 :   av2 = avma;
    3187      624495 :   s = gpowgs(leading_coeff(Q),delta);
    3188      624495 :   if (both_odd(dP, dQ)) sig = -sig;
    3189      624495 :   Z = Q;
    3190      624495 :   Q = RgX_pseudorem(P, Q);
    3191      624496 :   P = Z;
    3192     1459038 :   while(degpol(Q) > 0)
    3193             :   {
    3194      834539 :     delta = degpol(P) - degpol(Q); /* > 0 */
    3195      834543 :     Z = Lazard2(Q, leading_coeff(Q), s, delta);
    3196      834543 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    3197      834543 :     Q = nextSousResultant(P, Q, Z, s);
    3198      834542 :     P = Z;
    3199      834542 :     if (gc_needed(av,1))
    3200             :     {
    3201          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    3202          13 :       (void)gc_all(av2,2,&P,&Q);
    3203             :     }
    3204      834542 :     s = leading_coeff(P);
    3205             :   }
    3206      624494 :   if (!signe(Q)) { set_avma(av); return Rg_get_0(Q); }
    3207      624494 :   s = Lazard(leading_coeff(Q), s, degpol(P));
    3208      624495 :   if (sig == -1) s = gneg(s);
    3209      624495 :   if (DP) s = gdiv(s, gpowgs(DP,dQ));
    3210      624495 :   if (DQ) s = gdiv(s, gpowgs(DQ,dP));
    3211      624495 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    3212      624495 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    3213      624497 :   if (!sol) return gc_GEN(av, s);
    3214        2415 :   *sol = P; return gc_all(av, 2, &s, sol);
    3215             : }
    3216             : 
    3217             : static GEN
    3218          28 : RgX_resultant_FpX(GEN x, GEN y, GEN p)
    3219             : {
    3220          28 :   pari_sp av = avma;
    3221             :   GEN r;
    3222          28 :   if (lgefint(p) == 3)
    3223             :   {
    3224          14 :     ulong pp = uel(p, 2);
    3225          14 :     r = utoi(Flx_resultant(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), pp));
    3226             :   }
    3227             :   else
    3228          14 :     r = FpX_resultant(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3229          28 :   return gc_upto(av, Fp_to_mod(r, p));
    3230             : }
    3231             : 
    3232             : static GEN
    3233          21 : RgX_resultant_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3234             : {
    3235          21 :   pari_sp av = avma;
    3236          21 :   GEN r, T = RgX_to_FpX(pol, p);
    3237          21 :   r = FpXQX_resultant(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3238          21 :   return gc_upto(av, FpX_to_mod(r, p));
    3239             : }
    3240             : 
    3241             : static GEN
    3242     1151568 : resultant_fast(GEN x, GEN y)
    3243             : {
    3244             :   GEN p, pol;
    3245             :   long pa, t;
    3246     1151568 :   p = init_resultant(x,y);
    3247     1151570 :   if (p) return p;
    3248     1151472 :   t = RgX_type2(x,y, &p,&pol,&pa);
    3249     1151481 :   switch(t)
    3250             :   {
    3251        2688 :     case t_INT:    return ZX_resultant(x,y);
    3252          56 :     case t_FRAC:   return QX_resultant(x,y);
    3253          21 :     case t_FFELT:  return FFX_resultant(x,y,pol);
    3254          28 :     case t_INTMOD: return RgX_resultant_FpX(x, y, p);
    3255          21 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    3256          21 :                    return RgX_resultant_FpXQX(x, y, pol, p);
    3257     1148667 :     default:       return NULL;
    3258             :   }
    3259             : }
    3260             : 
    3261             : static GEN
    3262      169893 : RgX_resultant_sylvester(GEN x, GEN y)
    3263             : {
    3264      169893 :   pari_sp av = avma;
    3265      169893 :   return gc_upto(av, det(RgX_sylvestermatrix(x,y)));
    3266             : }
    3267             : 
    3268             : /* Return resultant(P,Q).
    3269             :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    3270             :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    3271             :  * in the "generic" case. */
    3272             : GEN
    3273     1151569 : resultant(GEN P, GEN Q)
    3274             : {
    3275     1151569 :   GEN z = resultant_fast(P,Q);
    3276     1151579 :   if (z) return z;
    3277     1148667 :   if (isinexact(P) || isinexact(Q)) return RgX_resultant_sylvester(P,Q);
    3278      978801 :   return RgX_resultant_all(P, Q, NULL);
    3279             : }
    3280             : 
    3281             : /*******************************************************************/
    3282             : /*                                                                 */
    3283             : /*               RESULTANT USING SYLVESTER MATRIX                  */
    3284             : /*                                                                 */
    3285             : /*******************************************************************/
    3286             : static GEN
    3287      371755 : syl_RgC(GEN x, long j, long d, long D, long cp)
    3288             : {
    3289      371755 :   GEN c = cgetg(d+1,t_COL);
    3290             :   long i;
    3291      990325 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    3292     2142732 :   for (   ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
    3293      990325 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    3294      371755 :   return c;
    3295             : }
    3296             : static GEN
    3297      169900 : syl_RgM(GEN x, GEN y, long cp)
    3298             : {
    3299      169900 :   long j, d, dx = degpol(x), dy = degpol(y);
    3300             :   GEN M;
    3301      169900 :   if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
    3302      169900 :   if (dy < 0) return zeromat(dx,dx);
    3303      169900 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    3304      442060 :   for (j=1; j<=dy; j++) gel(M,j)    = syl_RgC(x,j,d,j+dx, cp);
    3305      269495 :   for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
    3306      169900 :   return M;
    3307             : }
    3308             : GEN
    3309      169893 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
    3310             : GEN
    3311           7 : sylvestermatrix(GEN x, GEN y)
    3312             : {
    3313           7 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    3314           7 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    3315           7 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    3316           7 :   return syl_RgM(x,y,1);
    3317             : }
    3318             : 
    3319             : GEN
    3320          28 : resultant2(GEN x, GEN y)
    3321             : {
    3322          28 :   GEN r = init_resultant(x,y);
    3323          28 :   return r? r: RgX_resultant_sylvester(x,y);
    3324             : }
    3325             : 
    3326             : /* let vx = main variable of x, v0 a variable of highest priority;
    3327             :  * return a t_POL in variable v0:
    3328             :  * if vx <= v, return subst(x, v, pol_x(v0))
    3329             :  * if vx >  v, return scalarpol(x, v0) */
    3330             : static GEN
    3331         343 : fix_pol(GEN x, long v, long v0)
    3332             : {
    3333         343 :   long vx, tx = typ(x);
    3334         343 :   if (tx != t_POL)
    3335          42 :     vx = gvar(x);
    3336             :   else
    3337             :   { /* shortcut: almost nothing to do */
    3338         301 :     vx = varn(x);
    3339         301 :     if (v == vx)
    3340             :     {
    3341         119 :       if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
    3342         119 :       return x;
    3343             :     }
    3344             :   }
    3345         224 :   if (varncmp(v, vx) > 0)
    3346             :   {
    3347         217 :     x = gsubst(x, v, pol_x(v0));
    3348         217 :     if (typ(x) != t_POL) vx = gvar(x);
    3349             :     else
    3350             :     {
    3351         210 :       vx = varn(x);
    3352         210 :       if (vx == v0) return x;
    3353             :     }
    3354             :   }
    3355          49 :   if (varncmp(vx, v0) <= 0) pari_err_TYPE("polresultant", x);
    3356          42 :   return scalarpol_shallow(x, v0);
    3357             : }
    3358             : 
    3359             : /* resultant of x and y with respect to variable v, or with respect to their
    3360             :  * main variable if v < 0. */
    3361             : GEN
    3362         504 : polresultant0(GEN x, GEN y, long v, long flag)
    3363             : {
    3364         504 :   pari_sp av = avma;
    3365             : 
    3366         504 :   if (v >= 0)
    3367             :   {
    3368         147 :     long v0 = fetch_var_higher();
    3369         147 :     x = fix_pol(x,v, v0);
    3370         147 :     y = fix_pol(y,v, v0);
    3371             :   }
    3372         504 :   switch(flag)
    3373             :   {
    3374         497 :     case 2:
    3375         497 :     case 0: x=resultant(x,y); break;
    3376           7 :     case 1: x=resultant2(x,y); break;
    3377           0 :     default: pari_err_FLAG("polresultant");
    3378             :   }
    3379         504 :   if (v >= 0) (void)delete_var();
    3380         504 :   return gc_upto(av,x);
    3381             : }
    3382             : 
    3383             : static GEN
    3384          77 : RgX_extresultant_FpX(GEN x, GEN y, GEN p, GEN *u, GEN *v)
    3385             : {
    3386          77 :   pari_sp av = avma;
    3387          77 :   GEN r = FpX_extresultant(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, u, v);
    3388          77 :   if (signe(r) == 0) { *u = gen_0; *v = gen_0; return gc_const(av, gen_0); }
    3389          77 :   if (u) *u = FpX_to_mod(*u, p);
    3390          77 :   if (v) *v = FpX_to_mod(*v, p);
    3391          77 :   return gc_gcdext(av, Fp_to_mod(r, p), u, v);
    3392             : }
    3393             : 
    3394             : static GEN
    3395        1568 : RgX_extresultant_fast(GEN x, GEN y, GEN *U, GEN *V)
    3396             : {
    3397             :   GEN p, pol;
    3398             :   long pa;
    3399        1568 :   long t = RgX_type2(x, y, &p,&pol,&pa);
    3400        1568 :   switch(t)
    3401             :   {
    3402          77 :     case t_INTMOD: return RgX_extresultant_FpX(x, y, p, U, V);
    3403        1491 :     default:       return NULL;
    3404             :   }
    3405             : }
    3406             : 
    3407             : GEN
    3408        1575 : polresultantext0(GEN x, GEN y, long v)
    3409             : {
    3410        1575 :   GEN R = NULL, U, V;
    3411        1575 :   pari_sp av = avma;
    3412             : 
    3413        1575 :   if (v >= 0)
    3414             :   {
    3415          14 :     long v0 = fetch_var_higher();
    3416          14 :     x = fix_pol(x,v, v0);
    3417          14 :     y = fix_pol(y,v, v0);
    3418             :   }
    3419        1575 :   if (typ(x)==t_POL && typ(y)==t_POL)
    3420        1568 :     R = RgX_extresultant_fast(x, y, &U, &V);
    3421        1575 :   if (!R)
    3422        1498 :     R = subresext_i(x,y, &U,&V);
    3423        1575 :   if (v >= 0)
    3424             :   {
    3425          14 :     (void)delete_var();
    3426          14 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    3427          14 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    3428             :   }
    3429        1575 :   return gc_GEN(av, mkvec3(U,V,R));
    3430             : }
    3431             : GEN
    3432        1463 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    3433             : 
    3434             : /*******************************************************************/
    3435             : /*                                                                 */
    3436             : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    3437             : /*                                                                 */
    3438             : /*******************************************************************/
    3439             : 
    3440             : static GEN
    3441          14 : RgXQ_charpoly_FpXQ(GEN x, GEN T, GEN p, long v)
    3442             : {
    3443          14 :   pari_sp av = avma;
    3444             :   GEN r;
    3445          14 :   if (lgefint(p)==3)
    3446             :   {
    3447           0 :     ulong pp = p[2];
    3448           0 :     r = Flx_to_ZX(Flxq_charpoly(RgX_to_Flx(x, pp), RgX_to_Flx(T, pp), pp));
    3449             :   }
    3450             :   else
    3451          14 :     r = FpXQ_charpoly(RgX_to_FpX(x, p), RgX_to_FpX(T, p), p);
    3452          14 :   r = FpX_to_mod(r, p); setvarn(r, v);
    3453          14 :   return gc_upto(av, r);
    3454             : }
    3455             : 
    3456             : static GEN
    3457       12980 : RgXQ_charpoly_fast(GEN x, GEN T, long v)
    3458             : {
    3459             :   GEN p, pol;
    3460       12980 :   long pa, t = RgX_type2(x,T, &p,&pol,&pa);
    3461       12980 :   switch(t)
    3462             :   {
    3463        9431 :     case t_INT:    return ZXQ_charpoly(x, T, v);
    3464        2163 :     case t_FRAC:
    3465             :     {
    3466        2163 :       pari_sp av = avma;
    3467             :       GEN cT;
    3468        2163 :       T = Q_primitive_part(T, &cT);
    3469        2163 :       T = QXQ_charpoly(x, T, v);
    3470        2163 :       if (cT) T = gc_upto(av, T); /* silly rare case */
    3471        2163 :       return T;
    3472             :     }
    3473          14 :     case t_INTMOD: return RgXQ_charpoly_FpXQ(x, T, p, v);
    3474        1372 :     default:       return NULL;
    3475             :   }
    3476             : }
    3477             : 
    3478             : /* (v - x)^d */
    3479             : static GEN
    3480         126 : caract_const(pari_sp av, GEN x, long v, long d)
    3481         126 : { return gc_upto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    3482             : 
    3483             : GEN
    3484      977218 : RgXQ_charpoly_i(GEN x, GEN T, long v)
    3485             : {
    3486      977218 :   pari_sp av = avma;
    3487      977218 :   long d = degpol(T), dx = degpol(x), v0;
    3488             :   GEN ch, L;
    3489      977218 :   if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
    3490      977218 :   if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
    3491             : 
    3492      977148 :   v0 = fetch_var_higher();
    3493      977148 :   x = RgX_neg(x);
    3494      977148 :   gel(x,2) = gadd(gel(x,2), pol_x(v));
    3495      977144 :   setvarn(x, v0);
    3496      977144 :   T = leafcopy(T); setvarn(T, v0);
    3497      977141 :   ch = resultant(T, x);
    3498      977154 :   (void)delete_var();
    3499             :   /* test for silly input: x mod (deg 0 polynomial) */
    3500      977154 :   if (typ(ch) != t_POL)
    3501           7 :     pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
    3502      977147 :   L = leading_coeff(ch);
    3503      977147 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    3504      977147 :   return gc_upto(av, ch);
    3505             : }
    3506             : 
    3507             : /* return caract(Mod(x,T)) in variable v */
    3508             : GEN
    3509       12980 : RgXQ_charpoly(GEN x, GEN T, long v)
    3510             : {
    3511       12980 :   GEN ch = RgXQ_charpoly_fast(x, T, v);
    3512       12980 :   if (ch) return ch;
    3513        1372 :   return RgXQ_charpoly_i(x, T, v);
    3514             : }
    3515             : 
    3516             : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    3517             :  * algebra nf[t]/(Q(t)) */
    3518             : GEN
    3519         224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    3520             : {
    3521         224 :   const char *f = "rnfcharpoly";
    3522         224 :   long dQ = degpol(Q);
    3523         224 :   pari_sp av = avma;
    3524             :   GEN T;
    3525             : 
    3526         224 :   if (v < 0) v = 0;
    3527         224 :   nf = checknf(nf); T = nf_get_pol(nf);
    3528         224 :   Q = RgX_nffix(f, T,Q,0);
    3529         224 :   switch(typ(x))
    3530             :   {
    3531          28 :     case t_INT:
    3532          28 :     case t_FRAC: return caract_const(av, x, v, dQ);
    3533          91 :     case t_POLMOD:
    3534          91 :       x = polmod_nffix2(f,T,Q, x,0);
    3535          56 :       break;
    3536          56 :     case t_POL:
    3537          56 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    3538          42 :       break;
    3539          49 :     default: pari_err_TYPE(f,x);
    3540             :   }
    3541          98 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    3542             :   /* x a t_POL in variable vQ */
    3543          56 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    3544          56 :   if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
    3545          56 :   return gc_GEN(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    3546             : }
    3547             : 
    3548             : /*******************************************************************/
    3549             : /*                                                                 */
    3550             : /*                  GCD USING SUBRESULTANT                         */
    3551             : /*                                                                 */
    3552             : /*******************************************************************/
    3553             : static int inexact(GEN x, int *simple);
    3554             : static int
    3555        2198 : isinexactall(GEN x, int *simple)
    3556             : {
    3557        2198 :   long i, lx = lg(x);
    3558       13356 :   for (i=2; i<lx; i++)
    3559       11172 :     if (inexact(gel(x,i), simple)) return 1;
    3560        2184 :   return 0;
    3561             : }
    3562             : /* return 1 if coeff explosion is not possible */
    3563             : static int
    3564       11424 : inexact(GEN x, int *simple)
    3565             : {
    3566       11424 :   int junk = 0;
    3567       11424 :   switch(typ(x))
    3568             :   {
    3569        7721 :     case t_INT: case t_FRAC: return 0;
    3570             : 
    3571           7 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    3572             : 
    3573        2051 :     case t_INTMOD:
    3574             :     case t_FFELT:
    3575        2051 :       if (!*simple) *simple = 1;
    3576        2051 :       return 0;
    3577             : 
    3578          77 :     case t_COMPLEX:
    3579          77 :       return inexact(gel(x,1), simple)
    3580          77 :           || inexact(gel(x,2), simple);
    3581           0 :     case t_QUAD:
    3582           0 :       *simple = 0;
    3583           0 :       return inexact(gel(x,2), &junk)
    3584           0 :           || inexact(gel(x,3), &junk);
    3585             : 
    3586         819 :     case t_POLMOD:
    3587         819 :       return isinexactall(gel(x,1), simple);
    3588         700 :     case t_POL:
    3589         700 :       *simple = -1;
    3590         700 :       return isinexactall(x, &junk);
    3591          49 :     case t_RFRAC:
    3592          49 :       *simple = -1;
    3593          49 :       return inexact(gel(x,1), &junk)
    3594          49 :           || inexact(gel(x,2), &junk);
    3595             :   }
    3596           0 :   *simple = -1; return 0;
    3597             : }
    3598             : 
    3599             : /* x monomial, y t_POL in the same variable */
    3600             : static GEN
    3601        3717 : gcdmonome(GEN x, GEN y)
    3602             : {
    3603        3717 :   pari_sp av = avma;
    3604        3717 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    3605        3717 :   long i, l = lg(y);
    3606        3717 :   GEN t, v = cgetg(l, t_VEC);
    3607        3717 :   gel(v,1) = gel(x,dx+2);
    3608        7672 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    3609        3717 :   t = content(v); /* gcd(lc(x), cont(y)) */
    3610        3717 :   t = simplify_shallow(t);
    3611        3717 :   if (dx < e) e = dx;
    3612        3717 :   return gc_upto(av, monomialcopy(t, e, varn(x)));
    3613             : }
    3614             : 
    3615             : static GEN
    3616      111909 : RgX_gcd_FpX(GEN x, GEN y, GEN p)
    3617             : {
    3618      111909 :   pari_sp av = avma;
    3619             :   GEN r;
    3620      111909 :   if (lgefint(p) == 3)
    3621             :   {
    3622      111895 :     ulong pp = uel(p, 2);
    3623      111895 :     r = Flx_to_ZX_inplace(Flx_gcd(RgX_to_Flx(x, pp),
    3624             :                                   RgX_to_Flx(y, pp), pp));
    3625             :   }
    3626             :   else
    3627          14 :     r = FpX_gcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3628      111909 :   return gc_upto(av, FpX_to_mod(r, p));
    3629             : }
    3630             : 
    3631             : static GEN
    3632           7 : RgX_gcd_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3633             : {
    3634           7 :   pari_sp av = avma;
    3635           7 :   GEN r, T = RgX_to_FpX(pol, p);
    3636           7 :   if (signe(T)==0) pari_err_OP("gcd", x, y);
    3637           7 :   r = FpXQX_gcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3638           7 :   return gc_upto(av, FpXQX_to_mod(r, T, p));
    3639             : }
    3640             : 
    3641             : static GEN
    3642         595 : RgX_gcd_FpXk(GEN x, GEN y, GEN p)
    3643             : {
    3644         595 :   pari_sp av = avma;
    3645         595 :   GEN r = FpXk_gcd(Rg_to_FpXk(x, p), Rg_to_FpXk(y, p), p);
    3646         595 :   return gc_upto(av, gmul(r, gmodulsg(1,p)));
    3647             : }
    3648             : 
    3649             : static GEN
    3650       10878 : RgX_liftred(GEN x, GEN T)
    3651       10878 : { return RgXQX_red(liftpol_shallow(x), T); }
    3652             : 
    3653             : static GEN
    3654        2289 : RgX_gcd_ZXQX(GEN x, GEN y, GEN T)
    3655             : {
    3656        2289 :   pari_sp av = avma;
    3657        2289 :   GEN r = ZXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3658        2289 :   return gc_GEN(av, QXQX_to_mod_shallow(r, T));
    3659             : }
    3660             : 
    3661             : static GEN
    3662        3150 : RgX_gcd_QXQX(GEN x, GEN y, GEN T)
    3663             : {
    3664        3150 :   pari_sp av = avma;
    3665        3150 :   GEN r = QXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3666        3150 :   return gc_GEN(av, QXQX_to_mod_shallow(r, T));
    3667             : }
    3668             : 
    3669             : static GEN
    3670    10414600 : RgX_gcd_fast(GEN x, GEN y)
    3671             : {
    3672             :   GEN p, pol;
    3673             :   long pa;
    3674    10414600 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3675    10414600 :   switch(t)
    3676             :   {
    3677     8596223 :     case t_INT:    return ZX_gcd(x, y);
    3678        7525 :     case t_FRAC:   return QX_gcd(x, y);
    3679        2590 :     case t_FFELT:  return FFX_gcd(x, y, pol);
    3680      111909 :     case t_INTMOD: return RgX_gcd_FpX(x, y, p);
    3681           7 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    3682           7 :                    return RgX_gcd_FpXQX(x, y, pol, p);
    3683        2296 :     case RgX_type_code(t_POLMOD, t_INT):
    3684        2296 :                    return ZX_is_monic(pol)? RgX_gcd_ZXQX(x,y,pol): NULL;
    3685        3164 :     case RgX_type_code(t_POLMOD, t_FRAC):
    3686        6328 :                    return RgX_is_ZX(pol) && ZX_is_monic(pol) ?
    3687        6328 :                                             RgX_gcd_QXQX(x,y,pol): NULL;
    3688     1686063 :     case  RgX_type_code(t_POL, t_INT):
    3689     1686063 :                    return ZXk_gcd(x,y);
    3690         189 :     case  RgX_type_code(t_POL, t_FRAC):
    3691         189 :                    return QXk_gcd(x,y);
    3692         595 :     case  RgX_type_code(t_POL, t_INTMOD):
    3693         595 :                    return RgX_gcd_FpXk(x,y,p);
    3694        4039 :     default:       return NULL;
    3695             :   }
    3696             : }
    3697             : 
    3698             : /* x, y are t_POL in the same variable */
    3699             : GEN
    3700    10414600 : RgX_gcd(GEN x, GEN y)
    3701             : {
    3702             :   long dx, dy;
    3703             :   pari_sp av, av1;
    3704             :   GEN d, g, h, p1, p2, u, v;
    3705    10414600 :   int simple = 0;
    3706    10414600 :   GEN z = RgX_gcd_fast(x, y);
    3707    10414600 :   if (z) return z;
    3708        4060 :   if (isexactzero(y)) return RgX_copy(x);
    3709        4060 :   if (isexactzero(x)) return RgX_copy(y);
    3710        4060 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    3711         427 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    3712         343 :   if (isinexactall(x,&simple) || isinexactall(y,&simple))
    3713             :   {
    3714           7 :     av = avma; u = ggcd(content(x), content(y));
    3715           7 :     return gc_upto(av, scalarpol(u, varn(x)));
    3716             :   }
    3717             : 
    3718         336 :   av = avma;
    3719         336 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    3720             :   else
    3721             :   {
    3722         336 :     dx = lg(x); dy = lg(y);
    3723         336 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    3724         336 :     if (dy==3)
    3725             :     {
    3726           0 :       d = ggcd(gel(y,2), content(x));
    3727           0 :       return gc_upto(av, scalarpol(d, varn(x)));
    3728             :     }
    3729         336 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    3730         336 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    3731         336 :     d = ggcd(p1,p2);
    3732         336 :     av1 = avma;
    3733         336 :     g = h = gen_1;
    3734             :     for(;;)
    3735         301 :     {
    3736         637 :       GEN r = RgX_pseudorem(u,v);
    3737         637 :       long degq, du, dv, dr = lg(r);
    3738             : 
    3739         637 :       if (!signe(r)) break;
    3740         560 :       if (dr <= 3)
    3741             :       {
    3742         259 :         set_avma(av1);
    3743         259 :         return gc_upto(av, scalarpol(d, varn(x)));
    3744             :       }
    3745         301 :       du = lg(u); dv = lg(v); degq = du-dv;
    3746         301 :       u = v; p1 = g; g = leading_coeff(u);
    3747         301 :       switch(degq)
    3748             :       {
    3749          14 :         case 0: break;
    3750         280 :         case 1:
    3751         280 :           p1 = gmul(h,p1); h = g; break;
    3752           7 :         default:
    3753           7 :           p1 = gmul(gpowgs(h,degq), p1);
    3754           7 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    3755             :       }
    3756         301 :       v = RgX_Rg_div(r,p1);
    3757         301 :       if (gc_needed(av1,1))
    3758             :       {
    3759           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd, dr = %ld", degpol(r));
    3760           0 :         (void)gc_all(av1,4, &u,&v,&g,&h);
    3761             :       }
    3762             :     }
    3763          77 :     x = RgX_Rg_mul(primpart(v), d);
    3764             :   }
    3765          77 :   if (must_negate(x)) x = RgX_neg(x);
    3766          77 :   return gc_upto(av,x);
    3767             : }
    3768             : 
    3769             : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
    3770             : static GEN
    3771         399 : RgX_disc_i(GEN P)
    3772             : {
    3773         399 :   long n = degpol(P), dd;
    3774             :   GEN N, D, L, y;
    3775         399 :   if (!signe(P) || !n) return Rg_get_0(P);
    3776         392 :   if (n == 1) return Rg_get_1(P);
    3777         392 :   if (n == 2) {
    3778         126 :     GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
    3779         126 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    3780             :   }
    3781         266 :   y = RgX_deriv(P);
    3782         266 :   N = characteristic(P);
    3783         266 :   if (signe(N)) y = gmul(y, mkintmod(gen_1,N));
    3784         266 :   if (!signe(y)) return Rg_get_0(y);
    3785         266 :   dd = n - 2 - degpol(y);
    3786         266 :   if (isinexact(P))
    3787          21 :     D = resultant2(P,y);
    3788             :   else
    3789             :   {
    3790         245 :     D = RgX_resultant_all(P, y, NULL);
    3791         245 :     if (D == gen_0) return Rg_get_0(y);
    3792             :   }
    3793         266 :   L = leading_coeff(P);
    3794         266 :   if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
    3795         266 :   if (n & 2) D = gneg(D);
    3796         266 :   return D;
    3797             : }
    3798             : 
    3799             : static GEN
    3800          42 : RgX_disc_FpX(GEN x, GEN p)
    3801             : {
    3802          42 :   pari_sp av = avma;
    3803          42 :   GEN r = FpX_disc(RgX_to_FpX(x, p), p);
    3804          42 :   return gc_upto(av, Fp_to_mod(r, p));
    3805             : }
    3806             : 
    3807             : static GEN
    3808          28 : RgX_disc_FpXQX(GEN x, GEN pol, GEN p)
    3809             : {
    3810          28 :   pari_sp av = avma;
    3811          28 :   GEN r, T = RgX_to_FpX(pol, p);
    3812          28 :   r = FpXQX_disc(RgX_to_FpXQX(x, T, p), T, p);
    3813          28 :   return gc_upto(av, FpX_to_mod(r, p));
    3814             : }
    3815             : 
    3816             : static GEN
    3817      124140 : RgX_disc_fast(GEN x)
    3818             : {
    3819             :   GEN p, pol;
    3820             :   long pa;
    3821      124140 :   long t = RgX_type(x, &p,&pol,&pa);
    3822      124140 :   switch(t)
    3823             :   {
    3824      123629 :     case t_INT:    return ZX_disc(x);
    3825           7 :     case t_FRAC:   return QX_disc(x);
    3826          35 :     case t_FFELT:  return FFX_disc(x, pol);
    3827          42 :     case t_INTMOD: return RgX_disc_FpX(x, p);
    3828          28 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    3829          28 :                    return RgX_disc_FpXQX(x, pol, p);
    3830         399 :     default:       return NULL;
    3831             :   }
    3832             : }
    3833             : 
    3834             : GEN
    3835      124140 : RgX_disc(GEN x)
    3836             : {
    3837             :   pari_sp av;
    3838      124140 :   GEN z = RgX_disc_fast(x);
    3839      124140 :   if (z) return z;
    3840         399 :   av = avma;
    3841         399 :   return gc_upto(av, RgX_disc_i(x));
    3842             : }
    3843             : 
    3844             : GEN
    3845        4714 : poldisc0(GEN x, long v)
    3846             : {
    3847        4714 :   long v0, tx = typ(x);
    3848             :   pari_sp av;
    3849             :   GEN D;
    3850        4714 :   if (tx == t_POL && (v < 0 || v == varn(x))) return RgX_disc(x);
    3851          35 :   switch(tx)
    3852             :   {
    3853           0 :     case t_QUAD:
    3854           0 :       return quad_disc(x);
    3855           0 :     case t_POLMOD:
    3856           0 :       if (v >= 0 && varn(gel(x,1)) != v) break;
    3857           0 :       return RgX_disc(gel(x,1));
    3858          14 :     case t_QFB:
    3859          14 :       return icopy(qfb_disc(x));
    3860           0 :     case t_VEC: case t_COL: case t_MAT:
    3861           0 :       pari_APPLY_same(poldisc0(gel(x,i), v));
    3862             :   }
    3863          21 :   if (v < 0) pari_err_TYPE("poldisc",x);
    3864          21 :   av = avma; v0 = fetch_var_higher();
    3865          21 :   x = fix_pol(x,v, v0);
    3866          14 :   D = RgX_disc(x); (void)delete_var();
    3867          14 :   return gc_upto(av, D);
    3868             : }
    3869             : 
    3870             : GEN
    3871           7 : reduceddiscsmith(GEN x)
    3872             : {
    3873           7 :   long j, n = degpol(x);
    3874           7 :   pari_sp av = avma;
    3875             :   GEN xp, M;
    3876             : 
    3877           7 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    3878           7 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    3879           7 :   RgX_check_ZX(x,"poldiscreduced");
    3880           7 :   if (!gequal1(gel(x,n+2)))
    3881           0 :     pari_err_IMPL("nonmonic polynomial in poldiscreduced");
    3882           7 :   M = cgetg(n+1,t_MAT);
    3883           7 :   xp = ZX_deriv(x);
    3884          28 :   for (j=1; j<=n; j++)
    3885             :   {
    3886          21 :     gel(M,j) = RgX_to_RgC(xp, n);
    3887          21 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    3888             :   }
    3889           7 :   return gc_upto(av, ZM_snf(M));
    3890             : }
    3891             : 
    3892             : /***********************************************************************/
    3893             : /**                                                                   **/
    3894             : /**                       STURM ALGORITHM                             **/
    3895             : /**              (number of real roots of x in [a,b])                 **/
    3896             : /**                                                                   **/
    3897             : /***********************************************************************/
    3898             : static GEN
    3899         525 : R_to_Q_up(GEN x)
    3900             : {
    3901             :   long e;
    3902         525 :   switch(typ(x))
    3903             :   {
    3904         525 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3905           0 :     case t_REAL:
    3906           0 :       x = mantissa_real(x,&e);
    3907           0 :       return gmul2n(addiu(x,1), -e);
    3908           0 :     default: pari_err_TYPE("R_to_Q_up", x);
    3909             :       return NULL; /* LCOV_EXCL_LINE */
    3910             :   }
    3911             : }
    3912             : static GEN
    3913         525 : R_to_Q_down(GEN x)
    3914             : {
    3915             :   long e;
    3916         525 :   switch(typ(x))
    3917             :   {
    3918         525 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3919           0 :     case t_REAL:
    3920           0 :       x = mantissa_real(x,&e);
    3921           0 :       return gmul2n(subiu(x,1), -e);
    3922           0 :     default: pari_err_TYPE("R_to_Q_down", x);
    3923             :       return NULL; /* LCOV_EXCL_LINE */
    3924             :   }
    3925             : }
    3926             : 
    3927             : static long
    3928        1148 : sturmpart_i(GEN x, GEN ab)
    3929             : {
    3930        1148 :   long tx = typ(x);
    3931        1148 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    3932        1148 :   if (tx != t_POL)
    3933             :   {
    3934           0 :     if (is_real_t(tx)) return 0;
    3935           0 :     pari_err_TYPE("sturm",x);
    3936             :   }
    3937        1148 :   if (lg(x) == 3) return 0;
    3938        1148 :   if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
    3939        1148 :   (void)ZX_gcd_all(x, ZX_deriv(x), &x);
    3940        1148 :   if (ab)
    3941             :   {
    3942             :     GEN A, B;
    3943         525 :     if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
    3944         525 :     A = R_to_Q_down(gel(ab,1));
    3945         525 :     B = R_to_Q_up(gel(ab,2));
    3946         525 :     ab = mkvec2(A, B);
    3947             :   }
    3948        1148 :   return ZX_sturmpart(x, ab);
    3949             : }
    3950             : /* Deprecated: RgX_sturmpart() should be preferred  */
    3951             : long
    3952         385 : sturmpart(GEN x, GEN a, GEN b)
    3953             : {
    3954         385 :   pari_sp av = avma;
    3955         385 :   if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
    3956         385 :   if (!a) a = mkmoo();
    3957         385 :   if (!b) b = mkoo();
    3958         385 :   return gc_long(av, sturmpart_i(x, mkvec2(a,b)));
    3959             : }
    3960             : long
    3961         763 : RgX_sturmpart(GEN x, GEN ab)
    3962         763 : { pari_sp av = avma; return gc_long(av, sturmpart_i(x, ab)); }
    3963             : 
    3964             : /***********************************************************************/
    3965             : /**                                                                   **/
    3966             : /**                        GENERIC EXTENDED GCD                       **/
    3967             : /**                                                                   **/
    3968             : /***********************************************************************/
    3969             : /* assume typ(x) = typ(y) = t_POL */
    3970             : static GEN
    3971         848 : RgXQ_inv_i(GEN x, GEN y)
    3972             : {
    3973         848 :   long vx=varn(x), vy=varn(y);
    3974             :   pari_sp av;
    3975             :   GEN u, v, d;
    3976             : 
    3977         848 :   while (vx != vy)
    3978             :   {
    3979           0 :     if (varncmp(vx,vy) > 0)
    3980             :     {
    3981           0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3982           0 :       return scalarpol(d, vy);
    3983             :     }
    3984           0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3985           0 :     x = gel(x,2); vx = gvar(x);
    3986             :   }
    3987         848 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3988         848 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3989         848 :   d = gdiv(u,d);
    3990         848 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3991         848 :   return gc_upto(av, d);
    3992             : }
    3993             : 
    3994             : /*Assume x is a polynomial and y is not */
    3995             : static GEN
    3996         112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3997             : {
    3998         112 :   long vx = varn(x);
    3999         112 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    4000         112 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    4001          84 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    4002          56 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    4003             : }
    4004             : /* Assume x==0, y!=0 */
    4005             : static GEN
    4006          63 : zero_bezout(GEN y, GEN *U, GEN *V)
    4007             : {
    4008          63 :   *U=gen_0; *V = ginv(y); return gen_1;
    4009             : }
    4010             : 
    4011             : GEN
    4012         427 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    4013             : {
    4014         427 :   long tx=typ(x), ty=typ(y), vx;
    4015         427 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    4016         392 :   if (tx != t_POL)
    4017             :   {
    4018         140 :     if (ty == t_POL)
    4019          56 :       return scalar_bezout(y,x,v,u);
    4020             :     else
    4021             :     {
    4022          84 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    4023          84 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    4024          63 :       if (xis0) return zero_bezout(y,u,v);
    4025          42 :       else return zero_bezout(x,v,u);
    4026             :     }
    4027             :   }
    4028         252 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    4029         196 :   vx = varn(x);
    4030         196 :   if (vx != varn(y))
    4031           0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    4032           0 :                                    : scalar_bezout(y,x,v,u);
    4033         196 :   return RgX_extgcd(x,y,u,v);
    4034             : }
    4035             : 
    4036             : GEN
    4037         427 : gcdext0(GEN x, GEN y)
    4038             : {
    4039         427 :   GEN z=cgetg(4,t_VEC);
    4040         427 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    4041         427 :   return z;
    4042             : }
    4043             : 
    4044             : /*******************************************************************/
    4045             : /*                                                                 */
    4046             : /*                    GENERIC (modular) INVERSE                    */
    4047             : /*                                                                 */
    4048             : /*******************************************************************/
    4049             : 
    4050             : GEN
    4051       35063 : ginvmod(GEN x, GEN y)
    4052             : {
    4053       35063 :   long tx=typ(x);
    4054             : 
    4055       35063 :   switch(typ(y))
    4056             :   {
    4057       35063 :     case t_POL:
    4058       35063 :       if (tx==t_POL) return RgXQ_inv(x,y);
    4059       13629 :       if (is_scalar_t(tx)) return ginv(x);
    4060           0 :       break;
    4061           0 :     case t_INT:
    4062           0 :       if (tx==t_INT) return Fp_inv(x,y);
    4063           0 :       if (tx==t_POL) return gen_0;
    4064             :   }
    4065           0 :   pari_err_TYPE2("ginvmod",x,y);
    4066             :   return NULL; /* LCOV_EXCL_LINE */
    4067             : }
    4068             : 
    4069             : /***********************************************************************/
    4070             : /**                                                                   **/
    4071             : /**                         NEWTON POLYGON                            **/
    4072             : /**                                                                   **/
    4073             : /***********************************************************************/
    4074             : 
    4075             : /* assume leading coeff of x is nonzero */
    4076             : GEN
    4077          28 : newtonpoly(GEN x, GEN p)
    4078             : {
    4079          28 :   pari_sp av = avma;
    4080             :   long n, ind, a, b;
    4081             :   GEN y, vval;
    4082             : 
    4083          28 :   if (typ(x) != t_POL) pari_err_TYPE("newtonpoly",x);
    4084          28 :   n = degpol(x); if (n<=0) return cgetg(1,t_VEC);
    4085          28 :   vval = new_chunk(n+1);
    4086          28 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    4087         168 :   for (a = 0; a <= n; a++) vval[a] = gvaluation(gel(x,a),p);
    4088          42 :   for (a = 0, ind = 1; a < n; a++)
    4089             :   {
    4090          42 :     if (vval[a] != LONG_MAX) break;
    4091          14 :     gel(y,ind++) = mkoo();
    4092             :   }
    4093          84 :   for (b = a+1; b <= n; a = b, b = a+1)
    4094             :   {
    4095             :     long u1, u2, c;
    4096          70 :     while (vval[b] == LONG_MAX) b++;
    4097          56 :     u1 = vval[a] - vval[b];
    4098          56 :     u2 = b - a;
    4099         154 :     for (c = b+1; c <= n; c++)
    4100             :     {
    4101             :       long r1, r2;
    4102          98 :       if (vval[c] == LONG_MAX) continue;
    4103          70 :       r1 = vval[a] - vval[c];
    4104          70 :       r2 = c - a;
    4105          70 :       if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
    4106             :     }
    4107         154 :     while (ind <= b) gel(y,ind++) = sstoQ(u1,u2);
    4108             :   }
    4109          28 :   stackdummy((pari_sp)vval, av); return y;
    4110             : }
    4111             : 
    4112             : static GEN
    4113      274309 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
    4114             : {
    4115      274309 :   pari_sp av = avma;
    4116             :   GEN r;
    4117      274309 :   if (lgefint(p) == 3)
    4118             :   {
    4119      152402 :     ulong pp = uel(p, 2);
    4120      152402 :     r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
    4121             :                 RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
    4122             :   }
    4123             :   else
    4124      121907 :     r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
    4125      274309 :   return gc_upto(av, FpX_to_mod(r, p));
    4126             : }
    4127             : 
    4128             : static GEN
    4129          14 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
    4130             : {
    4131          14 :   pari_sp av = avma;
    4132             :   GEN r;
    4133          14 :   if (lgefint(p) == 3)
    4134             :   {
    4135           7 :     ulong pp = uel(p, 2);
    4136           7 :     r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
    4137             :                                    RgX_to_Flx(y, pp), pp));
    4138             :   }
    4139             :   else
    4140           7 :     r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    4141          14 :   return gc_upto(av, FpX_to_mod(r, p));
    4142             : }
    4143             : 
    4144             : static GEN
    4145       12054 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
    4146             : {
    4147       12054 :   pari_sp av = avma;
    4148             :   GEN r;
    4149       12054 :   if (lgefint(p) == 3)
    4150             :   {
    4151        6088 :     ulong pp = uel(p, 2);
    4152        6088 :     r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
    4153             :                                    RgX_to_Flx(y, pp), pp));
    4154             :   }
    4155             :   else
    4156        5966 :     r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    4157       12054 :   return gc_upto(av, FpX_to_mod(r, p));
    4158             : }
    4159             : 
    4160             : static GEN
    4161         385 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
    4162             : {
    4163         385 :   pari_sp av = avma;
    4164             :   GEN r;
    4165         385 :   GEN T = RgX_to_FpX(pol, p);
    4166         385 :   if (signe(T)==0) pari_err_OP("*",x,y);
    4167         385 :   if (lgefint(p) == 3)
    4168             :   {
    4169         241 :     ulong pp = uel(p, 2);
    4170         241 :     GEN Tp = ZX_to_Flx(T, pp);
    4171         241 :     r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
    4172             :                                RgX_to_FlxqX(y, Tp, pp),
    4173             :                                RgX_to_FlxqX(S, Tp, pp), Tp, pp));
    4174             :   }
    4175             :   else
    4176         144 :     r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
    4177             :                    RgX_to_FpXQX(S, T, p), T, p);
    4178         385 :   return gc_upto(av, FpXQX_to_mod(r, T, p));
    4179             : }
    4180             : 
    4181             : static GEN
    4182           0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    4183             : {
    4184           0 :   pari_sp av = avma;
    4185             :   GEN r;
    4186           0 :   GEN T = RgX_to_FpX(pol, p);
    4187           0 :   if (signe(T)==0) pari_err_OP("*",x,x);
    4188           0 :   if (lgefint(p) == 3)
    4189             :   {
    4190           0 :     ulong pp = uel(p, 2);
    4191           0 :     GEN Tp = ZX_to_Flx(T, pp);
    4192           0 :     r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
    4193             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    4194             :   }
    4195             :   else
    4196           0 :     r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    4197           0 :   return gc_upto(av, FpXQX_to_mod(r, T, p));
    4198             : }
    4199             : 
    4200             : static GEN
    4201           7 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    4202             : {
    4203           7 :   pari_sp av = avma;
    4204             :   GEN r;
    4205           7 :   GEN T = RgX_to_FpX(pol, p);
    4206           7 :   if (signe(T)==0) pari_err_OP("^",x,gen_m1);
    4207           7 :   if (lgefint(p) == 3)
    4208             :   {
    4209           7 :     ulong pp = uel(p, 2);
    4210           7 :     GEN Tp = ZX_to_Flx(T, pp);
    4211           7 :     r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
    4212             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    4213             :   }
    4214             :   else
    4215           0 :     r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    4216           7 :   return gc_upto(av, FpXQX_to_mod(r, T, p));
    4217             : }
    4218             : 
    4219             : static GEN
    4220     1527549 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
    4221             : {
    4222             :   GEN p, pol;
    4223             :   long pa;
    4224     1527549 :   long t = RgX_type3(x,y,T, &p,&pol,&pa);
    4225     1527550 :   switch(t)
    4226             :   {
    4227      582328 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
    4228      631592 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
    4229         105 :     case t_FFELT:  return FFXQ_mul(x, y, T, pol);
    4230      274309 :     case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
    4231         385 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    4232         385 :                    return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
    4233       38831 :     default:       return NULL;
    4234             :   }
    4235             : }
    4236             : 
    4237             : GEN
    4238     1527549 : RgXQ_mul(GEN x, GEN y, GEN T)
    4239             : {
    4240     1527549 :   GEN z = RgXQ_mul_fast(x, y, T);
    4241     1527545 :   if (!z) z = RgX_rem(RgX_mul(x, y), T);
    4242     1527545 :   return z;
    4243             : }
    4244             : 
    4245             : static GEN
    4246      458340 : RgXQ_sqr_fast(GEN x, GEN T)
    4247             : {
    4248             :   GEN p, pol;
    4249             :   long pa;
    4250      458340 :   long t = RgX_type2(x, T, &p,&pol,&pa);
    4251      458340 :   switch(t)
    4252             :   {
    4253      111565 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
    4254      339924 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
    4255           7 :     case t_FFELT:  return FFXQ_sqr(x, T, pol);
    4256          14 :     case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
    4257           0 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    4258           0 :                    return RgXQ_sqr_FpXQXQ(x, T, pol, p);
    4259        6830 :     default:       return NULL;
    4260             :   }
    4261             : }
    4262             : 
    4263             : GEN
    4264      458340 : RgXQ_sqr(GEN x, GEN T)
    4265             : {
    4266      458340 :   GEN z = RgXQ_sqr_fast(x, T);
    4267      458340 :   if (!z) z = RgX_rem(RgX_sqr(x), T);
    4268      458340 :   return z;
    4269             : }
    4270             : 
    4271             : static GEN
    4272      135069 : RgXQ_inv_fast(GEN x, GEN y)
    4273             : {
    4274             :   GEN p, pol;
    4275             :   long pa;
    4276      135069 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    4277      135069 :   switch(t)
    4278             :   {
    4279       90240 :     case t_INT:    return QXQ_inv(x,y);
    4280       31913 :     case t_FRAC:   return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
    4281          14 :     case t_FFELT:  return FFXQ_inv(x, y, pol);
    4282       12054 :     case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
    4283           7 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    4284           7 :                    return RgXQ_inv_FpXQXQ(x, y, pol, p);
    4285         841 :     default:       return NULL;
    4286             :   }
    4287             : }
    4288             : 
    4289             : GEN
    4290      135069 : RgXQ_inv(GEN x, GEN y)
    4291             : {
    4292      135069 :   GEN z = RgXQ_inv_fast(x, y);
    4293      135055 :   if (!z) z = RgXQ_inv_i(x, y);
    4294      135055 :   return z;
    4295             : }

Generated by: LCOV version 1.16