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

Generated by: LCOV version 1.13