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.16.2 lcov report (development 29419-8afb0ed749) Lines: 2245 2479 90.6 %
Date: 2024-07-02 09:03:41 Functions: 204 214 95.3 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.16