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 21196-f12677d) Lines: 1691 1867 90.6 %
Date: 2017-10-22 06:23:24 Functions: 145 152 95.4 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11