Code coverage tests

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

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

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

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

Generated by: LCOV version 1.16