===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/language/Attic/sumiter.c,v
retrieving revision 1.1.1.1
retrieving revision 1.2
diff -u -p -r1.1.1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/language/Attic/sumiter.c	2001/10/02 11:17:10	1.1.1.1
+++ OpenXM_contrib/pari-2.2/src/language/Attic/sumiter.c	2002/09/11 07:27:04	1.2
@@ -1,4 +1,4 @@
-/* $Id: sumiter.c,v 1.1.1.1 2001/10/02 11:17:10 noro Exp $
+/* $Id: sumiter.c,v 1.2 2002/09/11 07:27:04 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -21,6 +21,7 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #include "pari.h"
 #include "anal.h"
 extern void changevalue_p(entree *ep, GEN x);
+extern GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);
 
 /********************************************************************/
 /**                                                                **/
@@ -31,7 +32,7 @@ extern void changevalue_p(entree *ep, GEN x);
 void
 forpari(entree *ep, GEN a, GEN b, char *ch)
 {
-  ulong av,av0 = avma, lim;
+  gpmem_t av, av0 = avma, lim;
 
   b = gcopy(b); av=avma; lim = stack_lim(av,1);
  /* gcopy nedeed in case b gets overwritten in ch, as in
@@ -40,7 +41,7 @@ forpari(entree *ep, GEN a, GEN b, char *ch)
   push_val(ep, a);
   while (gcmp(a,b) <= 0)
   {
-    ulong av1=avma; (void)lisseq(ch); avma=av1;
+    gpmem_t av1=avma; (void)lisseq(ch); avma=av1;
     if (loop_break()) break;
     a = (GEN) ep->value; a = gadd(a,gun);
     if (low_stack(lim, stack_lim(av,1)))
@@ -58,7 +59,8 @@ static int negcmp(GEN x, GEN y) { return gcmp(y,x); }
 void
 forstep(entree *ep, GEN a, GEN b, GEN s, char *ch)
 {
-  long ss, av,av0 = avma, lim, i;
+  long ss, i;
+  gpmem_t av, av0 = avma, lim;
   GEN v = NULL;
   int (*cmp)(GEN,GEN);
 
@@ -75,7 +77,7 @@ forstep(entree *ep, GEN a, GEN b, GEN s, char *ch)
   i = 0;
   while (cmp(a,b) <= 0)
   {
-    long av1=avma; (void)lisseq(ch); avma=av1;
+    gpmem_t av1=avma; (void)lisseq(ch); avma=av1;
     if (loop_break()) break;
     if (v)
     {
@@ -103,7 +105,8 @@ sinitp(long a, long c, byteptr *ptr)
   byteptr p = *ptr;
   if (a <= 0) a = 2;
   if (maxprime() < (ulong)a) err(primer1);
-  while (a > c) c += *p++;
+  while (a > c)
+     NEXT_PRIME_VIADIFF(c,p);
   *ptr = p; return c;
 }
 
@@ -125,7 +128,7 @@ update_p(entree *ep, byteptr *ptr, long prime[])
     *ptr = diffptr;
     prime[2] = sinitp(a, 0, ptr);
   }
-  changevalue_p(ep, prime); 
+  changevalue_p(ep, prime);
 }
 
 static byteptr
@@ -153,8 +156,9 @@ prime_loop_init(GEN ga, GEN gb, long *a, long *b, long
 void
 forprime(entree *ep, GEN ga, GEN gb, char *ch)
 {
-  long prime[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
-  long av = avma, a,b;
+  long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
+  long a, b;
+  gpmem_t av = avma;
   byteptr p;
 
   p = prime_loop_init(ga,gb, &a,&b, prime);
@@ -179,13 +183,14 @@ forprime(entree *ep, GEN ga, GEN gb, char *ch)
 void
 fordiv(GEN a, entree *ep, char *ch)
 {
-  long i,av2,l, av = avma;
+  long i, l;
+  gpmem_t av2, av = avma;
   GEN t = divisors(a);
 
   push_val(ep, NULL); l=lg(t); av2 = avma;
   for (i=1; i<l; i++)
   {
-    ep->value = (void*) t[i]; 
+    ep->value = (void*) t[i];
     (void)lisseq(ch); if (loop_break()) break;
     avma = av2;
   }
@@ -198,111 +203,107 @@ fordiv(GEN a, entree *ep, char *ch)
  *   fl = 1: impose a1 <= ... <= an
  *   fl = 2:        a1 <  ... <  an
  */
-static GEN *fv_a, *fv_m, *fv_M;
-static long fv_n, fv_fl;
-static char *fv_ch;
+typedef struct {
+  GEN *a, *m, *M;
+  long n,fl;
+  char *ch;
+} fvdat;
 
 /* general case */
 static void
-fvloop(long i)
+fvloop(long i, fvdat *d)
 {
-  fv_a[i] = fv_m[i];
-  if (fv_fl && i > 1)
+  d->a[i] = d->m[i];
+  if (d->fl && i > 1)
   {
-    GEN p1 = gsub(fv_a[i], fv_a[i-1]);
+    GEN p1 = gsub(d->a[i], d->a[i-1]);
     if (gsigne(p1) < 0)
-      fv_a[i] = gadd(fv_a[i], gceil(gneg_i(p1)));
-    if (fv_fl == 2 && gegal(fv_a[i], fv_a[i-1]))
-      fv_a[i] = gadd(fv_a[i], gun);
+      d->a[i] = gadd(d->a[i], gceil(gneg_i(p1)));
+    if (d->fl == 2 && gegal(d->a[i], d->a[i-1]))
+      d->a[i] = gadd(d->a[i], gun);
   }
-  if (i+1 == fv_n)
-    while (gcmp(fv_a[i], fv_M[i]) <= 0)
+  if (i+1 == d->n)
+    while (gcmp(d->a[i], d->M[i]) <= 0)
     {
-      long av = avma; (void)lisseq(fv_ch); avma = av;
-      if (loop_break()) { fv_n = 0; return; }
-      fv_a[i] = gadd(fv_a[i], gun);
+      gpmem_t av = avma; (void)lisseq(d->ch); avma = av;
+      if (loop_break()) { d->n = 0; return; }
+      d->a[i] = gadd(d->a[i], gun);
     }
   else
-    while (gcmp(fv_a[i], fv_M[i]) <= 0)
+    while (gcmp(d->a[i], d->M[i]) <= 0)
     {
-      long av = avma; fvloop(i+1); avma = av;
-      if (!fv_n) return;
-      fv_a[i] = gadd(fv_a[i], gun);
+      gpmem_t av = avma; fvloop(i+1, d); avma = av;
+      if (!d->n) return;
+      d->a[i] = gadd(d->a[i], gun);
     }
 }
 
 /* we run through integers */
 static void
-fvloop_i(long i)
+fvloop_i(long i, fvdat *d)
 {
-  fv_a[i] = setloop(fv_m[i]);
-  if (fv_fl && i > 1)
+  d->a[i] = setloop(d->m[i]);
+  if (d->fl && i > 1)
   {
-    int c = cmpii(fv_a[i], fv_a[i-1]);
+    int c = cmpii(d->a[i], d->a[i-1]);
     if (c < 0)
     {
-      fv_a[i] = setloop(fv_a[i-1]);
+      d->a[i] = setloop(d->a[i-1]);
       c = 0;
     }
-    if (c == 0 && fv_fl == 2)
-      fv_a[i] = incloop(fv_a[i]);
+    if (c == 0 && d->fl == 2)
+      d->a[i] = incloop(d->a[i]);
   }
-  if (i+1 == fv_n)
-    while (gcmp(fv_a[i], fv_M[i]) <= 0)
+  if (i+1 == d->n)
+    while (gcmp(d->a[i], d->M[i]) <= 0)
     {
-      long av = avma; (void)lisseq(fv_ch); avma = av;
-      if (loop_break()) { fv_n = 0; return; }
-      fv_a[i] = incloop(fv_a[i]);
+      gpmem_t av = avma; (void)lisseq(d->ch); avma = av;
+      if (loop_break()) { d->n = 0; return; }
+      d->a[i] = incloop(d->a[i]);
     }
   else
-    while (gcmp(fv_a[i], fv_M[i]) <= 0)
+    while (gcmp(d->a[i], d->M[i]) <= 0)
     {
-      long av = avma; fvloop_i(i+1); avma = av;
-      if (!fv_n) return;
-      fv_a[i] = incloop(fv_a[i]);
+      gpmem_t av = avma; fvloop_i(i+1, d); avma = av;
+      if (!d->n) return;
+      d->a[i] = incloop(d->a[i]);
     }
 }
 
 void
 forvec(entree *ep, GEN x, char *c, long flag)
 {
-  long i, av = avma, tx = typ(x), n = fv_n, fl = fv_fl;
-  GEN *a = fv_a, *M = fv_M, *m = fv_m;
-  char *ch = fv_ch;
-  void (*fv_fun)(long) = fvloop_i;
+  gpmem_t av = avma;
+  long tx = typ(x);
+  fvdat D, *d = &D;
 
   if (!is_vec_t(tx)) err(talker,"not a vector in forvec");
-  if (fl<0 || fl>2) err(flagerr);
-  fv_n = lg(x);
-  fv_ch = c;
-  fv_fl = flag;
-  fv_a = (GEN*)cgetg(fv_n,t_VEC); push_val(ep, (GEN)fv_a);
-  fv_m = (GEN*)cgetg(fv_n,t_VEC);
-  fv_M = (GEN*)cgetg(fv_n,t_VEC);
-  if (fv_n == 1) lisseq(fv_ch);
+  if (flag<0 || flag>2) err(flagerr);
+  d->n = lg(x);
+  d->ch = c;
+  d->a = (GEN*)cgetg(d->n,t_VEC); push_val(ep, (GEN)d->a);
+  if (d->n == 1) (void)lisseq(d->ch);
   else
   {
-    for (i=1; i<fv_n; i++)
+    long i, t = t_INT;
+    d->fl = flag;
+    d->m = (GEN*)cgetg(d->n,t_VEC);
+    d->M = (GEN*)cgetg(d->n,t_VEC);
+    for (i=1; i<d->n; i++)
     {
-      GEN *c = (GEN*) x[i];
-      tx = typ(c);
-      if (! is_vec_t(tx) || lg(c)!=3)
+      GEN *e = (GEN*) x[i];
+      tx = typ(e);
+      if (! is_vec_t(tx) || lg(e)!=3)
 	err(talker,"not a vector of two-component vectors in forvec");
-      if (gcmp(c[1],c[2]) > 0) fv_n = 0;
-      /* in case x is an ep->value and c destroys it, we have to copy */
-      if (typ(c[1]) != t_INT) fv_fun = fvloop;
-      fv_m[i] = gcopy(c[1]);
-      fv_M[i] = gcopy(c[2]);
+      if (gcmp(e[1],e[2]) > 0) d->n = 0;
+      if (typ(e[1]) != t_INT) t = t_REAL;
+      /* in case x is an ep->value and lisexpr(d->ch) kills it, have to copy */
+      d->m[i] = gcopy(e[1]);
+      d->M[i] = gcopy(e[2]);
     }
-    fv_fun(1);
+    if (t == t_INT) fvloop_i(1, d); else fvloop(1, d);
   }
-  pop_val(ep);
-  fv_n = n;
-  fv_ch = ch;
-  fv_fl = fl;
-  fv_a = a;
-  fv_m = m;
-  fv_M = M; avma = av;
+  pop_val(ep); avma = av;
 }
 
 /********************************************************************/
@@ -314,14 +315,14 @@ forvec(entree *ep, GEN x, char *c, long flag)
 GEN
 somme(entree *ep, GEN a, GEN b, char *ch, GEN x)
 {
-  long av,av0 = avma, lim;
+  gpmem_t av, av0 = avma, lim;
   GEN p1;
 
   if (typ(a) != t_INT) err(talker,"non integral index in sum");
   if (!x) x = gzero;
   if (gcmp(b,a) < 0) return gcopy(x);
 
-  b = gfloor(b); 
+  b = gfloor(b);
   a = setloop(a);
   av=avma; lim = stack_lim(av,1);
   push_val(ep, a);
@@ -343,7 +344,8 @@ somme(entree *ep, GEN a, GEN b, char *ch, GEN x)
 GEN
 suminf(entree *ep, GEN a, char *ch, long prec)
 {
-  long fl,G,tetpil, av0 = avma, av,lim;
+  long fl, G;
+  gpmem_t tetpil, av0 = avma, av, lim;
   GEN p1,x = realun(prec);
 
   if (typ(a) != t_INT) err(talker,"non integral index in suminf");
@@ -373,14 +375,14 @@ suminf(entree *ep, GEN a, char *ch, long prec)
 GEN
 divsum(GEN num, entree *ep, char *ch)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN z, y = gzero, t = divisors(num);
   long i, l = lg(t);
 
   push_val(ep, NULL);
   for (i=1; i<l; i++)
   {
-    ep->value = (void*) t[i]; 
+    ep->value = (void*) t[i];
     z = lisseq(ch);
     if (did_break()) err(breaker,"divsum");
     y = gadd(y, z);
@@ -397,14 +399,14 @@ divsum(GEN num, entree *ep, char *ch)
 GEN
 produit(entree *ep, GEN a, GEN b, char *ch, GEN x)
 {
-  long av,av0 = avma, lim;
+  gpmem_t av, av0 = avma, lim;
   GEN p1;
 
   if (typ(a) != t_INT) err(talker,"non integral index in sum");
   if (!x) x = gun;
   if (gcmp(b,a) < 0) return gcopy(x);
 
-  b = gfloor(b); 
+  b = gfloor(b);
   a = setloop(a);
   av=avma; lim = stack_lim(av,1);
   push_val(ep, a);
@@ -438,7 +440,7 @@ prodinf0(entree *ep, GEN a, char *ch, long flag, long 
 GEN
 prodinf(entree *ep, GEN a, char *ch, long prec)
 {
-  ulong av0 = avma, av,lim;
+  gpmem_t av0 = avma, av, lim;
   long fl,G;
   GEN p1,x = realun(prec);
 
@@ -465,7 +467,7 @@ prodinf(entree *ep, GEN a, char *ch, long prec)
 GEN
 prodinf1(entree *ep, GEN a, char *ch, long prec)
 {
-  ulong av0 = avma, av,lim;
+  gpmem_t av0 = avma, av, lim;
   long fl,G;
   GEN p1,p2,x = realun(prec);
 
@@ -492,12 +494,12 @@ prodinf1(entree *ep, GEN a, char *ch, long prec)
 GEN
 prodeuler(entree *ep, GEN ga, GEN gb, char *ch, long prec)
 {
-  long prime[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
+  long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
   long a,b;
-  ulong av,av0 = avma, lim;
+  gpmem_t av, av0 = avma, lim;
   GEN p1,x = realun(prec);
   byteptr p;
-  
+
   av = avma;
   p = prime_loop_init(ga,gb, &a,&b, prime);
   if (!p) { avma = av; return x; }
@@ -530,15 +532,15 @@ prodeuler(entree *ep, GEN ga, GEN gb, char *ch, long p
 GEN
 direulerall(entree *ep, GEN ga, GEN gb, char *ch, GEN c)
 {
-  long prime[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
-  ulong av0 = avma,av, lim = (av0+bot)>>1;
+  long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
+  gpmem_t av0 = avma, av, lim = stack_lim(av0, 1);
   long p,n,i,j,k,tx,lx,a,b;
   GEN x,y,s,polnum,polden;
   byteptr d;
 
   d = prime_loop_init(ga,gb, &a,&b, prime);
   n = c? itos(c): b;
-  if (!d || b < 2 || n <= 0) { x=cgetg(2,t_VEC); x[1]=un; return x; }
+  if (!d || b < 2 || n <= 0) return _vec(gun);
   if (n < b) b = n;
   push_val(ep, prime);
 
@@ -636,7 +638,7 @@ vecteur(GEN nmax, entree *ep, char *ch)
 {
   GEN y,p1;
   long i,m;
-  long c[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
+  long c[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
 
   if (typ(nmax)!=t_INT || signe(nmax) < 0)
     err(talker,"bad number of components in vector");
@@ -657,6 +659,32 @@ vecteur(GEN nmax, entree *ep, char *ch)
 }
 
 GEN
+vecteursmall(GEN nmax, entree *ep, char *ch)
+{
+  GEN y,p1;
+  long i,m;
+  long c[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
+
+  if (typ(nmax)!=t_INT || signe(nmax) < 0)
+    err(talker,"bad number of components in vector");
+  m=itos(nmax); y=cgetg(m+1,t_VECSMALL);
+  if (!ep || !ch)
+  {
+    for (i=1; i<=m; i++) y[i] = 0;
+    return y;
+  }
+  push_val(ep, c);
+  for (i=1; i<=m; i++)
+  {
+    c[2]=i; p1 = lisseq(ch);
+    if (did_break()) err(breaker,"vector");
+    y[i] = itos(p1);
+  }
+  pop_val(ep); return y;
+}
+
+
+GEN
 vvecteur(GEN nmax, entree *ep, char *ch)
 {
   GEN y = vecteur(nmax,ep,ch);
@@ -668,11 +696,11 @@ matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, c
 {
   GEN y,z,p1;
   long i,j,m,n,s;
-  long c1[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 1};
-  long c2[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 1};
+  long c1[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 1};
+  long c2[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 1};
 
   s=signe(ncol);
-  if (ep1 == ep2) err(talker, "identical index variables in matrix");
+  if (ep1 == ep2 && ep1) err(talker, "identical index variables in matrix");
   if (typ(ncol)!=t_INT || s<0) err(talker,"bad number of columns in matrix");
   if (!s) return cgetg(1,t_MAT);
 
@@ -712,7 +740,7 @@ matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, c
 
 /********************************************************************/
 /**                                                                **/
-/**                    SOMMATION DE SERIES                         **/
+/**                         SUMMING SERIES                         **/
 /**                                                                **/
 /********************************************************************/
 
@@ -743,16 +771,19 @@ sumpos0(entree *ep, GEN a, char *ch, long flag, long p
 GEN
 sumalt(entree *ep, GEN a, char *ch, long prec)
 {
-  long av = avma, tetpil,k,N;
+  long k, N;
+  gpmem_t av = avma, tetpil;
   GEN s,az,c,x,e1,d;
 
   if (typ(a) != t_INT) err(talker,"non integral index in sumalt");
 
   push_val(ep, a);
-  e1=addsr(3,gsqrt(stoi(8),prec));
-  N=(long)(0.4*(bit_accuracy(prec) + 7));
-  d=gpuigs(e1,N); d=shiftr(addrr(d,divsr(1,d)),-1);
-  az=negi(gun); c=d; s=gzero;
+  e1 = addsr(3,gsqrt(stoi(8),prec));
+  N = (long)(0.4*(bit_accuracy(prec) + 7));
+  d = gpowgs(e1,N);
+  d = shiftr(addrr(d, ginv(d)),-1);
+  az = negi(gun); c = d;
+  s = gzero;
   for (k=0; ; k++)
   {
     x = lisexpr(ch); if (did_break()) err(breaker,"sumalt");
@@ -761,56 +792,62 @@ sumalt(entree *ep, GEN a, char *ch, long prec)
     if (k==N-1) break;
     a = addsi(1,a); ep->value = (void*) a;
   }
-  tetpil=avma; pop_val(ep);
+  tetpil = avma; pop_val(ep);
   return gerepile(av,tetpil,gdiv(s,d));
 }
 
 GEN
 sumalt2(entree *ep, GEN a, char *ch, long prec)
 {
-  long av = avma, tetpil,k,N;
+  long k, N;
+  gpmem_t av = avma, tetpil;
   GEN x,s,dn,pol;
 
   if (typ(a) != t_INT) err(talker,"non integral index in sumalt");
 
   push_val(ep, a);
-  N=(long)(0.31*(bit_accuracy(prec) + 5));
-  s=gzero; pol=polzagreel(N,N>>1,prec+1); dn=poleval(pol,gun);
-  pol[2]=lsub((GEN)pol[2],dn); pol=gdiv(pol,gsub(polx[0],gun));
-  N=lgef(pol)-2;
-  for (k=0; k<N; k++)
+  N = (long)(0.31*(bit_accuracy(prec) + 5));
+  pol = polzagreel(N,N>>1,prec+1);
+  dn = poleval(pol,gun);
+  pol[2] = lsub((GEN)pol[2],dn);
+  pol = gdiv(pol, gsub(polx[0],gun));
+  N = degpol(pol);
+  s = gzero;
+  for (k=0; k<=N; k++)
   {
     x = lisexpr(ch); if (did_break()) err(breaker,"sumalt2");
-    s=gadd(s,gmul(x,(GEN)pol[k+2]));
-    if (k==N-1) break;
+    s = gadd(s, gmul(x,(GEN)pol[k+2]));
+    if (k == N) break;
     a = addsi(1,a); ep->value = (void*) a;
   }
-  tetpil=avma; pop_val(ep);
-  return gerepile(av,tetpil,gdiv(s,dn));
+  tetpil = avma; pop_val(ep);
+  return gerepile(av,tetpil, gdiv(s,dn));
 }
 
 GEN
 sumpos(entree *ep, GEN a, char *ch, long prec)
 {
-  long av = avma,tetpil,k,kk,N,G;
+  long k, kk, N, G;
+  gpmem_t av = avma, tetpil;
   GEN p1,r,q1,reel,s,az,c,x,e1,d, *stock;
 
   if (typ(a) != t_INT) err(talker,"non integral index in sumpos");
 
   push_val(ep, NULL);
-  a=subis(a,1); reel=cgetr(prec);
-  e1=addsr(3,gsqrt(stoi(8),prec));
-  N=(long)(0.4*(bit_accuracy(prec) + 7));
-  d=gpuigs(e1,N); d=shiftr(addrr(d,divsr(1,d)),-1);
-  az=negi(gun); c=d; s=gzero;
+  a = subis(a,1); reel = cgetr(prec);
+  e1 = addsr(3,gsqrt(stoi(8),prec));
+  N = (long)(0.4*(bit_accuracy(prec) + 7));
+  d = gpowgs(e1,N); d = shiftr(addrr(d, ginv(d)),-1);
+  az = negi(gun); c = d; s = gzero;
+
   G = -bit_accuracy(prec) - 5;
-  stock=(GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k]=NULL;
+  stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;
   for (k=0; k<N; k++)
   {
-    if (k&1 && stock[k]) x=stock[k];
+    if (odd(k) && stock[k]) x = stock[k];
     else
     {
-      x=gzero; r=stoi(2*k+2);
+      x = gzero; r = stoi(2*k+2);
       for(kk=0;;kk++)
       {
         long ex;
@@ -818,14 +855,14 @@ sumpos(entree *ep, GEN a, char *ch, long prec)
         p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");
         gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);
 	x = gadd(x,reel); if (kk && ex < G) break;
-        r=shifti(r,1);
+        r = shifti(r,1);
       }
-      if (2*k<N) stock[2*k+1]=x;
+      if (2*k < N) stock[2*k+1] = x;
       q1 = addsi(k+1,a); ep->value = (void*) q1;
       p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");
       gaffect(p1,reel); x = gadd(reel, gmul2n(x,1));
     }
-    c=gadd(az,c); p1 = k&1? gneg_i(c): c;
+    c = gadd(az,c); p1 = k&1? gneg_i(c): c;
     s = gadd(s,gmul(x,p1));
     az = divii(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1));
   }
@@ -836,21 +873,23 @@ sumpos(entree *ep, GEN a, char *ch, long prec)
 GEN
 sumpos2(entree *ep, GEN a, char *ch, long prec)
 {
-  long av = avma,tetpil,k,kk,N,G;
+  long k, kk, N, G;
+  gpmem_t av = avma, tetpil;
   GEN p1,r,q1,reel,s,pol,dn,x, *stock;
 
   if (typ(a) != t_INT) err(talker,"non integral index in sumpos2");
 
   push_val(ep, a);
-  a=subis(a,1); reel=cgetr(prec);
-  N=(long)(0.31*(bit_accuracy(prec) + 5));
+  a = subis(a,1); reel = cgetr(prec);
+  N = (long)(0.31*(bit_accuracy(prec) + 5));
+
   G = -bit_accuracy(prec) - 5;
-  stock=(GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k]=NULL;
+  stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;
   for (k=1; k<=N; k++)
   {
     if (odd(k) || !stock[k])
     {
-      x=gzero; r=stoi(2*k);
+      x = gzero; r = stoi(2*k);
       for(kk=0;;kk++)
       {
         long ex;
@@ -858,21 +897,23 @@ sumpos2(entree *ep, GEN a, char *ch, long prec)
         p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");
         gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);
 	x=gadd(x,reel); if (kk && ex < G) break;
-        r=shifti(r,1);
+        r = shifti(r,1);
       }
-      if (2*k-1<N) stock[2*k]=x;
+      if (2*k-1 < N) stock[2*k] = x;
       q1 = addsi(k,a); ep->value = (void*) q1;
       p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");
       gaffect(p1,reel); stock[k] = gadd(reel, gmul2n(x,1));
     }
   }
   pop_val(ep); s = gzero;
-  pol=polzagreel(N,N>>1,prec+1); dn=poleval(pol,gun);
-  pol[2]=lsub((GEN)pol[2],dn); pol=gdiv(pol,gsub(gun,polx[0]));
+  pol = polzagreel(N,N>>1,prec+1);
+  dn = poleval(pol,gun);
+  pol[2] = lsub((GEN)pol[2],dn);
+  pol = gdiv(pol, gsub(gun,polx[0]));
   for (k=1; k<=lgef(pol)-2; k++)
   {
     p1 = gmul((GEN)pol[k+1],stock[k]);
-    if (k&1) p1 = gneg_i(p1);
+    if (odd(k)) p1 = gneg_i(p1);
     s = gadd(s,p1);
   }
   tetpil=avma; return gerepile(av,tetpil,gdiv(s,dn));
@@ -883,23 +924,15 @@ sumpos2(entree *ep, GEN a, char *ch, long prec)
 /**                NUMERICAL INTEGRATION (Romberg)                 **/
 /**                                                                **/
 /********************************************************************/
-GEN
-intnum0(entree *ep, GEN a, GEN b, char *ch, long flag, long prec)
-{
-  switch(flag)
-  {
-    case 0: return qromb  (ep,a,b,ch,prec);
-    case 1: return rombint(ep,a,b,ch,prec);
-    case 2: return qromi  (ep,a,b,ch,prec);
-    case 3: return qromo  (ep,a,b,ch,prec);
-    default: err(flagerr);
-  }
-  return NULL; /* not reached */
-}
+typedef struct {
+  entree *ep;
+  char *ch;
+} exprdat;
 
-#define JMAX 25
-#define JMAXP JMAX+3
-#define KLOC 4
+typedef struct {
+  GEN (*f)(void *, GEN);
+  void *E;
+} invfun;
 
 /* we need to make a copy in any case, cf forpari */
 static GEN
@@ -909,221 +942,242 @@ fix(GEN a, long prec)
   gaffect(a,p); return p;
 }
 
-extern GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);
-
+#if 0
 GEN
-qromb(entree *ep, GEN a, GEN b, char *ch, long prec)
+int_loop(entree *ep, char *ch)
 {
-  GEN ss,dss,s,h,p1,p2,qlint,del,x,sum;
-  long av = avma, av1, tetpil,j,j1,j2,lim,it,sig;
+  gpmem_t av = avma;
+  intstruct T;
 
+  T.in = NULL;
+  for(;;)
+  {
+    GEN x = Next(&T);
+    if (!x) return gerepileupto(av, T.out);
+    ep->value = (void*)x;
+    T.in = lisexpr(ch);
+  }
+}
+#endif
+
+/* f(x) */
+static GEN
+_gp_eval(void *dat, GEN x)
+{
+  exprdat *E = (exprdat*)dat;
+  E->ep->value = x;
+  return lisexpr(E->ch);
+}
+/* 1/x^2 f(1/x) */
+static GEN
+_invf(void *dat, GEN x)
+{
+  invfun *S = (invfun*)dat;
+  GEN y = ginv(x);
+  return gmul(S->f(S->E, y), gsqr(y));
+}
+
+#define swap(a,b) { GEN _x = a; a = b; b = _x; }
+static GEN
+interp(GEN h, GEN s, long j, long lim, long KLOC)
+{
+  gpmem_t av = avma;
+  long e1,e2;
+  GEN dss, ss = polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
+
+  e1 = gexpo(ss);
+  e2 = gexpo(dss);
+  if (e1-e2 <= lim && e1 >= -lim) { avma = av; return NULL; }
+  if (gcmp0(gimag(ss))) ss = greal(ss);
+  return ss;
+}
+
+static GEN
+qrom3(void *dat, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
+{
+  const long JMAX = 25, KLOC = 4;
+  GEN ss,s,h,p1,p2,qlint,del,x,sum;
+  long j, j1, it, sig;
+  gpmem_t av;
+
   a = fix(a,prec);
   b = fix(b,prec);
-  qlint=subrr(b,a); sig=signe(qlint);
-  if (!sig)  { avma=av; return gzero; }
-  if (sig<0) { setsigne(qlint,1); s=a; a=b; b=s; }
+  qlint = subrr(b,a); sig = signe(qlint);
+  if (!sig)  return gzero;
+  if (sig < 0) { setsigne(qlint,1); swap(a,b); }
 
-  s=new_chunk(JMAXP);
-  h=new_chunk(JMAXP);
+  s = new_chunk(JMAX+KLOC-1);
+  h = new_chunk(JMAX+KLOC-1);
   h[0] = (long)realun(prec);
 
-  push_val(ep, a);
-  p1=lisexpr(ch); if (p1 == a) p1=rcopy(p1);
-  ep->value = (void*)b;
-  p2=lisexpr(ch);
-  s[0]=lmul2n(gmul(qlint,gadd(p1,p2)),-1);
-  for (it=1,j=1; j<JMAX; j++,it=it<<1)
+  p1 = eval(dat, a); if (p1 == a) p1 = rcopy(p1);
+  p2 = eval(dat, b);
+  s[0] = lmul2n(gmul(qlint,gadd(p1,p2)),-1);
+  for (it=1,j=1; j<JMAX; j++, it<<=1)
   {
-    h[j]=lshiftr((GEN)h[j-1],-2);
-    av1=avma; del=divrs(qlint,it); x=addrr(a,shiftr(del,-1));
-    for (sum=gzero,j1=1; j1<=it; j1++,x=addrr(x,del))
-    {
-      ep->value = (void*) x; sum=gadd(sum,lisexpr(ch));
-    }
-    sum = gmul(sum,del); p1 = gadd((GEN)s[j-1],sum);
-    tetpil = avma;
-    s[j]=lpile(av1,tetpil,gmul2n(p1,-1));
+    h[j] = lshiftr((GEN)h[j-1],-2);
+    av = avma; del = divrs(qlint,it);
+    x = addrr(a, shiftr(del,-1));
+    for (sum = gzero, j1 = 1; j1 <= it; j1++, x = addrr(x,del))
+      sum = gadd(sum, eval(dat, x));
+    sum = gmul(sum,del); p1 = gadd((GEN)s[j-1], sum);
+    s[j] = lpileupto(av, gmul2n(p1,-1));
 
-    if (j>=KLOC)
-    {
-      tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
-      j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-j-6;
-      if (j1-j2 > lim || j1 < -lim)
-      {
-        pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);
-	tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));
-      }
-      avma=tetpil;
-    }
+    if (j >= KLOC && (ss = interp(h, s, j, bit_accuracy(prec)-j-6, KLOC)))
+      return gmulsg(sig,ss);
   }
-  err(intger2); return NULL; /* not reached */
+  return NULL;
 }
 
-GEN
-qromo(entree *ep, GEN a, GEN b, char *ch, long prec)
+static GEN
+qrom2(void *dat, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
 {
-  GEN ss,dss,s,h,p1,qlint,del,ddel,x,sum;
-  long av = avma,av1,tetpil,j,j1,j2,lim,it,sig;
+  const long JMAX = 16, KLOC = 4;
+  GEN ss,s,h,p1,qlint,del,ddel,x,sum;
+  long j, j1, it, sig;
+  gpmem_t av;
 
   a = fix(a, prec);
   b = fix(b, prec);
-  qlint=subrr(b,a); sig=signe(qlint);
-  if (!sig)  { avma=av; return gzero; }
-  if (sig<0) { setsigne(qlint,1); s=a; a=b; b=s; }
+  qlint = subrr(b,a); sig = signe(qlint);
+  if (!sig)  return gzero;
+  if (sig < 0) { setsigne(qlint,1); swap(a,b); }
 
-  s=new_chunk(JMAXP);
-  h=new_chunk(JMAXP);
+  s = new_chunk(JMAX+KLOC-1);
+  h = new_chunk(JMAX+KLOC-1);
   h[0] = (long)realun(prec);
 
-  p1 = shiftr(addrr(a,b),-1); push_val(ep, p1);
-  p1=lisexpr(ch); s[0]=lmul(qlint,p1);
-
+  p1 = shiftr(addrr(a,b),-1);
+  s[0] = lmul(qlint, eval(dat, p1));
   for (it=1, j=1; j<JMAX; j++, it*=3)
   {
     h[j] = ldivrs((GEN)h[j-1],9);
-    av1=avma; del=divrs(qlint,3*it); ddel=shiftr(del,1);
-    x=addrr(a,shiftr(del,-1)); sum=gzero;
-    for (j1=1; j1<=it; j1++)
+    av = avma; del = divrs(qlint,3*it); ddel = shiftr(del,1);
+    x = addrr(a, shiftr(del,-1));
+    for (sum = gzero, j1 = 1; j1 <= it; j1++)
     {
-      ep->value = (void*)x; sum=gadd(sum,lisexpr(ch)); x=addrr(x,ddel);
-      ep->value = (void*)x; sum=gadd(sum,lisexpr(ch)); x=addrr(x,del);
+      sum = gadd(sum, eval(dat, x)); x = addrr(x,ddel);
+      sum = gadd(sum, eval(dat, x)); x = addrr(x,del);
     }
     sum = gmul(sum,del); p1 = gdivgs((GEN)s[j-1],3);
-    tetpil = avma; s[j]=lpile(av1,tetpil,gadd(p1,sum));
+    s[j] = lpileupto(av, gadd(p1,sum));
 
-    if (j>=KLOC)
-    {
-      tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
-      j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6;
-      if ( j1-j2 > lim || j1 < -lim)
-      {
-        pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);
-	tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));
-      }
-      avma=tetpil;
-    }
+    if (j >= KLOC && (ss = interp(h, s, j, bit_accuracy(prec)-(3*j/2)-6, KLOC)))
+      return gmulsg(sig, ss);
   }
-  err(intger2); return NULL; /* not reached */
+  return NULL;
 }
 
-#undef JMAX
-#undef JMAXP
-#define JMAX 16
-#define JMAXP JMAX+3
+/* integrate after change of variables x --> 1/x */
+static GEN
+qromi(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
+{
+  GEN A = ginv(b), B = ginv(a);
+  invfun S;
+  S.f = eval;
+  S.E = E; return qrom2(&S, &_invf, A, B, prec);
+}
 
-GEN
-qromi(entree *ep, GEN a, GEN b, char *ch, long prec)
+/* a < b, assume b "small" (< 100 say) */
+static GEN
+rom_bsmall(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
 {
-  GEN ss,dss,s,h,q1,p1,p,qlint,del,ddel,x,sum;
-  long av = avma, av1,tetpil,j,j1,j2,lim,it,sig;
+  if (gcmpgs(a,-100) >= 0) return qrom2(E,eval,a,b,prec);
+  if (b == gun || gcmpgs(b, -1) >= 0)
+  { /* a < -100, b >= -1 */
+    GEN _1 = negi(gun); /* split at -1 */
+    return gadd(qromi(E,eval,a,_1,prec),
+                qrom2(E,eval,_1,b,prec));
+  }
+  /* a < -100, b < -1 */
+  return qromi(E,eval,a,b,prec);
+}
 
-  p=cgetr(prec); gaffect(ginv(a),p); a=p;
-  p=cgetr(prec); gaffect(ginv(b),p); b=p;
-  qlint=subrr(b,a); sig= -signe(qlint);
-  if (!sig)  { avma=av; return gzero; }
-  if (sig>0) { setsigne(qlint,1); s=a; a=b; b=s; }
+static GEN
+rombint(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
+{
+  long l = gcmp(b,a);
+  GEN z;
 
-  s=new_chunk(JMAXP);
-  h=new_chunk(JMAXP);
-  h[0] = (long)realun(prec);
-
-  x=divsr(2,addrr(a,b)); push_val(ep, x);
-  p1=gmul(lisexpr(ch),mulrr(x,x));
-  s[0]=lmul(qlint,p1);
-  for (it=1,j=1; j<JMAX; j++, it*=3)
+  if (!l) return gzero;
+  if (l < 0) swap(a,b);
+  if (gcmpgs(b,100) >= 0)
   {
-    h[j] = ldivrs((GEN)h[j-1],9);
-    av1=avma; del=divrs(qlint,3*it); ddel=shiftr(del,1);
-    x=addrr(a,shiftr(del,-1)); sum=gzero;
-    for (j1=1; j1<=it; j1++)
-    {
-      q1 = ginv(x); ep->value = (void*)q1;
-      p1=gmul(lisexpr(ch), gsqr(q1));
-      sum=gadd(sum,p1); x=addrr(x,ddel);
-
-      q1 = ginv(x); ep->value = (void*)q1;
-      p1=gmul(lisexpr(ch), gsqr(q1));
-      sum=gadd(sum,p1); x=addrr(x,del);
-    }
-    sum = gmul(sum,del); p1 = gdivgs((GEN)s[j-1],3);
-    tetpil=avma;
-    s[j]=lpile(av1,tetpil,gadd(p1,sum));
-
-    if (j>=KLOC)
-    {
-      tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
-      j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6;
-      if (j1-j2 > lim || j1 < -lim)
-      {
-        pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);
-	tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));
-      }
-    }
+    if (gcmpgs(a,1) >= 0)
+      z = qromi(E,eval,a,b,prec);
+    else /* split at 1 */
+      z = gadd(rom_bsmall(E,eval,a,gun,prec), qromi(E,eval,gun,b,prec));
   }
-  err(intger2); return NULL; /* not reached */
+  else
+    z = rom_bsmall(E,eval,a,b,prec);
+  if (l < 0) z = gneg(z);
+  return z;
 }
 
 GEN
-rombint(entree *ep, GEN a, GEN b, char *ch, long prec)
+intnum(void *E, GEN (*eval)(void*,GEN), GEN a, GEN b, long flag, long prec)
 {
-  GEN aa,bb,mun,p1,p2,p3;
-  long l,av,tetpil;
-
-  l=gcmp(b,a); if (!l) return gzero;
-  if (l<0) { bb=a; aa=b; } else { bb=b; aa=a; }
-  av=avma; mun = negi(gun);
-  if (gcmpgs(bb,100)>=0)
+  gpmem_t av = avma;
+  GEN z;
+  switch(flag)
   {
-    if (gcmpgs(aa,1)>=0) return qromi(ep,a,b,ch,prec);
-    p1=qromi(ep,gun,bb,ch,prec);
-    if (gcmpgs(aa,-100)>=0)
-    {
-      p2=qromo(ep,aa,gun,ch,prec); tetpil=avma;
-      return gerepile(av,tetpil,gmulsg(l,gadd(p1,p2)));
-    }
-    p2=qromo(ep,mun,gun,ch,prec); p3=gadd(p2,qromi(ep,aa,mun,ch,prec));
-    tetpil=avma; return gerepile(av,tetpil,gmulsg(l,gadd(p1,p3)));
+    case 0: z = qrom3  (E,eval,a,b,prec); break;
+    case 1: z = rombint(E,eval,a,b,prec); break;
+    case 2: z = qromi  (E,eval,a,b,prec); break;
+    case 3: z = qrom2  (E,eval,a,b,prec); break;
+    default: err(flagerr); z = NULL;
   }
-  if (gcmpgs(aa,-100)>=0) return qromo(ep,a,b,ch,prec);
-  if (gcmpgs(bb,-1)>=0)
-  {
-    p1=qromi(ep,aa,mun,ch,prec); p2=qromo(ep,mun,bb,ch,prec); tetpil=avma;
-    return gerepile(av,tetpil,gmulsg(l,gadd(p1,p2)));
-  }
-  return qromi(ep,a,b,ch,prec);
+  if (!z) err(intger2);
+  return gerepileupto(av, z);
 }
 
+GEN
+intnum0(entree *ep, GEN a, GEN b, char *ch, long flag, long prec)
+{
+  exprdat E;
+  GEN z;
+
+  E.ch = ch;
+  E.ep = ep; push_val(ep, NULL);
+  z = intnum(&E, &_gp_eval, a,b, flag, prec);
+  pop_val(ep); return z;
+}
 /********************************************************************/
 /**                                                                **/
-/**                  ZAGIER POLYNOMIALS (for sumiter)              **/
+/**                  ZAGIER POLYNOMIALS (for sumalt)               **/
 /**                                                                **/
 /********************************************************************/
 
 GEN
 polzag(long n, long m)
 {
-  long d1,d,r,k,av=avma,tetpil;
-  GEN p1,p2,pol1,g,s;
+  gpmem_t av = avma;
+  long d1, d, r, k;
+  GEN A, B, Bx, g, s;
 
-  if (m>=n || m<0)
+  if (m >= n || m < 0)
     err(talker,"first index must be greater than second in polzag");
-  d1=n-m; d=d1<<1; d1--; p1=gsub(gun,gmul2n(polx[0],1));
-  pol1=gsub(gun,polx[0]); p2=gmul(polx[0],pol1); r=(m+1)>>1;
-  g=gzero;
+  d1 = n-m; d = d1<<1; d1--;
+  A = gsub(gun, gmul2n(polx[0],1)); /* 1 - 2x */
+  B = gsub(gun, polx[0]); /* 1 - x */
+  Bx = gmul(polx[0],B); /* x - x^2 */
+  r = (m+1)>>1;
+  g = gzero;
   for (k=0; k<=d1; k++)
   {
-    s=binome(stoi(d),(k<<1)+1); if (k&1) s=negi(s);
-    g=gadd(g,gmul(s,gmul(gpuigs(polx[0],k),gpuigs(pol1,d1-k))));
+    s = binome(stoi(d), (k<<1)+1); if (k&1) s = negi(s);
+    g = gadd(g, gmul(s, gmul(gpowgs(polx[0],k), gpowgs(B,d1-k))));
   }
-  g=gmul(g,gpuigs(p2,r));
-  if (!(m&1)) g=gadd(gmul(p1,g),gmul2n(gmul(p2,derivpol(g)),1));
+  g = gmul(g, gpowgs(Bx,r));
+  if (!(m&1)) g = gadd(gmul(A,g), gmul2n(gmul(Bx,derivpol(g)),1));
   for (k=1; k<=r; k++)
   {
-    g=derivpol(g); g=gadd(gmul(p1,g),gmul2n(gmul(p2,derivpol(g)),1));
+    g = derivpol(g);
+    g = gadd(gmul(A,g), gmul2n(gmul(Bx,derivpol(g)),1));
   }
-  g = m ? gmul2n(g,(m-1)>>1):gmul2n(g,-1);
-  s=mulsi(n-m,mpfact(m+1));
-  tetpil=avma; return gerepile(av,tetpil,gdiv(g,s));
+  g = m? gmul2n(g,(m-1)>>1): gmul2n(g,-1);
+  s = mulsi(n-m, mpfact(m+1));
+  return gerepileupto(av, gdiv(g,s));
 }
 
 #ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */
@@ -1132,43 +1186,48 @@ polzag(long n, long m)
 GEN
 polzagreel(long n, long m, long prec)
 {
-  long d1,d,r,j,k,k2,av=avma,tetpil;
-  GEN p2,pol1,g,h,v,b,gend,s;
+  long d1, d, r, j, k, k2;
+  gpmem_t av = avma;
+  GEN Bx,B,g,h,v,b,s;
 
-  if (m>=n || m<0)
+  if (m >= n || m < 0)
     err(talker,"first index must be greater than second in polzag");
-  d1=n-m; d=d1<<1; d1--; pol1=gadd(gun,polx[0]);
-  p2=gmul(polx[0],pol1); r=(m+1)>>1; gend=stoi(d);
-  v=cgetg(d1+2,t_VEC); g=cgetg(d1+2,t_VEC);
-  v[d1+1]=un; b=mulri(realun(prec),gend); g[d1+1]=(long)b;
+  d1 = n-m; d = d1<<1; d1--;
+  B = gadd(gun,polx[0]);
+  Bx = gmul(polx[0],B);
+  r = (m+1)>>1;
+  v = cgetg(d1+2,t_VEC);
+  g = cgetg(d1+2,t_VEC);
+  v[d1+1] = un; b = stor(d, prec);
+  g[d1+1] = (long)b;
   for (k=1; k<=d1; k++)
   {
-    v[d1-k+1]=un;
+    v[d1-k+1] = un;
     for (j=1; j<k; j++)
-      v[d1-k+j+1]=(long)addii((GEN)v[d1-k+j+1],(GEN)v[d1-k+j+2]);
-    k2=k+k; b=divri(mulri(b,mulss(d-k2+1,d-k2)),mulss(k2,k2+1));
+      v[d1-k+j+1] = laddii((GEN)v[d1-k+j+1], (GEN)v[d1-k+j+2]);
+    k2 = k+k; b = divri(mulri(b,mulss(d-k2+1,d-k2)), mulss(k2,k2+1));
     for (j=1; j<=k; j++)
-      g[d1-k+j+1]=(long)mpadd((GEN)g[d1-k+j+1],mulri(b,(GEN)v[d1-k+j+1]));
-    g[d1-k+1]=(long)b;
+      g[d1-k+j+1] = lmpadd((GEN)g[d1-k+j+1], mulri(b,(GEN)v[d1-k+j+1]));
+    g[d1-k+1] = (long)b;
   }
-  h=cgetg(d1+3,t_POL); h[1]=evalsigne(1)+evallgef(d1+3);
-  for (k=0; k<=d1; k++) h[k+2]=g[k+1];
-  g=gmul(h,gpuigs(p2,r));
+  g = gmul(vec_to_pol(g,0), gpowgs(Bx,r));
   for (j=0; j<=r; j++)
   {
-    if (j) g=derivpol(g);
+    if (j) g = derivpol(g);
     if (j || !(m&1))
     {
-      h=cgetg(n+3,t_POL); h[1]=evalsigne(1)+evallgef(n+3);
-      h[2]=g[2];
+      h = cgetg(n+3,t_POL);
+      h[1] = evalsigne(1)|evallgef(n+3);
+      h[2] = g[2];
       for (k=1; k<n; k++)
-	h[k+2]=ladd(gmulsg(k+k+1,(GEN)g[k+2]),gmulsg(k<<1,(GEN)g[k+1]));
-      h[n+2]=lmulsg(n<<1,(GEN)g[n+1]); g=h;
+	h[k+2] = ladd(gmulsg(k+k+1,(GEN)g[k+2]), gmulsg(k<<1,(GEN)g[k+1]));
+      h[n+2] = lmulsg(n<<1, (GEN)g[n+1]);
+      g = h;
     }
   }
-  g = m ?gmul2n(g,(m-1)>>1):gmul2n(g,-1);
-  s=mulsi(n-m,mpfact(m+1));
-  tetpil=avma; return gerepile(av,tetpil,gdiv(g,s));
+  g = m? gmul2n(g,(m-1)>>1): gmul2n(g,-1);
+  s = mulsi(n-m, mpfact(m+1));
+  return gerepileupto(av, gdiv(g,s));
 }
 #ifdef _MSC_VER
 #pragma optimize("g",on)
@@ -1183,7 +1242,8 @@ polzagreel(long n, long m, long prec)
 GEN
 zbrent(entree *ep, GEN a, GEN b, char *ch, long prec)
 {
-  long av = avma,sig,iter,itmax;
+  long sig, iter, itmax;
+  gpmem_t av = avma;
   GEN c,d,e,tol,tol1,min1,min2,fa,fb,fc,p,q,r,s,xm;
 
   a = fix(a,prec);
@@ -1238,13 +1298,8 @@ zbrent(entree *ep, GEN a, GEN b, char *ch, long prec)
     else { d = xm; e = d; } /* bound decreasing too slowly, use bisection */
     a = b; fa = fb;
     if (gcmp(gabs(d,0),tol1) > 0) b = gadd(b,d);
-    else
-    {
-      if (gsigne(xm) > 0)
-        b = addrr(b,tol1);
-      else
-        b = subrr(b,tol1);
-    }
+    else if (gsigne(xm) > 0)      b = addrr(b,tol1);
+    else                          b = subrr(b,tol1);
     ep->value = (void*)b; fb = lisexpr(ch);
   }
   if (iter > itmax) err(talker,"too many iterations in solve");