/***********************************************************************/
/**                                                                   **/
/**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
/**                         (second part)                             **/
/**                                                                   **/
/***********************************************************************/
/* $Id: polarit2.c,v 1.2 1999/09/20 13:30:05 karim Exp $ */
#include "pari.h"

GEN
polsym(GEN x, long n)
{
  long av1,av2,dx=lgef(x)-3,i,k;
  GEN s,y,x_lead;

  if (n<0) err(impl,"polsym of a negative n");
  if (typ(x) != t_POL) err(typeer,"polsym");
  if (!signe(x)) err(zeropoler,"polsym");
  y=cgetg(n+2,t_COL); y[1]=lstoi(dx);
  x_lead=(GEN)x[dx+2]; if (gcmp1(x_lead)) x_lead=NULL;
  for (k=1; k<=n; k++)
  {
    av1=avma; s = (dx>=k)? gmulsg(k,(GEN)x[dx+2-k]): gzero;
    for (i=1; i<k && i<=dx; i++)
      s = gadd(s,gmul((GEN)y[k-i+1],(GEN)x[dx+2-i]));
    if (x_lead) s = gdiv(s,x_lead);
    av2=avma; y[k+1]=lpile(av1,av2,gneg(s));
  }
  return y;
}

static int (*polcmp_coeff_cmp)(GEN,GEN);

/* assume x and y are polynomials in the same variable whose coeffs can be
 * compared (used to sort polynomial factorizations)
 */
static int
polcmp(GEN x, GEN y)
{
  long i, lx = lgef(x), ly = lgef(y);
  int fl;
#if 0 /* for efficiency */
  if (typ(x) != t_POL || typ(y) != t_POL)
    err(talker,"not a polynomial in polcmp");
#endif
  if (lx > ly) return  1;
  if (lx < ly) return -1;
  for (i=lx-1; i>1; i--)
    if ((fl = polcmp_coeff_cmp((GEN)x[i], (GEN)y[i]))) return fl;
  return 0;
}

/* sort factorization so that factors appear by decreasing degree, in place */
GEN
sort_factor(GEN y, int (*cmp)(GEN,GEN))
{
  GEN w,c1,c2, p1 = (GEN)y[1], p2 = (GEN)y[2];
  long av=avma, n=lg(p1), i;
  int (*old)(GEN,GEN) = polcmp_coeff_cmp;

  c1 = new_chunk(n); c2 = new_chunk(n);
  polcmp_coeff_cmp = cmp;
  w = gen_sort(p1, cmp_IND | cmp_C, polcmp);
  for (i=1; i<n; i++) { c1[i]=p1[w[i]]; c2[i]=p2[w[i]]; }
  for (i=1; i<n; i++) { p1[i]=c1[i]; p2[i]=c2[i]; }
  polcmp_coeff_cmp = old; avma = av; return y;
}

/* for internal use */
GEN
centermod(GEN x, GEN p)
{
  long av,i,lx;
  GEN y,p1,ps2;

  ps2=shifti(p,-1);
  switch(typ(x))
  {
    case t_INT:
      y=modii(x,p);
      if (cmpii(y,ps2)>0) return subii(y,p);
      return y;

    case t_POL: lx=lgef(x);
      y=cgetg(lx,t_POL); y[1]=x[1];
      for (i=2; i<lx; i++)
      {
	av=avma; p1=modii((GEN)x[i],p);
	if (cmpii(p1,ps2)>0) p1=subii(p1,p);
	y[i]=lpileupto(av,p1);
      }
      return y;

    case t_COL: lx=lg(x);
      y=cgetg(lx,t_COL);
      for (i=1; i<lx; i++)
      {
	p1=modii((GEN)x[i],p);
	if (cmpii(p1,ps2)>0) p1=subii(p1,p);
	y[i]=(long)p1;
      }
      return y;
  }
  return x;
}

static GEN
decpol(GEN x, long klim)
{
  short int pos[200];
  long av=avma,av1,tete,k,kin,i,j,i1,i2,fl,d,nbfact;
  GEN res,p1,p2;

  kin=1; res=cgetg(lgef(x)-2,t_VEC); nbfact=0;
  p1=roots(x,DEFAULTPREC); d=lg(p1)-1; if (!klim) klim=d;
  do
  {
    fl=1;
    for (k=kin; k+k <= d && k<=klim; k++)
    {
      for (i=0; i<=k; i++) pos[i]=i;
      do
      {
	av1=avma; p2=gzero; j=0;
	for (i1=1; i1<=k; i1++) p2=gadd(p2,(GEN)p1[pos[i1]]);
	if (gexpo(gimag(p2))<-20 && gexpo(gsub(p2,ground(p2)))<-20)
	{
	  p2=gun;
	  for (i1=1; i1<=k; i1++)
	    p2=gmul(p2,gsub(polx[0],(GEN)p1[pos[i1]]));
          p2 = ground(p2);
	  if (gcmp0(gimag(p2)) && gcmp0(gmod(x,p2)))
	  {
	    res[++nbfact]=(long)p2; x=gdiv(x,p2);
            kin=k; p2=cgetg(d-k+1,t_COL);
            for (i=i1=i2=1; i<=d; i++)
	    {
	      if (i1<=k && i==pos[i1]) i1++;
	      else p2[i2++]=p1[i];
	    }
	    p1=p2; d-=k; fl=0; break;
	  }
	}
	avma=av1; pos[k]++;
	while (pos[k-j] > d-j) { j++; pos[k-j]++; }
	for (i=k-j+1; i<=k; i++) pos[i]=i+pos[k-j]-k+j;
      }
      while (j<k);
      if (!fl) break;
    }
    if (lgef(x)<=3) break;
  }
  while (!fl || (k+k <= d && k<=klim));
  if (lgef(x)>3) res[++nbfact]=(long)x;
  setlg(res,nbfact+1); tete=avma;
  return gerepile(av,tete,greal(res));
}

/* This code originally by Richard Schroeppel */

/* Note that PARI's idea of the maximum possible coefficient involves the
 * limit on the degree (klim).  Consider revising this.  If I don't respect
 * the degree limit when testing potential factors, there's the possibility
 * that I might identify a high degree factor that isn't irreducible, because
 * it's lower degree divisors were missed because they had a coefficient
 * outside the Borne limit for klim, but the higher degree factor had it's
 * coefficients within Borne.  This would still have the property that any
 * factors of degree <= klim were guaranteed irr, but higher degrees
 * (> 2*klim) might not be irr.
 *
 * The subroutine:
 *  fxn points at the first unconsidered factor for the current combination
 *  psf is the product-so-far, or 0 for a null product
 *  dlim is the degree limit remaining for unconsidered divisors
 *  other arguments are "global" and must already be setup
 *  as factors are found, they are put in cmbf_fax; the count is kept in
 *  cmbf_faxn; and they are divided out of cmbf_target; the degree and
 *  leading coefficient are updated; and the constituent modular factors
 *  are deleted from cmbf_modfax.
 *  exit value is 1 if any factors are found.
 *  If psf is 0, all factors made up from pieces at or after fxn will be
 *  found & removed.  If psf is not 0, only the factor which is a
 *  continuation of psf will be found.
 */
/* setup for calling cmbf = combine_factors */
static GEN cmbf_target;      /* target poly. Assume content removed */
static GEN cmbf_lc;          /* leading coefficient */
static GEN cmbf_abslc;       /* |lc| */
static GEN cmbf_abslcxtarget;/* abslc * target */
static GEN cmbf_mod;         /* modulus */
static GEN cmbf_fax;         /* result array + one cell for leftover cst. */
static GEN cmbf_modfax; /* array of modular factors.  Each has LC 1.1 based
                           indexing. Product should be congruent to a/lc(a). */
static long cmbf_degree;     /* degree(target) */
static long cmbf_modfaxn;    /* number of modular factors */
static long cmbf_faxn;       /* last used cell in fax. # of factors found */
static GEN cmbf_modhalf;

/* assume same variables, y normalized, non constant */
static GEN
polidivis(GEN x, GEN y, GEN bound)
{
  long dx,dy,dz,i,j,av, vx = varn(x);
  GEN z,p1,y_lead;

  dy=lgef(y)-3;
  dx=lgef(x)-3;
  dz=dx-dy; if (dz<0) return NULL;
  z=cgetg(dz+3,t_POL);
  x += 2; y += 2; z += 2;
  y_lead = (GEN)y[dy];
  if (gcmp1(y_lead)) y_lead = NULL;

  p1 = (GEN)x[dx];
  z[dz]=y_lead? ldivii(p1,y_lead): licopy(p1);
  for (i=dx-1; i>=dy; i--)
  {
    av=avma; p1=(GEN)x[i];
    for (j=i-dy+1; j<=i && j<=dz; j++)
      p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
    if (y_lead) { p1 = gdiv(p1,y_lead); if (typ(p1)!=t_INT) return NULL; }
    if (absi_cmp(p1, bound) > 0) return NULL;
    p1 = gerepileupto(av,p1);
    z[i-dy] = (long)p1;
  }
  av = avma;
  for (;; i--)
  {
    p1 = (GEN)x[i];
    /* we always enter this loop at least once */
    for (j=0; j<=i && j<=dz; j++)
      p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
    if (!gcmp0(p1)) return NULL;
    avma = av;
    if (!i) break;
  }
  z -= 2;
  z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
  return z;
}

static int
try_div(GEN x, GEN y, GEN *q)
{
  if (resii((GEN)x[2], (GEN)y[2]) != gzero)
  {
    if (DEBUGLEVEL>6) fprintferr(".");
    return 0;
  }
  *q = polidivis(x,y,cmbf_modhalf);
  if (*q) return 1;
  if (DEBUGLEVEL>6) fprintferr("*");
  return 0;
}

static int
cmbf(long fxn, GEN psf, long dlim, long hint)
{
  long newd,av,val2,val=0; /* assume failure */
  GEN newf,newpsf,quo,cont;

  if (dlim <= 0 || fxn > cmbf_modfaxn) return 0;

  /* first, try deeper factors without considering the current one */
  if (fxn < cmbf_modfaxn)
  {
    val = cmbf(fxn+1,psf,dlim,hint);
    if (val && psf) return val;
  }

  /* second, try including the current modular factor in the product */
  newf = (GEN)cmbf_modfax[fxn];
  if (!newf) return val; /* modular factor already used */

  newd = lgef(newf)-3; if (newd > dlim) return val;
  av = avma;
  if (!(newd%hint))
  {
    if (!psf) psf = cmbf_abslc;
    newpsf = centermod(gmul(psf,newf),cmbf_mod);

    /* try out the new combination */
    if (try_div(cmbf_abslcxtarget, newpsf, &quo))
    {
      /* found a factor */
      cont = content(newpsf);
      if (signe(leading_term(newpsf)) < 0) cont = negi(cont);
      newpsf = gdiv(newpsf,cont);
      cmbf_fax[++cmbf_faxn] = (long)newpsf; /* store factor */
      cmbf_modfax[fxn] = 0; /* remove used modular factor */
      if (DEBUGLEVEL > 2)
      { fprintferr("\n"); msgtimer("to find factor %Z",newpsf); }
      /* fix up target */
      cmbf_target = gdiv(quo,leading_term(newpsf));
      cmbf_degree = lgef(cmbf_target)-3;
      cmbf_lc = (GEN)cmbf_target[cmbf_degree+2];
      cmbf_abslc = absi(cmbf_lc);
      cmbf_abslcxtarget = gmul(cmbf_abslc,cmbf_target);
      return 1;
    }
  }
  /* degree too large, or no more modular factors? */
  if (newd == dlim || fxn == cmbf_modfaxn) return val;
  /* include more factors */
  val2 = cmbf(fxn+1,newpsf,dlim-newd,hint);
  if (val2) cmbf_modfax[fxn] = 0; /* remove used modular factor */
  else avma = av;
  return val || val2;
}

static long
min_deg(long jmax,ulong tabbit[])
{
  long j,k=1,j1=2;

  for (j=jmax; j>=0; j--)
  {
    for (  ; k<=15; k++)
    {
      if (tabbit[j]&j1) return ((jmax-j)<<4)+k;
      j1<<=1;
    }
    k=0; j1=1;
  }
  return (jmax<<4)+15;
}

static void
dotab(long d,long jmax,ulong *tabkbit,ulong *tmp)
{
  ulong rem,pro;
  long j,d1,d2;

  d1=d>>4; d2=d&15; rem=0;
  for (j=jmax-d1; j>=0; j--)
  {
    pro = tabkbit[j+d1]<<d2;
    tmp[j] = (pro&0xffff)+rem; rem=pro>>16;
  }
  for (j=jmax-d1; j>=0; j--) tabkbit[j] |= tmp[j];
}

/* lift factorisation mod p to mod p^e = pev */
GEN
hensel_lift_fact(GEN pol, GEN fact, GEN p, GEN pev, long e)
{
  long ev,i, nf = lg(fact);
  GEN aprov = pol, res = cgetg(nf, t_VEC);

  for (i=nf-1; i; i--)
  {
    GEN ae,be,u,v,ae2,be2,s,t,pe,pe2,z,g;
    long ltop = avma;

    ae = (GEN)fact[i]; /* lead coeff(ae) = 1 */
    be = Fp_poldivres(aprov,ae,p, NULL);
    g = (GEN)Fp_pol_extgcd(ae,be,p,&u,&v)[2]; /* deg g = 0 */
    if (!gcmp1(g))
    {
      g = mpinvmod(g, p);
      u = gmul(u, g);
      v = gmul(v, g);
    }
    for(pe=p,ev=1;;)
    {
      ev<<=1; pe2=(ev>=e)? pev: sqri(pe);
      g = gadd(aprov, gneg_i(gmul(ae,be)));

      g = Fp_pol_red(g, pe2); g = gdivexact(g, pe);
      z = Fp_pol_red(gmul(v,g), pe);
      t = Fp_poldivres(z,ae,pe, &s);
      t = gadd(gmul(u,g), gmul(t,be));
      t = Fp_pol_red(t, pe);
      be2 = gadd(be, gmul(t,pe));
      ae2 = gadd(ae, gmul(s,pe)); /* already reduced mod pe2 */
      if (ev>=e) break;

      g = gadd(gun, gneg_i(gadd(gmul(u,ae2),gmul(v,be2))));

      g = Fp_pol_red(g, pe2); g = gdivexact(g, pe);
      z = Fp_pol_red(gmul(v,g), pe);
      t = Fp_poldivres(z,ae,pe, &s);
      t = gadd(gmul(u,g), gmul(t,be));
      t = Fp_pol_red(t, pe);
      u = gadd(u, gmul(t,pe));
      v = gadd(v, gmul(s,pe));
      pe=pe2; ae=ae2; be=be2;
    }
    ae = gerepileupto(ltop, ae2);
    aprov = Fp_poldivres(aprov,ae,pev, NULL);
    ltop = avma;
    res[i] = (long)ae;
    if (DEBUGLEVEL > 4)
      fprintferr("...lifting factor of degree %3ld. Time = %ld\n",
                 lgef(ae)-3,timer2());
  }
  return res;
}

long split_berlekamp(GEN Q, GEN *t, GEN pp, GEN pps2);

/* assume degree(a) > 0, a(0) != 0, and a squarefree */
static GEN
squff(GEN a, long klim, long hint)
{
  GEN Q,prime,primes2,factorsmodp,p1,p2,y,g,z,w,pe,*tabd,*tabdnew;
  long av=avma,tetpil,va=varn(a),da=lgef(a)-3;
  long chosenp,p,nfacp,lbit,i,j,d,e,np,nmax,lgg,nf,nft;
  ulong *tabbit, *tabkbit, *tmp;
  byteptr pt=diffptr;
  const int NOMBDEP = 5;

  if (hint < 0) return decpol(a,klim);
  if (DEBUGLEVEL > 2) timer2();
  lbit=(da>>4)+1; nmax=da+1; i=da>>1;
  if (!klim || klim>i) klim=i;
  tabdnew = (GEN*)new_chunk(nmax);
  tabbit  = (ulong*)new_chunk(lbit);
  tabkbit = (ulong*)new_chunk(lbit);
  tmp     = (ulong*)new_chunk(lbit);
  p = np = 0; prime = icopy(gun);
  while (np < NOMBDEP)
  {
    GEN minuspolx;
    p += *pt++; if (!*pt) err(primer1);
    if (!smodis((GEN)a[da+2],p)) continue;
    prime[2] = p; z = Fp_pol(a, prime);
    if (gcmp0(discsr(z))) continue;
    z = lift_intern(z);

    for (j=0; j<lbit-1; j++) tabkbit[j]=0;
    tabkbit[j]=1;
    d=0; e=da; nfacp=0;
    w = polx[va]; minuspolx = gneg(w);
    while (d < (e>>1))
    {
      d++; w = Fp_pow_mod_pol(w, prime, z, prime);
      g = Fp_pol_gcd(z, gadd(w, minuspolx), prime);
      tabdnew[d]=g; lgg=lgef(g)-3;
      if (lgg>0)
      {
	z = Fp_poldivres(z, g, prime, NULL);
        e = lgef(z)-3;
        w = Fp_poldivres(w, z, prime, ONLY_REM);
        lgg /= d; nfacp += lgg;
        if (DEBUGLEVEL>5)
          fprintferr("   %3ld %s of degree %3ld\n",
                     lgg, lgg==1?"factor": "factors",d);
	for (j=1; j<=lgg; j++) dotab(d,lbit-1,tabkbit,tmp);
      }
    }
    if (e>0)
    {
      if (DEBUGLEVEL>5)
        fprintferr("   %3ld factor of degree %3ld\n",1,e);
      tabdnew[e]=z; nfacp++;
      dotab(e,lbit-1,tabkbit,tmp);
    }
    if (np)
      for (j=0; j<lbit; j++) tabbit[j] &= tabkbit[j];
    else
      for (j=0; j<lbit; j++) tabbit[j] = tabkbit[j];
    if (DEBUGLEVEL > 4)
      fprintferr("...tried prime %3ld (%-3ld %s). Time = %ld\n",
                  p,nfacp,nfacp==1?"factor": "factors",timer2());
    if (min_deg(lbit-1,tabbit) > klim)
    {
      avma=av; y=cgetg(2,t_COL); y[1]=lcopy(a); return y;
    }
    if (nfacp<nmax)
    {
      nmax=nfacp; tabd=tabdnew; chosenp=p;
      for (j=d+1; j<e; j++) tabd[j]=polun[va];
      tabdnew = (GEN*)new_chunk(da);
    }
    np++;
  }
  prime[2]=chosenp; primes2 = shifti(prime, -1);
  nf=nmax; nft=1;
  y=cgetg(nf+1,t_COL); factorsmodp=cgetg(nf+1,t_COL);
  Q = cgetg(da+1,t_MAT);
  for (i=1; i<=da; i++) Q[i] = lgetg(da+1, t_COL);
  p1 = (GEN)Q[1]; for (i=1; i<=da; i++) p1[i] = zero;
  for (d=1; nft<=nf; d++)
  {
    g=tabd[d]; lgg=lgef(g)-3;
    if (lgg)
    {
      long n = lgg/d;
      factorsmodp[nft] = (long)normalize_mod_p(g, prime);
      if (n==1) nft++;
      else
      {
        (void)split_berlekamp(Q, (GEN*)(factorsmodp+nft),prime,primes2);
        nft += n;
      }
    }
  }
  if (DEBUGLEVEL > 4) msgtimer("splitting mod p = %ld",chosenp);

  p1=gzero; for (i=2; i<=da+2; i++) p1=addii(p1,sqri((GEN)a[i]));
  p2=absi((GEN)a[da+2]);
  p1=addii(p2,addsi(1,racine(p1)));
  p1=mulii(p1,binome(stoi(da-1),klim));
  p1=shifti(mulii(p2,p1),1);
  for (e=1,pe=gun; ; e++)
  {
    pe = mulis(pe,chosenp);
    if (cmpii(pe,p1) >= 0) break;
  }
  if (DEBUGLEVEL > 4) fprintferr("coeff bound: %Z\n",p1);

  cmbf_target = a;
  cmbf_degree = lgef(a)-3;
  cmbf_lc = (GEN)a[lgef(a)-1];
  cmbf_abslc = absi(cmbf_lc);
  cmbf_abslcxtarget = gmul(cmbf_abslc,a);
  cmbf_mod = pe;
  cmbf_modhalf = shifti(pe,-1);
  cmbf_modfax = hensel_lift_fact(a,factorsmodp,prime,pe,e);
  cmbf_modfaxn = nf;
  cmbf_fax = cgetg(nf+2,t_COL);
  cmbf_faxn = 0;
  /* sorting factors decreasing by degree helps if klim is used
   * If not, can start with first arg of 2 instead of 1, saving some time.
   * Should be sending tabbit through for more efficiency ???
   */
  cmbf(1,NULL,klim,hint); /* The Call */
  if (cmbf_degree)
  {
    if (signe(cmbf_lc)<0) cmbf_target = gneg_i(cmbf_target);
    cmbf_fax[++cmbf_faxn] = (long)cmbf_target; /* leftover factor */
  }
  if (DEBUGLEVEL>6) fprintferr("\n");
  tetpil=avma; y=cgetg(cmbf_faxn+1,t_COL);
  for (i=1; i<=cmbf_faxn; i++) y[i]=lcopy((GEN)cmbf_fax[i]);
  return gerepile(av,tetpil,y);
}

/* klim=0 habituellement, sauf si l'on ne veut chercher que les facteurs
 * de degre <= klim
 */
GEN
factpol(GEN x, long klim, long hint)
{
  GEN fa,p1,d,t,v,w, y = cgetg(3,t_MAT);
  long av=avma,av2,lx,vx,i,j,k,ex,nbfac,zval;

  if (typ(x)!=t_POL) err(notpoler,"factpol");
  if (!signe(x)) err(zeropoler,"factpol");

  ex = nbfac = 0;
  p1 = x+2; while (gcmp0((GEN)*p1)) p1++;
  zval = p1 - (x + 2);
  lx = lgef(x) - zval;
  vx = varn(x);
  if (zval)
  {
    x = cgetg(lx, t_POL); p1 -= 2;
    for (i=2; i<lx; i++) x[i] = p1[i];
    x[1] = evalsigne(1)|evalvarn(vx)|evallgef(lx);
    nbfac++;
  }
  /* now x(0) != 0 */
  if (lx==3) goto END;
  p1 = cgetg(1,t_VEC); fa=cgetg(lx,t_VEC);
  for (i=1; i<lx; i++) fa[i] = (long)p1;
  d=content(x); if (gsigne(leading_term(x)) < 0) d = gneg_i(d);
  if (!gcmp1(d)) x=gdiv(x,d);
  if (lx==4) { nbfac++; ex++; fa[1] = (long)concatsp(p1,x); goto END; }

  w=derivpol(x); t=modulargcd(x,w);
  if (!gcmp1(t)) { x=gdeuc(x,t); w=gdeuc(w,t); }
  k=1;
  while (k)
  {
    ex++; w=gadd(w, gneg_i(derivpol(x))); k=signe(w);
    if (k) { t=modulargcd(x,w); x=gdeuc(x,t); w=gdeuc(w,t); } else t=x;
    if (lgef(t) > 3)
    {
      fa[ex] = (long)squff(t,klim,hint);
      nbfac += lg(fa[ex])-1;
    }
  }
END: av2=avma;
  v=cgetg(nbfac+1,t_COL); y[1]=(long)v;
  w=cgetg(nbfac+1,t_COL); y[2]=(long)w;
  if (zval) { v[1]=lpolx[vx]; w[1]=lstoi(zval); k=1; } else k=0;
  for (i=1; i<=ex; i++)
    for (j=1; j<lg(fa[i]); j++)
    {
      k++; v[k]=lcopy(gmael(fa,i,j)); w[k]=lstoi(i);
    }
  gerepilemanyvec(av,av2,y+1,2);
  return sort_factor(y, cmpii);
}

GEN
factpol2(GEN x, long klim)
{
  return factpol(x,klim,-1);
}

/***********************************************************************/
/**                                                                   **/
/**                          FACTORISATION                            **/
/**                                                                   **/
/***********************************************************************/
#define LT 17
#define assign_or_fail(x,y) {\
  if (y==NULL) y=x; else if (!gegal(x,y)) return 0;\
}
#define tsh 6
#define typs(x,y) ((x << tsh) | y)
#define typ1(x) (x >> tsh)
#define typ2(x) (x & ((1<<tsh)-1))

static long
poltype(GEN x, GEN *ptp, GEN *ptpol, long *ptpa)
{
  long t[LT]; /* code pour 0,1,2,3,61,62,63,67,7,81,82,83,86,87,91,93,97 */
  long tx = typ(x),lx,i,j,s,pa=BIGINT;
  GEN pcx=NULL, p=NULL,pol=NULL,p1,p2;

  if (is_scalar_t(tx))
  {
    if (tx == t_POLMOD) return 0;
    x = scalarpol(x,0);
  }
  for (i=2; i<LT; i++) t[i]=0; /* t[0..1] unused */
  lx = lgef(x);
  for (i=2; i<lx; i++)
  {
    p1=(GEN)x[i];
    switch(typ(p1))
    {
      case t_INT: case t_FRAC: case t_FRACN:
        break;
      case t_REAL:
        s = precision(p1); if (s < pa) pa = s;
        t[2]=1; break;
      case t_INTMOD:
	assign_or_fail((GEN)p1[1],p);
        t[3]=1; break;
      case t_COMPLEX:
        if (!pcx)
        {
          pcx = cgetg(5,t_POL); /* x^2 + 1 */
          pcx[1] = evalsigne(1)|evalvarn(0)|m_evallgef(5),
          pcx[2]=pcx[4]=un; pcx[3]=zero;
        }
	for (j=1; j<=2; j++)
	{
	  p2 = (GEN)p1[j];
	  switch(typ(p2))
	  {
	    case t_INT: case t_FRAC: case t_FRACN:
	      assign_or_fail(pcx,pol);
	      t[4]=1; break;
	    case t_REAL:
              s = precision(p2); if (s < pa) pa = s;
	      t[5]=1; break;
	    case t_INTMOD:
	      assign_or_fail((GEN)p2[1],p);
	      assign_or_fail(pcx,pol);
	      t[6]=1; break;
	    case t_PADIC:
	      s = precp(p2) + valp(p2); if (s < pa) pa = s;
	      assign_or_fail((GEN)p2[2],p);
	      assign_or_fail(pcx,pol);
	      t[7]=1; break;
	    default: return 0;
	  }
	}
	break;
      case t_PADIC:
        s = precp(p1) + valp(p1); if (s < pa) pa = s;
	assign_or_fail((GEN)p1[2],p);
        t[8]=1; break;
      case t_QUAD:
	for (j=2; j<=3; j++)
	{
	  p2 = (GEN)p1[j];
	  switch(typ(p2))
	  {
	    case t_INT: case t_FRAC: case t_FRACN:
	      assign_or_fail((GEN)p1[1],pol);
	      t[9]=1; break;
	    case t_REAL:
	      s = precision(p2); if (s < pa) pa = s;
	      if (gsigne(discsr((GEN)p1[1]))>0) t[10]=1; else t[12]=1;
	      break;
	    case t_INTMOD:
	      assign_or_fail((GEN)p2[1],p);
	      assign_or_fail((GEN)p1[1],pol);
	      t[11]=1; break;
	    case t_PADIC:
	      s = precp(p2) + valp(p2); if (s < pa) pa = s;
	      assign_or_fail((GEN)p2[2],p);
	      assign_or_fail((GEN)p1[1],pol);
	      t[13]=1; break;
	    default: return 0;
	  }
	}
	break;
      case t_POLMOD:
	assign_or_fail((GEN)p1[1],pol);
	for (j=1; j<=2; j++)
	{
	  GEN pbis = NULL, polbis = NULL;
	  long pabis;
	  switch(poltype((GEN)p1[j],&pbis,&polbis,&pabis))
	  {
	    case t_INT: t[14]=1; break;
	    case t_INTMOD: t[15]=1; break;
	    case t_PADIC: t[16]=1; if (pabis<pa) pa=pabis; break;
	    default: return 0;
	  }
	  if (pbis) assign_or_fail(pbis,p);
	  if (polbis) assign_or_fail(polbis,pol);
	}
	break;
      default: return 0;
    }
  }
  if (t[5]||t[12])
  {
    if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
    *ptpa=pa; return t_COMPLEX;
  }
  if (t[2]||t[10])
  {
    if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
    *ptpa=pa; return t[4]?t_COMPLEX:t_REAL;
  }
  if (t[6]||t[11]||t[15])
  {
    *ptpol=pol; *ptp=p;
    i = t[15]? t_POLMOD: (t[11]? t_QUAD: t_COMPLEX);
    return typs(i, t_INTMOD);
  }
  if (t[7]||t[13]||t[16])
  {
    *ptpol=pol; *ptp=p; *ptpa=pa;
    i = t[16]? t_POLMOD: (t[13]? t_QUAD: t_COMPLEX);
    return typs(i, t_PADIC);
  }
  if (t[4]||t[9]||t[14])
  {
    *ptpol=pol;
    i = t[14]? t_POLMOD: (t[9]? t_QUAD: t_COMPLEX);
    return typs(i, t_INT);
  }
  if (t[3]) { *ptp=p; return t_INTMOD; }
  if (t[8]) { *ptp=p; *ptpa=pa; return t_PADIC; }
  return t_INT;
}
#undef LT

GEN
factor0(GEN x,long flag)
{
  long tx=typ(x);

  if (flag<0) return factor(x);
  if (is_matvec_t(tx)) return gboundfact(x,flag);
  if (tx==t_INT || tx==t_FRAC || tx==t_FRACN) return boundfact(x,flag);
  err(talker,"partial factorization is not meaningful here");
  return NULL; /* not reached */
}

static GEN
poldeg1(long v, GEN x0, GEN x1)
{
  GEN x = cgetg(4,t_POL);
  x[1] = evalsigne(1) | evalvarn(v) | evallgef(4);
  x[2] = (long)x0; x[3] = (long)x1; return normalizepol(x);
}

GEN
factor(GEN x)
{
  long tx=typ(x),lx,av,tetpil,i,j,pa,v,r1;
  GEN  y,p,p1,p2,p3,p4,p5,pol;

  if (is_matvec_t(tx))
  {
    lx=lg(x); y=cgetg(lx,tx);
    for (i=1; i<lx; i++) y[i]=(long)factor((GEN)x[i]);
    return y;
  }
  if (gcmp0(x))
  {
    y=cgetg(3,t_MAT);
    p1=cgetg(2,t_COL); y[1]=(long)p1; p1[1]=zero;
    p2=cgetg(2,t_COL); y[2]=(long)p2; p2[1]=un;
    return y;
  }
  switch(tx)
  {
    case t_INT: return decomp(x);

    case t_FRACN:
      av=avma; x=gred(x); /* fall through */
    case t_FRAC:
      if (tx==t_FRAC) av=avma;
      p1 = decomp((GEN)x[1]);
      p2 = decomp((GEN)x[2]);
      p4 = concatsp((GEN)p1[1], (GEN)p2[1]);
      p5 = concatsp((GEN)p1[2], gneg_i((GEN)p2[2]));
      p3 = sindexsort(p4); lx = lg(p3);
      tetpil=avma; y=cgetg(3,t_MAT);
      p1 = cgetg(lx,t_COL); y[1]=(long)p1;
      p2 = cgetg(lx,t_COL); y[2]=(long)p2;
      for (i=1; i<lx; i++)
      {
      	p1[i] = licopy((GEN)p4[p3[i]]);
      	p2[i] = licopy((GEN)p5[p3[i]]);
      }
      return gerepile(av,tetpil,y);

    case t_POL:
      tx=poltype(x,&p,&pol,&pa);
      switch(tx)
      {
        case 0: err(impl,"factor for general polynomials");
	case t_INT: return factpol(x,0,1);
	case t_INTMOD: return factmod(x,p);

	case t_COMPLEX: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);
	  av=avma; p1=roots(x,pa); tetpil=avma;
          p2=cgetg(lx,t_COL);
	  for (i=1; i<lx; i++)
            p2[i] = (long)poldeg1(v, gneg((GEN)p1[i]), gun);
	  y[1]=lpile(av,tetpil,p2);
	  p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;
          y[2]=(long)p3; return y;

	case t_REAL: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);
	  av=avma; p1=roots(x,pa); tetpil=avma;
	  for(r1=1; r1<lx; r1++)
            if (signe(gmael(p1,r1,2))) break;
	  lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
	  for(i=1; i<r1; i++)
            p2[i] = (long)poldeg1(v, negr(gmael(p1,i,1)), gun);
	  for(   ; i<lx; i++)
	  {
	    GEN a = (GEN) p1[2*i-r1];
	    p = cgetg(5, t_POL); p2[i] = (long)p;
	    p[1] = evalsigne(1) | evalvarn(v) | evallgef(5);
	    p[2] = lnorm(a);
	    p[3] = lmul2n((GEN)a[1],1); setsigne(p[3],-signe(p[3]));
	    p[4] = un;
	  }
	  y[1]=lpile(av,tetpil,p2);
	  p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;
          y[2]=(long)p3; return y;

	case t_PADIC: return factorpadic4(x,p,pa);

        default:
        {
          long killv;
	  av=avma; x = dummycopy(x); lx=lgef(x);
          pol = dummycopy(pol);
          v = manage_var(4,NULL);
          for(i=2; i<lx; i++)
          {
            p1=(GEN)x[i];
            switch(typ(p1))
            {
              case t_QUAD: p1++;
              case t_COMPLEX:
                p2 = cgetg(3, t_POLMOD); x[i] = (long) p2;
                p2[1] = (long)pol;
                p2[2] = (long)poldeg1(v, (GEN)p1[1],(GEN)p1[2]);
            }
          }
          killv = (avma != (long)pol);
          if (killv) setvarn(pol, fetch_var());
          switch (typ2(tx))
          {
            case t_INT: p1 = polfnf(x,pol); break;
            case t_INTMOD: p1 = factmod9(x,p,pol); break;
	    default: err(impl,"factor of general polynomial");
          }
          switch (typ1(tx))
          {
            case t_POLMOD:
              if (killv) delete_var();
              return gerepileupto(av,p1);
            case t_COMPLEX: p5 = gi; break;
            case t_QUAD: p5=cgetg(4,t_QUAD);
              p5[1]=(long)pol; p5[2]=zero; p5[3]=un;
              break;
	    default: err(impl,"factor of general polynomial");
          }
          p2=(GEN)p1[1];
          for(i=1; i<lg(p2); i++)
          {
            p3=(GEN)p2[i];
            for(j=2; j<lgef(p3); j++)
            {
              p4=(GEN)p3[j];
              if(typ(p4)==t_POLMOD) p3[j]=lsubst((GEN)p4[2],v,p5);
            }
          }
          if (killv) delete_var();
          tetpil=avma; y=cgetg(3,t_MAT);
          y[1]=lcopy(p2);y[2]=lcopy((GEN)p1[2]);
          return gerepile(av,tetpil,y);
        }
      }

    case t_RFRACN:
      av=avma; x=gred_rfrac(x); /* fall through */
    case t_RFRAC:
      if (tx==t_RFRAC) av=avma;
      p1=factor((GEN)x[1]);
      p2=factor((GEN)x[2]); p3=gneg_i((GEN)p2[2]);
      tetpil=avma; y=cgetg(3,t_MAT);
      y[1]=lconcat((GEN)p1[1],(GEN)p2[1]);
      y[2]=lconcat((GEN)p1[2],p3);
      return gerepile(av,tetpil,y);
  }
  err(impl,"general factorization");
  return NULL; /* not reached */
}
#undef typ1
#undef typ2
#undef typs
#undef tsh

GEN
divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
{
  long i,k,lx = lg(x);

  if (lx == 1) return gun;
  if (lx == 2) return gcopy((GEN)x[1]);
  x = dummycopy(x); k = lx;
  while (k > 2)
  {
    if (DEBUGLEVEL>7)
      fprintferr("prod: remaining objects %ld\n",k-1);
    lx = k; k = 1;
    for (i=1; i<lx-1; i+=2)
      x[k++] = (long)mul((GEN)x[i],(GEN)x[i+1]);
    if (i < lx) x[k++] = x[i];
  }
  return (GEN)x[1];
}

static GEN static_nf;

static GEN
myidealmul(GEN x, GEN y) { return idealmul(static_nf, x, y); }

static GEN
myidealpow(GEN x, GEN n) { return idealpow(static_nf, x, n); }

GEN
factorback(GEN fa, GEN nf)
{
  long av=avma,k,l,lx;
  GEN ex, x;
  GEN (*_mul)(GEN,GEN);
  GEN (*_pow)(GEN,GEN);

  if (typ(fa)!=t_MAT || lg(fa)!=3)
    err(talker,"incorrect factorisation in factorback");
  ex=(GEN)fa[2]; fa=(GEN)fa[1];
  lx = lg(fa); if (lx == 1) return gun;
  x = cgetg(lx,t_VEC);
  if (nf)
  {
    static_nf = nf;
    _mul = &myidealmul;
    _pow = &myidealpow;
  }
  else
  {
    _mul = &gmul;
    _pow = &powgi;
  }
  for (l=1,k=1; k<lx; k++)
    if (signe(ex[k]))
      x[l++] = (long)_pow((GEN)fa[k],(GEN)ex[k]);
  setlg(x,l);
  return gerepileupto(av, divide_conquer_prod(x, _mul));
}

GEN
gisirreducible(GEN x)
{
  long av=avma, tx = typ(x),l,i;
  GEN y;

  if (is_matvec_t(tx))
  {
    l=lg(x); y=cgetg(l,tx);
    for (i=1; i<l; i++) y[i]=(long)gisirreducible((GEN)x[i]);
    return y;
  }
  if (is_intreal_t(tx) || is_frac_t(tx)) return gzero;
  if (tx!=t_POL) err(notpoler,"gisirreducible");
  l=lgef(x); if (l<=3) return gzero;
  y=factor(x); avma=av;
  return (lgef(gcoeff(y,1,1))==l)?gun:gzero;
}

/*******************************************************************/
/*                                                                 */
/*                         PGCD GENERAL                            */
/*                                                                 */
/*******************************************************************/

GEN
gcd0(GEN x, GEN y, long flag)
{
  switch(flag)
  {
    case 0: return ggcd(x,y);
    case 1: return modulargcd(x,y);
    case 2: return srgcd(x,y);
    default: err(flagerr,"gcd");
  }
  return NULL; /* not reached */
}

/* x is a COMPLEX or a QUAD */
static GEN
triv_cont_gcd(GEN x, GEN y)
{
  long av = avma, tetpil;
  GEN p1 = (typ(x)==t_COMPLEX)? ggcd((GEN)x[1],(GEN)x[2])
                              : ggcd((GEN)x[2],(GEN)x[3]);
  tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
}

static GEN
cont_gcd(GEN x, GEN y)
{
  long av = avma, tetpil;
  GEN p1 = content(x);

  tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
}

/* y is a PADIC, x a rational number or an INTMOD */
static GEN
padic_gcd(GEN x, GEN y)
{
  long v = ggval(x,(GEN)y[2]), w = valp(y);
  if (w < v) v = w;
  return gpuigs((GEN)y[2], v);
}

#define fix_frac(z) if (signe(z[2])<0)\
{\
  setsigne(z[1],-signe(z[1]));\
  setsigne(z[2],1);\
}

GEN
ggcd(GEN x, GEN y)
{
  long l,av,tetpil,i,vx,vy, tx = typ(x), ty = typ(y);
  GEN p1,z;

  if (tx>ty) { p1=x; x=y; y=p1; l=tx; tx=ty; ty=l; }
  if (is_matvec_t(ty))
  {
    l=lg(y); z=cgetg(l,ty);
    for (i=1; i<l; i++) z[i]=lgcd(x,(GEN)y[i]);
    return z;
  }
  if (gcmp0(x)) return gcopy(y);
  if (gcmp0(y)) return gcopy(x);
  if (is_const_t(tx))
  {
    if (ty == t_FRACN)
    {
      if (tx==t_INTMOD)
      {
        av=avma; y = gred(y); tetpil=avma;
        return gerepile(av,tetpil,ggcd(x,y));
      }
      ty = t_FRAC;
    }
    if (tx == t_FRACN) tx = t_FRAC;
    if (ty == tx) switch(tx)
    {
      case t_INT:
        return mppgcd(x,y);

      case t_INTMOD: z=cgetg(3,t_INTMOD);
        if (egalii((GEN)x[1],(GEN)y[1]))
          { copyifstack(x[1],z[1]); }
        else
          z[1]=lmppgcd((GEN)x[1],(GEN)y[1]);
        if (gcmp1((GEN)z[1])) z[2]=zero;
        else
        {
          av=avma; p1=mppgcd((GEN)z[1],(GEN)x[2]);
          if (!gcmp1(p1))
          {
            tetpil=avma;
            p1=gerepile(av,tetpil,mppgcd(p1,(GEN)y[2]));
          }
          z[2]=(long)p1;
        }
        return z;

      case t_FRAC: z=cgetg(3,t_FRAC);
        (void)new_chunk(lgefint(x[2])+lgefint(y[2]));
        p1 = divii((GEN)y[2], mppgcd((GEN)x[2], (GEN)y[2]));
        avma = (long)z;
        z[2] = lmulii(p1, (GEN)x[2]);
        z[1] = (long)mppgcd((GEN)x[1], (GEN)y[1]);
        fix_frac(z); return z;

      case t_COMPLEX:
        av=avma; p1=gdiv(x,y);
        if (gcmp0((GEN)p1[2]))
        {
          p1=(GEN)p1[1];
          switch(typ(p1))
          {
            case t_INT:
              avma=av; return gcopy(y);
            case t_FRAC: case t_FRACN:
              tetpil=avma; return gerepile(av,tetpil,gdiv(y,(GEN)p1[2]));
            default: avma=av; return gun;
          }
        }
        avma=av;
        if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) return gcopy(y);
        p1=gdiv(y,x); avma=av;
        if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) return gcopy(x);
        return triv_cont_gcd(y,x);

      case t_PADIC:
        if (!egalii((GEN)x[2],(GEN)y[2])) return gun;
        return gpuigs((GEN)y[2],min(valp(y),valp(x)));

      case t_QUAD:
        av=avma; p1=gdiv(x,y);
        if (gcmp0((GEN)p1[3]))
        {
          p1=(GEN)p1[2];
          if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
          tetpil=avma; return gerepile(av,tetpil, gdiv(y,(GEN)p1[2]));
        }
        avma=av;
        if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) return gcopy(y);
        p1=gdiv(y,x); avma=av;
        if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) return gcopy(x);
        return triv_cont_gcd(y,x);

      default: return gun; /* t_REAL */
    }
    if (is_const_t(ty)) switch(tx)
    {
      case t_INT:
        switch(ty)
        {
          case t_INTMOD: z=cgetg(3,t_INTMOD);
            copyifstack(y[1],z[1]); av=avma;
            p1=mppgcd((GEN)y[1],(GEN)y[2]);
            if (!gcmp1(p1))
              { tetpil=avma; p1=gerepile(av,tetpil,mppgcd(x,p1)); }
            z[2]=(long)p1; return z;

          case t_FRAC: z=cgetg(3,t_FRAC);
            z[1]=lmppgcd(x,(GEN)y[1]);
            z[2]=licopy((GEN)y[2]); return z;

          case t_COMPLEX: case t_QUAD:
            return triv_cont_gcd(y,x);

          case t_PADIC:
            return padic_gcd(x,y);

          default: /* t_REAL */
            return gun;
        }

      case t_INTMOD:
        switch(ty)
        {
          case t_FRAC:
            av=avma; p1=mppgcd((GEN)x[1],(GEN)y[2]); tetpil=avma;
            if (!gcmp1(p1))
              err(talker,"non invertible fraction in a gcd with INTMOD");
            return gerepile(av,tetpil, ggcd((GEN)y[1],x));

          case t_COMPLEX: case t_QUAD:
            return triv_cont_gcd(y,x);

          case t_PADIC:
            return padic_gcd(x,y);
        }

      case t_FRAC:
        switch(ty)
        {
          case t_COMPLEX: case t_QUAD:
            return triv_cont_gcd(y,x);

          case t_PADIC:
            return padic_gcd(x,y);
        }

      case t_COMPLEX: /* ty = PADIC or QUAD */
        return triv_cont_gcd(x,y);

      case t_PADIC: /* ty = QUAD */
        return triv_cont_gcd(y,x);

      default: /* tx = t_REAL */
        return gun;
    }
    return cont_gcd(y,x);
  }

  vx=gvar9(x); vy=gvar9(y);
  if (vy<vx) return cont_gcd(y,x);
  if (vx<vy) return cont_gcd(x,y);
  switch(tx)
  {
    case t_POLMOD:
      switch(ty)
      {
	case t_POLMOD: z=cgetg(3,t_POLMOD);
          if (gegal((GEN)x[1],(GEN)y[1]))
	    { copyifstack(x[1],z[1]); }
          else
            z[1] = lgcd((GEN)x[1],(GEN)y[1]);
	  if (lgef(z[1])<=3) z[2]=zero;
	  else
	  {
	    av=avma; p1=ggcd((GEN)z[1],(GEN)x[2]);
	    if (lgef(p1)>3)
	    {
	      tetpil=avma;
              p1=gerepile(av,tetpil,ggcd(p1,(GEN)y[2]));
	    }
	    z[2]=(long)p1;
	  }
	  return z;

	case t_POL: z=cgetg(3,t_POLMOD);
          copyifstack(x[1],z[1]); av=avma;
          p1=ggcd((GEN)x[1],(GEN)x[2]);
	  if (lgef(p1)>3)
            { tetpil=avma; p1=gerepile(av,tetpil,ggcd(y,p1)); }
	  z[2]=(long)p1; return z;

	case t_RFRAC:
	  av=avma; p1=ggcd((GEN)x[1],(GEN)y[2]); tetpil=avma;
          if (!gcmp1(p1))
            err(talker,"non invertible fraction in a gcd with POLMOD");
	  return gerepile(av,tetpil,ggcd((GEN)y[1],x));

	case t_RFRACN:
	  av=avma; p1=gred_rfrac(y); tetpil=avma;
	  return gerepile(av,tetpil,ggcd(p1,x));
      }
      break;

    case t_POL:
      switch(ty)
      {
	case t_POL:
	  return srgcd(x,y);

	case t_SER:
	  return gpuigs(polx[vx],min(valp(y),gval(x,vx)));

	case t_RFRAC: case t_RFRACN: av=avma; z=cgetg(3,ty);
          z[1]=lgcd(x,(GEN)y[1]);
          z[2]=lcopy((GEN)y[2]); return z;
      }
      break;

    case t_SER:
      switch(ty)
      {
	case t_SER:
	  return gpuigs(polx[vx],min(valp(x),valp(y)));

	case t_RFRAC: case t_RFRACN:
	  return gpuigs(polx[vx],min(valp(x),gval(y,vx)));
      }
      break;

    case t_RFRAC: case t_RFRACN: z=cgetg(3,ty);
      if (!is_rfrac_t(ty))
        err(talker,"forbidden gcd rational function with vector/matrix");
      p1 = gdiv((GEN)y[2], ggcd((GEN)x[2], (GEN)y[2]));
      tetpil = avma;
      z[2] = lpile((long)z,tetpil,gmul(p1, (GEN)x[2]));
      z[1] = lgcd((GEN)x[1], (GEN)y[1]); return z;
  }
  err(talker,"gcd vector/matrix with a forbidden type");
  return NULL; /* not reached */
}

GEN
glcm(GEN x, GEN y)
{
  long av=avma,tx,ty=typ(y),i,l;
  GEN p1,p2,z;

  if (is_matvec_t(ty))
  {
    l=lg(y); z=cgetg(l,ty);
    for (i=1; i<l; i++) z[i]=(long)glcm(x,(GEN)y[i]);
    return z;
  }
  tx=typ(x);
  if (is_matvec_t(tx))
  {
    l=lg(x); z=cgetg(l,tx);
    for (i=1; i<l; i++) z[i]=(long)glcm((GEN)x[i],y);
    return z;
  }
  if (gcmp0(x)) return gzero;
  if (tx==t_INT && ty==t_INT)
  {
    p1 = mppgcd(x,y); if (is_pm1(p1)) { avma = av; return mulii(x,y); }
    p2 = mulii(divii(y,p1), x);
    if (signe(p2)<0) setsigne(p2,1);
    return gerepileupto(av, p2);
  }
  p1 = ggcd(x,y); if (gcmp1(p1)) {avma = av; return gmul(x,y); }
  p2 = gmul(gdiv(y,p1), x);
  if (typ(p2)==t_INT && signe(p2)<0) setsigne(p2,1);
  if (typ(p2)==t_POL)
  {
    p1=leading_term(p2);
    if (typ(p1)==t_INT && signe(p1)<0) p2=gneg(p2);
  }
  return gerepileupto(av,p2);
}

static GEN
polgcdnun(GEN x, GEN y)
{
  long av1, av = avma, lim = stack_lim(av,1);
  GEN p1, yorig = y;

  for(;;)
  {
    av1 = avma; p1 = gres(x,y);
    if (gcmp0(p1))
    {
      avma = av1;
      return (y==yorig)? gcopy(y): gerepileupto(av,y);
    }
    x = y; y = p1;
    if (low_stack(lim,stack_lim(av,1)))
    {
      GEN *gptr[2]; x = gcopy(x); gptr[0]=&x; gptr[1]=&y;
      if(DEBUGMEM>1) err(warnmem,"polgcdnun");
      gerepilemanysp(av,av1,gptr,2);
    }
  }
}

/* return 1 if coeff explosion is not possible */
static int
issimplefield(GEN x)
{
  long lx,i;
  switch(typ(x))
  {
    case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
      return 1;
    case t_POL:
      lx=lgef(x);
      for (i=2; i<lx; i++)
	if (!issimplefield((GEN)x[i])) return 0;
      return 1;
    case t_COMPLEX: case t_POLMOD:
      return issimplefield((GEN)x[1]) || issimplefield((GEN)x[2]);
  }
  return 0;
}

static int
isrational(GEN x)
{
  long i;
  for (i=lgef(x)-1; i>1; i--)
    switch(typ(x[i]))
    {
      case t_INT:
      case t_FRAC: break;
      default: return 0;
    }
  return 1;
}

static int
isinexactfield(GEN x)
{
  long lx,i;
  switch(typ(x))
  {
    case t_REAL: case t_PADIC: case t_SER:
      return 1;
    case t_POL:
      lx=lgef(x);
      for (i=2; i<lx; i++)
	if (!isinexactfield((GEN)x[i])) return 0;
      return 1;
    case t_COMPLEX: case t_POLMOD:
      return isinexactfield((GEN)x[1]) || isinexactfield((GEN)x[2]);
  }
  return 0;
}

static GEN
gcdmonome(GEN x, GEN y)
{
  long tetpil,av=avma, lx=lgef(x), v=varn(x), e=gval(y,v);
  GEN p1,p2;

  if (lx-3<e) e=lx-3;
  p1=ggcd(leading_term(x),content(y)); p2=gpuigs(polx[v],e);
  tetpil=avma; return gerepile(av,tetpil,gmul(p1,p2));
}

/***********************************************************************/
/**                                                                   **/
/**                         BEZOUT GENERAL                            **/
/**                                                                   **/
/***********************************************************************/

static GEN
polinvinexact(GEN x, GEN y)
{
  long i,dx=lgef(x)-3,dy=lgef(y)-3,lz=dx+dy, av=avma, tetpil;
  GEN v,z;

  z=cgetg(dy+2,t_POL); z[1]=y[1];
  v=cgetg(lz+1,t_COL);
  for (i=1; i<lz; i++) v[i]=zero;
  v[lz]=un; v=gauss(sylvestermatrix(y,x),v);
  for (i=2; i<dy+2; i++) z[i]=v[lz-i+2];
  z = normalizepol_i(z, dy+2); tetpil = avma;
  return gerepile(av,tetpil,gcopy(z));
}

static GEN
polinvmod(GEN x, GEN y)
{
  long av,av1,tx,vx=varn(x),vy=varn(y);
  GEN u,v,d,p1;

  while (vx!=vy)
  {
    if (vx > vy)
    {
      d=cgetg(3,t_RFRAC);
      d[1]=(long)polun[vx];
      d[2]=lcopy(x); return d;
    }
    if (lgef(x)!=3) err(talker,"non-invertible polynomial in polinvmod");
    x=(GEN)x[2]; vx=gvar(x);
  }
  tx=typ(x);
  if (tx==t_POL)
  {
    if (isinexactfield(x) || isinexactfield(y))
      return polinvinexact(x,y);

    av=avma; d=subresext(x,y,&u,&v);
    if (gcmp0(d)) err(talker,"non-invertible polynomial in polinvmod");
    if (typ(d)==t_POL && varn(d)==vx)
    {
      if (lgef(d)>3) err(talker,"non-invertible polynomial in polinvmod");
      d=(GEN)d[2];
    }
    av1=avma; return gerepile(av,av1,gdiv(u,d));
  }
  if (!is_rfrac_t(tx)) err(typeer,"polinvmod");
  av=avma; p1=gmul((GEN)x[1],polinvmod((GEN)x[2],y));
  av1=avma; return gerepile(av,av1,gmod(p1,y));
}

GEN
gbezout(GEN x, GEN y, GEN *u, GEN *v)
{
  long tx=typ(x),ty=typ(y);

  if (tx==t_INT && ty==t_INT) return bezout(x,y,u,v);
  if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) err(typeer,"gbezout");
  return bezoutpol(x,y,u,v);
}

GEN
vecbezout(GEN x, GEN y)
{
  GEN z=cgetg(4,t_VEC);
  z[3]=(long)gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
  return z;
}

GEN
vecbezoutres(GEN x, GEN y)
{
  GEN z=cgetg(4,t_VEC);
  z[3]=(long)subresext(x,y,(GEN*)(z+1),(GEN*)(z+2));
  return z;
}

/*******************************************************************/
/*                                                                 */
/*                    CONTENU ET PARTIE PRIMITIVE                  */
/*                                                                 */
/*******************************************************************/

GEN
content(GEN x)
{
  long av,tetpil,lx,i,tx=typ(x);
  GEN p1,p2;

  if (is_scalar_t(tx))
    return tx==t_POLMOD? content((GEN)x[2]): gcopy(x);
  av = avma;
  switch(tx)
  {
    case t_RFRAC: case t_RFRACN:
      p1=content((GEN)x[1]);
      p2=content((GEN)x[2]); tetpil=avma;
      return gerepile(av,tetpil,gdiv(p1,p2));

    case t_VEC: case t_COL: case t_MAT:
      lx = lg(x); if (lx==1) return gun;
      p1=content((GEN)x[1]);
      for (i=2; i<lx; i++) p1 = ggcd(p1,content((GEN)x[i]));
      return gerepileupto(av,p1);

    case t_POL:
      if (!signe(x)) return gzero;
      lx = lgef(x); break;
    case t_SER:
      if (!signe(x)) return gzero;
      lx = lg(x); break;
    case t_QFR: case t_QFI:
      lx = 4; break;

    default: err(typeer,"content");
  }
  for (i=lontyp[tx]; i<lx; i++)
    if (typ(x[i]) != t_INT) break;
  lx--; p1=(GEN)x[lx];
  if (i > lx)
  { /* integer coeffs */
    while (lx>lontyp[tx])
    {
      lx--; p1=mppgcd(p1,(GEN)x[lx]);
      if (is_pm1(p1)) { avma=av; return gun; }
    }
  }
  else
  {
    while (lx>lontyp[tx])
    {
      lx--; p1=ggcd(p1,(GEN)x[lx]);
    }
    if (isinexactreal(p1)) { avma=av; return gun; }
  }
  return av==avma? gcopy(p1): gerepileupto(av,p1);
}

/*******************************************************************/
/*                                                                 */
/*                         SOUS RESULTANT                          */
/*                                                                 */
/*******************************************************************/
GEN
gdivexact(GEN x, GEN y)
{
  long i,lx;
  GEN z;
  if (gcmp1(y)) return x;
  switch(typ(x))
  {
    case t_INT:
      if (typ(y)==t_INT) return divii(x,y);
      if (!signe(x)) return gzero;
      break;
    case t_INTMOD:
    case t_POLMOD: return gmul(x,ginv(y));
    case t_POL:
      switch(typ(y))
      {
        case t_INTMOD:
        case t_POLMOD: return gmul(x,ginv(y));
        case t_POL:
          if (varn(x)==varn(y)) return poldivres(x,y, ONLY_DIVIDES_EXACT);
      }
      lx = lgef(x); z = cgetg(lx,t_POL);
      for (i=2; i<lx; i++) z[i]=(long)gdivexact((GEN)x[i],y);
      z[1]=x[1]; return z;
  }
  if (DEBUGLEVEL) err(warner,"missing case in gdivexact");
  return gdiv(x,y);
}

static GEN
init_resultant(GEN x, GEN y)
{
  long tx,ty;
  if (gcmp0(x) || gcmp0(y)) return gzero;
  tx = typ(x); ty = typ(y);
  if (is_scalar_t(tx) || is_scalar_t(ty))
  {
    if (tx==t_POL) return gpuigs(y,lgef(x)-3);
    if (ty==t_POL) return gpuigs(x,lgef(y)-3);
    return gun;
  }
  if (tx!=t_POL || ty!=t_POL) err(typeer,"subresall");
  if (varn(x)==varn(y)) return NULL;
  return (varn(x)<varn(y))? gpuigs(y,lgef(x)-3): gpuigs(x,lgef(y)-3);
}

/* return coefficients s.t x = x_0 X^n + ... + x_n */
static GEN 
revpol(GEN x)
{
  long i,n = lgef(x)-3;
  GEN y = cgetg(n+3,t_POL);
  y[1] = x[1]; x += 2; y += 2;
  for (i=0; i<=n; i++) y[i] = x[n-i];
  return y;
}

/* assume dx >= dy, y non constant */
GEN
pseudorem(GEN x, GEN y)
{
  long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,lx, p;

  if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");
  (void)new_chunk(2);
  dx=lgef(x)-3; x = revpol(x);
  dy=lgef(y)-3; y = revpol(y); dz=dx-dy; p = dz+1;
  av2 = avma; lim = stack_lim(av2,1);
  for (;;)
  {
    x[0] = lneg((GEN)x[0]); p--;
    for (i=1; i<=dy; i++)
      x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));
    for (   ; i<=dx; i++)
      x[i] = lmul((GEN)y[0], (GEN)x[i]);
    do { x++; dx--; } while (dx >= 0 && gcmp0((GEN)x[0]));
    if (dx < dy) break;
    if (low_stack(lim,stack_lim(av2,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"pseudorem dx = %ld >= %ld",dx,dy);
      gerepilemanycoeffs(av2,x,dx+1);
    }
  }
  if (dx < 0) return zeropol(vx);
  lx = dx+3; x -= 2;
  x[0]=evaltyp(t_POL) | evallg(lx);
  x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx);
  x = revpol(x) - 2;
  return gerepileupto(av, gmul(x, gpowgs((GEN)y[0], p)));
}

/* Si sol != NULL:
 *   met dans *sol le dernier polynome non nul de la polynomial remainder
 *   sequence si elle a ete effectuee, 0 sinon
 */
GEN
subresall(GEN u, GEN v, GEN *sol)
{
  long degq,av,av2,lim,tetpil,dx,dy,du,dv,dr,signh;
  GEN g,h,r,p1,p2,p3,p4;

  if (sol) *sol=gzero;
  if ((r = init_resultant(u,v))) return r;

  if (isinexactfield(u) || isinexactfield(v)) return resultant2(u,v);
  dx=lgef(u); dy=lgef(v); signh=1;
  if (dx<dy)
  {
    p1=u; u=v; v=p1; du=dx; dx=dy; dy=du;
    if ((dx&1) == 0 && (dy&1) == 0) signh = -signh;
  }
  if (dy==3) return gpowgs((GEN)v[2],dx-3);
  av=avma;
  p3=content(u); if (gcmp1(p3)) p3=NULL; else u=gdiv(u,p3);
  p4=content(v); if (gcmp1(p4)) p4=NULL; else v=gdiv(v,p4);
  g=gun; h=gun; av2 = avma; lim = stack_lim(av2,1);
  for(;;)
  {
    r = pseudorem(u,v); dr = lgef(r);
    if (dr==2)
    {
      if (sol) {avma = (long)(r+2); *sol=gerepileupto(av,v);} else avma = av;
      return gzero;
    }
    du=lgef(u); dv=lgef(v);
    degq=du-dv; u=v;
    p1 = g; g = leading_term(u);
    switch(degq)
    {
      case 0: break;
      case 1:
        p1 = gmul(h,p1); h = g; break;
      default:
        p1 = gmul(gpuigs(h,degq),p1);
        h = gdivexact(gpuigs(g,degq), gpuigs(h,degq-1));
    }
    if ((du&1) == 0 && (dv&1) == 0) signh = -signh;
    v = gdivexact(r,p1);
    if (dr==3) break;
    if (low_stack(lim,stack_lim(av2,1)))
    {
      GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
      if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr);
      gerepilemany(av2,gptr,4);
    }
  }
  if (dv==4) { tetpil=avma; p2=gcopy((GEN)v[2]); }
  else
  {
    if (dv == 3) err(bugparier,"subres");
    p1=gpuigs((GEN)v[2],dv-3); p2=gpuigs(h,dv-4);
    tetpil=avma; p2=gdiv(p1,p2);
  }
  if (p3) { p1=gpuigs(p3,dy-3); tetpil=avma; p2=gmul(p2,p1); }
  if (p4) { p1=gpuigs(p4,dx-3); tetpil=avma; p2=gmul(p2,p1); }
  if (signh<0) { tetpil=avma; p2=gneg(p2); }
  {
    GEN *gptr[2]; gptr[0]=&p2; if (sol) { *sol=gcopy(u); gptr[1]=sol; }
    gerepilemanysp(av,tetpil,gptr,sol?2:1);
    return p2;
  }
}

static GEN
scalar_res(GEN x, GEN y, GEN *U, GEN *V)
{
  long dx=lgef(x)-4;
  *V=gpuigs(y,dx); *U=gzero; return gmul(y,*V);
}

/* calcule U et V tel que Ux+By=resultant(x,y) */
GEN
subresext(GEN x, GEN y, GEN *U, GEN *V)
{
  long degq,av,tetpil,tx,ty,dx,dy,du,dv,dr,signh;
  GEN z,g,h,r,p1,p2,p3,p4,u,v,lpu,um1,uze, *gptr[2];

  if (gcmp0(x) || gcmp0(y)) { *U = *V = gzero; return gzero; }
  tx=typ(x); ty=typ(y);
  if (is_scalar_t(tx) || is_scalar_t(ty))
  {
    if (tx==t_POL) return scalar_res(x,y,U,V);
    if (ty==t_POL) return scalar_res(y,x,V,U);
    *U=ginv(x); *V=gzero; return gun;
  }
  if (tx!=t_POL || ty!=t_POL) err(typeer,"subresext");
  if (varn(x)!=varn(y))
    return (varn(x)<varn(y))? scalar_res(x,y,U,V)
                            : scalar_res(y,x,V,U);
  dx=lgef(x); dy=lgef(y); signh=1;
  if (dx<dy)
  {
    GEN *W = U; U=V; V=W;
    du=dx; dx=dy; dy=du; p1=x; x=y; y=p1;
    if ((dx&1) == 0 && (dy&1) == 0) signh = -signh;
  }
  if (dy==3)
  {
    *V=gpuigs((GEN)y[2],dx-4);
    *U=gzero; return gmul(*V,(GEN)y[2]);
  }
  av=avma;
  p3=content(x); if (gcmp1(p3)) {p3 = NULL; u=x; } else u=gdiv(x,p3);
  p4=content(y); if (gcmp1(p4)) {p4 = NULL; v=y; } else v=gdiv(y,p4);
  g=gun; h=gun; um1=gun; uze=gzero;
  for(;;)
  {
    du=lgef(u); dv=lgef(v); degq=du-dv;
    lpu=gpuigs((GEN)v[dv-1],degq+1);
    p1=gmul(lpu,u); p2=poldivres(p1,v,&r);
    dr=lgef(r); if (dr==2) { *U=gzero; *V=gzero; avma=av; return gzero; }

    p1=gsub(gmul(lpu,um1),gmul(p2,uze));
    um1=uze; uze=p1; u=v;
    p1 = g; g = leading_term(u);
    switch(degq)
    {
      case 0: break;
      case 1: p1 = gmul(h,p1); h = g; break;
      default:
        p1 = gmul(gpuigs(h,degq),p1);
        h = gdivexact(gpuigs(g,degq), gpuigs(h,degq-1));
    }
    if ((du & 1) == 0 && (dv & 1) == 0) signh= -signh;
    v=gdivexact(r,p1); uze=gdivexact(uze,p1);
    if (dr==3) break;
  }

  p2=(dv==4)?gun:gpuigs(gdiv((GEN)v[2],h),dv-4);
  if (p3) p2 = gmul(p2,gpuigs(p3,dy-3));
  if (p4) p2 = gmul(p2,gpuigs(p4,dx-3));
  if (signh<0) p2=gneg_i(p2);
  p3 = p3? gdiv(p2,p3): p2;

  tetpil=avma; z=gmul((GEN)v[2],p2); uze=gmul(uze,p3);
  gptr[0]=&z; gptr[1]=&uze; gerepilemanysp(av,tetpil,gptr,2);

  av=avma; p1 = gadd(z, gneg(gmul(uze,x))); tetpil = avma;
  if (!poldivis(p1,y,&p1)) err(bugparier,"subresext");
  *U=uze; *V=gerepile(av,tetpil,p1); return z;
}

static GEN
scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
{
  long v = varn(x);
  *U=gzero; *V=gdiv(polun[v],y); return polun[v];
}

static GEN
zero_bezout(GEN y, GEN *U, GEN *V)
{
  GEN x=content(y);
  *U=gzero; *V = gdiv(polun[varn(y)],x); return gmul(y,*V);
}

/* calcule U et V tel que Ux+Vy=GCD(x,y) par le sous-resultant */
GEN
bezoutpol(GEN x, GEN y, GEN *U, GEN *V)
{
  long degq,av,tetpil,tx,ty,dx,dy,du,dv,dr;
  GEN g,h,r,p1,p2,p3,p4,u,v,lpu,um1,uze,vze,xprim,yprim;
  GEN *gptr[3];

  if (gcmp0(x)) return zero_bezout(y,U,V);
  if (gcmp0(y)) return zero_bezout(x,V,U);
  tx=typ(x); ty=typ(y); av=avma;
  if (is_scalar_t(tx) || is_scalar_t(ty))
  {
    if (tx==t_POL) return scalar_bezout(x,y,U,V);
    if (ty==t_POL) return scalar_bezout(y,x,V,U);
    *U=ginv(x); *V=gzero; return polun[0];
  }
  if (tx!=t_POL || ty!=t_POL) err(typeer,"bezoutpol");
  if (varn(x)!=varn(y))
    return (varn(x)<varn(y))? scalar_bezout(x,y,U,V)
                            : scalar_bezout(y,x,V,U);
  dx=lgef(x); dy=lgef(y);
  if (dx<dy)
  {
    GEN *W=U; U=V; V=W;
    p1=x; x=y; y=p1; du=dx; dx=dy; dy=du;
  }
  if (dy==3) return scalar_bezout(x,y,U,V);

  p3=content(x); u=gdiv(x,p3);
  p4=content(y); v=gdiv(y,p4);
  xprim=u; yprim=v; g=gun; h=gun; um1=gun; uze=gzero;
  for(;;)
  {
    du=lgef(u); dv=lgef(v); degq=du-dv;
    lpu=gpuigs((GEN)v[dv-1],degq+1);
    p1=gmul(lpu,u); p2=poldivres(p1,v,&r);
    dr=lgef(r); if (dr<=2) break;
    p1=gsub(gmul(lpu,um1),gmul(p2,uze)); um1=uze; uze=p1;

    u=v; p1 = g; g  = leading_term(u);
    switch(degq)
    {
      case 0: break;
      case 1:
        p1 = gmul(h,p1); h = g; break;
      default:
        p1 = gmul(gpuigs(h,degq), p1);
        h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
    }
    v=gdiv(r,p1); uze=gdiv(uze,p1);
    if (dr==3) break;
  }
  if (!poldivis(gsub(v,gmul(uze,xprim)),yprim, &vze))
    err(warner,"non-exact computation in bezoutpol");
  uze=gdiv(uze,p3); vze=gdiv(vze,p4); p1=ginv(content(v));

  tetpil=avma; uze=gmul(uze,p1); vze=gmul(vze,p1); p1=gmul(v,p1);
  gptr[0]=&uze; gptr[1]=&vze; gptr[2]=&p1;
  gerepilemanysp(av,tetpil,gptr,3);
  *U=uze; *V=vze; return p1;
}



/*******************************************************************/
/*                                                                 */
/*               RESULTANT PAR L'ALGORITHME DE DUCOS               */
/*                                                                 */
/*******************************************************************/
GEN addshiftw(GEN x, GEN y, long d);

static GEN
reductum(GEN P)
{
  if (signe(P)==0) return P;
  return normalizepol_i(dummycopy(P),lgef(P)-1);
}

static GEN
Lazard(GEN x, GEN y, long n)
{
  long a, b;
  GEN c;
  
  if (n==1) return x;
  a=1; while (n >= (b=a+a)) a=b;
  c=x; n-=a;
  while (a>1)
  {
    a>>=1; c=gdivexact(gsqr(c),y);
    if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
  }
  return c;
}

static GEN
Lazard2(GEN F, GEN x, GEN y, long n)
{
  if (n<=1) return F;
  x = Lazard(x,y,n-1); return gdivexact(gmul(x,F),y);
}

static GEN
addshift(GEN x, GEN y)
{
  long v = varn(x);
  if (!signe(x)) return y;
  x = addshiftw(x,y,1);
  setvarn(x,v); return x;
}

static GEN
nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
{
  GEN p0, q0, z0, H, A;
  long p, q, j, av, lim, v = varn(P);

  z0 = leading_term(Z);
  p = degree(P); p0 = leading_term(P); P = reductum(P);
  q = degree(Q); q0 = leading_term(Q); Q = reductum(Q);

  av = avma; lim = stack_lim(av,1);
  H = gneg(reductum(Z));
  A = gmul((GEN)P[q+2],H);
  for (j = q+1; j < p; j++)
  {
    H = (degree(H) == q-1) ?
      addshift(reductum(H), gdivexact(gmul(gneg((GEN)H[q+1]),Q), q0)) :
      addshift(H, zeropol(v));
    A = gadd(A,gmul((GEN)P[j+2],H));
    if (low_stack(lim,stack_lim(av,1)))
    {
      GEN *gptr[2]; gptr[0]=&A; gptr[1]=&H;
      if(DEBUGMEM>1) err(warnmem,"nextSousResultant j = %ld/%ld",j,p);
      gerepilemany(av,gptr,2);
    }
  }
  P = normalizepol_i(P, q+2);
  A = gdivexact(gadd(A,gmul(z0,P)), p0);
  A = (degree(H) == q-1) ?
    gadd(gmul(q0,addshift(reductum(H),A)), gmul(gneg((GEN)H[q+1]),Q)) :
    gmul(q0, addshift(H,A));
  return gdivexact(A, ((p-q)&1)? s: gneg(s));
}

GEN
resultantducos(GEN P, GEN Q)
{
  long delta, av=avma, tetpil, lim = stack_lim(av,1);
  GEN Z, s;
  
  if ((Z = init_resultant(P,Q))) return Z;
  delta = degree(P) - degree(Q);
  if (delta < 0)
  {
    Z = ((degree(P)&1) && (degree(Q)&1)) ? gneg(Q) : Q;
    Q = P; P = Z; delta = -delta;
  }
  s = gun;
  if (degree(Q) > 0)
  {
    s = gpuigs(leading_term(Q),delta);
    Z = Q;
    Q = pseudorem(P, gneg(Q));
    P = Z;
    while(degree(Q) > 0)
    {
      if (low_stack(lim,stack_lim(av,1)))
      {
        GEN *gptr[2]; gptr[0]=&P; gptr[1]=&Q;
        if(DEBUGMEM>1) err(warnmem,"resultantducos, deg Q = %ld",degree(Q));
        gerepilemany(av,gptr,2); s = leading_term(P);
      }
      delta = degree(P) - degree(Q);
      Z = Lazard2(Q, leading_term(Q), s, delta);
      Q = nextSousResultant(P, Q, Z, s);
      P = Z;
      s = leading_term(P);
    }
  }
  if (!signe(Q)) { avma = av; return gzero; }
  s = Lazard(leading_term(Q), s, degree(P));
  tetpil = avma; return gerepile(av,tetpil,gcopy(s));
}

/*******************************************************************/
/*                                                                 */
/*               RESULTANT PAR MATRICE DE SYLVESTER                */
/*                                                                 */
/*******************************************************************/
GEN
sylvestermatrix_i(GEN x, GEN y)
{
  long d,dx,dy,i,j;
  GEN p1,p2;

  dx=lgef(x)-3; dy=lgef(y)-3; d=dx+dy;
  p1=cgetg(d+1,t_MAT);
  for (j=1; j<=dy; j++)
  {
    p2=cgetg(d+1,t_COL); p1[j]=(long)p2;
    for (i=1; i<j; i++) p2[i]=zero;
    for (   ; i<=j+dx; i++) p2[i]=x[dx-i+j+2];
    for (   ; i<=d; i++) p2[i]=zero;
  }
  for (j=1; j<=dx; j++)
  {
    p2=cgetg(d+1,t_COL); p1[j+dy]=(long)p2;
    for (i=1; i<j; i++) p2[i]=zero;
    for (   ; i<=j+dy; i++) p2[i]=y[dy-i+j+2];
    for (   ; i<=d; i++) p2[i]=zero;
  }
  return p1;
}

GEN 
sylvestermatrix(GEN x, GEN y)
{
  long i,j,lx;
  if (typ(x)!=t_POL || typ(y)!=t_POL) err(typeer,"sylvestermatrix");
  if (varn(x) != varn(y))
    err(talker,"not the same variables in sylvestermatrix");
  x = sylvestermatrix_i(x,y); lx = lg(x);
  for (i=1; i<lx; i++)
    for (j=1; j<lx; j++) coeff(x,i,j) = lcopy(gcoeff(x,i,j));
  return x;
}

GEN
resultant2(GEN x, GEN y)
{
  long av;
  GEN r;
  if ((r = init_resultant(x,y))) return r;
  av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
}

static GEN
fix_pol(GEN x, long v, long *mx)
{
  long vx;
  GEN p1;

  if (typ(x) == t_POL)
  {
    vx = varn(x);
    if (vx)
    {
      if (v>=vx) return gsubst(x,v,polx[0]);
      p1 = cgetg(3,t_POL);
      p1[1] = evalvarn(0)|evalsigne(signe(x))|evallgef(3);
      p1[2] = (long)x; return p1;
    }
    if (v)
    {
      *mx = 1;
      return gsubst(gsubst(x,0,polx[MAXVARN]),v,polx[0]);
    }
  }
  return x;
}

/* resultant of x and y with respect to variable v, or with respect to their
 * main variable if v < 0.
 */
GEN
polresultant0(GEN x, GEN y, long v, long flag)
{
  long av = avma, m = 0;

  if (v >= 0)
  {
    x = fix_pol(x,v, &m);
    y = fix_pol(y,v, &m);
  }
  switch(flag)
  {
    case 0: x=subresall(x,y,NULL); break;
    case 1: x=resultant2(x,y); break;
    case 2: x=resultantducos(x,y); break;
    default: err(flagerr,"polresultant");
  }
  if (m) x = gsubst(x,MAXVARN,polx[0]);
  return gerepileupto(av,x);
}

/*******************************************************************/
/*                                                                 */
/*           P.G.C.D PAR L'ALGORITHME DU SOUS RESULTANT            */
/*                                                                 */
/*******************************************************************/

GEN
srgcd(GEN x, GEN y)
{
  long av,av1,tetpil,tx=typ(x),ty=typ(y),dx,dy,vx,lim;
  GEN d,g,h,p1,p2,u,v;

  if (!signe(y)) return gcopy(x);
  if (!signe(x)) return gcopy(y);
  if (is_scalar_t(tx) || is_scalar_t(ty)) return gun;
  if (tx!=t_POL || ty!=t_POL) err(typeer,"srgcd");
  vx=varn(x); if (vx!=varn(y)) return gun;
  if (ismonome(x)) return gcdmonome(x,y);
  if (ismonome(y)) return gcdmonome(y,x);
  if (isrational(x) && isrational(y)) return modulargcd(x,y);

  av=avma;
  if (issimplefield(x) || issimplefield(y)) x = polgcdnun(x,y);
  else
  {
    dx=lgef(x); dy=lgef(y);
    if (dx<dy) { p1=x; x=y; y=p1; tx=dx; dx=dy; dy=tx; }
    p1=content(x); p2=content(y); d=ggcd(p1,p2);

    tetpil=avma; d=gmul(d,polun[vx]);
    if (dy==3) return gerepile(av,tetpil,d);

    av1=avma; lim=stack_lim(av1,1);
    u=gdiv(x,p1); v=gdiv(y,p2); g=h=gun;
    for(;;)
    {
      GEN r = pseudorem(u,v);
      long degq, du, dv, dr=lgef(r);

      if (dr<=3)
      {
        if (gcmp0(r)) break;
        avma=av1; return gerepile(av,tetpil,d);
      }
      if (DEBUGLEVEL > 9) fprintferr("srgcd: dr = %ld\n", dr);
      du=lgef(u); dv=lgef(v); degq=du-dv; u=v;
      switch(degq)
      {
        case 0:
          v = gdiv(r,g);
          g = leading_term(u);
          break;
        case 1:
          v = gdiv(r,gmul(h,g));
          h = g = leading_term(u);
          break;
        default:
          v = gdiv(r,gmul(gpuigs(h,degq),g));
          g = leading_term(u);
          h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
      }
      if (low_stack(lim, stack_lim(av1,1)))
      {
        GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
        if(DEBUGMEM>1) err(warnmem,"srgcd");
        gerepilemany(av1,gptr,4);
      }
    }
    p1 = content(v); if (!gcmp1(p1)) v = gdiv(v,p1);
    x = gmul(d,v);
  }

  if (typ(x)!=t_POL) x = gmul(polun[vx],x);
  else
  {
    p1=leading_term(x); ty=typ(p1);
    if ((is_frac_t(ty) || is_intreal_t(ty)) && gsigne(p1)<0) x = gneg(x);
  }
  return gerepileupto(av,x);
}

GEN qf_disc(GEN x, GEN y, GEN z);

GEN
poldisc0(GEN x, long v)
{
  long av,tx=typ(x),i;
  GEN z,p1,p2;

  switch(tx)
  {
    case t_POL:
      if (gcmp0(x)) return gzero;
      av = avma; i = 0;
      if (v >= 0 && v != varn(x)) x = fix_pol(x,v, &i);
      p1 = subres(x, derivpol(x));
      p2 = leading_term(x); if (!gcmp1(p2)) p1 = gdiv(p1,p2);
      if ((lgef(x)-3) & 2) p1 = gneg(p1);
      if (i) p1 = gsubst(p1, MAXVARN, polx[0]);
      return gerepileupto(av,p1);

    case t_COMPLEX:
      return stoi(-4);

    case t_QUAD: case t_POLMOD:
      return poldisc0((GEN)x[1], v);

    case t_QFR: case t_QFI:
      av = avma; return gerepileuptoint(av, qf_disc(x,NULL,NULL));

    case t_VEC: case t_COL: case t_MAT:
      i=lg(x); z=cgetg(i,tx);
      for (i--; i; i--) z[i]=(long)poldisc0((GEN)x[i], v);
      return z;
  }
  err(typeer,"discsr");
  return NULL; /* not reached */
}

GEN
discsr(GEN x)
{
  return poldisc0(x, -1);
}

GEN
reduceddiscsmith(GEN pol)
{
  long av=avma,tetpil,i,j,n;
  GEN polp,alpha,p1,m;

  if (typ(pol)!=t_POL) err(typeer,"reduceddiscsmith");
  n=lgef(pol)-3;
  if (n<=0) err(constpoler,"reduceddiscsmith");
  check_pol_int(pol);
  if (!gcmp1((GEN)pol[n+2]))
    err(talker,"non-monic polynomial in poldiscreduced");
  polp = derivpol(pol); alpha = polx[varn(pol)];
  m=cgetg(n+1,t_MAT);
  for (j=1; j<=n; j++)
  {
    p1=cgetg(n+1,t_COL); m[j]=(long)p1;
    for (i=1; i<=n; i++) p1[i]=(long)truecoeff(polp,i-1);
    if (j<n) polp = gres(gmul(alpha,polp), pol);
  }
  tetpil=avma; return gerepile(av,tetpil,smith(m));
}

/***********************************************************************/
/**							              **/
/**	                 ALGORITHME DE STURM                          **/
/**	         (number of real roots of x in ]a,b])		      **/
/**								      **/
/***********************************************************************/

/* if a (resp. b) is NULL, set it to -oo (resp. +oo) */
long
sturmpart(GEN x, GEN a, GEN b)
{
  long av = avma,sl,sr,s,t,r1;
  GEN g,h,u,v;

  if (typ(x) != t_POL) err(typeer,"sturm");
  if (gcmp0(x)) err(zeropoler,"sturm");
  s=lgef(x); if (s==3) return 0;

  sl = gsigne(leading_term(x));
  if (s==4)
  {
    s = b? gsigne(poleval(x,b)):  sl;
    t = a? gsigne(poleval(x,a)): -sl;
    avma = av; return (s == t)? 0: 1;
  }
  u=gdiv(x,content(x)); v=derivpol(x);
  v=gdiv(v,content(v));
  g=gun; h=gun;
  s = b? gsigne(poleval(u,b)): sl;
  t = a? gsigne(poleval(u,a)): ((lgef(u)&1)? sl: -sl);
  r1=0;
  sr = b? gsigne(poleval(v,b)): s;
  if (sr)
  {
    if (!s) s=sr;
    else if (sr!=s) { s= -s; r1--; }
  }
  sr = a? gsigne(poleval(v,a)): -t;
  if (sr)
  {
    if (!t) t=sr;
    else if (sr!=t) { t= -t; r1++; }
  }
  for(;;)
  {
    GEN p1, r = pseudorem(u,v);
    long du=lgef(u), dv=lgef(v), dr=lgef(r), degq=du-dv;

    if (dr<=2) err(talker,"not a squarefree polynomial in sturm");
    if (gsigne(leading_term(v)) > 0 || degq&1) r=gneg_i(r);
    sl = gsigne((GEN) r[dr-1]);
    sr = b? gsigne(poleval(r,b)): sl;
    if (sr)
    {
      if (!s) s=sr;
      else if (sr!=s) { s= -s; r1--; }
    }
    sr = a? gsigne(poleval(r,a)): ((dr&1)? sl: -sl);
    if (sr)
    {
      if (!t) t=sr;
      else if (sr!=t) { t= -t; r1++; }
    }
    if (dr==3) { avma=av; return r1; }

    u=v; p1 = g; g = gabs(leading_term(u),DEFAULTPREC);
    switch(degq)
    {
      case 0: break;
      case 1:
        p1 = gmul(h,p1); h = g; break;
      default:
        p1 = gmul(gpuigs(h,degq),p1);
        h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
    }
    v = gdiv(r,p1);
  }
}

/*******************************************************************/
/*                                                                 */
/*         POLYNOME QUADRATIQUE ASSOCIE A UN DISCRIMINANT          */
/*                                                                 */
/*******************************************************************/

GEN
quadpoly0(GEN x, long v)
{
  long res,l,tetpil,i,sx, tx = typ(x);
  GEN y,p1;

  if (is_matvec_t(tx))
  {
    l=lg(x); y=cgetg(l,tx);
    for (i=1; i<l; i++) y[i]=(long)quadpoly0((GEN)x[i],v);
    return y;
  }
  if (tx!=t_INT) err(arither1);
  if (v < 0) v = 0;
  sx=signe(x);
  if (!sx) err(talker,"zero discriminant in quadpoly");
  y=cgetg(5,t_POL);
  y[1]=evalsigne(1) | evalvarn(v) | evallgef(5); y[4]=un;
  res=mod4(x); if (sx<0 && res) res=4-res;
  if (res>1) err(funder2,"quadpoly");

  l=avma; p1=shifti(x,-2); setsigne(p1,-signe(p1));
  y[2] = (long) p1;
  if (!res) y[3] = zero;
  else
  {
    if (sx<0) { tetpil=avma; y[2] = lpile(l,tetpil,addsi(1,p1)); }
    y[3] = lnegi(gun);
  }
  return y;
}

GEN
quadpoly(GEN x)
{
  return quadpoly0(x,-1);
}

GEN
quadgen(GEN x)
{
  GEN y=cgetg(4,t_QUAD);
  y[1]=lquadpoly(x); y[2]=zero; y[3]=un; return y;
}

/*******************************************************************/
/*                                                                 */
/*                    INVERSE MODULO GENERAL                       */
/*                                                                 */
/*******************************************************************/

GEN
ginvmod(GEN x, GEN y)
{
  long tx=typ(x);

  switch(typ(y))
  {
    case t_POL:
      if (tx==t_POL) return polinvmod(x,y);
      if (is_scalar_t(tx)) return ginv(x);
      break;
    case t_INT:
      if (tx==t_INT) return mpinvmod(x,y);
      if (tx==t_POL) return gzero;
  }
  err(typeer,"ginvmod");
  return NULL; /* not reached */
}

/***********************************************************************/
/**							              **/
/**		            NEWTON POLYGON                            **/
/**								      **/
/***********************************************************************/

/* assume leading coeff of x is non-zero */
GEN
newtonpoly(GEN x, GEN p)
{
  GEN y;
  long n,ind,a,b,c,u1,u2,r1,r2;
  long *vval, num[] = {evaltyp(t_INT) | m_evallg(3), 0, 0};

  if (typ(x)!=t_POL) err(typeer,"newtonpoly");
  n=lgef(x)-3; if (n<=0) { y=cgetg(1,t_VEC); return y; }
  y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
  vval = (long *) gpmalloc(sizeof(long)*(n+1));
  for (a=0; a<=n; a++) vval[a] = ggval((GEN)x[a],p);
  for (a=0, ind=1; a<n; a++)
  {
    if (vval[a] != VERYBIGINT) break;
    y[ind++] = lstoi(VERYBIGINT);
  }
  for (b=a+1; b<=n; a=b, b=a+1)
  {
    while (vval[b] == VERYBIGINT) b++;
    u1=vval[a]-vval[b]; u2=b-a;
    for (c=b+1; c<=n; c++)
    {
      if (vval[c] == VERYBIGINT) continue;
      r1=vval[a]-vval[c]; r2=c-a;
      if (u1*r2<=u2*r1) { u1=r1; u2=r2; b=c; }
    }
    while (ind<=b) { affsi(u1,num); y[ind++] = ldivgs(num,u2); }
  }
  free(vval); return y;
}

int cmp_pol(GEN x, GEN y);

/* Factor polynomial a on the number field defined by polynomial t */
GEN
polfnf(GEN a, GEN t)
{
  GEN x0, y,p1,p2,u,g,fa,n,unt;
  long av=avma,tetpil,lx,v,i,k,vt;

  if (typ(a)!=t_POL || typ(t)!=t_POL) err(typeer,"polfnf");
  if (gcmp0(a)) return gcopy(a);
  vt=varn(t); v=varn(a);
  if (vt<=v)
    err(talker,"polynomial variable must be of higher priority than number field variable\nin factornf");
  u = gdiv(a, ggcd(a,derivpol(a)));
  unt = gmodulsg(1,t); u = gmul(unt,u);
  g = lift(u);
  for (k=-1; ; k++)
  {
    n = gsub(polx[MAXVARN], gmulsg(k,polx[vt]));
    n = subres(t, poleval(g,n));
    if (lgef(ggcd(n,derivpol(n))) == 3) break;
  }
  fa=factor(n); fa=(GEN)fa[1]; lx=lg(fa);
  y=cgetg(3,t_MAT);
  p1=cgetg(lx,t_COL); y[1]=(long)p1;
  p2=cgetg(lx,t_COL); y[2]=(long)p2;
  x0 = gadd(polx[v],gmulsg(k,gmodulcp(polx[vt],t)));
  for (i=1; i<lx; i++)
  {
    GEN b, pro = ggcd(u, gmul(unt, poleval((GEN)fa[i], x0)));
    long e;

    p1[i] = (typ(pro)==t_POL)? ldiv(pro,leading_term(pro)): (long)pro;
    if (gcmp1((GEN)p1[i])) err(talker,"reducible modulus in factornf");
    e=0; while (poldivis(a,(GEN)p1[i], &b)) { a = b; e++; }
    p2[i] = lstoi(e);
  }
  (void)sort_factor(y, cmp_pol);
  tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}