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 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.12.0 lcov report (development 22872-edcf83abb) Lines: 1831 2011 91.0 %
Date: 2018-07-20 05:36:03 Functions: 173 179 96.6 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /***********************************************************************/
      15             : /**                                                                   **/
      16             : /**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
      17             : /**                         (second part)                             **/
      18             : /**                                                                   **/
      19             : /***********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
      24             :  * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
      25             :  * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
      26             :  * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
      27             :  * Not memory clean in the latter case */
      28             : GEN
      29       34160 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      30             : {
      31       34160 :   long dP=degpol(P), i, k, m;
      32             :   pari_sp av1, av2;
      33             :   GEN s,y,P_lead;
      34             : 
      35       34160 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      36       34160 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      37       34160 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      38       34160 :   y = cgetg(n+2,t_COL);
      39       34160 :   if (y0)
      40             :   {
      41       11375 :     if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
      42       11375 :     m = lg(y0)-1;
      43       11375 :     for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
      44             :   }
      45             :   else
      46             :   {
      47       22785 :     m = 1;
      48       22785 :     gel(y,1) = stoi(dP);
      49             :   }
      50       34160 :   P += 2; /* strip codewords */
      51             : 
      52       34160 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      53       34160 :   if (P_lead)
      54             :   {
      55           7 :     if (N) P_lead = Fq_inv(P_lead,T,N);
      56           7 :     else if (T) P_lead = QXQ_inv(P_lead,T);
      57             :   }
      58       88760 :   for (k=m; k<=n; k++)
      59             :   {
      60       54600 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      61      237496 :     for (i=1; i<k && i<=dP; i++)
      62      182896 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      63       54600 :     if (N)
      64             :     {
      65       15960 :       s = Fq_red(s, T, N);
      66       15960 :       if (P_lead) s = Fq_mul(s, P_lead, T, N);
      67             :     }
      68       38640 :     else if (T)
      69             :     {
      70           0 :       s = grem(s, T);
      71           0 :       if (P_lead) s = grem(gmul(s, P_lead), T);
      72             :     }
      73             :     else
      74       38640 :       if (P_lead) s = gdiv(s, P_lead);
      75       54600 :     av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
      76             :   }
      77       34160 :   return y;
      78             : }
      79             : 
      80             : GEN
      81       19544 : polsym(GEN x, long n)
      82             : {
      83       19544 :   return polsym_gen(x, NULL, n, NULL,NULL);
      84             : }
      85             : 
      86             : /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
      87             : GEN
      88    41847626 : centermodii(GEN x, GEN p, GEN po2)
      89             : {
      90    41847626 :   GEN y = remii(x, p);
      91    41961460 :   switch(signe(y))
      92             :   {
      93     1986795 :     case 0: break;
      94    26146620 :     case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
      95    25986721 :       break;
      96    14155791 :     case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
      97    14088316 :       break;
      98             :   }
      99    41734086 :   return y;
     100             : }
     101             : 
     102             : static long
     103           0 : s_centermod(long x, ulong pp, ulong pps2)
     104             : {
     105           0 :   long y = x % (long)pp;
     106           0 :   if (y < 0) y += pp;
     107           0 :   return Fl_center(y, pp,pps2);
     108             : }
     109             : 
     110             : /* for internal use */
     111             : GEN
     112     6902834 : centermod_i(GEN x, GEN p, GEN ps2)
     113             : {
     114             :   long i, lx;
     115             :   pari_sp av;
     116             :   GEN y;
     117             : 
     118     6902834 :   if (!ps2) ps2 = shifti(p,-1);
     119     6902918 :   switch(typ(x))
     120             :   {
     121     3302810 :     case t_INT: return centermodii(x,p,ps2);
     122             : 
     123     2785541 :     case t_POL: lx = lg(x);
     124     2785541 :       y = cgetg(lx,t_POL); y[1] = x[1];
     125    20951103 :       for (i=2; i<lx; i++)
     126             :       {
     127    18164335 :         av = avma;
     128    18164335 :         gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
     129             :       }
     130     2786768 :       return normalizepol_lg(y, lx);
     131             : 
     132      813608 :     case t_COL: lx = lg(x);
     133      813608 :       y = cgetg(lx,t_COL);
     134      813607 :       for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
     135      813607 :       return y;
     136             : 
     137         959 :     case t_MAT: lx = lg(x);
     138         959 :       y = cgetg(lx,t_MAT);
     139         959 :       for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
     140         959 :       return y;
     141             : 
     142           0 :     case t_VECSMALL: lx = lg(x);
     143             :     {
     144           0 :       ulong pp = itou(p), pps2 = itou(ps2);
     145           0 :       y = cgetg(lx,t_VECSMALL);
     146           0 :       for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
     147           0 :       return y;
     148             :     }
     149             :   }
     150           0 :   return x;
     151             : }
     152             : 
     153             : GEN
     154     4349507 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
     155             : 
     156             : static GEN
     157          56 : RgX_Frobenius_deflate(GEN S, ulong p)
     158             : {
     159          56 :   GEN F = RgX_deflate(S, p);
     160          56 :   long i, l = lg(F);
     161         252 :   for (i=2; i<l; i++)
     162             :   {
     163          84 :     GEN Fi = gel(F,i), R;
     164          84 :     if (typ(Fi)==t_POL)
     165             :     {
     166          42 :       if (signe(RgX_deriv(Fi))==0)
     167          28 :         gel(F,i) = RgX_Frobenius_deflate(gel(F, i), p);
     168          28 :       else return NULL;
     169             :     }
     170          42 :     else if (ispower(Fi, utoi(p), &R))
     171          42 :       gel(F,i) = R;
     172           0 :     else return NULL;
     173             :   }
     174          42 :   return F;
     175             : }
     176             : 
     177             : 
     178             : static GEN
     179         154 : RgXY_squff(GEN f)
     180             : {
     181         154 :   long i, q, n = degpol(f);
     182         154 :   ulong p = itos_or_0(characteristic(f));
     183         154 :   GEN u = const_vec(n+1, pol_1(varn(f)));
     184         168 :   for(q = 1;;q *= p)
     185          14 :   {
     186         168 :     GEN t, v, tv, r = RgX_gcd(f, RgX_deriv(f));
     187         168 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
     188          35 :     t = RgX_div(f, r);
     189          35 :     if (degpol(t) > 0)
     190             :     {
     191             :       long j;
     192          14 :       for(j = 1;;j++)
     193             :       {
     194          21 :         v = RgX_gcd(r, t);
     195          14 :         tv = RgX_div(t, v);
     196          14 :         if (degpol(tv) > 0) gel(u, j*q) = tv;
     197          14 :         if (degpol(v) <= 0) break;
     198           7 :         r = RgX_div(r, v);
     199           7 :         t = v;
     200             :       }
     201           7 :       if (degpol(r) == 0) break;
     202             :     }
     203          28 :     if (!p) break;
     204          28 :     r = RgX_Frobenius_deflate(f, p);
     205          28 :     if (!r) { gel(u, q) = f; break; }
     206          14 :     f = r;
     207             :   }
     208         665 :   for (i = n; i; i--)
     209         665 :     if (degpol(gel(u,i))) break;
     210         154 :   setlg(u,i+1); return u;
     211             : }
     212             : 
     213             : /* Lmod contains modular factors of *F (NULL codes an empty slot: used factor)
     214             :  * Lfac accumulates irreducible factors as they are found.
     215             :  * p is a product of modular factors in Lmod[1..i-1] (NULL for p = 1), not
     216             :  * a rational factor of *F
     217             :  * Find an irreducible factor of *F divisible by p (by including
     218             :  * exhaustively further factors from Lmod[i..]); return 0 on failure, else 1.
     219             :  * Update Lmod, Lfac and *F */
     220             : static int
     221        8547 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
     222             : {
     223             :   GEN q;
     224        8547 :   if (i == lg(Lmod)) return 0;
     225        4466 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
     226        4235 :   if (!gel(Lmod,i)) return 0;
     227        4221 :   p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
     228        4221 :   q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
     229        4221 :   if (degpol(q))
     230             :   {
     231        3899 :     GEN R, Q = RgX_divrem(*F, q, &R);
     232        3899 :     if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
     233             :   }
     234        3920 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
     235        3668 :   return 0;
     236             : }
     237             : 
     238             : static GEN factor_domain(GEN x, GEN flag);
     239             : 
     240             : static GEN
     241         161 : RgXY_factor_squarefree(GEN f, GEN dom)
     242             : {
     243         161 :   pari_sp av = avma;
     244         161 :   GEN Lfac, Lmod, F = NULL, BLOC = NULL;
     245         161 :   long vy = gvar2(f), n = RgXY_degreex(f);
     246         161 :   ulong i, c = itou_or_0(residual_characteristic(f));
     247         161 :   long val = RgX_valrem(f, &f);
     248             :   for(;;)
     249             :   {
     250         245 :     for (i = 0; !c || i < c; i++)
     251             :     {
     252         245 :       BLOC = gpowgs(gaddgs(pol_x(vy), i), n+1);
     253         245 :       F = poleval(f, BLOC);
     254         245 :       if (issquarefree(c ? gmul(F,mkintmodu(1,c)): F)) break;
     255             :     }
     256         161 :     if (!c || i < c) break;
     257           0 :     n++;
     258             :   }
     259         161 :   if (DEBUGLEVEL >= 2)
     260           0 :     err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
     261         161 :   Lmod = gel(factor_domain(F,dom),1);
     262         161 :   if (DEBUGLEVEL >= 2)
     263           0 :     err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
     264         161 :   Lfac = vectrunc_init(lg(Lmod)+1);
     265         161 :   settyp(Lfac, t_COL);
     266         161 :   (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
     267         161 :   if (degpol(f)) vectrunc_append(Lfac, f);
     268         161 :   if (val)
     269             :   {
     270          21 :     GEN x = pol_x(varn(f)), c = NULL;
     271          21 :     if (dom) { c = Rg_get_1(dom); if (typ(c) == t_INT) c = NULL; }
     272          21 :     if (c) x = RgX_Rg_mul(x, c);
     273          21 :     vectrunc_append(Lfac, x);
     274             :   }
     275         161 :   return gerepilecopy(av, Lfac);
     276             : }
     277             : 
     278             : static GEN
     279         154 : FE_matconcat(GEN F, GEN E, long l)
     280             : {
     281         154 :   setlg(E,l); E = shallowconcat1(E);
     282         154 :   setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
     283             : }
     284             : 
     285             : static int
     286         315 : gen_cmp_RgXY(void *data, GEN x, GEN y)
     287             : {
     288         315 :   long vx = varn(x), vy = varn(y);
     289         315 :   return (vx == vy)? gen_cmp_RgX(data, x, y): -varncmp(vx, vy);
     290             : }
     291             : static GEN
     292         154 : RgXY_factor(GEN f, GEN dom)
     293             : {
     294         154 :   pari_sp av = avma;
     295             :   GEN C, F, E, cf, V;
     296             :   long i, j, l;
     297         154 :   if (dom) { GEN c = Rg_get_1(dom); if (typ(c) != t_INT) f = RgX_Rg_mul(f,c); }
     298         154 :   cf = content(f);
     299         154 :   V = RgXY_squff(gdiv(f, cf)); l = lg(V);
     300         154 :   C = factor_domain(cf, dom);
     301         154 :   F = cgetg(l+1, t_VEC); gel(F,1) = gel(C,1);
     302         154 :   E = cgetg(l+1, t_VEC); gel(E,1) = gel(C,2);
     303         399 :   for (i=1, j=2; i < l; i++)
     304             :   {
     305         245 :     GEN v = gel(V,i);
     306         245 :     if (degpol(v))
     307             :     {
     308         161 :       gel(F,j) = v = RgXY_factor_squarefree(v, dom);
     309         161 :       gel(E,j) = const_col(lg(v)-1, utoipos(i));
     310         161 :       j++;
     311             :     }
     312             :   }
     313         154 :   f = FE_matconcat(F,E,j);
     314         154 :   (void)sort_factor(f,(void*)cmp_universal, &gen_cmp_RgXY);
     315         154 :   return gerepilecopy(av, f);
     316             : }
     317             : 
     318             : /***********************************************************************/
     319             : /**                                                                   **/
     320             : /**                          FACTORIZATION                            **/
     321             : /**                                                                   **/
     322             : /***********************************************************************/
     323             : static long RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var);
     324             : #define assign_or_fail(x,y) { GEN __x = x;\
     325             :   if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
     326             : }
     327             : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
     328             : 
     329             : static const long tsh = 6;
     330             : static long
     331      312414 : code(long t1, long t2) { return (t1 << tsh) | t2; }
     332             : void
     333      390131 : RgX_type_decode(long x, long *t1, long *t2)
     334             : {
     335      390131 :   *t1 = x >> tsh;
     336      390131 :   *t2 = (x & ((1L<<tsh)-1));
     337      390131 : }
     338             : int
     339     2585282 : RgX_type_is_composite(long t) { return t >= tsh; }
     340             : 
     341             : static int
     342   565417532 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     343             : {
     344             :   long j;
     345   565417532 :   switch(typ(c))
     346             :   {
     347             :     case t_INT:
     348   460409014 :       break;
     349             :     case t_FRAC:
     350    14172966 :       t[1]=1; break;
     351             :       break;
     352             :     case t_REAL:
     353    29792007 :       update_prec(precision(c), pa);
     354    29792007 :       t[2]=1; break;
     355             :     case t_INTMOD:
     356    26014081 :       assign_or_fail(gel(c,1),p);
     357    26014081 :       t[3]=1; break;
     358             :     case t_FFELT:
     359     1358762 :       if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
     360     1358762 :       assign_or_fail(FF_p_i(c),p);
     361     1358762 :       t[5]=1; break;
     362             :     case t_COMPLEX:
     363    32050042 :       for (j=1; j<=2; j++)
     364             :       {
     365    21366697 :         GEN d = gel(c,j);
     366    21366697 :         switch(typ(d))
     367             :         {
     368             :           case t_INT: case t_FRAC:
     369     2355781 :             if (!*t2) *t2 = t_COMPLEX;
     370     2355781 :             t[1]=1; break;
     371             :           case t_REAL:
     372    19010895 :             update_prec(precision(d), pa);
     373    19010895 :             if (!*t2) *t2 = t_COMPLEX;
     374    19010895 :             t[2]=1; break;
     375             :           case t_INTMOD:
     376          14 :             assign_or_fail(gel(d,1),p);
     377          14 :             if (!signe(*p) || mod4(*p) != 3) return 0;
     378           7 :             if (!*t2) *t2 = t_COMPLEX;
     379           7 :             t[3]=1; break;
     380             :           case t_PADIC:
     381           7 :             update_prec(precp(d)+valp(d), pa);
     382           7 :             assign_or_fail(gel(d,2),p);
     383           7 :             if (!*t2) *t2 = t_COMPLEX;
     384           7 :             t[7]=1; break;
     385           0 :           default: return 0;
     386             :         }
     387             :       }
     388    10683345 :       if (!t[2]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
     389    10683345 :       break;
     390             :     case t_PADIC:
     391      208257 :       update_prec(precp(c)+valp(c), pa);
     392      208257 :       assign_or_fail(gel(c,2),p);
     393      208257 :       t[7]=1; break;
     394             :     case t_QUAD:
     395         343 :       assign_or_fail(gel(c,1),pol);
     396        1029 :       for (j=2; j<=3; j++)
     397             :       {
     398         686 :         GEN d = gel(c,j);
     399         686 :         switch(typ(d))
     400             :         {
     401             :           case t_INT: case t_FRAC:
     402         651 :             t[8]=1; break;
     403             :           case t_INTMOD:
     404          28 :             assign_or_fail(gel(d,1),p);
     405          28 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     406          28 :             t[3]=1; break;
     407             :           case t_PADIC:
     408           7 :             update_prec(precp(d)+valp(d), pa);
     409           7 :             assign_or_fail(gel(d,2),p);
     410           7 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     411           7 :             t[7]=1; break;
     412           0 :           default: return 0;
     413             :         }
     414             :       }
     415         343 :       break;
     416             :     case t_POLMOD:
     417     3339046 :       assign_or_fail(gel(c,1),pol);
     418     3339046 :       if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1))) return 0;
     419    20027220 :       for (j=1; j<=2; j++)
     420             :       {
     421             :         GEN pbis, polbis;
     422             :         long pabis;
     423     6676832 :         *t2 = t_POLMOD;
     424     6676832 :         switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
     425             :         {
     426     3682702 :           case t_INT: break;
     427      501080 :           case t_FRAC:   t[1]=1; break;
     428     2490957 :           case t_INTMOD: t[3]=1; break;
     429           7 :           case t_PADIC:  t[7]=1; update_prec(pabis,pa); break;
     430        4172 :           default: return 0;
     431             :         }
     432     6674746 :         if (pbis) assign_or_fail(pbis,p);
     433     6674746 :         if (polbis) assign_or_fail(polbis,pol);
     434             :       }
     435     3336778 :       break;
     436     5313598 :     case t_RFRAC: t[10] = 1;
     437     5313598 :       if (!settype(gel(c,1),t,p,pol,pa,ff,t2,var)) return 0;
     438     5313598 :       c = gel(c,2); /* fall through */
     439    19439571 :     case t_POL: t[10] = 1;
     440    19439571 :       if (!RgX_settype(c,t,p,pol,pa,ff,t2,var)) return 0;
     441    19439461 :       if (*var == NO_VARIABLE) { *var = varn(c); break; }
     442             :       /* if more than one free var, ensure varn() == *var fails. FIXME: should
     443             :        * keep the list of all variables, later t_POLMOD may cancel them */
     444    13358102 :       if (*var != varn(c)) *var = MAXVARN+1;
     445    13358102 :       break;
     446         133 :     default: return 0;
     447             :   }
     448   565415014 :   return 1;
     449             : }
     450             : /* t[0] unused. Other values, if set, indicate a coefficient of type
     451             :  * t[1] : t_FRAC
     452             :  * t[2] : t_REAL
     453             :  * t[3] : t_INTMOD
     454             :  * t[4] : Unused
     455             :  * t[5] : t_FFELT
     456             :  * t[6] : Unused
     457             :  * t[7] : t_PADIC
     458             :  * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
     459             :  * t[9]:  Unused
     460             :  * t[10]: t_POL (recursive factorisation) */
     461             : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
     462             :  * given by t) */
     463             : static long
     464    62369929 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
     465             : {
     466    62369929 :   if (t[10] && (!*pol || var!=varn(*pol))) return t_POL;
     467    56297705 :   if (t2) /* polmod/quad/complex of intmod/padic */
     468             :   {
     469     1482202 :     if (t[2] && (t[3]||t[7])) return 0;
     470     1482202 :     if (t[3]) return code(t2,t_INTMOD);
     471     1455455 :     if (t[7]) return code(t2,t_PADIC);
     472     1455420 :     if (t[2]) return t_COMPLEX;
     473      285415 :     if (t[1]) return code(t2,t_FRAC);
     474      178599 :     return code(t2,t_INT);
     475             :   }
     476    54815503 :   if (t[5]) /* ffelt */
     477             :   {
     478      215097 :     if (t[2]||t[8]||t[9]) return 0;
     479      215097 :     *pol=ff; return t_FFELT;
     480             :   }
     481    54600406 :   if (t[2]) /* inexact, real */
     482             :   {
     483     5071860 :     if (t[3]||t[7]||t[9]) return 0;
     484     5071860 :     return t_REAL;
     485             :   }
     486    49528546 :   if (t[10]) return t_POL;
     487    49528546 :   if (t[8]) return code(t_QUAD,t_INT);
     488    49528329 :   if (t[3]) return t_INTMOD;
     489    45558728 :   if (t[7]) return t_PADIC;
     490    45540276 :   if (t[1]) return t_FRAC;
     491    42565834 :   return t_INT;
     492             : }
     493             : 
     494             : static long
     495   125213886 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     496             : {
     497   125213886 :   long i, lx = lg(x);
     498   570762042 :   for (i=2; i<lx; i++)
     499   445548178 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     500   125213864 :   return 1;
     501             : }
     502             : 
     503             : static long
     504    20999309 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     505             : {
     506    20999309 :   long i, l = lg(x);
     507   135294745 :   for (i = 1; i<l; i++)
     508   114297410 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     509    20997335 :   return 1;
     510             : }
     511             : 
     512             : static long
     513     7836128 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     514             : {
     515     7836128 :   long i, l = lg(x);
     516    28570641 :   for (i = 1; i < l; i++)
     517    20736452 :     if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     518     7834189 :   return 1;
     519             : }
     520             : 
     521             : long
     522     9262114 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
     523             : {
     524     9262114 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0};
     525     9262114 :   long tx = typ(x), t2 = 0, var = NO_VARIABLE;
     526     9262114 :   GEN ff = NULL;
     527     9262114 :   *p = *pol = NULL; *pa = LONG_MAX;
     528     9262114 :   if (is_scalar_t(tx))
     529             :   {
     530      257063 :     if (tx == t_POLMOD) return 0;
     531      256671 :     if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     532             :   }
     533     9005051 :   else if (tx == t_MAT)
     534           7 :   { if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0; }
     535             :   else /* t_POL, t_SER */
     536     9005044 :     if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     537     9261722 :   return choosetype(t,t2,ff,pol,var);
     538             : }
     539             : 
     540             : long
     541     1951674 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     542             : {
     543     1951674 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     544     1951674 :   long t2 = 0, var = NO_VARIABLE;
     545     1951674 :   GEN ff = NULL;
     546     1951674 :   *p = *pol = NULL; *pa = LONG_MAX;
     547     1951674 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     548     1953160 :   return choosetype(t,t2,ff,pol,var);
     549             : }
     550             : 
     551             : long
     552         294 : RgX_Rg_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     553             : {
     554         294 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     555         294 :   long t2 = 0, var = NO_VARIABLE;
     556         294 :   GEN ff = NULL;
     557         294 :   *p = *pol = NULL; *pa = LONG_MAX;
     558         294 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     559         294 :   if (!settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     560         294 :   return choosetype(t,t2,ff,pol,var);
     561             : }
     562             : 
     563             : long
     564    46375732 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     565             : {
     566    46375732 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     567    46375732 :   long t2 = 0, var = NO_VARIABLE;
     568    46375732 :   GEN ff = NULL;
     569    46375732 :   *p = *pol = NULL; *pa = LONG_MAX;
     570    92752169 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     571    46376808 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     572    46376286 :   return choosetype(t,t2,ff,pol,var);
     573             : }
     574             : 
     575             : long
     576      689249 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
     577             : {
     578      689249 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     579      689249 :   long t2 = 0, var = NO_VARIABLE;
     580      689249 :   GEN ff = NULL;
     581      689249 :   *p = *pol = NULL; *pa = LONG_MAX;
     582     1378508 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     583     1378519 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
     584      689260 :       !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
     585      689261 :   return choosetype(t,t2,ff,pol,var);
     586             : }
     587             : 
     588             : long
     589       85002 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
     590             : {
     591       85002 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     592       85002 :   long t2 = 0, var = NO_VARIABLE;
     593       85002 :   GEN ff = NULL;
     594       85002 :   *p = *pol = NULL; *pa = LONG_MAX;
     595       85002 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     596       84036 :   return choosetype(t,t2,ff,pol,var);
     597             : }
     598             : 
     599             : long
     600      263249 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     601             : {
     602      263249 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     603      263249 :   long t2 = 0, var = NO_VARIABLE;
     604      263249 :   GEN ff = NULL;
     605      263249 :   *p = *pol = NULL; *pa = LONG_MAX;
     606      526106 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     607      263284 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     608      262822 :   return choosetype(t,t2,ff,pol,var);
     609             : }
     610             : 
     611             : long
     612     3744173 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     613             : {
     614     3744173 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     615     3744173 :   long t2 = 0, var = NO_VARIABLE;
     616     3744173 :   GEN ff = NULL;
     617     3744173 :   *p = *pol = NULL; *pa = LONG_MAX;
     618     7487870 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     619     3744278 :       !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     620     3743592 :   return choosetype(t,t2,ff,pol,var);
     621             : }
     622             : 
     623             : GEN
     624       37356 : factor0(GEN x, GEN flag)
     625             : {
     626             :   ulong B;
     627       37356 :   long tx = typ(x);
     628       37356 :   if (!flag) return factor(x);
     629         217 :   if ((tx != t_INT && tx!=t_FRAC) || typ(flag) != t_INT)
     630         175 :     return factor_domain(x, flag);
     631          42 :   if (signe(flag) < 0) pari_err_FLAG("factor");
     632          42 :   switch(lgefint(flag))
     633             :   {
     634           7 :     case 2: B = 0; break;
     635          35 :     case 3: B = flag[2]; break;
     636           0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     637             :              return NULL; /*LCOV_EXCL_LINE*/
     638             :   }
     639          42 :   return boundfact(x, B);
     640             : }
     641             : 
     642             : GEN
     643       83827 : deg1_from_roots(GEN L, long v)
     644             : {
     645       83827 :   long i, l = lg(L);
     646       83827 :   GEN z = cgetg(l,t_COL);
     647      191499 :   for (i=1; i<l; i++)
     648      107672 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     649       83827 :   return z;
     650             : }
     651             : GEN
     652         980 : roots_from_deg1(GEN x)
     653             : {
     654         980 :   long i,l = lg(x);
     655         980 :   GEN r = cgetg(l,t_VEC);
     656         980 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     657         980 :   return r;
     658             : }
     659             : 
     660             : static GEN
     661          42 : gauss_factor_p(GEN p)
     662             : {
     663          42 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     664          42 :   return mkcomplex(a, b);
     665             : }
     666             : 
     667             : static GEN
     668          49 : gauss_primpart(GEN x, GEN *c)
     669             : {
     670          49 :   GEN a = real_i(x), b = imag_i(x), n = gcdii(a, b);
     671          49 :   *c = n; if (n == gen_1) return x;
     672          49 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     673             : }
     674             : 
     675             : static GEN
     676          70 : gauss_primpart_try(GEN x, GEN c)
     677             : {
     678             :   GEN r, y;
     679          70 :   if (typ(x) == t_INT)
     680             :   {
     681          42 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     682             :   }
     683             :   else
     684             :   {
     685          28 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     686          28 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     687          14 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     688             :   }
     689          56 :   return y;
     690             : }
     691             : 
     692             : static int
     693          91 : gauss_cmp(GEN x, GEN y)
     694             : {
     695             :   int v;
     696          91 :   if (typ(x) != t_COMPLEX)
     697           0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     698          91 :   if (typ(y) != t_COMPLEX) return 1;
     699          63 :   v = cmpii(gel(x,2), gel(y,2));
     700          63 :   if (v) return v;
     701          28 :   return gcmp(gel(x,1), gel(y,1));
     702             : }
     703             : 
     704             : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
     705             : static GEN
     706         455 : gauss_normal(GEN x)
     707             : {
     708         455 :   if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
     709         385 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     710         385 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     711         385 :   return x;
     712             : }
     713             : 
     714             : static GEN
     715          49 : gauss_factor(GEN x)
     716             : {
     717          49 :   pari_sp av = avma;
     718          49 :   GEN a = real_i(x), b = imag_i(x), d = gen_1, n, y, fa, P, E, P2, E2;
     719          49 :   long t1 = typ(a);
     720          49 :   long t2 = typ(b), i, j, l, exp = 0;
     721          49 :   if (t1 == t_FRAC) d = gel(a,2);
     722          49 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     723          49 :   if (d == gen_1) y = x;
     724             :   else
     725             :   {
     726          21 :     y = gmul(x, d);
     727          21 :     a = real_i(y); t1 = typ(a);
     728          21 :     b = imag_i(y); t2 = typ(b);
     729             :   }
     730          49 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     731          49 :   y = gauss_primpart(y, &n);
     732          49 :   fa = factor(cxnorm(y));
     733          49 :   P = gel(fa,1);
     734          49 :   E = gel(fa,2); l = lg(P);
     735          49 :   P2 = cgetg(l, t_COL);
     736          49 :   E2 = cgetg(l, t_COL);
     737         105 :   for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
     738             :   { /* either p = 2 (ramified) or those factors split in Q(i) */
     739          56 :     GEN p = gel(P,i), w, w2, t, we, pe;
     740          56 :     long v, e = itos(gel(E,i));
     741          56 :     int is2 = absequaliu(p, 2);
     742          56 :     w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     743          56 :     w2 = gauss_normal( conj_i(w) );
     744             :     /* w * w2 * I^3 = p, w2 = conj(w) * I */
     745          56 :     pe = powiu(p, e);
     746          56 :     we = gpowgs(w, e);
     747          56 :     t = gauss_primpart_try( gmul(y, conj_i(we)), pe );
     748          56 :     if (t) y = t; /* y /= w^e */
     749             :     else {
     750             :       /* y /= conj(w)^e, should be y /= w2^e */
     751          14 :       y = gauss_primpart_try( gmul(y, we), pe );
     752          14 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     753             :     }
     754          56 :     gel(P,i) = w;
     755          56 :     v = Z_pvalrem(n, p, &n);
     756          56 :     if (v) {
     757           7 :       exp -= v; /* += 3*v mod 4 */
     758           7 :       if (is2) v <<= 1; /* 2 = w^2 I^3 */
     759             :       else {
     760           0 :         gel(P2,j) = w2;
     761           0 :         gel(E2,j) = utoipos(v); j++;
     762             :       }
     763           7 :       gel(E,i) = stoi(e + v);
     764             :     }
     765          56 :     v = Z_pvalrem(d, p, &d);
     766          56 :     if (v) {
     767           7 :       exp += v; /* -= 3*v mod 4 */
     768           7 :       if (is2) v <<= 1; /* 2 is ramified */
     769             :       else {
     770           7 :         gel(P2,j) = w2;
     771           7 :         gel(E2,j) = utoineg(v); j++;
     772             :       }
     773           7 :       gel(E,i) = stoi(e - v);
     774             :     }
     775          56 :     exp &= 3;
     776             :   }
     777          49 :   if (j > 1) {
     778           7 :     long k = 1;
     779           7 :     GEN P1 = cgetg(l, t_COL);
     780           7 :     GEN E1 = cgetg(l, t_COL);
     781             :     /* remove factors with exponent 0 */
     782          14 :     for (i = 1; i < l; i++)
     783           7 :       if (signe(gel(E,i)))
     784             :       {
     785           0 :         gel(P1,k) = gel(P,i);
     786           0 :         gel(E1,k) = gel(E,i);
     787           0 :         k++;
     788             :       }
     789           7 :     setlg(P1, k); setlg(E1, k);
     790           7 :     setlg(P2, j); setlg(E2, j);
     791           7 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     792             :   }
     793          49 :   if (!equali1(n) || !equali1(d))
     794             :   {
     795          28 :     GEN Fa = factor(Qdivii(n, d));
     796          28 :     P = gel(Fa,1); l = lg(P);
     797          28 :     E = gel(Fa,2);
     798          70 :     for (i = 1; i < l; i++)
     799             :     {
     800          42 :       GEN w, p = gel(P,i);
     801             :       long e;
     802             :       int is2;
     803          42 :       switch(mod4(p))
     804             :       {
     805          14 :         case 3: continue;
     806          14 :         case 2: is2 = 1; break;
     807          14 :         default:is2 = 0; break;
     808             :       }
     809          28 :       e = itos(gel(E,i));
     810          28 :       w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     811          28 :       gel(P,i) = w;
     812          28 :       if (is2)
     813          14 :         gel(E,i) = stoi(2*e);
     814             :       else
     815             :       {
     816          14 :         P = shallowconcat(P, gauss_normal( conj_i(w) ));
     817          14 :         E = shallowconcat(E, gel(E,i));
     818             :       }
     819          28 :       exp -= e; /* += 3*e mod 4 */
     820          28 :       exp &= 3;
     821             :     }
     822          28 :     gel(Fa,1) = P;
     823          28 :     gel(Fa,2) = E;
     824          28 :     fa = famat_mul_shallow(fa, Fa);
     825             :   }
     826          49 :   fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
     827             : 
     828          49 :   y = gmul(y, powIs(exp));
     829          49 :   if (!gequal1(y)) {
     830          35 :     gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
     831          35 :     gel(fa,2) = shallowconcat(gen_1,    gel(fa,2));
     832             :   }
     833          49 :   return gerepilecopy(av, fa);
     834             : }
     835             : 
     836             : GEN
     837          35 : Q_factor_limit(GEN x, ulong lim)
     838             : {
     839          35 :   pari_sp av = avma;
     840             :   GEN a, b;
     841          35 :   if (typ(x) == t_INT) return Z_factor_limit(x, lim);
     842          14 :   a = Z_factor_limit(gel(x,1), lim);
     843          14 :   b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
     844          14 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     845             : }
     846             : GEN
     847        2646 : Q_factor(GEN x)
     848             : {
     849        2646 :   pari_sp av = avma;
     850             :   GEN a, b;
     851        2646 :   if (typ(x) == t_INT) return Z_factor(x);
     852         189 :   a = Z_factor(gel(x,1));
     853         189 :   b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
     854         189 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     855             : }
     856             : /* replace t_QUAD/t_COMPLEX coeffs by t_POLMOD in T */
     857             : static GEN
     858         126 : RgX_fix_quad(GEN x, GEN T)
     859             : {
     860         126 :   long i, l, v = varn(T);
     861         126 :   GEN y = cgetg_copy(x,&l);
     862         630 :   for (i = 2; i < l; i++)
     863             :   {
     864         504 :     GEN c = gel(x,i);
     865         504 :     switch(typ(c))
     866             :     {
     867          56 :       case t_QUAD: c++;/* fall through */
     868          98 :       case t_COMPLEX: c = deg1pol_shallow(gel(c,2),gel(c,1),v);
     869             :     }
     870         504 :     gel(y,i) = c;
     871             :   }
     872         126 :   y[1] = x[1]; return y;
     873             : }
     874             : 
     875             : static GEN
     876        2282 : RgX_factor(GEN x, GEN dom)
     877             : {
     878             :   pari_sp av;
     879             :   long pa, v, lx, r1, i;
     880             :   GEN  p, pol, y, p1, p2;
     881        2282 :   long tx = dom ? RgX_Rg_type(x,dom,&p,&pol,&pa): RgX_type(x,&p,&pol,&pa);
     882        2282 :   switch(tx)
     883             :   {
     884           7 :     case 0: pari_err_IMPL("factor for general polynomials");
     885         154 :     case t_POL: return RgXY_factor(x, dom);
     886        1498 :     case t_INT: return ZX_factor(x);
     887           7 :     case t_FRAC: return QX_factor(x);
     888         175 :     case t_INTMOD: return factmod(x,p);
     889          42 :     case t_PADIC: return factorpadic(x,p,pa);
     890          98 :     case t_FFELT: return FFX_factor(x,pol);
     891             : 
     892          21 :     case t_COMPLEX: y = cgetg(3,t_MAT);
     893          21 :       av = avma; p1 = deg1_from_roots(roots(x,pa), varn(x));
     894          21 :       gel(y,1) = p1 = gerepileupto(av, p1);
     895          21 :       gel(y,2) = const_col(lg(p1)-1, gen_1); return y;
     896             : 
     897          28 :     case t_REAL: y=cgetg(3,t_MAT); v=varn(x);
     898          28 :       av=avma; p1=cleanroots(x,pa);
     899          28 :       lx = lg(p1);
     900          70 :       for (r1 = 1; r1 < lx; r1++)
     901          49 :         if (typ(gel(p1,r1)) == t_COMPLEX) break;
     902          28 :       lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
     903          70 :       for (i = 1; i < r1; i++)
     904          42 :         gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
     905          35 :       for (   ; i < lx; i++)
     906             :       {
     907           7 :         GEN a = gel(p1,2*i-r1);
     908           7 :         p = cgetg(5, t_POL); gel(p2,i) = p;
     909           7 :         p[1] = x[1];
     910           7 :         gel(p,2) = gnorm(a);
     911           7 :         gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
     912           7 :         gel(p,4) = gen_1;
     913             :       }
     914          28 :       gel(y,1) = gerepileupto(av,p2);
     915          28 :       gel(y,2) = const_col(lx-1, gen_1); return y;
     916             : 
     917             :     default:
     918             :     {
     919         252 :       GEN w = NULL, T = pol;
     920             :       long t1, t2;
     921         252 :       av = avma;
     922         252 :       RgX_type_decode(tx, &t1, &t2);
     923         252 :       if (t1 == t_COMPLEX) w = gen_I();
     924         203 :       else if (t1 == t_QUAD) w = mkquad(pol,gen_0,gen_1);
     925         252 :       if (w)
     926             :       { /* substitute I or w by t_POLMOD */
     927         126 :         T = leafcopy(pol); setvarn(T, fetch_var());
     928         126 :         x = RgX_fix_quad(x, T);
     929             :       }
     930         252 :       switch (t2)
     931             :       {
     932         161 :         case t_INT: case t_FRAC: p1 = nffactor(T,x); break;
     933             :         case t_INTMOD:
     934          56 :           T = RgX_to_FpX(T,p);
     935          56 :           if (FpX_is_irred(T,p)) { p1 = factmod(x,mkvec2(p,T)); break; }
     936             :         /*fall through*/
     937             :         default:
     938          56 :           if (w) (void)delete_var();
     939          56 :           pari_err_IMPL("factor for general polynomial");
     940             :           return NULL; /* LCOV_EXCL_LINE */
     941             :       }
     942         196 :       if (t1 == t_POLMOD) return gerepileupto(av, p1);
     943             :       /* substitute back I or w */
     944          98 :       gel(p1,1) = gsubst(liftpol_shallow(gel(p1,1)), varn(T), w);
     945          98 :       (void)delete_var(); return gerepilecopy(av, p1);
     946             :     }
     947             :   }
     948             : }
     949             : 
     950             : static GEN
     951       40801 : factor_domain(GEN x, GEN dom)
     952             : {
     953       40801 :   long tx = typ(x);
     954       40801 :   long tdom = dom ? typ(dom): 0;
     955             :   pari_sp av;
     956             : 
     957       40801 :   if (gequal0(x))
     958          63 :     switch(tx)
     959             :     {
     960             :       case t_INT:
     961             :       case t_COMPLEX:
     962             :       case t_POL:
     963          63 :       case t_RFRAC: return prime_fact(x);
     964           0 :       default: pari_err_TYPE("factor",x);
     965             :     }
     966       40738 :   av = avma;
     967       40738 :   switch(tx)
     968             :   {
     969        2233 :     case t_POL: return RgX_factor(x, dom);
     970             :     case t_RFRAC: {
     971          28 :       GEN a = gel(x,1), b = gel(x,2);
     972          28 :       GEN y = famat_inv_shallow(RgX_factor(b, dom));
     973          28 :       if (typ(a)==t_POL) y = famat_mul_shallow(RgX_factor(a, dom), y);
     974          28 :       return gerepilecopy(av, sort_factor_pol(y, cmp_universal));
     975             :     }
     976       38253 :     case t_INT:  if (tdom==0 || tdom==t_INT) return Z_factor(x);
     977         182 :     case t_FRAC: if (tdom==0 || tdom==t_INT) return Q_factor(x);
     978             :     case t_COMPLEX: /* fall through */
     979          49 :       if (tdom==0 || tdom==t_COMPLEX)
     980             :       {
     981          49 :         GEN y = gauss_factor(x); if (y) return y;
     982             :       }
     983             :       /* fall through */
     984             :   }
     985           0 :   pari_err_TYPE("factor",x);
     986             :   return NULL; /* LCOV_EXCL_LINE */
     987             : }
     988             : 
     989             : GEN
     990       40311 : factor(GEN x) { return factor_domain(x, NULL); }
     991             : 
     992             : /*******************************************************************/
     993             : /*                                                                 */
     994             : /*                     ROOTS --> MONIC POLYNOMIAL                  */
     995             : /*                                                                 */
     996             : /*******************************************************************/
     997             : static GEN
     998      689475 : normalized_mul(void *E, GEN x, GEN y)
     999             : {
    1000      689475 :   long a = gel(x,1)[1], b = gel(y,1)[1];
    1001             :   (void) E;
    1002     1378950 :   return mkvec2(mkvecsmall(a + b),
    1003     1378950 :                 RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
    1004             : }
    1005             : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
    1006             : static GEN
    1007      375492 : normalized_to_RgX(GEN L)
    1008             : {
    1009      375492 :   long i, a = gel(L,1)[1];
    1010      375492 :   GEN A = gel(L,2);
    1011      375492 :   GEN z = cgetg(a + 3, t_POL);
    1012      375492 :   z[1] = evalsigne(1) | evalvarn(varn(A));
    1013      375492 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
    1014      375492 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
    1015      375492 :   gel(z,i) = gen_1; return z;
    1016             : }
    1017             : 
    1018             : /* compute prod (x - a[i]) */
    1019             : GEN
    1020      331548 : roots_to_pol(GEN a, long v)
    1021             : {
    1022      331548 :   pari_sp av = avma;
    1023      331548 :   long i, k, lx = lg(a);
    1024             :   GEN L;
    1025      331548 :   if (lx == 1) return pol_1(v);
    1026      331513 :   L = cgetg(lx, t_VEC);
    1027      700422 :   for (k=1,i=1; i<lx-1; i+=2)
    1028             :   {
    1029      368909 :     GEN s = gel(a,i), t = gel(a,i+1);
    1030      368909 :     GEN x0 = gmul(s,t);
    1031      368909 :     GEN x1 = gneg(gadd(s,t));
    1032      368909 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1033             :   }
    1034      645321 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
    1035      313808 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
    1036      331513 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1037      331513 :   return gerepileupto(av, normalized_to_RgX(L));
    1038             : }
    1039             : 
    1040             : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
    1041             : GEN
    1042       43979 : roots_to_pol_r1(GEN a, long v, long r1)
    1043             : {
    1044       43979 :   pari_sp av = avma;
    1045       43979 :   long i, k, lx = lg(a);
    1046             :   GEN L;
    1047       43979 :   if (lx == 1) return pol_1(v);
    1048       43979 :   L = cgetg(lx, t_VEC);
    1049      300076 :   for (k=1,i=1; i<r1; i+=2)
    1050             :   {
    1051      256097 :     GEN s = gel(a,i), t = gel(a,i+1);
    1052      256097 :     GEN x0 = gmul(s,t);
    1053      256097 :     GEN x1 = gneg(gadd(s,t));
    1054      256097 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1055             :   }
    1056       64475 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
    1057       20496 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
    1058      149636 :   for (i=r1+1; i<lx; i++)
    1059             :   {
    1060      105657 :     GEN s = gel(a,i);
    1061      105657 :     GEN x0 = gnorm(s);
    1062      105657 :     GEN x1 = gneg(gtrace(s));
    1063      105657 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1064             :   }
    1065       43979 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1066       43979 :   return gerepileupto(av, normalized_to_RgX(L));
    1067             : }
    1068             : 
    1069             : /*******************************************************************/
    1070             : /*                                                                 */
    1071             : /*                          FACTORBACK                             */
    1072             : /*                                                                 */
    1073             : /*******************************************************************/
    1074             : static GEN
    1075           7 : idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
    1076             : static GEN
    1077         609 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
    1078             : static GEN
    1079        3339 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
    1080             : static GEN
    1081        5621 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
    1082             : static GEN
    1083        4234 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
    1084             : static GEN
    1085       17080 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
    1086             : static GEN
    1087    54169771 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
    1088             : static GEN
    1089    78793168 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
    1090             : static GEN
    1091    16398396 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
    1092             : static GEN
    1093      121716 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
    1094             : 
    1095             : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
    1096             : GEN
    1097    28182332 : gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
    1098             :                GEN (*_pow)(void*,GEN,GEN), void *data)
    1099             : {
    1100    28182332 :   pari_sp av = avma;
    1101             :   long k, l, lx;
    1102             :   GEN p,x;
    1103             : 
    1104    28182332 :   if (e) /* supplied vector of exponents */
    1105      860937 :     p = L;
    1106             :   else
    1107             :   {
    1108    27321395 :     switch(typ(L)) {
    1109             :       case t_VEC:
    1110             :       case t_COL: /* product of the L[i] */
    1111     3398764 :         return gerepileupto(av, gen_product(L, data, _mul));
    1112             :       case t_MAT: /* genuine factorization */
    1113    23922631 :         l = lg(L);
    1114    23922631 :         if (l == 3) break;
    1115             :         /*fall through*/
    1116             :       default:
    1117           7 :         pari_err_TYPE("factorback [not a factorization]", L);
    1118             :     }
    1119    23922624 :     p = gel(L,1);
    1120    23922624 :     e = gel(L,2);
    1121             :   }
    1122             :   /* p = elts, e = expo */
    1123    24783561 :   lx = lg(p);
    1124             :   /* check whether e is an integral vector of correct length */
    1125    24783561 :   switch(typ(e))
    1126             :   {
    1127             :     case t_VECSMALL:
    1128         343 :       if (lx != lg(e))
    1129           0 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1130         343 :       if (lx == 1) return gen_1;
    1131         266 :       x = cgetg(lx,t_VEC);
    1132         882 :       for (l=1,k=1; k<lx; k++)
    1133         616 :         if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
    1134         266 :       break;
    1135             :     case t_VEC: case t_COL:
    1136    24783218 :       if (lx != lg(e) || !RgV_is_ZV(e))
    1137           7 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1138    24783211 :       if (lx == 1) return gen_1;
    1139    24763469 :       x = cgetg(lx,t_VEC);
    1140   103875025 :       for (l=1,k=1; k<lx; k++)
    1141    79111556 :         if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
    1142    24763469 :       break;
    1143             :     default:
    1144           0 :       pari_err_TYPE("factorback [not an exponent vector]", e);
    1145           0 :       return NULL;
    1146             :   }
    1147    24763735 :   x[0] = evaltyp(t_VEC) | _evallg(l);
    1148    24763735 :   return gerepileupto(av, gen_product(x, data, _mul));
    1149             : }
    1150             : 
    1151             : GEN
    1152        3052 : idealfactorback(GEN nf, GEN L, GEN e, int red)
    1153             : {
    1154        3052 :   nf = checknf(nf);
    1155        3052 :   if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
    1156        2450 :   else     return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
    1157             : }
    1158             : 
    1159             : GEN
    1160       13426 : nffactorback(GEN nf, GEN L, GEN e)
    1161       13426 : { return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }
    1162             : 
    1163             : GEN
    1164     3496377 : FpV_factorback(GEN L, GEN e, GEN p)
    1165     3496377 : { return gen_factorback(L, e, &Fpmul, &Fppow, (void*)p); }
    1166             : 
    1167             : GEN
    1168    24653244 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
    1169             : GEN
    1170     1401223 : factorback(GEN fa) { return factorback2(fa, NULL); }
    1171             : 
    1172             : GEN
    1173          14 : vecprod(GEN v)
    1174             : {
    1175          14 :   pari_sp av = avma;
    1176          14 :   if (!is_vec_t(typ(v)))
    1177           0 :     pari_err_TYPE("vecprod", v);
    1178          14 :   if (lg(v) == 1) return gen_1;
    1179           7 :   return gerepilecopy(av, gen_product(v, NULL, mul));
    1180             : }
    1181             : 
    1182             : static int
    1183          70 : RgX_is_irred_i(GEN x)
    1184             : {
    1185             :   GEN y, p, pol;
    1186          70 :   long l = lg(x), pa;
    1187             : 
    1188          70 :   if (!signe(x) || l <= 3) return 0;
    1189          70 :   switch(RgX_type(x,&p,&pol,&pa))
    1190             :   {
    1191          21 :     case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
    1192           0 :     case t_COMPLEX: return l == 4;
    1193             :     case t_REAL:
    1194           0 :       if (l == 4) return 1;
    1195           0 :       if (l > 5) return 0;
    1196           0 :       return gsigne(RgX_disc(x)) > 0;
    1197             :   }
    1198          49 :   y = factor(x);
    1199          49 :   return (lg(gcoeff(y,1,1))==l);
    1200             : }
    1201             : static int
    1202          70 : RgX_is_irred(GEN x)
    1203             : {
    1204          70 :   pari_sp av = avma;
    1205          70 :   int r = RgX_is_irred_i(x);
    1206          70 :   avma = av; return r;
    1207             : }
    1208             : long
    1209          70 : isirreducible(GEN x)
    1210             : {
    1211          70 :   switch(typ(x))
    1212             :   {
    1213           0 :     case t_INT: case t_REAL: case t_FRAC: return 0;
    1214          70 :     case t_POL: return RgX_is_irred(x);
    1215             :   }
    1216           0 :   pari_err_TYPE("isirreducible",x);
    1217           0 :   return 0;
    1218             : }
    1219             : 
    1220             : /*******************************************************************/
    1221             : /*                                                                 */
    1222             : /*                         GENERIC GCD                             */
    1223             : /*                                                                 */
    1224             : /*******************************************************************/
    1225             : /* x is a COMPLEX or a QUAD */
    1226             : static GEN
    1227        2499 : triv_cont_gcd(GEN x, GEN y)
    1228             : {
    1229        2499 :   pari_sp av = avma;
    1230             :   GEN c;
    1231        2499 :   if (typ(x)==t_COMPLEX)
    1232             :   {
    1233        2177 :     GEN a = gel(x,1), b = gel(x,2);
    1234        2177 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
    1235          21 :     c = ggcd(a,b);
    1236             :   }
    1237             :   else
    1238         322 :     c = ggcd(gel(x,2),gel(x,3));
    1239         343 :   return gerepileupto(av, ggcd(c,y));
    1240             : }
    1241             : 
    1242             : /* y is a PADIC, x a rational number or an INTMOD */
    1243             : static GEN
    1244        2684 : padic_gcd(GEN x, GEN y)
    1245             : {
    1246        2684 :   GEN p = gel(y,2);
    1247        2684 :   long v = gvaluation(x,p), w = valp(y);
    1248        2684 :   if (w < v) v = w;
    1249        2684 :   return powis(p, v);
    1250             : }
    1251             : 
    1252             : /* x,y in Z[i], at least one of which is t_COMPLEX */
    1253             : static GEN
    1254         385 : gauss_gcd(GEN x, GEN y)
    1255             : {
    1256         385 :   pari_sp av = avma;
    1257             :   GEN dx, dy;
    1258         385 :   dx = denom_i(x); x = gmul(x, dx);
    1259         385 :   dy = denom_i(y); y = gmul(y, dy);
    1260        1197 :   while (!gequal0(y))
    1261             :   {
    1262         427 :     GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
    1263         427 :     x = y; y = z;
    1264             :   }
    1265         385 :   x = gauss_normal(x);
    1266         385 :   if (typ(x) == t_COMPLEX)
    1267             :   {
    1268         196 :     if      (gequal0(gel(x,2))) x = gel(x,1);
    1269         196 :     else if (gequal0(gel(x,1))) x = gel(x,2);
    1270             :   }
    1271         385 :   return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
    1272             : }
    1273             : 
    1274             : static int
    1275        2520 : c_is_rational(GEN x)
    1276        2520 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
    1277             : static GEN
    1278        1288 : c_zero_gcd(GEN c)
    1279             : {
    1280        1288 :   GEN x = gel(c,1), y = gel(c,2);
    1281        1288 :   long tx = typ(x), ty = typ(y);
    1282        1288 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
    1283          42 :   if (tx == t_PADIC || tx == t_INTMOD
    1284          35 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
    1285          35 :   return gauss_gcd(c, gen_0);
    1286             : }
    1287             : 
    1288             : /* gcd(x, 0) */
    1289             : static GEN
    1290     8227087 : zero_gcd(GEN x)
    1291             : {
    1292             :   pari_sp av;
    1293     8227087 :   switch(typ(x))
    1294             :   {
    1295       24436 :     case t_INT: return absi(x);
    1296        1302 :     case t_FRAC: return absfrac(x);
    1297        1288 :     case t_COMPLEX: return c_zero_gcd(x);
    1298        7000 :     case t_REAL: return gen_1;
    1299         742 :     case t_PADIC: return powis(gel(x,2), valp(x));
    1300         245 :     case t_SER: return pol_xnall(valp(x), varn(x));
    1301             :     case t_POLMOD: {
    1302        8338 :       GEN d = gel(x,2);
    1303        8338 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
    1304          28 :       return isinexact(d)? zero_gcd(d): gcopy(d);
    1305             :     }
    1306             :     case t_POL:
    1307     7921120 :       if (!isinexact(x)) break;
    1308          14 :       av = avma;
    1309          14 :       return gerepileupto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
    1310             : 
    1311             :     case t_RFRAC:
    1312      217144 :       if (!isinexact(x)) break;
    1313           0 :       av = avma;
    1314           0 :       return gerepileupto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
    1315             :   }
    1316     8183722 :   return gcopy(x);
    1317             : }
    1318             : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
    1319             : static GEN
    1320     9335490 : zero_gcd2(GEN y, GEN z)
    1321             : {
    1322             :   pari_sp av;
    1323     9335490 :   switch(typ(z))
    1324             :   {
    1325     8204640 :     case t_INT: return zero_gcd(y);
    1326             :     case t_INTMOD:
    1327     1129807 :       av = avma;
    1328     1129807 :       return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
    1329             :     case t_FFELT:
    1330        1043 :       av = avma;
    1331        1043 :       return gerepileupto(av, gmul(y, FF_1(z)));
    1332             :     default:
    1333           0 :       pari_err_TYPE("zero_gcd", z);
    1334             :   }
    1335           0 :   return NULL;
    1336             : }
    1337             : static GEN
    1338     2256947 : cont_gcd_pol_i(GEN x, GEN y) { return scalarpol(ggcd(content(x),y), varn(x));}
    1339             : /* tx = t_POL, y considered as constant */
    1340             : static GEN
    1341     2256940 : cont_gcd_pol(GEN x, GEN y)
    1342     2256940 : { pari_sp av = avma; return gerepileupto(av, cont_gcd_pol_i(x,y)); }
    1343             : /* tx = t_RFRAC, y considered as constant */
    1344             : static GEN
    1345      864904 : cont_gcd_rfrac(GEN x, GEN y)
    1346             : {
    1347      864904 :   pari_sp av = avma;
    1348      864904 :   GEN cx; x = primitive_part(x, &cx);
    1349             :   /* e.g. Mod(1,2) / (2*y+1) => primitive_part = Mod(1,2)*y^0 */
    1350      864904 :   if (typ(x) != t_RFRAC) x = cont_gcd_pol_i(x, y);
    1351      864897 :   else x = gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2));
    1352      864904 :   return gerepileupto(av, x);
    1353             : }
    1354             : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
    1355             : static GEN
    1356        6581 : cont_gcd_gen(GEN x, GEN y)
    1357             : {
    1358        6581 :   pari_sp av = avma;
    1359        6581 :   return gerepileupto(av, ggcd(content(x),y));
    1360             : }
    1361             : /* !is_const(tx), y considered as constant */
    1362             : static GEN
    1363     3071637 : cont_gcd(GEN x, long tx, GEN y)
    1364             : {
    1365     3071637 :   switch(tx)
    1366             :   {
    1367      808130 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1368     2256926 :     case t_POL: return cont_gcd_pol(x,y);
    1369        6581 :     default: return cont_gcd_gen(x,y);
    1370             :   }
    1371             : }
    1372             : static GEN
    1373    10076025 : gcdiq(GEN x, GEN y)
    1374             : {
    1375             :   GEN z;
    1376    10076025 :   if (!signe(x)) return Q_abs(y);
    1377     2532965 :   z = cgetg(3,t_FRAC);
    1378     2532965 :   gel(z,1) = gcdii(x,gel(y,1));
    1379     2532965 :   gel(z,2) = icopy(gel(y,2));
    1380     2532965 :   return z;
    1381             : }
    1382             : static GEN
    1383    11226837 : gcdqq(GEN x, GEN y)
    1384             : {
    1385    11226837 :   GEN z = cgetg(3,t_FRAC);
    1386    11226837 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1387    11226837 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1388    11226837 :   return z;
    1389             : }
    1390             : /* assume x,y t_INT or t_FRAC */
    1391             : GEN
    1392   104100235 : Q_gcd(GEN x, GEN y)
    1393             : {
    1394   104100235 :   long tx = typ(x), ty = typ(y);
    1395   104100235 :   if (tx == t_INT)
    1396    83650503 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1397             :   else
    1398    20449732 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1399             : }
    1400             : 
    1401             : GEN
    1402    26738760 : ggcd(GEN x, GEN y)
    1403             : {
    1404    26738760 :   long vx, vy, tx = typ(x), ty = typ(y);
    1405             :   pari_sp av, tetpil;
    1406             :   GEN p1,z;
    1407             : 
    1408    53477520 :   if (is_noncalc_t(tx) || is_matvec_t(tx) ||
    1409    53477520 :       is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
    1410    26738760 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1411             :   /* tx <= ty */
    1412    26738760 :   z = gisexactzero(x); if (z) return zero_gcd2(y,z);
    1413    23671161 :   z = gisexactzero(y); if (z) return zero_gcd2(x,z);
    1414    17403270 :   if (is_const_t(tx))
    1415             :   {
    1416    10466964 :     if (ty == tx) switch(tx)
    1417             :     {
    1418             :       case t_INT:
    1419     3210353 :         return gcdii(x,y);
    1420             : 
    1421     4147353 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1422     4147353 :         if (equalii(gel(x,1),gel(y,1)))
    1423     4147346 :           gel(z,1) = icopy(gel(x,1));
    1424             :         else
    1425           7 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1426     4147353 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1427             :         else
    1428             :         {
    1429     4147353 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1430     4147353 :           if (!equali1(p1))
    1431             :           {
    1432           7 :             p1 = gcdii(p1,gel(y,2));
    1433           7 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1434           7 :             else p1 = gerepileuptoint(av, p1);
    1435             :           }
    1436     4147353 :           gel(z,2) = p1;
    1437             :         }
    1438     4147353 :         return z;
    1439             : 
    1440             :       case t_FRAC:
    1441       32697 :         return gcdqq(x,y);
    1442             : 
    1443             :       case t_FFELT:
    1444        8841 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1445        8841 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1446             : 
    1447             :       case t_COMPLEX:
    1448          21 :         if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
    1449           7 :         return triv_cont_gcd(y,x);
    1450             : 
    1451             :       case t_PADIC:
    1452          14 :         if (!equalii(gel(x,2),gel(y,2))) return gen_1;
    1453           7 :         return powis(gel(y,2), minss(valp(x), valp(y)));
    1454             : 
    1455             :       case t_QUAD:
    1456         133 :         av=avma; p1=gdiv(x,y);
    1457         133 :         if (gequal0(gel(p1,3)))
    1458             :         {
    1459          14 :           p1=gel(p1,2);
    1460          14 :           if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
    1461           7 :           tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
    1462             :         }
    1463         119 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
    1464         112 :         p1 = ginv(p1); avma=av;
    1465         112 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1466         105 :         return triv_cont_gcd(y,x);
    1467             : 
    1468           0 :       default: return gen_1; /* t_REAL */
    1469             :     }
    1470     3067552 :     if (is_const_t(ty)) switch(tx)
    1471             :     {
    1472             :       case t_INT:
    1473       23005 :         switch(ty)
    1474             :         {
    1475        2737 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1476        2737 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1477        2737 :             p1 = gcdii(gel(y,1),gel(y,2));
    1478        2737 :             if (!equali1(p1)) {
    1479          14 :               p1 = gcdii(x,p1);
    1480          14 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1481             :               else
    1482          14 :                 p1 = gerepileuptoint(av, p1);
    1483             :             }
    1484        2737 :             gel(z,2) = p1; return z;
    1485             : 
    1486        6825 :           case t_REAL: return gen_1;
    1487             : 
    1488             :           case t_FRAC:
    1489        8022 :             return gcdiq(x,y);
    1490             : 
    1491             :           case t_COMPLEX:
    1492        2401 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1493        2072 :             return triv_cont_gcd(y,x);
    1494             : 
    1495             :           case t_FFELT:
    1496         224 :             if (!FF_equal0(y)) return FF_1(y);
    1497           0 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1498             : 
    1499             :           case t_PADIC:
    1500        2670 :             return padic_gcd(x,y);
    1501             : 
    1502             :           case t_QUAD:
    1503         126 :             return triv_cont_gcd(y,x);
    1504             :           default:
    1505           0 :             pari_err_TYPE2("gcd",x,y);
    1506             :         }
    1507             : 
    1508             :       case t_REAL:
    1509          14 :         switch(ty)
    1510             :         {
    1511             :           case t_INTMOD:
    1512             :           case t_FFELT:
    1513          14 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1514           0 :           default: return gen_1;
    1515             :         }
    1516             : 
    1517             :       case t_INTMOD:
    1518          49 :         switch(ty)
    1519             :         {
    1520             :           case t_FRAC:
    1521          14 :             av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
    1522          14 :             if (!equali1(p1)) pari_err_OP("gcd",x,y);
    1523           7 :             return ggcd(gel(y,1), x);
    1524             : 
    1525             :           case t_FFELT:
    1526             :           {
    1527          14 :             GEN p = gel(y,4);
    1528          14 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1529           7 :             if (!FF_equal0(y)) return FF_1(y);
    1530           0 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1531             :           }
    1532             : 
    1533             :           case t_COMPLEX: case t_QUAD:
    1534          14 :             return triv_cont_gcd(y,x);
    1535             : 
    1536             :           case t_PADIC:
    1537           7 :             return padic_gcd(x,y);
    1538             : 
    1539           0 :           default: pari_err_TYPE2("gcd",x,y);
    1540             :         }
    1541             : 
    1542             :       case t_FRAC:
    1543         210 :         switch(ty)
    1544             :         {
    1545             :           case t_COMPLEX:
    1546          84 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1547             :           case t_QUAD:
    1548         154 :             return triv_cont_gcd(y,x);
    1549             :           case t_FFELT:
    1550             :           {
    1551          42 :             GEN p = gel(y,4);
    1552          42 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1553          21 :             if (!FF_equal0(y)) return FF_1(y);
    1554           0 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1555             :           }
    1556             : 
    1557             :           case t_PADIC:
    1558           7 :             return padic_gcd(x,y);
    1559             : 
    1560           0 :           default: pari_err_TYPE2("gcd",x,y);
    1561             :         }
    1562             :       case t_FFELT:
    1563          70 :         switch(ty)
    1564             :         {
    1565             :           case t_PADIC:
    1566             :           {
    1567          42 :             GEN p = gel(y,2);
    1568          42 :             long v = valp(y);
    1569          42 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1570          14 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1571             :           }
    1572          28 :           default: pari_err_TYPE2("gcd",x,y);
    1573             :         }
    1574             : 
    1575             :       case t_COMPLEX:
    1576          14 :         switch(ty)
    1577             :         {
    1578             :           case t_PADIC:
    1579          14 :           case t_QUAD: return triv_cont_gcd(x,y);
    1580           0 :           default: pari_err_TYPE2("gcd",x,y);
    1581             :         }
    1582             : 
    1583             :       case t_PADIC:
    1584           7 :         switch(ty)
    1585             :         {
    1586           7 :           case t_QUAD: return triv_cont_gcd(y,x);
    1587           0 :           default: pari_err_TYPE2("gcd",x,y);
    1588             :         }
    1589             : 
    1590           0 :       default: return gen_1; /* tx = t_REAL */
    1591             :     }
    1592     3044183 :     return cont_gcd(y,ty, x);
    1593             :   }
    1594             : 
    1595     6936306 :   if (tx == t_POLMOD)
    1596             :   {
    1597       19372 :     if (ty == t_POLMOD)
    1598             :     {
    1599       19316 :       GEN T = gel(x,1);
    1600       19316 :       z = cgetg(3,t_POLMOD);
    1601       19316 :       T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
    1602       19316 :       gel(z,1) = T;
    1603       19316 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1604             :       else
    1605             :       {
    1606             :         GEN X, Y, d;
    1607       19316 :         av = avma; X = gel(x,2); Y = gel(y,2);
    1608       19316 :         d = ggcd(content(X), content(Y));
    1609       19316 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1610       19316 :         p1 = ggcd(T, X);
    1611       19316 :         gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
    1612             :       }
    1613       19316 :       return z;
    1614             :     }
    1615          56 :     vx = varn(gel(x,1));
    1616          56 :     switch(ty)
    1617             :     {
    1618             :       case t_POL:
    1619          28 :         vy = varn(y);
    1620          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1621          14 :         z = cgetg(3,t_POLMOD);
    1622          14 :         gel(z,1) = RgX_copy(gel(x,1));
    1623          14 :         av = avma; p1 = ggcd(gel(x,1),gel(x,2));
    1624          14 :         gel(z,2) = gerepileupto(av, ggcd(p1,y));
    1625          14 :         return z;
    1626             : 
    1627             :       case t_RFRAC:
    1628          28 :         vy = varn(gel(y,2));
    1629          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1630          28 :         av = avma;
    1631          28 :         p1 = ggcd(gel(x,1),gel(y,2));
    1632          28 :         if (degpol(p1)) pari_err_OP("gcd",x,y);
    1633          21 :         avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1634             :     }
    1635             :   }
    1636             : 
    1637     6916934 :   vx = gvar(x);
    1638     6916934 :   vy = gvar(y);
    1639     6916934 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1640     6904929 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1641             : 
    1642             :   /* vx = vy: same main variable */
    1643     6889480 :   switch(tx)
    1644             :   {
    1645             :     case t_POL:
    1646     6728574 :       switch(ty)
    1647             :       {
    1648     6671779 :         case t_POL: return RgX_gcd(x,y);
    1649             :         case t_SER:
    1650          21 :           z = ggcd(content(x), content(y));
    1651          21 :           return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
    1652       56774 :         case t_RFRAC: return cont_gcd_rfrac(y, x);
    1653             :       }
    1654           0 :       break;
    1655             : 
    1656             :     case t_SER:
    1657          14 :       z = ggcd(content(x), content(y));
    1658          14 :       switch(ty)
    1659             :       {
    1660           7 :         case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
    1661           7 :         case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
    1662             :       }
    1663           0 :       break;
    1664             : 
    1665             :     case t_RFRAC:
    1666             :     {
    1667      160892 :       GEN xd = gel(x,2), yd = gel(y,2);
    1668      160892 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1669      160892 :       z = cgetg(3,t_RFRAC); av = avma;
    1670      160892 :       gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1671      160892 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1672             :     }
    1673             :   }
    1674           0 :   pari_err_TYPE2("gcd",x,y);
    1675             :   return NULL; /* LCOV_EXCL_LINE */
    1676             : }
    1677             : GEN
    1678        3959 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
    1679             : 
    1680             : static GEN
    1681         329 : fix_lcm(GEN x)
    1682             : {
    1683             :   GEN t;
    1684         329 :   switch(typ(x))
    1685             :   {
    1686         231 :     case t_INT: if (signe(x)<0) x = negi(x);
    1687         231 :       break;
    1688             :     case t_POL:
    1689          98 :       if (lg(x) <= 2) break;
    1690          98 :       t = leading_coeff(x);
    1691          98 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1692             :   }
    1693         329 :   return x;
    1694             : }
    1695             : GEN
    1696        3101 : glcm0(GEN x, GEN y)
    1697             : {
    1698        3101 :   if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
    1699        2821 :   return glcm(x,y);
    1700             : }
    1701             : GEN
    1702        3619 : glcm(GEN x, GEN y)
    1703             : {
    1704             :   pari_sp av;
    1705             :   GEN z;
    1706        3619 :   if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
    1707          63 :   av = avma; z = ggcd(x,y);
    1708          63 :   if (!gequal1(z))
    1709             :   {
    1710          63 :     if (gequal0(z)) { avma = av; return gmul(x,y); }
    1711          49 :     y = gdiv(y,z);
    1712             :   }
    1713          49 :   return gerepileupto(av, fix_lcm(gmul(x,y)));
    1714             : }
    1715             : 
    1716             : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
    1717             : static int
    1718      603057 : pol_approx0(GEN r, GEN x, int exact)
    1719             : {
    1720             :   long i, lx,lr;
    1721      603057 :   if (exact) return gequal0(r);
    1722           0 :   lx = lg(x);
    1723           0 :   lr = lg(r); if (lr < lx) lx = lr;
    1724           0 :   for (i=2; i<lx; i++)
    1725           0 :     if (!approx_0(gel(r,i), gel(x,i))) return 0;
    1726           0 :   return 1;
    1727             : }
    1728             : 
    1729             : GEN
    1730      173593 : RgX_gcd_simple(GEN x, GEN y)
    1731             : {
    1732      173593 :   pari_sp av1, av = avma;
    1733      173593 :   GEN r, yorig = y;
    1734      173593 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1735             : 
    1736             :   for(;;)
    1737             :   {
    1738     1032521 :     av1 = avma; r = RgX_rem(x,y);
    1739      603057 :     if (pol_approx0(r, x, exact))
    1740             :     {
    1741      173593 :       avma = av1;
    1742      173593 :       if (y == yorig) return RgX_copy(y);
    1743      121198 :       y = normalizepol_approx(y, lg(y));
    1744      121198 :       if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
    1745        8561 :       return gerepileupto(av,y);
    1746             :     }
    1747      429464 :     x = y; y = r;
    1748      429464 :     if (gc_needed(av,1)) {
    1749           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1750           0 :       gerepileall(av,2, &x,&y);
    1751             :     }
    1752             :   }
    1753             : }
    1754             : GEN
    1755           0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
    1756             : {
    1757           0 :   pari_sp av = avma;
    1758             :   GEN q, r, d, d1, u, v, v1;
    1759           0 :   int exact = !(isinexactreal(a) || isinexactreal(b));
    1760             : 
    1761           0 :   d = a; d1 = b; v = gen_0; v1 = gen_1;
    1762             :   for(;;)
    1763             :   {
    1764           0 :     if (pol_approx0(d1, a, exact)) break;
    1765           0 :     q = poldivrem(d,d1, &r);
    1766           0 :     v = gsub(v, gmul(q,v1));
    1767           0 :     u=v; v=v1; v1=u;
    1768           0 :     u=r; d=d1; d1=u;
    1769             :   }
    1770           0 :   u = gsub(d, gmul(b,v));
    1771           0 :   u = RgX_div(u,a);
    1772             : 
    1773           0 :   gerepileall(av, 3, &u,&v,&d);
    1774           0 :   *pu = u;
    1775           0 :   *pv = v; return d;
    1776             : }
    1777             : /*******************************************************************/
    1778             : /*                                                                 */
    1779             : /*                    CONTENT / PRIMITIVE PART                     */
    1780             : /*                                                                 */
    1781             : /*******************************************************************/
    1782             : 
    1783             : GEN
    1784    78668391 : content(GEN x)
    1785             : {
    1786    78668391 :   long lx, i, t, tx = typ(x);
    1787    78668391 :   pari_sp av = avma;
    1788             :   GEN c;
    1789             : 
    1790    78668391 :   if (is_scalar_t(tx)) return zero_gcd(x);
    1791    78655163 :   switch(tx)
    1792             :   {
    1793             :     case t_RFRAC:
    1794             :     {
    1795      864939 :       GEN n = gel(x,1), d = gel(x,2);
    1796             :       /* -- varncmp(vn, vd) < 0 can't happen
    1797             :        * -- if n is POLMOD, its main variable (in the sense of gvar2)
    1798             :        *    has lower priority than denominator */
    1799      864939 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    1800      824601 :         n = isinexact(n)? zero_gcd(n): gcopy(n);
    1801             :       else
    1802       40338 :         n = content(n);
    1803      864939 :       return gerepileupto(av, gdiv(n, content(d)));
    1804             :     }
    1805             : 
    1806             :     case t_VEC: case t_COL:
    1807     7663719 :       lx = lg(x); if (lx==1) return gen_0;
    1808     7663712 :       break;
    1809             : 
    1810             :     case t_MAT:
    1811             :     {
    1812             :       long hx, j;
    1813         238 :       lx = lg(x);
    1814         238 :       if (lx == 1) return gen_0;
    1815         231 :       hx = lgcols(x);
    1816         231 :       if (hx == 1) return gen_0;
    1817         224 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    1818         224 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    1819         224 :       c = content(gel(x,1));
    1820         448 :       for (j=2; j<lx; j++)
    1821         224 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    1822         224 :       if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
    1823         224 :       return gerepileupto(av,c);
    1824             :     }
    1825             : 
    1826             :     case t_POL: case t_SER:
    1827    70126232 :       lx = lg(x); if (lx == 2) return gen_0;
    1828    70126218 :       break;
    1829          21 :     case t_VECSMALL: return utoi(zv_content(x));
    1830             :     case t_QFR: case t_QFI:
    1831          14 :       lx = 4; break;
    1832             : 
    1833           0 :     default: pari_err_TYPE("content",x);
    1834             :       return NULL; /* LCOV_EXCL_LINE */
    1835             :   }
    1836   230330864 :   for (i=lontyp[tx]; i<lx; i++)
    1837   165745807 :     if (typ(gel(x,i)) != t_INT) break;
    1838    77789944 :   lx--; c = gel(x,lx);
    1839    77789944 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    1840    77789944 :   if (i > lx)
    1841             :   { /* integer coeffs */
    1842   130877959 :     while (lx-- > lontyp[tx])
    1843             :     {
    1844    65450770 :       c = gcdii(c, gel(x,lx));
    1845    65450770 :       if (equali1(c)) { avma=av; return gen_1; }
    1846             :     }
    1847             :   }
    1848             :   else
    1849             :   {
    1850    13204887 :     if (isinexact(c)) c = zero_gcd(c);
    1851    49790989 :     while (lx-- > lontyp[tx])
    1852             :     {
    1853    23381215 :       GEN d = gel(x,lx);
    1854    23381215 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    1855    23381215 :       c = ggcd(c, d);
    1856             :     }
    1857    13204887 :     if (isinexact(c)) { avma=av; return gen_1; }
    1858             :   }
    1859    14047019 :   switch(typ(c))
    1860             :   {
    1861             :     case t_INT:
    1862      853773 :       if (signe(c) < 0) c = negi(c);
    1863      853773 :       break;
    1864             :     case t_VEC: case t_COL: case t_MAT:
    1865           0 :       pari_err_TYPE("content",x);
    1866             :   }
    1867             : 
    1868    14047019 :   return av==avma? gcopy(c): gerepileupto(av,c);
    1869             : }
    1870             : 
    1871             : GEN
    1872     1023221 : primitive_part(GEN x, GEN *ptc)
    1873             : {
    1874     1023221 :   pari_sp av = avma;
    1875     1023221 :   GEN c = content(x);
    1876     1023221 :   if (gequal1(c)) { avma = av; c = NULL; }
    1877       44223 :   else if (!gequal0(c)) x = gdiv(x,c);
    1878     1023221 :   if (ptc) *ptc = c;
    1879     1023221 :   return x;
    1880             : }
    1881             : GEN
    1882        2842 : primpart(GEN x) { return primitive_part(x, NULL); }
    1883             : 
    1884             : static GEN
    1885    23792249 : Q_content_v(GEN x, long i, long l)
    1886             : {
    1887    23792249 :   pari_sp av = avma;
    1888    23792249 :   GEN d = Q_content_safe(gel(x,i));
    1889    23792249 :   if (!d) return NULL;
    1890   124531591 :   for (i++; i < l; i++)
    1891             :   {
    1892   100739545 :     GEN c = Q_content_safe(gel(x,i));
    1893   100739545 :     if (!c) return NULL;
    1894   100739545 :     d = Q_gcd(d, c);
    1895             :   }
    1896    23792046 :   return gerepileupto(av, d);
    1897             : }
    1898             : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1899             :  * of Q(x2,...,xn)[x1] */
    1900             : GEN
    1901   139343631 : Q_content_safe(GEN x)
    1902             : {
    1903             :   long l;
    1904   139343631 :   switch(typ(x))
    1905             :   {
    1906   100863225 :     case t_INT:  return absi(x);
    1907    14384748 :     case t_FRAC: return absfrac(x);
    1908             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    1909    12284744 :       l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
    1910             :     case t_POL:
    1911    11785126 :       l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
    1912       25557 :     case t_POLMOD: return Q_content_safe(gel(x,2));
    1913             :     case t_RFRAC:
    1914             :     {
    1915             :       GEN a, b;
    1916          21 :       a = Q_content(gel(x,1)); if (!a) return NULL;
    1917          21 :       b = Q_content(gel(x,2)); if (!b) return NULL;
    1918          21 :       return gdiv(a, b);
    1919             :     }
    1920             :   }
    1921         210 :   return NULL;
    1922             : }
    1923             : GEN
    1924       92242 : Q_content(GEN x)
    1925             : {
    1926       92242 :   GEN c = Q_content_safe(x);
    1927       92242 :   if (!c) pari_err_TYPE("Q_content",x);
    1928       92242 :   return c;
    1929             : }
    1930             : 
    1931             : GEN
    1932       42094 : ZX_content(GEN x)
    1933             : {
    1934       42094 :   long i, l = lg(x);
    1935             :   GEN d;
    1936             :   pari_sp av;
    1937             : 
    1938       42094 :   if (l == 2) return gen_0;
    1939       42094 :   d = gel(x,2);
    1940       42094 :   if (l == 3) return absi(d);
    1941       35855 :   av = avma;
    1942       35855 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    1943       35855 :   if (signe(d) < 0) d = negi(d);
    1944       35855 :   return gerepileuptoint(av, d);
    1945             : }
    1946             : 
    1947             : static GEN
    1948       86499 : Z_content_v(GEN x, long i, long l)
    1949             : {
    1950       86499 :   pari_sp av = avma;
    1951       86499 :   GEN d = Z_content(gel(x,i));
    1952       86499 :   if (!d) return NULL;
    1953      381927 :   for (i++; i<l; i++)
    1954             :   {
    1955      298711 :     GEN c = Z_content(gel(x,i));
    1956      298711 :     if (!c) return NULL;
    1957      298585 :     d = gcdii(d, c); if (equali1(d)) return NULL;
    1958      298543 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    1959             :   }
    1960       83216 :   return gerepileuptoint(av, d);
    1961             : }
    1962             : /* return NULL for 1 */
    1963             : GEN
    1964      386582 : Z_content(GEN x)
    1965             : {
    1966             :   long l;
    1967      386582 :   switch(typ(x))
    1968             :   {
    1969             :     case t_INT:
    1970      283430 :       if (is_pm1(x)) return NULL;
    1971      282646 :       return absi(x);
    1972             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    1973        8967 :       l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
    1974             :     case t_POL:
    1975       94185 :       l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
    1976           0 :     case t_POLMOD: return Z_content(gel(x,2));
    1977             :   }
    1978           0 :   pari_err_TYPE("Z_content", x);
    1979             :   return NULL; /* LCOV_EXCL_LINE */
    1980             : }
    1981             : 
    1982             : static GEN
    1983    11275161 : Q_denom_v(GEN x, long i, long l)
    1984             : {
    1985    11275161 :   pari_sp av = avma;
    1986    11275161 :   GEN d = Q_denom_safe(gel(x,i));
    1987    11275161 :   if (!d) return NULL;
    1988    49219175 :   for (i++; i<l; i++)
    1989             :   {
    1990    37944014 :     GEN D = Q_denom_safe(gel(x,i));
    1991    37944014 :     if (!D) return NULL;
    1992    37944014 :     if (D != gen_1) d = lcmii(d, D);
    1993    37944014 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    1994             :   }
    1995    11275161 :   return gerepileuptoint(av, d);
    1996             : }
    1997             : /* NOT MEMORY CLEAN (because of t_FRAC).
    1998             :  * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1999             :  * of Q(x2,...,xn)[x1] */
    2000             : GEN
    2001    73021991 : Q_denom_safe(GEN x)
    2002             : {
    2003             :   long l;
    2004    73021991 :   switch(typ(x))
    2005             :   {
    2006    48128120 :     case t_INT: return gen_1;
    2007          28 :     case t_PADIC: l = valp(x); return l < 0? powiu(gel(x,2), -l): gen_1;
    2008    13491100 :     case t_FRAC: return gel(x,2);
    2009          14 :     case t_QUAD: return Q_denom_v(x, 2, 4);
    2010             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2011     8536562 :       l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
    2012             :     case t_POL: case t_SER:
    2013     2789685 :       l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
    2014       76356 :     case t_POLMOD: return Q_denom(gel(x,2));
    2015             :     case t_RFRAC:
    2016             :     {
    2017             :       GEN a, b;
    2018          56 :       a = Q_content(gel(x,1)); if (!a) return NULL;
    2019          56 :       b = Q_content(gel(x,2)); if (!b) return NULL;
    2020          56 :       return Q_denom(gdiv(a, b));
    2021             :     }
    2022             :   }
    2023          70 :   return NULL;
    2024             : }
    2025             : GEN
    2026    15268246 : Q_denom(GEN x)
    2027             : {
    2028    15268246 :   GEN d = Q_denom_safe(x);
    2029    15268246 :   if (!d) pari_err_TYPE("Q_denom",x);
    2030    15268246 :   return d;
    2031             : }
    2032             : 
    2033             : GEN
    2034     8534430 : Q_remove_denom(GEN x, GEN *ptd)
    2035             : {
    2036     8534430 :   GEN d = Q_denom_safe(x);
    2037     8534430 :   if (d) { if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d); }
    2038     8534430 :   if (ptd) *ptd = d;
    2039     8534430 :   return x;
    2040             : }
    2041             : 
    2042             : /* return y = x * d, assuming x rational, and d,y integral */
    2043             : GEN
    2044    47348544 : Q_muli_to_int(GEN x, GEN d)
    2045             : {
    2046             :   long i, l;
    2047             :   GEN y, xn, xd;
    2048             :   pari_sp av;
    2049             : 
    2050    47348544 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    2051    47348544 :   switch (typ(x))
    2052             :   {
    2053             :     case t_INT:
    2054    19986333 :       return mulii(x,d);
    2055             : 
    2056             :     case t_FRAC:
    2057    20076571 :       xn = gel(x,1);
    2058    20076571 :       xd = gel(x,2); av = avma;
    2059    20076571 :       y = mulii(xn, diviiexact(d, xd));
    2060    20076571 :       return gerepileuptoint(av, y);
    2061             :     case t_COMPLEX:
    2062           7 :       y = cgetg(3,t_COMPLEX);
    2063           7 :       gel(y,1) = Q_muli_to_int(gel(x,1),d);
    2064           7 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2065           7 :       return y;
    2066             :     case t_PADIC:
    2067          14 :       y = gcopy(x); if (!isint1(d)) setvalp(y, 0);
    2068          14 :       return y;
    2069             :     case t_QUAD:
    2070           7 :       y = cgetg(4,t_QUAD);
    2071           7 :       gel(y,1) = ZX_copy(gel(x,1));
    2072           7 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2073           7 :       gel(y,3) = Q_muli_to_int(gel(x,3),d); return y;
    2074             : 
    2075             :     case t_VEC: case t_COL: case t_MAT:
    2076     5112901 :       y = cgetg_copy(x, &l);
    2077     5112901 :       for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2078     5112901 :       return y;
    2079             : 
    2080             :     case t_POL: case t_SER:
    2081     2128716 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2082     2128716 :       for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2083     2128716 :       return y;
    2084             : 
    2085             :     case t_POLMOD:
    2086       43974 :       retmkpolmod(Q_muli_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2087             :     case t_RFRAC:
    2088          21 :       return gmul(x, d);
    2089             :   }
    2090           0 :   pari_err_TYPE("Q_muli_to_int",x);
    2091             :   return NULL; /* LCOV_EXCL_LINE */
    2092             : }
    2093             : 
    2094             : static void
    2095     1848126 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
    2096             : {
    2097             :   long e;
    2098     1848126 :   switch(typ(c))
    2099             :   {
    2100             :     case t_REAL:
    2101     1238439 :       *exact = 0;
    2102     1238439 :       if (!signe(c)) return;
    2103     1235892 :       e = expo(c) - bit_prec(c);
    2104     1235892 :       break;
    2105             :     case t_INT:
    2106      609589 :       if (!signe(c)) return;
    2107      250116 :       e = expi(c) + 32;
    2108      250116 :       break;
    2109             :     case t_FRAC:
    2110          70 :       e = expi(gel(c,1)) - expi(gel(c,2)) + 32;
    2111          70 :       if (exact) *D = lcmii(*D, gel(c,2));
    2112          70 :       break;
    2113             :     default:
    2114          28 :       pari_err_TYPE("rescale_to_int",c);
    2115             :       return; /* LCOV_EXCL_LINE */
    2116             :   }
    2117     1486078 :   if (e < *emin) *emin = e;
    2118             : }
    2119             : GEN
    2120      223439 : RgM_rescale_to_int(GEN x)
    2121             : {
    2122      223439 :   long lx = lg(x), i,j, hx, emin;
    2123             :   GEN D;
    2124             :   int exact;
    2125             : 
    2126      223439 :   if (lx == 1) return cgetg(1,t_MAT);
    2127      223439 :   hx = lgcols(x);
    2128      223439 :   exact = 1;
    2129      223439 :   emin = HIGHEXPOBIT;
    2130      223439 :   D = gen_1;
    2131      810689 :   for (j = 1; j < lx; j++)
    2132      587278 :     for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
    2133      223411 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2134      223355 :   return grndtoi(gmul2n(x, -emin), &i);
    2135             : }
    2136             : GEN
    2137         238 : RgX_rescale_to_int(GEN x)
    2138             : {
    2139         238 :   long lx = lg(x), i, emin;
    2140             :   GEN D;
    2141             :   int exact;
    2142         238 :   if (lx == 2) return gcopy(x); /* rare */
    2143         238 :   exact = 1;
    2144         238 :   emin = HIGHEXPOBIT;
    2145         238 :   D = gen_1;
    2146         238 :   for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
    2147         238 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2148         238 :   return grndtoi(gmul2n(x, -emin), &i);
    2149             : }
    2150             : 
    2151             : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    2152             : static GEN
    2153     8302225 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    2154             : {
    2155             :   long i, l;
    2156             :   GEN y, xn, xd;
    2157             :   pari_sp av;
    2158             : 
    2159     8302225 :   switch(typ(x))
    2160             :   {
    2161             :     case t_INT:
    2162     2503939 :       av = avma; y = diviiexact(x,d);
    2163     2503939 :       return gerepileuptoint(av, mulii(y,n));
    2164             : 
    2165             :     case t_FRAC:
    2166     3998152 :       xn = gel(x,1);
    2167     3998152 :       xd = gel(x,2); av = avma;
    2168     3998152 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    2169     3998152 :       return gerepileuptoint(av, y);
    2170             : 
    2171             :     case t_VEC: case t_COL: case t_MAT:
    2172      423992 :       y = cgetg_copy(x, &l);
    2173      423992 :       for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2174      423992 :       return y;
    2175             : 
    2176             :     case t_POL:
    2177     1376142 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2178     1376142 :       for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2179     1376142 :       return y;
    2180             : 
    2181             :     case t_POLMOD:
    2182           0 :       retmkpolmod(Q_divmuli_to_int(gel(x,2), d,n), RgX_copy(gel(x,1)));
    2183             :   }
    2184           0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    2185             :   return NULL; /* LCOV_EXCL_LINE */
    2186             : }
    2187             : 
    2188             : /* return x / d. x: rational; d,result: integral. */
    2189             : static GEN
    2190    16709146 : Q_divi_to_int(GEN x, GEN d)
    2191             : {
    2192             :   long i, l;
    2193             :   GEN y;
    2194             : 
    2195    16709146 :   switch(typ(x))
    2196             :   {
    2197             :     case t_INT:
    2198    13647996 :       return diviiexact(x,d);
    2199             : 
    2200             :     case t_VEC: case t_COL: case t_MAT:
    2201     1200293 :       y = cgetg_copy(x, &l);
    2202     1200293 :       for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2203     1200293 :       return y;
    2204             : 
    2205             :     case t_POL:
    2206     1857231 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2207     1857231 :       for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2208     1857231 :       return y;
    2209             : 
    2210             :     case t_POLMOD:
    2211        3626 :       retmkpolmod(Q_divi_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2212             :   }
    2213           0 :   pari_err_TYPE("Q_divi_to_int",x);
    2214             :   return NULL; /* LCOV_EXCL_LINE */
    2215             : }
    2216             : /* c t_FRAC */
    2217             : static GEN
    2218     3190433 : Q_divq_to_int(GEN x, GEN c)
    2219             : {
    2220     3190433 :   GEN n = gel(c,1), d = gel(c,2);
    2221     3190433 :   if (is_pm1(n)) {
    2222     1762798 :     GEN y = Q_muli_to_int(x,d);
    2223     1762798 :     if (signe(n) < 0) y = gneg(y);
    2224     1762798 :     return y;
    2225             :   }
    2226     1427635 :   return Q_divmuli_to_int(x, n,d);
    2227             : }
    2228             : 
    2229             : /* return y = x / c, assuming x,c rational, and y integral */
    2230             : GEN
    2231        9019 : Q_div_to_int(GEN x, GEN c)
    2232             : {
    2233        9019 :   switch(typ(c))
    2234             :   {
    2235        7787 :     case t_INT:  return Q_divi_to_int(x, c);
    2236        1232 :     case t_FRAC: return Q_divq_to_int(x, c);
    2237             :   }
    2238           0 :   pari_err_TYPE("Q_div_to_int",c);
    2239             :   return NULL; /* LCOV_EXCL_LINE */
    2240             : }
    2241             : /* return y = x * c, assuming x,c rational, and y integral */
    2242             : GEN
    2243           0 : Q_mul_to_int(GEN x, GEN c)
    2244             : {
    2245             :   GEN d, n;
    2246           0 :   switch(typ(c))
    2247             :   {
    2248           0 :     case t_INT: return Q_muli_to_int(x, c);
    2249             :     case t_FRAC:
    2250           0 :       n = gel(c,1);
    2251           0 :       d = gel(c,2);
    2252           0 :       return Q_divmuli_to_int(x, d,n);
    2253             :   }
    2254           0 :   pari_err_TYPE("Q_mul_to_int",c);
    2255             :   return NULL; /* LCOV_EXCL_LINE */
    2256             : }
    2257             : 
    2258             : GEN
    2259    14684978 : Q_primitive_part(GEN x, GEN *ptc)
    2260             : {
    2261    14684978 :   pari_sp av = avma;
    2262    14684978 :   GEN c = Q_content_safe(x);
    2263    14684978 :   if (c)
    2264             :   {
    2265    14684838 :     if (typ(c) == t_INT)
    2266             :     {
    2267    11495637 :       if (equali1(c)) { avma = av; c = NULL; }
    2268     2154008 :       else if (signe(c)) x = Q_divi_to_int(x, c);
    2269             :     }
    2270     3189201 :     else x = Q_divq_to_int(x, c);
    2271             :   }
    2272    14684978 :   if (ptc) *ptc = c;
    2273    14684978 :   return x;
    2274             : }
    2275             : GEN
    2276     1496608 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    2277             : 
    2278             : GEN
    2279       68989 : vec_Q_primpart(GEN x)
    2280       68989 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
    2281             : 
    2282             : /*******************************************************************/
    2283             : /*                                                                 */
    2284             : /*                           SUBRESULTANT                          */
    2285             : /*                                                                 */
    2286             : /*******************************************************************/
    2287             : /* for internal use */
    2288             : GEN
    2289     1812038 : gdivexact(GEN x, GEN y)
    2290             : {
    2291             :   long i,lx;
    2292             :   GEN z;
    2293     1812038 :   if (gequal1(y)) return x;
    2294     1808209 :   switch(typ(x))
    2295             :   {
    2296             :     case t_INT:
    2297     1442336 :       if (typ(y)==t_INT) return diviiexact(x,y);
    2298          35 :       if (!signe(x)) return gen_0;
    2299           0 :       break;
    2300             :     case t_INTMOD:
    2301             :     case t_FFELT:
    2302       14893 :     case t_POLMOD: return gmul(x,ginv(y));
    2303             :     case t_POL:
    2304      346990 :       switch(typ(y))
    2305             :       {
    2306             :         case t_INTMOD:
    2307             :         case t_FFELT:
    2308        5600 :         case t_POLMOD: return gmul(x,ginv(y));
    2309             :         case t_POL: { /* not stack-clean */
    2310             :           long v;
    2311       21224 :           if (varn(x)!=varn(y)) break;
    2312       20258 :           v = RgX_valrem(y,&y);
    2313       20258 :           if (v) x = RgX_shift_shallow(x,-v);
    2314       20258 :           if (!degpol(y)) { y = gel(y,2); break; }
    2315       18172 :           return RgX_div(x,y);
    2316             :         }
    2317             :       }
    2318      323218 :       return RgX_Rg_divexact(x, y);
    2319             :     case t_VEC: case t_COL: case t_MAT:
    2320        3990 :       lx = lg(x); z = new_chunk(lx);
    2321        3990 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    2322        3990 :       z[0] = x[0]; return z;
    2323             :   }
    2324           0 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    2325           0 :   return gdiv(x,y);
    2326             : }
    2327             : 
    2328             : static GEN
    2329       65013 : init_resultant(GEN x, GEN y)
    2330             : {
    2331       65013 :   long tx = typ(x), ty = typ(y), vx, vy;
    2332       65013 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    2333             :   {
    2334           0 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    2335           0 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    2336           0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    2337           0 :     return gen_1;
    2338             :   }
    2339       65013 :   if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
    2340       65013 :   if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
    2341       65013 :   if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
    2342       64936 :   vx = varn(x);
    2343       64936 :   vy = varn(y); if (vx == vy) return NULL;
    2344           0 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    2345             : }
    2346             : 
    2347             : static long
    2348      107759 : RgX_simpletype(GEN x)
    2349             : {
    2350      107759 :   long T = t_INT, i, lx = lg(x);
    2351      616415 :   for (i = 2; i < lx; i++)
    2352             :   {
    2353      517308 :     GEN c = gel(x,i);
    2354      517308 :     long tc = typ(c);
    2355      517308 :     switch(tc) {
    2356             :       case t_INT:
    2357      436626 :         break;
    2358             :       case t_FRAC:
    2359       37971 :         if (T == t_INT) T = t_FRAC;
    2360       37971 :         break;
    2361             :       default:
    2362       42711 :         if (isinexact(c)) return t_REAL;
    2363       34059 :         T = 0; break;
    2364             :     }
    2365             :   }
    2366       99107 :   return T;
    2367             : }
    2368             : 
    2369             : /* x an RgX, y a scalar */
    2370             : static GEN
    2371           0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    2372             : {
    2373           0 :   *V = gpowgs(y,degpol(x)-1);
    2374           0 :   *U = gen_0; return gmul(y, *V);
    2375             : }
    2376             : 
    2377             : /* return 0 if the subresultant chain can be interrupted.
    2378             :  * Set u = NULL if the resultant is 0. */
    2379             : static int
    2380        5165 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    2381             : {
    2382        5165 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    2383             :   long du, dv, dr, degq;
    2384             : 
    2385        5165 :   if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
    2386        5165 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    2387        5039 :   du = degpol(*u);
    2388        5039 :   dv = degpol(*v);
    2389        5039 :   degq = du - dv;
    2390        5039 :   if (*um1 == gen_1)
    2391        1847 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    2392        3192 :   else if (*um1 == gen_0)
    2393        1596 :     u0 = gen_0;
    2394             :   else /* except in those 2 cases, um1 is an RgX */
    2395        1596 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2396             : 
    2397        5039 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2398        1847 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2399             :   else
    2400        3192 :     u0 = gsub(u0, gmul(q,*uze));
    2401             : 
    2402        5039 :   *um1 = *uze;
    2403        5039 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2404             : 
    2405        5039 :   *u = *v; c = *g; *g  = leading_coeff(*u);
    2406        5039 :   switch(degq)
    2407             :   {
    2408          21 :     case 0: break;
    2409             :     case 1:
    2410        3912 :       c = gmul(*h,c); *h = *g; break;
    2411             :     default:
    2412        1106 :       c = gmul(gpowgs(*h,degq), c);
    2413        1106 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2414             :   }
    2415        5039 :   *v  = RgX_Rg_divexact(r,c);
    2416        5039 :   *uze= RgX_Rg_divexact(*uze,c);
    2417        5039 :   if (both_odd(du, dv)) *signh = -*signh;
    2418        5039 :   return (dr > 3);
    2419             : }
    2420             : 
    2421             : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2422             : static GEN
    2423        1742 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2424             : {
    2425             :   pari_sp av, av2;
    2426        1742 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2427             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2428             : 
    2429        1742 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2430        1742 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2431        1742 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2432        1742 :   if (tx != t_POL) {
    2433           0 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2434           0 :     return scalar_res(y,x,V,U);
    2435             :   }
    2436        1742 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2437        1742 :   if (varn(x) != varn(y))
    2438           0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2439           0 :                                         : scalar_res(y,x,V,U);
    2440        1742 :   if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
    2441        1742 :   if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
    2442        1742 :   dx = degpol(x);
    2443        1742 :   dy = degpol(y);
    2444        1742 :   signh = 1;
    2445        1742 :   if (dx < dy)
    2446             :   {
    2447         867 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2448         867 :     if (both_odd(dx, dy)) signh = -signh;
    2449             :   }
    2450        1742 :   if (dy == 0)
    2451             :   {
    2452           0 :     *V = gpowgs(gel(y,2),dx-1);
    2453           0 :     *U = gen_0; return gmul(*V,gel(y,2));
    2454             :   }
    2455        1742 :   av = avma;
    2456        1742 :   u = x = primitive_part(x, &cu);
    2457        1742 :   v = y = primitive_part(y, &cv);
    2458        1742 :   g = h = gen_1; av2 = avma;
    2459        1742 :   um1 = gen_1; uze = gen_0;
    2460             :   for(;;)
    2461             :   {
    2462        7650 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2463        2954 :     if (gc_needed(av2,1))
    2464             :     {
    2465           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2466           0 :       gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    2467             :     }
    2468             :   }
    2469             :   /* uze an RgX */
    2470        1742 :   if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
    2471        1735 :   z = gel(v,2); du = degpol(u);
    2472        1735 :   if (du > 1)
    2473             :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2474         167 :     p1 = gpowgs(gdiv(z,h),du-1);
    2475         167 :     z = gmul(z,p1);
    2476         167 :     uze = RgX_Rg_mul(uze, p1);
    2477             :   }
    2478        1735 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2479             : 
    2480        1735 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2481        1735 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2482             :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2483        1735 :   p1 = gen_1;
    2484        1735 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2485        1735 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2486        1735 :   cu = cu? gdiv(p1,cu): p1;
    2487        1735 :   cv = cv? gdiv(p1,cv): p1;
    2488        1735 :   z = gmul(z,p1);
    2489        1735 :   *U = RgX_Rg_mul(uze,cu);
    2490        1735 :   *V = RgX_Rg_mul(vze,cv);
    2491        1735 :   return z;
    2492             : }
    2493             : GEN
    2494           0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2495             : {
    2496           0 :   pari_sp av = avma;
    2497           0 :   GEN z = subresext_i(x, y, U, V);
    2498           0 :   gerepileall(av, 3, &z, U, V);
    2499           0 :   return z;
    2500             : }
    2501             : 
    2502             : static GEN
    2503          28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2504             : {
    2505          28 :   GEN x=content(y);
    2506          28 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2507             : }
    2508             : 
    2509             : static int
    2510      176624 : must_negate(GEN x)
    2511             : {
    2512      176624 :   GEN t = leading_coeff(x);
    2513      176624 :   switch(typ(t))
    2514             :   {
    2515             :     case t_INT: case t_REAL:
    2516      113309 :       return (signe(t) < 0);
    2517             :     case t_FRAC:
    2518         364 :       return (signe(gel(t,1)) < 0);
    2519             :   }
    2520       62951 :   return 0;
    2521             : }
    2522             : 
    2523             : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2524             : GEN
    2525         343 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2526             : {
    2527             :   pari_sp av, av2, tetpil;
    2528             :   long signh; /* junk */
    2529         343 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2530             :   GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
    2531             : 
    2532         343 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2533         343 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2534         343 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2535         343 :   vx=varn(x);
    2536         343 :   if (!signe(x))
    2537             :   {
    2538          14 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2539           7 :     *U = pol_0(vx); *V = pol_0(vx);
    2540           7 :     return pol_0(vx);
    2541             :   }
    2542         329 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2543         308 :   dx = degpol(x); dy = degpol(y);
    2544         308 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2545         308 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2546             : 
    2547         161 :   av = avma;
    2548         161 :   u = x = primitive_part(x, &cu);
    2549         161 :   v = y = primitive_part(y, &cv);
    2550         161 :   g = h = gen_1; av2 = avma;
    2551         161 :   um1 = gen_1; uze = gen_0;
    2552             :   for(;;)
    2553             :   {
    2554         203 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2555          21 :     if (gc_needed(av2,1))
    2556             :     {
    2557           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2558           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2559             :     }
    2560             :   }
    2561         161 :   if (uze != gen_0) {
    2562             :     GEN r;
    2563          42 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2564          42 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2565          42 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2566          42 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2567          42 :     p1 = ginv(content(v));
    2568             :   }
    2569             :   else /* y | x */
    2570             :   {
    2571         119 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2572         119 :     uze = pol_0(vx);
    2573         119 :     p1 = gen_1;
    2574             :   }
    2575         161 :   if (must_negate(v)) p1 = gneg(p1);
    2576         161 :   tetpil = avma;
    2577         161 :   z = RgX_Rg_mul(v,p1);
    2578         161 :   *U = RgX_Rg_mul(uze,p1);
    2579         161 :   *V = RgX_Rg_mul(vze,p1);
    2580         161 :   gptr[0] = &z;
    2581         161 :   gptr[1] = U;
    2582         161 :   gptr[2] = V;
    2583         161 :   gerepilemanysp(av,tetpil,gptr,3); return z;
    2584             : }
    2585             : 
    2586             : int
    2587          70 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2588             : {
    2589          70 :   pari_sp av = avma, av2, tetpil;
    2590             :   long signh; /* junk */
    2591             :   long vx;
    2592             :   GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
    2593             : 
    2594          70 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2595          70 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2596          70 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2597          70 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2598          70 :   if (!signe(T)) {
    2599           0 :     if (degpol(x) <= amax) {
    2600           0 :       *P = RgX_copy(x);
    2601           0 :       *Q = pol_1(varn(x));
    2602           0 :       return 1;
    2603             :     }
    2604           0 :     return 0;
    2605             :   }
    2606          70 :   if (amax+bmax >= degpol(T))
    2607           0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2608             :                     mkvec3(stoi(amax), stoi(bmax), T));
    2609          70 :   vx = varn(T);
    2610          70 :   u = x = primitive_part(x, &cu);
    2611          70 :   v = T = primitive_part(T, &cv);
    2612          70 :   g = h = gen_1; av2 = avma;
    2613          70 :   um1 = gen_1; uze = gen_0;
    2614             :   for(;;)
    2615             :   {
    2616         504 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2617         287 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
    2618         287 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2619         217 :     if (gc_needed(av2,1))
    2620             :     {
    2621           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2622           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2623             :     }
    2624             :   }
    2625          70 :   if (uze == gen_0)
    2626             :   {
    2627           0 :     avma = av; *P = pol_0(vx); *Q = pol_1(vx);
    2628           0 :     return 1;
    2629             :   }
    2630          70 :   if (cu) uze = RgX_Rg_div(uze,cu);
    2631          70 :   p1 = ginv(content(v));
    2632          70 :   if (must_negate(v)) p1 = gneg(p1);
    2633          70 :   tetpil = avma;
    2634          70 :   *P = RgX_Rg_mul(v,p1);
    2635          70 :   *Q = RgX_Rg_mul(uze,p1);
    2636          70 :   gptr[0] = P;
    2637          70 :   gptr[1] = Q;
    2638          70 :   gerepilemanysp(av,tetpil,gptr,2); return 1;
    2639             : }
    2640             : 
    2641             : /*******************************************************************/
    2642             : /*                                                                 */
    2643             : /*                RESULTANT USING DUCOS VARIANT                    */
    2644             : /*                                                                 */
    2645             : /*******************************************************************/
    2646             : /* x^n / y^(n-1), assume n > 0 */
    2647             : static GEN
    2648       20458 : Lazard(GEN x, GEN y, long n)
    2649             : {
    2650             :   long a;
    2651             :   GEN c;
    2652             : 
    2653       20458 :   if (n == 1) return x;
    2654        1369 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    2655        1369 :   c=x; n-=a;
    2656        4730 :   while (a>1)
    2657             :   {
    2658        1992 :     a>>=1; c=gdivexact(gsqr(c),y);
    2659        1992 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    2660             :   }
    2661        1369 :   return c;
    2662             : }
    2663             : 
    2664             : /* F (x/y)^(n-1), assume n >= 1 */
    2665             : static GEN
    2666       21196 : Lazard2(GEN F, GEN x, GEN y, long n)
    2667             : {
    2668       21196 :   if (n == 1) return F;
    2669         742 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    2670             : }
    2671             : 
    2672             : static GEN
    2673       21196 : RgX_neg_i(GEN x, long lx)
    2674             : {
    2675             :   long i;
    2676       21196 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    2677       21196 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    2678       21196 :   return y;
    2679             : }
    2680             : static GEN
    2681       63742 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    2682             : {
    2683             :   long i;
    2684             :   GEN z;
    2685       63742 :   if (isrationalzero(x)) return pol_0(varn(y));
    2686       63721 :   z = cgetg(ly,t_POL); z[1] = y[1];
    2687       63721 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2688       63721 :   return z;
    2689             : }
    2690             : static long
    2691       60837 : reductum_lg(GEN x, long lx)
    2692             : {
    2693       60837 :   long i = lx-2;
    2694       60837 :   while (i > 1 && gequal0(gel(x,i))) i--;
    2695       60837 :   return i+1;
    2696             : }
    2697             : 
    2698             : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
    2699             : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    2700             :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    2701             : static GEN
    2702       21196 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    2703             : {
    2704       21196 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
    2705             :   long p, q, j, lP, lQ;
    2706             :   pari_sp av;
    2707             : 
    2708       21196 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    2709       21196 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    2710             :   /* p > q. Very often p - 1 = q */
    2711       21196 :   av = avma;
    2712             :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    2713       21196 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    2714             : 
    2715       21196 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    2716       24843 :   for (j = q+1; j < p; j++)
    2717             :   {
    2718        3647 :     if (degpol(H) == q-1)
    2719             :     { /* h0 = coeff of degree q-1 = leading coeff */
    2720        3213 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    2721        3213 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    2722             :     }
    2723             :     else
    2724         434 :       H = RgX_shift_shallow(H, 1);
    2725        3647 :     if (j+2 < lP)
    2726             :     {
    2727        3073 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    2728        3073 :       A = A? RgX_add(A, TMP): TMP;
    2729             :     }
    2730        3647 :     if (gc_needed(av,1))
    2731             :     {
    2732         147 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    2733         147 :       gerepileall(av,A?2:1,&H,&A);
    2734             :     }
    2735             :   }
    2736       21196 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    2737       21196 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    2738       21196 :   A = A? RgX_add(A, TMP): TMP;
    2739       21196 :   A = RgX_Rg_divexact(A, p0);
    2740       21196 :   if (degpol(H) == q-1)
    2741             :   {
    2742       20888 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    2743       20888 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    2744             :   }
    2745             :   else
    2746         308 :     A = RgX_Rg_mul(addshift(H,A), q0);
    2747       21196 :   return RgX_Rg_divexact(A, s);
    2748             : }
    2749             : #undef addshift
    2750             : 
    2751             : /* Ducos's subresultant */
    2752             : GEN
    2753       21095 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    2754             : {
    2755             :   pari_sp av, av2;
    2756       21095 :   long dP, dQ, delta, sig = 1;
    2757             :   GEN cP, cQ, Z, s;
    2758             : 
    2759       21095 :   dP = degpol(P);
    2760       21095 :   dQ = degpol(Q); delta = dP - dQ;
    2761       21095 :   if (delta < 0)
    2762             :   {
    2763        1197 :     if (both_odd(dP, dQ)) sig = -1;
    2764        1197 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    2765             :   }
    2766       21095 :   if (sol) *sol = gen_0;
    2767       21095 :   av = avma;
    2768       21095 :   if (dQ <= 0)
    2769             :   {
    2770        1379 :     if (dQ < 0) return Rg_get_0(P);
    2771        1379 :     s = gpowgs(gel(Q,2), dP);
    2772        1379 :     if (sig == -1) s = gerepileupto(av, gneg(s));
    2773        1379 :     return s;
    2774             :   }
    2775             :   /* primitive_part is also possible here, but possibly very costly,
    2776             :    * and hardly ever worth it */
    2777       19716 :   P = Q_primitive_part(P, &cP);
    2778       19716 :   Q = Q_primitive_part(Q, &cQ);
    2779       19716 :   av2 = avma;
    2780       19716 :   s = gpowgs(leading_coeff(Q),delta);
    2781       19716 :   if (both_odd(dP, dQ)) sig = -sig;
    2782       19716 :   Z = Q;
    2783       19716 :   Q = RgX_pseudorem(P, Q);
    2784       19716 :   P = Z;
    2785       60628 :   while(degpol(Q) > 0)
    2786             :   {
    2787       21196 :     delta = degpol(P) - degpol(Q); /* > 0 */
    2788       21196 :     Z = Lazard2(Q, leading_coeff(Q), s, delta);
    2789       21196 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    2790       21196 :     Q = nextSousResultant(P, Q, Z, s);
    2791       21196 :     P = Z;
    2792       21196 :     if (gc_needed(av,1))
    2793             :     {
    2794          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    2795          13 :       gerepileall(av2,2,&P,&Q);
    2796             :     }
    2797       21196 :     s = leading_coeff(P);
    2798             :   }
    2799       19716 :   if (!signe(Q)) { avma = av; return Rg_get_0(Q); }
    2800       19716 :   s = Lazard(leading_coeff(Q), s, degpol(P));
    2801       19716 :   if (sig == -1) s = gneg(s);
    2802       19716 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    2803       19716 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    2804       19716 :   if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
    2805       18197 :   return gerepilecopy(av, s);
    2806             : }
    2807             : /* Return resultant(P,Q).
    2808             :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    2809             :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    2810             :  * in the "generic" case. */
    2811             : GEN
    2812       56354 : resultant(GEN P, GEN Q)
    2813             : {
    2814             :   long TP, TQ;
    2815       56354 :   GEN s, p = NULL;
    2816             : 
    2817       56354 :   if ((s = init_resultant(P,Q))) return s;
    2818       56277 :   if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
    2819        8638 :     return resultant2(P,Q); /* inexact */
    2820       47639 :   if (TP && TQ) /* rational */
    2821             :   {
    2822       28833 :     if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
    2823       13035 :     return QX_resultant(P,Q);
    2824             :   }
    2825       18806 :   if (RgX_is_FpX(P, &p) && RgX_is_FpX(Q, &p) && p)
    2826             :   {
    2827          14 :     pari_sp av = avma;
    2828          14 :     GEN r = FpX_resultant(RgX_to_FpX(P, p), RgX_to_FpX(Q, p), p);
    2829          14 :     return gerepileupto(av, Fp_to_mod(r, p));
    2830             :   }
    2831       18792 :   return RgX_resultant_all(P, Q, NULL);
    2832             : }
    2833             : 
    2834             : /*******************************************************************/
    2835             : /*                                                                 */
    2836             : /*               RESULTANT USING SYLVESTER MATRIX                  */
    2837             : /*                                                                 */
    2838             : /*******************************************************************/
    2839             : static GEN
    2840       72898 : syl_RgC(GEN x, long j, long d, long D, long cp)
    2841             : {
    2842       72898 :   GEN c = cgetg(d+1,t_COL);
    2843             :   long i;
    2844       72898 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    2845       72898 :   for (   ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
    2846       72898 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    2847       72898 :   return c;
    2848             : }
    2849             : static GEN
    2850        8666 : syl_RgM(GEN x, GEN y, long cp)
    2851             : {
    2852        8666 :   long j, d, dx = degpol(x), dy = degpol(y);
    2853             :   GEN M;
    2854        8666 :   if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
    2855        8666 :   if (dy < 0) return zeromat(dx,dx);
    2856        8666 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    2857        8666 :   for (j=1; j<=dy; j++) gel(M,j)    = syl_RgC(x,j,d,j+dx, cp);
    2858        8666 :   for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
    2859        8666 :   return M;
    2860             : }
    2861             : GEN
    2862        8659 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
    2863             : GEN
    2864           7 : sylvestermatrix(GEN x, GEN y)
    2865             : {
    2866           7 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    2867           7 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    2868           7 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    2869           7 :   return syl_RgM(x,y,1);
    2870             : }
    2871             : 
    2872             : GEN
    2873        8659 : resultant2(GEN x, GEN y)
    2874             : {
    2875        8659 :   pari_sp av = avma;
    2876        8659 :   GEN r = init_resultant(x,y);
    2877        8659 :   return r? r: gerepileupto(av, det(RgX_sylvestermatrix(x,y)));
    2878             : }
    2879             : 
    2880             : /* If x a t_POL, let vx = main variable of x; return a t_POL in variable v0:
    2881             :  * if vx <= v, return subst(x, v, pol_x(v0))
    2882             :  * if vx >  v, return scalarpol(x, v0) */
    2883             : static GEN
    2884        3626 : fix_pol(GEN x, long v, long v0)
    2885             : {
    2886             :   long vx;
    2887        3626 :   if (typ(x) != t_POL) return x;
    2888        3626 :   vx = varn(x);
    2889        3626 :   if (v == vx)
    2890             :   {
    2891        3479 :     if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
    2892        3479 :     return x;
    2893             :   }
    2894         147 :   if (varncmp(v, vx) > 0)
    2895             :   {
    2896         140 :     x = gsubst(x,v,pol_x(v0));
    2897         140 :     if (typ(x) == t_POL && varn(x) == v0) return x;
    2898             :   }
    2899          35 :   return scalarpol_shallow(x, v0);
    2900             : }
    2901             : 
    2902             : /* resultant of x and y with respect to variable v, or with respect to their
    2903             :  * main variable if v < 0. */
    2904             : GEN
    2905        2030 : polresultant0(GEN x, GEN y, long v, long flag)
    2906             : {
    2907        2030 :   long v0 = 0;
    2908        2030 :   pari_sp av = avma;
    2909             : 
    2910        2030 :   if (v >= 0)
    2911             :   {
    2912        1799 :     v0 = fetch_var_higher();
    2913        1799 :     x = fix_pol(x,v, v0);
    2914        1799 :     y = fix_pol(y,v, v0);
    2915             :   }
    2916        2030 :   switch(flag)
    2917             :   {
    2918             :     case 2:
    2919        2023 :     case 0: x=resultant(x,y); break;
    2920           7 :     case 1: x=resultant2(x,y); break;
    2921           0 :     default: pari_err_FLAG("polresultant");
    2922             :   }
    2923        2030 :   if (v >= 0) (void)delete_var();
    2924        2030 :   return gerepileupto(av,x);
    2925             : }
    2926             : GEN
    2927         875 : polresultantext0(GEN x, GEN y, long v)
    2928             : {
    2929             :   GEN R, U, V;
    2930         875 :   long v0 = 0;
    2931         875 :   pari_sp av = avma;
    2932             : 
    2933         875 :   if (v >= 0)
    2934             :   {
    2935          14 :     v0 = fetch_var_higher();
    2936          14 :     x = fix_pol(x,v, v0);
    2937          14 :     y = fix_pol(y,v, v0);
    2938             :   }
    2939         875 :   R = subresext_i(x,y, &U,&V);
    2940         875 :   if (v >= 0)
    2941             :   {
    2942          14 :     (void)delete_var();
    2943          14 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    2944          14 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    2945             :   }
    2946         875 :   return gerepilecopy(av, mkvec3(U,V,R));
    2947             : }
    2948             : GEN
    2949         847 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    2950             : 
    2951             : /*******************************************************************/
    2952             : /*                                                                 */
    2953             : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    2954             : /*                                                                 */
    2955             : /*******************************************************************/
    2956             : 
    2957             : /* (v - x)^d */
    2958             : static GEN
    2959         140 : caract_const(pari_sp av, GEN x, long v, long d)
    2960         140 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    2961             : 
    2962             : /* return caract(Mod(x,T)) in variable v */
    2963             : GEN
    2964       15635 : RgXQ_charpoly(GEN x, GEN T, long v)
    2965             : {
    2966       15635 :   pari_sp av = avma;
    2967       15635 :   long d = degpol(T), dx, vx, vp, v0;
    2968             :   GEN ch, L;
    2969             : 
    2970       15635 :   if (typ(x) != t_POL) return caract_const(av, x, v, d);
    2971       15621 :   vx = varn(x);
    2972       15621 :   vp = varn(T);
    2973       15621 :   if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
    2974       15621 :   if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
    2975       15621 :   dx = degpol(x);
    2976       15621 :   if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
    2977       15621 :   if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
    2978             : 
    2979       15565 :   v0 = fetch_var_higher();
    2980       15565 :   x = RgX_neg(x);
    2981       15565 :   gel(x,2) = gadd(gel(x,2), pol_x(v));
    2982       15565 :   setvarn(x, v0);
    2983       15565 :   T = leafcopy(T); setvarn(T, v0);
    2984       15565 :   ch = resultant(T, x);
    2985       15565 :   (void)delete_var();
    2986             :   /* test for silly input: x mod (deg 0 polynomial) */
    2987       15565 :   if (typ(ch) != t_POL)
    2988           7 :     pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
    2989       15558 :   L = leading_coeff(ch);
    2990       15558 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    2991       15558 :   return gerepileupto(av, ch);
    2992             : }
    2993             : 
    2994             : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    2995             :  * algebra nf[t]/(Q(t)) */
    2996             : GEN
    2997         224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    2998             : {
    2999         224 :   const char *f = "rnfcharpoly";
    3000         224 :   long dQ = degpol(Q);
    3001         224 :   pari_sp av = avma;
    3002             :   GEN T;
    3003             : 
    3004         224 :   if (v < 0) v = 0;
    3005         224 :   nf = checknf(nf); T = nf_get_pol(nf);
    3006         224 :   Q = RgX_nffix(f, T,Q,0);
    3007         224 :   switch(typ(x))
    3008             :   {
    3009             :     case t_INT:
    3010          28 :     case t_FRAC: return caract_const(av, x, v, dQ);
    3011             :     case t_POLMOD:
    3012          91 :       x = polmod_nffix2(f,T,Q, x,0);
    3013          56 :       break;
    3014             :     case t_POL:
    3015          56 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    3016          42 :       break;
    3017          49 :     default: pari_err_TYPE(f,x);
    3018             :   }
    3019          98 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    3020             :   /* x a t_POL in variable vQ */
    3021          56 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    3022          56 :   if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
    3023          56 :   return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    3024             : }
    3025             : 
    3026             : /*******************************************************************/
    3027             : /*                                                                 */
    3028             : /*                  GCD USING SUBRESULTANT                         */
    3029             : /*                                                                 */
    3030             : /*******************************************************************/
    3031             : static int inexact(GEN x, int *simple, int *rational);
    3032             : static int
    3033     7118581 : isinexactall(GEN x, int *simple, int *rational)
    3034             : {
    3035     7118581 :   long i, lx = lg(x);
    3036    35028821 :   for (i=2; i<lx; i++)
    3037    27910254 :     if (inexact(gel(x,i), simple, rational)) return 1;
    3038     7118567 :   return 0;
    3039             : }
    3040             : /* return 1 if coeff explosion is not possible */
    3041             : static int
    3042    27910464 : inexact(GEN x, int *simple, int *rational)
    3043             : {
    3044    27910464 :   int junk = 0;
    3045    27910464 :   switch(typ(x))
    3046             :   {
    3047    25109708 :     case t_INT: case t_FRAC: return 0;
    3048             : 
    3049           7 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    3050             : 
    3051             :     case t_INTMOD:
    3052             :     case t_FFELT:
    3053     2758798 :       *rational = 0;
    3054     2758798 :       if (!*simple) *simple = 1;
    3055     2758798 :       return 0;
    3056             : 
    3057             :     case t_COMPLEX:
    3058          77 :       *rational = 0;
    3059          77 :       return inexact(gel(x,1), simple, rational)
    3060          77 :           || inexact(gel(x,2), simple, rational);
    3061             :     case t_QUAD:
    3062           0 :       *rational = *simple = 0;
    3063           0 :       return inexact(gel(x,2), &junk, rational)
    3064           0 :           || inexact(gel(x,3), &junk, rational);
    3065             : 
    3066             :     case t_POLMOD:
    3067       11718 :       *rational = 0;
    3068       11718 :       return isinexactall(gel(x,1), simple, rational);
    3069             :     case t_POL:
    3070       30128 :       *rational = 0;
    3071       30128 :       *simple = -1;
    3072       30128 :       return isinexactall(x, &junk, rational);
    3073             :     case t_RFRAC:
    3074          28 :       *rational = 0;
    3075          28 :       *simple = -1;
    3076          28 :       return inexact(gel(x,1), &junk, rational)
    3077          28 :           || inexact(gel(x,2), &junk, rational);
    3078             :   }
    3079           0 :   *rational = 0;
    3080           0 :   *simple = -1; return 0;
    3081             : }
    3082             : 
    3083             : /* x monomial, y t_POL in the same variable */
    3084             : static GEN
    3085     7528465 : gcdmonome(GEN x, GEN y)
    3086             : {
    3087     7528465 :   pari_sp av = avma;
    3088     7528465 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    3089     7528465 :   long i, l = lg(y);
    3090     7528465 :   GEN t, v = cgetg(l, t_VEC);
    3091     7528465 :   gel(v,1) = gel(x,dx+2);
    3092     7528465 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    3093     7528465 :   t = content(v); /* gcd(lc(x), cont(y)) */
    3094     7528465 :   t = simplify_shallow(t);
    3095     7528465 :   if (dx < e) e = dx;
    3096     7528465 :   return gerepileupto(av, monomialcopy(t, e, varn(x)));
    3097             : }
    3098             : 
    3099             : /* x, y are t_POL in the same variable */
    3100             : GEN
    3101    11067249 : RgX_gcd(GEN x, GEN y)
    3102             : {
    3103             :   long dx, dy;
    3104             :   pari_sp av, av1;
    3105             :   GEN d, g, h, p1, p2, u, v;
    3106    11067249 :   int simple = 0, rational = 1;
    3107             : 
    3108    11067249 :   if (isexactzero(y)) return RgX_copy(x);
    3109    11066843 :   if (isexactzero(x)) return RgX_copy(y);
    3110    11066836 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    3111     4197196 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    3112     3538371 :   if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
    3113             :   {
    3114           7 :     av = avma; u = ggcd(content(x), content(y));
    3115           7 :     return gerepileupto(av, scalarpol(u, varn(x)));
    3116             :   }
    3117     3538364 :   if (rational) return QX_gcd(x,y); /* Q[X] */
    3118             : 
    3119      178556 :   av = avma;
    3120      178556 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    3121             :   else
    3122             :   {
    3123        4963 :     dx = lg(x); dy = lg(y);
    3124        4963 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    3125        4963 :     if (dy==3)
    3126             :     {
    3127           0 :       d = ggcd(gel(y,2), content(x));
    3128           0 :       return gerepileupto(av, scalarpol(d, varn(x)));
    3129             :     }
    3130        4963 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    3131        4963 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    3132        4963 :     d = ggcd(p1,p2);
    3133        4963 :     av1 = avma;
    3134        4963 :     g = h = gen_1;
    3135             :     for(;;)
    3136        1491 :     {
    3137        6454 :       GEN r = RgX_pseudorem(u,v);
    3138        6454 :       long degq, du, dv, dr = lg(r);
    3139             : 
    3140        6454 :       if (!signe(r)) break;
    3141        3654 :       if (dr <= 3)
    3142             :       {
    3143        2163 :         avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
    3144             :       }
    3145        1491 :       if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
    3146        1491 :       du = lg(u); dv = lg(v); degq = du-dv;
    3147        1491 :       u = v; p1 = g; g = leading_coeff(u);
    3148        1491 :       switch(degq)
    3149             :       {
    3150         189 :         case 0: break;
    3151             :         case 1:
    3152        1064 :           p1 = gmul(h,p1); h = g; break;
    3153             :         default:
    3154         238 :           p1 = gmul(gpowgs(h,degq), p1);
    3155         238 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    3156             :       }
    3157        1491 :       v = RgX_Rg_div(r,p1);
    3158        1491 :       if (gc_needed(av1,1))
    3159             :       {
    3160           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
    3161           0 :         gerepileall(av1,4, &u,&v,&g,&h);
    3162             :       }
    3163             :     }
    3164        2800 :     x = RgX_Rg_mul(primpart(v), d);
    3165             :   }
    3166      176393 :   if (must_negate(x)) x = RgX_neg(x);
    3167      176393 :   return gerepileupto(av,x);
    3168             : }
    3169             : 
    3170             : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
    3171             : static GEN
    3172        5432 : RgX_disc_aux(GEN P)
    3173             : {
    3174        5432 :   long n = degpol(P), TP, dd;
    3175             :   GEN D, L, y, p;
    3176        5432 :   if (!signe(P) || !n) return Rg_get_0(P);
    3177        5432 :   if (n == 1) return Rg_get_1(P);
    3178        5257 :   if (n == 2) {
    3179        1414 :     GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
    3180        1414 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    3181             :   }
    3182        3843 :   TP = RgX_simpletype(P);
    3183        3843 :   if (TP == t_INT) return ZX_disc(P);
    3184         462 :   if (TP == t_FRAC) return QX_disc(P);
    3185         462 :   p = NULL;
    3186         462 :   if (RgX_is_FpX(P, &p) && p)
    3187          28 :     return Fp_to_mod(FpX_disc(RgX_to_FpX(P,p), p), p);
    3188             : 
    3189         434 :   y = RgX_deriv(P);
    3190         434 :   if (!signe(y)) return Rg_get_0(y);
    3191         434 :   dd = degpol(P)-2 - degpol(y);
    3192         434 :   if (TP == t_REAL)
    3193          14 :     D = resultant2(P,y);
    3194             :   else
    3195             :   {
    3196         420 :     D = RgX_resultant_all(P, y, NULL);
    3197         420 :     if (D == gen_0) return Rg_get_0(y);
    3198             :   }
    3199         434 :   L = leading_coeff(P);
    3200         434 :   if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
    3201         434 :   if (n & 2) D = gneg(D);
    3202         434 :   return D;
    3203             : }
    3204             : GEN
    3205        2240 : RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }
    3206             : 
    3207             : GEN
    3208        3206 : poldisc0(GEN x, long v)
    3209             : {
    3210             :   pari_sp av;
    3211        3206 :   switch(typ(x))
    3212             :   {
    3213             :     case t_POL:
    3214             :     {
    3215             :       GEN D;
    3216        3192 :       long v0 = -1;
    3217        3192 :       av = avma;
    3218        3192 :       if (v >= 0 && v != varn(x))
    3219             :       {
    3220           0 :         v0 = fetch_var_higher();
    3221           0 :         x = fix_pol(x,v, v0);
    3222             :       }
    3223        3192 :       D = RgX_disc_aux(x);
    3224        3192 :       if (v0 >= 0) (void)delete_var();
    3225        3192 :       return gerepileupto(av, D);
    3226             :     }
    3227             : 
    3228             :     case t_COMPLEX:
    3229           0 :       return utoineg(4);
    3230             : 
    3231             :     case t_QUAD:
    3232           0 :       return quad_disc(x);
    3233             :     case t_POLMOD:
    3234           0 :       return poldisc0(gel(x,1), v);
    3235             : 
    3236             :     case t_QFR: case t_QFI:
    3237          14 :       av = avma; return gerepileuptoint(av, qfb_disc(x));
    3238             : 
    3239             :     case t_VEC: case t_COL: case t_MAT:
    3240             :     {
    3241             :       long i;
    3242           0 :       GEN z = cgetg_copy(x, &i);
    3243           0 :       for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
    3244           0 :       return z;
    3245             :     }
    3246             :   }
    3247           0 :   pari_err_TYPE("poldisc",x);
    3248             :   return NULL; /* LCOV_EXCL_LINE */
    3249             : }
    3250             : 
    3251             : GEN
    3252           7 : reduceddiscsmith(GEN x)
    3253             : {
    3254           7 :   long j, n = degpol(x);
    3255           7 :   pari_sp av = avma;
    3256             :   GEN xp, M;
    3257             : 
    3258           7 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    3259           7 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    3260           7 :   RgX_check_ZX(x,"poldiscreduced");
    3261           7 :   if (!gequal1(gel(x,n+2)))
    3262           0 :     pari_err_IMPL("non-monic polynomial in poldiscreduced");
    3263           7 :   M = cgetg(n+1,t_MAT);
    3264           7 :   xp = ZX_deriv(x);
    3265          28 :   for (j=1; j<=n; j++)
    3266             :   {
    3267          21 :     gel(M,j) = RgX_to_RgC(xp, n);
    3268          21 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    3269             :   }
    3270           7 :   return gerepileupto(av, ZM_snf(M));
    3271             : }
    3272             : 
    3273             : /***********************************************************************/
    3274             : /**                                                                   **/
    3275             : /**                       STURM ALGORITHM                             **/
    3276             : /**              (number of real roots of x in [a,b])                 **/
    3277             : /**                                                                   **/
    3278             : /***********************************************************************/
    3279             : static GEN
    3280         385 : R_to_Q_up(GEN x)
    3281             : {
    3282             :   long e;
    3283         385 :   switch(typ(x))
    3284             :   {
    3285         385 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3286             :     case t_REAL:
    3287           0 :       x = mantissa_real(x,&e);
    3288           0 :       return gmul2n(addiu(x,1), -e);
    3289           0 :     default: pari_err_TYPE("R_to_Q_up", x);
    3290             :       return NULL; /* LCOV_EXCL_LINE */
    3291             :   }
    3292             : }
    3293             : static GEN
    3294         385 : R_to_Q_down(GEN x)
    3295             : {
    3296             :   long e;
    3297         385 :   switch(typ(x))
    3298             :   {
    3299         371 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3300             :     case t_REAL:
    3301          14 :       x = mantissa_real(x,&e);
    3302          14 :       return gmul2n(subiu(x,1), -e);
    3303           0 :     default: pari_err_TYPE("R_to_Q_down", x);
    3304             :       return NULL; /* LCOV_EXCL_LINE */
    3305             :   }
    3306             : }
    3307             : 
    3308             : static long
    3309         385 : sturmpart_i(GEN x, GEN ab)
    3310             : {
    3311         385 :   long tx = typ(x);
    3312         385 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    3313         385 :   if (tx != t_POL)
    3314             :   {
    3315           0 :     if (is_real_t(tx)) return 0;
    3316           0 :     pari_err_TYPE("sturm",x);
    3317             :   }
    3318         385 :   if (lg(x) == 3) return 0;
    3319         385 :   if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
    3320         385 :   (void)ZX_gcd_all(x, ZX_deriv(x), &x);
    3321         385 :   if (ab)
    3322             :   {
    3323             :     GEN A, B;
    3324         385 :     if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
    3325         385 :     A = R_to_Q_down(gel(ab,1));
    3326         385 :     B = R_to_Q_up(gel(ab,2));
    3327         385 :     ab = mkvec2(A, B);
    3328             :   }
    3329         385 :   return ZX_sturmpart(x, ab);
    3330             : }
    3331             : /* Deprecated: RgX_sturmpart() should be preferred  */
    3332             : long
    3333         385 : sturmpart(GEN x, GEN a, GEN b)
    3334             : {
    3335         385 :   pari_sp av = avma;
    3336             :   long r;
    3337         385 :   if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
    3338         252 :   if (!a) a = mkmoo();
    3339         252 :   if (!b) b = mkoo();
    3340         252 :   r = sturmpart_i(x, mkvec2(a, b));
    3341         252 :   avma = av; return r;
    3342             : }
    3343             : long
    3344         133 : RgX_sturmpart(GEN x, GEN ab)
    3345             : {
    3346         133 :   pari_sp av = avma;
    3347         133 :   long r = sturmpart_i(x, ab);
    3348         133 :   avma = av; return r;
    3349             : }
    3350             : 
    3351             : /***********************************************************************/
    3352             : /**                                                                   **/
    3353             : /**                        GENERIC EXTENDED GCD                       **/
    3354             : /**                                                                   **/
    3355             : /***********************************************************************/
    3356             : /* assume typ(x) = typ(y) = t_POL */
    3357             : static GEN
    3358         867 : RgXQ_inv_i(GEN x, GEN y)
    3359             : {
    3360         867 :   long vx=varn(x), vy=varn(y);
    3361             :   pari_sp av;
    3362             :   GEN u, v, d;
    3363             : 
    3364        1734 :   while (vx != vy)
    3365             :   {
    3366           0 :     if (varncmp(vx,vy) > 0)
    3367             :     {
    3368           0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3369           0 :       return scalarpol(d, vy);
    3370             :     }
    3371           0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3372           0 :     x = gel(x,2); vx = gvar(x);
    3373             :   }
    3374         867 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3375         867 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3376         867 :   d = gdiv(u,d);
    3377         867 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3378         867 :   return gerepileupto(av, d);
    3379             : }
    3380             : 
    3381             : /*Assume x is a polynomial and y is not */
    3382             : static GEN
    3383         112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3384             : {
    3385         112 :   long vx = varn(x);
    3386         112 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    3387         112 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    3388          84 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    3389          56 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    3390             : }
    3391             : /* Assume x==0, y!=0 */
    3392             : static GEN
    3393          63 : zero_bezout(GEN y, GEN *U, GEN *V)
    3394             : {
    3395          63 :   *U=gen_0; *V = ginv(y); return gen_1;
    3396             : }
    3397             : 
    3398             : GEN
    3399         280 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    3400             : {
    3401         280 :   long tx=typ(x), ty=typ(y), vx;
    3402         280 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    3403         245 :   if (tx != t_POL)
    3404             :   {
    3405         140 :     if (ty == t_POL)
    3406          56 :       return scalar_bezout(y,x,v,u);
    3407             :     else
    3408             :     {
    3409          84 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    3410          84 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    3411          63 :       if (xis0) return zero_bezout(y,u,v);
    3412          42 :       else return zero_bezout(x,v,u);
    3413             :     }
    3414             :   }
    3415         105 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    3416          49 :   vx = varn(x);
    3417          49 :   if (vx != varn(y))
    3418           0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    3419           0 :                                    : scalar_bezout(y,x,v,u);
    3420          49 :   return RgX_extgcd(x,y,u,v);
    3421             : }
    3422             : 
    3423             : GEN
    3424         280 : gcdext0(GEN x, GEN y)
    3425             : {
    3426         280 :   GEN z=cgetg(4,t_VEC);
    3427         280 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    3428         280 :   return z;
    3429             : }
    3430             : 
    3431             : /*******************************************************************/
    3432             : /*                                                                 */
    3433             : /*                    GENERIC (modular) INVERSE                    */
    3434             : /*                                                                 */
    3435             : /*******************************************************************/
    3436             : 
    3437             : GEN
    3438       14897 : ginvmod(GEN x, GEN y)
    3439             : {
    3440       14897 :   long tx=typ(x);
    3441             : 
    3442       14897 :   switch(typ(y))
    3443             :   {
    3444             :     case t_POL:
    3445       14897 :       if (tx==t_POL) return RgXQ_inv(x,y);
    3446        2016 :       if (is_scalar_t(tx)) return ginv(x);
    3447           0 :       break;
    3448             :     case t_INT:
    3449           0 :       if (tx==t_INT) return Fp_inv(x,y);
    3450           0 :       if (tx==t_POL) return gen_0;
    3451             :   }
    3452           0 :   pari_err_TYPE2("ginvmod",x,y);
    3453             :   return NULL; /* LCOV_EXCL_LINE */
    3454             : }
    3455             : 
    3456             : /***********************************************************************/
    3457             : /**                                                                   **/
    3458             : /**                         NEWTON POLYGON                            **/
    3459             : /**                                                                   **/
    3460             : /***********************************************************************/
    3461             : 
    3462             : /* assume leading coeff of x is non-zero */
    3463             : GEN
    3464          28 : newtonpoly(GEN x, GEN p)
    3465             : {
    3466             :   GEN y;
    3467             :   long n,ind,a,b,c,u1,u2,r1,r2;
    3468          28 :   long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
    3469             : 
    3470          28 :   if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
    3471          28 :   n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
    3472          28 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    3473          28 :   vval = (long *) pari_malloc(sizeof(long)*(n+1));
    3474          28 :   for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
    3475          42 :   for (a=0, ind=1; a<n; a++)
    3476             :   {
    3477          42 :     if (vval[a] != LONG_MAX) break;
    3478          14 :     gel(y,ind++) = mkoo();
    3479             :   }
    3480          84 :   for (b=a+1; b<=n; a=b, b=a+1)
    3481             :   {
    3482          56 :     while (vval[b] == LONG_MAX) b++;
    3483          56 :     u1 = vval[a]-vval[b];
    3484          56 :     u2 = b-a;
    3485         154 :     for (c=b+1; c<=n; c++)
    3486             :     {
    3487          98 :       if (vval[c] == LONG_MAX) continue;
    3488          70 :       r1 = vval[a]-vval[c];
    3489          70 :       r2 = c-a;
    3490          70 :       if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
    3491             :     }
    3492          56 :     while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
    3493             :   }
    3494          28 :   pari_free(vval); return y;
    3495             : }
    3496             : 
    3497             : static GEN
    3498      262591 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
    3499             : {
    3500      262591 :   pari_sp av = avma;
    3501             :   GEN r;
    3502      262591 :   if (lgefint(p) == 3)
    3503             :   {
    3504      144256 :     ulong pp = uel(p, 2);
    3505      144256 :     r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
    3506             :                 RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
    3507             :   }
    3508             :   else
    3509      118335 :     r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
    3510      262591 :   return gerepileupto(av, FpX_to_mod(r, p));
    3511             : }
    3512             : 
    3513             : static GEN
    3514          14 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
    3515             : {
    3516          14 :   pari_sp av = avma;
    3517             :   GEN r;
    3518          14 :   if (lgefint(p) == 3)
    3519             :   {
    3520           7 :     ulong pp = uel(p, 2);
    3521           7 :     r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
    3522             :                                    RgX_to_Flx(y, pp), pp));
    3523             :   }
    3524             :   else
    3525           7 :     r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3526          14 :   return gerepileupto(av, FpX_to_mod(r, p));
    3527             : }
    3528             : 
    3529             : static GEN
    3530       11571 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
    3531             : {
    3532       11571 :   pari_sp av = avma;
    3533             :   GEN r;
    3534       11571 :   if (lgefint(p) == 3)
    3535             :   {
    3536        5789 :     ulong pp = uel(p, 2);
    3537        5789 :     r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
    3538             :                                    RgX_to_Flx(y, pp), pp));
    3539             :   }
    3540             :   else
    3541        5782 :     r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3542       11571 :   return gerepileupto(av, FpX_to_mod(r, p));
    3543             : }
    3544             : 
    3545             : static GEN
    3546           7 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
    3547             : {
    3548           7 :   pari_sp av = avma;
    3549             :   GEN r;
    3550           7 :   GEN T = RgX_to_FpX(pol, p);
    3551           7 :   if (signe(T)==0) pari_err_OP("*",x,y);
    3552           7 :   if (lgefint(p) == 3)
    3553             :   {
    3554           7 :     ulong pp = uel(p, 2);
    3555           7 :     GEN Tp = ZX_to_Flx(T, pp);
    3556           7 :     r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
    3557             :                                RgX_to_FlxqX(y, Tp, pp),
    3558             :                                RgX_to_FlxqX(S, Tp, pp), Tp, pp));
    3559             :   }
    3560             :   else
    3561           0 :     r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
    3562             :                    RgX_to_FpXQX(S, T, p), T, p);
    3563           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3564             : }
    3565             : 
    3566             : static GEN
    3567           0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    3568             : {
    3569           0 :   pari_sp av = avma;
    3570             :   GEN r;
    3571           0 :   GEN T = RgX_to_FpX(pol, p);
    3572           0 :   if (signe(T)==0) pari_err_OP("*",x,x);
    3573           0 :   if (lgefint(p) == 3)
    3574             :   {
    3575           0 :     ulong pp = uel(p, 2);
    3576           0 :     GEN Tp = ZX_to_Flx(T, pp);
    3577           0 :     r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
    3578             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3579             :   }
    3580             :   else
    3581           0 :     r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3582           0 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3583             : }
    3584             : 
    3585             : static GEN
    3586           7 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    3587             : {
    3588           7 :   pari_sp av = avma;
    3589             :   GEN r;
    3590           7 :   GEN T = RgX_to_FpX(pol, p);
    3591           7 :   if (signe(T)==0) pari_err_OP("^",x,gen_m1);
    3592           7 :   if (lgefint(p) == 3)
    3593             :   {
    3594           7 :     ulong pp = uel(p, 2);
    3595           7 :     GEN Tp = ZX_to_Flx(T, pp);
    3596           7 :     r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
    3597             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3598             :   }
    3599             :   else
    3600           0 :     r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3601           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3602             : }
    3603             : 
    3604             : #define code(t1,t2) ((t1 << 6) | t2)
    3605             : static GEN
    3606      689254 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
    3607             : {
    3608             :   GEN p, pol;
    3609             :   long pa;
    3610      689254 :   long t = RgX_type3(x,y,T, &p,&pol,&pa);
    3611      689261 :   switch(t)
    3612             :   {
    3613      367454 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
    3614       48515 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
    3615          14 :     case t_FFELT:  return FFXQ_mul(x, y, T, pol);
    3616      262591 :     case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
    3617             :     case code(t_POLMOD, t_INTMOD):
    3618           7 :                    return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
    3619       10680 :     default:       return NULL;
    3620             :   }
    3621             : }
    3622             : 
    3623             : GEN
    3624      689254 : RgXQ_mul(GEN x, GEN y, GEN T)
    3625             : {
    3626      689254 :   GEN z = RgXQ_mul_fast(x, y, T);
    3627      689231 :   if (!z) z = RgX_rem(RgX_mul(x, y), T);
    3628      689231 :   return z;
    3629             : }
    3630             : 
    3631             : static GEN
    3632      241199 : RgXQ_sqr_fast(GEN x, GEN T)
    3633             : {
    3634             :   GEN p, pol;
    3635             :   long pa;
    3636      241199 :   long t = RgX_type2(x, T, &p,&pol,&pa);
    3637      241222 :   switch(t)
    3638             :   {
    3639      228481 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
    3640       11564 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
    3641           0 :     case t_FFELT:  return FFXQ_sqr(x, T, pol);
    3642          14 :     case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
    3643             :     case code(t_POLMOD, t_INTMOD):
    3644           0 :                    return RgXQ_sqr_FpXQXQ(x, T, pol, p);
    3645        1163 :     default:       return NULL;
    3646             :   }
    3647             : }
    3648             : 
    3649             : GEN
    3650      241199 : RgXQ_sqr(GEN x, GEN T)
    3651             : {
    3652      241199 :   GEN z = RgXQ_sqr_fast(x, T);
    3653      241216 :   if (!z) z = RgX_rem(RgX_sqr(x), T);
    3654      241216 :   return z;
    3655             : }
    3656             : 
    3657             : static GEN
    3658       32174 : RgXQ_inv_fast(GEN x, GEN y)
    3659             : {
    3660             :   GEN p, pol;
    3661             :   long pa;
    3662       32174 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3663       32174 :   switch(t)
    3664             :   {
    3665       18980 :     case t_INT:    return QXQ_inv(x,y);
    3666         742 :     case t_FRAC:   return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
    3667          14 :     case t_FFELT:  return FFXQ_inv(x, y, pol);
    3668       11571 :     case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
    3669             :     case code(t_POLMOD, t_INTMOD):
    3670           7 :                    return RgXQ_inv_FpXQXQ(x, y, pol, p);
    3671         860 :     default:       return NULL;
    3672             :   }
    3673             : }
    3674             : #undef code
    3675             : 
    3676             : GEN
    3677       32174 : RgXQ_inv(GEN x, GEN y)
    3678             : {
    3679       32174 :   GEN z = RgXQ_inv_fast(x, y);
    3680       32160 :   if (!z) z = RgXQ_inv_i(x, y);
    3681       32160 :   return z;
    3682             : }

Generated by: LCOV version 1.13