Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - polarit2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24988-2584e74448) Lines: 1957 2160 90.6 %
Date: 2020-01-26 05:57:03 Functions: 192 201 95.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. 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       40929 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      30             : {
      31       40929 :   long dP=degpol(P), i, k, m;
      32             :   pari_sp av1, av2;
      33             :   GEN s,y,P_lead;
      34             : 
      35       40929 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      36       40929 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      37       40929 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      38       40929 :   y = cgetg(n+2,t_COL);
      39       40929 :   if (y0)
      40             :   {
      41       11893 :     if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
      42       11893 :     m = lg(y0)-1;
      43       11893 :     for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
      44             :   }
      45             :   else
      46             :   {
      47       29036 :     m = 1;
      48       29036 :     gel(y,1) = stoi(dP);
      49             :   }
      50       40929 :   P += 2; /* strip codewords */
      51             : 
      52       40929 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      53       40929 :   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      109851 :   for (k=m; k<=n; k++)
      59             :   {
      60       68922 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      61      305088 :     for (i=1; i<k && i<=dP; i++)
      62      236166 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      63       68922 :     if (N)
      64             :     {
      65       16954 :       s = Fq_red(s, T, N);
      66       16954 :       if (P_lead) s = Fq_mul(s, P_lead, T, N);
      67             :     }
      68       51968 :     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       51968 :       if (P_lead) s = gdiv(s, P_lead);
      75       68922 :     av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
      76             :   }
      77       40929 :   return y;
      78             : }
      79             : 
      80             : GEN
      81       25319 : polsym(GEN x, long n)
      82             : {
      83       25319 :   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    44614721 : centermodii(GEN x, GEN p, GEN po2)
      89             : {
      90    44614721 :   GEN y = remii(x, p);
      91    44600839 :   switch(signe(y))
      92             :   {
      93     3107984 :     case 0: break;
      94    26998430 :     case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
      95    26870151 :       break;
      96    14824893 :     case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
      97    14763368 :       break;
      98             :   }
      99    44411035 :   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     5987830 : centermod_i(GEN x, GEN p, GEN ps2)
     113             : {
     114             :   long i, lx;
     115             :   pari_sp av;
     116             :   GEN y;
     117             : 
     118     5987830 :   if (!ps2) ps2 = shifti(p,-1);
     119     5987820 :   switch(typ(x))
     120             :   {
     121     1357457 :     case t_INT: return centermodii(x,p,ps2);
     122             : 
     123     2840910 :     case t_POL: lx = lg(x);
     124     2840910 :       y = cgetg(lx,t_POL); y[1] = x[1];
     125    21454540 :       for (i=2; i<lx; i++)
     126             :       {
     127    18613474 :         av = avma;
     128    18613474 :         gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
     129             :       }
     130     2841066 :       return normalizepol_lg(y, lx);
     131             : 
     132     1788494 :     case t_COL: lx = lg(x);
     133     1788494 :       y = cgetg(lx,t_COL);
     134     1788630 :       for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
     135     1788547 :       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     3422477 : 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         161 : RgXY_squff(GEN f)
     180             : {
     181         161 :   long i, q, n = degpol(f);
     182         161 :   ulong p = itos_or_0(characteristic(f));
     183         161 :   GEN u = const_vec(n+1, pol_1(varn(f)));
     184         175 :   for(q = 1;;q *= p)
     185          14 :   {
     186         175 :     GEN t, v, tv, r = RgX_gcd(f, RgX_deriv(f));
     187         175 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
     188          42 :     t = RgX_div(f, r);
     189          42 :     if (degpol(t) > 0)
     190             :     {
     191             :       long j;
     192          35 :       for(j = 1;;j++)
     193             :       {
     194          56 :         v = RgX_gcd(r, t);
     195          35 :         tv = RgX_div(t, v);
     196          35 :         if (degpol(tv) > 0) gel(u, j*q) = tv;
     197          35 :         if (degpol(v) <= 0) break;
     198          21 :         r = RgX_div(r, v);
     199          21 :         t = v;
     200             :       }
     201          14 :       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         686 :   for (i = n; i; i--)
     209         686 :     if (degpol(gel(u,i))) break;
     210         161 :   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        8862 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
     222             : {
     223             :   GEN q;
     224        8862 :   if (i == lg(Lmod)) return 0;
     225        4648 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
     226        4410 :   if (!gel(Lmod,i)) return 0;
     227        4354 :   p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
     228        4354 :   q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
     229        4354 :   if (degpol(q))
     230             :   {
     231        3990 :     GEN R, Q = RgX_divrem(*F, q, &R);
     232        3990 :     if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
     233             :   }
     234        4046 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
     235        3780 :   return 0;
     236             : }
     237             : 
     238             : static GEN factor_domain(GEN x, GEN flag);
     239             : 
     240             : static GEN
     241         273 : ok_bloc(GEN f, GEN BLOC, ulong c)
     242             : {
     243         273 :   GEN F = poleval(f, BLOC);
     244         273 :   return issquarefree(c ? gmul(F,mkintmodu(1,c)): F)? F: NULL;
     245             : }
     246             : static GEN
     247         175 : RgXY_factor_squarefree(GEN f, GEN dom)
     248             : {
     249         175 :   pari_sp av = avma;
     250         175 :   ulong i, c = itou_or_0(residual_characteristic(f));
     251         175 :   long vy = gvar2(f), val = RgX_valrem(f, &f), n = RgXY_degreex(f);
     252         175 :   GEN y, Lmod, F = NULL, BLOC = NULL, Lfac = coltrunc_init(degpol(f)+2);
     253         175 :   if (val)
     254             :   {
     255          28 :     GEN x = pol_x(varn(f));
     256          28 :     if (dom)
     257             :     {
     258          14 :       GEN c = Rg_get_1(dom);
     259          14 :       if (typ(c) != t_INT) x = RgX_Rg_mul(x, c);
     260             :     }
     261          28 :     vectrunc_append(Lfac, x); if (!degpol(f)) return Lfac;
     262             :   }
     263         168 :   y = pol_x(vy);
     264             :   for(;;)
     265             :   {
     266         231 :     for (i = 0; !c || i < c; i++)
     267             :     {
     268         231 :       BLOC = gpowgs(gaddgs(y, i), n+1);
     269         231 :       if ((F = ok_bloc(f, BLOC, c))) break;
     270          77 :       if (c)
     271             :       {
     272          42 :         BLOC = random_FpX(n+1, vy, utoipos(c));
     273          42 :         gel(BLOC,lg(BLOC)-1) = gen_1;
     274          42 :         if ((F = ok_bloc(f, BLOC, c))) break;
     275             :       }
     276             :     }
     277         168 :     if (!c || i < c) break;
     278           0 :     n++;
     279             :   }
     280         168 :   if (DEBUGLEVEL >= 2)
     281           0 :     err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
     282         168 :   Lmod = gel(factor_domain(F,dom),1);
     283         168 :   if (DEBUGLEVEL >= 2)
     284           0 :     err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
     285         168 :   (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
     286         168 :   if (degpol(f)) vectrunc_append(Lfac, f);
     287         168 :   return gerepilecopy(av, Lfac);
     288             : }
     289             : 
     290             : static GEN
     291         161 : FE_matconcat(GEN F, GEN E, long l)
     292             : {
     293         161 :   setlg(E,l); E = shallowconcat1(E);
     294         161 :   setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
     295             : }
     296             : 
     297             : static int
     298         308 : gen_cmp_RgXY(void *data, GEN x, GEN y)
     299             : {
     300         308 :   long vx = varn(x), vy = varn(y);
     301         308 :   return (vx == vy)? gen_cmp_RgX(data, x, y): -varncmp(vx, vy);
     302             : }
     303             : static GEN
     304         161 : RgXY_factor(GEN f, GEN dom)
     305             : {
     306         161 :   pari_sp av = avma;
     307             :   GEN C, F, E, cf, V;
     308             :   long i, j, l;
     309         161 :   if (dom) { GEN c = Rg_get_1(dom); if (typ(c) != t_INT) f = RgX_Rg_mul(f,c); }
     310         161 :   cf = content(f);
     311         161 :   V = RgXY_squff(gdiv(f, cf)); l = lg(V);
     312         161 :   C = factor_domain(cf, dom);
     313         161 :   F = cgetg(l+1, t_VEC); gel(F,1) = gel(C,1);
     314         161 :   E = cgetg(l+1, t_VEC); gel(E,1) = gel(C,2);
     315         427 :   for (i=1, j=2; i < l; i++)
     316             :   {
     317         266 :     GEN v = gel(V,i);
     318         266 :     if (degpol(v))
     319             :     {
     320         175 :       gel(F,j) = v = RgXY_factor_squarefree(v, dom);
     321         175 :       gel(E,j) = const_col(lg(v)-1, utoipos(i));
     322         175 :       j++;
     323             :     }
     324             :   }
     325         161 :   f = FE_matconcat(F,E,j);
     326         161 :   (void)sort_factor(f,(void*)cmp_universal, &gen_cmp_RgXY);
     327         161 :   return gerepilecopy(av, f);
     328             : }
     329             : 
     330             : /***********************************************************************/
     331             : /**                                                                   **/
     332             : /**                          FACTORIZATION                            **/
     333             : /**                                                                   **/
     334             : /***********************************************************************/
     335             : static long RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var);
     336             : #define assign_or_fail(x,y) { GEN __x = x;\
     337             :   if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
     338             : }
     339             : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
     340             : 
     341             : static const long tsh = 6;
     342             : #define code(t1,t2) ((t1 << 6) | t2)
     343             : void
     344       75782 : RgX_type_decode(long x, long *t1, long *t2)
     345             : {
     346       75782 :   *t1 = x >> tsh;
     347       75782 :   *t2 = (x & ((1L<<tsh)-1));
     348       75782 : }
     349             : int
     350    10404699 : RgX_type_is_composite(long t) { return t >= tsh; }
     351             : 
     352             : static int
     353  1565427873 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     354             : {
     355             :   long j;
     356  1565427873 :   switch(typ(c))
     357             :   {
     358             :     case t_INT:
     359  1352627830 :       break;
     360             :     case t_FRAC:
     361    26495641 :       t[1]=1; break;
     362             :       break;
     363             :     case t_REAL:
     364    77262072 :       update_prec(precision(c), pa);
     365    77262072 :       t[2]=1; break;
     366             :     case t_INTMOD:
     367    23158459 :       assign_or_fail(gel(c,1),p);
     368    23158459 :       t[3]=1; break;
     369             :     case t_FFELT:
     370     1868810 :       if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
     371     1868810 :       assign_or_fail(FF_p_i(c),p);
     372     1868810 :       t[5]=1; break;
     373             :     case t_COMPLEX:
     374   161629294 :       for (j=1; j<=2; j++)
     375             :       {
     376   107752865 :         GEN d = gel(c,j);
     377   107752865 :         switch(typ(d))
     378             :         {
     379             :           case t_INT: case t_FRAC:
     380     2409743 :             if (!*t2) *t2 = t_COMPLEX;
     381     2409743 :             t[1]=1; break;
     382             :           case t_REAL:
     383   105343101 :             update_prec(precision(d), pa);
     384   105343101 :             if (!*t2) *t2 = t_COMPLEX;
     385   105343101 :             t[2]=1; break;
     386             :           case t_INTMOD:
     387          14 :             assign_or_fail(gel(d,1),p);
     388          14 :             if (!signe(*p) || mod4(*p) != 3) return 0;
     389           7 :             if (!*t2) *t2 = t_COMPLEX;
     390           7 :             t[3]=1; break;
     391             :           case t_PADIC:
     392           7 :             update_prec(precp(d)+valp(d), pa);
     393           7 :             assign_or_fail(gel(d,2),p);
     394           7 :             if (!*t2) *t2 = t_COMPLEX;
     395           7 :             t[7]=1; break;
     396           0 :           default: return 0;
     397             :         }
     398             :       }
     399    53876429 :       if (!t[2]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
     400    53876429 :       break;
     401             :     case t_PADIC:
     402     1214142 :       update_prec(precp(c)+valp(c), pa);
     403     1214142 :       assign_or_fail(gel(c,2),p);
     404     1214142 :       t[7]=1; break;
     405             :     case t_QUAD:
     406        1512 :       assign_or_fail(gel(c,1),pol);
     407        4536 :       for (j=2; j<=3; j++)
     408             :       {
     409        3024 :         GEN d = gel(c,j);
     410        3024 :         switch(typ(d))
     411             :         {
     412             :           case t_INT: case t_FRAC:
     413        2989 :             t[8]=1; break;
     414             :           case t_INTMOD:
     415          28 :             assign_or_fail(gel(d,1),p);
     416          28 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     417          28 :             t[3]=1; break;
     418             :           case t_PADIC:
     419           7 :             update_prec(precp(d)+valp(d), pa);
     420           7 :             assign_or_fail(gel(d,2),p);
     421           7 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     422           7 :             t[7]=1; break;
     423           0 :           default: return 0;
     424             :         }
     425             :       }
     426        1512 :       break;
     427             :     case t_POLMOD:
     428     4512527 :       assign_or_fail(gel(c,1),pol);
     429     4512340 :       if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1))) return 0;
     430    27055778 :       for (j=1; j<=2; j++)
     431             :       {
     432             :         GEN pbis, polbis;
     433             :         long pabis;
     434     9020790 :         *t2 = t_POLMOD;
     435     9020790 :         switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
     436             :         {
     437     5682233 :           case t_INT: break;
     438      756441 :           case t_FRAC:   t[1]=1; break;
     439     2576875 :           case t_INTMOD: t[3]=1; break;
     440           7 :           case t_PADIC:  t[7]=1; update_prec(pabis,pa); break;
     441       10500 :           default: return 0;
     442             :         }
     443     9015556 :         if (pbis) assign_or_fail(pbis,p);
     444     9015556 :         if (polbis) assign_or_fail(polbis,pol);
     445             :       }
     446     4507099 :       break;
     447     5320773 :     case t_RFRAC: t[10] = 1;
     448     5320773 :       if (!settype(gel(c,1),t,p,pol,pa,ff,t2,var)) return 0;
     449     5320773 :       c = gel(c,2); /* fall through */
     450    24409044 :     case t_POL: t[10] = 1;
     451    24409044 :       if (!RgX_settype(c,t,p,pol,pa,ff,t2,var)) return 0;
     452    24409031 :       if (*var == NO_VARIABLE) { *var = varn(c); break; }
     453             :       /* if more than one free var, ensure varn() == *var fails. FIXME: should
     454             :        * keep the list of all variables, later t_POLMOD may cancel them */
     455    16881101 :       if (*var != varn(c)) *var = MAXVARN+1;
     456    16881101 :       break;
     457        1400 :     default: return 0;
     458             :   }
     459  1565421025 :   return 1;
     460             : }
     461             : /* t[0] unused. Other values, if set, indicate a coefficient of type
     462             :  * t[1] : t_FRAC
     463             :  * t[2] : t_REAL
     464             :  * t[3] : t_INTMOD
     465             :  * t[4] : Unused
     466             :  * t[5] : t_FFELT
     467             :  * t[6] : Unused
     468             :  * t[7] : t_PADIC
     469             :  * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
     470             :  * t[9]:  Unused
     471             :  * t[10]: t_POL (recursive factorisation) */
     472             : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
     473             :  * given by t) */
     474             : static long
     475   108059836 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
     476             : {
     477   108059836 :   if (t[10] && (!*pol || var!=varn(*pol))) return t_POL;
     478   100542980 :   if (t2) /* polmod/quad/complex of intmod/padic */
     479             :   {
     480     3448424 :     if (t[2] && (t[3]||t[7])) return 0;
     481     3448424 :     if (t[3]) return code(t2,t_INTMOD);
     482     3419171 :     if (t[7]) return code(t2,t_PADIC);
     483     3419136 :     if (t[2]) return t_COMPLEX;
     484      420704 :     if (t[1]) return code(t2,t_FRAC);
     485      246228 :     return code(t2,t_INT);
     486             :   }
     487    97094556 :   if (t[5]) /* ffelt */
     488             :   {
     489      219682 :     if (t[2]||t[8]||t[9]) return 0;
     490      219682 :     *pol=ff; return t_FFELT;
     491             :   }
     492    96874874 :   if (t[2]) /* inexact, real */
     493             :   {
     494    12894208 :     if (t[3]||t[7]||t[9]) return 0;
     495    12894208 :     return t_REAL;
     496             :   }
     497    83980666 :   if (t[10]) return t_POL;
     498    83980666 :   if (t[8]) return code(t_QUAD,t_INT);
     499    83980022 :   if (t[3]) return t_INTMOD;
     500    80290337 :   if (t[7]) return t_PADIC;
     501    79953045 :   if (t[1]) return t_FRAC;
     502    74582898 :   return t_INT;
     503             : }
     504             : 
     505             : static long
     506   188115084 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     507             : {
     508   188115084 :   long i, lx = lg(x);
     509   762277019 :   for (i=2; i<lx; i++)
     510   574166021 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     511   188110998 :   return 1;
     512             : }
     513             : 
     514             : static long
     515   116408886 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     516             : {
     517   116408886 :   long i, l = lg(x);
     518  1101650161 :   for (i = 1; i<l; i++)
     519   985244040 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     520   116406121 :   return 1;
     521             : }
     522             : 
     523             : static long
     524    21102430 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     525             : {
     526    21102430 :   long i, l = lg(x);
     527   125627058 :   for (i = 1; i < l; i++)
     528   104526609 :     if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     529    21100449 :   return 1;
     530             : }
     531             : 
     532             : long
     533    19425517 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
     534             : {
     535    19425517 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0};
     536    19425517 :   long t2 = 0, var = NO_VARIABLE;
     537    19425517 :   GEN ff = NULL;
     538    19425517 :   *p = *pol = NULL; *pa = LONG_MAX;
     539    19425517 :   switch(typ(x))
     540             :   {
     541             :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
     542             :     case t_COMPLEX: case t_PADIC: case t_QUAD:
     543      696770 :       if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     544      696770 :       break;
     545             :     case t_POL: case t_SER:
     546    18164547 :       if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     547    18164550 :       break;
     548             :     case t_VEC: case t_COL:
     549          14 :       if(!RgC_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
     550          14 :       break;
     551             :     case t_MAT:
     552         126 :       if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
     553         126 :       break;
     554      564060 :     default: return 0;
     555             :   }
     556    18861460 :   return choosetype(t,t2,ff,pol,var);
     557             : }
     558             : 
     559             : long
     560      328773 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     561             : {
     562      328773 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     563      328773 :   long t2 = 0, var = NO_VARIABLE;
     564      328773 :   GEN ff = NULL;
     565      328773 :   *p = *pol = NULL; *pa = LONG_MAX;
     566      328773 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     567      328726 :   return choosetype(t,t2,ff,pol,var);
     568             : }
     569             : 
     570             : long
     571         294 : RgX_Rg_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     572             : {
     573         294 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     574         294 :   long t2 = 0, var = NO_VARIABLE;
     575         294 :   GEN ff = NULL;
     576         294 :   *p = *pol = NULL; *pa = LONG_MAX;
     577         294 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     578         294 :   if (!settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     579         294 :   return choosetype(t,t2,ff,pol,var);
     580             : }
     581             : 
     582             : long
     583    71737117 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     584             : {
     585    71737117 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     586    71737117 :   long t2 = 0, var = NO_VARIABLE;
     587    71737117 :   GEN ff = NULL;
     588    71737117 :   *p = *pol = NULL; *pa = LONG_MAX;
     589   143470313 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     590    71737228 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     591    71733084 :   return choosetype(t,t2,ff,pol,var);
     592             : }
     593             : 
     594             : long
     595      580718 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
     596             : {
     597      580718 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     598      580718 :   long t2 = 0, var = NO_VARIABLE;
     599      580718 :   GEN ff = NULL;
     600      580718 :   *p = *pol = NULL; *pa = LONG_MAX;
     601     1161434 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     602     1161434 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
     603      580716 :       !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
     604      580717 :   return choosetype(t,t2,ff,pol,var);
     605             : }
     606             : 
     607             : long
     608      131107 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
     609             : {
     610      131107 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     611      131107 :   long t2 = 0, var = NO_VARIABLE;
     612      131107 :   GEN ff = NULL;
     613      131107 :   *p = *pol = NULL; *pa = LONG_MAX;
     614      131107 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     615      130141 :   return choosetype(t,t2,ff,pol,var);
     616             : }
     617             : 
     618             : long
     619          70 : RgV_type(GEN x, GEN *p, GEN *pol, long *pa)
     620             : {
     621          70 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     622          70 :   long t2 = 0, var = NO_VARIABLE;
     623          70 :   GEN ff = NULL;
     624          70 :   *p = *pol = NULL; *pa = LONG_MAX;
     625          70 :   if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     626          70 :   return choosetype(t,t2,ff,pol,var);
     627             : }
     628             : 
     629             : long
     630         196 : RgV_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     631             : {
     632         196 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     633         196 :   long t2 = 0, var = NO_VARIABLE;
     634         196 :   GEN ff = NULL;
     635         196 :   *p = *pol = NULL; *pa = LONG_MAX;
     636         392 :   if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     637         196 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     638         196 :   return choosetype(t,t2,ff,pol,var);
     639             : }
     640             : 
     641             : long
     642    11882235 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     643             : {
     644    11882235 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     645    11882235 :   long t2 = 0, var = NO_VARIABLE;
     646    11882235 :   GEN ff = NULL;
     647    11882235 :   *p = *pol = NULL; *pa = LONG_MAX;
     648    23764036 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     649    11883019 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     650    11881017 :   return choosetype(t,t2,ff,pol,var);
     651             : }
     652             : 
     653             : long
     654     4544719 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     655             : {
     656     4544719 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     657     4544719 :   long t2 = 0, var = NO_VARIABLE;
     658     4544719 :   GEN ff = NULL;
     659     4544719 :   *p = *pol = NULL; *pa = LONG_MAX;
     660     9088962 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     661     4544824 :       !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     662     4544138 :   return choosetype(t,t2,ff,pol,var);
     663             : }
     664             : 
     665             : GEN
     666       37448 : factor0(GEN x, GEN flag)
     667             : {
     668             :   ulong B;
     669       37448 :   long tx = typ(x);
     670       37448 :   if (!flag) return factor(x);
     671         217 :   if ((tx != t_INT && tx!=t_FRAC) || typ(flag) != t_INT)
     672         175 :     return factor_domain(x, flag);
     673          42 :   if (signe(flag) < 0) pari_err_FLAG("factor");
     674          42 :   switch(lgefint(flag))
     675             :   {
     676           7 :     case 2: B = 0; break;
     677          35 :     case 3: B = flag[2]; break;
     678           0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     679             :              return NULL; /*LCOV_EXCL_LINE*/
     680             :   }
     681          42 :   return boundfact(x, B);
     682             : }
     683             : 
     684             : GEN
     685       91523 : deg1_from_roots(GEN L, long v)
     686             : {
     687       91523 :   long i, l = lg(L);
     688       91523 :   GEN z = cgetg(l,t_COL);
     689      211873 :   for (i=1; i<l; i++)
     690      120350 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     691       91523 :   return z;
     692             : }
     693             : GEN
     694        1078 : roots_from_deg1(GEN x)
     695             : {
     696        1078 :   long i,l = lg(x);
     697        1078 :   GEN r = cgetg(l,t_VEC);
     698        1078 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     699        1078 :   return r;
     700             : }
     701             : 
     702             : static GEN
     703          42 : gauss_factor_p(GEN p)
     704             : {
     705          42 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     706          42 :   return mkcomplex(a, b);
     707             : }
     708             : 
     709             : static GEN
     710          49 : gauss_primpart(GEN x, GEN *c)
     711             : {
     712          49 :   GEN a = real_i(x), b = imag_i(x), n = gcdii(a, b);
     713          49 :   *c = n; if (n == gen_1) return x;
     714          49 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     715             : }
     716             : 
     717             : static GEN
     718          70 : gauss_primpart_try(GEN x, GEN c)
     719             : {
     720             :   GEN r, y;
     721          70 :   if (typ(x) == t_INT)
     722             :   {
     723          42 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     724             :   }
     725             :   else
     726             :   {
     727          28 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     728          28 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     729          14 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     730             :   }
     731          56 :   return y;
     732             : }
     733             : 
     734             : static int
     735          91 : gauss_cmp(GEN x, GEN y)
     736             : {
     737             :   int v;
     738          91 :   if (typ(x) != t_COMPLEX)
     739           0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     740          91 :   if (typ(y) != t_COMPLEX) return 1;
     741          63 :   v = cmpii(gel(x,2), gel(y,2));
     742          63 :   if (v) return v;
     743          28 :   return gcmp(gel(x,1), gel(y,1));
     744             : }
     745             : 
     746             : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
     747             : static GEN
     748         455 : gauss_normal(GEN x)
     749             : {
     750         455 :   if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
     751         385 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     752         385 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     753         385 :   return x;
     754             : }
     755             : 
     756             : static GEN
     757          49 : gauss_factor(GEN x)
     758             : {
     759          49 :   pari_sp av = avma;
     760          49 :   GEN a = real_i(x), b = imag_i(x), d = gen_1, n, y, fa, P, E, P2, E2;
     761          49 :   long t1 = typ(a);
     762          49 :   long t2 = typ(b), i, j, l, exp = 0;
     763          49 :   if (t1 == t_FRAC) d = gel(a,2);
     764          49 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     765          49 :   if (d == gen_1) y = x;
     766             :   else
     767             :   {
     768          21 :     y = gmul(x, d);
     769          21 :     a = real_i(y); t1 = typ(a);
     770          21 :     b = imag_i(y); t2 = typ(b);
     771             :   }
     772          49 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     773          49 :   y = gauss_primpart(y, &n);
     774          49 :   fa = factor(cxnorm(y));
     775          49 :   P = gel(fa,1);
     776          49 :   E = gel(fa,2); l = lg(P);
     777          49 :   P2 = cgetg(l, t_COL);
     778          49 :   E2 = cgetg(l, t_COL);
     779         105 :   for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
     780             :   { /* either p = 2 (ramified) or those factors split in Q(i) */
     781          56 :     GEN p = gel(P,i), w, w2, t, we, pe;
     782          56 :     long v, e = itos(gel(E,i));
     783          56 :     int is2 = absequaliu(p, 2);
     784          56 :     w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     785          56 :     w2 = gauss_normal( conj_i(w) );
     786             :     /* w * w2 * I^3 = p, w2 = conj(w) * I */
     787          56 :     pe = powiu(p, e);
     788          56 :     we = gpowgs(w, e);
     789          56 :     t = gauss_primpart_try( gmul(y, conj_i(we)), pe );
     790          56 :     if (t) y = t; /* y /= w^e */
     791             :     else {
     792             :       /* y /= conj(w)^e, should be y /= w2^e */
     793          14 :       y = gauss_primpart_try( gmul(y, we), pe );
     794          14 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     795             :     }
     796          56 :     gel(P,i) = w;
     797          56 :     v = Z_pvalrem(n, p, &n);
     798          56 :     if (v) {
     799           7 :       exp -= v; /* += 3*v mod 4 */
     800           7 :       if (is2) v <<= 1; /* 2 = w^2 I^3 */
     801             :       else {
     802           0 :         gel(P2,j) = w2;
     803           0 :         gel(E2,j) = utoipos(v); j++;
     804             :       }
     805           7 :       gel(E,i) = stoi(e + v);
     806             :     }
     807          56 :     v = Z_pvalrem(d, p, &d);
     808          56 :     if (v) {
     809           7 :       exp += v; /* -= 3*v mod 4 */
     810           7 :       if (is2) v <<= 1; /* 2 is ramified */
     811             :       else {
     812           7 :         gel(P2,j) = w2;
     813           7 :         gel(E2,j) = utoineg(v); j++;
     814             :       }
     815           7 :       gel(E,i) = stoi(e - v);
     816             :     }
     817          56 :     exp &= 3;
     818             :   }
     819          49 :   if (j > 1) {
     820           7 :     long k = 1;
     821           7 :     GEN P1 = cgetg(l, t_COL);
     822           7 :     GEN E1 = cgetg(l, t_COL);
     823             :     /* remove factors with exponent 0 */
     824          14 :     for (i = 1; i < l; i++)
     825           7 :       if (signe(gel(E,i)))
     826             :       {
     827           0 :         gel(P1,k) = gel(P,i);
     828           0 :         gel(E1,k) = gel(E,i);
     829           0 :         k++;
     830             :       }
     831           7 :     setlg(P1, k); setlg(E1, k);
     832           7 :     setlg(P2, j); setlg(E2, j);
     833           7 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     834             :   }
     835          49 :   if (!equali1(n) || !equali1(d))
     836             :   {
     837          28 :     GEN Fa = factor(Qdivii(n, d));
     838          28 :     P = gel(Fa,1); l = lg(P);
     839          28 :     E = gel(Fa,2);
     840          70 :     for (i = 1; i < l; i++)
     841             :     {
     842          42 :       GEN w, p = gel(P,i);
     843             :       long e;
     844             :       int is2;
     845          42 :       switch(mod4(p))
     846             :       {
     847          14 :         case 3: continue;
     848          14 :         case 2: is2 = 1; break;
     849          14 :         default:is2 = 0; break;
     850             :       }
     851          28 :       e = itos(gel(E,i));
     852          28 :       w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     853          28 :       gel(P,i) = w;
     854          28 :       if (is2)
     855          14 :         gel(E,i) = stoi(2*e);
     856             :       else
     857             :       {
     858          14 :         P = shallowconcat(P, gauss_normal( conj_i(w) ));
     859          14 :         E = shallowconcat(E, gel(E,i));
     860             :       }
     861          28 :       exp -= e; /* += 3*e mod 4 */
     862          28 :       exp &= 3;
     863             :     }
     864          28 :     gel(Fa,1) = P;
     865          28 :     gel(Fa,2) = E;
     866          28 :     fa = famat_mul_shallow(fa, Fa);
     867             :   }
     868          49 :   fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
     869             : 
     870          49 :   y = gmul(y, powIs(exp));
     871          49 :   if (!gequal1(y)) {
     872          35 :     gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
     873          35 :     gel(fa,2) = shallowconcat(gen_1,    gel(fa,2));
     874             :   }
     875          49 :   return gerepilecopy(av, fa);
     876             : }
     877             : 
     878             : GEN
     879          35 : Q_factor_limit(GEN x, ulong lim)
     880             : {
     881          35 :   pari_sp av = avma;
     882             :   GEN a, b;
     883          35 :   if (typ(x) == t_INT) return Z_factor_limit(x, lim);
     884          14 :   a = Z_factor_limit(gel(x,1), lim);
     885          14 :   b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
     886          14 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     887             : }
     888             : GEN
     889        4249 : Q_factor(GEN x)
     890             : {
     891        4249 :   pari_sp av = avma;
     892             :   GEN a, b;
     893        4249 :   if (typ(x) == t_INT) return Z_factor(x);
     894         497 :   a = Z_factor(gel(x,1));
     895         497 :   b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
     896         497 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     897             : }
     898             : /* replace t_QUAD/t_COMPLEX coeffs by t_POLMOD in T */
     899             : static GEN
     900         126 : RgX_fix_quad(GEN x, GEN T)
     901             : {
     902         126 :   long i, l, v = varn(T);
     903         126 :   GEN y = cgetg_copy(x,&l);
     904         630 :   for (i = 2; i < l; i++)
     905             :   {
     906         504 :     GEN c = gel(x,i);
     907         504 :     switch(typ(c))
     908             :     {
     909          56 :       case t_QUAD: c++;/* fall through */
     910          98 :       case t_COMPLEX: c = deg1pol_shallow(gel(c,2),gel(c,1),v);
     911             :     }
     912         504 :     gel(y,i) = c;
     913             :   }
     914         126 :   y[1] = x[1]; return y;
     915             : }
     916             : 
     917             : static GEN
     918        2345 : RgX_factor(GEN x, GEN dom)
     919             : {
     920             :   pari_sp av;
     921             :   long pa, v, lx, r1, i;
     922             :   GEN  p, pol, y, p1, p2;
     923        2345 :   long tx = dom ? RgX_Rg_type(x,dom,&p,&pol,&pa): RgX_type(x,&p,&pol,&pa);
     924        2345 :   switch(tx)
     925             :   {
     926           7 :     case 0: pari_err_IMPL("factor for general polynomials");
     927         161 :     case t_POL: return RgXY_factor(x, dom);
     928        1554 :     case t_INT: return ZX_factor(x);
     929           7 :     case t_FRAC: return QX_factor(x);
     930         175 :     case t_INTMOD: return factmod(x,p);
     931          42 :     case t_PADIC: return factorpadic(x,p,pa);
     932          98 :     case t_FFELT: return FFX_factor(x,pol);
     933             : 
     934          21 :     case t_COMPLEX: y = cgetg(3,t_MAT);
     935          21 :       av = avma; p1 = deg1_from_roots(roots(x,pa), varn(x));
     936          21 :       gel(y,1) = p1 = gerepileupto(av, p1);
     937          21 :       gel(y,2) = const_col(lg(p1)-1, gen_1); return y;
     938             : 
     939          28 :     case t_REAL: y=cgetg(3,t_MAT); v=varn(x);
     940          28 :       av=avma; p1=cleanroots(x,pa);
     941          28 :       lx = lg(p1);
     942          70 :       for (r1 = 1; r1 < lx; r1++)
     943          49 :         if (typ(gel(p1,r1)) == t_COMPLEX) break;
     944          28 :       lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
     945          70 :       for (i = 1; i < r1; i++)
     946          42 :         gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
     947          35 :       for (   ; i < lx; i++)
     948             :       {
     949           7 :         GEN a = gel(p1,2*i-r1);
     950           7 :         p = cgetg(5, t_POL); gel(p2,i) = p;
     951           7 :         p[1] = x[1];
     952           7 :         gel(p,2) = gnorm(a);
     953           7 :         gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
     954           7 :         gel(p,4) = gen_1;
     955             :       }
     956          28 :       gel(y,1) = gerepileupto(av,p2);
     957          28 :       gel(y,2) = const_col(lx-1, gen_1); return y;
     958             : 
     959             :     default:
     960             :     {
     961         252 :       GEN w = NULL, T = pol;
     962             :       long t1, t2;
     963         252 :       av = avma;
     964         252 :       RgX_type_decode(tx, &t1, &t2);
     965         252 :       if (t1 == t_COMPLEX) w = gen_I();
     966         203 :       else if (t1 == t_QUAD) w = mkquad(pol,gen_0,gen_1);
     967         252 :       if (w)
     968             :       { /* substitute I or w by t_POLMOD */
     969         126 :         T = leafcopy(pol); setvarn(T, fetch_var());
     970         126 :         x = RgX_fix_quad(x, T);
     971             :       }
     972         252 :       switch (t2)
     973             :       {
     974         161 :         case t_INT: case t_FRAC: p1 = nffactor(T,x); break;
     975             :         case t_INTMOD:
     976          56 :           T = RgX_to_FpX(T,p);
     977          56 :           if (FpX_is_irred(T,p)) { p1 = factmod(x,mkvec2(p,T)); break; }
     978             :         /*fall through*/
     979             :         default:
     980          56 :           if (w) (void)delete_var();
     981          56 :           pari_err_IMPL("factor for general polynomial");
     982             :           return NULL; /* LCOV_EXCL_LINE */
     983             :       }
     984         196 :       if (t1 == t_POLMOD) return gerepileupto(av, p1);
     985             :       /* substitute back I or w */
     986          98 :       gel(p1,1) = gsubst(liftpol_shallow(gel(p1,1)), varn(T), w);
     987          98 :       (void)delete_var(); return gerepilecopy(av, p1);
     988             :     }
     989             :   }
     990             : }
     991             : 
     992             : static GEN
     993       41501 : factor_domain(GEN x, GEN dom)
     994             : {
     995       41501 :   long tx = typ(x);
     996       41501 :   long tdom = dom ? typ(dom): 0;
     997             :   pari_sp av;
     998             : 
     999       41501 :   if (gequal0(x))
    1000          63 :     switch(tx)
    1001             :     {
    1002             :       case t_INT:
    1003             :       case t_COMPLEX:
    1004             :       case t_POL:
    1005          63 :       case t_RFRAC: return prime_fact(x);
    1006           0 :       default: pari_err_TYPE("factor",x);
    1007             :     }
    1008       41438 :   av = avma;
    1009       41438 :   switch(tx)
    1010             :   {
    1011        2240 :     case t_POL: return RgX_factor(x, dom);
    1012             :     case t_RFRAC: {
    1013          35 :       GEN a = gel(x,1), b = gel(x,2);
    1014          35 :       GEN y = famat_inv_shallow(RgX_factor(b, dom));
    1015          35 :       if (typ(a)==t_POL) y = famat_mul_shallow(RgX_factor(a, dom), y);
    1016          35 :       return gerepilecopy(av, sort_factor_pol(y, cmp_universal));
    1017             :     }
    1018       38631 :     case t_INT:  if (tdom==0 || tdom==t_INT) return Z_factor(x);
    1019         490 :     case t_FRAC: if (tdom==0 || tdom==t_INT) return Q_factor(x);
    1020             :     case t_COMPLEX: /* fall through */
    1021          49 :       if (tdom==0 || tdom==t_COMPLEX)
    1022             :       {
    1023          49 :         GEN y = gauss_factor(x); if (y) return y;
    1024             :       }
    1025             :       /* fall through */
    1026             :   }
    1027           0 :   pari_err_TYPE("factor",x);
    1028             :   return NULL; /* LCOV_EXCL_LINE */
    1029             : }
    1030             : 
    1031             : GEN
    1032       40997 : factor(GEN x) { return factor_domain(x, NULL); }
    1033             : 
    1034             : /*******************************************************************/
    1035             : /*                                                                 */
    1036             : /*                     ROOTS --> MONIC POLYNOMIAL                  */
    1037             : /*                                                                 */
    1038             : /*******************************************************************/
    1039             : static GEN
    1040      731428 : normalized_mul(void *E, GEN x, GEN y)
    1041             : {
    1042      731428 :   long a = gel(x,1)[1], b = gel(y,1)[1];
    1043             :   (void) E;
    1044     1462856 :   return mkvec2(mkvecsmall(a + b),
    1045     1462856 :                 RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
    1046             : }
    1047             : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
    1048             : static GEN
    1049      382418 : normalized_to_RgX(GEN L)
    1050             : {
    1051      382418 :   long i, a = gel(L,1)[1];
    1052      382418 :   GEN A = gel(L,2);
    1053      382418 :   GEN z = cgetg(a + 3, t_POL);
    1054      382418 :   z[1] = evalsigne(1) | evalvarn(varn(A));
    1055      382418 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
    1056      382418 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
    1057      382418 :   gel(z,i) = gen_1; return z;
    1058             : }
    1059             : 
    1060             : /* compute prod (x - a[i]) */
    1061             : GEN
    1062      324878 : roots_to_pol(GEN a, long v)
    1063             : {
    1064      324878 :   pari_sp av = avma;
    1065      324878 :   long i, k, lx = lg(a);
    1066             :   GEN L;
    1067      324878 :   if (lx == 1) return pol_1(v);
    1068      324830 :   L = cgetg(lx, t_VEC);
    1069      670362 :   for (k=1,i=1; i<lx-1; i+=2)
    1070             :   {
    1071      345532 :     GEN s = gel(a,i), t = gel(a,i+1);
    1072      345532 :     GEN x0 = gmul(s,t);
    1073      345532 :     GEN x1 = gneg(gadd(s,t));
    1074      345532 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1075             :   }
    1076      643388 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
    1077      318558 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
    1078      324830 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1079      324830 :   return gerepileupto(av, normalized_to_RgX(L));
    1080             : }
    1081             : 
    1082             : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
    1083             : GEN
    1084       57588 : roots_to_pol_r1(GEN a, long v, long r1)
    1085             : {
    1086       57588 :   pari_sp av = avma;
    1087       57588 :   long i, k, lx = lg(a);
    1088             :   GEN L;
    1089       57588 :   if (lx == 1) return pol_1(v);
    1090       57588 :   L = cgetg(lx, t_VEC);
    1091      320209 :   for (k=1,i=1; i<r1; i+=2)
    1092             :   {
    1093      262621 :     GEN s = gel(a,i), t = gel(a,i+1);
    1094      262621 :     GEN x0 = gmul(s,t);
    1095      262621 :     GEN x1 = gneg(gadd(s,t));
    1096      262621 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1097             :   }
    1098       78469 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
    1099       20881 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
    1100      223842 :   for (i=r1+1; i<lx; i++)
    1101             :   {
    1102      166254 :     GEN s = gel(a,i);
    1103      166254 :     GEN x0 = gnorm(s);
    1104      166254 :     GEN x1 = gneg(gtrace(s));
    1105      166254 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1106             :   }
    1107       57588 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1108       57588 :   return gerepileupto(av, normalized_to_RgX(L));
    1109             : }
    1110             : 
    1111             : /*******************************************************************/
    1112             : /*                                                                 */
    1113             : /*                          FACTORBACK                             */
    1114             : /*                                                                 */
    1115             : /*******************************************************************/
    1116             : static GEN
    1117          62 : idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
    1118             : static GEN
    1119         419 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
    1120             : static GEN
    1121        4130 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
    1122             : static GEN
    1123        7057 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
    1124             : static GEN
    1125       41928 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
    1126             : static GEN
    1127       57336 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
    1128             : static GEN
    1129    54196910 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
    1130             : static GEN
    1131    78847422 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
    1132             : static GEN
    1133    24179586 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
    1134             : static GEN
    1135      203889 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
    1136             : 
    1137             : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
    1138             : GEN
    1139    30036704 : gen_factorback(GEN L, GEN e, void *data, GEN (*_mul)(void*,GEN,GEN),
    1140             :                GEN (*_pow)(void*,GEN,GEN))
    1141             : {
    1142    30036704 :   pari_sp av = avma;
    1143             :   long k, l, lx;
    1144             :   GEN p,x;
    1145             : 
    1146    30036704 :   if (e) /* supplied vector of exponents */
    1147      929874 :     p = L;
    1148             :   else
    1149             :   {
    1150    29106830 :     switch(typ(L)) {
    1151             :       case t_VEC:
    1152             :       case t_COL: /* product of the L[i] */
    1153     5158112 :         return gerepileupto(av, gen_product(L, data, _mul));
    1154             :       case t_MAT: /* genuine factorization */
    1155    23948718 :         l = lg(L);
    1156    23948718 :         if (l == 3) break;
    1157             :         /*fall through*/
    1158             :       default:
    1159           7 :         pari_err_TYPE("factorback [not a factorization]", L);
    1160             :     }
    1161    23948711 :     p = gel(L,1);
    1162    23948711 :     e = gel(L,2);
    1163             :   }
    1164             :   /* p = elts, e = expo */
    1165    24878585 :   lx = lg(p);
    1166             :   /* check whether e is an integral vector of correct length */
    1167    24878585 :   switch(typ(e))
    1168             :   {
    1169             :     case t_VECSMALL:
    1170        4662 :       if (lx != lg(e))
    1171           0 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1172        4662 :       if (lx == 1) return gen_1;
    1173        3655 :       x = cgetg(lx,t_VEC);
    1174       42828 :       for (l=1,k=1; k<lx; k++)
    1175       39173 :         if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
    1176        3655 :       break;
    1177             :     case t_VEC: case t_COL:
    1178    24873923 :       if (lx != lg(e) || !RgV_is_ZV(e))
    1179           7 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1180    24873916 :       if (lx == 1) return gen_1;
    1181    24853131 :       x = cgetg(lx,t_VEC);
    1182   104218004 :       for (l=1,k=1; k<lx; k++)
    1183    79364873 :         if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
    1184    24853131 :       break;
    1185             :     default:
    1186           0 :       pari_err_TYPE("factorback [not an exponent vector]", e);
    1187           0 :       return NULL;
    1188             :   }
    1189    24856786 :   x[0] = evaltyp(t_VEC) | _evallg(l);
    1190    24856786 :   return gerepileupto(av, gen_product(x, data, _mul));
    1191             : }
    1192             : 
    1193             : GEN
    1194        4284 : idealfactorback(GEN nf, GEN L, GEN e, int red)
    1195             : {
    1196        4284 :   nf = checknf(nf);
    1197        4284 :   if (red) return gen_factorback(L, e, (void*)nf, &idmulred, &idpowred);
    1198        3927 :   else     return gen_factorback(L, e, (void*)nf, &idmul, &idpow);
    1199             : }
    1200             : 
    1201             : GEN
    1202       15961 : nffactorback(GEN nf, GEN L, GEN e)
    1203       15961 : { return gen_factorback(L, e, (void*)checknf(nf), &eltmul, &eltpow); }
    1204             : 
    1205             : GEN
    1206     5316772 : FpV_factorback(GEN L, GEN e, GEN p)
    1207     5316772 : { return gen_factorback(L, e, (void*)p, &Fpmul, &Fppow); }
    1208             : 
    1209             : ulong
    1210           0 : Flv_factorback(GEN L, GEN e, ulong p)
    1211             : {
    1212           0 :   long i, l = lg(e);
    1213           0 :   ulong r = 1UL, ri = 1UL;
    1214           0 :   for (i = 1; i < l; i++)
    1215             :   {
    1216           0 :     long c = e[i];
    1217           0 :     if (!c) continue;
    1218           0 :     if (c < 0)
    1219           0 :       ri = Fl_mul(ri, Fl_powu(L[i],-c,p), p);
    1220             :     else
    1221           0 :       r = Fl_mul(r, Fl_powu(L[i],c,p), p);
    1222             :   }
    1223           0 :   if (ri != 1UL) r = Fl_div(r, ri, p);
    1224           0 :   return r;
    1225             : }
    1226             : GEN
    1227        3135 : FlxqV_factorback(GEN L, GEN e, GEN Tp, ulong p)
    1228             : {
    1229        3135 :   pari_sp av = avma;
    1230        3135 :   GEN Hi = NULL, H = NULL;
    1231        3135 :   long i, l = lg(L), v = get_Flx_var(Tp);
    1232      212705 :   for (i = 1; i < l; i++)
    1233             :   {
    1234      209570 :     GEN x, ei = gel(e,i);
    1235      209570 :     long s = signe(ei);
    1236      209570 :     if (!s) continue;
    1237      183328 :     x = Flxq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
    1238      183328 :     if (s > 0)
    1239       93186 :       H = H? Flxq_mul(H, x, Tp, p): x;
    1240             :     else
    1241       90142 :       Hi = Hi? Flxq_mul(Hi, x, Tp, p): x;
    1242             :   }
    1243        3135 :   if (!Hi)
    1244             :   {
    1245          30 :     if (!H) { set_avma(av); return mkvecsmall2(v,1); }
    1246          30 :     return gerepileuptoleaf(av, H);
    1247             :   }
    1248        3105 :   Hi = Flxq_inv(Hi, Tp, p);
    1249        3105 :   return gerepileuptoleaf(av, H? Flxq_mul(H,Hi,Tp,p): Hi);
    1250             : }
    1251             : GEN
    1252           7 : FqV_factorback(GEN L, GEN e, GEN Tp, GEN p)
    1253             : {
    1254           7 :   pari_sp av = avma;
    1255           7 :   GEN Hi = NULL, H = NULL;
    1256           7 :   long i, l = lg(L);
    1257          21 :   for (i = 1; i < l; i++)
    1258             :   {
    1259          14 :     GEN x, ei = gel(e,i);
    1260          14 :     long s = signe(ei);
    1261          14 :     if (!s) continue;
    1262          14 :     x = Fq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
    1263          14 :     if (s > 0)
    1264          14 :       H = H? Fq_mul(H, x, Tp, p): x;
    1265             :     else
    1266           0 :       Hi = Hi? Fq_mul(Hi, x, Tp, p): x;
    1267             :   }
    1268           7 :   if (!Hi)
    1269             :   {
    1270           7 :     if (!H) { set_avma(av); return gen_1; }
    1271           7 :     return gerepileupto(av, H);
    1272             :   }
    1273           0 :   Hi = Fq_inv(Hi, Tp, p);
    1274           0 :   return gerepileupto(av, H? Fq_mul(H,Hi,Tp,p): Hi);
    1275             : }
    1276             : 
    1277             : GEN
    1278    24683468 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, NULL, &mul, &powi); }
    1279             : GEN
    1280     1426260 : factorback(GEN fa) { return factorback2(fa, NULL); }
    1281             : 
    1282             : GEN
    1283          63 : vecprod(GEN v)
    1284             : {
    1285          63 :   pari_sp av = avma;
    1286          63 :   if (!is_vec_t(typ(v)))
    1287           0 :     pari_err_TYPE("vecprod", v);
    1288          63 :   if (lg(v) == 1) return gen_1;
    1289          56 :   return gerepilecopy(av, gen_product(v, NULL, mul));
    1290             : }
    1291             : 
    1292             : static int
    1293          70 : RgX_is_irred_i(GEN x)
    1294             : {
    1295             :   GEN y, p, pol;
    1296          70 :   long l = lg(x), pa;
    1297             : 
    1298          70 :   if (!signe(x) || l <= 3) return 0;
    1299          70 :   switch(RgX_type(x,&p,&pol,&pa))
    1300             :   {
    1301          21 :     case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
    1302           0 :     case t_COMPLEX: return l == 4;
    1303             :     case t_REAL:
    1304           0 :       if (l == 4) return 1;
    1305           0 :       if (l > 5) return 0;
    1306           0 :       return gsigne(RgX_disc(x)) > 0;
    1307             :   }
    1308          49 :   y = RgX_factor(x, NULL);
    1309          49 :   return (lg(gcoeff(y,1,1))==l);
    1310             : }
    1311             : static int
    1312          70 : RgX_is_irred(GEN x)
    1313          70 : { pari_sp av = avma; return gc_bool(av, RgX_is_irred_i(x)); }
    1314             : long
    1315          70 : polisirreducible(GEN x)
    1316             : {
    1317          70 :   long tx = typ(x);
    1318          70 :   if (tx == t_POL) return RgX_is_irred(x);
    1319           0 :   if (!is_scalar_t(tx)) pari_err_TYPE("polisirreducible",x);
    1320           0 :   return 0;
    1321             : }
    1322             : 
    1323             : /*******************************************************************/
    1324             : /*                                                                 */
    1325             : /*                         GENERIC GCD                             */
    1326             : /*                                                                 */
    1327             : /*******************************************************************/
    1328             : /* x is a COMPLEX or a QUAD */
    1329             : static GEN
    1330        2499 : triv_cont_gcd(GEN x, GEN y)
    1331             : {
    1332        2499 :   pari_sp av = avma;
    1333             :   GEN c;
    1334        2499 :   if (typ(x)==t_COMPLEX)
    1335             :   {
    1336        2177 :     GEN a = gel(x,1), b = gel(x,2);
    1337        2177 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
    1338          21 :     c = ggcd(a,b);
    1339             :   }
    1340             :   else
    1341         322 :     c = ggcd(gel(x,2),gel(x,3));
    1342         343 :   return gerepileupto(av, ggcd(c,y));
    1343             : }
    1344             : 
    1345             : /* y is a PADIC, x a rational number or an INTMOD */
    1346             : static GEN
    1347        2506 : padic_gcd(GEN x, GEN y)
    1348             : {
    1349        2506 :   GEN p = gel(y,2);
    1350        2506 :   long v = gvaluation(x,p), w = valp(y);
    1351        2506 :   if (w < v) v = w;
    1352        2506 :   return powis(p, v);
    1353             : }
    1354             : 
    1355             : /* x,y in Z[i], at least one of which is t_COMPLEX */
    1356             : static GEN
    1357         385 : gauss_gcd(GEN x, GEN y)
    1358             : {
    1359         385 :   pari_sp av = avma;
    1360             :   GEN dx, dy;
    1361         385 :   dx = denom_i(x); x = gmul(x, dx);
    1362         385 :   dy = denom_i(y); y = gmul(y, dy);
    1363        1197 :   while (!gequal0(y))
    1364             :   {
    1365         427 :     GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
    1366         427 :     x = y; y = z;
    1367             :   }
    1368         385 :   x = gauss_normal(x);
    1369         385 :   if (typ(x) == t_COMPLEX)
    1370             :   {
    1371         196 :     if      (gequal0(gel(x,2))) x = gel(x,1);
    1372         196 :     else if (gequal0(gel(x,1))) x = gel(x,2);
    1373             :   }
    1374         385 :   return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
    1375             : }
    1376             : 
    1377             : static int
    1378        2520 : c_is_rational(GEN x)
    1379        2520 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
    1380             : static GEN
    1381        1281 : c_zero_gcd(GEN c)
    1382             : {
    1383        1281 :   GEN x = gel(c,1), y = gel(c,2);
    1384        1281 :   long tx = typ(x), ty = typ(y);
    1385        1281 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
    1386          42 :   if (tx == t_PADIC || tx == t_INTMOD
    1387          35 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
    1388          35 :   return gauss_gcd(c, gen_0);
    1389             : }
    1390             : 
    1391             : /* gcd(x, 0) */
    1392             : static GEN
    1393     8225307 : zero_gcd(GEN x)
    1394             : {
    1395             :   pari_sp av;
    1396     8225307 :   switch(typ(x))
    1397             :   {
    1398       27957 :     case t_INT: return absi(x);
    1399        5001 :     case t_FRAC: return absfrac(x);
    1400        1281 :     case t_COMPLEX: return c_zero_gcd(x);
    1401        3913 :     case t_REAL: return gen_1;
    1402        1043 :     case t_PADIC: return powis(gel(x,2), valp(x));
    1403         245 :     case t_SER: return pol_xnall(valp(x), varn(x));
    1404             :     case t_POLMOD: {
    1405        3314 :       GEN d = gel(x,2);
    1406        3314 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
    1407         406 :       return isinexact(d)? zero_gcd(d): gcopy(d);
    1408             :     }
    1409             :     case t_POL:
    1410     7920035 :       if (!isinexact(x)) break;
    1411          14 :       av = avma;
    1412          14 :       return gerepileupto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
    1413             : 
    1414             :     case t_RFRAC:
    1415      217067 :       if (!isinexact(x)) break;
    1416           0 :       av = avma;
    1417           0 :       return gerepileupto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
    1418             :   }
    1419     8182539 :   return gcopy(x);
    1420             : }
    1421             : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
    1422             : static GEN
    1423     9337533 : zero_gcd2(GEN y, GEN z)
    1424             : {
    1425             :   pari_sp av;
    1426     9337533 :   switch(typ(z))
    1427             :   {
    1428     8211128 :     case t_INT: return zero_gcd(y);
    1429             :     case t_INTMOD:
    1430     1125355 :       av = avma;
    1431     1125355 :       return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
    1432             :     case t_FFELT:
    1433        1050 :       av = avma;
    1434        1050 :       return gerepileupto(av, gmul(y, FF_1(z)));
    1435             :     default:
    1436           0 :       pari_err_TYPE("zero_gcd", z);
    1437             :   }
    1438           0 :   return NULL;
    1439             : }
    1440             : static GEN
    1441     2785979 : cont_gcd_pol_i(GEN x, GEN y) { return scalarpol(ggcd(content(x),y), varn(x));}
    1442             : /* tx = t_POL, y considered as constant */
    1443             : static GEN
    1444     2785972 : cont_gcd_pol(GEN x, GEN y)
    1445     2785972 : { pari_sp av = avma; return gerepileupto(av, cont_gcd_pol_i(x,y)); }
    1446             : /* tx = t_RFRAC, y considered as constant */
    1447             : static GEN
    1448      864918 : cont_gcd_rfrac(GEN x, GEN y)
    1449             : {
    1450      864918 :   pari_sp av = avma;
    1451      864918 :   GEN cx; x = primitive_part(x, &cx);
    1452             :   /* e.g. Mod(1,2) / (2*y+1) => primitive_part = Mod(1,2)*y^0 */
    1453      864918 :   if (typ(x) != t_RFRAC) x = cont_gcd_pol_i(x, y);
    1454      864911 :   else x = gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2));
    1455      864918 :   return gerepileupto(av, x);
    1456             : }
    1457             : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
    1458             : static GEN
    1459        2572 : cont_gcd_gen(GEN x, GEN y)
    1460             : {
    1461        2572 :   pari_sp av = avma;
    1462        2572 :   return gerepileupto(av, ggcd(content(x),y));
    1463             : }
    1464             : /* !is_const(tx), y considered as constant */
    1465             : static GEN
    1466     2799464 : cont_gcd(GEN x, long tx, GEN y)
    1467             : {
    1468     2799464 :   switch(tx)
    1469             :   {
    1470       10934 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1471     2785958 :     case t_POL: return cont_gcd_pol(x,y);
    1472        2572 :     default: return cont_gcd_gen(x,y);
    1473             :   }
    1474             : }
    1475             : static GEN
    1476    13338378 : gcdiq(GEN x, GEN y)
    1477             : {
    1478             :   GEN z;
    1479    13338378 :   if (!signe(x)) return Q_abs(y);
    1480     2809943 :   z = cgetg(3,t_FRAC);
    1481     2809943 :   gel(z,1) = gcdii(x,gel(y,1));
    1482     2809943 :   gel(z,2) = icopy(gel(y,2));
    1483     2809943 :   return z;
    1484             : }
    1485             : static GEN
    1486    21147050 : gcdqq(GEN x, GEN y)
    1487             : {
    1488    21147050 :   GEN z = cgetg(3,t_FRAC);
    1489    21147050 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1490    21147050 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1491    21147050 :   return z;
    1492             : }
    1493             : /* assume x,y t_INT or t_FRAC */
    1494             : GEN
    1495   153026380 : Q_gcd(GEN x, GEN y)
    1496             : {
    1497   153026380 :   long tx = typ(x), ty = typ(y);
    1498   153026380 :   if (tx == t_INT)
    1499   119673618 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1500             :   else
    1501    33352762 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1502             : }
    1503             : 
    1504             : GEN
    1505    27208442 : ggcd(GEN x, GEN y)
    1506             : {
    1507    27208442 :   long vx, vy, tx = typ(x), ty = typ(y);
    1508             :   pari_sp av, tetpil;
    1509             :   GEN p1,z;
    1510             : 
    1511    54416884 :   if (is_noncalc_t(tx) || is_matvec_t(tx) ||
    1512    54416884 :       is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
    1513    27208442 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1514             :   /* tx <= ty */
    1515    27208442 :   z = gisexactzero(x); if (z) return zero_gcd2(y,z);
    1516    24147467 :   z = gisexactzero(y); if (z) return zero_gcd2(x,z);
    1517    17870909 :   if (is_const_t(tx))
    1518             :   {
    1519     9880567 :     if (ty == tx) switch(tx)
    1520             :     {
    1521             :       case t_INT:
    1522     2934831 :         return gcdii(x,y);
    1523             : 
    1524     4083310 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1525     4083310 :         if (equalii(gel(x,1),gel(y,1)))
    1526     4083303 :           gel(z,1) = icopy(gel(x,1));
    1527             :         else
    1528           7 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1529     4083310 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1530             :         else
    1531             :         {
    1532     4083310 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1533     4083310 :           if (!equali1(p1))
    1534             :           {
    1535           7 :             p1 = gcdii(p1,gel(y,2));
    1536           7 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1537           7 :             else p1 = gerepileuptoint(av, p1);
    1538             :           }
    1539     4083310 :           gel(z,2) = p1;
    1540             :         }
    1541     4083310 :         return z;
    1542             : 
    1543             :       case t_FRAC:
    1544       67528 :         return gcdqq(x,y);
    1545             : 
    1546             :       case t_FFELT:
    1547        5397 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1548        5397 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1549             : 
    1550             :       case t_COMPLEX:
    1551          21 :         if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
    1552           7 :         return triv_cont_gcd(y,x);
    1553             : 
    1554             :       case t_PADIC:
    1555          14 :         if (!equalii(gel(x,2),gel(y,2))) return gen_1;
    1556           7 :         return powis(gel(y,2), minss(valp(x), valp(y)));
    1557             : 
    1558             :       case t_QUAD:
    1559         133 :         av=avma; p1=gdiv(x,y);
    1560         133 :         if (gequal0(gel(p1,3)))
    1561             :         {
    1562          14 :           p1=gel(p1,2);
    1563          14 :           if (typ(p1)==t_INT) { set_avma(av); return gcopy(y); }
    1564           7 :           tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
    1565             :         }
    1566         119 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {set_avma(av); return gcopy(y);}
    1567         112 :         p1 = ginv(p1); set_avma(av);
    1568         112 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1569         105 :         return triv_cont_gcd(y,x);
    1570             : 
    1571           0 :       default: return gen_1; /* t_REAL */
    1572             :     }
    1573     2789333 :     if (is_const_t(ty)) switch(tx)
    1574             :     {
    1575             :       case t_INT:
    1576       16959 :         switch(ty)
    1577             :         {
    1578          56 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1579          56 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1580          56 :             p1 = gcdii(gel(y,1),gel(y,2));
    1581          56 :             if (!equali1(p1)) {
    1582          14 :               p1 = gcdii(x,p1);
    1583          14 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1584             :               else
    1585          14 :                 p1 = gerepileuptoint(av, p1);
    1586             :             }
    1587          56 :             gel(z,2) = p1; return z;
    1588             : 
    1589        8246 :           case t_REAL: return gen_1;
    1590             : 
    1591             :           case t_FRAC:
    1592        3526 :             return gcdiq(x,y);
    1593             : 
    1594             :           case t_COMPLEX:
    1595        2401 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1596        2072 :             return triv_cont_gcd(y,x);
    1597             : 
    1598             :           case t_FFELT:
    1599         112 :             if (!FF_equal0(y)) return FF_1(y);
    1600           0 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1601             : 
    1602             :           case t_PADIC:
    1603        2492 :             return padic_gcd(x,y);
    1604             : 
    1605             :           case t_QUAD:
    1606         126 :             return triv_cont_gcd(y,x);
    1607             :           default:
    1608           0 :             pari_err_TYPE2("gcd",x,y);
    1609             :         }
    1610             : 
    1611             :       case t_REAL:
    1612          14 :         switch(ty)
    1613             :         {
    1614             :           case t_INTMOD:
    1615             :           case t_FFELT:
    1616          14 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1617           0 :           default: return gen_1;
    1618             :         }
    1619             : 
    1620             :       case t_INTMOD:
    1621          49 :         switch(ty)
    1622             :         {
    1623             :           case t_FRAC:
    1624          14 :             av = avma; p1=gcdii(gel(x,1),gel(y,2)); set_avma(av);
    1625          14 :             if (!equali1(p1)) pari_err_OP("gcd",x,y);
    1626           7 :             return ggcd(gel(y,1), x);
    1627             : 
    1628             :           case t_FFELT:
    1629             :           {
    1630          14 :             GEN p = gel(y,4);
    1631          14 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1632           7 :             if (!FF_equal0(y)) return FF_1(y);
    1633           0 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1634             :           }
    1635             : 
    1636             :           case t_COMPLEX: case t_QUAD:
    1637          14 :             return triv_cont_gcd(y,x);
    1638             : 
    1639             :           case t_PADIC:
    1640           7 :             return padic_gcd(x,y);
    1641             : 
    1642           0 :           default: pari_err_TYPE2("gcd",x,y);
    1643             :         }
    1644             : 
    1645             :       case t_FRAC:
    1646         210 :         switch(ty)
    1647             :         {
    1648             :           case t_COMPLEX:
    1649          84 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1650             :           case t_QUAD:
    1651         154 :             return triv_cont_gcd(y,x);
    1652             :           case t_FFELT:
    1653             :           {
    1654          42 :             GEN p = gel(y,4);
    1655          42 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1656          21 :             if (!FF_equal0(y)) return FF_1(y);
    1657           0 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1658             :           }
    1659             : 
    1660             :           case t_PADIC:
    1661           7 :             return padic_gcd(x,y);
    1662             : 
    1663           0 :           default: pari_err_TYPE2("gcd",x,y);
    1664             :         }
    1665             :       case t_FFELT:
    1666          70 :         switch(ty)
    1667             :         {
    1668             :           case t_PADIC:
    1669             :           {
    1670          42 :             GEN p = gel(y,2);
    1671          42 :             long v = valp(y);
    1672          42 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1673          14 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1674             :           }
    1675          28 :           default: pari_err_TYPE2("gcd",x,y);
    1676             :         }
    1677             : 
    1678             :       case t_COMPLEX:
    1679          14 :         switch(ty)
    1680             :         {
    1681             :           case t_PADIC:
    1682          14 :           case t_QUAD: return triv_cont_gcd(x,y);
    1683           0 :           default: pari_err_TYPE2("gcd",x,y);
    1684             :         }
    1685             : 
    1686             :       case t_PADIC:
    1687           7 :         switch(ty)
    1688             :         {
    1689           7 :           case t_QUAD: return triv_cont_gcd(y,x);
    1690           0 :           default: pari_err_TYPE2("gcd",x,y);
    1691             :         }
    1692             : 
    1693           0 :       default: return gen_1; /* tx = t_REAL */
    1694             :     }
    1695     2772010 :     return cont_gcd(y,ty, x);
    1696             :   }
    1697             : 
    1698     7990342 :   if (tx == t_POLMOD)
    1699             :   {
    1700       15308 :     if (ty == t_POLMOD)
    1701             :     {
    1702       15252 :       GEN T = gel(x,1);
    1703       15252 :       z = cgetg(3,t_POLMOD);
    1704       15252 :       T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
    1705       15252 :       gel(z,1) = T;
    1706       15252 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1707             :       else
    1708             :       {
    1709             :         GEN X, Y, d;
    1710       15252 :         av = avma; X = gel(x,2); Y = gel(y,2);
    1711       15252 :         d = ggcd(content(X), content(Y));
    1712       15252 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1713       15252 :         p1 = ggcd(T, X);
    1714       15252 :         gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
    1715             :       }
    1716       15252 :       return z;
    1717             :     }
    1718          56 :     vx = varn(gel(x,1));
    1719          56 :     switch(ty)
    1720             :     {
    1721             :       case t_POL:
    1722          28 :         vy = varn(y);
    1723          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1724          14 :         z = cgetg(3,t_POLMOD);
    1725          14 :         gel(z,1) = RgX_copy(gel(x,1));
    1726          14 :         av = avma; p1 = ggcd(gel(x,1),gel(x,2));
    1727          14 :         gel(z,2) = gerepileupto(av, ggcd(p1,y));
    1728          14 :         return z;
    1729             : 
    1730             :       case t_RFRAC:
    1731          28 :         vy = varn(gel(y,2));
    1732          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1733          28 :         av = avma;
    1734          28 :         p1 = ggcd(gel(x,1),gel(y,2));
    1735          28 :         if (degpol(p1)) pari_err_OP("gcd",x,y);
    1736          21 :         set_avma(av); return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1737             :     }
    1738             :   }
    1739             : 
    1740     7975034 :   vx = gvar(x);
    1741     7975034 :   vy = gvar(y);
    1742     7975034 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1743     7963029 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1744             : 
    1745             :   /* vx = vy: same main variable */
    1746     7947580 :   switch(tx)
    1747             :   {
    1748             :     case t_POL:
    1749     7786660 :       switch(ty)
    1750             :       {
    1751     6932655 :         case t_POL: return RgX_gcd(x,y);
    1752             :         case t_SER:
    1753          21 :           z = ggcd(content(x), content(y));
    1754          21 :           return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
    1755      853984 :         case t_RFRAC: return cont_gcd_rfrac(y, x);
    1756             :       }
    1757           0 :       break;
    1758             : 
    1759             :     case t_SER:
    1760          14 :       z = ggcd(content(x), content(y));
    1761          14 :       switch(ty)
    1762             :       {
    1763           7 :         case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
    1764           7 :         case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
    1765             :       }
    1766           0 :       break;
    1767             : 
    1768             :     case t_RFRAC:
    1769             :     {
    1770      160906 :       GEN xd = gel(x,2), yd = gel(y,2);
    1771      160906 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1772      160906 :       z = cgetg(3,t_RFRAC); av = avma;
    1773      160906 :       gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1774      160906 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1775             :     }
    1776             :   }
    1777           0 :   pari_err_TYPE2("gcd",x,y);
    1778             :   return NULL; /* LCOV_EXCL_LINE */
    1779             : }
    1780             : GEN
    1781        4071 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
    1782             : 
    1783             : static GEN
    1784        1715 : fix_lcm(GEN x)
    1785             : {
    1786             :   GEN t;
    1787        1715 :   switch(typ(x))
    1788             :   {
    1789        1617 :     case t_INT: if (signe(x)<0) x = negi(x);
    1790        1617 :       break;
    1791             :     case t_POL:
    1792          98 :       if (lg(x) <= 2) break;
    1793          98 :       t = leading_coeff(x);
    1794          98 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1795             :   }
    1796        1715 :   return x;
    1797             : }
    1798             : GEN
    1799        2877 : glcm0(GEN x, GEN y)
    1800             : {
    1801        2877 :   if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
    1802        2828 :   return glcm(x,y);
    1803             : }
    1804             : GEN
    1805        1617 : ZV_lcm(GEN x) { return fix_lcm(gassoc_proto(lcmii,x,NULL)); }
    1806             : GEN
    1807        3171 : glcm(GEN x, GEN y)
    1808             : {
    1809             :   pari_sp av;
    1810             :   GEN z;
    1811        3171 :   if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
    1812          63 :   av = avma; z = ggcd(x,y);
    1813          63 :   if (!gequal1(z))
    1814             :   {
    1815          63 :     if (gequal0(z)) { set_avma(av); return gmul(x,y); }
    1816          49 :     y = gdiv(y,z);
    1817             :   }
    1818          49 :   return gerepileupto(av, fix_lcm(gmul(x,y)));
    1819             : }
    1820             : 
    1821             : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
    1822             : static int
    1823           0 : pol_approx0(GEN r, GEN x, int exact)
    1824             : {
    1825             :   long i, lx,lr;
    1826           0 :   if (exact) return gequal0(r);
    1827           0 :   lx = lg(x);
    1828           0 :   lr = lg(r); if (lr < lx) lx = lr;
    1829           0 :   for (i=2; i<lx; i++)
    1830           0 :     if (!approx_0(gel(r,i), gel(x,i))) return 0;
    1831           0 :   return 1;
    1832             : }
    1833             : 
    1834             : GEN
    1835           0 : RgX_gcd_simple(GEN x, GEN y)
    1836             : {
    1837           0 :   pari_sp av1, av = avma;
    1838           0 :   GEN r, yorig = y;
    1839           0 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1840             : 
    1841             :   for(;;)
    1842             :   {
    1843           0 :     av1 = avma; r = RgX_rem(x,y);
    1844           0 :     if (pol_approx0(r, x, exact))
    1845             :     {
    1846           0 :       set_avma(av1);
    1847           0 :       if (y == yorig) return RgX_copy(y);
    1848           0 :       y = normalizepol_approx(y, lg(y));
    1849           0 :       if (lg(y) == 3) { set_avma(av); return pol_1(varn(x)); }
    1850           0 :       return gerepileupto(av,y);
    1851             :     }
    1852           0 :     x = y; y = r;
    1853           0 :     if (gc_needed(av,1)) {
    1854           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1855           0 :       gerepileall(av,2, &x,&y);
    1856             :     }
    1857             :   }
    1858             : }
    1859             : GEN
    1860           0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
    1861             : {
    1862           0 :   pari_sp av = avma;
    1863             :   GEN q, r, d, d1, u, v, v1;
    1864           0 :   int exact = !(isinexactreal(a) || isinexactreal(b));
    1865             : 
    1866           0 :   d = a; d1 = b; v = gen_0; v1 = gen_1;
    1867             :   for(;;)
    1868             :   {
    1869           0 :     if (pol_approx0(d1, a, exact)) break;
    1870           0 :     q = poldivrem(d,d1, &r);
    1871           0 :     v = gsub(v, gmul(q,v1));
    1872           0 :     u=v; v=v1; v1=u;
    1873           0 :     u=r; d=d1; d1=u;
    1874             :   }
    1875           0 :   u = gsub(d, gmul(b,v));
    1876           0 :   u = RgX_div(u,a);
    1877             : 
    1878           0 :   gerepileall(av, 3, &u,&v,&d);
    1879           0 :   *pu = u;
    1880           0 :   *pv = v; return d;
    1881             : }
    1882             : /*******************************************************************/
    1883             : /*                                                                 */
    1884             : /*                    CONTENT / PRIMITIVE PART                     */
    1885             : /*                                                                 */
    1886             : /*******************************************************************/
    1887             : 
    1888             : GEN
    1889    73252153 : content(GEN x)
    1890             : {
    1891    73252153 :   long lx, i, t, tx = typ(x);
    1892    73252153 :   pari_sp av = avma;
    1893             :   GEN c;
    1894             : 
    1895    73252153 :   if (is_scalar_t(tx)) return zero_gcd(x);
    1896    73244400 :   switch(tx)
    1897             :   {
    1898             :     case t_RFRAC:
    1899             :     {
    1900      864953 :       GEN n = gel(x,1), d = gel(x,2);
    1901             :       /* -- varncmp(vn, vd) < 0 can't happen
    1902             :        * -- if n is POLMOD, its main variable (in the sense of gvar2)
    1903             :        *    has lower priority than denominator */
    1904      864953 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    1905      824608 :         n = isinexact(n)? zero_gcd(n): gcopy(n);
    1906             :       else
    1907       40345 :         n = content(n);
    1908      864953 :       return gerepileupto(av, gdiv(n, content(d)));
    1909             :     }
    1910             : 
    1911             :     case t_VEC: case t_COL:
    1912     1981870 :       lx = lg(x); if (lx==1) return gen_0;
    1913     1981863 :       break;
    1914             : 
    1915             :     case t_MAT:
    1916             :     {
    1917             :       long hx, j;
    1918         245 :       lx = lg(x);
    1919         245 :       if (lx == 1) return gen_0;
    1920         238 :       hx = lgcols(x);
    1921         238 :       if (hx == 1) return gen_0;
    1922         231 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    1923         231 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    1924         231 :       c = content(gel(x,1));
    1925         462 :       for (j=2; j<lx; j++)
    1926         231 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    1927         231 :       if (typ(c) == t_INTMOD || isinexact(c)) { set_avma(av); return gen_1; }
    1928         231 :       return gerepileupto(av,c);
    1929             :     }
    1930             : 
    1931             :     case t_POL: case t_SER:
    1932    70397297 :       lx = lg(x); if (lx == 2) return gen_0;
    1933    70397283 :       break;
    1934          21 :     case t_VECSMALL: return utoi(zv_content(x));
    1935             :     case t_QFR: case t_QFI:
    1936          14 :       lx = 4; break;
    1937             : 
    1938           0 :     default: pari_err_TYPE("content",x);
    1939             :       return NULL; /* LCOV_EXCL_LINE */
    1940             :   }
    1941   210877930 :   for (i=lontyp[tx]; i<lx; i++)
    1942   151674002 :     if (typ(gel(x,i)) != t_INT) break;
    1943    72379160 :   lx--; c = gel(x,lx);
    1944    72379160 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    1945    72379160 :   if (i > lx)
    1946             :   { /* integer coeffs */
    1947   120154987 :     while (lx-- > lontyp[tx])
    1948             :     {
    1949    59558979 :       c = gcdii(c, gel(x,lx));
    1950    59558979 :       if (equali1(c)) { set_avma(av); return gen_1; }
    1951             :     }
    1952             :   }
    1953             :   else
    1954             :   {
    1955    13175232 :     if (isinexact(c)) c = zero_gcd(c);
    1956    49690252 :     while (lx-- > lontyp[tx])
    1957             :     {
    1958    23339788 :       GEN d = gel(x,lx);
    1959    23339788 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    1960    23339788 :       c = ggcd(c, d);
    1961             :     }
    1962    13175232 :     if (isinexact(c)) { set_avma(av); return gen_1; }
    1963             :   }
    1964    14567312 :   switch(typ(c))
    1965             :   {
    1966             :     case t_INT:
    1967     1400623 :       if (signe(c) < 0) c = negi(c);
    1968     1400623 :       break;
    1969             :     case t_VEC: case t_COL: case t_MAT:
    1970           0 :       pari_err_TYPE("content",x);
    1971             :   }
    1972             : 
    1973    14567312 :   return av==avma? gcopy(c): gerepileupto(av,c);
    1974             : }
    1975             : 
    1976             : GEN
    1977     1157471 : primitive_part(GEN x, GEN *ptc)
    1978             : {
    1979     1157471 :   pari_sp av = avma;
    1980     1157471 :   GEN c = content(x);
    1981     1157471 :   if (gequal1(c)) { set_avma(av); c = NULL; }
    1982       80206 :   else if (!gequal0(c)) x = gdiv(x,c);
    1983     1157471 :   if (ptc) *ptc = c;
    1984     1157471 :   return x;
    1985             : }
    1986             : GEN
    1987        2191 : primpart(GEN x) { return primitive_part(x, NULL); }
    1988             : 
    1989             : static GEN
    1990    45494491 : Q_content_v(GEN x, long i, long l)
    1991             : {
    1992    45494491 :   pari_sp av = avma;
    1993    45494491 :   GEN d = Q_content_safe(gel(x,i));
    1994    45494518 :   if (!d) return NULL;
    1995   198516112 :   for (i++; i < l; i++)
    1996             :   {
    1997   153022401 :     GEN c = Q_content_safe(gel(x,i));
    1998   153022491 :     if (!c) return NULL;
    1999   153022491 :     d = Q_gcd(d, c);
    2000             :   }
    2001    45493711 :   return gerepileupto(av, d);
    2002             : }
    2003             : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    2004             :  * of Q(x2,...,xn)[x1] */
    2005             : GEN
    2006   232720794 : Q_content_safe(GEN x)
    2007             : {
    2008             :   long l;
    2009   232720794 :   switch(typ(x))
    2010             :   {
    2011   159400764 :     case t_INT:  return absi(x);
    2012    27395812 :     case t_FRAC: return absfrac(x);
    2013             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2014    15401441 :       l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
    2015             :     case t_POL:
    2016    30493338 :       l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
    2017       28970 :     case t_POLMOD: return Q_content_safe(gel(x,2));
    2018             :     case t_RFRAC:
    2019             :     {
    2020             :       GEN a, b;
    2021          49 :       a = Q_content(gel(x,1)); if (!a) return NULL;
    2022          49 :       b = Q_content(gel(x,2)); if (!b) return NULL;
    2023          49 :       return gdiv(a, b);
    2024             :     }
    2025             :   }
    2026         420 :   return NULL;
    2027             : }
    2028             : GEN
    2029       84965 : Q_content(GEN x)
    2030             : {
    2031       84965 :   GEN c = Q_content_safe(x);
    2032       84965 :   if (!c) pari_err_TYPE("Q_content",x);
    2033       84965 :   return c;
    2034             : }
    2035             : 
    2036             : GEN
    2037       13146 : ZX_content(GEN x)
    2038             : {
    2039       13146 :   long i, l = lg(x);
    2040             :   GEN d;
    2041             :   pari_sp av;
    2042             : 
    2043       13146 :   if (l == 2) return gen_0;
    2044       13146 :   d = gel(x,2);
    2045       13146 :   if (l == 3) return absi(d);
    2046        9170 :   av = avma;
    2047        9170 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    2048        9170 :   if (signe(d) < 0) d = negi(d);
    2049        9170 :   return gerepileuptoint(av, d);
    2050             : }
    2051             : 
    2052             : static GEN
    2053      116578 : Z_content_v(GEN x, long i, long l)
    2054             : {
    2055      116578 :   pari_sp av = avma;
    2056      116578 :   GEN d = Z_content(gel(x,i));
    2057      116578 :   if (!d) return NULL;
    2058      487389 :   for (i++; i<l; i++)
    2059             :   {
    2060      375361 :     GEN c = Z_content(gel(x,i));
    2061      375361 :     if (!c) return NULL;
    2062      375186 :     d = gcdii(d, c); if (equali1(d)) return NULL;
    2063      375067 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    2064             :   }
    2065      112028 :   return gerepileuptoint(av, d);
    2066             : }
    2067             : /* return NULL for 1 */
    2068             : GEN
    2069      493752 : Z_content(GEN x)
    2070             : {
    2071             :   long l;
    2072      493752 :   switch(typ(x))
    2073             :   {
    2074             :     case t_INT:
    2075      357847 :       if (is_pm1(x)) return NULL;
    2076      356811 :       return absi(x);
    2077             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2078       11886 :       l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
    2079             :     case t_POL:
    2080      124019 :       l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
    2081           0 :     case t_POLMOD: return Z_content(gel(x,2));
    2082             :   }
    2083           0 :   pari_err_TYPE("Z_content", x);
    2084             :   return NULL; /* LCOV_EXCL_LINE */
    2085             : }
    2086             : 
    2087             : static GEN
    2088    13408387 : Q_denom_v(GEN x, long i, long l)
    2089             : {
    2090    13408387 :   pari_sp av = avma;
    2091    13408387 :   GEN d = Q_denom_safe(gel(x,i));
    2092    13408394 :   if (!d) return NULL;
    2093    53876182 :   for (i++; i<l; i++)
    2094             :   {
    2095    40467821 :     GEN D = Q_denom_safe(gel(x,i));
    2096    40467788 :     if (!D) return NULL;
    2097    40467788 :     if (D != gen_1) d = lcmii(d, D);
    2098    40467788 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    2099             :   }
    2100    13408361 :   return gerepileuptoint(av, d);
    2101             : }
    2102             : /* NOT MEMORY CLEAN (because of t_FRAC).
    2103             :  * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    2104             :  * of Q(x2,...,xn)[x1] */
    2105             : GEN
    2106    66924853 : Q_denom_safe(GEN x)
    2107             : {
    2108             :   long l;
    2109    66924853 :   switch(typ(x))
    2110             :   {
    2111    41998880 :     case t_INT: return gen_1;
    2112          28 :     case t_PADIC: l = valp(x); return l < 0? powiu(gel(x,2), -l): gen_1;
    2113    11351766 :     case t_FRAC: return gel(x,2);
    2114         448 :     case t_QUAD: return Q_denom_v(x, 2, 4);
    2115             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    2116    10074664 :       l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
    2117             :     case t_POL: case t_SER:
    2118     3395152 :       l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
    2119      103789 :     case t_POLMOD: return Q_denom(gel(x,2));
    2120             :     case t_RFRAC:
    2121             :     {
    2122             :       GEN a, b;
    2123          56 :       a = Q_content(gel(x,1)); if (!a) return NULL;
    2124          56 :       b = Q_content(gel(x,2)); if (!b) return NULL;
    2125          56 :       return Q_denom(gdiv(a, b));
    2126             :     }
    2127             :   }
    2128          70 :   return NULL;
    2129             : }
    2130             : GEN
    2131     1673176 : Q_denom(GEN x)
    2132             : {
    2133     1673176 :   GEN d = Q_denom_safe(x);
    2134     1673176 :   if (!d) pari_err_TYPE("Q_denom",x);
    2135     1673176 :   return d;
    2136             : }
    2137             : 
    2138             : GEN
    2139    11375341 : Q_remove_denom(GEN x, GEN *ptd)
    2140             : {
    2141    11375341 :   GEN d = Q_denom_safe(x);
    2142    11375209 :   if (d) { if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d); }
    2143    11375216 :   if (ptd) *ptd = d;
    2144    11375216 :   return x;
    2145             : }
    2146             : 
    2147             : /* return y = x * d, assuming x rational, and d,y integral */
    2148             : GEN
    2149    66281667 : Q_muli_to_int(GEN x, GEN d)
    2150             : {
    2151             :   long i, l;
    2152             :   GEN y, xn, xd;
    2153             :   pari_sp av;
    2154             : 
    2155    66281667 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    2156    66281667 :   switch (typ(x))
    2157             :   {
    2158             :     case t_INT:
    2159    22772795 :       return mulii(x,d);
    2160             : 
    2161             :     case t_FRAC:
    2162    32819392 :       xn = gel(x,1);
    2163    32819392 :       xd = gel(x,2); av = avma;
    2164    32819392 :       y = mulii(xn, diviiexact(d, xd));
    2165    32819392 :       return gerepileuptoint(av, y);
    2166             :     case t_COMPLEX:
    2167           7 :       y = cgetg(3,t_COMPLEX);
    2168           7 :       gel(y,1) = Q_muli_to_int(gel(x,1),d);
    2169           7 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2170           7 :       return y;
    2171             :     case t_PADIC:
    2172          14 :       y = gcopy(x); if (!isint1(d)) setvalp(y, 0);
    2173          14 :       return y;
    2174             :     case t_QUAD:
    2175         175 :       y = cgetg(4,t_QUAD);
    2176         175 :       gel(y,1) = ZX_copy(gel(x,1));
    2177         175 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2178         175 :       gel(y,3) = Q_muli_to_int(gel(x,3),d); return y;
    2179             : 
    2180             :     case t_VEC: case t_COL: case t_MAT:
    2181     5856287 :       y = cgetg_copy(x, &l);
    2182     5856287 :       for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2183     5856287 :       return y;
    2184             : 
    2185             :     case t_POL: case t_SER:
    2186     4771130 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2187     4771130 :       for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2188     4771130 :       return y;
    2189             : 
    2190             :     case t_POLMOD:
    2191       61846 :       retmkpolmod(Q_muli_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2192             :     case t_RFRAC:
    2193          21 :       return gmul(x, d);
    2194             :   }
    2195           0 :   pari_err_TYPE("Q_muli_to_int",x);
    2196             :   return NULL; /* LCOV_EXCL_LINE */
    2197             : }
    2198             : 
    2199             : static void
    2200     2209719 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
    2201             : {
    2202             :   long e;
    2203     2209719 :   switch(typ(c))
    2204             :   {
    2205             :     case t_REAL:
    2206     1559295 :       *exact = 0;
    2207     1559295 :       e = signe(c)? expo(c) - bit_prec(c): expo(c);
    2208     1559295 :       break;
    2209             :     case t_INT:
    2210      647533 :       if (!signe(c)) return;
    2211      275454 :       e = expi(c) + 32;
    2212      275454 :       break;
    2213             :     case t_FRAC:
    2214        2863 :       e = expi(gel(c,1)) - expi(gel(c,2)) + 32;
    2215        2863 :       if (exact) *D = lcmii(*D, gel(c,2));
    2216        2863 :       break;
    2217             :     default:
    2218          28 :       pari_err_TYPE("rescale_to_int",c);
    2219             :       return; /* LCOV_EXCL_LINE */
    2220             :   }
    2221     1837612 :   if (e < *emin) *emin = e;
    2222             : }
    2223             : GEN
    2224      227259 : RgM_rescale_to_int(GEN x)
    2225             : {
    2226      227259 :   long lx = lg(x), i,j, hx, emin;
    2227             :   GEN D;
    2228             :   int exact;
    2229             : 
    2230      227259 :   if (lx == 1) return cgetg(1,t_MAT);
    2231      227259 :   hx = lgcols(x);
    2232      227259 :   exact = 1;
    2233      227259 :   emin = HIGHEXPOBIT;
    2234      227259 :   D = gen_1;
    2235      834055 :   for (j = 1; j < lx; j++)
    2236      606824 :     for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
    2237      227231 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2238      227175 :   return grndtoi(gmul2n(x, -emin), &i);
    2239             : }
    2240             : GEN
    2241       36960 : RgX_rescale_to_int(GEN x)
    2242             : {
    2243       36960 :   long lx = lg(x), i, emin;
    2244             :   GEN D;
    2245             :   int exact;
    2246       36960 :   if (lx == 2) return gcopy(x); /* rare */
    2247       36960 :   exact = 1;
    2248       36960 :   emin = HIGHEXPOBIT;
    2249       36960 :   D = gen_1;
    2250       36960 :   for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
    2251       36960 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2252       36225 :   return grndtoi(gmul2n(x, -emin), &i);
    2253             : }
    2254             : 
    2255             : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    2256             : static GEN
    2257    10128740 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    2258             : {
    2259             :   long i, l;
    2260             :   GEN y, xn, xd;
    2261             :   pari_sp av;
    2262             : 
    2263    10128740 :   switch(typ(x))
    2264             :   {
    2265             :     case t_INT:
    2266     2693473 :       av = avma; y = diviiexact(x,d);
    2267     2693473 :       return gerepileuptoint(av, mulii(y,n));
    2268             : 
    2269             :     case t_FRAC:
    2270     5072392 :       xn = gel(x,1);
    2271     5072392 :       xd = gel(x,2); av = avma;
    2272     5072392 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    2273     5072392 :       return gerepileuptoint(av, y);
    2274             : 
    2275             :     case t_VEC: case t_COL: case t_MAT:
    2276      446342 :       y = cgetg_copy(x, &l);
    2277      446342 :       for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2278      446342 :       return y;
    2279             : 
    2280             :     case t_POL:
    2281     1916533 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2282     1916533 :       for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2283     1916533 :       return y;
    2284             : 
    2285             :     case t_RFRAC:
    2286           0 :       av = avma;
    2287           0 :       return gerepileupto(av, gmul(x,mkfrac(n,d)));
    2288             : 
    2289             :     case t_POLMOD:
    2290           0 :       retmkpolmod(Q_divmuli_to_int(gel(x,2), d,n), RgX_copy(gel(x,1)));
    2291             :   }
    2292           0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    2293             :   return NULL; /* LCOV_EXCL_LINE */
    2294             : }
    2295             : 
    2296             : /* return x / d. x: rational; d,result: integral. */
    2297             : static GEN
    2298    22990832 : Q_divi_to_int(GEN x, GEN d)
    2299             : {
    2300             :   long i, l;
    2301             :   GEN y;
    2302             : 
    2303    22990832 :   switch(typ(x))
    2304             :   {
    2305             :     case t_INT:
    2306    18414829 :       return diviiexact(x,d);
    2307             : 
    2308             :     case t_VEC: case t_COL: case t_MAT:
    2309     1325268 :       y = cgetg_copy(x, &l);
    2310     1325303 :       for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2311     1325252 :       return y;
    2312             : 
    2313             :     case t_POL:
    2314     3244268 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2315     3244268 :       for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2316     3244268 :       return y;
    2317             : 
    2318             :     case t_RFRAC:
    2319           7 :       return gdiv(x,d);
    2320             : 
    2321             :     case t_POLMOD:
    2322        6460 :       retmkpolmod(Q_divi_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
    2323             :   }
    2324           0 :   pari_err_TYPE("Q_divi_to_int",x);
    2325             :   return NULL; /* LCOV_EXCL_LINE */
    2326             : }
    2327             : /* c t_FRAC */
    2328             : static GEN
    2329     6316591 : Q_divq_to_int(GEN x, GEN c)
    2330             : {
    2331     6316591 :   GEN n = gel(c,1), d = gel(c,2);
    2332     6316591 :   if (is_pm1(n)) {
    2333     4325175 :     GEN y = Q_muli_to_int(x,d);
    2334     4325175 :     if (signe(n) < 0) y = gneg(y);
    2335     4325175 :     return y;
    2336             :   }
    2337     1991416 :   return Q_divmuli_to_int(x, n,d);
    2338             : }
    2339             : 
    2340             : /* return y = x / c, assuming x,c rational, and y integral */
    2341             : GEN
    2342        4375 : Q_div_to_int(GEN x, GEN c)
    2343             : {
    2344        4375 :   switch(typ(c))
    2345             :   {
    2346        2954 :     case t_INT:  return Q_divi_to_int(x, c);
    2347        1421 :     case t_FRAC: return Q_divq_to_int(x, c);
    2348             :   }
    2349           0 :   pari_err_TYPE("Q_div_to_int",c);
    2350             :   return NULL; /* LCOV_EXCL_LINE */
    2351             : }
    2352             : /* return y = x * c, assuming x,c rational, and y integral */
    2353             : GEN
    2354           0 : Q_mul_to_int(GEN x, GEN c)
    2355             : {
    2356             :   GEN d, n;
    2357           0 :   switch(typ(c))
    2358             :   {
    2359           0 :     case t_INT: return Q_muli_to_int(x, c);
    2360             :     case t_FRAC:
    2361           0 :       n = gel(c,1);
    2362           0 :       d = gel(c,2);
    2363           0 :       return Q_divmuli_to_int(x, d,n);
    2364             :   }
    2365           0 :   pari_err_TYPE("Q_mul_to_int",c);
    2366             :   return NULL; /* LCOV_EXCL_LINE */
    2367             : }
    2368             : 
    2369             : GEN
    2370    34079156 : Q_primitive_part(GEN x, GEN *ptc)
    2371             : {
    2372    34079156 :   pari_sp av = avma;
    2373    34079156 :   GEN c = Q_content_safe(x);
    2374    34079171 :   if (c)
    2375             :   {
    2376    34078822 :     if (typ(c) == t_INT)
    2377             :     {
    2378    27763652 :       if (equali1(c)) { set_avma(av); c = NULL; }
    2379     3682602 :       else if (signe(c)) x = Q_divi_to_int(x, c);
    2380             :     }
    2381     6315170 :     else x = Q_divq_to_int(x, c);
    2382             :   }
    2383    34079144 :   if (ptc) *ptc = c;
    2384    34079144 :   return x;
    2385             : }
    2386             : GEN
    2387     1889397 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    2388             : 
    2389             : GEN
    2390      105393 : vec_Q_primpart(GEN x)
    2391      105393 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
    2392             : 
    2393             : /*******************************************************************/
    2394             : /*                                                                 */
    2395             : /*                           SUBRESULTANT                          */
    2396             : /*                                                                 */
    2397             : /*******************************************************************/
    2398             : /* for internal use */
    2399             : GEN
    2400     4470496 : gdivexact(GEN x, GEN y)
    2401             : {
    2402             :   long i,lx;
    2403             :   GEN z;
    2404     4470496 :   if (gequal1(y)) return x;
    2405     4466617 :   if (typ(y) == t_POLMOD) return gmul(x, ginv(y));
    2406     4466540 :   switch(typ(x))
    2407             :   {
    2408             :     case t_INT:
    2409     3983962 :       if (typ(y)==t_INT) return diviiexact(x,y);
    2410          35 :       if (!signe(x)) return gen_0;
    2411           0 :       break;
    2412             :     case t_INTMOD:
    2413             :     case t_FFELT:
    2414        9639 :     case t_POLMOD: return gmul(x,ginv(y));
    2415             :     case t_POL:
    2416      468557 :       switch(typ(y))
    2417             :       {
    2418             :         case t_INTMOD:
    2419             :         case t_FFELT:
    2420         714 :         case t_POLMOD: return gmul(x,ginv(y));
    2421             :         case t_POL: { /* not stack-clean */
    2422             :           long v;
    2423       25845 :           if (varn(x)!=varn(y)) break;
    2424       24879 :           v = RgX_valrem(y,&y);
    2425       24879 :           if (v) x = RgX_shift_shallow(x,-v);
    2426       24879 :           if (!degpol(y)) { y = gel(y,2); break; }
    2427       22156 :           return RgX_div(x,y);
    2428             :         }
    2429             :       }
    2430      445687 :       return RgX_Rg_divexact(x, y);
    2431             :     case t_VEC: case t_COL: case t_MAT:
    2432        4382 :       lx = lg(x); z = new_chunk(lx);
    2433        4382 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    2434        4382 :       z[0] = x[0]; return z;
    2435             :   }
    2436           0 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    2437           0 :   return gdiv(x,y);
    2438             : }
    2439             : 
    2440             : static GEN
    2441      213517 : init_resultant(GEN x, GEN y)
    2442             : {
    2443      213517 :   long tx = typ(x), ty = typ(y), vx, vy;
    2444      213517 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    2445             :   {
    2446          14 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    2447          14 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    2448           0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    2449           0 :     return gen_1;
    2450             :   }
    2451      213503 :   if (tx!=t_POL) pari_err_TYPE("resultant",x);
    2452      213503 :   if (ty!=t_POL) pari_err_TYPE("resultant",y);
    2453      213503 :   if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
    2454      213426 :   vx = varn(x);
    2455      213426 :   vy = varn(y); if (vx == vy) return NULL;
    2456           7 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    2457             : }
    2458             : 
    2459             : /* x an RgX, y a scalar */
    2460             : static GEN
    2461           0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    2462             : {
    2463           0 :   *V = gpowgs(y,degpol(x)-1);
    2464           0 :   *U = gen_0; return gmul(y, *V);
    2465             : }
    2466             : 
    2467             : /* return 0 if the subresultant chain can be interrupted.
    2468             :  * Set u = NULL if the resultant is 0. */
    2469             : static int
    2470        7598 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    2471             : {
    2472        7598 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    2473             :   long du, dv, dr, degq;
    2474             : 
    2475        7598 :   if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
    2476        7598 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    2477        7479 :   du = degpol(*u);
    2478        7479 :   dv = degpol(*v);
    2479        7479 :   degq = du - dv;
    2480        7479 :   if (*um1 == gen_1)
    2481        2496 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    2482        4983 :   else if (*um1 == gen_0)
    2483        2239 :     u0 = gen_0;
    2484             :   else /* except in those 2 cases, um1 is an RgX */
    2485        2744 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2486             : 
    2487        7479 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2488        2496 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2489             :   else
    2490        4983 :     u0 = gsub(u0, gmul(q,*uze));
    2491             : 
    2492        7479 :   *um1 = *uze;
    2493        7479 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2494             : 
    2495        7479 :   *u = *v; c = *g; *g  = leading_coeff(*u);
    2496        7479 :   switch(degq)
    2497             :   {
    2498          63 :     case 0: break;
    2499             :     case 1:
    2500        5625 :       c = gmul(*h,c); *h = *g; break;
    2501             :     default:
    2502        1791 :       c = gmul(gpowgs(*h,degq), c);
    2503        1791 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2504             :   }
    2505        7479 :   if (typ(c) == t_POLMOD)
    2506             :   {
    2507         895 :     c = ginv(c);
    2508         895 :     *v = RgX_Rg_mul(r,c);
    2509         895 :     *uze = RgX_Rg_mul(*uze,c);
    2510             :   }
    2511             :   else
    2512             :   {
    2513        6584 :     *v  = RgX_Rg_divexact(r,c);
    2514        6584 :     *uze= RgX_Rg_divexact(*uze,c);
    2515             :   }
    2516        7479 :   if (both_odd(du, dv)) *signh = -*signh;
    2517        7479 :   return (dr > 3);
    2518             : }
    2519             : 
    2520             : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2521             : static GEN
    2522        2377 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2523             : {
    2524             :   pari_sp av, av2;
    2525        2377 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2526             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2527             : 
    2528        2377 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2529        2377 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2530        2377 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2531        2377 :   if (tx != t_POL) {
    2532           0 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2533           0 :     return scalar_res(y,x,V,U);
    2534             :   }
    2535        2377 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2536        2377 :   if (varn(x) != varn(y))
    2537           0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2538           0 :                                         : scalar_res(y,x,V,U);
    2539        2377 :   if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
    2540        2377 :   if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
    2541        2377 :   dx = degpol(x);
    2542        2377 :   dy = degpol(y);
    2543        2377 :   signh = 1;
    2544        2377 :   if (dx < dy)
    2545             :   {
    2546         886 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2547         886 :     if (both_odd(dx, dy)) signh = -signh;
    2548             :   }
    2549        2377 :   if (dy == 0)
    2550             :   {
    2551           0 :     *V = gpowgs(gel(y,2),dx-1);
    2552           0 :     *U = gen_0; return gmul(*V,gel(y,2));
    2553             :   }
    2554        2377 :   av = avma;
    2555        2377 :   u = x = primitive_part(x, &cu);
    2556        2377 :   v = y = primitive_part(y, &cv);
    2557        2377 :   g = h = gen_1; av2 = avma;
    2558        2377 :   um1 = gen_1; uze = gen_0;
    2559             :   for(;;)
    2560             :   {
    2561       11657 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2562        4640 :     if (gc_needed(av2,1))
    2563             :     {
    2564           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2565           0 :       gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    2566             :     }
    2567             :   }
    2568             :   /* uze an RgX */
    2569        2377 :   if (!u) { *U = *V = gen_0; set_avma(av); return gen_0; }
    2570        2370 :   z = gel(v,2); du = degpol(u);
    2571        2370 :   if (du > 1)
    2572             :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2573         259 :     p1 = gpowgs(gdiv(z,h),du-1);
    2574         259 :     z = gmul(z,p1);
    2575         259 :     uze = RgX_Rg_mul(uze, p1);
    2576             :   }
    2577        2370 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2578             : 
    2579        2370 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2580        2370 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2581             :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2582        2370 :   p1 = gen_1;
    2583        2370 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2584        2370 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2585        2370 :   cu = cu? gdiv(p1,cu): p1;
    2586        2370 :   cv = cv? gdiv(p1,cv): p1;
    2587        2370 :   z = gmul(z,p1);
    2588        2370 :   *U = RgX_Rg_mul(uze,cu);
    2589        2370 :   *V = RgX_Rg_mul(vze,cv);
    2590        2370 :   return z;
    2591             : }
    2592             : GEN
    2593           0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2594             : {
    2595           0 :   pari_sp av = avma;
    2596           0 :   GEN z = subresext_i(x, y, U, V);
    2597           0 :   gerepileall(av, 3, &z, U, V);
    2598           0 :   return z;
    2599             : }
    2600             : 
    2601             : static GEN
    2602          28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2603             : {
    2604          28 :   GEN x=content(y);
    2605          28 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2606             : }
    2607             : 
    2608             : static int
    2609        2345 : must_negate(GEN x)
    2610             : {
    2611        2345 :   GEN t = leading_coeff(x);
    2612        2345 :   switch(typ(t))
    2613             :   {
    2614             :     case t_INT: case t_REAL:
    2615         469 :       return (signe(t) < 0);
    2616             :     case t_FRAC:
    2617           0 :       return (signe(gel(t,1)) < 0);
    2618             :   }
    2619        1876 :   return 0;
    2620             : }
    2621             : 
    2622             : static GEN
    2623          63 : gc_gcdext(pari_sp av, GEN r, GEN *u, GEN *v)
    2624             : {
    2625          63 :   if (!u && !v) return gerepileupto(av, r);
    2626          63 :   if (u  &&  v) gerepileall(av, 3, &r, u, v);
    2627           0 :   else          gerepileall(av, 2, &r, u ? u: v);
    2628          63 :   return r;
    2629             : }
    2630             : 
    2631             : static GEN
    2632          56 : RgX_extgcd_FpX(GEN x, GEN y, GEN p, GEN *u, GEN *v)
    2633             : {
    2634          56 :   pari_sp av = avma;
    2635          56 :   GEN r = FpX_extgcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, u, v);
    2636          56 :   if (u) *u = FpX_to_mod(*u, p);
    2637          56 :   if (v) *v = FpX_to_mod(*v, p);
    2638          56 :   return gc_gcdext(av, FpX_to_mod(r, p), u, v);
    2639             : }
    2640             : 
    2641             : static GEN
    2642           7 : RgX_extgcd_FpXQX(GEN x, GEN y, GEN pol, GEN p, GEN *U, GEN *V)
    2643             : {
    2644           7 :   pari_sp av = avma;
    2645           7 :   GEN r, T = RgX_to_FpX(pol, p);
    2646           7 :   r = FpXQX_extgcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p, U, V);
    2647           7 :   return gc_gcdext(av, FpXQX_to_mod(r, T, p), U, V);
    2648             : }
    2649             : 
    2650             : static GEN
    2651         378 : RgX_extgcd_fast(GEN x, GEN y, GEN *U, GEN *V)
    2652             : {
    2653             :   GEN p, pol;
    2654             :   long pa;
    2655         378 :   long t = RgX_type(x, &p,&pol,&pa);
    2656         378 :   switch(t)
    2657             :   {
    2658          21 :     case t_FFELT:  return FFX_extgcd(x, y, pol, U, V);
    2659          56 :     case t_INTMOD: return RgX_extgcd_FpX(x, y, p, U, V);
    2660             :     case code(t_POLMOD, t_INTMOD):
    2661           7 :                    return RgX_extgcd_FpXQX(x, y, pol, p, U, V);
    2662         294 :     default:       return NULL;
    2663             :   }
    2664             : }
    2665             : 
    2666             : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2667             : GEN
    2668         413 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2669             : {
    2670             :   pari_sp av, av2, tetpil;
    2671             :   long signh; /* junk */
    2672         413 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2673             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
    2674             : 
    2675         413 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2676         413 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2677         413 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2678         413 :   vx=varn(x);
    2679         413 :   if (!signe(x))
    2680             :   {
    2681          14 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2682           7 :     *U = pol_0(vx); *V = pol_0(vx);
    2683           7 :     return pol_0(vx);
    2684             :   }
    2685         399 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2686         378 :   r = RgX_extgcd_fast(x, y, U, V);
    2687         378 :   if (r) return r;
    2688         294 :   dx = degpol(x); dy = degpol(y);
    2689         294 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2690         294 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2691             : 
    2692         168 :   av = avma;
    2693         168 :   u = x = primitive_part(x, &cu);
    2694         168 :   v = y = primitive_part(y, &cv);
    2695         168 :   g = h = gen_1; av2 = avma;
    2696         168 :   um1 = gen_1; uze = gen_0;
    2697             :   for(;;)
    2698             :   {
    2699         420 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2700         126 :     if (gc_needed(av2,1))
    2701             :     {
    2702           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2703           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2704             :     }
    2705             :   }
    2706         168 :   if (uze != gen_0) {
    2707             :     GEN r;
    2708          56 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2709          56 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2710          56 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2711          56 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2712          56 :     p1 = ginv(content(v));
    2713             :   }
    2714             :   else /* y | x */
    2715             :   {
    2716         112 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2717         112 :     uze = pol_0(vx);
    2718         112 :     p1 = gen_1;
    2719             :   }
    2720         168 :   if (must_negate(v)) p1 = gneg(p1);
    2721         168 :   tetpil = avma;
    2722         168 :   z = RgX_Rg_mul(v,p1);
    2723         168 :   *U = RgX_Rg_mul(uze,p1);
    2724         168 :   *V = RgX_Rg_mul(vze,p1);
    2725         168 :   gptr[0] = &z;
    2726         168 :   gptr[1] = U;
    2727         168 :   gptr[2] = V;
    2728         168 :   gerepilemanysp(av,tetpil,gptr,3); return z;
    2729             : }
    2730             : 
    2731             : int
    2732          70 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2733             : {
    2734          70 :   pari_sp av = avma, av2, tetpil;
    2735             :   long signh; /* junk */
    2736             :   long vx;
    2737             :   GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
    2738             : 
    2739          70 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2740          70 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2741          70 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2742          70 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2743          70 :   if (!signe(T)) {
    2744           0 :     if (degpol(x) <= amax) {
    2745           0 :       *P = RgX_copy(x);
    2746           0 :       *Q = pol_1(varn(x));
    2747           0 :       return 1;
    2748             :     }
    2749           0 :     return 0;
    2750             :   }
    2751          70 :   if (amax+bmax >= degpol(T))
    2752           0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2753             :                     mkvec3(stoi(amax), stoi(bmax), T));
    2754          70 :   vx = varn(T);
    2755          70 :   u = x = primitive_part(x, &cu);
    2756          70 :   v = T = primitive_part(T, &cv);
    2757          70 :   g = h = gen_1; av2 = avma;
    2758          70 :   um1 = gen_1; uze = gen_0;
    2759             :   for(;;)
    2760             :   {
    2761         504 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2762         287 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) return gc_bool(av,0);
    2763         287 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2764         217 :     if (gc_needed(av2,1))
    2765             :     {
    2766           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2767           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2768             :     }
    2769             :   }
    2770          70 :   if (uze == gen_0)
    2771             :   {
    2772           0 :     set_avma(av); *P = pol_0(vx); *Q = pol_1(vx);
    2773           0 :     return 1;
    2774             :   }
    2775          70 :   if (cu) uze = RgX_Rg_div(uze,cu);
    2776          70 :   p1 = ginv(content(v));
    2777          70 :   if (must_negate(v)) p1 = gneg(p1);
    2778          70 :   tetpil = avma;
    2779          70 :   *P = RgX_Rg_mul(v,p1);
    2780          70 :   *Q = RgX_Rg_mul(uze,p1);
    2781          70 :   gptr[0] = P;
    2782          70 :   gptr[1] = Q;
    2783          70 :   gerepilemanysp(av,tetpil,gptr,2); return 1;
    2784             : }
    2785             : 
    2786             : /*******************************************************************/
    2787             : /*                                                                 */
    2788             : /*                RESULTANT USING DUCOS VARIANT                    */
    2789             : /*                                                                 */
    2790             : /*******************************************************************/
    2791             : /* x^n / y^(n-1), assume n > 0 */
    2792             : static GEN
    2793       23440 : Lazard(GEN x, GEN y, long n)
    2794             : {
    2795             :   long a;
    2796             :   GEN c;
    2797             : 
    2798       23440 :   if (n == 1) return x;
    2799        1654 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    2800        1654 :   c=x; n-=a;
    2801        5606 :   while (a>1)
    2802             :   {
    2803        2298 :     a>>=1; c=gdivexact(gsqr(c),y);
    2804        2298 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    2805             :   }
    2806        1654 :   return c;
    2807             : }
    2808             : 
    2809             : /* F (x/y)^(n-1), assume n >= 1 */
    2810             : static GEN
    2811       25075 : Lazard2(GEN F, GEN x, GEN y, long n)
    2812             : {
    2813       25075 :   if (n == 1) return F;
    2814         840 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    2815             : }
    2816             : 
    2817             : static GEN
    2818       25075 : RgX_neg_i(GEN x, long lx)
    2819             : {
    2820             :   long i;
    2821       25075 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    2822       25075 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    2823       25075 :   return y;
    2824             : }
    2825             : static GEN
    2826       75082 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    2827             : {
    2828             :   long i;
    2829             :   GEN z;
    2830       75082 :   if (isrationalzero(x)) return pol_0(varn(y));
    2831       75061 :   z = cgetg(ly,t_POL); z[1] = y[1];
    2832       75061 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2833       75061 :   return z;
    2834             : }
    2835             : static long
    2836       71925 : reductum_lg(GEN x, long lx)
    2837             : {
    2838       71925 :   long i = lx-2;
    2839       71925 :   while (i > 1 && gequal0(gel(x,i))) i--;
    2840       71925 :   return i+1;
    2841             : }
    2842             : 
    2843             : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
    2844             : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    2845             :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    2846             : static GEN
    2847       25075 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    2848             : {
    2849       25075 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
    2850             :   long p, q, j, lP, lQ;
    2851             :   pari_sp av;
    2852             : 
    2853       25075 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    2854       25075 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    2855             :   /* p > q. Very often p - 1 = q */
    2856       25075 :   av = avma;
    2857             :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    2858       25075 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    2859             : 
    2860       25075 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    2861       29016 :   for (j = q+1; j < p; j++)
    2862             :   {
    2863        3941 :     if (degpol(H) == q-1)
    2864             :     { /* h0 = coeff of degree q-1 = leading coeff */
    2865        3486 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    2866        3486 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    2867             :     }
    2868             :     else
    2869         455 :       H = RgX_shift_shallow(H, 1);
    2870        3941 :     if (j+2 < lP)
    2871             :     {
    2872        3150 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    2873        3150 :       A = A? RgX_add(A, TMP): TMP;
    2874             :     }
    2875        3941 :     if (gc_needed(av,1))
    2876             :     {
    2877         147 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    2878         147 :       gerepileall(av,A?2:1,&H,&A);
    2879             :     }
    2880             :   }
    2881       25075 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    2882       25075 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    2883       25075 :   A = A? RgX_add(A, TMP): TMP;
    2884       25075 :   A = RgX_Rg_divexact(A, p0);
    2885       25075 :   if (degpol(H) == q-1)
    2886             :   {
    2887       24746 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    2888       24746 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    2889             :   }
    2890             :   else
    2891         329 :     A = RgX_Rg_mul(addshift(H,A), q0);
    2892       25075 :   return RgX_Rg_divexact(A, s);
    2893             : }
    2894             : #undef addshift
    2895             : 
    2896             : /* Ducos's subresultant */
    2897             : GEN
    2898       24119 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    2899             : {
    2900             :   pari_sp av, av2;
    2901       24119 :   long dP, dQ, delta, sig = 1;
    2902             :   GEN cP, cQ, Z, s;
    2903             : 
    2904       24119 :   dP = degpol(P);
    2905       24119 :   dQ = degpol(Q); delta = dP - dQ;
    2906       24119 :   if (delta < 0)
    2907             :   {
    2908        1288 :     if (both_odd(dP, dQ)) sig = -1;
    2909        1288 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    2910             :   }
    2911       24119 :   if (sol) *sol = gen_0;
    2912       24119 :   av = avma;
    2913       24119 :   if (dQ <= 0)
    2914             :   {
    2915        1519 :     if (dQ < 0) return Rg_get_0(P);
    2916        1519 :     s = gpowgs(gel(Q,2), dP);
    2917        1519 :     if (sig == -1) s = gerepileupto(av, gneg(s));
    2918        1519 :     return s;
    2919             :   }
    2920             :   /* primitive_part is also possible here, but possibly very costly,
    2921             :    * and hardly ever worth it */
    2922       22600 :   P = Q_primitive_part(P, &cP);
    2923       22600 :   Q = Q_primitive_part(Q, &cQ);
    2924       22600 :   av2 = avma;
    2925       22600 :   s = gpowgs(leading_coeff(Q),delta);
    2926       22600 :   if (both_odd(dP, dQ)) sig = -sig;
    2927       22600 :   Z = Q;
    2928       22600 :   Q = RgX_pseudorem(P, Q);
    2929       22600 :   P = Z;
    2930       70275 :   while(degpol(Q) > 0)
    2931             :   {
    2932       25075 :     delta = degpol(P) - degpol(Q); /* > 0 */
    2933       25075 :     Z = Lazard2(Q, leading_coeff(Q), s, delta);
    2934       25075 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    2935       25075 :     Q = nextSousResultant(P, Q, Z, s);
    2936       25075 :     P = Z;
    2937       25075 :     if (gc_needed(av,1))
    2938             :     {
    2939          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    2940          13 :       gerepileall(av2,2,&P,&Q);
    2941             :     }
    2942       25075 :     s = leading_coeff(P);
    2943             :   }
    2944       22600 :   if (!signe(Q)) { set_avma(av); return Rg_get_0(Q); }
    2945       22600 :   s = Lazard(leading_coeff(Q), s, degpol(P));
    2946       22600 :   if (sig == -1) s = gneg(s);
    2947       22600 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    2948       22600 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    2949       22600 :   if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
    2950       20542 :   return gerepilecopy(av, s);
    2951             : }
    2952             : 
    2953             : static GEN
    2954          14 : RgX_resultant_FpX(GEN x, GEN y, GEN p)
    2955             : {
    2956          14 :   pari_sp av = avma;
    2957             :   GEN r;
    2958          14 :   if (lgefint(p) == 3)
    2959             :   {
    2960           7 :     ulong pp = uel(p, 2);
    2961           7 :     r = utoi(Flx_resultant(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), pp));
    2962             :   }
    2963             :   else
    2964           7 :     r = FpX_resultant(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    2965          14 :   return gerepileupto(av, Fp_to_mod(r, p));
    2966             : }
    2967             : 
    2968             : static GEN
    2969          21 : RgX_resultant_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    2970             : {
    2971          21 :   pari_sp av = avma;
    2972          21 :   GEN r, T = RgX_to_FpX(pol, p);
    2973          21 :   r = FpXQX_resultant(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    2974          21 :   return gerepileupto(av, FpX_to_mod(r, p));
    2975             : }
    2976             : 
    2977             : static GEN
    2978      213496 : resultant_fast(GEN x, GEN y)
    2979             : {
    2980             :   GEN p, pol;
    2981             :   long pa, t;
    2982      213496 :   p = init_resultant(x,y);
    2983      213496 :   if (p) return p;
    2984      213398 :   t = RgX_type2(x,y, &p,&pol,&pa);
    2985      213398 :   switch(t)
    2986             :   {
    2987       21144 :     case t_INT:    return ZX_resultant(x,y);
    2988       12641 :     case t_FRAC:   return QX_resultant(x,y);
    2989          21 :     case t_FFELT:  return FFX_resultant(x,y,pol);
    2990          14 :     case t_INTMOD: return RgX_resultant_FpX(x, y, p);
    2991             :     case code(t_POLMOD, t_INTMOD):
    2992          21 :                    return RgX_resultant_FpXQX(x, y, pol, p);
    2993      179557 :     default:       return NULL;
    2994             :   }
    2995             : }
    2996             : 
    2997             : static GEN
    2998      158651 : RgX_resultant_sylvester(GEN x, GEN y)
    2999             : {
    3000      158651 :   pari_sp av = avma;
    3001      158651 :   return gerepileupto(av, det(RgX_sylvestermatrix(x,y)));
    3002             : }
    3003             : 
    3004             : /* Return resultant(P,Q).
    3005             :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    3006             :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    3007             :  * in the "generic" case. */
    3008             : GEN
    3009      213496 : resultant(GEN P, GEN Q)
    3010             : {
    3011      213496 :   GEN z = resultant_fast(P,Q);
    3012      213496 :   if (z) return z;
    3013      179557 :   if (isinexact(P) || isinexact(Q)) return RgX_resultant_sylvester(P,Q);
    3014       20927 :   return RgX_resultant_all(P, Q, NULL);
    3015             : }
    3016             : 
    3017             : /*******************************************************************/
    3018             : /*                                                                 */
    3019             : /*               RESULTANT USING SYLVESTER MATRIX                  */
    3020             : /*                                                                 */
    3021             : /*******************************************************************/
    3022             : static GEN
    3023      216189 : syl_RgC(GEN x, long j, long d, long D, long cp)
    3024             : {
    3025      216189 :   GEN c = cgetg(d+1,t_COL);
    3026             :   long i;
    3027      216189 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    3028      216189 :   for (   ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
    3029      216189 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    3030      216189 :   return c;
    3031             : }
    3032             : static GEN
    3033      158658 : syl_RgM(GEN x, GEN y, long cp)
    3034             : {
    3035      158658 :   long j, d, dx = degpol(x), dy = degpol(y);
    3036             :   GEN M;
    3037      158658 :   if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
    3038      158658 :   if (dy < 0) return zeromat(dx,dx);
    3039      158658 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    3040      158658 :   for (j=1; j<=dy; j++) gel(M,j)    = syl_RgC(x,j,d,j+dx, cp);
    3041      158658 :   for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
    3042      158658 :   return M;
    3043             : }
    3044             : GEN
    3045      158651 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
    3046             : GEN
    3047           7 : sylvestermatrix(GEN x, GEN y)
    3048             : {
    3049           7 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    3050           7 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    3051           7 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    3052           7 :   return syl_RgM(x,y,1);
    3053             : }
    3054             : 
    3055             : GEN
    3056          21 : resultant2(GEN x, GEN y)
    3057             : {
    3058          21 :   GEN r = init_resultant(x,y);
    3059          21 :   return r? r: RgX_resultant_sylvester(x,y);
    3060             : }
    3061             : 
    3062             : /* let vx = main variable of x, v0 a variable of highest priority;
    3063             :  * return a t_POL in variable v0:
    3064             :  * if vx <= v, return subst(x, v, pol_x(v0))
    3065             :  * if vx >  v, return scalarpol(x, v0) */
    3066             : static GEN
    3067        3717 : fix_pol(GEN x, long v, long v0)
    3068             : {
    3069        3717 :   long vx, tx = typ(x);
    3070        3717 :   if (tx != t_POL)
    3071          21 :     vx = gvar(x);
    3072             :   else
    3073             :   { /* shortcut: almost nothing to do */
    3074        3696 :     vx = varn(x);
    3075        3696 :     if (v == vx)
    3076             :     {
    3077        3535 :       if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
    3078        3535 :       return x;
    3079             :     }
    3080             :   }
    3081         182 :   if (varncmp(v, vx) > 0)
    3082             :   {
    3083         175 :     x = gsubst(x, v, pol_x(v0));
    3084         175 :     if (typ(x) != t_POL) vx = gvar(x);
    3085             :     else
    3086             :     {
    3087         168 :       vx = varn(x);
    3088         168 :       if (vx == v0) return x;
    3089             :     }
    3090             :   }
    3091          49 :   if (varncmp(vx, v0) <= 0) pari_err_TYPE("polresultant", x);
    3092          42 :   return scalarpol_shallow(x, v0);
    3093             : }
    3094             : 
    3095             : /* resultant of x and y with respect to variable v, or with respect to their
    3096             :  * main variable if v < 0. */
    3097             : GEN
    3098        2170 : polresultant0(GEN x, GEN y, long v, long flag)
    3099             : {
    3100        2170 :   pari_sp av = avma;
    3101             : 
    3102        2170 :   if (v >= 0)
    3103             :   {
    3104        1834 :     long v0 = fetch_var_higher();
    3105        1834 :     x = fix_pol(x,v, v0);
    3106        1834 :     y = fix_pol(y,v, v0);
    3107             :   }
    3108        2170 :   switch(flag)
    3109             :   {
    3110             :     case 2:
    3111        2163 :     case 0: x=resultant(x,y); break;
    3112           7 :     case 1: x=resultant2(x,y); break;
    3113           0 :     default: pari_err_FLAG("polresultant");
    3114             :   }
    3115        2170 :   if (v >= 0) (void)delete_var();
    3116        2170 :   return gerepileupto(av,x);
    3117             : }
    3118             : GEN
    3119        1491 : polresultantext0(GEN x, GEN y, long v)
    3120             : {
    3121             :   GEN R, U, V;
    3122        1491 :   pari_sp av = avma;
    3123             : 
    3124        1491 :   if (v >= 0)
    3125             :   {
    3126          14 :     long v0 = fetch_var_higher();
    3127          14 :     x = fix_pol(x,v, v0);
    3128          14 :     y = fix_pol(y,v, v0);
    3129             :   }
    3130        1491 :   R = subresext_i(x,y, &U,&V);
    3131        1491 :   if (v >= 0)
    3132             :   {
    3133          14 :     (void)delete_var();
    3134          14 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    3135          14 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    3136             :   }
    3137        1491 :   return gerepilecopy(av, mkvec3(U,V,R));
    3138             : }
    3139             : GEN
    3140        1463 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    3141             : 
    3142             : /*******************************************************************/
    3143             : /*                                                                 */
    3144             : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    3145             : /*                                                                 */
    3146             : /*******************************************************************/
    3147             : 
    3148             : /* (v - x)^d */
    3149             : static GEN
    3150         126 : caract_const(pari_sp av, GEN x, long v, long d)
    3151         126 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    3152             : 
    3153             : /* return caract(Mod(x,T)) in variable v */
    3154             : GEN
    3155       17637 : RgXQ_charpoly(GEN x, GEN T, long v)
    3156             : {
    3157       17637 :   pari_sp av = avma;
    3158       17637 :   long d = degpol(T), dx, vx, vp, v0;
    3159             :   GEN ch, L;
    3160             : 
    3161       17637 :   if (typ(x) != t_POL) return caract_const(av, x, v, d);
    3162       17637 :   vx = varn(x);
    3163       17637 :   vp = varn(T);
    3164       17637 :   if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
    3165       17637 :   if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
    3166       17637 :   dx = degpol(x);
    3167       17637 :   if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
    3168       17637 :   if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
    3169             : 
    3170       17581 :   v0 = fetch_var_higher();
    3171       17581 :   x = RgX_neg(x);
    3172       17581 :   gel(x,2) = gadd(gel(x,2), pol_x(v));
    3173       17581 :   setvarn(x, v0);
    3174       17581 :   T = leafcopy(T); setvarn(T, v0);
    3175       17581 :   ch = resultant(T, x);
    3176       17581 :   (void)delete_var();
    3177             :   /* test for silly input: x mod (deg 0 polynomial) */
    3178       17581 :   if (typ(ch) != t_POL)
    3179           7 :     pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
    3180       17574 :   L = leading_coeff(ch);
    3181       17574 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    3182       17574 :   return gerepileupto(av, ch);
    3183             : }
    3184             : 
    3185             : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    3186             :  * algebra nf[t]/(Q(t)) */
    3187             : GEN
    3188         224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    3189             : {
    3190         224 :   const char *f = "rnfcharpoly";
    3191         224 :   long dQ = degpol(Q);
    3192         224 :   pari_sp av = avma;
    3193             :   GEN T;
    3194             : 
    3195         224 :   if (v < 0) v = 0;
    3196         224 :   nf = checknf(nf); T = nf_get_pol(nf);
    3197         224 :   Q = RgX_nffix(f, T,Q,0);
    3198         224 :   switch(typ(x))
    3199             :   {
    3200             :     case t_INT:
    3201          28 :     case t_FRAC: return caract_const(av, x, v, dQ);
    3202             :     case t_POLMOD:
    3203          91 :       x = polmod_nffix2(f,T,Q, x,0);
    3204          56 :       break;
    3205             :     case t_POL:
    3206          56 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    3207          42 :       break;
    3208          49 :     default: pari_err_TYPE(f,x);
    3209             :   }
    3210          98 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    3211             :   /* x a t_POL in variable vQ */
    3212          56 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    3213          56 :   if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
    3214          56 :   return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    3215             : }
    3216             : 
    3217             : /*******************************************************************/
    3218             : /*                                                                 */
    3219             : /*                  GCD USING SUBRESULTANT                         */
    3220             : /*                                                                 */
    3221             : /*******************************************************************/
    3222             : static int inexact(GEN x, int *simple);
    3223             : static int
    3224       37821 : isinexactall(GEN x, int *simple)
    3225             : {
    3226       37821 :   long i, lx = lg(x);
    3227      136850 :   for (i=2; i<lx; i++)
    3228       99043 :     if (inexact(gel(x,i), simple)) return 1;
    3229       37807 :   return 0;
    3230             : }
    3231             : /* return 1 if coeff explosion is not possible */
    3232             : static int
    3233       99253 : inexact(GEN x, int *simple)
    3234             : {
    3235       99253 :   int junk = 0;
    3236       99253 :   switch(typ(x))
    3237             :   {
    3238       63889 :     case t_INT: case t_FRAC: return 0;
    3239             : 
    3240           7 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    3241             : 
    3242             :     case t_INTMOD:
    3243             :     case t_FFELT:
    3244        3654 :       if (!*simple) *simple = 1;
    3245        3654 :       return 0;
    3246             : 
    3247             :     case t_COMPLEX:
    3248          77 :       return inexact(gel(x,1), simple)
    3249          77 :           || inexact(gel(x,2), simple);
    3250             :     case t_QUAD:
    3251           0 :       *simple = 0;
    3252           0 :       return inexact(gel(x,2), &junk)
    3253           0 :           || inexact(gel(x,3), &junk);
    3254             : 
    3255             :     case t_POLMOD:
    3256         819 :       return isinexactall(gel(x,1), simple);
    3257             :     case t_POL:
    3258       30779 :       *simple = -1;
    3259       30779 :       return isinexactall(x, &junk);
    3260             :     case t_RFRAC:
    3261          28 :       *simple = -1;
    3262          28 :       return inexact(gel(x,1), &junk)
    3263          28 :           || inexact(gel(x,2), &junk);
    3264             :   }
    3265           0 :   *simple = -1; return 0;
    3266             : }
    3267             : 
    3268             : /* x monomial, y t_POL in the same variable */
    3269             : static GEN
    3270     1706251 : gcdmonome(GEN x, GEN y)
    3271             : {
    3272     1706251 :   pari_sp av = avma;
    3273     1706251 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    3274     1706251 :   long i, l = lg(y);
    3275     1706251 :   GEN t, v = cgetg(l, t_VEC);
    3276     1706251 :   gel(v,1) = gel(x,dx+2);
    3277     1706251 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    3278     1706251 :   t = content(v); /* gcd(lc(x), cont(y)) */
    3279     1706251 :   t = simplify_shallow(t);
    3280     1706251 :   if (dx < e) e = dx;
    3281     1706251 :   return gerepileupto(av, monomialcopy(t, e, varn(x)));
    3282             : }
    3283             : 
    3284             : static GEN
    3285      198226 : RgX_gcd_FpX(GEN x, GEN y, GEN p)
    3286             : {
    3287      198226 :   pari_sp av = avma;
    3288             :   GEN r;
    3289      198226 :   if (lgefint(p) == 3)
    3290             :   {
    3291      198212 :     ulong pp = uel(p, 2);
    3292      198212 :     r = Flx_to_ZX_inplace(Flx_gcd(RgX_to_Flx(x, pp),
    3293             :                                   RgX_to_Flx(y, pp), pp));
    3294             :   }
    3295             :   else
    3296          14 :     r = FpX_gcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3297      198226 :   return gerepileupto(av, FpX_to_mod(r, p));
    3298             : }
    3299             : 
    3300             : static GEN
    3301           7 : RgX_gcd_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3302             : {
    3303           7 :   pari_sp av = avma;
    3304           7 :   GEN r, T = RgX_to_FpX(pol, p);
    3305           7 :   if (signe(T)==0) pari_err_OP("gcd", x, y);
    3306           7 :   r = FpXQX_gcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3307           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3308             : }
    3309             : 
    3310             : static GEN
    3311        7112 : RgX_liftred(GEN x, GEN T)
    3312        7112 : { return RgXQX_red(liftpol_shallow(x), T); }
    3313             : 
    3314             : static GEN
    3315        2268 : RgX_gcd_ZXQX(GEN x, GEN y, GEN T)
    3316             : {
    3317        2268 :   pari_sp av = avma;
    3318        2268 :   GEN r = ZXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3319        2268 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    3320             : }
    3321             : 
    3322             : static GEN
    3323        1288 : RgX_gcd_QXQX(GEN x, GEN y, GEN T)
    3324             : {
    3325        1288 :   pari_sp av = avma;
    3326        1288 :   GEN r = QXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3327        1288 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    3328             : }
    3329             : 
    3330             : static GEN
    3331    11328433 : RgX_gcd_fast(GEN x, GEN y)
    3332             : {
    3333             :   GEN p, pol;
    3334             :   long pa;
    3335    11328433 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3336    11328433 :   switch(t)
    3337             :   {
    3338     9412301 :     case t_INT:    return ZX_gcd(x, y);
    3339        2429 :     case t_FRAC:   return QX_gcd(x, y);
    3340        2520 :     case t_FFELT:  return FFX_gcd(x, y, pol);
    3341      198226 :     case t_INTMOD: return RgX_gcd_FpX(x, y, p);
    3342             :     case code(t_POLMOD, t_INTMOD):
    3343           7 :                    return RgX_gcd_FpXQX(x, y, pol, p);
    3344             :     case code(t_POLMOD, t_INT):
    3345        2275 :                    return ZX_is_monic(pol)? RgX_gcd_ZXQX(x,y,pol): NULL;
    3346             :     case code(t_POLMOD, t_FRAC):
    3347        2604 :                    return RgX_is_ZX(pol) && ZX_is_monic(pol) ?
    3348        2590 :                                             RgX_gcd_QXQX(x,y,pol): NULL;
    3349     1709373 :     default:       return NULL;
    3350             :   }
    3351             : }
    3352             : 
    3353             : /* x, y are t_POL in the same variable */
    3354             : GEN
    3355    11328433 : RgX_gcd(GEN x, GEN y)
    3356             : {
    3357             :   long dx, dy;
    3358             :   pari_sp av, av1;
    3359             :   GEN d, g, h, p1, p2, u, v;
    3360    11328433 :   int simple = 0;
    3361    11328433 :   GEN z = RgX_gcd_fast(x, y);
    3362    11328433 :   if (z) return z;
    3363     1709394 :   if (isexactzero(y)) return RgX_copy(x);
    3364     1709366 :   if (isexactzero(x)) return RgX_copy(y);
    3365     1709366 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    3366        4053 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    3367        3115 :   if (isinexactall(x,&simple) || isinexactall(y,&simple))
    3368             :   {
    3369           7 :     av = avma; u = ggcd(content(x), content(y));
    3370           7 :     return gerepileupto(av, scalarpol(u, varn(x)));
    3371             :   }
    3372             : 
    3373        3108 :   av = avma;
    3374        3108 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    3375             :   else
    3376             :   {
    3377        3108 :     dx = lg(x); dy = lg(y);
    3378        3108 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    3379        3108 :     if (dy==3)
    3380             :     {
    3381           0 :       d = ggcd(gel(y,2), content(x));
    3382           0 :       return gerepileupto(av, scalarpol(d, varn(x)));
    3383             :     }
    3384        3108 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    3385        3108 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    3386        3108 :     d = ggcd(p1,p2);
    3387        3108 :     av1 = avma;
    3388        3108 :     g = h = gen_1;
    3389             :     for(;;)
    3390        1190 :     {
    3391        4298 :       GEN r = RgX_pseudorem(u,v);
    3392        4298 :       long degq, du, dv, dr = lg(r);
    3393             : 
    3394        4298 :       if (!signe(r)) break;
    3395        2191 :       if (dr <= 3)
    3396             :       {
    3397        1001 :         set_avma(av1);
    3398        1001 :         return gerepileupto(av, scalarpol(d, varn(x)));
    3399             :       }
    3400        1190 :       du = lg(u); dv = lg(v); degq = du-dv;
    3401        1190 :       u = v; p1 = g; g = leading_coeff(u);
    3402        1190 :       switch(degq)
    3403             :       {
    3404         189 :         case 0: break;
    3405             :         case 1:
    3406         910 :           p1 = gmul(h,p1); h = g; break;
    3407             :         default:
    3408          91 :           p1 = gmul(gpowgs(h,degq), p1);
    3409          91 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    3410             :       }
    3411        1190 :       v = RgX_Rg_div(r,p1);
    3412        1190 :       if (gc_needed(av1,1))
    3413             :       {
    3414           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd, dr = %ld", degpol(r));
    3415           0 :         gerepileall(av1,4, &u,&v,&g,&h);
    3416             :       }
    3417             :     }
    3418        2107 :     x = RgX_Rg_mul(primpart(v), d);
    3419             :   }
    3420        2107 :   if (must_negate(x)) x = RgX_neg(x);
    3421        2107 :   return gerepileupto(av,x);
    3422             : }
    3423             : 
    3424             : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
    3425             : static GEN
    3426        1337 : RgX_disc_i(GEN P)
    3427             : {
    3428        1337 :   long n = degpol(P), dd;
    3429             :   GEN N, D, L, y;
    3430        1337 :   if (!signe(P) || !n) return Rg_get_0(P);
    3431        1330 :   if (n == 1) return Rg_get_1(P);
    3432        1323 :   if (n == 2) {
    3433         630 :     GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
    3434         630 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    3435             :   }
    3436         693 :   y = RgX_deriv(P);
    3437         693 :   N = characteristic(P);
    3438         693 :   if (signe(N)) y = gmul(y, mkintmod(gen_1,N));
    3439         693 :   if (!signe(y)) return Rg_get_0(y);
    3440         693 :   dd = n - 2 - degpol(y);
    3441         693 :   if (isinexact(P))
    3442          14 :     D = resultant2(P,y);
    3443             :   else
    3444             :   {
    3445         679 :     D = RgX_resultant_all(P, y, NULL);
    3446         679 :     if (D == gen_0) return Rg_get_0(y);
    3447             :   }
    3448         693 :   L = leading_coeff(P);
    3449         693 :   if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
    3450         693 :   if (n & 2) D = gneg(D);
    3451         693 :   return D;
    3452             : }
    3453             : 
    3454             : static GEN
    3455          42 : RgX_disc_FpX(GEN x, GEN p)
    3456             : {
    3457          42 :   pari_sp av = avma;
    3458          42 :   GEN r = FpX_disc(RgX_to_FpX(x, p), p);
    3459          42 :   return gerepileupto(av, Fp_to_mod(r, p));
    3460             : }
    3461             : 
    3462             : static GEN
    3463          28 : RgX_disc_FpXQX(GEN x, GEN pol, GEN p)
    3464             : {
    3465          28 :   pari_sp av = avma;
    3466          28 :   GEN r, T = RgX_to_FpX(pol, p);
    3467          28 :   r = FpXQX_disc(RgX_to_FpXQX(x, T, p), T, p);
    3468          28 :   return gerepileupto(av, FpX_to_mod(r, p));
    3469             : }
    3470             : 
    3471             : static GEN
    3472        6608 : RgX_disc_fast(GEN x)
    3473             : {
    3474             :   GEN p, pol;
    3475             :   long pa;
    3476        6608 :   long t = RgX_type(x, &p,&pol,&pa);
    3477        6609 :   switch(t)
    3478             :   {
    3479        5146 :     case t_INT:    return ZX_disc(x);
    3480          21 :     case t_FRAC:   return QX_disc(x);
    3481          35 :     case t_FFELT:  return FFX_disc(x, pol);
    3482          42 :     case t_INTMOD: return RgX_disc_FpX(x, p);
    3483             :     case code(t_POLMOD, t_INTMOD):
    3484          28 :                    return RgX_disc_FpXQX(x, pol, p);
    3485        1337 :     default:       return NULL;
    3486             :   }
    3487             : }
    3488             : 
    3489             : GEN
    3490        6608 : RgX_disc(GEN x)
    3491             : {
    3492             :   pari_sp av;
    3493        6608 :   GEN z = RgX_disc_fast(x);
    3494        6609 :   if (z) return z;
    3495        1337 :   av = avma;
    3496        1337 :   return gerepileupto(av, RgX_disc_i(x));
    3497             : }
    3498             : 
    3499             : GEN
    3500        3871 : poldisc0(GEN x, long v)
    3501             : {
    3502        3871 :   long v0, tx = typ(x);
    3503             :   pari_sp av;
    3504             :   GEN D;
    3505        3871 :   if (tx == t_POL && (v < 0 || v == varn(x))) return RgX_disc(x);
    3506          35 :   switch(tx)
    3507             :   {
    3508             :     case t_QUAD:
    3509           0 :       return quad_disc(x);
    3510             :     case t_POLMOD:
    3511           0 :       if (v >= 0 && varn(gel(x,1)) != v) break;
    3512           0 :       return RgX_disc(gel(x,1));
    3513             :     case t_QFR: case t_QFI:
    3514          14 :       av = avma; return gerepileuptoint(av, qfb_disc(x));
    3515             :     case t_VEC: case t_COL: case t_MAT:
    3516             :     {
    3517             :       long i;
    3518           0 :       GEN z = cgetg_copy(x, &i);
    3519           0 :       for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
    3520           0 :       return z;
    3521             :     }
    3522             :   }
    3523          21 :   if (v < 0) pari_err_TYPE("poldisc",x);
    3524          21 :   av = avma; v0 = fetch_var_higher();
    3525          21 :   x = fix_pol(x,v, v0);
    3526          14 :   D = RgX_disc(x); (void)delete_var();
    3527          14 :   return gerepileupto(av, D);
    3528             : }
    3529             : 
    3530             : GEN
    3531           7 : reduceddiscsmith(GEN x)
    3532             : {
    3533           7 :   long j, n = degpol(x);
    3534           7 :   pari_sp av = avma;
    3535             :   GEN xp, M;
    3536             : 
    3537           7 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    3538           7 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    3539           7 :   RgX_check_ZX(x,"poldiscreduced");
    3540           7 :   if (!gequal1(gel(x,n+2)))
    3541           0 :     pari_err_IMPL("non-monic polynomial in poldiscreduced");
    3542           7 :   M = cgetg(n+1,t_MAT);
    3543           7 :   xp = ZX_deriv(x);
    3544          28 :   for (j=1; j<=n; j++)
    3545             :   {
    3546          21 :     gel(M,j) = RgX_to_RgC(xp, n);
    3547          21 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    3548             :   }
    3549           7 :   return gerepileupto(av, ZM_snf(M));
    3550             : }
    3551             : 
    3552             : /***********************************************************************/
    3553             : /**                                                                   **/
    3554             : /**                       STURM ALGORITHM                             **/
    3555             : /**              (number of real roots of x in [a,b])                 **/
    3556             : /**                                                                   **/
    3557             : /***********************************************************************/
    3558             : static GEN
    3559         882 : R_to_Q_up(GEN x)
    3560             : {
    3561             :   long e;
    3562         882 :   switch(typ(x))
    3563             :   {
    3564         882 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3565             :     case t_REAL:
    3566           0 :       x = mantissa_real(x,&e);
    3567           0 :       return gmul2n(addiu(x,1), -e);
    3568           0 :     default: pari_err_TYPE("R_to_Q_up", x);
    3569             :       return NULL; /* LCOV_EXCL_LINE */
    3570             :   }
    3571             : }
    3572             : static GEN
    3573         882 : R_to_Q_down(GEN x)
    3574             : {
    3575             :   long e;
    3576         882 :   switch(typ(x))
    3577             :   {
    3578         868 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3579             :     case t_REAL:
    3580          14 :       x = mantissa_real(x,&e);
    3581          14 :       return gmul2n(subiu(x,1), -e);
    3582           0 :     default: pari_err_TYPE("R_to_Q_down", x);
    3583             :       return NULL; /* LCOV_EXCL_LINE */
    3584             :   }
    3585             : }
    3586             : 
    3587             : static long
    3588         882 : sturmpart_i(GEN x, GEN ab)
    3589             : {
    3590         882 :   long tx = typ(x);
    3591         882 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    3592         882 :   if (tx != t_POL)
    3593             :   {
    3594           0 :     if (is_real_t(tx)) return 0;
    3595           0 :     pari_err_TYPE("sturm",x);
    3596             :   }
    3597         882 :   if (lg(x) == 3) return 0;
    3598         882 :   if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
    3599         882 :   (void)ZX_gcd_all(x, ZX_deriv(x), &x);
    3600         882 :   if (ab)
    3601             :   {
    3602             :     GEN A, B;
    3603         882 :     if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
    3604         882 :     A = R_to_Q_down(gel(ab,1));
    3605         882 :     B = R_to_Q_up(gel(ab,2));
    3606         882 :     ab = mkvec2(A, B);
    3607             :   }
    3608         882 :   return ZX_sturmpart(x, ab);
    3609             : }
    3610             : /* Deprecated: RgX_sturmpart() should be preferred  */
    3611             : long
    3612         882 : sturmpart(GEN x, GEN a, GEN b)
    3613             : {
    3614         882 :   pari_sp av = avma;
    3615         882 :   if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
    3616         749 :   if (!a) a = mkmoo();
    3617         749 :   if (!b) b = mkoo();
    3618         749 :   return gc_long(av, sturmpart_i(x, mkvec2(a,b)));
    3619             : }
    3620             : long
    3621         133 : RgX_sturmpart(GEN x, GEN ab)
    3622         133 : { pari_sp av = avma; return gc_long(av, sturmpart_i(x, ab)); }
    3623             : 
    3624             : /***********************************************************************/
    3625             : /**                                                                   **/
    3626             : /**                        GENERIC EXTENDED GCD                       **/
    3627             : /**                                                                   **/
    3628             : /***********************************************************************/
    3629             : /* assume typ(x) = typ(y) = t_POL */
    3630             : static GEN
    3631         886 : RgXQ_inv_i(GEN x, GEN y)
    3632             : {
    3633         886 :   long vx=varn(x), vy=varn(y);
    3634             :   pari_sp av;
    3635             :   GEN u, v, d;
    3636             : 
    3637        1772 :   while (vx != vy)
    3638             :   {
    3639           0 :     if (varncmp(vx,vy) > 0)
    3640             :     {
    3641           0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3642           0 :       return scalarpol(d, vy);
    3643             :     }
    3644           0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3645           0 :     x = gel(x,2); vx = gvar(x);
    3646             :   }
    3647         886 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3648         886 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3649         886 :   d = gdiv(u,d);
    3650         886 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3651         886 :   return gerepileupto(av, d);
    3652             : }
    3653             : 
    3654             : /*Assume x is a polynomial and y is not */
    3655             : static GEN
    3656         112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3657             : {
    3658         112 :   long vx = varn(x);
    3659         112 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    3660         112 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    3661          84 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    3662          56 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    3663             : }
    3664             : /* Assume x==0, y!=0 */
    3665             : static GEN
    3666          63 : zero_bezout(GEN y, GEN *U, GEN *V)
    3667             : {
    3668          63 :   *U=gen_0; *V = ginv(y); return gen_1;
    3669             : }
    3670             : 
    3671             : GEN
    3672         350 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    3673             : {
    3674         350 :   long tx=typ(x), ty=typ(y), vx;
    3675         350 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    3676         315 :   if (tx != t_POL)
    3677             :   {
    3678         140 :     if (ty == t_POL)
    3679          56 :       return scalar_bezout(y,x,v,u);
    3680             :     else
    3681             :     {
    3682          84 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    3683          84 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    3684          63 :       if (xis0) return zero_bezout(y,u,v);
    3685          42 :       else return zero_bezout(x,v,u);
    3686             :     }
    3687             :   }
    3688         175 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    3689         119 :   vx = varn(x);
    3690         119 :   if (vx != varn(y))
    3691           0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    3692           0 :                                    : scalar_bezout(y,x,v,u);
    3693         119 :   return RgX_extgcd(x,y,u,v);
    3694             : }
    3695             : 
    3696             : GEN
    3697         350 : gcdext0(GEN x, GEN y)
    3698             : {
    3699         350 :   GEN z=cgetg(4,t_VEC);
    3700         350 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    3701         350 :   return z;
    3702             : }
    3703             : 
    3704             : /*******************************************************************/
    3705             : /*                                                                 */
    3706             : /*                    GENERIC (modular) INVERSE                    */
    3707             : /*                                                                 */
    3708             : /*******************************************************************/
    3709             : 
    3710             : GEN
    3711       16779 : ginvmod(GEN x, GEN y)
    3712             : {
    3713       16779 :   long tx=typ(x);
    3714             : 
    3715       16779 :   switch(typ(y))
    3716             :   {
    3717             :     case t_POL:
    3718       16779 :       if (tx==t_POL) return RgXQ_inv(x,y);
    3719        3813 :       if (is_scalar_t(tx)) return ginv(x);
    3720           0 :       break;
    3721             :     case t_INT:
    3722           0 :       if (tx==t_INT) return Fp_inv(x,y);
    3723           0 :       if (tx==t_POL) return gen_0;
    3724             :   }
    3725           0 :   pari_err_TYPE2("ginvmod",x,y);
    3726             :   return NULL; /* LCOV_EXCL_LINE */
    3727             : }
    3728             : 
    3729             : /***********************************************************************/
    3730             : /**                                                                   **/
    3731             : /**                         NEWTON POLYGON                            **/
    3732             : /**                                                                   **/
    3733             : /***********************************************************************/
    3734             : 
    3735             : /* assume leading coeff of x is non-zero */
    3736             : GEN
    3737          28 : newtonpoly(GEN x, GEN p)
    3738             : {
    3739          28 :   pari_sp av = avma;
    3740             :   long n, ind, a, b;
    3741             :   GEN y, vval;
    3742             : 
    3743          28 :   if (typ(x) != t_POL) pari_err_TYPE("newtonpoly",x);
    3744          28 :   n = degpol(x); if (n<=0) return cgetg(1,t_VEC);
    3745          28 :   vval = new_chunk(n+1);
    3746          28 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    3747          28 :   for (a = 0; a <= n; a++) vval[a] = gvaluation(gel(x,a),p);
    3748          42 :   for (a = 0, ind = 1; a < n; a++)
    3749             :   {
    3750          42 :     if (vval[a] != LONG_MAX) break;
    3751          14 :     gel(y,ind++) = mkoo();
    3752             :   }
    3753          84 :   for (b = a+1; b <= n; a = b, b = a+1)
    3754             :   {
    3755             :     long u1, u2, c;
    3756          56 :     while (vval[b] == LONG_MAX) b++;
    3757          56 :     u1 = vval[a] - vval[b];
    3758          56 :     u2 = b - a;
    3759         154 :     for (c = b+1; c <= n; c++)
    3760             :     {
    3761             :       long r1, r2;
    3762          98 :       if (vval[c] == LONG_MAX) continue;
    3763          70 :       r1 = vval[a] - vval[c];
    3764          70 :       r2 = c - a;
    3765          70 :       if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
    3766             :     }
    3767          56 :     while (ind <= b) gel(y,ind++) = sstoQ(u1,u2);
    3768             :   }
    3769          28 :   stackdummy((pari_sp)vval, av); return y;
    3770             : }
    3771             : 
    3772             : static GEN
    3773      274288 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
    3774             : {
    3775      274288 :   pari_sp av = avma;
    3776             :   GEN r;
    3777      274288 :   if (lgefint(p) == 3)
    3778             :   {
    3779      152381 :     ulong pp = uel(p, 2);
    3780      152381 :     r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
    3781             :                 RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
    3782             :   }
    3783             :   else
    3784      121907 :     r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
    3785      274288 :   return gerepileupto(av, FpX_to_mod(r, p));
    3786             : }
    3787             : 
    3788             : static GEN
    3789          14 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
    3790             : {
    3791          14 :   pari_sp av = avma;
    3792             :   GEN r;
    3793          14 :   if (lgefint(p) == 3)
    3794             :   {
    3795           7 :     ulong pp = uel(p, 2);
    3796           7 :     r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
    3797             :                                    RgX_to_Flx(y, pp), pp));
    3798             :   }
    3799             :   else
    3800           7 :     r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3801          14 :   return gerepileupto(av, FpX_to_mod(r, p));
    3802             : }
    3803             : 
    3804             : static GEN
    3805       12054 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
    3806             : {
    3807       12054 :   pari_sp av = avma;
    3808             :   GEN r;
    3809       12054 :   if (lgefint(p) == 3)
    3810             :   {
    3811        6088 :     ulong pp = uel(p, 2);
    3812        6088 :     r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
    3813             :                                    RgX_to_Flx(y, pp), pp));
    3814             :   }
    3815             :   else
    3816        5966 :     r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3817       12054 :   return gerepileupto(av, FpX_to_mod(r, p));
    3818             : }
    3819             : 
    3820             : static GEN
    3821         385 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
    3822             : {
    3823         385 :   pari_sp av = avma;
    3824             :   GEN r;
    3825         385 :   GEN T = RgX_to_FpX(pol, p);
    3826         385 :   if (signe(T)==0) pari_err_OP("*",x,y);
    3827         385 :   if (lgefint(p) == 3)
    3828             :   {
    3829         241 :     ulong pp = uel(p, 2);
    3830         241 :     GEN Tp = ZX_to_Flx(T, pp);
    3831         241 :     r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
    3832             :                                RgX_to_FlxqX(y, Tp, pp),
    3833             :                                RgX_to_FlxqX(S, Tp, pp), Tp, pp));
    3834             :   }
    3835             :   else
    3836         144 :     r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
    3837             :                    RgX_to_FpXQX(S, T, p), T, p);
    3838         385 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3839             : }
    3840             : 
    3841             : static GEN
    3842           0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    3843             : {
    3844           0 :   pari_sp av = avma;
    3845             :   GEN r;
    3846           0 :   GEN T = RgX_to_FpX(pol, p);
    3847           0 :   if (signe(T)==0) pari_err_OP("*",x,x);
    3848           0 :   if (lgefint(p) == 3)
    3849             :   {
    3850           0 :     ulong pp = uel(p, 2);
    3851           0 :     GEN Tp = ZX_to_Flx(T, pp);
    3852           0 :     r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
    3853             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3854             :   }
    3855             :   else
    3856           0 :     r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3857           0 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3858             : }
    3859             : 
    3860             : static GEN
    3861           7 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    3862             : {
    3863           7 :   pari_sp av = avma;
    3864             :   GEN r;
    3865           7 :   GEN T = RgX_to_FpX(pol, p);
    3866           7 :   if (signe(T)==0) pari_err_OP("^",x,gen_m1);
    3867           7 :   if (lgefint(p) == 3)
    3868             :   {
    3869           7 :     ulong pp = uel(p, 2);
    3870           7 :     GEN Tp = ZX_to_Flx(T, pp);
    3871           7 :     r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
    3872             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3873             :   }
    3874             :   else
    3875           0 :     r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3876           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3877             : }
    3878             : 
    3879             : static GEN
    3880      580718 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
    3881             : {
    3882             :   GEN p, pol;
    3883             :   long pa;
    3884      580718 :   long t = RgX_type3(x,y,T, &p,&pol,&pa);
    3885      580717 :   switch(t)
    3886             :   {
    3887      236285 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
    3888       49915 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
    3889         105 :     case t_FFELT:  return FFXQ_mul(x, y, T, pol);
    3890      274288 :     case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
    3891             :     case code(t_POLMOD, t_INTMOD):
    3892         385 :                    return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
    3893       19739 :     default:       return NULL;
    3894             :   }
    3895             : }
    3896             : 
    3897             : GEN
    3898      580718 : RgXQ_mul(GEN x, GEN y, GEN T)
    3899             : {
    3900      580718 :   GEN z = RgXQ_mul_fast(x, y, T);
    3901      580712 :   if (!z) z = RgX_rem(RgX_mul(x, y), T);
    3902      580712 :   return z;
    3903             : }
    3904             : 
    3905             : static GEN
    3906       86599 : RgXQ_sqr_fast(GEN x, GEN T)
    3907             : {
    3908             :   GEN p, pol;
    3909             :   long pa;
    3910       86599 :   long t = RgX_type2(x, T, &p,&pol,&pa);
    3911       86599 :   switch(t)
    3912             :   {
    3913       71381 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
    3914       11598 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
    3915           7 :     case t_FFELT:  return FFXQ_sqr(x, T, pol);
    3916          14 :     case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
    3917             :     case code(t_POLMOD, t_INTMOD):
    3918           0 :                    return RgXQ_sqr_FpXQXQ(x, T, pol, p);
    3919        3599 :     default:       return NULL;
    3920             :   }
    3921             : }
    3922             : 
    3923             : GEN
    3924       86599 : RgXQ_sqr(GEN x, GEN T)
    3925             : {
    3926       86599 :   GEN z = RgXQ_sqr_fast(x, T);
    3927       86599 :   if (!z) z = RgX_rem(RgX_sqr(x), T);
    3928       86599 :   return z;
    3929             : }
    3930             : 
    3931             : static GEN
    3932       31191 : RgXQ_inv_fast(GEN x, GEN y)
    3933             : {
    3934             :   GEN p, pol;
    3935             :   long pa;
    3936       31191 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3937       31191 :   switch(t)
    3938             :   {
    3939       17448 :     case t_INT:    return QXQ_inv(x,y);
    3940         789 :     case t_FRAC:   return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
    3941          14 :     case t_FFELT:  return FFXQ_inv(x, y, pol);
    3942       12054 :     case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
    3943             :     case code(t_POLMOD, t_INTMOD):
    3944           7 :                    return RgXQ_inv_FpXQXQ(x, y, pol, p);
    3945         879 :     default:       return NULL;
    3946             :   }
    3947             : }
    3948             : 
    3949             : GEN
    3950       31191 : RgXQ_inv(GEN x, GEN y)
    3951             : {
    3952       31191 :   GEN z = RgXQ_inv_fast(x, y);
    3953       31177 :   if (!z) z = RgXQ_inv_i(x, y);
    3954       31177 :   return z;
    3955             : }

Generated by: LCOV version 1.13