===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/gen3.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/gen3.c	2001/10/02 11:17:04	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/gen3.c	2002/09/11 07:26:51	1.2
@@ -1,4 +1,4 @@
-/* $Id: gen3.c,v 1.1 2001/10/02 11:17:04 noro Exp $
+/* $Id: gen3.c,v 1.2 2002/09/11 07:26:51 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -197,14 +197,13 @@ padicprec(GEN x, GEN p)
   return 0; /* not reached */
 }
 
-/* Degre de x par rapport a la variable v si v>=0, par rapport a la variable
- * principale si v<0. On suppose x fraction rationnelle ou polynome.
- * Convention deg(0)=-1.
+/* Degree of x (scalar, t_POL, t_RFRAC) wrt variable v if v >= 0,
+ * wrt to main variable if v < 0.  Convention: deg(0) = -1.
  */
 long
 poldegree(GEN x, long v)
 {
-  long tx=typ(x), av, w, d;
+  long tx = typ(x), lx,w,i,d;
 
   if (is_scalar_t(tx)) return gcmp0(x)? -1: 0;
   switch(tx)
@@ -213,9 +212,13 @@ poldegree(GEN x, long v)
       w = varn(x);
       if (v < 0 || v == w) return degpol(x);
       if (v < w) return signe(x)? 0: -1;
-      av = avma; x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
-      if (gvar(x)) { d = gcmp0(x)? -1: 0; } else d = degpol(x);
-      avma = av; return d;
+      lx = lgef(x); d = -1;
+      for (i=2; i<lx; i++)
+      {
+        long e = poldegree((GEN)x[i], v);
+        if (e > d) d = e;
+      }
+      return d;
 
     case t_RFRAC: case t_RFRACN:
       if (gcmp0((GEN)x[1])) return -1;
@@ -237,7 +240,8 @@ degree(GEN x)
 GEN
 pollead(GEN x, long v)
 {
-  long l,tx = typ(x),av,tetpil,w;
+  long l, tx = typ(x), w;
+  gpmem_t av, tetpil;
   GEN xinit;
 
   if (is_scalar_t(tx)) return gcopy(x);
@@ -375,12 +379,13 @@ ismonome(GEN x)
 
 /* assume z[1] was created last */
 #define fix_frac_if_int(z) if (is_pm1(z[2]))\
-  z = gerepileuptoint((long)(z+3), (GEN)z[1]);
+  z = gerepileuptoint((gpmem_t)(z+3), (GEN)z[1]);
 
 GEN
 gmulsg(long s, GEN y)
 {
-  long ty=typ(y),ly=lg(y),i,av,tetpil;
+  long ty=typ(y), ly=lg(y), i;
+  gpmem_t av, tetpil;
   GEN z,p1,p2;
 
   switch(ty)
@@ -393,7 +398,7 @@ gmulsg(long s, GEN y)
 
     case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
       (void)new_chunk(lgefint(p2)<<2); /* HACK */
-      p1=mulsi(s,(GEN)y[2]); avma=(long)z;
+      p1=mulsi(s,(GEN)y[2]); avma=(gpmem_t)z;
       z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
 
     case t_FRAC:
@@ -453,7 +458,7 @@ gmulsg(long s, GEN y)
       if (!s) return zeropol(gvar(y));
       z = cgetg(3, t_RFRAC);
       i = ggcd(stoi(s),(GEN)y[2])[2];
-      avma = (long)z;
+      avma = (gpmem_t)z;
       if (i == 1)
       {
         z[1]=lmulgs((GEN)y[1], s);
@@ -490,8 +495,9 @@ gmulsg(long s, GEN y)
 GEN
 gdivgs(GEN x, long s)
 {
-  static long court[] = { evaltyp(t_INT) | m_evallg(3),0,0 };
-  long tx=typ(x),lx=lg(x),i,av;
+  static long court[] = { evaltyp(t_INT) | _evallg(3),0,0 };
+  long tx=typ(x), lx=lg(x), i;
+  gpmem_t av;
   GEN z,p1;
 
   if (!s) err(gdiver2);
@@ -560,20 +566,36 @@ gdivgs(GEN x, long s)
       z[1]=x[1]; return z;
 
     case t_RFRAC:
-      z = cgetg(3, t_RFRAC);
-      i = ggcd(stoi(s),(GEN)x[1])[2];
-      avma = (long)z;
-      if (i == 1)
+      av = avma;
+      p1 = ggcd(stoi(s),(GEN)x[1]);
+      if (typ(p1) == t_INT)
       {
-        z[2]=lmulsg(s,(GEN)x[2]);
-        z[1]=lcopy((GEN)x[1]); return z;
+        avma = av;
+        z = cgetg(3, t_RFRAC);
+        i = p1[2];
+        if (i == 1)
+        {
+          z[1] = lcopy((GEN)x[1]);
+          z[2] = lmulsg(s,(GEN)x[2]);
+        }
+        else
+        {
+          z[1] = ldivgs((GEN)x[1], i);
+          z[2] = lmulgs((GEN)x[2], s/i);
+        }
       }
-      z[1] = ldivgs((GEN)x[1], i);
-      z[2] = lmulgs((GEN)x[2], s/i); return z;
+      else /* t_FRAC */
+      {
+        z = cgetg(3, t_RFRAC);
+        z[1] = ldiv((GEN)x[1], p1);
+        z[2] = lmul((GEN)x[2], gdivsg(s,p1));
+        z = gerepilecopy(av, z);
+      }
+      return z;
 
     case t_RFRACN: z=cgetg(3,t_RFRACN);
-      z[2]=lmulsg(s,(GEN)x[2]);
-      z[1]=lcopy((GEN)x[1]); return z;
+      z[1]=lcopy((GEN)x[1]);
+      z[2]=lmulsg(s,(GEN)x[2]); return z;
 
     case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
       for (i=1; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
@@ -584,22 +606,38 @@ gdivgs(GEN x, long s)
 
 /*******************************************************************/
 /*                                                                 */
-/*                    MODULO GENERAL                               */
+/*                    GENERIC REMAINDER                            */
 /*                                                                 */
 /*******************************************************************/
+/* euclidean quotient for scalars of admissible types */
+static GEN
+_quot(GEN x, GEN y)
+{
+  GEN q = gdiv(x,y), f = gfloor(q);
+  if (gsigne(y) < 0 && !gegal(f,q)) f = gadd(f,gun);
+  return f;
+}
+static GEN
+quot(GEN x, GEN y)
+{
+  gpmem_t av = avma;
+  return gerepileupto(av, _quot(x, y));
+}
 
 GEN
 gmod(GEN x, GEN y)
 {
-  long av,tetpil,i, tx=typ(x), ty = typ(y);
+  gpmem_t av, tetpil;
+  long i,lx,ty, tx = typ(x);
   GEN z,p1;
 
   if (is_matvec_t(tx))
   {
-    i=lg(x); z=cgetg(i,tx);
-    for (i--; i; i--) z[i]=lmod((GEN)x[i],y);
+    lx=lg(x); z=cgetg(lx,tx);
+    for (i=1; i<lx; i++) z[i] = lmod((GEN)x[i],y);
     return z;
   }
+  ty = typ(y);
   switch(ty)
   {
     case t_INT:
@@ -612,8 +650,11 @@ gmod(GEN x, GEN y)
           z[1]=lmppgcd((GEN)x[1],y);
 	  z[2]=lmodii((GEN)x[2],(GEN)z[1]); return z;
 
-	case t_FRAC: case t_FRACN:
-	  av=avma; if (tx==t_FRACN) x=gred(x);
+        case t_FRACN:
+	  av=avma; x = gred(x);
+          return gerepileupto(av, gmod(x,y));
+	case t_FRAC:
+	  av=avma;
 	  p1=mulii((GEN)x[1],mpinvmod((GEN)x[2],y));
 	  tetpil=avma; return gerepile(av,tetpil,modii(p1,y));
 
@@ -624,35 +665,37 @@ gmod(GEN x, GEN y)
 
 	case t_PADIC:
         {
-          long p1[] = {evaltyp(t_INTMOD)|m_evallg(3),0,0};
+          long p1[] = {evaltyp(t_INTMOD)|_evallg(3),0,0};
           p1[1] = (long)y;
           p1[2] = lgeti(lgefint(y));
           gaffect(x,p1); return (GEN)p1[2];
         }
 	case t_POLMOD: case t_POL:
 	  return gzero;
+      /* case t_REAL could be defined as below, but conlicting semantic
+       * with lift(x * Mod(1,y)). */
 
-	default: err(operf,"%",tx,ty);
+	default: err(operf,"%",x,y);
       }
 
     case t_REAL: case t_FRAC: case t_FRACN:
       switch(tx)
       {
 	case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
-	  av=avma; p1 = gfloor(gdiv(x,y)); p1 = gneg_i(gmul(p1,y));
-          tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
+	  av = avma;
+          return gerepileupto(av, gsub(x, gmul(_quot(x,y),y)));
 
 	case t_POLMOD: case t_POL:
 	  return gzero;
 
-	default: err(operf,"%",tx,ty);
+	default: err(operf,"%",x,y);
       }
 
     case t_POL:
       if (is_scalar_t(tx))
       {
         if (tx!=t_POLMOD || varn(x[1]) > varn(y))
-          return (lgef(y) < 4)? gzero: gcopy(x);
+          return degpol(y)? gcopy(x): gzero;
 	if (varn(x[1])!=varn(y)) return gzero;
         z=cgetg(3,t_POLMOD);
         z[1]=lgcd((GEN)x[1],y);
@@ -663,8 +706,10 @@ gmod(GEN x, GEN y)
 	case t_POL:
 	  return gres(x,y);
 
-	case t_RFRAC: case t_RFRACN:
-	  av=avma; if (tx==t_RFRACN) x=gred_rfrac(x);
+        case t_RFRACN:
+	  av=avma; return gerepileupto(av, gmod(gred_rfrac(x), y));
+	case t_RFRAC:
+	  av=avma;
 	  p1=gmul((GEN)x[1],ginvmod((GEN)x[2],y)); tetpil=avma;
           return gerepile(av,tetpil,gres(p1,y));
 
@@ -672,14 +717,14 @@ gmod(GEN x, GEN y)
           if (ismonome(y) && varn(x) == varn(y))
           {
             long d = degpol(y);
-            if (lg(x)-2 + valp(x) < d) err(operi,"%",tx,ty);
+            if (lg(x)-2 + valp(x) < d) err(operi,"%",x,y);
             av = avma; 
             return gerepileupto(av, gmod(gtrunc(x), y));
           }
-	default: err(operf,"%",tx,ty);
+	default: err(operf,"%",x,y);
       }
   }
-  err(operf,"%",tx,ty);
+  err(operf,"%",x,y);
   return NULL; /* not reached */
 }
 
@@ -696,7 +741,7 @@ gmodulsg(long x, GEN y)
     case t_POL: z=cgetg(3,t_POLMOD);
       z[1]=lcopy(y); z[2]=lstoi(x); return z;
   }
-  err(operf,"%",t_INT,typ(y)); return NULL; /* not reached */
+  err(operf,"%",stoi(x),y); return NULL; /* not reached */
 }
 
 GEN
@@ -743,7 +788,7 @@ gmodulo(GEN x,GEN y)
       if (tx!=t_POL && !is_rfrac_t(tx) && tx!=t_SER) break;
       z[2]=(long)specialmod(x,y); return z;
   }
-  err(operf,"%",tx,typ(y)); return NULL; /* not reached */
+  err(operf,"%",x,y); return NULL; /* not reached */
 }
 
 GEN
@@ -776,7 +821,7 @@ gmodulcp(GEN x,GEN y)
       if (tx!=t_POL && !is_rfrac_t(tx) && tx!=t_SER) break;
       z[2]=(long)specialmod(x,y); return z;
   }
-  err(operf,"%",tx,typ(y)); return NULL; /* not reached */
+  err(operf,"%",x,y); return NULL; /* not reached */
 }
 
 GEN
@@ -793,52 +838,202 @@ Mod0(GEN x,GEN y,long flag)
 
 /*******************************************************************/
 /*                                                                 */
-/*                 DIVISION ENTIERE GENERALE                       */
-/*            DIVISION ENTIERE AVEC RESTE GENERALE                 */
+/*                 GENERIC EUCLIDEAN DIVISION                      */
 /*                                                                 */
 /*******************************************************************/
 
 GEN
 gdivent(GEN x, GEN y)
 {
-  long tx=typ(x),ty=typ(y);
+  long tx = typ(x);
 
-  if (tx == t_INT)
+  if (is_matvec_t(tx))
   {
-    if (ty == t_INT) return truedvmdii(x,y,NULL);
-    if (ty!=t_POL) err(typeer,"gdivent");
-    return gzero;
+    long i, lx = lg(x);
+    GEN z = cgetg(lx,tx);
+    for (i=1; i<lx; i++) z[i] = (long)gdivent((GEN)x[i],y);
+    return z;
   }
-  if (ty != t_POL) err(typeer,"gdivent");
-  if (tx == t_POL) return gdeuc(x,y);
-  if (! is_scalar_t(tx)) err(typeer,"gdivent");
-  return gzero;
+  switch(typ(y))
+  {
+    case t_INT:
+      switch(tx)
+      { /* equal to, but more efficient than, quot(x,y) */
+        case t_INT: return truedvmdii(x,y,NULL);
+        case t_REAL:
+        case t_FRAC:
+        case t_FRACN: return quot(x,y);
+        case t_POL: return gdiv(x,y);
+      }
+      break;
+    case t_REAL:
+    case t_FRAC:
+    case t_FRACN: return quot(x,y);
+    case t_POL:
+      if (is_scalar_t(tx))
+      {
+        if (tx == t_POLMOD) break;
+        return degpol(y)? gzero: gdiv(x,y);
+      }
+      if (tx == t_POL) return gdeuc(x,y);
+  }
+  err(operf,"\\",x,y);
+  return NULL; /* not reached */
 }
 
+/* with remainder */
+static GEN
+quotrem(GEN x, GEN y, GEN *r)
+{
+  gpmem_t av;
+  GEN q = quot(x,y);
+  av = avma;
+  *r = gerepileupto(av, gsub(x, gmul(q,y)));
+  return q;
+}
+
 GEN
 gdiventres(GEN x, GEN y)
 {
-  long tx=typ(x),ty=typ(y);
-  GEN z = cgetg(3,t_COL);
+  long tx = typ(x);
+  GEN z,q,r;
 
-  if (tx==t_INT)
+  if (is_matvec_t(tx))
   {
-    if (ty==t_INT)
+    long i, lx = lg(x);
+    z = cgetg(lx,tx);
+    for (i=1; i<lx; i++) z[i] = (long)gdiventres((GEN)x[i],y);
+    return z;
+  }
+  z = cgetg(3,t_COL);
+  switch(typ(y))
+  {
+    case t_INT:
+      switch(tx)
+      { /* equal to, but more efficient than next case */
+        case t_INT:
+          z[1] = (long)truedvmdii(x,y,(GEN*)(z+2));
+          return z;
+        case t_REAL:
+        case t_FRAC:
+        case t_FRACN:
+          q = quotrem(x,y,&r);
+          z[1] = (long)q;
+          z[2] = (long)r; return z;
+        case t_POL: return gdiv(x,y);
+      }
+      break;
+    case t_REAL:
+    case t_FRAC:
+    case t_FRACN:
+          q = quotrem(x,y,&r);
+          z[1] = (long)q;
+          z[2] = (long)r; return z;
+    case t_POL:
+      if (is_scalar_t(tx))
+      {
+        if (tx == t_POLMOD) break;
+        if (degpol(y))
+        {
+          q = gzero;
+          r = gcopy(x);
+        }
+        else
+        {
+          q = gdiv(x,y);
+          r = gzero;
+        }
+        return q;
+      }
+      if (tx == t_POL)
+      {
+        z[1]=ldivres(x,y,(GEN *)(z+2));
+        return z;
+      }
+  }
+  err(operf,"\\",x,y);
+  return NULL; /* not reached */
+}
+
+extern GEN swap_vars(GEN b0, long v);
+
+GEN
+divrem(GEN x, GEN y, long v)
+{
+  gpmem_t av = avma;
+  long vx,vy;
+  GEN z,q,r;
+  if (v < 0 || typ(y) != t_POL || typ(x) != t_POL) return gdiventres(x,y);
+  vx = varn(x); if (vx != v) x = swap_vars(x,v);
+  vy = varn(y); if (vy != v) y = swap_vars(y,v);
+  q = poldivres(x,y, &r);
+  if (v && (vx != v || vy != v))
+  {
+    q = poleval(q, polx[v]);
+    r = poleval(r, polx[v]);
+  }
+  z = cgetg(3,t_COL);
+  z[1] = (long)q;
+  z[2] = (long)r; return gerepilecopy(av, z);
+}
+
+static int
+is_scal(long t) { return t==t_INT || t==t_FRAC || t==t_FRACN; }
+
+GEN
+diviiround(GEN x, GEN y)
+{
+  gpmem_t av1, av = avma;
+  GEN q,r;
+  int fl;
+
+  q = dvmdii(x,y,&r); /* q = x/y rounded towards 0, sgn(r)=sgn(x) */
+  if (r==gzero) return q;
+  av1 = avma;
+  fl = absi_cmp(shifti(r,1),y);
+  avma = av1; cgiv(r);
+  if (fl >= 0) /* If 2*|r| >= |y| */
+  {
+    long sz = signe(x)*signe(y);
+    if (fl || sz > 0) q = gerepileuptoint(av, addis(q,sz));
+  }
+  return q;
+}
+
+/* If x and y are not both scalars, same as gdivent.
+ * Otherwise, compute the quotient x/y, rounded to the nearest integer
+ * (towards +oo in case of tie). */
+GEN
+gdivround(GEN x, GEN y)
+{
+  gpmem_t av1, av;
+  long tx=typ(x),ty=typ(y);
+  GEN q,r;
+  int fl;
+
+  if (tx==t_INT && ty==t_INT) return diviiround(x,y);
+  av = avma;
+  if (is_scal(tx) && is_scal(ty))
+  { /* same as diviiround but less efficient */
+    q = quotrem(x,y,&r);
+    av1 = avma;
+    fl = gcmp(gmul2n(gabs(r,0),1), gabs(y,0));
+    avma = av1; cgiv(r);
+    if (fl >= 0) /* If 2*|r| >= |y| */
     {
-      z[1]=(long)truedvmdii(x,y,(GEN*)(z+2));
-      return z;
+      long sz = gsigne(y);
+      if (fl || sz > 0) q = gerepileupto(av, gaddgs(q, sz));
     }
-    if (ty!=t_POL) err(typeer,"gdiventres");
-    z[1]=zero; z[2]=licopy(x); return z;
+    return q;
   }
-  if (ty != t_POL) err(typeer,"gdiventres");
-  if (tx == t_POL)
+  if (is_matvec_t(tx))
   {
-    z[1]=ldivres(x,y,(GEN *)(z+2));
+    long i,lx = lg(x);
+    GEN z = cgetg(lx,tx);
+    for (i=1; i<lx; i++) z[i] = (long)gdivround((GEN)x[i],y);
     return z;
   }
-  if (! is_scalar_t(tx)) err(typeer,"gdiventres");
-  z[1]=zero; z[2]=lcopy(x); return z;
+  return gdivent(x,y);
 }
 
 GEN
@@ -857,43 +1052,6 @@ gdivmod(GEN x, GEN y, GEN *pr)
   return poldivres(x,y,pr);
 }
 
-/* When x and y are integers, compute the quotient x/y, rounded to the
- * nearest integer. If there is a tie, the quotient is rounded towards zero.
- * If x and y are not both integers, same as gdivent.
- */
-GEN
-gdivround(GEN x, GEN y)
-{
-  long tx=typ(x),ty=typ(y);
-
-  if (tx==t_INT)
-  {
-    if (ty==t_INT)
-    {
-      long av = avma, av1,fl;
-      GEN r, q=dvmdii(x,y,&r); /* q = x/y rounded towards 0, sqn(r)=sgn(x) */
-
-      if (r==gzero) return q;
-      av1 = avma;
-      fl = absi_cmp(shifti(r,1),y);
-      avma = av1; cgiv(r);
-      if (fl >= 0) /* If 2*|r| >= |y| */
-      {
-        long sz = signe(x)*signe(y);
-	if (fl || sz > 0)
-          { av1=avma; q = gerepile(av,av1,addis(q,sz)); }
-      }
-      return q;
-    }
-    if (ty!=t_POL) err(typeer,"gdivround");
-    return gzero;
-  }
-  if (ty != t_POL) err(typeer,"gdivround");
-  if (tx == t_POL) return gdeuc(x,y);
-  if (! is_scalar_t(tx)) err(typeer,"gdivround");
-  return gzero;
-}
-
 /*******************************************************************/
 /*                                                                 */
 /*                       SHIFT D'UN GEN                            */
@@ -904,13 +1062,19 @@ gdivround(GEN x, GEN y)
 GEN
 gshift(GEN x, long n)
 {
+  return gshift3(x, n, 0);
+}
+
+GEN
+gshift3(GEN x, long n, long flag)
+{
   long i,l,lx, tx = typ(x);
   GEN y;
 
   switch(tx)
   {
     case t_INT:
-      return shifti(x,n);
+      return shifti3(x,n,flag);
     case t_REAL:
       return shiftr(x,n);
 
@@ -929,7 +1093,8 @@ extern GEN mulscalrfrac(GEN x, GEN y);
 GEN
 gmul2n(GEN x, long n)
 {
-  long tx,lx,i,k,l,av,tetpil;
+  long tx, lx, i, k, l;
+  gpmem_t av, tetpil;
   GEN p2,p1,y;
 
   tx=typ(x);
@@ -951,12 +1116,12 @@ gmul2n(GEN x, long n)
       if (n > 0)
       {
         y=cgetg(3,t_INTMOD); p2=(GEN)x[1];
-        av=avma; new_chunk(2 * lgefint(p2) + n); /* HACK */
+        av=avma; (void)new_chunk(2 * lgefint(p2) + n); /* HACK */
         p1 = shifti((GEN)x[2],n); avma=av;
         y[2]=lmodii(p1,p2); icopyifstack(p2,y[1]); return y;
       }
-      l=avma; y=gmul2n(gun,n); tetpil=avma;
-      return gerepile(l,tetpil,gmul(y,x));
+      av=avma; y=gmul2n(gun,n); tetpil=avma;
+      return gerepile(av,tetpil,gmul(y,x));
 
     case t_FRAC: case t_FRACN:
       l = vali((GEN)x[1]);
@@ -1007,13 +1172,16 @@ gmul2n(GEN x, long n)
       return y;
 
     case t_PADIC:
-      l=avma; y=gmul2n(gun,n); tetpil=avma;
-      return gerepile(l,tetpil,gmul(y,x));
+      av=avma; y=gmul2n(gun,n); tetpil=avma;
+      return gerepile(av,tetpil,gmul(y,x));
   }
   err(typeer,"gmul2n");
   return NULL; /* not reached */
 }
 
+GEN
+real2n(long n, long prec) { GEN z = realun(prec); setexpo(z, n); return z; }
+
 /*******************************************************************/
 /*                                                                 */
 /*                      INVERSE D'UN GEN                           */
@@ -1024,7 +1192,8 @@ extern GEN fix_rfrac_if_pol(GEN x, GEN y);
 GEN
 ginv(GEN x)
 {
-  long tx=typ(x),av,tetpil,s;
+  long tx=typ(x), s;
+  gpmem_t av, tetpil;
   GEN z,y,p1,p2;
 
   switch(tx)
@@ -1098,6 +1267,18 @@ ginv(GEN x)
       return y;
     case t_MAT:
       return (lg(x)==1)? cgetg(1,t_MAT): invmat(x);
+    case t_VECSMALL:
+    {
+      long i,lx = lg(x);
+      y = cgetg(lx,t_VECSMALL);
+      for (i=1; i<lx; i++)
+      {
+        long xi=x[i];
+	  if (xi<1 || xi>=lx) err(talker,"incorrect permtuation to inverse");
+        y[xi] = i;
+      }
+      return y;
+    }
   }
   err(typeer,"inverse");
   return NULL; /* not reached */
@@ -1113,14 +1294,15 @@ ginv(GEN x)
 static GEN
 gconvsp(GEN x, int flpile)
 {
-  long v = varn(x), av,tetpil,i;
+  long v = varn(x), i;
+  gpmem_t av, tetpil;
   GEN p1,y;
 
   if (gcmp0(x)) return zeropol(v);
   av=avma; y=dummycopy(x); settyp(y,t_POL);
   i=lg(x)-1; while (i>1 && gcmp0((GEN)y[i])) i--;
   setlgef(y,i+1);
-  p1=gpuigs(polx[v],valp(x));
+  p1=gpowgs(polx[v],valp(x));
   tetpil=avma; p1=gmul(p1,y);
   return flpile? gerepile(av,tetpil,p1): p1;
 }
@@ -1128,7 +1310,7 @@ gconvsp(GEN x, int flpile)
 GEN
 gsubst0(GEN x, GEN T, GEN y)
 {
-  ulong av;
+  gpmem_t av;
   long d, v;
   if (typ(T) != t_POL || !ismonome(T) || !gcmp1(leading_term(T)))
     err(talker,"variable number expected in subst");
@@ -1142,7 +1324,8 @@ GEN
 gsubst(GEN x, long v, GEN y)
 {
   long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
-  long l,vx,vy,e,ex,ey,tetpil,av,i,j,k,jb,limite;
+  long l, vx, vy, e, ex, ey, i, j, k, jb;
+  gpmem_t tetpil, av, limite;
   GEN t,p1,p2,z;
 
   if (ty==t_MAT)
@@ -1178,21 +1361,29 @@ gsubst(GEN x, long v, GEN y)
   switch(tx)
   {
     case t_POL:
-      l=lgef(x);
+      l = lgef(x);
       if (l==2)
         return (ty==t_MAT)? gscalmat(gzero,ly-1): gzero;
 
-      vx=varn(x);
-      if (vx<v)
+      vx = varn(x);
+      if (vx < v)
       {
+        if (gvar(y) > vx)
+        { /* easy special case */
+          z = cgetg(l, t_POL); z[1] = x[1];
+          for (i=2; i<l; i++) z[i] = lsubst((GEN)x[i],v,y);
+          return normalizepol_i(z,l);
+        }
+        /* general case */
 	av=avma; p1=polx[vx]; z= gsubst((GEN)x[l-1],v,y);
 	for (i=l-1; i>2; i--) z=gadd(gmul(z,p1),gsubst((GEN)x[i-1],v,y));
 	return gerepileupto(av,z);
       }
+      /* v <= vx */
       if (ty!=t_MAT)
-        return (vx>v)? gcopy(x): poleval(x,y);
+        return (vx > v)? gcopy(x): poleval(x,y);
 
-      if (vx>v) return gscalmat(x,ly-1);
+      if (vx > v) return gscalmat(x,ly-1);
       if (l==3) return gscalmat((GEN)x[2],ly-1);
       av=avma; z=(GEN)x[l-1];
       for (i=l-1; i>2; i--) z=gaddmat((GEN)x[i-1],gmul(z,y));
@@ -1211,7 +1402,7 @@ gsubst(GEN x, long v, GEN y)
         p2=gsubst(p1,v,y); tetpil=avma; z=tayl(p2,vx,lx-2);
         if (ex)
         {
-          p1=gpuigs(polx[vx],ex); tetpil=avma; z=gmul(z,p1);
+          p1=gpowgs(polx[vx],ex); tetpil=avma; z=gmul(z,p1);
         }
         return gerepile(av,tetpil,z);
       }
@@ -1232,7 +1423,7 @@ gsubst(GEN x, long v, GEN y)
 	    av=avma; z = zeroser(vy,0);
 	    for (i=lx-1; i>=2; i--)
               z = gadd((GEN)x[i], gmul(y,z));
-	    if (ex) z = gmul(z, gpuigs(y,ex));
+	    if (ex) z = gmul(z, gpowgs(y,ex));
 	    return gerepileupto(av,z);
 	  }
 
@@ -1265,8 +1456,8 @@ gsubst(GEN x, long v, GEN y)
 	  }
 	  if (!ex) return gerepilecopy(av,z);
 
-          if (l<ly) { setlg(y,l); p1=gpuigs(y,ex); setlg(y,ly); }
-          else p1=gpuigs(y,ex);
+          if (l<ly) { setlg(y,l); p1=gpowgs(y,ex); setlg(y,ly); }
+          else p1=gpowgs(y,ex);
           tetpil=avma; return gerepile(av,tetpil,gmul(z,p1));
 
         case t_POL: case t_RFRAC: case t_RFRACN:
@@ -1289,8 +1480,9 @@ gsubst(GEN x, long v, GEN y)
 
     case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
       for (i=1; i<lx; i++) z[i]=lsubst((GEN)x[i],v,y);
+      return z;
   }
-  return z;
+  return gcopy(x);
 }
 
 /*******************************************************************/
@@ -1302,37 +1494,48 @@ gsubst(GEN x, long v, GEN y)
 GEN
 recip(GEN x)
 {
-  long tetpil, av=avma, v=varn(x);
+  long v=varn(x), lx = lg(x);
+  gpmem_t tetpil, av=avma;
   GEN p1,p2,a,y,u;
 
   if (typ(x)!=t_SER) err(talker,"not a series in serreverse");
-  if (valp(x)!=1) err(talker,"valuation not equal to 1 in serreverse");
+  if (valp(x)!=1 || lx < 3)
+    err(talker,"valuation not equal to 1 in serreverse");
 
   a=(GEN)x[2];
   if (gcmp1(a))
   {
-    long i,j,k, lx=lg(x), lim=stack_lim(av,2);
+    long i, j, k, mi;
+    gpmem_t lim=stack_lim(av, 2);
 
-    u=cgetg(lx,t_SER); y=cgetg(lx,t_SER);
-    u[1]=y[1]=evalsigne(1) | evalvalp(1) | evalvarn(v);
-    u[2]=un; u[3]=lmulsg(-2,(GEN)x[3]);
-    y[2]=un; y[3]=lneg((GEN)x[3]);
+    mi = lx-1; while (mi>=3 && gcmp0((GEN)x[mi])) mi--;
+    u = cgetg(lx,t_SER);
+    y = cgetg(lx,t_SER);
+    u[1] = y[1] = evalsigne(1) | evalvalp(1) | evalvarn(v);
+    u[2] = y[2] = un;
+    if (lx > 3)
+    {
+      u[3] = lmulsg(-2,(GEN)x[3]);
+      y[3] = lneg((GEN)x[3]);
+    }
     for (i=3; i<lx-1; )
     {
       for (j=3; j<i+1; j++)
       {
-        p1=(GEN)u[j];
-        for (k=j-1; k>2; k--)
-          p1=gsub(p1,gmul((GEN)u[k],(GEN)x[j-k+2]));
-        u[j]=lsub(p1,(GEN)x[j]);
+        p1 = (GEN)x[j];
+        for (k=max(3,j+2-mi); k<j; k++)
+          p1 = gadd(p1, gmul((GEN)u[k],(GEN)x[j-k+2]));
+        u[j] = lsub((GEN)u[j], p1);
       }
-      p1=gmulsg(i,(GEN)x[i+1]);
-      for (k=2; k<i; k++)
+      p1 = gmulsg(i,(GEN)x[i+1]);
+      for (k=2; k<min(i,mi); k++)
       {
-        p2=gmul((GEN)x[k+1],(GEN)u[i-k+2]);
-        p1=gadd(p1,gmulsg(k,p2));
+        p2 = gmul((GEN)x[k+1],(GEN)u[i-k+2]);
+        p1 = gadd(p1, gmulsg(k,p2));
       }
-      i++; u[i]=lneg(p1); y[i]=ldivgs((GEN)u[i],i-1);
+      i++;
+      u[i] = lneg(p1);
+      y[i] = ldivgs((GEN)u[i],i-1);
       if (low_stack(lim, stack_lim(av,2)))
       {
 	GEN *gptr[2];
@@ -1343,9 +1546,9 @@ recip(GEN x)
     }
     return gerepilecopy(av,y);
   }
-  y=gdiv(x,a); y[2]=un; y=recip(y);
-  a=gdiv(polx[v],a); tetpil=avma;
-  return gerepile(av,tetpil,gsubst(y,v,a));
+  y = gdiv(x,a); y[2] = un; y = recip(y);
+  a = gdiv(polx[v],a); tetpil = avma;
+  return gerepile(av,tetpil, gsubst(y,v,a));
 }
 
 /*******************************************************************/
@@ -1359,7 +1562,7 @@ derivpol(GEN x)
   long i,lx = lgef(x)-1;
   GEN y;
 
-  if (lx<3) return gzero;
+  if (lx<3) return zeropol(varn(x));
   y = cgetg(lx,t_POL);
   for (i=2; i<lx ; i++) y[i] = lmulsg(i-1,(GEN)x[i+1]);
   y[1] = x[1]; return normalizepol_i(y,i);
@@ -1370,7 +1573,7 @@ derivser(GEN x)
 {
   long i,j,vx = varn(x), e = valp(x), lx = lg(x);
   GEN y;
-  if (gcmp0(x)) return zeroser(vx,e-1);
+  if (gcmp0(x)) return zeroser(vx,e? e-1: 0);
   if (e)
   {
     y=cgetg(lx,t_SER); y[1] = evalsigne(1) | evalvalp(e-1) | evalvarn(vx);
@@ -1389,7 +1592,8 @@ derivser(GEN x)
 GEN
 deriv(GEN x, long v)
 {
-  long lx,vx,tx,e,i,j,l,av,tetpil;
+  long lx, vx, tx, e, i, j;
+  gpmem_t av, tetpil;
   GEN y,p1,p2;
 
   tx=typ(x); if (is_const_t(tx)) return gzero;
@@ -1419,11 +1623,11 @@ deriv(GEN x, long v)
       {
         if (!signe(x)) return gcopy(x);
         lx=lg(x); e=valp(x);
-	l=avma;
+	av=avma;
 	for (i=2; i<lx; i++)
         {
           if (!gcmp0(deriv((GEN)x[i],v))) break;
-          avma=l;
+          avma=av;
         }
 	if (i==lx) return ggrando(polx[vx],e+lx-2);
 	y=cgetg(lx-i+2,t_SER);
@@ -1434,11 +1638,11 @@ deriv(GEN x, long v)
       return derivser(x);
 
     case t_RFRAC: case t_RFRACN: av=avma; y=cgetg(3,tx);
-      y[2]=lsqr((GEN)x[2]); l=avma;
+      y[2]=lsqr((GEN)x[2]); av=avma;
       p1=gmul((GEN)x[2],deriv((GEN)x[1],v));
       p2=gmul(gneg_i((GEN)x[1]),deriv((GEN)x[2],v));
       tetpil=avma; p1=gadd(p1,p2);
-      if (tx==t_RFRACN) { y[1]=lpile(l,tetpil,p1); return y; }
+      if (tx==t_RFRACN) { y[1]=lpile(av,tetpil,p1); return y; }
       y[1]=(long)p1; return gerepileupto(av,gred_rfrac(y));
 
     case t_VEC: case t_COL: case t_MAT: lx=lg(x); y=cgetg(lx,tx);
@@ -1469,7 +1673,8 @@ triv_integ(GEN x, long v, long tx, long lx)
 GEN
 integ(GEN x, long v)
 {
-  long lx,tx,e,i,j,vx,n,av=avma,tetpil;
+  long lx, tx, e, i, j, vx, n;
+  gpmem_t av=avma, tetpil;
   GEN y,p1;
 
   tx = typ(x);
@@ -1614,7 +1819,7 @@ gfloor(GEN x)
 GEN
 gfrac(GEN x)
 {
-  long av = avma,tetpil;
+  gpmem_t av = avma, tetpil;
   GEN p1 = gneg_i(gfloor(x));
 
   tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
@@ -1624,7 +1829,8 @@ GEN
 gceil(GEN x)
 {
   GEN y,p1;
-  long i,lx,tx=typ(x),av,tetpil;
+  long i, lx, tx=typ(x);
+  gpmem_t av, tetpil;
 
   switch(tx)
   {
@@ -1671,7 +1877,8 @@ GEN
 ground(GEN x)
 {
   GEN y,p1;
-  long i,lx,tx=typ(x),av,tetpil;
+  long i, lx, tx=typ(x);
+  gpmem_t av;
 
   switch(tx)
   {
@@ -1683,14 +1890,13 @@ ground(GEN x)
       long ex, s = signe(x);
       if (s==0 || (ex=expo(x)) < -1) return gzero;
       if (ex < 0) return s>0? gun: negi(gun);
-      av=avma; p1 = realun(3 + (ex>>TWOPOTBITS_IN_LONG));
-      setexpo(p1,-1); /* p1 = 0.5 */
-      p1 = addrr(x,p1); tetpil = avma;
-      return gerepile(av,tetpil,mpent(p1));
+      av = avma;
+      p1 = real2n(-1, 3 + (ex>>TWOPOTBITS_IN_LONG)); /* = 0.5 */
+      p1 = addrr(x,p1);
+      return gerepileuptoint(av, mpent(p1));
     }
     case t_FRAC: case t_FRACN:
-      av=avma; p1 = addii(shifti((GEN)x[2], -1), (GEN)x[1]);
-      return gerepileuptoint(av, truedvmdii(p1, (GEN)x[2], NULL));
+      return diviiround((GEN)x[1], (GEN)x[2]);
 
     case t_POLMOD: y=cgetg(3,t_POLMOD);
       copyifstack(x[1],y[1]);
@@ -1715,7 +1921,8 @@ GEN
 grndtoi(GEN x, long *e)
 {
   GEN y,p1;
-  long i,tx=typ(x), lx,av,ex,e1;
+  long i, tx=typ(x), lx, ex, e1;
+  gpmem_t av;
 
   *e = -HIGHEXPOBIT;
   switch(tx)
@@ -1734,7 +1941,7 @@ grndtoi(GEN x, long *e)
       lx=lg(x); e1 = ex - bit_accuracy(lx) + 1;
       settyp(p1,t_INT); setlgefint(p1,lx);
       y=shifti(p1,e1); if (signe(x)<0) y=addsi(-1,y);
-      y = gerepileupto(av,y);
+      y = gerepileuptoint(av,y);
 
       if (e1<=0) { av=avma; e1=expo(subri(x,y)); avma=av; }
       *e=e1; return y;
@@ -1763,7 +1970,8 @@ grndtoi(GEN x, long *e)
 GEN
 gcvtoi(GEN x, long *e)
 {
-  long tx=typ(x), lx,i,ex,av,e1;
+  long tx=typ(x), lx, i, ex, e1;
+  gpmem_t av;
   GEN y;
 
   *e = -HIGHEXPOBIT;
@@ -1795,7 +2003,7 @@ gcvtoi(GEN x, long *e)
 GEN
 ceil_safe(GEN x)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long e;
   GEN y = gcvtoi(x,&e);
   if (e < 0) e = 0;
@@ -1806,7 +2014,8 @@ ceil_safe(GEN x)
 GEN
 gtrunc(GEN x)
 {
-  long tx=typ(x),av,tetpil,i,v;
+  long tx=typ(x), i, v;
+  gpmem_t av;
   GEN y;
 
   switch(tx)
@@ -1826,12 +2035,12 @@ gtrunc(GEN x)
       if (!v) return gcopy((GEN)x[4]);
       if (v>0)
       { /* here p^v is an integer */
-        av=avma; y=gpuigs((GEN)x[2],v); tetpil=avma;
-        return gerepile(av,tetpil, mulii(y,(GEN)x[4]));
+        av = avma; y = gpowgs((GEN)x[2],v);
+        return gerepileuptoint(av, mulii(y,(GEN)x[4]));
       }
       y=cgetg(3,t_FRAC);
-      y[1]=licopy((GEN)x[4]);
-      y[2]=lpuigs((GEN)x[2],-v);
+      y[1] = licopy((GEN)x[4]);
+      y[2] = lpowgs((GEN)x[2],-v);
       return y;
 
     case t_RFRAC: case t_RFRACN:
@@ -2032,7 +2241,8 @@ scalarser(GEN x, long v, long prec)
 GEN
 gtoser(GEN x, long v)
 {
-  long tx=typ(x),lx,i,j,l,av,tetpil;
+  long tx=typ(x), lx, i, j, l;
+  gpmem_t av, tetpil;
   GEN y,p1,p2;
 
   if (v<0) v = 0;
@@ -2112,12 +2322,7 @@ gtovec(GEN x)
     for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
     return y;
   }
-  if (tx == t_VECSMALL)
-  {
-    lx=lg(x); y=cgetg(lx,t_VEC);
-    for (i=1; i<lx; i++) y[i]=lstoi(x[i]);
-    return y;
-  }
+  if (tx == t_VECSMALL) return small_to_vec(x);
   if (!signe(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
   lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
   for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
@@ -2160,74 +2365,104 @@ compo(GEN x, long n)
   return gcopy((GEN)x[l]);
 }
 
+/* assume v > varn(x), extract coeff of polx[v]^n */
+static GEN
+multi_coeff(GEN x, long n, long v, long dx)
+{
+  long i;
+  GEN z;
+  if (dx == 0) return polcoeff_i((GEN)x[2], n, v);
+  z = cgetg(dx+3, t_POL);
+  z[1] = x[1];
+  for (i=2; i<dx+3; i++) z[i] = (long)polcoeff_i((GEN)x[i], n, v);
+  return normalizepol(z);
+}
+
+/* assume x a t_POL */
+static GEN
+_polcoeff(GEN x, long n, long v)
+{
+  long w, dx;
+  dx = degpol(x);
+  if (dx < 0) return gzero;
+  if (v < 0 || v == (w=varn(x)))
+    return (n < 0 || n > dx)? gzero: (GEN)x[n+2];
+  if (w > v) return n? gzero: x;
+  /* w < v */
+  return multi_coeff(x, n, v, dx);
+}
+
+/* assume x a t_SER */
+static GEN
+_sercoeff(GEN x, long n, long v)
+{
+  long w, dx, ex;
+  if (!signe(x)) return gzero;
+  dx = lg(x)-3; ex = valp(x); n -= ex;
+  if (v < 0 || v == (w=varn(x)))
+  {
+    if (n > dx) err(talker,"non existent component in truecoeff");
+    return (n < 0)? gzero: (GEN)x[n+2];
+  }
+  if (w > v) return n? gzero: x;
+  /* w < v */
+  return multi_coeff(x, n, v, dx);
+}
+
+/* assume x a t_RFRAC(n) */
+static GEN
+_rfraccoeff(GEN x, long n, long v)
+{
+  GEN P,Q, p = (GEN)x[1], q = (GEN)x[2];
+  long vp = gvar(p), vq = gvar(q);
+  if (v < 0) v = min(vp, vq);
+  P = (vp == v)? p: swap_vars(p, v);
+  Q = (vq == v)? q: swap_vars(q, v);
+  if (!ismonome(Q)) err(typeer, "polcoeff");
+  n += degpol(Q);
+  return gdiv(_polcoeff(P, n, v), leading_term(Q));
+}
+
+GEN
+polcoeff_i(GEN x, long n, long v)
+{
+  switch(typ(x))
+  {
+    case t_POL: return _polcoeff(x,n,v);
+    case t_SER: return _sercoeff(x,n,v);
+    case t_RFRACN: err(typeer, "polcoeff"); /* too expensive to reduce */
+    case t_RFRAC: return _rfraccoeff(x,n,v);
+    default: return n? gzero: x;
+  }
+}
+
 /* with respect to the main variable if v<0, with respect to the variable v
    otherwise. v ignored if x is not a polynomial/series. */
-
 GEN
 polcoeff0(GEN x, long n, long v)
 {
-  long tx=typ(x),lx,ex,w,av,tetpil;
-  GEN xinit;
+  long tx=typ(x);
+  gpmem_t av;
 
   if (is_scalar_t(tx)) return n? gzero: gcopy(x);
 
+  av = avma;
   switch(tx)
   {
+    case t_POL: x = _polcoeff(x,n,v); break;
+    case t_SER: x = _sercoeff(x,n,v); break;
+    case t_RFRAC:
+    case t_RFRACN: x = _rfraccoeff(x,n,v); break;
+   
     case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
-      if (n<1 || n>=lg(x)) break;
-      return gcopy((GEN)x[n]);
+      if (n>=1 && n<lg(x)) return gcopy((GEN)x[n]);
+    /* fall through */
 
-    case t_POL:
-      if (n<0) return gzero;
-      w=varn(x);
-      if (v < 0 || v == w)
-	return (n>=lgef(x)-2)? gzero: gcopy((GEN)x[n+2]);
-      if (v < w) return n? gzero: gcopy(x);
-      av=avma; xinit=x;
-      x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
-      if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
-      if (typ(x) == t_POL)
-      {
-        if (n>=lgef(x)-2) { avma=av; return gzero; }
-        x = (GEN)x[n+2];
-      }
-      else
-        x = polcoeff0(x, n, 0);
-      tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
-
-    case t_SER:
-      w=varn(x);
-      if (v < 0 || v == w)
-      {
-	if (!signe(x)) return gzero;
-	lx=lg(x); ex=valp(x); if (n<ex) return gzero;
-	if (n>=ex+lx-2) break;
-	return gcopy((GEN)x[n-ex+2]);
-      }
-      if (v < w) return n?  gzero: gcopy(x);
-      av=avma; xinit=x;
-      x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
-      if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
-      if (gcmp0(x)) { avma=av; return gzero; }
-      if (typ(x) == t_SER)
-      {
-        lx=lg(x); ex=valp(x); if (n<ex) return gzero;
-        if (n>=ex+lx-2) break;
-        x = (GEN)x[n-ex+2];
-      }
-      else
-        x = polcoeff0(x, n, 0);
-      tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
-
-    case t_RFRAC: case t_RFRACN:
-      w = precdl; av = avma;
-      if (v<0) v = gvar(x);
-      ex = ggval((GEN)x[2], polx[v]);
-      precdl = n + ex + 1; x = gtoser(x, v); precdl = w;
-      return gerepileupto(av, polcoeff0(x, n, v));
+    default: err(talker,"nonexistent component in truecoeff");
   }
-  err(talker,"nonexistent component in truecoeff");
-  return NULL; /* not reached */
+  if (x == gzero) return x;
+  if (avma == av) return gcopy(x);
+  return gerepilecopy(av, x);
 }
 
 GEN
@@ -2239,7 +2474,8 @@ truecoeff(GEN x, long n)
 GEN
 denom(GEN x)
 {
-  long lx,i,av,tetpil;
+  long lx, i;
+  gpmem_t av, tetpil;
   GEN s,t;
 
   switch(typ(x))
@@ -2285,7 +2521,7 @@ denom(GEN x)
 GEN
 numer(GEN x)
 {
-  long av,tetpil;
+  gpmem_t av, tetpil;
   GEN s;
 
   switch(typ(x))
@@ -2407,7 +2643,8 @@ lift_intern0(GEN x, long v)
 GEN
 centerlift0(GEN x, long v)
 {
-  long lx,tx=typ(x),i,av;
+  long lx, tx=typ(x), i;
+  gpmem_t av;
   GEN y;
 
   switch(tx)
@@ -2463,7 +2700,8 @@ centerlift(GEN x)
 static GEN
 op_ReIm(GEN f(GEN), GEN x)
 {
-  long lx,i,j,av, tx = typ(x);
+  long lx, i, j, tx = typ(x);
+  gpmem_t av;
   GEN z;
 
   switch(tx)
@@ -2549,7 +2787,7 @@ gimag(GEN x)
 static long
 _egal(GEN x, GEN y)
 {
-  long av = avma;
+  gpmem_t av = avma;
   long r = gegal(simplify_i(x), simplify_i(y));
   avma = av; return r;
 }
@@ -2590,7 +2828,8 @@ gnot(GEN x) { return gcmp0(x)? gun: gzero; }
 GEN
 geval(GEN x)
 {
-  long av,tetpil,lx,i, tx = typ(x);
+  long lx, i, tx = typ(x);
+  gpmem_t av, tetpil;
   GEN y,z;
 
   if (is_const_t(tx)) return gcopy(x);
@@ -2710,7 +2949,7 @@ simplify_i(GEN x)
 GEN
 simplify(GEN x)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   return gerepilecopy(av, simplify_i(x));
 }
 
@@ -2722,7 +2961,8 @@ simplify(GEN x)
 static GEN
 qfeval0_i(GEN q, GEN x, long n)
 {
-  long i,j,av=avma;
+  long i, j;
+  gpmem_t av=avma;
   GEN res=gzero;
 
   for (i=2;i<n;i++)
@@ -2738,7 +2978,8 @@ qfeval0_i(GEN q, GEN x, long n)
 static GEN
 qfeval0(GEN q, GEN x, long n)
 {
-  long i,j,av=avma;
+  long i, j;
+  gpmem_t av=avma;
   GEN res=gzero;
 
   for (i=2;i<n;i++)
@@ -2753,7 +2994,8 @@ qfeval0(GEN q, GEN x, long n)
 static GEN
 qfeval0(GEN q, GEN x, long n)
 {
-  long i,j,av=avma;
+  long i, j;
+  gpmem_t av=avma;
   GEN res = gmul(gcoeff(q,1,1), gsqr((GEN)x[1]));
 
   for (i=2; i<n; i++)
@@ -2796,7 +3038,8 @@ qfeval(GEN q, GEN x)
 static GEN
 qfbeval0_i(GEN q, GEN x, GEN y, long n)
 {
-  long i,j,av=avma;
+  long i, j;
+  gpmem_t av=avma;
   GEN res = gmul(gcoeff(q,1,1), mulii((GEN)x[1],(GEN)y[1]));
 
   for (i=2;i<n;i++)
@@ -2815,7 +3058,8 @@ qfbeval0_i(GEN q, GEN x, GEN y, long n)
 static GEN
 qfbeval0(GEN q, GEN x, GEN y, long n)
 {
-  long i,j,av=avma;
+  long i, j;
+  gpmem_t av=avma;
   GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1],(GEN)y[1]));
 
   for (i=2;i<n;i++)
@@ -2833,7 +3077,8 @@ qfbeval0(GEN q, GEN x, GEN y, long n)
 static GEN
 qfbeval0(GEN q, GEN x, GEN y, long n)
 {
-  long i,j,av=avma;
+  long i, j;
+  gpmem_t av=avma;
   GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1], (GEN)y[1]));
 
   for (i=2; i<n; i++)
@@ -2908,44 +3153,6 @@ qf_base_change(GEN q, GEN M, int flag)
   return res;
 }
 
-/* compute M'.M */
-GEN
-gram_matrix(GEN M)
-{
-  long n=lg(M),i,j,k, av;
-  GEN res = cgetg(n,t_MAT),p1;
-
-  if (n==1)
-  {
-    if (typ(M) != t_MAT)
-      err(talker,"invalid data in gram_matrix");
-    return res;
-  }
-  if (typ(M) != t_MAT || lg(M[1]) != n)
-    err(talker,"not a square matrix in gram_matrix");
-
-  for (i=1;i<n;i++) res[i]=lgetg(n,t_COL);
-  av=avma;
-  for (i=1;i<n;i++)
-  {
-    p1 = gzero;
-    for (k=1;k<n;k++)
-      p1 = gadd(p1, gsqr(gcoeff(M,k,i)));
-    coeff(res,i,i) = (long) gerepileupto(av,p1);
-    av=avma;
-  }
-  for (i=2;i<n;i++)
-    for (j=1;j<i;j++)
-    {
-      p1=gzero;
-      for (k=1;k<n;k++)
-	p1 = gadd(p1, gmul(gcoeff(M,k,i),gcoeff(M,k,j)));
-      coeff(res,i,j)=coeff(res,j,i) = lpileupto(av,p1);
-      av=avma;
-    }
-  return res;
-}
-
 /* return Re(x * y), x and y scalars */
 GEN
 mul_real(GEN x, GEN y)
@@ -2954,7 +3161,7 @@ mul_real(GEN x, GEN y)
   {
     if (typ(y) == t_COMPLEX)
     {
-      long av=avma, tetpil;
+      gpmem_t av=avma, tetpil;
       GEN p1 = gmul((GEN)x[1], (GEN) y[1]);
       GEN p2 = gneg(gmul((GEN)x[2], (GEN) y[2]));
       tetpil=avma; return gerepile(av,tetpil,gadd(p1,p2));
@@ -2970,7 +3177,8 @@ mul_real(GEN x, GEN y)
 GEN
 mulmat_real(GEN x, GEN y)
 {
-  long av,i,j,k, lx = lg(x), ly = lg(y), l = lg(x[1]);
+  long i, j, k, lx = lg(x), ly = lg(y), l = lg(x[1]);
+  gpmem_t av;
   GEN p1, z = cgetg(ly,t_MAT);
 
   for (j=1; j<ly; j++)
@@ -2990,7 +3198,8 @@ mulmat_real(GEN x, GEN y)
 static GEN
 hqfeval0(GEN q, GEN x, long n)
 {
-  long i,j, av=avma;
+  long i, j;
+  gpmem_t av=avma;
   GEN res=gzero;
 
   for (i=2;i<n;i++)
@@ -3028,29 +3237,31 @@ hqfeval(GEN q, GEN x)
 GEN
 poleval(GEN x, GEN y)
 {
-  long av,tetpil,i,j,imin,tx=typ(x);
-  GEN p1,p2,p3,r,s;
+  long i, j, imin, tx = typ(x);
+  gpmem_t av0 = avma, av, lim;
+  GEN p1, p2, r, s;
 
   if (is_scalar_t(tx)) return gcopy(x);
   switch(tx)
   {
     case t_POL:
-      i=lgef(x)-1; imin=2; break;
+      i = lgef(x)-1; imin = 2; break;
 
-    case t_RFRAC: case t_RFRACN: av=avma;
-      p1=poleval((GEN)x[1],y);
-      p2=poleval((GEN)x[2],y); tetpil=avma;
-      return gerepile(av,tetpil,gdiv(p1,p2));
+    case t_RFRAC: case t_RFRACN:
+      p1 = poleval((GEN)x[1],y);
+      p2 = poleval((GEN)x[2],y);
+      return gerepileupto(av0, gdiv(p1,p2));
 
     case t_VEC: case t_COL:
-      i=lg(x)-1; imin=1; break;
+      i = lg(x)-1; imin = 1; break;
     default: err(typeer,"poleval");
       return NULL; /* not reached */
   }
   if (i<=imin)
     return (i==imin)? gcopy((GEN)x[imin]): gzero;
 
-  av=avma; p1=(GEN)x[i]; i--;
+  lim = stack_lim(av0,2);
+  p1 = (GEN)x[i]; i--;
   if (typ(y)!=t_COMPLEX)
   {
 #if 0 /* standard Horner's rule */
@@ -3063,21 +3274,31 @@ poleval(GEN x, GEN y)
       for (j=i; gcmp0((GEN)x[j]); j--)
         if (j==imin)
         {
-          if (i!=j) y = gpuigs(y,i-j+1);
-          tetpil=avma; return gerepile(av,tetpil,gmul(p1,y));
+          if (i!=j) y = gpowgs(y, i-j+1);
+          return gerepileupto(av0, gmul(p1,y));
         }
-      r = (i==j)? y: gpuigs(y,i-j+1);
+      r = (i==j)? y: gpowgs(y, i-j+1);
       p1 = gadd(gmul(p1,r), (GEN)x[j]);
+      if (low_stack(lim, stack_lim(av0,2)))
+      {
+        if (DEBUGMEM>1) err(warnmem,"poleval: i = %ld",i);
+        p1 = gerepileupto(av0, p1);
+      }
     }
-    return gerepileupto(av,p1);
+    return gerepileupto(av0,p1);
   }
 
-  p2=(GEN)x[i]; i--; r=gtrace(y); s=gneg_i(gnorm(y));
+  p2 = (GEN)x[i]; i--; r = gtrace(y); s = gneg_i(gnorm(y));
+  av = avma;
   for ( ; i>=imin; i--)
   {
-    p3=gadd(p2,gmul(r,p1));
-    p2=gadd((GEN)x[i],gmul(s,p1)); p1=p3;
+    GEN p3 = gadd(p2, gmul(r, p1));
+    p2 = gadd((GEN)x[i], gmul(s, p1)); p1 = p3;
+    if (low_stack(lim, stack_lim(av0,2)))
+    {
+      if (DEBUGMEM>1) err(warnmem,"poleval: i = %ld",i);
+      gerepileall(av, 2, &p1, &p2);
+    }
   }
-  p1=gmul(y,p1); tetpil=avma;
-  return gerepile(av,tetpil,gadd(p1,p2));
+  return gerepileupto(av0, gadd(p2, gmul(y,p1)));
 }