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 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.10.0 lcov report (development 22303-eb3e11d) Lines: 1781 1979 90.0 %
Date: 2018-04-21 06:16:28 Functions: 166 174 95.4 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. 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       33292 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      30             : {
      31       33292 :   long dP=degpol(P), i, k, m;
      32             :   pari_sp av1, av2;
      33             :   GEN s,y,P_lead;
      34             : 
      35       33292 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      36       33292 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      37       33292 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      38       33292 :   y = cgetg(n+2,t_COL);
      39       33292 :   if (y0)
      40             :   {
      41       11375 :     if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
      42       11375 :     m = lg(y0)-1;
      43       11375 :     for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
      44             :   }
      45             :   else
      46             :   {
      47       21917 :     m = 1;
      48       21917 :     gel(y,1) = stoi(dP);
      49             :   }
      50       33292 :   P += 2; /* strip codewords */
      51             : 
      52       33292 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      53       33292 :   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       86317 :   for (k=m; k<=n; k++)
      59             :   {
      60       53025 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      61      233772 :     for (i=1; i<k && i<=dP; i++)
      62      180747 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      63       53025 :     if (N)
      64             :     {
      65       15960 :       s = Fq_red(s, T, N);
      66       15960 :       if (P_lead) s = Fq_mul(s, P_lead, T, N);
      67             :     }
      68       37065 :     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       37065 :       if (P_lead) s = gdiv(s, P_lead);
      75       53025 :     av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
      76             :   }
      77       33292 :   return y;
      78             : }
      79             : 
      80             : GEN
      81       18676 : polsym(GEN x, long n)
      82             : {
      83       18676 :   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    55653776 : centermodii(GEN x, GEN p, GEN po2)
      89             : {
      90    55653776 :   GEN y = remii(x, p);
      91    55668107 :   switch(signe(y))
      92             :   {
      93     2060103 :     case 0: break;
      94    33442583 :     case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
      95    33401680 :       break;
      96    20276234 :     case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
      97    20255839 :       break;
      98             :   }
      99    55606809 :   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     8375085 : centermod_i(GEN x, GEN p, GEN ps2)
     113             : {
     114             :   long i, lx;
     115             :   pari_sp av;
     116             :   GEN y;
     117             : 
     118     8375085 :   if (!ps2) ps2 = shifti(p,-1);
     119     8375752 :   switch(typ(x))
     120             :   {
     121     3672886 :     case t_INT: return centermodii(x,p,ps2);
     122             : 
     123     3904411 :     case t_POL: lx = lg(x);
     124     3904411 :       y = cgetg(lx,t_POL); y[1] = x[1];
     125    32099766 :       for (i=2; i<lx; i++)
     126             :       {
     127    28194744 :         av = avma;
     128    28194744 :         gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
     129             :       }
     130     3905022 :       return normalizepol_lg(y, lx);
     131             : 
     132      797496 :     case t_COL: lx = lg(x);
     133      797496 :       y = cgetg(lx,t_COL);
     134      797497 :       for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
     135      797495 :       return y;
     136             : 
     137         959 :     case t_MAT: lx = lg(x);
     138         959 :       y = cgetg(lx,t_MAT);
     139         959 :       for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
     140         959 :       return y;
     141             : 
     142           0 :     case t_VECSMALL: lx = lg(x);
     143             :     {
     144           0 :       ulong pp = itou(p), pps2 = itou(ps2);
     145           0 :       y = cgetg(lx,t_VECSMALL);
     146           0 :       for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
     147           0 :       return y;
     148             :     }
     149             :   }
     150           0 :   return x;
     151             : }
     152             : 
     153             : GEN
     154     4805745 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
     155             : 
     156             : static GEN
     157          28 : RgX_Frobenius_deflate(GEN S, ulong p)
     158             : {
     159          28 :   GEN F = RgX_deflate(S, p);
     160          28 :   long i, l = lg(F);
     161         126 :   for (i=2; i<l; i++)
     162             :   {
     163          42 :     GEN Fi = gel(F,i), R;
     164          42 :     if (typ(Fi)==t_POL)
     165             :     {
     166          21 :       if (signe(RgX_deriv(Fi))==0)
     167          14 :         gel(F,i) = RgX_Frobenius_deflate(gel(F, i), p);
     168          14 :       else return NULL;
     169             :     }
     170          21 :     else if (ispower(Fi, utoi(p), &R))
     171          21 :       gel(F,i) = R;
     172           0 :     else return NULL;
     173             :   }
     174          21 :   return F;
     175             : }
     176             : 
     177             : 
     178             : static GEN
     179          84 : RgXY_squff(GEN f)
     180             : {
     181          84 :   long i, q, n = degpol(f);
     182          84 :   ulong p = itos_or_0(characteristic(f));
     183          84 :   GEN u = const_vec(n+1, pol_1(varn(f)));
     184          91 :   for(q = 1;;q *= p)
     185             :   {
     186          91 :     GEN t, v, tv, r = RgX_gcd(f, RgX_deriv(f));
     187          91 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
     188          21 :     t = RgX_div(f, r);
     189          21 :     if (degpol(t) > 0)
     190             :     {
     191             :       long j;
     192          14 :       for(j = 1;;j++)
     193             :       {
     194          14 :         v = RgX_gcd(r, t);
     195          14 :         tv = RgX_div(t, v);
     196          14 :         if (degpol(tv) > 0) gel(u, j*q) = tv;
     197          14 :         if (degpol(v) <= 0) break;
     198           7 :         r = RgX_div(r, v);
     199           7 :         t = v;
     200           7 :       }
     201           7 :       if (degpol(r) == 0) break;
     202             :     }
     203          14 :     if (!p) break;
     204          14 :     r = RgX_Frobenius_deflate(f, p);
     205          14 :     if (!r) { gel(u, q) = f; break; }
     206           7 :     f = r;
     207           7 :   }
     208         413 :   for (i = n; i; i--)
     209         413 :     if (degpol(gel(u,i))) break;
     210          84 :   setlg(u,i+1); return u;
     211             : }
     212             : 
     213             : /* Lmod contains modular factors of *F (NULL codes an empty slot: used factor)
     214             :  * Lfac accumulates irreducible factors as they are found.
     215             :  * p is a product of modular factors in Lmod[1..i-1] (NULL for p = 1), not
     216             :  * a rational factor of *F
     217             :  * Find an irreducible factor of *F divisible by p (by including
     218             :  * exhaustively further factors from Lmod[i..]); return 0 on failure, else 1.
     219             :  * Update Lmod, Lfac and *F */
     220             : static int
     221        4487 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
     222             : {
     223             :   GEN q;
     224        4487 :   if (i == lg(Lmod)) return 0;
     225        2352 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
     226        2226 :   if (!gel(Lmod,i)) return 0;
     227        2219 :   p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
     228        2219 :   q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
     229        2219 :   if (degpol(q))
     230             :   {
     231        2023 :     GEN R, Q = RgX_divrem(*F, q, &R);
     232        2023 :     if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
     233             :   }
     234        2044 :   if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
     235        1904 :   return 0;
     236             : }
     237             : 
     238             : static GEN
     239          91 : RgXY_factor_squarefree(GEN f)
     240             : {
     241          91 :   pari_sp av = avma;
     242          91 :   GEN Lfac, Lmod, F = NULL, BLOC = NULL;
     243          91 :   long vy = gvar2(f), n = RgXY_degreex(f);
     244          91 :   ulong i, c = itou_or_0(residual_characteristic(f));
     245          91 :   long val = RgX_valrem(f, &f);
     246             :   for(;;)
     247             :   {
     248         154 :     for (i = 0; !c || i < c; i++)
     249             :     {
     250         154 :       BLOC = gpowgs(gaddgs(pol_x(vy), i), n+1);
     251         154 :       F = poleval(f, BLOC);
     252         154 :       if (issquarefree(c ? gmul(F,gmodulss(1,c)): F)) break;
     253             :     }
     254          91 :     if (!c || i < c) break;
     255           0 :     n++;
     256           0 :   }
     257          91 :   if (DEBUGLEVEL >= 2)
     258           0 :     err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
     259          91 :   Lmod = gel(factor(F),1);
     260          91 :   if (DEBUGLEVEL >= 2)
     261           0 :     err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
     262          91 :   Lfac = vectrunc_init(lg(Lmod)+1);
     263          91 :   settyp(Lfac, t_COL);
     264          91 :   (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
     265          91 :   if (degpol(f)) vectrunc_append(Lfac, f);
     266          91 :   if (val) vectrunc_append(Lfac, pol_x(varn(f)));
     267          91 :   return gerepilecopy(av, Lfac);
     268             : }
     269             : 
     270             : static GEN
     271          84 : FE_matconcat(GEN F, GEN E, long l)
     272             : {
     273          84 :   setlg(E,l); E = shallowconcat1(E);
     274          84 :   setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
     275             : }
     276             : 
     277             : static GEN
     278          84 : RgXY_factor(GEN f)
     279             : {
     280          84 :   pari_sp av = avma;
     281             :   GEN F, E;
     282          84 :   GEN cnt = content(f);
     283          84 :   GEN V = RgXY_squff(gdiv(f, cnt));
     284          84 :   long i, j, l = lg(V);
     285          84 :   GEN C = factor(cnt);
     286          84 :   F = cgetg(l+1, t_VEC);
     287          84 :   E = cgetg(l+1, t_VEC);
     288          84 :   gel(F,1) = gel(C,1);
     289          84 :   gel(E,1) = gel(C,2);
     290         217 :   for (i=1, j=2; i < l; i++)
     291         133 :     if (degpol(gel(V,i)))
     292             :     {
     293          91 :       GEN Fj = RgXY_factor_squarefree(gel(V,i));
     294          91 :       gel(F, j) = Fj;
     295          91 :       gel(E, j) = const_vec(lg(Fj)-1, stoi(i));
     296          91 :       j++;
     297             :     }
     298          84 :   return gerepilecopy(av, sort_factor_pol(FE_matconcat(F,E,j), cmp_universal));
     299             : }
     300             : 
     301             : /***********************************************************************/
     302             : /**                                                                   **/
     303             : /**                          FACTORIZATION                            **/
     304             : /**                                                                   **/
     305             : /***********************************************************************/
     306             : #define assign_or_fail(x,y) { GEN __x = x;\
     307             :   if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
     308             : }
     309             : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
     310             : 
     311             : static const long tsh = 6;
     312             : static long
     313      222004 : code(long t1, long t2) { return (t1 << tsh) | t2; }
     314             : void
     315      401982 : RgX_type_decode(long x, long *t1, long *t2)
     316             : {
     317      401982 :   *t1 = x >> tsh;
     318      401982 :   *t2 = (x & ((1L<<tsh)-1));
     319      401982 : }
     320             : int
     321     2455071 : RgX_type_is_composite(long t) { return t >= tsh; }
     322             : 
     323             : static int
     324   544216103 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     325             : {
     326             :   long j;
     327   544216103 :   switch(typ(c))
     328             :   {
     329             :     case t_INT:
     330   456477584 :       break;
     331             :     case t_FRAC:
     332    13350597 :       t[1]=1; break;
     333             :       break;
     334             :     case t_REAL:
     335    26373772 :       update_prec(precision(c), pa);
     336    26373772 :       t[2]=1; break;
     337             :     case t_INTMOD:
     338    20197277 :       assign_or_fail(gel(c,1),p);
     339    20197277 :       t[3]=1; break;
     340             :     case t_FFELT:
     341     1358762 :       if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
     342     1358762 :       assign_or_fail(FF_p_i(c),p);
     343     1358762 :       t[5]=1; break;
     344             :     case t_COMPLEX:
     345    23673673 :       for (j=1; j<=2; j++)
     346             :       {
     347    15782451 :         GEN d = gel(c,j);
     348    15782451 :         switch(typ(d))
     349             :         {
     350             :           case t_INT: case t_FRAC:
     351      644267 :             t[4]=1; break;
     352             :           case t_REAL:
     353    15138163 :             update_prec(precision(d), pa);
     354    15138163 :             t[6]=1; break;
     355             :           case t_INTMOD:
     356          14 :             assign_or_fail(gel(d,1),p);
     357          14 :             if (!signe(*p) || mod4(*p) != 3) return 0;
     358           7 :             if (!*t2) *t2 = t_COMPLEX;
     359           7 :             t[3]=1; break;
     360             :           case t_PADIC:
     361           7 :             update_prec(precp(d)+valp(d), pa);
     362           7 :             assign_or_fail(gel(d,2),p);
     363           7 :             if (!*t2) *t2 = t_COMPLEX;
     364           7 :             t[7]=1; break;
     365           0 :           default: return 0;
     366             :         }
     367             :       }
     368     7891222 :       if (!t[6]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
     369     7891222 :       break;
     370             :     case t_PADIC:
     371      208411 :       update_prec(precp(c)+valp(c), pa);
     372      208411 :       assign_or_fail(gel(c,2),p);
     373      208411 :       t[7]=1; break;
     374             :     case t_QUAD:
     375         210 :       assign_or_fail(gel(c,1),pol);
     376         630 :       for (j=2; j<=3; j++)
     377             :       {
     378         420 :         GEN d = gel(c,j);
     379         420 :         switch(typ(d))
     380             :         {
     381             :           case t_INT: case t_FRAC:
     382         385 :             t[8]=1; break;
     383             :           case t_INTMOD:
     384          28 :             assign_or_fail(gel(d,1),p);
     385          28 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     386          28 :             t[3]=1; break;
     387             :           case t_PADIC:
     388           7 :             update_prec(precp(d)+valp(d), pa);
     389           7 :             assign_or_fail(gel(d,2),p);
     390           7 :             if (*t2 != t_POLMOD) *t2 = t_QUAD;
     391           7 :             t[7]=1; break;
     392           0 :           default: return 0;
     393             :         }
     394             :       }
     395         210 :       break;
     396             :     case t_POLMOD:
     397     3132752 :       assign_or_fail(gel(c,1),pol);
     398     3132752 :       if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1)))
     399         182 :         return 0;
     400    18790646 :       for (j=1; j<=2; j++)
     401             :       {
     402             :         GEN pbis, polbis;
     403             :         long pabis;
     404     6264545 :         *t2 = t_POLMOD;
     405     6264545 :         switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
     406             :         {
     407     3284912 :           case t_INT:  break;
     408      486891 :           case t_FRAC: t[1]=1; break;
     409     2490943 :           case t_INTMOD: t[3]=1; break;
     410           7 :           case t_PADIC: t[7]=1; update_prec(pabis,pa); break;
     411        3584 :           default: return 0;
     412             :         }
     413     6262753 :         if (pbis) assign_or_fail(pbis,p);
     414     6262753 :         if (polbis) assign_or_fail(polbis,pol);
     415             :       }
     416     3130778 :       break;
     417    14007975 :     case t_POL: t[10] = 1;
     418             :     {
     419             :       GEN pbis, polbis;
     420             :       long tbis, pabis;
     421    14007975 :       tbis = RgX_type(c,&pbis,&polbis,&pabis);
     422    14007975 :       if (tbis && tbis!=t_POL && !polbis)
     423             :       {
     424    13527753 :         if (*var == NO_VARIABLE) { *var = varn(c); break; }
     425     8324622 :         if (*var == varn(c)) break;
     426             :       }
     427      670946 :       *var = MAXVARN+1; break; /* ensure varn() == *var fails in choosetype */
     428             :     }
     429     1217534 :     default: return 0;
     430             :   }
     431   542996588 :   return 1;
     432             : }
     433             : /* t[0] unused. Other values, if set, indicate a coefficient of type
     434             :  * t[1] : t_FRAC
     435             :  * t[2] : t_REAL
     436             :  * t[3] : t_INTMOD
     437             :  * t[4] : t_COMPLEX of rationals (t_INT/t_FRAC)
     438             :  * t[5] : t_FFELT
     439             :  * t[6] : t_COMPLEX of t_REAL
     440             :  * t[7] : t_PADIC
     441             :  * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
     442             :  * t[9]:  Unused
     443             :  * t[10]: t_POL (recursive factorisation) */
     444             : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
     445             :  * given by t) */
     446             : static long
     447    75377088 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
     448             : {
     449    75377088 :   if (t[10] && (!*pol || var!=varn(*pol))) return t_POL;
     450    70046761 :   if (t[5]) /* ffelt */
     451             :   {
     452      222832 :     if (t2 ||t[2]||t[4]||t[6]||t[8]||t[9]) return 0;
     453      222832 :     *pol=ff; return t_FFELT;
     454             :   }
     455    69823929 :   if (t[6]) /* inexact, complex */
     456             :   {
     457     1145525 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
     458     1145525 :     return t_COMPLEX;
     459             :   }
     460    68678404 :   if (t[2]) /* inexact, real */
     461             :   {
     462     4040397 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
     463     4040397 :     return t[4]?t_COMPLEX:t_REAL;
     464             :   }
     465    64638007 :   if (t2) /* polmod/quad/complex of intmod/padic */
     466             :   {
     467      216558 :     if (t[3]) return code(t2,t_INTMOD);
     468      189888 :     if (t[7]) return code(t2,t_PADIC);
     469      189867 :     if (t[1]) return code(t2,t_FRAC);
     470      123393 :     return code(t2,t_INT);
     471             :   }
     472    64421449 :   if (t[10]) return t_POL;
     473    64421449 :   if (t[8]) return code(t_QUAD,t_INT);
     474    64421337 :   if (t[4]) return code(t_COMPLEX,t_INT);
     475    64416003 :   if (t[3]) return t_INTMOD;
     476    60821979 :   if (t[7]) return t_PADIC;
     477    60803506 :   if (t[1]) return t_FRAC;
     478    57805988 :   return t_INT;
     479             : }
     480             : 
     481             : static long
     482   121412476 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     483             : {
     484   121412476 :   long i, lx = lg(x);
     485   568716321 :   for (i=2; i<lx; i++)
     486   448526035 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     487   120190286 :   return 1;
     488             : }
     489             : 
     490             : static long
     491    15956883 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     492             : {
     493    15956883 :   long i, l = lg(x);
     494   111420426 :   for (i = 1; i<l; i++)
     495    95465216 :     if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     496    15955210 :   return 1;
     497             : }
     498             : 
     499             : static long
     500     5355522 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
     501             : {
     502     5355522 :   long i, l = lg(x);
     503    21048127 :   for (i = 1; i < l; i++)
     504    15694243 :     if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
     505     5353884 :   return 1;
     506             : }
     507             : 
     508             : long
     509     8719616 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
     510             : {
     511     8719616 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0};
     512     8719616 :   long tx = typ(x), t2 = 0, var = NO_VARIABLE;
     513     8719616 :   GEN ff = NULL;
     514     8719616 :   *p = *pol = NULL; *pa = LONG_MAX;
     515     8719616 :   if (is_scalar_t(tx))
     516             :   {
     517      226578 :     if (tx == t_POLMOD) return 0;
     518      226193 :     if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     519             :   }
     520     8493038 :   else if (tx == t_MAT)
     521           7 :   { if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0; }
     522             :   else /* t_POL, t_SER */
     523     8493031 :     if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     524     8719231 :   return choosetype(t,t2,ff,pol,var);
     525             : }
     526             : 
     527             : long
     528    16667361 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     529             : {
     530    16667361 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     531    16667361 :   long t2 = 0, var = NO_VARIABLE;
     532    16667361 :   GEN ff = NULL;
     533    16667361 :   *p = *pol = NULL; *pa = LONG_MAX;
     534    16667361 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     535    16648535 :   return choosetype(t,t2,ff,pol,var);
     536             : }
     537             : 
     538             : long
     539    47640502 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     540             : {
     541    47640502 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     542    47640502 :   long t2 = 0, var = NO_VARIABLE;
     543    47640502 :   GEN ff = NULL;
     544    47640502 :   *p = *pol = NULL; *pa = LONG_MAX;
     545    94093533 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     546    47651509 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     547    46442703 :   return choosetype(t,t2,ff,pol,var);
     548             : }
     549             : 
     550             : long
     551      718390 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
     552             : {
     553      718390 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     554      718390 :   long t2 = 0, var = NO_VARIABLE;
     555      718390 :   GEN ff = NULL;
     556      718390 :   *p = *pol = NULL; *pa = LONG_MAX;
     557     1436784 :   if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     558     1436787 :       !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
     559      718394 :       !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
     560      718397 :   return choosetype(t,t2,ff,pol,var);
     561             : }
     562             : 
     563             : long
     564       83538 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
     565             : {
     566       83538 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     567       83538 :   long t2 = 0, var = NO_VARIABLE;
     568       83538 :   GEN ff = NULL;
     569       83538 :   *p = *pol = NULL; *pa = LONG_MAX;
     570       83538 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
     571       82838 :   return choosetype(t,t2,ff,pol,var);
     572             : }
     573             : 
     574             : long
     575      263032 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     576             : {
     577      263032 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     578      263032 :   long t2 = 0, var = NO_VARIABLE;
     579      263032 :   GEN ff = NULL;
     580      263032 :   *p = *pol = NULL; *pa = LONG_MAX;
     581      525672 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     582      263067 :       !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     583      262605 :   return choosetype(t,t2,ff,pol,var);
     584             : }
     585             : 
     586             : long
     587     2504693 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
     588             : {
     589     2504693 :   long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
     590     2504693 :   long t2 = 0, var = NO_VARIABLE;
     591     2504693 :   GEN ff = NULL;
     592     2504693 :   *p = *pol = NULL; *pa = LONG_MAX;
     593     5008945 :   if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
     594     2504798 :       !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
     595     2504147 :   return choosetype(t,t2,ff,pol,var);
     596             : }
     597             : 
     598             : GEN
     599           0 : factor0(GEN x,long flag) { return (flag<0)? factor(x): boundfact(x,flag); }
     600             : 
     601             : /* only present for interface with GP */
     602             : GEN
     603       37164 : gp_factor0(GEN x, GEN flag)
     604             : {
     605             :   ulong B;
     606       37164 :   if (!flag) return factor(x);
     607          35 :   if (typ(flag) != t_INT || signe(flag) < 0) pari_err_FLAG("factor");
     608          35 :   switch(lgefint(flag))
     609             :   {
     610           7 :     case 2: B = 0; break;
     611          28 :     case 3: B = flag[2]; break;
     612           0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     613             :              return NULL; /*LCOV_EXCL_LINE*/
     614             :   }
     615          35 :   return boundfact(x, B);
     616             : }
     617             : 
     618             : GEN
     619       83403 : deg1_from_roots(GEN L, long v)
     620             : {
     621       83403 :   long i, l = lg(L);
     622       83403 :   GEN z = cgetg(l,t_COL);
     623      190457 :   for (i=1; i<l; i++)
     624      107054 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     625       83403 :   return z;
     626             : }
     627             : GEN
     628         980 : roots_from_deg1(GEN x)
     629             : {
     630         980 :   long i,l = lg(x);
     631         980 :   GEN r = cgetg(l,t_VEC);
     632         980 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     633         980 :   return r;
     634             : }
     635             : 
     636             : static GEN
     637          35 : gauss_factor_p(GEN p)
     638             : {
     639          35 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     640          35 :   return mkcomplex(a, b);
     641             : }
     642             : 
     643             : static GEN
     644          42 : gauss_primpart(GEN x, GEN *c)
     645             : {
     646          42 :   GEN a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
     647          42 :   *c = n; if (n == gen_1) return x;
     648          42 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     649             : }
     650             : 
     651             : static GEN
     652          70 : gauss_primpart_try(GEN x, GEN c)
     653             : {
     654             :   GEN r, y;
     655          70 :   if (typ(x) == t_INT)
     656             :   {
     657          42 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     658             :   }
     659             :   else
     660             :   {
     661          28 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     662          28 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     663          14 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     664             :   }
     665          56 :   return y;
     666             : }
     667             : 
     668             : static int
     669          77 : gauss_cmp(GEN x, GEN y)
     670             : {
     671             :   int v;
     672          77 :   if (typ(x) != t_COMPLEX)
     673           0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     674          77 :   if (typ(y) != t_COMPLEX) return 1;
     675          49 :   v = cmpii(gel(x,2), gel(y,2));
     676          49 :   if (v) return v;
     677          21 :   return gcmp(gel(x,1), gel(y,1));
     678             : }
     679             : 
     680             : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
     681             : static GEN
     682         448 : gauss_normal(GEN x)
     683             : {
     684         448 :   if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
     685         378 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     686         378 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     687         378 :   return x;
     688             : }
     689             : 
     690             : static GEN
     691          42 : gauss_factor(GEN x)
     692             : {
     693          42 :   pari_sp av = avma;
     694          42 :   GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
     695          42 :   long t1 = typ(a);
     696          42 :   long t2 = typ(b), i, j, l, exp = 0;
     697          42 :   if (t1 == t_FRAC) d = gel(a,2);
     698          42 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     699          42 :   if (d == gen_1) y = x;
     700             :   else
     701             :   {
     702          14 :     y = gmul(x, d);
     703          14 :     a = gel(y,1); t1 = typ(a);
     704          14 :     b = gel(y,2); t2 = typ(b);
     705             :   }
     706          42 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     707          42 :   y = gauss_primpart(y, &n);
     708          42 :   fa = factor(cxnorm(y));
     709          42 :   P = gel(fa,1);
     710          42 :   E = gel(fa,2); l = lg(P);
     711          42 :   P2 = cgetg(l, t_COL);
     712          42 :   E2 = cgetg(l, t_COL);
     713          98 :   for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
     714             :   { /* either p = 2 (ramified) or those factors split in Q(i) */
     715          56 :     GEN p = gel(P,i), w, w2, t, we, pe;
     716          56 :     long v, e = itos(gel(E,i));
     717          56 :     int is2 = absequaliu(p, 2);
     718          56 :     w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     719          56 :     w2 = gauss_normal( gconj(w) );
     720             :     /* w * w2 * I^3 = p, w2 = gconj(w) * I */
     721          56 :     pe = powiu(p, e);
     722          56 :     we = gpowgs(w, e);
     723          56 :     t = gauss_primpart_try( gmul(y, gconj(we)), pe );
     724          56 :     if (t) y = t; /* y /= w^e */
     725             :     else {
     726             :       /* y /= conj(w)^e, should be y /= w2^e */
     727          14 :       y = gauss_primpart_try( gmul(y, we), pe );
     728          14 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     729             :     }
     730          56 :     gel(P,i) = w;
     731          56 :     v = Z_pvalrem(n, p, &n);
     732          56 :     if (v) {
     733           7 :       exp -= v; /* += 3*v mod 4 */
     734           7 :       if (is2) v <<= 1; /* 2 = w^2 I^3 */
     735             :       else {
     736           0 :         gel(P2,j) = w2;
     737           0 :         gel(E2,j) = utoipos(v); j++;
     738             :       }
     739           7 :       gel(E,i) = stoi(e + v);
     740             :     }
     741          56 :     v = Z_pvalrem(d, p, &d);
     742          56 :     if (v) {
     743           7 :       exp += v; /* -= 3*v mod 4 */
     744           7 :       if (is2) v <<= 1; /* 2 is ramified */
     745             :       else {
     746           7 :         gel(P2,j) = w2;
     747           7 :         gel(E2,j) = utoineg(v); j++;
     748             :       }
     749           7 :       gel(E,i) = stoi(e - v);
     750             :     }
     751          56 :     exp &= 3;
     752             :   }
     753          42 :   if (j > 1) {
     754           7 :     long k = 1;
     755           7 :     GEN P1 = cgetg(l, t_COL);
     756           7 :     GEN E1 = cgetg(l, t_COL);
     757             :     /* remove factors with exponent 0 */
     758          14 :     for (i = 1; i < l; i++)
     759           7 :       if (signe(gel(E,i)))
     760             :       {
     761           0 :         gel(P1,k) = gel(P,i);
     762           0 :         gel(E1,k) = gel(E,i);
     763           0 :         k++;
     764             :       }
     765           7 :     setlg(P1, k); setlg(E1, k);
     766           7 :     setlg(P2, j); setlg(E2, j);
     767           7 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     768             :   }
     769          42 :   if (!is_pm1(n) || !is_pm1(d))
     770             :   {
     771          21 :     GEN Fa = factor(gdiv(n, d));
     772          21 :     P = gel(Fa,1); l = lg(P);
     773          21 :     E = gel(Fa,2);
     774          49 :     for (i = 1; i < l; i++)
     775             :     {
     776          28 :       GEN w, p = gel(P,i);
     777             :       long e;
     778             :       int is2;
     779          28 :       switch(mod4(p))
     780             :       {
     781          14 :         case 3: continue;
     782           7 :         case 2: is2 = 1; break;
     783           7 :         default:is2 = 0; break;
     784             :       }
     785          14 :       e = itos(gel(E,i));
     786          14 :       w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     787          14 :       gel(P,i) = w;
     788          14 :       if (is2)
     789           7 :         gel(E,i) = stoi(2*e);
     790             :       else
     791             :       {
     792           7 :         P = shallowconcat(P, gauss_normal( gconj(w) ));
     793           7 :         E = shallowconcat(E, gel(E,i));
     794             :       }
     795          14 :       exp -= e; /* += 3*e mod 4 */
     796          14 :       exp &= 3;
     797             :     }
     798          21 :     gel(Fa,1) = P;
     799          21 :     gel(Fa,2) = E;
     800          21 :     fa = famat_mul_shallow(fa, Fa);
     801             :   }
     802          42 :   fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
     803             : 
     804          42 :   y = gmul(y, powIs(exp));
     805          42 :   if (!gequal1(y)) {
     806          35 :     gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
     807          35 :     gel(fa,2) = shallowconcat(gen_1,    gel(fa,2));
     808             :   }
     809          42 :   return gerepilecopy(av, fa);
     810             : }
     811             : 
     812             : GEN
     813           0 : Q_factor_limit(GEN x, ulong lim)
     814             : {
     815           0 :   pari_sp av = avma;
     816             :   GEN a, b;
     817           0 :   if (typ(x) == t_INT) return Z_factor_limit(x, lim);
     818           0 :   a = Z_factor_limit(gel(x,1), lim);
     819           0 :   b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
     820           0 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     821             : }
     822             : GEN
     823        2527 : Q_factor(GEN x)
     824             : {
     825        2527 :   pari_sp av = avma;
     826             :   GEN a, b;
     827        2527 :   if (typ(x) == t_INT) return Z_factor(x);
     828         182 :   a = Z_factor(gel(x,1));
     829         182 :   b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
     830         182 :   return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
     831             : }
     832             : /* replace t_QUAD/t_COMPLEX coeffs by t_POLMOD in T */
     833             : static GEN
     834          56 : RgX_fix_quad(GEN x, GEN T)
     835             : {
     836          56 :   long i, l, v = varn(T);
     837          56 :   GEN y = cgetg_copy(x,&l);
     838         294 :   for (i = 2; i < l; i++)
     839             :   {
     840         238 :     GEN c = gel(x,i);
     841         238 :     switch(typ(c))
     842             :     {
     843          28 :       case t_QUAD: c++;/* fall through */
     844          70 :       case t_COMPLEX: c = deg1pol_shallow(gel(c,2),gel(c,1),v);
     845             :     }
     846         238 :     gel(y,i) = c;
     847             :   }
     848          56 :   y[1] = x[1]; return y;
     849             : }
     850             : 
     851             : GEN
     852       40384 : factor(GEN x)
     853             : {
     854       40384 :   long tx=typ(x), lx, i, pa, v, r1;
     855             :   pari_sp av, tetpil;
     856             :   GEN  y, p, p1, p2, pol;
     857             : 
     858       40384 :   if (gequal0(x))
     859          63 :     switch(tx)
     860             :     {
     861             :       case t_INT:
     862             :       case t_COMPLEX:
     863             :       case t_POL:
     864          63 :       case t_RFRAC: return prime_fact(x);
     865           0 :       default: pari_err_TYPE("factor",x);
     866             :     }
     867       40321 :   av = avma;
     868       40321 :   switch(tx)
     869             :   {
     870       38228 :     case t_INT: return Z_factor(x);
     871         168 :     case t_FRAC: return Q_factor(x);
     872             : 
     873             :     case t_POL:
     874        1876 :       tx=RgX_type(x,&p,&pol,&pa);
     875        1876 :       switch(tx)
     876             :       {
     877           7 :         case 0: pari_err_IMPL("factor for general polynomials");
     878          84 :         case t_POL: return RgXY_factor(x);
     879        1400 :         case t_INT: return ZX_factor(x);
     880           7 :         case t_FRAC: return QX_factor(x);
     881          70 :         case t_INTMOD: return factmod(x,p);
     882          28 :         case t_PADIC: return factorpadic(x,p,pa);
     883          98 :         case t_FFELT: return FFX_factor(x,pol);
     884             : 
     885          14 :         case t_COMPLEX: y = cgetg(3,t_MAT);
     886          14 :           av = avma; p1 = deg1_from_roots(roots(x,pa), varn(x));
     887          14 :           gel(y,1) = p1 = gerepileupto(av, p1);
     888          14 :           gel(y,2) = const_col(lg(p1)-1, gen_1); return y;
     889             : 
     890          21 :         case t_REAL: y=cgetg(3,t_MAT); v=varn(x);
     891          21 :           av=avma; p1=cleanroots(x,pa); tetpil=avma;
     892          21 :           lx = lg(p1);
     893          49 :           for (r1 = 1; r1 < lx; r1++)
     894          35 :             if (typ(gel(p1,r1)) == t_COMPLEX) break;
     895          21 :           lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
     896          49 :           for (i = 1; i < r1; i++)
     897          28 :             gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
     898          28 :           for (   ; i < lx; i++)
     899             :           {
     900           7 :             GEN a = gel(p1,2*i-r1);
     901           7 :             p = cgetg(5, t_POL); gel(p2,i) = p;
     902           7 :             p[1] = x[1];
     903           7 :             gel(p,2) = gnorm(a);
     904           7 :             gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
     905           7 :             gel(p,4) = gen_1;
     906             :           }
     907          21 :           gel(y,1) = gerepile(av,tetpil,p2);
     908          21 :           gel(y,2) = const_col(lx-1, gen_1); return y;
     909             : 
     910             :         default:
     911             :         {
     912         147 :           GEN w = NULL, T = pol;
     913             :           long t1, t2;
     914         147 :           RgX_type_decode(tx, &t1, &t2);
     915         147 :           if (t1 == t_COMPLEX) w = gen_I();
     916         119 :           else if (t1 == t_QUAD) w = mkquad(pol,gen_0,gen_1);
     917         147 :           if (w)
     918             :           { /* substitute I or w by t_POLMOD */
     919          56 :             T = leafcopy(pol); setvarn(T, fetch_var());
     920          56 :             x = RgX_fix_quad(x, T);
     921             :           }
     922         147 :           switch (t2)
     923             :           {
     924          91 :             case t_INT: case t_FRAC: p1 = nffactor(T,x); break;
     925             :             case t_INTMOD:
     926          35 :               T = RgX_to_FpX(T,p);
     927          35 :               if (FpX_is_irred(T,p)) { p1 = factmod(x,mkvec2(p,T)); break; }
     928             :             /*fall through*/
     929             :             default:
     930          35 :               if (w) (void)delete_var();
     931          35 :               pari_err_IMPL("factor for general polynomial");
     932             :               return NULL; /* LCOV_EXCL_LINE */
     933             :           }
     934         112 :           if (t1 == t_POLMOD) return gerepileupto(av, p1);
     935             :           /* substitute back I or w */
     936          35 :           gel(p1,1) = gsubst(liftpol_shallow(gel(p1,1)), varn(T), w);
     937          35 :           (void)delete_var(); return gerepilecopy(av, p1);
     938             :         }
     939             :       }
     940             :     case t_RFRAC: {
     941           7 :       GEN a = gel(x,1), b = gel(x,2);
     942           7 :       y = famat_inv_shallow(factor(b));
     943           7 :       if (typ(a)==t_POL && varn(a)==varn(b)) y = famat_mul_shallow(factor(a),y);
     944           7 :       return gerepilecopy(av, y);
     945             :     }
     946             :     case t_COMPLEX:
     947          42 :       y = gauss_factor(x); if (y) return y;
     948             :       /* fall through */
     949             :   }
     950           0 :   pari_err_TYPE("factor",x);
     951             :   return NULL; /* LCOV_EXCL_LINE */
     952             : }
     953             : 
     954             : /*******************************************************************/
     955             : /*                                                                 */
     956             : /*                     ROOTS --> MONIC POLYNOMIAL                  */
     957             : /*                                                                 */
     958             : /*******************************************************************/
     959             : static GEN
     960      634070 : normalized_mul(void *E, GEN x, GEN y)
     961             : {
     962      634070 :   long a = gel(x,1)[1], b = gel(y,1)[1];
     963             :   (void) E;
     964     1268140 :   return mkvec2(mkvecsmall(a + b),
     965     1268140 :                 RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
     966             : }
     967             : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
     968             : static GEN
     969      320499 : normalized_to_RgX(GEN L)
     970             : {
     971      320499 :   long i, a = gel(L,1)[1];
     972      320499 :   GEN A = gel(L,2);
     973      320499 :   GEN z = cgetg(a + 3, t_POL);
     974      320499 :   z[1] = evalsigne(1) | evalvarn(varn(A));
     975      320499 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
     976      320499 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
     977      320499 :   gel(z,i) = gen_1; return z;
     978             : }
     979             : 
     980             : /* compute prod (x - a[i]) */
     981             : GEN
     982      276716 : roots_to_pol(GEN a, long v)
     983             : {
     984      276716 :   pari_sp av = avma;
     985      276716 :   long i, k, lx = lg(a);
     986             :   GEN L;
     987      276716 :   if (lx == 1) return pol_1(v);
     988      276681 :   L = cgetg(lx, t_VEC);
     989      589928 :   for (k=1,i=1; i<lx-1; i+=2)
     990             :   {
     991      313247 :     GEN s = gel(a,i), t = gel(a,i+1);
     992      313247 :     GEN x0 = gmul(s,t);
     993      313247 :     GEN x1 = gneg(gadd(s,t));
     994      313247 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     995             :   }
     996      536299 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
     997      259618 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
     998      276681 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
     999      276681 :   return gerepileupto(av, normalized_to_RgX(L));
    1000             : }
    1001             : 
    1002             : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
    1003             : GEN
    1004       43818 : roots_to_pol_r1(GEN a, long v, long r1)
    1005             : {
    1006       43818 :   pari_sp av = avma;
    1007       43818 :   long i, k, lx = lg(a);
    1008             :   GEN L;
    1009       43818 :   if (lx == 1) return pol_1(v);
    1010       43818 :   L = cgetg(lx, t_VEC);
    1011      299873 :   for (k=1,i=1; i<r1; i+=2)
    1012             :   {
    1013      256055 :     GEN s = gel(a,i), t = gel(a,i+1);
    1014      256055 :     GEN x0 = gmul(s,t);
    1015      256055 :     GEN x1 = gneg(gadd(s,t));
    1016      256055 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1017             :   }
    1018       64272 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
    1019       20454 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
    1020      149013 :   for (i=r1+1; i<lx; i++)
    1021             :   {
    1022      105195 :     GEN s = gel(a,i);
    1023      105195 :     GEN x0 = gnorm(s);
    1024      105195 :     GEN x1 = gneg(gtrace(s));
    1025      105195 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
    1026             :   }
    1027       43818 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
    1028       43818 :   return gerepileupto(av, normalized_to_RgX(L));
    1029             : }
    1030             : 
    1031             : /*******************************************************************/
    1032             : /*                                                                 */
    1033             : /*                          FACTORBACK                             */
    1034             : /*                                                                 */
    1035             : /*******************************************************************/
    1036             : static GEN
    1037           7 : idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
    1038             : static GEN
    1039         609 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
    1040             : static GEN
    1041        2198 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
    1042             : static GEN
    1043        3423 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
    1044             : static GEN
    1045        4164 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
    1046             : static GEN
    1047       16940 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
    1048             : static GEN
    1049    54162803 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
    1050             : static GEN
    1051    78776658 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
    1052             : static GEN
    1053    16094057 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
    1054             : static GEN
    1055      119504 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
    1056             : 
    1057             : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
    1058             : GEN
    1059    28043739 : gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
    1060             :                GEN (*_pow)(void*,GEN,GEN), void *data)
    1061             : {
    1062    28043739 :   pari_sp av = avma;
    1063             :   long k, l, lx;
    1064             :   GEN p,x;
    1065             : 
    1066    28043739 :   if (e) /* supplied vector of exponents */
    1067      852894 :     p = L;
    1068             :   else
    1069             :   {
    1070    27190845 :     switch(typ(L)) {
    1071             :       case t_VEC:
    1072             :       case t_COL: /* product of the L[i] */
    1073     3274378 :         return gerepileupto(av, gen_product(L, data, _mul));
    1074             :       case t_MAT: /* genuine factorization */
    1075    23916467 :         l = lg(L);
    1076    23916467 :         if (l == 1) return gen_1;
    1077    23916460 :         if (l == 3) break;
    1078             :         /*fall through*/
    1079             :       default:
    1080           7 :         pari_err_TYPE("factorback [not a factorization]", L);
    1081             :     }
    1082    23916453 :     p = gel(L,1);
    1083    23916453 :     e = gel(L,2);
    1084             :   }
    1085             :   /* p = elts, e = expo */
    1086    24769347 :   lx = lg(p);
    1087             :   /* check whether e is an integral vector of correct length */
    1088    24769347 :   switch(typ(e))
    1089             :   {
    1090             :     case t_VECSMALL:
    1091         294 :       if (lx != lg(e))
    1092           0 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1093         294 :       if (lx == 1) return gen_1;
    1094         266 :       x = cgetg(lx,t_VEC);
    1095         882 :       for (l=1,k=1; k<lx; k++)
    1096         616 :         if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
    1097         266 :       break;
    1098             :     case t_VEC: case t_COL:
    1099    24769053 :       if (lx != lg(e) || !RgV_is_ZV(e))
    1100           7 :         pari_err_TYPE("factorback [not an exponent vector]", e);
    1101    24769046 :       if (lx == 1) return gen_1;
    1102    24749570 :       x = cgetg(lx,t_VEC);
    1103   103836160 :       for (l=1,k=1; k<lx; k++)
    1104    79086590 :         if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
    1105    24749570 :       break;
    1106             :     default:
    1107           0 :       pari_err_TYPE("factorback [not an exponent vector]", e);
    1108           0 :       return NULL;
    1109             :   }
    1110    24749836 :   x[0] = evaltyp(t_VEC) | _evallg(l);
    1111    24749836 :   return gerepileupto(av, gen_product(x, data, _mul));
    1112             : }
    1113             : 
    1114             : GEN
    1115        1911 : idealfactorback(GEN nf, GEN L, GEN e, int red)
    1116             : {
    1117        1911 :   nf = checknf(nf);
    1118        1911 :   if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
    1119        1309 :   else     return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
    1120             : }
    1121             : 
    1122             : GEN
    1123       13286 : nffactorback(GEN nf, GEN L, GEN e)
    1124       13286 : { return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }
    1125             : 
    1126             : GEN
    1127     3369663 : FpV_factorback(GEN L, GEN e, GEN p)
    1128     3369663 : { return gen_factorback(L, e, &Fpmul, &Fppow, (void*)p); }
    1129             : 
    1130             : GEN
    1131    24643682 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
    1132             : GEN
    1133     1395129 : factorback(GEN fa) { return factorback2(fa, NULL); }
    1134             : 
    1135             : GEN
    1136          14 : vecprod(GEN v)
    1137             : {
    1138          14 :   pari_sp av = avma;
    1139          14 :   if (!is_vec_t(typ(v)))
    1140           0 :     pari_err_TYPE("vecprod", v);
    1141          14 :   if (lg(v) == 1) return gen_1;
    1142           7 :   return gerepilecopy(av, gen_product(v, NULL, mul));
    1143             : }
    1144             : 
    1145             : static int
    1146          63 : RgX_is_irred_i(GEN x)
    1147             : {
    1148             :   GEN y, p, pol;
    1149          63 :   long l = lg(x), pa;
    1150             : 
    1151          63 :   if (!signe(x) || l <= 3) return 0;
    1152          63 :   switch(RgX_type(x,&p,&pol,&pa))
    1153             :   {
    1154          14 :     case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
    1155           0 :     case t_COMPLEX: return l == 4;
    1156             :     case t_REAL:
    1157           0 :       if (l == 4) return 1;
    1158           0 :       if (l > 5) return 0;
    1159           0 :       return gsigne(RgX_disc(x)) > 0;
    1160             :   }
    1161          49 :   y = factor(x);
    1162          49 :   return (lg(gcoeff(y,1,1))==l);
    1163             : }
    1164             : static int
    1165          63 : RgX_is_irred(GEN x)
    1166             : {
    1167          63 :   pari_sp av = avma;
    1168          63 :   int r = RgX_is_irred_i(x);
    1169          63 :   avma = av; return r;
    1170             : }
    1171             : long
    1172          63 : isirreducible(GEN x)
    1173             : {
    1174          63 :   switch(typ(x))
    1175             :   {
    1176           0 :     case t_INT: case t_REAL: case t_FRAC: return 0;
    1177          63 :     case t_POL: return RgX_is_irred(x);
    1178             :   }
    1179           0 :   pari_err_TYPE("isirreducible",x);
    1180           0 :   return 0;
    1181             : }
    1182             : 
    1183             : /*******************************************************************/
    1184             : /*                                                                 */
    1185             : /*                         GENERIC GCD                             */
    1186             : /*                                                                 */
    1187             : /*******************************************************************/
    1188             : /* x is a COMPLEX or a QUAD */
    1189             : static GEN
    1190        2499 : triv_cont_gcd(GEN x, GEN y)
    1191             : {
    1192        2499 :   pari_sp av = avma;
    1193             :   GEN c;
    1194        2499 :   if (typ(x)==t_COMPLEX)
    1195             :   {
    1196        2177 :     GEN a = gel(x,1), b = gel(x,2);
    1197        2177 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
    1198          21 :     c = ggcd(a,b);
    1199             :   }
    1200             :   else
    1201         322 :     c = ggcd(gel(x,2),gel(x,3));
    1202         343 :   return gerepileupto(av, ggcd(c,y));
    1203             : }
    1204             : 
    1205             : /* y is a PADIC, x a rational number or an INTMOD */
    1206             : static GEN
    1207        2684 : padic_gcd(GEN x, GEN y)
    1208             : {
    1209        2684 :   GEN p = gel(y,2);
    1210        2684 :   long v = gvaluation(x,p), w = valp(y);
    1211        2684 :   if (w < v) v = w;
    1212        2684 :   return powis(p, v);
    1213             : }
    1214             : 
    1215             : /* x,y in Z[i], at least one of which is t_COMPLEX */
    1216             : static GEN
    1217         385 : gauss_gcd(GEN x, GEN y)
    1218             : {
    1219         385 :   pari_sp av = avma;
    1220             :   GEN dx, dy;
    1221         385 :   dx = denom(x); x = gmul(x, dx);
    1222         385 :   dy = denom(y); y = gmul(y, dy);
    1223        1197 :   while (!gequal0(y))
    1224             :   {
    1225         427 :     GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
    1226         427 :     x = y; y = z;
    1227             :   }
    1228         385 :   x = gauss_normal(x);
    1229         385 :   if (typ(x) == t_COMPLEX)
    1230             :   {
    1231         196 :     if      (gequal0(gel(x,2))) x = gel(x,1);
    1232         196 :     else if (gequal0(gel(x,1))) x = gel(x,2);
    1233             :   }
    1234         385 :   return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
    1235             : }
    1236             : 
    1237             : static int
    1238        2520 : c_is_rational(GEN x)
    1239        2520 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
    1240             : static GEN
    1241        1288 : c_zero_gcd(GEN c)
    1242             : {
    1243        1288 :   GEN x = gel(c,1), y = gel(c,2);
    1244        1288 :   long tx = typ(x), ty = typ(y);
    1245        1288 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
    1246          42 :   if (tx == t_PADIC || tx == t_INTMOD
    1247          35 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
    1248          35 :   return gauss_gcd(c, gen_0);
    1249             : }
    1250             : 
    1251             : /* gcd(x, 0) */
    1252             : static GEN
    1253     8204365 : zero_gcd(GEN x)
    1254             : {
    1255             :   pari_sp av;
    1256     8204365 :   switch(typ(x))
    1257             :   {
    1258       24268 :     case t_INT: return absi(x);
    1259        1295 :     case t_FRAC: return absfrac(x);
    1260        1288 :     case t_COMPLEX: return c_zero_gcd(x);
    1261        7210 :     case t_REAL: return gen_1;
    1262         735 :     case t_PADIC: return powis(gel(x,2), valp(x));
    1263         210 :     case t_SER: return pol_xn(valp(x), varn(x));
    1264             :     case t_POLMOD: {
    1265        8338 :       GEN d = gel(x,2);
    1266        8338 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
    1267          28 :       return isinexact(d)? zero_gcd(d): gcopy(d);
    1268             :     }
    1269             :     case t_POL:
    1270     7921099 :       if (!isinexact(x)) break;
    1271          14 :       av = avma;
    1272          14 :       return gerepileupto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
    1273             : 
    1274             :     case t_RFRAC:
    1275      217137 :       if (!isinexact(x)) break;
    1276           0 :       av = avma;
    1277           0 :       return gerepileupto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
    1278             :   }
    1279     8161007 :   return gcopy(x);
    1280             : }
    1281             : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
    1282             : static GEN
    1283     8749464 : zero_gcd2(GEN y, GEN z)
    1284             : {
    1285             :   pari_sp av;
    1286     8749464 :   switch(typ(z))
    1287             :   {
    1288     8181890 :     case t_INT: return zero_gcd(y);
    1289             :     case t_INTMOD:
    1290      566531 :       av = avma;
    1291      566531 :       return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
    1292             :     case t_FFELT:
    1293        1043 :       av = avma;
    1294        1043 :       return gerepileupto(av, gmul(y, FF_1(z)));
    1295             :     default:
    1296           0 :       pari_err_TYPE("zero_gcd", z);
    1297             :   }
    1298           0 :   return NULL;
    1299             : }
    1300             : /* tx = t_RFRAC, y considered as constant */
    1301             : static GEN
    1302      864876 : cont_gcd_rfrac(GEN x, GEN y)
    1303             : {
    1304      864876 :   pari_sp av = avma;
    1305      864876 :   GEN cx; x = primitive_part(x, &cx);
    1306      864876 :   return gerepileupto(av, gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2)));
    1307             : }
    1308             : /* tx = t_POL, y considered as constant */
    1309             : static GEN
    1310     2256744 : cont_gcd_pol(GEN x, GEN y)
    1311             : {
    1312     2256744 :   pari_sp av = avma;
    1313     2256744 :   return gerepileupto(av, scalarpol(ggcd(content(x),y), varn(x)));
    1314             : }
    1315             : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
    1316             : static GEN
    1317        6588 : cont_gcd_gen(GEN x, GEN y)
    1318             : {
    1319        6588 :   pari_sp av = avma;
    1320        6588 :   return gerepileupto(av, ggcd(content(x),y));
    1321             : }
    1322             : /* !is_const(tx), y considered as constant */
    1323             : static GEN
    1324     3071427 : cont_gcd(GEN x, long tx, GEN y)
    1325             : {
    1326     3071427 :   switch(tx)
    1327             :   {
    1328      808109 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1329     2256730 :     case t_POL: return cont_gcd_pol(x,y);
    1330        6588 :     default: return cont_gcd_gen(x,y);
    1331             :   }
    1332             : }
    1333             : static GEN
    1334     9780703 : gcdiq(GEN x, GEN y)
    1335             : {
    1336             :   GEN z;
    1337     9780703 :   if (!signe(x)) return Q_abs(y);
    1338     2458920 :   z = cgetg(3,t_FRAC);
    1339     2458920 :   gel(z,1) = gcdii(x,gel(y,1));
    1340     2458920 :   gel(z,2) = icopy(gel(y,2));
    1341     2458920 :   return z;
    1342             : }
    1343             : static GEN
    1344    10946197 : gcdqq(GEN x, GEN y)
    1345             : {
    1346    10946197 :   GEN z = cgetg(3,t_FRAC);
    1347    10946197 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1348    10946197 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1349    10946197 :   return z;
    1350             : }
    1351             : /* assume x,y t_INT or t_FRAC */
    1352             : GEN
    1353   104851665 : Q_gcd(GEN x, GEN y)
    1354             : {
    1355   104851665 :   long tx = typ(x), ty = typ(y);
    1356   104851665 :   if (tx == t_INT)
    1357    84909721 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1358             :   else
    1359    19941944 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1360             : }
    1361             : 
    1362             : GEN
    1363    24077207 : ggcd(GEN x, GEN y)
    1364             : {
    1365    24077207 :   long vx, vy, tx = typ(x), ty = typ(y);
    1366             :   pari_sp av, tetpil;
    1367             :   GEN p1,z;
    1368             : 
    1369    48154414 :   if (is_noncalc_t(tx) || is_matvec_t(tx) ||
    1370    48154414 :       is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
    1371    24077207 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1372             :   /* tx <= ty */
    1373    24077207 :   z = gisexactzero(x); if (z) return zero_gcd2(y,z);
    1374    21032365 :   z = gisexactzero(y); if (z) return zero_gcd2(x,z);
    1375    15327743 :   if (is_const_t(tx))
    1376             :   {
    1377     8391780 :     if (ty == tx) switch(tx)
    1378             :     {
    1379             :       case t_INT:
    1380     3210277 :         return gcdii(x,y);
    1381             : 
    1382     2074065 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1383     2074065 :         if (equalii(gel(x,1),gel(y,1)))
    1384     2074065 :           gel(z,1) = icopy(gel(x,1));
    1385             :         else
    1386           0 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1387     2074065 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1388             :         else
    1389             :         {
    1390     2074065 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1391     2074065 :           if (!is_pm1(p1))
    1392             :           {
    1393           0 :             p1 = gcdii(p1,gel(y,2));
    1394           0 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1395           0 :             else p1 = gerepileuptoint(av, p1);
    1396             :           }
    1397     2074065 :           gel(z,2) = p1;
    1398             :         }
    1399     2074065 :         return z;
    1400             : 
    1401             :       case t_FRAC:
    1402       32571 :         return gcdqq(x,y);
    1403             : 
    1404             :       case t_FFELT:
    1405        8841 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1406        8841 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1407             : 
    1408             :       case t_COMPLEX:
    1409          21 :         if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
    1410           7 :         return triv_cont_gcd(y,x);
    1411             : 
    1412             :       case t_PADIC:
    1413          14 :         if (!equalii(gel(x,2),gel(y,2))) return gen_1;
    1414           7 :         return powis(gel(y,2), minss(valp(x), valp(y)));
    1415             : 
    1416             :       case t_QUAD:
    1417         133 :         av=avma; p1=gdiv(x,y);
    1418         133 :         if (gequal0(gel(p1,3)))
    1419             :         {
    1420          14 :           p1=gel(p1,2);
    1421          14 :           if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
    1422           7 :           tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
    1423             :         }
    1424         119 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
    1425         112 :         p1 = ginv(p1); avma=av;
    1426         112 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1427         105 :         return triv_cont_gcd(y,x);
    1428             : 
    1429           0 :       default: return gen_1; /* t_REAL */
    1430             :     }
    1431     3065858 :     if (is_const_t(ty)) switch(tx)
    1432             :     {
    1433             :       case t_INT:
    1434       21521 :         switch(ty)
    1435             :         {
    1436        1386 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1437        1386 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1438        1386 :             p1 = gcdii(gel(y,1),gel(y,2));
    1439        1386 :             if (!is_pm1(p1)) {
    1440          14 :               p1 = gcdii(x,p1);
    1441          14 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1442             :               else
    1443          14 :                 p1 = gerepileuptoint(av, p1);
    1444             :             }
    1445        1386 :             gel(z,2) = p1; return z;
    1446             : 
    1447        6825 :           case t_REAL: return gen_1;
    1448             : 
    1449             :           case t_FRAC:
    1450        7889 :             return gcdiq(x,y);
    1451             : 
    1452             :           case t_COMPLEX:
    1453        2401 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1454        2072 :             return triv_cont_gcd(y,x);
    1455             : 
    1456             :           case t_FFELT:
    1457         224 :             if (!FF_equal0(y)) return FF_1(y);
    1458           0 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1459             : 
    1460             :           case t_PADIC:
    1461        2670 :             return padic_gcd(x,y);
    1462             : 
    1463             :           case t_QUAD:
    1464         126 :             return triv_cont_gcd(y,x);
    1465             :           default:
    1466           0 :             pari_err_TYPE2("gcd",x,y);
    1467             :         }
    1468             : 
    1469             :       case t_REAL:
    1470          14 :         switch(ty)
    1471             :         {
    1472             :           case t_INTMOD:
    1473             :           case t_FFELT:
    1474          14 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1475           0 :           default: return gen_1;
    1476             :         }
    1477             : 
    1478             :       case t_INTMOD:
    1479          49 :         switch(ty)
    1480             :         {
    1481             :           case t_FRAC:
    1482          14 :             av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
    1483          14 :             if (!is_pm1(p1)) pari_err_OP("gcd",x,y);
    1484           7 :             return ggcd(gel(y,1), x);
    1485             : 
    1486             :           case t_FFELT:
    1487             :           {
    1488          14 :             GEN p = gel(y,4);
    1489          14 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1490           7 :             if (!FF_equal0(y)) return FF_1(y);
    1491           0 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1492             :           }
    1493             : 
    1494             :           case t_COMPLEX: case t_QUAD:
    1495          14 :             return triv_cont_gcd(y,x);
    1496             : 
    1497             :           case t_PADIC:
    1498           7 :             return padic_gcd(x,y);
    1499             : 
    1500           0 :           default: pari_err_TYPE2("gcd",x,y);
    1501             :         }
    1502             : 
    1503             :       case t_FRAC:
    1504         210 :         switch(ty)
    1505             :         {
    1506             :           case t_COMPLEX:
    1507          84 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1508             :           case t_QUAD:
    1509         154 :             return triv_cont_gcd(y,x);
    1510             :           case t_FFELT:
    1511             :           {
    1512          42 :             GEN p = gel(y,4);
    1513          42 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1514          21 :             if (!FF_equal0(y)) return FF_1(y);
    1515           0 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1516             :           }
    1517             : 
    1518             :           case t_PADIC:
    1519           7 :             return padic_gcd(x,y);
    1520             : 
    1521           0 :           default: pari_err_TYPE2("gcd",x,y);
    1522             :         }
    1523             :       case t_FFELT:
    1524          70 :         switch(ty)
    1525             :         {
    1526             :           case t_PADIC:
    1527             :           {
    1528          42 :             GEN p = gel(y,2);
    1529          42 :             long v = valp(y);
    1530          42 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1531          14 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1532             :           }
    1533          28 :           default: pari_err_TYPE2("gcd",x,y);
    1534             :         }
    1535             : 
    1536             :       case t_COMPLEX:
    1537          14 :         switch(ty)
    1538             :         {
    1539             :           case t_PADIC:
    1540          14 :           case t_QUAD: return triv_cont_gcd(x,y);
    1541           0 :           default: pari_err_TYPE2("gcd",x,y);
    1542             :         }
    1543             : 
    1544             :       case t_PADIC:
    1545           7 :         switch(ty)
    1546             :         {
    1547           7 :           case t_QUAD: return triv_cont_gcd(y,x);
    1548           0 :           default: pari_err_TYPE2("gcd",x,y);
    1549             :         }
    1550             : 
    1551           0 :       default: return gen_1; /* tx = t_REAL */
    1552             :     }
    1553     3043973 :     return cont_gcd(y,ty, x);
    1554             :   }
    1555             : 
    1556     6935963 :   if (tx == t_POLMOD)
    1557             :   {
    1558       19358 :     if (ty == t_POLMOD)
    1559             :     {
    1560       19302 :       GEN T = gel(x,1);
    1561       19302 :       z = cgetg(3,t_POLMOD);
    1562       19302 :       T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
    1563       19302 :       gel(z,1) = T;
    1564       19302 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1565             :       else
    1566             :       {
    1567             :         GEN X, Y, d;
    1568       19302 :         av = avma; X = gel(x,2); Y = gel(y,2);
    1569       19302 :         d = ggcd(content(X), content(Y));
    1570       19302 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1571       19302 :         p1 = ggcd(T, X);
    1572       19302 :         gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
    1573             :       }
    1574       19302 :       return z;
    1575             :     }
    1576          56 :     vx = varn(gel(x,1));
    1577          56 :     switch(ty)
    1578             :     {
    1579             :       case t_POL:
    1580          28 :         vy = varn(y);
    1581          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1582          14 :         z = cgetg(3,t_POLMOD);
    1583          14 :         gel(z,1) = RgX_copy(gel(x,1));
    1584          14 :         av = avma; p1 = ggcd(gel(x,1),gel(x,2));
    1585          14 :         gel(z,2) = gerepileupto(av, ggcd(p1,y));
    1586          14 :         return z;
    1587             : 
    1588             :       case t_RFRAC:
    1589          28 :         vy = varn(gel(y,2));
    1590          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1591          28 :         av = avma;
    1592          28 :         p1 = ggcd(gel(x,1),gel(y,2));
    1593          28 :         if (degpol(p1)) pari_err_OP("gcd",x,y);
    1594          21 :         avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1595             :     }
    1596             :   }
    1597             : 
    1598     6916605 :   vx = gvar(x);
    1599     6916605 :   vy = gvar(y);
    1600     6916605 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1601     6904600 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1602             : 
    1603             :   /* vx = vy: same main variable */
    1604     6889151 :   switch(tx)
    1605             :   {
    1606             :     case t_POL:
    1607     6728245 :       switch(ty)
    1608             :       {
    1609     6671464 :         case t_POL: return RgX_gcd(x,y);
    1610             :         case t_SER:
    1611          14 :           z = ggcd(content(x), content(y));
    1612          14 :           return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
    1613       56767 :         case t_RFRAC: return cont_gcd_rfrac(y, x);
    1614             :       }
    1615           0 :       break;
    1616             : 
    1617             :     case t_SER:
    1618          14 :       z = ggcd(content(x), content(y));
    1619          14 :       switch(ty)
    1620             :       {
    1621           7 :         case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
    1622           7 :         case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
    1623             :       }
    1624           0 :       break;
    1625             : 
    1626             :     case t_RFRAC:
    1627             :     {
    1628      160892 :       GEN xd = gel(x,2), yd = gel(y,2);
    1629      160892 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1630      160892 :       z = cgetg(3,t_RFRAC); av = avma;
    1631      160892 :       gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1632      160892 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1633             :     }
    1634             :   }
    1635           0 :   pari_err_TYPE2("gcd",x,y);
    1636             :   return NULL; /* LCOV_EXCL_LINE */
    1637             : }
    1638             : GEN
    1639        3960 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
    1640             : 
    1641             : static GEN
    1642         133 : fix_lcm(GEN x)
    1643             : {
    1644             :   GEN t;
    1645         133 :   switch(typ(x))
    1646             :   {
    1647          35 :     case t_INT: if (signe(x)<0) x = negi(x);
    1648          35 :       break;
    1649             :     case t_POL:
    1650          98 :       if (lg(x) <= 2) break;
    1651          98 :       t = leading_coeff(x);
    1652          98 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1653             :   }
    1654         133 :   return x;
    1655             : }
    1656             : GEN
    1657        2905 : glcm0(GEN x, GEN y)
    1658             : {
    1659        2905 :   if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
    1660        2821 :   return glcm(x,y);
    1661             : }
    1662             : GEN
    1663        3199 : glcm(GEN x, GEN y)
    1664             : {
    1665             :   pari_sp av;
    1666             :   GEN z;
    1667        3199 :   if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
    1668          63 :   av = avma; z = ggcd(x,y);
    1669          63 :   if (!gequal1(z))
    1670             :   {
    1671          63 :     if (gequal0(z)) { avma = av; return gmul(x,y); }
    1672          49 :     y = gdiv(y,z);
    1673             :   }
    1674          49 :   return gerepileupto(av, fix_lcm(gmul(x,y)));
    1675             : }
    1676             : 
    1677             : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
    1678             : static int
    1679      302589 : pol_approx0(GEN r, GEN x, int exact)
    1680             : {
    1681             :   long i, lx,lr;
    1682      302589 :   if (exact) return gequal0(r);
    1683           0 :   lx = lg(x);
    1684           0 :   lr = lg(r); if (lr < lx) lx = lr;
    1685           0 :   for (i=2; i<lx; i++)
    1686           0 :     if (!approx_0(gel(r,i), gel(x,i))) return 0;
    1687           0 :   return 1;
    1688             : }
    1689             : 
    1690             : GEN
    1691       87164 : RgX_gcd_simple(GEN x, GEN y)
    1692             : {
    1693       87164 :   pari_sp av1, av = avma;
    1694       87164 :   GEN r, yorig = y;
    1695       87164 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1696             : 
    1697             :   for(;;)
    1698             :   {
    1699      302589 :     av1 = avma; r = RgX_rem(x,y);
    1700      302589 :     if (pol_approx0(r, x, exact))
    1701             :     {
    1702       87164 :       avma = av1;
    1703       87164 :       if (y == yorig) return RgX_copy(y);
    1704       60942 :       y = normalizepol_approx(y, lg(y));
    1705       60942 :       if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
    1706        4340 :       return gerepileupto(av,y);
    1707             :     }
    1708      215425 :     x = y; y = r;
    1709      215425 :     if (gc_needed(av,1)) {
    1710           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1711           0 :       gerepileall(av,2, &x,&y);
    1712             :     }
    1713      215425 :   }
    1714             : }
    1715             : GEN
    1716           0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
    1717             : {
    1718           0 :   pari_sp av = avma;
    1719             :   GEN q, r, d, d1, u, v, v1;
    1720           0 :   int exact = !(isinexactreal(a) || isinexactreal(b));
    1721             : 
    1722           0 :   d = a; d1 = b; v = gen_0; v1 = gen_1;
    1723             :   for(;;)
    1724             :   {
    1725           0 :     if (pol_approx0(d1, a, exact)) break;
    1726           0 :     q = poldivrem(d,d1, &r);
    1727           0 :     v = gsub(v, gmul(q,v1));
    1728           0 :     u=v; v=v1; v1=u;
    1729           0 :     u=r; d=d1; d1=u;
    1730           0 :   }
    1731           0 :   u = gsub(d, gmul(b,v));
    1732           0 :   u = RgX_div(u,a);
    1733             : 
    1734           0 :   gerepileall(av, 3, &u,&v,&d);
    1735           0 :   *pu = u;
    1736           0 :   *pv = v; return d;
    1737             : }
    1738             : /*******************************************************************/
    1739             : /*                                                                 */
    1740             : /*                    CONTENT / PRIMITIVE PART                     */
    1741             : /*                                                                 */
    1742             : /*******************************************************************/
    1743             : 
    1744             : GEN
    1745    77308784 : content(GEN x)
    1746             : {
    1747    77308784 :   long lx, i, t, tx = typ(x);
    1748    77308784 :   pari_sp av = avma;
    1749             :   GEN c;
    1750             : 
    1751    77308784 :   if (is_scalar_t(tx)) return zero_gcd(x);
    1752    77295703 :   switch(tx)
    1753             :   {
    1754             :     case t_RFRAC:
    1755             :     {
    1756      864890 :       GEN n = gel(x,1), d = gel(x,2);
    1757             :       /* -- varncmp(vn, vd) < 0 can't happen
    1758             :        * -- if n is POLMOD, its main variable (in the sense of gvar2)
    1759             :        *    has lower priority than denominator */
    1760      864890 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    1761      824566 :         n = isinexact(n)? zero_gcd(n): gcopy(n);
    1762             :       else
    1763       40324 :         n = content(n);
    1764      864890 :       return gerepileupto(av, gdiv(n, content(d)));
    1765             :     }
    1766             : 
    1767             :     case t_VEC: case t_COL:
    1768     7566135 :       lx = lg(x); if (lx==1) return gen_0;
    1769     7566128 :       break;
    1770             : 
    1771             :     case t_MAT:
    1772             :     {
    1773             :       long hx, j;
    1774         238 :       lx = lg(x);
    1775         238 :       if (lx == 1) return gen_0;
    1776         231 :       hx = lgcols(x);
    1777         231 :       if (hx == 1) return gen_0;
    1778         224 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    1779         224 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    1780         224 :       c = content(gel(x,1));
    1781         448 :       for (j=2; j<lx; j++)
    1782         224 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    1783         224 :       if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
    1784         224 :       return gerepileupto(av,c);
    1785             :     }
    1786             : 
    1787             :     case t_POL: case t_SER:
    1788    68864405 :       lx = lg(x); if (lx == 2) return gen_0;
    1789    68864391 :       break;
    1790          21 :     case t_VECSMALL: return utoi(zv_content(x));
    1791             :     case t_QFR: case t_QFI:
    1792          14 :       lx = 4; break;
    1793             : 
    1794           0 :     default: pari_err_TYPE("content",x);
    1795             :       return NULL; /* LCOV_EXCL_LINE */
    1796             :   }
    1797   224952391 :   for (i=lontyp[tx]; i<lx; i++)
    1798   161441355 :     if (typ(gel(x,i)) != t_INT) break;
    1799    76430533 :   lx--; c = gel(x,lx);
    1800    76430533 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    1801    76430533 :   if (i > lx)
    1802             :   { /* integer coeffs */
    1803   128282365 :     while (lx-- > lontyp[tx])
    1804             :     {
    1805    64013813 :       c = gcdii(c, gel(x,lx));
    1806    64013813 :       if (is_pm1(c)) { avma=av; return gen_1; }
    1807             :     }
    1808             :   }
    1809             :   else
    1810             :   {
    1811    12919497 :     if (isinexact(c)) c = zero_gcd(c);
    1812    46558949 :     while (lx-- > lontyp[tx])
    1813             :     {
    1814    20719955 :       GEN d = gel(x,lx);
    1815    20719955 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    1816    20719955 :       c = ggcd(c, d);
    1817             :     }
    1818    12919497 :     if (isinexact(c)) { avma=av; return gen_1; }
    1819             :   }
    1820    13677013 :   switch(typ(c))
    1821             :   {
    1822             :     case t_INT:
    1823      769374 :       if (signe(c) < 0) c = negi(c);
    1824      769374 :       break;
    1825             :     case t_VEC: case t_COL: case t_MAT:
    1826           0 :       pari_err_TYPE("content",x);
    1827             :   }
    1828             : 
    1829    13677013 :   return av==avma? gcopy(c): gerepileupto(av,c);
    1830             : }
    1831             : 
    1832             : GEN
    1833     1007796 : primitive_part(GEN x, GEN *ptc)
    1834             : {
    1835     1007796 :   pari_sp av = avma;
    1836     1007796 :   GEN c = content(x);
    1837     1007796 :   if (gequal1(c)) { avma = av; c = NULL; }
    1838       33135 :   else if (!gequal0(c)) x = gdiv(x,c);
    1839     1007796 :   if (ptc) *ptc = c;
    1840     1007796 :   return x;
    1841             : }
    1842             : GEN
    1843        2835 : primpart(GEN x) { return primitive_part(x, NULL); }
    1844             : 
    1845             : static GEN
    1846    24079382 : Q_content_v(GEN x, long i, long l)
    1847             : {
    1848    24079382 :   pari_sp av = avma;
    1849    24079382 :   GEN d = Q_content_safe(gel(x,i));
    1850    24079382 :   if (!d) return NULL;
    1851   125570630 :   for (i++; i < l; i++)
    1852             :   {
    1853   101491451 :     GEN c = Q_content_safe(gel(x,i));
    1854   101491451 :     if (!c) return NULL;
    1855   101491451 :     d = Q_gcd(d, c);
    1856             :   }
    1857    24079179 :   return gerepileupto(av, d);
    1858             : }
    1859             : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1860             :  * of Q(x2,...,xn)[x1] */
    1861             : GEN
    1862   140209274 : Q_content_safe(GEN x)
    1863             : {
    1864             :   long l;
    1865   140209274 :   switch(typ(x))
    1866             :   {
    1867   101892561 :     case t_INT:  return absi(x);
    1868    13939592 :     case t_FRAC: return absfrac(x);
    1869             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    1870    12893790 :       l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
    1871             :     case t_POL:
    1872    11461099 :       l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
    1873       22092 :     case t_POLMOD: return Q_content_safe(gel(x,2));
    1874             :   }
    1875         140 :   return NULL;
    1876             : }
    1877             : GEN
    1878       99951 : Q_content(GEN x)
    1879             : {
    1880       99951 :   GEN c = Q_content_safe(x);
    1881       99951 :   if (!c) pari_err_TYPE("Q_content",x);
    1882       99951 :   return c;
    1883             : }
    1884             : 
    1885             : GEN
    1886       38678 : ZX_content(GEN x)
    1887             : {
    1888       38678 :   long i, l = lg(x);
    1889             :   GEN d;
    1890             :   pari_sp av;
    1891             : 
    1892       38678 :   if (l == 2) return gen_0;
    1893       38678 :   d = gel(x,2);
    1894       38678 :   if (l == 3) return absi(d);
    1895       32775 :   av = avma;
    1896       32775 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    1897       32775 :   if (signe(d) < 0) d = absi(d);
    1898       32775 :   return gerepileuptoint(av, d);
    1899             : }
    1900             : 
    1901             : static GEN
    1902       82110 : Z_content_v(GEN x, long i, long l)
    1903             : {
    1904       82110 :   pari_sp av = avma;
    1905       82110 :   GEN d = Z_content(gel(x,i));
    1906       82110 :   if (!d) return NULL;
    1907      356867 :   for (i++; i<l; i++)
    1908             :   {
    1909      277984 :     GEN c = Z_content(gel(x,i));
    1910      277984 :     if (!c) return NULL;
    1911      277858 :     d = gcdii(d, c); if (is_pm1(d)) return NULL;
    1912      277816 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    1913             :   }
    1914       78883 :   return gerepileuptoint(av, d);
    1915             : }
    1916             : /* return NULL for 1 */
    1917             : GEN
    1918      361326 : Z_content(GEN x)
    1919             : {
    1920             :   long l;
    1921      361326 :   switch(typ(x))
    1922             :   {
    1923             :     case t_INT:
    1924      264110 :       if (is_pm1(x)) return NULL;
    1925      263340 :       return absi(x);
    1926             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    1927        7903 :       l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
    1928             :     case t_POL:
    1929       89313 :       l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
    1930           0 :     case t_POLMOD: return Z_content(gel(x,2));
    1931             :   }
    1932           0 :   pari_err_TYPE("Z_content", x);
    1933             :   return NULL; /* LCOV_EXCL_LINE */
    1934             : }
    1935             : 
    1936             : static GEN
    1937    10761295 : Q_denom_v(GEN x, long i, long l)
    1938             : {
    1939    10761295 :   pari_sp av = avma;
    1940    10761295 :   GEN d = Q_denom(gel(x,i));
    1941    46591719 :   for (i++; i<l; i++)
    1942             :   {
    1943    35830424 :     GEN D = Q_denom(gel(x,i));
    1944    35830424 :     if (D != gen_1) d = lcmii(d, D);
    1945    35830424 :     if ((i & 255) == 0) d = gerepileuptoint(av, d);
    1946             :   }
    1947    10761295 :   return gerepileuptoint(av, d);
    1948             : }
    1949             : /* NOT MEMORY CLEAN (because of t_FRAC).
    1950             :  * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1951             :  * of Q(x2,...,xn)[x1] */
    1952             : GEN
    1953    69170418 : Q_denom(GEN x)
    1954             : {
    1955             :   long l;
    1956    69170418 :   switch(typ(x))
    1957             :   {
    1958    45455912 :     case t_INT: return gen_1;
    1959    12806358 :     case t_FRAC: return gel(x,2);
    1960             :     case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
    1961     8191698 :       l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
    1962             :     case t_POL:
    1963     2640514 :       l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
    1964       75936 :     case t_POLMOD: return Q_denom(gel(x,2));
    1965             :   }
    1966           0 :   pari_err_TYPE("Q_denom",x);
    1967             :   return NULL; /* LCOV_EXCL_LINE */
    1968             : }
    1969             : 
    1970             : GEN
    1971     8184865 : Q_remove_denom(GEN x, GEN *ptd)
    1972             : {
    1973     8184865 :   GEN d = Q_denom(x);
    1974     8184865 :   if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d);
    1975     8184865 :   if (ptd) *ptd = d;
    1976     8184865 :   return x;
    1977             : }
    1978             : 
    1979             : /* return y = x * d, assuming x rational, and d,y integral */
    1980             : GEN
    1981    44009800 : Q_muli_to_int(GEN x, GEN d)
    1982             : {
    1983             :   long i, l;
    1984             :   GEN y, xn, xd;
    1985             :   pari_sp av;
    1986             : 
    1987    44009800 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    1988    44009800 :   switch (typ(x))
    1989             :   {
    1990             :     case t_INT:
    1991    18224557 :       return mulii(x,d);
    1992             : 
    1993             :     case t_FRAC:
    1994    18980597 :       xn = gel(x,1);
    1995    18980597 :       xd = gel(x,2); av = avma;
    1996    18980597 :       y = mulii(xn, diviiexact(d, xd));
    1997    18980597 :       return gerepileuptoint(av, y);
    1998             :     case t_COMPLEX:
    1999           0 :       y = cgetg(3,t_COMPLEX);
    2000           0 :       gel(y,1) = Q_muli_to_int(gel(x,1),d);
    2001           0 :       gel(y,2) = Q_muli_to_int(gel(x,2),d);
    2002           0 :       return y;
    2003             : 
    2004             :     case t_VEC: case t_COL: case t_MAT:
    2005     4918135 :       y = cgetg_copy(x, &l);
    2006     4918135 :       for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2007     4918135 :       return y;
    2008             : 
    2009             :     case t_POL:
    2010     1844840 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2011     1844840 :       for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    2012     1844840 :       return y;
    2013             : 
    2014             :     case t_POLMOD:
    2015       41671 :       y = cgetg(3, t_POLMOD);
    2016       41671 :       gel(y,1) = RgX_copy(gel(x,1));
    2017       41671 :       gel(y,2) = Q_muli_to_int(gel(x,2), d);
    2018       41671 :       return y;
    2019             :   }
    2020           0 :   pari_err_TYPE("Q_muli_to_int",x);
    2021             :   return NULL; /* LCOV_EXCL_LINE */
    2022             : }
    2023             : 
    2024             : static void
    2025     1840356 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
    2026             : {
    2027             :   long e;
    2028     1840356 :   switch(typ(c))
    2029             :   {
    2030             :     case t_REAL:
    2031     1232321 :       *exact = 0;
    2032     1232321 :       if (!signe(c)) return;
    2033     1229942 :       e = expo(c) - bit_prec(c);
    2034     1229942 :       break;
    2035             :     case t_INT:
    2036      607937 :       if (!signe(c)) return;
    2037      248674 :       e = expi(c) + 32;
    2038      248674 :       break;
    2039             :     case t_FRAC:
    2040          70 :       e = expi(gel(c,1)) - expi(gel(c,2)) + 32;
    2041          70 :       if (exact) *D = lcmii(*D, gel(c,2));
    2042          70 :       break;
    2043             :     default:
    2044          28 :       pari_err_TYPE("rescale_to_int",c);
    2045             :       return; /* LCOV_EXCL_LINE */
    2046             :   }
    2047     1478686 :   if (e < *emin) *emin = e;
    2048             : }
    2049             : GEN
    2050      223019 : RgM_rescale_to_int(GEN x)
    2051             : {
    2052      223019 :   long lx = lg(x), i,j, hx, emin;
    2053             :   GEN D;
    2054             :   int exact;
    2055             : 
    2056      223019 :   if (lx == 1) return cgetg(1,t_MAT);
    2057      223019 :   hx = lgcols(x);
    2058      223019 :   exact = 1;
    2059      223019 :   emin = HIGHEXPOBIT;
    2060      223019 :   D = gen_1;
    2061      808750 :   for (j = 1; j < lx; j++)
    2062      585759 :     for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
    2063      222991 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2064      222935 :   return grndtoi(gmul2n(x, -emin), &i);
    2065             : }
    2066             : GEN
    2067         238 : RgX_rescale_to_int(GEN x)
    2068             : {
    2069         238 :   long lx = lg(x), i, emin;
    2070             :   GEN D;
    2071             :   int exact;
    2072         238 :   if (lx == 2) return gcopy(x); /* rare */
    2073         238 :   exact = 1;
    2074         238 :   emin = HIGHEXPOBIT;
    2075         238 :   D = gen_1;
    2076         238 :   for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
    2077         238 :   if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
    2078         238 :   return grndtoi(gmul2n(x, -emin), &i);
    2079             : }
    2080             : 
    2081             : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    2082             : static GEN
    2083     8259014 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    2084             : {
    2085             :   long i, l;
    2086             :   GEN y, xn, xd;
    2087             :   pari_sp av;
    2088             : 
    2089     8259014 :   switch(typ(x))
    2090             :   {
    2091             :     case t_INT:
    2092     2493594 :       av = avma; y = diviiexact(x,d);
    2093     2493594 :       return gerepileuptoint(av, mulii(y,n));
    2094             : 
    2095             :     case t_FRAC:
    2096     3974819 :       xn = gel(x,1);
    2097     3974819 :       xd = gel(x,2); av = avma;
    2098     3974819 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    2099     3974819 :       return gerepileuptoint(av, y);
    2100             : 
    2101             :     case t_VEC: case t_COL: case t_MAT:
    2102      427444 :       y = cgetg_copy(x, &l);
    2103      427444 :       for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2104      427444 :       return y;
    2105             : 
    2106             :     case t_POL:
    2107     1363157 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2108     1363157 :       for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    2109     1363157 :       return y;
    2110             : 
    2111             :     case t_POLMOD:
    2112           0 :       y = cgetg(3, t_POLMOD);
    2113           0 :       gel(y,1) = RgX_copy(gel(x,1));
    2114           0 :       gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
    2115           0 :       return y;
    2116             :   }
    2117           0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    2118             :   return NULL; /* LCOV_EXCL_LINE */
    2119             : }
    2120             : 
    2121             : /* return x / d. x: rational; d,result: integral. */
    2122             : static GEN
    2123    16975303 : Q_divi_to_int(GEN x, GEN d)
    2124             : {
    2125             :   long i, l;
    2126             :   GEN y;
    2127             : 
    2128    16975303 :   switch(typ(x))
    2129             :   {
    2130             :     case t_INT:
    2131    13835118 :       return diviiexact(x,d);
    2132             : 
    2133             :     case t_VEC: case t_COL: case t_MAT:
    2134     1308474 :       y = cgetg_copy(x, &l);
    2135     1308474 :       for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2136     1308474 :       return y;
    2137             : 
    2138             :     case t_POL:
    2139     1828085 :       y = cgetg_copy(x, &l); y[1] = x[1];
    2140     1828085 :       for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    2141     1828085 :       return y;
    2142             : 
    2143             :     case t_POLMOD:
    2144        3626 :       y = cgetg(3, t_POLMOD);
    2145        3626 :       gel(y,1) = RgX_copy(gel(x,1));
    2146        3626 :       gel(y,2) = Q_divi_to_int(gel(x,2), d);
    2147        3626 :       return y;
    2148             :   }
    2149           0 :   pari_err_TYPE("Q_divi_to_int",x);
    2150             :   return NULL; /* LCOV_EXCL_LINE */
    2151             : }
    2152             : /* c t_FRAC */
    2153             : static GEN
    2154     3025546 : Q_divq_to_int(GEN x, GEN c)
    2155             : {
    2156     3025546 :   GEN n = gel(c,1), d = gel(c,2);
    2157     3025546 :   if (is_pm1(n)) {
    2158     1608470 :     GEN y = Q_muli_to_int(x,d);
    2159     1608470 :     if (signe(n) < 0) y = gneg(y);
    2160     1608470 :     return y;
    2161             :   }
    2162     1417076 :   return Q_divmuli_to_int(x, n,d);
    2163             : }
    2164             : 
    2165             : /* return y = x / c, assuming x,c rational, and y integral */
    2166             : GEN
    2167        7486 : Q_div_to_int(GEN x, GEN c)
    2168             : {
    2169        7486 :   switch(typ(c))
    2170             :   {
    2171        6275 :     case t_INT:  return Q_divi_to_int(x, c);
    2172        1211 :     case t_FRAC: return Q_divq_to_int(x, c);
    2173             :   }
    2174           0 :   pari_err_TYPE("Q_div_to_int",c);
    2175             :   return NULL; /* LCOV_EXCL_LINE */
    2176             : }
    2177             : /* return y = x * c, assuming x,c rational, and y integral */
    2178             : GEN
    2179           0 : Q_mul_to_int(GEN x, GEN c)
    2180             : {
    2181             :   GEN d, n;
    2182           0 :   switch(typ(c))
    2183             :   {
    2184           0 :     case t_INT: return Q_muli_to_int(x, c);
    2185             :     case t_FRAC:
    2186           0 :       n = gel(c,1);
    2187           0 :       d = gel(c,2);
    2188           0 :       return Q_divmuli_to_int(x, d,n);
    2189             :   }
    2190           0 :   pari_err_TYPE("Q_mul_to_int",c);
    2191             :   return NULL; /* LCOV_EXCL_LINE */
    2192             : }
    2193             : 
    2194             : GEN
    2195    14507483 : Q_primitive_part(GEN x, GEN *ptc)
    2196             : {
    2197    14507483 :   pari_sp av = avma;
    2198    14507483 :   GEN c = Q_content_safe(x);
    2199    14507483 :   if (c)
    2200             :   {
    2201    14507343 :     if (typ(c) == t_INT)
    2202             :     {
    2203    11483008 :       if (is_pm1(c)) { avma = av; c = NULL; }
    2204     2148312 :       else if (signe(c)) x = Q_divi_to_int(x, c);
    2205             :     }
    2206     3024335 :     else x = Q_divq_to_int(x, c);
    2207             :   }
    2208    14507483 :   if (ptc) *ptc = c;
    2209    14507483 :   return x;
    2210             : }
    2211             : GEN
    2212     1445084 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    2213             : 
    2214             : GEN
    2215       68618 : vec_Q_primpart(GEN x)
    2216       68618 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
    2217             : 
    2218             : /*******************************************************************/
    2219             : /*                                                                 */
    2220             : /*                           SUBRESULTANT                          */
    2221             : /*                                                                 */
    2222             : /*******************************************************************/
    2223             : /* for internal use */
    2224             : GEN
    2225     1791241 : gdivexact(GEN x, GEN y)
    2226             : {
    2227             :   long i,lx;
    2228             :   GEN z;
    2229     1791241 :   if (gequal1(y)) return x;
    2230     1787384 :   switch(typ(x))
    2231             :   {
    2232             :     case t_INT:
    2233     1431206 :       if (typ(y)==t_INT) return diviiexact(x,y);
    2234          35 :       if (!signe(x)) return gen_0;
    2235           0 :       break;
    2236             :     case t_INTMOD:
    2237             :     case t_FFELT:
    2238       13696 :     case t_POLMOD: return gmul(x,ginv(y));
    2239             :     case t_POL:
    2240      338506 :       switch(typ(y))
    2241             :       {
    2242             :         case t_INTMOD:
    2243             :         case t_FFELT:
    2244        1344 :         case t_POLMOD: return gmul(x,ginv(y));
    2245             :         case t_POL: { /* not stack-clean */
    2246             :           long v;
    2247       20209 :           if (varn(x)!=varn(y)) break;
    2248       19243 :           v = RgX_valrem(y,&y);
    2249       19243 :           if (v) x = RgX_shift_shallow(x,-v);
    2250       19243 :           if (!degpol(y)) { y = gel(y,2); break; }
    2251       17164 :           return RgX_div(x,y);
    2252             :         }
    2253             :       }
    2254      319998 :       return RgX_Rg_divexact(x, y);
    2255             :     case t_VEC: case t_COL: case t_MAT:
    2256        3976 :       lx = lg(x); z = new_chunk(lx);
    2257        3976 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    2258        3976 :       z[0] = x[0]; return z;
    2259             :   }
    2260           0 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    2261           0 :   return gdiv(x,y);
    2262             : }
    2263             : 
    2264             : static GEN
    2265       59856 : init_resultant(GEN x, GEN y)
    2266             : {
    2267       59856 :   long tx = typ(x), ty = typ(y), vx, vy;
    2268       59856 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    2269             :   {
    2270         210 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    2271         210 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    2272           0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    2273           0 :     return gen_1;
    2274             :   }
    2275       59646 :   if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
    2276       59646 :   if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
    2277       59646 :   if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
    2278       59569 :   vx = varn(x);
    2279       59569 :   vy = varn(y); if (vx == vy) return NULL;
    2280           0 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    2281             : }
    2282             : 
    2283             : static long
    2284       96955 : RgX_simpletype(GEN x)
    2285             : {
    2286       96955 :   long T = t_INT, i, lx = lg(x);
    2287      574891 :   for (i = 2; i < lx; i++)
    2288             :   {
    2289      486609 :     GEN c = gel(x,i);
    2290      486609 :     long tc = typ(c);
    2291      486609 :     switch(tc) {
    2292             :       case t_INT:
    2293      415501 :         break;
    2294             :       case t_FRAC:
    2295       32693 :         if (T == t_INT) T = t_FRAC;
    2296       32693 :         break;
    2297             :       default:
    2298       38415 :         if (isinexact(c)) return t_REAL;
    2299       29742 :         T = 0; break;
    2300             :     }
    2301             :   }
    2302       88282 :   return T;
    2303             : }
    2304             : 
    2305             : /* x an RgX, y a scalar */
    2306             : static GEN
    2307           0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    2308             : {
    2309           0 :   *V = gpowgs(y,degpol(x)-1);
    2310           0 :   *U = gen_0; return gmul(y, *V);
    2311             : }
    2312             : 
    2313             : /* return 0 if the subresultant chain can be interrupted.
    2314             :  * Set u = NULL if the resultant is 0. */
    2315             : static int
    2316        5172 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    2317             : {
    2318        5172 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    2319             :   long du, dv, dr, degq;
    2320             : 
    2321        5172 :   if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
    2322        5172 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    2323        5046 :   du = degpol(*u);
    2324        5046 :   dv = degpol(*v);
    2325        5046 :   degq = du - dv;
    2326        5046 :   if (*um1 == gen_1)
    2327        1854 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    2328        3192 :   else if (*um1 == gen_0)
    2329        1596 :     u0 = gen_0;
    2330             :   else /* except in those 2 cases, um1 is an RgX */
    2331        1596 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2332             : 
    2333        5046 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2334        1854 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2335             :   else
    2336        3192 :     u0 = gsub(u0, gmul(q,*uze));
    2337             : 
    2338        5046 :   *um1 = *uze;
    2339        5046 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2340             : 
    2341        5046 :   *u = *v; c = *g; *g  = leading_coeff(*u);
    2342        5046 :   switch(degq)
    2343             :   {
    2344          21 :     case 0: break;
    2345             :     case 1:
    2346        3919 :       c = gmul(*h,c); *h = *g; break;
    2347             :     default:
    2348        1106 :       c = gmul(gpowgs(*h,degq), c);
    2349        1106 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2350             :   }
    2351        5046 :   *v  = RgX_Rg_divexact(r,c);
    2352        5046 :   *uze= RgX_Rg_divexact(*uze,c);
    2353        5046 :   if (both_odd(du, dv)) *signh = -*signh;
    2354        5046 :   return (dr > 3);
    2355             : }
    2356             : 
    2357             : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2358             : static GEN
    2359        1749 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2360             : {
    2361             :   pari_sp av, av2;
    2362        1749 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2363             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2364             : 
    2365        1749 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2366        1749 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2367        1749 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2368        1749 :   if (tx != t_POL) {
    2369           0 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2370           0 :     return scalar_res(y,x,V,U);
    2371             :   }
    2372        1749 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2373        1749 :   if (varn(x) != varn(y))
    2374           0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2375           0 :                                         : scalar_res(y,x,V,U);
    2376        1749 :   if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
    2377        1749 :   if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
    2378        1749 :   dx = degpol(x);
    2379        1749 :   dy = degpol(y);
    2380        1749 :   signh = 1;
    2381        1749 :   if (dx < dy)
    2382             :   {
    2383         874 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2384         874 :     if (both_odd(dx, dy)) signh = -signh;
    2385             :   }
    2386        1749 :   if (dy == 0)
    2387             :   {
    2388           0 :     *V = gpowgs(gel(y,2),dx-1);
    2389           0 :     *U = gen_0; return gmul(*V,gel(y,2));
    2390             :   }
    2391        1749 :   av = avma;
    2392        1749 :   u = x = primitive_part(x, &cu);
    2393        1749 :   v = y = primitive_part(y, &cv);
    2394        1749 :   g = h = gen_1; av2 = avma;
    2395        1749 :   um1 = gen_1; uze = gen_0;
    2396             :   for(;;)
    2397             :   {
    2398        4703 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2399        2954 :     if (gc_needed(av2,1))
    2400             :     {
    2401           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2402           0 :       gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    2403             :     }
    2404        2954 :   }
    2405             :   /* uze an RgX */
    2406        1749 :   if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
    2407        1742 :   z = gel(v,2); du = degpol(u);
    2408        1742 :   if (du > 1)
    2409             :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2410         167 :     p1 = gpowgs(gdiv(z,h),du-1);
    2411         167 :     z = gmul(z,p1);
    2412         167 :     uze = RgX_Rg_mul(uze, p1);
    2413             :   }
    2414        1742 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2415             : 
    2416        1742 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2417        1742 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2418             :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2419        1742 :   p1 = gen_1;
    2420        1742 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2421        1742 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2422        1742 :   cu = cu? gdiv(p1,cu): p1;
    2423        1742 :   cv = cv? gdiv(p1,cv): p1;
    2424        1742 :   z = gmul(z,p1);
    2425        1742 :   *U = RgX_Rg_mul(uze,cu);
    2426        1742 :   *V = RgX_Rg_mul(vze,cv);
    2427        1742 :   return z;
    2428             : }
    2429             : GEN
    2430           0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2431             : {
    2432           0 :   pari_sp av = avma;
    2433           0 :   GEN z = subresext_i(x, y, U, V);
    2434           0 :   gerepileall(av, 3, &z, U, V);
    2435           0 :   return z;
    2436             : }
    2437             : 
    2438             : static GEN
    2439          28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2440             : {
    2441          28 :   GEN x=content(y);
    2442          28 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2443             : }
    2444             : 
    2445             : static int
    2446       90188 : must_negate(GEN x)
    2447             : {
    2448       90188 :   GEN t = leading_coeff(x);
    2449       90188 :   switch(typ(t))
    2450             :   {
    2451             :     case t_INT: case t_REAL:
    2452       57274 :       return (signe(t) < 0);
    2453             :     case t_FRAC:
    2454         364 :       return (signe(gel(t,1)) < 0);
    2455             :   }
    2456       32550 :   return 0;
    2457             : }
    2458             : 
    2459             : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2460             : GEN
    2461         343 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2462             : {
    2463             :   pari_sp av, av2, tetpil;
    2464             :   long signh; /* junk */
    2465         343 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2466             :   GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
    2467             : 
    2468         343 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2469         343 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2470         343 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2471         343 :   vx=varn(x);
    2472         343 :   if (!signe(x))
    2473             :   {
    2474          14 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2475           7 :     *U = pol_0(vx); *V = pol_0(vx);
    2476           7 :     return pol_0(vx);
    2477             :   }
    2478         329 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2479         308 :   dx = degpol(x); dy = degpol(y);
    2480         308 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2481         308 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2482             : 
    2483         161 :   av = avma;
    2484         161 :   u = x = primitive_part(x, &cu);
    2485         161 :   v = y = primitive_part(y, &cv);
    2486         161 :   g = h = gen_1; av2 = avma;
    2487         161 :   um1 = gen_1; uze = gen_0;
    2488             :   for(;;)
    2489             :   {
    2490         182 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2491          21 :     if (gc_needed(av2,1))
    2492             :     {
    2493           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2494           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2495             :     }
    2496          21 :   }
    2497         161 :   if (uze != gen_0) {
    2498             :     GEN r;
    2499          42 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2500          42 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2501          42 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2502          42 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2503          42 :     p1 = ginv(content(v));
    2504             :   }
    2505             :   else /* y | x */
    2506             :   {
    2507         119 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2508         119 :     uze = pol_0(vx);
    2509         119 :     p1 = gen_1;
    2510             :   }
    2511         161 :   if (must_negate(v)) p1 = gneg(p1);
    2512         161 :   tetpil = avma;
    2513         161 :   z = RgX_Rg_mul(v,p1);
    2514         161 :   *U = RgX_Rg_mul(uze,p1);
    2515         161 :   *V = RgX_Rg_mul(vze,p1);
    2516         161 :   gptr[0] = &z;
    2517         161 :   gptr[1] = U;
    2518         161 :   gptr[2] = V;
    2519         161 :   gerepilemanysp(av,tetpil,gptr,3); return z;
    2520             : }
    2521             : 
    2522             : int
    2523          70 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2524             : {
    2525          70 :   pari_sp av = avma, av2, tetpil;
    2526             :   long signh; /* junk */
    2527             :   long vx;
    2528             :   GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
    2529             : 
    2530          70 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2531          70 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2532          70 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2533          70 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2534          70 :   if (!signe(T)) {
    2535           0 :     if (degpol(x) <= amax) {
    2536           0 :       *P = RgX_copy(x);
    2537           0 :       *Q = pol_1(varn(x));
    2538           0 :       return 1;
    2539             :     }
    2540           0 :     return 0;
    2541             :   }
    2542          70 :   if (amax+bmax >= degpol(T))
    2543           0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2544             :                     mkvec3(stoi(amax), stoi(bmax), T));
    2545          70 :   vx = varn(T);
    2546          70 :   u = x = primitive_part(x, &cu);
    2547          70 :   v = T = primitive_part(T, &cv);
    2548          70 :   g = h = gen_1; av2 = avma;
    2549          70 :   um1 = gen_1; uze = gen_0;
    2550             :   for(;;)
    2551             :   {
    2552         287 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2553         287 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
    2554         287 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2555         217 :     if (gc_needed(av2,1))
    2556             :     {
    2557           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2558           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2559             :     }
    2560         217 :   }
    2561          70 :   if (uze == gen_0)
    2562             :   {
    2563           0 :     avma = av; *P = pol_0(vx); *Q = pol_1(vx);
    2564           0 :     return 1;
    2565             :   }
    2566          70 :   if (cu) uze = RgX_Rg_div(uze,cu);
    2567          70 :   p1 = ginv(content(v));
    2568          70 :   if (must_negate(v)) p1 = gneg(p1);
    2569          70 :   tetpil = avma;
    2570          70 :   *P = RgX_Rg_mul(v,p1);
    2571          70 :   *Q = RgX_Rg_mul(uze,p1);
    2572          70 :   gptr[0] = P;
    2573          70 :   gptr[1] = Q;
    2574          70 :   gerepilemanysp(av,tetpil,gptr,2); return 1;
    2575             : }
    2576             : 
    2577             : /*******************************************************************/
    2578             : /*                                                                 */
    2579             : /*                RESULTANT USING DUCOS VARIANT                    */
    2580             : /*                                                                 */
    2581             : /*******************************************************************/
    2582             : /* x^n / y^(n-1), assume n > 0 */
    2583             : static GEN
    2584       19473 : Lazard(GEN x, GEN y, long n)
    2585             : {
    2586             :   long a;
    2587             :   GEN c;
    2588             : 
    2589       19473 :   if (n == 1) return x;
    2590        1397 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    2591        1397 :   c=x; n-=a;
    2592        4814 :   while (a>1)
    2593             :   {
    2594        2020 :     a>>=1; c=gdivexact(gsqr(c),y);
    2595        2020 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    2596             :   }
    2597        1397 :   return c;
    2598             : }
    2599             : 
    2600             : /* F (x/y)^(n-1), assume n >= 1 */
    2601             : static GEN
    2602       19761 : Lazard2(GEN F, GEN x, GEN y, long n)
    2603             : {
    2604       19761 :   if (n == 1) return F;
    2605         735 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    2606             : }
    2607             : 
    2608             : static GEN
    2609       19761 : RgX_neg_i(GEN x, long lx)
    2610             : {
    2611             :   long i;
    2612       19761 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    2613       19761 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    2614       19761 :   return y;
    2615             : }
    2616             : static GEN
    2617       59472 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    2618             : {
    2619             :   long i;
    2620             :   GEN z;
    2621       59472 :   if (isrationalzero(x)) return pol_0(varn(y));
    2622       59451 :   z = cgetg(ly,t_POL); z[1] = y[1];
    2623       59451 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2624       59451 :   return z;
    2625             : }
    2626             : static long
    2627       56560 : reductum_lg(GEN x, long lx)
    2628             : {
    2629       56560 :   long i = lx-2;
    2630       56560 :   while (i > 1 && gequal0(gel(x,i))) i--;
    2631       56560 :   return i+1;
    2632             : }
    2633             : 
    2634             : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
    2635             : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    2636             :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    2637             : static GEN
    2638       19761 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    2639             : {
    2640       19761 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
    2641             :   long p, q, j, lP, lQ;
    2642             :   pari_sp av;
    2643             : 
    2644       19761 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    2645       19761 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    2646             :   /* p > q. Very often p - 1 = q */
    2647       19761 :   av = avma;
    2648             :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    2649       19761 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    2650             : 
    2651       19761 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    2652       23401 :   for (j = q+1; j < p; j++)
    2653             :   {
    2654        3640 :     if (degpol(H) == q-1)
    2655             :     { /* h0 = coeff of degree q-1 = leading coeff */
    2656        3206 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    2657        3206 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    2658             :     }
    2659             :     else
    2660         434 :       H = RgX_shift_shallow(H, 1);
    2661        3640 :     if (j+2 < lP)
    2662             :     {
    2663        3066 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    2664        3066 :       A = A? RgX_add(A, TMP): TMP;
    2665             :     }
    2666        3640 :     if (gc_needed(av,1))
    2667             :     {
    2668         147 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    2669         147 :       gerepileall(av,A?2:1,&H,&A);
    2670             :     }
    2671             :   }
    2672       19761 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    2673       19761 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    2674       19761 :   A = A? RgX_add(A, TMP): TMP;
    2675       19761 :   A = RgX_Rg_divexact(A, p0);
    2676       19761 :   if (degpol(H) == q-1)
    2677             :   {
    2678       19467 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    2679       19467 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    2680             :   }
    2681             :   else
    2682         294 :     A = RgX_Rg_mul(addshift(H,A), q0);
    2683       19761 :   return RgX_Rg_divexact(A, s);
    2684             : }
    2685             : #undef addshift
    2686             : 
    2687             : /* Ducos's subresultant */
    2688             : GEN
    2689       20432 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    2690             : {
    2691             :   pari_sp av, av2;
    2692       20432 :   long dP, dQ, delta, sig = 1;
    2693             :   GEN cP, cQ, Z, s;
    2694             : 
    2695       20432 :   dP = degpol(P);
    2696       20432 :   dQ = degpol(Q); delta = dP - dQ;
    2697       20432 :   if (delta < 0)
    2698             :   {
    2699        1463 :     if (both_odd(dP, dQ)) sig = -1;
    2700        1463 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    2701             :   }
    2702       20432 :   if (sol) *sol = gen_0;
    2703       20432 :   av = avma;
    2704       20432 :   if (dQ <= 0)
    2705             :   {
    2706        1694 :     if (dQ < 0) return Rg_get_0(P);
    2707        1694 :     s = gpowgs(gel(Q,2), dP);
    2708        1694 :     if (sig == -1) s = gerepileupto(av, gneg(s));
    2709        1694 :     return s;
    2710             :   }
    2711             :   /* primitive_part is also possible here, but possibly very costly,
    2712             :    * and hardly ever worth it */
    2713       18738 :   P = Q_primitive_part(P, &cP);
    2714       18738 :   Q = Q_primitive_part(Q, &cQ);
    2715       18738 :   av2 = avma;
    2716       18738 :   s = gpowgs(leading_coeff(Q),delta);
    2717       18738 :   if (both_odd(dP, dQ)) sig = -sig;
    2718       18738 :   Z = Q;
    2719       18738 :   Q = RgX_pseudorem(P, Q);
    2720       18738 :   P = Z;
    2721       57237 :   while(degpol(Q) > 0)
    2722             :   {
    2723       19761 :     delta = degpol(P) - degpol(Q); /* > 0 */
    2724       19761 :     Z = Lazard2(Q, leading_coeff(Q), s, delta);
    2725       19761 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    2726       19761 :     Q = nextSousResultant(P, Q, Z, s);
    2727       19761 :     P = Z;
    2728       19761 :     if (gc_needed(av,1))
    2729             :     {
    2730          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    2731          13 :       gerepileall(av2,2,&P,&Q);
    2732             :     }
    2733       19761 :     s = leading_coeff(P);
    2734             :   }
    2735       18738 :   if (!signe(Q)) { avma = av; return Rg_get_0(Q); }
    2736       18738 :   s = Lazard(leading_coeff(Q), s, degpol(P));
    2737       18738 :   if (sig == -1) s = gneg(s);
    2738       18738 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    2739       18738 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    2740       18738 :   if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
    2741       16995 :   return gerepilecopy(av, s);
    2742             : }
    2743             : /* Return resultant(P,Q).
    2744             :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    2745             :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    2746             :  * in the "generic" case. */
    2747             : GEN
    2748       51176 : resultant(GEN P, GEN Q)
    2749             : {
    2750             :   long TP, TQ;
    2751       51176 :   GEN s, p = NULL;
    2752             : 
    2753       51176 :   if ((s = init_resultant(P,Q))) return s;
    2754       50889 :   if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
    2755        8659 :     return resultant2(P,Q); /* inexact */
    2756       42230 :   if (TP && TQ) /* rational */
    2757             :   {
    2758       24682 :     if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
    2759       10441 :     return QX_resultant(P,Q);
    2760             :   }
    2761       17548 :   if (RgX_is_FpX(P, &p) && RgX_is_FpX(Q, &p) && p)
    2762             :   {
    2763          14 :     pari_sp av = avma;
    2764          14 :     GEN r = FpX_resultant(RgX_to_FpX(P, p), RgX_to_FpX(Q, p), p);
    2765          14 :     return gerepileupto(av, Fp_to_mod(r, p));
    2766             :   }
    2767       17534 :   return RgX_resultant_all(P, Q, NULL);
    2768             : }
    2769             : 
    2770             : /*******************************************************************/
    2771             : /*                                                                 */
    2772             : /*               RESULTANT USING SYLVESTER MATRIX                  */
    2773             : /*                                                                 */
    2774             : /*******************************************************************/
    2775             : static GEN
    2776       72632 : syl_RgC(GEN x, long j, long d, long D, long cp)
    2777             : {
    2778       72632 :   GEN c = cgetg(d+1,t_COL);
    2779             :   long i;
    2780       72632 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    2781       72632 :   for (   ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
    2782       72632 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    2783       72632 :   return c;
    2784             : }
    2785             : static GEN
    2786        8687 : syl_RgM(GEN x, GEN y, long cp)
    2787             : {
    2788        8687 :   long j, d, dx = degpol(x), dy = degpol(y);
    2789             :   GEN M;
    2790        8687 :   if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
    2791        8687 :   if (dy < 0) return zeromat(dx,dx);
    2792        8687 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    2793        8687 :   for (j=1; j<=dy; j++) gel(M,j)    = syl_RgC(x,j,d,j+dx, cp);
    2794        8687 :   for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
    2795        8687 :   return M;
    2796             : }
    2797             : GEN
    2798        8680 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
    2799             : GEN
    2800           7 : sylvestermatrix(GEN x, GEN y)
    2801             : {
    2802           7 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    2803           7 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    2804           7 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    2805           7 :   return syl_RgM(x,y,1);
    2806             : }
    2807             : 
    2808             : GEN
    2809        8680 : resultant2(GEN x, GEN y)
    2810             : {
    2811        8680 :   pari_sp av = avma;
    2812        8680 :   GEN r = init_resultant(x,y);
    2813        8680 :   return r? r: gerepileupto(av, det(RgX_sylvestermatrix(x,y)));
    2814             : }
    2815             : 
    2816             : /* If x a t_POL, let vx = main variable of x; return a t_POL in variable v0:
    2817             :  * if vx <= v, return subst(x, v, pol_x(v0))
    2818             :  * if vx >  v, return scalarpol(x, v0) */
    2819             : static GEN
    2820        3458 : fix_pol(GEN x, long v, long v0)
    2821             : {
    2822             :   long vx;
    2823        3458 :   if (typ(x) != t_POL) return x;
    2824        3458 :   vx = varn(x);
    2825        3458 :   if (v == vx)
    2826             :   {
    2827        3311 :     if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
    2828        3311 :     return x;
    2829             :   }
    2830         147 :   if (varncmp(v, vx) > 0)
    2831             :   {
    2832         140 :     x = gsubst(x,v,pol_x(v0));
    2833         140 :     if (typ(x) == t_POL && varn(x) == v0) return x;
    2834             :   }
    2835          35 :   return scalarpol_shallow(x, v0);
    2836             : }
    2837             : 
    2838             : /* resultant of x and y with respect to variable v, or with respect to their
    2839             :  * main variable if v < 0. */
    2840             : GEN
    2841        1946 : polresultant0(GEN x, GEN y, long v, long flag)
    2842             : {
    2843        1946 :   long v0 = 0;
    2844        1946 :   pari_sp av = avma;
    2845             : 
    2846        1946 :   if (v >= 0)
    2847             :   {
    2848        1715 :     v0 = fetch_var_higher();
    2849        1715 :     x = fix_pol(x,v, v0);
    2850        1715 :     y = fix_pol(y,v, v0);
    2851             :   }
    2852        1946 :   switch(flag)
    2853             :   {
    2854             :     case 2:
    2855        1939 :     case 0: x=resultant(x,y); break;
    2856           7 :     case 1: x=resultant2(x,y); break;
    2857           0 :     default: pari_err_FLAG("polresultant");
    2858             :   }
    2859        1946 :   if (v >= 0) (void)delete_var();
    2860        1946 :   return gerepileupto(av,x);
    2861             : }
    2862             : GEN
    2863         875 : polresultantext0(GEN x, GEN y, long v)
    2864             : {
    2865             :   GEN R, U, V;
    2866         875 :   long v0 = 0;
    2867         875 :   pari_sp av = avma;
    2868             : 
    2869         875 :   if (v >= 0)
    2870             :   {
    2871          14 :     v0 = fetch_var_higher();
    2872          14 :     x = fix_pol(x,v, v0);
    2873          14 :     y = fix_pol(y,v, v0);
    2874             :   }
    2875         875 :   R = subresext_i(x,y, &U,&V);
    2876         875 :   if (v >= 0)
    2877             :   {
    2878          14 :     (void)delete_var();
    2879          14 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    2880          14 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    2881             :   }
    2882         875 :   return gerepilecopy(av, mkvec3(U,V,R));
    2883             : }
    2884             : GEN
    2885         847 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    2886             : 
    2887             : /*******************************************************************/
    2888             : /*                                                                 */
    2889             : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    2890             : /*                                                                 */
    2891             : /*******************************************************************/
    2892             : 
    2893             : /* (v - x)^d */
    2894             : static GEN
    2895         140 : caract_const(pari_sp av, GEN x, long v, long d)
    2896         140 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    2897             : 
    2898             : /* return caract(Mod(x,T)) in variable v */
    2899             : GEN
    2900       14461 : RgXQ_charpoly(GEN x, GEN T, long v)
    2901             : {
    2902       14461 :   pari_sp av = avma;
    2903       14461 :   long d = degpol(T), dx, vx, vp, v0;
    2904             :   GEN ch, L;
    2905             : 
    2906       14461 :   if (typ(x) != t_POL) return caract_const(av, x, v, d);
    2907       14447 :   vx = varn(x);
    2908       14447 :   vp = varn(T);
    2909       14447 :   if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
    2910       14447 :   if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
    2911       14447 :   dx = degpol(x);
    2912       14447 :   if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
    2913       14447 :   if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
    2914             : 
    2915       14391 :   v0 = fetch_var_higher();
    2916       14391 :   x = RgX_neg(x);
    2917       14391 :   gel(x,2) = gadd(gel(x,2), pol_x(v));
    2918       14391 :   setvarn(x, v0);
    2919       14391 :   T = leafcopy(T); setvarn(T, v0);
    2920       14391 :   ch = resultant(T, x);
    2921       14391 :   (void)delete_var();
    2922             :   /* test for silly input: x mod (deg 0 polynomial) */
    2923       14391 :   if (typ(ch) != t_POL)
    2924           7 :     pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
    2925       14384 :   L = leading_coeff(ch);
    2926       14384 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    2927       14384 :   return gerepileupto(av, ch);
    2928             : }
    2929             : 
    2930             : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    2931             :  * algebra nf[t]/(Q(t)) */
    2932             : GEN
    2933         224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    2934             : {
    2935         224 :   const char *f = "rnfcharpoly";
    2936         224 :   long dQ = degpol(Q);
    2937         224 :   pari_sp av = avma;
    2938             :   GEN T;
    2939             : 
    2940         224 :   if (v < 0) v = 0;
    2941         224 :   nf = checknf(nf); T = nf_get_pol(nf);
    2942         224 :   Q = RgX_nffix(f, T,Q,0);
    2943         224 :   switch(typ(x))
    2944             :   {
    2945             :     case t_INT:
    2946          28 :     case t_FRAC: return caract_const(av, x, v, dQ);
    2947             :     case t_POLMOD:
    2948          91 :       x = polmod_nffix2(f,T,Q, x,0);
    2949          49 :       break;
    2950             :     case t_POL:
    2951          56 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    2952          42 :       break;
    2953          49 :     default: pari_err_TYPE(f,x);
    2954             :   }
    2955          91 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    2956             :   /* x a t_POL in variable vQ */
    2957          49 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    2958          49 :   if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
    2959          49 :   return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    2960             : }
    2961             : 
    2962             : /*******************************************************************/
    2963             : /*                                                                 */
    2964             : /*                  GCD USING SUBRESULTANT                         */
    2965             : /*                                                                 */
    2966             : /*******************************************************************/
    2967             : static int inexact(GEN x, int *simple, int *rational);
    2968             : static int
    2969     6944596 : isinexactall(GEN x, int *simple, int *rational)
    2970             : {
    2971     6944596 :   long i, lx = lg(x);
    2972    33458343 :   for (i=2; i<lx; i++)
    2973    26513761 :     if (inexact(gel(x,i), simple, rational)) return 1;
    2974     6944582 :   return 0;
    2975             : }
    2976             : /* return 1 if coeff explosion is not possible */
    2977             : static int
    2978    26513971 : inexact(GEN x, int *simple, int *rational)
    2979             : {
    2980    26513971 :   int junk = 0;
    2981    26513971 :   switch(typ(x))
    2982             :   {
    2983    25088855 :     case t_INT: case t_FRAC: return 0;
    2984             : 
    2985           7 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    2986             : 
    2987             :     case t_INTMOD:
    2988             :     case t_FFELT:
    2989     1383501 :       *rational = 0;
    2990     1383501 :       if (!*simple) *simple = 1;
    2991     1383501 :       return 0;
    2992             : 
    2993             :     case t_COMPLEX:
    2994          77 :       *rational = 0;
    2995         154 :       return inexact(gel(x,1), simple, rational)
    2996          77 :           || inexact(gel(x,2), simple, rational);
    2997             :     case t_QUAD:
    2998           0 :       *rational = *simple = 0;
    2999           0 :       return inexact(gel(x,2), &junk, rational)
    3000           0 :           || inexact(gel(x,3), &junk, rational);
    3001             : 
    3002             :     case t_POLMOD:
    3003       11718 :       *rational = 0;
    3004       11718 :       return isinexactall(gel(x,1), simple, rational);
    3005             :     case t_POL:
    3006       29785 :       *rational = 0;
    3007       29785 :       *simple = -1;
    3008       29785 :       return isinexactall(x, &junk, rational);
    3009             :     case t_RFRAC:
    3010          28 :       *rational = 0;
    3011          28 :       *simple = -1;
    3012          56 :       return inexact(gel(x,1), &junk, rational)
    3013          28 :           || inexact(gel(x,2), &junk, rational);
    3014             :   }
    3015           0 :   *rational = 0;
    3016           0 :   *simple = -1; return 0;
    3017             : }
    3018             : 
    3019             : /* x monomial, y t_POL in the same variable */
    3020             : static GEN
    3021     7446453 : gcdmonome(GEN x, GEN y)
    3022             : {
    3023     7446453 :   pari_sp av = avma;
    3024     7446453 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    3025     7446453 :   long i, l = lg(y);
    3026     7446453 :   GEN t, v = cgetg(l, t_VEC);
    3027     7446453 :   gel(v,1) = gel(x,dx+2);
    3028     7446453 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    3029     7446453 :   t = content(v); /* gcd(lc(x), cont(y)) */
    3030     7446453 :   t = simplify_shallow(t);
    3031     7446453 :   if (dx < e) e = dx;
    3032     7446453 :   return gerepileupto(av, monomialcopy(t, e, varn(x)));
    3033             : }
    3034             : 
    3035             : /* x, y are t_POL in the same variable */
    3036             : GEN
    3037    10898388 : RgX_gcd(GEN x, GEN y)
    3038             : {
    3039             :   long dx, dy;
    3040             :   pari_sp av, av1;
    3041             :   GEN d, g, h, p1, p2, u, v;
    3042    10898388 :   int simple = 0, rational = 1;
    3043             : 
    3044    10898388 :   if (isexactzero(y)) return RgX_copy(x);
    3045    10898010 :   if (isexactzero(x)) return RgX_copy(y);
    3046    10898003 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    3047     4099077 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    3048     3451550 :   if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
    3049             :   {
    3050           7 :     av = avma; u = ggcd(content(x), content(y));
    3051           7 :     return gerepileupto(av, scalarpol(u, varn(x)));
    3052             :   }
    3053     3451543 :   if (rational) return QX_gcd(x,y); /* Q[X] */
    3054             : 
    3055       92099 :   av = avma;
    3056       92099 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    3057             :   else
    3058             :   {
    3059        4935 :     dx = lg(x); dy = lg(y);
    3060        4935 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    3061        4935 :     if (dy==3)
    3062             :     {
    3063           0 :       d = ggcd(gel(y,2), content(x));
    3064           0 :       return gerepileupto(av, scalarpol(d, varn(x)));
    3065             :     }
    3066        4935 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    3067        4935 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    3068        4935 :     d = ggcd(p1,p2);
    3069        4935 :     av1 = avma;
    3070        4935 :     g = h = gen_1;
    3071             :     for(;;)
    3072             :     {
    3073        6377 :       GEN r = RgX_pseudorem(u,v);
    3074        6377 :       long degq, du, dv, dr = lg(r);
    3075             : 
    3076        6377 :       if (!signe(r)) break;
    3077        3584 :       if (dr <= 3)
    3078             :       {
    3079        2142 :         avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
    3080             :       }
    3081        1442 :       if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
    3082        1442 :       du = lg(u); dv = lg(v); degq = du-dv;
    3083        1442 :       u = v; p1 = g; g = leading_coeff(u);
    3084        1442 :       switch(degq)
    3085             :       {
    3086         189 :         case 0: break;
    3087             :         case 1:
    3088        1029 :           p1 = gmul(h,p1); h = g; break;
    3089             :         default:
    3090         224 :           p1 = gmul(gpowgs(h,degq), p1);
    3091         224 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    3092             :       }
    3093        1442 :       v = RgX_Rg_div(r,p1);
    3094        1442 :       if (gc_needed(av1,1))
    3095             :       {
    3096           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
    3097           0 :         gerepileall(av1,4, &u,&v,&g,&h);
    3098             :       }
    3099        1442 :     }
    3100        2793 :     x = RgX_Rg_mul(primpart(v), d);
    3101             :   }
    3102       89957 :   if (must_negate(x)) x = RgX_neg(x);
    3103       89957 :   return gerepileupto(av,x);
    3104             : }
    3105             : 
    3106             : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
    3107             : static GEN
    3108        5145 : RgX_disc_aux(GEN P)
    3109             : {
    3110        5145 :   long n = degpol(P), TP, dd;
    3111             :   GEN D, L, y, p;
    3112        5145 :   if (!signe(P) || !n) return Rg_get_0(P);
    3113        5145 :   if (n == 1) return Rg_get_1(P);
    3114        4970 :   if (n == 2) {
    3115        1134 :     GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
    3116        1134 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    3117             :   }
    3118        3836 :   TP = RgX_simpletype(P);
    3119        3836 :   if (TP == t_INT) return ZX_disc(P);
    3120         483 :   if (TP == t_FRAC) return QX_disc(P);
    3121         483 :   p = NULL;
    3122         483 :   if (RgX_is_FpX(P, &p) && p)
    3123          28 :     return Fp_to_mod(FpX_disc(RgX_to_FpX(P,p), p), p);
    3124             : 
    3125         455 :   y = RgX_deriv(P);
    3126         455 :   if (!signe(y)) return Rg_get_0(y);
    3127         455 :   dd = degpol(P)-2 - degpol(y);
    3128         455 :   if (TP == t_REAL)
    3129          14 :     D = resultant2(P,y);
    3130             :   else
    3131             :   {
    3132         441 :     D = RgX_resultant_all(P, y, NULL);
    3133         441 :     if (D == gen_0) return Rg_get_0(y);
    3134             :   }
    3135         455 :   L = leading_coeff(P);
    3136         455 :   if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
    3137         455 :   if (n & 2) D = gneg(D);
    3138         455 :   return D;
    3139             : }
    3140             : GEN
    3141        1953 : RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }
    3142             : 
    3143             : GEN
    3144        3206 : poldisc0(GEN x, long v)
    3145             : {
    3146             :   pari_sp av;
    3147        3206 :   switch(typ(x))
    3148             :   {
    3149             :     case t_POL:
    3150             :     {
    3151             :       GEN D;
    3152        3192 :       long v0 = -1;
    3153        3192 :       av = avma;
    3154        3192 :       if (v >= 0 && v != varn(x))
    3155             :       {
    3156           0 :         v0 = fetch_var_higher();
    3157           0 :         x = fix_pol(x,v, v0);
    3158             :       }
    3159        3192 :       D = RgX_disc_aux(x);
    3160        3192 :       if (v0 >= 0) (void)delete_var();
    3161        3192 :       return gerepileupto(av, D);
    3162             :     }
    3163             : 
    3164             :     case t_COMPLEX:
    3165           0 :       return utoineg(4);
    3166             : 
    3167             :     case t_QUAD:
    3168           0 :       return quad_disc(x);
    3169             :     case t_POLMOD:
    3170           0 :       return poldisc0(gel(x,1), v);
    3171             : 
    3172             :     case t_QFR: case t_QFI:
    3173          14 :       av = avma; return gerepileuptoint(av, qfb_disc(x));
    3174             : 
    3175             :     case t_VEC: case t_COL: case t_MAT:
    3176             :     {
    3177             :       long i;
    3178           0 :       GEN z = cgetg_copy(x, &i);
    3179           0 :       for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
    3180           0 :       return z;
    3181             :     }
    3182             :   }
    3183           0 :   pari_err_TYPE("poldisc",x);
    3184             :   return NULL; /* LCOV_EXCL_LINE */
    3185             : }
    3186             : 
    3187             : GEN
    3188           7 : reduceddiscsmith(GEN x)
    3189             : {
    3190           7 :   long j, n = degpol(x);
    3191           7 :   pari_sp av = avma;
    3192             :   GEN xp, M;
    3193             : 
    3194           7 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    3195           7 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    3196           7 :   RgX_check_ZX(x,"poldiscreduced");
    3197           7 :   if (!gequal1(gel(x,n+2)))
    3198           0 :     pari_err_IMPL("non-monic polynomial in poldiscreduced");
    3199           7 :   M = cgetg(n+1,t_MAT);
    3200           7 :   xp = ZX_deriv(x);
    3201          28 :   for (j=1; j<=n; j++)
    3202             :   {
    3203          21 :     gel(M,j) = RgX_to_RgC(xp, n);
    3204          21 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    3205             :   }
    3206           7 :   return gerepileupto(av, ZM_snf(M));
    3207             : }
    3208             : 
    3209             : /***********************************************************************/
    3210             : /**                                                                   **/
    3211             : /**                       STURM ALGORITHM                             **/
    3212             : /**              (number of real roots of x in [a,b])                 **/
    3213             : /**                                                                   **/
    3214             : /***********************************************************************/
    3215             : static GEN
    3216         833 : R_to_Q_up(GEN x)
    3217             : {
    3218             :   long e;
    3219         833 :   switch(typ(x))
    3220             :   {
    3221         833 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3222             :     case t_REAL:
    3223           0 :       x = mantissa_real(x,&e);
    3224           0 :       return gmul2n(addiu(x,1), -e);
    3225           0 :     default: pari_err_TYPE("R_to_Q_up", x);
    3226             :       return NULL; /* LCOV_EXCL_LINE */
    3227             :   }
    3228             : }
    3229             : static GEN
    3230         833 : R_to_Q_down(GEN x)
    3231             : {
    3232             :   long e;
    3233         833 :   switch(typ(x))
    3234             :   {
    3235         819 :     case t_INT: case t_FRAC: case t_INFINITY: return x;
    3236             :     case t_REAL:
    3237          14 :       x = mantissa_real(x,&e);
    3238          14 :       return gmul2n(subiu(x,1), -e);
    3239           0 :     default: pari_err_TYPE("R_to_Q_down", x);
    3240             :       return NULL; /* LCOV_EXCL_LINE */
    3241             :   }
    3242             : }
    3243             : 
    3244             : static long
    3245         847 : sturmpart_i(GEN x, GEN ab)
    3246             : {
    3247         847 :   long tx = typ(x);
    3248         847 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    3249         847 :   if (tx != t_POL)
    3250             :   {
    3251           0 :     if (is_real_t(tx)) return 0;
    3252           0 :     pari_err_TYPE("sturm",x);
    3253             :   }
    3254         847 :   if (lg(x) == 3) return 0;
    3255         847 :   if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
    3256         847 :   if (!ZX_is_squarefree(x))
    3257          14 :     pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,x);
    3258         833 :   if (ab)
    3259             :   {
    3260             :     GEN A, B;
    3261         833 :     if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
    3262         833 :     A = R_to_Q_down(gel(ab,1));
    3263         833 :     B = R_to_Q_up(gel(ab,2));
    3264         833 :     ab = mkvec2(A, B);
    3265             :   }
    3266         833 :   return ZX_sturmpart(x, ab);
    3267             : }
    3268             : /* Deprecated: RgX_sturmpart() should be preferred  */
    3269             : long
    3270         847 : sturmpart(GEN x, GEN a, GEN b)
    3271             : {
    3272         847 :   pari_sp av = avma;
    3273             :   long r;
    3274         847 :   if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
    3275         714 :   if (!a) a = mkmoo();
    3276         714 :   if (!b) b = mkoo();
    3277         714 :   r = sturmpart_i(x, mkvec2(a, b));
    3278         700 :   avma = av; return r;
    3279             : }
    3280             : long
    3281         133 : RgX_sturmpart(GEN x, GEN ab)
    3282             : {
    3283         133 :   pari_sp av = avma;
    3284         133 :   long r = sturmpart_i(x, ab);
    3285         133 :   avma = av; return r;
    3286             : }
    3287             : 
    3288             : /***********************************************************************/
    3289             : /**                                                                   **/
    3290             : /**                        GENERIC EXTENDED GCD                       **/
    3291             : /**                                                                   **/
    3292             : /***********************************************************************/
    3293             : /* assume typ(x) = typ(y) = t_POL */
    3294             : static GEN
    3295         874 : RgXQ_inv_i(GEN x, GEN y)
    3296             : {
    3297         874 :   long vx=varn(x), vy=varn(y);
    3298             :   pari_sp av;
    3299             :   GEN u, v, d;
    3300             : 
    3301        1748 :   while (vx != vy)
    3302             :   {
    3303           0 :     if (varncmp(vx,vy) > 0)
    3304             :     {
    3305           0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3306           0 :       return scalarpol(d, vy);
    3307             :     }
    3308           0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3309           0 :     x = gel(x,2); vx = gvar(x);
    3310             :   }
    3311         874 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3312         874 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3313         874 :   d = gdiv(u,d);
    3314         874 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3315         874 :   return gerepileupto(av, d);
    3316             : }
    3317             : 
    3318             : /*Assume x is a polynomial and y is not */
    3319             : static GEN
    3320         112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3321             : {
    3322         112 :   long vx = varn(x);
    3323         112 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    3324         112 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    3325          84 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    3326          56 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    3327             : }
    3328             : /* Assume x==0, y!=0 */
    3329             : static GEN
    3330          63 : zero_bezout(GEN y, GEN *U, GEN *V)
    3331             : {
    3332          63 :   *U=gen_0; *V = ginv(y); return gen_1;
    3333             : }
    3334             : 
    3335             : GEN
    3336         280 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    3337             : {
    3338         280 :   long tx=typ(x), ty=typ(y), vx;
    3339         280 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    3340         245 :   if (tx != t_POL)
    3341             :   {
    3342         140 :     if (ty == t_POL)
    3343          56 :       return scalar_bezout(y,x,v,u);
    3344             :     else
    3345             :     {
    3346          84 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    3347          84 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    3348          63 :       if (xis0) return zero_bezout(y,u,v);
    3349          42 :       else return zero_bezout(x,v,u);
    3350             :     }
    3351             :   }
    3352         105 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    3353          49 :   vx = varn(x);
    3354          49 :   if (vx != varn(y))
    3355           0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    3356           0 :                                    : scalar_bezout(y,x,v,u);
    3357          49 :   return RgX_extgcd(x,y,u,v);
    3358             : }
    3359             : 
    3360             : GEN
    3361         280 : gcdext0(GEN x, GEN y)
    3362             : {
    3363         280 :   GEN z=cgetg(4,t_VEC);
    3364         280 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    3365         280 :   return z;
    3366             : }
    3367             : 
    3368             : /*******************************************************************/
    3369             : /*                                                                 */
    3370             : /*                    GENERIC (modular) INVERSE                    */
    3371             : /*                                                                 */
    3372             : /*******************************************************************/
    3373             : 
    3374             : GEN
    3375       14904 : ginvmod(GEN x, GEN y)
    3376             : {
    3377       14904 :   long tx=typ(x);
    3378             : 
    3379       14904 :   switch(typ(y))
    3380             :   {
    3381             :     case t_POL:
    3382       14904 :       if (tx==t_POL) return RgXQ_inv(x,y);
    3383        2023 :       if (is_scalar_t(tx)) return ginv(x);
    3384           0 :       break;
    3385             :     case t_INT:
    3386           0 :       if (tx==t_INT) return Fp_inv(x,y);
    3387           0 :       if (tx==t_POL) return gen_0;
    3388             :   }
    3389           0 :   pari_err_TYPE2("ginvmod",x,y);
    3390             :   return NULL; /* LCOV_EXCL_LINE */
    3391             : }
    3392             : 
    3393             : /***********************************************************************/
    3394             : /**                                                                   **/
    3395             : /**                         NEWTON POLYGON                            **/
    3396             : /**                                                                   **/
    3397             : /***********************************************************************/
    3398             : 
    3399             : /* assume leading coeff of x is non-zero */
    3400             : GEN
    3401          28 : newtonpoly(GEN x, GEN p)
    3402             : {
    3403             :   GEN y;
    3404             :   long n,ind,a,b,c,u1,u2,r1,r2;
    3405          28 :   long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
    3406             : 
    3407          28 :   if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
    3408          28 :   n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
    3409          28 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    3410          28 :   vval = (long *) pari_malloc(sizeof(long)*(n+1));
    3411          28 :   for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
    3412          42 :   for (a=0, ind=1; a<n; a++)
    3413             :   {
    3414          42 :     if (vval[a] != LONG_MAX) break;
    3415          14 :     gel(y,ind++) = mkoo();
    3416             :   }
    3417          84 :   for (b=a+1; b<=n; a=b, b=a+1)
    3418             :   {
    3419          56 :     while (vval[b] == LONG_MAX) b++;
    3420          56 :     u1 = vval[a]-vval[b];
    3421          56 :     u2 = b-a;
    3422         154 :     for (c=b+1; c<=n; c++)
    3423             :     {
    3424          98 :       if (vval[c] == LONG_MAX) continue;
    3425          70 :       r1 = vval[a]-vval[c];
    3426          70 :       r2 = c-a;
    3427          70 :       if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
    3428             :     }
    3429          56 :     while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
    3430             :   }
    3431          28 :   pari_free(vval); return y;
    3432             : }
    3433             : 
    3434             : static GEN
    3435      262591 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
    3436             : {
    3437      262591 :   pari_sp av = avma;
    3438             :   GEN r;
    3439      262591 :   if (lgefint(p) == 3)
    3440             :   {
    3441      144256 :     ulong pp = uel(p, 2);
    3442      144256 :     r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
    3443             :                 RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
    3444             :   }
    3445             :   else
    3446      118335 :     r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
    3447      262591 :   return gerepileupto(av, FpX_to_mod(r, p));
    3448             : }
    3449             : 
    3450             : static GEN
    3451          14 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
    3452             : {
    3453          14 :   pari_sp av = avma;
    3454             :   GEN r;
    3455          14 :   if (lgefint(p) == 3)
    3456             :   {
    3457           7 :     ulong pp = uel(p, 2);
    3458           7 :     r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
    3459             :                                    RgX_to_Flx(y, pp), pp));
    3460             :   }
    3461             :   else
    3462           7 :     r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3463          14 :   return gerepileupto(av, FpX_to_mod(r, p));
    3464             : }
    3465             : 
    3466             : static GEN
    3467       11571 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
    3468             : {
    3469       11571 :   pari_sp av = avma;
    3470             :   GEN r;
    3471       11571 :   if (lgefint(p) == 3)
    3472             :   {
    3473        5789 :     ulong pp = uel(p, 2);
    3474        5789 :     r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
    3475             :                                    RgX_to_Flx(y, pp), pp));
    3476             :   }
    3477             :   else
    3478        5782 :     r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3479       11571 :   return gerepileupto(av, FpX_to_mod(r, p));
    3480             : }
    3481             : 
    3482             : static GEN
    3483           7 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
    3484             : {
    3485           7 :   pari_sp av = avma;
    3486             :   GEN r;
    3487           7 :   GEN T = RgX_to_FpX(pol, p);
    3488           7 :   if (lgefint(p) == 3)
    3489             :   {
    3490           7 :     ulong pp = uel(p, 2);
    3491           7 :     GEN Tp = ZX_to_Flx(T, pp);
    3492           7 :     r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
    3493             :                                RgX_to_FlxqX(y, Tp, pp),
    3494             :                                RgX_to_FlxqX(S, Tp, pp), Tp, pp));
    3495             :   }
    3496             :   else
    3497           0 :     r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
    3498             :                    RgX_to_FpXQX(S, T, p), T, p);
    3499           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3500             : }
    3501             : 
    3502             : static GEN
    3503           0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    3504             : {
    3505           0 :   pari_sp av = avma;
    3506             :   GEN r;
    3507           0 :   GEN T = RgX_to_FpX(pol, p);
    3508           0 :   if (lgefint(p) == 3)
    3509             :   {
    3510           0 :     ulong pp = uel(p, 2);
    3511           0 :     GEN Tp = ZX_to_Flx(T, pp);
    3512           0 :     r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
    3513             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3514             :   }
    3515             :   else
    3516           0 :     r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3517           0 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3518             : }
    3519             : 
    3520             : static GEN
    3521           7 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
    3522             : {
    3523           7 :   pari_sp av = avma;
    3524             :   GEN r;
    3525           7 :   GEN T = RgX_to_FpX(pol, p);
    3526           7 :   if (lgefint(p) == 3)
    3527             :   {
    3528           7 :     ulong pp = uel(p, 2);
    3529           7 :     GEN Tp = ZX_to_Flx(T, pp);
    3530           7 :     r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
    3531             :                                RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3532             :   }
    3533             :   else
    3534           0 :     r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3535           7 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3536             : }
    3537             : 
    3538             : #define code(t1,t2) ((t1 << 6) | t2)
    3539             : static GEN
    3540      718387 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
    3541             : {
    3542             :   GEN p, pol;
    3543             :   long pa;
    3544      718387 :   long t = RgX_type3(x,y,T, &p,&pol,&pa);
    3545      718389 :   switch(t)
    3546             :   {
    3547      402378 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
    3548       42530 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
    3549          14 :     case t_FFELT:  return FFXQ_mul(x, y, T, pol);
    3550      262591 :     case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
    3551             :     case code(t_POLMOD, t_INTMOD):
    3552           7 :                    return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
    3553       10869 :     default:       return NULL;
    3554             :   }
    3555             : }
    3556             : 
    3557             : GEN
    3558      718383 : RgXQ_mul(GEN x, GEN y, GEN T)
    3559             : {
    3560      718383 :   GEN z = RgXQ_mul_fast(x, y, T);
    3561      718389 :   if (!z) z = RgX_rem(RgX_mul(x, y), T);
    3562      718389 :   return z;
    3563             : }
    3564             : 
    3565             : static GEN
    3566      277409 : RgXQ_sqr_fast(GEN x, GEN T)
    3567             : {
    3568             :   GEN p, pol;
    3569             :   long pa;
    3570      277409 :   long t = RgX_type2(x, T, &p,&pol,&pa);
    3571      277412 :   switch(t)
    3572             :   {
    3573      265379 :     case t_INT:    return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
    3574       10843 :     case t_FRAC:   return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
    3575           0 :     case t_FFELT:  return FFXQ_sqr(x, T, pol);
    3576          14 :     case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
    3577             :     case code(t_POLMOD, t_INTMOD):
    3578           0 :                    return RgXQ_sqr_FpXQXQ(x, T, pol, p);
    3579        1176 :     default:       return NULL;
    3580             :   }
    3581             : }
    3582             : 
    3583             : GEN
    3584      277411 : RgXQ_sqr(GEN x, GEN T)
    3585             : {
    3586      277411 :   GEN z = RgXQ_sqr_fast(x, T);
    3587      277400 :   if (!z) z = RgX_rem(RgX_sqr(x), T);
    3588      277400 :   return z;
    3589             : }
    3590             : 
    3591             : static GEN
    3592       29269 : RgXQ_inv_fast(GEN x, GEN y)
    3593             : {
    3594             :   GEN p, pol;
    3595             :   long pa;
    3596       29269 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3597       29269 :   switch(t)
    3598             :   {
    3599       16061 :     case t_INT:    return QXQ_inv(x,y);
    3600         749 :     case t_FRAC:   return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
    3601          14 :     case t_FFELT:  return FFXQ_inv(x, y, pol);
    3602       11571 :     case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
    3603             :     case code(t_POLMOD, t_INTMOD):
    3604           7 :                    return RgXQ_inv_FpXQXQ(x, y, pol, p);
    3605         867 :     default:       return NULL;
    3606             :   }
    3607             : }
    3608             : #undef code
    3609             : 
    3610             : GEN
    3611       29269 : RgXQ_inv(GEN x, GEN y)
    3612             : {
    3613       29269 :   GEN z = RgXQ_inv_fast(x, y);
    3614       29255 :   if (!z) z = RgXQ_inv_i(x, y);
    3615       29255 :   return z;
    3616             : }

Generated by: LCOV version 1.11