/********************************************************************/
/**                                                                **/
/**                      GENERIC OPERATIONS                        **/
/**                        (second part)                           **/
/**                                                                **/
/********************************************************************/
/* $Id: gen2.c,v 1.3 1999/09/23 17:50:56 karim Exp $ */
#include "pari.h"

/*******************************************************************/
/*                                                                 */
/*                    OPERATIONS BY VALUE                          */
/*            f is a pointer to the function called.               */
/*            result is gaffect-ed to last parameter               */
/*                                                                 */
/*******************************************************************/

void
gop1z(GEN (*f)(GEN), GEN x, GEN y)
{
  long av=avma; gaffect(f(x),y); avma=av;
}

void
gop2z(GEN (*f)(GEN, GEN), GEN x, GEN y, GEN z)
{
  long av=avma; gaffect(f(x,y),z); avma=av;
}

void
gops2gsz(GEN (*f)(GEN, long), GEN x, long s, GEN z)
{
  long av=avma; gaffect(f(x,s),z); avma=av;
}

void
gops2sgz(GEN (*f)(long, GEN), long s, GEN y, GEN z)
{
  long av=avma; gaffect(f(s,y),z); avma=av;
}

void
gops2ssz(GEN (*f)(long, long), long s, long y, GEN z)
{
  long av=avma; gaffect(f(s,y),z); avma=av;
}

/*******************************************************************/
/*                                                                 */
/*                OPERATIONS USING SMALL INTEGERS                  */
/*                                                                 */
/*******************************************************************/

/* small int prototype */
static long court_p[] = { evaltyp(t_INT) | m_evallg(3),0,0 };

GEN
gopsg2(GEN (*f)(GEN, GEN), long s, GEN y)
{
  affsi(s,court_p); return f(court_p,y);
}

GEN
gopgs2(GEN (*f)(GEN, GEN), GEN y, long s)
{
  affsi(s,court_p); return f(y,court_p);
}

long
opgs2(int (*f)(GEN, GEN), GEN y, long s)
{
  affsi(s,court_p); return f(y,court_p);
}

void
gopsg2z(GEN (*f)(GEN, GEN), long s, GEN y, GEN z)
{
  long av=avma;
  affsi(s,court_p); gaffect(f(court_p,y),z); avma=av;
}

void
gopgs2z(GEN (*f)(GEN, GEN), GEN y, long s, GEN z)
{
  long av=avma;
  affsi(s,court_p); gaffect(f(y,court_p),z); avma=av;
}

/*******************************************************************/
/*                                                                 */
/*                    CREATION OF A P-ADIC GEN                     */
/*                                                                 */
/*******************************************************************/

GEN
cgetp(GEN x)
{
  GEN y = cgetg(5,t_PADIC);
  y[1] = evalprecp(precp(x)) | evalvalp(0);
  icopyifstack(x[2], y[2]);
  y[3] = licopy((GEN)x[3]);
  y[4] = lgeti(lgefint(x[3])); return y;
}

/* y[4] not filled */
GEN
cgetp2(GEN x, long v)
{
  GEN y = cgetg(5,t_PADIC);
  y[1] = evalprecp(precp(x)) | evalvalp(v);
  icopyifstack(x[2], y[2]);
  y[3] = licopy((GEN)x[3]); return y;
}

/*******************************************************************/
/*                                                                 */
/*                       CLONING & COPY                            */
/*                  Replicate an existing GEN                      */
/*                                                                 */
/*******************************************************************/
/* lontyp = 0 means non recursive type
 * otherwise:
 *   lontyp = number of codewords
 *   if not in stack, we don't copy the words in [lontyp,lontyp2[
 */
const  long  lontyp[] = { 0,0,0,1,1,1,1,2,1,1, 2,2,0,1,1,1,1,1,1,1, 2,0,0 };
static long lontyp2[] = { 0,0,0,2,1,1,1,3,2,2, 2,2,0,1,1,1,1,1,1,1, 2,0,0 };

/* can't do a memcpy there: avma and x may overlap. memmove is slower */
GEN
gcopy(GEN x)
{
  long tx=typ(x),lx,i;
  GEN y;

  if (tx == t_SMALL) return x;
  lx = lg(x); y=new_chunk(lx);
  if (! is_recursive_t(tx))
    for (i=lx-1; i>=0; i--) y[i]=x[i];
  else
  {
    if (tx==t_POL || tx==t_LIST) lx=lgef(x);
    for (i=0; i<lontyp[tx];  i++) y[i]=x[i];
    for (   ; i<lontyp2[tx]; i++) copyifstack(x[i],y[i]);
    for (   ; i<lx;          i++) y[i]=lcopy((GEN)x[i]);
  }
  /* unsetisclone(y); useless because gunclone checks isonstack */
  return y;
}

GEN
gcopy_i(GEN x, long lx)
{
  long tx=typ(x),i;
  GEN y;

  if (tx == t_SMALL) return x;
  y=cgetg(lx,tx);
  if (! is_recursive_t(tx))
    for (i=lx-1; i>0; i--) y[i]=x[i];
  else
  {
    for (i=1; i<lontyp[tx];  i++) y[i]=x[i];
    for (   ; i<lontyp2[tx]; i++) copyifstack(x[i],y[i]);
    for (   ; i<lx;          i++) y[i]=lcopy((GEN)x[i]);
  }
  return y;
}

GEN
forcecopy(GEN x)
{
  long tx=typ(x),lx,i;
  GEN y;

  if (tx == t_SMALL) return x;
  lx=lg(x); y=new_chunk(lx);
  if (! is_recursive_t(tx))
    for (i=lx-1; i>=0; i--) y[i]=x[i];
  else
  {
    if (tx==t_POL || tx==t_LIST) lx=lgef(x);
    for (i=0; i<lontyp[tx]; i++) y[i]=x[i];
    for (   ; i<lx;         i++) y[i]=(long)forcecopy((GEN)x[i]);
  }
  unsetisclone(y); return y;
}

GEN
dummycopy(GEN x)
{
  long tx=typ(x), lx=lg(x),i;
  GEN y=new_chunk(lx);

  switch(tx)
  {
    case t_POLMOD:
      y[1]=x[1]; y[2]=(long) dummycopy((GEN)x[2]);
      break;
    case t_MAT:
      for (i=lx-1;i;i--) y[i]=(long)dummycopy((GEN)x[i]);
      break;
    default:
      for (i=lx-1;i;i--) y[i]=x[i];
  }
  y[0]=x[0]; return y;
}

/* assume x is a t_POL. FOR INTERNAL USE! */
GEN
dummyclone(GEN x)
{
  long lx = lgef(x);
  GEN z = (GEN)gpmalloc(lx*sizeof(long));
  while (--lx >= 0) z[lx] = x[lx];
  return z;
}

long
taille(GEN x)
{
  long i,n,lx=lg(x),tx=typ(x);

  n = lx;
  if (is_recursive_t(tx))
  {
    if (tx==t_POL || tx==t_LIST) lx = lgef(x);
    for (i=lontyp[tx]; i<lx; i++) n += taille((GEN)x[i]);
  }
  return n;
}

long
glength(GEN x)
{
  switch(typ(x))
  {
    case t_INT: return lgefint(x)-2;
    case t_POL: case t_LIST: return lgef(x)-2;
    case t_REAL: return signe(x)? lg(x)-2: 0;
    case t_STR: return strlen(GSTR(x));
  }
  return lg(x)-lontyp[typ(x)];
}

GEN
matsize(GEN x)
{
  GEN y=cgetg(3,t_VEC);

  switch(typ(x))
  {
    case t_VEC:
      y[1]=un; y[2]=lstoi(lg(x)-1); break;
    case t_COL:
      y[1]=lstoi(lg(x)-1); y[2]=un; break;
    case t_MAT:
      y[2]=lstoi(lg(x)-1);
      y[1]=(lg(x)==1)? zero: lstoi(lg(x[1])-1); break;
    default: err(typeer,"matsize");
  }
  return y;
}

long
taille2(GEN x) { return taille(x)<<TWOPOTBYTES_IN_LONG; }

GEN
brutcopy(GEN x, GEN y)
{
  long i,lx,tx=typ(x);
  GEN z;

  if (! is_recursive_t(tx))
  {
    lx = (tx==t_INT)? lgefint(x): lg(x);
    for (i=0; i<lx; i++) y[i] = x[i];
  }
  else
  {
    lx=lg(x); z = y + lx;
    if (tx==t_POL || tx==t_LIST) lx = lgef(x);
    for (i=0; i<lontyp[tx]; i++) y[i] = x[i];
    for (   ; i<lx; i++)
    {
      y[i] = (long)brutcopy((GEN)x[i], z);
      z += taille((GEN)x[i]);
    }
  }
  unsetisclone(y); return y;
}

GEN
gclone(GEN x)
{
  x = brutcopy(x, newbloc(taille(x)));
  setisclone(x); return x;
}

/*******************************************************************/
/*                                                                 */
/*                            GREFFE                               */
/*              Greffe d'une serie sur un polynome                 */
/*                                                                 */
/*******************************************************************/

GEN
greffe(GEN x, long l, long use_stack)
{
  long i,e,k,vx;
  GEN y;

  if (typ(x)!=t_POL) err(notpoler,"greffe");
  if (use_stack) y = cgetg(l,t_SER);
  else
  {
    y = (GEN) gpmalloc(l*sizeof(long));
    y[0] = evaltyp(t_SER) | evallg(l);
  }
  if (gcmp0(x))
  {
    y[1]=evalvalp(l-2) | evalvarn(varn(x));
    i=2; while(i<l) { y[i]=x[2]; i++; }
    return y;
  }
  vx=varn(x); e=gval(x,vx);
  y[1]=evalsigne(1) | evalvalp(e) | evalvarn(vx);
  k=lgef(x)-e-1; i=l-1;
  if (k<l)
    while (i>k) { y[i]=zero; i--; }
  while (i>=2) { y[i]=x[i+e]; i--; }
  return y;
}

/*******************************************************************/
/*                                                                 */
/*                 CONVERSION GEN --> long                         */
/*                                                                 */
/*******************************************************************/

long
gtolong(GEN x)
{
  long y,tx=typ(x),av=avma;

  switch(tx)
  {
    case t_INT:
      return itos(x);
    case t_REAL: case t_FRAC: case t_FRACN:
      y=itos(ground(x)); avma=av; return y;
    case t_COMPLEX:
      if (gcmp0((GEN)x[2])) return gtolong((GEN)x[1]); break;
    case t_QUAD:
      if (gcmp0((GEN)x[3])) return gtolong((GEN)x[2]); break;
  }
  err(typeer,"gtolong");
  return 0; /* not reached */
}

/*******************************************************************/
/*                                                                 */
/*                      COMPARE TO ZERO                            */
/*        returns 1 whenever the GEN x is 0, and 0 otherwise       */
/*                                                                 */
/*******************************************************************/

int
gcmp0(GEN x)
{
  switch(typ(x))
  {
    case t_INT: case t_REAL: case t_POL: case t_SER:
      return !signe(x);

    case t_INTMOD: case t_POLMOD:
      return gcmp0((GEN)x[2]);

    case t_FRAC: case t_FRACN:
      return !signe(x[1]);

    case t_COMPLEX:
     /* is 0 iff norm(x) would be 0 (can happen with Re(x) and Im(x) != 0
      * only if Re(x) and Im(x) are of type t_REAL). See mp.c:addrr().
      */
      if (gcmp0((GEN)x[1]))
      {
	if (gcmp0((GEN)x[2])) return 1;
	if (typ(x[1])!=t_REAL || typ(x[2])!=t_REAL) return 0;
	return (expo(x[1])>expo(x[2]));
      }
      if (gcmp0((GEN)x[2]))
      {
	if (typ(x[1])!=t_REAL || typ(x[2])!=t_REAL) return 0;
	return (expo(x[2])>expo(x[1]));
      }
      return 0;

    case t_PADIC:
      return !signe(x[4]);

    case t_QUAD:
      return gcmp0((GEN)x[2]) && gcmp0((GEN)x[3]);

    case t_RFRAC: case t_RFRACN:
      return gcmp0((GEN)x[1]);

    case t_VEC: case t_COL: case t_MAT:
    {
      long i;
      for (i=lg(x)-1; i; i--)
        if (!gcmp0((GEN)x[i])) return 0;
      return 1;
    }
  }
  return 0;
}

/*******************************************************************/
/*                                                                 */
/*                      COMPARE TO 1 and -1                        */
/*     returns 1 whenever the GEN x is 1 (resp. -1), 0 otherwise   */
/*                                                                 */
/*******************************************************************/

int
gcmp1(GEN x)
{
  switch(typ(x))
  {
    case t_INT:
      return is_pm1(x) && signe(x)==1;

    case t_REAL:
      if (signe(x) > 0 && expo(x)==0 && x[2]==HIGHBIT)
      {
        long i,lx = lg(x);
        for (i=3; i<lx; i++)
          if (x[i]) return 0;
        return 1;
      }
      return 0;

    case t_INTMOD: case t_POLMOD:
      return gcmp1((GEN)x[2]);

    case t_FRAC:
      return gcmp1((GEN)x[1]) && gcmp1((GEN)x[2]);

    case t_FRACN:
      return egalii((GEN)x[1],(GEN)x[2]);

    case t_COMPLEX:
      return gcmp1((GEN)x[1]) && gcmp0((GEN)x[2]);

    case t_PADIC:
      return !valp(x) && gcmp1((GEN)x[4]);

    case t_QUAD:
      return gcmp1((GEN)x[2]) && gcmp0((GEN)x[3]);

    case t_POL:
      return lgef(x)==3 && gcmp1((GEN)x[2]);
  }
  return 0;
}

int
gcmp_1(GEN x)
{
  long l,y;
  GEN p1;

  switch(typ(x))
  {
    case t_INT:
      return is_pm1(x) && signe(x)== -1;

    case t_REAL:
      if (signe(x) < 0 && expo(x)==0 && x[2]==HIGHBIT)
      {
        long i,lx = lg(x);
        for (i=3; i<lx; i++)
          if (x[i]) return 0;
        return 1;
      }
      return 0;

    case t_INTMOD:
      l=avma; y=egalii(addsi(1,(GEN)x[2]), (GEN)x[1]); avma=l; return y;

    case t_FRAC: case t_FRACN:
      l = signe(x[1]);
      return l && l == -signe(x[2]) && !absi_cmp((GEN)x[1],(GEN)x[2]);

    case t_COMPLEX:
      return gcmp_1((GEN)x[1]) && gcmp0((GEN)x[2]);

    case t_QUAD:
      return gcmp_1((GEN)x[2]) && gcmp0((GEN)x[3]);

    case t_PADIC:
      l=avma; y=gegal(addsi(1,(GEN)x[4]), (GEN)x[3]); avma=l; return y;

    case t_POLMOD:
      l=avma; p1=gadd(gun,(GEN)x[2]);
      y = signe(p1) && !gegal(p1,(GEN)x[1]); avma=l; return !y;

    case t_POL:
      return lgef(x)==3 && gcmp_1((GEN)x[2]);
  }
  return 0;
}

/*******************************************************************/
/*                                                                 */
/*                       SIGNED COMPARISON                         */
/*     returns the sign of x - y when it makes sense. 0 otherwise  */
/*                                                                 */
/*******************************************************************/

int
gcmp(GEN x, GEN y)
{
  long tx,ty,f,av;

  tx=typ(x); ty=typ(y);
  if (is_intreal_t(tx))
    { if (is_intreal_t(ty)) return mpcmp(x,y); }
  else
  {
    if (tx==t_STR)
    {
      if (ty != t_STR) return 1;
      return strcmp(GSTR(x),GSTR(y));
    }
    if (!is_frac_t(tx)) err(typeer,"comparison");
  }
  if (ty == t_STR) return -1;
  if (!is_intreal_t(ty) && !is_frac_t(ty)) err(typeer,"comparison");
  av=avma; y=gneg_i(y); f=gsigne(gadd(x,y)); avma=av; return f;
}

int
lexcmp(GEN x, GEN y)
{
  const long tx=typ(x), ty=typ(y);
  long lx,ly,l,fl,i;

  ly=lg(y);
  if (!is_matvec_t(tx))
  {
    if (!is_matvec_t(ty)) return gcmp(x,y);
    if (ly==1) return 1;
    fl = lexcmp(x,(GEN)y[1]);
    if (fl) return fl;
    return (ly>2)? -1:0;
  }

  lx=lg(x);
  if (!is_matvec_t(ty))
  {
    if (lx==1) return -1;
    fl = lexcmp(y,(GEN)x[1]);
    if (fl) return -fl;
    return (lx>2)? 1:0;
  }

  /* x and y are matvec_t */
  if (ly==1) return (lx==1)?0:1;
  if (lx==1) return -1;
  if (ty==t_MAT)
  {
    if (tx != t_MAT)
    {
      fl = lexcmp(x,(GEN)y[1]);
      if (fl) return fl;
      return (ly>2)?-1:0;
    }
  }
  else if (tx==t_MAT)
  {
    fl = lexcmp(y,(GEN)x[1]);
    if (fl) return -fl;
    return (ly>2)?1:0;
  }

  /* tx = ty = t_MAT, or x and y are both vect_t */
  l=min(lx,ly);
  for (i=1; i<l; i++)
  {
    fl = lexcmp((GEN)x[i],(GEN)y[i]);
    if (fl) return fl;
  }
  return (ly != lx)? -1 : 0;
}

/*****************************************************************/
/*                                                               */
/*                          EQUALITY                             */
/*                returns 1 if x == y, 0 otherwise               */
/*                                                               */
/*****************************************************************/

static int
polegal(GEN x, GEN y)
{
  long i, lx = lgef(x);

  if (x[1] != y[1] && (lgef(y) != lx || lx > 3)) return 0;
  for (i = 2; i < lx; i++)
    if (!gegal((GEN)x[i],(GEN)y[i])) return 0;
  return 1;
}

#define MASK(x) (((ulong)(x)) & (TYPBITS | LGBITS))
static int
vecegal(GEN x, GEN y)
{
  long i, tx = typ(x);

  if (!is_matvec_t(tx)) return gegal(x,y);
  if (MASK(x[0]) != MASK(y[0])) return 0;

  i = lg(x)-1;
  if (tx != t_MAT)
  {
    for ( ; i; i--)
      if (! gegal((GEN)x[i],(GEN)y[i]) ) return 0;
    return 1;
  }
  for ( ; i; i--)
    if (! vecegal((GEN)x[i],(GEN)y[i]) ) return 0;
  return 1;
}
#undef MASK

#define MASK(x) (((ulong)(x)) & (LGEFINTBITS | SIGNBITS))
int
egalii(GEN x, GEN y)
{
  long i;
  if (MASK(x[1]) != MASK(y[1])) return 0;
  i = lgefint(x)-1; while (i>1 && x[i]==y[i]) i--;
  return i==1;
}
#undef MASK

int
gegal(GEN x, GEN y)
{
  long av,i,tx,ty;
  
  if (x == y) return 1;
  tx = typ(x); ty = typ(y);
  if (tx!=ty)
    { if (tx == t_STR || ty == t_STR) return 0; }
  else
    switch(tx)
    {
      case t_INT:
        return egalii(x,y);

      case t_POL:
        return polegal(x,y);

      case t_COMPLEX:
	return gegal((GEN)x[1],(GEN)y[1]) && gegal((GEN)x[2],(GEN)y[2]);

      case t_INTMOD: case t_POLMOD:
	return gegal((GEN)x[2],(GEN)y[2])
            && (x[1]==y[1] || gegal((GEN)x[1],(GEN)y[1]));

      case t_QFR:
	    if (!gegal((GEN)x[4],(GEN)y[4])) return 0; /* fall through */
      case t_QUAD: case t_QFI:
	return gegal((GEN)x[1],(GEN)y[1])
	    && gegal((GEN)x[2],(GEN)y[2])
	    && gegal((GEN)x[3],(GEN)y[3]);

      case t_FRAC:
	return gegal((GEN)x[1], (GEN)y[1])
            && gegal((GEN)x[2], (GEN)y[2]);

      case t_FRACN: case t_RFRAC: case t_RFRACN:
	av=avma; i=gegal(gmul((GEN)x[1],(GEN)y[2]),gmul((GEN)x[2],(GEN)y[1]));
	avma=av; return i;

      case t_STR:
        return !strcmp(GSTR(x),GSTR(y));

      case t_VEC: case t_COL: case t_MAT:
        return vecegal(x,y);
    }
  av=avma; y=gneg_i(y); i=gcmp0(gadd(x,y));
  avma=av; return i;
}

/*******************************************************************/
/*                                                                 */
/*                          VALUATION                              */
/*             p is either an int or a polynomial.                 */
/*  returns the largest exponent of p dividing x when this makes   */
/*  sense : error for types real, integermod and polymod if p does */
/*  not divide the modulus, q-adic if q!=p.                        */
/*                                                                 */
/*******************************************************************/

static long
minval(GEN x, GEN p, long first, long lx)
{
  long i,k, val = VERYBIGINT;
  for (i=first; i<lx; i++)
  {
    k = ggval((GEN)x[i],p);
    if (k < val) val = k;
  }
  return val;
}

long
ggval(GEN x, GEN p)
{
  long tx=typ(x), tp=typ(p), av, limit,vx,v,i,val;
  GEN p1,p2;

  if (isexactzero(x)) return VERYBIGINT;
  if (is_const_t(tx) && tp==t_POL) return 0;
  switch(tx)
  {
    case t_INT:
      if (tp != t_INT) break;
      av=avma; val = pvaluation(x,p, &p1);
      avma=av; return val;

    case t_INTMOD:
      av=avma;
      p1=cgeti(lgefint(x[1]));
      p2=cgeti(lgefint(x[2]));
      if (tp!=t_INT || !mpdivis((GEN)x[1],p,p1)) break;
      if (!mpdivis((GEN)x[2],p,p2)) { avma=av; return 0; }
      val=1; while (mpdivis(p1,p,p1) && mpdivis(p2,p,p2)) val++;
      avma=av; return val;

    case t_PADIC:
      if (tp!=t_INT || !gegal(p,(GEN)x[2])) break;
      return valp(x);

    case t_POLMOD:
      if (tp==t_INT) return ggval((GEN)x[2],p);
      av = avma;
      if (tp!=t_POL || !poldivis((GEN)x[1],p,&p1)) break;
      if (!poldivis((GEN)x[2],p,&p2)) { avma=av; return 0; }
      val=1; while (poldivis(p1,p,&p1)&&poldivis(p2,p,&p2)) val++;
      avma = av; return val;

    case t_POL:
      if (tp==t_POL)
      {
	v=varn(p); vx=varn(x);
	if (vx == v)
	{
	  if ((p>=(GEN)polx && p <= (GEN)(polx+MAXVARN)) || ismonome(p))
          {
            i=2; while (isexactzero((GEN)x[i])) i++;
            return i-2;
          }
	  av = avma; limit=stack_lim(av,1);
	  for (val=0; ; val++)
	  {
	    if (!poldivis(x,p,&x)) break;
            if (low_stack(limit, stack_lim(av,1)))
	    {
	      if(DEBUGMEM>1) err(warnmem,"ggval");
	      x = gerepileupto(av, gcopy(x));
	    }
	  }
	  avma = av; return val;
	}
	if (vx > v) return 0;
      }
      else
      {
        if (tp!=t_INT) break;
        i=2; while (isexactzero((GEN)x[i])) i++;
      }
      return minval(x,p,i,lgef(x));

    case t_SER:
      if (tp!=t_POL && tp!=t_SER && tp!=t_INT) break;
      v=gvar(p); vx=varn(x);
      if (vx==v) return (long)(valp(x)/ggval(p,polx[v]));
      if (vx>v) return 0;
      return minval(x,p,2,lg(x));

    case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
      return ggval((GEN)x[1],p) - ggval((GEN)x[2],p);

    case t_COMPLEX: case t_QUAD: case t_VEC: case t_COL: case t_MAT:
      return minval(x,p,1,lg(x));
  }
  err(talker,"forbidden or conflicting type in gval");
  return 0; /* not reached */
}

/* x is non-zero */
long
svaluation(ulong x, ulong p, long *py)
{
  ulong v = 0;
  for(;;)
  {
    if (x%p) { *py = x; return v; }
    v++; x/=p;
  }
}

/* x is a non-zero integer */
long
pvaluation(GEN x, GEN p, GEN *py)
{
  long av,v;
  GEN p1,p2;

  if (!is_bigint(x))
  {
    long y;
    if (!is_bigint(p))
    {
      v = svaluation(x[2],p[2], &y);
      if (signe(x) < 0) y = -y;
      *py = stoi(y);
    }
    else
    {
      v = 0;
      *py = icopy(x);
    }
    return v;
  }
  av = avma; v = 0; (void)new_chunk(lgefint(x));
  for(;;)
  {
    p1 = dvmdii(x,p,&p2);
    if (p2 != gzero) { avma=av; *py = icopy(x); return v; }
    v++; x = p1;
  }
}

/*******************************************************************/
/*                                                                 */
/*                       NEGATION: Create -x                       */
/*                                                                 */
/*******************************************************************/

GEN
gneg(GEN x)
{
  long tx=typ(x),lx,i;
  GEN y;

  if (gcmp0(x)) return gcopy(x);
  switch(tx)
  {
    case t_INT: case t_REAL:
      return mpneg(x);

    case t_INTMOD: y=cgetg(3,t_INTMOD);
      icopyifstack(x[1],y[1]);
      y[2] = lsubii((GEN)y[1],(GEN)x[2]);
      break;

    case t_POLMOD: y=cgetg(3,t_POLMOD);
      copyifstack(x[1],y[1]);
      y[2]=lneg((GEN)x[2]); break;

    case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
      y=cgetg(3,tx);
      y[1]=lneg((GEN)x[1]);
      y[2]=lcopy((GEN)x[2]); break;

    case t_PADIC:
      y=cgetp2(x,valp(x));
      y[4]=lsubii((GEN)x[3],(GEN)x[4]);
      break;

    case t_QUAD:
      y=cgetg(4,t_QUAD); copyifstack(x[1],y[1]);
      y[2]=lneg((GEN)x[2]);
      y[3]=lneg((GEN)x[3]); break;

    case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
      lx=lg(x); y=cgetg(lx,tx);
      for (i=1; i<lx; i++) y[i]=lneg((GEN)x[i]);
      break;

    case t_POL:
      lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=lneg((GEN)x[i]);
      break;

    case t_SER:
      lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=lneg((GEN)x[i]);
      break;

    default:
      err(typeer,"negation");
  }
  return y;
}

GEN
gneg_i(GEN x)
{
  long tx=typ(x),lx,i;
  GEN y;

  if (gcmp0(x)) return x;
  switch(tx)
  {
    case t_INT: case t_REAL:
      return mpneg(x);

    case t_INTMOD: y=cgetg(3,t_INTMOD); y[1]=x[1];
      y[2] = lsubii((GEN)y[1],(GEN)x[2]);
      break;

    case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
      y=cgetg(3,tx); y[2]=x[2];
      y[1]=(long)gneg_i((GEN)x[1]); break;

    case t_PADIC: y = cgetg(5,t_PADIC); y[2]=x[2]; y[3]=x[3];
      y[1] = evalprecp(precp(x)) | evalvalp(valp(x));
      y[4]=lsubii((GEN)x[3],(GEN)x[4]); break;

    case t_POLMOD: y=cgetg(3,t_POLMOD); y[1]=x[1];
      y[2]=(long)gneg_i((GEN)x[2]); break;

    case t_QUAD: y=cgetg(4,t_QUAD); y[1]=x[1];
      y[2]=(long)gneg_i((GEN)x[2]);
      y[3]=(long)gneg_i((GEN)x[3]); break;

    case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
      lx=lg(x); y=cgetg(lx,tx);
      for (i=1; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
      break;

    case t_POL: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
      break;

    case t_SER: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
      break;

    default:
      err(typeer,"negation");
  }
  return y;
}

/******************************************************************/
/*                                                                */
/*                       ABSOLUTE VALUE                           */
/*    Create abs(x) if x is integer, real, fraction or complex.   */
/*                       Error otherwise.                         */
/*                                                                */
/******************************************************************/

GEN
gabs(GEN x, long prec)
{
  long tx=typ(x),lx,i,l,tetpil;
  GEN y,p1;

  switch(tx)
  {
    case t_INT: case t_REAL:
      return mpabs(x);

    case t_FRAC: case t_FRACN: y=cgetg(lg(x),tx);
      y[1]=labsi((GEN)x[1]);
      y[2]=labsi((GEN)x[2]); return y;

    case t_COMPLEX:
      l=avma; p1=gnorm(x); tetpil=avma;
      return gerepile(l,tetpil,gsqrt(p1,prec));

    case t_QUAD:
      l=avma; p1=gmul(x, realun(prec)); tetpil=avma;
      return gerepile(l,tetpil,gabs(p1,prec));

    case t_POL:
      lx=lgef(x); if (lx<=2) return gcopy(x);
      p1=(GEN)x[lx-1];
      switch(typ(p1))
      {
	case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
	  if (gsigne(p1)<0) return gneg(x);
      }
      return gcopy(x);

    case t_VEC: case t_COL: case t_MAT:
      lx=lg(x); y=cgetg(lx,tx);
      for (i=1; i<lx; i++)
        y[i]=(long)gabs((GEN)x[i],prec);
      return y;
  }
  err(typeer,"gabs");
  return NULL; /* not reached */
}

GEN
gmax(GEN x, GEN y)
{
  if (gcmp(x,y)>0) y=x; return gcopy(y);
}

GEN
gmin(GEN x, GEN y)
{
  if (gcmp(x,y)<0) y=x; return gcopy(y);
}

GEN
vecmax(GEN x)
{
  long tx=typ(x),lx,lx2,i,j;
  GEN *p1,s;

  if (!is_matvec_t(tx)) return gcopy(x);
  lx=lg(x); if (lx==1) return stoi(-VERYBIGINT);
  if (tx!=t_MAT)
  {
    s=(GEN)x[1];
    for (i=2; i<lx; i++)
      if (gcmp((GEN)x[i],s) > 0) s=(GEN)x[i];
  }
  else
  {
    lx2 = lg(x[1]);
    if (lx2==1) return stoi(-VERYBIGINT);
    s=gmael(x,1,1); i=2;
    for (j=1; j<lx; j++)
    {
      p1 = (GEN *) x[j];
      for (; i<lx2; i++)
	if (gcmp(p1[i],s) > 0) s=p1[i];
      i=1;
    }
  }
  return gcopy(s);
}

GEN
vecmin(GEN x)
{
  long tx=typ(x),lx,lx2,i,j;
  GEN *p1,s;

  if (!is_matvec_t(tx)) return gcopy(x);
  lx=lg(x); if (lx==1) return stoi(VERYBIGINT);
  if (tx!=t_MAT)
  {
    s=(GEN)x[1];
    for (i=2; i<lx; i++)
      if (gcmp((GEN)x[i],s) < 0) s=(GEN)x[i];
  }
  else
  {
    lx2 = lg(x[1]);
    if (lx2==1) return stoi(VERYBIGINT);
    s=gmael(x,1,1); i=2;
    for (j=1; j<lx; j++)
    {
      p1 = (GEN *) x[j];
      for (; i<lx2; i++)
	if (gcmp(p1[i],s) < 0) s=p1[i];
      i=1;
    }
  }
  return gcopy(s);
}

/*******************************************************************/
/*                                                                 */
/*                      AFFECT long --> GEN                        */
/*         affect long s to GEN x. Useful for initialization.      */
/*                                                                 */
/*******************************************************************/

static void
padicaff0(GEN x)
{
  if (signe(x[4]))
  {
    setvalp(x, valp(x)|precp(x));
    affsi(0,(GEN)x[4]);
  }
}

void
gaffsg(long s, GEN x)
{
  long l,i,v;
  GEN p;

  switch(typ(x))
  {
    case t_INT:
      affsi(s,x); break;

    case t_REAL:
      affsr(s,x); break;

    case t_INTMOD:
      modsiz(s,(GEN)x[1],(GEN)x[2]); break;

    case t_FRAC: case t_FRACN:
      affsi(s,(GEN)x[1]); affsi(1,(GEN)x[2]); break;

    case t_COMPLEX:
      gaffsg(s,(GEN)x[1]); gaffsg(0,(GEN)x[2]); break;

    case t_PADIC:
      if (!s) { padicaff0(x); break; }
      p = (GEN)x[2];
      v = (is_bigint(p))? 0: svaluation(s,p[2], &i);
      setvalp(x,v); modsiz(i,(GEN)x[3],(GEN)x[4]);
      break;

    case t_QUAD:
      gaffsg(s,(GEN)x[2]); gaffsg(0,(GEN)x[3]); break;

    case t_POLMOD:
      gaffsg(s,(GEN)x[2]); break;

    case t_POL:
      v=varn(x);
      if (!s) x[1]=evallgef(2) | evalvarn(v);
      else
      {
        x[1]=evalsigne(1) | evallgef(3) | evalvarn(v);
	gaffsg(s,(GEN)x[2]);
      }
      break;

    case t_SER:
      v=varn(x); gaffsg(s,(GEN)x[2]); l=lg(x);
      if (!s)
        x[1] = evalvalp(l-2) | evalvarn(v);
      else
        x[1] = evalsigne(1) | evalvalp(0) | evalvarn(v);
      for (i=3; i<l; i++) gaffsg(0,(GEN)x[i]);
      break;

    case t_RFRAC: case t_RFRACN:
      gaffsg(s,(GEN)x[1]); gaffsg(1,(GEN)x[2]); break;

    case t_VEC: case t_COL: case t_MAT:
      if (lg(x)!=2) err(assigneri,t_INT,typ(x));
      gaffsg(s,(GEN)x[1]); break;

    default: err(typeer,"gaffsg");
  }
}

/*******************************************************************/
/*                                                                 */
/*                     GENERIC AFFECTATION                         */
/*         Affect the content of x to y, whenever possible         */
/*                                                                 */
/*******************************************************************/

GEN
realzero(long prec)
{
  GEN x=cgetr(3);
  x[1]=evalexpo(-bit_accuracy(prec));
  x[2]=0; return x;
}
GEN
realun(long prec)
{
  GEN x=cgetr(prec); affsr(1,x);
  return x;
}

void
gaffect(GEN x, GEN y)
{
  long i,j,k,l,v,vy,lx,ly,tx,ty,av;
  GEN p1,num,den;


  tx=typ(x); ty=typ(y); lx=lg(x); ly=lg(y);
  if (is_scalar_t(tx))
  {
    if (is_scalar_t(ty))
    {
      switch(tx)
      {
	case t_INT:
	  switch(ty)
	  {
	    case t_INT:
	      /* y = gzero, gnil, gun or gdeux */
	      if (is_universal_constant(y))
	      {
		if (y==gzero) err(overwriter,"integer: 0 (gzero)");
		if (y==gun)   err(overwriter,"integer: 1 (gun)");
		if (y==gdeux) err(overwriter,"integer: 2 (gdeux)");
		err(overwriter,"integer: nil (gnil)");
	      }
	      affii(x,y); break;

	    case t_REAL:
	      if (y == gpi || y==geuler) err(overwriter,"real");
	      affir(x,y); break;

	    case t_INTMOD:
	      modiiz(x,(GEN)y[1],(GEN)y[2]); break;

	    case t_FRAC:
	      if (y == ghalf) err(overwriter,"fraction: 1/2 (ghalf)");
	    case t_FRACN:
	      affii(x,(GEN)y[1]); affsi(1,(GEN)y[2]); break;

	    case t_COMPLEX:
	      if (y == gi) err(overwriter,"complex number: I (gi)");
	      gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;

	    case t_PADIC:
              if (!signe(x)) { padicaff0(y); break; }
	      l=avma;
	      setvalp(y, pvaluation(x,(GEN)y[2],&p1));
	      modiiz(p1,(GEN)y[3],(GEN)y[4]);
	      avma=l; break;

	    case t_QUAD:
	      gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;

	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;

	    default: err(typeer,"gaffect");
	  }
	  break;
	
	case t_REAL:
	  switch(ty)
	  {
	    case t_REAL:
	      affrr(x,y); break;

	    case t_COMPLEX:
	      gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;

	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;

	    case t_INT: case t_INTMOD: case t_FRAC:
	    case t_FRACN: case t_PADIC: case t_QUAD:
              err(assignerf,tx,ty);

	    default: err(typeer,"gaffect");
	  }
	  break;
	
	case t_INTMOD:
	  switch(ty)
	  {
	    case t_INTMOD:
	      if (!divise((GEN)x[1],(GEN)y[1]))
                err(assigneri,tx,ty);
	      modiiz((GEN)x[2],(GEN)y[1],(GEN)y[2]); break;

	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;

	    case t_INT: case t_REAL: case t_FRAC:
	    case t_FRACN: case t_COMPLEX: case t_QUAD:
	      err(assignerf,tx,ty);

	    case t_PADIC:
	      err(assignerf,tx,ty);
	    default: err(typeer,"gaffect");
	  }
	  break;

	case t_FRAC:
	  switch(ty)
	  {
	    case t_INT:
	      if (! mpdivis((GEN)x[1],(GEN)x[2],y))
                err(assigneri,tx,ty);
	      break;

	    case t_REAL:
	      diviiz((GEN)x[1],(GEN)x[2],y); break;
	
	    case t_INTMOD:
	      av=avma; p1=mpinvmod((GEN)x[2],(GEN)y[1]);
	      modiiz(mulii((GEN)x[1],p1),(GEN)y[1],(GEN)y[2]);
	      avma=av; break;

	    case t_FRAC: case t_FRACN:
	      affii((GEN)x[1],(GEN)y[1]);
	      affii((GEN)x[2],(GEN)y[2]);
	      break;
	
	    case t_COMPLEX:
	      gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;

	    case t_PADIC:
	      if (!signe(x[1])) { padicaff0(y); break; }
	      av=avma; num=(GEN)x[1]; den=(GEN)x[2];
	      v = pvaluation(num, (GEN) y[2], &num);
	      if (!v) v = -pvaluation(den,(GEN)y[2],&den);
	      p1=mulii(num,mpinvmod(den,(GEN)y[3]));
	      modiiz(p1,(GEN)y[3],(GEN)y[4]); avma=av;
	      setvalp(y,v); break;

	    case t_QUAD:
	      gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;
	    default: err(typeer,"gaffect");
	  }
	  break;

	case t_FRACN:
	  switch(ty)
	  {
	    case t_INT:
	      if (! mpdivis((GEN)x[1],(GEN)x[2],y)) err(assigneri,tx,ty);
	      break;

	    case t_REAL:
	      diviiz((GEN)x[1],(GEN)x[2],y); break;
	    case t_INTMOD:
	      av=avma; gaffect(gred(x),y); avma=av; break;
	    case t_FRAC:
	      gredz(x,y); break;
	    case t_FRACN:
	      affii((GEN)x[1],(GEN)y[1]);
	      affii((GEN)x[2],(GEN)y[2]); break;
	    case t_COMPLEX:
	      gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
	    case t_PADIC:
	      if (!signe(x[1])) { padicaff0(y); break; }
	      av=avma; num=(GEN)x[1]; den=(GEN)x[2];
	      v =  pvaluation(num,(GEN)y[2],&num)
	         - pvaluation(den,(GEN)y[2],&den);
	      p1=mulii(num,mpinvmod(den,(GEN)y[3]));
	      modiiz(p1,(GEN)y[3],(GEN)y[4]); avma=av;
	      setvalp(y,v); break;

	    case t_QUAD:
	      gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;
	    default: err(typeer,"gaffect");
	  }
	  break;
	
	case t_COMPLEX:
	  switch(ty)
	  {
	    case t_INT: case t_REAL: case t_INTMOD:
	    case t_FRAC: case t_FRACN: case t_PADIC: case t_QUAD:
	      if (!gcmp0((GEN)x[2])) err(assigneri,tx,ty);
	      gaffect((GEN)x[1],y); break;
	
	    case t_COMPLEX:
	      gaffect((GEN)x[1],(GEN)y[1]);
	      gaffect((GEN)x[2],(GEN)y[2]);
	      break;

	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;

	    default: err(typeer,"gaffect");
	  }
	  break;
	
	case t_PADIC:
	  switch(ty)
	  {
	    case t_INTMOD:
	      if (valp(x)<0) err(assigneri,tx,ty);
	      av=avma;
	      v = pvaluation((GEN)y[1],(GEN)x[2],&p1);
              k = signe(x[4])? (precp(x)+valp(x)): valp(x)+1;
	      if (!gcmp1(p1) || v > k) err(assigneri,tx,ty);
	      modiiz(gtrunc(x),(GEN)y[1],(GEN)y[2]); avma=av; break;
	
	    case t_PADIC:
	      if (!egalii((GEN)x[2],(GEN)y[2])) err(assigneri,tx,ty);
	      modiiz((GEN)x[4],(GEN)y[3],(GEN)y[4]);
	      setvalp(y,valp(x)); break;
	
	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;

	    case t_INT: case t_REAL: case t_FRAC:
	    case t_FRACN: case t_COMPLEX: case t_QUAD:
              err(assignerf,tx,ty);

	    default: err(typeer,"gaffect");
	  }
	  break;
	
	case t_QUAD:
	  switch(ty)
	  {
	    case t_INT: case t_INTMOD: case t_FRAC:
	    case t_FRACN: case t_PADIC:
	      if (!gcmp0((GEN)x[3])) err(assigneri,tx,ty);
	      gaffect((GEN)x[2],y); break;

	    case t_REAL:
	      av=avma; p1=co8(x,ly); gaffect(p1,y); avma=av; break;
	    case t_COMPLEX:
	      ly=precision(y);
	      if (ly)
	      {
		av=avma; p1=co8(x,ly); gaffect(p1,y); avma=av;
	      }
	      else
	      {
                if (!gcmp0((GEN)x[3])) err(assigneri,tx,ty);
		gaffect((GEN)x[2],y);
	      }
	      break;

	    case t_QUAD:
	      if (! gegal((GEN)x[1],(GEN)y[1])) err(assigneri,tx,ty);
	      affii((GEN)x[2],(GEN)y[2]);
	      affii((GEN)x[3],(GEN)y[3]);
	      break;
	    case t_POLMOD:
	      gaffect(x,(GEN)y[2]); break;
	    default: err(typeer,"gaffect");
	  }
	  break;

	case t_POLMOD:
	  if (ty!=t_POLMOD) err(assignerf,tx,ty);
	  if (! gdivise((GEN)x[1],(GEN)y[1])) err(assigneri,tx,ty);
	  gmodz((GEN)x[2],(GEN)y[1],(GEN)y[2]); break;
	default: err(typeer,"gaffect");
      }
      return;
    }

    /* here y is not scalar */
    switch(ty)
    {
      case t_POL:
	v=varn(y);
	if (y==polun[v] || y==polx[v])
	  err(talker,"trying to overwrite a universal polynomial");
	gaffect(x,(GEN)y[2]);
	for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
	if (gcmp0(x)) y[1]=evallgef(2)+evalvarn(v);
	else y[1]=evalsigne(1)+evallgef(3)+evalvarn(v);
	break;

      case t_SER:
	v=varn(y); gaffect(x,(GEN)y[2]);
	if (gcmp0(x))
          y[1] = evalvalp(ly-2) | evalvarn(v);
	else
	  y[1] = evalsigne(1) | evalvalp(0) | evalvarn(v);
        for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
	break;

      case t_RFRAC: case t_RFRACN:
	gaffect(x,(GEN)y[1]); gaffsg(1,(GEN)y[2]); break;

      case t_VEC: case t_COL: case t_MAT:
	err(assignerf,tx,ty);

      default: err(typeer,"gaffect");
    }
    return;
  }

  if (is_const_t(ty)) {
      /* initial_value macro is not defined now... */
#define initial_value(ep) ((ep)+1)
      if (tx == t_POL) {
	  entree *var = varentries[ordvar[varn(x)]];
	  /* Is a bound expression, thus should not loop: */
	  if (var && var->value != (void*)initial_value(var)) {
	      gaffect(geval(x),y);
	      return;
	  }	  
      } else if (is_rfrac_t(tx)) {
	  GEN num = (GEN)x[1], den = (GEN)x[2];
	  entree *varnum = varentries[ordvar[varn(num)]];
	  entree *varden = varentries[ordvar[varn(den)]];

	  /* Are bound expressions, thus should not loop: */
	  if (varnum && varnum->value != (void*)initial_value(varnum)
	   && varden && varden->value != (void*)initial_value(varden)) {
	      gaffect(geval(x),y);
	      return;
	  }
      }
#undef initial_value
      err(assignerf,tx,ty);
  }
  
  switch(tx)
  {
    case t_POL:
      v=varn(x);
      switch(ty)
      {
	case t_POL:
	  vy=varn(y); if (vy>v) err(assignerf,tx,ty);
	  if (vy==v)
	  {
	    l=lgef(x); if (l>ly) err(assigneri,tx,ty);
	    y[1]=x[1]; for (i=2; i<l; i++) gaffect((GEN)x[i],(GEN)y[i]);
	  }
	  else
	  {
	    gaffect(x,(GEN)y[2]);
	    for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
	    if (signe(x)) y[1]=evalsigne(1)+evallgef(3)+evalvarn(vy);
	    else y[1]=evallgef(2)+evalvarn(vy);
	  }
	  break;
	
	case t_SER:
	  vy=varn(y); if (vy>v) err(assignerf,tx,ty);
	  if (!signe(x)) { gaffsg(0,y); return; }
	  if (vy==v)
	  {
	    i=gval(x,v); y[1]=evalvarn(v) | evalvalp(i) | evalsigne(1);
	    k=lgef(x)-i; if (k>ly) k=ly;
	    for (j=2; j<k; j++) gaffect((GEN)x[i+j],(GEN)y[j]);
	    for (   ; j<ly; j++) gaffsg(0,(GEN)y[j]);
	  }
	  else
	  {
	    gaffect(x,(GEN)y[2]);
	    if (!signe(x))
              y[1] = evalvalp(ly-2) | evalvarn(vy);
	    else
	      y[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
            for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
	  }
	  break;
	
	case t_POLMOD:
	  gmodz(x,(GEN)y[1],(GEN)y[2]); break;

	case t_RFRAC: case t_RFRACN:
	  gaffect(x,(GEN)y[1]); gaffsg(1,(GEN)y[2]); break;

	case t_VEC: case t_COL: case t_MAT:
	  err(assignerf,tx,ty);
	
	default: err(typeer,"gaffect");
      }
      break;

    case t_SER:
      if (ty!=t_SER) err(assignerf,tx,ty);
      v=varn(x); vy=varn(y); if (vy>v) err(assignerf,tx,ty);
      if (vy==v)
      {
	y[1]=x[1]; k=lx; if (k>ly) k=ly;
        for (i=2; i< k; i++) gaffect((GEN)x[i],(GEN)y[i]);
        for (   ; i<ly; i++) gaffsg(0,(GEN)y[i]);
      }
      else
      {
	gaffect(x,(GEN)y[2]);
	if (!signe(x))
          y[1] = evalvalp(ly-2) | evalvarn(vy);
	else
	  y[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
        for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
      }
      break;

    case t_RFRAC: case t_RFRACN:
      switch(ty)
      {
	case t_POL: case t_VEC: case t_COL: case t_MAT:
	  err(assignerf,tx,ty);

	case t_POLMOD:
	  av=avma; p1=ginvmod((GEN)x[2],(GEN)y[1]);
	  gmodz(gmul((GEN)x[1],p1),(GEN)y[1],(GEN)y[2]);
	  avma=av; break;

	case t_SER:
	  gdivz((GEN)x[1],(GEN)x[2],y); break;
	
	case t_RFRAC:
	  if (tx==t_RFRACN) gredz(x,y);
	  else
	  {
	    gaffect((GEN)x[1],(GEN)y[1]);
	    gaffect((GEN)x[2],(GEN)y[2]);
	  }
	  break;
	
	case t_RFRACN:
	  gaffect((GEN)x[1],(GEN)y[1]);
	  gaffect((GEN)x[2],(GEN)y[2]); break;
	
	default: err(typeer,"gaffect");
      }
      break;

    case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
      if (ty != tx || ly != lx) err(assigneri,tx,ty);
      for (i=1; i<lx; i++) gaffect((GEN)x[i],(GEN)y[i]);
      break;

    default: err(typeer,"gaffect");
  }
}

/*******************************************************************/
/*                                                                 */
/*           CONVERSION QUAD --> REAL, COMPLEX OR P-ADIC           */
/*                                                                 */
/*******************************************************************/

GEN
co8(GEN x, long prec)
{
  long av=avma,tetpil;
  GEN p1, p = (GEN) x[1];

  p1 = subii(sqri((GEN)p[3]), shifti((GEN)p[2],2));
  if (signe(p1) > 0)
  {
    p1 = subri(gsqrt(p1,prec), (GEN)p[3]);
    setexpo(p1, expo(p1)-1);
  }
  else
  {
    p1 = gsub(gsqrt(p1,prec), (GEN)p[3]);
    p1[1] = lmul2n((GEN)p1[1],-1);
    setexpo(p1[2], expo(p1[2])-1);
  }/* p1 = (-b + sqrt(D)) / 2 */
  p1 = gmul((GEN)x[3],p1); tetpil=avma;
  return gerepile(av,tetpil,gadd((GEN)x[2],p1));
}

GEN
cvtop(GEN x, GEN p, long l)
{
  GEN p1,p2,p3;
  long av,tetpil,n;

  if (typ(p)!=t_INT)
    err(talker,"not an integer modulus in cvtop or gcvtop");
  if (gcmp0(x)) return ggrandocp(p,l);
  switch(typ(x))
  {
    case t_INT:
      return gadd(x,ggrandocp(p,ggval(x,p)+l));

    case t_INTMOD:
      n=ggval((GEN)x[1],p); if (n>l) n=l;
      return gadd((GEN)x[2],ggrandocp(p,n));

    case t_FRAC: case t_FRACN:
      n = ggval((GEN)x[1],p) - ggval((GEN)x[2],p);
      return gadd(x,ggrandocp(p,n+l));

    case t_COMPLEX:
      av=avma; p1=gsqrt(gaddgs(ggrandocp(p,l),-1),0);
      p1=gmul(p1,(GEN)x[2]); tetpil=avma;
      return gerepile(av,tetpil,gadd(p1,(GEN)x[1]));

    case t_PADIC:
      return gprec(x,l);

    case t_QUAD:
      av=avma; p1=(GEN)x[1]; p3=gmul2n((GEN)p1[3],-1);
      p2=gsub(gsqr(p3),(GEN)p1[2]);
      switch(typ(p2))
      {
	case t_INT:
	  n=ggval(p2,p); p2=gadd(p2,ggrandocp(p,n+l)); break;

	case t_INTMOD: case t_PADIC:
	  break;

	case t_FRAC: case t_FRACN:
	  n = ggval((GEN)p2[1],p) - ggval((GEN)p2[2],p);
	  p2=gadd(p2,ggrandocp(p,n+l)); break;

	default: err(assigneri,t_QUAD,t_QUAD);
      }
      p2=gsqrt(p2,0); p1=gmul((GEN)x[3],gsub(p2,p3)); tetpil=avma;
      return gerepile(av,tetpil,gadd((GEN)x[2],p1));

  }
  err(typeer,"cvtop");
  return NULL; /* not reached */
}

GEN
gcvtop(GEN x, GEN p, long r)
{
  long i,lx, tx=typ(x);
  GEN y;

  if (is_const_t(tx)) return cvtop(x,p,r);
  switch(tx)
  {
    case t_POL:
      lx=lgef(x); y=cgetg(lx,t_POL); y[1]=x[1];
      for (i=2; i<lx; i++)
	y[i]=(long)gcvtop((GEN)x[i],p,r);
      break;

    case t_SER:
      lx=lg(x); y=cgetg(lx,t_SER); y[1]=x[1];
      for (i=2; i<lx; i++)
	y[i]=(long)gcvtop((GEN)x[i],p,r);
      break;

    case t_POLMOD: case t_RFRAC: case t_RFRACN:
    case t_VEC: case t_COL: case t_MAT:
      lx=lg(x); y=cgetg(lx,tx);
      for (i=1; i<lx; i++)
	y[i]=(long)gcvtop((GEN)x[i],p,r);
      break;

    default: err(typeer,"gcvtop");
  }
  return y;
}

long
gexpo(GEN x)
{
  long tx=typ(x),lx,e,i,y,av;

  switch(tx)
  {
    case t_INT:
      return expi(x);

    case t_FRAC: case t_FRACN:
      if (!signe(x[1])) return -HIGHEXPOBIT;
      return expi((GEN)x[1]) - expi((GEN)x[2]);

    case t_REAL:
      return expo(x);

    case t_COMPLEX:
      e = gexpo((GEN)x[1]);
      y = gexpo((GEN)x[2]); return max(e,y);

    case t_QUAD:
      av=avma; x = co8(x,3); avma=av; return gexpo(x);

    case t_POL: case t_SER: case t_VEC: case t_COL: case t_MAT:
      lx=(tx==t_POL)? lgef(x): lg(x);
      y = -HIGHEXPOBIT;
      for (i=lontyp[tx]; i<lx; i++) { e=gexpo((GEN)x[i]); if (e>y) y=e; }
      return y;
  }
  err(typeer,"gexpo");
  return 0; /* not reached */
}

long
gsize(GEN x)
{
  return gcmp0(x)? 0: (long) ((gexpo(x)+1) * L2SL10) + 1;
}

/* Normalize series x in place.
 * Assumption: x,x[2],...,x[lg(x)-1] have been created in that order.
 * All intermediate objects will be destroyed.
 */
GEN
normalize(GEN x)
{
  long i,j, lx = lg(x);

  if (typ(x)!=t_SER) err(typeer,"normalize");
  if (lx==2) { setsigne(x,0); avma = (long) x; return x; }
  if (! isexactzero((GEN)x[2])) { setsigne(x,1); return x; }

  for (i=3; i<lx; i++)
    if (! isexactzero((GEN)x[i]))
    {
      long tetpil = avma;
      GEN p1 = cgetg(lx-i+2,t_SER);
      p1[1] = evalsigne(1) | evalvalp(valp(x)+i-2) | evalvarn(varn(x));
      j=i; i=2; while (j<lx) p1[i++] = lcopy((GEN)x[j++]);
      return gerepile((long) (x+lx),tetpil,p1);
    }
  avma = (long) (x+lx); return zeroser(varn(x),lx-2);
}

GEN
normalizepol_i(GEN x, long lx)
{
  long i;
  for (i=lx-1; i>1; i--)
    if (! isexactzero((GEN)x[i])) break;
  setlgef(x,i+1);
  for (; i>1; i--)
    if (! gcmp0((GEN)x[i]) ) { setsigne(x,1); return x; }
  setsigne(x,0); return x;
}

/* Normalize polynomial x in place. See preceding comment */
GEN
normalizepol(GEN x)
{
  if (typ(x)!=t_POL) err(typeer,"normalizepol");
  return normalizepol_i(x, lgef(x));
}

int
gsigne(GEN x)
{
  switch(typ(x))
  {
    case t_INT: case t_REAL:
      return signe(x);

    case t_FRAC: case t_FRACN:
      return (signe(x[2])>0) ? signe(x[1]) : -signe(x[1]);
  }
  err(typeer,"gsigne");
  return 0; /* not reached */
}

int ff_poltype(GEN *x, GEN *p, GEN *pol);

GEN
gsqr(GEN x)
{
  long tx=typ(x),lx,i,j,k,l,av,tetpil;
  GEN z,p1,p2,p3,p4;

  if (is_scalar_t(tx))
    switch(tx)
    {
      case t_INT:
	return sqri(x);
	
      case t_REAL:
	return mulrr(x,x);
	
      case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
        (void)new_chunk(lgefint(p2)<<2);
        p1=sqri((GEN)x[2]); avma=(long)z;
        z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;

      case t_FRAC: case t_FRACN:
	z=cgetg(3,tx);
	z[1]=lsqri((GEN)x[1]);
	z[2]=lsqri((GEN)x[2]);
	return z; /* reduction is useless ! */

      case t_COMPLEX:
	z=cgetg(lg(x),tx); l=avma;
	p1=gadd((GEN)x[1],(GEN)x[2]);
	p2=gadd((GEN)x[1],gneg_i((GEN)x[2]));
	p3=gmul((GEN)x[1],(GEN)x[2]);
	tetpil=avma;
	z[1]=lmul(p1,p2); z[2]=lshift(p3,1);
	gerepilemanyvec(l,tetpil,z+1,2);
	return z;
	
      case t_PADIC:
	z=cgetg(5,t_PADIC);
	z[2] = lcopy((GEN)x[2]);
	if (!cmpsi(2,(GEN)x[2]))
	{
	  i=precp(x)+1; av=avma;
	  p1=addii((GEN)x[3],shifti((GEN)x[4],1));
	  if (!gcmp0(p1)) { j=vali(p1); if (j<i) i=j; }
	  avma=av;
          z[3]=lshifti((GEN)x[3],i);
	  z[1]=evalprecp(precp(x)+i) | evalvalp(2*valp(x));
	}
	else
	{
	  z[3]=licopy((GEN)x[3]);
	  z[1]=evalprecp(precp(x)) | evalvalp(2*valp(x));
	}
	z[4]=lgeti(lg(z[3])); av=avma;
	modiiz(sqri((GEN)x[4]),(GEN)z[3],(GEN)z[4]);
	avma=av; return z;
	
      case t_QUAD:
	p1=(GEN)x[1]; z=cgetg(lg(x),tx); l=avma;
	p2=gsqr((GEN)x[2]); p3=gsqr((GEN)x[3]);
	p4=gmul(gneg_i((GEN)p1[2]),p3);

	if (gcmp0((GEN)p1[3]))
	{
	  tetpil=avma;
	  z[2]=lpile(l,tetpil,gadd(p4,p2));
	  l=avma; p2=gmul((GEN)x[2],(GEN)x[3]); tetpil=avma;
	  z[3]=lpile(l,tetpil,gmul2n(p2,1));
	  copyifstack(p1,z[1]); return z;
	}

	p1=gmul((GEN)x[2],(GEN)x[3]);
	p1=gmul2n(p1,1); tetpil=avma;
	z[2]=ladd(p2,p4); z[3]=ladd(p1,p3);
	gerepilemanyvec(l,tetpil,z+2,2);
	copyifstack(x[1],z[1]); return z;

      case t_POLMOD:
        z=cgetg(lg(x),tx); copyifstack(x[1],z[1]);
	l=avma; p1=gsqr((GEN)x[2]); tetpil=avma;
        z[2]=lpile(l,tetpil, gres(p1,(GEN)z[1]));
	return z;
    }

  switch(tx)
  {
    case t_POL:
    {
      GEN a = x, p = NULL, pol = NULL;
      long vx = varn(x);
      av = avma;
      if (ff_poltype(&x,&p,&pol))
      {
        z = quicksqr(x+2, lgef(x)-2);
        if (p) z = Fp_pol(z,p);
        if (pol) z = from_Kronecker(z,pol);
        z = gerepileupto(av, z);
      }
      else
      {
        avma = av;
        z = quicksqr(a+2, lgef(a)-2);
      }
      setvarn(z, vx); return z;
    }

    case t_SER:
      if (gcmp0(x)) return zeroser(varn(x), 2*valp(x));
      lx = lg(x); z = cgetg(lx,tx);
      z[1] = evalsigne(1) | evalvalp(2*valp(x)) | evalvarn(varn(x));
      x += 2; z += 2; lx -= 3;
      p2 = (GEN)gpmalloc((lx+1)*sizeof(long));
      for (i=0; i<=lx; i++)
      {
	p2[i] = !isexactzero((GEN)x[i]);
        p1=gzero; av=avma; l=(i+1)>>1;
        for (j=0; j<l; j++)
          if (p2[j] && p2[i-j])
            p1 = gadd(p1, gmul((GEN)x[j],(GEN)x[i-j]));
        p1 = gshift(p1,1);
        if ((i&1) == 0 && p2[i>>1])
          p1 = gadd(p1, gsqr((GEN)x[i>>1]));
        z[i] = lpileupto(av,p1);
      }
      z -= 2; free(p2); return normalize(z);

    case t_RFRAC: case t_RFRACN:
      z=cgetg(3,tx);
      z[1]=lsqr((GEN)x[1]);
      z[2]=lsqr((GEN)x[2]); return z;

    case t_QFR: return sqcompreal(x);
    case t_QFI: return sqcompimag(x);
	
    case t_MAT:
      lx=lg(x);
      if (lx==1) return cgetg(1,tx);
      if (lx != lg(x[1])) err(gmuleri,tx,tx);
      z=cgetg(lx,tx);
      for (j=1; j<lx; j++)
      {
        z[j]=lgetg(lx,t_COL);
        for (i=1; i<lx; i++)
        {
          p1=gzero; l=avma;
          for (k=1; k<lx; k++)
          {
            p2=gmul(gcoeff(x,i,k),gcoeff(x,k,j));
            tetpil=avma; p1=gadd(p1,p2);
          }
          coeff(z,i,j)=lpile(l,tetpil,p1);
        }
      }
      return z;

    case t_VEC: case t_COL:
      err(gmuleri,tx,tx);
  }
  err(typeer,"gsqr");
  return NULL; /* not reached */
}

/*******************************************************************/
/*                                                                 */
/*                              LISTS                              */
/*                                                                 */
/*******************************************************************/

GEN
listcreate(long n)
{
  GEN list;

  if (n<0) err(talker,"negative length in listcreate");
  list=cgetg(n+2,t_LIST); list[1]=evallgef(2);
  return list;
}

static void
listaffect(GEN list, long index, GEN object)
{
  if (index<lgef(list) && isclone(list[index]))
    gunclone((GEN)list[index]);
  list[index]=lclone(object);
}

void
listkill(GEN list)
{
  long lx,i;

  if (typ(list)!=t_LIST) err(typeer,"listkill");
  lx=lgef(list);
  for (i=2; i<lx; i++)
    if (isclone(list[i])) gunclone((GEN)list[i]);
  list[1]=evallgef(2); return;
}

GEN
listput(GEN list, GEN object, long index)
{
  long lx=lgef(list);

  if (typ(list)!=t_LIST) err(typeer,"listput");
  if (index < 0) err(talker,"negative index (%ld) in listput", index);
  if (!index || index >= lx-1)
  {
    index = lx-1; lx++;
    if (lx > lg(list))
      err(talker,"no more room in this list (size %ld)", lg(list)-2);
  }
  listaffect(list,index+1,object);
  list[1]=evallgef(lx);
  return (GEN)list[index+1];
}

GEN
listinsert(GEN list, GEN object, long index)
{
  long lx = lgef(list), i;

  if (typ(list)!=t_LIST) err(typeer,"listinsert");
  if (index <= 0 || index > lx-1) err(talker,"bad index in listinsert");
  lx++; if (lx > lg(list)) err(talker,"no more room in this list");

  for (i=lx-2; i > index; i--) list[i+1]=list[i];
  list[index+1] = lclone(object);
  list[1] = evallgef(lx); return (GEN)list[index+1];
}

GEN
gtolist(GEN x)
{
  long tx,lx,i;
  GEN list;

  if (!x) { list = cgetg(2, t_LIST); list[1] = evallgef(2); return list; }

  tx = typ(x);
  lx = (tx==t_LIST)? lgef(x): lg(x);
  switch(tx)
  {
    case t_VEC: case t_COL:
      lx++; x--; /* fall through */
    case t_LIST:
      list = cgetg(lx,t_LIST);
      for (i=2; i<lx; i++) list[i] = lclone((GEN)x[i]);
      break;
    default: err(typeer,"gtolist");
  }
  list[1] = evallgef(lx); return list;
}

GEN
listconcat(GEN list1, GEN list2)
{
  long i,l1,lx;
  GEN list;

  if (typ(list1)!=t_LIST || typ(list2)!=t_LIST)
    err(typeer,"listconcat");
  l1=lgef(list1)-2; lx=l1+lgef(list2);
  list = listcreate(lx-2);
  for (i=2; i<=l1+1; i++) listaffect(list,i,(GEN)list1[i]);
  for (   ; i<lx; i++) listaffect(list,i,(GEN)list2[i-l1]);
  list[1]=evallgef(lx); return list;
}

GEN
listsort(GEN list, long flag)
{
  long i,av=avma, c=list[1], lx = lgef(list)-1;
  GEN perm,vec,l;

  if (typ(list) != t_LIST) err(typeer,"listsort");
  vec = list+1;
  vec[0] = evaltyp(t_VEC) | evallg(lx);
  perm = sindexsort(vec);
  list[1] = c; l = cgetg(lx,t_VEC);
  for (i=1; i<lx; i++) l[i] = vec[perm[i]];
  if (flag)
  {
    c=1; vec[1] = l[1];
    for (i=2; i<lx; i++)
      if (!gegal((GEN)l[i], (GEN)vec[c]))
        vec[++c] = l[i];
      else
        if (isclone(l[i])) gunclone((GEN)l[i]);
    setlgef(list,c+2);
  }
  else
    for (i=1; i<lx; i++) vec[i] = l[i];

  avma=av; return list;
}