===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/rootpol.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/rootpol.c	2001/10/02 11:17:05	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/rootpol.c	2002/09/11 07:26:52	1.2
@@ -1,4 +1,4 @@
-/* $Id: rootpol.c,v 1.1 2001/10/02 11:17:05 noro Exp $
+/* $Id: rootpol.c,v 1.2 2002/09/11 07:26:52 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -22,8 +22,6 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #include "pari.h"
 
 extern GEN polrecip_i(GEN x);
-extern GEN poldeflate(GEN x0, long *m);
-extern GEN roots_to_pol(GEN a, long v);
 #define pariINFINITY 100000
 #define NEWTON_MAX 10
 
@@ -65,7 +63,7 @@ quickmulcc(GEN x, GEN y)
     }
     if (ty==t_COMPLEX)
     {
-      long av,tetpil;
+      gpmem_t av, tetpil;
       GEN p1,p2;
 
       z=cgetg(3,t_COMPLEX); av=avma;
@@ -97,7 +95,8 @@ static GEN
 mysquare(GEN p)
 {
   GEN s,aux1,aux2;
-  long i,j,n=degpol(p),nn,ltop,lbot;
+  long i, j, n=degpol(p), nn;
+  gpmem_t ltop, lbot;
 
   if (n==-1) return gcopy(p);
   nn=n<<1; s=cgetg(nn+3,t_POL);
@@ -142,7 +141,7 @@ karasquare(GEN p)
 {
   GEN p1,s0,s1,s2,aux;
   long n=degpol(p),n0,n1,i,var,nn0;
-  ulong ltop;
+  gpmem_t ltop;
 
   if (n<=KARASQUARE_LIMIT) return mysquare(p);
   ltop=avma;
@@ -174,7 +173,7 @@ cook_square(GEN p)
 {
   GEN p0,p1,p2,p3,q,aux0,aux1,r,aux,plus,moins;
   long n=degpol(p),n0,n3,i,j,var;
-  ulong ltop = avma;
+  gpmem_t ltop = avma;
 
   if (n<=COOK_SQUARE_LIMIT) return karasquare(p);
 
@@ -347,7 +346,7 @@ log2ir(GEN x)
   return 1.+ (double) expo(x)+log2( (double)(ulong) x[2]) - (double) BITS_IN_LONG;
 }
 
-static double
+double
 mylog2(GEN z)
 {
   double x,y;
@@ -360,58 +359,59 @@ mylog2(GEN z)
   return x+0.5*log2( 1 + exp2(2*(y-x)));
 }
 
+/* find s such that  A_h <= 2^s <= 2 A_i  for one h and all i < n = deg(p),
+ * with  A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and  p = sum p_i X^i */
 static long
 findpower(GEN p)
 {
-  double x, logbinomial,pente,pentemax=-pariINFINITY;
-  long n=degpol(p),i;
+  double x, logbinomial, mins = pariINFINITY;
+  long n = degpol(p),i;
 
-  logbinomial = mylog2((GEN) p[n+2]);
+  logbinomial = mylog2((GEN)p[n+2]); /* log2(lc * binom(n,i)) */
   for (i=n-1; i>=0; i--)
   {
     logbinomial += log2((double) (i+1) / (double) (n-i));
-    x = mylog2((GEN) p[2+i])-logbinomial;
-    if (x>-pariINFINITY)
+    x = mylog2((GEN)p[i+2]);
+    if (x != -pariINFINITY)
     {
-      pente = x/ (double) (n-i);
-      if (pente > pentemax) pentemax = pente;
+      double s = (logbinomial - x) / (double) (n-i);
+      if (s < mins) mins = s;
     }
   }
-  return (long) -floor(pentemax);
+  i = (long)ceil(mins);
+  if (i - mins > 1 - 1e-12) i--;
+  return i;
 }
 
 /* returns the exponent for the procedure modulus, from the newton diagram */
 static long
-polygone_newton(GEN p, long k)
+newton_polygon(GEN p, long k)
 {
   double *logcoef,pente;
-  long n=degpol(p),i,j,h,l,*sommet,pentelong;
+  long n=degpol(p),i,j,h,l,*vertex,pentelong;
 
-  logcoef=(double*) gpmalloc((n+1)*sizeof(double));
-  sommet=(long*) gpmalloc((n+1)*sizeof(long));
+  logcoef = (double*) gpmalloc((n+1)*sizeof(double));
+  vertex  = (long*)   gpmalloc((n+1)*sizeof(long));
 
-  /* sommet[i]=1 si i est un sommet, =0 sinon */
-  for (i=0; i<=n; i++) { logcoef[i]=mylog2((GEN)p[2+i]); sommet[i]=0; }
-  sommet[0]=1; i=0;
-  while (i<n)
+  /* vertex[i]=1 if i a vertex of convex hull, 0 otherwise */
+  for (i=0; i<=n; i++) { logcoef[i] = mylog2((GEN)p[2+i]); vertex[i] = 0; }
+  vertex[0]=1;
+  for (i=0; i < n; i=h)
   {
-    pente=logcoef[i+1]-logcoef[i];
-    h=i+1;
+    pente = logcoef[i+1]-logcoef[i];
+    h = i+1;
     for (j=i+1; j<=n; j++)
-    {
-      if (pente<(logcoef[j]-logcoef[i])/(double)(j-i))
+      if (pente < (logcoef[j]-logcoef[i])/(double)(j-i))
       {
-	h=j;
-	pente=(logcoef[j]-logcoef[i])/(double)(j-i);
+	h = j;
+	pente = (logcoef[j]-logcoef[i])/(double)(j-i);
       }
-    }
-    i=h;
-    sommet[h]=1;
+    vertex[h] = 1;
   }
-  h=k; while (!sommet[h]) h++;
-  l=k-1; while (!sommet[l]) l--;
-  pentelong=(long) floor((logcoef[h]-logcoef[l])/(double)(h-l)+0.5);
-  free(logcoef); free(sommet); return pentelong;
+  h = k;   while (!vertex[h]) h++;
+  l = k-1; while (!vertex[l]) l--;
+  pentelong = (long) floor((logcoef[h]-logcoef[l])/(double)(h-l) + 0.5);
+  free(logcoef); free(vertex); return pentelong;
 }
 
 /* change z into z*2^e, where z is real or complex of real */
@@ -444,10 +444,8 @@ myshiftic(GEN z, long e)
 static GEN
 myrealun(long bitprec)
 {
-  GEN x;
   if (bitprec < 0) bitprec = 0;
-  x = cgetr(bitprec/BITS_IN_LONG+3);
-  affsr(1,x); return x;
+  return realun(bitprec/BITS_IN_LONG+3);
 }
 
 static GEN
@@ -547,7 +545,7 @@ homothetie(GEN p, GEN R, long bitprec)
   return r;
 }
 
-/* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) */
+/* change q in 2^(n*e) p(x*2^(-e)), n=deg(q)  [ ~as above with R = 2^-e ]*/
 static void
 homothetie2n(GEN p, long e)
 {
@@ -569,98 +567,92 @@ homothetie_gauss(GEN p, long e,long f)
   }
 }
 
-static long
-valuation(GEN p)
-{
-  long j=0,n=degpol(p);
-
-  while ((j<=n) && isexactzero((GEN)p[j+2])) j++;
-  return j;
-}
-
 /* provides usually a good lower bound on the largest modulus of the roots,
 puts in k an upper bound of the number of roots near the largest root
 at a distance eps */
 static double
 lower_bound(GEN p, long *k, double eps)
 {
-  long n=degpol(p),i,j,ltop=avma;
-  GEN a,s,icd;
+  long n=degpol(p), i, j;
+  gpmem_t ltop=avma;
+  GEN a,s,S,icd;
   double r,*rho;
 
   if (n<4) { *k=n; return 0.; }
-  a=cgetg(6,t_POL); s=cgetg(6,t_POL);
+  S = cgetg(5,t_VEC);
+  a = cgetg(5,t_VEC);
   rho=(double *) gpmalloc(4*sizeof(double));
   icd = gdiv(realun(DEFAULTPREC), (GEN) p[n+2]);
-  for (i=1; i<=4; i++) a[i+1]=lmul(icd,(GEN)p[n+2-i]);
+  for (i=1; i<=4; i++) a[i] = lmul(icd,(GEN)p[n+2-i]);
   for (i=1; i<=4; i++)
   {
-    s[i+1]=lmulsg(i,(GEN)a[i+1]);
+    s = gmulsg(i,(GEN)a[i]);
     for (j=1; j<i; j++)
-      s[i+1]=ladd((GEN)s[i+1],gmul((GEN)s[j+1],(GEN)a[i+1-j]));
-    s[i+1]=lneg((GEN)s[i+1]);
-    r=gtodouble(gabs((GEN) s[i+1],3));
-    if (r<=0.)  /* should not be strictly negative */
-      rho[i-1]=0.;
+      s = gadd(s, gmul((GEN)S[j],(GEN)a[i-j]));
+    S[i] = lneg(s); /* Newton sum S_i */
+    r = gtodouble(gabs(s,3));
+    if (r == 0.)
+      rho[i-1] = 0.;
     else
-      rho[i-1]=exp(log(r/(double)n)/(double) i);
+      rho[i-1] = exp(log(r/(double)n)/(double) i);
   }
-  r=0.;
+  r = 0.;
   for (i=0; i<4; i++) if (r<rho[i]) r=rho[i];
   if (r>0. && eps<1.2)
-    *k=(long) floor((n*rho[0]/r+n)/(1+exp(-eps)*cos(eps)));
+    *k = (long) floor((n*rho[0]/r+n)/(1+exp(-eps)*cos(eps)));
   else
-    *k=n;
-  free(rho); avma=ltop; return r;
+    *k = n;
+  free(rho); avma = ltop; return r;
 }
 
 /* returns the maximum of the modulus of p with a relative error tau */
-static GEN
+GEN
 max_modulus(GEN p, double tau)
 {
-  GEN q,aux,gunr;
-  long i,j,k,valuat,n=degpol(p),nn,ltop=avma,bitprec,imax,e;
-  double r,rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
+  GEN r,q,aux,gunr;
+  gpmem_t av, ltop = avma;
+  long i,k,n=degpol(p),nn,bitprec,M,e;
+  double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
 
+  r = cgeti(BIGDEFAULTPREC);
+  av = avma;
+
   eps = - 1/log(1.5*tau2); /* > 0 */
-  bitprec=(long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
-  gunr=myrealun(bitprec+2*n);
-  aux=gdiv(gunr,(GEN) p[2+n]);
-  q=gmul(aux,p); q[2+n]=lcopy(gunr);
-  k=nn=n;
-  e=findpower(q); homothetie2n(q,e); r=-(double) e;
-  q=mygprec(q,bitprec+(n<<1));
+  bitprec = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
+  gunr = myrealun(bitprec+2*n);
+  aux = gdiv(gunr,(GEN) p[2+n]);
+  q = gmul(aux,p); q[2+n] = (long)gunr;
+  e = findpower(q);
+  homothetie2n(q,e);
+  affsi(e, r);
+  q = mygprec(q,bitprec+(n<<1));
   pol_to_gaussint(q,bitprec);
-  imax=(long) ((log(log(4.*n)/(2*tau2))) / log(2.)) + 2;
+  M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
+  nn = n;
   for (i=0,e=0;;)
-  {
-    rho=lower_bound(q,&k,eps);
-    if (rho>exp2(-(double) e)) e = (long) -floor(log2(rho));
-    r -= e / exp2((double)i);
-    if (++i == imax) {
-      avma=ltop; 
-      return gpui(dbltor(2.),dbltor(r),DEFAULTPREC);
-    }
+  { /* nn = deg(q) */
+    rho = lower_bound(q,&k,eps);
+    if (rho > exp2(-(double) e)) e = (long) -floor(log2(rho));
+    affii(shifti(addis(r,e), 1), r);
+    if (++i == M) break;
 
-    if (k<nn)
-      bitprec=(long) ((double) k* log2(1./tau2)+
-                      (double) (nn-k)*log2(1./eps)+
-                      3*log2((double) nn))+1;
-    else
-      bitprec=(long) ((double) nn* log2(1./tau2)+
-                      3.*log2((double) nn))+1;
+    bitprec = (long) ((double) k* log2(1./tau2) +
+                     (double) (nn-k)*log2(1./eps) + 3*log2((double) nn)) + 1;
     homothetie_gauss(q,e,bitprec-(long)floor(mylog2((GEN) q[2+nn])+0.5));
-    valuat=valuation(q);
-    if (valuat>0)
-    {
-      nn-=valuat; setlgef(q,nn+3);
-      for (j=0; j<=nn; j++) q[2+j]=q[2+valuat+j];
-    }
+    nn -= polvaluation(q, &q);
     set_karasquare_limit(gexpo(q));
-    q = gerepileupto(ltop, graeffe(q));
-    tau2=1.5*tau2; eps=1/log(1./tau2);
+    q = gerepileupto(av, graeffe(q));
+    tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
+    eps = -1/log(tau2); /* > 0 */
     e = findpower(q);
   }
+  if (!signe(r)) { avma = ltop; return realun(DEFAULTPREC); } 
+  r = itor(r, DEFAULTPREC);
+  setexpo(r, expo(r) - M);
+  rho = rtodbl(r);
+  /* rho = sum e_i 2^-i */
+  avma = ltop;
+  return mpexp(dbltor(-rho * LOG2)); /* 2^-r */
 }
 
 /* return the k-th modulus (in ascending order) of p, rel. error tau*/
@@ -668,41 +660,35 @@ static GEN
 modulus(GEN p, long k, double tau)
 {
   GEN q,gunr;
-  long i,j,kk=k,imax,n=degpol(p),nn,nnn,valuat,av,ltop=avma,bitprec,decprec,e;
+  long i, kk=k, imax, n=degpol(p), nn, bitprec, decprec, e;
+  gpmem_t av, ltop=avma;
   double tau2,r;
 
-  tau2=tau/6; nn=n;
+  tau2 = tau/6;
   bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2)));
   decprec=(long) ((double) bitprec * L2SL10)+1;
   gunr=myrealun(bitprec);
   av = avma;
-  q=gprec(p,decprec);
-  q=gmul(gunr,q);
-  e=polygone_newton(q,k);
+  q = gprec(p,decprec);
+  q = gmul(gunr,q);
+  e = newton_polygon(q,k);
   homothetie2n(q,e);
-  r=(double) e;
+  r = (double) e;
   imax=(long) ((log2(3./tau)+log2(log(4.*(double) n)) ))+1;
   for (i=1; i<imax; i++)
   {
-    q=eval_rel_pol(q,bitprec);
+    q = eval_rel_pol(q,bitprec);
+    kk -= polvaluation(q, &q);
+    nn = degpol(q);
 
-    nnn=degpol(q); valuat=valuation(q);
-    if (valuat>0)
-    {
-      kk-=valuat;
-      for (j=0; j<=nnn-valuat; j++) q[2+j]=q[2+valuat+j];
-      setlgef(q,nnn-valuat+3);
-    }
-    nn=nnn-valuat;
-
     set_karasquare_limit(bitprec);
     q = gerepileupto(av, graeffe(q));
-    e=polygone_newton(q,kk);
+    e = newton_polygon(q,kk);
     r += e / exp2((double)i);
     q=gmul(gunr,q);
     homothetie2n(q,e);
 
-    tau2=1.5*tau2; if (tau2>1.) tau2=1.;
+    tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
     bitprec= 1+(long) ((double) nn*(2.+log2(3.*(double) nn)+log2(1./tau2)));
   }
   avma=ltop; return mpexp(dbltor(-r * LOG2));
@@ -715,7 +701,8 @@ static GEN
 pre_modulus(GEN p, long k, double tau, GEN rmin, GEN rmax)
 {
   GEN R, q, aux;
-  long n=degpol(p),i,imax,imax2,bitprec,ltop=avma, av;
+  long n=degpol(p), i, imax, imax2, bitprec;
+  gpmem_t ltop=avma, av;
   double tau2, aux2;
 
   tau2=tau/6.;
@@ -750,7 +737,7 @@ pre_modulus(GEN p, long k, double tau, GEN rmin, GEN r
 static GEN
 min_modulus(GEN p, double tau)
 {
-  long av=avma;
+  gpmem_t av=avma;
   GEN r;
 
   if (isexactzero((GEN)p[2])) return realzero(3);
@@ -764,7 +751,8 @@ static long
 dual_modulus(GEN p, GEN R, double tau, long l)
 {
   GEN q;
-  long i,j,imax,k,delta_k=0,n=degpol(p),nn,nnn,valuat,ltop=avma,bitprec,ll=l;
+  long i, imax, k, delta_k=0, n=degpol(p), nn, nnn, valuat, bitprec, ll=l;
+  gpmem_t ltop=avma;
   double logmax,aux,tau2;
 
   tau2=7.*tau/8.;
@@ -778,24 +766,18 @@ dual_modulus(GEN p, GEN R, double tau, long l)
     bitprec=6*nn-5*ll+(long) ((double) nn*(log2(1/tau2)+8.*tau2/7.));
 
     q=eval_rel_pol(q,bitprec);
-    nnn=degpol(q); valuat=valuation(q);
-    if (valuat>0)
-    {
-      delta_k+=valuat;
-      for (j=0; j<=nnn-valuat; j++) q[2+j]=q[2+valuat+j];
-      setlgef(q,nnn-valuat+3);
-    }
+    nnn = degpol(q); valuat = polvaluation(q, &q);
     ll= (-valuat<nnn-n)? ll-valuat: ll+nnn-n; /* min(ll-valuat,ll+nnn-n) */
     if (ll<0) ll=0;
 
-    nn=nnn-valuat;
+    nn = nnn-valuat; delta_k += valuat;
     if (nn==0) return delta_k;
 
     set_karasquare_limit(bitprec);
     q = gerepileupto(ltop, graeffe(q));
     tau2=tau2*7./4.;
   }
-  k=-1; logmax=- (double) pariINFINITY;
+  k = -1; logmax = - (double)pariINFINITY;
   for (i=0; i<=degpol(q); i++)
   {
     aux=mylog2((GEN)q[2+i]);
@@ -831,7 +813,7 @@ gmulbyi(GEN z)
 static void
 fft(GEN Omega, GEN p, GEN f, long Step, long l)
 {
-  ulong ltop;
+  gpmem_t ltop;
   long i,l1,l2,l3,rap=Lmax/l,rapi,Step4;
   GEN f1,f2,f3,f02,f13,g02,g13,ff;
 
@@ -885,20 +867,6 @@ fft(GEN Omega, GEN p, GEN f, long Step, long l)
   for (i=0; i<l; i++) f[i]=ff[i+1];
 }
 
-extern void mpsincos(GEN x, GEN *s, GEN *c);
-
-/* return exp(ix), x a t_REAL */
-static GEN
-exp_i(GEN x)
-{
-  GEN v;
-  
-  if (!signe(x)) return realun(lg(x)); /* should not happen */
-  v = cgetg(3,t_COMPLEX);
-  mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
-  return v;
-}
-
 /* e(1/N) */
 static GEN
 RUgen(long N, long bitprec)
@@ -907,7 +875,7 @@ RUgen(long N, long bitprec)
   if (N == 2) return mpneg(realun(bitprec));
   if (N == 4) return gi;
   pi2 = gmul2n(mppi(bitprec/BITS_IN_LONG+3), 1);
-  return exp_i(gdivgs(pi2,N));
+  return exp_Ir(gdivgs(pi2,N));
 }
 
 /* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
@@ -947,7 +915,7 @@ initRUgen(long N, long bitprec)
 }
 
 /* returns 1 if p has only real coefficients, 0 else */
-static long
+static int
 isreal(GEN p)
 {
   long n=degpol(p),i=0;
@@ -958,10 +926,11 @@ isreal(GEN p)
 
 static void
 parameters(GEN p, double *mu, double *gamma,
-           long polreal, double param, double param2)
+           int polreal, double param, double param2)
 {
   GEN q,pc,Omega,coef,RU,prim,aux,aux0,ggamma,gx,mygpi;
-  long ltop=avma,limite=stack_lim(ltop,1),n=degpol(p),bitprec,NN,K,i,j,ltop2;
+  long n=degpol(p), bitprec, NN, K, i, j;
+  gpmem_t ltop2, ltop=avma, limite=stack_lim(ltop, 1);
   double lx;
 
   bitprec=gexpo(p)+(long)param2+8;
@@ -969,7 +938,7 @@ parameters(GEN p, double *mu, double *gamma,
   K=NN/Lmax; if (K%2==1) K++; NN=Lmax*K;
   mygpi=mppi(bitprec/BITS_IN_LONG+3);
   aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */
-  prim = exp_i(aux);
+  prim = exp_Ir(aux);
   aux = gmulbyi(aux);
   RU = myrealun(bitprec);
 
@@ -1026,12 +995,13 @@ parameters(GEN p, double *mu, double *gamma,
 static void
 dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H, long polreal)
 {
-  GEN Omega,q,qd,pc,pdc,alpha,beta,gamma,RU,aux,U,W,mygpi,prim,prim2;
-  long limite,n=degpol(p),i,j,K,ltop;
+  GEN Omega,q,qd,pc,pdc,alpha,beta,gamma,RU,aux,U,W,mygpi,prim,prim2,*gptr[3];
+  long n=degpol(p),i,j,K;
+  gpmem_t ltop;
 
   mygpi=mppi(bitprec/BITS_IN_LONG+3);
   aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */
-  prim = exp_i(aux);
+  prim = exp_Ir(aux);
   aux = gmulbyi(aux);
   prim2 = myrealun(bitprec);
   RU=cgetg(n+2,t_VEC); RU++;
@@ -1051,7 +1021,7 @@ dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H
 
   if (polreal) K=K/2+1;
 
-  ltop=avma; limite = stack_lim(ltop,1);
+  ltop=avma;
   W=cgetg(k+1,t_VEC); U=cgetg(k+1,t_VEC);
   for (i=1; i<=k; i++) W[i]=U[i]=zero;
 
@@ -1102,13 +1072,8 @@ dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H
       }
     }
     prim2=gmul(prim2,prim);
-    if (low_stack(limite, stack_lim(ltop,1)))
-    {
-      GEN *gptr[3];
-      if(DEBUGMEM>1) err(warnmem,"dft");
-      gptr[0]=&W; gptr[1]=&U; gptr[2]=&prim2;
-      gerepilemany(ltop,gptr,3);
-    }
+    gptr[0]=&W; gptr[1]=&U; gptr[2]=&prim2;
+    gerepilemany(ltop,gptr,3);
   }
 
   for (i=1; i<=k; i++)
@@ -1129,7 +1094,7 @@ static GEN
 refine_H(GEN F, GEN G, GEN HH, long bitprec, long shiftbitprec)
 {
   GEN H=HH,D,aux;
-  ulong ltop=avma, limite=stack_lim(ltop,1);
+  gpmem_t ltop=avma, limite=stack_lim(ltop, 1);
   long error=0,i,bitprec1,bitprec2;
 
   D=gsub(gun,gres(gmul(HH,G),F)); error=gexpo(D);
@@ -1139,9 +1104,9 @@ refine_H(GEN F, GEN G, GEN HH, long bitprec, long shif
   {
     if (low_stack(limite, stack_lim(ltop,1)))
     {
-      GEN *gptr[2];
+      GEN *gptr[2]; gptr[0]=&D; gptr[1]=&H;
       if(DEBUGMEM>1) err(warnmem,"refine_H");
-      gptr[0]=&D; gptr[1]=&H; gerepilemany(ltop,gptr,2);
+      gerepilemany(ltop,gptr,2);
     }
     bitprec1=-error+shiftbitprec;
     aux=gmul(mygprec(H,bitprec1),mygprec(D,bitprec1));
@@ -1154,8 +1119,8 @@ refine_H(GEN F, GEN G, GEN HH, long bitprec, long shif
     D=gsub(gun,gres(gmul(H,G),F));
     error=gexpo(D); if (error<-bitprec1) error=-bitprec1;
   }
-  if (error<=-bitprec/2) return gerepilecopy(ltop,H);
-  avma=ltop; return gzero; /* procedure failed */
+  if (error>-bitprec/2) return NULL; /* FAIL */
+  return gerepilecopy(ltop,H);
 }
 
 /* return 0 if fails, 1 else */
@@ -1163,8 +1128,9 @@ static long
 refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, double gamma)
 {
   GEN pp,FF,GG,r,HH,f0;
-  long error,i,bitprec1=0,bitprec2,ltop=avma,shiftbitprec;
-  long shiftbitprec2,n=degpol(p),enh,normF,normG,limite=stack_lim(ltop,1);
+  long error, i, bitprec1=0, bitprec2, shiftbitprec;
+  long shiftbitprec2, n=degpol(p), enh, normF, normG;
+  gpmem_t ltop=avma, limite=stack_lim(ltop, 1);
 
   FF=*F; HH=H;
   GG=poldivres(p,*F,&r);
@@ -1184,16 +1150,15 @@ refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, d
     }
     if (low_stack(limite, stack_lim(ltop,1)))
     {
-      GEN *gptr[4];
+      GEN *gptr[4]; gptr[0]=&FF; gptr[1]=&GG; gptr[2]=&r; gptr[3]=&HH;
       if(DEBUGMEM>1) err(warnmem,"refine_F");
-      gptr[0]=&FF; gptr[1]=&GG; gptr[2]=&r; gptr[3]=&HH;
       gerepilemany(ltop,gptr,4);
     }
 
     bitprec1=-error+shiftbitprec2;
     HH=refine_H(mygprec(FF,bitprec1),mygprec(GG,bitprec1),
 		mygprec(HH,bitprec1),1-error,shiftbitprec2);
-    if (HH==gzero) return 0; /* procedure failed */
+    if (!HH) return 0; /* procedure failed */
 
     bitprec1=-error+shiftbitprec;
     r=gmul(mygprec(HH,bitprec1),mygprec(r,bitprec1));
@@ -1210,12 +1175,8 @@ refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, d
     GG=poldivres(pp,mygprec(FF,bitprec1),&r);
     error=gexpo(r); if (error<-bitprec1) error=-bitprec1;
   }
-  if (error<=-bitprec)
-  {
-    *F=FF; *G=GG;
-    return 1; /* procedure succeeds */
-  }
-  return 0; /* procedure failed */
+  if (error>-bitprec) return 0; /* FAIL */
+  *F=FF; *G=GG; return 1;
 }
 
 /* returns F and G from the unit circle U such that |p-FG|<2^(-bitprec) |cd|,
@@ -1225,8 +1186,9 @@ split_fromU(GEN p, long k, double delta, long bitprec,
             GEN *F, GEN *G, double param, double param2)
 {
   GEN pp,FF,GG,H;
-  long n=degpol(p),NN,bitprec2,
-  ltop=avma,polreal=isreal(p);
+  long n=degpol(p),NN,bitprec2;
+  int polreal=isreal(p);
+  gpmem_t ltop;
   double mu,gamma;
 
   pp=gdiv(p,(GEN)p[2+n]);
@@ -1244,7 +1206,7 @@ split_fromU(GEN p, long k, double delta, long bitprec,
     bitprec2=(long) (((double) NN*delta-mu)/log(2.))+gexpo(pp)+8;
     dft(pp,k,NN,bitprec2,FF,H,polreal);
     if (refine_F(pp,&FF,&GG,H,bitprec,gamma)) break;
-    NN=(NN<<1); avma=ltop;
+    NN <<= 1; avma=ltop;
   }
   *G=gmul(GG,(GEN)p[2+n]); *F=FF;
 }
@@ -1301,7 +1263,8 @@ extern GEN addshiftpol(GEN x, GEN y, long d);
 static GEN
 shiftpol(GEN p, GEN b)
 {
-  long av = avma,i, limit = stack_lim(av,1);
+  long i;
+  gpmem_t av = avma, limit = stack_lim(av, 1);
   GEN q = gzero;
 
   if (gcmp0(b)) return p;
@@ -1326,7 +1289,7 @@ conformal_pol(GEN p, GEN a, long bitprec)
 {
   GEN r,pui,num,aux, unr = myrealun(bitprec);
   long n=degpol(p), i;
-  ulong av, limit;
+  gpmem_t av, limit;
 
   aux = pui = cgetg(4,t_POL);
   pui[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(4);
@@ -1388,7 +1351,7 @@ compute_radius(GEN* radii, GEN p, long k, double aux, 
   {
     if (!signe(radii[i]) || cmprr(radii[i], p1) < 0)
       affrr(p1, radii[i]);
-    else      
+    else
       p1 = radii[i];
   }
   *delta = rtodbl(gmul2n(mplog(divrr(rmax,rmin)), -1));
@@ -1417,7 +1380,8 @@ static void
 conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, long bitprec,
                   double aux, GEN *F,GEN *G)
 {
-  long bitprec2,n=degpol(p),decprec,i,ltop = avma, av;
+  long bitprec2, n=degpol(p), decprec, i;
+  gpmem_t ltop = avma, av;
   GEN q,FF,GG,a,R, *gptr[2];
   GEN rho,invrho;
   double delta,param,param2;
@@ -1432,7 +1396,7 @@ conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, 
   for (i=1; i<=n; i++)
     if (signe(radii[i])) /* updating array radii */
     {
-      long a = avma;
+      gpmem_t a = avma;
       GEN p1 = gsqr(radii[i]);
       /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
       p1 = divrr(gmul2n((subrs(p1,1)),1),
@@ -1443,7 +1407,7 @@ conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, 
 
   rho = compute_radius(radii, q,k,aux/10., &delta);
   invrho = update_radius(radii, rho, &param, &param2);
-  
+
   bitprec2 += (long) (((double)n) * fabs(log2ir(rho)) + 1.);
   R = mygprec(invrho,bitprec2);
   q = scalepol(q,R,bitprec2);
@@ -1470,6 +1434,14 @@ conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, 
   gptr[0]=F; gptr[1]=G; gerepilemany(ltop,gptr,2);
 }
 
+static GEN
+myrealzero(void)
+{
+  GEN x = cgetr(3);
+  x[1] = evalexpo(-bit_accuracy(3));
+  return x;
+}
+
 /* split p, this time with no scaling. returns in F and G two polynomials
 such that |p-FG|< 2^(-bitprec)|p| */
 static void
@@ -1480,7 +1452,7 @@ split_2(GEN p, long bitprec, GEN ctr, double thickness
   long n=degpol(p),i,j,k,bitprec2;
   GEN q,FF,GG,R;
   GEN *radii = (GEN*) cgetg(n+1, t_VEC);
-  for (i=2; i<n; i++) radii[i]=realzero(3);
+  for (i=2; i<n; i++) radii[i] = myrealzero();
   aux = thickness/(double) n/4.;
   radii[1] = rmin = min_modulus(p, aux);
   radii[n] = rmax = max_modulus(p, aux);
@@ -1752,8 +1724,10 @@ mygprec_absolute(GEN x, long bitprec)
   {
     case t_REAL:
       e=gexpo(x);
-      if (e<-bitprec || !signe(x)) { y=dbltor(0.); setexpo(y,-bitprec); }
-      else y=mygprec(x,bitprec+e);
+      if (e<-bitprec || !signe(x))
+        y=realzero_bit(-bitprec);
+      else
+        y=mygprec(x,bitprec+e);
       break;
     case t_COMPLEX:
       if (gexpo((GEN)x[2])<-bitprec)
@@ -1781,7 +1755,7 @@ a_posteriori_errors(GEN p, GEN roots_pol, long err)
   setexpo(sigma, err + (long)log2((double)n) + 1);
   overn=dbltor(1./n);
   shatzle=gdiv(gpui(sigma,overn,0),
-	       gsub(gpui(gsub(gun,sigma),overn,0),
+	       gsub(gpui(gsub(realun(DEFAULTPREC),sigma),overn,0),
 		    gpui(sigma,overn,0)));
   shatzle=gmul2n(shatzle,1); e_max=-pariINFINITY;
   for (i=1; i<=n; i++)
@@ -1799,24 +1773,21 @@ a_posteriori_errors(GEN p, GEN roots_pol, long err)
 /**                                                                **/
 /********************************************************************/
 static GEN
-append_root(GEN roots_pol, GEN a)
-{
-  long l = lg(roots_pol);
-  setlg(roots_pol, l+1); return (GEN)(roots_pol[l] = lclone(a));
-}
+append_clone(GEN r, GEN a) { a = gclone(a); appendL(r, a); return a; }
 
 /* put roots in placeholder roots_pol so that |P-L_1...L_n|<2^(-bitprec)|P|
  *  and returns  prod (x-roots_pol[i]) for i=1..degree(p) */
 static GEN
 split_complete(GEN p, long bitprec, GEN roots_pol)
 {
-  long n=degpol(p),decprec,ltop;
+  long n=degpol(p), decprec;
+  gpmem_t ltop;
   GEN p1,F,G,a,b,m1,m2,m;
 
   if (n==1)
   {
     a=gneg_i(gdiv((GEN)p[2],(GEN)p[3]));
-    append_root(roots_pol,a); return p;
+    (void)append_clone(roots_pol,a); return p;
   }
   ltop = avma;
   if (n==2)
@@ -1827,8 +1798,8 @@ split_complete(GEN p, long bitprec, GEN roots_pol)
     p1 = gmul2n((GEN)p[4],1);
     a = gneg_i(gdiv(gadd(F,(GEN)p[3]), p1));
     b =        gdiv(gsub(F,(GEN)p[3]), p1);
-    a = append_root(roots_pol,a);
-    b = append_root(roots_pol,b); avma = ltop;
+    a = append_clone(roots_pol,a);
+    b = append_clone(roots_pol,b); avma = ltop;
     m=gmul(gsub(polx[varn(p)],mygprec(a,3*bitprec)),
 	   gsub(polx[varn(p)],mygprec(b,3*bitprec)));
     return gmul(m,(GEN)p[4]);
@@ -1840,21 +1811,24 @@ split_complete(GEN p, long bitprec, GEN roots_pol)
 }
 
 /* compute a bound on the maximum modulus of roots of p */
-static GEN
+GEN
 cauchy_bound(GEN p)
 {
-  long i,n=degpol(p);
-  GEN x=gzero,y,lc;
+  long i, n = degpol(p), prec = DEFAULTPREC;
+  GEN lc, y, q = gmul(p, realun(prec)), x = gzero;
 
-  lc=gabs((GEN)p[n+2],DEFAULTPREC); /* leading coefficient */
-  lc=gdiv(dbltor(1.),lc);
+  if (n <= 0) err(constpoler,"cauchy_bound");
+
+  lc = gabs((GEN)q[n+2],prec); /* leading coefficient */
+  lc = ginv(lc);
   for (i=0; i<n; i++)
   {
-    y=gmul(gabs((GEN) p[i+2],DEFAULTPREC),lc);
-    y=gpui(y,dbltor(1./(n-i)),DEFAULTPREC);
-    if (gcmp(y,x) > 0) x=y;
+    y = (GEN)q[i+2]; if (gcmp0(y)) continue;
+    y = gmul(gabs(y,prec), lc);
+    y = divrs(mplog(y), n-i);
+    if (gcmp(y,x) > 0) x = y;
   }
-  return x;
+  return gexp(x, prec);
 }
 
 static GEN
@@ -1925,32 +1899,30 @@ fix_roots(GEN r, GEN *m, long h, long bitprec)
 static GEN
 all_roots(GEN p, long bitprec)
 {
-  GEN pd,q,roots_pol,m;
+  GEN lc,pd,q,roots_pol,m;
   long bitprec0, bitprec2,n=degpol(p),i,e,h;
-  ulong av;
+  gpmem_t av;
 
-#if 0
-  pd = poldeflate(p, &h);
-#else
-  pd = p; h = 1;
-#endif
+  pd = poldeflate(p, &h); lc = leading_term(pd);
   e = 2*gexpo(cauchy_bound(pd)); if (e<0) e=0;
-  bitprec0=bitprec + gexpo(pd) - gexpo(leading_term(pd)) + (long)log2(n/h)+1+e;
+  bitprec0=bitprec + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
+  bitprec2 = bitprec0; e = 0;
   for (av=avma,i=1;; i++,avma=av)
   {
-    roots_pol = cgetg(n+1,t_VEC); setlg(roots_pol,1); 
-    bitprec2 = bitprec0 + (1<<i)*n;
+    roots_pol = cget1(n+1,t_VEC);
+    bitprec2 += e + (1<<i)*n;
     q = gmul(myrealun(bitprec2), mygprec(pd,bitprec2));
     m = split_complete(q,bitprec2,roots_pol);
     roots_pol = fix_roots(roots_pol, &m, h, bitprec2);
+    q = mygprec_special(p,bitprec2); lc = leading_term(q);
+    if (h > 1) m = gmul(m,lc);
 
-    e = gexpo(gsub(mygprec_special(p,bitprec2), m))
-      - gexpo(leading_term(q)) + (long)log2((double)n) + 1;
+    e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
     if (e<-2*bitprec2) e=-2*bitprec2; /* to avoid e=-pariINFINITY */
     if (e < 0)
     {
-      e = a_posteriori_errors(q,roots_pol,e);
-      if (e < -bitprec) return roots_pol;
+      e = bitprec + a_posteriori_errors(p,roots_pol,e);
+      if (e < 0) return roots_pol;
     }
     if (DEBUGLEVEL > 7)
       fprintferr("all_roots: restarting, i = %ld, e = %ld\n", i,e);
@@ -2079,7 +2051,7 @@ isrealappr(GEN x, long e)
 static int
 isconj(GEN x, GEN y, long e)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long i= (gexpo( gsub((GEN)x[1],(GEN)y[1]) ) < e
         && gexpo( gadd((GEN)x[2],(GEN)y[2]) ) < e);
   avma = av; return i;
@@ -2091,7 +2063,7 @@ isconj(GEN x, GEN y, long e)
 GEN
 roots(GEN p, long l)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long n,i,k,s,t,e;
   GEN c,L,p1,res,rea,com;