Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - polarit2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 20777-d2a9243) Lines: 1523 1705 89.3 %
Date: 2017-06-25 05:59:24 Functions: 132 139 95.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11