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.13.0 lcov report (development 25825-edc109b529) Lines: 2151 2389 90.0 %
Date: 2020-09-21 06:08:33 Functions: 198 208 95.2 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13