/*******************************************************************/
/*                                                                 */
/*                ROOTS OF COMPLEX POLYNOMIALS                     */
/*                (written by Xavier Gourdon)                      */
/*                                                                 */
/*******************************************************************/
/* $Id: rootpol.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
#include "pari.h"

GEN polrecip_i(GEN x);
#define pariINFINITY 100000
#define NEWTON_MAX 10

static GEN gunr, globalu;
static long KARASQUARE_LIMIT, COOK_SQUARE_LIMIT, Lmax, step4;
static double *radius;

/********************************************************************/
/**                                                                **/
/**                      ARITHMETIQUE RAPIDE                       **/
/**                                                                **/
/********************************************************************/

/* fast product of x,y which must be integer or complex of integer */
static GEN
quickmulcc(GEN x, GEN y)
{
  long tx=typ(x),ty=typ(y);
  GEN z;

  if (tx==t_INT)
  {
    if (ty==t_INT) return mulii(x,y);
    if (ty==t_COMPLEX)
    {
      z=cgetg(3,t_COMPLEX);
      z[1]=(long) mulii(x,(GEN) y[1]);
      z[2]=(long) mulii(x,(GEN) y[2]);
      return z;
    }
  }

  if (tx==t_COMPLEX)
  {
    if (ty==t_INT)
    {
      z=cgetg(3,t_COMPLEX);
      z[1]=(long) mulii((GEN)x[1],y);
      z[2]=(long) mulii((GEN)x[2],y);
      return z;
    }
    if (ty==t_COMPLEX)
    {
      long av,tetpil;
      GEN p1,p2;

      z=cgetg(3,t_COMPLEX); av=avma;
      p1=mulii((GEN)x[1],(GEN)y[1]); p2=mulii((GEN)x[2],(GEN)y[2]);
      x=addii((GEN)x[1],(GEN)x[2]); y=addii((GEN)y[1],(GEN)y[2]);
      y=mulii(x,y); x=addii(p1,p2);
      tetpil=avma; z[1]=lsubii(p1,p2); z[2]=lsubii(y,x);
      gerepilemanyvec(av,tetpil,z+1,2);
      return z;
    }
  }
  err(talker,"bug in quickmulcc");
  return NULL; /* not reached */
}

static void
set_karasquare_limit(long bitprec)
{
  if (bitprec<600) { KARASQUARE_LIMIT=8; COOK_SQUARE_LIMIT=400; return; }
  if (bitprec<2000) { KARASQUARE_LIMIT=4; COOK_SQUARE_LIMIT=200; return; }
  if (bitprec<3000) { KARASQUARE_LIMIT=4; COOK_SQUARE_LIMIT=125; return; }
  if (bitprec<5000) { KARASQUARE_LIMIT=2; COOK_SQUARE_LIMIT=75; return; }
  KARASQUARE_LIMIT=1; COOK_SQUARE_LIMIT=50;
}

/* the pari library does not have specific procedure for the square of
polynomials. This one is twice faster than gmul */
static GEN
mysquare(GEN p)
{
  GEN s,aux1,aux2;
  long i,j,n=lgef(p)-3,nn,ltop,lbot;

  if (n==-1) return gcopy(p);
  nn=n<<1; s=cgetg(nn+3,t_POL);
  s[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(nn+3);
  for (i=0; i<=n; i++)
  {
    aux1=gzero; ltop=avma;
    for (j=0; j<(i+1)>>1; j++)
    {
      aux2=quickmulcc((GEN) p[j+2],(GEN)p[i-j+2]);
      aux1=gadd(aux1,aux2);
    }
    if (i%2==1) { lbot=avma; s[i+2]=lpile(ltop,lbot,gshift(aux1,1)); }
    else
    {
      aux1=gshift(aux1,1);
      aux2=quickmulcc((GEN)p[2+(i>>1)],(GEN)p[2+(i>>1)]);
      lbot=avma; s[i+2]=lpile(ltop,lbot,gadd(aux1,aux2));
    }
  }
  for (i=n+1; i<=nn; i++)
  {
    aux1=gzero; ltop=avma;
    for (j=i-n; j<(i+1)>>1; j++)
    {
      aux2=quickmulcc((GEN)p[j+2],(GEN)p[i-j+2]);
      aux1=gadd(aux1,aux2);
    }
    if (i%2==1) { lbot=avma; s[i+2]=lpile(ltop,lbot,gshift(aux1,1)); }
    else
    {
      aux1=gshift(aux1,1);
      aux2=quickmulcc((GEN)p[2+(i>>1)],(GEN)p[2+(i>>1)]);
      lbot=avma; s[i+2]=lpile(ltop,lbot,gadd(aux1,aux2));
    }
  }
  return s;
}

static GEN
karasquare(GEN p)
{
  GEN p1,s0,s1,s2,aux;
  long n=lgef(p)-3,n0,n1,i,var,nn0,ltop,lbot;

  if (n<=KARASQUARE_LIMIT) return mysquare(p);
  ltop=avma;
  var=evalsigne(1)+evalvarn(varn(p)); n0=n>>1; n1=n-n0-1;
  setlgef(p,n0+3); /* hack to have the first half of p */
  s0=karasquare(p);
  p1=cgetg(n1+3,t_POL); p1[1]=var+evallgef(n1+3);
  for (i=2; i<=n1+2; i++) p1[i]=p[1+i+n0];
  s2=karasquare(p1);
  s1=karasquare(gadd(p,p1));
  s1=gsub(s1,gadd(s0,s2));
  nn0=n0<<1;
  aux=cgetg((n<<1)+3,t_POL); aux[1]=var+evallgef(2*n+3);
  var=lgef(s0);
  for (i=2; i<var; i++) aux[i]=s0[i];
  for (   ; i<=nn0+2; i++) aux[i]=zero;
  var=lgef(s2);
  for (i=2; i<var; i++) aux[2+i+nn0]=s2[i];
  for (i=var-2; i<=(n1<<1); i++) aux[4+i+nn0]=zero;
  aux[3+nn0]=zero;
  for (i=3; i<=lgef(s1); i++)
    aux[i+n0]=ladd((GEN) aux[i+n0],(GEN) s1[i-1]);
  setlgef(p,n+3); /* recover all the polynomial p */
  lbot=avma; return gerepile(ltop,lbot,gcopy(aux));
}

static GEN
cook_square(GEN p)
{
  GEN p0,p1,p2,p3,q,aux0,aux1,r,aux,plus,moins;
  long n=lgef(p)-3,n0,n3,i,j,ltop=avma,lbot,var;

  if (n<=COOK_SQUARE_LIMIT) return karasquare(p);

  n0=(n+1)/4; n3=n+1-3*n0;
  p0=cgetg(n0+2,t_POL); p1=cgetg(n0+2,t_POL); p2=cgetg(n0+2,t_POL);
  p3=cgetg(n3+2,t_POL);
  var=evalsigne(1)|evalvarn(varn(p));
  p0[1]=p1[1]=p2[1]=var|evallgef(n0+2); p3[1]=var|evallgef(n3+2);

  for (i=0; i<n0; i++)
  {
    p0[i+2]=p[i+2]; p1[i+2]=p[i+n0+2]; p2[i+2]=p[i+2*n0+2];
  }
  for (i=0; i<n3; i++) p3[i+2]=p[i+3*n0+2];

  q=cgetg(8,t_VEC); q=q+4;

  q[0]=(long) p0;
  aux0=gadd(p0,p2); aux1=gadd(p1,p3);
  q[-1]=lsub(aux0,aux1); q[1]=ladd(aux0,aux1);
  aux0=gadd(p0,gmulgs(p2,4)); aux1=gmulgs(gadd(p1,gmulgs(p3,4)),2);
  q[-2]=lsub(aux0,aux1); q[2]=ladd(aux0,aux1);
  aux0=gadd(p0,gmulgs(p2,9)); aux1=gmulgs(gadd(p1,gmulgs(p3,9)),3);
  q[-3]=lsub(aux0,aux1); q[3]=ladd(aux0,aux1);
  for (i=-3; i<=3; i++) q[i]=(long) cook_square((GEN)q[i]);
  r=new_chunk(7);
  plus=cgetg(4,t_VEC); moins=cgetg(4,t_VEC);
  for (i=1; i<=3; i++)
  {
    plus[i]=ladd((GEN)q[-i],(GEN)q[i]);
    moins[i]=lsub((GEN)q[-i],(GEN)q[i]);
  }
  r[0]=q[0];
  r[1]=ldivgs(
	      gsub(
		   gsub(gmulgs((GEN)moins[2],9),(GEN)moins[3]),
		   gmulgs((GEN)moins[1],45)),
	      60);
  r[2]=ldivgs(
	      gadd(
		   gadd(gmulgs((GEN)plus[1],270),gmulgs((GEN)q[0],-490)),
		   gadd(gmulgs((GEN)plus[2],-27),gmulgs((GEN)plus[3],2))),
	      360);
  r[3]=ldivgs(
	      gadd(
		   gadd(gmulgs((GEN)moins[1],13),gmulgs((GEN)moins[2],-8)),
		   (GEN)moins[3]),
	      48);
  r[4]=ldivgs(
	      gadd(
		   gadd(gmulgs((GEN)q[0],56),gmulgs((GEN)plus[1],-39)),
		   gsub(gmulgs((GEN)plus[2],12),(GEN)plus[3])),
	      144);
  r[5]=ldivgs(
	      gsub(
		   gadd(gmulgs((GEN)moins[1],-5),gmulgs((GEN)moins[2],4)),
		   (GEN)moins[3]),
	      240);
  r[6]=ldivgs(
	      gadd(
		   gadd(gmulgs((GEN)q[0],-20),gmulgs((GEN)plus[1],15)),
		   gadd(gmulgs((GEN)plus[2],-6),(GEN)plus[3])),
	      720);
  q=cgetg(2*n+3,t_POL); q[1]=var|evallgef(2*n+3);
  for (i=0; i<=2*n; i++) q[i+2]=zero;
  for (i=0; i<=6; i++)
  {
    aux=(GEN) r[i];
    for (j=0; j<=lgef(aux)-3; j++)
      q[n0*i+j+2]=ladd((GEN)q[n0*i+j+2],(GEN)aux[j+2]);
  }
  lbot=avma; return gerepile(ltop,lbot,gcopy(q));
}

static GEN
graeffe(GEN p)
{
  GEN p0,p1,s0,s1,ss1;
  long n=lgef(p)-3,n0,n1,i,auxi,ns1;

  if (n==0) return gcopy(p);
  n0=n>>1; n1=(n-1)>>1;
  auxi=evalsigne(1)|evalvarn(varn(p));
  p0=cgetg(n0+3,t_POL); p0[1]=auxi|evallgef(n0+3);
  p1=cgetg(n1+3,t_POL); p1[1]=auxi|evallgef(n1+3);
  for (i=0; i<=n0; i++) p0[i+2]=p[2+(i<<1)];
  for (i=0; i<=n1; i++) p1[i+2]=p[3+(i<<1)];

  s0=cook_square(p0);
  s1=cook_square(p1); ns1 = lgef(s1)-3;
  ss1 = cgetg(ns1+4, t_POL);
  ss1[1] = auxi | evallgef(ns1+4);
  ss1[2]=zero;
  for (i=0; i<=ns1; i++) ss1[3+i]=lneg((GEN)s1[2+i]);
  /* now ss1 contains -x * s1 */
  return gadd(s0,ss1);
}

/********************************************************************/
/**                                                                **/
/**        FACTORISATION SQUAREFREE AVEC LE GCD MODULAIRE          **/
/**                                                                **/
/********************************************************************/

/* return a n x 2 matrix:
 *   col 1 contains the i's such that A_i non constant
 *   col 2 the A_i's, s.t. pol = A_i1^i1.A_i2^i2...A_in^in.
 * if pol is constant return [;]
 */
GEN
square_free_factorization(GEN pol)
{
  long deg,i,j,m;
  GEN p1,x,t1,v1,t,v,A;

  if (typ(pol)!=t_POL) err(typeer,"square_free_factorization");
  deg=lgef(pol)-3; if (deg<1) return cgetg(1,t_MAT);
  p1 = content(pol); if (!gcmp1(p1)) pol = gdiv(pol,p1);

  x=cgetg(3,t_MAT);
  if (deg > 1)
  {
    t1 = modulargcd(pol,derivpol(pol));
    if (isscalar(t1)) deg = 1;
  }
  if (deg==1)
  {
    x[1]=lgetg(2,t_COL); p1=(GEN)x[1]; p1[1]=un;
    x[2]=lgetg(2,t_COL); p1=(GEN)x[2]; p1[1]=(long)pol; return x;
  }
  A=new_chunk(deg+1); v1=gdivexact(pol,t1); v=v1; i=0;
  while (lgef(v)>3)
  {
    v=modulargcd(t1,v1); i++;
    A[i]=(long)gdivexact(v1,v);
    t=gdivexact(t1,v); v1=v; t1=t;
  }
  m=1; x[1]=lgetg(deg+1,t_COL); x[2]=lgetg(deg+1,t_COL);
  for (j=1; j<=i; j++)
    if (isnonscalar(A[j]))
    {
      p1=(GEN)x[1]; p1[m] = lstoi(j);
      p1=(GEN)x[2]; p1[m] = A[j];
      m++;
    }
  setlg(x[1],m); setlg(x[2],m); return x;
}

/********************************************************************/
/**                                                                **/
/**                 CALCUL DU MODULE DES RACINES                   **/
/**                                                                **/
/********************************************************************/

static double
log2ir(GEN x)
{
  double l;

  if (signe(x)==0) return (double) -pariINFINITY;
  if (typ(x)==t_INT)
  {
    if (lgef(x)==3) return (double) log2( (double)(ulong) x[2]);
    l=(double)(ulong) x[2]+
	(double)(ulong) x[3] / exp2((double) BITS_IN_LONG);
    return log2(l)+ (double) BITS_IN_LONG * (lgef(x)-3.);
  }
  /* else x is real */
  return 1.+ (double) expo(x)+log2( (double)(ulong) x[2]) - (double) BITS_IN_LONG;
}

static double
mylog2(GEN z)
{
  double x,y;

  if (typ(z)!=t_COMPLEX) return log2ir(z);

  x=log2ir((GEN) z[1]); y=log2ir((GEN) z[2]);
  if (fabs(x-y)>10) return (x>y)? x: y;
  return x+0.5*log2( 1 + exp2(2*(y-x)));
}

static long
findpower(GEN p)
{
  double x, logbinomial,pente,pentemax=-pariINFINITY;
  long n=lgef(p)-3,i;

  logbinomial=mylog2((GEN) p[n+2]);
  for (i=n-1; i>=0; i--)
  {
    logbinomial=logbinomial+log2((double) (i+1) / (double) (n-i));
    x=mylog2((GEN) p[2+i])-logbinomial;
    if (x>-pariINFINITY)
    {
      pente=x/ (double) (n-i);
      if (pente>pentemax) pentemax=pente;
    }
  }
  return (long) -floor(pentemax);
}

/* returns the exponent for the procedure modulus, from the newton diagram */
static long
polygone_newton(GEN p, long k)
{
  double *logcoef,pente;
  long n=lgef(p)-3,i,j,h,l,*sommet,pentelong;

  logcoef=(double*) gpmalloc((n+1)*sizeof(double));
  sommet=(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)
  {
    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))
      {
	h=j;
	pente=(logcoef[j]-logcoef[i])/(double)(j-i);
      }
    }
    i=h;
    sommet[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;
}

/* change z into z*2^e, where z is real or complex of real */
static void
myshiftrc(GEN z, long e)
{
  if (typ(z)==t_COMPLEX)
  {
    if (signe(z[1])!=0) setexpo(z[1], expo(z[1])+e);
    if (signe(z[2])!=0) setexpo(z[2], expo(z[2])+e);
  }
  else
    if (signe(z)!=0) setexpo(z,expo(z)+e);
}

/* return z*2^e, where z is integer or complex of integer (destroy z) */
static GEN
myshiftic(GEN z, long e)
{
  if (typ(z)==t_COMPLEX)
  {
    z[1]=signe(z[1])? lmpshift((GEN) z[1],e): zero;
    z[2]=lmpshift((GEN) z[2],e);
    return z;
  }
  return signe(z)? mpshift(z,e): gzero;
}

static GEN
mygprecrc(GEN x, long bitprec, long e)
{
  long tx=typ(x);
  GEN y;

  if (bitprec<0) bitprec=0; /* should rarely happen */
  switch(tx)
  {
    case t_REAL:
      y=cgetr(bitprec/BITS_IN_LONG+3); affrr(x,y);
      if (!signe(x)) setexpo(y,-bitprec+e);
      break;
    case t_COMPLEX:
      y=cgetg(3,t_COMPLEX);
      y[1]=(long) mygprecrc((GEN)x[1],bitprec,e);
      y[2]=(long) mygprecrc((GEN)x[2],bitprec,e);
      break;
    default: y=gcopy(x);
  }
  return y;
}

/* gprec behaves badly with the zero for polynomials.
The second parameter in mygprec is the precision in base 2 */
static GEN
mygprec(GEN x, long bitprec)
{
  long tx=typ(x),lx,i,e;
  GEN y;

  switch(tx)
  {
    case t_POL:
      lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1]; e=gexpo(x);
      for (i=2; i<lx; i++) y[i]=(long) mygprecrc((GEN)x[i],bitprec,e);
      break;

    default: y=mygprecrc(x,bitprec,0);
  }
  return y;
}

/* the round fonction has a bug in pari. Thus I create mygfloor, using gfloor
which has no bug (destroy z)*/
static GEN
mygfloor(GEN z)
{
  if (typ(z)!=t_COMPLEX) return gfloor(z);
  z[1]=lfloor((GEN)z[1]); z[2]=lfloor((GEN)z[2]); return z;
}

/* returns a polynomial q in (Z[i])[x] keeping bitprec bits of p */
static GEN
eval_rel_pol(GEN p,long bitprec)
{
  long e=gexpo(p),n=lgef(p),i,shift;
  GEN q = gprec(p,(long) ((double) bitprec * L2SL10)+2);

  shift=bitprec-e+1;
  for (i=2; i<n; i++)
    q[i]=(long) mygfloor(myshiftic((GEN)q[i],shift));
  return q;
}

/* normalize a polynomial p, that is change it with coefficients in Z[i],
after making product by 2^bitprec */
static void
pol_to_gaussint(GEN p, long bitprec)
{
  long i,n=lgef(p);
  for (i=2; i<n; i++)
  {
    myshiftrc((GEN) p[i],bitprec);
    p[i]=(long) mygfloor((GEN) p[i]);
  }
}

/* returns p(R*x)/R^n (in R or R[i]), n=deg(p), bitprec bits of precision */
static GEN
homothetie(GEN p, double R, long bitprec)
{
  GEN q,r,gR,aux;
  long n=lgef(p)-3,i;

  gR=mygprec(dbltor(1/R),bitprec);
  q=mygprec(p,bitprec);
  r=cgetg(n+3,t_POL); r[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(n+3);
  aux=gR; r[n+2]=q[n+2];
  for (i=n-1; i>=0; i--)
  {
    r[i+2]=lmul(aux,(GEN)q[i+2]);
    aux=gmul(aux,gR);
  }
  return r;
}

/* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) */
static void
homothetie2n(GEN p, long e)
{
  if (e)
  {
    long i,n=lgef(p)-1;
    for (i=2; i<=n; i++) myshiftrc((GEN) p[i],(n-i)*e);
  }
}

/* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
static void
homothetie_gauss(GEN p, long e,long f)
{
  if (e||f)
  {
    long i, n=lgef(p)-1;
    for (i=2; i<=n; i++) p[i]=(long) myshiftic((GEN) p[i],f+(n-i)*e);
  }
}

static long
valuation(GEN p)
{
  long j=0,n=lgef(p)-3;

  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=lgef(p)-3,i,j,ltop=avma;
  GEN a,s,icd;
  double r,*rho;

  if (n<4) { *k=n; return 0.; }
  a=cgetg(6,t_POL); s=cgetg(6,t_POL);
  rho=(double *) gpmalloc(4*sizeof(double));
  icd=gdiv(gunr,(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++)
  {
    s[i+1]=lmulsg(i,(GEN)a[i+1]);
    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.;
    else
      rho[i-1]=exp(log(r/(double)n)/(double) i);
  }
  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)));
  else
    *k=n;
  free(rho); avma=ltop; return r;
}

/* returns the maximum of the modulus of p with a relative error tau */
static GEN
gmax_modulus(GEN p, double tau)
{
  GEN q,aux,new_gunr;
  long i,j,k,valuat,n=lgef(p)-3,nn,ltop=avma,bitprec,imax,e;
  double r=0.,rho,tau2,eps;

  if (tau>3.0) tau=3.0; /* fix PZ 7.2.98: ensures eps is positive */
  eps=1/log(4./tau); tau2=tau/6.;
  bitprec=(long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
  new_gunr=mygprec(gunr,bitprec+2*n);
  aux=gdiv(new_gunr,(GEN) p[2+n]);
  q=gmul(aux,p); q[2+n]=lcopy(new_gunr);
  k=nn=n;
  e=findpower(q); homothetie2n(q,e); r=-(double) e;
  q=mygprec(q,bitprec+(n<<1));
  pol_to_gaussint(q,bitprec);
  imax=(long) ((log(log(4.*n))+log(3./tau))/log(2.))+2;
  for (i=0,e=0;;)
  {
    rho=lower_bound(q,&k,eps);
    if (rho>exp2(-(double) e)) e = (long) -floor(log2(rho));
    r = r-(double) e/ exp2( (double) i);
    if (++i == imax)
    {
      avma=ltop;
      return gpui(dbltor(2.),dbltor(r),DEFAULTPREC);
    }

    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;
    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];
    }
    set_karasquare_limit(gexpo(q));
    q = gerepileupto(ltop, graeffe(q));
    tau2=1.5*tau2; eps=1/log(1./tau2);
    e = findpower(q);
  }
}

static double
max_modulus(GEN p, double tau)
{
  return gtodouble(gmax_modulus(p,tau));
}

/* return the k-th modulus (in ascending order) of p, rel. error tau*/
static double
modulus(GEN p, long k, double tau)
{
  GEN q,new_gunr;
  long i,j,kk=k,imax,n=lgef(p)-3,nn,nnn,valuat,av,ltop=avma,bitprec,decprec,e;
  double tau2,r;

  tau2=tau/6; nn=n;
  bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2)));
  decprec=(long) ((double) bitprec * L2SL10)+1;
  new_gunr=gprec(gunr,decprec);
  av = avma;
  q=gprec(p,decprec);
  q=gmul(new_gunr,q);
  e=polygone_newton(q,k);
  homothetie2n(q,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);

    nnn=lgef(q)-3; 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);
    r=r+e/exp2((double)i);
    q=gmul(new_gunr,q);
    homothetie2n(q,e);

    tau2=1.5*tau2; if (tau2>1.) tau2=1.;
    bitprec= 1+(long) ((double) nn*(2.+log2(3.*(double) nn)+log2(1./tau2)));
  }
  avma=ltop; return exp2(-r);
}

/* return the k-th modulus r_k of p, rel. error tau, knowing that
rmin < r_k < rmax. This helps because the information enable us to use
less precision... quicker ! */
static double
pre_modulus(GEN p, long k, double tau, double rmin, double rmax)
{
  GEN q;
  long n=lgef(p)-3,i,imax,imax2,bitprec,ltop=avma;
  double aux,tau2,R;

  tau2=tau/6.; aux=sqrt(rmax/rmin)*exp(4*tau2);
  imax=(long) log2(log((double)n)/log(aux));
  if (imax<=0) return modulus(p,k,tau);

  R=sqrt(rmin*rmax);
  bitprec=(long) ((double) n*(2.+log2(aux)+log2(1./tau2)));
  q=homothetie(p,R,bitprec);
  imax2=(long) ((log2(3./tau)+log2(log(4.*(double) n)) ))+1;
  if (imax>imax2) imax=imax2;

  for (i=0; i<imax; i++)
  {
    q=eval_rel_pol(q,bitprec);
    set_karasquare_limit(bitprec);
    q = gerepileupto(ltop, graeffe(q));
    aux=aux*aux*exp(2*tau2);
    tau2=1.5*tau2;
    bitprec= (long) ((double) n*(2.+log2(aux)+log2(1./(1-exp(-tau2)))));
    q=gmul(mygprec(gunr,bitprec),q);
  }

  aux=modulus(q,k,exp2((double)imax)*tau/3.);
  R=R*exp(log(aux)*exp2(-(double)imax));
  avma=ltop; return R;
}

/* returns the minimum of the modulus of p with a relative error tau */
static double
min_modulus(GEN p, double tau)
{
  long ltop=avma;
  double r;

  if (isexactzero((GEN)p[2])) return 0.;
  r=1./max_modulus(polrecip_i(p),tau); avma=ltop; return r;
}

/* returns k such that r_k e^(-tau) < R < r_{ k+1 } e^tau.
l is such that you know in advance that l<= k <= n-l */
static long
dual_modulus(GEN p, double R, double tau, long l)
{
  GEN q;
  long i,j,imax,k,delta_k=0,n=lgef(p)-3,nn,nnn,valuat,ltop=avma,bitprec,ll=l;
  double logmax,aux,tau2;

  tau2=7.*tau/8.;
  bitprec=6*n-5*l+(long) ((double) n*(log2(1/tau2)+8.*tau2/7.));
  q=homothetie(p,R,bitprec);
  nn=n;
  imax=(long)(log(log(2.*(double)n)/tau2)/log(7./4.)+1);

  for (i=0; i<imax; i++)
  {
    bitprec=6*nn-5*ll+(long) ((double) nn*(log2(1/tau2)+8.*tau2/7.));

    q=eval_rel_pol(q,bitprec);
    nnn=lgef(q)-3; 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);
    }
    ll= (-valuat<nnn-n)? ll-valuat: ll+nnn-n; /* min(ll-valuat,ll+nnn-n) */
    if (ll<0) ll=0;

    nn=nnn-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;
  for (i=0; i<=lgef(q)-3; i++)
  {
    aux=mylog2((GEN)q[2+i]);
    if (aux>logmax) { logmax=aux; k=i; }
  }
  avma=ltop; return delta_k+k;
}

/********************************************************************/
/**                                                                **/
/**       CALCUL D'UN FACTEUR PAR INTEGRATION SUR LE CERCLE        **/
/**                                                                **/
/********************************************************************/

static GEN
gmulbyi(GEN z)
{
  GEN aux = cgetg(3,t_COMPLEX);

  if (typ(z)!=t_COMPLEX)
  {
    aux[1]=zero;
    aux[2]=(long) z;
  }
  else
  {
    aux[1]=lneg((GEN)z[2]);
    aux[2]=z[1];
  }
  return aux;
}

static void
fft(GEN Omega, GEN p, GEN f, long step, long l)
{
  long i,l1,l2,l3,rap=Lmax/l,rapi,step4,ltop,lbot;
  GEN f1,f2,f3,f02,f13,g02,g13,ff;

  if (l==2)
  {
    f[0]=ladd((GEN)p[0],(GEN)p[step]);
    f[1]=lsub((GEN)p[0],(GEN)p[step]); return;
  }
  if (l==4)
  {
    f1=gadd((GEN)p[0],(GEN)p[(step<<1)]);
    f3=gadd((GEN)p[step],(GEN)p[3*step]);
    f[0]=ladd(f1,f3);
    f[2]=lsub(f1,f3);

    f2=gsub((GEN)p[0],(GEN)p[(step<<1)]);
    f02=gsub((GEN)p[step],(GEN)p[3*step]);
    f02=gmulbyi(f02);
    f[1]=ladd(f2,f02);
    f[3]=lsub(f2,f02);
    return;
  }

  l1=(l>>2); l2=(l1<<1); l3=l1+l2; step4=(step<<2);

  ltop=avma;
  fft(Omega,p,f,step4,l1);
  fft(Omega,p+step,f+l1,step4,l1);
  fft(Omega,p+(step<<1),f+l2,step4,l1);
  fft(Omega,p+3*step,f+l3,step4,l1);

  ff=cgetg(l+1,t_VEC);
  for (i=0; i<l1; i++)
  {
    rapi=rap*i;
    f1=gmul((GEN)Omega[rapi],(GEN) f[i+l1]);
    f2=gmul((GEN)Omega[(rapi<<1)],(GEN) f[i+l2]);
    f3=gmul((GEN)Omega[3*rapi],(GEN) f[i+l3]);

    f02=gadd((GEN)f[i],f2);
    g02=gsub((GEN)f[i],f2);
    f13=gadd(f1,f3);
    g13=gmulbyi(gsub(f1,f3));

    ff[i+1]=ladd(f02,f13);
    ff[i+l1+1]=ladd(g02,g13);
    ff[i+l2+1]=lsub(f02,f13);
    ff[i+l3+1]=lsub(g02,g13);
  }
  lbot=avma; ff=gerepile(ltop,lbot,gcopy(ff));
  for (i=0; i<l; i++) f[i]=ff[i+1];
}

/* returns a vector RU which contains exp(2*i*k*Pi/NN), k=0..NN-1 */
static GEN
initRU(long NN, long bitprec)
{
  GEN prim,aux,RU,mygpi;
  long i,N2=(NN>>1),N4=(NN>>2),N8=(NN>>3),decprec;

  RU=cgetg(NN+1,t_VEC); RU++;
  mygpi=mppi(bitprec/BITS_IN_LONG+3);
  aux=gmul(gi,gdivgs(mygpi,NN/2)); /* 2i Pi/NN */
  decprec=(long) ((double) bitprec * L2SL10)+1;
  prim=gexp(aux,decprec);
  RU[0]=lprec(gunr,decprec);

  for (i=1; i<=N8; i++) RU[i]=lmul(prim,(GEN)RU[i-1]);
  for (i=1; i<N8; i++)
  {
    aux=cgetg(3,t_COMPLEX);
    aux[1]=((GEN)RU[i])[2]; aux[2]=((GEN)RU[i])[1];
    RU[N4-i]=(long) aux;
  }
  for (i=0; i<N4; i++) RU[i+N4]=(long)gmulbyi((GEN)RU[i]);
  for (i=0; i<N2; i++) RU[i+N2]=lneg((GEN)RU[i]);
  return RU;
}

/* returns 1 if p has only real coefficients, 0 else */
static long
isreal(GEN p)
{
  long n=lgef(p)-3,i=0;

  while (i<=n && typ(p[i+2])!=t_COMPLEX) i++;
  return (i>n);
}

static void
parameters(GEN p, double *mu, double *gamma,
           long polreal, double param, double param2)
{
  GEN q,pc,Omega,coef,RU,prim,aux,ggamma,gx,mygpi;
  long ltop=avma,limite=stack_lim(ltop,1),n=lgef(p)-3,bitprec,NN,K,i,j,ltop2;
  double lx;

  bitprec=gexpo(p)+(long)param2+8;
  NN=(long) (param*3.14)+1; if (NN<Lmax) NN=Lmax;
  K=NN/Lmax; if (K%2==1) K++; NN=Lmax*K;
  mygpi=mppi(bitprec/BITS_IN_LONG+3);
  aux=gmul(gi,gdivgs(mygpi,NN/2)); /* 2i Pi/NN */
  prim=gexp(aux,(long) ((double) bitprec * L2SL10)+1);
  RU=mygprec(gunr,bitprec);

  Omega=initRU(Lmax,bitprec);

  q=mygprec(p,bitprec);
  pc=cgetg(Lmax+1,t_VEC); pc++;
  for (i=n+1; i<Lmax; i++) pc[i]=zero;

  coef=cgetg(Lmax+1,t_VEC); coef++;
  *mu=(double)pariINFINITY; *gamma=0.;
  ggamma=gmul(gzero,gunr);
  if (polreal) K=K/2+1;
  ltop2=avma;
  for (i=0; i<K; i++)
  {
    aux=mygprec(gunr,bitprec);
    for (j=0; j<=n; j++)
    {
      pc[j]=lmul((GEN)q[j+2],aux);
      aux=gmul(aux,RU); /* RU = prim^i, aux=prim^(ij) */
    }

    fft(Omega,pc,coef,1,Lmax);
    for (j=0; j<Lmax; j++)
    {
      aux=gprec((GEN)coef[j],DEFAULTPREC);
      gx=gabs(aux,DEFAULTPREC);
      lx=gtodouble(mplog(gx));
      if (lx<*mu) *mu=lx;
      if (polreal && (i>0 && i<K-1))
      {
	gx=gdiv(gdeux,gx);
	ggamma=gadd(ggamma,gx);
      }
      else
	ggamma=gadd(ggamma,ginv(gx));
    }
    RU=gmul(RU,prim);
    if (low_stack(limite, stack_lim(ltop,1)))
    {
      GEN *gptr[2];
      if(DEBUGMEM>1) err(warnmem,"parameters");
      gptr[0]=&ggamma; gptr[1]=&RU; gerepilemany(ltop2,gptr,2);
    }
  }
  ggamma=gdivgs(ggamma,NN);
  *gamma=gtodouble(glog(ggamma,DEFAULTPREC))/log(2.);
  avma=ltop;
}

/* NN is a multiple of Lmax */
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=lgef(p)-3,i,j,K,ltop;

  mygpi=mppi(bitprec/BITS_IN_LONG+3);
  aux=gmul(gi,gdivgs(mygpi,NN/2)); /* 2i Pi/NN */
  prim=gexp(aux,(long) ((double) bitprec * L2SL10)+1);
  prim2=mygprec(gunr,bitprec);
  RU=cgetg(n+2,t_VEC); RU++;

  Omega=initRU(Lmax,bitprec);
  K=NN/Lmax; q=mygprec(p,bitprec);
  qd=derivpol(q);

  pc=cgetg(Lmax+1,t_VEC); pc++;
  for (i=n+1; i<Lmax; i++) pc[i]=zero;
  pdc=cgetg(Lmax+1,t_VEC); pdc++;
  for (i=n; i<Lmax; i++) pdc[i]=zero;

  alpha=cgetg(Lmax+1,t_VEC); alpha++;
  beta=cgetg(Lmax+1,t_VEC); beta++;
  gamma=cgetg(Lmax+1,t_VEC); gamma++;

  if (polreal) K=K/2+1;

  ltop=avma; limite = stack_lim(ltop,1);
  W=cgetg(k+1,t_VEC); U=cgetg(k+1,t_VEC);
  for (i=1; i<=k; i++) W[i]=U[i]=zero;

  for (i=0; i<K; i++)
  {
    RU[0]=(long) gun;
    for (j=1; j<=n; j++) RU[j]=lmul((GEN)RU[j-1],prim2);
    /* RU[j]=prim^{ ij }=prim2^j */

    for (j=0; j<n; j++) pdc[j]=lmul((GEN)qd[j+2],(GEN)RU[j]);
    fft(Omega,pdc,alpha,1,Lmax);
    for (j=0; j<=n; j++) pc[j]=lmul((GEN)q[j+2],(GEN)RU[j]);
    fft(Omega,pc,beta,1,Lmax);
    for (j=0; j<Lmax; j++) gamma[j]=linv((GEN)beta[j]);
    for (j=0; j<Lmax; j++) beta[j]=lmul((GEN)alpha[j],(GEN)gamma[j]);
    fft(Omega,beta,alpha,1,Lmax);
    fft(Omega,gamma,beta,1,Lmax);

    if (polreal) /* p has real coefficients */
    {
      if (i>0 && i<K-1)
      {
	for (j=1; j<=k; j++)
	{
	  aux=gmul((GEN)alpha[j+1],(GEN)RU[j+1]);
	  W[j]=ladd((GEN)W[j],gshift(greal(aux),1));
	  aux=gmul((GEN)beta[j],(GEN)RU[j]);
	  U[j]=ladd((GEN)U[j],gshift(greal(aux),1));
	}
      }
      else
      {
	for (j=1; j<=k; j++)
	{
	  aux=gmul((GEN)alpha[j+1],(GEN)RU[j+1]);
	  W[j]=ladd((GEN)W[j],greal(aux));
	  aux=gmul((GEN)beta[j],(GEN)RU[j]);
	  U[j]=ladd((GEN)U[j],greal(aux));
	}
      }
    }
    else
    {
      for (j=1; j<=k; j++)
      {
	W[j]=ladd((GEN)W[j],gmul((GEN)alpha[j+1],(GEN)RU[j+1]));
	U[j]=ladd((GEN)U[j],gmul((GEN)beta[j],(GEN)RU[j]));
      }
    }
    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);
    }
  }

  for (i=1; i<=k; i++)
  {
    aux=(GEN)W[i];
    for (j=1; j<i; j++) aux=gadd(aux,gmul((GEN)W[i-j],(GEN)F[k+2-j]));
    F[k+2-i] = ldivgs(aux,-i*NN);
  }
  for (i=0; i<k; i++)
  {
    aux=(GEN)U[k-i];
    for (j=1+i; j<k; j++) aux=gadd(aux,gmul((GEN)F[2+j],(GEN)U[j-i]));
    H[i+2] = ldivgs(aux,NN);
  }
}

static GEN
refine_H(GEN F, GEN G, GEN HH, long bitprec, long shiftbitprec)
{
  GEN H=HH,D,aux;
  long ltop=avma, limite=stack_lim(ltop,1),error=0,i,bitprec1,bitprec2,lbot;

  D=gsub(gun,gres(gmul(HH,G),F)); error=gexpo(D);
  bitprec2=bitprec+shiftbitprec;

  for (i=0; (error>-bitprec && i<NEWTON_MAX) && error<=0; i++)
  {
    if (low_stack(limite, stack_lim(ltop,1)))
    {
      GEN *gptr[2];
      if(DEBUGMEM>1) err(warnmem,"refine_H");
      gptr[0]=&D; gptr[1]=&H; gerepilemany(ltop,gptr,2);
    }
    bitprec1=-error+shiftbitprec;
    aux=gmul(mygprec(H,bitprec1),mygprec(D,bitprec1));
    aux=mygprec(aux,bitprec1);
    aux=gres(aux,mygprec(F,bitprec1));

    bitprec1=-error*2+shiftbitprec;
    if (bitprec1>bitprec2) bitprec1=bitprec2;
    H=gadd(mygprec(H,bitprec1),aux);
    D=gsub(gun,gres(gmul(H,G),F));
    error=gexpo(D); if (error<-bitprec1) error=-bitprec1;
  }
  if (error<=-bitprec/2) { lbot=avma; return gerepile(ltop,lbot,gcopy(H)); }
  avma=ltop; return gzero; /* procedure failed */
}

/* return 0 if fails, 1 else */
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=lgef(p)-3,enh,normF,normG,limite=stack_lim(ltop,1);

  FF=*F; HH=H;
  GG=poldivres(p,*F,&r);
  normF=gexpo(FF);
  normG=gexpo(GG);
  enh=gexpo(H); if (enh<0) enh=0;
  shiftbitprec=normF+2*normG+enh+(long) (4.*log2((double)n)+gamma) +1;
  shiftbitprec2=enh+2*(normF+normG)+(long) (2.*gamma+5.*log2((double)n))+1;
  bitprec2=bitprec+shiftbitprec;
  error=gexpo(r);
  if (error<-bitprec) error=1-bitprec;
  for (i=0; (error>-bitprec && i<NEWTON_MAX) && error<=0; i++)
  {
    if ((bitprec1==bitprec2) && (i>=2))
    {
      shiftbitprec+=n; shiftbitprec2+=n; bitprec2+=n;
    }
    if (low_stack(limite, stack_lim(ltop,1)))
    {
      GEN *gptr[4];
      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 */

    bitprec1=-error+shiftbitprec;
    r=gmul(mygprec(HH,bitprec1),mygprec(r,bitprec1));
    r=mygprec(r,bitprec1);
    f0=gres(r,mygprec(FF,bitprec1));

    bitprec1=-2*error+shiftbitprec;
    if (bitprec1>bitprec2) bitprec1=bitprec2;
    FF=gadd(mygprec(FF,bitprec1),f0);

    bitprec1=-3*error+shiftbitprec;
    if (bitprec1>bitprec2) bitprec1=bitprec2;
    pp=mygprec(p,bitprec1);
    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 */
}

/* returns F and G from the unit circle U such that |p-FG|<2^(-bitprec) |cd|,
where cd is the leading coefficient of p */
static void
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=lgef(p)-3,NN,bitprec2,
  ltop=avma,polreal=isreal(p);
  double mu,gamma;

  pp=gdiv(p,(GEN)p[2+n]);
  Lmax=4; while (Lmax<=n) Lmax=(Lmax<<1);
  parameters(pp,&mu,&gamma,polreal,param,param2);

  H =cgetg(k+2,t_POL); H[1] =evalsigne(1) | evalvarn(varn(p)) | evallgef(k+2);
  FF=cgetg(k+3,t_POL); FF[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(k+3);
  FF[k+2]=un;

  NN=(long) (0.5/delta); NN+=(NN%2); if (NN<2) NN=2;
  NN=NN*Lmax; ltop=avma;
  for(;;)
  {
    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;
  }
  *G=gmul(GG,(GEN)p[2+n]); *F=FF;
}

static void
optimize_split(GEN p, long k, double delta, long bitprec,
            GEN *F, GEN *G, double param, double param2)
{
  long n=lgef(p)-3;
  GEN FF,GG;

  if (k<=n/2)
    split_fromU(p,k,delta,bitprec,F,G,param,param2);
  else
  { /* start from the reciprocal of p */
    split_fromU(polrecip_i(p),n-k,delta,bitprec,&FF,&GG,param,param2);
    *F=polrecip(GG); *G=polrecip(FF);
  }
}

/********************************************************************/
/**                                                                **/
/**             RECHERCHE DU CERCLE DE SEPARATION                  **/
/**                                                                **/
/********************************************************************/

/* return p(2^e*x) *2^(-n*e) */
static void
scalepol2n(GEN p, long e)
{
  long i,n=lgef(p)-1;
  for (i=2; i<=n; i++) p[i]=lmul2n((GEN)p[i],(i-n)*e);
}

/* returns p(x/R)*R^n */
static GEN
scalepol(GEN p, GEN R, long bitprec)
{
  GEN q,aux,gR;
  long i;

  aux = gR = mygprec(R,bitprec); q = mygprec(p,bitprec);
  for (i=lgef(p)-2; i>=2; i--)
  {
    q[i]=lmul(aux,(GEN)q[i]);
    aux = gmul(aux,gR);
  }
  return q;
}

/* returns q(x)=p(b), where b is (usually) a polynomial  */
static GEN
shiftpol(GEN p, GEN b)
{
  long av = avma,i, limit = stack_lim(av,1);
  GEN q = gzero;

  for (i=lgef(p)-1; i>=2; i--)
  {
    q = gadd((GEN)p[i], gmul(q,b));
    if (low_stack(limit, stack_lim(av,1))) q = gerepileupto(av,q);
  }
  return gerepileupto(av,q);
}

/* return (aX-1)^n * p[ (X-a) / (aX-1) ] */
static GEN
conformal_pol(GEN p, GEN a, long bitprec)
{
  GEN r,pui,num,aux;
  long n=lgef(p)-3, i;

  aux = pui = cgetg(4,t_POL);
  pui[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(4);
  pui[2] = (long) mygprec(gneg_i(gunr),bitprec);
  pui[3] = lconj(a);
  num = cgetg(4,t_POL);
  num[1] = pui[1];
  num[2] = lneg(a);
  num[3] = (long) mygprec(gunr,bitprec);
  r=(GEN)p[2+n];
  for (i=n-1; ; i--)
  {
    r=gadd(gmul(r,num),gmul(aux,(GEN) p[2+i]));
    if (i==0) return r;
    aux=gmul(pui,aux);
  }
}

/* apply the conformal mapping then split from U */
static void
conformal_mapping(GEN p, long k, long bitprec, double aux,GEN *F,GEN *G)
{
  long bitprec2,bitprec3,n=lgef(p)-3,decprec,i,ltop = avma, av;
  GEN q,FF,GG,a,R, *gptr[2];
  double rmin,rmax,rho,delta,aux2,param,param2,r1,r2;

  bitprec2=bitprec+(long) (n*(2.*log2(2.732)+log2(1.5)))+1;
  a=gsqrt(stoi(3),10);
  a=gmul(mygprec(a,bitprec2),mygprec(globalu,bitprec2));
  a=gdivgs(a,-6); /* a=-globalu/2/sqrt(3) */

  av = avma; q=mygprec(p,bitprec2);
  q = conformal_pol(q,a,bitprec2);
  for (i=1; i<=n; i++)
    if (radius[i]!=0.) /* updating array radius */
    {
      aux2=radius[i]*radius[i];
      aux2=2.*(aux2-1)/(aux2-3.*radius[i]+3.);
      radius[i]=sqrt(1+aux2);
    }
  if (k>1)
  {
    i=k-1; while (i>0 && radius[i]==0.) i--;
    r1=radius[i]; r2=radius[k];
    rmin=pre_modulus(q,k,aux/10.,r1,r2);
  }
  else
    rmin=min_modulus(q,aux/10.);
  radius[k]=rmin;

  if (k+1<n)
  {
    i=k+2; while (i<=n && radius[i]==0.) i++;
    r2=radius[i]; r1=radius[k+1];
    rmax=pre_modulus(q,k+1,aux/10.,r1,r2);
  }
  else /* k+1=n */
    rmax=max_modulus(q,aux/10.);
  radius[k+1]=rmax;

  rho=sqrt(rmin*rmax); delta=0.5*log(rmax/rmin);
  if (delta>1.) delta=1.;

  if (rho<1.) bitprec3=(long) ((double)n*log2(1./rho))+1;
  else bitprec3=(long) ((double)n*log2(rho))+1;
  bitprec2=bitprec2+bitprec3;

  R=mygprec(dbltor(1/rho),bitprec2);
  q=scalepol(q,R,bitprec2);
  gptr[0]=&q; gptr[1]=&R; gerepilemany(av,gptr,2);

  aux2=radius[k];
  for (i=k-1; i>=1; i--)
  {
    if (radius[i]==0.) radius[i]=aux2;
    else aux2=radius[i];
  }
  aux2=radius[k+1];
  for (i=k+1; i<=n; i++)
  {
    if (radius[i]==0.) radius[i]=aux2;
    else aux2=radius[i];
  }
  param=0.; param2=0.;
  for (i=1; i<=n; i++)
  {
    radius[i]=radius[i]/rho;
    aux2=fabs(1-radius[i]);
    param+=1./aux2;
    if (aux2<1.) param2-=log2(aux2);
  }
  optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
  bitprec2=bitprec2+n; R = ginv(R);
  FF=scalepol(FF,R,bitprec2);
  GG=scalepol(GG,R,bitprec2);

  a=mygprec(a,bitprec2);
  FF=conformal_pol(FF,a,bitprec2);
  GG=conformal_pol(GG,a,bitprec2);
  a=ginv(gsub(gun,gmul(a,gconj(a))));
  a=glog(a,(long) (bitprec2 * L2SL10)+1);

  decprec=(long) ((bitprec+n) * L2SL10)+1;
  FF=gmul(FF,gexp(gmulgs(a,k),decprec));
  GG=gmul(GG,gexp(gmulgs(a,n-k),decprec));

  *F=mygprec(FF,bitprec+n); *G=mygprec(GG,bitprec+n);
  gptr[0]=F; gptr[1]=G; gerepilemany(ltop,gptr,2);
}

/* split p, this time with no scaling. returns in F and G two polynomials
such that |p-FG|< 2^(-bitprec)|p| */
static void
split_2(GEN p, long bitprec, double thickness, GEN *F, GEN *G)
{
  double rmin,rmax,rho,kappa,aux,delta,param,param2,r1,r2;
  long n=lgef(p)-3,i,j,k,disti,diste,bitprec2;
  GEN q,FF,GG,R;

  radius=(double *) gpmalloc((n+1)*sizeof(double));
  for (i=2; i<n; i++) radius[i]=0.;
  rmin=min_modulus(p,thickness/(double) n/4.);
  rmax=max_modulus(p,thickness/(double) n/4.);
  i=1; j=n; radius[1]=rmin; radius[n]=rmax;
  rho=sqrt(rmin*rmax);
  k=dual_modulus(p,rho,thickness/(double) n/4.,1);
  if (k<n/5. || (n/2.<k && k<(4*n)/5.)) { rmax=rho; j=k+1; radius[j]=rho; }
  else { rmin=rho; i=k; radius[i]=rho; }
  while (j>i+1)
  {
    disti= (i<n-j)? i: n-j; /* min(i,n-j) */
    diste= (j<n-i)? j: n-i;
    kappa=1.-log(1.+(double)disti)/log(1.+(double)diste);
    if (i+j<n+1)
      rho=exp( (log(rmin)+log(rmax)*(1+kappa))/(2+kappa));
    else if (i+j==n+1) rho=sqrt(rmin*rmax);
    else
      rho=exp( (log(rmin)*(1+kappa)+log(rmax))/(2+kappa));
    /* use log(rmax) - log(rmin) since rmax / rmin can overflow */
    aux=(log(rmax)-log(rmin))/(j-i)/4.;

    k=(i<n+1-j)? i: n+1-j; /* min(i,n+1-j) */
    k=dual_modulus(p,rho,aux,k);
    if (k-i < j-k-1) { rmax=rho; j=k+1; radius[j]=rho*exp(-aux); }
    else
    {
      if (k-i > j-k-1) { rmin=rho; i=k; radius[i]=rho*exp(aux); }
      else
      {
	if (2*k>n) { rmax=rho; j=k+1; radius[j]=rho*exp(-aux); }
	else { rmin=rho; i=k; radius[i]=rho*exp(aux); }
      }
    }
  }
  aux=log(rmax)-log(rmin);

  if (step4==0)
  {
    if (k>1)
    {
      i=k-1; while ((i>0) && (radius[i]==0.)) i--;
      r1=radius[i]; r2=radius[k];
      rmin=pre_modulus(p,k,aux/10.,r1,r2);
    }
    else /* k=1 */
      rmin=min_modulus(p,aux/10.);
    radius[k]=rmin;

    if (k+1<n)
    {
      i=k+2; while ((i<=n) && (radius[i]==0.)) i++;
      r2=radius[i]; r1=radius[k+1];
      rmax=pre_modulus(p,k+1,aux/10.,r1,r2);
    }
    else /* k+1=n */
      rmax=max_modulus(p,aux/10.);
    radius[k+1]=rmax;
    rho=sqrt(rmin*rmax); delta=0.5*(log(rmax)-log(rmin));

    aux=radius[k];
    for (i=k-1; i>=1; i--)
    {
      if (radius[i]==0.) radius[i]=aux;
      else aux=radius[i];
    }
    aux=radius[k+1];
    for (i=k+1; i<=n; i++)
    {
      if (radius[i]==0.) radius[i]=aux;
      else aux=radius[i];
    }
    param=0.; param2=0.;
    for (i=1; i<=n; i++)
    {
      radius[i]=radius[i]/rho;
      aux=fabs(1.-radius[i]);
      param+=1./aux;
      if (aux<1) param2-=log2(aux);
    }
    if (rho<1.) bitprec2=(long) ((double)n*log2(1./rho))+1;
    else bitprec2=(long) ((double)n*log2(rho))+1;
    bitprec2 += bitprec;

    R=mygprec(dbltor(1./rho),bitprec2);
    q=scalepol(p,R,bitprec2);
    optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
    bitprec2 += n; R=ginv(R);
  }
  else
  {
    rho=sqrt(rmax*rmin);
    if (rho<1.) bitprec2=(long) ((double)n*log2(1./rho))+1;
    else bitprec2=(long) ((double)n*log2(rho))+1;
    bitprec2=bitprec+bitprec2;

    R=mygprec(dbltor(1/rho),bitprec2);
    q=scalepol(p,R,bitprec2);
    for (i=1; i<=n; i++)
      if (radius[i]!=0.) radius[i]=radius[i]/rho;

    conformal_mapping(q, k, bitprec2, aux, &FF, &GG);
    bitprec2 += n; R = ginv(mygprec(R,bitprec2));
  }
  free(radius);
  FF=scalepol(FF,R,bitprec2); GG=scalepol(GG,R,bitprec2);
  *F=mygprec(FF,bitprec+n); *G=mygprec(GG,bitprec+n);
}

/* procedure corresponding to steps 5,6,.. page 44 in the RR n. 1852 */
/* put in F and G two polynomial such that |p-FG|<2^(-bitprec)|p|
where the maximum modulus of the roots of p is <=1 and the sum of roots
is zero */

static void
split_1(GEN p, long bitprec, GEN *F, GEN *G)
{
  long bitprec2,bitprec3,i,imax,n=lgef(p)-3, polreal = isreal(p);
  double rmax,rmin,thickness,quo;
  GEN q,qq,newq,FF,GG,v,gr,r;

  r=gmax_modulus(p,0.01);
  bitprec2=bitprec+n;
  gr=mygprec(ginv(r),bitprec2);
  q=scalepol(p,gr,bitprec2);

  bitprec2=bitprec+gexpo(q)-gexpo(p);

  bitprec3=bitprec2+(long)((double)n*2.*log2(3.)+1);
  v=cgetg(5,t_VEC); v++;
  v[0]=lmulgs(mygprec(gunr,bitprec3),2); v[1]=lneg((GEN)v[0]);
  v[2]=lmul((GEN)v[0],gi); v[3]=lneg((GEN)v[2]);
  imax = polreal? 3: 4;
  q=mygprec(q,bitprec3); thickness=1.;
  for (i=0; i<imax; i++)
  {
    qq=shiftpol(q,gadd(polx[varn(p)],(GEN)v[i]));
    rmin=min_modulus(qq,0.05);
    if (3./rmin > thickness)
    {
      rmax=max_modulus(qq,0.05); quo = rmax/rmin;
      if ((float)quo > (float)thickness)
      {
	thickness=quo; newq=qq; globalu=(GEN)v[i];
      }
    }
    if (thickness>2.) break;
    if (polreal && (i==1 && thickness>1.5)) break;
  }
  bitprec3=bitprec2+(long)((double) n*log2(3.)+1)+gexpo(newq)-gexpo(q);
  split_2(newq,bitprec3,log(thickness),&FF,&GG);
  globalu=gsub(polx[varn(p)],mygprec(globalu,bitprec3));
  FF=shiftpol(FF,globalu); GG=shiftpol(GG,globalu);

  gr=ginv(gr);
  bitprec2=bitprec2+gexpo(FF)+gexpo(GG)-gexpo(q);
  *F=scalepol(FF,gr,bitprec2); *G=scalepol(GG,gr,bitprec2);
}

/* put in F and G two polynomials such that |P-FG|<2^(-bitprec)|P|,
where the maximum modulus of the roots of p is <0.5 */
static void
split_0_2(GEN p, long bitprec, GEN *F, GEN *G)
{
  GEN q,b,FF,GG;
  long n=lgef(p)-3,k,bitprec2,i;
  double auxnorm;

  step4=1;
  auxnorm=(double) n*log2(1+exp2(mylog2((GEN)p[n+1])
				 -mylog2((GEN)p[n+2]))/(double)n);
  bitprec2=bitprec+1+(long) (log2((double)n)+auxnorm);

  q=mygprec(p,bitprec2);
  b=gdivgs(gdiv((GEN)q[n+1],(GEN)q[n+2]),-n);
  q=shiftpol(q,gadd(polx[varn(p)],b));

  k=0; while
      (gexpo((GEN)q[k+2])<-(bitprec2+2*(n-k)+gexpo(q))
       || gcmp0((GEN)q[k+2])) k++;
  if (k>0)
  {
    if (k>n/2) k=n/2;
    bitprec2+=(k<<1);
    FF=cgetg(k+3,t_POL); FF[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(k+3);
    for (i=0; i<k; i++) FF[i+2]=zero;
    FF[k+2]=(long) mygprec(gunr,bitprec2);
    GG=cgetg(n-k+3,t_POL); GG[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(n-k+3);
    for (i=0; i<=n-k; i++) GG[i+2]=q[i+k+2];
    b=gsub(polx[varn(p)],mygprec(b,bitprec2));
  }
  else
  {
    split_1(q,bitprec2,&FF,&GG);
    bitprec2=bitprec+gexpo(FF)+gexpo(GG)-gexpo(p)+(long)auxnorm+1;
    b=gsub(polx[varn(p)],mygprec(b,bitprec2));
    FF=mygprec(FF,bitprec2);
  }
  GG=mygprec(GG,bitprec2);
  *F=shiftpol(FF,b); *G=shiftpol(GG,b);
}

/* put in F and G two polynomials such that |P-FG|<2^(-bitprec)|P|,
where the maximum modulus of the roots of p is <2 */
static void
split_0_1(GEN p, long bitprec, GEN *F, GEN *G)
{
  GEN q,FF,GG;
  long n=lgef(p)-3,bitprec2,normp,pow;
  double aux;

  normp=gexpo(p);
  aux=exp2(mylog2((GEN)p[n+1])-mylog2((GEN)p[n+2]))/(double)n;
  pow=2;
  if (aux<=2.5){ split_0_2(p,bitprec,F,G); return; }

  if (aux<=1.) pow=1;
  scalepol2n(p,pow); /* p <- 4^(-n) p(4*x) */
  bitprec2=bitprec+pow*n+gexpo(p)-normp;
  q=mygprec(p,bitprec2);
  split_1(q,bitprec2,&FF,&GG);
  scalepol2n(FF,-pow); scalepol2n(GG,-pow);
  bitprec2=bitprec+gexpo(FF)+gexpo(GG)-normp;
  *F=mygprec(FF,bitprec2); *G=mygprec(GG,bitprec2);
}

/* put in F and G two polynomials such that |P-FG|<2^(-bitprec)|P| */
static void
split_0(GEN p, long bitprec, GEN *F, GEN *G)
{
  GEN FF,GG,q,R;
  long n=lgef(p)-3,k=0,i;

  while (gexpo((GEN)p[k+2])<-bitprec) k++;
  if (k>0)
  {
    if (k>n/2) k=n/2;
    FF=cgetg(k+3,t_POL);
    FF[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(k+3);
    for (i=0; i<k; i++) FF[i+2]=lcopy(gzero);
    FF[k+2]=(long) mygprec(gunr,bitprec);
    GG=cgetg(n-k+3,t_POL);
    GG[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(n-k+3);
    for (i=0; i<=n-k; i++) GG[i+2]=p[i+k+2];
  }
  else
  {
    R = gmax_modulus(p,0.05);
    if (gexpo(R)<1 && gtodouble(R)<1.9) split_0_1(p,bitprec,&FF,&GG);
    else
    {
      q=cgetg(n+3,t_POL); q[1]=p[1];
      for (i=0; i<=n; i++) q[i+2]=p[n-i+2]; /* p^* with copy of ptr */
      R = gmax_modulus(q,0.05);
      if (gexpo(R)<1 && gtodouble(R)<1.9)
      {
	split_0_1(q,bitprec,&FF,&GG);
	FF=polrecip(FF); GG=polrecip(GG);
      }
      else
      {
	step4=0;
	split_2(p,bitprec,1.2837,&FF,&GG);
      }
    }
  }
  *F=FF; *G=GG;
}

/********************************************************************/
/**                                                                **/
/**     CALCUL A POSTERIORI DE L'ERREUR ABSOLUE SUR LES RACINES    **/
/**                                                                **/
/********************************************************************/

static GEN
root_error(long n, long k, GEN roots_pol, GEN sigma, GEN shatzle)
{
  GEN rho,d,eps,epsbis,eps2,prod,aux,rap=NULL;
  long i,j,m;

  d=cgetg(n+1,t_VEC);
  for (i=1; i<=n; i++)
  {
    if (i!=k)
    {
      aux=gsub((GEN)roots_pol[i],(GEN)roots_pol[k]);
      d[i]=(long) gabs(mygprec(aux,31),4);
    }
  }
  rho=gabs(mygprec((GEN)roots_pol[k],31),4);
  if (gcmp(rho,dbltor(1.))==-1) rho=gun;
  eps=gmul(rho,shatzle);
  aux=gmul(gpui(rho,stoi(n),4),sigma);

  for (j=1; j<=2 || (j<=5 && gcmp(rap,dbltor(1.2))==1); j++)
  {
    m=n; prod=gun;
    epsbis=gdivgs(gmulgs(eps,5),4);
    for (i=1; i<=n; i++)
    {
      if (i!=k && gcmp((GEN)d[i],epsbis)==1)
      {
	m--;
	prod=gmul(prod,gsub((GEN)d[i],eps));
      }
    }
    eps2=gdiv(gmul2n(aux,2*m-2),prod);
    eps2=gpui(eps2,dbltor(1./m),4);
    rap=gdiv(eps,eps2); eps=eps2;
  }
  return eps;
}

/* round a complex or real number x to an absolute value of 2^(-e) */
static GEN
mygprec_absolute(GEN x, long bitprec)
{
  long tx=typ(x),e;
  GEN y;

  switch(tx)
  {
    case t_REAL:
      e=gexpo(x);
      if (e<-bitprec || !signe(x)) { y=dbltor(0.); setexpo(y,-bitprec); }
      else y=mygprec(x,bitprec+e);
      break;
    case t_COMPLEX:
      if (gexpo((GEN)x[2])<-bitprec)
	y=mygprec_absolute((GEN)x[1],bitprec);
      else
      {
	y=cgetg(3,t_COMPLEX);
	y[1]=(long) mygprec_absolute((GEN)x[1],bitprec);
	y[2]=(long) mygprec_absolute((GEN)x[2],bitprec);
      }
      break;

    default: y=mygprec(x,bitprec);
  }
  return y;
}

static long
a_posteriori_errors(GEN p, GEN roots_pol, long err)
{
  GEN sigma,overn,shatzle,x;
  long i,n=lgef(p)-3,e,e_max;

  sigma=mygprec(dbltor(1.),10);
  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),
		    gpui(sigma,overn,0)));
  shatzle=gmul2n(shatzle,1); e_max=-pariINFINITY;
  for (i=1; i<=n; i++)
  {
    x=root_error(n,i,roots_pol,sigma,shatzle);
    e=gexpo(x); if (e>e_max) e_max=e;
    roots_pol[i]=(long) mygprec_absolute((GEN)roots_pol[i],-e);
  }
  return e_max;
}

/********************************************************************/
/**                                                                **/
/**                           MAIN                                 **/
/**                                                                **/
/********************************************************************/

/* compute roots in 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 *iroots)
{
  long n=lgef(p)-3,decprec,ltop,lbot;
  GEN F,G,a,b,m1,m2,m;
  GEN *gptr[2];

  if (n==1)
  {
    a=gneg_i(gdiv((GEN)p[2],(GEN)p[3]));
    (*iroots)++; (*roots_pol)[*iroots]=(long) a;
    return p;
  }
  ltop=avma;
  if (n==2)
  {
    F=gsub(gsqr((GEN)p[3]),gmul2n(gmul((GEN)p[2],(GEN)p[4]),2));
    decprec=(long) ((double) bitprec * L2SL10)+1;
    F=gsqrt(F,decprec);
    a=gneg_i(gdiv(gadd((GEN)p[3],F),gmul2n((GEN)p[4],1)));
    b=gdiv(gsub(F,(GEN)p[3]),gmul2n((GEN)p[4],1));

    gptr[0]=&a; gptr[1]=&b;
    gerepilemany(ltop,gptr,2);
    (*iroots)++; (*roots_pol)[*iroots]=(long) a;
    (*iroots)++; (*roots_pol)[*iroots]=(long) b;
    m=gmul(gsub(polx[varn(p)],mygprec(a,3*bitprec)),
	   gsub(polx[varn(p)],mygprec(b,3*bitprec)));
    return gmul(m,(GEN)p[4]);
  }
  split_0(p,bitprec,&F,&G);
  m1=split_complete(F,bitprec,roots_pol,iroots);
  m2=split_complete(G,bitprec,roots_pol,iroots);
  lbot=avma;
  m=gmul(m1,m2); *roots_pol=gcopy(*roots_pol);
  gptr[0]=roots_pol; gptr[1]=&m;
  gerepilemanysp(ltop,lbot,gptr,2); return m;
}

/* compute a bound on the maximum modulus of roots of p */
static GEN
cauchy_bound(GEN p)
{
  long i,n=lgef(p)-3;
  GEN x=gzero,y,lc;

  lc=gabs((GEN)p[n+2],DEFAULTPREC); /* leading coefficient */
  lc=gdiv(dbltor(1.),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;
  }
  return x;
}

static GEN
mygprecrc_special(GEN x, long bitprec, long e)
{
  long tx=typ(x),lx,ex;
  GEN y;

  if (bitprec<=0) bitprec=0; /* should not happen */
  switch(tx)
  {
    case t_REAL:
      lx=bitprec/BITS_IN_LONG+3;
      if (lx<lg(x)) lx=lg(x);
      y=cgetr(lx); affrr(x,y); ex=-bitprec+e;
      if (!signe(x) && expo(x)>ex) setexpo(y,ex);
      break;
    case t_COMPLEX:
      y=cgetg(3,t_COMPLEX);
      y[1]=(long) mygprecrc_special((GEN)x[1],bitprec,e);
      y[2]=(long) mygprecrc_special((GEN)x[2],bitprec,e);
      break;
    default: y=gcopy(x);
  }
  return y;
}

/* like mygprec but keep at least the same precision as before */
static GEN
mygprec_special(GEN x, long bitprec)
{
  long tx=typ(x),lx,i,e;
  GEN y;

  switch(tx)
  {
    case t_POL:
      lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1]; e=gexpo(x);
      for (i=2; i<lx; i++) y[i]=(long) mygprecrc_special((GEN)x[i],bitprec,e);
      break;

    default: y=mygprecrc_special(x,bitprec,0);
  }
  return y;
}

static GEN
all_roots(GEN p, long bitprec)
{
  GEN q,roots_pol,m;
  long bitprec2,n=lgef(p)-3,iroots,i,j,e,av,tetpil;

  roots_pol=cgetg(n+1,t_VEC); av=avma;
  for (i=1; i<=n; i++) roots_pol[i]=zero;
  for (j=1; j<=10; j++)
  {
    iroots=0; e = 2*gexpo(cauchy_bound(p)); if (e<0) e=0;
    bitprec2=bitprec+(1<<j)*n+gexpo(p)-gexpo((GEN)p[2+n])
      +(long) log2(n)+1+e;
    q=gmul(mygprec(gunr,bitprec2),mygprec(p,bitprec2));
    m=split_complete(q,bitprec2,&roots_pol,&iroots);
    e = gexpo(gsub(mygprec_special(p,bitprec2),m))
      - gexpo((GEN)q[2+n])+(long) log2((double)n)+1;
    if (e<-2*bitprec2) e=-2*bitprec2; /* to avoid e=-pariINFINITY */
    if (a_posteriori_errors(q,roots_pol,e) < -bitprec) return roots_pol;

    tetpil=avma; roots_pol=gerepile(av,tetpil,gcopy(roots_pol));
  }
  err(bugparier,"all_roots");
  return NULL; /* not reached */
}

/* true if x is an exact scalar, that is integer or rational */
static int
isexactscalar(GEN x)
{
  long tx=typ(x);
  return (tx==t_INT || is_frac_t(tx));
}

static int
isexactpol(GEN p)
{
  long i,n=lgef(p)-3;

  for (i=0; i<=n; i++)
    if (isexactscalar((GEN)p[i+2])==0) return 0;
  return 1;
}

static long
isvalidcoeff(GEN x)
{
  long tx=typ(x);

  switch(tx)
  {
    case t_INT: case t_REAL: case t_FRAC: case t_FRACN: return 1;
    case t_COMPLEX:
      if (isvalidcoeff((GEN)x[1]) && isvalidcoeff((GEN)x[2])) return 1;
  }
  return 0;
}

static long
isvalidpol(GEN p)
{
  long i,n = lgef(p);
  for (i=2; i<n; i++)
    if (!isvalidcoeff((GEN)p[i])) return 0;
  return 1;
}

static GEN
solve_exact_pol(GEN p, long bitprec)
{
  GEN S,ex,factors,roots_pol,roots_fact;
  long i,j,k,m,n,iroots;

  n=lgef(p)-3;

  iroots=0;
  roots_pol=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) roots_pol[i]=zero;

  S=square_free_factorization(p);
  ex=(GEN) S[1]; factors=(GEN) S[2];
  for (i=1; i<lg(factors); i++)
  {
    roots_fact=all_roots((GEN)factors[i],bitprec);
    n=lgef(factors[i])-3; m=itos((GEN)ex[i]);
    for (j=1; j<=n; j++)
      for (k=1; k<=m; k++) roots_pol[++iroots] = roots_fact[j];
  }
  return roots_pol;
}

/* return the roots of p with absolute error bitprec */
static GEN
roots_com(GEN p, long l)
{
  long bitprec;

  if (gcmp0(p)) err(zeropoler,"roots");
  if (typ(p)!=t_POL)
  {
    if (!isvalidcoeff(p)) err(typeer,"roots");
    return cgetg(1,t_VEC); /* constant polynomial */
  }
  if (!isvalidpol(p)) err(talker,"invalid coefficients in roots");
  if (lgef(p) == 3) return cgetg(1,t_VEC); /* constant polynomial */
  if (l<3) l=3;
  bitprec=bit_accuracy(l); gunr=realun(l);
  return isexactpol(p)? solve_exact_pol(p,bitprec): all_roots(p,bitprec);
}

static GEN
tocomplex(GEN x, long l)
{
  GEN y=cgetg(3,t_COMPLEX);

  y[1]=lgetr(l);
  if (typ(x) == t_COMPLEX)
    { y[2]=lgetr(l); gaffect(x,y); }
  else
    { gaffect(x,(GEN)y[1]); y[2]=(long)realzero(l); }
 return y;
}

/* Check if x is approximately real with precision e */
int
isrealappr(GEN x, long e)
{
  long tx=typ(x),lx,i;
  switch(tx)
  {
    case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
      return 1;
    case t_COMPLEX:
      return (gexpo((GEN)x[2]) < e);
    case t_QUAD:
      err(impl,"isrealappr for type t_QUAD");
    case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
    case t_VEC: case t_COL: case t_MAT:
      lx = (tx==t_POL)?lgef(x): lg(x);
      for (i=lontyp[tx]; i<lx; i++)
        if (! isrealappr((GEN)x[i],e)) return 0;
      return 1;
    default: err(typeer,"isrealappr"); return 0;
  }
}

/* x,y sont de type t_COMPLEX */
static int
isconj(GEN x, GEN y, long e)
{
  return (gexpo( gsub((GEN)x[1],(GEN)y[1]) ) < e
       && gexpo( gadd((GEN)x[2],(GEN)y[2]) ) < e);
}

/* returns the vector of roots of p, with guaranteed absolute error
 * BASE^(-(l-2)), where BASE=2^{ BITS_IN_LONG }.
 */
GEN
roots(GEN p, long l)
{
  long av,av1,tetpil,n,j,k,s, e = 5 - bit_accuracy(l);
  GEN p1,p2,p3,p22,res;

  av=avma; p1=roots_com(p,l); n=lg(p1);
  if (isrealappr(p,e))
  {
    p3=cgetg(n,t_COL); s=0;
    for (j=1; j<n; j++)
    {
      p2=(GEN)p1[j];
      if (typ(p2) != t_COMPLEX) { p3[++s]=(long)p2; p1[j]=zero; }
      else if (isrealappr(p2,e)) { p3[++s]=p2[1]; p1[j]=zero; }
    }
    setlg(p3,s+1); p2=sort(p3); setlg(p3,n);
    tetpil=avma; res=cgetg(n,t_COL);
    for (j=1; j<=s; j++) res[j]=(long)tocomplex((GEN)p2[j],l);
    for (j=1; j<n; j++)
    {
      p2=(GEN)p1[j];
      if (typ(p2) == t_COMPLEX)
      {
	res[++s]=(long)tocomplex(p2,l);
        av1=avma;
	for (k=j+1; k<n; k++)
	{
	  p22=(GEN)p1[k];
	  if (typ(p22) == t_COMPLEX && isconj(p2,p22,e))
          {
            avma=av1; res[++s]=(long)tocomplex(p22,l);
            p1[k]=zero; break;
          }
	}
	if (k==n) err(bugparier,"roots (conjugates)");
      }
    }
    return gerepile(av,tetpil,res);
  }
  tetpil=avma; res=cgetg(n,t_COL);
  for (j=1; j<n; j++) res[j]=(long)tocomplex((GEN)p1[j],l);
  return gerepile(av,tetpil,res);
}

GEN
roots0(GEN p, long flag,long l)
{
  switch(flag)
  {
    case 0: return roots(p,l);
    case 1: return rootsold(p,l);
    default: err(flagerr,"polroots");
  }
  return NULL; /* not reached */
}