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.8.0 lcov report (development 19369-efd6c3d) Lines: 1494 1738 86.0 %
Date: 2016-08-29 06:11:50 Functions: 122 133 91.7 %
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             : #define addshift(x,y) addshiftpol((x),(y),1)
      24             : 
      25             : /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
      26             :  * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
      27             :  * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
      28             :  * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
      29             :  * Not memory clean in the latter case */
      30             : GEN
      31       23611 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      32             : {
      33       23611 :   long dP=degpol(P), i, k, m;
      34             :   pari_sp av1, av2;
      35             :   GEN s,y,P_lead;
      36             : 
      37       23611 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      38       23611 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      39       23611 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      40       23611 :   y = cgetg(n+2,t_COL);
      41       23611 :   if (y0)
      42             :   {
      43       10675 :     if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
      44       10675 :     m = lg(y0)-1;
      45       10675 :     for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
      46             :   }
      47             :   else
      48             :   {
      49       12936 :     m = 1;
      50       12936 :     gel(y,1) = stoi(dP);
      51             :   }
      52       23611 :   P += 2; /* strip codewords */
      53             : 
      54       23611 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      55       23611 :   if (P_lead)
      56             :   {
      57           7 :     if (N) P_lead = Fq_inv(P_lead,T,N);
      58           7 :     else if (T) P_lead = QXQ_inv(P_lead,T);
      59             :   }
      60       61705 :   for (k=m; k<=n; k++)
      61             :   {
      62       38094 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      63      186277 :     for (i=1; i<k && i<=dP; i++)
      64      148183 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      65       38094 :     if (N)
      66             :     {
      67       15764 :       s = Fq_red(s, T, N);
      68       15764 :       if (P_lead) s = Fq_mul(s, P_lead, T, N);
      69             :     }
      70       22330 :     else if (T)
      71             :     {
      72        2044 :       s = grem(s, T);
      73        2044 :       if (P_lead) s = grem(gmul(s, P_lead), T);
      74             :     }
      75             :     else
      76       20286 :       if (P_lead) s = gdiv(s, P_lead);
      77       38094 :     av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
      78             :   }
      79       23611 :   return y;
      80             : }
      81             : 
      82             : GEN
      83        8120 : polsym(GEN x, long n)
      84             : {
      85        8120 :   return polsym_gen(x, NULL, n, NULL,NULL);
      86             : }
      87             : 
      88             : /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
      89             : GEN
      90   117425691 : centermodii(GEN x, GEN p, GEN po2)
      91             : {
      92   117425691 :   GEN y = remii(x, p);
      93   117415220 :   switch(signe(y))
      94             :   {
      95    25963648 :     case 0: break;
      96    60606912 :     case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
      97    60591165 :       break;
      98    30886384 :     case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
      99    30872774 :       break;
     100             :   }
     101   117385863 :   return y;
     102             : }
     103             : 
     104             : static long
     105           0 : s_centermod(long x, ulong pp, ulong pps2)
     106             : {
     107           0 :   long y = x % (long)pp;
     108           0 :   if (y < 0) y += pp;
     109           0 :   return Fl_center(y, pp,pps2);
     110             : }
     111             : 
     112             : /* for internal use */
     113             : GEN
     114     7386038 : centermod_i(GEN x, GEN p, GEN ps2)
     115             : {
     116             :   long i, lx;
     117             :   pari_sp av;
     118             :   GEN y;
     119             : 
     120     7386038 :   if (!ps2) ps2 = shifti(p,-1);
     121     7386034 :   switch(typ(x))
     122             :   {
     123     3663037 :     case t_INT: return centermodii(x,p,ps2);
     124             : 
     125     2875013 :     case t_POL: lx = lg(x);
     126     2875013 :       y = cgetg(lx,t_POL); y[1] = x[1];
     127    21915262 :       for (i=2; i<lx; i++)
     128             :       {
     129    19039979 :         av = avma;
     130    19039979 :         gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
     131             :       }
     132     2875283 :       return normalizepol_lg(y, lx);
     133             : 
     134      847144 :     case t_COL: lx = lg(x);
     135      847144 :       y = cgetg(lx,t_COL);
     136      847144 :       for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
     137      847144 :       return y;
     138             : 
     139         840 :     case t_MAT: lx = lg(x);
     140         840 :       y = cgetg(lx,t_MAT);
     141         840 :       for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
     142         840 :       return y;
     143             : 
     144           0 :     case t_VECSMALL: lx = lg(x);
     145             :     {
     146           0 :       ulong pp = itou(p), pps2 = itou(ps2);
     147           0 :       y = cgetg(lx,t_VECSMALL);
     148           0 :       for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
     149           0 :       return y;
     150             :     }
     151             :   }
     152           0 :   return x;
     153             : }
     154             : 
     155             : GEN
     156     4832226 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
     157             : 
     158             : /***********************************************************************/
     159             : /**                                                                   **/
     160             : /**                          FACTORIZATION                            **/
     161             : /**                                                                   **/
     162             : /***********************************************************************/
     163             : #define assign_or_fail(x,y) { GEN __x = x;\
     164             :   if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
     165             : }
     166             : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
     167             : 
     168             : static const long tsh = 6;
     169             : static long
     170         175 : code(long t1, long t2) { return (t1 << tsh) | t2; }
     171             : void
     172         140 : RgX_type_decode(long x, long *t1, long *t2)
     173             : {
     174         140 :   *t1 = x >> tsh;
     175         140 :   *t2 = (x & ((1L<<tsh)-1));
     176         140 : }
     177             : int
     178      281643 : RgX_type_is_composite(long t) { return t >= tsh; }
     179             : 
     180             : long
     181      387252 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     182             : {
     183      387252 :   long t[] = {0,0,0,0,0,0,0,0,0,0};
     184      387252 :   long tx = typ(x), lx, i, j, t2 = 0;
     185      387252 :   GEN ff = NULL;
     186      387252 :   *p = *pol = NULL; *pa = LONG_MAX;
     187      387252 :   if (is_scalar_t(tx))
     188             :   {
     189           0 :     if (tx == t_POLMOD) return 0;
     190           0 :     x = scalarpol(x,0);
     191             :   }
     192             :   /* t[0..1] unused. Other values, if set, indicate a coefficient of type
     193             :    * t[2] : t_REAL
     194             :    * t[3] : t_INTMOD
     195             :    * t[4] : t_COMPLEX of rationals (t_INT/t_FRAC)
     196             :    * t[5] : t_FFELT
     197             :    * t[6] : t_COMPLEX of t_REAL
     198             :    * t[7] : t_PADIC
     199             :    * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
     200             :    * t[9]: t_POLMOD of rationals (t_INT/t_FRAC) */
     201             :   /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
     202             :    * given by t) */
     203      387252 :   lx = lg(x);
     204     2243735 :   for (i=2; i<lx; i++)
     205             :   {
     206     1863210 :     GEN c = gel(x,i);
     207     1863210 :     switch(typ(c))
     208             :     {
     209             :       case t_INT: case t_FRAC:
     210     1391550 :         break;
     211             :       case t_REAL:
     212        1946 :         update_prec(precision(c), pa);
     213        1946 :         t[2]=1; break;
     214             :       case t_INTMOD:
     215         189 :         assign_or_fail(gel(c,1),p);
     216         189 :         t[3]=1; break;
     217             :       case t_FFELT:
     218      461832 :         if (!ff) ff=c; else if (!FF_samefield(c,ff)) return 0;
     219      461832 :         assign_or_fail(FF_p_i(c),p);
     220      461832 :         t[5]=1; break;
     221             :       case t_COMPLEX:
     222          91 :         for (j=1; j<=2; j++)
     223             :         {
     224          63 :           GEN d = gel(c,j);
     225          63 :           switch(typ(d))
     226             :           {
     227             :             case t_INT: case t_FRAC:
     228          35 :               t[4]=1; break;
     229             :             case t_REAL:
     230           7 :               update_prec(precision(d), pa);
     231           7 :               t[6]=1; break;
     232             :             case t_INTMOD:
     233          14 :               assign_or_fail(gel(d,1),p);
     234          14 :               if (!signe(*p) || mod4(*p) != 3) return 0;
     235           7 :               if (!t2) t2 = t_COMPLEX;
     236           7 :               t[3]=1; break;
     237             :             case t_PADIC:
     238           7 :               update_prec(precp(d)+valp(d), pa);
     239           7 :               assign_or_fail(gel(d,2),p);
     240           7 :               if (!t2) t2 = t_COMPLEX;
     241           7 :               t[7]=1; break;
     242           0 :             default: return 0;
     243             :           }
     244             :         }
     245          28 :         if (!t[6]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
     246          28 :         break;
     247             :       case t_PADIC:
     248          70 :         update_prec(precp(c)+valp(c), pa);
     249          70 :         assign_or_fail(gel(c,2),p);
     250          70 :         t[7]=1; break;
     251             :       case t_QUAD:
     252          28 :         assign_or_fail(gel(c,1),pol);
     253          84 :         for (j=2; j<=3; j++)
     254             :         {
     255          56 :           GEN d = gel(c,j);
     256          56 :           switch(typ(d))
     257             :           {
     258             :             case t_INT: case t_FRAC:
     259          21 :               t[8]=1; break;
     260             :             case t_INTMOD:
     261          28 :               assign_or_fail(gel(d,1),p);
     262          28 :               if (t2 != t_POLMOD) t2 = t_QUAD;
     263          28 :               t[3]=1; break;
     264             :             case t_PADIC:
     265           7 :               update_prec(precp(d)+valp(d), pa);
     266           7 :               assign_or_fail(gel(d,2),p);
     267           7 :               if (t2 != t_POLMOD) t2 = t_QUAD;
     268           7 :               t[7]=1; break;
     269           0 :             default: return 0;
     270             :           }
     271             :         }
     272          28 :         break;
     273             :       case t_POLMOD:
     274         840 :         assign_or_fail(gel(c,1),pol);
     275        5040 :         for (j=1; j<=2; j++)
     276             :         {
     277        1680 :           GEN pbis = NULL, polbis = NULL;
     278             :           long pabis;
     279        1680 :           switch(RgX_type(gel(c,j),&pbis,&polbis,&pabis))
     280             :           {
     281        1659 :             case t_INT: t[9]=1; break;
     282          14 :             case t_INTMOD: t[3]=1; t2 = t_POLMOD; break;
     283           7 :             case t_PADIC: t[7]=1; t2 = t_POLMOD; update_prec(pabis,pa); break;
     284           0 :             default: return 0;
     285             :           }
     286        1680 :           if (pbis) assign_or_fail(pbis,p);
     287        1680 :           if (polbis) assign_or_fail(polbis,pol);
     288             :         }
     289         840 :         break;
     290        6720 :       default: return 0;
     291             :     }
     292             :   }
     293      380525 :   if (t[5]) /* ffelt */
     294             :   {
     295      103453 :     if (t2 ||t[2]||t[4]||t[6]||t[8]||t[9]) return 0;
     296      103453 :     *pol=ff; return t_FFELT;
     297             :   }
     298      277072 :   if (t[6]) /* inexact, complex */
     299             :   {
     300           7 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
     301           7 :     return t_COMPLEX;
     302             :   }
     303      277065 :   if (t[2]) /* inexact, real */
     304             :   {
     305         280 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
     306         280 :     return t[4]?t_COMPLEX:t_REAL;
     307             :   }
     308      276785 :   if (t2) /* polmod/quad/complex of intmod/padic */
     309             :   {
     310          56 :     if (t[3]) return code(t2,t_INTMOD);
     311          21 :     if (t[7]) return code(t2,t_PADIC);
     312             :   }
     313      276729 :   if (t[9]) return code(t_POLMOD,t_INT);
     314      276624 :   if (t[8]) return code(t_QUAD,t_INT);
     315      276617 :   if (t[4]) return code(t_COMPLEX,t_INT);
     316      276610 :   if (t[3]) return t_INTMOD;
     317      276533 :   if (t[7]) return t_PADIC;
     318      276505 :   return t_INT;
     319             : }
     320             : 
     321             : GEN
     322           0 : factor0(GEN x,long flag)
     323             : {
     324           0 :   return (flag<0)? factor(x): boundfact(x,flag);
     325             : }
     326             : 
     327             : /* only present for interface with GP */
     328             : GEN
     329       36951 : gp_factor0(GEN x, GEN flag)
     330             : {
     331             :   ulong B;
     332       36951 :   if (!flag) return factor(x);
     333          35 :   if (typ(flag) != t_INT || signe(flag) < 0) pari_err_FLAG("factor");
     334          35 :   switch(lgefint(flag))
     335             :   {
     336           7 :     case 2: B = 0; break;
     337          28 :     case 3: B = flag[2]; break;
     338           0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     339           0 :              return NULL; /*not reached*/
     340             :   }
     341          35 :   return boundfact(x, B);
     342             : }
     343             : 
     344             : GEN
     345       39263 : deg1_from_roots(GEN L, long v)
     346             : {
     347       39263 :   long i, l = lg(L);
     348       39263 :   GEN z = cgetg(l,t_COL);
     349       79366 :   for (i=1; i<l; i++)
     350       40103 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     351       39263 :   return z;
     352             : }
     353             : GEN
     354         959 : roots_from_deg1(GEN x)
     355             : {
     356         959 :   long i,l = lg(x);
     357         959 :   GEN r = cgetg(l,t_VEC);
     358         959 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     359         959 :   return r;
     360             : }
     361             : 
     362             : static GEN
     363          35 : gauss_factor_p(GEN p)
     364             : {
     365          35 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     366          35 :   return mkcomplex(a, b);
     367             : }
     368             : 
     369             : static GEN
     370          42 : gauss_primpart(GEN x, GEN *c)
     371             : {
     372          42 :   GEN a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
     373          42 :   *c = n; if (n == gen_1) return x;
     374          42 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     375             : }
     376             : 
     377             : static GEN
     378          70 : gauss_primpart_try(GEN x, GEN c)
     379             : {
     380             :   GEN r, y;
     381          70 :   if (typ(x) == t_INT)
     382             :   {
     383          42 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     384             :   }
     385             :   else
     386             :   {
     387          28 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     388          28 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     389          14 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     390             :   }
     391          56 :   return y;
     392             : }
     393             : 
     394             : static int
     395          77 : gauss_cmp(GEN x, GEN y)
     396             : {
     397             :   int v;
     398          77 :   if (typ(x) != t_COMPLEX)
     399           0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     400          77 :   if (typ(y) != t_COMPLEX) return 1;
     401          49 :   v = cmpii(gel(x,2), gel(y,2));
     402          49 :   if (v) return v;
     403          21 :   return gcmp(gel(x,1), gel(y,1));
     404             : }
     405             : 
     406             : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
     407             : static GEN
     408         112 : gauss_normal(GEN x)
     409             : {
     410         112 :   if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
     411         105 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     412         105 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     413         105 :   return x;
     414             : }
     415             : 
     416             : static GEN
     417          42 : gauss_factor(GEN x)
     418             : {
     419          42 :   pari_sp av = avma;
     420          42 :   GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
     421          42 :   long t1 = typ(a);
     422          42 :   long t2 = typ(b), i, j, l, exp = 0;
     423          42 :   if (t1 == t_FRAC) d = gel(a,2);
     424          42 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     425          42 :   if (d == gen_1) y = x;
     426             :   else
     427             :   {
     428          14 :     y = gmul(x, d);
     429          14 :     a = gel(y,1); t1 = typ(a);
     430          14 :     b = gel(y,2); t2 = typ(b);
     431             :   }
     432          42 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     433          42 :   y = gauss_primpart(y, &n);
     434          42 :   fa = factor(cxnorm(y));
     435          42 :   P = gel(fa,1);
     436          42 :   E = gel(fa,2); l = lg(P);
     437          42 :   P2 = cgetg(l, t_COL);
     438          42 :   E2 = cgetg(l, t_COL);
     439          98 :   for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
     440             :   { /* either p = 2 (ramified) or those factors split in Q(i) */
     441          56 :     GEN p = gel(P,i), w, w2, t, we, pe;
     442          56 :     long v, e = itos(gel(E,i));
     443          56 :     int is2 = absequaliu(p, 2);
     444          56 :     w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     445          56 :     w2 = gauss_normal( gconj(w) );
     446             :     /* w * w2 * I^3 = p, w2 = gconj(w) * I */
     447          56 :     pe = powiu(p, e);
     448          56 :     we = gpowgs(w, e);
     449          56 :     t = gauss_primpart_try( gmul(y, gconj(we)), pe );
     450          56 :     if (t) y = t; /* y /= w^e */
     451             :     else {
     452             :       /* y /= conj(w)^e, should be y /= w2^e */
     453          14 :       y = gauss_primpart_try( gmul(y, we), pe );
     454          14 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     455             :     }
     456          56 :     gel(P,i) = w;
     457          56 :     v = Z_pvalrem(n, p, &n);
     458          56 :     if (v) {
     459           7 :       exp -= v; /* += 3*v mod 4 */
     460           7 :       if (is2) v <<= 1; /* 2 = w^2 I^3 */
     461             :       else {
     462           0 :         gel(P2,j) = w2;
     463           0 :         gel(E2,j) = utoipos(v); j++;
     464             :       }
     465           7 :       gel(E,i) = stoi(e + v);
     466             :     }
     467          56 :     v = Z_pvalrem(d, p, &d);
     468          56 :     if (v) {
     469           7 :       exp += v; /* -= 3*v mod 4 */
     470           7 :       if (is2) v <<= 1; /* 2 is ramified */
     471             :       else {
     472           7 :         gel(P2,j) = w2;
     473           7 :         gel(E2,j) = utoineg(v); j++;
     474             :       }
     475           7 :       gel(E,i) = stoi(e - v);
     476             :     }
     477          56 :     exp &= 3;
     478             :   }
     479          42 :   if (j > 1) {
     480           7 :     long k = 1;
     481           7 :     GEN P1 = cgetg(l, t_COL);
     482           7 :     GEN E1 = cgetg(l, t_COL);
     483             :     /* remove factors with exponent 0 */
     484          14 :     for (i = 1; i < l; i++)
     485           7 :       if (signe(gel(E,i)))
     486             :       {
     487           0 :         gel(P1,k) = gel(P,i);
     488           0 :         gel(E1,k) = gel(E,i);
     489           0 :         k++;
     490             :       }
     491           7 :     setlg(P1, k); setlg(E1, k);
     492           7 :     setlg(P2, j); setlg(E2, j);
     493           7 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     494             :   }
     495          42 :   if (!is_pm1(n) || !is_pm1(d))
     496             :   {
     497          21 :     GEN Fa = factor(gdiv(n, d));
     498          21 :     P = gel(Fa,1); l = lg(P);
     499          21 :     E = gel(Fa,2);
     500          49 :     for (i = 1; i < l; i++)
     501             :     {
     502          28 :       GEN w, p = gel(P,i);
     503             :       long e;
     504             :       int is2;
     505          28 :       switch(mod4(p))
     506             :       {
     507          14 :         case 3: continue;
     508           7 :         case 2: is2 = 1; break;
     509           7 :         default:is2 = 0; break;
     510             :       }
     511          14 :       e = itos(gel(E,i));
     512          14 :       w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     513          14 :       gel(P,i) = w;
     514          14 :       if (is2)
     515           7 :         gel(E,i) = stoi(2*e);
     516             :       else
     517             :       {
     518           7 :         P = shallowconcat(P, gauss_normal( gconj(w) ));
     519           7 :         E = shallowconcat(E, gel(E,i));
     520             :       }
     521          14 :       exp -= e; /* += 3*e mod 4 */
     522          14 :       exp &= 3;
     523             :     }
     524          21 :     gel(Fa,1) = P;
     525          21 :     gel(Fa,2) = E;
     526          21 :     fa = famat_mul_shallow(fa, Fa);
     527             :   }
     528          42 :   fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
     529             : 
     530          42 :   y = gmul(y, powIs(exp));
     531          42 :   if (!gequal1(y)) {
     532          35 :     gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
     533          35 :     gel(fa,2) = shallowconcat(gen_1,    gel(fa,2));
     534             :   }
     535          42 :   return gerepilecopy(av, fa);
     536             : }
     537             : 
     538             : GEN
     539       39946 : factor(GEN x)
     540             : {
     541       39946 :   long tx=typ(x), lx, i, pa, v, r1;
     542             :   pari_sp av, tetpil;
     543             :   GEN  y, p, p1, p2, pol;
     544             : 
     545       39946 :   if (gequal0(x))
     546          56 :     switch(tx)
     547             :     {
     548             :       case t_INT:
     549             :       case t_COMPLEX:
     550             :       case t_POL:
     551          56 :       case t_RFRAC: return prime_fact(x);
     552           0 :       default: pari_err_TYPE("factor",x);
     553             :     }
     554       39891 :   av = avma;
     555       39891 :   switch(tx)
     556             :   {
     557       39072 :     case t_INT: return Z_factor(x);
     558             : 
     559             :     case t_FRAC:
     560         308 :       p1 = Z_factor(gel(x,1));
     561         308 :       p2 = Z_factor(gel(x,2)); gel(p2,2) = gneg_i(gel(p2,2));
     562         308 :       return gerepilecopy(av, merge_factor_i(p1,p2));
     563             : 
     564             :     case t_POL:
     565         462 :       tx=RgX_type(x,&p,&pol,&pa);
     566         462 :       switch(tx)
     567             :       {
     568           7 :         case 0: pari_err_IMPL("factor for general polynomials");
     569         252 :         case t_INT: return QX_factor(x);
     570           7 :         case t_INTMOD: return factmod(x,p);
     571             : 
     572           7 :         case t_COMPLEX: y=cgetg(3,t_MAT); lx=lg(x)-2;
     573           7 :           av = avma; p1 = roots(x,pa); tetpil = avma;
     574           7 :           p1 = deg1_from_roots(p1, varn(x));
     575           7 :           gel(y,1) = gerepile(av,tetpil,p1);
     576           7 :           gel(y,2) = const_col(lx-1, gen_1); return y;
     577             : 
     578          14 :         case t_REAL: y=cgetg(3,t_MAT); lx=lg(x)-2; v=varn(x);
     579          14 :           av=avma; p1=cleanroots(x,pa); tetpil=avma;
     580          35 :           for(r1=1; r1<lx; r1++)
     581          28 :             if (typ(gel(p1,r1)) == t_COMPLEX) break;
     582          14 :           lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
     583          35 :           for(i=1; i<r1; i++)
     584          21 :             gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
     585          21 :           for(   ; i<lx; i++)
     586             :           {
     587           7 :             GEN a = gel(p1,2*i-r1);
     588           7 :             p = cgetg(5, t_POL); gel(p2,i) = p;
     589           7 :             p[1] = x[1];
     590           7 :             gel(p,2) = gnorm(a);
     591           7 :             gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
     592           7 :             gel(p,4) = gen_1;
     593             :           }
     594          14 :           gel(y,1) = gerepile(av,tetpil,p2);
     595          14 :           gel(y,2) = const_col(lx-1, gen_1); return y;
     596             : 
     597          14 :         case t_PADIC: return factorpadic(x,p,pa);
     598             : 
     599          35 :         case t_FFELT: return FFX_factor(x,pol);
     600             : 
     601             :         default:
     602             :         {
     603             :           GEN w;
     604             :           long killv, t1, t2;
     605         126 :           x = leafcopy(x); lx=lg(x);
     606         126 :           pol = leafcopy(pol);
     607         126 :           v = pari_var_next_temp();
     608        1001 :           for(i=2; i<lx; i++)
     609             :           {
     610         875 :             p1 = gel(x,i);
     611         875 :             switch(typ(p1))
     612             :             {
     613          28 :               case t_QUAD: p1++;
     614             :               case t_COMPLEX:
     615          49 :                 gel(x,i) = mkpolmod(deg1pol_shallow(gel(p1,2), gel(p1,1), v), pol);
     616             :             }
     617             :           }
     618         126 :           killv = (avma != (pari_sp)pol);
     619         126 :           if (killv) setvarn(pol, fetch_var());
     620         126 :           RgX_type_decode(tx, &t1, &t2);
     621         126 :           switch (t2)
     622             :           {
     623          70 :             case t_INT: p1 = polfnf(x,pol); break;
     624             :             case t_INTMOD:
     625          35 :               pol = RgX_to_FpX(pol, p);
     626          35 :               if (FpX_is_squarefree(pol,p) && FpX_nbfact(pol, p) == 1)
     627             :               {
     628          21 :                 p1 = factorff(x,p,pol); break;
     629             :               }
     630             :             /*fall through*/
     631          35 :             default: pari_err_IMPL("factor for general polynomial");
     632           0 :               return NULL; /* not reached */
     633             :           }
     634          91 :           switch (t1)
     635             :           {
     636             :             case t_POLMOD:
     637          63 :               if (killv) (void)delete_var();
     638          63 :               return gerepileupto(av,p1);
     639          14 :             case t_COMPLEX: w = gen_I(); break;
     640          14 :             case t_QUAD: w = mkquad(pol,gen_0,gen_1);
     641          14 :               break;
     642           0 :             default: pari_err_IMPL("factor for general polynomial");
     643           0 :               return NULL; /* not reached */
     644             :           }
     645          28 :           p2=gel(p1,1);
     646          70 :           for(i=1; i<lg(p2); i++)
     647             :           {
     648          42 :             GEN P = gel(p2,i);
     649             :             long j;
     650         140 :             for(j=2; j<lg(P); j++)
     651             :             {
     652          98 :               GEN p = gel(P,j);
     653          98 :               if(typ(p)==t_POLMOD) gel(P,j) = gsubst(gel(p,2),v,w);
     654             :             }
     655             :           }
     656          28 :           if (killv) (void)delete_var();
     657          28 :           return gerepilecopy(av, p1);
     658             :         }
     659             :       }
     660             :     case t_RFRAC: {
     661           7 :       GEN a = gel(x,1), b = gel(x,2);
     662           7 :       y = factor(b); gel(y,2) = gneg_i(gel(y,2));
     663           7 :       if (typ(a)==t_POL && varn(a)==varn(b)) y = famat_mul_shallow(factor(a), y);
     664           7 :       return gerepilecopy(av, y);
     665             :     }
     666             : 
     667             :     case t_COMPLEX:
     668          42 :       y = gauss_factor(x);
     669          42 :       if (y) return y;
     670             :       /* fall through */
     671             :   }
     672           0 :   pari_err_TYPE("factor",x);
     673           0 :   return NULL; /* not reached */
     674             : }
     675             : 
     676             : /*******************************************************************/
     677             : /*                                                                 */
     678             : /*                     ROOTS --> MONIC POLYNOMIAL                  */
     679             : /*                                                                 */
     680             : /*******************************************************************/
     681             : static GEN
     682      412981 : normalized_mul(void *E, GEN x, GEN y)
     683             : {
     684      412981 :   long a = gel(x,1)[1], b = gel(y,1)[1];
     685             :   (void) E;
     686      825962 :   return mkvec2(mkvecsmall(a + b),
     687      825962 :                 RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
     688             : }
     689             : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
     690             : static GEN
     691      260513 : normalized_to_RgX(GEN L)
     692             : {
     693      260513 :   long i, a = gel(L,1)[1];
     694      260513 :   GEN A = gel(L,2);
     695      260513 :   GEN z = cgetg(a + 3, t_POL);
     696      260513 :   z[1] = evalsigne(1) | evalvarn(varn(A));
     697      260513 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
     698      260513 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
     699      260513 :   gel(z,i) = gen_1; return z;
     700             : }
     701             : 
     702             : /* compute prod (x - a[i]) */
     703             : GEN
     704      235672 : roots_to_pol(GEN a, long v)
     705             : {
     706      235672 :   pari_sp av = avma;
     707      235672 :   long i, k, lx = lg(a);
     708             :   GEN L;
     709      235672 :   if (lx == 1) return pol_1(v);
     710      235672 :   L = cgetg(lx, t_VEC);
     711      502405 :   for (k=1,i=1; i<lx-1; i+=2)
     712             :   {
     713      266733 :     GEN s = gel(a,i), t = gel(a,i+1);
     714      266733 :     GEN x0 = gmul(s,t);
     715      266733 :     GEN x1 = gneg(gadd(s,t));
     716      266733 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     717             :   }
     718      463372 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
     719      227700 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
     720      235672 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
     721      235672 :   return gerepileupto(av, normalized_to_RgX(L));
     722             : }
     723             : 
     724             : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
     725             : GEN
     726       24841 : roots_to_pol_r1(GEN a, long v, long r1)
     727             : {
     728       24841 :   pari_sp av = avma;
     729       24841 :   long i, k, lx = lg(a);
     730             :   GEN L;
     731       24841 :   if (lx == 1) return pol_1(v);
     732       24841 :   L = cgetg(lx, t_VEC);
     733      122472 :   for (k=1,i=1; i<r1; i+=2)
     734             :   {
     735       97631 :     GEN s = gel(a,i), t = gel(a,i+1);
     736       97631 :     GEN x0 = gmul(s,t);
     737       97631 :     GEN x1 = gneg(gadd(s,t));
     738       97631 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     739             :   }
     740       30791 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
     741        5950 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
     742      100321 :   for (i=r1+1; i<lx; i++)
     743             :   {
     744       75480 :     GEN s = gel(a,i);
     745       75480 :     GEN x0 = gnorm(s);
     746       75480 :     GEN x1 = gneg(gtrace(s));
     747       75480 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     748             :   }
     749       24841 :   setlg(L, k); L = gen_product(L, NULL, normalized_mul);
     750       24841 :   return gerepileupto(av, normalized_to_RgX(L));
     751             : }
     752             : 
     753             : /*******************************************************************/
     754             : /*                                                                 */
     755             : /*                          FACTORBACK                             */
     756             : /*                                                                 */
     757             : /*******************************************************************/
     758             : static GEN
     759           0 : idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
     760             : static GEN
     761          56 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
     762             : static GEN
     763           7 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
     764             : static GEN
     765        6048 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
     766             : static GEN
     767        2779 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
     768             : static GEN
     769        4494 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
     770             : static GEN
     771     2646153 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
     772             : static GEN
     773     4247352 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
     774             : static GEN
     775           0 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
     776             : static GEN
     777          49 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
     778             : 
     779             : #if 0
     780             : static GEN
     781             : _ellmul(void *ell, GEN x, GEN y) { return elladd((GEN) ell, x, y); }
     782             : static GEN
     783             : _ellpow(void *ell, GEN x, GEN n) { return ellmul((GEN) ell, x, n); }
     784             : #endif
     785             : 
     786             : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
     787             : GEN
     788     1623311 : gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
     789             :                GEN (*_pow)(void*,GEN,GEN), void *data)
     790             : {
     791     1623311 :   pari_sp av = avma;
     792             :   long k, l, lx;
     793             :   GEN p,x;
     794             : 
     795     1623311 :   if (e) /* supplied vector of exponents */
     796      465094 :     p = L;
     797             :   else
     798             :   {
     799     1158217 :     switch(typ(L)) {
     800             :       case t_VEC:
     801             :       case t_COL: /* product of the L[i] */
     802         128 :         return gerepileupto(av, gen_product(L, data, _mul));
     803             :       case t_MAT: /* genuine factorization */
     804     1158089 :         l = lg(L);
     805     1158089 :         if (l == 1) return gen_1;
     806     1158089 :         if (l == 3) break;
     807             :         /*fall through*/
     808             :       default:
     809           7 :         pari_err_TYPE("factorback [not a factorization]", L);
     810             :     }
     811     1158082 :     p = gel(L,1);
     812     1158082 :     e = gel(L,2);
     813             :   }
     814             :   /* p = elts, e = expo */
     815     1623176 :   lx = lg(p);
     816             :   /* check whether e is an integral vector of correct length */
     817     1623176 :   switch(typ(e))
     818             :   {
     819             :     case t_VECSMALL:
     820         252 :       if (lx != lg(e))
     821           0 :         pari_err_TYPE("factorback [not an exponent vector]", e);
     822         252 :       if (lx == 1) return gen_1;
     823         252 :       x = cgetg(lx,t_VEC);
     824         854 :       for (l=1,k=1; k<lx; k++)
     825         602 :         if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
     826         252 :       break;
     827             :     case t_VEC: case t_COL:
     828     1622924 :       if (lx != lg(e) || !RgV_is_ZV(e))
     829           7 :         pari_err_TYPE("factorback [not an exponent vector]", e);
     830     1622917 :       if (lx == 1) return gen_1;
     831     1622586 :       x = cgetg(lx,t_VEC);
     832     5887515 :       for (l=1,k=1; k<lx; k++)
     833     4264929 :         if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
     834     1622586 :       break;
     835             :     default:
     836           0 :       pari_err_TYPE("factorback [not an exponent vector]", e);
     837           0 :       return NULL;
     838             :   }
     839     1622838 :   x[0] = evaltyp(t_VEC) | _evallg(l);
     840     1622838 :   return gerepileupto(av, gen_product(x, data, _mul));
     841             : }
     842             : 
     843             : GEN
     844        6097 : idealfactorback(GEN nf, GEN L, GEN e, int red)
     845             : {
     846        6097 :   nf = checknf(nf);
     847        6097 :   if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
     848        6041 :   else     return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
     849             : }
     850             : 
     851             : GEN
     852        2191 : nffactorback(GEN nf, GEN L, GEN e)
     853        2191 : { return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }
     854             : 
     855             : GEN
     856          49 : FpV_factorback(GEN L, GEN e, GEN p)
     857          49 : { return gen_factorback(L, e, &Fpmul, &Fppow, (void*)p); }
     858             : 
     859             : GEN
     860     1611726 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
     861             : GEN
     862     1157179 : factorback(GEN fa) { return factorback2(fa, NULL); }
     863             : 
     864             : static int
     865          63 : RgX_is_irred_i(GEN x)
     866             : {
     867             :   GEN y, p, pol;
     868          63 :   long l = lg(x), pa;
     869             : 
     870          63 :   if (!signe(x) || l <= 3) return 0;
     871          63 :   switch(RgX_type(x,&p,&pol,&pa))
     872             :   {
     873          14 :     case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
     874           0 :     case t_COMPLEX: return l == 4;
     875             :     case t_REAL:
     876           0 :       if (l == 4) return 1;
     877           0 :       if (l > 5) return 0;
     878           0 :       return gsigne(RgX_disc(x)) > 0;
     879             :   }
     880          49 :   y = factor(x);
     881          49 :   return (lg(gcoeff(y,1,1))==l);
     882             : }
     883             : static int
     884          63 : RgX_is_irred(GEN x)
     885             : {
     886          63 :   pari_sp av = avma;
     887          63 :   int r = RgX_is_irred_i(x);
     888          63 :   avma = av; return r;
     889             : }
     890             : long
     891          63 : isirreducible(GEN x)
     892             : {
     893          63 :   switch(typ(x))
     894             :   {
     895           0 :     case t_INT: case t_REAL: case t_FRAC: return 0;
     896          63 :     case t_POL: return RgX_is_irred(x);
     897             :   }
     898           0 :   pari_err_TYPE("isirreducible",x);
     899           0 :   return 0;
     900             : }
     901             : 
     902             : /*******************************************************************/
     903             : /*                                                                 */
     904             : /*                         GENERIC GCD                             */
     905             : /*                                                                 */
     906             : /*******************************************************************/
     907             : /* x is a COMPLEX or a QUAD */
     908             : static GEN
     909         357 : triv_cont_gcd(GEN x, GEN y)
     910             : {
     911         357 :   pari_sp av = avma;
     912             :   GEN c;
     913         357 :   if (typ(x)==t_COMPLEX)
     914             :   {
     915          35 :     GEN a = gel(x,1), b = gel(x,2);
     916          35 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
     917          21 :     c = ggcd(a,b);
     918             :   }
     919             :   else
     920         322 :     c = ggcd(gel(x,2),gel(x,3));
     921         343 :   return gerepileupto(av, ggcd(c,y));
     922             : }
     923             : 
     924             : /* y is a PADIC, x a rational number or an INTMOD */
     925             : static GEN
     926         147 : padic_gcd(GEN x, GEN y)
     927             : {
     928         147 :   GEN p = gel(y,2);
     929         147 :   long v = gvaluation(x,p), w = valp(y);
     930         147 :   if (w < v) v = w;
     931         147 :   return powis(p, v);
     932             : }
     933             : 
     934             : /* x,y in Z[i], at least one of which is t_COMPLEX */
     935             : static GEN
     936          49 : gauss_gcd(GEN x, GEN y)
     937             : {
     938          49 :   pari_sp av = avma;
     939             :   GEN dx, dy;
     940          49 :   dx = denom(x); x = gmul(x, dx);
     941          49 :   dy = denom(y); y = gmul(y, dy);
     942         140 :   while (!gequal0(y))
     943             :   {
     944          42 :     GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
     945          42 :     x = y; y = z;
     946             :   }
     947          49 :   x = gauss_normal(x);
     948          49 :   if (typ(x) == t_COMPLEX)
     949             :   {
     950          35 :     if      (gequal0(gel(x,2))) x = gel(x,1);
     951          35 :     else if (gequal0(gel(x,1))) x = gel(x,2);
     952             :   }
     953          49 :   return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
     954             : }
     955             : 
     956             : static int
     957          56 : c_is_rational(GEN x)
     958          56 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
     959             : static GEN
     960          28 : c_zero_gcd(GEN c)
     961             : {
     962          28 :   GEN x = gel(c,1), y = gel(c,2);
     963          28 :   long tx = typ(x), ty = typ(y);
     964          28 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
     965          21 :   if (tx == t_PADIC || tx == t_INTMOD
     966          14 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
     967          14 :   return gauss_gcd(c, gen_0);
     968             : }
     969             : 
     970             : /* gcd(x, 0) */
     971             : static GEN
     972     8166132 : zero_gcd(GEN x)
     973             : {
     974             :   pari_sp av;
     975     8166132 :   switch(typ(x))
     976             :   {
     977       35855 :     case t_INT: return absi(x);
     978        1197 :     case t_FRAC: return absfrac(x);
     979          28 :     case t_COMPLEX: return c_zero_gcd(x);
     980        1757 :     case t_REAL: return gen_1;
     981          42 :     case t_PADIC: return powis(gel(x,2), valp(x));
     982         210 :     case t_SER: return monomial(gen_1, valp(x), varn(x));
     983             :     case t_POLMOD: {
     984       13426 :       GEN d = gel(x,2);
     985       13426 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
     986        1386 :       return isinexact(d)? zero_gcd(d): gcopy(d);
     987             :     }
     988             :     case t_POL:
     989     7896445 :       if (!isinexact(x)) break;
     990           0 :       av = avma;
     991           0 :       return gerepileupto(av,
     992           0 :         monomialcopy(content(x), RgX_val(x), varn(x))
     993             :       );
     994             : 
     995             :     case t_RFRAC:
     996      216969 :       if (!isinexact(x)) break;
     997           0 :       av = avma;
     998           0 :       return gerepileupto(av,
     999           0 :         gdiv(zero_gcd(gel(x,1)), gel(x,2))
    1000             :       );
    1001             :   }
    1002     8113617 :   return gcopy(x);
    1003             : }
    1004             : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
    1005             : static GEN
    1006     8155282 : zero_gcd2(GEN y, GEN z)
    1007             : {
    1008             :   pari_sp av;
    1009     8155282 :   switch(typ(z))
    1010             :   {
    1011     8151005 :     case t_INT: return zero_gcd(y);
    1012             :     case t_INTMOD:
    1013        3367 :       av = avma;
    1014        3367 :       return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
    1015             :     case t_FFELT:
    1016         910 :       av = avma;
    1017         910 :       return gerepileupto(av, gmul(y, FF_1(z)));
    1018             :     default:
    1019           0 :       pari_err_TYPE("zero_gcd", z);
    1020             :   }
    1021           0 :   return NULL;
    1022             : }
    1023             : /* tx = t_RFRAC, y considered as constant */
    1024             : static GEN
    1025      864892 : cont_gcd_rfrac(GEN x, GEN y)
    1026             : {
    1027      864892 :   pari_sp av = avma;
    1028      864892 :   GEN cx; x = primitive_part(x, &cx);
    1029      864892 :   return gerepileupto(av, gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2)));
    1030             : }
    1031             : /* tx = t_POL, y considered as constant */
    1032             : static GEN
    1033     2316430 : cont_gcd_pol(GEN x, GEN y)
    1034             : {
    1035     2316430 :   pari_sp av = avma;
    1036     2316430 :   return gerepileupto(av, scalarpol(ggcd(content(x),y), varn(x)));
    1037             : }
    1038             : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
    1039             : static GEN
    1040       10962 : cont_gcd_gen(GEN x, GEN y)
    1041             : {
    1042       10962 :   pari_sp av = avma;
    1043       10962 :   return gerepileupto(av, ggcd(content(x),y));
    1044             : }
    1045             : /* !is_const(tx), y considered as constant */
    1046             : static GEN
    1047     3135788 : cont_gcd(GEN x, long tx, GEN y)
    1048             : {
    1049     3135788 :   switch(tx)
    1050             :   {
    1051      808445 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1052     2316381 :     case t_POL: return cont_gcd_pol(x,y);
    1053       10962 :     default: return cont_gcd_gen(x,y);
    1054             :   }
    1055             : }
    1056             : static GEN
    1057      792139 : gcdiq(GEN x, GEN y)
    1058             : {
    1059             :   GEN z;
    1060      792139 :   if (!signe(x)) return Q_abs(y);
    1061      198022 :   z = cgetg(3,t_FRAC);
    1062      198022 :   gel(z,1) = gcdii(x,gel(y,1));
    1063      198022 :   gel(z,2) = icopy(gel(y,2));
    1064      198022 :   return z;
    1065             : }
    1066             : static GEN
    1067     1122460 : gcdqq(GEN x, GEN y)
    1068             : {
    1069     1122460 :   GEN z = cgetg(3,t_FRAC);
    1070     1122460 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1071     1122460 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1072     1122460 :   return z;
    1073             : }
    1074             : /* assume x,y t_INT or t_FRAC */
    1075             : GEN
    1076    58660826 : Q_gcd(GEN x, GEN y)
    1077             : {
    1078    58660826 :   long tx = typ(x), ty = typ(y);
    1079    58660826 :   if (tx == t_INT)
    1080    56877442 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1081             :   else
    1082     1783384 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1083             : }
    1084             : 
    1085             : GEN
    1086    21596993 : ggcd(GEN x, GEN y)
    1087             : {
    1088    21596993 :   long l, i, vx, vy, tx = typ(x), ty = typ(y);
    1089             :   pari_sp av, tetpil;
    1090             :   GEN p1,z;
    1091             : 
    1092    21596993 :   if (is_noncalc_t(tx) || is_noncalc_t(ty)) pari_err_TYPE2("gcd",x,y);
    1093    21596993 :   if (is_matvec_t(ty))
    1094             :   {
    1095           0 :     z = cgetg_copy(y, &l);
    1096           0 :     for (i=1; i<l; i++) gel(z,i) = ggcd(x,gel(y,i));
    1097           0 :     return z;
    1098             :   }
    1099    21596993 :   if (is_matvec_t(tx))
    1100             :   {
    1101           0 :     z = cgetg_copy(x, &l);
    1102           0 :     for (i=1; i<l; i++) gel(z,i) = ggcd(gel(x,i),y);
    1103           0 :     return z;
    1104             :   }
    1105    21596993 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1106             :   /* tx <= ty */
    1107    21596993 :   z = gisexactzero(x); if (z) return zero_gcd2(y,z);
    1108    18574300 :   z = gisexactzero(y); if (z) return zero_gcd2(x,z);
    1109    13441711 :   if (is_const_t(tx))
    1110             :   {
    1111     6526300 :     if (ty == tx) switch(tx)
    1112             :     {
    1113             :       case t_INT:
    1114     3383435 :         return gcdii(x,y);
    1115             : 
    1116        1211 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1117        1211 :         if (equalii(gel(x,1),gel(y,1)))
    1118        1211 :           gel(z,1) = icopy(gel(x,1));
    1119             :         else
    1120           0 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1121        1211 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1122             :         else
    1123             :         {
    1124        1211 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1125        1211 :           if (!is_pm1(p1))
    1126             :           {
    1127           0 :             p1 = gcdii(p1,gel(y,2));
    1128           0 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1129           0 :             else p1 = gerepileuptoint(av, p1);
    1130             :           }
    1131        1211 :           gel(z,2) = p1;
    1132             :         }
    1133        1211 :         return z;
    1134             : 
    1135             :       case t_FRAC:
    1136       14847 :         return gcdqq(x,y);
    1137             : 
    1138             :       case t_FFELT:
    1139        7931 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1140        7931 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1141             : 
    1142             :       case t_COMPLEX:
    1143          14 :         if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
    1144           7 :         return triv_cont_gcd(y,x);
    1145             : 
    1146             :       case t_PADIC:
    1147          14 :         if (!equalii(gel(x,2),gel(y,2))) return gen_1;
    1148           7 :         return powis(gel(y,2), minss(valp(x), valp(y)));
    1149             : 
    1150             :       case t_QUAD:
    1151         133 :         av=avma; p1=gdiv(x,y);
    1152         133 :         if (gequal0(gel(p1,3)))
    1153             :         {
    1154          14 :           p1=gel(p1,2);
    1155          14 :           if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
    1156           7 :           tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
    1157             :         }
    1158         119 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
    1159         112 :         p1 = ginv(p1); avma=av;
    1160         112 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1161         105 :         return triv_cont_gcd(y,x);
    1162             : 
    1163           0 :       default: return gen_1; /* t_REAL */
    1164             :     }
    1165     3118715 :     if (is_const_t(ty)) switch(tx)
    1166             :     {
    1167             :       case t_INT:
    1168       10724 :         switch(ty)
    1169             :         {
    1170          42 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1171          42 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1172          42 :             p1 = gcdii(gel(y,1),gel(y,2));
    1173          42 :             if (!is_pm1(p1)) {
    1174          14 :               p1 = gcdii(x,p1);
    1175          14 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1176             :               else
    1177          14 :                 p1 = gerepileuptoint(av, p1);
    1178             :             }
    1179          42 :             gel(z,2) = p1; return z;
    1180             : 
    1181        7210 :           case t_REAL: return gen_1;
    1182             : 
    1183             :           case t_FRAC:
    1184        2065 :             return gcdiq(x,y);
    1185             : 
    1186             :           case t_COMPLEX:
    1187          21 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1188           0 :             return triv_cont_gcd(y,x);
    1189             : 
    1190             :           case t_FFELT:
    1191        1127 :             if (!FF_equal0(y)) return FF_1(y);
    1192           0 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1193             : 
    1194             :           case t_PADIC:
    1195         133 :             return padic_gcd(x,y);
    1196             : 
    1197             :           case t_QUAD:
    1198         126 :             return triv_cont_gcd(y,x);
    1199             :           default:
    1200           0 :             pari_err_TYPE2("gcd",x,y);
    1201             :         }
    1202             : 
    1203             :       case t_REAL:
    1204          28 :         switch(ty)
    1205             :         {
    1206             :           case t_INTMOD:
    1207             :           case t_FFELT:
    1208          14 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1209          14 :           default: return gen_1;
    1210             :         }
    1211             : 
    1212             :       case t_INTMOD:
    1213          63 :         switch(ty)
    1214             :         {
    1215             :           case t_FRAC:
    1216          14 :             av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
    1217          14 :             if (!is_pm1(p1)) pari_err_OP("gcd",x,y);
    1218           7 :             return ggcd(gel(y,1), x);
    1219             : 
    1220             :           case t_FFELT:
    1221             :           {
    1222          28 :             GEN p = gel(y,4);
    1223          28 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1224          21 :             if (!FF_equal0(y)) return FF_1(y);
    1225           0 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1226             :           }
    1227             : 
    1228             :           case t_COMPLEX: case t_QUAD:
    1229          14 :             return triv_cont_gcd(y,x);
    1230             : 
    1231             :           case t_PADIC:
    1232           7 :             return padic_gcd(x,y);
    1233             : 
    1234           0 :           default: pari_err_TYPE2("gcd",x,y);
    1235             :         }
    1236             : 
    1237             :       case t_FRAC:
    1238         140 :         switch(ty)
    1239             :         {
    1240             :           case t_COMPLEX:
    1241          14 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1242             :           case t_QUAD:
    1243          84 :             return triv_cont_gcd(y,x);
    1244             :           case t_FFELT:
    1245             :           {
    1246          42 :             GEN p = gel(y,4);
    1247          42 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1248          21 :             if (!FF_equal0(y)) return FF_1(y);
    1249           0 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1250             :           }
    1251             : 
    1252             :           case t_PADIC:
    1253           7 :             return padic_gcd(x,y);
    1254             : 
    1255           0 :           default: pari_err_TYPE2("gcd",x,y);
    1256             :         }
    1257             :       case t_FFELT:
    1258          70 :         switch(ty)
    1259             :         {
    1260             :           case t_PADIC:
    1261             :           {
    1262          42 :             GEN p = gel(y,2);
    1263          42 :             long v = valp(y);
    1264          42 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1265          14 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1266             :           }
    1267          28 :           default: pari_err_TYPE2("gcd",x,y);
    1268             :         }
    1269             : 
    1270             :       case t_COMPLEX:
    1271          14 :         switch(ty)
    1272             :         {
    1273             :           case t_PADIC:
    1274          14 :           case t_QUAD: return triv_cont_gcd(x,y);
    1275           0 :           default: pari_err_TYPE2("gcd",x,y);
    1276             :         }
    1277             : 
    1278             :       case t_PADIC:
    1279           7 :         switch(ty)
    1280             :         {
    1281           7 :           case t_QUAD: return triv_cont_gcd(y,x);
    1282           0 :           default: pari_err_TYPE2("gcd",x,y);
    1283             :         }
    1284             : 
    1285           0 :       default: return gen_1; /* tx = t_REAL */
    1286             :     }
    1287     3107669 :     return cont_gcd(y,ty, x);
    1288             :   }
    1289             : 
    1290     6915411 :   if (tx == t_POLMOD)
    1291             :   {
    1292        9779 :     if (ty == t_POLMOD)
    1293             :     {
    1294        9443 :       GEN T = gel(x,1);
    1295        9443 :       z = cgetg(3,t_POLMOD);
    1296        9443 :       T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
    1297        9443 :       gel(z,1) = T;
    1298        9443 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1299             :       else
    1300             :       {
    1301             :         GEN X, Y, d;
    1302        9443 :         av = avma; X = gel(x,2); Y = gel(y,2);
    1303        9443 :         d = ggcd(content(X), content(Y));
    1304        9443 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1305        9443 :         p1 = ggcd(T, X);
    1306        9443 :         gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
    1307             :       }
    1308        9443 :       return z;
    1309             :     }
    1310         336 :     vx = varn(gel(x,1));
    1311         336 :     switch(ty)
    1312             :     {
    1313             :       case t_POL:
    1314         308 :         vy = varn(y);
    1315         308 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1316         259 :         z = cgetg(3,t_POLMOD);
    1317         259 :         gel(z,1) = RgX_copy(gel(x,1));
    1318         259 :         av = avma; p1 = ggcd(gel(x,1),gel(x,2));
    1319         259 :         gel(z,2) = gerepileupto(av, ggcd(p1,y));
    1320         259 :         return z;
    1321             : 
    1322             :       case t_RFRAC:
    1323          28 :         vy = varn(gel(y,2));
    1324          28 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1325          28 :         av = avma;
    1326          28 :         p1 = ggcd(gel(x,1),gel(y,2));
    1327          28 :         if (degpol(p1)) pari_err_OP("gcd",x,y);
    1328          21 :         avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1329             :     }
    1330             :   }
    1331             : 
    1332     6905632 :   vx = gvar(x);
    1333     6905632 :   vy = gvar(y);
    1334     6905632 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1335     6893536 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1336             : 
    1337             :   /* vx = vy: same main variable */
    1338     6877513 :   switch(tx)
    1339             :   {
    1340             :     case t_POL:
    1341     6716607 :       switch(ty)
    1342             :       {
    1343     6660146 :         case t_POL: return RgX_gcd(x,y);
    1344             :         case t_SER:
    1345          14 :           z = ggcd(content(x), content(y));
    1346          14 :           return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
    1347       56447 :         case t_RFRAC: return cont_gcd_rfrac(y, x);
    1348             :       }
    1349           0 :       break;
    1350             : 
    1351             :     case t_SER:
    1352          14 :       z = ggcd(content(x), content(y));
    1353          14 :       switch(ty)
    1354             :       {
    1355           7 :         case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
    1356           7 :         case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
    1357             :       }
    1358           0 :       break;
    1359             : 
    1360             :     case t_RFRAC:
    1361             :     {
    1362      160892 :       GEN xd = gel(x,2), yd = gel(y,2);
    1363      160892 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1364      160892 :       z = cgetg(3,t_RFRAC); av = avma;
    1365      160892 :       gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1366      160892 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1367             :     }
    1368             :   }
    1369           0 :   pari_err_TYPE2("gcd",x,y);
    1370           0 :   return NULL; /* not reached */
    1371             : }
    1372             : GEN
    1373        3798 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
    1374             : 
    1375             : /* x a t_VEC,t_COL or t_MAT */
    1376             : static GEN
    1377           0 : vec_lcm(GEN x)
    1378             : {
    1379           0 :   if (typ(x) == t_MAT)
    1380             :   {
    1381           0 :     long i, l = lg(x);
    1382           0 :     GEN z = cgetg(l, t_VEC);
    1383           0 :     for (i = 1; i < l; i++) gel(z,i) = glcm0(gel(x,i), NULL);
    1384           0 :     x = z;
    1385             :   }
    1386           0 :   return glcm0(x, NULL);
    1387             : }
    1388             : static GEN
    1389        2905 : scal_lcm(GEN x, GEN y)
    1390             : {
    1391        2905 :   pari_sp av = avma;
    1392        2905 :   long tx = typ(x), ty = typ(y);
    1393        2905 :   if (is_matvec_t(tx)) x = vec_lcm(x);
    1394        2905 :   if (is_matvec_t(ty)) y = vec_lcm(y);
    1395        2905 :   return gerepileupto(av, glcm(x, y));
    1396             : }
    1397             : 
    1398             : static GEN
    1399          63 : fix_lcm(GEN x)
    1400             : {
    1401             :   GEN t;
    1402          63 :   switch(typ(x))
    1403             :   {
    1404           7 :     case t_INT: if (signe(x)<0) x = negi(x);
    1405           7 :       break;
    1406             :     case t_POL:
    1407          56 :       if (lg(x) <= 2) break;
    1408          56 :       t = leading_coeff(x);
    1409          56 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1410             :   }
    1411          63 :   return x;
    1412             : }
    1413             : 
    1414             : GEN
    1415        2884 : glcm0(GEN x, GEN y) {
    1416        2884 :   if (!y && lg(x) == 2)
    1417             :   {
    1418          14 :     long tx = typ(x);
    1419          14 :     if (is_vec_t(tx))
    1420             :     {
    1421          14 :       x = gel(x,1);
    1422          14 :       tx = typ(x);
    1423          14 :       return is_matvec_t(tx)? vec_lcm(x): fix_lcm(x);
    1424             :     }
    1425             :   }
    1426        2870 :   return gassoc_proto(scal_lcm,x,y);
    1427             : }
    1428             : 
    1429             : GEN
    1430        3150 : glcm(GEN x, GEN y)
    1431             : {
    1432             :   long tx, ty, i, l;
    1433             :   pari_sp av;
    1434             :   GEN p1, z;
    1435             : 
    1436        3150 :   ty = typ(y);
    1437        3150 :   if (is_matvec_t(ty))
    1438             :   {
    1439           0 :     z = cgetg_copy(y, &l);
    1440           0 :     for (i=1; i<l; i++) gel(z,i) = glcm(x,gel(y,i));
    1441           0 :     return z;
    1442             :   }
    1443        3150 :   tx = typ(x);
    1444        3150 :   if (is_matvec_t(tx))
    1445             :   {
    1446           0 :     z = cgetg_copy(x, &l);
    1447           0 :     for (i=1; i<l; i++) gel(z,i) = glcm(gel(x,i),y);
    1448           0 :     return z;
    1449             :   }
    1450        3150 :   if (tx==t_INT && ty==t_INT) return lcmii(x,y);
    1451          49 :   if (gequal0(x)) return gen_0;
    1452             : 
    1453          49 :   av = avma;
    1454          49 :   p1 = ggcd(x,y); if (!gequal1(p1)) y = gdiv(y,p1);
    1455          49 :   return gerepileupto(av, fix_lcm(gmul(x,y)));
    1456             : }
    1457             : 
    1458             : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
    1459             : static int
    1460        2107 : pol_approx0(GEN r, GEN x, int exact)
    1461             : {
    1462             :   long i, lx,lr;
    1463        2107 :   if (exact) return gequal0(r);
    1464           0 :   lx = lg(x);
    1465           0 :   lr = lg(r); if (lr < lx) lx = lr;
    1466           0 :   for (i=2; i<lx; i++)
    1467           0 :     if (!approx_0(gel(r,i), gel(x,i))) return 0;
    1468           0 :   return 1;
    1469             : }
    1470             : 
    1471             : GEN
    1472         777 : RgX_gcd_simple(GEN x, GEN y)
    1473             : {
    1474         777 :   pari_sp av1, av = avma;
    1475         777 :   GEN r, yorig = y;
    1476         777 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1477             : 
    1478             :   for(;;)
    1479             :   {
    1480        2107 :     av1 = avma; r = RgX_rem(x,y);
    1481        2107 :     if (pol_approx0(r, x, exact))
    1482             :     {
    1483         777 :       avma = av1;
    1484         777 :       if (y == yorig) return RgX_copy(y);
    1485         665 :       y = normalizepol_approx(y, lg(y));
    1486         665 :       if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
    1487          98 :       return gerepileupto(av,y);
    1488             :     }
    1489        1330 :     x = y; y = r;
    1490        1330 :     if (gc_needed(av,1)) {
    1491           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1492           0 :       gerepileall(av,2, &x,&y);
    1493             :     }
    1494        1330 :   }
    1495             : }
    1496             : GEN
    1497           0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
    1498             : {
    1499           0 :   pari_sp av = avma;
    1500             :   GEN q, r, d, d1, u, v, v1;
    1501           0 :   int exact = !(isinexactreal(a) || isinexactreal(b));
    1502             : 
    1503           0 :   d = a; d1 = b; v = gen_0; v1 = gen_1;
    1504             :   for(;;)
    1505             :   {
    1506           0 :     if (pol_approx0(d1, a, exact)) break;
    1507           0 :     q = poldivrem(d,d1, &r);
    1508           0 :     v = gsub(v, gmul(q,v1));
    1509           0 :     u=v; v=v1; v1=u;
    1510           0 :     u=r; d=d1; d1=u;
    1511           0 :   }
    1512           0 :   u = gsub(d, gmul(b,v));
    1513           0 :   u = RgX_div(u,a);
    1514             : 
    1515           0 :   gerepileall(av, 3, &u,&v,&d);
    1516           0 :   *pu = u;
    1517           0 :   *pv = v; return d;
    1518             : }
    1519             : /*******************************************************************/
    1520             : /*                                                                 */
    1521             : /*                    CONTENT / PRIMITIVE PART                     */
    1522             : /*                                                                 */
    1523             : /*******************************************************************/
    1524             : 
    1525             : GEN
    1526    73367792 : content(GEN x)
    1527             : {
    1528    73367792 :   long lx, i, t, tx = typ(x);
    1529    73367792 :   pari_sp av = avma;
    1530             :   GEN c;
    1531             : 
    1532    73367792 :   if (is_scalar_t(tx)) return zero_gcd(x);
    1533    73354660 :   switch(tx)
    1534             :   {
    1535             :     case t_RFRAC:
    1536             :     {
    1537      864906 :       GEN n = gel(x,1), d = gel(x,2);
    1538             :       /* -- varncmp(vn, vd) < 0 can't happen
    1539             :        * -- if n is POLMOD, its main variable (in the sense of gvar2)
    1540             :        *    has lower priority than denominator */
    1541      864906 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    1542      824650 :         n = isinexact(n)? zero_gcd(n): gcopy(n);
    1543             :       else
    1544       40256 :         n = content(n);
    1545      864906 :       return gerepileupto(av, gdiv(n, content(d)));
    1546             :     }
    1547             : 
    1548             :     case t_VEC: case t_COL:
    1549     7290816 :       lx = lg(x); if (lx==1) return gen_0;
    1550     7290809 :       break;
    1551             : 
    1552             :     case t_MAT:
    1553             :     {
    1554             :       long hx, j;
    1555         294 :       lx = lg(x);
    1556         294 :       if (lx == 1) return gen_0;
    1557         287 :       hx = lgcols(x);
    1558         287 :       if (hx == 1) return gen_0;
    1559         280 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    1560         280 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    1561         280 :       c = content(gel(x,1));
    1562         798 :       for (j=2; j<lx; j++)
    1563         518 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    1564         280 :       if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
    1565         280 :       return gerepileupto(av,c);
    1566             :     }
    1567             : 
    1568             :     case t_POL: case t_SER:
    1569    65198609 :       lx = lg(x); if (lx == 2) return gen_0;
    1570    65198595 :       break;
    1571          21 :     case t_VECSMALL: return utoi(zv_content(x));
    1572             :     case t_QFR: case t_QFI:
    1573          14 :       lx = 4; break;
    1574             : 
    1575           0 :     default: pari_err_TYPE("content",x);
    1576           0 :       return NULL; /* not reached */
    1577             :   }
    1578   207717197 :   for (i=lontyp[tx]; i<lx; i++)
    1579   147883352 :     if (typ(gel(x,i)) != t_INT) break;
    1580    72489418 :   lx--; c = gel(x,lx);
    1581    72489418 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    1582    72489418 :   if (i > lx)
    1583             :   { /* integer coeffs */
    1584   120520926 :     while (lx-- > lontyp[tx])
    1585             :     {
    1586    59989089 :       c = gcdii(c, gel(x,lx));
    1587    59989089 :       if (is_pm1(c)) { avma=av; return gen_1; }
    1588             :     }
    1589             :   }
    1590             :   else
    1591             :   {
    1592    12655573 :     if (isinexact(c)) c = zero_gcd(c);
    1593    43515430 :     while (lx-- > lontyp[tx])
    1594             :     {
    1595    18204284 :       GEN d = gel(x,lx);
    1596    18204284 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    1597    18204284 :       c = ggcd(c, d);
    1598             :     }
    1599    12655573 :     if (isinexact(c)) { avma=av; return gen_1; }
    1600             :   }
    1601    13353565 :   switch(typ(c))
    1602             :   {
    1603             :     case t_INT:
    1604      708170 :       if (signe(c) < 0) c = negi(c);
    1605      708170 :       break;
    1606             :     case t_VEC: case t_COL: case t_MAT:
    1607           0 :       pari_err_TYPE("content",x);
    1608             :   }
    1609             : 
    1610    13353565 :   return av==avma? gcopy(c): gerepileupto(av,c);
    1611             : }
    1612             : 
    1613             : GEN
    1614     1163504 : primitive_part(GEN x, GEN *ptc)
    1615             : {
    1616     1163504 :   pari_sp av = avma;
    1617     1163504 :   GEN c = content(x);
    1618     1163504 :   if (gequal1(c)) { avma = av; c = NULL; }
    1619       32017 :   else if (!gequal0(c)) x = gdiv(x,c);
    1620     1163504 :   if (ptc) *ptc = c;
    1621     1163504 :   return x;
    1622             : }
    1623             : GEN
    1624        3101 : primpart(GEN x) { return primitive_part(x, NULL); }
    1625             : 
    1626             : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1627             :  * of Q(x2,...,xn)[x1] */
    1628             : GEN
    1629    81795464 : Q_content(GEN x)
    1630             : {
    1631             :   long i, l;
    1632             :   GEN d;
    1633             :   pari_sp av;
    1634             : 
    1635    81795464 :   switch(typ(x))
    1636             :   {
    1637    64696000 :     case t_INT:  return absi(x);
    1638     1271635 :     case t_FRAC: return absfrac(x);
    1639             : 
    1640             :     case t_VEC: case t_COL: case t_MAT:
    1641     7230693 :       l = lg(x); if (l == 1) return gen_1;
    1642     7230658 :       av = avma; d = Q_content(gel(x,1));
    1643    37914441 :       for (i=2; i<l; i++)
    1644             :       {
    1645    30683783 :         d = Q_gcd(d, Q_content(gel(x,i)));
    1646    30683783 :         if ((i & 255) == 0) d = gerepileupto(av, d);
    1647             :       }
    1648     7230658 :       return gerepileupto(av, d);
    1649             : 
    1650             :     case t_POL:
    1651     8597115 :       l = lg(x); if (l == 2) return gen_0;
    1652     8594609 :       av = avma; d = Q_content(gel(x,2));
    1653     8594609 :       for (i=3; i<l; i++) d = Q_gcd(d, Q_content(gel(x,i)));
    1654     8594609 :       return gerepileupto(av, d);
    1655           0 :     case t_POLMOD: return Q_content(gel(x,2));
    1656          21 :     case t_COMPLEX: return Q_gcd(Q_content(gel(x,1)), Q_content(gel(x,2)));
    1657             :   }
    1658           0 :   pari_err_TYPE("Q_content",x);
    1659           0 :   return NULL; /* not reached */
    1660             : }
    1661             : 
    1662             : GEN
    1663       13153 : ZX_content(GEN x)
    1664             : {
    1665       13153 :   long i, l = lg(x);
    1666             :   GEN d;
    1667             :   pari_sp av;
    1668             : 
    1669       13153 :   if (l == 2) return gen_0;
    1670       13153 :   d = gel(x,2);
    1671       13153 :   if (l == 3) return absi(d);
    1672        9177 :   av = avma;
    1673        9177 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    1674        9177 :   if (signe(d) < 0) d = absi(d);
    1675        9177 :   return gerepileuptoint(av, d);
    1676             : }
    1677             : 
    1678             : /* NOT MEMORY CLEAN (because of t_FRAC).
    1679             :  * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1680             :  * of Q(x2,...,xn)[x1] */
    1681             : GEN
    1682    20980594 : Q_denom(GEN x)
    1683             : {
    1684             :   long i, l;
    1685             :   GEN d, D;
    1686             :   pari_sp av;
    1687             : 
    1688    20980594 :   switch(typ(x))
    1689             :   {
    1690    13402730 :     case t_INT: return gen_1;
    1691     4174882 :     case t_FRAC: return gel(x,2);
    1692             : 
    1693             :     case t_VEC: case t_COL: case t_MAT:
    1694     2955566 :       l = lg(x); if (l == 1) return gen_1;
    1695     2954068 :       av = avma; d = Q_denom(gel(x,1));
    1696    15879008 :       for (i=2; i<l; i++)
    1697             :       {
    1698    12924940 :         D = Q_denom(gel(x,i));
    1699    12924940 :         if (D != gen_1) d = lcmii(d, D);
    1700    12924940 :         if ((i & 255) == 0) d = gerepileuptoint(av, d);
    1701             :       }
    1702     2954068 :       return gerepileuptoint(av, d);
    1703             : 
    1704             :     case t_POL:
    1705      446562 :       l = lg(x); if (l == 2) return gen_1;
    1706      446212 :       av = avma; d = Q_denom(gel(x,2));
    1707     1800109 :       for (i=3; i<l; i++)
    1708             :       {
    1709     1353897 :         D = Q_denom(gel(x,i));
    1710     1353897 :         if (D != gen_1) d = lcmii(d, D);
    1711             :       }
    1712      446212 :       return gerepileuptoint(av, d);
    1713             : 
    1714         854 :     case t_POLMOD: return Q_denom(gel(x,2));
    1715             :   }
    1716           0 :   pari_err_TYPE("Q_denom",x);
    1717           0 :   return NULL; /* not reached */
    1718             : }
    1719             : 
    1720             : GEN
    1721     2027373 : Q_remove_denom(GEN x, GEN *ptd)
    1722             : {
    1723     2027373 :   GEN d = Q_denom(x);
    1724     2027373 :   if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d);
    1725     2027373 :   if (ptd) *ptd = d;
    1726     2027373 :   return x;
    1727             : }
    1728             : 
    1729             : /* return y = x * d, assuming x rational, and d,y integral */
    1730             : GEN
    1731    11742484 : Q_muli_to_int(GEN x, GEN d)
    1732             : {
    1733             :   long i, l;
    1734             :   GEN y, xn, xd;
    1735             :   pari_sp av;
    1736             : 
    1737    11742484 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    1738    11742484 :   switch (typ(x))
    1739             :   {
    1740             :     case t_INT:
    1741     5461141 :       return mulii(x,d);
    1742             : 
    1743             :     case t_FRAC:
    1744     4554105 :       xn = gel(x,1);
    1745     4554105 :       xd = gel(x,2); av = avma;
    1746     4554105 :       y = mulii(xn, diviiexact(d, xd));
    1747     4554105 :       return gerepileuptoint(av, y);
    1748             : 
    1749             :     case t_VEC: case t_COL: case t_MAT:
    1750     1274381 :       y = cgetg_copy(x, &l);
    1751     1274381 :       for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    1752     1274381 :       return y;
    1753             : 
    1754             :     case t_POL:
    1755      452857 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1756      452857 :       for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    1757      452857 :       return y;
    1758             : 
    1759             :     case t_POLMOD:
    1760           0 :       y = cgetg(3, t_POLMOD);
    1761           0 :       gel(y,1) = RgX_copy(gel(x,1));
    1762           0 :       gel(y,2) = Q_muli_to_int(gel(x,2), d);
    1763           0 :       return y;
    1764             :   }
    1765           0 :   pari_err_TYPE("Q_muli_to_int",x);
    1766           0 :   return NULL; /* not reached */
    1767             : }
    1768             : 
    1769             : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    1770             : static GEN
    1771      146629 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    1772             : {
    1773             :   long i, l;
    1774             :   GEN y, xn, xd;
    1775             :   pari_sp av;
    1776             : 
    1777      146629 :   switch(typ(x))
    1778             :   {
    1779             :     case t_INT:
    1780       50176 :       av = avma; y = diviiexact(x,d);
    1781       50176 :       return gerepileuptoint(av, mulii(y,n));
    1782             : 
    1783             :     case t_FRAC:
    1784       83767 :       xn = gel(x,1);
    1785       83767 :       xd = gel(x,2); av = avma;
    1786       83767 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    1787       83767 :       return gerepileuptoint(av, y);
    1788             : 
    1789             :     case t_VEC: case t_COL: case t_MAT:
    1790       11013 :       y = cgetg_copy(x, &l);
    1791       11013 :       for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    1792       11013 :       return y;
    1793             : 
    1794             :     case t_POL:
    1795        1673 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1796        1673 :       for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    1797        1673 :       return y;
    1798             : 
    1799             :     case t_POLMOD:
    1800           0 :       y = cgetg(3, t_POLMOD);
    1801           0 :       gel(y,1) = RgX_copy(gel(x,1));
    1802           0 :       gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
    1803           0 :       return y;
    1804             :   }
    1805           0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    1806           0 :   return NULL; /* not reached */
    1807             : }
    1808             : 
    1809             : /* return x / d. x: rational; d,result: integral. */
    1810             : static GEN
    1811    10136040 : Q_divi_to_int(GEN x, GEN d)
    1812             : {
    1813             :   long i, l;
    1814             :   GEN y;
    1815             : 
    1816    10136040 :   switch(typ(x))
    1817             :   {
    1818             :     case t_INT:
    1819     8328130 :       return diviiexact(x,d);
    1820             : 
    1821             :     case t_VEC: case t_COL: case t_MAT:
    1822      646844 :       y = cgetg_copy(x, &l);
    1823      646844 :       for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    1824      646844 :       return y;
    1825             : 
    1826             :     case t_POL:
    1827     1161066 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1828     1161066 :       for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    1829     1161066 :       return y;
    1830             : 
    1831             :     case t_POLMOD:
    1832           0 :       y = cgetg(3, t_POLMOD);
    1833           0 :       gel(y,1) = RgX_copy(gel(x,1));
    1834           0 :       gel(y,2) = Q_divi_to_int(gel(x,2), d);
    1835           0 :       return y;
    1836             :   }
    1837           0 :   pari_err_TYPE("Q_divi_to_int",x);
    1838           0 :   return NULL; /* not reached */
    1839             : }
    1840             : /* c t_FRAC */
    1841             : static GEN
    1842      164659 : Q_divq_to_int(GEN x, GEN c)
    1843             : {
    1844      164659 :   GEN n = gel(c,1), d = gel(c,2);
    1845      164659 :   if (is_pm1(n)) {
    1846      159170 :     GEN y = Q_muli_to_int(x,d);
    1847      159170 :     if (signe(n) < 0) y = gneg(y);
    1848      159170 :     return y;
    1849             :   }
    1850        5489 :   return Q_divmuli_to_int(x, n,d);
    1851             : }
    1852             : 
    1853             : /* return y = x / c, assuming x,c rational, and y integral */
    1854             : GEN
    1855        2086 : Q_div_to_int(GEN x, GEN c)
    1856             : {
    1857        2086 :   switch(typ(c))
    1858             :   {
    1859        2044 :     case t_INT:  return Q_divi_to_int(x, c);
    1860          42 :     case t_FRAC: return Q_divq_to_int(x, c);
    1861             :   }
    1862           0 :   pari_err_TYPE("Q_div_to_int",c);
    1863           0 :   return NULL; /* not reached */
    1864             : }
    1865             : /* return y = x * c, assuming x,c rational, and y integral */
    1866             : GEN
    1867           0 : Q_mul_to_int(GEN x, GEN c)
    1868             : {
    1869             :   GEN d, n;
    1870           0 :   switch(typ(c))
    1871             :   {
    1872           0 :     case t_INT: return Q_muli_to_int(x, c);
    1873             :     case t_FRAC:
    1874           0 :       n = gel(c,1);
    1875           0 :       d = gel(c,2);
    1876           0 :       return Q_divmuli_to_int(x, d,n);
    1877             :   }
    1878           0 :   pari_err_TYPE("Q_mul_to_int",c);
    1879           0 :   return NULL; /* not reached */
    1880             : }
    1881             : 
    1882             : GEN
    1883    10538224 : Q_primitive_part(GEN x, GEN *ptc)
    1884             : {
    1885    10538224 :   pari_sp av = avma;
    1886    10538224 :   GEN c = Q_content(x);
    1887    10538224 :   if (typ(c) == t_INT)
    1888             :   {
    1889    10373607 :     if (is_pm1(c)) { avma = av; c = NULL; }
    1890     1532310 :     else if (signe(c)) x = Q_divi_to_int(x, c);
    1891             :   }
    1892      164617 :   else x = Q_divq_to_int(x, c);
    1893    10538224 :   if (ptc) *ptc = c;
    1894    10538224 :   return x;
    1895             : }
    1896             : GEN
    1897     1239919 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    1898             : 
    1899             : /*******************************************************************/
    1900             : /*                                                                 */
    1901             : /*                           SUBRESULTANT                          */
    1902             : /*                                                                 */
    1903             : /*******************************************************************/
    1904             : /* for internal use */
    1905             : GEN
    1906     3253632 : gdivexact(GEN x, GEN y)
    1907             : {
    1908             :   long i,lx;
    1909             :   GEN z;
    1910     3253632 :   if (gequal1(y)) return x;
    1911     3243419 :   switch(typ(x))
    1912             :   {
    1913             :     case t_INT:
    1914     2452044 :       if (typ(y)==t_INT) return diviiexact(x,y);
    1915         238 :       if (!signe(x)) return gen_0;
    1916           0 :       break;
    1917             :     case t_INTMOD:
    1918       12779 :     case t_POLMOD: return gmul(x,ginv(y));
    1919             :     case t_POL:
    1920      778232 :       switch(typ(y))
    1921             :       {
    1922             :         case t_INTMOD:
    1923           0 :         case t_POLMOD: return gmul(x,ginv(y));
    1924             :         case t_POL: { /* not stack-clean */
    1925             :           long v;
    1926       25011 :           if (varn(x)!=varn(y)) break;
    1927       24164 :           v = RgX_valrem(y,&y);
    1928       24164 :           if (v) x = RgX_shift_shallow(x,-v);
    1929       24164 :           if (!degpol(y)) { y = gel(y,2); break; }
    1930       20783 :           return RgX_div(x,y);
    1931             :         }
    1932             :       }
    1933      757449 :       return RgX_Rg_divexact(x, y);
    1934             :     case t_VEC: case t_COL: case t_MAT:
    1935          98 :       lx = lg(x); z = new_chunk(lx);
    1936          98 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1937          98 :       z[0] = x[0]; return z;
    1938             :   }
    1939         266 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    1940         266 :   return gdiv(x,y);
    1941             : }
    1942             : 
    1943             : static GEN
    1944       87702 : init_resultant(GEN x, GEN y)
    1945             : {
    1946       87702 :   long tx = typ(x), ty = typ(y), vx, vy;
    1947       87702 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    1948             :   {
    1949           0 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    1950           0 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    1951           0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    1952           0 :     return gen_1;
    1953             :   }
    1954       87702 :   if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
    1955       87702 :   if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
    1956       87702 :   if (!signe(x) || !signe(y)) return gmul(RgX_get_0(x),RgX_get_0(y)); /*type*/
    1957       87625 :   vx = varn(x);
    1958       87625 :   vy = varn(y); if (vx == vy) return NULL;
    1959           0 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    1960             : }
    1961             : 
    1962             : static long
    1963      178687 : RgX_simpletype(GEN x)
    1964             : {
    1965      178687 :   long T = t_INT, i, lx = lg(x);
    1966     1063460 :   for (i = 2; i < lx; i++)
    1967             :   {
    1968      884787 :     GEN c = gel(x,i);
    1969      884787 :     long tc = typ(c);
    1970      884787 :     switch(tc) {
    1971             :       case t_INT:
    1972      779113 :         break;
    1973             :       case t_FRAC:
    1974       30865 :         if (T == t_INT) T = t_FRAC;
    1975       30865 :         break;
    1976             :       default:
    1977       74809 :         if (isinexact(c)) return t_REAL;
    1978       74795 :         T = 0; break;
    1979             :     }
    1980             :   }
    1981      178673 :   return T;
    1982             : }
    1983             : 
    1984             : /* x an RgX, y a scalar */
    1985             : static GEN
    1986           0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    1987             : {
    1988           0 :   *V = gpowgs(y,degpol(x)-1);
    1989           0 :   *U = gen_0; return gmul(y, *V);
    1990             : }
    1991             : 
    1992             : /* return 0 if the subresultant chain can be interrupted.
    1993             :  * Set u = NULL if the resultant is 0. */
    1994             : static int
    1995       24028 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    1996             : {
    1997       24028 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    1998             :   long du, dv, dr, degq;
    1999             : 
    2000       24028 :   if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
    2001       24028 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    2002       23902 :   du = degpol(*u);
    2003       23902 :   dv = degpol(*v);
    2004       23902 :   degq = du - dv;
    2005       23902 :   if (*um1 == gen_1)
    2006       11043 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    2007       12859 :   else if (*um1 == gen_0)
    2008        5054 :     u0 = gen_0;
    2009             :   else /* except in those 2 cases, um1 is an RgX */
    2010        7805 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2011             : 
    2012       23902 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2013       11043 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2014             :   else
    2015       12859 :     u0 = gsub(u0, gmul(q,*uze));
    2016             : 
    2017       23902 :   *um1 = *uze;
    2018       23902 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2019             : 
    2020       23902 :   *u = *v; c = *g; *g  = leading_coeff(*u);
    2021       23902 :   switch(degq)
    2022             :   {
    2023          70 :     case 0: break;
    2024             :     case 1:
    2025       16041 :       c = gmul(*h,c); *h = *g; break;
    2026             :     default:
    2027        7791 :       c = gmul(gpowgs(*h,degq), c);
    2028        7791 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2029             :   }
    2030       23902 :   *v  = RgX_Rg_divexact(r,c);
    2031       23902 :   *uze= RgX_Rg_divexact(*uze,c);
    2032       23902 :   if (both_odd(du, dv)) *signh = -*signh;
    2033       23902 :   return (dr > 3);
    2034             : }
    2035             : 
    2036             : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2037             : static GEN
    2038       11673 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2039             : {
    2040             :   pari_sp av, av2;
    2041       11673 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2042             :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2043             : 
    2044       11673 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2045       11673 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2046       11673 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2047       11666 :   if (tx != t_POL) {
    2048           0 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2049           0 :     return scalar_res(y,x,V,U);
    2050             :   }
    2051       11666 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2052       11666 :   if (varn(x) != varn(y))
    2053           0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2054           0 :                                         : scalar_res(y,x,V,U);
    2055       11666 :   if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
    2056       11666 :   if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
    2057       11666 :   dx = degpol(x);
    2058       11666 :   dy = degpol(y);
    2059       11666 :   signh = 1;
    2060       11666 :   if (dx < dy)
    2061             :   {
    2062       10721 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2063       10721 :     if (both_odd(dx, dy)) signh = -signh;
    2064             :   }
    2065       11666 :   if (dy == 0)
    2066             :   {
    2067         714 :     *V = gpowgs(gel(y,2),dx-1);
    2068         714 :     *U = gen_0; return gmul(*V,gel(y,2));
    2069             :   }
    2070       10952 :   av = avma;
    2071       10952 :   u = x = primitive_part(x, &cu);
    2072       10952 :   v = y = primitive_part(y, &cv);
    2073       10952 :   g = h = gen_1; av2 = avma;
    2074       10952 :   um1 = gen_1; uze = gen_0;
    2075             :   for(;;)
    2076             :   {
    2077       23734 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2078       12782 :     if (gc_needed(av2,1))
    2079             :     {
    2080           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2081           0 :       gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    2082             :     }
    2083       12782 :   }
    2084             :   /* uze an RgX */
    2085       10952 :   if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
    2086       10945 :   z = gel(v,2); du = degpol(u);
    2087       10945 :   if (du > 1)
    2088             :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2089         683 :     p1 = gpowgs(gdiv(z,h),du-1);
    2090         683 :     z = gmul(z,p1);
    2091         683 :     uze = RgX_Rg_mul(uze, p1);
    2092             :   }
    2093       10945 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2094             : 
    2095       10945 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2096       10945 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2097             :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2098       10945 :   p1 = gen_1;
    2099       10945 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2100       10945 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2101       10945 :   cu = cu? gdiv(p1,cu): p1;
    2102       10945 :   cv = cv? gdiv(p1,cv): p1;
    2103       10945 :   z = gmul(z,p1);
    2104       10945 :   *U = RgX_Rg_mul(uze,cu);
    2105       10945 :   *V = RgX_Rg_mul(vze,cv);
    2106       10945 :   return z;
    2107             : }
    2108             : GEN
    2109           0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2110             : {
    2111           0 :   pari_sp av = avma;
    2112           0 :   GEN z = subresext_i(x, y, U, V);
    2113           0 :   gerepileall(av, 3, &z, U, V);
    2114           0 :   return z;
    2115             : }
    2116             : 
    2117             : static GEN
    2118          28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2119             : {
    2120          28 :   GEN x=content(y);
    2121          28 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2122             : }
    2123             : 
    2124             : static int
    2125        3178 : must_negate(GEN x)
    2126             : {
    2127        3178 :   GEN t = leading_coeff(x);
    2128        3178 :   switch(typ(t))
    2129             :   {
    2130             :     case t_INT: case t_REAL:
    2131         994 :       return (signe(t) < 0);
    2132             :     case t_FRAC:
    2133           0 :       return (signe(gel(t,1)) < 0);
    2134             :   }
    2135        2184 :   return 0;
    2136             : }
    2137             : 
    2138             : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2139             : GEN
    2140         343 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2141             : {
    2142             :   pari_sp av, av2, tetpil;
    2143             :   long signh; /* junk */
    2144         343 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2145             :   GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
    2146             : 
    2147         343 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2148         343 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2149         343 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2150         343 :   vx=varn(x);
    2151         343 :   if (!signe(x))
    2152             :   {
    2153          14 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2154           7 :     *U = pol_0(vx); *V = pol_0(vx);
    2155           7 :     return pol_0(vx);
    2156             :   }
    2157         329 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2158         308 :   dx = degpol(x); dy = degpol(y);
    2159         308 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2160         308 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2161             : 
    2162         161 :   av = avma;
    2163         161 :   u = x = primitive_part(x, &cu);
    2164         161 :   v = y = primitive_part(y, &cv);
    2165         161 :   g = h = gen_1; av2 = avma;
    2166         161 :   um1 = gen_1; uze = gen_0;
    2167             :   for(;;)
    2168             :   {
    2169         182 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2170          21 :     if (gc_needed(av2,1))
    2171             :     {
    2172           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2173           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2174             :     }
    2175          21 :   }
    2176         161 :   if (uze != gen_0) {
    2177             :     GEN r;
    2178          42 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2179          42 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2180          42 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2181          42 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2182          42 :     p1 = ginv(content(v));
    2183             :   }
    2184             :   else /* y | x */
    2185             :   {
    2186         119 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2187         119 :     uze = pol_0(vx);
    2188         119 :     p1 = gen_1;
    2189             :   }
    2190         161 :   if (must_negate(v)) p1 = gneg(p1);
    2191         161 :   tetpil = avma;
    2192         161 :   z = RgX_Rg_mul(v,p1);
    2193         161 :   *U = RgX_Rg_mul(uze,p1);
    2194         161 :   *V = RgX_Rg_mul(vze,p1);
    2195         161 :   gptr[0] = &z;
    2196         161 :   gptr[1] = U;
    2197         161 :   gptr[2] = V;
    2198         161 :   gerepilemanysp(av,tetpil,gptr,3); return z;
    2199             : }
    2200             : 
    2201             : int
    2202          56 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2203             : {
    2204          56 :   pari_sp av = avma, av2, tetpil;
    2205             :   long signh; /* junk */
    2206             :   long vx;
    2207             :   GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
    2208             : 
    2209          56 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2210          56 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2211          56 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2212          56 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2213          56 :   if (!signe(T)) {
    2214           0 :     if (degpol(x) <= amax) {
    2215           0 :       *P = RgX_copy(x);
    2216           0 :       *Q = pol_1(varn(x));
    2217           0 :       return 1;
    2218             :     }
    2219           0 :     return 0;
    2220             :   }
    2221          56 :   if (amax+bmax >= degpol(T))
    2222           0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2223             :                     mkvec3(stoi(amax), stoi(bmax), T));
    2224          56 :   vx = varn(T);
    2225          56 :   u = x = primitive_part(x, &cu);
    2226          56 :   v = T = primitive_part(T, &cv);
    2227          56 :   g = h = gen_1; av2 = avma;
    2228          56 :   um1 = gen_1; uze = gen_0;
    2229             :   for(;;)
    2230             :   {
    2231         112 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2232         112 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
    2233         112 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2234          56 :     if (gc_needed(av2,1))
    2235             :     {
    2236           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2237           0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2238             :     }
    2239          56 :   }
    2240          56 :   if (uze == gen_0)
    2241             :   {
    2242           0 :     avma = av; *P = pol_0(vx); *Q = pol_1(vx);
    2243           0 :     return 1;
    2244             :   }
    2245          56 :   if (cu) uze = RgX_Rg_div(uze,cu);
    2246          56 :   p1 = ginv(content(v));
    2247          56 :   if (must_negate(v)) p1 = gneg(p1);
    2248          56 :   tetpil = avma;
    2249          56 :   *P = RgX_Rg_mul(v,p1);
    2250          56 :   *Q = RgX_Rg_mul(uze,p1);
    2251          56 :   gptr[0] = P;
    2252          56 :   gptr[1] = Q;
    2253          56 :   gerepilemanysp(av,tetpil,gptr,2); return 1;
    2254             : }
    2255             : 
    2256             : /*******************************************************************/
    2257             : /*                                                                 */
    2258             : /*                RESULTANT USING DUCOS VARIANT                    */
    2259             : /*                                                                 */
    2260             : /*******************************************************************/
    2261             : /* x^n / y^(n-1), assume n > 0 */
    2262             : static GEN
    2263       64792 : Lazard(GEN x, GEN y, long n)
    2264             : {
    2265             :   long a;
    2266             :   GEN c;
    2267             : 
    2268       64792 :   if (n == 1) return x;
    2269        1320 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    2270        1320 :   c=x; n-=a;
    2271        4576 :   while (a>1)
    2272             :   {
    2273        1936 :     a>>=1; c=gdivexact(gsqr(c),y);
    2274        1936 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    2275             :   }
    2276        1320 :   return c;
    2277             : }
    2278             : 
    2279             : /* F (x/y)^(n-1), assume n >= 1 */
    2280             : static GEN
    2281      113701 : Lazard2(GEN F, GEN x, GEN y, long n)
    2282             : {
    2283      113701 :   if (n == 1) return F;
    2284         693 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    2285             : }
    2286             : 
    2287             : static GEN
    2288      113701 : RgX_neg_i(GEN x, long lx)
    2289             : {
    2290             :   long i;
    2291      113701 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    2292      113701 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    2293      113701 :   return y;
    2294             : }
    2295             : static GEN
    2296      340844 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    2297             : {
    2298             :   long i;
    2299             :   GEN z;
    2300      340844 :   if (isrationalzero(x)) return pol_0(varn(y));
    2301      340830 :   z = cgetg(ly,t_POL); z[1] = y[1];
    2302      340830 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2303      340830 :   return z;
    2304             : }
    2305             : static long
    2306      337974 : reductum_lg(GEN x, long lx)
    2307             : {
    2308      337974 :   long i = lx-2;
    2309      337974 :   while (i > 1 && gequal0(gel(x,i))) i--;
    2310      337974 :   return i+1;
    2311             : }
    2312             : 
    2313             : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    2314             :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    2315             : static GEN
    2316      113701 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    2317             : {
    2318      113701 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
    2319             :   long p, q, j, lP, lQ;
    2320             :   pari_sp av;
    2321             : 
    2322      113701 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    2323      113701 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    2324             :   /* p > q. Very often p - 1 = q */
    2325      113701 :   av = avma;
    2326             :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    2327      113701 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    2328             : 
    2329      113701 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    2330      117306 :   for (j = q+1; j < p; j++)
    2331             :   {
    2332        3605 :     if (degpol(H) == q-1)
    2333             :     { /* h0 = coeff of degree q-1 = leading coeff */
    2334        3185 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    2335        3185 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    2336             :     }
    2337             :     else
    2338         420 :       H = RgX_shift_shallow(H, 1);
    2339        3605 :     if (j+2 < lP)
    2340             :     {
    2341        3038 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    2342        3038 :       A = A? RgX_add(A, TMP): TMP;
    2343             :     }
    2344        3605 :     if (gc_needed(av,1))
    2345             :     {
    2346         147 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    2347         147 :       gerepileall(av,A?2:1,&H,&A);
    2348             :     }
    2349             :   }
    2350      113701 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    2351      113701 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    2352      113701 :   A = A? RgX_add(A, TMP): TMP;
    2353      113701 :   A = RgX_Rg_divexact(A, p0);
    2354      113701 :   if (degpol(H) == q-1)
    2355             :   {
    2356      113386 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    2357      113386 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    2358             :   }
    2359             :   else
    2360         315 :     A = RgX_Rg_mul(addshift(H,A), q0);
    2361      113701 :   return RgX_Rg_divexact(A, s);
    2362             : }
    2363             : 
    2364             : /* Ducos's subresultant */
    2365             : GEN
    2366       65499 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    2367             : {
    2368             :   pari_sp av, av2;
    2369       65499 :   long dP, dQ, delta, sig = 1;
    2370             :   GEN cP, cQ, Z, s;
    2371             : 
    2372       65499 :   dP = degpol(P);
    2373       65499 :   dQ = degpol(Q); delta = dP - dQ;
    2374       65499 :   if (delta < 0)
    2375             :   {
    2376        2408 :     if (both_odd(dP, dQ)) sig = -1;
    2377        2408 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    2378             :   }
    2379       65499 :   if (sol) *sol = gen_0;
    2380       65499 :   av = avma;
    2381       65499 :   if (dQ <= 0)
    2382             :   {
    2383        1400 :     if (dQ < 0) return RgX_get_0(P);
    2384        1400 :     s = gpowgs(gel(Q,2), dP);
    2385        1400 :     if (sig == -1) s = gerepileupto(av, gneg(s));
    2386        1400 :     return s;
    2387             :   }
    2388       64099 :   P = primitive_part(P, &cP);
    2389       64099 :   Q = primitive_part(Q, &cQ);
    2390       64099 :   av2 = avma;
    2391       64099 :   s = gpowgs(leading_coeff(Q),delta);
    2392       64099 :   if (both_odd(dP, dQ)) sig = -sig;
    2393       64099 :   Z = Q;
    2394       64099 :   Q = RgX_pseudorem(P, Q);
    2395       64099 :   P = Z;
    2396      241899 :   while(degpol(Q) > 0)
    2397             :   {
    2398      113701 :     delta = degpol(P) - degpol(Q); /* > 0 */
    2399      113701 :     Z = Lazard2(Q, leading_coeff(Q), s, delta);
    2400      113701 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    2401      113701 :     Q = nextSousResultant(P, Q, Z, s);
    2402      113701 :     P = Z;
    2403      113701 :     if (gc_needed(av,1))
    2404             :     {
    2405          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    2406          13 :       gerepileall(av2,2,&P,&Q);
    2407             :     }
    2408      113701 :     s = leading_coeff(P);
    2409             :   }
    2410       64099 :   if (!signe(Q)) { avma = av; return RgX_get_0(Q); }
    2411       64099 :   s = Lazard(leading_coeff(Q), s, degpol(P));
    2412       64099 :   if (sig == -1) s = gneg(s);
    2413       64099 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    2414       64099 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    2415       64099 :   if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
    2416       63035 :   return gerepilecopy(av, s);
    2417             : }
    2418             : /* Return resultant(P,Q). If sol != NULL: set *sol to the last non-constant
    2419             :  * polynomial in the prs IF the sequence was computed, and gen_0 otherwise.
    2420             :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    2421             :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    2422             :  * in the "generic" case. */
    2423             : GEN
    2424       87681 : resultant_all(GEN P, GEN Q, GEN *sol)
    2425             : {
    2426             :   long TP, TQ;
    2427             :   GEN s;
    2428             : 
    2429       87681 :   if (sol) *sol = gen_0;
    2430       87681 :   if ((s = init_resultant(P,Q))) return s;
    2431       87604 :   if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
    2432           0 :     return resultant2(P,Q); /* inexact */
    2433       87604 :   if (TP && TQ) /* rational */
    2434             :   {
    2435       22413 :     if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
    2436       10482 :     return QX_resultant(P,Q);
    2437             :   }
    2438       65191 :   return RgX_resultant_all(P, Q, sol);
    2439             : }
    2440             : 
    2441             : /*******************************************************************/
    2442             : /*                                                                 */
    2443             : /*               RESULTANT USING SYLVESTER MATRIX                  */
    2444             : /*                                                                 */
    2445             : /*******************************************************************/
    2446             : static GEN
    2447           0 : _pol_0(void)
    2448             : {
    2449           0 :   GEN x = cgetg(3,t_POL);
    2450           0 :   x[1] = 0;
    2451           0 :   gel(x,2) = gen_0; return x;
    2452             : }
    2453             : 
    2454             : static GEN
    2455         133 : sylvester_col(GEN x, long j, long d, long D)
    2456             : {
    2457         133 :   GEN c = cgetg(d+1,t_COL);
    2458             :   long i;
    2459         133 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    2460         133 :   for (   ; i<=D; i++) gel(c,i) = gel(x,D-i+2);
    2461         133 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    2462         133 :   return c;
    2463             : }
    2464             : 
    2465             : GEN
    2466          28 : sylvestermatrix_i(GEN x, GEN y)
    2467             : {
    2468             :   long j,d,dx,dy;
    2469             :   GEN M;
    2470             : 
    2471          28 :   dx = degpol(x); if (dx < 0) { dx = 0; x = _pol_0(); }
    2472          28 :   dy = degpol(y); if (dy < 0) { dy = 0; y = _pol_0(); }
    2473          28 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    2474          28 :   for (j=1; j<=dy; j++) gel(M,j)    = sylvester_col(x,j,d,j+dx);
    2475          28 :   for (j=1; j<=dx; j++) gel(M,j+dy) = sylvester_col(y,j,d,j+dy);
    2476          28 :   return M;
    2477             : }
    2478             : 
    2479             : GEN
    2480           7 : sylvestermatrix(GEN x, GEN y)
    2481             : {
    2482             :   long i,j,lx;
    2483           7 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    2484           7 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    2485           7 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    2486           7 :   x = sylvestermatrix_i(x,y); lx = lg(x);
    2487          28 :   for (i=1; i<lx; i++)
    2488          21 :     for (j=1; j<lx; j++) gcoeff(x,i,j) = gcopy(gcoeff(x,i,j));
    2489           7 :   return x;
    2490             : }
    2491             : 
    2492             : GEN
    2493          21 : resultant2(GEN x, GEN y)
    2494             : {
    2495             :   pari_sp av;
    2496             :   GEN r;
    2497          21 :   if ((r = init_resultant(x,y))) return r;
    2498          21 :   av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
    2499             : }
    2500             : 
    2501             : /* If x a t_POL, let vx = main variable of x; return a t_POL in variable v0:
    2502             :  * if vx <= v, return subst(x, v, pol_x(v0))
    2503             :  * if vx >  v, return scalarpol(x, v0) */
    2504             : static GEN
    2505        2254 : fix_pol(GEN x, long v, long v0)
    2506             : {
    2507             :   long vx;
    2508        2254 :   if (typ(x) != t_POL) return x;
    2509        2254 :   vx = varn(x);
    2510        2254 :   if (v == vx)
    2511             :   {
    2512        2191 :     if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
    2513        2191 :     return x;
    2514             :   }
    2515          63 :   if (varncmp(v, vx) > 0)
    2516             :   {
    2517          56 :     x = gsubst(x,v,pol_x(v0));
    2518          56 :     if (typ(x) == t_POL && varn(x) == v0) return x;
    2519             :   }
    2520          35 :   return scalarpol_shallow(x, v0);
    2521             : }
    2522             : 
    2523             : /* resultant of x and y with respect to variable v, or with respect to their
    2524             :  * main variable if v < 0. */
    2525             : GEN
    2526        1323 : polresultant0(GEN x, GEN y, long v, long flag)
    2527             : {
    2528        1323 :   long v0 = 0;
    2529        1323 :   pari_sp av = avma;
    2530             : 
    2531        1323 :   if (v >= 0)
    2532             :   {
    2533        1113 :     v0 = fetch_var_higher();
    2534        1113 :     x = fix_pol(x,v, v0);
    2535        1113 :     y = fix_pol(y,v, v0);
    2536             :   }
    2537        1323 :   switch(flag)
    2538             :   {
    2539             :     case 2:
    2540        1316 :     case 0: x=resultant_all(x,y,NULL); break;
    2541           7 :     case 1: x=resultant2(x,y); break;
    2542           0 :     default: pari_err_FLAG("polresultant");
    2543             :   }
    2544        1323 :   if (v >= 0) (void)delete_var();
    2545        1323 :   return gerepileupto(av,x);
    2546             : }
    2547             : GEN
    2548         868 : polresultantext0(GEN x, GEN y, long v)
    2549             : {
    2550             :   GEN R, U, V;
    2551         868 :   long v0 = 0;
    2552         868 :   pari_sp av = avma;
    2553             : 
    2554         868 :   if (v >= 0)
    2555             :   {
    2556          14 :     v0 = fetch_var_higher();
    2557          14 :     x = fix_pol(x,v, v0);
    2558          14 :     y = fix_pol(y,v, v0);
    2559             :   }
    2560         868 :   R = subresext_i(x,y, &U,&V);
    2561         868 :   if (v >= 0)
    2562             :   {
    2563          14 :     (void)delete_var();
    2564          14 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    2565          14 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    2566             :   }
    2567         868 :   return gerepilecopy(av, mkvec3(U,V,R));
    2568             : }
    2569             : GEN
    2570         840 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    2571             : 
    2572             : /*******************************************************************/
    2573             : /*                                                                 */
    2574             : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    2575             : /*                                                                 */
    2576             : /*******************************************************************/
    2577             : 
    2578             : /* (v - x)^d */
    2579             : static GEN
    2580         133 : caract_const(pari_sp av, GEN x, long v, long d)
    2581         133 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    2582             : 
    2583             : /* return caract(Mod(x,T)) in variable v */
    2584             : GEN
    2585       59409 : RgXQ_charpoly(GEN x, GEN T, long v)
    2586             : {
    2587       59409 :   pari_sp av = avma;
    2588       59409 :   long d = degpol(T), dx, vx, vp, v0;
    2589             :   GEN ch, L;
    2590             : 
    2591       59409 :   if (typ(x) != t_POL) return caract_const(av, x, v, d);
    2592       59395 :   vx = varn(x);
    2593       59395 :   vp = varn(T);
    2594       59395 :   if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
    2595       59395 :   if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
    2596       59395 :   dx = degpol(x);
    2597       59395 :   if (dx <= 0)
    2598          49 :     return dx? monomial(gen_1, d, v): caract_const(av, gel(x,2), v, d);
    2599             : 
    2600       59346 :   v0 = fetch_var_higher();
    2601       59346 :   x = RgX_neg(x);
    2602       59346 :   gel(x,2) = gadd(gel(x,2), pol_x(v));
    2603       59346 :   setvarn(x, v0);
    2604       59346 :   T = leafcopy(T); setvarn(T, v0);
    2605       59346 :   ch = resultant_all(T, x, NULL);
    2606       59346 :   (void)delete_var();
    2607             :   /* test for silly input: x mod (deg 0 polynomial) */
    2608       59346 :   if (typ(ch) != t_POL) { avma = av; return pol_1(v); }
    2609             : 
    2610       59346 :   L = leading_coeff(ch);
    2611       59346 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    2612       59346 :   return gerepileupto(av, ch);
    2613             : }
    2614             : 
    2615             : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    2616             :  * algebra nf[t]/(Q(t)) */
    2617             : GEN
    2618         224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    2619             : {
    2620         224 :   const char *f = "rnfcharpoly";
    2621         224 :   long dQ = degpol(Q);
    2622         224 :   pari_sp av = avma;
    2623             :   GEN T;
    2624             : 
    2625         224 :   if (v < 0) v = 0;
    2626         224 :   nf = checknf(nf); T = nf_get_pol(nf);
    2627         224 :   Q = RgX_nffix(f, T,Q,0);
    2628         224 :   switch(typ(x))
    2629             :   {
    2630             :     case t_INT:
    2631          28 :     case t_FRAC: return caract_const(av, x, v, dQ);
    2632             :     case t_POLMOD:
    2633          91 :       x = polmod_nffix2(f,T,Q, x,0);
    2634          49 :       break;
    2635             :     case t_POL:
    2636          56 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    2637          42 :       break;
    2638          49 :     default: pari_err_TYPE(f,x);
    2639             :   }
    2640          91 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    2641             :   /* x a t_POL in variable vQ */
    2642          49 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    2643          49 :   if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
    2644          49 :   return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    2645             : }
    2646             : 
    2647             : /*******************************************************************/
    2648             : /*                                                                 */
    2649             : /*                  GCD USING SUBRESULTANT                         */
    2650             : /*                                                                 */
    2651             : /*******************************************************************/
    2652             : static int inexact(GEN x, int *simple, int *rational);
    2653             : static int
    2654     6946858 : isinexactall(GEN x, int *simple, int *rational)
    2655             : {
    2656     6946858 :   long i, lx = lg(x);
    2657    32851176 :   for (i=2; i<lx; i++)
    2658    25904318 :     if (inexact(gel(x,i), simple, rational)) return 1;
    2659     6946858 :   return 0;
    2660             : }
    2661             : /* return 1 if coeff explosion is not possible */
    2662             : static int
    2663    25904374 : inexact(GEN x, int *simple, int *rational)
    2664             : {
    2665    25904374 :   int junk = 0;
    2666    25904374 :   switch(typ(x))
    2667             :   {
    2668    25864075 :     case t_INT: case t_FRAC: return 0;
    2669             : 
    2670           0 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    2671             : 
    2672             :     case t_INTMOD:
    2673             :     case t_FFELT:
    2674        7287 :       *rational = 0;
    2675        7287 :       if (!*simple) *simple = 1;
    2676        7287 :       return 0;
    2677             : 
    2678             :     case t_COMPLEX:
    2679           0 :       *rational = 0;
    2680           0 :       return inexact(gel(x,1), simple, rational)
    2681           0 :           || inexact(gel(x,2), simple, rational);
    2682             :     case t_QUAD:
    2683           0 :       *rational = *simple = 0;
    2684           0 :       return inexact(gel(x,2), &junk, rational)
    2685           0 :           || inexact(gel(x,3), &junk, rational);
    2686             : 
    2687             :     case t_POLMOD:
    2688        3787 :       *rational = 0;
    2689        3787 :       return isinexactall(gel(x,1), simple, rational);
    2690             :     case t_POL:
    2691       29197 :       *rational = 0;
    2692       29197 :       *simple = -1;
    2693       29197 :       return isinexactall(x, &junk, rational);
    2694             :     case t_RFRAC:
    2695          28 :       *rational = 0;
    2696          28 :       *simple = -1;
    2697          56 :       return inexact(gel(x,1), &junk, rational)
    2698          28 :           || inexact(gel(x,2), &junk, rational);
    2699             :   }
    2700           0 :   *rational = 0;
    2701           0 :   *simple = -1; return 0;
    2702             : }
    2703             : 
    2704             : /* x monomial, y t_POL in the same variable */
    2705             : static GEN
    2706     7153646 : gcdmonome(GEN x, GEN y)
    2707             : {
    2708     7153646 :   pari_sp av = avma;
    2709     7153646 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    2710     7153646 :   long i, l = lg(y);
    2711     7153646 :   GEN t, v = cgetg(l, t_VEC);
    2712     7153646 :   gel(v,1) = gel(x,dx+2);
    2713     7153646 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    2714     7153646 :   t = content(v); /* gcd(lc(x), cont(y)) */
    2715     7153646 :   t = simplify_shallow(t);
    2716     7153646 :   if (dx < e) e = dx;
    2717     7153646 :   return gerepileupto(av, monomialcopy(t, e, varn(x)));
    2718             : }
    2719             : 
    2720             : /* x, y are t_POL in the same variable */
    2721             : GEN
    2722    10610926 : RgX_gcd(GEN x, GEN y)
    2723             : {
    2724             :   long dx, dy;
    2725             :   pari_sp av, av1;
    2726             :   GEN d, g, h, p1, p2, u, v;
    2727    10610926 :   int simple = 0, rational = 1;
    2728             : 
    2729    10610926 :   if (isexactzero(y)) return RgX_copy(x);
    2730    10610590 :   if (isexactzero(x)) return RgX_copy(y);
    2731    10610583 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    2732     4083174 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    2733     3456937 :   if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
    2734             :   {
    2735           0 :     av = avma; u = ggcd(content(x), content(y));
    2736           0 :     return gerepileupto(av, scalarpol(u, varn(x)));
    2737             :   }
    2738     3456937 :   if (rational) return QX_gcd(x,y); /* Q[X] */
    2739             : 
    2740        4137 :   av = avma;
    2741        4137 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    2742             :   else
    2743             :   {
    2744        3360 :     dx = lg(x); dy = lg(y);
    2745        3360 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    2746        3360 :     if (dy==3)
    2747             :     {
    2748           0 :       d = ggcd(gel(y,2), content(x));
    2749           0 :       return gerepileupto(av, scalarpol(d, varn(x)));
    2750             :     }
    2751        3360 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    2752        3360 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    2753        3360 :     d = ggcd(p1,p2);
    2754        3360 :     av1 = avma;
    2755        3360 :     g = h = gen_1;
    2756             :     for(;;)
    2757             :     {
    2758        4494 :       GEN r = RgX_pseudorem(u,v);
    2759        4494 :       long degq, du, dv, dr = lg(r);
    2760             : 
    2761        4494 :       if (!signe(r)) break;
    2762        2310 :       if (dr <= 3)
    2763             :       {
    2764        1176 :         avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
    2765             :       }
    2766        1134 :       if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
    2767        1134 :       du = lg(u); dv = lg(v); degq = du-dv;
    2768        1134 :       u = v; p1 = g; g = leading_coeff(u);
    2769        1134 :       switch(degq)
    2770             :       {
    2771         189 :         case 0: break;
    2772             :         case 1:
    2773         763 :           p1 = gmul(h,p1); h = g; break;
    2774             :         default:
    2775         182 :           p1 = gmul(gpowgs(h,degq), p1);
    2776         182 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    2777             :       }
    2778        1134 :       v = RgX_Rg_div(r,p1);
    2779        1134 :       if (gc_needed(av1,1))
    2780             :       {
    2781           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
    2782           0 :         gerepileall(av1,4, &u,&v,&g,&h);
    2783             :       }
    2784        1134 :     }
    2785        2184 :     x = RgX_Rg_mul(primpart(v), d);
    2786             :   }
    2787        2961 :   if (must_negate(x)) x = RgX_neg(x);
    2788        2961 :   return gerepileupto(av,x);
    2789             : }
    2790             : 
    2791             : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
    2792             : static GEN
    2793        4235 : RgX_disc_aux(GEN P)
    2794             : {
    2795        4235 :   long n = degpol(P), TP, dd;
    2796             :   GEN D, L, y, p;
    2797        4235 :   if (!signe(P) || !n) return RgX_get_0(P);
    2798        4235 :   if (n == 1) return RgX_get_1(P);
    2799        4130 :   if (n == 2) {
    2800         651 :     GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
    2801         651 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    2802             :   }
    2803        3479 :   TP = RgX_simpletype(P);
    2804        3479 :   if (TP == t_INT) return ZX_disc(P);
    2805         217 :   if (TP == t_FRAC) return QX_disc(P);
    2806         217 :   p = NULL;
    2807         217 :   if (RgX_is_FpX(P, &p) && p)
    2808          28 :     return Fp_to_mod(FpX_disc(RgX_to_FpX(P,p), p), p);
    2809             : 
    2810         189 :   y = RgX_deriv(P);
    2811         189 :   if (!signe(y)) return RgX_get_0(y);
    2812         189 :   dd = degpol(P)-2 - degpol(y);
    2813         189 :   if (TP == t_REAL)
    2814          14 :     D = resultant2(P,y);
    2815             :   else
    2816             :   {
    2817         175 :     D = RgX_resultant_all(P, y, NULL);
    2818         175 :     if (D == gen_0) return RgX_get_0(y);
    2819             :   }
    2820         189 :   L = leading_coeff(P);
    2821         189 :   if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
    2822         189 :   if (n & 2) D = gneg(D);
    2823         189 :   return D;
    2824             : }
    2825             : GEN
    2826        1176 : RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }
    2827             : 
    2828             : GEN
    2829        3073 : poldisc0(GEN x, long v)
    2830             : {
    2831             :   pari_sp av;
    2832        3073 :   switch(typ(x))
    2833             :   {
    2834             :     case t_POL:
    2835             :     {
    2836             :       GEN D;
    2837        3059 :       long v0 = -1;
    2838        3059 :       av = avma;
    2839        3059 :       if (v >= 0 && v != varn(x))
    2840             :       {
    2841           0 :         v0 = fetch_var_higher();
    2842           0 :         x = fix_pol(x,v, v0);
    2843             :       }
    2844        3059 :       D = RgX_disc_aux(x);
    2845        3059 :       if (v0 >= 0) (void)delete_var();
    2846        3059 :       return gerepileupto(av, D);
    2847             :     }
    2848             : 
    2849             :     case t_COMPLEX:
    2850           0 :       return utoineg(4);
    2851             : 
    2852             :     case t_QUAD:
    2853           0 :       return quad_disc(x);
    2854             :     case t_POLMOD:
    2855           0 :       return poldisc0(gel(x,1), v);
    2856             : 
    2857             :     case t_QFR: case t_QFI:
    2858          14 :       av = avma; return gerepileuptoint(av, qfb_disc(x));
    2859             : 
    2860             :     case t_VEC: case t_COL: case t_MAT:
    2861             :     {
    2862             :       long i;
    2863           0 :       GEN z = cgetg_copy(x, &i);
    2864           0 :       for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
    2865           0 :       return z;
    2866             :     }
    2867             :   }
    2868           0 :   pari_err_TYPE("poldisc",x);
    2869           0 :   return NULL; /* not reached */
    2870             : }
    2871             : 
    2872             : GEN
    2873           7 : reduceddiscsmith(GEN x)
    2874             : {
    2875           7 :   long j, n = degpol(x);
    2876           7 :   pari_sp av = avma;
    2877             :   GEN xp, M;
    2878             : 
    2879           7 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    2880           7 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    2881           7 :   RgX_check_ZX(x,"poldiscreduced");
    2882           7 :   if (!gequal1(gel(x,n+2)))
    2883           0 :     pari_err_IMPL("non-monic polynomial in poldiscreduced");
    2884           7 :   M = cgetg(n+1,t_MAT);
    2885           7 :   xp = ZX_deriv(x);
    2886          28 :   for (j=1; j<=n; j++)
    2887             :   {
    2888          21 :     gel(M,j) = RgX_to_RgC(xp, n);
    2889          21 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    2890             :   }
    2891           7 :   return gerepileupto(av, ZM_snf(M));
    2892             : }
    2893             : 
    2894             : /***********************************************************************/
    2895             : /**                                                                   **/
    2896             : /**                       STURM ALGORITHM                             **/
    2897             : /**              (number of real roots of x in [a,b])                 **/
    2898             : /**                                                                   **/
    2899             : /***********************************************************************/
    2900             : static int
    2901        1603 : exact_sturm(GEN a)
    2902             : {
    2903        1603 :   switch(typ(a))
    2904             :   {
    2905        1596 :     case t_INT: case t_FRAC: case t_INFINITY: return 1;
    2906           7 :     default: return 0;
    2907             :   }
    2908             : }
    2909             : 
    2910             : /* Deprecated: support the old format: if a (resp. b) is NULL, set it
    2911             :  * to -oo resp. +oo). ZX_sturmpart() should be preferred  */
    2912             : static long
    2913         833 : sturmpart_i(GEN x, GEN a, GEN b)
    2914             : {
    2915             :   long sl, sr, s, t, r1;
    2916         833 :   pari_sp av = avma;
    2917             :   int integral;
    2918             :   GEN g,h,u,v;
    2919             : 
    2920         833 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    2921         833 :   t = typ(x);
    2922         833 :   if (t != t_POL)
    2923             :   {
    2924           0 :     if (is_real_t(t)) return 0;
    2925           0 :     pari_err_TYPE("sturm",x);
    2926             :   }
    2927         833 :   s = lg(x); if (s==3) return 0;
    2928         833 :   u = primpart(x);
    2929         833 :   integral = RgX_is_ZX(u);
    2930         833 :   if (integral && !ZX_is_squarefree(u))
    2931           7 :     pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,u);
    2932         826 :   if (!b && a && typ(a) == t_VEC && lg(a) == 3)
    2933             :   { /* new format */
    2934         126 :     if (integral && exact_sturm(gel(a,1)) && exact_sturm(gel(a,2)))
    2935         126 :       return ZX_sturmpart(u, a);
    2936             :     /* but can't use new function; convert to old form */
    2937           0 :     b = gel(a,2);
    2938           0 :     if (typ(b) == t_INFINITY)
    2939             :     {
    2940           0 :       if (inf_get_sign(b) < 0) return 0;
    2941           0 :       b = NULL;
    2942             :     }
    2943           0 :     a = gel(a,1);
    2944           0 :     if (typ(a) == t_INFINITY)
    2945             :     {
    2946           0 :       if (inf_get_sign(a) > 0) return 0;
    2947           0 :       a = NULL;
    2948             :     }
    2949             :   }
    2950         700 :   else if (integral)
    2951             :   {
    2952         679 :     GEN A = a? a: mkmoo();
    2953         679 :     GEN B = b? b: mkoo();
    2954         679 :     if (exact_sturm(A) && exact_sturm(B)) return ZX_sturmpart(u, mkvec2(A,B));
    2955             :   }
    2956             :   /* legacy code: should only be used if we have a t_REAL somewhere; and even
    2957             :    * then, the calling program should be changed */
    2958          28 :   sl = gsigne(leading_coeff(u));
    2959          28 :   t = a? gsigne(poleval(u,a)): (odd(s)? sl: -sl);
    2960          28 :   if (s==4)
    2961             :   {
    2962           0 :     if (t == 0) return 1;
    2963           0 :     s = b? gsigne(poleval(u,b)):  sl;
    2964           0 :     return (s == t)? 0: 1;
    2965             :   }
    2966          28 :   s = b? gsigne(poleval(u,b)): sl;
    2967          28 :   r1= (t == 0)? 1: 0;
    2968          28 :   v = primpart(RgX_deriv(x));
    2969          28 :   sr = b? gsigne(poleval(v,b)): s;
    2970          28 :   if (sr)
    2971             :   {
    2972          28 :     if (!s) s=sr;
    2973          21 :     else if (sr!=s) { s= -s; r1--; }
    2974             :   }
    2975          28 :   sr = a? gsigne(poleval(v,a)): -t;
    2976          28 :   if (sr)
    2977             :   {
    2978          28 :     if (!t) t=sr;
    2979          21 :     else if (sr!=t) { t= -t; r1++; }
    2980             :   }
    2981          28 :   g=gen_1; h=gen_1;
    2982             :   for(;;)
    2983             :   {
    2984         147 :     GEN p1, r = RgX_pseudorem(u,v);
    2985         147 :     long du=lg(u), dv=lg(v), dr=lg(r), degq=du-dv;
    2986             : 
    2987         147 :     if (dr<=2) pari_err_PREC("polsturm");
    2988         140 :     if (gsigne(leading_coeff(v)) > 0 || degq&1) r=gneg_i(r);
    2989         140 :     sl = gsigne(gel(r,dr-1));
    2990         140 :     sr = b? gsigne(poleval(r,b)): sl;
    2991         140 :     if (sr)
    2992             :     {
    2993         140 :       if (!s) s=sr;
    2994         140 :       else if (sr!=s) { s= -s; r1--; }
    2995             :     }
    2996         140 :     sr = a? gsigne(poleval(r,a)): ((dr&1)? sl: -sl);
    2997         140 :     if (sr)
    2998             :     {
    2999         133 :       if (!t) t=sr;
    3000         133 :       else if (sr!=t) { t= -t; r1++; }
    3001             :     }
    3002         140 :     if (dr==3) return r1;
    3003             : 
    3004         119 :     u=v; p1 = g; g = gabs(leading_coeff(u),DEFAULTPREC);
    3005         119 :     switch(degq)
    3006             :     {
    3007           0 :       case 0: break;
    3008             :       case 1:
    3009         119 :         p1 = gmul(h,p1); h = g; break;
    3010             :       default:
    3011           0 :         p1 = gmul(gpowgs(h,degq),p1);
    3012           0 :         h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
    3013             :     }
    3014         119 :     v = RgX_Rg_divexact(r,p1);
    3015         119 :     if (gc_needed(av,1))
    3016             :     {
    3017           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"polsturm, dr = %ld",dr);
    3018           0 :       gerepileall(av,4,&u,&v,&g,&h);
    3019             :     }
    3020         119 :   }
    3021             : }
    3022             : long
    3023         833 : sturmpart(GEN x, GEN a, GEN b)
    3024             : {
    3025         833 :   pari_sp av = avma;
    3026         833 :   long r = sturmpart_i(x,a,b);
    3027         819 :   avma = av; return r;
    3028             : }
    3029             : long
    3030           0 : RgX_sturmpart(GEN x, GEN ab) { return sturmpart(x, ab, NULL); }
    3031             : 
    3032             : /***********************************************************************/
    3033             : /**                                                                   **/
    3034             : /**                        GENERIC EXTENDED GCD                       **/
    3035             : /**                                                                   **/
    3036             : /***********************************************************************/
    3037             : /* assume typ(x) = typ(y) = t_POL */
    3038             : GEN
    3039       10805 : RgXQ_inv(GEN x, GEN y)
    3040             : {
    3041       10805 :   long vx=varn(x), vy=varn(y);
    3042             :   pari_sp av;
    3043             :   GEN u, v, d;
    3044             : 
    3045       21610 :   while (vx != vy)
    3046             :   {
    3047           0 :     if (varncmp(vx,vy) > 0)
    3048             :     {
    3049           0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3050           0 :       return scalarpol(d, vy);
    3051             :     }
    3052           0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3053           0 :     x = gel(x,2); vx = gvar(x);
    3054             :   }
    3055       10805 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3056       10805 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3057       10798 :   d = gdiv(u,d);
    3058       10798 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3059       10798 :   return gerepileupto(av, d);
    3060             : }
    3061             : 
    3062             : /*Assume x is a polynomial and y is not */
    3063             : static GEN
    3064         112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3065             : {
    3066         112 :   long vx = varn(x);
    3067         112 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    3068         112 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    3069          84 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    3070          56 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    3071             : }
    3072             : /* Assume x==0, y!=0 */
    3073             : static GEN
    3074          63 : zero_bezout(GEN y, GEN *U, GEN *V)
    3075             : {
    3076          63 :   *U=gen_0; *V = ginv(y); return gen_1;
    3077             : }
    3078             : 
    3079             : GEN
    3080         280 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    3081             : {
    3082         280 :   long tx=typ(x), ty=typ(y), vx;
    3083         280 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    3084         245 :   if (tx != t_POL)
    3085             :   {
    3086         140 :     if (ty == t_POL)
    3087          56 :       return scalar_bezout(y,x,v,u);
    3088             :     else
    3089             :     {
    3090          84 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    3091          84 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    3092          63 :       if (xis0) return zero_bezout(y,u,v);
    3093          42 :       else return zero_bezout(x,v,u);
    3094             :     }
    3095             :   }
    3096         105 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    3097          49 :   vx = varn(x);
    3098          49 :   if (vx != varn(y))
    3099           0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    3100           0 :                                    : scalar_bezout(y,x,v,u);
    3101          49 :   return RgX_extgcd(x,y,u,v);
    3102             : }
    3103             : 
    3104             : GEN
    3105         280 : gcdext0(GEN x, GEN y)
    3106             : {
    3107         280 :   GEN z=cgetg(4,t_VEC);
    3108         280 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    3109         280 :   return z;
    3110             : }
    3111             : 
    3112             : /*******************************************************************/
    3113             : /*                                                                 */
    3114             : /*                    GENERIC (modular) INVERSE                    */
    3115             : /*                                                                 */
    3116             : /*******************************************************************/
    3117             : 
    3118             : GEN
    3119        4669 : ginvmod(GEN x, GEN y)
    3120             : {
    3121        4669 :   long tx=typ(x);
    3122             : 
    3123        4669 :   switch(typ(y))
    3124             :   {
    3125             :     case t_POL:
    3126        4669 :       if (tx==t_POL) return RgXQ_inv(x,y);
    3127        3941 :       if (is_scalar_t(tx)) return ginv(x);
    3128           0 :       break;
    3129             :     case t_INT:
    3130           0 :       if (tx==t_INT) return Fp_inv(x,y);
    3131           0 :       if (tx==t_POL) return gen_0;
    3132             :   }
    3133           0 :   pari_err_TYPE2("ginvmod",x,y);
    3134           0 :   return NULL; /* not reached */
    3135             : }
    3136             : 
    3137             : /***********************************************************************/
    3138             : /**                                                                   **/
    3139             : /**                         NEWTON POLYGON                            **/
    3140             : /**                                                                   **/
    3141             : /***********************************************************************/
    3142             : 
    3143             : /* assume leading coeff of x is non-zero */
    3144             : GEN
    3145          28 : newtonpoly(GEN x, GEN p)
    3146             : {
    3147             :   GEN y;
    3148             :   long n,ind,a,b,c,u1,u2,r1,r2;
    3149          28 :   long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
    3150             : 
    3151          28 :   if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
    3152          28 :   n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
    3153          28 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    3154          28 :   vval = (long *) pari_malloc(sizeof(long)*(n+1));
    3155          28 :   for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
    3156          42 :   for (a=0, ind=1; a<n; a++)
    3157             :   {
    3158          42 :     if (vval[a] != LONG_MAX) break;
    3159          14 :     gel(y,ind++) = mkoo();
    3160             :   }
    3161          84 :   for (b=a+1; b<=n; a=b, b=a+1)
    3162             :   {
    3163          56 :     while (vval[b] == LONG_MAX) b++;
    3164          56 :     u1 = vval[a]-vval[b];
    3165          56 :     u2 = b-a;
    3166         154 :     for (c=b+1; c<=n; c++)
    3167             :     {
    3168          98 :       if (vval[c] == LONG_MAX) continue;
    3169          70 :       r1 = vval[a]-vval[c];
    3170          70 :       r2 = c-a;
    3171          70 :       if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
    3172             :     }
    3173          56 :     while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
    3174             :   }
    3175          28 :   pari_free(vval); return y;
    3176             : }

Generated by: LCOV version 1.11