// -*- mode:C++ ; compile-command: "g++ -I.. -I../include -g -c -Wall -DHAVE_CONFIG_H -DIN_GIAC csturm.cc" -*-
#include "giacPCH.h"

/*
 *  Copyright (C) 2000,14 B. Parisse, Institut Fourier, 38402 St Martin d'Heres
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
using namespace std;
#include <cmath>
#include <stdexcept>
#include <map>
#include "gen.h"
#include "csturm.h"
#include "vecteur.h"
#include "modpoly.h"
#include "unary.h"
#include "symbolic.h"
#include "usual.h"
#include "sym2poly.h"
#include "solve.h"
#include "prog.h"
#include "subst.h"
#include "permu.h"
#include "series.h"
#include "alg_ext.h"
#include "ti89.h"
#include "plot.h"
#include "modfactor.h"
#include"giacintl.h"

#ifndef NO_NAMESPACE_GIAC
namespace giac {
#endif // ndef NO_NAMESPACE_GIAC

  // compute Sturm sequence of r0 and r1,
  // returns gcd (without content)
  // and compute list of quotients, coeffP, coeffR
  // such that coeffR*r_(k+2) = Q_k*r_(k+1) - coeffP_k*r_k
  gen csturm_seq(modpoly & r0,modpoly & r1,vecteur & listquo,vecteur & coeffP, vecteur & coeffR,GIAC_CONTEXT){
    listquo.clear();
    coeffP.clear();
    coeffR.clear();
    if (r0.empty())
      return r1;
    if (r1.empty())
      return r0;
    gen tmp;
    lcmdeno(r0,tmp,contextptr);
    if (ck_is_positive(-tmp,contextptr))
      r0=-r0;
    r0=r0/abs(lgcd(r0),contextptr);
    lcmdeno(r1,tmp,contextptr);
    if (ck_is_positive(-tmp,contextptr))
      r1=-r1;
    r1=r1/abs(lgcd(r0),contextptr);
    // set auxiliary constants g and h to 1
    gen g(1),h(1);
    modpoly a(r0),b(r1),quo,r;
    gen b0(1);
    for (int loop_counter=0;;++loop_counter){
      int m=int(a.size())-1;
      int n=int(b.size())-1;
      int ddeg=m-n; // should be 1 generically
      if (!n) { // if b is constant, gcd=1
	return 1;
      }
      b0=b.front();
      if (b.front().type==_VECT) {
	// ddeg should be even if b0 is a _POLY1
	if (ddeg%2==0)
	  *logptr(contextptr) << gettext("Singular parametric Sturm sequence ") << a << "/" << b << endl;
      }
      else
	b0=abs(b.front(),contextptr); 
      coeffP.push_back(pow(b0,ddeg+1));
      DivRem(coeffP.back()*a,b,0,quo,r);
      listquo.push_back(quo);
      coeffR.push_back(g*pow(h,ddeg));
      if (r.empty()){
	return b/abs(lgcd(b),contextptr);
      }
      // remainder is non 0, loop continue: a <- b
      a=b;
      // now divides r by g*h^(m-n) and change sign, result is the new b
      b= -r/coeffR.back();
      g=b0;
      h=pow(b0,ddeg)/pow(h,ddeg-1);
    } // end while loop
  }

  static gen csturm_horner(const modpoly & p,const gen & a){
    if (p.size()==1 && p.front().type==_POLY && p.front()._POLYptr->dim==1){
      // patch for "sparse modpoly"
      vector< monomial<gen> >::const_iterator it=p.front()._POLYptr->coord.begin(),itend=p.front()._POLYptr->coord.end();
      gen res=0,anum=a,aden=1,den=1;
      int oldpui=0,pui;
      if (a.type==_FRAC){
	anum=a._FRACptr->num;
	aden=a._FRACptr->den;
      }
      for (;it!=itend;++it){
	pui=it->index.front();
	if (oldpui){
	  res = res * pow(anum,oldpui-pui,context0);
	  den = den * pow(aden,oldpui-pui,context0);
	}
	res += it->value*den;
	oldpui=pui;
      }
      if (oldpui){
	res = res *pow(a,oldpui,context0);
	den = den *pow(a,oldpui,context0);
      }
      return res/den;
    }
    else
      return horner(p,a);
  }

  static int csturm_vertex_ab(const modpoly & r0,const modpoly & r1,const vecteur & listquo,const vecteur & coeffP, const vecteur & coeffR,const gen & a,int start,GIAC_CONTEXT){
    int n=int(listquo.size()),j,k,res=0;
    vecteur R(n+2);
    R[0]=csturm_horner(r0,a);
    R[1]=csturm_horner(r1,a);
    for (j=0;j<n;j++)
      R[j+2]=(-coeffP[j]*R[j]+R[j+1]*horner(listquo[j],a))/coeffR[j];
    // signes
    for (j=start;j<n+2;j++){
      if (R[j]!=0) break;
    }
    for (k=j+1;k<n+2;k++){
      if (is_zero(R[k])) continue;
      if (fastsign(R[j],context0)*fastsign(R[k],context0)<0 
	  // is_positive(-R[j]*R[k],contextptr)
	  ){
	res++;
	j=k;
      }
    }
    return res;
  }

#if 0
  // compute vertex index at a (==0 unless s(a)==0) 
  static int csturm_vertex_a(const modpoly & s,const modpoly & r,const gen & a,int direction,GIAC_CONTEXT){
    int j;
    modpoly s1,s2;
    gen sa=horner(s,a,0,s1);
    if (!is_zero(sa)) return 0;
    for (j=1;;j++){
      sa=horner(s1,a,0,s2);
      if (!is_zero(sa))
	break;
      s1=s2;
    }
    if (direction==1) j=0;
    gen tmp=sign(sa,contextptr)*sign(horner(r,a),contextptr);
    return tmp.val*((j%2)?-1:1);
  }
#endif

  void change_scale(modpoly & p,const gen & l){
    int n=int(p.size());
    gen lton(l);
    for (int i=n-2;i>=0;--i){
      p[i] = p[i] * lton;
      lton = lton * l;
    }
  }

  void back_change_scale(modpoly & p,const gen & l){
    int n=int(p.size());
    gen lton(l);
    for (int i=n-2;i>=0;--i){
      p[i] = p[i] / lton;
      lton = lton * l;
    }
  }

  // p(x)->p(a*x+b)
  modpoly linear_changevar(const modpoly & p,const gen & a,const gen & b){
    modpoly res(taylor(p,b));
    change_scale(res,a);
    return res;
  }

  // p(a*x+b)->p(x)
  // t=a*x+b -> pgcd(t)=g((t-b)/a)
  modpoly inv_linear_changevar(const modpoly & p,const gen & a,const gen & b){
    gen A=inv(a,context0); 
    gen B=-b/a;
    modpoly res(taylor(p,B));
    change_scale(res,A);
    return res;
  }

  // Find roots of R, S=R' at precision eps, returns number of roots
  // if eps==0 does not compute intervals for roots
  static int csturm_realroots(const modpoly & S,const modpoly & R,const vecteur & listquo,const vecteur & coeffP, const vecteur & coeffR,const gen & a,const gen & b,const gen & t0, const gen & t1,vecteur & realroots,double eps,GIAC_CONTEXT){
    if (is_inf(t0)) // replace with max(R)
      return csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,-linfnorm(R,contextptr),t1,realroots,eps,contextptr);
    if (is_inf(t1)) // replace with max(R)
      return csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,t0,linfnorm(R,contextptr),realroots,eps,contextptr);    
    int n1=csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t0,1,contextptr);
    int n2=csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t1,1,contextptr);
    int n=(n2-n1);
    if (!eps || !n)
      return n;
    /* disabled localization of roots, do isolation of roots instead
    if (is_strictly_greater(eps,(t1-t0)*abs(b,contextptr),contextptr)){
      realroots.push_back(makevecteur(makevecteur(a+t0*b,a+t1*b),n));
      return n;
    }
    */
    if (n==1){
      gen T0=t0,T1=t1,T2;
      int s0=fastsign(csturm_horner(R,T0),contextptr);
      // int s1=fastsign(csturm_horner(R,T1),contextptr);
      int s2;
      gen delta=evalf_double(log((T1-T0)*abs(b,contextptr)/eps,contextptr)/log(2.,contextptr),1,contextptr);
      if (delta.type!=_DOUBLE_){
	realroots=vecteur(1,gentypeerr(contextptr));
	return -2;
      }
      int nstep=int(delta._DOUBLE_val+1);
      for (int step=0;step<nstep;++step){
	T2=(T0+T1)/2;
	s2=fastsign(csturm_horner(R,T2),contextptr);
	if (!s2){
	  realroots.push_back(makevecteur(a+T2*b,n));
	  return n;
	}
	if (s2==s0)
	  T0=T2;
	else
	  T1=T2;
      }
      realroots.push_back(makevecteur(makevecteur(a+T0*b,a+T1*b),n));
      return n;
    }
    gen T0=t0,T1=t1,t01;
    for (;;){
      t01=(T0+T1)/2;
      int n01=csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t01,1,contextptr);
      if (n01!=n1 && n01!=n2)
	break;
      if (n01==n1)
	T0=t01;
      else
	T1=t01;
    }
    if (csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,T0,t01,realroots,eps,contextptr)==-2)
      return -2;
    if (csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,t01,T1,realroots,eps,contextptr)==-2)
      return -2;
    return n;
  }

  // Find complex sturm sequence for P(a+(b-a)*x)
  // If P is "pseudo"-real on [a,b] and eps>0 put roots in [a,b]
  // at precision eps inside realroots
  // returns a,b,R,S,g,listquo,coeffP,coeffR,typeseq
  // with typeseq=0 (complex Sturm) or 1 (limit)
  // If b-a is real and horiz_sturm is not empty, it tries to replace
  // the variable by im(a)*i in horiz_sturm and if no quotient in horiz_sturm
  // has a leading 0 coefficient, 
  // it returns im(a)*i,im(a)*i+1,R,S,g,listquo,coeffP,coeffR,typeseq
  // If b-a is pure imaginary and vert_sturm is not empty, it tries to replace
  // the variable by re(a) and returns re(a),re(a)+i,R,S,g,listquo,coeffP,coeffR,typeseq
  static vecteur csturm_segment_seq(const modpoly & P,const gen & a,const gen & b,vecteur & realroots,double eps,vecteur & horiz_sturm,vecteur & vert_sturm,GIAC_CONTEXT){
    // try with horiz_sturm and vert_sturm
    gen ab(b-a);
    /* // Optimization fails for sturmab(x^3-1,-1-i,1+i)
    if (is_zero(re(ab,contextptr))){ // b-a is pure imaginary
      if (vert_sturm.empty()){
	gen A=gen(makevecteur(1,0),_POLY1__VECT);
	vert_sturm.push_back(undef);
	vecteur tmp;
	vert_sturm=csturm_segment_seq(P,A,A+cst_i,tmp,eps,horiz_sturm,vert_sturm,contextptr);
	if (is_undef(vert_sturm))
	  return vert_sturm;
      }
      if (vert_sturm.size()==9){
	vecteur res(vert_sturm);
	gen A=re(a,contextptr);
	res[0]=A; // re(a)
	res[1]=A+cst_i; // re(a)+i
	res[2]=apply1st(res[2],A,horner); // R 
	res[3]=apply1st(res[3],A,horner); // S
	res[4]=horner(res[4],A); // g
	vecteur tmp(*res[5]._VECTptr);
	int tmps=tmp.size();
	for (int j=0;j<tmps;++j)
	  tmp[j]=apply1st(tmp[j],A,horner); 
	res[5]=tmp; // listquo
	res[6]=apply1st(res[6],A,horner); // coeffP
	res[7]=apply1st(res[7],A,horner); // coeffR
	if (res[6].type==_VECT && !equalposcomp(*res[6]._VECTptr,0))
	  return res;
	else
	  CERR << "list of quotients is not regular" << endl;
      }
    }
    */
    modpoly Q(taylor(P,a));
    change_scale(Q,b-a);
    // now x is in 0..1
    gen gtmp=apply(Q,re,contextptr);
    if (gtmp.type!=_VECT)
      return vecteur(1,gensizeerr(contextptr));
    modpoly R=trim(*gtmp._VECTptr,0);
    gtmp=apply(Q,im,contextptr);
    if (gtmp.type!=_VECT)
      return vecteur(1,gensizeerr(contextptr));
    modpoly S=trim(*gtmp._VECTptr,0);
    modpoly listquo,coeffP,coeffR;
    gen g=csturm_seq(S,R,listquo,coeffP,coeffR,contextptr);
    int typeseq=-1;
    if (debug_infolevel)
      *logptr(contextptr)  << "segment " << a << ".." << b << ", im/re:" << S << "|" << R << ", gcd:" << g << endl;
    if (g.type==_VECT && g._VECTptr->size()==P.size()){
      // if g==P (up to a constant), use real Sturm sequences
      if (debug_infolevel)
	*logptr(contextptr)  << "Real-kind roots: " << g << endl;
      R=*g._VECTptr;
      S=derivative(R);
      g=csturm_seq(S,R,listquo,coeffP,coeffR,contextptr);
      typeseq=csturm_realroots(S,R,listquo,coeffP,coeffR,a,b-a,0,1,realroots,eps,contextptr);
      if (typeseq==-2)
	return realroots;
    }
    if (g.type==_VECT)
      g=inv_linear_changevar(*g._VECTptr,b-a,a);
    vecteur res= makevecteur(a,b,R,S,g,listquo,coeffP,coeffR,typeseq);
    return res;
  }

  // index for segment a,b (2* number of roots when summed over a closed
  // polygon). Note that if S=ImP along the segment is 0 we remove
  // the roots on [a,b] using real Sturm sequences
  // If S=0 at a or b, this is simply ignored
  // Indeed the computed index is then the same as if S was of the
  // sign of R, and since R!=0 if S is 0 this is a property of the vertex
  // not of the segment (note that contrary to counting real roots
  // on an interval, S can vanish as many times as long as R keeps
  // the same sign, without modifying the algebraic number of Im=0
  // cuts if S has the same sign on both end)
  static int csturm_segment(const vecteur & seq,const gen & a,const gen & b,GIAC_CONTEXT){
    gen t0,t1;
    if (seq.size()!=9)
      return -(RAND_MAX/2);
    gen aseq=seq[0];
    gen bseq=seq[1];
    gen directeur=(b-a)/(bseq-aseq);
    t0=(a-aseq)/(bseq-aseq);
    if ( !is_zero(im(directeur,contextptr)) || !is_zero(im(t0,contextptr)) )
      return -(RAND_MAX/2);
    t0=re(t0,contextptr); // t0=normal(t0);
    t1=re(t0+directeur,contextptr); // t1=normal(t0+directeur);
    int signe=1;
    if (is_strictly_greater(t0,t1,contextptr)){
      signe=-1;
      swapgen(t0,t1);
    } 
    const modpoly & R=*seq[2]._VECTptr;
    const modpoly & S=*seq[3]._VECTptr;
    gen g=seq[4];
    const modpoly & listquo=*seq[5]._VECTptr;
    const modpoly & coeffP=*seq[6]._VECTptr;
    const modpoly & coeffR=*seq[7]._VECTptr;
    int debut=(seq[8].val==-1)?0:1;
    int tmp = csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t0,debut,contextptr);
    int res = tmp;
    tmp = csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t1,debut,contextptr);
    res -= tmp;
    // tmp = (-csturm_vertex_a(S,R,t0,1,contextptr)+csturm_vertex_a(S,R,t1,-1,contextptr));
    // res += tmp;
    res=(debut?1:signe)*res;
    if (debug_infolevel)
      *logptr(contextptr)  << "segment " << a << ".." << b << " index contribution " << res << endl;
    return res;
  }

  static bool csturm_square_seq(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,gen & pgcd,vecteur & realroots,double eps,vecteur & seq1,vecteur & seq2,vecteur & seq3,vecteur & seq4,vecteur & horiz_sturm,vecteur & vert_sturm,GIAC_CONTEXT){
    gen A=a0+cst_i*b0,B=a1+cst_i*b0;
    vecteur rroots;
    seq1=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr);
    if (is_undef(seq1))
      return false;
    pgcd=seq1[4];
    if (!is_one(pgcd)){
      return false;
    }
    A=a1+cst_i*b0; B=a1+cst_i*b1;
    seq2=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr);
    if (is_undef(seq2))
      return false;
    pgcd=seq2[4];
    if (!is_one(pgcd)){
      return false;
    }
    A=a1+cst_i*b1; B=a0+cst_i*b1;
    seq3=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr);
    if (is_undef(seq3))
      return false;
    pgcd=seq3[4];
    if (!is_one(pgcd)){
      return false;
    }
    A=a0+cst_i*b1; B=a0+cst_i*b0;
    seq4=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr);
    if (is_undef(seq4))
      return false;
    pgcd=seq4[4];
    if (!is_one(pgcd)){
      return false;
    }
    realroots=mergevecteur(realroots,rroots);
    return true;
  }

  // find 2* number of roots of P inside the square of vertex of affixes a,b
  // roots on the square are not counted. P must not vanish at the vertices.
  // The complex Sturm sequences must be known
  // returns -1 on error
  static int csturm_square(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,const vecteur & seq1,const vecteur & seq2,const vecteur & seq3,const vecteur & seq4,GIAC_CONTEXT){
    int ind,tmp;
    ind = 0;
    gen A=a0+cst_i*b0,B=a1+cst_i*b0;
    tmp = csturm_segment(seq1,A,B,contextptr);
    if (tmp==-(RAND_MAX/2))
      return -1;
    ind += tmp;
    A=a1+cst_i*b0; B=a1+cst_i*b1;
    tmp = csturm_segment(seq2,A,B,contextptr);
    if (tmp==-(RAND_MAX/2))
      return -1;
    ind += tmp;
    A=a1+cst_i*b1; B=a0+cst_i*b1;
    tmp = csturm_segment(seq3,A,B,contextptr);
    if (tmp==-(RAND_MAX/2))
      return -1;
    ind += tmp;
    A=a0+cst_i*b1; B=a0+cst_i*b0;
    tmp = csturm_segment(seq4,A,B,contextptr);
    if (tmp==-(RAND_MAX/2))
      return -1;
    ind += tmp;
    return ind;
  }

  static void csturm_normalize(modpoly & p,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & roots){
    int n=int(p.size())-1;
    // Make sure that x->a+i*x does not return a multiple 
    // of a real polynomial with the multiple non real
    // If degree of p is even the multiple will be a real (because of lcoeff)
    if (n%2){
      // If degree is odd then look at q=p(x-a_{n-1}/n*an)
      // it has the same property
      // if its cst coeff is zero remove
      gen an=p.front();
      gen b=p[1];
      gen shift=-b/n/an;
      modpoly q(taylor(p,shift));
      gen q0;
      // remove valuation
      int qs=int(q.size());
      int n1=0;
      for (;qs>0;--qs,++n1){
	if (!is_zero(q0=q[qs-1]))
	  break;
      }
      if (is_zero(re(q0,context0))){
	q=cst_i*q;
	p=cst_i*p;
      }
      if (n1){
	q=modpoly(q.begin(),q.begin()+qs);
	gen a=re(shift,context0),b=im(shift,context0);
	if (is_greater(a,a0,context0) && is_greater(b,b0,context0) && is_greater(a1,a,context0) && is_greater(b1,b,context0))
	  roots.push_back(makevecteur(shift,n1));
	p=taylor(q,-shift);
      }
    }
  }

  void ab2a0b0a1b1(const gen & a,const gen & b,gen & a0,gen & b0,gen & a1,gen & b1,GIAC_CONTEXT){
    a0=re(a,contextptr); b0=im(a,contextptr);
    a1=re(b,contextptr); b1=im(b,contextptr);
    if (ck_is_greater(a0,a1,contextptr)) swapgen(a0,a1);
    if (ck_is_greater(b0,b1,contextptr)) swapgen(b0,b1);
  }

  // find 2* number of roots of P inside the square of vertex of affixes a,b
  // excluding those on the square
  // returns -1 on error
  int csturm_square(const gen & p,const gen & a,const gen & b,gen& pgcd,GIAC_CONTEXT){
    if (p.type==_POLY){
      int res=0;
      factorization f(sqff(*p._POLYptr));
      factorization::const_iterator it=f.begin(),itend=f.end();
      for (;it!=itend;++it){
	int tmp=csturm_square(polynome2poly1(it->fact),a,b,pgcd,contextptr);
	if (tmp==-1)
	  return -1;
	res += it->mult*tmp;
      }
      return res;
    }
    if (p.type!=_VECT)
      return 0;
    modpoly P=*p._VECTptr;
    vecteur realroots;
    gen a0,b0,a1,b1;
    ab2a0b0a1b1(a,b,a0,b0,a1,b1,contextptr);
    csturm_normalize(P,a0,b0,a1,b1,realroots);
    int evident=0;
    if (!realroots.empty()){
      gen r=realroots.front();
      if (r.type==_VECT && r._VECTptr->size()==2)
	r=r._VECTptr->front();
      gen rx=re(r,contextptr),ry=im(r,contextptr);
      if ( ( is_zero(ry) && (rx==a0 || rx==a1) ) ||
	   ( is_zero(rx) && (ry==b0 || ry==b1) ) )
	;
      else
	evident=1;
    }
    if (P.size()<2)
      return evident;
    vecteur seq1,seq2,seq3,seq4,horiz_seq,vert_seq;    
    if (!csturm_square_seq(P,a0,b0,a1,b1,pgcd,realroots,0.0,seq1,seq2,seq3,seq4,horiz_seq,vert_seq,contextptr)){
      if (pgcd.type!=_VECT)	      
	return -1;
      modpoly g=(*pgcd._VECTptr)/pgcd[0];
      // true factorization found, restart with each factor
      modpoly p1=P/g;
      int n1=csturm_square(p1,a,b,pgcd,contextptr);
      if (n1==-1)
	return -1;
      int n2=csturm_square(g,a,b,pgcd,contextptr);
      if (n2==-1)
	return -1;
      return evident+n1+n2;
    }
    int tmp=csturm_square(P,a0,b0,a1,b1,seq1,seq2,seq3,seq4,contextptr);
    if (tmp==-1)
      return tmp;
    return evident+tmp;
  }

  static void complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps);

  static bool complex_roots_split(const modpoly & P,const gen & pgcd,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps){
    if (pgcd.type!=_VECT)
      return false;
    modpoly g=(*pgcd._VECTptr)/pgcd[0];
    // true factorization found, restart with each factor
    modpoly p1=P/g;
    csturm_normalize(p1,a0,b0,a1,b1,realroots);
    csturm_normalize(g,a0,b0,a1,b1,realroots);
    complex_roots(p1,a0,b0,a1,b1,realroots,complexroots,eps);
    complex_roots(g,a0,b0,a1,b1,realroots,complexroots,eps);
    return true;
  }

#if 0
  // check that arg is >=pi/8 (assumes im(g)>=0)
  static bool arg_geq_pi_8(const gen & g){
    gen gr=re(g,context0),gi=im(g,context0);
    if (is_positive(-gr,context0))
      return true;
    // ? gi/gr>=sqrt(2)-1
    gen r=gi/gr+1;
    if (is_positive(r*r-2,context0))
      return true;
    return false;
  }

  // is im(b/a)>=0, tested without quotient
  static bool arg_in_0_pi(const gen & a,const gen & b){
    gen A(a),B(b);
    if (A.type==_FRAC && is_integer(A._FRACptr->den) && is_positive(A._FRACptr->den,context0))
      A=A._FRACptr->num;
    if (B.type==_FRAC && is_integer(B._FRACptr->den) && is_positive(B._FRACptr->den,context0))
      B=B._FRACptr->num;
    gen c=re(A,context0)*im(B,context0)+re(B,context0)*im(A,context0);
    return is_positive(c,context0);
  }

  static gen hornerarg(const modpoly & p,const gen & x){
    if (p.empty())
      return 0;
    if (x.type!=_FRAC || !is_integer(x._FRACptr->den))
      return horner(p,x);
    fraction & f =*x._FRACptr;
    gen num=f.num,den=f.den,d=den;
    if (is_positive(-f.den,context0)){
      num=-num; den=-den; d=den;
    }
    modpoly::const_iterator it=p.begin(),itend=p.end();
    gen res(*it);
    ++it;
    if (it==itend)
      return res;
    for (;;){
      res=res*num+(*it)*d;
      ++it;
      if (it==itend)
	break;
      d=d*den;   
    }
    return res;
  }  

  // Find one complex root inside a0,b0->a1,b1, return false if not found
  static bool complex_1root(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & complexroots,double eps){
    return false; // disabled since it is not faster!!
    int step,nstep =int(evalf_double(log(max(a1-a0,b1-b0,context0)/eps,context0)/log(2.0,context0),1,context0)._DOUBLE_val+0.5);
    if (nstep<4)
      return false;
    // First compute P at the 4 vertex and check whether P[vertex_n+1]/P[vertex_n] is in C^+
    gen P0=hornerarg(P,a0+cst_i*b0),P2=hornerarg(P,a1+cst_i*b0),
      P4=hornerarg(P,a1+cst_i*b1),P6=hornerarg(P,a0+cst_i*b1);
    if (!(arg_in_0_pi(P0,P2) && arg_in_0_pi(P2,P4) && arg_in_0_pi(P4,P6) && arg_in_0_pi(P6,P0)))
      return false;
    gen A0(a0),A2(a1),B0(b0),B2(b1),A1,B1;
    for (step=0;step<2*nstep;step++){
      A1=(A0+A2)/2;
      B1=(B0+B2)/2;
      gen P1=hornerarg(P,A1+cst_i*B0),P7=hornerarg(P,A0+cst_i*B1),P8=hornerarg(P,A1+cst_i*B1),P3,P5;
      bool found=false;
    /*
      P6(A0,B2) - P5(A1,B2) - P4(A2,B2)  
         |            |           |      
      P7(A0,B1) - P8(A1,B1) - P3(A2,B1) 
         |            |            |      
      P0(A0,B0) - P1(A1,B0) - P2(A2,B0)  
     */
      // ? P0, P1, P8, P7
      if (arg_in_0_pi(P0,P1) && arg_in_0_pi(P1,P8) && arg_in_0_pi(P8,P7) && arg_in_0_pi(P7,P0)){
	A2=A1;
	B2=B1;
	P2=P1;
	P4=P8;
	P6=P7;
	if (step<nstep)
	  continue;
	found=true;
      }
      if (!found){
	P3=hornerarg(P,A2+cst_i*B1);
	// ? P1, P2, P3, P8
	if (arg_in_0_pi(P1,P2) && arg_in_0_pi(P2,P3) && arg_in_0_pi(P3,P8) && arg_in_0_pi(P8,P1)){
	  A0=A1;
	  B2=B1;
	  P0=P1;
	  P4=P3;
	  P6=P8;
	  if (step<nstep)
	    continue;
	  found=true;
	}
      }
      if (!found){
	// P8, P3, P4, P5
	P5=hornerarg(P,A1+cst_i*B2);
	if (arg_in_0_pi(P8,P3) && arg_in_0_pi(P3,P4) && arg_in_0_pi(P4,P5) && arg_in_0_pi(P5,P8)){
	  A0=A1;
	  B0=B1;
	  P0=P8;
	  P2=P3;
	  P6=P5;
	  if (step<nstep)
	    continue;
	  found=true;
	}
      }
      if (!found){
	// P7 P8 P5 P6
	if (arg_in_0_pi(P7,P8) && arg_in_0_pi(P8,P5) && arg_in_0_pi(P5,P6) && arg_in_0_pi(P6,P7)){
	  A2=A1;
	  B0=B1;
	  P0=P7;
	  P2=P8;
	  P4=P5;
	  if (step<nstep)
	    continue;
	  found=true;
	}
      }
      if (!found)
	return false;
      // Final check that there is indeed a root inside rectangle
      // args must be >= pi/8 and degree of (P)*max square length/distance to original square <= pi/8
      gen dist=min(min(A0-a0,a1-A2,context0),min(B0-b0,b1-B2,context0),context0);
      if (is_zero(dist))
	continue;
      gen max_sq=max(A2-A0,B2-B0,context0);
      if (is_greater((int(P.size())-2)*max_sq/dist,cst_pi/8,context0))
	continue;
      gen r1=P2/P0, r2=P4/P2, r3=P6/P4, r4=P0/P6;
      if (arg_geq_pi_8(r1) && arg_geq_pi_8(r2) && arg_geq_pi_8(r3) && arg_geq_pi_8(r4)){
	complexroots.push_back(makevecteur(makevecteur(A0+cst_i*B0,A2+cst_i*B2),1));
	return true;
      }
    }
    return false;
  }
#endif

  static gen round2util(const gen & num,const gen & den,int n){
    if (num.type==_CPLX){
      gen r=round2util(*num._CPLXptr,den,n);
      gen i=round2util(*(num._CPLXptr+1),den,n);
      return r+cst_i*i;
    }
    // num must be a _ZINT
    mpz_t tmp1,tmp2;
    mpz_init_set(tmp1,*num._ZINTptr);
    mpz_mul_2exp(tmp1,tmp1,n+1); // tmp1=2^(n+1)*num
    mpz_add(tmp1,tmp1,*den._ZINTptr); //      + den
    mpz_init_set(tmp2,*den._ZINTptr);
    mpz_mul_ui(tmp2,tmp2,2); // tmp2=2*den
    mpz_fdiv_q(tmp1,tmp1,tmp2);
    gen res=tmp1;
    mpz_clear(tmp1); mpz_clear(tmp2);
    return res;
  }

  void in_round2(gen & x,const gen & deuxn, int n){
    if (x.type==_INT_ || x.type==_ZINT)
      return ;
    if (x.type==_FRAC && x._FRACptr->den.type==_CPLX)
      x=fraction(x._FRACptr->num*conj(x._FRACptr->den,context0),x._FRACptr->den.squarenorm(context0));
    if (x.type==_FRAC && x._FRACptr->den.type==_ZINT && 
	(x._FRACptr->num.type==_ZINT || 
	 (x._FRACptr->num.type==_CPLX && x._FRACptr->num._CPLXptr->type==_ZINT && (x._FRACptr->num._CPLXptr+1)->type==_ZINT)) 
	){
      gen num=x._FRACptr->num,d=x._FRACptr->den;
      x=round2util(num,d,n);
      x=x/deuxn;
      return;
    }
    x=_floor(x*deuxn+plus_one_half,context0)/deuxn;
  }

  void round2(gen & x,int n){
    if (x.type==_INT_ || x.type==_ZINT)
      return ;
    gen deuxn;
    if (n<30)
      deuxn = (1<<n);
    else {
      mpz_t tmp;
      mpz_init_set_si(tmp,1);
      mpz_mul_2exp(tmp,tmp,n);
      deuxn=tmp;
      mpz_clear(tmp);
    }
    in_round2(x,deuxn,n);
  }

  void round2(gen & x,const gen & deuxn,GIAC_CONTEXT){
    if (x.type==_INT_ || x.type==_ZINT)
      return;
    if (x.type!=_FRAC)
      x=_floor(x*deuxn+plus_one_half,context0)/deuxn;
    else {
      gen n=x._FRACptr->num,d=x._FRACptr->den;
      if (d.type==_INT_){
	int di=d.val,ni=1;
	while (di>1){ di=di>>1; ni=ni<<1;}
	if (ni==d.val)
	  return;
      }
      n=2*n*deuxn+d;
      x=iquo(n,2*d)/deuxn;
    }
  }

  // Find one complex root inside a0,b0->a1,b1, return false if not found
  // algo: Newton method in exact mode starting from center
  bool newton_complex_1root(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & complexroots,double eps){
    if (is_positive(a1-a0-0.01,context0) ||
	is_positive(b1-b0-0.01,context0))
      return false;
    gen x0=(a0+a1)/2+cst_i*(b0+b1)/2;
    modpoly Pprime=derivative(P);
    int n=int(-std::log(eps)/std::log(2.0)+.5); // for rounding
    gen eps2=pow(2,-(n+1),context0);
    for (int ii=0;ii<n;ii++){
      gen Pprimex0=horner(Pprime,x0,0,false);
      if (is_zero(Pprimex0,context0))
	return false;
      gen dx=horner(P,x0,0,false)/Pprimex0;
      gen absdx=dx*conj(dx,context0);
      x0=x0-dx;
      gen r=re(x0,context0),i=im(x0,context0);
      if (is_positive(a0-r,context0) || is_positive(r-a1,context0) || 
	  is_positive(b0-i,context0) || is_positive(i-b1,context0))
	return false;
      round2(r,n+4);
      round2(i,n+4);
      x0=r+cst_i*i;
      if (is_positive(absdx-eps2*eps2,context0))
	continue;
      // make a small square around x0 
      // and check that there is indeed a root inside
      gen A0=r-eps2;
      gen A1=r+eps2;
      gen B0=i-eps2;
      gen B1=i+eps2;
      gen tmp;
      if (csturm_square(P,A0+cst_i*B0,A1+cst_i*B1,tmp,context0)==2){
	complexroots.push_back(makevecteur(makevecteur(A0+cst_i*B0,A1+cst_i*B1),1));
	return true;
      }
    }
    return false;
  }

  // Find complex roots of P in a0,b0 -> a1,b1
  static int complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,const vecteur & seq1,const vecteur & seq2,const vecteur & seq3,const vecteur & seq4,vecteur & realroots,vecteur & complexroots,double eps,vecteur & horiz_sturm,vecteur & vert_sturm){
    int n=csturm_square(P,a0,b0,a1,b1,seq1,seq2,seq3,seq4,context0);
    if (debug_infolevel && n)
      CERR << a0 << "," << b0 << ".." << a1 << "," << b1 << ":" << n/2 << endl;
    if (n<=0)
      return 2*n;
    if (eps<=0){
      *logptr(context0) << gettext("Bad precision, using 1e-12 instead of ")+print_DOUBLE_(eps,14) << endl;
      eps=1e-12;
    }
    if (is_strictly_greater(eps,a1-a0,context0) && is_strictly_greater(eps,b1-b0,context0)){
      gen r(makevecteur(a0+cst_i*b0,a1+cst_i*b1));
      complexroots.push_back(makevecteur(r,gen(n)/2));
      return n;
    }
    if (n==2 && newton_complex_1root(P,a0,b0,a1,b1,complexroots,eps))
      return n;
    gen a01=(a0+a1)/2,b01=(b0+b1)/2,pgcd;
    vecteur seqvert,seqhoriz;
    gen A=a0+cst_i*b01,B=a1+cst_i*b01;
    seqhoriz=csturm_segment_seq(P,A,B,realroots,eps,horiz_sturm,vert_sturm,context0);
    if (is_undef(seqhoriz)){
      realroots=seqhoriz;
      return -2;
    }
    pgcd=seqhoriz[4];
    if (is_one(pgcd)){
      A=a01+cst_i*b0; B=a01+cst_i*b1;
      seqvert=csturm_segment_seq(P,A,B,realroots,eps,horiz_sturm,vert_sturm,context0);
      if (is_undef(seqvert)){
	realroots=seqvert;
	return -2;
      }
      pgcd=seqvert[4];
    }
    if (!is_one(pgcd)){
      complex_roots_split(P,pgcd,a0,b0,a1,b1,realroots,complexroots,eps);
      return n;
    }
    /*
      (a0,b1)  - (a01,b1)  -  (a1,b1)         seq3                seq3
         |     n4   |     n3    |        seq4   n4      seqvert    n3     seq2
      (a0,b01) - (a01,b01) -  (a1,b01)        seqhoriz           seqhoriz
         |     n1   |     n2    |        seq4   n1      seqvert   n2      seq2
      (a0,b0)  - (a01,b0)  -  (a1,b0)         seq1                seq1
     */
    int n1=complex_roots(P,a0,b0,a01,b01,seq1,seqvert,seqhoriz,seq4,realroots,complexroots,eps,horiz_sturm,vert_sturm),nadd;
    if (n1==-2)
      return -2;
    if (n1==n)
      return n;
    n1 += (nadd=complex_roots(P,a01,b0,a1,b01,seq1,seq2,seqhoriz,seqvert,realroots,complexroots,eps,horiz_sturm,vert_sturm));
    if (nadd==-2)
      return -2;
    if (n1==n)
      return n;
    n1 += (nadd=complex_roots(P,a01,b01,a1,b1,seqhoriz,seq2,seq3,seqvert,realroots,complexroots,eps,horiz_sturm,vert_sturm));
    if (nadd==-2)
      return -2;
    if (n1==n)
      return n;
    n1 += (nadd=complex_roots(P,a0,b01,a01,b1,seqhoriz,seqvert,seq3,seq4,realroots,complexroots,eps,horiz_sturm,vert_sturm));
    if (nadd==-2)
      return -2;
    return n;
  }

  // Find complex roots of P in a0,b0 -> a1,b1
  static void complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps){
    if (P.size()<2)
      return;
    vecteur Seq1,Seq2,Seq3,Seq4,horiz_sturm,vert_sturm;
    gen pgcd;
    if (!csturm_square_seq(P,a0,b0,a1,b1,pgcd,realroots,eps,Seq1,Seq2,Seq3,Seq4,horiz_sturm,vert_sturm,context0))
      complex_roots_split(P,pgcd,a0,b0,a1,b1,realroots,complexroots,eps);
    else
      complex_roots(P,a0,b0,a1,b1,Seq1,Seq2,Seq3,Seq4,realroots,complexroots,eps,horiz_sturm,vert_sturm);
  }

  // Find complex roots of P in a0,b0 -> a1,b1
  bool complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,gen & pgcd,vecteur & roots,double eps){
    vecteur realroots,complexroots;
    complex_roots(P,a0,b0,a1,b1,realroots,complexroots,eps);
    if (is_undef(realroots))
      return false;
    roots=mergevecteur(roots,mergevecteur(realroots,complexroots));
    return true;
  }

  vecteur crationalroot(polynome & p,bool complexe){
    vectpoly v;
    int i=1;
    polynome qrem;
    environment * env= new environment;
    env->complexe=complexe || !is_zero(im(p,context0));
    vecteur w;
    if (!do_linearfind(p,env,qrem,v,w,i))
      w.clear();
    delete env;
    p=qrem;
    return w;
  }

  vecteur keep_in_rectangle(const vecteur & croots,const gen A0,const gen & B0,const gen & A1,const gen & B1,bool embed,GIAC_CONTEXT){
    vecteur roots;
    const_iterateur it=croots.begin(),itend=croots.end();
    for (;it!=itend;++it){
      gen a=re(*it,contextptr),b=im(*it,contextptr);
      if (is_greater(a,A0,contextptr)&&is_greater(A1,a,contextptr)&&is_greater(b,B0,contextptr)&&is_greater(B1,b,contextptr))
	roots.push_back(embed?makevecteur(*it,1):*it);
    }
    return roots;
  }

  gen square_modulus(const gen & g,GIAC_CONTEXT){
    return g.squarenorm(contextptr);
  }

  // P is the polynomial, P1 derivative, v list of approx roots 
  // (initially should have at least n bits precision), 
  // epsn is the target number of bits precision int(std::log(eps)/std::log(2.)-.5);
  // epsg2surdeg2 is eps^2/degree(P)^2 as a gen, epsg is the target precision
  // v[i] is set by newton_improve to be at distance at most vradius[i] of a root
  // kmax is the maximal number of Newton iterations
  bool newton_improve(const vecteur & P,const vecteur & P1,bool Preal,vecteur & v,vecteur & vradius,int i,int kmax,int n,int epsn,const gen & epsg2surdeg2,const gen & epsg){
    gen r=v[i];
    bool nextconj=false;
    if (Preal && i+1<int(v.size()))
      nextconj=is_exactly_zero(r-conj(v[i+1],context0));
    if (r.type==_FRAC || is_cinteger(r))
      return true;
    // find nearest root from v
    gen deltar=plus_inf,delta;
    for (unsigned j=0;j<v.size();++j){
      if (int(j)==i)
	continue;
      gen tmp=abs(r-v[j],context0);
      if (is_strictly_greater(deltar,tmp,context0))
	deltar=tmp;
    }
    if (is_zero(deltar))
      return false;
    deltar=deltar/3;
    gen sumdr2=0;
    int N=n; // effective value of number of bits for computation
#ifdef HAVE_LIBMPFR
    if (r.type==_REAL && mpfr_get_prec(r._REALptr->inf)>N)
      N=mpfr_get_prec(r._REALptr->inf);
    if (r.type==_CPLX && r._CPLXptr->type==_REAL && mpfr_get_prec(r._CPLXptr->_REALptr->inf)>N)
      N=mpfr_get_prec(r._CPLXptr->_REALptr->inf);
#endif
#if 0 // def HAVE_LIBMPFI
    gen deuxN=pow(2,N,context0);
    gen rr,ri,dr,di;
    reim(r,rr,ri,context0);
    if (Preal && !nextconj)
      r=eval(gen(makevecteur(rr-plus_one/deuxN,rr+plus_one/deuxN),_INTERVAL__VECT),1,context0);
    else
      r=eval(gen(makevecteur(rr-plus_one/deuxN,rr+plus_one/deuxN),_INTERVAL__VECT),1,context0)+cst_i*eval(gen(makevecteur(ri-plus_one/deuxN,ri+plus_one/deuxN),_INTERVAL__VECT),1,context0);
    for (int k=0;k<kmax;++k){
      // check if root precision is ok
      // otherwise try to improve root precision with Newton method
      gen P1r=horner(P1,r,0,false);
      if (is_exactly_zero(P1r)){
	delta=plus_inf;
	break;
      }
      gen Pr=horner(P,r,0,false);
      delta=Pr/P1r;
      bool test;
      if (Preal && ! nextconj){
	test=contains(delta,dr);
	delta=abs(delta,context0);
      }
      else {
	reim(delta,dr,di,context0);
	test= contains(rr,dr) && contains(ri,di);
	delta=square_modulus(delta,context0);
      }
      if (test){
	v[i]=r; // we can certify there is a root in r by Brouwer fixed thm
	vradius[i]=-1;
	if (nextconj){
	  v[i+1]=conj(r,context0);
	  vradius[i+1]=vradius[i];
	  ++i;
	}
	break;
      }
      if (delta.type==_REAL){
	if (real_interval * ptr=dynamic_cast<real_interval *>(delta._REALptr)){
	  mpfr_t tmp; mpfr_init(tmp);
	  mpfi_get_right(tmp,ptr->infsup);
	  delta=real_object(tmp);
	  mpfr_clear(tmp);
	}
      }
      sumdr2 += delta;
      if (!is_greater(deltar*deltar,sumdr2,context0)){
	CERR << "Unable to certify " << v[i] << endl ;
	return false;
      }
      if (N<P.size()-epsn){
	// add 10 bits of precision or double it
	if (N<-epsn/2){
	  deuxN=deuxN*deuxN;
	  N*=2;
	}
	else {
	  deuxN=1024*deuxN;
	  N+=10;
	}
      }
      r -= Pr/P1r;
      // change precision to N
      reim(r,rr,ri,context0);
      if (rr.type==_REAL){
	if (real_interval * ptr=dynamic_cast<real_interval *>(rr._REALptr))
	  mpfi_set_prec(ptr->infsup,N);
      }
      if (ri.type==_REAL){
	if (real_interval * ptr=dynamic_cast<real_interval *>(ri._REALptr))
	  mpfi_set_prec(ptr->infsup,N);
      }
      r=rr+cst_i*ri;
    } // end for k
#else
    if (N>int(P.size())/4-epsn/2)
      N=int(P.size())/4-epsn/2;
    gen deuxN=pow(2,N,context0);
    for (int k=0;k<kmax;++k){
      in_round2(r,deuxN,N); // round2(r,deuxN,context0);
      // check if root precision is ok
      // otherwise try to improve root precision with Newton method
      gen P1r=horner(P1,r,0,false);
      if (is_exactly_zero(P1r)){
	delta=plus_inf;
	break;
      }
      gen Pr=horner(P,r,0,false);
      delta=square_modulus(Pr,context0)/square_modulus(P1r,context0);
      if (is_greater(epsg2surdeg2,delta,context0)){
	v[i]=r; // we can certify there is a root at distance <= eps from r
	if (is_exactly_zero(Pr))
	  vradius[i]=0;
	else
	  vradius[i]=n*sqrt(accurate_evalf(delta,100),context0);
	// problem with double underflow
	if (!is_exactly_zero(vradius[i]))
	  vradius[i]=min(epsg,pow(plus_two,int(evalf_double(ln(vradius[i],context0),1,context0)._DOUBLE_val/std::log(2.))+1),context0);
	if (debug_infolevel)
	  CERR << CLOCK() << " isolated " << r << " radius " << vradius[i] << endl;
	if (nextconj){
	  v[i+1]=conj(r,context0);
	  vradius[i+1]=vradius[i];
	  ++i;
	}
	break;
      }
      sumdr2 += delta;
      if (!is_greater(deltar*deltar,sumdr2,context0)){
	CERR << "Unable to certify " << v[i] << endl ;
	return false;
      }
      if (N<int(P.size())-epsn){
	// add 10 bits of precision or double it
	if (N<-epsn){
	  deuxN=deuxN*deuxN;
	  N*=2;
	}
	else {
	  deuxN=1024*deuxN;
	  N+=10;
	}
      }
      in_round2(Pr,deuxN,N); in_round2(P1r,deuxN,N); // round2(Pr,deuxN,context0); round2(P1r,deuxN,context0);
      r -= Pr/P1r;
    } // end for k
#endif
    if (!is_greater(epsg*epsg,delta,context0))
      return false;
    return true;
  }

  // find roots of polynomial P at precision eps using proot or 
  // complex Sturm sequences
  // P must have numeric coefficients, in Q[i]
  vecteur complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,bool complexe,double eps,bool use_proot){
    if (P.empty())
      return P;
    eps=std::abs(eps);
    if (eps>1e-6)
      eps=1e-6;
    if (use_proot){
      int epsn=int(std::log(eps)/std::log(2.)-.5);
      gen epsg=pow(2,epsn,context0);
      gen epsg2surdeg2=(epsg*epsg)/int((P.size()+1)*(P.size()+1));	
      // first try proot with double precision, if it's not enough increase
      int n=45;
      bool Preal=is_zero(im(P,context0));
      modpoly P1=derivative(P);
      for (;n<400;n*=2){
	double cureps=std::pow(2.0,-n);
	if (debug_infolevel)
	  CERR << CLOCK() << " proot at precision " << cureps << endl;
	vecteur v=proot(P,cureps,n);
	if (debug_infolevel)
	  CERR << CLOCK() << " proot end at precision " << cureps << endl;
	vecteur vradius(v.size());
	unsigned i=0;
	int kmax=int(std::log(eps)/std::log(cureps))+4;
	for (;i<v.size();++i){
	  newton_improve(P,P1,Preal,v,vradius,i,kmax,n,epsn,epsg2surdeg2,epsg);
	} // end for i
	if (i==v.size()){
	  vecteur res;
	  for (unsigned j=0;j<v.size();++j){
	    if (Preal && is_zero(im(v[j],context0))){
	      if (is_exactly_zero(vradius[j]) || vradius[j]==-1){
		if (is_greater(v[j],a0,context0) && is_greater(a1,v[j],context0) && is_greater(0,b0,context0) && is_greater(b1,0,context0))
		  res.push_back(makevecteur(v[j],1));
		continue;
	      }
	      gen P1=horner(P,v[j]-vradius[j],0,false),P2=horner(P,v[j]+vradius[j],0,false);
	      if (P1.type==_FRAC) P1=P1._FRACptr->num;
	      if (P2.type==_FRAC) P2=P2._FRACptr->num;
	      if (is_strictly_positive(-P1*P2,context0)){
		if (is_greater(v[j],a0,context0) && is_greater(a1,v[j],context0) && is_greater(0,b0,context0) && is_greater(b1,0,context0))
		  res.push_back(makevecteur(eval(change_subtype(makevecteur(v[j]-vradius[j],v[j]+vradius[j]),_INTERVAL__VECT),1,context0),1));
		continue;
	      }
	    }
	    gen R,I;
	    reim(v[j],R,I,context0);
	    if (is_greater(R,a0,context0) && is_greater(a1,R,context0) && is_greater(I,b0,context0) && is_greater(b1,I,context0)){
	      if (is_exactly_zero(vradius[j]))
		res.push_back(makevecteur(v[j],1));
	      else {
#ifdef HAVE_LIBMPFI
		gen a,b;
		reim(v[j],a,b,context0);
		res.push_back(makevecteur(eval(change_subtype(makevecteur(a-vradius[j],a+vradius[j]),_INTERVAL__VECT),1,context0)+cst_i*eval(change_subtype(makevecteur(b-vradius[j],b+vradius[j]),_INTERVAL__VECT),1,context0),1));
#else
		res.push_back(makevecteur(makevecteur(ratnormal(v[j]-vradius[j]*(1+cst_i)),ratnormal(v[j]+vradius[j]*(1+cst_i))),1));
#endif
	      }
	    }
	  }
	  return res;
	} // end if i==v.size()
      } // end for n
      CERR << "proot isolation did not work, trying complex Sturm sequences" << endl;
    }
    bool aplati=(a0==a1) && (b0==b1);
    if (!aplati && complexe && (a0==a1 || b0==b1) )
      return vecteur(1,gensizeerr(gettext("Square is flat!")));
    gen A0(a0),B0(b0),A1(a1),B1(b1);
    {
      // initial rectangle: |roots|< 1+ max(|a_i|)/|a_n|
      gen maxai=_max(*apply(P,abs,context0)._VECTptr,context0);
      gen tmp=1+maxai/abs(P.front(),context0);
      if (aplati){
	A0=-tmp;
	B0=-tmp;
	A1=tmp;
	B1=tmp;
      }
      if (is_inf(A0)) A0=-tmp;
      if (is_inf(B0)) B0=-tmp;
      if (is_inf(A1)) A1=tmp;
      if (is_inf(B1)) B1=tmp;
    }
    gen tmp;
    modpoly p(*apply(P,exact,context0)._VECTptr);
    lcmdeno(p,tmp,context0);
    polynome pp(poly12polynome(p));
    if (!complexe){
      gen tmp=gcd(re(pp,context0),im(pp,context0));
      if (tmp.type!=_POLY)
	return vecteur(0);
      pp=*tmp._POLYptr;
    }
    vecteur croots=crationalroot(pp,complexe);
    vecteur roots=keep_in_rectangle(croots,A0,B0,A1,B1,true,context0);
    p=polynome2poly1(pp);
    gen an=p.front();
    if (!is_zero(im(an,context0)))
      p=conj(p.front(),context0)*p;
    if (!complexe){ // real root isolation
      modpoly R=p;
      modpoly S=derivative(R);
      vecteur listquo,coeffP,coeffR;
      csturm_seq(S,R,listquo,coeffP,coeffR,context0);
      // sparse polynomial patch
      if (pp.coord.size()<p.size()/10.){
	R=vecteur(1,poly12polynome(R,1,1));
	S=vecteur(1,poly12polynome(S,1,1));
      }
      csturm_realroots(S,R,listquo,coeffP,coeffR,0,1,A0,A1,roots,eps,context0);
      return roots;
    }
    csturm_normalize(p,A0,B0,A1,B1,roots);
    gen pgcd;
    if (!complex_roots(p,A0,B0,A1,B1,pgcd,roots,eps))
      return vecteur(1,gensizeerr(context0));
    return roots;
  }


  gen complexroot(const gen & g,bool complexe,GIAC_CONTEXT){
    vecteur v=gen2vecteur(g);
    bool use_vas=!complexe ,use_proot=true;
#ifndef HAVE_LIBMPFR
    use_proot=false;
#endif
    bool isolation=false;
    if (!v.empty() && v[0]==at_sturm){
      use_vas=false;
      use_proot=false;
      v.erase(v.begin());
    }
    if (v.empty())
      return gensizeerr(contextptr);
    if (v.size()<2){
      isolation=true;
      v.push_back(epsilon(contextptr));
    }
    if (v.size()==3)
      v.insert(v.begin()+1,epsilon(contextptr));
    gen p=v.front(),prec=evalf_double(v[1],1,contextptr);
    if (prec.type!=_DOUBLE_)
      return gentypeerr(contextptr);
    double eps=prec._DOUBLE_val;
    if (eps>=1.0)
      eps=std::pow(10.,-eps);
    if (eps<=0)
      eps=epsilon(contextptr);
    if (eps<=0)
      eps=1e-12;
    if (v[0].type==_VECT && has_num_coeff(v[0])){
      v=proot(*v[0]._VECTptr,eps);
      vecteur w;
      for (unsigned i=0;i<v.size();++i){
	if (is_zero(im(v[i],contextptr)))
	  w.push_back(makevecteur(v[i],1));
      }
      return w;
    }
    unsigned vs=unsigned(v.size());
    gen A(0),B(0),a0(minus_inf),b0(minus_inf),a1(plus_inf),b1(plus_inf);
    if (vs>3){
      A=v[2];
      B=v[3];
      a0=re(A,contextptr); b0=im(A,contextptr);
      a1=re(B,contextptr);b1=im(B,contextptr);
    }
    if (is_greater(a0,a1,contextptr))
      swapgen(a0,a1);
    if (is_greater(b0,b1,contextptr))
      swapgen(b0,b1);
    vecteur vas_res;
    if (p.type==_VECT){
      if (use_vas && vas(*p._VECTptr,a0,a1,isolation?1e300:eps,vas_res,true,contextptr))
	return vas_res;
      return complex_roots(*p._VECTptr,a0,b0,a1,b1,complexe,eps,use_proot);
    }
    if (use_vas && vas(symb2poly_num(v[0],contextptr),a0,a1,isolation?1e300:eps,vas_res,true,contextptr))
      return vas_res;
    vecteur l,l0;
    lidnt(p,l0,false);
    if (l0.size()!=1)
      return gentypeerr(contextptr);
    l=alg_lvar(p);
    gen px=_e2r(makesequence(p,l),contextptr);
    if (px.type==_FRAC)
      px=px._FRACptr->num;
    if (px.type!=_POLY)
      return vecteur(0);
    factorization f(sqff(*px._POLYptr));
    factorization::const_iterator it=f.begin(),itend=f.end();
    vecteur res;
    for (;it!=itend;++it){
      gen P=_poly2symb(makesequence(it->fact,l),contextptr);
      P=_e2r(makesequence(P,l0.front()),contextptr);
      if (P.type!=_VECT)
	continue;
      vecteur tmp=complex_roots(*P._VECTptr,a0,b0,a1,b1,complexe,eps,use_proot);
      if (is_undef(tmp))
	return tmp;
      iterateur jt=tmp.begin(),jtend=tmp.end();
      for (;jt!=jtend;++jt){
	if (jt->type==_VECT && jt->_VECTptr->size()==2)
	  jt->_VECTptr->back()=it->mult*jt->_VECTptr->back();
      }
      res=mergevecteur(res,tmp);
    }
    return res;
  }

  gen _complexroot(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    gen res=complexroot(g,true,contextptr);
    if (res.type==_VECT)
      gen_sort_f_context(res._VECTptr->begin(),res._VECTptr->end(),complex_sort,contextptr);
    return res;
    // return _sorta(complexroot(g,true,contextptr),contextptr);
  }
  static const char _complexroot_s []="complexroot";
  static define_unary_function_eval (__complexroot,&giac::_complexroot,_complexroot_s);
  define_unary_function_ptr5( at_complexroot ,alias_at_complexroot,&__complexroot,0,true);

  gen _realroot(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    gen res; bool evalf_after=false;
    if (g.type==_VECT && !g._VECTptr->empty() && g._VECTptr->back()==at_evalf){
      res=complexroot(gen(vecteur(g._VECTptr->begin(),g._VECTptr->end()-1),g.subtype),false,contextptr);
      evalf_after=true;
    }
    else
      res=complexroot(g,false,contextptr);
    if (res.type!=_VECT)
      return res;
    vecteur v=*res._VECTptr;
    for (unsigned i=0;i<v.size();++i){
      if (v[i].type==_VECT && v[i]._VECTptr->size()==2){
	gen a=v[i]._VECTptr->front(),b=v[i]._VECTptr->back();
	if (a.type==_VECT && a.subtype==_INTERVAL__VECT){
	  if (evalf_after)
	    v[i]=evalf((a._VECTptr->front()+a._VECTptr->back())/2,1,contextptr);
	  else {
	    a=eval(a,1,contextptr);
	    v[i]=makevecteur(a,b);
	  }
	}
	else {
	  if (evalf_after)
	    v[i]=evalf(a,1,contextptr);
	}
      }
    }
    return v;
  }
  static const char _realroot_s []="realroot";
  static define_unary_function_eval (__realroot,&giac::_realroot,_realroot_s);
  define_unary_function_ptr5( at_realroot ,alias_at_realroot,&__realroot,0,true);

  static vecteur crationalroot(const gen & g0,bool complexe){
    gen g(g0),a,b;
    if (g.type==_VECT){
      if (g.subtype==_SEQ__VECT){
	vecteur & tmp=*g._VECTptr;
	if (tmp.size()!=3)
	  return vecteur(1,gendimerr(context0));
	g=tmp[0];
	a=tmp[1];
	b=tmp[2];
      }
      else {
	g=poly12polynome(*g._VECTptr);
      }
    }
    gen a0,b0,a1,b1;
    ab2a0b0a1b1(a,b,a0,b0,a1,b1,context0);
    vecteur l;
    lvar(g,l);
    if (l.empty())
      return vecteur(0);
    if (l.size()!=1)
      return vecteur(1,gentypeerr(context0));
    gen px=_e2r(makevecteur(g,l),context0);
    if (px.type==_FRAC)
      px=px._FRACptr->num;
    if (px.type!=_POLY)
      return vecteur(0);
    factorization f(sqff(*px._POLYptr));
    factorization::const_iterator it=f.begin(),itend=f.end();
    vecteur res;
    for (;it!=itend;++it){
      polynome p=it->fact;
      vecteur tmp=crationalroot(p,complexe);
      res=mergevecteur(res,tmp);
    }
    if (a0!=a1 || b0!=b1)
      res=keep_in_rectangle(res,a0,b0,a1,b1,false,context0);
    return res;
  }
  gen _crationalroot(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    return crationalroot(g,true);
  }
  static const char _crationalroot_s []="crationalroot";
  static define_unary_function_eval (__crationalroot,&giac::_crationalroot,_crationalroot_s);
  define_unary_function_ptr5( at_crationalroot ,alias_at_crationalroot,&__crationalroot,0,true);

  gen _rationalroot(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    return crationalroot(g,false);
  }
  static const char _rationalroot_s []="rationalroot";
  static define_unary_function_eval (__rationalroot,&giac::_rationalroot,_rationalroot_s);
  define_unary_function_ptr5( at_rationalroot ,alias_at_rationalroot,&__rationalroot,0,true);

  // convert numerator of g to a list
  vecteur symb2poly_num(const gen & g_,GIAC_CONTEXT){
    gen g(g_);
    if (g.type!=_VECT)
      g=makesequence(g,ggb_var(g));
    gen tmp=_symb2poly(g,contextptr);
    if (tmp.type==_FRAC)
      tmp=tmp._FRACptr->num;
    if (tmp.type!=_VECT)
      return vecteur(1,gensizeerr(contextptr));
    return *tmp._VECTptr;
  }
  // VAS implementation. Based on Xcas implementation by  Alkiviadis G. Akritas,
  // A first C++ implementation was written by Spyros Kehagias and others
  // but it was too close to the Xcas code 
  // This implementation is much faster, using basic data structures of giac
  // number of sign changes of the coefficients of P, returns -1 on error
  int variations(const modpoly & P,GIAC_CONTEXT){
    int res=0,n=int(P.size());
    if (!n)
      return -1;
    int previous=fastsign(P.front(),contextptr);
    if (previous==0)
      return -1;
    for (int i=1;i<n;i++){
      if (is_exactly_zero(P[i]))
	continue;
      int current=fastsign(P[i],contextptr);
      if (!current)
	return -1;
      if (current!=previous){
	++res;
	previous=current;
      }
    }
    return res;
  }

#ifndef M_LN2
#define M_LN2 0.6931471805599454
#endif

  // like (ln(n/d)+shift*ln2)/expo, but faster for large integers
  gen LMQ_evalf(const gen & n,const gen & d,double shift,int expo,GIAC_CONTEXT){
#ifndef USE_GMP_REPLACEMENTS
    if (is_integer(n) && is_integer(d)){
      long int nexp=0,dexp=0;
      double nmant,dmant;
      if (n.type==_INT_)
	nmant=n.val;
      else
	nmant=mpz_get_d_2exp (&nexp,*n._ZINTptr);
      if (d.type==_INT_)
	dmant=d.val;
      else
	dmant=mpz_get_d_2exp (&dexp,*d._ZINTptr);
      return ( std::log(-nmant/dmant) + (nexp-dexp+shift)*M_LN2 )/expo;
    }
#endif
    return ( ln(evalf(-n/d,1,contextptr),contextptr) + shift*M_LN2 )/gen(expo);
  }

  static bool compute_lnabsmantexpo(const vecteur & cl,vector<double> & cllnabsmant,vector<long int> & clexpo,vector<short int> & clsign,GIAC_CONTEXT){
    int k=int(cl.size());
    cllnabsmant.resize(k);
    clexpo.resize(k);
    clsign.resize(k);
    for (int i=0;i<k;++i){
      gen tmp=sign(cl[i],contextptr);
      if (tmp.type!=_INT_)
	return false;
      clsign[i]=tmp.val;
      double mant;
      long int expo=0;
      if (is_integer(cl[i])){
	if (cl[i].type==_ZINT){
#ifdef USE_GMP_REPLACEMENTS
	  mant=evalf_double(cl[i],1,contextptr)._DOUBLE_val;
#else
	  mant=mpz_get_d_2exp (&expo,*cl[i]._ZINTptr);
#endif
	}
	else mant=cl[i].val;
      }
      else {
	if (cl[i].type==_DOUBLE_){
	  mant=cl[i]._DOUBLE_val;
	}
	else {
	  mant=evalf_double(cl[i],1,contextptr)._DOUBLE_val;
	}
      }
      mant=std::log(std::abs(mant));
      cllnabsmant[i]=mant;
      clexpo[i]=expo;
    }
    return true;
  }


  gen posubLMQ(const modpoly & P,GIAC_CONTEXT){
    //---implements the Local_Max_Quadratic method (LMQ) to compute an
    //---upper bound on the values of the POSITIVE roots of p(x).
    
    //---Reference:"Linear and Quadratic Complexity Bounds on the Values of the 
    //---Positive Roots of Polynomials" by Alkiviadis G. Akritas. 
    //---Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009. 
    int k=int(P.size());
    if (k<=1)
      return 0;
    vecteur cl;
    if (is_strictly_positive(P.front(),contextptr))
      cl=P;
    else
      cl=-P;
    reverse(cl.begin(),cl.end());
    vector<double> cllnabsmant; 
    vector<long int> clexpo;
    vector<short int> clsign;
    if (!compute_lnabsmantexpo(cl,cllnabsmant,clexpo,clsign,contextptr))
      return gensizeerr(contextptr);
    gen tempmax=minus_inf;
    vector<int> timesused(k+1,1);
    for (int m=k-1;m>=1;--m){
      if (clsign[m-1]==-1){ // is_strictly_positive(-cl[m-1],contextptr)
	gen tempmin=plus_inf;
	for (int n=k;n>m;--n){
	  if (clsign[n-1]==1){ // is_strictly_positive(cl[n-1],contextptr)
	    gen temp= (cllnabsmant[m-1]-cllnabsmant[n-1] + (clexpo[m-1]-clexpo[n-1]+timesused[n-1])*M_LN2)/(n-m);// LMQ_evalf(cl[m-1],cl[n-1],timesused[n-1],n-m,contextptr);
	    // gen temp=pow(-cl[m-1]/cl[n-1]*pow(plus_two,timesused[n-1]),inv(n-m,contextptr),contextptr);
	    // temp=evalf(temp,1,contextptr);
	    ++timesused[n-1];
	    if (is_strictly_greater(tempmin,temp,contextptr))
	      tempmin=temp;
	  }
	}
	if (is_strictly_greater(tempmin,tempmax,contextptr))
	  tempmax=tempmin;
      }
    }
    return _ceil(65*exp(tempmax,contextptr)/64,contextptr); // small margin
  }

  gen _posubLMQ(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    vecteur v;
    if (g.type==_VECT && g.subtype!=_SEQ__VECT)
      v=*g._VECTptr;
    else
      v=symb2poly_num(g,contextptr);
    return posubLMQ(v,contextptr);
  }
  static const char _posubLMQ_s []="posubLMQ";
  static define_unary_function_eval (__posubLMQ,&giac::_posubLMQ,_posubLMQ_s);
  define_unary_function_ptr5( at_posubLMQ ,alias_at_posubLMQ,&__posubLMQ,0,true);

  gen poslbdLMQ(const modpoly & P,GIAC_CONTEXT){
    //---implements the Local_Max_Quadratic method (LMQ) to compute a
    //---lower bound on the values of the POSITIVE roots of p(x).
    
    //---Reference:"Linear and Quadratic Complexity Bounds on the Values of the 
    //---Positive Roots of Polynomials" by Alkiviadis G. Akritas. 
    //---Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009. 
    int k=int(P.size());
    if (k<=1)
      return 0;
    vecteur cl(P);
    reverse(cl.begin(),cl.end());
    if (is_strictly_positive(-cl.front(),contextptr))
      cl=-cl;
    vector<double> cllnabsmant; 
    vector<long int> clexpo;
    vector<short int> clsign;
    if (!compute_lnabsmantexpo(cl,cllnabsmant,clexpo,clsign,contextptr))
      return gensizeerr(contextptr);
    gen tempmax=minus_inf;
    vector<int> timesused(k,1);
    for (int m=1;m<k;++m){
      if (clsign[m]==-1){ // is_strictly_positive(-cl[m],contextptr)
	gen tempmin=plus_inf;
	for (int n=0;n<m;++n){
	  if (clsign[n]==1){ // is_strictly_positive(cl[n],contextptr)
	    // gen temp=pow(-cl[m]/cl[n]*pow(plus_two,timesused[n]),inv(m-n,contextptr),contextptr);
	    // temp=evalf(temp,1,contextptr);
	    gen temp= (cllnabsmant[m]-cllnabsmant[n] + (clexpo[m]-clexpo[n]+timesused[n])*M_LN2)/(m-n);// LMQ_evalf(cl[m],cl[n],timesused[n],m-n,contextptr);
	    ++timesused[n];
	    if (is_strictly_greater(tempmin,temp,contextptr))
	      tempmin=temp;
	  }
	}
	if (is_strictly_greater(tempmin,tempmax,contextptr))
	  tempmax=tempmin;
      }
    }
    tempmax=tempmax/M_LN2;
    tempmax=-_ceil(tempmax,contextptr);
    tempmax=pow(plus_two,tempmax,contextptr);
    return tempmax; 
  }
  
  gen _poslbdLMQ(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    vecteur v;
    if (g.type==_VECT && g.subtype!=_SEQ__VECT)
      v=*g._VECTptr;
    else
      v=symb2poly_num(g,contextptr);
    return poslbdLMQ(v,contextptr);
  }
  static const char _poslbdLMQ_s []="poslbdLMQ";
  static define_unary_function_eval (__poslbdLMQ,&giac::_poslbdLMQ,_poslbdLMQ_s);
  define_unary_function_ptr5( at_poslbdLMQ ,alias_at_poslbdLMQ,&__poslbdLMQ,0,true);

  vecteur makeinterval(const gen & a,const gen & b){
    if (is_strictly_greater(a,b,context0))
      return makevecteur(b,a);
    else
      return makevecteur(a,b);
  }

  bool vas_sort(const gen & a,const gen &b){
    gen a1(a),b1(b);
    if (a.type==_VECT && a._VECTptr->size()==2)
      a1=a._VECTptr->front();
    if (b.type==_VECT && b._VECTptr->size()==2)
      b1=b._VECTptr->front();
    return is_strictly_greater(b1,a1,context0);
  }

  // P is assumed to be squarefree and without rational roots
  // find roots of P((ax+b)/(cx+d))
  vecteur VAS_positive_roots(const modpoly & P,const gen & ap,const gen & bp,const gen & cp,const gen & dp,GIAC_CONTEXT){
    //---The steps below correspond to the steps described in the reference below.
    
    //---Reference:	"A Comparative Study of Two Real Root Isolation Methods" 
    //---by Alkiviadis G. Akritas and Adam W. Strzebonski. 
    //---Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005.
    vecteur res; // root isolation intervals
    vecteur intervals_to_process;
    // STEP 1
    int v0=variations(P,contextptr);
    if (!v0)
      return res;
    gen ub=posubLMQ(P,contextptr);
    if (v0==1)
      return vecteur(1,makeinterval(0,ap*ub));
    intervals_to_process.push_back(makevecteur(ap, bp, cp, dp, P,v0));

    // STEP 2
    while (!intervals_to_process.empty()){
      gen tmp=intervals_to_process.back();
      intervals_to_process.pop_back();
      if (tmp.type!=_VECT || tmp._VECTptr->size()!=6)
	return vecteur(1,gensizeerr("VAS interval"+tmp.print()));
      vecteur & tmpv=*tmp._VECTptr;
      gen a=tmpv[0],b=tmpv[1],c=tmpv[2],d=tmpv[3], genf=tmpv[4],genv=tmpv[5];
      if (genf.type!=_VECT || genv.type!=_INT_)
	return vecteur(1,gensizeerr("VAS interval"+tmp.print()));
      int v=genv.val;
      modpoly f = *genf._VECTptr;

      // STEP 3
      gen lb=poslbdLMQ(f,contextptr);

      // STEP 4
      if (is_strictly_greater(lb,16,contextptr)){
	change_scale(f,lb);
	a=lb*a; c=lb*c; lb=1;
      }
      
      // STEP 5
      if (is_greater(lb,1,contextptr)){
	// f=taylor(f,lb);
	change_scale(f,lb);
	f=taylor(f,1);
	back_change_scale(f,lb);
	b = lb*a + b; d = lb*c + d;
	if (is_zero(f.back())){
	  res.push_back(b/d);
	  f.pop_back();
	}
	v=variations(f,contextptr);
	if (!v)
	  continue;
	if (v==1){
	  if (!is_zero(c))
	    res.push_back(makeinterval(a/c,b/d));
	  else
	    res.push_back(makeinterval(b,b+a*posubLMQ(f,contextptr)));
	  continue;
	}
      }

      // STEP 6
      modpoly f1=taylor(f,1),f2;
      gen a1=a, b1=a+b, c1=c, d1=c+d;
      int r=0;
      if (is_zero(f1.back())){
	f1.pop_back();
	res.push_back(b1/d1);
	r=1;
      }
      int v1=variations(f1,contextptr);
      int v2=v-v1-r; 
      gen a2=b, b2=a+b, c2=d, d2=c+d;
      
      // STEP 7
      if (v2>1){
	f2=f;
	reverse(f2.begin(),f2.end());
	f2=taylor(f2,1);
	if (is_zero(f2.back()))
	  f2.pop_back();
	v2=variations(f2,contextptr);
      }

      // STEP 8
      if (v1<v2){
	swapgen(a1,a2);
	swapgen(b1,b2);
	swapgen(c1,c2);
	swapgen(d1,d2);
	swap(f1,f2);
	int i=v1; v1=v2; v2=i;
      }

      // STEP 9
      if (v1==0) continue;
      if (v1==1){
	if (!is_zero(c1))
	  res.push_back(makeinterval(a1/c1,b1/d1));
	else 
	  res.push_back(makeinterval(b1,b1 + a1*posubLMQ(f1,contextptr)));
      }
      else 
	intervals_to_process.push_back(makevecteur(a1,b1,c1,d1,f1,v1));

      // STEP 10
      if (v2==0) continue;
      if (v2==1){
	if (!is_zero(c2))
	  res.push_back(makeinterval(a2/c2,b2/d2));
	else 
	  res.push_back(makeinterval(b2,b2 + a2*posubLMQ(f2,contextptr)));
      }
      else 
	intervals_to_process.push_back(makevecteur(a2,b2,c2,d2,f2,v2));
    }
    gen_sort_f(res.begin(),res.end(),vas_sort);
    return res;
  }

  gen _VAS_positive(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    vecteur v;
    if (g.type==_VECT && g.subtype!=_SEQ__VECT)
      v=*g._VECTptr;
    else
      v=symb2poly_num(g,contextptr);
    return VAS_positive_roots(v,1,0,0,1,contextptr);
  }
  static const char _VAS_positive_s []="VAS_positive";
  static define_unary_function_eval (__VAS_positive,&giac::_VAS_positive,_VAS_positive_s);
  define_unary_function_ptr5( at_VAS_positive ,alias_at_VAS_positive,&__VAS_positive,0,true);

  // square-free factorization of p, then remove all exponents
  // optionally remove factors with even multiplicities
  modpoly remove_multiplicities(const modpoly & p,factorization & f,bool odd_only,GIAC_CONTEXT){
    vecteur res(1,1),tmp;
    polynome P;
    poly12polynome(p,1,P,1);
    P=P/lgcd(P);
    f=sqff(P);
    factorization::const_iterator it=f.begin(),itend=f.end();
    for (;it!=itend;++it){
      if (odd_only && it->mult%2==0)
	continue;
      polynome2poly1(it->fact,1,tmp);
      res=operator_times(res,tmp,0);
    }
    return res;
  }

  gen vas(const modpoly & p,GIAC_CONTEXT){
    vecteur v(p);
    vecteur res;
    bool has_zero=false;
    if (is_zero(v.back())){
      has_zero=true;
      v.pop_back();
    }
    res=VAS_positive_roots(v,1,0,0,1,contextptr);
    vecteur w(v);
    change_scale(w,-1);
    if (w.size()%2==0)
      w=-w;
    if (w==v){
      v=-res;
      reverse(v.begin(),v.end());
      iterateur it=v.begin(),itend=v.end();
      for (;it!=itend;++it){
	if (it->type==_VECT)
	  reverse(it->_VECTptr->begin(),it->_VECTptr->end());
      }
      if (has_zero)
	v.push_back(0);
      res=mergevecteur(v,res);
    }
    else {
      v=VAS_positive_roots(w,-1,0,0,1,contextptr);
      if (has_zero)
	v.push_back(0);
      res=mergevecteur(v,res);
    }
    return res;
  }
  gen _VAS(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    vecteur v;
    if (g.type==_VECT && g.subtype!=_SEQ__VECT)
      v=*g._VECTptr;
    else
      v=symb2poly_num(g,contextptr);
    factorization f;
    v=remove_multiplicities(v,f,false,contextptr);
    return vas(v,contextptr);
  }
  static const char _VAS_s []="VAS";
  static define_unary_function_eval (__VAS,&giac::_VAS,_VAS_s);
  define_unary_function_ptr5( at_VAS ,alias_at_VAS,&__VAS,0,true);

  static const double bisection_newton_eps=1e-3;

  // returns true if a root of p is found by Newton method, such that res-eps>a0
  // res+eps<b0 and sign of p changes between res-eps and res+eps
  static bool bisection_newton(const modpoly & P,const modpoly & Pprime,const gen & a0,const gen & a1,gen & x0,double eps,gen & eps2,GIAC_CONTEXT){
    if (is_greater(bisection_newton_eps,x0-a0,contextptr) || is_greater(bisection_newton_eps,a1-x0,contextptr))
      return false; // bisection is faster if the root is too close to the isolation interval boundaries
    int n=int(-std::log(eps)/std::log(2.0)+.5); // for rounding
    eps2=pow(2,-(n+1),contextptr);
    gen deuxn4(pow(2,n+4,contextptr));
    in_round2(x0,deuxn4,n+4); // round2(x0,deuxn4,contextptr);
    for (int ii=0;ii<n;ii++){
      gen Pprimex0=horner(Pprime,x0,0,false);
      if (is_zero(Pprimex0,contextptr))
	return false;
      gen Px0=horner(P,x0,0,false);
      in_round2(Px0,deuxn4,n+4); // round2(Px0,deuxn4,contextptr);
      in_round2(Pprimex0,deuxn4,n+4); // round2(Pprimex0,deuxn4,contextptr);
      gen dx=Px0/Pprimex0;
      in_round2(dx,deuxn4,n+4); // round2(dx,deuxn4,contextptr);
      x0=x0-dx;
      if (is_positive(a0-x0,contextptr) || is_positive(x0-a1,contextptr))
	return false;
      if (is_greater(abs(dx,contextptr),eps2,contextptr))
	continue;
      if (is_positive(-horner(P,x0-eps2,0,false)*horner(P,x0+eps2,0,false),contextptr))
	return true;
    }
    return false;
    
  }

  gen bisection(const modpoly & p,const gen & a0,const gen &b0,double eps,GIAC_CONTEXT){
    int nsteps=int(std::ceil(std::log(evalf_double(b0-a0,1,contextptr)._DOUBLE_val/eps)/M_LN2));
    int trynewtonstep=int(nsteps-std::log(bisection_newton_eps/eps)/M_LN2);
    gen a(a0),b(b0),m,eps2;
    modpoly dp=derivative(p);
    int s1=fastsign(horner(p,a,0,false),contextptr);
    if (s1==0)
      s1=fastsign(horner(dp,a,0,false),contextptr);
    int s2=fastsign(horner(p,b,0,false),contextptr);
    if (s2==0)
      s2=-fastsign(horner(dp,b,0,false),contextptr);
    if (s1==s2)
      return undef;
    for (int i=0;i<nsteps;++i){
      m=(a+b)/2;
      if (i==trynewtonstep && bisection_newton(p,dp,a0,b0,m,eps,eps2,contextptr))
	return makevecteur(m-eps2,m+eps2);
      int s=fastsign(horner(p,m,0,false),contextptr);
      if (s==0)
	return m;
      if (s==s1)
	a=m;
      else
	b=m;
    }
    return makevecteur(a,b);
  }

  int multiplicity(const factorization & f,const gen & interval,GIAC_CONTEXT){
    factorization::const_iterator it=f.begin(),itend=f.end();
    if (interval.type==_VECT && interval._VECTptr->size()==2){
      for (;it!=itend;++it){
	if (is_strictly_positive(-it->fact(interval._VECTptr->front())*it->fact(interval._VECTptr->back()),contextptr))
	  return it->mult;
      }
      for (it=f.begin();it!=itend;++it){
	if (is_positive(-it->fact(interval._VECTptr->front())*it->fact(interval._VECTptr->back()),contextptr))
	  return it->mult;
      }
    }
    else {
      for (;it!=itend;++it){
	if (is_zero(it->fact(interval)))
	  return it->mult;
      }
    }
    return 0;
  }

  static void add_vasres(vecteur & vasres,const gen & a,const gen & a0,const gen & b0,int mult,bool with_mult,GIAC_CONTEXT){
    if (a0==b0 || (is_greater(a,a0,contextptr) && is_greater(b0,a,contextptr)) )
      vasres.push_back(with_mult?gen(makevecteur(a,mult)):a);
  }

  // isolate and find real roots of P at precision eps between a and b
  // returns a list of intervals or of rationals
  bool vas(const modpoly & P,const gen & a0,const gen &b0,double eps,vecteur & vasres,bool with_mult,GIAC_CONTEXT){
    if (P.size()<=3){
      if (P.size()<2)
	return true;
      gen a(P[0]),b(P[1]);
      if (P.size()==2){
	a=-b/a;
	add_vasres(vasres,a,a0,b0,1,with_mult,contextptr);
	return true;
      }
      gen c(P[2]);
      gen delta=b*b-4*a*c;
      if (is_zero(delta)){
	a=-b/a/2;
	add_vasres(vasres,a,a0,b0,2,with_mult,contextptr);
	return true;	
      }
      if (is_positive(delta,contextptr)){
	delta=sqrt(delta,contextptr)/a/2;
	c=-b/a/2;
	add_vasres(vasres,c-delta,a0,b0,1,with_mult,contextptr);
	add_vasres(vasres,c+delta,a0,b0,1,with_mult,contextptr);
      }
      return true;
    }
    gen a(a0),b(b0);
    if (a==b){
      a=minus_inf;
      b=plus_inf;
    }
    // check and convert coeffs of P
    modpoly p(P);
    iterateur it=p.begin(),itend=p.end();
    for (;it!=itend;++it){
      *it=exact(*it,contextptr);
    }
    gen tmp;
    lcmdeno(p,tmp,contextptr);
    for (it=p.begin();it!=itend;++it){
      if (!is_integer(*it))
	return false;
    }
    p=divvecteur(p,lgcd(p));
    factorization f;
    p=remove_multiplicities(p,f,false,contextptr);
    tmp=vas(p,contextptr);
    if (tmp.type!=_VECT)
      return false;
    vecteur v=*tmp._VECTptr;
    // now improve precision by bisection 
    it=v.begin(),itend=v.end();
    for (;it!=itend;++it){
      if (it->type!=_VECT){
	if (is_greater(*it,a,contextptr) && is_greater(b,*it,contextptr)){
	  if (with_mult){
	    int n=multiplicity(f,*it,contextptr);
	    vasres.push_back(makevecteur(*it,n));
	  }
	  else
	    vasres.push_back(*it);
	}
	continue;
      }
      if (it->_VECTptr->size()!=2)
	return false;
      gen A=it->_VECTptr->front(),B=it->_VECTptr->back();
      if (is_strictly_greater(a,B,contextptr) || is_strictly_greater(A,b,contextptr))
	continue;
      gen interval=bisection(p,max(A,a,contextptr),min(B,b,contextptr),eps,contextptr);
      if (is_undef(interval))
	continue;
      if (interval.type==_VECT)
	interval.subtype=_INTERVAL__VECT;
      if (with_mult)
	vasres.push_back(makevecteur(interval,multiplicity(f,interval,contextptr)));
      else
	vasres.push_back( (interval.type==_VECT && interval._VECTptr->size()==2)?evalf((interval._VECTptr->front()+interval._VECTptr->back())/2,1,contextptr):interval);
    }
    return true;
  }

#ifndef NO_NAMESPACE_GIAC
} // namespace giac
#endif // ndef NO_NAMESPACE_GIAC