/* Copyright (C) 2000  The PARI group.

This file is part of the PARI/GP package.

PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.

Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */

/***********************************************************************/
/**                                                                   **/
/**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
/**                         (second part)                             **/
/**                                                                   **/
/***********************************************************************/
#include "pari.h"
#include "paripriv.h"

#define addshift(x,y) addshiftpol((x),(y),1)

/* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
 * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
 * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
 * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
 * Not memory clean in the latter case */
GEN
polsym_gen(GEN P, GEN y0, pari_long n, GEN T, GEN N)
{
  pari_long dP=degpol(P), i, k, m;
  pari_sp av1, av2;
  GEN s,y,P_lead;

  if (n<0) pari_err_IMPL("polsym of a negative n");
  if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
  if (!signe(P)) pari_err_ROOTS0("polsym");
  y = cgetg(n+2,t_COL);
  if (y0)
  {
    if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
    m = lg(y0)-1;
    for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
  }
  else
  {
    m = 1;
    gel(y,1) = stoi(dP);
  }
  P += 2; /* strip codewords */

  P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
  if (P_lead)
  {
    if (N) P_lead = Fq_inv(P_lead,T,N);
    else if (T) P_lead = QXQ_inv(P_lead,T);
  }
  for (k=m; k<=n; k++)
  {
    av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
    for (i=1; i<k && i<=dP; i++)
      s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
    if (N)
    {
      s = Fq_red(s, T, N);
      if (P_lead) s = Fq_mul(s, P_lead, T, N);
    }
    else if (T)
    {
      s = grem(s, T);
      if (P_lead) s = grem(gmul(s, P_lead), T);
    }
    else
      if (P_lead) s = gdiv(s, P_lead);
    av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
  }
  return y;
}

GEN
polsym(GEN x, pari_long n)
{
  return polsym_gen(x, NULL, n, NULL,NULL);
}

/* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
GEN
centermodii(GEN x, GEN p, GEN po2)
{
  GEN y = remii(x, p);
  switch(signe(y))
  {
    case 0: break;
    case 1: if (po2 && absi_cmp(y,po2) > 0) y = subii(y, p);
      break;
    case -1: if (!po2 || absi_cmp(y,po2) > 0) y = addii(y, p);
      break;
  }
  return y;
}

static pari_long
s_centermod(pari_long x, ulong pp, ulong pps2)
{
  pari_long y = x % (pari_long)pp;
  if (y < 0) y += pp;
  return Fl_center(y, pp,pps2);
}

/* for internal use */
GEN
centermod_i(GEN x, GEN p, GEN ps2)
{
  pari_long i, lx;
  pari_sp av;
  GEN y;

  if (!ps2) ps2 = shifti(p,-1);
  switch(typ(x))
  {
    case t_INT: return centermodii(x,p,ps2);

    case t_POL: lx = lg(x);
      y = cgetg(lx,t_POL); y[1] = x[1];
      for (i=2; i<lx; i++)
      {
        av = avma;
        gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
      }
      return normalizepol_lg(y, lx);

    case t_COL: lx = lg(x);
      y = cgetg(lx,t_COL);
      for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
      return y;

    case t_MAT: lx = lg(x);
      y = cgetg(lx,t_MAT);
      for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
      return y;

    case t_VECSMALL: lx = lg(x);
    {
      ulong pp = itou(p), pps2 = itou(ps2);
      y = cgetg(lx,t_VECSMALL);
      for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
      return y;
    }
  }
  return x;
}

GEN
centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }

/***********************************************************************/
/**                                                                   **/
/**                          FACTORIZATION                            **/
/**                                                                   **/
/***********************************************************************/
#define assign_or_fail(x,y) { GEN __x = x;\
  if (y==NULL) y=__x; else if (!gequal(__x,y)) return 0;\
}

static const pari_long tsh = 6;
static pari_long
RgX_type_code(pari_long t1, pari_long t2) { return (t1 << tsh) | t2; }
void
RgX_type_decode(pari_long x, pari_long *t1, pari_long *t2)
{
  *t1 = x >> tsh;
  *t2 = (x & ((1LL<<tsh)-1));
}
int
RgX_type_is_composite(pari_long t) { return t >= tsh; }

pari_long
RgX_type(GEN x, GEN *ptp, GEN *ptpol, pari_long *ptpa)
{
  pari_long t[16];
  pari_long tx = typ(x), lx, i, j, s, pa = LONG_MAX;
  GEN pcx=NULL, p=NULL, pol=NULL, ff=NULL;

  if (is_scalar_t(tx))
  {
    if (tx == t_POLMOD) return 0;
    x = scalarpol(x,0);
  }
  for (i=2; i<16; i++) t[i]=0;
  /* t[0..1] unused. Other values, if set, indicate a coefficient of type
   * t[2] : t_REAL
   * t[3] : t_INTMOD
   * t[4] : t_COMPLEX of rationals (t_INT/t_FRAC)
   * t[5] : t_COMPLEX of t_REAL
   * t[6] : t_COMPLEX of t_INTMOD
   * t[7] : t_COMPLEX of t_PADIC
   * t[8] : t_PADIC
   * t[9] : t_QUAD of rationals (t_INT/t_FRAC)
   * t[10]: t_QUAD of t_INTMOD
   * t[11]: t_QUAD of t_PADIC
   * t[12]: t_POLMOD of rationals (t_INT/t_FRAC)
   * t[13]: t_POLMOD of t_INTMOD
   * t[14]: t_POLMOD of t_PADIC
   * t[15]: t_FFELT */
  lx = lg(x);
  for (i=2; i<lx; i++)
  {
    GEN c = gel(x,i);
    switch(typ(c))
    {
      case t_INT: case t_FRAC:
        break;
      case t_REAL:
        s = precision(c); if (s < pa) pa = s;
        t[2]=1; break;
      case t_INTMOD:
        assign_or_fail(gel(c,1),p);
        t[3]=1; break;
      case t_FFELT:
        if (ff==NULL) ff=c;
        else if (!FF_samefield(c,ff)) return 0;
        assign_or_fail(FF_p_i(c),p);
        t[15]=1; break;
      case t_COMPLEX:
        if (!pcx) pcx = mkpoln(3, gen_1,gen_0,gen_1); /* x^2 + 1 */
        for (j=1; j<=2; j++)
        {
          GEN d = gel(c,j);
          switch(typ(d))
          {
            case t_INT: case t_FRAC:
              assign_or_fail(pcx,pol);
              t[4]=1; break;
            case t_REAL:
              s = precision(d); if (s < pa) pa = s;
              t[5]=1; break;
            case t_INTMOD:
              assign_or_fail(gel(d,1),p);
              if (!signe(p) || mod4(p) != 3) return 0;
              assign_or_fail(pcx,pol);
              t[6]=1; break;
            case t_PADIC:
              s = precp(d) + valp(d); if (s < pa) pa = s;
              assign_or_fail(gel(d,2),p);
              assign_or_fail(pcx,pol);
              t[7]=1; break;
            default: return 0;
          }
        }
        break;
      case t_PADIC:
        s = precp(c) + valp(c); if (s < pa) pa = s;
        assign_or_fail(gel(c,2),p);
        t[8]=1; break;
      case t_QUAD:
        for (j=2; j<=3; j++)
        {
          GEN d = gel(c,j);
          switch(typ(d))
          {
            case t_INT: case t_FRAC:
              assign_or_fail(gel(c,1),pol);
              t[9]=1; break;
            case t_INTMOD:
              assign_or_fail(gel(d,1),p);
              assign_or_fail(gel(c,1),pol);
              t[10]=1; break;
            case t_PADIC:
              s = precp(d) + valp(d); if (s < pa) pa = s;
              assign_or_fail(gel(d,2),p);
              assign_or_fail(gel(c,1),pol);
              t[11]=1; break;
            default: return 0;
          }
        }
        break;
      case t_POLMOD:
        assign_or_fail(gel(c,1),pol);
        for (j=1; j<=2; j++)
        {
          GEN pbis = NULL, polbis = NULL;
          pari_long pabis;
          switch(RgX_type(gel(c,j),&pbis,&polbis,&pabis))
          {
            case t_INT: t[12]=1; break;
            case t_INTMOD: t[13]=1; break;
            case t_PADIC: t[14]=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])
  {
    if (t[3]||t[6]||t[7]||t[8]||t[10]||t[11]||t[12]||t[13]||t[14]) return 0;
    *ptpa=pa; return t_COMPLEX;
  }
  if (t[2])
  {
    if (t[3]||t[6]||t[7]||t[8]||t[10]||t[11]||t[12]||t[13]||t[14]) return 0;
    *ptpa=pa; return t[4]?t_COMPLEX:t_REAL;
  }
  if (t[6]||t[10]||t[13])
  {
    *ptpol=pol; *ptp=p;
    i = t[13]? t_POLMOD: (t[10]? t_QUAD: t_COMPLEX);
    return RgX_type_code(i, t_INTMOD);
  }
  if (t[7]||t[11]||t[14])
  {
    *ptpol=pol; *ptp=p; *ptpa=pa;
    i = t[14]? t_POLMOD: (t[11]? t_QUAD: t_COMPLEX);
    return RgX_type_code(i, t_PADIC);
  }
  if (t[4]||t[9]||t[12])
  {
    *ptpol=pol;
    i = t[12]? t_POLMOD: (t[9]? t_QUAD: t_COMPLEX);
    return RgX_type_code(i, t_INT);
  }
  if (t[15])
  {
    if (t[8]) return 0;
    *ptp=p; *ptpol=ff;
    return t_FFELT;
  }
  if (t[3]) { *ptp=p; return t_INTMOD; }
  if (t[8]) { *ptp=p; *ptpa=pa; return t_PADIC; }
  return t_INT;
}

GEN
factor0(GEN x, pari_long flag)
{
  return (flag<0)? factor(x): boundfact(x,flag);
}

/* only present for interface with GP */
GEN
gp_factor0(GEN x, GEN flag)
{
  ulong B;
  if (!flag) return factor(x);
  if (typ(flag) != t_INT || signe(flag) < 0) pari_err_FLAG("factor");
  switch(lgefint(flag))
  {
    case 2: B = 0; break;
    case 3: B = flag[2]; break;
    default: pari_err_OVERFLOW("factor [large prime bound]");
             return NULL; /*not reached*/
  }
  return boundfact(x, B);
}

GEN
deg1_from_roots(GEN L, pari_long v)
{
  pari_long i, l = lg(L);
  GEN z = cgetg(l,t_COL);
  for (i=1; i<l; i++)
    gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
  return z;
}
GEN
roots_from_deg1(GEN x)
{
  pari_long i,l = lg(x);
  GEN r = cgetg(l,t_VEC);
  for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
  return r;
}

static GEN
gauss_factor_p(GEN p)
{
  GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
  return mkcomplex(a, b);
}

static GEN
gauss_primpart(GEN x, GEN *c)
{
  GEN a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
  *c = n; if (n == gen_1) return x;
  retmkcomplex(diviiexact(a,n), diviiexact(b,n));
}

static GEN
gauss_primpart_try(GEN x, GEN c)
{
  GEN r, y;
  if (typ(x) == t_INT)
  {
    y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
  }
  else
  {
    GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
    gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
    gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
  }
  return y;
}

static int
gauss_cmp(GEN x, GEN y)
{
  int v;
  if (typ(x) != t_COMPLEX)
    return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
  if (typ(y) != t_COMPLEX) return 1;
  v = cmpii(gel(x,2), gel(y,2));
  if (v) return v;
  return gcmp(gel(x,1), gel(y,1));
}

/* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
static GEN
gauss_normal(GEN x)
{
  if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
  if (signe(gel(x,1)) < 0) x = gneg(x);
  if (signe(gel(x,2)) < 0) x = mulcxI(x);
  return x;
}

static GEN
gauss_factor(GEN x)
{
  pari_sp av = avma;
  GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
  pari_long t1 = typ(a);
  pari_long t2 = typ(b), i, j, l, exp = 0;
  if (t1 == t_FRAC) d = gel(a,2);
  if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
  if (d == gen_1) y = x;
  else
  {
    y = gmul(x, d);
    a = gel(y,1); t1 = typ(a);
    b = gel(y,2); t2 = typ(b);
  }
  if (t1 != t_INT || t2 != t_INT) return NULL;
  y = gauss_primpart(y, &n);
  fa = factor(cxnorm(y));
  P = gel(fa,1);
  E = gel(fa,2); l = lg(P);
  P2 = cgetg(l, t_COL);
  E2 = cgetg(l, t_COL);
  for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
  { /* either p = 2 (ramified) or those factors split in Q(i) */
    GEN p = gel(P,i), w, w2, t, we, pe;
    pari_long v, e = itos(gel(E,i));
    int is2 = equaliu(p, 2);
    w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
    w2 = gauss_normal( gconj(w) );
    /* w * w2 * I^3 = p, w2 = gconj(w) * I */
    pe = powiu(p, e);
    we = gpowgs(w, e);
    t = gauss_primpart_try( gmul(y, gconj(we)), pe );
    if (t) y = t; /* y /= w^e */
    else {
      /* y /= conj(w)^e, should be y /= w2^e */
      y = gauss_primpart_try( gmul(y, we), pe );
      swap(w, w2); exp -= e; /* += 3*e mod 4 */
    }
    gel(P,i) = w;
    v = Z_pvalrem(n, p, &n);
    if (v) {
      exp -= v; /* += 3*v mod 4 */
      if (is2) v <<= 1; /* 2 = w^2 I^3 */
      else {
        gel(P2,j) = w2;
        gel(E2,j) = utoipos(v); j++;
      }
      gel(E,i) = stoi(e + v);
    }
    v = Z_pvalrem(d, p, &d);
    if (v) {
      exp += v; /* -= 3*v mod 4 */
      if (is2) v <<= 1; /* 2 is ramified */
      else {
        gel(P2,j) = w2;
        gel(E2,j) = utoineg(v); j++;
      }
      gel(E,i) = stoi(e - v);
    }
    exp &= 3;
  }
  if (j > 1) {
    pari_long k = 1;;
    GEN P1 = cgetg(l, t_COL);
    GEN E1 = cgetg(l, t_COL);
    /* remove factors with exponent 0 */
    for (i = 1; i < l; i++)
      if (signe(gel(E,i)))
      {
        gel(P1,k) = gel(P,i);
        gel(E1,k) = gel(E,i);
        k++;
      }
    setlg(P1, k); setlg(E1, k);
    setlg(P2, j); setlg(E2, j);
    fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
  }
  if (!is_pm1(n) || !is_pm1(d))
  {
    GEN Fa = factor(gdiv(n, d));
    P = gel(Fa,1); l = lg(P);
    E = gel(Fa,2);
    for (i = 1; i < l; i++)
    {
      GEN w, p = gel(P,i);
      pari_long e;
      int is2;
      switch(mod4(p))
      {
        case 3: continue;
        case 2: is2 = 1; break;
        default:is2 = 0; break;
      }
      e = itos(gel(E,i));
      w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
      gel(P,i) = w;
      if (is2)
        gel(E,i) = stoi(e << 1);
      else
      {
        P = shallowconcat(P, gauss_normal( gconj(w) ));
        E = shallowconcat(E, gel(E,i));
      }
      exp -= e; /* += 3*e mod 4 */
      exp &= 3;
    }
    gel(Fa,1) = P;
    gel(Fa,2) = E;
    fa = famat_mul_shallow(fa, Fa);
  }
  fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);

  y = gmul(y, powIs(exp));
  if (!gequal1(y)) {
    gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
    gel(fa,2) = shallowconcat(gen_1,    gel(fa,2));
  }
  return gerepilecopy(av, fa);
}

GEN
factor(GEN x)
{
  pari_long tx=typ(x), lx, i, pa, v, r1;
  pari_sp av, tetpil;
  GEN  y, p, p1, p2, pol;

  if (gequal0(x))
    switch(tx)
    {
      case t_INT:
      case t_COMPLEX:
      case t_POL:
      case t_RFRAC: return prime_fact(x);
      default: pari_err_TYPE("factor",x);
    }
  av = avma;
  switch(tx)
  {
    case t_INT: return Z_factor(x);

    case t_FRAC:
      p1 = Z_factor(gel(x,1));
      p2 = Z_factor(gel(x,2)); gel(p2,2) = gneg_i(gel(p2,2));
      return gerepilecopy(av, merge_factor_i(p1,p2));

    case t_POL:
      tx=RgX_type(x,&p,&pol,&pa);
      switch(tx)
      {
        case 0: pari_err_IMPL("factor for general polynomials");
        case t_INT: return QX_factor(x);
        case t_INTMOD: return factmod(x,p);

        case t_COMPLEX: y=cgetg(3,t_MAT); lx=lg(x)-2;
          av = avma; p1 = roots(x,pa); tetpil = avma;
          p1 = deg1_from_roots(p1, varn(x));
          gel(y,1) = gerepile(av,tetpil,p1);
          gel(y,2) = const_col(lx-1, gen_1); return y;

        case t_REAL: y=cgetg(3,t_MAT); lx=lg(x)-2; v=varn(x);
          av=avma; p1=cleanroots(x,pa); tetpil=avma;
          for(r1=1; r1<lx; r1++)
            if (typ(gel(p1,r1)) == t_COMPLEX) break;
          lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
          for(i=1; i<r1; i++)
            gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
          for(   ; i<lx; i++)
          {
            GEN a = gel(p1,2*i-r1);
            p = cgetg(5, t_POL); gel(p2,i) = p;
            p[1] = x[1];
            gel(p,2) = gnorm(a);
            gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
            gel(p,4) = gen_1;
          }
          gel(y,1) = gerepile(av,tetpil,p2);
          gel(y,2) = const_col(lx-1, gen_1); return y;

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

        case t_FFELT: return FFX_factor(x,pol);

        default:
        {
          GEN w;
          pari_long killv, t1, t2;
          x = leafcopy(x); lx=lg(x);
          pol = leafcopy(pol);
          v = pari_var_next_temp();
          for(i=2; i<lx; i++)
          {
            p1 = gel(x,i);
            switch(typ(p1))
            {
              case t_QUAD: p1++;
              case t_COMPLEX:
                gel(x,i) = mkpolmod(deg1pol_shallow(gel(p1,2), gel(p1,1), v), pol);
            }
          }
          killv = (avma != (pari_sp)pol);
          if (killv) setvarn(pol, fetch_var());
          RgX_type_decode(tx, &t1, &t2);
          switch (t2)
          {
            case t_INT: p1 = polfnf(x,pol); break;
            case t_INTMOD:
              pol = RgX_to_FpX(pol, p);
              if (FpX_is_squarefree(pol,p) && FpX_nbfact(pol, p) == 1)
              {
                p1 = factorff(x,p,pol); break;
              }
            /*fall through*/
            default: pari_err_IMPL("factor for general polynomial");
              return NULL; /* not reached */
          }
          switch (t1)
          {
            case t_POLMOD:
              if (killv) (void)delete_var();
              return gerepileupto(av,p1);
            case t_COMPLEX: w = gen_I(); break;
            case t_QUAD: w = mkquad(pol,gen_0,gen_1);
              break;
            default: pari_err_IMPL("factor for general polynomial");
              return NULL; /* not reached */
          }
          p2=gel(p1,1);
          for(i=1; i<lg(p2); i++)
          {
            GEN P = gel(p2,i);
            pari_long j;
            for(j=2; j<lg(P); j++)
            {
              GEN p = gel(P,j);
              if(typ(p)==t_POLMOD) gel(P,j) = gsubst(gel(p,2),v,w);
            }
          }
          if (killv) (void)delete_var();
          return gerepilecopy(av, p1);
        }
      }
    case t_RFRAC: {
      GEN a = gel(x,1), b = gel(x,2);
      y = factor(b); gel(y,2) = gneg_i(gel(y,2));
      if (typ(a)==t_POL && varn(a)==varn(b)) y = famat_mul_shallow(factor(a), y);
      return gerepilecopy(av, y);
    }

    case t_COMPLEX:
      y = gauss_factor(x);
      if (y) return y;
      /* fall through */
  }
  pari_err_TYPE("factor",x);
  return NULL; /* not reached */
}

/*******************************************************************/
/*                                                                 */
/*                     ROOTS --> MONIC POLYNOMIAL                  */
/*                                                                 */
/*******************************************************************/
static GEN
normalized_mul(GEN x, GEN y)
{
  pari_long a = gel(x,1)[1], b = gel(y,1)[1];
  return mkvec2(mkvecsmall(a + b),
                RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
}
/* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
static GEN
normalized_to_RgX(GEN L)
{
  pari_long i, a = gel(L,1)[1];
  GEN A = gel(L,2);
  GEN z = cgetg(a + 3, t_POL);
  z[1] = evalsigne(1) | evalvarn(varn(A));
  for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
  for (     ; i < a+2;   i++) gel(z,i) = gen_0;
  gel(z,i) = gen_1; return z;
}

/* compute prod (x - a[i]) */
GEN
roots_to_pol(GEN a, pari_long v)
{
  pari_sp av = avma;
  pari_long i, k, lx = lg(a);
  GEN L;
  if (lx == 1) return pol_1(v);
  L = cgetg(lx, t_VEC);
  for (k=1,i=1; i<lx-1; i+=2)
  {
    GEN s = gel(a,i), t = gel(a,i+1);
    GEN x0 = gmul(s,t);
    GEN x1 = gneg(gadd(s,t));
    gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
  }
  if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
                                  scalarpol_shallow(gneg(gel(a,i)), v));
  setlg(L, k); L = divide_conquer_prod(L, normalized_mul);
  return gerepileupto(av, normalized_to_RgX(L));
}

/* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
GEN
roots_to_pol_r1(GEN a, pari_long v, pari_long r1)
{
  pari_sp av = avma;
  pari_long i, k, lx = lg(a);
  GEN L;
  if (lx == 1) return pol_1(v);
  L = cgetg(lx, t_VEC);
  for (k=1,i=1; i<r1; i+=2)
  {
    GEN s = gel(a,i), t = gel(a,i+1);
    GEN x0 = gmul(s,t);
    GEN x1 = gneg(gadd(s,t));
    gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
  }
  if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
                                    scalarpol_shallow(gneg(gel(a,i)), v));
  for (i=r1+1; i<lx; i++)
  {
    GEN s = gel(a,i);
    GEN x0 = gnorm(s);
    GEN x1 = gneg(gtrace(s));
    gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
  }
  setlg(L, k); L = divide_conquer_prod(L, normalized_mul);
  return gerepileupto(av, normalized_to_RgX(L));
}

GEN
divide_conquer_assoc(GEN x, void *data, GEN (*mul)(void *,GEN,GEN))
{
  pari_sp ltop, lim;
  pari_long i,k,lx = lg(x);

  if (lx == 1) return gen_1;
  if (lx == 2) return gcopy(gel(x,1));
  x = leafcopy(x); k = lx;
  ltop=avma; lim = stack_lim(ltop,1);
  while (k > 2)
  {
    if (DEBUGLEVEL>7)
      err_printf("prod: remaining objects %ld\n",k-1);
    lx = k; k = 1;
    for (i=1; i<lx-1; i+=2)
      gel(x,k++) = mul(data,gel(x,i),gel(x,i+1));
    if (i < lx) gel(x,k++) = gel(x,i);
    if (low_stack(lim,stack_lim(ltop,1)))
      gerepilecoeffs(ltop,x+1,(int)k-1);
  }
  return gel(x,1);
}

static GEN
_domul(void *data, GEN x, GEN y)
{
  GEN (*mul)(GEN,GEN)=(GEN (*)(GEN,GEN)) data;
  return mul(x,y);
}
GEN
divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
{ return divide_conquer_assoc(x, (void *)mul, _domul); }

/*******************************************************************/
/*                                                                 */
/*                          FACTORBACK                             */
/*                                                                 */
/*******************************************************************/
static GEN
idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
static GEN
idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
static GEN
idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
static GEN
idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
static GEN
eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
static GEN
eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
static GEN
mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
static GEN
powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}

#if 0
static GEN
_ellmul(void *ell, GEN x, GEN y) { return elladd((GEN) ell, x, y); }
static GEN
_ellpow(void *ell, GEN x, GEN n) { return ellmul((GEN) ell, x, n); }
#endif

/* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
GEN
gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
               GEN (*_pow)(void*,GEN,GEN), void *data)
{
  pari_sp av = avma;
  pari_long k, l, lx;
  GEN p,x;

  if (e) /* supplied vector of exponents */
    p = L;
  else
  {
    switch(typ(L)) {
      case t_VEC:
      case t_COL: /* product of the L[i] */
        return gerepileupto(av, divide_conquer_assoc(L, data, _mul));
      case t_MAT: /* genuine factorization */
        l = lg(L);
        if (l == 1) return gen_1;
        if (l == 3) break;
        /*fall through*/
      default:
        pari_err_TYPE("factorback [not a factorization]", L);
    }
    p = gel(L,1);
    e = gel(L,2);
  }
  /* p = elts, e = expo */
  lx = lg(p);
  /* check whether e is an integral vector of correct length */
  if (!is_vec_t(typ(e)) || lx != lg(e) || !RgV_is_ZV(e))
    pari_err_TYPE("factorback [not an exponent vector]", e);
  if (lx == 1) return gen_1;
  x = cgetg(lx,t_VEC);
  for (l=1,k=1; k<lx; k++)
    if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
  x[0] = evaltyp(t_VEC) | _evallg(l);
  return gerepileupto(av, divide_conquer_assoc(x, data, _mul));
}

GEN
idealfactorback(GEN nf, GEN L, GEN e, int red)
{
  nf = checknf(nf);
  if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
  else     return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
}

GEN
nffactorback(GEN nf, GEN L, GEN e)
{ return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }

GEN
factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
GEN
factorback(GEN fa) { return factorback2(fa, NULL); }

static int
RgX_is_irred_i(GEN x)
{
  GEN y, p, pol;
  pari_long l = lg(x), pa;

  if (!signe(x) || l <= 3) return 0;
  switch(RgX_type(x,&p,&pol,&pa))
  {
    case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
    case t_COMPLEX: return l == 4;
    case t_REAL:
      if (l == 4) return 1;
      if (l > 5) return 0;
      return gsigne(RgX_disc(x)) > 0;
  }
  y = factor(x);
  return (lg(gcoeff(y,1,1))==l);
}
static int
RgX_is_irred(GEN x)
{
  pari_sp av = avma;
  int r = RgX_is_irred_i(x);
  avma = av; return r;
}
pari_long
isirreducible(GEN x)
{
  switch(typ(x))
  {
    case t_INT: case t_REAL: case t_FRAC: return 0;
    case t_POL: return RgX_is_irred(x);
  }
  pari_err_TYPE("isirreducible",x);
  return 0;
}

/*******************************************************************/
/*                                                                 */
/*                         GENERIC GCD                             */
/*                                                                 */
/*******************************************************************/
/* x is a COMPLEX or a QUAD */
static GEN
triv_cont_gcd(GEN x, GEN y)
{
  pari_sp av = avma;
  GEN c;
  if (typ(x)==t_COMPLEX)
  {
    GEN a = gel(x,1), b = gel(x,2);
    if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
    c = ggcd(a,b);
  }
  else
    c = ggcd(gel(x,2),gel(x,3));
  return gerepileupto(av, ggcd(c,y));
}

/* y is a PADIC, x a rational number or an INTMOD */
static GEN
padic_gcd(GEN x, GEN y)
{
  GEN p = gel(y,2);
  pari_long v = gvaluation(x,p), w = valp(y);
  if (w < v) v = w;
  return powis(p, v);
}

/* x,y in Z[i], at least one of which is t_COMPLEX */
static GEN
gauss_gcd(GEN x, GEN y)
{
  pari_sp av = avma;
  GEN dx, dy;
  dx = denom(x); x = gmul(x, dx);
  dy = denom(y); y = gmul(y, dy);
  while (!gequal0(y))
  {
    GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
    x = y; y = z;
  }
  x = gauss_normal(x);
  if (typ(x) == t_COMPLEX)
  {
    if      (gequal0(gel(x,2))) x = gel(x,1);
    else if (gequal0(gel(x,1))) x = gel(x,2);
  }
  return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
}

static int
is_rational(GEN x) { pari_long t = typ(x); return is_rational_t(t); }
static int
c_is_rational(GEN x)
{
  return (is_rational(gel(x,1)) && is_rational(gel(x,2)));
}
static GEN
c_zero_gcd(GEN c)
{
  GEN x = gel(c,1), y = gel(c,2);
  pari_long tx = typ(x), ty = typ(y);
  if (tx == t_REAL || ty == t_REAL) return gen_1;
  if (tx == t_PADIC || tx == t_INTMOD
   || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
  return gauss_gcd(c, gen_0);
}

/* y == 0 */
static GEN
zero_gcd(GEN x, pari_long tx)
{
  pari_sp av;
  switch(tx)
  {
    case t_INT: return absi(x);
    case t_FRAC: return absfrac(x);
    case t_COMPLEX: return c_zero_gcd(x);
    case t_REAL: return gen_1;
    case t_PADIC: return powis(gel(x,2), valp(x));
    case t_SER: return monomial(gen_1, valp(x), varn(x));
    case t_POLMOD: {
      GEN d = gel(x,2);
      if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
      return isinexact(d)? zero_gcd(d, typ(d)): gcopy(d);
    }
    case t_POL:
      if (!isinexact(x)) break;
      av = avma;
      return gerepileupto(av,
        monomialcopy(content(x), RgX_val(x), varn(x))
      );

    case t_RFRAC:
      if (!isinexact(x)) break;
      av = avma;
      return gerepileupto(av,
        gdiv(zero_gcd(gel(x,1), typ(gel(x,1))), gel(x,2))
      );
  }
  return gcopy(x);
}

/* tx = t_RFRAC, y considered as constant */
static GEN
cont_gcd_rfrac(GEN x, GEN y)
{
  pari_sp av = avma;
  GEN cx; x = primitive_part(x, &cx);
  return gerepileupto(av, gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2)));
}
/* tx = t_POL, y considered as constant */
static GEN
cont_gcd_pol(GEN x, GEN y)
{
  pari_sp av = avma;
  return gerepileupto(av, scalarpol(ggcd(content(x),y), varn(x)));
}
/* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
static GEN
cont_gcd_gen(GEN x, GEN y)
{
  pari_sp av = avma;
  return gerepileupto(av, ggcd(content(x),y));
}
/* !is_const(tx), y considered as constant */
static GEN
cont_gcd(GEN x, pari_long tx, GEN y)
{
  switch(tx)
  {
    case t_RFRAC: return cont_gcd_rfrac(x,y);
    case t_POL: return cont_gcd_pol(x,y);
    default: return cont_gcd_gen(x,y);
  }
}
static GEN
gcdiq(GEN x, GEN y)
{
  GEN z;
  if (!signe(x)) return Q_abs(y);
  z = cgetg(3,t_FRAC);
  gel(z,1) = gcdii(x,gel(y,1));
  gel(z,2) = icopy(gel(y,2));
  return z;
}
static GEN
gcdqq(GEN x, GEN y)
{
  GEN z = cgetg(3,t_FRAC);
  gel(z,1) = gcdii(gel(x,1), gel(y,1));
  gel(z,2) = lcmii(gel(x,2), gel(y,2));
  return z;
}
/* assume x,y t_INT or t_FRAC */
GEN
Q_gcd(GEN x, GEN y)
{
  pari_long tx = typ(x), ty = typ(y);
  if (tx == t_INT)
  { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
  else
  { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
}

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

  if (is_noncalc_t(tx) || is_noncalc_t(ty)) pari_err_TYPE2("gcd",x,y);
  if (is_matvec_t(ty))
  {
    z = cgetg_copy(y, &l);
    for (i=1; i<l; i++) gel(z,i) = ggcd(x,gel(y,i));
    return z;
  }
  if (is_matvec_t(tx))
  {
    z = cgetg_copy(x, &l);
    for (i=1; i<l; i++) gel(z,i) = ggcd(gel(x,i),y);
    return z;
  }
  if (tx>ty) { swap(x,y); lswap(tx,ty); }
  /* tx <= ty */
  if (isrationalzero(x)) return zero_gcd(y, ty);
  if (isrationalzero(y)) return zero_gcd(x, tx);
  if (is_const_t(tx))
  {
    if (ty == tx) switch(tx)
    {
      case t_INT:
        return gcdii(x,y);

      case t_INTMOD: z=cgetg(3,t_INTMOD);
        if (equalii(gel(x,1),gel(y,1)))
          gel(z,1) = icopy(gel(x,1));
        else
          gel(z,1) = gcdii(gel(x,1),gel(y,1));
        if (gequal1(gel(z,1))) gel(z,2) = gen_0;
        else
        {
          av = avma; p1 = gcdii(gel(z,1),gel(x,2));
          if (!is_pm1(p1))
          {
            p1 = gcdii(p1,gel(y,2));
            if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
            else p1 = gerepileuptoint(av, p1);
          }
          gel(z,2) = p1;
        }
        return z;

      case t_FRAC:
        return gcdqq(x,y);

      case t_FFELT:
        if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
        return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);

      case t_COMPLEX:
        if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
        return triv_cont_gcd(y,x);

      case t_PADIC:
        if (!equalii(gel(x,2),gel(y,2))) return gen_1;
        return powis(gel(y,2), minss(valp(x), valp(y)));

      case t_QUAD:
        av=avma; p1=gdiv(x,y);
        if (gequal0(gel(p1,3)))
        {
          p1=gel(p1,2);
          if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
          tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
        }
        if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
        p1 = ginv(p1); avma=av;
        if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
        return triv_cont_gcd(y,x);

      default: return gen_1; /* t_REAL */
    }
    if (is_const_t(ty)) switch(tx)
    {
      case t_INT:
        switch(ty)
        {
          case t_INTMOD: z = cgetg(3,t_INTMOD);
            gel(z,1) = icopy(gel(y,1)); av = avma;
            p1 = gcdii(gel(y,1),gel(y,2));
            if (!is_pm1(p1)) {
              p1 = gcdii(x,p1);
              if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
              else
                p1 = gerepileuptoint(av, p1);
            }
            gel(z,2) = p1; return z;

          case t_REAL: return gen_1;

          case t_FRAC:
            return gcdiq(x,y);

          case t_COMPLEX:
            if (c_is_rational(y)) return gauss_gcd(x,y);
            return triv_cont_gcd(y,x);

          case t_FFELT:
            if (!FF_equal0(y)) return FF_1(y);
            return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);

          case t_PADIC:
            return padic_gcd(x,y);

          case t_QUAD:
            return triv_cont_gcd(y,x);
          default:
            pari_err_TYPE2("gcd",x,y);
        }

      case t_REAL:
        switch(ty)
        {
          case t_INTMOD:
          case t_FFELT:
          case t_PADIC: pari_err_TYPE2("gcd",x,y);
          default: return gen_1;
        }

      case t_INTMOD:
        switch(ty)
        {
          case t_FRAC:
            av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
            if (!is_pm1(p1)) pari_err_OP("gcd",x,y);
            return ggcd(gel(y,1), x);

          case t_FFELT:
          {
            GEN p = gel(y,4);
            if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
            if (!FF_equal0(y)) return FF_1(y);
            return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
          }

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

          case t_PADIC:
            return padic_gcd(x,y);

          default: pari_err_TYPE2("gcd",x,y);
        }

      case t_FRAC:
        switch(ty)
        {
          case t_COMPLEX:
            if (c_is_rational(y)) return gauss_gcd(x,y);
          case t_QUAD:
            return triv_cont_gcd(y,x);
          case t_FFELT:
          {
            GEN p = gel(y,4);
            if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
            if (!FF_equal0(y)) return FF_1(y);
            return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
          }

          case t_PADIC:
            return padic_gcd(x,y);

          default: pari_err_TYPE2("gcd",x,y);
        }
      case t_FFELT:
        switch(ty)
        {
          case t_PADIC:
          {
            GEN p = gel(y,2);
            pari_long v = valp(y);
            if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
            return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
          }
          default: pari_err_TYPE2("gcd",x,y);
        }

      case t_COMPLEX:
        switch(ty)
        {
          case t_PADIC:
          case t_QUAD: return triv_cont_gcd(x,y);
          default: pari_err_TYPE2("gcd",x,y);
        }

      case t_PADIC:
        switch(ty)
        {
          case t_QUAD: return triv_cont_gcd(y,x);
          default: pari_err_TYPE2("gcd",x,y);
        }

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

  if (tx == t_POLMOD)
  {
    if (ty == t_POLMOD)
    {
      GEN T = gel(x,1);
      z = cgetg(3,t_POLMOD);
      T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
      gel(z,1) = T;
      if (degpol(T) <= 0) gel(z,2) = gen_0;
      else
      {
        GEN X, Y, d;
        av = avma; X = gel(x,2); Y = gel(y,2);
        d = ggcd(content(X), content(Y));
        if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
        p1 = ggcd(T, X);
        gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
      }
      return z;
    }
    vx = varn(gel(x,1));
    switch(ty)
    {
      case t_POL:
        vy = varn(y);
        if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
        z = cgetg(3,t_POLMOD);
        gel(z,1) = RgX_copy(gel(x,1));
        av = avma; p1 = ggcd(gel(x,1),gel(x,2));
        gel(z,2) = gerepileupto(av, ggcd(p1,y));
        return z;

      case t_RFRAC:
        vy = varn(gel(y,2));
        if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
        av = avma;
        p1 = ggcd(gel(x,1),gel(y,2));
        if (degpol(p1)) pari_err_OP("gcd",x,y);
        avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    }
  }

  vx = gvar(x);
  vy = gvar(y);
  if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
  if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);

  /* vx = vy: same main variable */
  switch(tx)
  {
    case t_POL:
      switch(ty)
      {
        case t_POL: return RgX_gcd(x,y);
        case t_SER:
          z = ggcd(content(x), content(y));
          return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
        case t_RFRAC: return cont_gcd_rfrac(y, x);
      }
      break;

    case t_SER:
      z = ggcd(content(x), content(y));
      switch(ty)
      {
        case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
        case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
      }
      break;

    case t_RFRAC:
    {
      GEN xd = gel(x,2), yd = gel(y,2);
      if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
      z = cgetg(3,t_RFRAC); av = avma;
      gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
      gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    }
  }
  pari_err_TYPE2("gcd",x,y);
  return NULL; /* not reached */
}
GEN
ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }

/* x a t_VEC,t_COL or t_MAT */
static GEN
vec_lcm(GEN x)
{
  if (typ(x) == t_MAT)
  {
    pari_long i, l = lg(x);
    GEN z = cgetg(l, t_VEC);
    for (i = 1; i < l; i++) gel(z,i) = glcm0(gel(x,i), NULL);
    x = z;
  }
  return glcm0(x, NULL);
}
static GEN
scal_lcm(GEN x, GEN y)
{
  pari_sp av = avma;
  pari_long tx = typ(x), ty = typ(y);
  if (is_matvec_t(tx)) x = vec_lcm(x);
  if (is_matvec_t(ty)) y = vec_lcm(y);
  return gerepileupto(av, glcm(x, y));
}

static GEN
fix_lcm(GEN x)
{
  GEN t;
  switch(typ(x))
  {
    case t_INT: if (signe(x)<0) x = negi(x);
      break;
    case t_POL:
      if (lg(x) <= 2) break;
      t = leading_term(x);
      if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
  }
  return x;
}

GEN
glcm0(GEN x, GEN y) {
  if (!y && lg(x) == 2)
  {
    pari_long tx = typ(x);
    if (is_vec_t(tx))
    {
      x = gel(x,1);
      tx = typ(x);
      return is_matvec_t(tx)? vec_lcm(x): fix_lcm(x);
    }
  }
  return gassoc_proto(scal_lcm,x,y);
}

GEN
glcm(GEN x, GEN y)
{
  pari_long tx, ty, i, l;
  pari_sp av;
  GEN p1, z;

  ty = typ(y);
  if (is_matvec_t(ty))
  {
    z = cgetg_copy(y, &l);
    for (i=1; i<l; i++) gel(z,i) = glcm(x,gel(y,i));
    return z;
  }
  tx = typ(x);
  if (is_matvec_t(tx))
  {
    z = cgetg_copy(x, &l);
    for (i=1; i<l; i++) gel(z,i) = glcm(gel(x,i),y);
    return z;
  }
  if (tx==t_INT && ty==t_INT) return lcmii(x,y);
  if (gequal0(x)) return gen_0;

  av = avma;
  p1 = ggcd(x,y); if (!gequal1(p1)) y = gdiv(y,p1);
  return gerepileupto(av, fix_lcm(gmul(x,y)));
}

/* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
static int
pol_approx0(GEN r, GEN x, int exact)
{
  pari_long i, lx,lr;
  if (exact) return gequal0(r);
  lx = lg(x);
  lr = lg(r); if (lr < lx) lx = lr;
  for (i=2; i<lx; i++)
    if (!approx_0(gel(r,i), gel(x,i))) return 0;
  return 1;
}

GEN
RgX_gcd_simple(GEN x, GEN y)
{
  pari_sp av1, av = avma, lim = stack_lim(av, 1);
  GEN r, yorig = y;
  int exact = !(isinexactreal(x) || isinexactreal(y));

  for(;;)
  {
    av1 = avma; r = RgX_rem(x,y);
    if (pol_approx0(r, x, exact))
    {
      avma = av1;
      if (y == yorig) return RgX_copy(y);
      y = normalizepol_approx(y, lg(y));
      if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
      return gerepileupto(av,y);
    }
    x = y; y = r;
    if (low_stack(lim,stack_lim(av,1))) {
      if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
      gerepileall(av,2, &x,&y);
    }
  }
}
GEN
RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
{
  pari_sp av = avma;
  GEN q, r, d, d1, u, v, v1;
  int exact = !(isinexactreal(a) || isinexactreal(b));

  d = a; d1 = b; v = gen_0; v1 = gen_1;
  for(;;)
  {
    if (pol_approx0(d1, a, exact)) break;
    q = poldivrem(d,d1, &r);
    v = gsub(v, gmul(q,v1));
    u=v; v=v1; v1=u;
    u=r; d=d1; d1=u;
  }
  u = gsub(d, gmul(b,v));
  u = RgX_div(u,a);

  gerepileall(av, 3, &u,&v,&d);
  *pu = u;
  *pv = v; return d;
}
/*******************************************************************/
/*                                                                 */
/*                    CONTENT / PRIMITIVE PART                     */
/*                                                                 */
/*******************************************************************/

GEN
content(GEN x)
{
  pari_long lx, i, t, tx = typ(x);
  pari_sp av = avma;
  GEN c;

  if (is_scalar_t(tx)) return zero_gcd(x, tx);
  switch(tx)
  {
    case t_RFRAC:
    {
      GEN n = gel(x,1), d = gel(x,2);
      /* -- varncmp(vn, vd) < 0 can't happen
       * -- if n is POLMOD, its main variable (in the sense of gvar2)
       *    has lower priority than denominator */
      if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
        n = isinexact(n)? zero_gcd(n, typ(n)): gcopy(n);
      else
        n = content(n);
      return gerepileupto(av, gdiv(n, content(d)));
    }

    case t_VEC: case t_COL:
      lx = lg(x); if (lx==1) return gen_1;
      break;

    case t_MAT:
    {
      pari_long hx, j;
      lx = lg(x);
      if (lx == 1) return gen_1;
      hx = lgcols(x);
      if (hx == 1) return gen_1;
      if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
      if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
      c = content(gel(x,1));
      for (j=2; j<lx; j++)
        for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
      if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
      return gerepileupto(av,c);
    }

    case t_POL: case t_SER:
      lx = lg(x); if (lx == 2) return gen_0;
      break;
    case t_QFR: case t_QFI:
      lx = 4; break;

    default: pari_err_TYPE("content",x);
      return NULL; /* not reached */
  }
  for (i=lontyp[tx]; i<lx; i++)
    if (typ(gel(x,i)) != t_INT) break;
  lx--; c = gel(x,lx);
  t = typ(c); if (is_matvec_t(t)) c = content(c);
  if (i > lx)
  { /* integer coeffs */
    while (lx-- > lontyp[tx])
    {
      c = gcdii(c, gel(x,lx));
      if (is_pm1(c)) { avma=av; return gen_1; }
    }
  }
  else
  {
    if (isinexact(c)) c = zero_gcd(c, typ(c));
    while (lx-- > lontyp[tx])
    {
      GEN d = gel(x,lx);
      t = typ(d); if (is_matvec_t(t)) d = content(d);
      c = ggcd(c, d);
    }
    if (isinexact(c)) { avma=av; return gen_1; }
  }
  switch(typ(c))
  {
    case t_INT:
      if (signe(c) < 0) c = negi(c);
      break;
    case t_VEC: case t_COL: case t_MAT:
      pari_err_TYPE("content",x);
  }

  return av==avma? gcopy(c): gerepileupto(av,c);
}

GEN
primitive_part(GEN x, GEN *ptc)
{
  pari_sp av = avma;
  GEN c = content(x);
  if (gequal1(c)) { avma = av; c = NULL; }
  else if (!gequal0(c)) x = gdiv(x,c);
  if (ptc) *ptc = c;
  return x;
}
GEN
primpart(GEN x) { return primitive_part(x, NULL); }

/* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
 * of Q(x2,...,xn)[x1] */
GEN
Q_content(GEN x)
{
  pari_long i, l;
  GEN d;
  pari_sp av;

  switch(typ(x))
  {
    case t_INT:  return absi(x);
    case t_FRAC: return absfrac(x);

    case t_VEC: case t_COL: case t_MAT:
      l = lg(x); if (l == 1) return gen_1;
      av = avma; d = Q_content(gel(x,1));
      for (i=2; i<l; i++)
      {
        d = Q_gcd(d, Q_content(gel(x,i)));
        if ((i & 255) == 0) d = gerepileupto(av, d);
      }
      return gerepileupto(av, d);

    case t_POL:
      l = lg(x); if (l == 2) return gen_0;
      av = avma; d = Q_content(gel(x,2));
      for (i=3; i<l; i++) d = Q_gcd(d, Q_content(gel(x,i)));
      return gerepileupto(av, d);
    case t_POLMOD: return Q_content(gel(x,2));
    case t_COMPLEX: return Q_gcd(Q_content(gel(x,1)), Q_content(gel(x,2)));
  }
  pari_err_TYPE("Q_content",x);
  return NULL; /* not reached */
}

GEN
ZX_content(GEN x)
{
  pari_long i, l = lg(x);
  GEN d;
  pari_sp av;

  if (l == 2) return gen_0;
  d = gel(x,2);
  if (l == 3) return absi(d);
  av = avma;
  for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
  if (signe(d) < 0) d = absi(d);
  return gerepileuptoint(av, d);
}

/* NOT MEMORY CLEAN (because of t_FRAC).
 * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
 * of Q(x2,...,xn)[x1] */
GEN
Q_denom(GEN x)
{
  pari_long i, l;
  GEN d, D;
  pari_sp av;

  switch(typ(x))
  {
    case t_INT: return gen_1;
    case t_FRAC: return gel(x,2);

    case t_VEC: case t_COL: case t_MAT:
      l = lg(x); if (l == 1) return gen_1;
      av = avma; d = Q_denom(gel(x,1));
      for (i=2; i<l; i++)
      {
        D = Q_denom(gel(x,i));
        if (D != gen_1) d = lcmii(d, D);
        if ((i & 255) == 0) d = gerepileuptoint(av, d);
      }
      return gerepileuptoint(av, d);

    case t_POL:
      l = lg(x); if (l == 2) return gen_1;
      av = avma; d = Q_denom(gel(x,2));
      for (i=3; i<l; i++)
      {
        D = Q_denom(gel(x,i));
        if (D != gen_1) d = lcmii(d, D);
      }
      return gerepileuptoint(av, d);
  }
  pari_err_TYPE("Q_denom",x);
  return NULL; /* not reached */
}

GEN
Q_remove_denom(GEN x, GEN *ptd)
{
  GEN d = Q_denom(x);
  if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d);
  if (ptd) *ptd = d;
  return x;
}

/* return y = x * d, assuming x rational, and d,y integral */
GEN
Q_muli_to_int(GEN x, GEN d)
{
  pari_long i, l;
  GEN y, xn, xd;
  pari_sp av;

  if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
  switch (typ(x))
  {
    case t_INT:
      return mulii(x,d);

    case t_FRAC:
      xn = gel(x,1);
      xd = gel(x,2); av = avma;
      y = mulii(xn, diviiexact(d, xd));
      return gerepileuptoint(av, y);

    case t_VEC: case t_COL: case t_MAT:
      y = cgetg_copy(x, &l);
      for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
      return y;

    case t_POL:
      y = cgetg_copy(x, &l); y[1] = x[1];
      for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
      return y;

    case t_POLMOD:
      y = cgetg(3, t_POLMOD);
      gel(y,1) = RgX_copy(gel(x,1));
      gel(y,2) = Q_muli_to_int(gel(x,2), d);
      return y;
  }
  pari_err_TYPE("Q_muli_to_int",x);
  return NULL; /* not reached */
}

/* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
static GEN
Q_divmuli_to_int(GEN x, GEN d, GEN n)
{
  pari_long i, l;
  GEN y, xn, xd;
  pari_sp av;

  switch(typ(x))
  {
    case t_INT:
      av = avma; y = diviiexact(x,d);
      return gerepileuptoint(av, mulii(y,n));

    case t_FRAC:
      xn = gel(x,1);
      xd = gel(x,2); av = avma;
      y = mulii(diviiexact(xn, d), diviiexact(n, xd));
      return gerepileuptoint(av, y);

    case t_VEC: case t_COL: case t_MAT:
      y = cgetg_copy(x, &l);
      for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
      return y;

    case t_POL:
      y = cgetg_copy(x, &l); y[1] = x[1];
      for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
      return y;

    case t_POLMOD:
      y = cgetg(3, t_POLMOD);
      gel(y,1) = RgX_copy(gel(x,1));
      gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
      return y;
  }
  pari_err_TYPE("Q_divmuli_to_int",x);
  return NULL; /* not reached */
}

/* return x / d. x: rational; d,result: integral. */
static GEN
Q_divi_to_int(GEN x, GEN d)
{
  pari_long i, l;
  GEN y;

  switch(typ(x))
  {
    case t_INT:
      return diviiexact(x,d);

    case t_VEC: case t_COL: case t_MAT:
      y = cgetg_copy(x, &l);
      for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
      return y;

    case t_POL:
      y = cgetg_copy(x, &l); y[1] = x[1];
      for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
      return y;

    case t_POLMOD:
      y = cgetg(3, t_POLMOD);
      gel(y,1) = RgX_copy(gel(x,1));
      gel(y,2) = Q_divi_to_int(gel(x,2), d);
      return y;
  }
  pari_err_TYPE("Q_divi_to_int",x);
  return NULL; /* not reached */
}
/* c t_FRAC */
static GEN
Q_divq_to_int(GEN x, GEN c)
{
  GEN n = gel(c,1), d = gel(c,2);
  if (is_pm1(n)) {
    GEN y = Q_muli_to_int(x,d);
    if (signe(n) < 0) y = gneg(y);
    return y;
  }
  return Q_divmuli_to_int(x, n,d);
}

/* return y = x / c, assuming x,c rational, and y integral */
GEN
Q_div_to_int(GEN x, GEN c)
{
  switch(typ(c))
  {
    case t_INT:  return Q_divi_to_int(x, c);
    case t_FRAC: return Q_divq_to_int(x, c);
  }
  pari_err_TYPE("Q_div_to_int",c);
  return NULL; /* not reached */
}
/* return y = x * c, assuming x,c rational, and y integral */
GEN
Q_mul_to_int(GEN x, GEN c)
{
  GEN d, n;
  switch(typ(c))
  {
    case t_INT: return Q_muli_to_int(x, c);
    case t_FRAC:
      n = gel(c,1);
      d = gel(c,2);
      return Q_divmuli_to_int(x, d,n);
  }
  pari_err_TYPE("Q_mul_to_int",c);
  return NULL; /* not reached */
}

GEN
Q_primitive_part(GEN x, GEN *ptc)
{
  pari_sp av = avma;
  GEN c = Q_content(x);
  if (typ(c) == t_INT)
  {
    if (is_pm1(c)) { avma = av; c = NULL; }
    else if (signe(c)) x = Q_divi_to_int(x, c);
  }
  else x = Q_divq_to_int(x, c);
  if (ptc) *ptc = c;
  return x;
}
GEN
Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }

/*******************************************************************/
/*                                                                 */
/*                           SUBRESULTANT                          */
/*                                                                 */
/*******************************************************************/
/* for internal use */
GEN
gdivexact(GEN x, GEN y)
{
  pari_long i,lx;
  GEN z;
  if (gequal1(y)) return x;
  switch(typ(x))
  {
    case t_INT:
      if (typ(y)==t_INT) return diviiexact(x,y);
      if (!signe(x)) return gen_0;
      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: { /* not stack-clean */
          pari_long v;
          if (varn(x)!=varn(y)) break;
          v = RgX_valrem(y,&y);
          if (v) x = RgX_shift_shallow(x,-v);
          if (!degpol(y)) { y = gel(y,2); break; }
          return RgX_div(x,y);
        }
      }
      return RgX_Rg_divexact(x, y);
    case t_VEC: case t_COL: case t_MAT:
      lx = lg(x); z = new_chunk(lx);
      for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
      z[0] = x[0]; return z;
  }
  if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
  return gdiv(x,y);
}

static GEN
init_resultant(GEN x, GEN y)
{
  pari_long tx = typ(x), ty = typ(y), vx, vy;
  if (is_scalar_t(tx) || is_scalar_t(ty))
  {
    if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    if (tx==t_POL) return gpowgs(y, degpol(x));
    if (ty==t_POL) return gpowgs(x, degpol(y));
    return gen_1;
  }
  if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
  if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
  if (!signe(x) || !signe(y)) return gmul(RgX_get_0(x),RgX_get_0(y)); /*type*/
  vx = varn(x);
  vy = varn(y); if (vx == vy) return NULL;
  return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
}

static pari_long
RgX_simpletype(GEN x)
{
  pari_long T = t_INT, i, lx = lg(x);
  for (i = 2; i < lx; i++)
  {
    GEN c = gel(x,i);
    pari_long tc = typ(c);
    switch(tc) {
      case t_INT:
        break;
      case t_FRAC:
        if (T == t_INT) T = t_FRAC;
        break;
      default:
        if (isinexact(c)) return t_REAL;
        T = 0; break;
    }
  }
  return T;
}

/* x an RgX, y a scalar */
static GEN
scalar_res(GEN x, GEN y, GEN *U, GEN *V)
{
  *V = gpowgs(y,degpol(x)-1);
  *U = gen_0; return gmul(y, *V);
}

/* return 0 if the subresultant chain can be interrupted.
 * Set u = NULL if the resultant is 0. */
static int
subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, pari_long *signh)
{
  GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
  pari_long du, dv, dr, degq;

  dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
  du = degpol(*u);
  dv = degpol(*v);
  degq = du - dv;
  if (*um1 == gen_1)
    u0 = gpowgs(gel(*v,dv+2),degq+1);
  else if (*um1 == gen_0)
    u0 = gen_0;
  else /* except in those 2 cases, um1 is an RgX */
    u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));

  if (*uze == gen_0) /* except in that case, uze is an RgX */
    u0 = scalarpol(u0, varn(*u)); /* now an RgX */
  else
    u0 = gsub(u0, gmul(q,*uze));

  *um1 = *uze;
  *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */

  *u = *v; c = *g; *g  = leading_term(*u);
  switch(degq)
  {
    case 0: break;
    case 1:
      c = gmul(*h,c); *h = *g; break;
    default:
      c = gmul(gpowgs(*h,degq), c);
      *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
  }
  *v  = RgX_Rg_divexact(r,c);
  *uze= RgX_Rg_divexact(*uze,c);
  if (both_odd(du, dv)) *signh = -*signh;
  return (dr > 3);
}

/* compute U, V s.t Ux + Vy = resultant(x,y) */
static GEN
subresext_i(GEN x, GEN y, GEN *U, GEN *V)
{
  pari_sp av, av2, lim;
  pari_long dx, dy, du, signh, tx = typ(x), ty = typ(y);
  GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;

  if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
  if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
  if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
  if (tx != t_POL) {
    if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    return scalar_res(y,x,V,U);
  }
  if (ty != t_POL) return scalar_res(x,y,U,V);
  if (varn(x) != varn(y))
    return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
                                        : scalar_res(y,x,V,U);
  dx = degpol(x); dy = degpol(y); signh = 1;
  if (dx < dy)
  {
    pswap(U,V); lswap(dx,dy); swap(x,y);
    if (both_odd(dx, dy)) signh = -signh;
  }
  if (dy == 0)
  {
    *V = gpowgs(gel(y,2),dx-1);
    *U = gen_0; return gmul(*V,gel(y,2));
  }
  av = avma;
  u = x = primitive_part(x, &cu);
  v = y = primitive_part(y, &cv);
  g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
  um1 = gen_1; uze = gen_0;
  for(;;)
  {
    if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    if (low_stack(lim,stack_lim(av2,1)))
    {
      if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
      gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    }
  }
  /* uze an RgX */
  if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
  z = gel(v,2); du = degpol(u);
  if (du > 1)
  { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    p1 = gpowgs(gdiv(z,h),du-1);
    z = gmul(z,p1);
    uze = RgX_Rg_mul(uze, p1);
  }
  if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }

  vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
  if (signe(r)) pari_warn(warner,"inexact computation in subresext");
  /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
  p1 = gen_1;
  if (cu) p1 = gmul(p1, gpowgs(cu,dy));
  if (cv) p1 = gmul(p1, gpowgs(cv,dx));
  cu = cu? gdiv(p1,cu): p1;
  cv = cv? gdiv(p1,cv): p1;
  z = gmul(z,p1);
  *U = RgX_Rg_mul(uze,cu);
  *V = RgX_Rg_mul(vze,cv);
  return z;
}
GEN
subresext(GEN x, GEN y, GEN *U, GEN *V)
{
  pari_sp av = avma;
  GEN z = subresext_i(x, y, U, V);
  gerepileall(av, 3, &z, U, V);
  return z;
}

static GEN
zero_extgcd(GEN y, GEN *U, GEN *V, pari_long vx)
{
  GEN x=content(y);
  *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
}

static int
must_negate(GEN x)
{
  GEN t = leading_term(x);
  switch(typ(t))
  {
    case t_INT: case t_REAL:
      return (signe(t) < 0);
    case t_FRAC:
      return (signe(gel(t,1)) < 0);
  }
  return 0;
}

/* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
GEN
RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
{
  pari_sp av, av2, tetpil, lim;
  pari_long signh; /* junk */
  pari_long dx, dy, vx, tx = typ(x), ty = typ(y);
  GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];

  if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
  if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
  if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
  vx=varn(x);
  if (!signe(x))
  {
    if (signe(y)) return zero_extgcd(y,U,V,vx);
    *U = pol_0(vx); *V = pol_0(vx);
    return pol_0(vx);
  }
  if (!signe(y)) return zero_extgcd(x,V,U,vx);
  dx = degpol(x); dy = degpol(y);
  if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
  if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }

  av = avma;
  u = x = primitive_part(x, &cu);
  v = y = primitive_part(y, &cv);
  g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
  um1 = gen_1; uze = gen_0;
  for(;;)
  {
    if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    if (low_stack(lim,stack_lim(av2,1)))
    {
      if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
      gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    }
  }
  if (uze != gen_0) {
    GEN r;
    vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    if (cu) uze = RgX_Rg_div(uze,cu);
    if (cv) vze = RgX_Rg_div(vze,cv);
    p1 = ginv(content(v));
  }
  else /* y | x */
  {
    vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    uze = pol_0(vx);
    p1 = gen_1;
  }
  if (must_negate(v)) p1 = gneg(p1);
  tetpil = avma;
  z = RgX_Rg_mul(v,p1);
  *U = RgX_Rg_mul(uze,p1);
  *V = RgX_Rg_mul(vze,p1);
  gptr[0] = &z;
  gptr[1] = U;
  gptr[2] = V;
  gerepilemanysp(av,tetpil,gptr,3); return z;
}

int
RgXQ_ratlift(GEN x, GEN T, pari_long amax, pari_long bmax, GEN *P, GEN *Q)
{
  pari_sp av = avma, av2, tetpil, lim;
  pari_long signh; /* junk */
  pari_long vx;
  GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];

  if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
  if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
  if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
  if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
  if (!signe(T)) {
    if (degpol(x) <= amax) {
      *P = RgX_copy(x);
      *Q = pol_1(varn(x));
      return 1;
    }
    return 0;
  }
  if (amax+bmax >= degpol(T))
    pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
                    mkvec3(stoi(amax), stoi(bmax), T));
  vx = varn(T);
  u = x = primitive_part(x, &cu);
  v = T = primitive_part(T, &cv);
  g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
  um1 = gen_1; uze = gen_0;
  for(;;)
  {
    (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
    if (typ(v)!=t_POL || degpol(v)<=amax) break;
    if (low_stack(lim,stack_lim(av2,1)))
    {
      if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
      gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    }
  }
  if (uze == gen_0)
  {
    avma = av; *P = pol_0(vx); *Q = pol_1(vx);
    return 1;
  }
  if (cu) uze = RgX_Rg_div(uze,cu);
  p1 = ginv(content(v));
  if (must_negate(v)) p1 = gneg(p1);
  tetpil = avma;
  *P = RgX_Rg_mul(v,p1);
  *Q = RgX_Rg_mul(uze,p1);
  gptr[0] = P;
  gptr[1] = Q;
  gerepilemanysp(av,tetpil,gptr,2); return 1;
}

/*******************************************************************/
/*                                                                 */
/*                RESULTANT USING DUCOS VARIANT                    */
/*                                                                 */
/*******************************************************************/
/* x^n / y^(n-1), assume n > 0 */
static GEN
Lazard(GEN x, GEN y, pari_long n)
{
  pari_long a;
  GEN c;

  if (n == 1) return x;
  a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
  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;
}

/* F (x/y)^(n-1), assume n >= 1 */
static GEN
Lazard2(GEN F, GEN x, GEN y, pari_long n)
{
  if (n == 1) return F;
  return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
}

static GEN
RgX_neg_i(GEN x, pari_long lx)
{
  pari_long i;
  GEN y = cgetg(lx, t_POL); y[1] = x[1];
  for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
  return y;
}
static GEN
RgX_Rg_mul_i(GEN y, GEN x, pari_long ly)
{
  pari_long i;
  GEN z;
  if (isrationalzero(x)) return pol_0(varn(y));
  z = cgetg(ly,t_POL); z[1] = y[1];
  for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
  return z;
}
static pari_long
reductum_lg(GEN x, pari_long lx)
{
  pari_long i = lx-2;
  while (i > 1 && gequal0(gel(x,i))) i--;
  return i+1;
}

/* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
 * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
static GEN
nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
{
  GEN p0, q0, h0, TMP, H, A, z0 = leading_term(Z);
  pari_long p, q, j, lP, lQ;
  pari_sp av, lim;

  p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
  q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
  /* p > q. Very often p - 1 = q */
  av = avma; lim = stack_lim(av,1);
  /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
  H = RgX_neg_i(Z, lQ); /* deg H < q */

  A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
  for (j = q+1; j < p; j++)
  {
    if (degpol(H) == q-1)
    { /* h0 = coeff of degree q-1 = leading coeff */
      h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
      H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    }
    else
      H = RgX_shift_shallow(H, 1);
    if (j+2 < lP)
    {
      TMP = RgX_Rg_mul(H, gel(P,j+2));
      A = A? RgX_add(A, TMP): TMP;
    }
    if (low_stack(lim,stack_lim(av,1)))
    {
      if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
      gerepileall(av,A?2:1,&H,&A);
    }
  }
  if (q+2 < lP) lP = reductum_lg(P, q+3);
  TMP = RgX_Rg_mul_i(P, z0, lP);
  A = A? RgX_add(A, TMP): TMP;
  A = RgX_Rg_divexact(A, p0);
  if (degpol(H) == q-1)
  {
    h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
  }
  else
    A = RgX_Rg_mul(addshift(H,A), q0);
  return RgX_Rg_divexact(A, s);
}

/* Ducos's subresultant */
GEN
RgX_resultant_all(GEN P, GEN Q, GEN *sol)
{
  pari_sp av, av2, lim;
  pari_long dP, dQ, delta, sig = 1;
  GEN cP, cQ, Z, s;

  dP = degpol(P);
  dQ = degpol(Q); delta = dP - dQ;
  if (delta < 0)
  {
    if (both_odd(dP, dQ)) sig = -1;
    swap(P,Q); lswap(dP, dQ); delta = -delta;
  }
  if (sol) *sol = gen_0;
  av = avma;
  if (dQ <= 0)
  {
    if (dQ < 0) return RgX_get_0(P);
    s = gpowgs(gel(Q,2), dP);
    if (sig == -1) s = gerepileupto(av, gneg(s));
    return s;
  }
  P = primitive_part(P, &cP);
  Q = primitive_part(Q, &cQ);
  av2 = avma; lim = stack_lim(av2,1);
  s = gpowgs(leading_term(Q),delta);
  if (both_odd(dP, dQ)) sig = -sig;
  Z = Q;
  Q = RgX_pseudorem(P, Q);
  P = Z;
  while(degpol(Q) > 0)
  {
    delta = degpol(P) - degpol(Q); /* > 0 */
    Z = Lazard2(Q, leading_term(Q), s, delta);
    if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    Q = nextSousResultant(P, Q, Z, s);
    P = Z;
    if (low_stack(lim,stack_lim(av,1)))
    {
      if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
      gerepileall(av2,2,&P,&Q);
    }
    s = leading_term(P);
  }
  if (!signe(Q)) { avma = av; return RgX_get_0(Q); }
  av2 = avma;
  s = Lazard(leading_term(Q), s, degpol(P));
  if (sig == -1) s = gneg(s);
  if (cP) s = gmul(s, gpowgs(cP,dQ));
  if (cQ) s = gmul(s, gpowgs(cQ,dP));
  if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
  return (avma == av2)? gerepilecopy(av, s): gerepileupto(av, s);
}
/* Return resultant(P,Q). If sol != NULL: set *sol to the last non-constant
 * polynomial in the prs IF the sequence was computed, and gen_0 otherwise.
 * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
 * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
 * in the "generic" case. */
GEN
resultant_all(GEN P, GEN Q, GEN *sol)
{
  pari_long TP, TQ;
  GEN s;

  if (sol) *sol = gen_0;
  if ((s = init_resultant(P,Q))) return s;
  if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
    return resultant2(P,Q); /* inexact */
  if (TP && TQ) /* rational */
  {
    if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
    return QX_resultant(P,Q);
  }
  return RgX_resultant_all(P, Q, sol);
}

/*******************************************************************/
/*                                                                 */
/*               RESULTANT USING SYLVESTER MATRIX                  */
/*                                                                 */
/*******************************************************************/
static GEN
_pol_0(void)
{
  GEN x = cgetg(3,t_POL);
  x[1] = 0;
  gel(x,2) = gen_0; return x;
}

static GEN
sylvester_col(GEN x, pari_long j, pari_long d, pari_long D)
{
  GEN c = cgetg(d+1,t_COL);
  pari_long i;
  for (i=1; i< j; i++) gel(c,i) = gen_0;
  for (   ; i<=D; i++) gel(c,i) = gel(x,D-i+2);
  for (   ; i<=d; i++) gel(c,i) = gen_0;
  return c;
}

GEN
sylvestermatrix_i(GEN x, GEN y)
{
  pari_long j,d,dx,dy;
  GEN M;

  dx = degpol(x); if (dx < 0) { dx = 0; x = _pol_0(); }
  dy = degpol(y); if (dy < 0) { dy = 0; y = _pol_0(); }
  d = dx+dy; M = cgetg(d+1,t_MAT);
  for (j=1; j<=dy; j++) gel(M,j)    = sylvester_col(x,j,d,j+dx);
  for (j=1; j<=dx; j++) gel(M,j+dy) = sylvester_col(y,j,d,j+dy);
  return M;
}

GEN
sylvestermatrix(GEN x, GEN y)
{
  pari_long i,j,lx;
  if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
  if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
  if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
  x = sylvestermatrix_i(x,y); lx = lg(x);
  for (i=1; i<lx; i++)
    for (j=1; j<lx; j++) gcoeff(x,i,j) = gcopy(gcoeff(x,i,j));
  return x;
}

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

/* Let vx = main variable of x. Return a polynomial in variable 0:
 * if vx is 0 and v != 0, set *mx = 1 and replace vx by pol_x(MAXVARN)
 * if vx = v, copy x, set its main variable to 0 and return
 * if vx < v, return subst(x, v, pol_x(0))
 * if vx > v, return scalarpol(x, 0) */
static GEN
fix_pol(GEN x, pari_long v, pari_long *mx)
{
  pari_long vx;
  if (typ(x) != t_POL) return x;
  vx = varn(x);
  if (v == vx)
  {
    if (v) { x = leafcopy(x); setvarn(x, 0); }
    return x;
  }
  if (!vx)
  {
    *mx = 1;
    x = poleval(x, pol_x(MAXVARN));
    vx = varn(x);
    if (v == vx) { setvarn(x, 0); return x; }
  }
  if (varncmp(v, vx) > 0)
  {
    x = gsubst(x,v,pol_x(0));
    if (typ(x) == t_POL && varn(x) == 0) return x;
  }
  return scalarpol_shallow(x, 0);
}

/* 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, pari_long v, pari_long flag)
{
  pari_long m = 0;
  pari_sp av = avma;

  if (v >= 0)
  {
    x = fix_pol(x,v, &m);
    y = fix_pol(y,v, &m);
  }
  switch(flag)
  {
    case 2:
    case 0: x=resultant_all(x,y,NULL); break;
    case 1: x=resultant2(x,y); break;
    default: pari_err_FLAG("polresultant");
  }
  if (m) x = gsubst(x,MAXVARN,pol_x(0));
  return gerepileupto(av,x);
}
GEN
polresultantext0(GEN x, GEN y, pari_long v)
{
  GEN R, U, V;
  pari_long m = 0;
  pari_sp av = avma;

  if (v >= 0)
  {
    x = fix_pol(x,v, &m);
    y = fix_pol(y,v, &m);
  }
  R = subresext_i(x,y, &U,&V);
  if (m)
  {
    U = gsubst(gsubst(U, 0, pol_x(v)), MAXVARN, pol_x(0));
    V = gsubst(gsubst(V, 0, pol_x(v)), MAXVARN, pol_x(0));
    R = gsubst(R,MAXVARN,pol_x(0));
  }
  else if (v >= 0)
  {
    if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
  }
  return gerepilecopy(av, mkvec3(U,V,R));
}
GEN
polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }

/*******************************************************************/
/*                                                                 */
/*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
/*                                                                 */
/*******************************************************************/

/* (v - x)^d */
static GEN
caract_const(pari_sp av, GEN x, pari_long v, pari_long d)
{ return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }

/* return caract(Mod(x,T)) in variable v */
GEN
RgXQ_charpoly(GEN x, GEN T, pari_long v)
{
  pari_sp av = avma;
  pari_long d = degpol(T), dx, vx, vp;
  GEN ch, L;

  if (typ(x) != t_POL) return caract_const(av, x, v, d);
  vx = varn(x);
  vp = varn(T);
  if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
  if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
  dx = degpol(x);
  if (dx <= 0)
    return dx? monomial(gen_1, d, v): caract_const(av, gel(x,2), v, d);

  x = RgX_neg(x);
  if (varn(x) == MAXVARN) { setvarn(x, 0); T = leafcopy(T); setvarn(T, 0); }
  gel(x,2) = gadd(gel(x,2), pol_x(MAXVARN));
  ch = resultant_all(T, x, NULL);
  if (v != MAXVARN)
  {
    if (typ(ch) == t_POL && varn(ch) == MAXVARN)
      setvarn(ch, v);
    else
      ch = gsubst(ch, MAXVARN, pol_x(v));
  }
  /* test for silly input: x mod (deg 0 polynomial) */
  if (typ(ch) != t_POL) { avma = av; return pol_1(v); }

  L = leading_term(ch);
  if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
  return gerepileupto(av, ch);
}

/* characteristic polynomial (in v) of x over nf, where x is an element of the
 * algebra nf[t]/(Q(t)) */
GEN
rnfcharpoly(GEN nf, GEN Q, GEN x, pari_long v)
{
  const char *f = "rnfcharpoly";
  pari_long dQ = degpol(Q);
  pari_sp av = avma;
  GEN T;

  if (v < 0) v = 0;
  nf = checknf(nf); T = nf_get_pol(nf);
  Q = RgX_nffix(f, T,Q,0);
  switch(typ(x))
  {
    case t_INT:
    case t_FRAC: return caract_const(av, x, v, dQ);
    case t_POLMOD:
      x = polmod_nffix2(f,T,Q, x,0);
      break;
    case t_POL:
      x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
      break;
    default: pari_err_TYPE(f,x);
  }
  if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
  /* x a t_POL in variable vQ */
  if (degpol(x) >= dQ) x = RgX_rem(x, Q);
  if (dQ <= 1) return caract_const(av, constant_term(x), v, 1);
  return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
}

/*******************************************************************/
/*                                                                 */
/*                  GCD USING SUBRESULTANT                         */
/*                                                                 */
/*******************************************************************/
static int inexact(GEN x, int *simple, int *rational);
static int
isinexactall(GEN x, int *simple, int *rational)
{
  pari_long i, lx = lg(x);
  for (i=2; i<lx; i++)
    if (inexact(gel(x,i), simple, rational)) return 1;
  return 0;
}
/* return 1 if coeff explosion is not possible */
static int
inexact(GEN x, int *simple, int *rational)
{
  int junk = 0;
  switch(typ(x))
  {
    case t_INT: case t_FRAC: return 0;

    case t_REAL: case t_PADIC: case t_SER: return 1;

    case t_INTMOD:
    case t_FFELT:
      *rational = 0;
      if (!*simple) *simple = 1;
      return 0;

    case t_COMPLEX:
      *rational = 0;
      return inexact(gel(x,1), simple, rational)
          || inexact(gel(x,2), simple, rational);
    case t_QUAD:
      *rational = *simple = 0;
      return inexact(gel(x,2), &junk, rational)
          || inexact(gel(x,3), &junk, rational);

    case t_POLMOD:
      *rational = 0;
      return isinexactall(gel(x,1), simple, rational);
    case t_POL:
      *rational = 0;
      *simple = -1;
      return isinexactall(x, &junk, rational);
    case t_RFRAC:
      *rational = 0;
      *simple = -1;
      return inexact(gel(x,1), &junk, rational)
          || inexact(gel(x,2), &junk, rational);
  }
  *rational = 0;
  *simple = -1; return 0;
}

/* x monomial, y t_POL in the same variable */
static GEN
gcdmonome(GEN x, GEN y)
{
  pari_sp av = avma;
  pari_long dx = degpol(x), e = RgX_valrem(y, &y);
  pari_long i, l = lg(y);
  GEN t, v = cgetg(l, t_VEC);
  gel(v,1) = gel(x,dx+2);
  for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
  t = content(v); /* gcd(lc(x), cont(y)) */
  t = simplify_shallow(t);
  if (dx < e) e = dx;
  return gerepileupto(av, monomialcopy(t, e, varn(x)));
}

/* x, y are t_POL in the same variable */
GEN
RgX_gcd(GEN x, GEN y)
{
  pari_long dx, dy;
  pari_sp av, av1, lim;
  GEN d, g, h, p1, p2, u, v;
  int simple = 0, rational = 1;

  if (isexactzero(y)) return RgX_copy(x);
  if (isexactzero(x)) return RgX_copy(y);
  if (RgX_is_monomial(x)) return gcdmonome(x,y);
  if (RgX_is_monomial(y)) return gcdmonome(y,x);
  if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
  {
    av = avma; u = ggcd(content(x), content(y));
    return gerepileupto(av, scalarpol(u, varn(x)));
  }
  if (rational) return QX_gcd(x,y); /* Q[X] */

  av = avma;
  if (simple > 0) x = RgX_gcd_simple(x,y);
  else
  {
    dx = lg(x); dy = lg(y);
    if (dx < dy) { swap(x,y); lswap(dx,dy); }
    if (dy==3)
    {
      d = ggcd(gel(y,2), content(x));
      return gerepileupto(av, scalarpol(d, varn(x)));
    }
    u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    d = ggcd(p1,p2);
    av1 = avma; lim = stack_lim(av1,1);
    g = h = gen_1;
    for(;;)
    {
      GEN r = RgX_pseudorem(u,v);
      pari_long degq, du, dv, dr = lg(r);

      if (!signe(r)) break;
      if (dr <= 3)
      {
        avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
      }
      if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
      du = lg(u); dv = lg(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(gpowgs(h,degq), p1);
          h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
      }
      v = RgX_Rg_div(r,p1);
      if (low_stack(lim, stack_lim(av1,1)))
      {
        if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
        gerepileall(av1,4, &u,&v,&g,&h);
      }
    }
    x = RgX_Rg_mul(primpart(v), d);
  }
  if (must_negate(x)) x = RgX_neg(x);
  return gerepileupto(av,x);
}

static GEN
RgX_disc_aux(GEN x)
{
  pari_long dx = degpol(x), Tx;
  GEN D, L, y, p;
  if (!signe(x) || !dx) return RgX_get_0(x);
  if (dx == 1) return RgX_get_1(x);
  if (dx == 2) {
    GEN a = gel(x,4), b = gel(x,3), c = gel(x,2);
    return gsub(gsqr(b), gmul2n(gmul(a,c),2));
  }
  Tx = RgX_simpletype(x);
  if (Tx == t_INT) return ZX_disc(x);
  if (Tx == t_FRAC) return QX_disc(x);
  p = NULL;
  if (RgX_is_FpX(x, &p) && p)
    return Fp_to_mod(FpX_disc(RgX_to_FpX(x,p), p), p);

  y = RgX_deriv(x);
  if (!signe(y)) return RgX_get_0(y);
  if (Tx == t_REAL)
    D = resultant2(x,y);
  else
  {
    D = RgX_resultant_all(x, y, NULL);
    if (D == gen_0) return RgX_get_0(y);
  }
  L = leading_term(x); if (!gequal1(L)) D = gdiv(D,L);
  if (dx & 2) D = gneg(D);
  return D;
}
GEN
RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }

GEN
poldisc0(GEN x, pari_long v)
{
  pari_long i;
  pari_sp av;
  GEN z, D;

  switch(typ(x))
  {
    case t_POL:
      av = avma; i = 0;
      if (v >= 0 && v != varn(x)) x = fix_pol(x,v, &i);
      D = RgX_disc_aux(x);
      if (i) D = gsubst(D, MAXVARN, pol_x(0));
      return gerepileupto(av, D);

    case t_COMPLEX:
      return utoineg(4);

    case t_QUAD:
      return quad_disc(x);
    case t_POLMOD:
      return poldisc0(gel(x,1), v);

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

    case t_VEC: case t_COL: case t_MAT:
      z = cgetg_copy(x, &i);
      for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
      return z;
  }
  pari_err_TYPE("poldisc",x);
  return NULL; /* not reached */
}

GEN
reduceddiscsmith(GEN x)
{
  pari_long j, n = degpol(x);
  pari_sp av = avma;
  GEN xp, M;

  if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
  if (n<=0) pari_err_CONSTPOL("poldiscreduced");
  RgX_check_ZX(x,"poldiscreduced");
  if (!gequal1(gel(x,n+2)))
    pari_err_IMPL("non-monic polynomial in poldiscreduced");
  M = cgetg(n+1,t_MAT);
  xp = ZX_deriv(x);
  for (j=1; j<=n; j++)
  {
    gel(M,j) = RgX_to_RgV(xp, n);
    if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
  }
  return gerepileupto(av, ZM_snf(M));
}

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

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

  if (gequal0(x)) pari_err_ROOTS0("sturm");
  t = typ(x);
  if (t != t_POL)
  {
    if (t == t_INT || t == t_REAL || t == t_FRAC) return 0;
    pari_err_TYPE("sturm",x);
  }
  s=lg(x); if (s==3) return 0;

  sl = gsigne(leading_term(x));
  if (s==4)
  {
    t = a? gsigne(poleval(x,a)): -sl;
    if (t == 0) { avma = av; return 0; }
    s = b? gsigne(poleval(x,b)):  sl;
    avma = av; return (s == t)? 0: 1;
  }
  u = primpart(x);
  v = primpart(RgX_deriv(x));
  g=gen_1; h=gen_1;
  s = b? gsigne(poleval(u,b)): sl;
  t = a? gsigne(poleval(u,a)): ((lg(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 = RgX_pseudorem(u,v);
    pari_long du=lg(u), dv=lg(v), dr=lg(r), degq=du-dv;

    if (dr<=2) pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,x);
    if (gsigne(leading_term(v)) > 0 || degq&1) r=gneg_i(r);
    sl = gsigne(gel(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(gpowgs(h,degq),p1);
        h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
    }
    v = RgX_Rg_divexact(r,p1);
    if (low_stack(lim,stack_lim(av,1)))
    {
      if(DEBUGMEM>1) pari_warn(warnmem,"polsturm, dr = %ld",dr);
      gerepileall(av,4,&u,&v,&g,&h);
    }
  }
}

/***********************************************************************/
/**                                                                   **/
/**                        GENERIC EXTENDED GCD                       **/
/**                                                                   **/
/***********************************************************************/
/* assume typ(x) = typ(y) = t_POL */
GEN
RgXQ_inv(GEN x, GEN y)
{
  pari_long vx=varn(x), vy=varn(y);
  pari_sp av;
  GEN u, v, d;

  while (vx != vy)
  {
    if (varncmp(vx,vy) > 0)
    {
      d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
      return scalarpol(d, vy);
    }
    if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    x = gel(x,2); vx = gvar(x);
  }
  av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
  if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
  d = gdiv(u,d);
  if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
  return gerepileupto(av, d);
}

/*Assume x is a polynomial and y is not */
static GEN
scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
{
  pari_long vx = varn(x);
  int xis0 = signe(x)==0, yis0 = gequal0(y);
  if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
  if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
  *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
}
/* Assume x==0, y!=0 */
static GEN
zero_bezout(GEN y, GEN *U, GEN *V)
{
  *U=gen_0; *V = ginv(y); return gen_1;
}

GEN
gbezout(GEN x, GEN y, GEN *u, GEN *v)
{
  pari_long tx=typ(x), ty=typ(y), vx;
  if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
  if (tx != t_POL)
  {
    if (ty == t_POL)
      return scalar_bezout(y,x,v,u);
    else
    {
      int xis0 = gequal0(x), yis0 = gequal0(y);
      if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
      if (xis0) return zero_bezout(y,u,v);
      else return zero_bezout(x,v,u);
    }
  }
  else if (ty != t_POL) return scalar_bezout(x,y,u,v);
  vx = varn(x);
  if (vx != varn(y))
    return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
                                   : scalar_bezout(y,x,v,u);
  return RgX_extgcd(x,y,u,v);
}

GEN
gcdext0(GEN x, GEN y)
{
  GEN z=cgetg(4,t_VEC);
  gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
  return z;
}

/*******************************************************************/
/*                                                                 */
/*                    GENERIC (modular) INVERSE                    */
/*                                                                 */
/*******************************************************************/

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

  switch(typ(y))
  {
    case t_POL:
      if (tx==t_POL) return RgXQ_inv(x,y);
      if (is_scalar_t(tx)) return ginv(x);
      break;
    case t_INT:
      if (tx==t_INT) return Fp_inv(x,y);
      if (tx==t_POL) return gen_0;
  }
  pari_err_TYPE2("ginvmod",x,y);
  return NULL; /* not reached */
}

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

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

  if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
  n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
  y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
  vval = (pari_long *) pari_malloc(sizeof(pari_long)*(n+1));
  for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
  for (a=0, ind=1; a<n; a++)
  {
    if (vval[a] != LONG_MAX) break;
    gel(y,ind++) = utoipos(LONG_MAX);
  }
  for (b=a+1; b<=n; a=b, b=a+1)
  {
    while (vval[b] == LONG_MAX) b++;
    u1=vval[a]-vval[b]; u2=b-a;
    for (c=b+1; c<=n; c++)
    {
      if (vval[c] == LONG_MAX) 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); gel(y,ind++) = gdivgs(num,u2); }
  }
  pari_free(vval); return y;
}